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 11 PetscErrorCode MatSeqDenseSymmetrize_Private(Mat A, PetscBool hermitian) 12 { 13 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 14 PetscInt j, k, n = A->rmap->n; 15 PetscScalar *v; 16 PetscErrorCode ierr; 17 18 PetscFunctionBegin; 19 if (A->rmap->n != A->cmap->n) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot symmetrize a rectangular matrix"); 20 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 21 if (!hermitian) { 22 for (k=0;k<n;k++) { 23 for (j=k;j<n;j++) { 24 v[j*mat->lda + k] = v[k*mat->lda + j]; 25 } 26 } 27 } else { 28 for (k=0;k<n;k++) { 29 for (j=k;j<n;j++) { 30 v[j*mat->lda + k] = PetscConj(v[k*mat->lda + j]); 31 } 32 } 33 } 34 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 35 PetscFunctionReturn(0); 36 } 37 38 PETSC_EXTERN PetscErrorCode MatSeqDenseInvertFactors_Private(Mat A) 39 { 40 #if defined(PETSC_MISSING_LAPACK_POTRF) 41 PetscFunctionBegin; 42 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable."); 43 #else 44 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 45 PetscErrorCode ierr; 46 PetscBLASInt info,n; 47 48 PetscFunctionBegin; 49 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 50 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 51 if (A->factortype == MAT_FACTOR_LU) { 52 if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present"); 53 if (!mat->fwork) { 54 mat->lfwork = n; 55 ierr = PetscMalloc1(mat->lfwork,&mat->fwork);CHKERRQ(ierr); 56 ierr = PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));CHKERRQ(ierr); 57 } 58 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 59 PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info)); 60 ierr = PetscFPTrapPop();CHKERRQ(ierr); 61 ierr = PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);CHKERRQ(ierr); 62 } else if (A->factortype == MAT_FACTOR_CHOLESKY) { 63 if (A->spd) { 64 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 65 PetscStackCallBLAS("LAPACKpotri",LAPACKpotri_("L",&n,mat->v,&mat->lda,&info)); 66 ierr = PetscFPTrapPop();CHKERRQ(ierr); 67 ierr = MatSeqDenseSymmetrize_Private(A,PETSC_TRUE);CHKERRQ(ierr); 68 #if defined(PETSC_USE_COMPLEX) 69 } else if (A->hermitian) { 70 if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present"); 71 if (!mat->fwork) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Fwork not present"); 72 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 73 PetscStackCallBLAS("LAPACKhetri",LAPACKhetri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info)); 74 ierr = PetscFPTrapPop();CHKERRQ(ierr); 75 ierr = MatSeqDenseSymmetrize_Private(A,PETSC_TRUE);CHKERRQ(ierr); 76 #endif 77 } else { /* symmetric case */ 78 if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present"); 79 if (!mat->fwork) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Fwork not present"); 80 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 81 PetscStackCallBLAS("LAPACKsytri",LAPACKsytri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info)); 82 ierr = PetscFPTrapPop();CHKERRQ(ierr); 83 ierr = MatSeqDenseSymmetrize_Private(A,PETSC_FALSE);CHKERRQ(ierr); 84 } 85 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad Inversion: zero pivot in row %D",(PetscInt)info-1); 86 ierr = PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);CHKERRQ(ierr); 87 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 88 #endif 89 90 A->ops->solve = NULL; 91 A->ops->matsolve = NULL; 92 A->ops->solvetranspose = NULL; 93 A->ops->matsolvetranspose = NULL; 94 A->ops->solveadd = NULL; 95 A->ops->solvetransposeadd = NULL; 96 A->factortype = MAT_FACTOR_NONE; 97 ierr = PetscFree(A->solvertype);CHKERRQ(ierr); 98 PetscFunctionReturn(0); 99 } 100 101 PetscErrorCode MatZeroRowsColumns_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 102 { 103 PetscErrorCode ierr; 104 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 105 PetscInt m = l->lda, n = A->cmap->n,r = A->rmap->n, i,j; 106 PetscScalar *slot,*bb,*v; 107 const PetscScalar *xx; 108 109 PetscFunctionBegin; 110 #if defined(PETSC_USE_DEBUG) 111 for (i=0; i<N; i++) { 112 if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 113 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); 114 if (rows[i] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Col %D requested to be zeroed greater than or equal number of cols %D",rows[i],A->cmap->n); 115 } 116 #endif 117 if (!N) PetscFunctionReturn(0); 118 119 /* fix right hand side if needed */ 120 if (x && b) { 121 Vec xt; 122 123 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 124 ierr = VecDuplicate(x,&xt);CHKERRQ(ierr); 125 ierr = VecCopy(x,xt);CHKERRQ(ierr); 126 ierr = VecScale(xt,-1.0);CHKERRQ(ierr); 127 ierr = MatMultAdd(A,xt,b,b);CHKERRQ(ierr); 128 ierr = VecDestroy(&xt);CHKERRQ(ierr); 129 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 130 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 131 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]]; 132 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 133 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 134 } 135 136 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 137 for (i=0; i<N; i++) { 138 slot = v + rows[i]*m; 139 ierr = PetscArrayzero(slot,r);CHKERRQ(ierr); 140 } 141 for (i=0; i<N; i++) { 142 slot = v + rows[i]; 143 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 144 } 145 if (diag != 0.0) { 146 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 147 for (i=0; i<N; i++) { 148 slot = v + (m+1)*rows[i]; 149 *slot = diag; 150 } 151 } 152 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 153 PetscFunctionReturn(0); 154 } 155 156 PetscErrorCode MatPtAPNumeric_SeqDense_SeqDense(Mat A,Mat P,Mat C) 157 { 158 Mat_SeqDense *c = (Mat_SeqDense*)(C->data); 159 PetscErrorCode ierr; 160 161 PetscFunctionBegin; 162 if (c->ptapwork) { 163 ierr = (*C->ops->matmultnumeric)(A,P,c->ptapwork);CHKERRQ(ierr); 164 ierr = (*C->ops->transposematmultnumeric)(P,c->ptapwork,C);CHKERRQ(ierr); 165 } else { /* first time went trough the basic. Should we add better dispatching for subclasses? */ 166 ierr = MatPtAP_Basic(A,P,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr); 167 } 168 PetscFunctionReturn(0); 169 } 170 171 PetscErrorCode MatPtAPSymbolic_SeqDense_SeqDense(Mat A,Mat P,PetscReal fill,Mat *C) 172 { 173 Mat_SeqDense *c; 174 PetscBool flg; 175 PetscErrorCode ierr; 176 177 PetscFunctionBegin; 178 ierr = PetscObjectTypeCompare((PetscObject)P,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 179 ierr = MatCreate(PetscObjectComm((PetscObject)A),C);CHKERRQ(ierr); 180 ierr = MatSetSizes(*C,P->cmap->n,P->cmap->n,P->cmap->N,P->cmap->N);CHKERRQ(ierr); 181 ierr = MatSetType(*C,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 182 ierr = MatSeqDenseSetPreallocation(*C,NULL);CHKERRQ(ierr); 183 c = (Mat_SeqDense*)((*C)->data); 184 ierr = MatCreate(PetscObjectComm((PetscObject)A),&c->ptapwork);CHKERRQ(ierr); 185 ierr = MatSetSizes(c->ptapwork,A->rmap->n,P->cmap->n,A->rmap->N,P->cmap->N);CHKERRQ(ierr); 186 ierr = MatSetType(c->ptapwork,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 187 ierr = MatSeqDenseSetPreallocation(c->ptapwork,NULL);CHKERRQ(ierr); 188 PetscFunctionReturn(0); 189 } 190 191 PETSC_INTERN PetscErrorCode MatPtAP_SeqDense_SeqDense(Mat A,Mat P,MatReuse reuse,PetscReal fill,Mat *C) 192 { 193 PetscErrorCode ierr; 194 195 PetscFunctionBegin; 196 if (reuse == MAT_INITIAL_MATRIX) { 197 ierr = MatPtAPSymbolic_SeqDense_SeqDense(A,P,fill,C);CHKERRQ(ierr); 198 } 199 ierr = PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);CHKERRQ(ierr); 200 ierr = (*(*C)->ops->ptapnumeric)(A,P,*C);CHKERRQ(ierr); 201 ierr = PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);CHKERRQ(ierr); 202 PetscFunctionReturn(0); 203 } 204 205 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 206 { 207 Mat B = NULL; 208 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 209 Mat_SeqDense *b; 210 PetscErrorCode ierr; 211 PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->N,i; 212 MatScalar *av=a->a; 213 PetscBool isseqdense; 214 215 PetscFunctionBegin; 216 if (reuse == MAT_REUSE_MATRIX) { 217 ierr = PetscObjectTypeCompare((PetscObject)*newmat,MATSEQDENSE,&isseqdense);CHKERRQ(ierr); 218 if (!isseqdense) SETERRQ1(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_USER,"Cannot reuse matrix of type %s",((PetscObject)(*newmat))->type_name); 219 } 220 if (reuse != MAT_REUSE_MATRIX) { 221 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 222 ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr); 223 ierr = MatSetType(B,MATSEQDENSE);CHKERRQ(ierr); 224 ierr = MatSeqDenseSetPreallocation(B,NULL);CHKERRQ(ierr); 225 b = (Mat_SeqDense*)(B->data); 226 } else { 227 b = (Mat_SeqDense*)((*newmat)->data); 228 ierr = PetscArrayzero(b->v,m*n);CHKERRQ(ierr); 229 } 230 for (i=0; i<m; i++) { 231 PetscInt j; 232 for (j=0;j<ai[1]-ai[0];j++) { 233 b->v[*aj*m+i] = *av; 234 aj++; 235 av++; 236 } 237 ai++; 238 } 239 240 if (reuse == MAT_INPLACE_MATRIX) { 241 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 242 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 243 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 244 } else { 245 if (B) *newmat = B; 246 ierr = MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 247 ierr = MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 248 } 249 PetscFunctionReturn(0); 250 } 251 252 PETSC_INTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 253 { 254 Mat B; 255 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 256 PetscErrorCode ierr; 257 PetscInt i, j; 258 PetscInt *rows, *nnz; 259 MatScalar *aa = a->v, *vals; 260 261 PetscFunctionBegin; 262 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 263 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 264 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 265 ierr = PetscCalloc3(A->rmap->n,&rows,A->rmap->n,&nnz,A->rmap->n,&vals);CHKERRQ(ierr); 266 for (j=0; j<A->cmap->n; j++) { 267 for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) ++nnz[i]; 268 aa += a->lda; 269 } 270 ierr = MatSeqAIJSetPreallocation(B,PETSC_DETERMINE,nnz);CHKERRQ(ierr); 271 aa = a->v; 272 for (j=0; j<A->cmap->n; j++) { 273 PetscInt numRows = 0; 274 for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) {rows[numRows] = i; vals[numRows++] = aa[i];} 275 ierr = MatSetValues(B,numRows,rows,1,&j,vals,INSERT_VALUES);CHKERRQ(ierr); 276 aa += a->lda; 277 } 278 ierr = PetscFree3(rows,nnz,vals);CHKERRQ(ierr); 279 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 280 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 281 282 if (reuse == MAT_INPLACE_MATRIX) { 283 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 284 } else { 285 *newmat = B; 286 } 287 PetscFunctionReturn(0); 288 } 289 290 PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str) 291 { 292 Mat_SeqDense *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data; 293 const PetscScalar *xv; 294 PetscScalar *yv; 295 PetscBLASInt N,m,ldax,lday,one = 1; 296 PetscErrorCode ierr; 297 298 PetscFunctionBegin; 299 ierr = MatDenseGetArrayRead(X,&xv);CHKERRQ(ierr); 300 ierr = MatDenseGetArray(Y,&yv);CHKERRQ(ierr); 301 ierr = PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);CHKERRQ(ierr); 302 ierr = PetscBLASIntCast(X->rmap->n,&m);CHKERRQ(ierr); 303 ierr = PetscBLASIntCast(x->lda,&ldax);CHKERRQ(ierr); 304 ierr = PetscBLASIntCast(y->lda,&lday);CHKERRQ(ierr); 305 if (ldax>m || lday>m) { 306 PetscInt j; 307 308 for (j=0; j<X->cmap->n; j++) { 309 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&alpha,xv+j*ldax,&one,yv+j*lday,&one)); 310 } 311 } else { 312 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&alpha,xv,&one,yv,&one)); 313 } 314 ierr = MatDenseRestoreArrayRead(X,&xv);CHKERRQ(ierr); 315 ierr = MatDenseRestoreArray(Y,&yv);CHKERRQ(ierr); 316 ierr = PetscLogFlops(PetscMax(2*N-1,0));CHKERRQ(ierr); 317 PetscFunctionReturn(0); 318 } 319 320 static PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info) 321 { 322 PetscLogDouble N = A->rmap->n*A->cmap->n; 323 324 PetscFunctionBegin; 325 info->block_size = 1.0; 326 info->nz_allocated = N; 327 info->nz_used = N; 328 info->nz_unneeded = 0; 329 info->assemblies = A->num_ass; 330 info->mallocs = 0; 331 info->memory = ((PetscObject)A)->mem; 332 info->fill_ratio_given = 0; 333 info->fill_ratio_needed = 0; 334 info->factor_mallocs = 0; 335 PetscFunctionReturn(0); 336 } 337 338 static PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha) 339 { 340 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 341 PetscScalar *v; 342 PetscErrorCode ierr; 343 PetscBLASInt one = 1,j,nz,lda; 344 345 PetscFunctionBegin; 346 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 347 ierr = PetscBLASIntCast(a->lda,&lda);CHKERRQ(ierr); 348 if (lda>A->rmap->n) { 349 ierr = PetscBLASIntCast(A->rmap->n,&nz);CHKERRQ(ierr); 350 for (j=0; j<A->cmap->n; j++) { 351 PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&alpha,v+j*lda,&one)); 352 } 353 } else { 354 ierr = PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);CHKERRQ(ierr); 355 PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&alpha,v,&one)); 356 } 357 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 358 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 359 PetscFunctionReturn(0); 360 } 361 362 static PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool *fl) 363 { 364 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 365 PetscInt i,j,m = A->rmap->n,N = a->lda; 366 const PetscScalar *v; 367 PetscErrorCode ierr; 368 369 PetscFunctionBegin; 370 *fl = PETSC_FALSE; 371 if (A->rmap->n != A->cmap->n) PetscFunctionReturn(0); 372 ierr = MatDenseGetArrayRead(A,&v);CHKERRQ(ierr); 373 for (i=0; i<m; i++) { 374 for (j=i; j<m; j++) { 375 if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) PetscFunctionReturn(0); 376 } 377 } 378 ierr = MatDenseRestoreArrayRead(A,&v);CHKERRQ(ierr); 379 *fl = PETSC_TRUE; 380 PetscFunctionReturn(0); 381 } 382 383 PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues) 384 { 385 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 386 PetscErrorCode ierr; 387 PetscInt lda = (PetscInt)mat->lda,j,m; 388 389 PetscFunctionBegin; 390 ierr = PetscLayoutReference(A->rmap,&newi->rmap);CHKERRQ(ierr); 391 ierr = PetscLayoutReference(A->cmap,&newi->cmap);CHKERRQ(ierr); 392 ierr = MatSeqDenseSetPreallocation(newi,NULL);CHKERRQ(ierr); 393 if (cpvalues == MAT_COPY_VALUES) { 394 const PetscScalar *av; 395 PetscScalar *v; 396 397 ierr = MatDenseGetArrayRead(A,&av);CHKERRQ(ierr); 398 ierr = MatDenseGetArray(newi,&v);CHKERRQ(ierr); 399 if (lda>A->rmap->n) { 400 m = A->rmap->n; 401 for (j=0; j<A->cmap->n; j++) { 402 ierr = PetscArraycpy(v+j*m,av+j*lda,m);CHKERRQ(ierr); 403 } 404 } else { 405 ierr = PetscArraycpy(v,av,A->rmap->n*A->cmap->n);CHKERRQ(ierr); 406 } 407 ierr = MatDenseRestoreArray(newi,&v);CHKERRQ(ierr); 408 ierr = MatDenseRestoreArrayRead(A,&av);CHKERRQ(ierr); 409 } 410 PetscFunctionReturn(0); 411 } 412 413 PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat) 414 { 415 PetscErrorCode ierr; 416 417 PetscFunctionBegin; 418 ierr = MatCreate(PetscObjectComm((PetscObject)A),newmat);CHKERRQ(ierr); 419 ierr = MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 420 ierr = MatSetType(*newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 421 ierr = MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);CHKERRQ(ierr); 422 PetscFunctionReturn(0); 423 } 424 425 static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy) 426 { 427 MatFactorInfo info; 428 PetscErrorCode ierr; 429 430 PetscFunctionBegin; 431 ierr = MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);CHKERRQ(ierr); 432 ierr = (*fact->ops->lufactor)(fact,0,0,&info);CHKERRQ(ierr); 433 PetscFunctionReturn(0); 434 } 435 436 static PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy) 437 { 438 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 439 PetscErrorCode ierr; 440 const PetscScalar *x; 441 PetscScalar *y; 442 PetscBLASInt one = 1,info,m; 443 444 PetscFunctionBegin; 445 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 446 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 447 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 448 ierr = PetscArraycpy(y,x,A->rmap->n);CHKERRQ(ierr); 449 if (A->factortype == MAT_FACTOR_LU) { 450 #if defined(PETSC_MISSING_LAPACK_GETRS) 451 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 452 #else 453 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 454 PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 455 ierr = PetscFPTrapPop();CHKERRQ(ierr); 456 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve"); 457 #endif 458 } else if (A->factortype == MAT_FACTOR_CHOLESKY) { 459 #if defined(PETSC_MISSING_LAPACK_POTRS) 460 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 461 #else 462 if (A->spd) { 463 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 464 PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info)); 465 ierr = PetscFPTrapPop();CHKERRQ(ierr); 466 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve"); 467 #if defined(PETSC_USE_COMPLEX) 468 } else if (A->hermitian) { 469 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 470 PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 471 ierr = PetscFPTrapPop();CHKERRQ(ierr); 472 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve"); 473 #endif 474 } else { /* symmetric case */ 475 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 476 PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 477 ierr = PetscFPTrapPop();CHKERRQ(ierr); 478 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve"); 479 } 480 #endif 481 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 482 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 483 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 484 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 485 PetscFunctionReturn(0); 486 } 487 488 static PetscErrorCode MatMatSolve_SeqDense(Mat A,Mat B,Mat X) 489 { 490 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 491 PetscErrorCode ierr; 492 const PetscScalar *b; 493 PetscScalar *x; 494 PetscInt n; 495 PetscBLASInt nrhs,info,m; 496 497 PetscFunctionBegin; 498 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 499 ierr = MatGetSize(B,NULL,&n);CHKERRQ(ierr); 500 ierr = PetscBLASIntCast(n,&nrhs);CHKERRQ(ierr); 501 ierr = MatDenseGetArrayRead(B,&b);CHKERRQ(ierr); 502 ierr = MatDenseGetArray(X,&x);CHKERRQ(ierr); 503 504 ierr = PetscArraycpy(x,b,m*nrhs);CHKERRQ(ierr); 505 506 if (A->factortype == MAT_FACTOR_LU) { 507 #if defined(PETSC_MISSING_LAPACK_GETRS) 508 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 509 #else 510 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 511 PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info)); 512 ierr = PetscFPTrapPop();CHKERRQ(ierr); 513 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve"); 514 #endif 515 } else if (A->factortype == MAT_FACTOR_CHOLESKY) { 516 #if defined(PETSC_MISSING_LAPACK_POTRS) 517 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 518 #else 519 if (A->spd) { 520 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 521 PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info)); 522 ierr = PetscFPTrapPop();CHKERRQ(ierr); 523 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve"); 524 #if defined(PETSC_USE_COMPLEX) 525 } else if (A->hermitian) { 526 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 527 PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info)); 528 ierr = PetscFPTrapPop();CHKERRQ(ierr); 529 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve"); 530 #endif 531 } else { /* symmetric case */ 532 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 533 PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info)); 534 ierr = PetscFPTrapPop();CHKERRQ(ierr); 535 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve"); 536 } 537 #endif 538 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 539 540 ierr = MatDenseRestoreArrayRead(B,&b);CHKERRQ(ierr); 541 ierr = MatDenseRestoreArray(X,&x);CHKERRQ(ierr); 542 ierr = PetscLogFlops(nrhs*(2.0*m*m - m));CHKERRQ(ierr); 543 PetscFunctionReturn(0); 544 } 545 546 static PetscErrorCode MatConjugate_SeqDense(Mat); 547 548 static PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy) 549 { 550 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 551 PetscErrorCode ierr; 552 const PetscScalar *x; 553 PetscScalar *y; 554 PetscBLASInt one = 1,info,m; 555 556 PetscFunctionBegin; 557 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 558 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 559 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 560 ierr = PetscArraycpy(y,x,A->rmap->n);CHKERRQ(ierr); 561 if (A->factortype == MAT_FACTOR_LU) { 562 #if defined(PETSC_MISSING_LAPACK_GETRS) 563 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 564 #else 565 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 566 PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 567 ierr = PetscFPTrapPop();CHKERRQ(ierr); 568 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve"); 569 #endif 570 } else if (A->factortype == MAT_FACTOR_CHOLESKY) { 571 #if defined(PETSC_MISSING_LAPACK_POTRS) 572 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 573 #else 574 if (A->spd) { 575 #if defined(PETSC_USE_COMPLEX) 576 ierr = MatConjugate_SeqDense(A);CHKERRQ(ierr); 577 #endif 578 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 579 PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info)); 580 ierr = PetscFPTrapPop();CHKERRQ(ierr); 581 #if defined(PETSC_USE_COMPLEX) 582 ierr = MatConjugate_SeqDense(A);CHKERRQ(ierr); 583 #endif 584 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve"); 585 #if defined(PETSC_USE_COMPLEX) 586 } else if (A->hermitian) { 587 ierr = MatConjugate_SeqDense(A);CHKERRQ(ierr); 588 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 589 PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 590 ierr = PetscFPTrapPop();CHKERRQ(ierr); 591 ierr = MatConjugate_SeqDense(A);CHKERRQ(ierr); 592 #endif 593 } else { /* symmetric case */ 594 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 595 PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 596 ierr = PetscFPTrapPop();CHKERRQ(ierr); 597 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve"); 598 } 599 #endif 600 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 601 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 602 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 603 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 604 PetscFunctionReturn(0); 605 } 606 607 /* ---------------------------------------------------------------*/ 608 /* COMMENT: I have chosen to hide row permutation in the pivots, 609 rather than put it in the Mat->row slot.*/ 610 PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo) 611 { 612 #if defined(PETSC_MISSING_LAPACK_GETRF) 613 PetscFunctionBegin; 614 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable."); 615 #else 616 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 617 PetscErrorCode ierr; 618 PetscBLASInt n,m,info; 619 620 PetscFunctionBegin; 621 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 622 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 623 if (!mat->pivots) { 624 ierr = PetscMalloc1(A->rmap->n,&mat->pivots);CHKERRQ(ierr); 625 ierr = PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));CHKERRQ(ierr); 626 } 627 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 628 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 629 PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info)); 630 ierr = PetscFPTrapPop();CHKERRQ(ierr); 631 632 if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization"); 633 if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization"); 634 635 A->ops->solve = MatSolve_SeqDense; 636 A->ops->matsolve = MatMatSolve_SeqDense; 637 A->ops->solvetranspose = MatSolveTranspose_SeqDense; 638 A->factortype = MAT_FACTOR_LU; 639 640 ierr = PetscFree(A->solvertype);CHKERRQ(ierr); 641 ierr = PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);CHKERRQ(ierr); 642 643 ierr = PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);CHKERRQ(ierr); 644 #endif 645 PetscFunctionReturn(0); 646 } 647 648 /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */ 649 PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo) 650 { 651 #if defined(PETSC_MISSING_LAPACK_POTRF) 652 PetscFunctionBegin; 653 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable."); 654 #else 655 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 656 PetscErrorCode ierr; 657 PetscBLASInt info,n; 658 659 PetscFunctionBegin; 660 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 661 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 662 if (A->spd) { 663 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 664 PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info)); 665 ierr = PetscFPTrapPop();CHKERRQ(ierr); 666 #if defined(PETSC_USE_COMPLEX) 667 } else if (A->hermitian) { 668 if (!mat->pivots) { 669 ierr = PetscMalloc1(A->rmap->n,&mat->pivots);CHKERRQ(ierr); 670 ierr = PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));CHKERRQ(ierr); 671 } 672 if (!mat->fwork) { 673 PetscScalar dummy; 674 675 mat->lfwork = -1; 676 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 677 PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info)); 678 ierr = PetscFPTrapPop();CHKERRQ(ierr); 679 mat->lfwork = (PetscInt)PetscRealPart(dummy); 680 ierr = PetscMalloc1(mat->lfwork,&mat->fwork);CHKERRQ(ierr); 681 ierr = PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));CHKERRQ(ierr); 682 } 683 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 684 PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info)); 685 ierr = PetscFPTrapPop();CHKERRQ(ierr); 686 #endif 687 } else { /* symmetric case */ 688 if (!mat->pivots) { 689 ierr = PetscMalloc1(A->rmap->n,&mat->pivots);CHKERRQ(ierr); 690 ierr = PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));CHKERRQ(ierr); 691 } 692 if (!mat->fwork) { 693 PetscScalar dummy; 694 695 mat->lfwork = -1; 696 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 697 PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info)); 698 ierr = PetscFPTrapPop();CHKERRQ(ierr); 699 mat->lfwork = (PetscInt)PetscRealPart(dummy); 700 ierr = PetscMalloc1(mat->lfwork,&mat->fwork);CHKERRQ(ierr); 701 ierr = PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));CHKERRQ(ierr); 702 } 703 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 704 PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info)); 705 ierr = PetscFPTrapPop();CHKERRQ(ierr); 706 } 707 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1); 708 709 A->ops->solve = MatSolve_SeqDense; 710 A->ops->matsolve = MatMatSolve_SeqDense; 711 A->ops->solvetranspose = MatSolveTranspose_SeqDense; 712 A->factortype = MAT_FACTOR_CHOLESKY; 713 714 ierr = PetscFree(A->solvertype);CHKERRQ(ierr); 715 ierr = PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);CHKERRQ(ierr); 716 717 ierr = PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);CHKERRQ(ierr); 718 #endif 719 PetscFunctionReturn(0); 720 } 721 722 723 PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy) 724 { 725 PetscErrorCode ierr; 726 MatFactorInfo info; 727 728 PetscFunctionBegin; 729 info.fill = 1.0; 730 731 ierr = MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);CHKERRQ(ierr); 732 ierr = (*fact->ops->choleskyfactor)(fact,0,&info);CHKERRQ(ierr); 733 PetscFunctionReturn(0); 734 } 735 736 PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info) 737 { 738 PetscFunctionBegin; 739 fact->assembled = PETSC_TRUE; 740 fact->preallocated = PETSC_TRUE; 741 fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense; 742 fact->ops->solve = MatSolve_SeqDense; 743 fact->ops->matsolve = MatMatSolve_SeqDense; 744 fact->ops->solvetranspose = MatSolveTranspose_SeqDense; 745 PetscFunctionReturn(0); 746 } 747 748 PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info) 749 { 750 PetscFunctionBegin; 751 fact->preallocated = PETSC_TRUE; 752 fact->assembled = PETSC_TRUE; 753 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense; 754 fact->ops->solve = MatSolve_SeqDense; 755 fact->ops->matsolve = MatMatSolve_SeqDense; 756 fact->ops->solvetranspose = MatSolveTranspose_SeqDense; 757 PetscFunctionReturn(0); 758 } 759 760 /* uses LAPACK */ 761 PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact) 762 { 763 PetscErrorCode ierr; 764 765 PetscFunctionBegin; 766 ierr = MatCreate(PetscObjectComm((PetscObject)A),fact);CHKERRQ(ierr); 767 ierr = MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 768 ierr = MatSetType(*fact,MATDENSE);CHKERRQ(ierr); 769 if (ftype == MAT_FACTOR_LU) { 770 (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense; 771 } else { 772 (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense; 773 } 774 (*fact)->factortype = ftype; 775 776 ierr = PetscFree((*fact)->solvertype);CHKERRQ(ierr); 777 ierr = PetscStrallocpy(MATSOLVERPETSC,&(*fact)->solvertype);CHKERRQ(ierr); 778 PetscFunctionReturn(0); 779 } 780 781 /* ------------------------------------------------------------------*/ 782 static PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx) 783 { 784 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 785 PetscScalar *x,*v = mat->v,zero = 0.0,xt; 786 const PetscScalar *b; 787 PetscErrorCode ierr; 788 PetscInt m = A->rmap->n,i; 789 PetscBLASInt o = 1,bm; 790 791 PetscFunctionBegin; 792 #if defined(PETSC_HAVE_CUDA) 793 if (A->valid_GPU_matrix == PETSC_OFFLOAD_GPU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented"); 794 #endif 795 if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */ 796 ierr = PetscBLASIntCast(m,&bm);CHKERRQ(ierr); 797 if (flag & SOR_ZERO_INITIAL_GUESS) { 798 /* this is a hack fix, should have another version without the second BLASdotu */ 799 ierr = VecSet(xx,zero);CHKERRQ(ierr); 800 } 801 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 802 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 803 its = its*lits; 804 if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 805 while (its--) { 806 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 807 for (i=0; i<m; i++) { 808 PetscStackCallBLAS("BLASdotu",xt = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o)); 809 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 810 } 811 } 812 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 813 for (i=m-1; i>=0; i--) { 814 PetscStackCallBLAS("BLASdotu",xt = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o)); 815 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 816 } 817 } 818 } 819 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 820 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 821 PetscFunctionReturn(0); 822 } 823 824 /* -----------------------------------------------------------------*/ 825 PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy) 826 { 827 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 828 const PetscScalar *v = mat->v,*x; 829 PetscScalar *y; 830 PetscErrorCode ierr; 831 PetscBLASInt m, n,_One=1; 832 PetscScalar _DOne=1.0,_DZero=0.0; 833 834 PetscFunctionBegin; 835 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 836 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 837 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 838 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 839 if (!A->rmap->n || !A->cmap->n) { 840 PetscBLASInt i; 841 for (i=0; i<n; i++) y[i] = 0.0; 842 } else { 843 PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One)); 844 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 845 } 846 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 847 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 848 PetscFunctionReturn(0); 849 } 850 851 PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy) 852 { 853 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 854 PetscScalar *y,_DOne=1.0,_DZero=0.0; 855 PetscErrorCode ierr; 856 PetscBLASInt m, n, _One=1; 857 const PetscScalar *v = mat->v,*x; 858 859 PetscFunctionBegin; 860 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 861 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 862 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 863 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 864 if (!A->rmap->n || !A->cmap->n) { 865 PetscBLASInt i; 866 for (i=0; i<m; i++) y[i] = 0.0; 867 } else { 868 PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One)); 869 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);CHKERRQ(ierr); 870 } 871 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 872 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 873 PetscFunctionReturn(0); 874 } 875 876 PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 877 { 878 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 879 const PetscScalar *v = mat->v,*x; 880 PetscScalar *y,_DOne=1.0; 881 PetscErrorCode ierr; 882 PetscBLASInt m, n, _One=1; 883 884 PetscFunctionBegin; 885 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 886 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 887 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 888 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 889 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 890 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 891 PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One)); 892 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 893 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 894 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n);CHKERRQ(ierr); 895 PetscFunctionReturn(0); 896 } 897 898 PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 899 { 900 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 901 const PetscScalar *v = mat->v,*x; 902 PetscScalar *y; 903 PetscErrorCode ierr; 904 PetscBLASInt m, n, _One=1; 905 PetscScalar _DOne=1.0; 906 907 PetscFunctionBegin; 908 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 909 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 910 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 911 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 912 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 913 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 914 PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One)); 915 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 916 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 917 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n);CHKERRQ(ierr); 918 PetscFunctionReturn(0); 919 } 920 921 /* -----------------------------------------------------------------*/ 922 static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 923 { 924 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 925 PetscErrorCode ierr; 926 PetscInt i; 927 928 PetscFunctionBegin; 929 *ncols = A->cmap->n; 930 if (cols) { 931 ierr = PetscMalloc1(A->cmap->n+1,cols);CHKERRQ(ierr); 932 for (i=0; i<A->cmap->n; i++) (*cols)[i] = i; 933 } 934 if (vals) { 935 const PetscScalar *v; 936 937 ierr = MatDenseGetArrayRead(A,&v);CHKERRQ(ierr); 938 ierr = PetscMalloc1(A->cmap->n+1,vals);CHKERRQ(ierr); 939 v += row; 940 for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;} 941 ierr = MatDenseRestoreArrayRead(A,&v);CHKERRQ(ierr); 942 } 943 PetscFunctionReturn(0); 944 } 945 946 static PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 947 { 948 PetscErrorCode ierr; 949 950 PetscFunctionBegin; 951 if (cols) {ierr = PetscFree(*cols);CHKERRQ(ierr);} 952 if (vals) {ierr = PetscFree(*vals);CHKERRQ(ierr); } 953 PetscFunctionReturn(0); 954 } 955 /* ----------------------------------------------------------------*/ 956 static PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv) 957 { 958 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 959 PetscScalar *av; 960 PetscInt i,j,idx=0; 961 #if defined(PETSC_HAVE_CUDA) 962 PetscOffloadFlag oldf; 963 #endif 964 PetscErrorCode ierr; 965 966 PetscFunctionBegin; 967 ierr = MatDenseGetArray(A,&av);CHKERRQ(ierr); 968 if (!mat->roworiented) { 969 if (addv == INSERT_VALUES) { 970 for (j=0; j<n; j++) { 971 if (indexn[j] < 0) {idx += m; continue;} 972 #if defined(PETSC_USE_DEBUG) 973 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); 974 #endif 975 for (i=0; i<m; i++) { 976 if (indexm[i] < 0) {idx++; continue;} 977 #if defined(PETSC_USE_DEBUG) 978 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); 979 #endif 980 av[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 981 } 982 } 983 } else { 984 for (j=0; j<n; j++) { 985 if (indexn[j] < 0) {idx += m; continue;} 986 #if defined(PETSC_USE_DEBUG) 987 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); 988 #endif 989 for (i=0; i<m; i++) { 990 if (indexm[i] < 0) {idx++; continue;} 991 #if defined(PETSC_USE_DEBUG) 992 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); 993 #endif 994 av[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 995 } 996 } 997 } 998 } else { 999 if (addv == INSERT_VALUES) { 1000 for (i=0; i<m; i++) { 1001 if (indexm[i] < 0) { idx += n; continue;} 1002 #if defined(PETSC_USE_DEBUG) 1003 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); 1004 #endif 1005 for (j=0; j<n; j++) { 1006 if (indexn[j] < 0) { idx++; continue;} 1007 #if defined(PETSC_USE_DEBUG) 1008 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); 1009 #endif 1010 av[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 1011 } 1012 } 1013 } else { 1014 for (i=0; i<m; i++) { 1015 if (indexm[i] < 0) { idx += n; continue;} 1016 #if defined(PETSC_USE_DEBUG) 1017 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); 1018 #endif 1019 for (j=0; j<n; j++) { 1020 if (indexn[j] < 0) { idx++; continue;} 1021 #if defined(PETSC_USE_DEBUG) 1022 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); 1023 #endif 1024 av[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 1025 } 1026 } 1027 } 1028 } 1029 /* hack to prevent unneeded copy to the GPU while returning the array */ 1030 #if defined(PETSC_HAVE_CUDA) 1031 oldf = A->valid_GPU_matrix; 1032 A->valid_GPU_matrix = PETSC_OFFLOAD_GPU; 1033 #endif 1034 ierr = MatDenseRestoreArray(A,&av);CHKERRQ(ierr); 1035 #if defined(PETSC_HAVE_CUDA) 1036 A->valid_GPU_matrix = (oldf == PETSC_OFFLOAD_UNALLOCATED ? PETSC_OFFLOAD_UNALLOCATED : PETSC_OFFLOAD_CPU); 1037 #endif 1038 PetscFunctionReturn(0); 1039 } 1040 1041 static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[]) 1042 { 1043 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1044 const PetscScalar *vv; 1045 PetscInt i,j; 1046 PetscErrorCode ierr; 1047 1048 PetscFunctionBegin; 1049 ierr = MatDenseGetArrayRead(A,&vv);CHKERRQ(ierr); 1050 /* row-oriented output */ 1051 for (i=0; i<m; i++) { 1052 if (indexm[i] < 0) {v += n;continue;} 1053 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); 1054 for (j=0; j<n; j++) { 1055 if (indexn[j] < 0) {v++; continue;} 1056 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); 1057 *v++ = vv[indexn[j]*mat->lda + indexm[i]]; 1058 } 1059 } 1060 ierr = MatDenseRestoreArrayRead(A,&vv);CHKERRQ(ierr); 1061 PetscFunctionReturn(0); 1062 } 1063 1064 /* -----------------------------------------------------------------*/ 1065 1066 static PetscErrorCode MatLoad_SeqDense_Binary(Mat newmat,PetscViewer viewer) 1067 { 1068 Mat_SeqDense *a; 1069 PetscErrorCode ierr; 1070 PetscInt *scols,i,j,nz,header[4]; 1071 int fd; 1072 PetscMPIInt size; 1073 PetscInt *rowlengths = 0,M,N,*cols,grows,gcols; 1074 PetscScalar *vals,*svals,*v,*w; 1075 MPI_Comm comm; 1076 1077 PetscFunctionBegin; 1078 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 1079 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1080 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor"); 1081 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 1082 ierr = PetscBinaryRead(fd,header,4,NULL,PETSC_INT);CHKERRQ(ierr); 1083 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not matrix object"); 1084 M = header[1]; N = header[2]; nz = header[3]; 1085 1086 /* set global size if not set already*/ 1087 if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) { 1088 ierr = MatSetSizes(newmat,M,N,M,N);CHKERRQ(ierr); 1089 } else { 1090 /* if sizes and type are already set, check if the vector global sizes are correct */ 1091 ierr = MatGetSize(newmat,&grows,&gcols);CHKERRQ(ierr); 1092 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); 1093 } 1094 a = (Mat_SeqDense*)newmat->data; 1095 if (!a->user_alloc) { 1096 ierr = MatSeqDenseSetPreallocation(newmat,NULL);CHKERRQ(ierr); 1097 } 1098 1099 if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */ 1100 a = (Mat_SeqDense*)newmat->data; 1101 v = a->v; 1102 /* Allocate some temp space to read in the values and then flip them 1103 from row major to column major */ 1104 ierr = PetscMalloc1(M*N > 0 ? M*N : 1,&w);CHKERRQ(ierr); 1105 /* read in nonzero values */ 1106 ierr = PetscBinaryRead(fd,w,M*N,NULL,PETSC_SCALAR);CHKERRQ(ierr); 1107 /* now flip the values and store them in the matrix*/ 1108 for (j=0; j<N; j++) { 1109 for (i=0; i<M; i++) { 1110 *v++ =w[i*N+j]; 1111 } 1112 } 1113 ierr = PetscFree(w);CHKERRQ(ierr); 1114 } else { 1115 /* read row lengths */ 1116 ierr = PetscMalloc1(M+1,&rowlengths);CHKERRQ(ierr); 1117 ierr = PetscBinaryRead(fd,rowlengths,M,NULL,PETSC_INT);CHKERRQ(ierr); 1118 1119 a = (Mat_SeqDense*)newmat->data; 1120 v = a->v; 1121 1122 /* read column indices and nonzeros */ 1123 ierr = PetscMalloc1(nz+1,&scols);CHKERRQ(ierr); 1124 cols = scols; 1125 ierr = PetscBinaryRead(fd,cols,nz,NULL,PETSC_INT);CHKERRQ(ierr); 1126 ierr = PetscMalloc1(nz+1,&svals);CHKERRQ(ierr); 1127 vals = svals; 1128 ierr = PetscBinaryRead(fd,vals,nz,NULL,PETSC_SCALAR);CHKERRQ(ierr); 1129 1130 /* insert into matrix */ 1131 for (i=0; i<M; i++) { 1132 for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j]; 1133 svals += rowlengths[i]; scols += rowlengths[i]; 1134 } 1135 ierr = PetscFree(vals);CHKERRQ(ierr); 1136 ierr = PetscFree(cols);CHKERRQ(ierr); 1137 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 1138 } 1139 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1140 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1141 PetscFunctionReturn(0); 1142 } 1143 1144 PetscErrorCode MatLoad_SeqDense(Mat newMat, PetscViewer viewer) 1145 { 1146 PetscBool isbinary, ishdf5; 1147 PetscErrorCode ierr; 1148 1149 PetscFunctionBegin; 1150 PetscValidHeaderSpecific(newMat,MAT_CLASSID,1); 1151 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 1152 /* force binary viewer to load .info file if it has not yet done so */ 1153 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 1154 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1155 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5, &ishdf5);CHKERRQ(ierr); 1156 if (isbinary) { 1157 ierr = MatLoad_SeqDense_Binary(newMat,viewer);CHKERRQ(ierr); 1158 } else if (ishdf5) { 1159 #if defined(PETSC_HAVE_HDF5) 1160 ierr = MatLoad_Dense_HDF5(newMat,viewer);CHKERRQ(ierr); 1161 #else 1162 SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5"); 1163 #endif 1164 } else { 1165 SETERRQ2(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"Viewer type %s not yet supported for reading %s matrices",((PetscObject)viewer)->type_name,((PetscObject)newMat)->type_name); 1166 } 1167 PetscFunctionReturn(0); 1168 } 1169 1170 static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer) 1171 { 1172 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1173 PetscErrorCode ierr; 1174 PetscInt i,j; 1175 const char *name; 1176 PetscScalar *v,*av; 1177 PetscViewerFormat format; 1178 #if defined(PETSC_USE_COMPLEX) 1179 PetscBool allreal = PETSC_TRUE; 1180 #endif 1181 1182 PetscFunctionBegin; 1183 ierr = MatDenseGetArrayRead(A,(const PetscScalar**)&av);CHKERRQ(ierr); 1184 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1185 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 1186 PetscFunctionReturn(0); /* do nothing for now */ 1187 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 1188 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 1189 for (i=0; i<A->rmap->n; i++) { 1190 v = av + i; 1191 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 1192 for (j=0; j<A->cmap->n; j++) { 1193 #if defined(PETSC_USE_COMPLEX) 1194 if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) { 1195 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));CHKERRQ(ierr); 1196 } else if (PetscRealPart(*v)) { 1197 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));CHKERRQ(ierr); 1198 } 1199 #else 1200 if (*v) { 1201 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);CHKERRQ(ierr); 1202 } 1203 #endif 1204 v += a->lda; 1205 } 1206 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 1207 } 1208 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 1209 } else { 1210 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 1211 #if defined(PETSC_USE_COMPLEX) 1212 /* determine if matrix has all real values */ 1213 v = av; 1214 for (i=0; i<A->rmap->n*A->cmap->n; i++) { 1215 if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;} 1216 } 1217 #endif 1218 if (format == PETSC_VIEWER_ASCII_MATLAB) { 1219 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 1220 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);CHKERRQ(ierr); 1221 ierr = PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 1222 ierr = PetscViewerASCIIPrintf(viewer,"%s = [\n",name);CHKERRQ(ierr); 1223 } 1224 1225 for (i=0; i<A->rmap->n; i++) { 1226 v = av + i; 1227 for (j=0; j<A->cmap->n; j++) { 1228 #if defined(PETSC_USE_COMPLEX) 1229 if (allreal) { 1230 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));CHKERRQ(ierr); 1231 } else { 1232 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));CHKERRQ(ierr); 1233 } 1234 #else 1235 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);CHKERRQ(ierr); 1236 #endif 1237 v += a->lda; 1238 } 1239 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 1240 } 1241 if (format == PETSC_VIEWER_ASCII_MATLAB) { 1242 ierr = PetscViewerASCIIPrintf(viewer,"];\n");CHKERRQ(ierr); 1243 } 1244 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 1245 } 1246 ierr = MatDenseRestoreArrayRead(A,(const PetscScalar**)&av);CHKERRQ(ierr); 1247 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1248 PetscFunctionReturn(0); 1249 } 1250 1251 static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer) 1252 { 1253 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1254 PetscErrorCode ierr; 1255 int fd; 1256 PetscInt ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n; 1257 PetscScalar *av,*v,*anonz,*vals; 1258 PetscViewerFormat format; 1259 1260 PetscFunctionBegin; 1261 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 1262 ierr = MatDenseGetArrayRead(A,(const PetscScalar**)&av);CHKERRQ(ierr); 1263 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1264 if (format == PETSC_VIEWER_NATIVE) { 1265 /* store the matrix as a dense matrix */ 1266 ierr = PetscMalloc1(4,&col_lens);CHKERRQ(ierr); 1267 1268 col_lens[0] = MAT_FILE_CLASSID; 1269 col_lens[1] = m; 1270 col_lens[2] = n; 1271 col_lens[3] = MATRIX_BINARY_FORMAT_DENSE; 1272 1273 ierr = PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 1274 ierr = PetscFree(col_lens);CHKERRQ(ierr); 1275 1276 /* write out matrix, by rows */ 1277 ierr = PetscMalloc1(m*n+1,&vals);CHKERRQ(ierr); 1278 v = av; 1279 for (j=0; j<n; j++) { 1280 for (i=0; i<m; i++) { 1281 vals[j + i*n] = *v++; 1282 } 1283 } 1284 ierr = PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 1285 ierr = PetscFree(vals);CHKERRQ(ierr); 1286 } else { 1287 ierr = PetscMalloc1(4+nz,&col_lens);CHKERRQ(ierr); 1288 1289 col_lens[0] = MAT_FILE_CLASSID; 1290 col_lens[1] = m; 1291 col_lens[2] = n; 1292 col_lens[3] = nz; 1293 1294 /* store lengths of each row and write (including header) to file */ 1295 for (i=0; i<m; i++) col_lens[4+i] = n; 1296 ierr = PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 1297 1298 /* Possibly should write in smaller increments, not whole matrix at once? */ 1299 /* store column indices (zero start index) */ 1300 ict = 0; 1301 for (i=0; i<m; i++) { 1302 for (j=0; j<n; j++) col_lens[ict++] = j; 1303 } 1304 ierr = PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 1305 ierr = PetscFree(col_lens);CHKERRQ(ierr); 1306 1307 /* store nonzero values */ 1308 ierr = PetscMalloc1(nz+1,&anonz);CHKERRQ(ierr); 1309 ict = 0; 1310 for (i=0; i<m; i++) { 1311 v = av + i; 1312 for (j=0; j<n; j++) { 1313 anonz[ict++] = *v; v += a->lda; 1314 } 1315 } 1316 ierr = PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 1317 ierr = PetscFree(anonz);CHKERRQ(ierr); 1318 } 1319 ierr = MatDenseRestoreArrayRead(A,(const PetscScalar**)&av);CHKERRQ(ierr); 1320 PetscFunctionReturn(0); 1321 } 1322 1323 #include <petscdraw.h> 1324 static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa) 1325 { 1326 Mat A = (Mat) Aa; 1327 PetscErrorCode ierr; 1328 PetscInt m = A->rmap->n,n = A->cmap->n,i,j; 1329 int color = PETSC_DRAW_WHITE; 1330 const PetscScalar *v; 1331 PetscViewer viewer; 1332 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 1333 PetscViewerFormat format; 1334 1335 PetscFunctionBegin; 1336 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 1337 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1338 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 1339 1340 /* Loop over matrix elements drawing boxes */ 1341 ierr = MatDenseGetArrayRead(A,&v);CHKERRQ(ierr); 1342 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 1343 ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr); 1344 /* Blue for negative and Red for positive */ 1345 for (j = 0; j < n; j++) { 1346 x_l = j; x_r = x_l + 1.0; 1347 for (i = 0; i < m; i++) { 1348 y_l = m - i - 1.0; 1349 y_r = y_l + 1.0; 1350 if (PetscRealPart(v[j*m+i]) > 0.) color = PETSC_DRAW_RED; 1351 else if (PetscRealPart(v[j*m+i]) < 0.) color = PETSC_DRAW_BLUE; 1352 else continue; 1353 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1354 } 1355 } 1356 ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr); 1357 } else { 1358 /* use contour shading to indicate magnitude of values */ 1359 /* first determine max of all nonzero values */ 1360 PetscReal minv = 0.0, maxv = 0.0; 1361 PetscDraw popup; 1362 1363 for (i=0; i < m*n; i++) { 1364 if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]); 1365 } 1366 if (minv >= maxv) maxv = minv + PETSC_SMALL; 1367 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 1368 ierr = PetscDrawScalePopup(popup,minv,maxv);CHKERRQ(ierr); 1369 1370 ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr); 1371 for (j=0; j<n; j++) { 1372 x_l = j; 1373 x_r = x_l + 1.0; 1374 for (i=0; i<m; i++) { 1375 y_l = m - i - 1.0; 1376 y_r = y_l + 1.0; 1377 color = PetscDrawRealToColor(PetscAbsScalar(v[j*m+i]),minv,maxv); 1378 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1379 } 1380 } 1381 ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr); 1382 } 1383 ierr = MatDenseRestoreArrayRead(A,&v);CHKERRQ(ierr); 1384 PetscFunctionReturn(0); 1385 } 1386 1387 static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer) 1388 { 1389 PetscDraw draw; 1390 PetscBool isnull; 1391 PetscReal xr,yr,xl,yl,h,w; 1392 PetscErrorCode ierr; 1393 1394 PetscFunctionBegin; 1395 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 1396 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 1397 if (isnull) PetscFunctionReturn(0); 1398 1399 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 1400 xr += w; yr += h; xl = -w; yl = -h; 1401 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 1402 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 1403 ierr = PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);CHKERRQ(ierr); 1404 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);CHKERRQ(ierr); 1405 ierr = PetscDrawSave(draw);CHKERRQ(ierr); 1406 PetscFunctionReturn(0); 1407 } 1408 1409 PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer) 1410 { 1411 PetscErrorCode ierr; 1412 PetscBool iascii,isbinary,isdraw; 1413 1414 PetscFunctionBegin; 1415 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1416 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1417 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 1418 1419 if (iascii) { 1420 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 1421 } else if (isbinary) { 1422 ierr = MatView_SeqDense_Binary(A,viewer);CHKERRQ(ierr); 1423 } else if (isdraw) { 1424 ierr = MatView_SeqDense_Draw(A,viewer);CHKERRQ(ierr); 1425 } 1426 PetscFunctionReturn(0); 1427 } 1428 1429 static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A,const PetscScalar array[]) 1430 { 1431 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1432 1433 PetscFunctionBegin; 1434 a->unplacedarray = a->v; 1435 a->unplaced_user_alloc = a->user_alloc; 1436 a->v = (PetscScalar*) array; 1437 #if defined(PETSC_HAVE_CUDA) 1438 A->valid_GPU_matrix = PETSC_OFFLOAD_CPU; 1439 #endif 1440 PetscFunctionReturn(0); 1441 } 1442 1443 static PetscErrorCode MatDenseResetArray_SeqDense(Mat A) 1444 { 1445 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1446 1447 PetscFunctionBegin; 1448 a->v = a->unplacedarray; 1449 a->user_alloc = a->unplaced_user_alloc; 1450 a->unplacedarray = NULL; 1451 #if defined(PETSC_HAVE_CUDA) 1452 A->valid_GPU_matrix = PETSC_OFFLOAD_CPU; 1453 #endif 1454 PetscFunctionReturn(0); 1455 } 1456 1457 PetscErrorCode MatDestroy_SeqDense(Mat mat) 1458 { 1459 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 1460 PetscErrorCode ierr; 1461 1462 PetscFunctionBegin; 1463 #if defined(PETSC_USE_LOG) 1464 PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n); 1465 #endif 1466 ierr = PetscFree(l->pivots);CHKERRQ(ierr); 1467 ierr = PetscFree(l->fwork);CHKERRQ(ierr); 1468 ierr = MatDestroy(&l->ptapwork);CHKERRQ(ierr); 1469 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 1470 ierr = PetscFree(mat->data);CHKERRQ(ierr); 1471 1472 ierr = PetscObjectChangeTypeName((PetscObject)mat,0);CHKERRQ(ierr); 1473 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetLDA_C",NULL);CHKERRQ(ierr); 1474 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);CHKERRQ(ierr); 1475 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL);CHKERRQ(ierr); 1476 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL);CHKERRQ(ierr); 1477 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);CHKERRQ(ierr); 1478 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);CHKERRQ(ierr); 1479 #if defined(PETSC_HAVE_ELEMENTAL) 1480 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);CHKERRQ(ierr); 1481 #endif 1482 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);CHKERRQ(ierr); 1483 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1484 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1485 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1486 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqbaij_seqdense_C",NULL);CHKERRQ(ierr); 1487 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqbaij_seqdense_C",NULL);CHKERRQ(ierr); 1488 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqbaij_seqdense_C",NULL);CHKERRQ(ierr); 1489 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatPtAP_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1490 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1491 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1492 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1493 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumn_C",NULL);CHKERRQ(ierr); 1494 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumn_C",NULL);CHKERRQ(ierr); 1495 PetscFunctionReturn(0); 1496 } 1497 1498 static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout) 1499 { 1500 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1501 PetscErrorCode ierr; 1502 PetscInt k,j,m,n,M; 1503 PetscScalar *v,tmp; 1504 1505 PetscFunctionBegin; 1506 m = A->rmap->n; M = mat->lda; n = A->cmap->n; 1507 if (reuse == MAT_INPLACE_MATRIX && m == n) { /* in place transpose */ 1508 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 1509 for (j=0; j<m; j++) { 1510 for (k=0; k<j; k++) { 1511 tmp = v[j + k*M]; 1512 v[j + k*M] = v[k + j*M]; 1513 v[k + j*M] = tmp; 1514 } 1515 } 1516 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 1517 } else { /* out-of-place transpose */ 1518 Mat tmat; 1519 Mat_SeqDense *tmatd; 1520 PetscScalar *v2; 1521 PetscInt M2; 1522 1523 if (reuse != MAT_REUSE_MATRIX) { 1524 ierr = MatCreate(PetscObjectComm((PetscObject)A),&tmat);CHKERRQ(ierr); 1525 ierr = MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);CHKERRQ(ierr); 1526 ierr = MatSetType(tmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1527 ierr = MatSeqDenseSetPreallocation(tmat,NULL);CHKERRQ(ierr); 1528 } else tmat = *matout; 1529 1530 ierr = MatDenseGetArrayRead(A,(const PetscScalar**)&v);CHKERRQ(ierr); 1531 ierr = MatDenseGetArray(tmat,&v2);CHKERRQ(ierr); 1532 tmatd = (Mat_SeqDense*)tmat->data; 1533 M2 = tmatd->lda; 1534 for (j=0; j<n; j++) { 1535 for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M]; 1536 } 1537 ierr = MatDenseRestoreArray(tmat,&v2);CHKERRQ(ierr); 1538 ierr = MatDenseRestoreArrayRead(A,(const PetscScalar**)&v);CHKERRQ(ierr); 1539 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1540 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1541 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) *matout = tmat; 1542 else { 1543 ierr = MatHeaderMerge(A,&tmat);CHKERRQ(ierr); 1544 } 1545 } 1546 PetscFunctionReturn(0); 1547 } 1548 1549 static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg) 1550 { 1551 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1552 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1553 PetscInt i; 1554 const PetscScalar *v1,*v2; 1555 PetscErrorCode ierr; 1556 1557 PetscFunctionBegin; 1558 if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1559 if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1560 ierr = MatDenseGetArrayRead(A1,&v1);CHKERRQ(ierr); 1561 ierr = MatDenseGetArrayRead(A2,&v2);CHKERRQ(ierr); 1562 for (i=0; i<A1->cmap->n; i++) { 1563 ierr = PetscArraycmp(v1,v2,A1->rmap->n,flg);CHKERRQ(ierr); 1564 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 1565 v1 += mat1->lda; 1566 v2 += mat2->lda; 1567 } 1568 ierr = MatDenseRestoreArrayRead(A1,&v1);CHKERRQ(ierr); 1569 ierr = MatDenseRestoreArrayRead(A2,&v2);CHKERRQ(ierr); 1570 *flg = PETSC_TRUE; 1571 PetscFunctionReturn(0); 1572 } 1573 1574 static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1575 { 1576 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1577 PetscInt i,n,len; 1578 PetscScalar *x; 1579 const PetscScalar *vv; 1580 PetscErrorCode ierr; 1581 1582 PetscFunctionBegin; 1583 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1584 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1585 len = PetscMin(A->rmap->n,A->cmap->n); 1586 ierr = MatDenseGetArrayRead(A,&vv);CHKERRQ(ierr); 1587 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1588 for (i=0; i<len; i++) { 1589 x[i] = vv[i*mat->lda + i]; 1590 } 1591 ierr = MatDenseRestoreArrayRead(A,&vv);CHKERRQ(ierr); 1592 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1593 PetscFunctionReturn(0); 1594 } 1595 1596 static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1597 { 1598 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1599 const PetscScalar *l,*r; 1600 PetscScalar x,*v,*vv; 1601 PetscErrorCode ierr; 1602 PetscInt i,j,m = A->rmap->n,n = A->cmap->n; 1603 1604 PetscFunctionBegin; 1605 ierr = MatDenseGetArray(A,&vv);CHKERRQ(ierr); 1606 if (ll) { 1607 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1608 ierr = VecGetArrayRead(ll,&l);CHKERRQ(ierr); 1609 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1610 for (i=0; i<m; i++) { 1611 x = l[i]; 1612 v = vv + i; 1613 for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;} 1614 } 1615 ierr = VecRestoreArrayRead(ll,&l);CHKERRQ(ierr); 1616 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1617 } 1618 if (rr) { 1619 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1620 ierr = VecGetArrayRead(rr,&r);CHKERRQ(ierr); 1621 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1622 for (i=0; i<n; i++) { 1623 x = r[i]; 1624 v = vv + i*mat->lda; 1625 for (j=0; j<m; j++) (*v++) *= x; 1626 } 1627 ierr = VecRestoreArrayRead(rr,&r);CHKERRQ(ierr); 1628 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1629 } 1630 ierr = MatDenseRestoreArray(A,&vv);CHKERRQ(ierr); 1631 PetscFunctionReturn(0); 1632 } 1633 1634 PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1635 { 1636 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1637 PetscScalar *v,*vv; 1638 PetscReal sum = 0.0; 1639 PetscInt lda =mat->lda,m=A->rmap->n,i,j; 1640 PetscErrorCode ierr; 1641 1642 PetscFunctionBegin; 1643 ierr = MatDenseGetArrayRead(A,(const PetscScalar**)&vv);CHKERRQ(ierr); 1644 v = vv; 1645 if (type == NORM_FROBENIUS) { 1646 if (lda>m) { 1647 for (j=0; j<A->cmap->n; j++) { 1648 v = vv+j*lda; 1649 for (i=0; i<m; i++) { 1650 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1651 } 1652 } 1653 } else { 1654 #if defined(PETSC_USE_REAL___FP16) 1655 PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n; 1656 *nrm = BLASnrm2_(&cnt,v,&one); 1657 } 1658 #else 1659 for (i=0; i<A->cmap->n*A->rmap->n; i++) { 1660 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1661 } 1662 } 1663 *nrm = PetscSqrtReal(sum); 1664 #endif 1665 ierr = PetscLogFlops(2.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1666 } else if (type == NORM_1) { 1667 *nrm = 0.0; 1668 for (j=0; j<A->cmap->n; j++) { 1669 v = vv + j*mat->lda; 1670 sum = 0.0; 1671 for (i=0; i<A->rmap->n; i++) { 1672 sum += PetscAbsScalar(*v); v++; 1673 } 1674 if (sum > *nrm) *nrm = sum; 1675 } 1676 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1677 } else if (type == NORM_INFINITY) { 1678 *nrm = 0.0; 1679 for (j=0; j<A->rmap->n; j++) { 1680 v = vv + j; 1681 sum = 0.0; 1682 for (i=0; i<A->cmap->n; i++) { 1683 sum += PetscAbsScalar(*v); v += mat->lda; 1684 } 1685 if (sum > *nrm) *nrm = sum; 1686 } 1687 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1688 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 1689 ierr = MatDenseRestoreArrayRead(A,(const PetscScalar**)&vv);CHKERRQ(ierr); 1690 PetscFunctionReturn(0); 1691 } 1692 1693 static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg) 1694 { 1695 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1696 PetscErrorCode ierr; 1697 1698 PetscFunctionBegin; 1699 switch (op) { 1700 case MAT_ROW_ORIENTED: 1701 aij->roworiented = flg; 1702 break; 1703 case MAT_NEW_NONZERO_LOCATIONS: 1704 case MAT_NEW_NONZERO_LOCATION_ERR: 1705 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1706 case MAT_NEW_DIAGONALS: 1707 case MAT_KEEP_NONZERO_PATTERN: 1708 case MAT_IGNORE_OFF_PROC_ENTRIES: 1709 case MAT_USE_HASH_TABLE: 1710 case MAT_IGNORE_ZERO_ENTRIES: 1711 case MAT_IGNORE_LOWER_TRIANGULAR: 1712 case MAT_SORTED_FULL: 1713 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1714 break; 1715 case MAT_SPD: 1716 case MAT_SYMMETRIC: 1717 case MAT_STRUCTURALLY_SYMMETRIC: 1718 case MAT_HERMITIAN: 1719 case MAT_SYMMETRY_ETERNAL: 1720 /* These options are handled directly by MatSetOption() */ 1721 break; 1722 default: 1723 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 1724 } 1725 PetscFunctionReturn(0); 1726 } 1727 1728 static PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1729 { 1730 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1731 PetscErrorCode ierr; 1732 PetscInt lda=l->lda,m=A->rmap->n,j; 1733 PetscScalar *v; 1734 1735 PetscFunctionBegin; 1736 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 1737 if (lda>m) { 1738 for (j=0; j<A->cmap->n; j++) { 1739 ierr = PetscArrayzero(v+j*lda,m);CHKERRQ(ierr); 1740 } 1741 } else { 1742 ierr = PetscArrayzero(v,A->rmap->n*A->cmap->n);CHKERRQ(ierr); 1743 } 1744 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 1745 PetscFunctionReturn(0); 1746 } 1747 1748 static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1749 { 1750 PetscErrorCode ierr; 1751 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1752 PetscInt m = l->lda, n = A->cmap->n, i,j; 1753 PetscScalar *slot,*bb,*v; 1754 const PetscScalar *xx; 1755 1756 PetscFunctionBegin; 1757 #if defined(PETSC_USE_DEBUG) 1758 for (i=0; i<N; i++) { 1759 if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 1760 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); 1761 } 1762 #endif 1763 if (!N) PetscFunctionReturn(0); 1764 1765 /* fix right hand side if needed */ 1766 if (x && b) { 1767 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1768 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1769 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]]; 1770 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1771 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1772 } 1773 1774 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 1775 for (i=0; i<N; i++) { 1776 slot = v + rows[i]; 1777 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 1778 } 1779 if (diag != 0.0) { 1780 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 1781 for (i=0; i<N; i++) { 1782 slot = v + (m+1)*rows[i]; 1783 *slot = diag; 1784 } 1785 } 1786 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 1787 PetscFunctionReturn(0); 1788 } 1789 1790 static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A,PetscInt *lda) 1791 { 1792 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1793 1794 PetscFunctionBegin; 1795 *lda = mat->lda; 1796 PetscFunctionReturn(0); 1797 } 1798 1799 PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[]) 1800 { 1801 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1802 1803 PetscFunctionBegin; 1804 *array = mat->v; 1805 PetscFunctionReturn(0); 1806 } 1807 1808 PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1809 { 1810 PetscFunctionBegin; 1811 PetscFunctionReturn(0); 1812 } 1813 1814 /*@C 1815 MatDenseGetLDA - gets the leading dimension of the array returned from MatDenseGetArray() 1816 1817 Logically Collective on Mat 1818 1819 Input Parameter: 1820 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1821 1822 Output Parameter: 1823 . lda - the leading dimension 1824 1825 Level: intermediate 1826 1827 .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatSeqDenseSetLDA() 1828 @*/ 1829 PetscErrorCode MatDenseGetLDA(Mat A,PetscInt *lda) 1830 { 1831 PetscErrorCode ierr; 1832 1833 PetscFunctionBegin; 1834 ierr = PetscUseMethod(A,"MatDenseGetLDA_C",(Mat,PetscInt*),(A,lda));CHKERRQ(ierr); 1835 PetscFunctionReturn(0); 1836 } 1837 1838 /*@C 1839 MatDenseGetArray - gives access to the array where the data for a SeqDense matrix is stored 1840 1841 Logically Collective on Mat 1842 1843 Input Parameter: 1844 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1845 1846 Output Parameter: 1847 . array - pointer to the data 1848 1849 Level: intermediate 1850 1851 .seealso: MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead() 1852 @*/ 1853 PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array) 1854 { 1855 PetscErrorCode ierr; 1856 1857 PetscFunctionBegin; 1858 ierr = PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1859 PetscFunctionReturn(0); 1860 } 1861 1862 /*@C 1863 MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray() 1864 1865 Logically Collective on Mat 1866 1867 Input Parameters: 1868 + mat - a MATSEQDENSE or MATMPIDENSE matrix 1869 - array - pointer to the data 1870 1871 Level: intermediate 1872 1873 .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead() 1874 @*/ 1875 PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array) 1876 { 1877 PetscErrorCode ierr; 1878 1879 PetscFunctionBegin; 1880 ierr = PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1881 if (array) *array = NULL; 1882 ierr = PetscObjectStateIncrease((PetscObject)A);CHKERRQ(ierr); 1883 PetscFunctionReturn(0); 1884 } 1885 1886 /*@C 1887 MatDenseGetArrayRead - gives access to the array where the data for a SeqDense matrix is stored 1888 1889 Not Collective 1890 1891 Input Parameter: 1892 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1893 1894 Output Parameter: 1895 . array - pointer to the data 1896 1897 Level: intermediate 1898 1899 .seealso: MatDenseRestoreArray(), MatDenseGetArray(), MatDenseRestoreArrayRead() 1900 @*/ 1901 PetscErrorCode MatDenseGetArrayRead(Mat A,const PetscScalar **array) 1902 { 1903 PetscErrorCode ierr; 1904 1905 PetscFunctionBegin; 1906 ierr = PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array));CHKERRQ(ierr); 1907 PetscFunctionReturn(0); 1908 } 1909 1910 /*@C 1911 MatDenseRestoreArrayRead - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray() 1912 1913 Not Collective 1914 1915 Input Parameters: 1916 + mat - a MATSEQDENSE or MATMPIDENSE matrix 1917 - array - pointer to the data 1918 1919 Level: intermediate 1920 1921 .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArray() 1922 @*/ 1923 PetscErrorCode MatDenseRestoreArrayRead(Mat A,const PetscScalar **array) 1924 { 1925 PetscErrorCode ierr; 1926 1927 PetscFunctionBegin; 1928 ierr = PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array));CHKERRQ(ierr); 1929 if (array) *array = NULL; 1930 PetscFunctionReturn(0); 1931 } 1932 1933 static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 1934 { 1935 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1936 PetscErrorCode ierr; 1937 PetscInt i,j,nrows,ncols,blda; 1938 const PetscInt *irow,*icol; 1939 PetscScalar *av,*bv,*v = mat->v; 1940 Mat newmat; 1941 1942 PetscFunctionBegin; 1943 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1944 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1945 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1946 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1947 1948 /* Check submatrixcall */ 1949 if (scall == MAT_REUSE_MATRIX) { 1950 PetscInt n_cols,n_rows; 1951 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1952 if (n_rows != nrows || n_cols != ncols) { 1953 /* resize the result matrix to match number of requested rows/columns */ 1954 ierr = MatSetSizes(*B,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1955 } 1956 newmat = *B; 1957 } else { 1958 /* Create and fill new matrix */ 1959 ierr = MatCreate(PetscObjectComm((PetscObject)A),&newmat);CHKERRQ(ierr); 1960 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1961 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1962 ierr = MatSeqDenseSetPreallocation(newmat,NULL);CHKERRQ(ierr); 1963 } 1964 1965 /* Now extract the data pointers and do the copy,column at a time */ 1966 ierr = MatDenseGetArray(newmat,&bv);CHKERRQ(ierr); 1967 ierr = MatDenseGetLDA(newmat,&blda);CHKERRQ(ierr); 1968 for (i=0; i<ncols; i++) { 1969 av = v + mat->lda*icol[i]; 1970 for (j=0; j<nrows; j++) bv[j] = av[irow[j]]; 1971 bv += blda; 1972 } 1973 ierr = MatDenseRestoreArray(newmat,&bv);CHKERRQ(ierr); 1974 1975 /* Assemble the matrices so that the correct flags are set */ 1976 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1977 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1978 1979 /* Free work space */ 1980 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1981 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1982 *B = newmat; 1983 PetscFunctionReturn(0); 1984 } 1985 1986 static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1987 { 1988 PetscErrorCode ierr; 1989 PetscInt i; 1990 1991 PetscFunctionBegin; 1992 if (scall == MAT_INITIAL_MATRIX) { 1993 ierr = PetscCalloc1(n+1,B);CHKERRQ(ierr); 1994 } 1995 1996 for (i=0; i<n; i++) { 1997 ierr = MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1998 } 1999 PetscFunctionReturn(0); 2000 } 2001 2002 static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 2003 { 2004 PetscFunctionBegin; 2005 PetscFunctionReturn(0); 2006 } 2007 2008 static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 2009 { 2010 PetscFunctionBegin; 2011 PetscFunctionReturn(0); 2012 } 2013 2014 static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 2015 { 2016 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data; 2017 PetscErrorCode ierr; 2018 const PetscScalar *va; 2019 PetscScalar *vb; 2020 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 2021 2022 PetscFunctionBegin; 2023 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 2024 if (A->ops->copy != B->ops->copy) { 2025 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2026 PetscFunctionReturn(0); 2027 } 2028 if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 2029 ierr = MatDenseGetArrayRead(A,&va);CHKERRQ(ierr); 2030 ierr = MatDenseGetArray(B,&vb);CHKERRQ(ierr); 2031 if (lda1>m || lda2>m) { 2032 for (j=0; j<n; j++) { 2033 ierr = PetscArraycpy(vb+j*lda2,va+j*lda1,m);CHKERRQ(ierr); 2034 } 2035 } else { 2036 ierr = PetscArraycpy(vb,va,A->rmap->n*A->cmap->n);CHKERRQ(ierr); 2037 } 2038 ierr = MatDenseRestoreArray(B,&vb);CHKERRQ(ierr); 2039 ierr = MatDenseRestoreArrayRead(A,&va);CHKERRQ(ierr); 2040 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2041 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2042 PetscFunctionReturn(0); 2043 } 2044 2045 static PetscErrorCode MatSetUp_SeqDense(Mat A) 2046 { 2047 PetscErrorCode ierr; 2048 2049 PetscFunctionBegin; 2050 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 2051 PetscFunctionReturn(0); 2052 } 2053 2054 static PetscErrorCode MatConjugate_SeqDense(Mat A) 2055 { 2056 PetscInt i,nz = A->rmap->n*A->cmap->n; 2057 PetscScalar *aa; 2058 PetscErrorCode ierr; 2059 2060 PetscFunctionBegin; 2061 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr); 2062 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 2063 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr); 2064 PetscFunctionReturn(0); 2065 } 2066 2067 static PetscErrorCode MatRealPart_SeqDense(Mat A) 2068 { 2069 PetscInt i,nz = A->rmap->n*A->cmap->n; 2070 PetscScalar *aa; 2071 PetscErrorCode ierr; 2072 2073 PetscFunctionBegin; 2074 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr); 2075 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 2076 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr); 2077 PetscFunctionReturn(0); 2078 } 2079 2080 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 2081 { 2082 PetscInt i,nz = A->rmap->n*A->cmap->n; 2083 PetscScalar *aa; 2084 PetscErrorCode ierr; 2085 2086 PetscFunctionBegin; 2087 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr); 2088 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 2089 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr); 2090 PetscFunctionReturn(0); 2091 } 2092 2093 /* ----------------------------------------------------------------*/ 2094 PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 2095 { 2096 PetscErrorCode ierr; 2097 2098 PetscFunctionBegin; 2099 if (scall == MAT_INITIAL_MATRIX) { 2100 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2101 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 2102 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2103 } 2104 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 2105 if ((*C)->ops->matmultnumeric) { 2106 ierr = (*(*C)->ops->matmultnumeric)(A,B,*C);CHKERRQ(ierr); 2107 } else { 2108 ierr = MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 2109 } 2110 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 2111 PetscFunctionReturn(0); 2112 } 2113 2114 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 2115 { 2116 PetscErrorCode ierr; 2117 PetscInt m=A->rmap->n,n=B->cmap->n; 2118 Mat Cmat; 2119 PetscBool flg; 2120 2121 PetscFunctionBegin; 2122 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 2123 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 2124 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 2125 ierr = MatSetType(Cmat,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 2126 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 2127 *C = Cmat; 2128 PetscFunctionReturn(0); 2129 } 2130 2131 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2132 { 2133 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2134 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2135 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2136 PetscBLASInt m,n,k; 2137 const PetscScalar *av,*bv; 2138 PetscScalar *cv; 2139 PetscScalar _DOne=1.0,_DZero=0.0; 2140 PetscErrorCode ierr; 2141 PetscBool flg; 2142 2143 PetscFunctionBegin; 2144 /* Handle case where where user provided the final C matrix rather than calling MatMatMult() with MAT_INITIAL_MATRIX*/ 2145 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);CHKERRQ(ierr); 2146 if (flg) { 2147 C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqDense; 2148 ierr = (*C->ops->matmultnumeric)(A,B,C);CHKERRQ(ierr); 2149 PetscFunctionReturn(0); 2150 } 2151 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&flg);CHKERRQ(ierr); 2152 if (flg) { 2153 C->ops->matmultnumeric = MatMatMultNumeric_SeqBAIJ_SeqDense; 2154 ierr = (*C->ops->matmultnumeric)(A,B,C);CHKERRQ(ierr); 2155 PetscFunctionReturn(0); 2156 } 2157 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2158 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2159 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 2160 if (!m || !n || !k) PetscFunctionReturn(0); 2161 ierr = MatDenseGetArrayRead(A,&av);CHKERRQ(ierr); 2162 ierr = MatDenseGetArrayRead(B,&bv);CHKERRQ(ierr); 2163 ierr = MatDenseGetArray(C,&cv);CHKERRQ(ierr); 2164 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda)); 2165 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr); 2166 ierr = MatDenseRestoreArrayRead(A,&av);CHKERRQ(ierr); 2167 ierr = MatDenseRestoreArrayRead(B,&bv);CHKERRQ(ierr); 2168 ierr = MatDenseRestoreArray(C,&cv);CHKERRQ(ierr); 2169 PetscFunctionReturn(0); 2170 } 2171 2172 PetscErrorCode MatMatTransposeMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 2173 { 2174 PetscErrorCode ierr; 2175 2176 PetscFunctionBegin; 2177 if (scall == MAT_INITIAL_MATRIX) { 2178 ierr = PetscLogEventBegin(MAT_MatTransposeMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2179 ierr = MatMatTransposeMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 2180 ierr = PetscLogEventEnd(MAT_MatTransposeMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2181 } 2182 ierr = PetscLogEventBegin(MAT_MatTransposeMultNumeric,A,B,0,0);CHKERRQ(ierr); 2183 ierr = MatMatTransposeMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 2184 ierr = PetscLogEventEnd(MAT_MatTransposeMultNumeric,A,B,0,0);CHKERRQ(ierr); 2185 PetscFunctionReturn(0); 2186 } 2187 2188 PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 2189 { 2190 PetscErrorCode ierr; 2191 PetscInt m=A->rmap->n,n=B->rmap->n; 2192 Mat Cmat; 2193 PetscBool flg; 2194 2195 PetscFunctionBegin; 2196 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 2197 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 2198 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 2199 ierr = MatSetType(Cmat,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 2200 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 2201 *C = Cmat; 2202 PetscFunctionReturn(0); 2203 } 2204 2205 PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2206 { 2207 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2208 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2209 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2210 PetscBLASInt m,n,k; 2211 PetscScalar _DOne=1.0,_DZero=0.0; 2212 PetscErrorCode ierr; 2213 2214 PetscFunctionBegin; 2215 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2216 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2217 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 2218 if (!m || !n || !k) PetscFunctionReturn(0); 2219 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 2220 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr); 2221 PetscFunctionReturn(0); 2222 } 2223 2224 PetscErrorCode MatTransposeMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 2225 { 2226 PetscErrorCode ierr; 2227 2228 PetscFunctionBegin; 2229 if (scall == MAT_INITIAL_MATRIX) { 2230 ierr = PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2231 ierr = MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 2232 ierr = PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2233 } 2234 ierr = PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 2235 ierr = MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 2236 ierr = PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 2237 PetscFunctionReturn(0); 2238 } 2239 2240 PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 2241 { 2242 PetscErrorCode ierr; 2243 PetscInt m=A->cmap->n,n=B->cmap->n; 2244 Mat Cmat; 2245 PetscBool flg; 2246 2247 PetscFunctionBegin; 2248 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 2249 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 2250 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 2251 ierr = MatSetType(Cmat,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 2252 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 2253 *C = Cmat; 2254 PetscFunctionReturn(0); 2255 } 2256 2257 PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2258 { 2259 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2260 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2261 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2262 PetscBLASInt m,n,k; 2263 PetscScalar _DOne=1.0,_DZero=0.0; 2264 PetscErrorCode ierr; 2265 2266 PetscFunctionBegin; 2267 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2268 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2269 ierr = PetscBLASIntCast(A->rmap->n,&k);CHKERRQ(ierr); 2270 if (!m || !n || !k) PetscFunctionReturn(0); 2271 PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 2272 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr); 2273 PetscFunctionReturn(0); 2274 } 2275 2276 static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 2277 { 2278 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2279 PetscErrorCode ierr; 2280 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2281 PetscScalar *x; 2282 const PetscScalar *aa; 2283 2284 PetscFunctionBegin; 2285 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2286 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2287 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2288 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2289 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2290 for (i=0; i<m; i++) { 2291 x[i] = aa[i]; if (idx) idx[i] = 0; 2292 for (j=1; j<n; j++) { 2293 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;} 2294 } 2295 } 2296 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2297 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2298 PetscFunctionReturn(0); 2299 } 2300 2301 static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 2302 { 2303 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2304 PetscErrorCode ierr; 2305 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2306 PetscScalar *x; 2307 PetscReal atmp; 2308 const PetscScalar *aa; 2309 2310 PetscFunctionBegin; 2311 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2312 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2313 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2314 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2315 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2316 for (i=0; i<m; i++) { 2317 x[i] = PetscAbsScalar(aa[i]); 2318 for (j=1; j<n; j++) { 2319 atmp = PetscAbsScalar(aa[i+a->lda*j]); 2320 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 2321 } 2322 } 2323 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2324 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2325 PetscFunctionReturn(0); 2326 } 2327 2328 static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 2329 { 2330 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2331 PetscErrorCode ierr; 2332 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2333 PetscScalar *x; 2334 const PetscScalar *aa; 2335 2336 PetscFunctionBegin; 2337 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2338 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2339 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2340 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2341 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2342 for (i=0; i<m; i++) { 2343 x[i] = aa[i]; if (idx) idx[i] = 0; 2344 for (j=1; j<n; j++) { 2345 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;} 2346 } 2347 } 2348 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2349 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2350 PetscFunctionReturn(0); 2351 } 2352 2353 static PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 2354 { 2355 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2356 PetscErrorCode ierr; 2357 PetscScalar *x; 2358 const PetscScalar *aa; 2359 2360 PetscFunctionBegin; 2361 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2362 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2363 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2364 ierr = PetscArraycpy(x,aa+col*a->lda,A->rmap->n);CHKERRQ(ierr); 2365 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2366 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2367 PetscFunctionReturn(0); 2368 } 2369 2370 PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms) 2371 { 2372 PetscErrorCode ierr; 2373 PetscInt i,j,m,n; 2374 const PetscScalar *a; 2375 2376 PetscFunctionBegin; 2377 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 2378 ierr = PetscArrayzero(norms,n);CHKERRQ(ierr); 2379 ierr = MatDenseGetArrayRead(A,&a);CHKERRQ(ierr); 2380 if (type == NORM_2) { 2381 for (i=0; i<n; i++) { 2382 for (j=0; j<m; j++) { 2383 norms[i] += PetscAbsScalar(a[j]*a[j]); 2384 } 2385 a += m; 2386 } 2387 } else if (type == NORM_1) { 2388 for (i=0; i<n; i++) { 2389 for (j=0; j<m; j++) { 2390 norms[i] += PetscAbsScalar(a[j]); 2391 } 2392 a += m; 2393 } 2394 } else if (type == NORM_INFINITY) { 2395 for (i=0; i<n; i++) { 2396 for (j=0; j<m; j++) { 2397 norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]); 2398 } 2399 a += m; 2400 } 2401 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType"); 2402 ierr = MatDenseRestoreArrayRead(A,&a);CHKERRQ(ierr); 2403 if (type == NORM_2) { 2404 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 2405 } 2406 PetscFunctionReturn(0); 2407 } 2408 2409 static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx) 2410 { 2411 PetscErrorCode ierr; 2412 PetscScalar *a; 2413 PetscInt m,n,i; 2414 2415 PetscFunctionBegin; 2416 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 2417 ierr = MatDenseGetArray(x,&a);CHKERRQ(ierr); 2418 for (i=0; i<m*n; i++) { 2419 ierr = PetscRandomGetValue(rctx,a+i);CHKERRQ(ierr); 2420 } 2421 ierr = MatDenseRestoreArray(x,&a);CHKERRQ(ierr); 2422 PetscFunctionReturn(0); 2423 } 2424 2425 static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool *missing,PetscInt *d) 2426 { 2427 PetscFunctionBegin; 2428 *missing = PETSC_FALSE; 2429 PetscFunctionReturn(0); 2430 } 2431 2432 /* vals is not const */ 2433 static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals) 2434 { 2435 PetscErrorCode ierr; 2436 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2437 PetscScalar *v; 2438 2439 PetscFunctionBegin; 2440 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2441 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 2442 *vals = v+col*a->lda; 2443 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 2444 PetscFunctionReturn(0); 2445 } 2446 2447 static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals) 2448 { 2449 PetscFunctionBegin; 2450 *vals = 0; /* user cannot accidently use the array later */ 2451 PetscFunctionReturn(0); 2452 } 2453 2454 /* -------------------------------------------------------------------*/ 2455 static struct _MatOps MatOps_Values = { MatSetValues_SeqDense, 2456 MatGetRow_SeqDense, 2457 MatRestoreRow_SeqDense, 2458 MatMult_SeqDense, 2459 /* 4*/ MatMultAdd_SeqDense, 2460 MatMultTranspose_SeqDense, 2461 MatMultTransposeAdd_SeqDense, 2462 0, 2463 0, 2464 0, 2465 /* 10*/ 0, 2466 MatLUFactor_SeqDense, 2467 MatCholeskyFactor_SeqDense, 2468 MatSOR_SeqDense, 2469 MatTranspose_SeqDense, 2470 /* 15*/ MatGetInfo_SeqDense, 2471 MatEqual_SeqDense, 2472 MatGetDiagonal_SeqDense, 2473 MatDiagonalScale_SeqDense, 2474 MatNorm_SeqDense, 2475 /* 20*/ MatAssemblyBegin_SeqDense, 2476 MatAssemblyEnd_SeqDense, 2477 MatSetOption_SeqDense, 2478 MatZeroEntries_SeqDense, 2479 /* 24*/ MatZeroRows_SeqDense, 2480 0, 2481 0, 2482 0, 2483 0, 2484 /* 29*/ MatSetUp_SeqDense, 2485 0, 2486 0, 2487 0, 2488 0, 2489 /* 34*/ MatDuplicate_SeqDense, 2490 0, 2491 0, 2492 0, 2493 0, 2494 /* 39*/ MatAXPY_SeqDense, 2495 MatCreateSubMatrices_SeqDense, 2496 0, 2497 MatGetValues_SeqDense, 2498 MatCopy_SeqDense, 2499 /* 44*/ MatGetRowMax_SeqDense, 2500 MatScale_SeqDense, 2501 MatShift_Basic, 2502 0, 2503 MatZeroRowsColumns_SeqDense, 2504 /* 49*/ MatSetRandom_SeqDense, 2505 0, 2506 0, 2507 0, 2508 0, 2509 /* 54*/ 0, 2510 0, 2511 0, 2512 0, 2513 0, 2514 /* 59*/ 0, 2515 MatDestroy_SeqDense, 2516 MatView_SeqDense, 2517 0, 2518 0, 2519 /* 64*/ 0, 2520 0, 2521 0, 2522 0, 2523 0, 2524 /* 69*/ MatGetRowMaxAbs_SeqDense, 2525 0, 2526 0, 2527 0, 2528 0, 2529 /* 74*/ 0, 2530 0, 2531 0, 2532 0, 2533 0, 2534 /* 79*/ 0, 2535 0, 2536 0, 2537 0, 2538 /* 83*/ MatLoad_SeqDense, 2539 0, 2540 MatIsHermitian_SeqDense, 2541 0, 2542 0, 2543 0, 2544 /* 89*/ MatMatMult_SeqDense_SeqDense, 2545 MatMatMultSymbolic_SeqDense_SeqDense, 2546 MatMatMultNumeric_SeqDense_SeqDense, 2547 MatPtAP_SeqDense_SeqDense, 2548 MatPtAPSymbolic_SeqDense_SeqDense, 2549 /* 94*/ MatPtAPNumeric_SeqDense_SeqDense, 2550 MatMatTransposeMult_SeqDense_SeqDense, 2551 MatMatTransposeMultSymbolic_SeqDense_SeqDense, 2552 MatMatTransposeMultNumeric_SeqDense_SeqDense, 2553 0, 2554 /* 99*/ 0, 2555 0, 2556 0, 2557 MatConjugate_SeqDense, 2558 0, 2559 /*104*/ 0, 2560 MatRealPart_SeqDense, 2561 MatImaginaryPart_SeqDense, 2562 0, 2563 0, 2564 /*109*/ 0, 2565 0, 2566 MatGetRowMin_SeqDense, 2567 MatGetColumnVector_SeqDense, 2568 MatMissingDiagonal_SeqDense, 2569 /*114*/ 0, 2570 0, 2571 0, 2572 0, 2573 0, 2574 /*119*/ 0, 2575 0, 2576 0, 2577 0, 2578 0, 2579 /*124*/ 0, 2580 MatGetColumnNorms_SeqDense, 2581 0, 2582 0, 2583 0, 2584 /*129*/ 0, 2585 MatTransposeMatMult_SeqDense_SeqDense, 2586 MatTransposeMatMultSymbolic_SeqDense_SeqDense, 2587 MatTransposeMatMultNumeric_SeqDense_SeqDense, 2588 0, 2589 /*134*/ 0, 2590 0, 2591 0, 2592 0, 2593 0, 2594 /*139*/ 0, 2595 0, 2596 0, 2597 0, 2598 0, 2599 /*144*/ MatCreateMPIMatConcatenateSeqMat_SeqDense 2600 }; 2601 2602 /*@C 2603 MatCreateSeqDense - Creates a sequential dense matrix that 2604 is stored in column major order (the usual Fortran 77 manner). Many 2605 of the matrix operations use the BLAS and LAPACK routines. 2606 2607 Collective 2608 2609 Input Parameters: 2610 + comm - MPI communicator, set to PETSC_COMM_SELF 2611 . m - number of rows 2612 . n - number of columns 2613 - data - optional location of matrix data in column major order. Set data=NULL for PETSc 2614 to control all matrix memory allocation. 2615 2616 Output Parameter: 2617 . A - the matrix 2618 2619 Notes: 2620 The data input variable is intended primarily for Fortran programmers 2621 who wish to allocate their own matrix memory space. Most users should 2622 set data=NULL. 2623 2624 Level: intermediate 2625 2626 .seealso: MatCreate(), MatCreateDense(), MatSetValues() 2627 @*/ 2628 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 2629 { 2630 PetscErrorCode ierr; 2631 2632 PetscFunctionBegin; 2633 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2634 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2635 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 2636 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 2637 PetscFunctionReturn(0); 2638 } 2639 2640 /*@C 2641 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 2642 2643 Collective 2644 2645 Input Parameters: 2646 + B - the matrix 2647 - data - the array (or NULL) 2648 2649 Notes: 2650 The data input variable is intended primarily for Fortran programmers 2651 who wish to allocate their own matrix memory space. Most users should 2652 need not call this routine. 2653 2654 Level: intermediate 2655 2656 .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA() 2657 2658 @*/ 2659 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 2660 { 2661 PetscErrorCode ierr; 2662 2663 PetscFunctionBegin; 2664 ierr = PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));CHKERRQ(ierr); 2665 PetscFunctionReturn(0); 2666 } 2667 2668 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 2669 { 2670 Mat_SeqDense *b; 2671 PetscErrorCode ierr; 2672 2673 PetscFunctionBegin; 2674 B->preallocated = PETSC_TRUE; 2675 2676 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2677 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2678 2679 b = (Mat_SeqDense*)B->data; 2680 b->Mmax = B->rmap->n; 2681 b->Nmax = B->cmap->n; 2682 if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n; 2683 2684 ierr = PetscIntMultError(b->lda,b->Nmax,NULL);CHKERRQ(ierr); 2685 if (!data) { /* petsc-allocated storage */ 2686 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2687 ierr = PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);CHKERRQ(ierr); 2688 ierr = PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 2689 2690 b->user_alloc = PETSC_FALSE; 2691 } else { /* user-allocated storage */ 2692 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2693 b->v = data; 2694 b->user_alloc = PETSC_TRUE; 2695 } 2696 B->assembled = PETSC_TRUE; 2697 PetscFunctionReturn(0); 2698 } 2699 2700 #if defined(PETSC_HAVE_ELEMENTAL) 2701 PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 2702 { 2703 Mat mat_elemental; 2704 PetscErrorCode ierr; 2705 const PetscScalar *array; 2706 PetscScalar *v_colwise; 2707 PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols; 2708 2709 PetscFunctionBegin; 2710 ierr = PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);CHKERRQ(ierr); 2711 ierr = MatDenseGetArrayRead(A,&array);CHKERRQ(ierr); 2712 /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */ 2713 k = 0; 2714 for (j=0; j<N; j++) { 2715 cols[j] = j; 2716 for (i=0; i<M; i++) { 2717 v_colwise[j*M+i] = array[k++]; 2718 } 2719 } 2720 for (i=0; i<M; i++) { 2721 rows[i] = i; 2722 } 2723 ierr = MatDenseRestoreArrayRead(A,&array);CHKERRQ(ierr); 2724 2725 ierr = MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);CHKERRQ(ierr); 2726 ierr = MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 2727 ierr = MatSetType(mat_elemental,MATELEMENTAL);CHKERRQ(ierr); 2728 ierr = MatSetUp(mat_elemental);CHKERRQ(ierr); 2729 2730 /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */ 2731 ierr = MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);CHKERRQ(ierr); 2732 ierr = MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2733 ierr = MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2734 ierr = PetscFree3(v_colwise,rows,cols);CHKERRQ(ierr); 2735 2736 if (reuse == MAT_INPLACE_MATRIX) { 2737 ierr = MatHeaderReplace(A,&mat_elemental);CHKERRQ(ierr); 2738 } else { 2739 *newmat = mat_elemental; 2740 } 2741 PetscFunctionReturn(0); 2742 } 2743 #endif 2744 2745 /*@C 2746 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 2747 2748 Input parameter: 2749 + A - the matrix 2750 - lda - the leading dimension 2751 2752 Notes: 2753 This routine is to be used in conjunction with MatSeqDenseSetPreallocation(); 2754 it asserts that the preallocation has a leading dimension (the LDA parameter 2755 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 2756 2757 Level: intermediate 2758 2759 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 2760 2761 @*/ 2762 PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda) 2763 { 2764 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2765 2766 PetscFunctionBegin; 2767 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); 2768 b->lda = lda; 2769 b->changelda = PETSC_FALSE; 2770 b->Mmax = PetscMax(b->Mmax,lda); 2771 PetscFunctionReturn(0); 2772 } 2773 2774 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 2775 { 2776 PetscErrorCode ierr; 2777 PetscMPIInt size; 2778 2779 PetscFunctionBegin; 2780 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2781 if (size == 1) { 2782 if (scall == MAT_INITIAL_MATRIX) { 2783 ierr = MatDuplicate(inmat,MAT_COPY_VALUES,outmat);CHKERRQ(ierr); 2784 } else { 2785 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2786 } 2787 } else { 2788 ierr = MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 2789 } 2790 PetscFunctionReturn(0); 2791 } 2792 2793 /*MC 2794 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 2795 2796 Options Database Keys: 2797 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 2798 2799 Level: beginner 2800 2801 .seealso: MatCreateSeqDense() 2802 2803 M*/ 2804 PetscErrorCode MatCreate_SeqDense(Mat B) 2805 { 2806 Mat_SeqDense *b; 2807 PetscErrorCode ierr; 2808 PetscMPIInt size; 2809 2810 PetscFunctionBegin; 2811 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 2812 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 2813 2814 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 2815 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2816 B->data = (void*)b; 2817 2818 b->roworiented = PETSC_TRUE; 2819 2820 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetLDA_C",MatDenseGetLDA_SeqDense);CHKERRQ(ierr); 2821 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2822 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2823 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);CHKERRQ(ierr); 2824 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);CHKERRQ(ierr); 2825 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2826 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2827 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);CHKERRQ(ierr); 2828 #if defined(PETSC_HAVE_ELEMENTAL) 2829 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);CHKERRQ(ierr); 2830 #endif 2831 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 2832 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2833 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2834 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2835 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqbaij_seqdense_C",MatMatMult_SeqBAIJ_SeqDense);CHKERRQ(ierr); 2836 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqbaij_seqdense_C",MatMatMultSymbolic_SeqBAIJ_SeqDense);CHKERRQ(ierr); 2837 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqbaij_seqdense_C",MatMatMultNumeric_SeqBAIJ_SeqDense);CHKERRQ(ierr); 2838 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaij_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2839 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijperm_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2840 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijperm_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2841 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijperm_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2842 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijperm_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2843 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijsell_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2844 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijsell_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2845 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijsell_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2846 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijsell_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2847 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijmkl_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2848 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijmkl_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2849 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijmkl_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2850 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijmkl_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2851 2852 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2853 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2854 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2855 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijperm_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2856 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijperm_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2857 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijperm_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2858 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijsell_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2859 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijsell_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2860 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijsell_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2861 2862 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijmkl_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2863 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijmkl_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2864 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijmkl_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2865 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);CHKERRQ(ierr); 2866 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);CHKERRQ(ierr); 2867 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 2868 PetscFunctionReturn(0); 2869 } 2870 2871 /*@C 2872 MatDenseGetColumn - gives access to a column of a dense matrix. This is only the local part of the column. You MUST call MatDenseRestoreColumn() to avoid memory bleeding. 2873 2874 Not Collective 2875 2876 Input Parameter: 2877 + mat - a MATSEQDENSE or MATMPIDENSE matrix 2878 - col - column index 2879 2880 Output Parameter: 2881 . vals - pointer to the data 2882 2883 Level: intermediate 2884 2885 .seealso: MatDenseRestoreColumn() 2886 @*/ 2887 PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals) 2888 { 2889 PetscErrorCode ierr; 2890 2891 PetscFunctionBegin; 2892 ierr = PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));CHKERRQ(ierr); 2893 PetscFunctionReturn(0); 2894 } 2895 2896 /*@C 2897 MatDenseRestoreColumn - returns access to a column of a dense matrix which is returned by MatDenseGetColumn(). 2898 2899 Not Collective 2900 2901 Input Parameter: 2902 . mat - a MATSEQDENSE or MATMPIDENSE matrix 2903 2904 Output Parameter: 2905 . vals - pointer to the data 2906 2907 Level: intermediate 2908 2909 .seealso: MatDenseGetColumn() 2910 @*/ 2911 PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals) 2912 { 2913 PetscErrorCode ierr; 2914 2915 PetscFunctionBegin; 2916 ierr = PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));CHKERRQ(ierr); 2917 PetscFunctionReturn(0); 2918 } 2919