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 = VecGetArrayWrite(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 = VecRestoreArrayWrite(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 = VecGetArrayWrite(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 = VecRestoreArrayWrite(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 #if defined(PETSC_HAVE_CUDA) 1483 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqdensecuda_C",NULL);CHKERRQ(ierr); 1484 #endif 1485 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);CHKERRQ(ierr); 1486 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1487 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1488 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1489 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqbaij_seqdense_C",NULL);CHKERRQ(ierr); 1490 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqbaij_seqdense_C",NULL);CHKERRQ(ierr); 1491 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqbaij_seqdense_C",NULL);CHKERRQ(ierr); 1492 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatPtAP_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1493 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1494 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1495 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1496 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumn_C",NULL);CHKERRQ(ierr); 1497 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumn_C",NULL);CHKERRQ(ierr); 1498 PetscFunctionReturn(0); 1499 } 1500 1501 static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout) 1502 { 1503 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1504 PetscErrorCode ierr; 1505 PetscInt k,j,m,n,M; 1506 PetscScalar *v,tmp; 1507 1508 PetscFunctionBegin; 1509 m = A->rmap->n; M = mat->lda; n = A->cmap->n; 1510 if (reuse == MAT_INPLACE_MATRIX && m == n) { /* in place transpose */ 1511 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 1512 for (j=0; j<m; j++) { 1513 for (k=0; k<j; k++) { 1514 tmp = v[j + k*M]; 1515 v[j + k*M] = v[k + j*M]; 1516 v[k + j*M] = tmp; 1517 } 1518 } 1519 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 1520 } else { /* out-of-place transpose */ 1521 Mat tmat; 1522 Mat_SeqDense *tmatd; 1523 PetscScalar *v2; 1524 PetscInt M2; 1525 1526 if (reuse != MAT_REUSE_MATRIX) { 1527 ierr = MatCreate(PetscObjectComm((PetscObject)A),&tmat);CHKERRQ(ierr); 1528 ierr = MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);CHKERRQ(ierr); 1529 ierr = MatSetType(tmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1530 ierr = MatSeqDenseSetPreallocation(tmat,NULL);CHKERRQ(ierr); 1531 } else tmat = *matout; 1532 1533 ierr = MatDenseGetArrayRead(A,(const PetscScalar**)&v);CHKERRQ(ierr); 1534 ierr = MatDenseGetArray(tmat,&v2);CHKERRQ(ierr); 1535 tmatd = (Mat_SeqDense*)tmat->data; 1536 M2 = tmatd->lda; 1537 for (j=0; j<n; j++) { 1538 for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M]; 1539 } 1540 ierr = MatDenseRestoreArray(tmat,&v2);CHKERRQ(ierr); 1541 ierr = MatDenseRestoreArrayRead(A,(const PetscScalar**)&v);CHKERRQ(ierr); 1542 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1543 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1544 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) *matout = tmat; 1545 else { 1546 ierr = MatHeaderMerge(A,&tmat);CHKERRQ(ierr); 1547 } 1548 } 1549 PetscFunctionReturn(0); 1550 } 1551 1552 static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg) 1553 { 1554 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1555 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1556 PetscInt i; 1557 const PetscScalar *v1,*v2; 1558 PetscErrorCode ierr; 1559 1560 PetscFunctionBegin; 1561 if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1562 if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1563 ierr = MatDenseGetArrayRead(A1,&v1);CHKERRQ(ierr); 1564 ierr = MatDenseGetArrayRead(A2,&v2);CHKERRQ(ierr); 1565 for (i=0; i<A1->cmap->n; i++) { 1566 ierr = PetscArraycmp(v1,v2,A1->rmap->n,flg);CHKERRQ(ierr); 1567 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 1568 v1 += mat1->lda; 1569 v2 += mat2->lda; 1570 } 1571 ierr = MatDenseRestoreArrayRead(A1,&v1);CHKERRQ(ierr); 1572 ierr = MatDenseRestoreArrayRead(A2,&v2);CHKERRQ(ierr); 1573 *flg = PETSC_TRUE; 1574 PetscFunctionReturn(0); 1575 } 1576 1577 static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1578 { 1579 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1580 PetscInt i,n,len; 1581 PetscScalar *x; 1582 const PetscScalar *vv; 1583 PetscErrorCode ierr; 1584 1585 PetscFunctionBegin; 1586 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1587 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1588 len = PetscMin(A->rmap->n,A->cmap->n); 1589 ierr = MatDenseGetArrayRead(A,&vv);CHKERRQ(ierr); 1590 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1591 for (i=0; i<len; i++) { 1592 x[i] = vv[i*mat->lda + i]; 1593 } 1594 ierr = MatDenseRestoreArrayRead(A,&vv);CHKERRQ(ierr); 1595 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1596 PetscFunctionReturn(0); 1597 } 1598 1599 static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1600 { 1601 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1602 const PetscScalar *l,*r; 1603 PetscScalar x,*v,*vv; 1604 PetscErrorCode ierr; 1605 PetscInt i,j,m = A->rmap->n,n = A->cmap->n; 1606 1607 PetscFunctionBegin; 1608 ierr = MatDenseGetArray(A,&vv);CHKERRQ(ierr); 1609 if (ll) { 1610 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1611 ierr = VecGetArrayRead(ll,&l);CHKERRQ(ierr); 1612 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1613 for (i=0; i<m; i++) { 1614 x = l[i]; 1615 v = vv + i; 1616 for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;} 1617 } 1618 ierr = VecRestoreArrayRead(ll,&l);CHKERRQ(ierr); 1619 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1620 } 1621 if (rr) { 1622 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1623 ierr = VecGetArrayRead(rr,&r);CHKERRQ(ierr); 1624 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1625 for (i=0; i<n; i++) { 1626 x = r[i]; 1627 v = vv + i*mat->lda; 1628 for (j=0; j<m; j++) (*v++) *= x; 1629 } 1630 ierr = VecRestoreArrayRead(rr,&r);CHKERRQ(ierr); 1631 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1632 } 1633 ierr = MatDenseRestoreArray(A,&vv);CHKERRQ(ierr); 1634 PetscFunctionReturn(0); 1635 } 1636 1637 PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1638 { 1639 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1640 PetscScalar *v,*vv; 1641 PetscReal sum = 0.0; 1642 PetscInt lda =mat->lda,m=A->rmap->n,i,j; 1643 PetscErrorCode ierr; 1644 1645 PetscFunctionBegin; 1646 ierr = MatDenseGetArrayRead(A,(const PetscScalar**)&vv);CHKERRQ(ierr); 1647 v = vv; 1648 if (type == NORM_FROBENIUS) { 1649 if (lda>m) { 1650 for (j=0; j<A->cmap->n; j++) { 1651 v = vv+j*lda; 1652 for (i=0; i<m; i++) { 1653 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1654 } 1655 } 1656 } else { 1657 #if defined(PETSC_USE_REAL___FP16) 1658 PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n; 1659 *nrm = BLASnrm2_(&cnt,v,&one); 1660 } 1661 #else 1662 for (i=0; i<A->cmap->n*A->rmap->n; i++) { 1663 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1664 } 1665 } 1666 *nrm = PetscSqrtReal(sum); 1667 #endif 1668 ierr = PetscLogFlops(2.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1669 } else if (type == NORM_1) { 1670 *nrm = 0.0; 1671 for (j=0; j<A->cmap->n; j++) { 1672 v = vv + j*mat->lda; 1673 sum = 0.0; 1674 for (i=0; i<A->rmap->n; i++) { 1675 sum += PetscAbsScalar(*v); v++; 1676 } 1677 if (sum > *nrm) *nrm = sum; 1678 } 1679 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1680 } else if (type == NORM_INFINITY) { 1681 *nrm = 0.0; 1682 for (j=0; j<A->rmap->n; j++) { 1683 v = vv + j; 1684 sum = 0.0; 1685 for (i=0; i<A->cmap->n; i++) { 1686 sum += PetscAbsScalar(*v); v += mat->lda; 1687 } 1688 if (sum > *nrm) *nrm = sum; 1689 } 1690 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1691 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 1692 ierr = MatDenseRestoreArrayRead(A,(const PetscScalar**)&vv);CHKERRQ(ierr); 1693 PetscFunctionReturn(0); 1694 } 1695 1696 static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg) 1697 { 1698 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1699 PetscErrorCode ierr; 1700 1701 PetscFunctionBegin; 1702 switch (op) { 1703 case MAT_ROW_ORIENTED: 1704 aij->roworiented = flg; 1705 break; 1706 case MAT_NEW_NONZERO_LOCATIONS: 1707 case MAT_NEW_NONZERO_LOCATION_ERR: 1708 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1709 case MAT_NEW_DIAGONALS: 1710 case MAT_KEEP_NONZERO_PATTERN: 1711 case MAT_IGNORE_OFF_PROC_ENTRIES: 1712 case MAT_USE_HASH_TABLE: 1713 case MAT_IGNORE_ZERO_ENTRIES: 1714 case MAT_IGNORE_LOWER_TRIANGULAR: 1715 case MAT_SORTED_FULL: 1716 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1717 break; 1718 case MAT_SPD: 1719 case MAT_SYMMETRIC: 1720 case MAT_STRUCTURALLY_SYMMETRIC: 1721 case MAT_HERMITIAN: 1722 case MAT_SYMMETRY_ETERNAL: 1723 /* These options are handled directly by MatSetOption() */ 1724 break; 1725 default: 1726 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 1727 } 1728 PetscFunctionReturn(0); 1729 } 1730 1731 static PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1732 { 1733 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1734 PetscErrorCode ierr; 1735 PetscInt lda=l->lda,m=A->rmap->n,j; 1736 PetscScalar *v; 1737 1738 PetscFunctionBegin; 1739 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 1740 if (lda>m) { 1741 for (j=0; j<A->cmap->n; j++) { 1742 ierr = PetscArrayzero(v+j*lda,m);CHKERRQ(ierr); 1743 } 1744 } else { 1745 ierr = PetscArrayzero(v,A->rmap->n*A->cmap->n);CHKERRQ(ierr); 1746 } 1747 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 1748 PetscFunctionReturn(0); 1749 } 1750 1751 static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1752 { 1753 PetscErrorCode ierr; 1754 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1755 PetscInt m = l->lda, n = A->cmap->n, i,j; 1756 PetscScalar *slot,*bb,*v; 1757 const PetscScalar *xx; 1758 1759 PetscFunctionBegin; 1760 #if defined(PETSC_USE_DEBUG) 1761 for (i=0; i<N; i++) { 1762 if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 1763 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); 1764 } 1765 #endif 1766 if (!N) PetscFunctionReturn(0); 1767 1768 /* fix right hand side if needed */ 1769 if (x && b) { 1770 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1771 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1772 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]]; 1773 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1774 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1775 } 1776 1777 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 1778 for (i=0; i<N; i++) { 1779 slot = v + rows[i]; 1780 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 1781 } 1782 if (diag != 0.0) { 1783 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 1784 for (i=0; i<N; i++) { 1785 slot = v + (m+1)*rows[i]; 1786 *slot = diag; 1787 } 1788 } 1789 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 1790 PetscFunctionReturn(0); 1791 } 1792 1793 static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A,PetscInt *lda) 1794 { 1795 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1796 1797 PetscFunctionBegin; 1798 *lda = mat->lda; 1799 PetscFunctionReturn(0); 1800 } 1801 1802 PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[]) 1803 { 1804 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1805 1806 PetscFunctionBegin; 1807 *array = mat->v; 1808 PetscFunctionReturn(0); 1809 } 1810 1811 PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1812 { 1813 PetscFunctionBegin; 1814 PetscFunctionReturn(0); 1815 } 1816 1817 /*@C 1818 MatDenseGetLDA - gets the leading dimension of the array returned from MatDenseGetArray() 1819 1820 Logically Collective on Mat 1821 1822 Input Parameter: 1823 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1824 1825 Output Parameter: 1826 . lda - the leading dimension 1827 1828 Level: intermediate 1829 1830 .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatSeqDenseSetLDA() 1831 @*/ 1832 PetscErrorCode MatDenseGetLDA(Mat A,PetscInt *lda) 1833 { 1834 PetscErrorCode ierr; 1835 1836 PetscFunctionBegin; 1837 ierr = PetscUseMethod(A,"MatDenseGetLDA_C",(Mat,PetscInt*),(A,lda));CHKERRQ(ierr); 1838 PetscFunctionReturn(0); 1839 } 1840 1841 /*@C 1842 MatDenseGetArray - gives access to the array where the data for a SeqDense matrix is stored 1843 1844 Logically Collective on Mat 1845 1846 Input Parameter: 1847 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1848 1849 Output Parameter: 1850 . array - pointer to the data 1851 1852 Level: intermediate 1853 1854 .seealso: MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead() 1855 @*/ 1856 PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array) 1857 { 1858 PetscErrorCode ierr; 1859 1860 PetscFunctionBegin; 1861 ierr = PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1862 PetscFunctionReturn(0); 1863 } 1864 1865 /*@C 1866 MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray() 1867 1868 Logically Collective on Mat 1869 1870 Input Parameters: 1871 + mat - a MATSEQDENSE or MATMPIDENSE matrix 1872 - array - pointer to the data 1873 1874 Level: intermediate 1875 1876 .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead() 1877 @*/ 1878 PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array) 1879 { 1880 PetscErrorCode ierr; 1881 1882 PetscFunctionBegin; 1883 ierr = PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1884 if (array) *array = NULL; 1885 ierr = PetscObjectStateIncrease((PetscObject)A);CHKERRQ(ierr); 1886 PetscFunctionReturn(0); 1887 } 1888 1889 /*@C 1890 MatDenseGetArrayRead - gives access to the array where the data for a SeqDense matrix is stored 1891 1892 Not Collective 1893 1894 Input Parameter: 1895 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1896 1897 Output Parameter: 1898 . array - pointer to the data 1899 1900 Level: intermediate 1901 1902 .seealso: MatDenseRestoreArray(), MatDenseGetArray(), MatDenseRestoreArrayRead() 1903 @*/ 1904 PetscErrorCode MatDenseGetArrayRead(Mat A,const PetscScalar **array) 1905 { 1906 PetscErrorCode ierr; 1907 1908 PetscFunctionBegin; 1909 ierr = PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array));CHKERRQ(ierr); 1910 PetscFunctionReturn(0); 1911 } 1912 1913 /*@C 1914 MatDenseRestoreArrayRead - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray() 1915 1916 Not Collective 1917 1918 Input Parameters: 1919 + mat - a MATSEQDENSE or MATMPIDENSE matrix 1920 - array - pointer to the data 1921 1922 Level: intermediate 1923 1924 .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArray() 1925 @*/ 1926 PetscErrorCode MatDenseRestoreArrayRead(Mat A,const PetscScalar **array) 1927 { 1928 PetscErrorCode ierr; 1929 1930 PetscFunctionBegin; 1931 ierr = PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array));CHKERRQ(ierr); 1932 if (array) *array = NULL; 1933 PetscFunctionReturn(0); 1934 } 1935 1936 static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 1937 { 1938 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1939 PetscErrorCode ierr; 1940 PetscInt i,j,nrows,ncols,blda; 1941 const PetscInt *irow,*icol; 1942 PetscScalar *av,*bv,*v = mat->v; 1943 Mat newmat; 1944 1945 PetscFunctionBegin; 1946 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1947 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1948 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1949 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1950 1951 /* Check submatrixcall */ 1952 if (scall == MAT_REUSE_MATRIX) { 1953 PetscInt n_cols,n_rows; 1954 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1955 if (n_rows != nrows || n_cols != ncols) { 1956 /* resize the result matrix to match number of requested rows/columns */ 1957 ierr = MatSetSizes(*B,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1958 } 1959 newmat = *B; 1960 } else { 1961 /* Create and fill new matrix */ 1962 ierr = MatCreate(PetscObjectComm((PetscObject)A),&newmat);CHKERRQ(ierr); 1963 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1964 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1965 ierr = MatSeqDenseSetPreallocation(newmat,NULL);CHKERRQ(ierr); 1966 } 1967 1968 /* Now extract the data pointers and do the copy,column at a time */ 1969 ierr = MatDenseGetArray(newmat,&bv);CHKERRQ(ierr); 1970 ierr = MatDenseGetLDA(newmat,&blda);CHKERRQ(ierr); 1971 for (i=0; i<ncols; i++) { 1972 av = v + mat->lda*icol[i]; 1973 for (j=0; j<nrows; j++) bv[j] = av[irow[j]]; 1974 bv += blda; 1975 } 1976 ierr = MatDenseRestoreArray(newmat,&bv);CHKERRQ(ierr); 1977 1978 /* Assemble the matrices so that the correct flags are set */ 1979 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1980 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1981 1982 /* Free work space */ 1983 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1984 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1985 *B = newmat; 1986 PetscFunctionReturn(0); 1987 } 1988 1989 static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1990 { 1991 PetscErrorCode ierr; 1992 PetscInt i; 1993 1994 PetscFunctionBegin; 1995 if (scall == MAT_INITIAL_MATRIX) { 1996 ierr = PetscCalloc1(n+1,B);CHKERRQ(ierr); 1997 } 1998 1999 for (i=0; i<n; i++) { 2000 ierr = MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 2001 } 2002 PetscFunctionReturn(0); 2003 } 2004 2005 static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 2006 { 2007 PetscFunctionBegin; 2008 PetscFunctionReturn(0); 2009 } 2010 2011 static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 2012 { 2013 PetscFunctionBegin; 2014 PetscFunctionReturn(0); 2015 } 2016 2017 static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 2018 { 2019 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data; 2020 PetscErrorCode ierr; 2021 const PetscScalar *va; 2022 PetscScalar *vb; 2023 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 2024 2025 PetscFunctionBegin; 2026 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 2027 if (A->ops->copy != B->ops->copy) { 2028 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2029 PetscFunctionReturn(0); 2030 } 2031 if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 2032 ierr = MatDenseGetArrayRead(A,&va);CHKERRQ(ierr); 2033 ierr = MatDenseGetArray(B,&vb);CHKERRQ(ierr); 2034 if (lda1>m || lda2>m) { 2035 for (j=0; j<n; j++) { 2036 ierr = PetscArraycpy(vb+j*lda2,va+j*lda1,m);CHKERRQ(ierr); 2037 } 2038 } else { 2039 ierr = PetscArraycpy(vb,va,A->rmap->n*A->cmap->n);CHKERRQ(ierr); 2040 } 2041 ierr = MatDenseRestoreArray(B,&vb);CHKERRQ(ierr); 2042 ierr = MatDenseRestoreArrayRead(A,&va);CHKERRQ(ierr); 2043 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2044 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2045 PetscFunctionReturn(0); 2046 } 2047 2048 static PetscErrorCode MatSetUp_SeqDense(Mat A) 2049 { 2050 PetscErrorCode ierr; 2051 2052 PetscFunctionBegin; 2053 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 2054 PetscFunctionReturn(0); 2055 } 2056 2057 static PetscErrorCode MatConjugate_SeqDense(Mat A) 2058 { 2059 PetscInt i,nz = A->rmap->n*A->cmap->n; 2060 PetscScalar *aa; 2061 PetscErrorCode ierr; 2062 2063 PetscFunctionBegin; 2064 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr); 2065 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 2066 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr); 2067 PetscFunctionReturn(0); 2068 } 2069 2070 static PetscErrorCode MatRealPart_SeqDense(Mat A) 2071 { 2072 PetscInt i,nz = A->rmap->n*A->cmap->n; 2073 PetscScalar *aa; 2074 PetscErrorCode ierr; 2075 2076 PetscFunctionBegin; 2077 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr); 2078 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 2079 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr); 2080 PetscFunctionReturn(0); 2081 } 2082 2083 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 2084 { 2085 PetscInt i,nz = A->rmap->n*A->cmap->n; 2086 PetscScalar *aa; 2087 PetscErrorCode ierr; 2088 2089 PetscFunctionBegin; 2090 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr); 2091 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 2092 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr); 2093 PetscFunctionReturn(0); 2094 } 2095 2096 /* ----------------------------------------------------------------*/ 2097 PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 2098 { 2099 PetscErrorCode ierr; 2100 2101 PetscFunctionBegin; 2102 if (scall == MAT_INITIAL_MATRIX) { 2103 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2104 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 2105 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2106 } 2107 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 2108 if ((*C)->ops->matmultnumeric) { 2109 ierr = (*(*C)->ops->matmultnumeric)(A,B,*C);CHKERRQ(ierr); 2110 } else { 2111 ierr = MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 2112 } 2113 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 2114 PetscFunctionReturn(0); 2115 } 2116 2117 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 2118 { 2119 PetscErrorCode ierr; 2120 PetscInt m=A->rmap->n,n=B->cmap->n; 2121 Mat Cmat; 2122 PetscBool flg; 2123 2124 PetscFunctionBegin; 2125 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 2126 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 2127 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 2128 ierr = MatSetType(Cmat,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 2129 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 2130 *C = Cmat; 2131 PetscFunctionReturn(0); 2132 } 2133 2134 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2135 { 2136 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2137 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2138 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2139 PetscBLASInt m,n,k; 2140 const PetscScalar *av,*bv; 2141 PetscScalar *cv; 2142 PetscScalar _DOne=1.0,_DZero=0.0; 2143 PetscErrorCode ierr; 2144 PetscBool flg; 2145 2146 PetscFunctionBegin; 2147 /* Handle case where where user provided the final C matrix rather than calling MatMatMult() with MAT_INITIAL_MATRIX*/ 2148 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);CHKERRQ(ierr); 2149 if (flg) { 2150 C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqDense; 2151 ierr = (*C->ops->matmultnumeric)(A,B,C);CHKERRQ(ierr); 2152 PetscFunctionReturn(0); 2153 } 2154 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&flg);CHKERRQ(ierr); 2155 if (flg) { 2156 C->ops->matmultnumeric = MatMatMultNumeric_SeqBAIJ_SeqDense; 2157 ierr = (*C->ops->matmultnumeric)(A,B,C);CHKERRQ(ierr); 2158 PetscFunctionReturn(0); 2159 } 2160 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2161 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2162 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 2163 if (!m || !n || !k) PetscFunctionReturn(0); 2164 ierr = MatDenseGetArrayRead(A,&av);CHKERRQ(ierr); 2165 ierr = MatDenseGetArrayRead(B,&bv);CHKERRQ(ierr); 2166 ierr = MatDenseGetArray(C,&cv);CHKERRQ(ierr); 2167 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda)); 2168 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr); 2169 ierr = MatDenseRestoreArrayRead(A,&av);CHKERRQ(ierr); 2170 ierr = MatDenseRestoreArrayRead(B,&bv);CHKERRQ(ierr); 2171 ierr = MatDenseRestoreArray(C,&cv);CHKERRQ(ierr); 2172 PetscFunctionReturn(0); 2173 } 2174 2175 PetscErrorCode MatMatTransposeMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 2176 { 2177 PetscErrorCode ierr; 2178 2179 PetscFunctionBegin; 2180 if (scall == MAT_INITIAL_MATRIX) { 2181 ierr = PetscLogEventBegin(MAT_MatTransposeMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2182 ierr = MatMatTransposeMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 2183 ierr = PetscLogEventEnd(MAT_MatTransposeMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2184 } 2185 ierr = PetscLogEventBegin(MAT_MatTransposeMultNumeric,A,B,0,0);CHKERRQ(ierr); 2186 ierr = MatMatTransposeMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 2187 ierr = PetscLogEventEnd(MAT_MatTransposeMultNumeric,A,B,0,0);CHKERRQ(ierr); 2188 PetscFunctionReturn(0); 2189 } 2190 2191 PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 2192 { 2193 PetscErrorCode ierr; 2194 PetscInt m=A->rmap->n,n=B->rmap->n; 2195 Mat Cmat; 2196 PetscBool flg; 2197 2198 PetscFunctionBegin; 2199 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 2200 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 2201 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 2202 ierr = MatSetType(Cmat,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 2203 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 2204 *C = Cmat; 2205 PetscFunctionReturn(0); 2206 } 2207 2208 PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2209 { 2210 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2211 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2212 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2213 PetscBLASInt m,n,k; 2214 PetscScalar _DOne=1.0,_DZero=0.0; 2215 PetscErrorCode ierr; 2216 2217 PetscFunctionBegin; 2218 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2219 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2220 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 2221 if (!m || !n || !k) PetscFunctionReturn(0); 2222 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 2223 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr); 2224 PetscFunctionReturn(0); 2225 } 2226 2227 PetscErrorCode MatTransposeMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 2228 { 2229 PetscErrorCode ierr; 2230 2231 PetscFunctionBegin; 2232 if (scall == MAT_INITIAL_MATRIX) { 2233 ierr = PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2234 ierr = MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 2235 ierr = PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2236 } 2237 ierr = PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 2238 ierr = MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 2239 ierr = PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 2240 PetscFunctionReturn(0); 2241 } 2242 2243 PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 2244 { 2245 PetscErrorCode ierr; 2246 PetscInt m=A->cmap->n,n=B->cmap->n; 2247 Mat Cmat; 2248 PetscBool flg; 2249 2250 PetscFunctionBegin; 2251 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 2252 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 2253 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 2254 ierr = MatSetType(Cmat,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 2255 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 2256 *C = Cmat; 2257 PetscFunctionReturn(0); 2258 } 2259 2260 PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2261 { 2262 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2263 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2264 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2265 PetscBLASInt m,n,k; 2266 PetscScalar _DOne=1.0,_DZero=0.0; 2267 PetscErrorCode ierr; 2268 2269 PetscFunctionBegin; 2270 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2271 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2272 ierr = PetscBLASIntCast(A->rmap->n,&k);CHKERRQ(ierr); 2273 if (!m || !n || !k) PetscFunctionReturn(0); 2274 PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 2275 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr); 2276 PetscFunctionReturn(0); 2277 } 2278 2279 static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 2280 { 2281 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2282 PetscErrorCode ierr; 2283 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2284 PetscScalar *x; 2285 const PetscScalar *aa; 2286 2287 PetscFunctionBegin; 2288 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2289 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2290 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2291 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2292 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2293 for (i=0; i<m; i++) { 2294 x[i] = aa[i]; if (idx) idx[i] = 0; 2295 for (j=1; j<n; j++) { 2296 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;} 2297 } 2298 } 2299 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2300 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2301 PetscFunctionReturn(0); 2302 } 2303 2304 static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 2305 { 2306 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2307 PetscErrorCode ierr; 2308 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2309 PetscScalar *x; 2310 PetscReal atmp; 2311 const PetscScalar *aa; 2312 2313 PetscFunctionBegin; 2314 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2315 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2316 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2317 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2318 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2319 for (i=0; i<m; i++) { 2320 x[i] = PetscAbsScalar(aa[i]); 2321 for (j=1; j<n; j++) { 2322 atmp = PetscAbsScalar(aa[i+a->lda*j]); 2323 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 2324 } 2325 } 2326 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2327 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2328 PetscFunctionReturn(0); 2329 } 2330 2331 static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 2332 { 2333 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2334 PetscErrorCode ierr; 2335 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2336 PetscScalar *x; 2337 const PetscScalar *aa; 2338 2339 PetscFunctionBegin; 2340 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2341 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2342 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2343 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2344 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2345 for (i=0; i<m; i++) { 2346 x[i] = aa[i]; if (idx) idx[i] = 0; 2347 for (j=1; j<n; j++) { 2348 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;} 2349 } 2350 } 2351 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2352 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2353 PetscFunctionReturn(0); 2354 } 2355 2356 static PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 2357 { 2358 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2359 PetscErrorCode ierr; 2360 PetscScalar *x; 2361 const PetscScalar *aa; 2362 2363 PetscFunctionBegin; 2364 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2365 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2366 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2367 ierr = PetscArraycpy(x,aa+col*a->lda,A->rmap->n);CHKERRQ(ierr); 2368 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2369 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2370 PetscFunctionReturn(0); 2371 } 2372 2373 PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms) 2374 { 2375 PetscErrorCode ierr; 2376 PetscInt i,j,m,n; 2377 const PetscScalar *a; 2378 2379 PetscFunctionBegin; 2380 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 2381 ierr = PetscArrayzero(norms,n);CHKERRQ(ierr); 2382 ierr = MatDenseGetArrayRead(A,&a);CHKERRQ(ierr); 2383 if (type == NORM_2) { 2384 for (i=0; i<n; i++) { 2385 for (j=0; j<m; j++) { 2386 norms[i] += PetscAbsScalar(a[j]*a[j]); 2387 } 2388 a += m; 2389 } 2390 } else if (type == NORM_1) { 2391 for (i=0; i<n; i++) { 2392 for (j=0; j<m; j++) { 2393 norms[i] += PetscAbsScalar(a[j]); 2394 } 2395 a += m; 2396 } 2397 } else if (type == NORM_INFINITY) { 2398 for (i=0; i<n; i++) { 2399 for (j=0; j<m; j++) { 2400 norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]); 2401 } 2402 a += m; 2403 } 2404 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType"); 2405 ierr = MatDenseRestoreArrayRead(A,&a);CHKERRQ(ierr); 2406 if (type == NORM_2) { 2407 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 2408 } 2409 PetscFunctionReturn(0); 2410 } 2411 2412 static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx) 2413 { 2414 PetscErrorCode ierr; 2415 PetscScalar *a; 2416 PetscInt m,n,i; 2417 2418 PetscFunctionBegin; 2419 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 2420 ierr = MatDenseGetArray(x,&a);CHKERRQ(ierr); 2421 for (i=0; i<m*n; i++) { 2422 ierr = PetscRandomGetValue(rctx,a+i);CHKERRQ(ierr); 2423 } 2424 ierr = MatDenseRestoreArray(x,&a);CHKERRQ(ierr); 2425 PetscFunctionReturn(0); 2426 } 2427 2428 static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool *missing,PetscInt *d) 2429 { 2430 PetscFunctionBegin; 2431 *missing = PETSC_FALSE; 2432 PetscFunctionReturn(0); 2433 } 2434 2435 /* vals is not const */ 2436 static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals) 2437 { 2438 PetscErrorCode ierr; 2439 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2440 PetscScalar *v; 2441 2442 PetscFunctionBegin; 2443 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2444 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 2445 *vals = v+col*a->lda; 2446 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 2447 PetscFunctionReturn(0); 2448 } 2449 2450 static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals) 2451 { 2452 PetscFunctionBegin; 2453 *vals = 0; /* user cannot accidently use the array later */ 2454 PetscFunctionReturn(0); 2455 } 2456 2457 /* -------------------------------------------------------------------*/ 2458 static struct _MatOps MatOps_Values = { MatSetValues_SeqDense, 2459 MatGetRow_SeqDense, 2460 MatRestoreRow_SeqDense, 2461 MatMult_SeqDense, 2462 /* 4*/ MatMultAdd_SeqDense, 2463 MatMultTranspose_SeqDense, 2464 MatMultTransposeAdd_SeqDense, 2465 0, 2466 0, 2467 0, 2468 /* 10*/ 0, 2469 MatLUFactor_SeqDense, 2470 MatCholeskyFactor_SeqDense, 2471 MatSOR_SeqDense, 2472 MatTranspose_SeqDense, 2473 /* 15*/ MatGetInfo_SeqDense, 2474 MatEqual_SeqDense, 2475 MatGetDiagonal_SeqDense, 2476 MatDiagonalScale_SeqDense, 2477 MatNorm_SeqDense, 2478 /* 20*/ MatAssemblyBegin_SeqDense, 2479 MatAssemblyEnd_SeqDense, 2480 MatSetOption_SeqDense, 2481 MatZeroEntries_SeqDense, 2482 /* 24*/ MatZeroRows_SeqDense, 2483 0, 2484 0, 2485 0, 2486 0, 2487 /* 29*/ MatSetUp_SeqDense, 2488 0, 2489 0, 2490 0, 2491 0, 2492 /* 34*/ MatDuplicate_SeqDense, 2493 0, 2494 0, 2495 0, 2496 0, 2497 /* 39*/ MatAXPY_SeqDense, 2498 MatCreateSubMatrices_SeqDense, 2499 0, 2500 MatGetValues_SeqDense, 2501 MatCopy_SeqDense, 2502 /* 44*/ MatGetRowMax_SeqDense, 2503 MatScale_SeqDense, 2504 MatShift_Basic, 2505 0, 2506 MatZeroRowsColumns_SeqDense, 2507 /* 49*/ MatSetRandom_SeqDense, 2508 0, 2509 0, 2510 0, 2511 0, 2512 /* 54*/ 0, 2513 0, 2514 0, 2515 0, 2516 0, 2517 /* 59*/ 0, 2518 MatDestroy_SeqDense, 2519 MatView_SeqDense, 2520 0, 2521 0, 2522 /* 64*/ 0, 2523 0, 2524 0, 2525 0, 2526 0, 2527 /* 69*/ MatGetRowMaxAbs_SeqDense, 2528 0, 2529 0, 2530 0, 2531 0, 2532 /* 74*/ 0, 2533 0, 2534 0, 2535 0, 2536 0, 2537 /* 79*/ 0, 2538 0, 2539 0, 2540 0, 2541 /* 83*/ MatLoad_SeqDense, 2542 0, 2543 MatIsHermitian_SeqDense, 2544 0, 2545 0, 2546 0, 2547 /* 89*/ MatMatMult_SeqDense_SeqDense, 2548 MatMatMultSymbolic_SeqDense_SeqDense, 2549 MatMatMultNumeric_SeqDense_SeqDense, 2550 MatPtAP_SeqDense_SeqDense, 2551 MatPtAPSymbolic_SeqDense_SeqDense, 2552 /* 94*/ MatPtAPNumeric_SeqDense_SeqDense, 2553 MatMatTransposeMult_SeqDense_SeqDense, 2554 MatMatTransposeMultSymbolic_SeqDense_SeqDense, 2555 MatMatTransposeMultNumeric_SeqDense_SeqDense, 2556 0, 2557 /* 99*/ 0, 2558 0, 2559 0, 2560 MatConjugate_SeqDense, 2561 0, 2562 /*104*/ 0, 2563 MatRealPart_SeqDense, 2564 MatImaginaryPart_SeqDense, 2565 0, 2566 0, 2567 /*109*/ 0, 2568 0, 2569 MatGetRowMin_SeqDense, 2570 MatGetColumnVector_SeqDense, 2571 MatMissingDiagonal_SeqDense, 2572 /*114*/ 0, 2573 0, 2574 0, 2575 0, 2576 0, 2577 /*119*/ 0, 2578 0, 2579 0, 2580 0, 2581 0, 2582 /*124*/ 0, 2583 MatGetColumnNorms_SeqDense, 2584 0, 2585 0, 2586 0, 2587 /*129*/ 0, 2588 MatTransposeMatMult_SeqDense_SeqDense, 2589 MatTransposeMatMultSymbolic_SeqDense_SeqDense, 2590 MatTransposeMatMultNumeric_SeqDense_SeqDense, 2591 0, 2592 /*134*/ 0, 2593 0, 2594 0, 2595 0, 2596 0, 2597 /*139*/ 0, 2598 0, 2599 0, 2600 0, 2601 0, 2602 /*144*/ MatCreateMPIMatConcatenateSeqMat_SeqDense 2603 }; 2604 2605 /*@C 2606 MatCreateSeqDense - Creates a sequential dense matrix that 2607 is stored in column major order (the usual Fortran 77 manner). Many 2608 of the matrix operations use the BLAS and LAPACK routines. 2609 2610 Collective 2611 2612 Input Parameters: 2613 + comm - MPI communicator, set to PETSC_COMM_SELF 2614 . m - number of rows 2615 . n - number of columns 2616 - data - optional location of matrix data in column major order. Set data=NULL for PETSc 2617 to control all matrix memory allocation. 2618 2619 Output Parameter: 2620 . A - the matrix 2621 2622 Notes: 2623 The data input variable is intended primarily for Fortran programmers 2624 who wish to allocate their own matrix memory space. Most users should 2625 set data=NULL. 2626 2627 Level: intermediate 2628 2629 .seealso: MatCreate(), MatCreateDense(), MatSetValues() 2630 @*/ 2631 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 2632 { 2633 PetscErrorCode ierr; 2634 2635 PetscFunctionBegin; 2636 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2637 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2638 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 2639 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 2640 PetscFunctionReturn(0); 2641 } 2642 2643 /*@C 2644 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 2645 2646 Collective 2647 2648 Input Parameters: 2649 + B - the matrix 2650 - data - the array (or NULL) 2651 2652 Notes: 2653 The data input variable is intended primarily for Fortran programmers 2654 who wish to allocate their own matrix memory space. Most users should 2655 need not call this routine. 2656 2657 Level: intermediate 2658 2659 .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA() 2660 2661 @*/ 2662 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 2663 { 2664 PetscErrorCode ierr; 2665 2666 PetscFunctionBegin; 2667 ierr = PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));CHKERRQ(ierr); 2668 PetscFunctionReturn(0); 2669 } 2670 2671 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 2672 { 2673 Mat_SeqDense *b; 2674 PetscErrorCode ierr; 2675 2676 PetscFunctionBegin; 2677 B->preallocated = PETSC_TRUE; 2678 2679 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2680 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2681 2682 b = (Mat_SeqDense*)B->data; 2683 b->Mmax = B->rmap->n; 2684 b->Nmax = B->cmap->n; 2685 if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n; 2686 2687 ierr = PetscIntMultError(b->lda,b->Nmax,NULL);CHKERRQ(ierr); 2688 if (!data) { /* petsc-allocated storage */ 2689 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2690 ierr = PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);CHKERRQ(ierr); 2691 ierr = PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 2692 2693 b->user_alloc = PETSC_FALSE; 2694 } else { /* user-allocated storage */ 2695 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2696 b->v = data; 2697 b->user_alloc = PETSC_TRUE; 2698 } 2699 B->assembled = PETSC_TRUE; 2700 PetscFunctionReturn(0); 2701 } 2702 2703 #if defined(PETSC_HAVE_ELEMENTAL) 2704 PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 2705 { 2706 Mat mat_elemental; 2707 PetscErrorCode ierr; 2708 const PetscScalar *array; 2709 PetscScalar *v_colwise; 2710 PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols; 2711 2712 PetscFunctionBegin; 2713 ierr = PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);CHKERRQ(ierr); 2714 ierr = MatDenseGetArrayRead(A,&array);CHKERRQ(ierr); 2715 /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */ 2716 k = 0; 2717 for (j=0; j<N; j++) { 2718 cols[j] = j; 2719 for (i=0; i<M; i++) { 2720 v_colwise[j*M+i] = array[k++]; 2721 } 2722 } 2723 for (i=0; i<M; i++) { 2724 rows[i] = i; 2725 } 2726 ierr = MatDenseRestoreArrayRead(A,&array);CHKERRQ(ierr); 2727 2728 ierr = MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);CHKERRQ(ierr); 2729 ierr = MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 2730 ierr = MatSetType(mat_elemental,MATELEMENTAL);CHKERRQ(ierr); 2731 ierr = MatSetUp(mat_elemental);CHKERRQ(ierr); 2732 2733 /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */ 2734 ierr = MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);CHKERRQ(ierr); 2735 ierr = MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2736 ierr = MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2737 ierr = PetscFree3(v_colwise,rows,cols);CHKERRQ(ierr); 2738 2739 if (reuse == MAT_INPLACE_MATRIX) { 2740 ierr = MatHeaderReplace(A,&mat_elemental);CHKERRQ(ierr); 2741 } else { 2742 *newmat = mat_elemental; 2743 } 2744 PetscFunctionReturn(0); 2745 } 2746 #endif 2747 2748 /*@C 2749 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 2750 2751 Input parameter: 2752 + A - the matrix 2753 - lda - the leading dimension 2754 2755 Notes: 2756 This routine is to be used in conjunction with MatSeqDenseSetPreallocation(); 2757 it asserts that the preallocation has a leading dimension (the LDA parameter 2758 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 2759 2760 Level: intermediate 2761 2762 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 2763 2764 @*/ 2765 PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda) 2766 { 2767 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2768 2769 PetscFunctionBegin; 2770 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); 2771 b->lda = lda; 2772 b->changelda = PETSC_FALSE; 2773 b->Mmax = PetscMax(b->Mmax,lda); 2774 PetscFunctionReturn(0); 2775 } 2776 2777 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 2778 { 2779 PetscErrorCode ierr; 2780 PetscMPIInt size; 2781 2782 PetscFunctionBegin; 2783 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2784 if (size == 1) { 2785 if (scall == MAT_INITIAL_MATRIX) { 2786 ierr = MatDuplicate(inmat,MAT_COPY_VALUES,outmat);CHKERRQ(ierr); 2787 } else { 2788 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2789 } 2790 } else { 2791 ierr = MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 2792 } 2793 PetscFunctionReturn(0); 2794 } 2795 2796 /*MC 2797 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 2798 2799 Options Database Keys: 2800 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 2801 2802 Level: beginner 2803 2804 .seealso: MatCreateSeqDense() 2805 2806 M*/ 2807 PetscErrorCode MatCreate_SeqDense(Mat B) 2808 { 2809 Mat_SeqDense *b; 2810 PetscErrorCode ierr; 2811 PetscMPIInt size; 2812 2813 PetscFunctionBegin; 2814 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 2815 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 2816 2817 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 2818 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2819 B->data = (void*)b; 2820 2821 b->roworiented = PETSC_TRUE; 2822 2823 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetLDA_C",MatDenseGetLDA_SeqDense);CHKERRQ(ierr); 2824 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2825 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2826 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);CHKERRQ(ierr); 2827 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);CHKERRQ(ierr); 2828 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2829 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2830 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);CHKERRQ(ierr); 2831 #if defined(PETSC_HAVE_ELEMENTAL) 2832 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);CHKERRQ(ierr); 2833 #endif 2834 #if defined(PETSC_HAVE_CUDA) 2835 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqdensecuda_C",MatConvert_SeqDense_SeqDenseCUDA);CHKERRQ(ierr); 2836 #endif 2837 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 2838 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2839 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2840 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2841 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqbaij_seqdense_C",MatMatMult_SeqBAIJ_SeqDense);CHKERRQ(ierr); 2842 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqbaij_seqdense_C",MatMatMultSymbolic_SeqBAIJ_SeqDense);CHKERRQ(ierr); 2843 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqbaij_seqdense_C",MatMatMultNumeric_SeqBAIJ_SeqDense);CHKERRQ(ierr); 2844 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaij_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2845 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijperm_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2846 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijperm_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2847 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijperm_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2848 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijperm_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2849 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijsell_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2850 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijsell_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2851 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijsell_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2852 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijsell_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2853 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijmkl_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2854 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijmkl_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2855 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijmkl_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2856 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijmkl_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2857 2858 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2859 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2860 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2861 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijperm_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2862 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijperm_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2863 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijperm_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2864 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijsell_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2865 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijsell_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2866 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijsell_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2867 2868 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijmkl_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2869 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijmkl_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2870 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijmkl_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2871 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);CHKERRQ(ierr); 2872 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);CHKERRQ(ierr); 2873 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 2874 PetscFunctionReturn(0); 2875 } 2876 2877 /*@C 2878 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. 2879 2880 Not Collective 2881 2882 Input Parameter: 2883 + mat - a MATSEQDENSE or MATMPIDENSE matrix 2884 - col - column index 2885 2886 Output Parameter: 2887 . vals - pointer to the data 2888 2889 Level: intermediate 2890 2891 .seealso: MatDenseRestoreColumn() 2892 @*/ 2893 PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals) 2894 { 2895 PetscErrorCode ierr; 2896 2897 PetscFunctionBegin; 2898 ierr = PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));CHKERRQ(ierr); 2899 PetscFunctionReturn(0); 2900 } 2901 2902 /*@C 2903 MatDenseRestoreColumn - returns access to a column of a dense matrix which is returned by MatDenseGetColumn(). 2904 2905 Not Collective 2906 2907 Input Parameter: 2908 . mat - a MATSEQDENSE or MATMPIDENSE matrix 2909 2910 Output Parameter: 2911 . vals - pointer to the data 2912 2913 Level: intermediate 2914 2915 .seealso: MatDenseGetColumn() 2916 @*/ 2917 PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals) 2918 { 2919 PetscErrorCode ierr; 2920 2921 PetscFunctionBegin; 2922 ierr = PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));CHKERRQ(ierr); 2923 PetscFunctionReturn(0); 2924 } 2925