1 2 /* 3 Defines the basic matrix operations for sequential dense. 4 */ 5 6 #include <../src/mat/impls/dense/seq/dense.h> /*I "petscmat.h" I*/ 7 #include <petscblaslapack.h> 8 9 #include <../src/mat/impls/aij/seq/aij.h> 10 EXTERN_C_BEGIN 11 #undef __FUNCT__ 12 #define __FUNCT__ "MatConvert_SeqDense_SeqAIJ" 13 PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 14 { 15 Mat B; 16 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 17 PetscErrorCode ierr; 18 PetscInt i; 19 PetscInt *rows; 20 MatScalar *aa = a->v; 21 22 PetscFunctionBegin; 23 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 24 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 25 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 26 ierr = MatSeqAIJSetPreallocation(B,A->cmap->n,PETSC_NULL);CHKERRQ(ierr); 27 28 ierr = PetscMalloc(A->rmap->n*sizeof(PetscInt),&rows);CHKERRQ(ierr); 29 for (i=0; i<A->rmap->n; i++) rows[i] = i; 30 31 for (i=0; i<A->cmap->n; i++) { 32 ierr = MatSetValues(B,A->rmap->n,rows,1,&i,aa,INSERT_VALUES);CHKERRQ(ierr); 33 aa += a->lda; 34 } 35 ierr = PetscFree(rows);CHKERRQ(ierr); 36 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 37 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 38 39 if (reuse == MAT_REUSE_MATRIX) { 40 ierr = MatHeaderReplace(A,B);CHKERRQ(ierr); 41 } else { 42 *newmat = B; 43 } 44 PetscFunctionReturn(0); 45 } 46 EXTERN_C_END 47 48 #undef __FUNCT__ 49 #define __FUNCT__ "MatAXPY_SeqDense" 50 PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str) 51 { 52 Mat_SeqDense *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data; 53 PetscScalar oalpha = alpha; 54 PetscInt j; 55 PetscBLASInt N,m,ldax,lday,one = 1; 56 PetscErrorCode ierr; 57 58 PetscFunctionBegin; 59 N = PetscBLASIntCast(X->rmap->n*X->cmap->n); 60 m = PetscBLASIntCast(X->rmap->n); 61 ldax = PetscBLASIntCast(x->lda); 62 lday = PetscBLASIntCast(y->lda); 63 if (ldax>m || lday>m) { 64 for (j=0; j<X->cmap->n; j++) { 65 BLASaxpy_(&m,&oalpha,x->v+j*ldax,&one,y->v+j*lday,&one); 66 } 67 } else { 68 BLASaxpy_(&N,&oalpha,x->v,&one,y->v,&one); 69 } 70 ierr = PetscLogFlops(PetscMax(2*N-1,0));CHKERRQ(ierr); 71 PetscFunctionReturn(0); 72 } 73 74 #undef __FUNCT__ 75 #define __FUNCT__ "MatGetInfo_SeqDense" 76 PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info) 77 { 78 PetscInt N = A->rmap->n*A->cmap->n; 79 80 PetscFunctionBegin; 81 info->block_size = 1.0; 82 info->nz_allocated = (double)N; 83 info->nz_used = (double)N; 84 info->nz_unneeded = (double)0; 85 info->assemblies = (double)A->num_ass; 86 info->mallocs = 0; 87 info->memory = ((PetscObject)A)->mem; 88 info->fill_ratio_given = 0; 89 info->fill_ratio_needed = 0; 90 info->factor_mallocs = 0; 91 PetscFunctionReturn(0); 92 } 93 94 #undef __FUNCT__ 95 #define __FUNCT__ "MatScale_SeqDense" 96 PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha) 97 { 98 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 99 PetscScalar oalpha = alpha; 100 PetscErrorCode ierr; 101 PetscBLASInt one = 1,j,nz,lda = PetscBLASIntCast(a->lda); 102 103 PetscFunctionBegin; 104 if (lda>A->rmap->n) { 105 nz = PetscBLASIntCast(A->rmap->n); 106 for (j=0; j<A->cmap->n; j++) { 107 BLASscal_(&nz,&oalpha,a->v+j*lda,&one); 108 } 109 } else { 110 nz = PetscBLASIntCast(A->rmap->n*A->cmap->n); 111 BLASscal_(&nz,&oalpha,a->v,&one); 112 } 113 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 114 PetscFunctionReturn(0); 115 } 116 117 #undef __FUNCT__ 118 #define __FUNCT__ "MatIsHermitian_SeqDense" 119 PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool *fl) 120 { 121 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 122 PetscInt i,j,m = A->rmap->n,N; 123 PetscScalar *v = a->v; 124 125 PetscFunctionBegin; 126 *fl = PETSC_FALSE; 127 if (A->rmap->n != A->cmap->n) PetscFunctionReturn(0); 128 N = a->lda; 129 130 for (i=0; i<m; i++) { 131 for (j=i+1; j<m; j++) { 132 if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) PetscFunctionReturn(0); 133 } 134 } 135 *fl = PETSC_TRUE; 136 PetscFunctionReturn(0); 137 } 138 139 #undef __FUNCT__ 140 #define __FUNCT__ "MatDuplicateNoCreate_SeqDense" 141 PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues) 142 { 143 Mat_SeqDense *mat = (Mat_SeqDense*)A->data,*l; 144 PetscErrorCode ierr; 145 PetscInt lda = (PetscInt)mat->lda,j,m; 146 147 PetscFunctionBegin; 148 ierr = PetscLayoutReference(A->rmap,&newi->rmap);CHKERRQ(ierr); 149 ierr = PetscLayoutReference(A->cmap,&newi->cmap);CHKERRQ(ierr); 150 ierr = MatSeqDenseSetPreallocation(newi,PETSC_NULL);CHKERRQ(ierr); 151 if (cpvalues == MAT_COPY_VALUES) { 152 l = (Mat_SeqDense*)newi->data; 153 if (lda>A->rmap->n) { 154 m = A->rmap->n; 155 for (j=0; j<A->cmap->n; j++) { 156 ierr = PetscMemcpy(l->v+j*m,mat->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 157 } 158 } else { 159 ierr = PetscMemcpy(l->v,mat->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 160 } 161 } 162 newi->assembled = PETSC_TRUE; 163 PetscFunctionReturn(0); 164 } 165 166 #undef __FUNCT__ 167 #define __FUNCT__ "MatDuplicate_SeqDense" 168 PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat) 169 { 170 PetscErrorCode ierr; 171 172 PetscFunctionBegin; 173 ierr = MatCreate(((PetscObject)A)->comm,newmat);CHKERRQ(ierr); 174 ierr = MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 175 ierr = MatSetType(*newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 176 ierr = MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);CHKERRQ(ierr); 177 PetscFunctionReturn(0); 178 } 179 180 181 extern PetscErrorCode MatLUFactor_SeqDense(Mat,IS,IS,const MatFactorInfo*); 182 183 #undef __FUNCT__ 184 #define __FUNCT__ "MatLUFactorNumeric_SeqDense" 185 PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy) 186 { 187 MatFactorInfo info; 188 PetscErrorCode ierr; 189 190 PetscFunctionBegin; 191 ierr = MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);CHKERRQ(ierr); 192 ierr = MatLUFactor_SeqDense(fact,0,0,&info);CHKERRQ(ierr); 193 PetscFunctionReturn(0); 194 } 195 196 #undef __FUNCT__ 197 #define __FUNCT__ "MatSolve_SeqDense" 198 PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy) 199 { 200 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 201 PetscErrorCode ierr; 202 PetscScalar *x,*y; 203 PetscBLASInt one = 1,info,m = PetscBLASIntCast(A->rmap->n); 204 205 PetscFunctionBegin; 206 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 207 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 208 ierr = PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 209 if (A->factortype == MAT_FACTOR_LU) { 210 #if defined(PETSC_MISSING_LAPACK_GETRS) 211 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 212 #else 213 LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info); 214 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve"); 215 #endif 216 } else if (A->factortype == MAT_FACTOR_CHOLESKY){ 217 #if defined(PETSC_MISSING_LAPACK_POTRS) 218 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 219 #else 220 LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info); 221 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve"); 222 #endif 223 } 224 else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 225 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 226 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 227 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 228 PetscFunctionReturn(0); 229 } 230 231 #undef __FUNCT__ 232 #define __FUNCT__ "MatMatSolve_SeqDense" 233 PetscErrorCode MatMatSolve_SeqDense(Mat A,Mat B,Mat X) 234 { 235 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 236 PetscErrorCode ierr; 237 PetscScalar *b,*x; 238 PetscInt n; 239 PetscBLASInt nrhs,info,m=PetscBLASIntCast(A->rmap->n); 240 PetscBool flg; 241 242 PetscFunctionBegin; 243 ierr = PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,PETSC_NULL);CHKERRQ(ierr); 244 if (!flg) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix"); 245 ierr = PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,PETSC_NULL);CHKERRQ(ierr); 246 if (!flg) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix"); 247 248 ierr = MatGetSize(B,PETSC_NULL,&n);CHKERRQ(ierr); 249 nrhs = PetscBLASIntCast(n); 250 ierr = MatDenseGetArray(B,&b);CHKERRQ(ierr); 251 ierr = MatDenseGetArray(X,&x);CHKERRQ(ierr); 252 253 ierr = PetscMemcpy(x,b,m*nrhs*sizeof(PetscScalar));CHKERRQ(ierr); 254 255 if (A->factortype == MAT_FACTOR_LU) { 256 #if defined(PETSC_MISSING_LAPACK_GETRS) 257 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 258 #else 259 LAPACKgetrs_("N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info); 260 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve"); 261 #endif 262 } else if (A->factortype == MAT_FACTOR_CHOLESKY){ 263 #if defined(PETSC_MISSING_LAPACK_POTRS) 264 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 265 #else 266 LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info); 267 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve"); 268 #endif 269 } 270 else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 271 272 ierr = MatDenseRestoreArray(B,&b);CHKERRQ(ierr); 273 ierr = MatDenseRestoreArray(X,&x);CHKERRQ(ierr); 274 ierr = PetscLogFlops(nrhs*(2.0*m*m - m));CHKERRQ(ierr); 275 PetscFunctionReturn(0); 276 } 277 278 #undef __FUNCT__ 279 #define __FUNCT__ "MatSolveTranspose_SeqDense" 280 PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy) 281 { 282 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 283 PetscErrorCode ierr; 284 PetscScalar *x,*y; 285 PetscBLASInt one = 1,info,m = PetscBLASIntCast(A->rmap->n); 286 287 PetscFunctionBegin; 288 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 289 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 290 ierr = PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 291 /* assume if pivots exist then use LU; else Cholesky */ 292 if (mat->pivots) { 293 #if defined(PETSC_MISSING_LAPACK_GETRS) 294 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 295 #else 296 LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info); 297 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve"); 298 #endif 299 } else { 300 #if defined(PETSC_MISSING_LAPACK_POTRS) 301 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 302 #else 303 LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info); 304 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve"); 305 #endif 306 } 307 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 308 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 309 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 310 PetscFunctionReturn(0); 311 } 312 313 #undef __FUNCT__ 314 #define __FUNCT__ "MatSolveAdd_SeqDense" 315 PetscErrorCode MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 316 { 317 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 318 PetscErrorCode ierr; 319 PetscScalar *x,*y,sone = 1.0; 320 Vec tmp = 0; 321 PetscBLASInt one = 1,info,m = PetscBLASIntCast(A->rmap->n); 322 323 PetscFunctionBegin; 324 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 325 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 326 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 327 if (yy == zz) { 328 ierr = VecDuplicate(yy,&tmp);CHKERRQ(ierr); 329 ierr = PetscLogObjectParent(A,tmp);CHKERRQ(ierr); 330 ierr = VecCopy(yy,tmp);CHKERRQ(ierr); 331 } 332 ierr = PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 333 /* assume if pivots exist then use LU; else Cholesky */ 334 if (mat->pivots) { 335 #if defined(PETSC_MISSING_LAPACK_GETRS) 336 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 337 #else 338 LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info); 339 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve"); 340 #endif 341 } else { 342 #if defined(PETSC_MISSING_LAPACK_POTRS) 343 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 344 #else 345 LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info); 346 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve"); 347 #endif 348 } 349 if (tmp) { 350 ierr = VecAXPY(yy,sone,tmp);CHKERRQ(ierr); 351 ierr = VecDestroy(&tmp);CHKERRQ(ierr); 352 } else { 353 ierr = VecAXPY(yy,sone,zz);CHKERRQ(ierr); 354 } 355 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 356 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 357 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n);CHKERRQ(ierr); 358 PetscFunctionReturn(0); 359 } 360 361 #undef __FUNCT__ 362 #define __FUNCT__ "MatSolveTransposeAdd_SeqDense" 363 PetscErrorCode MatSolveTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 364 { 365 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 366 PetscErrorCode ierr; 367 PetscScalar *x,*y,sone = 1.0; 368 Vec tmp; 369 PetscBLASInt one = 1,info,m = PetscBLASIntCast(A->rmap->n); 370 371 PetscFunctionBegin; 372 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 373 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 374 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 375 if (yy == zz) { 376 ierr = VecDuplicate(yy,&tmp);CHKERRQ(ierr); 377 ierr = PetscLogObjectParent(A,tmp);CHKERRQ(ierr); 378 ierr = VecCopy(yy,tmp);CHKERRQ(ierr); 379 } 380 ierr = PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 381 /* assume if pivots exist then use LU; else Cholesky */ 382 if (mat->pivots) { 383 #if defined(PETSC_MISSING_LAPACK_GETRS) 384 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 385 #else 386 LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info); 387 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve"); 388 #endif 389 } else { 390 #if defined(PETSC_MISSING_LAPACK_POTRS) 391 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 392 #else 393 LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info); 394 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve"); 395 #endif 396 } 397 if (tmp) { 398 ierr = VecAXPY(yy,sone,tmp);CHKERRQ(ierr); 399 ierr = VecDestroy(&tmp);CHKERRQ(ierr); 400 } else { 401 ierr = VecAXPY(yy,sone,zz);CHKERRQ(ierr); 402 } 403 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 404 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 405 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n);CHKERRQ(ierr); 406 PetscFunctionReturn(0); 407 } 408 409 /* ---------------------------------------------------------------*/ 410 /* COMMENT: I have chosen to hide row permutation in the pivots, 411 rather than put it in the Mat->row slot.*/ 412 #undef __FUNCT__ 413 #define __FUNCT__ "MatLUFactor_SeqDense" 414 PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo) 415 { 416 #if defined(PETSC_MISSING_LAPACK_GETRF) 417 PetscFunctionBegin; 418 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable."); 419 #else 420 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 421 PetscErrorCode ierr; 422 PetscBLASInt n,m,info; 423 424 PetscFunctionBegin; 425 n = PetscBLASIntCast(A->cmap->n); 426 m = PetscBLASIntCast(A->rmap->n); 427 if (!mat->pivots) { 428 ierr = PetscMalloc((A->rmap->n+1)*sizeof(PetscBLASInt),&mat->pivots);CHKERRQ(ierr); 429 ierr = PetscLogObjectMemory(A,A->rmap->n*sizeof(PetscBLASInt));CHKERRQ(ierr); 430 } 431 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 432 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 433 LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info); 434 ierr = PetscFPTrapPop();CHKERRQ(ierr); 435 436 if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization"); 437 if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization"); 438 A->ops->solve = MatSolve_SeqDense; 439 A->ops->solvetranspose = MatSolveTranspose_SeqDense; 440 A->ops->solveadd = MatSolveAdd_SeqDense; 441 A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense; 442 A->factortype = MAT_FACTOR_LU; 443 444 ierr = PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);CHKERRQ(ierr); 445 #endif 446 PetscFunctionReturn(0); 447 } 448 449 #undef __FUNCT__ 450 #define __FUNCT__ "MatCholeskyFactor_SeqDense" 451 PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo) 452 { 453 #if defined(PETSC_MISSING_LAPACK_POTRF) 454 PetscFunctionBegin; 455 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable."); 456 #else 457 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 458 PetscErrorCode ierr; 459 PetscBLASInt info,n = PetscBLASIntCast(A->cmap->n); 460 461 PetscFunctionBegin; 462 ierr = PetscFree(mat->pivots);CHKERRQ(ierr); 463 464 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 465 LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info); 466 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1); 467 A->ops->solve = MatSolve_SeqDense; 468 A->ops->solvetranspose = MatSolveTranspose_SeqDense; 469 A->ops->solveadd = MatSolveAdd_SeqDense; 470 A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense; 471 A->factortype = MAT_FACTOR_CHOLESKY; 472 ierr = PetscLogFlops((A->cmap->n*A->cmap->n*A->cmap->n)/3.0);CHKERRQ(ierr); 473 #endif 474 PetscFunctionReturn(0); 475 } 476 477 478 #undef __FUNCT__ 479 #define __FUNCT__ "MatCholeskyFactorNumeric_SeqDense" 480 PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy) 481 { 482 PetscErrorCode ierr; 483 MatFactorInfo info; 484 485 PetscFunctionBegin; 486 info.fill = 1.0; 487 ierr = MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);CHKERRQ(ierr); 488 ierr = MatCholeskyFactor_SeqDense(fact,0,&info);CHKERRQ(ierr); 489 PetscFunctionReturn(0); 490 } 491 492 #undef __FUNCT__ 493 #define __FUNCT__ "MatCholeskyFactorSymbolic_SeqDense" 494 PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info) 495 { 496 PetscFunctionBegin; 497 fact->assembled = PETSC_TRUE; 498 fact->preallocated = PETSC_TRUE; 499 fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense; 500 PetscFunctionReturn(0); 501 } 502 503 #undef __FUNCT__ 504 #define __FUNCT__ "MatLUFactorSymbolic_SeqDense" 505 PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info) 506 { 507 PetscFunctionBegin; 508 fact->preallocated = PETSC_TRUE; 509 fact->assembled = PETSC_TRUE; 510 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense; 511 PetscFunctionReturn(0); 512 } 513 514 EXTERN_C_BEGIN 515 #undef __FUNCT__ 516 #define __FUNCT__ "MatGetFactor_seqdense_petsc" 517 PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact) 518 { 519 PetscErrorCode ierr; 520 521 PetscFunctionBegin; 522 ierr = MatCreate(((PetscObject)A)->comm,fact);CHKERRQ(ierr); 523 ierr = MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 524 ierr = MatSetType(*fact,((PetscObject)A)->type_name);CHKERRQ(ierr); 525 if (ftype == MAT_FACTOR_LU){ 526 (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense; 527 } else { 528 (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense; 529 } 530 (*fact)->factortype = ftype; 531 PetscFunctionReturn(0); 532 } 533 EXTERN_C_END 534 535 /* ------------------------------------------------------------------*/ 536 #undef __FUNCT__ 537 #define __FUNCT__ "MatSOR_SeqDense" 538 PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx) 539 { 540 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 541 PetscScalar *x,*b,*v = mat->v,zero = 0.0,xt; 542 PetscErrorCode ierr; 543 PetscInt m = A->rmap->n,i; 544 PetscBLASInt o = 1,bm = PetscBLASIntCast(m); 545 546 PetscFunctionBegin; 547 if (flag & SOR_ZERO_INITIAL_GUESS) { 548 /* this is a hack fix, should have another version without the second BLASdot */ 549 ierr = VecSet(xx,zero);CHKERRQ(ierr); 550 } 551 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 552 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 553 its = its*lits; 554 if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 555 while (its--) { 556 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 557 for (i=0; i<m; i++) { 558 xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o); 559 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 560 } 561 } 562 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 563 for (i=m-1; i>=0; i--) { 564 xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o); 565 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 566 } 567 } 568 } 569 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 570 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 571 PetscFunctionReturn(0); 572 } 573 574 /* -----------------------------------------------------------------*/ 575 #undef __FUNCT__ 576 #define __FUNCT__ "MatMultTranspose_SeqDense" 577 PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy) 578 { 579 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 580 PetscScalar *v = mat->v,*x,*y; 581 PetscErrorCode ierr; 582 PetscBLASInt m, n,_One=1; 583 PetscScalar _DOne=1.0,_DZero=0.0; 584 585 PetscFunctionBegin; 586 m = PetscBLASIntCast(A->rmap->n); 587 n = PetscBLASIntCast(A->cmap->n); 588 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 589 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 590 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 591 BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One); 592 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 593 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 594 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 595 PetscFunctionReturn(0); 596 } 597 598 #undef __FUNCT__ 599 #define __FUNCT__ "MatMult_SeqDense" 600 PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy) 601 { 602 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 603 PetscScalar *v = mat->v,*x,*y,_DOne=1.0,_DZero=0.0; 604 PetscErrorCode ierr; 605 PetscBLASInt m, n, _One=1; 606 607 PetscFunctionBegin; 608 m = PetscBLASIntCast(A->rmap->n); 609 n = PetscBLASIntCast(A->cmap->n); 610 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 611 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 612 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 613 BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One); 614 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 615 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 616 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);CHKERRQ(ierr); 617 PetscFunctionReturn(0); 618 } 619 620 #undef __FUNCT__ 621 #define __FUNCT__ "MatMultAdd_SeqDense" 622 PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 623 { 624 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 625 PetscScalar *v = mat->v,*x,*y,_DOne=1.0; 626 PetscErrorCode ierr; 627 PetscBLASInt m, n, _One=1; 628 629 PetscFunctionBegin; 630 m = PetscBLASIntCast(A->rmap->n); 631 n = PetscBLASIntCast(A->cmap->n); 632 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 633 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 634 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 635 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 636 BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One); 637 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 638 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 639 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n);CHKERRQ(ierr); 640 PetscFunctionReturn(0); 641 } 642 643 #undef __FUNCT__ 644 #define __FUNCT__ "MatMultTransposeAdd_SeqDense" 645 PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 646 { 647 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 648 PetscScalar *v = mat->v,*x,*y; 649 PetscErrorCode ierr; 650 PetscBLASInt m, n, _One=1; 651 PetscScalar _DOne=1.0; 652 653 PetscFunctionBegin; 654 m = PetscBLASIntCast(A->rmap->n); 655 n = PetscBLASIntCast(A->cmap->n); 656 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 657 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 658 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 659 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 660 BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One); 661 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 662 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 663 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n);CHKERRQ(ierr); 664 PetscFunctionReturn(0); 665 } 666 667 /* -----------------------------------------------------------------*/ 668 #undef __FUNCT__ 669 #define __FUNCT__ "MatGetRow_SeqDense" 670 PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 671 { 672 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 673 PetscScalar *v; 674 PetscErrorCode ierr; 675 PetscInt i; 676 677 PetscFunctionBegin; 678 *ncols = A->cmap->n; 679 if (cols) { 680 ierr = PetscMalloc((A->cmap->n+1)*sizeof(PetscInt),cols);CHKERRQ(ierr); 681 for (i=0; i<A->cmap->n; i++) (*cols)[i] = i; 682 } 683 if (vals) { 684 ierr = PetscMalloc((A->cmap->n+1)*sizeof(PetscScalar),vals);CHKERRQ(ierr); 685 v = mat->v + row; 686 for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;} 687 } 688 PetscFunctionReturn(0); 689 } 690 691 #undef __FUNCT__ 692 #define __FUNCT__ "MatRestoreRow_SeqDense" 693 PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 694 { 695 PetscErrorCode ierr; 696 PetscFunctionBegin; 697 if (cols) {ierr = PetscFree(*cols);CHKERRQ(ierr);} 698 if (vals) {ierr = PetscFree(*vals);CHKERRQ(ierr); } 699 PetscFunctionReturn(0); 700 } 701 /* ----------------------------------------------------------------*/ 702 #undef __FUNCT__ 703 #define __FUNCT__ "MatSetValues_SeqDense" 704 PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv) 705 { 706 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 707 PetscInt i,j,idx=0; 708 709 PetscFunctionBegin; 710 if (v) PetscValidScalarPointer(v,6); 711 if (!mat->roworiented) { 712 if (addv == INSERT_VALUES) { 713 for (j=0; j<n; j++) { 714 if (indexn[j] < 0) {idx += m; continue;} 715 #if defined(PETSC_USE_DEBUG) 716 if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 717 #endif 718 for (i=0; i<m; i++) { 719 if (indexm[i] < 0) {idx++; continue;} 720 #if defined(PETSC_USE_DEBUG) 721 if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 722 #endif 723 mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 724 } 725 } 726 } else { 727 for (j=0; j<n; j++) { 728 if (indexn[j] < 0) {idx += m; continue;} 729 #if defined(PETSC_USE_DEBUG) 730 if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 731 #endif 732 for (i=0; i<m; i++) { 733 if (indexm[i] < 0) {idx++; continue;} 734 #if defined(PETSC_USE_DEBUG) 735 if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 736 #endif 737 mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 738 } 739 } 740 } 741 } else { 742 if (addv == INSERT_VALUES) { 743 for (i=0; i<m; i++) { 744 if (indexm[i] < 0) { idx += n; continue;} 745 #if defined(PETSC_USE_DEBUG) 746 if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 747 #endif 748 for (j=0; j<n; j++) { 749 if (indexn[j] < 0) { idx++; continue;} 750 #if defined(PETSC_USE_DEBUG) 751 if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 752 #endif 753 mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 754 } 755 } 756 } else { 757 for (i=0; i<m; i++) { 758 if (indexm[i] < 0) { idx += n; continue;} 759 #if defined(PETSC_USE_DEBUG) 760 if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 761 #endif 762 for (j=0; j<n; j++) { 763 if (indexn[j] < 0) { idx++; continue;} 764 #if defined(PETSC_USE_DEBUG) 765 if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 766 #endif 767 mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 768 } 769 } 770 } 771 } 772 PetscFunctionReturn(0); 773 } 774 775 #undef __FUNCT__ 776 #define __FUNCT__ "MatGetValues_SeqDense" 777 PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[]) 778 { 779 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 780 PetscInt i,j; 781 782 PetscFunctionBegin; 783 /* row-oriented output */ 784 for (i=0; i<m; i++) { 785 if (indexm[i] < 0) {v += n;continue;} 786 if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D requested larger than number rows %D",indexm[i],A->rmap->n); 787 for (j=0; j<n; j++) { 788 if (indexn[j] < 0) {v++; continue;} 789 if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column %D requested larger than number columns %D",indexn[j],A->cmap->n); 790 *v++ = mat->v[indexn[j]*mat->lda + indexm[i]]; 791 } 792 } 793 PetscFunctionReturn(0); 794 } 795 796 /* -----------------------------------------------------------------*/ 797 798 #undef __FUNCT__ 799 #define __FUNCT__ "MatLoad_SeqDense" 800 PetscErrorCode MatLoad_SeqDense(Mat newmat,PetscViewer viewer) 801 { 802 Mat_SeqDense *a; 803 PetscErrorCode ierr; 804 PetscInt *scols,i,j,nz,header[4]; 805 int fd; 806 PetscMPIInt size; 807 PetscInt *rowlengths = 0,M,N,*cols,grows,gcols; 808 PetscScalar *vals,*svals,*v,*w; 809 MPI_Comm comm = ((PetscObject)viewer)->comm; 810 811 PetscFunctionBegin; 812 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 813 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor"); 814 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 815 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 816 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not matrix object"); 817 M = header[1]; N = header[2]; nz = header[3]; 818 819 /* set global size if not set already*/ 820 if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) { 821 ierr = MatSetSizes(newmat,M,N,M,N);CHKERRQ(ierr); 822 } else { 823 /* if sizes and type are already set, check if the vector global sizes are correct */ 824 ierr = MatGetSize(newmat,&grows,&gcols);CHKERRQ(ierr); 825 if (M != grows || N != gcols) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different length (%d, %d) than the input matrix (%d, %d)",M,N,grows,gcols); 826 } 827 ierr = MatSeqDenseSetPreallocation(newmat,PETSC_NULL);CHKERRQ(ierr); 828 829 if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */ 830 a = (Mat_SeqDense*)newmat->data; 831 v = a->v; 832 /* Allocate some temp space to read in the values and then flip them 833 from row major to column major */ 834 ierr = PetscMalloc((M*N > 0 ? M*N : 1)*sizeof(PetscScalar),&w);CHKERRQ(ierr); 835 /* read in nonzero values */ 836 ierr = PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);CHKERRQ(ierr); 837 /* now flip the values and store them in the matrix*/ 838 for (j=0; j<N; j++) { 839 for (i=0; i<M; i++) { 840 *v++ =w[i*N+j]; 841 } 842 } 843 ierr = PetscFree(w);CHKERRQ(ierr); 844 } else { 845 /* read row lengths */ 846 ierr = PetscMalloc((M+1)*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 847 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 848 849 a = (Mat_SeqDense*)newmat->data; 850 v = a->v; 851 852 /* read column indices and nonzeros */ 853 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),&scols);CHKERRQ(ierr); 854 cols = scols; 855 ierr = PetscBinaryRead(fd,cols,nz,PETSC_INT);CHKERRQ(ierr); 856 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&svals);CHKERRQ(ierr); 857 vals = svals; 858 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 859 860 /* insert into matrix */ 861 for (i=0; i<M; i++) { 862 for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j]; 863 svals += rowlengths[i]; scols += rowlengths[i]; 864 } 865 ierr = PetscFree(vals);CHKERRQ(ierr); 866 ierr = PetscFree(cols);CHKERRQ(ierr); 867 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 868 } 869 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 870 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 871 872 PetscFunctionReturn(0); 873 } 874 875 #undef __FUNCT__ 876 #define __FUNCT__ "MatView_SeqDense_ASCII" 877 static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer) 878 { 879 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 880 PetscErrorCode ierr; 881 PetscInt i,j; 882 const char *name; 883 PetscScalar *v; 884 PetscViewerFormat format; 885 #if defined(PETSC_USE_COMPLEX) 886 PetscBool allreal = PETSC_TRUE; 887 #endif 888 889 PetscFunctionBegin; 890 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 891 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 892 PetscFunctionReturn(0); /* do nothing for now */ 893 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 894 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 895 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer,"Matrix Object");CHKERRQ(ierr); 896 for (i=0; i<A->rmap->n; i++) { 897 v = a->v + i; 898 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 899 for (j=0; j<A->cmap->n; j++) { 900 #if defined(PETSC_USE_COMPLEX) 901 if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) { 902 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i) ",j,PetscRealPart(*v),PetscImaginaryPart(*v));CHKERRQ(ierr); 903 } else if (PetscRealPart(*v)) { 904 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",j,PetscRealPart(*v));CHKERRQ(ierr); 905 } 906 #else 907 if (*v) { 908 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",j,*v);CHKERRQ(ierr); 909 } 910 #endif 911 v += a->lda; 912 } 913 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 914 } 915 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 916 } else { 917 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 918 #if defined(PETSC_USE_COMPLEX) 919 /* determine if matrix has all real values */ 920 v = a->v; 921 for (i=0; i<A->rmap->n*A->cmap->n; i++) { 922 if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break ;} 923 } 924 #endif 925 if (format == PETSC_VIEWER_ASCII_MATLAB) { 926 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 927 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);CHKERRQ(ierr); 928 ierr = PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 929 ierr = PetscViewerASCIIPrintf(viewer,"%s = [\n",name);CHKERRQ(ierr); 930 } else { 931 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer,"Matrix Object");CHKERRQ(ierr); 932 } 933 934 for (i=0; i<A->rmap->n; i++) { 935 v = a->v + i; 936 for (j=0; j<A->cmap->n; j++) { 937 #if defined(PETSC_USE_COMPLEX) 938 if (allreal) { 939 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));CHKERRQ(ierr); 940 } else { 941 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));CHKERRQ(ierr); 942 } 943 #else 944 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);CHKERRQ(ierr); 945 #endif 946 v += a->lda; 947 } 948 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 949 } 950 if (format == PETSC_VIEWER_ASCII_MATLAB) { 951 ierr = PetscViewerASCIIPrintf(viewer,"];\n");CHKERRQ(ierr); 952 } 953 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 954 } 955 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 956 PetscFunctionReturn(0); 957 } 958 959 #undef __FUNCT__ 960 #define __FUNCT__ "MatView_SeqDense_Binary" 961 static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer) 962 { 963 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 964 PetscErrorCode ierr; 965 int fd; 966 PetscInt ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n; 967 PetscScalar *v,*anonz,*vals; 968 PetscViewerFormat format; 969 970 PetscFunctionBegin; 971 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 972 973 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 974 if (format == PETSC_VIEWER_NATIVE) { 975 /* store the matrix as a dense matrix */ 976 ierr = PetscMalloc(4*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 977 col_lens[0] = MAT_FILE_CLASSID; 978 col_lens[1] = m; 979 col_lens[2] = n; 980 col_lens[3] = MATRIX_BINARY_FORMAT_DENSE; 981 ierr = PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 982 ierr = PetscFree(col_lens);CHKERRQ(ierr); 983 984 /* write out matrix, by rows */ 985 ierr = PetscMalloc((m*n+1)*sizeof(PetscScalar),&vals);CHKERRQ(ierr); 986 v = a->v; 987 for (j=0; j<n; j++) { 988 for (i=0; i<m; i++) { 989 vals[j + i*n] = *v++; 990 } 991 } 992 ierr = PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 993 ierr = PetscFree(vals);CHKERRQ(ierr); 994 } else { 995 ierr = PetscMalloc((4+nz)*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 996 col_lens[0] = MAT_FILE_CLASSID; 997 col_lens[1] = m; 998 col_lens[2] = n; 999 col_lens[3] = nz; 1000 1001 /* store lengths of each row and write (including header) to file */ 1002 for (i=0; i<m; i++) col_lens[4+i] = n; 1003 ierr = PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 1004 1005 /* Possibly should write in smaller increments, not whole matrix at once? */ 1006 /* store column indices (zero start index) */ 1007 ict = 0; 1008 for (i=0; i<m; i++) { 1009 for (j=0; j<n; j++) col_lens[ict++] = j; 1010 } 1011 ierr = PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 1012 ierr = PetscFree(col_lens);CHKERRQ(ierr); 1013 1014 /* store nonzero values */ 1015 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&anonz);CHKERRQ(ierr); 1016 ict = 0; 1017 for (i=0; i<m; i++) { 1018 v = a->v + i; 1019 for (j=0; j<n; j++) { 1020 anonz[ict++] = *v; v += a->lda; 1021 } 1022 } 1023 ierr = PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 1024 ierr = PetscFree(anonz);CHKERRQ(ierr); 1025 } 1026 PetscFunctionReturn(0); 1027 } 1028 1029 #undef __FUNCT__ 1030 #define __FUNCT__ "MatView_SeqDense_Draw_Zoom" 1031 PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa) 1032 { 1033 Mat A = (Mat) Aa; 1034 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1035 PetscErrorCode ierr; 1036 PetscInt m = A->rmap->n,n = A->cmap->n,color,i,j; 1037 PetscScalar *v = a->v; 1038 PetscViewer viewer; 1039 PetscDraw popup; 1040 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,scale,maxv = 0.0; 1041 PetscViewerFormat format; 1042 1043 PetscFunctionBegin; 1044 1045 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 1046 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1047 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 1048 1049 /* Loop over matrix elements drawing boxes */ 1050 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 1051 /* Blue for negative and Red for positive */ 1052 color = PETSC_DRAW_BLUE; 1053 for(j = 0; j < n; j++) { 1054 x_l = j; 1055 x_r = x_l + 1.0; 1056 for(i = 0; i < m; i++) { 1057 y_l = m - i - 1.0; 1058 y_r = y_l + 1.0; 1059 #if defined(PETSC_USE_COMPLEX) 1060 if (PetscRealPart(v[j*m+i]) > 0.) { 1061 color = PETSC_DRAW_RED; 1062 } else if (PetscRealPart(v[j*m+i]) < 0.) { 1063 color = PETSC_DRAW_BLUE; 1064 } else { 1065 continue; 1066 } 1067 #else 1068 if (v[j*m+i] > 0.) { 1069 color = PETSC_DRAW_RED; 1070 } else if (v[j*m+i] < 0.) { 1071 color = PETSC_DRAW_BLUE; 1072 } else { 1073 continue; 1074 } 1075 #endif 1076 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1077 } 1078 } 1079 } else { 1080 /* use contour shading to indicate magnitude of values */ 1081 /* first determine max of all nonzero values */ 1082 for(i = 0; i < m*n; i++) { 1083 if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]); 1084 } 1085 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 1086 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 1087 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 1088 for(j = 0; j < n; j++) { 1089 x_l = j; 1090 x_r = x_l + 1.0; 1091 for(i = 0; i < m; i++) { 1092 y_l = m - i - 1.0; 1093 y_r = y_l + 1.0; 1094 color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(v[j*m+i])); 1095 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1096 } 1097 } 1098 } 1099 PetscFunctionReturn(0); 1100 } 1101 1102 #undef __FUNCT__ 1103 #define __FUNCT__ "MatView_SeqDense_Draw" 1104 PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer) 1105 { 1106 PetscDraw draw; 1107 PetscBool isnull; 1108 PetscReal xr,yr,xl,yl,h,w; 1109 PetscErrorCode ierr; 1110 1111 PetscFunctionBegin; 1112 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 1113 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 1114 if (isnull) PetscFunctionReturn(0); 1115 1116 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 1117 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 1118 xr += w; yr += h; xl = -w; yl = -h; 1119 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 1120 ierr = PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);CHKERRQ(ierr); 1121 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 1122 PetscFunctionReturn(0); 1123 } 1124 1125 #undef __FUNCT__ 1126 #define __FUNCT__ "MatView_SeqDense" 1127 PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer) 1128 { 1129 PetscErrorCode ierr; 1130 PetscBool iascii,isbinary,isdraw; 1131 1132 PetscFunctionBegin; 1133 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1134 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1135 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 1136 1137 if (iascii) { 1138 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 1139 } else if (isbinary) { 1140 ierr = MatView_SeqDense_Binary(A,viewer);CHKERRQ(ierr); 1141 } else if (isdraw) { 1142 ierr = MatView_SeqDense_Draw(A,viewer);CHKERRQ(ierr); 1143 } else { 1144 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported by dense matrix",((PetscObject)viewer)->type_name); 1145 } 1146 PetscFunctionReturn(0); 1147 } 1148 1149 #undef __FUNCT__ 1150 #define __FUNCT__ "MatDestroy_SeqDense" 1151 PetscErrorCode MatDestroy_SeqDense(Mat mat) 1152 { 1153 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 1154 PetscErrorCode ierr; 1155 1156 PetscFunctionBegin; 1157 #if defined(PETSC_USE_LOG) 1158 PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n); 1159 #endif 1160 ierr = PetscFree(l->pivots);CHKERRQ(ierr); 1161 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 1162 ierr = PetscFree(mat->data);CHKERRQ(ierr); 1163 1164 ierr = PetscObjectChangeTypeName((PetscObject)mat,0);CHKERRQ(ierr); 1165 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatDenseGetArray_C","",PETSC_NULL);CHKERRQ(ierr); 1166 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatDenseRestoreArray_C","",PETSC_NULL);CHKERRQ(ierr); 1167 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatSeqDenseSetPreallocation_C","",PETSC_NULL);CHKERRQ(ierr); 1168 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMult_seqaij_seqdense_C","",PETSC_NULL);CHKERRQ(ierr); 1169 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C","",PETSC_NULL);CHKERRQ(ierr); 1170 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C","",PETSC_NULL);CHKERRQ(ierr); 1171 PetscFunctionReturn(0); 1172 } 1173 1174 #undef __FUNCT__ 1175 #define __FUNCT__ "MatTranspose_SeqDense" 1176 PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout) 1177 { 1178 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1179 PetscErrorCode ierr; 1180 PetscInt k,j,m,n,M; 1181 PetscScalar *v,tmp; 1182 1183 PetscFunctionBegin; 1184 v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n; 1185 if (reuse == MAT_REUSE_MATRIX && *matout == A) { /* in place transpose */ 1186 if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can not transpose non-square matrix in place"); 1187 else { 1188 for (j=0; j<m; j++) { 1189 for (k=0; k<j; k++) { 1190 tmp = v[j + k*M]; 1191 v[j + k*M] = v[k + j*M]; 1192 v[k + j*M] = tmp; 1193 } 1194 } 1195 } 1196 } else { /* out-of-place transpose */ 1197 Mat tmat; 1198 Mat_SeqDense *tmatd; 1199 PetscScalar *v2; 1200 1201 if (reuse == MAT_INITIAL_MATRIX) { 1202 ierr = MatCreate(((PetscObject)A)->comm,&tmat);CHKERRQ(ierr); 1203 ierr = MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);CHKERRQ(ierr); 1204 ierr = MatSetType(tmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1205 ierr = MatSeqDenseSetPreallocation(tmat,PETSC_NULL);CHKERRQ(ierr); 1206 } else { 1207 tmat = *matout; 1208 } 1209 tmatd = (Mat_SeqDense*)tmat->data; 1210 v = mat->v; v2 = tmatd->v; 1211 for (j=0; j<n; j++) { 1212 for (k=0; k<m; k++) v2[j + k*n] = v[k + j*M]; 1213 } 1214 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1215 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1216 *matout = tmat; 1217 } 1218 PetscFunctionReturn(0); 1219 } 1220 1221 #undef __FUNCT__ 1222 #define __FUNCT__ "MatEqual_SeqDense" 1223 PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg) 1224 { 1225 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1226 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1227 PetscInt i,j; 1228 PetscScalar *v1 = mat1->v,*v2 = mat2->v; 1229 1230 PetscFunctionBegin; 1231 if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1232 if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1233 for (i=0; i<A1->rmap->n; i++) { 1234 v1 = mat1->v+i; v2 = mat2->v+i; 1235 for (j=0; j<A1->cmap->n; j++) { 1236 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1237 v1 += mat1->lda; v2 += mat2->lda; 1238 } 1239 } 1240 *flg = PETSC_TRUE; 1241 PetscFunctionReturn(0); 1242 } 1243 1244 #undef __FUNCT__ 1245 #define __FUNCT__ "MatGetDiagonal_SeqDense" 1246 PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1247 { 1248 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1249 PetscErrorCode ierr; 1250 PetscInt i,n,len; 1251 PetscScalar *x,zero = 0.0; 1252 1253 PetscFunctionBegin; 1254 ierr = VecSet(v,zero);CHKERRQ(ierr); 1255 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1256 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1257 len = PetscMin(A->rmap->n,A->cmap->n); 1258 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1259 for (i=0; i<len; i++) { 1260 x[i] = mat->v[i*mat->lda + i]; 1261 } 1262 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1263 PetscFunctionReturn(0); 1264 } 1265 1266 #undef __FUNCT__ 1267 #define __FUNCT__ "MatDiagonalScale_SeqDense" 1268 PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1269 { 1270 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1271 PetscScalar *l,*r,x,*v; 1272 PetscErrorCode ierr; 1273 PetscInt i,j,m = A->rmap->n,n = A->cmap->n; 1274 1275 PetscFunctionBegin; 1276 if (ll) { 1277 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1278 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1279 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1280 for (i=0; i<m; i++) { 1281 x = l[i]; 1282 v = mat->v + i; 1283 for (j=0; j<n; j++) { (*v) *= x; v+= m;} 1284 } 1285 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1286 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 1287 } 1288 if (rr) { 1289 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1290 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1291 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1292 for (i=0; i<n; i++) { 1293 x = r[i]; 1294 v = mat->v + i*m; 1295 for (j=0; j<m; j++) { (*v++) *= x;} 1296 } 1297 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1298 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 1299 } 1300 PetscFunctionReturn(0); 1301 } 1302 1303 #undef __FUNCT__ 1304 #define __FUNCT__ "MatNorm_SeqDense" 1305 PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1306 { 1307 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1308 PetscScalar *v = mat->v; 1309 PetscReal sum = 0.0; 1310 PetscInt lda=mat->lda,m=A->rmap->n,i,j; 1311 PetscErrorCode ierr; 1312 1313 PetscFunctionBegin; 1314 if (type == NORM_FROBENIUS) { 1315 if (lda>m) { 1316 for (j=0; j<A->cmap->n; j++) { 1317 v = mat->v+j*lda; 1318 for (i=0; i<m; i++) { 1319 #if defined(PETSC_USE_COMPLEX) 1320 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1321 #else 1322 sum += (*v)*(*v); v++; 1323 #endif 1324 } 1325 } 1326 } else { 1327 for (i=0; i<A->cmap->n*A->rmap->n; i++) { 1328 #if defined(PETSC_USE_COMPLEX) 1329 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1330 #else 1331 sum += (*v)*(*v); v++; 1332 #endif 1333 } 1334 } 1335 *nrm = PetscSqrtReal(sum); 1336 ierr = PetscLogFlops(2.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1337 } else if (type == NORM_1) { 1338 *nrm = 0.0; 1339 for (j=0; j<A->cmap->n; j++) { 1340 v = mat->v + j*mat->lda; 1341 sum = 0.0; 1342 for (i=0; i<A->rmap->n; i++) { 1343 sum += PetscAbsScalar(*v); v++; 1344 } 1345 if (sum > *nrm) *nrm = sum; 1346 } 1347 ierr = PetscLogFlops(A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1348 } else if (type == NORM_INFINITY) { 1349 *nrm = 0.0; 1350 for (j=0; j<A->rmap->n; j++) { 1351 v = mat->v + j; 1352 sum = 0.0; 1353 for (i=0; i<A->cmap->n; i++) { 1354 sum += PetscAbsScalar(*v); v += mat->lda; 1355 } 1356 if (sum > *nrm) *nrm = sum; 1357 } 1358 ierr = PetscLogFlops(A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1359 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 1360 PetscFunctionReturn(0); 1361 } 1362 1363 #undef __FUNCT__ 1364 #define __FUNCT__ "MatSetOption_SeqDense" 1365 PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg) 1366 { 1367 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1368 PetscErrorCode ierr; 1369 1370 PetscFunctionBegin; 1371 switch (op) { 1372 case MAT_ROW_ORIENTED: 1373 aij->roworiented = flg; 1374 break; 1375 case MAT_NEW_NONZERO_LOCATIONS: 1376 case MAT_NEW_NONZERO_LOCATION_ERR: 1377 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1378 case MAT_NEW_DIAGONALS: 1379 case MAT_KEEP_NONZERO_PATTERN: 1380 case MAT_IGNORE_OFF_PROC_ENTRIES: 1381 case MAT_USE_HASH_TABLE: 1382 case MAT_SYMMETRIC: 1383 case MAT_STRUCTURALLY_SYMMETRIC: 1384 case MAT_HERMITIAN: 1385 case MAT_SYMMETRY_ETERNAL: 1386 case MAT_IGNORE_LOWER_TRIANGULAR: 1387 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1388 break; 1389 default: 1390 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 1391 } 1392 PetscFunctionReturn(0); 1393 } 1394 1395 #undef __FUNCT__ 1396 #define __FUNCT__ "MatZeroEntries_SeqDense" 1397 PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1398 { 1399 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1400 PetscErrorCode ierr; 1401 PetscInt lda=l->lda,m=A->rmap->n,j; 1402 1403 PetscFunctionBegin; 1404 if (lda>m) { 1405 for (j=0; j<A->cmap->n; j++) { 1406 ierr = PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 1407 } 1408 } else { 1409 ierr = PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1410 } 1411 PetscFunctionReturn(0); 1412 } 1413 1414 #undef __FUNCT__ 1415 #define __FUNCT__ "MatZeroRows_SeqDense" 1416 PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1417 { 1418 PetscErrorCode ierr; 1419 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1420 PetscInt m = l->lda, n = A->cmap->n, i,j; 1421 PetscScalar *slot,*bb; 1422 const PetscScalar *xx; 1423 1424 PetscFunctionBegin; 1425 #if defined(PETSC_USE_DEBUG) 1426 for (i=0; i<N; i++) { 1427 if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 1428 if (rows[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D requested to be zeroed greater than or equal number of rows %D",rows[i],A->rmap->n); 1429 } 1430 #endif 1431 1432 /* fix right hand side if needed */ 1433 if (x && b) { 1434 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1435 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1436 for (i=0; i<N; i++) { 1437 bb[rows[i]] = diag*xx[rows[i]]; 1438 } 1439 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1440 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1441 } 1442 1443 for (i=0; i<N; i++) { 1444 slot = l->v + rows[i]; 1445 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 1446 } 1447 if (diag != 0.0) { 1448 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 1449 for (i=0; i<N; i++) { 1450 slot = l->v + (m+1)*rows[i]; 1451 *slot = diag; 1452 } 1453 } 1454 PetscFunctionReturn(0); 1455 } 1456 1457 #undef __FUNCT__ 1458 #define __FUNCT__ "MatDenseGetArray_SeqDense" 1459 PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[]) 1460 { 1461 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1462 1463 PetscFunctionBegin; 1464 if (mat->lda != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get array for Dense matrices with LDA different from number of rows"); 1465 *array = mat->v; 1466 PetscFunctionReturn(0); 1467 } 1468 1469 #undef __FUNCT__ 1470 #define __FUNCT__ "MatDenseRestoreArray_SeqDense" 1471 PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1472 { 1473 PetscFunctionBegin; 1474 *array = 0; /* user cannot accidently use the array later */ 1475 PetscFunctionReturn(0); 1476 } 1477 1478 #undef __FUNCT__ 1479 #define __FUNCT__ "MatDenseGetArray" 1480 /*@C 1481 MatDenseGetArray - gives access to the array where the data for a SeqDense matrix is stored 1482 1483 Not Collective 1484 1485 Input Parameter: 1486 . mat - a MATSEQDENSE matrix 1487 1488 Output Parameter: 1489 . array - pointer to the data 1490 1491 Level: intermediate 1492 1493 .seealso: MatDenseRestoreArray() 1494 @*/ 1495 PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array) 1496 { 1497 PetscErrorCode ierr; 1498 1499 PetscFunctionBegin; 1500 ierr = PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1501 PetscFunctionReturn(0); 1502 } 1503 1504 #undef __FUNCT__ 1505 #define __FUNCT__ "MatDenseRestoreArray" 1506 /*@C 1507 MatDenseRestoreArray - returns access to the array where the data for a SeqDense matrix is stored obtained by MatDenseGetArray() 1508 1509 Not Collective 1510 1511 Input Parameters: 1512 . mat - a MATSEQDENSE matrix 1513 . array - pointer to the data 1514 1515 Level: intermediate 1516 1517 .seealso: MatDenseGetArray() 1518 @*/ 1519 PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array) 1520 { 1521 PetscErrorCode ierr; 1522 1523 PetscFunctionBegin; 1524 ierr = PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1525 PetscFunctionReturn(0); 1526 } 1527 1528 #undef __FUNCT__ 1529 #define __FUNCT__ "MatGetSubMatrix_SeqDense" 1530 static PetscErrorCode MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 1531 { 1532 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1533 PetscErrorCode ierr; 1534 PetscInt i,j,nrows,ncols; 1535 const PetscInt *irow,*icol; 1536 PetscScalar *av,*bv,*v = mat->v; 1537 Mat newmat; 1538 1539 PetscFunctionBegin; 1540 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1541 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1542 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1543 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1544 1545 /* Check submatrixcall */ 1546 if (scall == MAT_REUSE_MATRIX) { 1547 PetscInt n_cols,n_rows; 1548 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1549 if (n_rows != nrows || n_cols != ncols) { 1550 /* resize the result matrix to match number of requested rows/columns */ 1551 ierr = MatSetSizes(*B,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1552 } 1553 newmat = *B; 1554 } else { 1555 /* Create and fill new matrix */ 1556 ierr = MatCreate(((PetscObject)A)->comm,&newmat);CHKERRQ(ierr); 1557 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1558 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1559 ierr = MatSeqDenseSetPreallocation(newmat,PETSC_NULL);CHKERRQ(ierr); 1560 } 1561 1562 /* Now extract the data pointers and do the copy,column at a time */ 1563 bv = ((Mat_SeqDense*)newmat->data)->v; 1564 1565 for (i=0; i<ncols; i++) { 1566 av = v + mat->lda*icol[i]; 1567 for (j=0; j<nrows; j++) { 1568 *bv++ = av[irow[j]]; 1569 } 1570 } 1571 1572 /* Assemble the matrices so that the correct flags are set */ 1573 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1574 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1575 1576 /* Free work space */ 1577 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1578 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1579 *B = newmat; 1580 PetscFunctionReturn(0); 1581 } 1582 1583 #undef __FUNCT__ 1584 #define __FUNCT__ "MatGetSubMatrices_SeqDense" 1585 PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1586 { 1587 PetscErrorCode ierr; 1588 PetscInt i; 1589 1590 PetscFunctionBegin; 1591 if (scall == MAT_INITIAL_MATRIX) { 1592 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1593 } 1594 1595 for (i=0; i<n; i++) { 1596 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1597 } 1598 PetscFunctionReturn(0); 1599 } 1600 1601 #undef __FUNCT__ 1602 #define __FUNCT__ "MatAssemblyBegin_SeqDense" 1603 PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 1604 { 1605 PetscFunctionBegin; 1606 PetscFunctionReturn(0); 1607 } 1608 1609 #undef __FUNCT__ 1610 #define __FUNCT__ "MatAssemblyEnd_SeqDense" 1611 PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 1612 { 1613 PetscFunctionBegin; 1614 PetscFunctionReturn(0); 1615 } 1616 1617 #undef __FUNCT__ 1618 #define __FUNCT__ "MatCopy_SeqDense" 1619 PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1620 { 1621 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense *)B->data; 1622 PetscErrorCode ierr; 1623 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 1624 1625 PetscFunctionBegin; 1626 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 1627 if (A->ops->copy != B->ops->copy) { 1628 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1629 PetscFunctionReturn(0); 1630 } 1631 if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1632 if (lda1>m || lda2>m) { 1633 for (j=0; j<n; j++) { 1634 ierr = PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));CHKERRQ(ierr); 1635 } 1636 } else { 1637 ierr = PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1638 } 1639 PetscFunctionReturn(0); 1640 } 1641 1642 #undef __FUNCT__ 1643 #define __FUNCT__ "MatSetUp_SeqDense" 1644 PetscErrorCode MatSetUp_SeqDense(Mat A) 1645 { 1646 PetscErrorCode ierr; 1647 1648 PetscFunctionBegin; 1649 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1650 PetscFunctionReturn(0); 1651 } 1652 1653 #undef __FUNCT__ 1654 #define __FUNCT__ "MatConjugate_SeqDense" 1655 static PetscErrorCode MatConjugate_SeqDense(Mat A) 1656 { 1657 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1658 PetscInt i,nz = A->rmap->n*A->cmap->n; 1659 PetscScalar *aa = a->v; 1660 1661 PetscFunctionBegin; 1662 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 1663 PetscFunctionReturn(0); 1664 } 1665 1666 #undef __FUNCT__ 1667 #define __FUNCT__ "MatRealPart_SeqDense" 1668 static PetscErrorCode MatRealPart_SeqDense(Mat A) 1669 { 1670 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1671 PetscInt i,nz = A->rmap->n*A->cmap->n; 1672 PetscScalar *aa = a->v; 1673 1674 PetscFunctionBegin; 1675 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1676 PetscFunctionReturn(0); 1677 } 1678 1679 #undef __FUNCT__ 1680 #define __FUNCT__ "MatImaginaryPart_SeqDense" 1681 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 1682 { 1683 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1684 PetscInt i,nz = A->rmap->n*A->cmap->n; 1685 PetscScalar *aa = a->v; 1686 1687 PetscFunctionBegin; 1688 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1689 PetscFunctionReturn(0); 1690 } 1691 1692 /* ----------------------------------------------------------------*/ 1693 #undef __FUNCT__ 1694 #define __FUNCT__ "MatMatMult_SeqDense_SeqDense" 1695 PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1696 { 1697 PetscErrorCode ierr; 1698 1699 PetscFunctionBegin; 1700 if (scall == MAT_INITIAL_MATRIX){ 1701 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1702 } 1703 ierr = MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1704 PetscFunctionReturn(0); 1705 } 1706 1707 #undef __FUNCT__ 1708 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqDense" 1709 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1710 { 1711 PetscErrorCode ierr; 1712 PetscInt m=A->rmap->n,n=B->cmap->n; 1713 Mat Cmat; 1714 1715 PetscFunctionBegin; 1716 if (A->cmap->n != B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"A->cmap->n %d != B->rmap->n %d\n",A->cmap->n,B->rmap->n); 1717 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1718 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1719 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1720 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 1721 Cmat->assembled = PETSC_TRUE; 1722 Cmat->ops->matmult = MatMatMult_SeqDense_SeqDense; 1723 *C = Cmat; 1724 PetscFunctionReturn(0); 1725 } 1726 1727 #undef __FUNCT__ 1728 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqDense" 1729 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1730 { 1731 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1732 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1733 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1734 PetscBLASInt m,n,k; 1735 PetscScalar _DOne=1.0,_DZero=0.0; 1736 1737 PetscFunctionBegin; 1738 m = PetscBLASIntCast(A->rmap->n); 1739 n = PetscBLASIntCast(B->cmap->n); 1740 k = PetscBLASIntCast(A->cmap->n); 1741 BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda); 1742 PetscFunctionReturn(0); 1743 } 1744 1745 #undef __FUNCT__ 1746 #define __FUNCT__ "MatTransposeMatMult_SeqDense_SeqDense" 1747 PetscErrorCode MatTransposeMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1748 { 1749 PetscErrorCode ierr; 1750 1751 PetscFunctionBegin; 1752 if (scall == MAT_INITIAL_MATRIX){ 1753 ierr = MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1754 } 1755 ierr = MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1756 PetscFunctionReturn(0); 1757 } 1758 1759 #undef __FUNCT__ 1760 #define __FUNCT__ "MatTransposeMatMultSymbolic_SeqDense_SeqDense" 1761 PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1762 { 1763 PetscErrorCode ierr; 1764 PetscInt m=A->cmap->n,n=B->cmap->n; 1765 Mat Cmat; 1766 1767 PetscFunctionBegin; 1768 if (A->rmap->n != B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"A->rmap->n %d != B->rmap->n %d\n",A->rmap->n,B->rmap->n); 1769 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1770 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1771 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1772 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 1773 Cmat->assembled = PETSC_TRUE; 1774 *C = Cmat; 1775 PetscFunctionReturn(0); 1776 } 1777 1778 #undef __FUNCT__ 1779 #define __FUNCT__ "MatTransposeMatMultNumeric_SeqDense_SeqDense" 1780 PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1781 { 1782 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1783 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1784 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1785 PetscBLASInt m,n,k; 1786 PetscScalar _DOne=1.0,_DZero=0.0; 1787 1788 PetscFunctionBegin; 1789 m = PetscBLASIntCast(A->cmap->n); 1790 n = PetscBLASIntCast(B->cmap->n); 1791 k = PetscBLASIntCast(A->rmap->n); 1792 /* 1793 Note the m and n arguments below are the number rows and columns of A', not A! 1794 */ 1795 BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda); 1796 PetscFunctionReturn(0); 1797 } 1798 1799 #undef __FUNCT__ 1800 #define __FUNCT__ "MatGetRowMax_SeqDense" 1801 PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 1802 { 1803 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1804 PetscErrorCode ierr; 1805 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1806 PetscScalar *x; 1807 MatScalar *aa = a->v; 1808 1809 PetscFunctionBegin; 1810 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1811 1812 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1813 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1814 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1815 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1816 for (i=0; i<m; i++) { 1817 x[i] = aa[i]; if (idx) idx[i] = 0; 1818 for (j=1; j<n; j++){ 1819 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1820 } 1821 } 1822 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1823 PetscFunctionReturn(0); 1824 } 1825 1826 #undef __FUNCT__ 1827 #define __FUNCT__ "MatGetRowMaxAbs_SeqDense" 1828 PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 1829 { 1830 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1831 PetscErrorCode ierr; 1832 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1833 PetscScalar *x; 1834 PetscReal atmp; 1835 MatScalar *aa = a->v; 1836 1837 PetscFunctionBegin; 1838 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1839 1840 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1841 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1842 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1843 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1844 for (i=0; i<m; i++) { 1845 x[i] = PetscAbsScalar(aa[i]); 1846 for (j=1; j<n; j++){ 1847 atmp = PetscAbsScalar(aa[i+m*j]); 1848 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 1849 } 1850 } 1851 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1852 PetscFunctionReturn(0); 1853 } 1854 1855 #undef __FUNCT__ 1856 #define __FUNCT__ "MatGetRowMin_SeqDense" 1857 PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 1858 { 1859 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1860 PetscErrorCode ierr; 1861 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1862 PetscScalar *x; 1863 MatScalar *aa = a->v; 1864 1865 PetscFunctionBegin; 1866 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1867 1868 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1869 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1870 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1871 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1872 for (i=0; i<m; i++) { 1873 x[i] = aa[i]; if (idx) idx[i] = 0; 1874 for (j=1; j<n; j++){ 1875 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1876 } 1877 } 1878 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1879 PetscFunctionReturn(0); 1880 } 1881 1882 #undef __FUNCT__ 1883 #define __FUNCT__ "MatGetColumnVector_SeqDense" 1884 PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 1885 { 1886 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1887 PetscErrorCode ierr; 1888 PetscScalar *x; 1889 1890 PetscFunctionBegin; 1891 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1892 1893 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1894 ierr = PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1895 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1896 PetscFunctionReturn(0); 1897 } 1898 1899 1900 #undef __FUNCT__ 1901 #define __FUNCT__ "MatGetColumnNorms_SeqDense" 1902 PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms) 1903 { 1904 PetscErrorCode ierr; 1905 PetscInt i,j,m,n; 1906 PetscScalar *a; 1907 1908 PetscFunctionBegin; 1909 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 1910 ierr = PetscMemzero(norms,n*sizeof(PetscReal));CHKERRQ(ierr); 1911 ierr = MatDenseGetArray(A,&a);CHKERRQ(ierr); 1912 if (type == NORM_2) { 1913 for (i=0; i<n; i++ ){ 1914 for (j=0; j<m; j++) { 1915 norms[i] += PetscAbsScalar(a[j]*a[j]); 1916 } 1917 a += m; 1918 } 1919 } else if (type == NORM_1) { 1920 for (i=0; i<n; i++ ){ 1921 for (j=0; j<m; j++) { 1922 norms[i] += PetscAbsScalar(a[j]); 1923 } 1924 a += m; 1925 } 1926 } else if (type == NORM_INFINITY) { 1927 for (i=0; i<n; i++ ){ 1928 for (j=0; j<m; j++) { 1929 norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]); 1930 } 1931 a += m; 1932 } 1933 } else SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"Unknown NormType"); 1934 ierr = MatDenseRestoreArray(A,&a);CHKERRQ(ierr); 1935 if (type == NORM_2) { 1936 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 1937 } 1938 PetscFunctionReturn(0); 1939 } 1940 1941 #undef __FUNCT__ 1942 #define __FUNCT__ "MatSetRandom_SeqDense" 1943 static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx) 1944 { 1945 PetscErrorCode ierr; 1946 PetscScalar *a; 1947 PetscInt m,n,i; 1948 1949 PetscFunctionBegin; 1950 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 1951 ierr = MatDenseGetArray(x,&a);CHKERRQ(ierr); 1952 for (i=0; i<m*n; i++) { 1953 ierr = PetscRandomGetValue(rctx,a+i);CHKERRQ(ierr); 1954 } 1955 ierr = MatDenseRestoreArray(x,&a);CHKERRQ(ierr); 1956 PetscFunctionReturn(0); 1957 } 1958 1959 1960 /* -------------------------------------------------------------------*/ 1961 static struct _MatOps MatOps_Values = {MatSetValues_SeqDense, 1962 MatGetRow_SeqDense, 1963 MatRestoreRow_SeqDense, 1964 MatMult_SeqDense, 1965 /* 4*/ MatMultAdd_SeqDense, 1966 MatMultTranspose_SeqDense, 1967 MatMultTransposeAdd_SeqDense, 1968 0, 1969 0, 1970 0, 1971 /*10*/ 0, 1972 MatLUFactor_SeqDense, 1973 MatCholeskyFactor_SeqDense, 1974 MatSOR_SeqDense, 1975 MatTranspose_SeqDense, 1976 /*15*/ MatGetInfo_SeqDense, 1977 MatEqual_SeqDense, 1978 MatGetDiagonal_SeqDense, 1979 MatDiagonalScale_SeqDense, 1980 MatNorm_SeqDense, 1981 /*20*/ MatAssemblyBegin_SeqDense, 1982 MatAssemblyEnd_SeqDense, 1983 MatSetOption_SeqDense, 1984 MatZeroEntries_SeqDense, 1985 /*24*/ MatZeroRows_SeqDense, 1986 0, 1987 0, 1988 0, 1989 0, 1990 /*29*/ MatSetUp_SeqDense, 1991 0, 1992 0, 1993 0, 1994 0, 1995 /*34*/ MatDuplicate_SeqDense, 1996 0, 1997 0, 1998 0, 1999 0, 2000 /*39*/ MatAXPY_SeqDense, 2001 MatGetSubMatrices_SeqDense, 2002 0, 2003 MatGetValues_SeqDense, 2004 MatCopy_SeqDense, 2005 /*44*/ MatGetRowMax_SeqDense, 2006 MatScale_SeqDense, 2007 0, 2008 0, 2009 0, 2010 /*49*/ MatSetRandom_SeqDense, 2011 0, 2012 0, 2013 0, 2014 0, 2015 /*54*/ 0, 2016 0, 2017 0, 2018 0, 2019 0, 2020 /*59*/ 0, 2021 MatDestroy_SeqDense, 2022 MatView_SeqDense, 2023 0, 2024 0, 2025 /*64*/ 0, 2026 0, 2027 0, 2028 0, 2029 0, 2030 /*69*/ MatGetRowMaxAbs_SeqDense, 2031 0, 2032 0, 2033 0, 2034 0, 2035 /*74*/ 0, 2036 0, 2037 0, 2038 0, 2039 0, 2040 /*79*/ 0, 2041 0, 2042 0, 2043 0, 2044 /*83*/ MatLoad_SeqDense, 2045 0, 2046 MatIsHermitian_SeqDense, 2047 0, 2048 0, 2049 0, 2050 /*89*/ MatMatMult_SeqDense_SeqDense, 2051 MatMatMultSymbolic_SeqDense_SeqDense, 2052 MatMatMultNumeric_SeqDense_SeqDense, 2053 0, 2054 0, 2055 /*94*/ 0, 2056 0, 2057 0, 2058 0, 2059 0, 2060 /*99*/ 0, 2061 0, 2062 0, 2063 MatConjugate_SeqDense, 2064 0, 2065 /*104*/0, 2066 MatRealPart_SeqDense, 2067 MatImaginaryPart_SeqDense, 2068 0, 2069 0, 2070 /*109*/MatMatSolve_SeqDense, 2071 0, 2072 MatGetRowMin_SeqDense, 2073 MatGetColumnVector_SeqDense, 2074 0, 2075 /*114*/0, 2076 0, 2077 0, 2078 0, 2079 0, 2080 /*119*/0, 2081 0, 2082 0, 2083 0, 2084 0, 2085 /*124*/0, 2086 MatGetColumnNorms_SeqDense, 2087 0, 2088 0, 2089 0, 2090 /*129*/0, 2091 MatTransposeMatMult_SeqDense_SeqDense, 2092 MatTransposeMatMultSymbolic_SeqDense_SeqDense, 2093 MatTransposeMatMultNumeric_SeqDense_SeqDense, 2094 }; 2095 2096 #undef __FUNCT__ 2097 #define __FUNCT__ "MatCreateSeqDense" 2098 /*@C 2099 MatCreateSeqDense - Creates a sequential dense matrix that 2100 is stored in column major order (the usual Fortran 77 manner). Many 2101 of the matrix operations use the BLAS and LAPACK routines. 2102 2103 Collective on MPI_Comm 2104 2105 Input Parameters: 2106 + comm - MPI communicator, set to PETSC_COMM_SELF 2107 . m - number of rows 2108 . n - number of columns 2109 - data - optional location of matrix data in column major order. Set data=PETSC_NULL for PETSc 2110 to control all matrix memory allocation. 2111 2112 Output Parameter: 2113 . A - the matrix 2114 2115 Notes: 2116 The data input variable is intended primarily for Fortran programmers 2117 who wish to allocate their own matrix memory space. Most users should 2118 set data=PETSC_NULL. 2119 2120 Level: intermediate 2121 2122 .keywords: dense, matrix, LAPACK, BLAS 2123 2124 .seealso: MatCreate(), MatCreateDense(), MatSetValues() 2125 @*/ 2126 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 2127 { 2128 PetscErrorCode ierr; 2129 2130 PetscFunctionBegin; 2131 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2132 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2133 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 2134 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 2135 PetscFunctionReturn(0); 2136 } 2137 2138 #undef __FUNCT__ 2139 #define __FUNCT__ "MatSeqDenseSetPreallocation" 2140 /*@C 2141 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 2142 2143 Collective on MPI_Comm 2144 2145 Input Parameters: 2146 + A - the matrix 2147 - data - the array (or PETSC_NULL) 2148 2149 Notes: 2150 The data input variable is intended primarily for Fortran programmers 2151 who wish to allocate their own matrix memory space. Most users should 2152 need not call this routine. 2153 2154 Level: intermediate 2155 2156 .keywords: dense, matrix, LAPACK, BLAS 2157 2158 .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA() 2159 2160 @*/ 2161 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 2162 { 2163 PetscErrorCode ierr; 2164 2165 PetscFunctionBegin; 2166 ierr = PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));CHKERRQ(ierr); 2167 PetscFunctionReturn(0); 2168 } 2169 2170 EXTERN_C_BEGIN 2171 #undef __FUNCT__ 2172 #define __FUNCT__ "MatSeqDenseSetPreallocation_SeqDense" 2173 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 2174 { 2175 Mat_SeqDense *b; 2176 PetscErrorCode ierr; 2177 2178 PetscFunctionBegin; 2179 B->preallocated = PETSC_TRUE; 2180 2181 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2182 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2183 2184 b = (Mat_SeqDense*)B->data; 2185 b->Mmax = B->rmap->n; 2186 b->Nmax = B->cmap->n; 2187 if(b->lda <= 0 || b->changelda) b->lda = B->rmap->n; 2188 2189 if (!data) { /* petsc-allocated storage */ 2190 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2191 ierr = PetscMalloc(b->lda*b->Nmax*sizeof(PetscScalar),&b->v);CHKERRQ(ierr); 2192 ierr = PetscMemzero(b->v,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 2193 ierr = PetscLogObjectMemory(B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 2194 b->user_alloc = PETSC_FALSE; 2195 } else { /* user-allocated storage */ 2196 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2197 b->v = data; 2198 b->user_alloc = PETSC_TRUE; 2199 } 2200 B->assembled = PETSC_TRUE; 2201 PetscFunctionReturn(0); 2202 } 2203 EXTERN_C_END 2204 2205 #undef __FUNCT__ 2206 #define __FUNCT__ "MatSeqDenseSetLDA" 2207 /*@C 2208 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 2209 2210 Input parameter: 2211 + A - the matrix 2212 - lda - the leading dimension 2213 2214 Notes: 2215 This routine is to be used in conjunction with MatSeqDenseSetPreallocation(); 2216 it asserts that the preallocation has a leading dimension (the LDA parameter 2217 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 2218 2219 Level: intermediate 2220 2221 .keywords: dense, matrix, LAPACK, BLAS 2222 2223 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 2224 2225 @*/ 2226 PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda) 2227 { 2228 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2229 2230 PetscFunctionBegin; 2231 if (lda < B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"LDA %D must be at least matrix dimension %D",lda,B->rmap->n); 2232 b->lda = lda; 2233 b->changelda = PETSC_FALSE; 2234 b->Mmax = PetscMax(b->Mmax,lda); 2235 PetscFunctionReturn(0); 2236 } 2237 2238 /*MC 2239 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 2240 2241 Options Database Keys: 2242 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 2243 2244 Level: beginner 2245 2246 .seealso: MatCreateSeqDense() 2247 2248 M*/ 2249 2250 EXTERN_C_BEGIN 2251 #undef __FUNCT__ 2252 #define __FUNCT__ "MatCreate_SeqDense" 2253 PetscErrorCode MatCreate_SeqDense(Mat B) 2254 { 2255 Mat_SeqDense *b; 2256 PetscErrorCode ierr; 2257 PetscMPIInt size; 2258 2259 PetscFunctionBegin; 2260 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 2261 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 2262 2263 ierr = PetscNewLog(B,Mat_SeqDense,&b);CHKERRQ(ierr); 2264 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2265 B->data = (void*)b; 2266 2267 b->pivots = 0; 2268 b->roworiented = PETSC_TRUE; 2269 b->v = 0; 2270 b->changelda = PETSC_FALSE; 2271 2272 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatDenseGetArray_C","MatDenseGetArray_SeqDense",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2273 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatDenseRestoreArray_C","MatDenseRestoreArray_SeqDense",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2274 2275 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqdense_seqaij_C","MatConvert_SeqDense_SeqAIJ",MatConvert_SeqDense_SeqAIJ);CHKERRQ(ierr); 2276 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C", 2277 "MatGetFactor_seqdense_petsc", 2278 MatGetFactor_seqdense_petsc);CHKERRQ(ierr); 2279 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqDenseSetPreallocation_C", 2280 "MatSeqDenseSetPreallocation_SeqDense", 2281 MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 2282 2283 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqaij_seqdense_C", 2284 "MatMatMult_SeqAIJ_SeqDense", 2285 MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2286 2287 2288 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C", 2289 "MatMatMultSymbolic_SeqAIJ_SeqDense", 2290 MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2291 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C", 2292 "MatMatMultNumeric_SeqAIJ_SeqDense", 2293 MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2294 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 2295 PetscFunctionReturn(0); 2296 } 2297 EXTERN_C_END 2298