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 17 PetscFunctionBegin; 18 PetscCheck(A->rmap->n == A->cmap->n,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot symmetrize a rectangular matrix"); 19 PetscCall(MatDenseGetArray(A,&v)); 20 if (!hermitian) { 21 for (k=0;k<n;k++) { 22 for (j=k;j<n;j++) { 23 v[j*mat->lda + k] = v[k*mat->lda + j]; 24 } 25 } 26 } else { 27 for (k=0;k<n;k++) { 28 for (j=k;j<n;j++) { 29 v[j*mat->lda + k] = PetscConj(v[k*mat->lda + j]); 30 } 31 } 32 } 33 PetscCall(MatDenseRestoreArray(A,&v)); 34 PetscFunctionReturn(0); 35 } 36 37 PETSC_EXTERN PetscErrorCode MatSeqDenseInvertFactors_Private(Mat A) 38 { 39 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 40 PetscBLASInt info,n; 41 42 PetscFunctionBegin; 43 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 44 PetscCall(PetscBLASIntCast(A->cmap->n,&n)); 45 if (A->factortype == MAT_FACTOR_LU) { 46 PetscCheck(mat->pivots,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present"); 47 if (!mat->fwork) { 48 mat->lfwork = n; 49 PetscCall(PetscMalloc1(mat->lfwork,&mat->fwork)); 50 PetscCall(PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt))); 51 } 52 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 53 PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info)); 54 PetscCall(PetscFPTrapPop()); 55 PetscCall(PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0)); 56 } else if (A->factortype == MAT_FACTOR_CHOLESKY) { 57 if (A->spd) { 58 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 59 PetscStackCallBLAS("LAPACKpotri",LAPACKpotri_("L",&n,mat->v,&mat->lda,&info)); 60 PetscCall(PetscFPTrapPop()); 61 PetscCall(MatSeqDenseSymmetrize_Private(A,PETSC_TRUE)); 62 #if defined(PETSC_USE_COMPLEX) 63 } else if (A->hermitian) { 64 PetscCheck(mat->pivots,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present"); 65 PetscCheck(mat->fwork,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Fwork not present"); 66 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 67 PetscStackCallBLAS("LAPACKhetri",LAPACKhetri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info)); 68 PetscCall(PetscFPTrapPop()); 69 PetscCall(MatSeqDenseSymmetrize_Private(A,PETSC_TRUE)); 70 #endif 71 } else { /* symmetric case */ 72 PetscCheck(mat->pivots,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present"); 73 PetscCheck(mat->fwork,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Fwork not present"); 74 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 75 PetscStackCallBLAS("LAPACKsytri",LAPACKsytri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info)); 76 PetscCall(PetscFPTrapPop()); 77 PetscCall(MatSeqDenseSymmetrize_Private(A,PETSC_FALSE)); 78 } 79 PetscCheck(!info,PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad Inversion: zero pivot in row %" PetscInt_FMT,(PetscInt)info-1); 80 PetscCall(PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0)); 81 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 82 83 A->ops->solve = NULL; 84 A->ops->matsolve = NULL; 85 A->ops->solvetranspose = NULL; 86 A->ops->matsolvetranspose = NULL; 87 A->ops->solveadd = NULL; 88 A->ops->solvetransposeadd = NULL; 89 A->factortype = MAT_FACTOR_NONE; 90 PetscCall(PetscFree(A->solvertype)); 91 PetscFunctionReturn(0); 92 } 93 94 PetscErrorCode MatZeroRowsColumns_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 95 { 96 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 97 PetscInt m = l->lda, n = A->cmap->n,r = A->rmap->n, i,j; 98 PetscScalar *slot,*bb,*v; 99 const PetscScalar *xx; 100 101 PetscFunctionBegin; 102 if (PetscDefined(USE_DEBUG)) { 103 for (i=0; i<N; i++) { 104 PetscCheck(rows[i] >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 105 PetscCheck(rows[i] < A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %" PetscInt_FMT " requested to be zeroed greater than or equal number of rows %" PetscInt_FMT,rows[i],A->rmap->n); 106 PetscCheck(rows[i] < A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Col %" PetscInt_FMT " requested to be zeroed greater than or equal number of cols %" PetscInt_FMT,rows[i],A->cmap->n); 107 } 108 } 109 if (!N) PetscFunctionReturn(0); 110 111 /* fix right hand side if needed */ 112 if (x && b) { 113 Vec xt; 114 115 PetscCheck(A->rmap->n == A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 116 PetscCall(VecDuplicate(x,&xt)); 117 PetscCall(VecCopy(x,xt)); 118 PetscCall(VecScale(xt,-1.0)); 119 PetscCall(MatMultAdd(A,xt,b,b)); 120 PetscCall(VecDestroy(&xt)); 121 PetscCall(VecGetArrayRead(x,&xx)); 122 PetscCall(VecGetArray(b,&bb)); 123 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]]; 124 PetscCall(VecRestoreArrayRead(x,&xx)); 125 PetscCall(VecRestoreArray(b,&bb)); 126 } 127 128 PetscCall(MatDenseGetArray(A,&v)); 129 for (i=0; i<N; i++) { 130 slot = v + rows[i]*m; 131 PetscCall(PetscArrayzero(slot,r)); 132 } 133 for (i=0; i<N; i++) { 134 slot = v + rows[i]; 135 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 136 } 137 if (diag != 0.0) { 138 PetscCheck(A->rmap->n == A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 139 for (i=0; i<N; i++) { 140 slot = v + (m+1)*rows[i]; 141 *slot = diag; 142 } 143 } 144 PetscCall(MatDenseRestoreArray(A,&v)); 145 PetscFunctionReturn(0); 146 } 147 148 PetscErrorCode MatPtAPNumeric_SeqDense_SeqDense(Mat A,Mat P,Mat C) 149 { 150 Mat_SeqDense *c = (Mat_SeqDense*)(C->data); 151 152 PetscFunctionBegin; 153 if (c->ptapwork) { 154 PetscCall((*C->ops->matmultnumeric)(A,P,c->ptapwork)); 155 PetscCall((*C->ops->transposematmultnumeric)(P,c->ptapwork,C)); 156 } else SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"Must call MatPtAPSymbolic_SeqDense_SeqDense() first"); 157 PetscFunctionReturn(0); 158 } 159 160 PetscErrorCode MatPtAPSymbolic_SeqDense_SeqDense(Mat A,Mat P,PetscReal fill,Mat C) 161 { 162 Mat_SeqDense *c; 163 PetscBool cisdense; 164 165 PetscFunctionBegin; 166 PetscCall(MatSetSizes(C,P->cmap->n,P->cmap->n,P->cmap->N,P->cmap->N)); 167 PetscCall(PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"")); 168 if (!cisdense) { 169 PetscBool flg; 170 171 PetscCall(PetscObjectTypeCompare((PetscObject)P,((PetscObject)A)->type_name,&flg)); 172 PetscCall(MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE)); 173 } 174 PetscCall(MatSetUp(C)); 175 c = (Mat_SeqDense*)C->data; 176 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&c->ptapwork)); 177 PetscCall(MatSetSizes(c->ptapwork,A->rmap->n,P->cmap->n,A->rmap->N,P->cmap->N)); 178 PetscCall(MatSetType(c->ptapwork,((PetscObject)C)->type_name)); 179 PetscCall(MatSetUp(c->ptapwork)); 180 PetscFunctionReturn(0); 181 } 182 183 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 184 { 185 Mat B = NULL; 186 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 187 Mat_SeqDense *b; 188 PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->N,i; 189 const MatScalar *av; 190 PetscBool isseqdense; 191 192 PetscFunctionBegin; 193 if (reuse == MAT_REUSE_MATRIX) { 194 PetscCall(PetscObjectTypeCompare((PetscObject)*newmat,MATSEQDENSE,&isseqdense)); 195 PetscCheck(isseqdense,PetscObjectComm((PetscObject)*newmat),PETSC_ERR_USER,"Cannot reuse matrix of type %s",((PetscObject)(*newmat))->type_name); 196 } 197 if (reuse != MAT_REUSE_MATRIX) { 198 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&B)); 199 PetscCall(MatSetSizes(B,m,n,m,n)); 200 PetscCall(MatSetType(B,MATSEQDENSE)); 201 PetscCall(MatSeqDenseSetPreallocation(B,NULL)); 202 b = (Mat_SeqDense*)(B->data); 203 } else { 204 b = (Mat_SeqDense*)((*newmat)->data); 205 PetscCall(PetscArrayzero(b->v,m*n)); 206 } 207 PetscCall(MatSeqAIJGetArrayRead(A,&av)); 208 for (i=0; i<m; i++) { 209 PetscInt j; 210 for (j=0;j<ai[1]-ai[0];j++) { 211 b->v[*aj*m+i] = *av; 212 aj++; 213 av++; 214 } 215 ai++; 216 } 217 PetscCall(MatSeqAIJRestoreArrayRead(A,&av)); 218 219 if (reuse == MAT_INPLACE_MATRIX) { 220 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 221 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 222 PetscCall(MatHeaderReplace(A,&B)); 223 } else { 224 if (B) *newmat = B; 225 PetscCall(MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY)); 226 PetscCall(MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY)); 227 } 228 PetscFunctionReturn(0); 229 } 230 231 PETSC_INTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 232 { 233 Mat B = NULL; 234 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 235 PetscInt i, j; 236 PetscInt *rows, *nnz; 237 MatScalar *aa = a->v, *vals; 238 239 PetscFunctionBegin; 240 PetscCall(PetscCalloc3(A->rmap->n,&rows,A->rmap->n,&nnz,A->rmap->n,&vals)); 241 if (reuse != MAT_REUSE_MATRIX) { 242 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&B)); 243 PetscCall(MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N)); 244 PetscCall(MatSetType(B,MATSEQAIJ)); 245 for (j=0; j<A->cmap->n; j++) { 246 for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || (i == j && A->cmap->n == A->rmap->n)) ++nnz[i]; 247 aa += a->lda; 248 } 249 PetscCall(MatSeqAIJSetPreallocation(B,PETSC_DETERMINE,nnz)); 250 } else B = *newmat; 251 aa = a->v; 252 for (j=0; j<A->cmap->n; j++) { 253 PetscInt numRows = 0; 254 for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || (i == j && A->cmap->n == A->rmap->n)) {rows[numRows] = i; vals[numRows++] = aa[i];} 255 PetscCall(MatSetValues(B,numRows,rows,1,&j,vals,INSERT_VALUES)); 256 aa += a->lda; 257 } 258 PetscCall(PetscFree3(rows,nnz,vals)); 259 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 260 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 261 262 if (reuse == MAT_INPLACE_MATRIX) { 263 PetscCall(MatHeaderReplace(A,&B)); 264 } else if (reuse != MAT_REUSE_MATRIX) *newmat = B; 265 PetscFunctionReturn(0); 266 } 267 268 PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str) 269 { 270 Mat_SeqDense *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data; 271 const PetscScalar *xv; 272 PetscScalar *yv; 273 PetscBLASInt N,m,ldax = 0,lday = 0,one = 1; 274 275 PetscFunctionBegin; 276 PetscCall(MatDenseGetArrayRead(X,&xv)); 277 PetscCall(MatDenseGetArray(Y,&yv)); 278 PetscCall(PetscBLASIntCast(X->rmap->n*X->cmap->n,&N)); 279 PetscCall(PetscBLASIntCast(X->rmap->n,&m)); 280 PetscCall(PetscBLASIntCast(x->lda,&ldax)); 281 PetscCall(PetscBLASIntCast(y->lda,&lday)); 282 if (ldax>m || lday>m) { 283 PetscInt j; 284 285 for (j=0; j<X->cmap->n; j++) { 286 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&alpha,xv+j*ldax,&one,yv+j*lday,&one)); 287 } 288 } else { 289 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&alpha,xv,&one,yv,&one)); 290 } 291 PetscCall(MatDenseRestoreArrayRead(X,&xv)); 292 PetscCall(MatDenseRestoreArray(Y,&yv)); 293 PetscCall(PetscLogFlops(PetscMax(2.0*N-1,0))); 294 PetscFunctionReturn(0); 295 } 296 297 static PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info) 298 { 299 PetscLogDouble N = A->rmap->n*A->cmap->n; 300 301 PetscFunctionBegin; 302 info->block_size = 1.0; 303 info->nz_allocated = N; 304 info->nz_used = N; 305 info->nz_unneeded = 0; 306 info->assemblies = A->num_ass; 307 info->mallocs = 0; 308 info->memory = ((PetscObject)A)->mem; 309 info->fill_ratio_given = 0; 310 info->fill_ratio_needed = 0; 311 info->factor_mallocs = 0; 312 PetscFunctionReturn(0); 313 } 314 315 PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha) 316 { 317 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 318 PetscScalar *v; 319 PetscBLASInt one = 1,j,nz,lda = 0; 320 321 PetscFunctionBegin; 322 PetscCall(MatDenseGetArray(A,&v)); 323 PetscCall(PetscBLASIntCast(a->lda,&lda)); 324 if (lda>A->rmap->n) { 325 PetscCall(PetscBLASIntCast(A->rmap->n,&nz)); 326 for (j=0; j<A->cmap->n; j++) { 327 PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&alpha,v+j*lda,&one)); 328 } 329 } else { 330 PetscCall(PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz)); 331 PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&alpha,v,&one)); 332 } 333 PetscCall(PetscLogFlops(A->rmap->n*A->cmap->n)); 334 PetscCall(MatDenseRestoreArray(A,&v)); 335 PetscFunctionReturn(0); 336 } 337 338 PetscErrorCode MatShift_SeqDense(Mat A,PetscScalar alpha) 339 { 340 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 341 PetscScalar *v; 342 PetscInt j,k; 343 344 PetscFunctionBegin; 345 PetscCall(MatDenseGetArray(A,&v)); 346 k = PetscMin(A->rmap->n,A->cmap->n); 347 for (j=0; j<k; j++) v[j+j*a->lda] += alpha; 348 PetscCall(PetscLogFlops(k)); 349 PetscCall(MatDenseRestoreArray(A,&v)); 350 PetscFunctionReturn(0); 351 } 352 353 static PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool *fl) 354 { 355 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 356 PetscInt i,j,m = A->rmap->n,N = a->lda; 357 const PetscScalar *v; 358 359 PetscFunctionBegin; 360 *fl = PETSC_FALSE; 361 if (A->rmap->n != A->cmap->n) PetscFunctionReturn(0); 362 PetscCall(MatDenseGetArrayRead(A,&v)); 363 for (i=0; i<m; i++) { 364 for (j=i; j<m; j++) { 365 if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) { 366 goto restore; 367 } 368 } 369 } 370 *fl = PETSC_TRUE; 371 restore: 372 PetscCall(MatDenseRestoreArrayRead(A,&v)); 373 PetscFunctionReturn(0); 374 } 375 376 static PetscErrorCode MatIsSymmetric_SeqDense(Mat A,PetscReal rtol,PetscBool *fl) 377 { 378 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 379 PetscInt i,j,m = A->rmap->n,N = a->lda; 380 const PetscScalar *v; 381 382 PetscFunctionBegin; 383 *fl = PETSC_FALSE; 384 if (A->rmap->n != A->cmap->n) PetscFunctionReturn(0); 385 PetscCall(MatDenseGetArrayRead(A,&v)); 386 for (i=0; i<m; i++) { 387 for (j=i; j<m; j++) { 388 if (PetscAbsScalar(v[i+j*N] - v[j+i*N]) > rtol) { 389 goto restore; 390 } 391 } 392 } 393 *fl = PETSC_TRUE; 394 restore: 395 PetscCall(MatDenseRestoreArrayRead(A,&v)); 396 PetscFunctionReturn(0); 397 } 398 399 PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues) 400 { 401 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 402 PetscInt lda = (PetscInt)mat->lda,j,m,nlda = lda; 403 PetscBool isdensecpu; 404 405 PetscFunctionBegin; 406 PetscCall(PetscLayoutReference(A->rmap,&newi->rmap)); 407 PetscCall(PetscLayoutReference(A->cmap,&newi->cmap)); 408 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { /* propagate LDA */ 409 PetscCall(MatDenseSetLDA(newi,lda)); 410 } 411 PetscCall(PetscObjectTypeCompare((PetscObject)newi,MATSEQDENSE,&isdensecpu)); 412 if (isdensecpu) PetscCall(MatSeqDenseSetPreallocation(newi,NULL)); 413 if (cpvalues == MAT_COPY_VALUES) { 414 const PetscScalar *av; 415 PetscScalar *v; 416 417 PetscCall(MatDenseGetArrayRead(A,&av)); 418 PetscCall(MatDenseGetArrayWrite(newi,&v)); 419 PetscCall(MatDenseGetLDA(newi,&nlda)); 420 m = A->rmap->n; 421 if (lda>m || nlda>m) { 422 for (j=0; j<A->cmap->n; j++) { 423 PetscCall(PetscArraycpy(v+j*nlda,av+j*lda,m)); 424 } 425 } else { 426 PetscCall(PetscArraycpy(v,av,A->rmap->n*A->cmap->n)); 427 } 428 PetscCall(MatDenseRestoreArrayWrite(newi,&v)); 429 PetscCall(MatDenseRestoreArrayRead(A,&av)); 430 } 431 PetscFunctionReturn(0); 432 } 433 434 PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat) 435 { 436 PetscFunctionBegin; 437 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),newmat)); 438 PetscCall(MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n)); 439 PetscCall(MatSetType(*newmat,((PetscObject)A)->type_name)); 440 PetscCall(MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues)); 441 PetscFunctionReturn(0); 442 } 443 444 static PetscErrorCode MatSolve_SeqDense_Internal_LU(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k, PetscBool T) 445 { 446 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 447 PetscBLASInt info; 448 449 PetscFunctionBegin; 450 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 451 PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_(T ? "T" : "N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info)); 452 PetscCall(PetscFPTrapPop()); 453 PetscCheck(!info,PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve %d",(int)info); 454 PetscCall(PetscLogFlops(nrhs*(2.0*m*m - m))); 455 PetscFunctionReturn(0); 456 } 457 458 static PetscErrorCode MatConjugate_SeqDense(Mat); 459 460 static PetscErrorCode MatSolve_SeqDense_Internal_Cholesky(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k, PetscBool T) 461 { 462 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 463 PetscBLASInt info; 464 465 PetscFunctionBegin; 466 if (A->spd) { 467 if (PetscDefined(USE_COMPLEX) && T) PetscCall(MatConjugate_SeqDense(A)); 468 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 469 PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info)); 470 PetscCall(PetscFPTrapPop()); 471 PetscCheck(!info,PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve %d",(int)info); 472 if (PetscDefined(USE_COMPLEX) && T) PetscCall(MatConjugate_SeqDense(A)); 473 #if defined(PETSC_USE_COMPLEX) 474 } else if (A->hermitian) { 475 if (T) PetscCall(MatConjugate_SeqDense(A)); 476 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 477 PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info)); 478 PetscCall(PetscFPTrapPop()); 479 PetscCheck(!info,PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve %d",(int)info); 480 if (T) PetscCall(MatConjugate_SeqDense(A)); 481 #endif 482 } else { /* symmetric case */ 483 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 484 PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info)); 485 PetscCall(PetscFPTrapPop()); 486 PetscCheck(!info,PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve %d",(int)info); 487 } 488 PetscCall(PetscLogFlops(nrhs*(2.0*m*m - m))); 489 PetscFunctionReturn(0); 490 } 491 492 static PetscErrorCode MatSolve_SeqDense_Internal_QR(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k) 493 { 494 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 495 PetscBLASInt info; 496 char trans; 497 498 PetscFunctionBegin; 499 if (PetscDefined(USE_COMPLEX)) { 500 trans = 'C'; 501 } else { 502 trans = 'T'; 503 } 504 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 505 { /* lwork depends on the number of right-hand sides */ 506 PetscBLASInt nlfwork, lfwork = -1; 507 PetscScalar fwork; 508 509 PetscStackCallBLAS("LAPACKormqr",LAPACKormqr_("L", &trans, &m,&nrhs,&mat->rank,mat->v,&mat->lda,mat->tau,x,&ldx,&fwork,&lfwork,&info)); 510 nlfwork = (PetscBLASInt)PetscRealPart(fwork); 511 if (nlfwork > mat->lfwork) { 512 mat->lfwork = nlfwork; 513 PetscCall(PetscFree(mat->fwork)); 514 PetscCall(PetscMalloc1(mat->lfwork,&mat->fwork)); 515 } 516 } 517 PetscStackCallBLAS("LAPACKormqr",LAPACKormqr_("L", &trans, &m,&nrhs,&mat->rank,mat->v,&mat->lda,mat->tau,x,&ldx,mat->fwork,&mat->lfwork,&info)); 518 PetscCall(PetscFPTrapPop()); 519 PetscCheck(!info,PETSC_COMM_SELF,PETSC_ERR_LIB,"ORMQR - Bad orthogonal transform %d",(int)info); 520 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 521 PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U", "N", "N", &mat->rank,&nrhs,mat->v,&mat->lda,x,&ldx,&info)); 522 PetscCall(PetscFPTrapPop()); 523 PetscCheck(!info,PETSC_COMM_SELF,PETSC_ERR_LIB,"TRTRS - Bad triangular solve %d",(int)info); 524 for (PetscInt j = 0; j < nrhs; j++) { 525 for (PetscInt i = mat->rank; i < k; i++) { 526 x[j*ldx + i] = 0.; 527 } 528 } 529 PetscCall(PetscLogFlops(nrhs*(4.0*m*mat->rank - PetscSqr(mat->rank)))); 530 PetscFunctionReturn(0); 531 } 532 533 static PetscErrorCode MatSolveTranspose_SeqDense_Internal_QR(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k) 534 { 535 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 536 PetscBLASInt info; 537 538 PetscFunctionBegin; 539 if (A->rmap->n == A->cmap->n && mat->rank == A->rmap->n) { 540 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 541 PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U", "T", "N", &m,&nrhs,mat->v,&mat->lda,x,&ldx,&info)); 542 PetscCall(PetscFPTrapPop()); 543 PetscCheck(!info,PETSC_COMM_SELF,PETSC_ERR_LIB,"TRTRS - Bad triangular solve %d",(int)info); 544 if (PetscDefined(USE_COMPLEX)) PetscCall(MatConjugate_SeqDense(A)); 545 { /* lwork depends on the number of right-hand sides */ 546 PetscBLASInt nlfwork, lfwork = -1; 547 PetscScalar fwork; 548 549 PetscStackCallBLAS("LAPACKormqr",LAPACKormqr_("L", "N", &m,&nrhs,&mat->rank,mat->v,&mat->lda,mat->tau,x,&ldx,&fwork,&lfwork,&info)); 550 nlfwork = (PetscBLASInt)PetscRealPart(fwork); 551 if (nlfwork > mat->lfwork) { 552 mat->lfwork = nlfwork; 553 PetscCall(PetscFree(mat->fwork)); 554 PetscCall(PetscMalloc1(mat->lfwork,&mat->fwork)); 555 } 556 } 557 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 558 PetscStackCallBLAS("LAPACKormqr",LAPACKormqr_("L", "N", &m,&nrhs,&mat->rank,mat->v,&mat->lda,mat->tau,x,&ldx,mat->fwork,&mat->lfwork,&info)); 559 PetscCall(PetscFPTrapPop()); 560 PetscCheck(!info,PETSC_COMM_SELF,PETSC_ERR_LIB,"ORMQR - Bad orthogonal transform %d",(int)info); 561 if (PetscDefined(USE_COMPLEX)) PetscCall(MatConjugate_SeqDense(A)); 562 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"QR factored matrix cannot be used for transpose solve"); 563 PetscCall(PetscLogFlops(nrhs*(4.0*m*mat->rank - PetscSqr(mat->rank)))); 564 PetscFunctionReturn(0); 565 } 566 567 static PetscErrorCode MatSolve_SeqDense_SetUp(Mat A, Vec xx, Vec yy, PetscScalar **_y, PetscBLASInt *_m, PetscBLASInt *_k) 568 { 569 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 570 PetscScalar *y; 571 PetscBLASInt m=0, k=0; 572 573 PetscFunctionBegin; 574 PetscCall(PetscBLASIntCast(A->rmap->n,&m)); 575 PetscCall(PetscBLASIntCast(A->cmap->n,&k)); 576 if (k < m) { 577 PetscCall(VecCopy(xx, mat->qrrhs)); 578 PetscCall(VecGetArray(mat->qrrhs,&y)); 579 } else { 580 PetscCall(VecCopy(xx, yy)); 581 PetscCall(VecGetArray(yy,&y)); 582 } 583 *_y = y; 584 *_k = k; 585 *_m = m; 586 PetscFunctionReturn(0); 587 } 588 589 static PetscErrorCode MatSolve_SeqDense_TearDown(Mat A, Vec xx, Vec yy, PetscScalar **_y, PetscBLASInt *_m, PetscBLASInt *_k) 590 { 591 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 592 PetscScalar *y = NULL; 593 PetscBLASInt m, k; 594 595 PetscFunctionBegin; 596 y = *_y; 597 *_y = NULL; 598 k = *_k; 599 m = *_m; 600 if (k < m) { 601 PetscScalar *yv; 602 PetscCall(VecGetArray(yy,&yv)); 603 PetscCall(PetscArraycpy(yv, y, k)); 604 PetscCall(VecRestoreArray(yy,&yv)); 605 PetscCall(VecRestoreArray(mat->qrrhs, &y)); 606 } else { 607 PetscCall(VecRestoreArray(yy,&y)); 608 } 609 PetscFunctionReturn(0); 610 } 611 612 static PetscErrorCode MatSolve_SeqDense_LU(Mat A, Vec xx, Vec yy) 613 { 614 PetscScalar *y = NULL; 615 PetscBLASInt m = 0, k = 0; 616 617 PetscFunctionBegin; 618 PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k)); 619 PetscCall(MatSolve_SeqDense_Internal_LU(A, y, m, m, 1, k, PETSC_FALSE)); 620 PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k)); 621 PetscFunctionReturn(0); 622 } 623 624 static PetscErrorCode MatSolveTranspose_SeqDense_LU(Mat A, Vec xx, Vec yy) 625 { 626 PetscScalar *y = NULL; 627 PetscBLASInt m = 0, k = 0; 628 629 PetscFunctionBegin; 630 PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k)); 631 PetscCall(MatSolve_SeqDense_Internal_LU(A, y, m, m, 1, k, PETSC_TRUE)); 632 PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k)); 633 PetscFunctionReturn(0); 634 } 635 636 static PetscErrorCode MatSolve_SeqDense_Cholesky(Mat A, Vec xx, Vec yy) 637 { 638 PetscScalar *y = NULL; 639 PetscBLASInt m = 0, k = 0; 640 641 PetscFunctionBegin; 642 PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k)); 643 PetscCall(MatSolve_SeqDense_Internal_Cholesky(A, y, m, m, 1, k, PETSC_FALSE)); 644 PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k)); 645 PetscFunctionReturn(0); 646 } 647 648 static PetscErrorCode MatSolveTranspose_SeqDense_Cholesky(Mat A, Vec xx, Vec yy) 649 { 650 PetscScalar *y = NULL; 651 PetscBLASInt m = 0, k = 0; 652 653 PetscFunctionBegin; 654 PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k)); 655 PetscCall(MatSolve_SeqDense_Internal_Cholesky(A, y, m, m, 1, k, PETSC_TRUE)); 656 PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k)); 657 PetscFunctionReturn(0); 658 } 659 660 static PetscErrorCode MatSolve_SeqDense_QR(Mat A, Vec xx, Vec yy) 661 { 662 PetscScalar *y = NULL; 663 PetscBLASInt m = 0, k = 0; 664 665 PetscFunctionBegin; 666 PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k)); 667 PetscCall(MatSolve_SeqDense_Internal_QR(A, y, PetscMax(m,k), m, 1, k)); 668 PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k)); 669 PetscFunctionReturn(0); 670 } 671 672 static PetscErrorCode MatSolveTranspose_SeqDense_QR(Mat A, Vec xx, Vec yy) 673 { 674 PetscScalar *y = NULL; 675 PetscBLASInt m = 0, k = 0; 676 677 PetscFunctionBegin; 678 PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k)); 679 PetscCall(MatSolveTranspose_SeqDense_Internal_QR(A, y, PetscMax(m,k), m, 1, k)); 680 PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k)); 681 PetscFunctionReturn(0); 682 } 683 684 static PetscErrorCode MatMatSolve_SeqDense_SetUp(Mat A, Mat B, Mat X, PetscScalar **_y, PetscBLASInt *_ldy, PetscBLASInt *_m, PetscBLASInt *_nrhs, PetscBLASInt *_k) 685 { 686 const PetscScalar *b; 687 PetscScalar *y; 688 PetscInt n, _ldb, _ldx; 689 PetscBLASInt nrhs=0,m=0,k=0,ldb=0,ldx=0,ldy=0; 690 691 PetscFunctionBegin; 692 *_ldy=0; *_m=0; *_nrhs=0; *_k=0; *_y = NULL; 693 PetscCall(PetscBLASIntCast(A->rmap->n,&m)); 694 PetscCall(PetscBLASIntCast(A->cmap->n,&k)); 695 PetscCall(MatGetSize(B,NULL,&n)); 696 PetscCall(PetscBLASIntCast(n,&nrhs)); 697 PetscCall(MatDenseGetLDA(B,&_ldb)); 698 PetscCall(PetscBLASIntCast(_ldb, &ldb)); 699 PetscCall(MatDenseGetLDA(X,&_ldx)); 700 PetscCall(PetscBLASIntCast(_ldx, &ldx)); 701 if (ldx < m) { 702 PetscCall(MatDenseGetArrayRead(B,&b)); 703 PetscCall(PetscMalloc1(nrhs * m, &y)); 704 if (ldb == m) { 705 PetscCall(PetscArraycpy(y,b,ldb*nrhs)); 706 } else { 707 for (PetscInt j = 0; j < nrhs; j++) { 708 PetscCall(PetscArraycpy(&y[j*m],&b[j*ldb],m)); 709 } 710 } 711 ldy = m; 712 PetscCall(MatDenseRestoreArrayRead(B,&b)); 713 } else { 714 if (ldb == ldx) { 715 PetscCall(MatCopy(B, X, SAME_NONZERO_PATTERN)); 716 PetscCall(MatDenseGetArray(X,&y)); 717 } else { 718 PetscCall(MatDenseGetArray(X,&y)); 719 PetscCall(MatDenseGetArrayRead(B,&b)); 720 for (PetscInt j = 0; j < nrhs; j++) { 721 PetscCall(PetscArraycpy(&y[j*ldx],&b[j*ldb],m)); 722 } 723 PetscCall(MatDenseRestoreArrayRead(B,&b)); 724 } 725 ldy = ldx; 726 } 727 *_y = y; 728 *_ldy = ldy; 729 *_k = k; 730 *_m = m; 731 *_nrhs = nrhs; 732 PetscFunctionReturn(0); 733 } 734 735 static PetscErrorCode MatMatSolve_SeqDense_TearDown(Mat A, Mat B, Mat X, PetscScalar **_y, PetscBLASInt *_ldy, PetscBLASInt *_m, PetscBLASInt *_nrhs, PetscBLASInt *_k) 736 { 737 PetscScalar *y; 738 PetscInt _ldx; 739 PetscBLASInt k,ldy,nrhs,ldx=0; 740 741 PetscFunctionBegin; 742 y = *_y; 743 *_y = NULL; 744 k = *_k; 745 ldy = *_ldy; 746 nrhs = *_nrhs; 747 PetscCall(MatDenseGetLDA(X,&_ldx)); 748 PetscCall(PetscBLASIntCast(_ldx, &ldx)); 749 if (ldx != ldy) { 750 PetscScalar *xv; 751 PetscCall(MatDenseGetArray(X,&xv)); 752 for (PetscInt j = 0; j < nrhs; j++) { 753 PetscCall(PetscArraycpy(&xv[j*ldx],&y[j*ldy],k)); 754 } 755 PetscCall(MatDenseRestoreArray(X,&xv)); 756 PetscCall(PetscFree(y)); 757 } else { 758 PetscCall(MatDenseRestoreArray(X,&y)); 759 } 760 PetscFunctionReturn(0); 761 } 762 763 static PetscErrorCode MatMatSolve_SeqDense_LU(Mat A, Mat B, Mat X) 764 { 765 PetscScalar *y; 766 PetscBLASInt m, k, ldy, nrhs; 767 768 PetscFunctionBegin; 769 PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k)); 770 PetscCall(MatSolve_SeqDense_Internal_LU(A, y, ldy, m, nrhs, k, PETSC_FALSE)); 771 PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k)); 772 PetscFunctionReturn(0); 773 } 774 775 static PetscErrorCode MatMatSolveTranspose_SeqDense_LU(Mat A, Mat B, Mat X) 776 { 777 PetscScalar *y; 778 PetscBLASInt m, k, ldy, nrhs; 779 780 PetscFunctionBegin; 781 PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k)); 782 PetscCall(MatSolve_SeqDense_Internal_LU(A, y, ldy, m, nrhs, k, PETSC_TRUE)); 783 PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k)); 784 PetscFunctionReturn(0); 785 } 786 787 static PetscErrorCode MatMatSolve_SeqDense_Cholesky(Mat A, Mat B, Mat X) 788 { 789 PetscScalar *y; 790 PetscBLASInt m, k, ldy, nrhs; 791 792 PetscFunctionBegin; 793 PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k)); 794 PetscCall(MatSolve_SeqDense_Internal_Cholesky(A, y, ldy, m, nrhs, k, PETSC_FALSE)); 795 PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k)); 796 PetscFunctionReturn(0); 797 } 798 799 static PetscErrorCode MatMatSolveTranspose_SeqDense_Cholesky(Mat A, Mat B, Mat X) 800 { 801 PetscScalar *y; 802 PetscBLASInt m, k, ldy, nrhs; 803 804 PetscFunctionBegin; 805 PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k)); 806 PetscCall(MatSolve_SeqDense_Internal_Cholesky(A, y, ldy, m, nrhs, k, PETSC_TRUE)); 807 PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k)); 808 PetscFunctionReturn(0); 809 } 810 811 static PetscErrorCode MatMatSolve_SeqDense_QR(Mat A, Mat B, Mat X) 812 { 813 PetscScalar *y; 814 PetscBLASInt m, k, ldy, nrhs; 815 816 PetscFunctionBegin; 817 PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k)); 818 PetscCall(MatSolve_SeqDense_Internal_QR(A, y, ldy, m, nrhs, k)); 819 PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k)); 820 PetscFunctionReturn(0); 821 } 822 823 static PetscErrorCode MatMatSolveTranspose_SeqDense_QR(Mat A, Mat B, Mat X) 824 { 825 PetscScalar *y; 826 PetscBLASInt m, k, ldy, nrhs; 827 828 PetscFunctionBegin; 829 PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k)); 830 PetscCall(MatSolveTranspose_SeqDense_Internal_QR(A, y, ldy, m, nrhs, k)); 831 PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k)); 832 PetscFunctionReturn(0); 833 } 834 835 /* ---------------------------------------------------------------*/ 836 /* COMMENT: I have chosen to hide row permutation in the pivots, 837 rather than put it in the Mat->row slot.*/ 838 PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo) 839 { 840 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 841 PetscBLASInt n,m,info; 842 843 PetscFunctionBegin; 844 PetscCall(PetscBLASIntCast(A->cmap->n,&n)); 845 PetscCall(PetscBLASIntCast(A->rmap->n,&m)); 846 if (!mat->pivots) { 847 PetscCall(PetscMalloc1(A->rmap->n,&mat->pivots)); 848 PetscCall(PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt))); 849 } 850 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 851 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 852 PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info)); 853 PetscCall(PetscFPTrapPop()); 854 855 PetscCheck(info>=0,PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization %d",(int)info); 856 PetscCheck(info<=0,PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization %d",(int)info); 857 858 A->ops->solve = MatSolve_SeqDense_LU; 859 A->ops->matsolve = MatMatSolve_SeqDense_LU; 860 A->ops->solvetranspose = MatSolveTranspose_SeqDense_LU; 861 A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_LU; 862 A->factortype = MAT_FACTOR_LU; 863 864 PetscCall(PetscFree(A->solvertype)); 865 PetscCall(PetscStrallocpy(MATSOLVERPETSC,&A->solvertype)); 866 867 PetscCall(PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3)); 868 PetscFunctionReturn(0); 869 } 870 871 static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy) 872 { 873 MatFactorInfo info; 874 875 PetscFunctionBegin; 876 PetscCall(MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES)); 877 PetscCall((*fact->ops->lufactor)(fact,NULL,NULL,&info)); 878 PetscFunctionReturn(0); 879 } 880 881 PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info) 882 { 883 PetscFunctionBegin; 884 fact->preallocated = PETSC_TRUE; 885 fact->assembled = PETSC_TRUE; 886 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense; 887 PetscFunctionReturn(0); 888 } 889 890 /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */ 891 PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo) 892 { 893 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 894 PetscBLASInt info,n; 895 896 PetscFunctionBegin; 897 PetscCall(PetscBLASIntCast(A->cmap->n,&n)); 898 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 899 if (A->spd) { 900 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 901 PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info)); 902 PetscCall(PetscFPTrapPop()); 903 #if defined(PETSC_USE_COMPLEX) 904 } else if (A->hermitian) { 905 if (!mat->pivots) { 906 PetscCall(PetscMalloc1(A->rmap->n,&mat->pivots)); 907 PetscCall(PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt))); 908 } 909 if (!mat->fwork) { 910 PetscScalar dummy; 911 912 mat->lfwork = -1; 913 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 914 PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info)); 915 PetscCall(PetscFPTrapPop()); 916 mat->lfwork = (PetscInt)PetscRealPart(dummy); 917 PetscCall(PetscMalloc1(mat->lfwork,&mat->fwork)); 918 PetscCall(PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt))); 919 } 920 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 921 PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info)); 922 PetscCall(PetscFPTrapPop()); 923 #endif 924 } else { /* symmetric case */ 925 if (!mat->pivots) { 926 PetscCall(PetscMalloc1(A->rmap->n,&mat->pivots)); 927 PetscCall(PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt))); 928 } 929 if (!mat->fwork) { 930 PetscScalar dummy; 931 932 mat->lfwork = -1; 933 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 934 PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info)); 935 PetscCall(PetscFPTrapPop()); 936 mat->lfwork = (PetscInt)PetscRealPart(dummy); 937 PetscCall(PetscMalloc1(mat->lfwork,&mat->fwork)); 938 PetscCall(PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt))); 939 } 940 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 941 PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info)); 942 PetscCall(PetscFPTrapPop()); 943 } 944 PetscCheck(!info,PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %" PetscInt_FMT,(PetscInt)info-1); 945 946 A->ops->solve = MatSolve_SeqDense_Cholesky; 947 A->ops->matsolve = MatMatSolve_SeqDense_Cholesky; 948 A->ops->solvetranspose = MatSolveTranspose_SeqDense_Cholesky; 949 A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_Cholesky; 950 A->factortype = MAT_FACTOR_CHOLESKY; 951 952 PetscCall(PetscFree(A->solvertype)); 953 PetscCall(PetscStrallocpy(MATSOLVERPETSC,&A->solvertype)); 954 955 PetscCall(PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0)); 956 PetscFunctionReturn(0); 957 } 958 959 static PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy) 960 { 961 MatFactorInfo info; 962 963 PetscFunctionBegin; 964 info.fill = 1.0; 965 966 PetscCall(MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES)); 967 PetscCall((*fact->ops->choleskyfactor)(fact,NULL,&info)); 968 PetscFunctionReturn(0); 969 } 970 971 PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info) 972 { 973 PetscFunctionBegin; 974 fact->assembled = PETSC_TRUE; 975 fact->preallocated = PETSC_TRUE; 976 fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense; 977 PetscFunctionReturn(0); 978 } 979 980 PetscErrorCode MatQRFactor_SeqDense(Mat A,IS col,const MatFactorInfo *minfo) 981 { 982 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 983 PetscBLASInt n,m,info, min, max; 984 985 PetscFunctionBegin; 986 PetscCall(PetscBLASIntCast(A->cmap->n,&n)); 987 PetscCall(PetscBLASIntCast(A->rmap->n,&m)); 988 max = PetscMax(m, n); 989 min = PetscMin(m, n); 990 if (!mat->tau) { 991 PetscCall(PetscMalloc1(min,&mat->tau)); 992 PetscCall(PetscLogObjectMemory((PetscObject)A,min*sizeof(PetscScalar))); 993 } 994 if (!mat->pivots) { 995 PetscCall(PetscMalloc1(n,&mat->pivots)); 996 PetscCall(PetscLogObjectMemory((PetscObject)A,n*sizeof(PetscScalar))); 997 } 998 if (!mat->qrrhs) { 999 PetscCall(MatCreateVecs(A, NULL, &(mat->qrrhs))); 1000 } 1001 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 1002 if (!mat->fwork) { 1003 PetscScalar dummy; 1004 1005 mat->lfwork = -1; 1006 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 1007 PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&m,&n,mat->v,&mat->lda,mat->tau,&dummy,&mat->lfwork,&info)); 1008 PetscCall(PetscFPTrapPop()); 1009 mat->lfwork = (PetscInt)PetscRealPart(dummy); 1010 PetscCall(PetscMalloc1(mat->lfwork,&mat->fwork)); 1011 PetscCall(PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt))); 1012 } 1013 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 1014 PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&m,&n,mat->v,&mat->lda,mat->tau,mat->fwork,&mat->lfwork,&info)); 1015 PetscCall(PetscFPTrapPop()); 1016 PetscCheck(!info,PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to QR factorization %d",(int)info); 1017 // TODO: try to estimate rank or test for and use geqp3 for rank revealing QR. For now just say rank is min of m and n 1018 mat->rank = min; 1019 1020 A->ops->solve = MatSolve_SeqDense_QR; 1021 A->ops->matsolve = MatMatSolve_SeqDense_QR; 1022 A->factortype = MAT_FACTOR_QR; 1023 if (m == n) { 1024 A->ops->solvetranspose = MatSolveTranspose_SeqDense_QR; 1025 A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_QR; 1026 } 1027 1028 PetscCall(PetscFree(A->solvertype)); 1029 PetscCall(PetscStrallocpy(MATSOLVERPETSC,&A->solvertype)); 1030 1031 PetscCall(PetscLogFlops(2.0*min*min*(max-min/3.0))); 1032 PetscFunctionReturn(0); 1033 } 1034 1035 static PetscErrorCode MatQRFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy) 1036 { 1037 MatFactorInfo info; 1038 1039 PetscFunctionBegin; 1040 info.fill = 1.0; 1041 1042 PetscCall(MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES)); 1043 PetscUseMethod(fact,"MatQRFactor_C",(Mat,IS,const MatFactorInfo *),(fact,NULL,&info)); 1044 PetscFunctionReturn(0); 1045 } 1046 1047 PetscErrorCode MatQRFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info) 1048 { 1049 PetscFunctionBegin; 1050 fact->assembled = PETSC_TRUE; 1051 fact->preallocated = PETSC_TRUE; 1052 PetscCall(PetscObjectComposeFunction((PetscObject)fact,"MatQRFactorNumeric_C",MatQRFactorNumeric_SeqDense)); 1053 PetscFunctionReturn(0); 1054 } 1055 1056 /* uses LAPACK */ 1057 PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact) 1058 { 1059 PetscFunctionBegin; 1060 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),fact)); 1061 PetscCall(MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n)); 1062 PetscCall(MatSetType(*fact,MATDENSE)); 1063 (*fact)->trivialsymbolic = PETSC_TRUE; 1064 if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) { 1065 (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense; 1066 (*fact)->ops->ilufactorsymbolic = MatLUFactorSymbolic_SeqDense; 1067 } else if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1068 (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense; 1069 } else if (ftype == MAT_FACTOR_QR) { 1070 PetscCall(PetscObjectComposeFunction((PetscObject)(*fact),"MatQRFactorSymbolic_C",MatQRFactorSymbolic_SeqDense)); 1071 } 1072 (*fact)->factortype = ftype; 1073 1074 PetscCall(PetscFree((*fact)->solvertype)); 1075 PetscCall(PetscStrallocpy(MATSOLVERPETSC,&(*fact)->solvertype)); 1076 PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&(*fact)->preferredordering[MAT_FACTOR_LU])); 1077 PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&(*fact)->preferredordering[MAT_FACTOR_ILU])); 1078 PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&(*fact)->preferredordering[MAT_FACTOR_CHOLESKY])); 1079 PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&(*fact)->preferredordering[MAT_FACTOR_ICC])); 1080 PetscFunctionReturn(0); 1081 } 1082 1083 /* ------------------------------------------------------------------*/ 1084 static PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx) 1085 { 1086 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1087 PetscScalar *x,*v = mat->v,zero = 0.0,xt; 1088 const PetscScalar *b; 1089 PetscInt m = A->rmap->n,i; 1090 PetscBLASInt o = 1,bm = 0; 1091 1092 PetscFunctionBegin; 1093 #if defined(PETSC_HAVE_CUDA) 1094 PetscCheck(A->offloadmask != PETSC_OFFLOAD_GPU,PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented"); 1095 #endif 1096 if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */ 1097 PetscCall(PetscBLASIntCast(m,&bm)); 1098 if (flag & SOR_ZERO_INITIAL_GUESS) { 1099 /* this is a hack fix, should have another version without the second BLASdotu */ 1100 PetscCall(VecSet(xx,zero)); 1101 } 1102 PetscCall(VecGetArray(xx,&x)); 1103 PetscCall(VecGetArrayRead(bb,&b)); 1104 its = its*lits; 1105 PetscCheck(its > 0,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %" PetscInt_FMT " and local its %" PetscInt_FMT " both positive",its,lits); 1106 while (its--) { 1107 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 1108 for (i=0; i<m; i++) { 1109 PetscStackCallBLAS("BLASdotu",xt = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o)); 1110 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 1111 } 1112 } 1113 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 1114 for (i=m-1; i>=0; i--) { 1115 PetscStackCallBLAS("BLASdotu",xt = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o)); 1116 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 1117 } 1118 } 1119 } 1120 PetscCall(VecRestoreArrayRead(bb,&b)); 1121 PetscCall(VecRestoreArray(xx,&x)); 1122 PetscFunctionReturn(0); 1123 } 1124 1125 /* -----------------------------------------------------------------*/ 1126 PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy) 1127 { 1128 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1129 const PetscScalar *v = mat->v,*x; 1130 PetscScalar *y; 1131 PetscBLASInt m, n,_One=1; 1132 PetscScalar _DOne=1.0,_DZero=0.0; 1133 1134 PetscFunctionBegin; 1135 PetscCall(PetscBLASIntCast(A->rmap->n,&m)); 1136 PetscCall(PetscBLASIntCast(A->cmap->n,&n)); 1137 PetscCall(VecGetArrayRead(xx,&x)); 1138 PetscCall(VecGetArrayWrite(yy,&y)); 1139 if (!A->rmap->n || !A->cmap->n) { 1140 PetscBLASInt i; 1141 for (i=0; i<n; i++) y[i] = 0.0; 1142 } else { 1143 PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One)); 1144 PetscCall(PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n)); 1145 } 1146 PetscCall(VecRestoreArrayRead(xx,&x)); 1147 PetscCall(VecRestoreArrayWrite(yy,&y)); 1148 PetscFunctionReturn(0); 1149 } 1150 1151 PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy) 1152 { 1153 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1154 PetscScalar *y,_DOne=1.0,_DZero=0.0; 1155 PetscBLASInt m, n, _One=1; 1156 const PetscScalar *v = mat->v,*x; 1157 1158 PetscFunctionBegin; 1159 PetscCall(PetscBLASIntCast(A->rmap->n,&m)); 1160 PetscCall(PetscBLASIntCast(A->cmap->n,&n)); 1161 PetscCall(VecGetArrayRead(xx,&x)); 1162 PetscCall(VecGetArrayWrite(yy,&y)); 1163 if (!A->rmap->n || !A->cmap->n) { 1164 PetscBLASInt i; 1165 for (i=0; i<m; i++) y[i] = 0.0; 1166 } else { 1167 PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One)); 1168 PetscCall(PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n)); 1169 } 1170 PetscCall(VecRestoreArrayRead(xx,&x)); 1171 PetscCall(VecRestoreArrayWrite(yy,&y)); 1172 PetscFunctionReturn(0); 1173 } 1174 1175 PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 1176 { 1177 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1178 const PetscScalar *v = mat->v,*x; 1179 PetscScalar *y,_DOne=1.0; 1180 PetscBLASInt m, n, _One=1; 1181 1182 PetscFunctionBegin; 1183 PetscCall(PetscBLASIntCast(A->rmap->n,&m)); 1184 PetscCall(PetscBLASIntCast(A->cmap->n,&n)); 1185 PetscCall(VecCopy(zz,yy)); 1186 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 1187 PetscCall(VecGetArrayRead(xx,&x)); 1188 PetscCall(VecGetArray(yy,&y)); 1189 PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One)); 1190 PetscCall(VecRestoreArrayRead(xx,&x)); 1191 PetscCall(VecRestoreArray(yy,&y)); 1192 PetscCall(PetscLogFlops(2.0*A->rmap->n*A->cmap->n)); 1193 PetscFunctionReturn(0); 1194 } 1195 1196 PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 1197 { 1198 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1199 const PetscScalar *v = mat->v,*x; 1200 PetscScalar *y; 1201 PetscBLASInt m, n, _One=1; 1202 PetscScalar _DOne=1.0; 1203 1204 PetscFunctionBegin; 1205 PetscCall(PetscBLASIntCast(A->rmap->n,&m)); 1206 PetscCall(PetscBLASIntCast(A->cmap->n,&n)); 1207 PetscCall(VecCopy(zz,yy)); 1208 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 1209 PetscCall(VecGetArrayRead(xx,&x)); 1210 PetscCall(VecGetArray(yy,&y)); 1211 PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One)); 1212 PetscCall(VecRestoreArrayRead(xx,&x)); 1213 PetscCall(VecRestoreArray(yy,&y)); 1214 PetscCall(PetscLogFlops(2.0*A->rmap->n*A->cmap->n)); 1215 PetscFunctionReturn(0); 1216 } 1217 1218 /* -----------------------------------------------------------------*/ 1219 static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 1220 { 1221 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1222 PetscInt i; 1223 1224 PetscFunctionBegin; 1225 *ncols = A->cmap->n; 1226 if (cols) { 1227 PetscCall(PetscMalloc1(A->cmap->n,cols)); 1228 for (i=0; i<A->cmap->n; i++) (*cols)[i] = i; 1229 } 1230 if (vals) { 1231 const PetscScalar *v; 1232 1233 PetscCall(MatDenseGetArrayRead(A,&v)); 1234 PetscCall(PetscMalloc1(A->cmap->n,vals)); 1235 v += row; 1236 for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;} 1237 PetscCall(MatDenseRestoreArrayRead(A,&v)); 1238 } 1239 PetscFunctionReturn(0); 1240 } 1241 1242 static PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 1243 { 1244 PetscFunctionBegin; 1245 if (ncols) *ncols = 0; 1246 if (cols) PetscCall(PetscFree(*cols)); 1247 if (vals) PetscCall(PetscFree(*vals)); 1248 PetscFunctionReturn(0); 1249 } 1250 /* ----------------------------------------------------------------*/ 1251 static PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv) 1252 { 1253 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1254 PetscScalar *av; 1255 PetscInt i,j,idx=0; 1256 #if defined(PETSC_HAVE_CUDA) 1257 PetscOffloadMask oldf; 1258 #endif 1259 1260 PetscFunctionBegin; 1261 PetscCall(MatDenseGetArray(A,&av)); 1262 if (!mat->roworiented) { 1263 if (addv == INSERT_VALUES) { 1264 for (j=0; j<n; j++) { 1265 if (indexn[j] < 0) {idx += m; continue;} 1266 PetscCheck(indexn[j] < A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT,indexn[j],A->cmap->n-1); 1267 for (i=0; i<m; i++) { 1268 if (indexm[i] < 0) {idx++; continue;} 1269 PetscCheck(indexm[i] < A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT,indexm[i],A->rmap->n-1); 1270 av[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 1271 } 1272 } 1273 } else { 1274 for (j=0; j<n; j++) { 1275 if (indexn[j] < 0) {idx += m; continue;} 1276 PetscCheck(indexn[j] < A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT,indexn[j],A->cmap->n-1); 1277 for (i=0; i<m; i++) { 1278 if (indexm[i] < 0) {idx++; continue;} 1279 PetscCheck(indexm[i] < A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT,indexm[i],A->rmap->n-1); 1280 av[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 1281 } 1282 } 1283 } 1284 } else { 1285 if (addv == INSERT_VALUES) { 1286 for (i=0; i<m; i++) { 1287 if (indexm[i] < 0) { idx += n; continue;} 1288 PetscCheck(indexm[i] < A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT,indexm[i],A->rmap->n-1); 1289 for (j=0; j<n; j++) { 1290 if (indexn[j] < 0) { idx++; continue;} 1291 PetscCheck(indexn[j] < A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT,indexn[j],A->cmap->n-1); 1292 av[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 1293 } 1294 } 1295 } else { 1296 for (i=0; i<m; i++) { 1297 if (indexm[i] < 0) { idx += n; continue;} 1298 PetscCheck(indexm[i] < A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT,indexm[i],A->rmap->n-1); 1299 for (j=0; j<n; j++) { 1300 if (indexn[j] < 0) { idx++; continue;} 1301 PetscCheck(indexn[j] < A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT,indexn[j],A->cmap->n-1); 1302 av[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 1303 } 1304 } 1305 } 1306 } 1307 /* hack to prevent unneeded copy to the GPU while returning the array */ 1308 #if defined(PETSC_HAVE_CUDA) 1309 oldf = A->offloadmask; 1310 A->offloadmask = PETSC_OFFLOAD_GPU; 1311 #endif 1312 PetscCall(MatDenseRestoreArray(A,&av)); 1313 #if defined(PETSC_HAVE_CUDA) 1314 A->offloadmask = (oldf == PETSC_OFFLOAD_UNALLOCATED ? PETSC_OFFLOAD_UNALLOCATED : PETSC_OFFLOAD_CPU); 1315 #endif 1316 PetscFunctionReturn(0); 1317 } 1318 1319 static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[]) 1320 { 1321 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1322 const PetscScalar *vv; 1323 PetscInt i,j; 1324 1325 PetscFunctionBegin; 1326 PetscCall(MatDenseGetArrayRead(A,&vv)); 1327 /* row-oriented output */ 1328 for (i=0; i<m; i++) { 1329 if (indexm[i] < 0) {v += n;continue;} 1330 PetscCheck(indexm[i] < A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %" PetscInt_FMT " requested larger than number rows %" PetscInt_FMT,indexm[i],A->rmap->n); 1331 for (j=0; j<n; j++) { 1332 if (indexn[j] < 0) {v++; continue;} 1333 PetscCheck(indexn[j] < A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column %" PetscInt_FMT " requested larger than number columns %" PetscInt_FMT,indexn[j],A->cmap->n); 1334 *v++ = vv[indexn[j]*mat->lda + indexm[i]]; 1335 } 1336 } 1337 PetscCall(MatDenseRestoreArrayRead(A,&vv)); 1338 PetscFunctionReturn(0); 1339 } 1340 1341 /* -----------------------------------------------------------------*/ 1342 1343 PetscErrorCode MatView_Dense_Binary(Mat mat,PetscViewer viewer) 1344 { 1345 PetscBool skipHeader; 1346 PetscViewerFormat format; 1347 PetscInt header[4],M,N,m,lda,i,j,k; 1348 const PetscScalar *v; 1349 PetscScalar *vwork; 1350 1351 PetscFunctionBegin; 1352 PetscCall(PetscViewerSetUp(viewer)); 1353 PetscCall(PetscViewerBinaryGetSkipHeader(viewer,&skipHeader)); 1354 PetscCall(PetscViewerGetFormat(viewer,&format)); 1355 if (skipHeader) format = PETSC_VIEWER_NATIVE; 1356 1357 PetscCall(MatGetSize(mat,&M,&N)); 1358 1359 /* write matrix header */ 1360 header[0] = MAT_FILE_CLASSID; header[1] = M; header[2] = N; 1361 header[3] = (format == PETSC_VIEWER_NATIVE) ? MATRIX_BINARY_FORMAT_DENSE : M*N; 1362 if (!skipHeader) PetscCall(PetscViewerBinaryWrite(viewer,header,4,PETSC_INT)); 1363 1364 PetscCall(MatGetLocalSize(mat,&m,NULL)); 1365 if (format != PETSC_VIEWER_NATIVE) { 1366 PetscInt nnz = m*N, *iwork; 1367 /* store row lengths for each row */ 1368 PetscCall(PetscMalloc1(nnz,&iwork)); 1369 for (i=0; i<m; i++) iwork[i] = N; 1370 PetscCall(PetscViewerBinaryWriteAll(viewer,iwork,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT)); 1371 /* store column indices (zero start index) */ 1372 for (k=0, i=0; i<m; i++) 1373 for (j=0; j<N; j++, k++) 1374 iwork[k] = j; 1375 PetscCall(PetscViewerBinaryWriteAll(viewer,iwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT)); 1376 PetscCall(PetscFree(iwork)); 1377 } 1378 /* store matrix values as a dense matrix in row major order */ 1379 PetscCall(PetscMalloc1(m*N,&vwork)); 1380 PetscCall(MatDenseGetArrayRead(mat,&v)); 1381 PetscCall(MatDenseGetLDA(mat,&lda)); 1382 for (k=0, i=0; i<m; i++) 1383 for (j=0; j<N; j++, k++) 1384 vwork[k] = v[i+lda*j]; 1385 PetscCall(MatDenseRestoreArrayRead(mat,&v)); 1386 PetscCall(PetscViewerBinaryWriteAll(viewer,vwork,m*N,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR)); 1387 PetscCall(PetscFree(vwork)); 1388 PetscFunctionReturn(0); 1389 } 1390 1391 PetscErrorCode MatLoad_Dense_Binary(Mat mat,PetscViewer viewer) 1392 { 1393 PetscBool skipHeader; 1394 PetscInt header[4],M,N,m,nz,lda,i,j,k; 1395 PetscInt rows,cols; 1396 PetscScalar *v,*vwork; 1397 1398 PetscFunctionBegin; 1399 PetscCall(PetscViewerSetUp(viewer)); 1400 PetscCall(PetscViewerBinaryGetSkipHeader(viewer,&skipHeader)); 1401 1402 if (!skipHeader) { 1403 PetscCall(PetscViewerBinaryRead(viewer,header,4,NULL,PETSC_INT)); 1404 PetscCheck(header[0] == MAT_FILE_CLASSID,PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Not a matrix object in file"); 1405 M = header[1]; N = header[2]; 1406 PetscCheck(M >= 0,PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix row size (%" PetscInt_FMT ") in file is negative",M); 1407 PetscCheck(N >= 0,PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix column size (%" PetscInt_FMT ") in file is negative",N); 1408 nz = header[3]; 1409 PetscCheck(nz == MATRIX_BINARY_FORMAT_DENSE || nz >= 0,PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Unknown matrix format %" PetscInt_FMT " in file",nz); 1410 } else { 1411 PetscCall(MatGetSize(mat,&M,&N)); 1412 PetscCheck(M >= 0 && N >= 0,PETSC_COMM_SELF,PETSC_ERR_USER,"Matrix binary file header was skipped, thus the user must specify the global sizes of input matrix"); 1413 nz = MATRIX_BINARY_FORMAT_DENSE; 1414 } 1415 1416 /* setup global sizes if not set */ 1417 if (mat->rmap->N < 0) mat->rmap->N = M; 1418 if (mat->cmap->N < 0) mat->cmap->N = N; 1419 PetscCall(MatSetUp(mat)); 1420 /* check if global sizes are correct */ 1421 PetscCall(MatGetSize(mat,&rows,&cols)); 1422 PetscCheck(M == rows && N == cols,PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different sizes (%" PetscInt_FMT ", %" PetscInt_FMT ") than the input matrix (%" PetscInt_FMT ", %" PetscInt_FMT ")",M,N,rows,cols); 1423 1424 PetscCall(MatGetSize(mat,NULL,&N)); 1425 PetscCall(MatGetLocalSize(mat,&m,NULL)); 1426 PetscCall(MatDenseGetArray(mat,&v)); 1427 PetscCall(MatDenseGetLDA(mat,&lda)); 1428 if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense format */ 1429 PetscInt nnz = m*N; 1430 /* read in matrix values */ 1431 PetscCall(PetscMalloc1(nnz,&vwork)); 1432 PetscCall(PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR)); 1433 /* store values in column major order */ 1434 for (j=0; j<N; j++) 1435 for (i=0; i<m; i++) 1436 v[i+lda*j] = vwork[i*N+j]; 1437 PetscCall(PetscFree(vwork)); 1438 } else { /* matrix in file is sparse format */ 1439 PetscInt nnz = 0, *rlens, *icols; 1440 /* read in row lengths */ 1441 PetscCall(PetscMalloc1(m,&rlens)); 1442 PetscCall(PetscViewerBinaryReadAll(viewer,rlens,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT)); 1443 for (i=0; i<m; i++) nnz += rlens[i]; 1444 /* read in column indices and values */ 1445 PetscCall(PetscMalloc2(nnz,&icols,nnz,&vwork)); 1446 PetscCall(PetscViewerBinaryReadAll(viewer,icols,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT)); 1447 PetscCall(PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR)); 1448 /* store values in column major order */ 1449 for (k=0, i=0; i<m; i++) 1450 for (j=0; j<rlens[i]; j++, k++) 1451 v[i+lda*icols[k]] = vwork[k]; 1452 PetscCall(PetscFree(rlens)); 1453 PetscCall(PetscFree2(icols,vwork)); 1454 } 1455 PetscCall(MatDenseRestoreArray(mat,&v)); 1456 PetscCall(MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY)); 1457 PetscCall(MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY)); 1458 PetscFunctionReturn(0); 1459 } 1460 1461 PetscErrorCode MatLoad_SeqDense(Mat newMat, PetscViewer viewer) 1462 { 1463 PetscBool isbinary, ishdf5; 1464 1465 PetscFunctionBegin; 1466 PetscValidHeaderSpecific(newMat,MAT_CLASSID,1); 1467 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 1468 /* force binary viewer to load .info file if it has not yet done so */ 1469 PetscCall(PetscViewerSetUp(viewer)); 1470 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary)); 1471 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5, &ishdf5)); 1472 if (isbinary) { 1473 PetscCall(MatLoad_Dense_Binary(newMat,viewer)); 1474 } else if (ishdf5) { 1475 #if defined(PETSC_HAVE_HDF5) 1476 PetscCall(MatLoad_Dense_HDF5(newMat,viewer)); 1477 #else 1478 SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5"); 1479 #endif 1480 } else { 1481 SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"Viewer type %s not yet supported for reading %s matrices",((PetscObject)viewer)->type_name,((PetscObject)newMat)->type_name); 1482 } 1483 PetscFunctionReturn(0); 1484 } 1485 1486 static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer) 1487 { 1488 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1489 PetscInt i,j; 1490 const char *name; 1491 PetscScalar *v,*av; 1492 PetscViewerFormat format; 1493 #if defined(PETSC_USE_COMPLEX) 1494 PetscBool allreal = PETSC_TRUE; 1495 #endif 1496 1497 PetscFunctionBegin; 1498 PetscCall(MatDenseGetArrayRead(A,(const PetscScalar**)&av)); 1499 PetscCall(PetscViewerGetFormat(viewer,&format)); 1500 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 1501 PetscFunctionReturn(0); /* do nothing for now */ 1502 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 1503 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 1504 for (i=0; i<A->rmap->n; i++) { 1505 v = av + i; 1506 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i)); 1507 for (j=0; j<A->cmap->n; j++) { 1508 #if defined(PETSC_USE_COMPLEX) 1509 if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) { 1510 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v))); 1511 } else if (PetscRealPart(*v)) { 1512 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",j,(double)PetscRealPart(*v))); 1513 } 1514 #else 1515 if (*v) { 1516 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",j,(double)*v)); 1517 } 1518 #endif 1519 v += a->lda; 1520 } 1521 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 1522 } 1523 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 1524 } else { 1525 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 1526 #if defined(PETSC_USE_COMPLEX) 1527 /* determine if matrix has all real values */ 1528 for (j=0; j<A->cmap->n; j++) { 1529 v = av + j*a->lda; 1530 for (i=0; i<A->rmap->n; i++) { 1531 if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;} 1532 } 1533 } 1534 #endif 1535 if (format == PETSC_VIEWER_ASCII_MATLAB) { 1536 PetscCall(PetscObjectGetName((PetscObject)A,&name)); 1537 PetscCall(PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT " \n",A->rmap->n,A->cmap->n)); 1538 PetscCall(PetscViewerASCIIPrintf(viewer,"%s = zeros(%" PetscInt_FMT ",%" PetscInt_FMT ");\n",name,A->rmap->n,A->cmap->n)); 1539 PetscCall(PetscViewerASCIIPrintf(viewer,"%s = [\n",name)); 1540 } 1541 1542 for (i=0; i<A->rmap->n; i++) { 1543 v = av + i; 1544 for (j=0; j<A->cmap->n; j++) { 1545 #if defined(PETSC_USE_COMPLEX) 1546 if (allreal) { 1547 PetscCall(PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v))); 1548 } else { 1549 PetscCall(PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v))); 1550 } 1551 #else 1552 PetscCall(PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v)); 1553 #endif 1554 v += a->lda; 1555 } 1556 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 1557 } 1558 if (format == PETSC_VIEWER_ASCII_MATLAB) { 1559 PetscCall(PetscViewerASCIIPrintf(viewer,"];\n")); 1560 } 1561 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 1562 } 1563 PetscCall(MatDenseRestoreArrayRead(A,(const PetscScalar**)&av)); 1564 PetscCall(PetscViewerFlush(viewer)); 1565 PetscFunctionReturn(0); 1566 } 1567 1568 #include <petscdraw.h> 1569 static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa) 1570 { 1571 Mat A = (Mat) Aa; 1572 PetscInt m = A->rmap->n,n = A->cmap->n,i,j; 1573 int color = PETSC_DRAW_WHITE; 1574 const PetscScalar *v; 1575 PetscViewer viewer; 1576 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 1577 PetscViewerFormat format; 1578 1579 PetscFunctionBegin; 1580 PetscCall(PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer)); 1581 PetscCall(PetscViewerGetFormat(viewer,&format)); 1582 PetscCall(PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr)); 1583 1584 /* Loop over matrix elements drawing boxes */ 1585 PetscCall(MatDenseGetArrayRead(A,&v)); 1586 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 1587 PetscDrawCollectiveBegin(draw); 1588 /* Blue for negative and Red for positive */ 1589 for (j = 0; j < n; j++) { 1590 x_l = j; x_r = x_l + 1.0; 1591 for (i = 0; i < m; i++) { 1592 y_l = m - i - 1.0; 1593 y_r = y_l + 1.0; 1594 if (PetscRealPart(v[j*m+i]) > 0.) color = PETSC_DRAW_RED; 1595 else if (PetscRealPart(v[j*m+i]) < 0.) color = PETSC_DRAW_BLUE; 1596 else continue; 1597 PetscCall(PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color)); 1598 } 1599 } 1600 PetscDrawCollectiveEnd(draw); 1601 } else { 1602 /* use contour shading to indicate magnitude of values */ 1603 /* first determine max of all nonzero values */ 1604 PetscReal minv = 0.0, maxv = 0.0; 1605 PetscDraw popup; 1606 1607 for (i=0; i < m*n; i++) { 1608 if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]); 1609 } 1610 if (minv >= maxv) maxv = minv + PETSC_SMALL; 1611 PetscCall(PetscDrawGetPopup(draw,&popup)); 1612 PetscCall(PetscDrawScalePopup(popup,minv,maxv)); 1613 1614 PetscDrawCollectiveBegin(draw); 1615 for (j=0; j<n; j++) { 1616 x_l = j; 1617 x_r = x_l + 1.0; 1618 for (i=0; i<m; i++) { 1619 y_l = m - i - 1.0; 1620 y_r = y_l + 1.0; 1621 color = PetscDrawRealToColor(PetscAbsScalar(v[j*m+i]),minv,maxv); 1622 PetscCall(PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color)); 1623 } 1624 } 1625 PetscDrawCollectiveEnd(draw); 1626 } 1627 PetscCall(MatDenseRestoreArrayRead(A,&v)); 1628 PetscFunctionReturn(0); 1629 } 1630 1631 static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer) 1632 { 1633 PetscDraw draw; 1634 PetscBool isnull; 1635 PetscReal xr,yr,xl,yl,h,w; 1636 1637 PetscFunctionBegin; 1638 PetscCall(PetscViewerDrawGetDraw(viewer,0,&draw)); 1639 PetscCall(PetscDrawIsNull(draw,&isnull)); 1640 if (isnull) PetscFunctionReturn(0); 1641 1642 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 1643 xr += w; yr += h; xl = -w; yl = -h; 1644 PetscCall(PetscDrawSetCoordinates(draw,xl,yl,xr,yr)); 1645 PetscCall(PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer)); 1646 PetscCall(PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A)); 1647 PetscCall(PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL)); 1648 PetscCall(PetscDrawSave(draw)); 1649 PetscFunctionReturn(0); 1650 } 1651 1652 PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer) 1653 { 1654 PetscBool iascii,isbinary,isdraw; 1655 1656 PetscFunctionBegin; 1657 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 1658 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary)); 1659 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw)); 1660 if (iascii) PetscCall(MatView_SeqDense_ASCII(A,viewer)); 1661 else if (isbinary) PetscCall(MatView_Dense_Binary(A,viewer)); 1662 else if (isdraw) PetscCall(MatView_SeqDense_Draw(A,viewer)); 1663 PetscFunctionReturn(0); 1664 } 1665 1666 static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A,const PetscScalar *array) 1667 { 1668 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1669 1670 PetscFunctionBegin; 1671 PetscCheck(!a->vecinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 1672 PetscCheck(!a->matinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 1673 PetscCheck(!a->unplacedarray,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreArray() first"); 1674 a->unplacedarray = a->v; 1675 a->unplaced_user_alloc = a->user_alloc; 1676 a->v = (PetscScalar*) array; 1677 a->user_alloc = PETSC_TRUE; 1678 #if defined(PETSC_HAVE_CUDA) 1679 A->offloadmask = PETSC_OFFLOAD_CPU; 1680 #endif 1681 PetscFunctionReturn(0); 1682 } 1683 1684 static PetscErrorCode MatDenseResetArray_SeqDense(Mat A) 1685 { 1686 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1687 1688 PetscFunctionBegin; 1689 PetscCheck(!a->vecinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 1690 PetscCheck(!a->matinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 1691 a->v = a->unplacedarray; 1692 a->user_alloc = a->unplaced_user_alloc; 1693 a->unplacedarray = NULL; 1694 #if defined(PETSC_HAVE_CUDA) 1695 A->offloadmask = PETSC_OFFLOAD_CPU; 1696 #endif 1697 PetscFunctionReturn(0); 1698 } 1699 1700 static PetscErrorCode MatDenseReplaceArray_SeqDense(Mat A,const PetscScalar *array) 1701 { 1702 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1703 1704 PetscFunctionBegin; 1705 PetscCheck(!a->vecinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 1706 PetscCheck(!a->matinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 1707 if (!a->user_alloc) PetscCall(PetscFree(a->v)); 1708 a->v = (PetscScalar*) array; 1709 a->user_alloc = PETSC_FALSE; 1710 #if defined(PETSC_HAVE_CUDA) 1711 A->offloadmask = PETSC_OFFLOAD_CPU; 1712 #endif 1713 PetscFunctionReturn(0); 1714 } 1715 1716 PetscErrorCode MatDestroy_SeqDense(Mat mat) 1717 { 1718 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 1719 1720 PetscFunctionBegin; 1721 #if defined(PETSC_USE_LOG) 1722 PetscLogObjectState((PetscObject)mat,"Rows %" PetscInt_FMT " Cols %" PetscInt_FMT,mat->rmap->n,mat->cmap->n); 1723 #endif 1724 PetscCall(VecDestroy(&(l->qrrhs))); 1725 PetscCall(PetscFree(l->tau)); 1726 PetscCall(PetscFree(l->pivots)); 1727 PetscCall(PetscFree(l->fwork)); 1728 PetscCall(MatDestroy(&l->ptapwork)); 1729 if (!l->user_alloc) PetscCall(PetscFree(l->v)); 1730 if (!l->unplaced_user_alloc) PetscCall(PetscFree(l->unplacedarray)); 1731 PetscCheck(!l->vecinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 1732 PetscCheck(!l->matinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 1733 PetscCall(VecDestroy(&l->cvec)); 1734 PetscCall(MatDestroy(&l->cmat)); 1735 PetscCall(PetscFree(mat->data)); 1736 1737 PetscCall(PetscObjectChangeTypeName((PetscObject)mat,NULL)); 1738 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatQRFactor_C",NULL)); 1739 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatQRFactorSymbolic_C",NULL)); 1740 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatQRFactorNumeric_C",NULL)); 1741 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetLDA_C",NULL)); 1742 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseSetLDA_C",NULL)); 1743 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL)); 1744 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL)); 1745 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL)); 1746 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL)); 1747 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseReplaceArray_C",NULL)); 1748 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayRead_C",NULL)); 1749 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayRead_C",NULL)); 1750 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayWrite_C",NULL)); 1751 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayWrite_C",NULL)); 1752 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL)); 1753 #if defined(PETSC_HAVE_ELEMENTAL) 1754 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL)); 1755 #endif 1756 #if defined(PETSC_HAVE_SCALAPACK) 1757 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_scalapack_C",NULL)); 1758 #endif 1759 #if defined(PETSC_HAVE_CUDA) 1760 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqdensecuda_C",NULL)); 1761 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",NULL)); 1762 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdense_C",NULL)); 1763 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdense_seqdensecuda_C",NULL)); 1764 #endif 1765 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL)); 1766 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqaij_seqdense_C",NULL)); 1767 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdense_seqdense_C",NULL)); 1768 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqbaij_seqdense_C",NULL)); 1769 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqsbaij_seqdense_C",NULL)); 1770 1771 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumn_C",NULL)); 1772 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumn_C",NULL)); 1773 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVec_C",NULL)); 1774 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVec_C",NULL)); 1775 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVecRead_C",NULL)); 1776 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVecRead_C",NULL)); 1777 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVecWrite_C",NULL)); 1778 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVecWrite_C",NULL)); 1779 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetSubMatrix_C",NULL)); 1780 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreSubMatrix_C",NULL)); 1781 PetscFunctionReturn(0); 1782 } 1783 1784 static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout) 1785 { 1786 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1787 PetscInt k,j,m = A->rmap->n, M = mat->lda, n = A->cmap->n; 1788 PetscScalar *v,tmp; 1789 1790 PetscFunctionBegin; 1791 if (reuse == MAT_INPLACE_MATRIX) { 1792 if (m == n) { /* in place transpose */ 1793 PetscCall(MatDenseGetArray(A,&v)); 1794 for (j=0; j<m; j++) { 1795 for (k=0; k<j; k++) { 1796 tmp = v[j + k*M]; 1797 v[j + k*M] = v[k + j*M]; 1798 v[k + j*M] = tmp; 1799 } 1800 } 1801 PetscCall(MatDenseRestoreArray(A,&v)); 1802 } else { /* reuse memory, temporary allocates new memory */ 1803 PetscScalar *v2; 1804 PetscLayout tmplayout; 1805 1806 PetscCall(PetscMalloc1((size_t)m*n,&v2)); 1807 PetscCall(MatDenseGetArray(A,&v)); 1808 for (j=0; j<n; j++) { 1809 for (k=0; k<m; k++) v2[j + (size_t)k*n] = v[k + (size_t)j*M]; 1810 } 1811 PetscCall(PetscArraycpy(v,v2,(size_t)m*n)); 1812 PetscCall(PetscFree(v2)); 1813 PetscCall(MatDenseRestoreArray(A,&v)); 1814 /* cleanup size dependent quantities */ 1815 PetscCall(VecDestroy(&mat->cvec)); 1816 PetscCall(MatDestroy(&mat->cmat)); 1817 PetscCall(PetscFree(mat->pivots)); 1818 PetscCall(PetscFree(mat->fwork)); 1819 PetscCall(MatDestroy(&mat->ptapwork)); 1820 /* swap row/col layouts */ 1821 mat->lda = n; 1822 tmplayout = A->rmap; 1823 A->rmap = A->cmap; 1824 A->cmap = tmplayout; 1825 } 1826 } else { /* out-of-place transpose */ 1827 Mat tmat; 1828 Mat_SeqDense *tmatd; 1829 PetscScalar *v2; 1830 PetscInt M2; 1831 1832 if (reuse == MAT_INITIAL_MATRIX) { 1833 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&tmat)); 1834 PetscCall(MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n)); 1835 PetscCall(MatSetType(tmat,((PetscObject)A)->type_name)); 1836 PetscCall(MatSeqDenseSetPreallocation(tmat,NULL)); 1837 } else tmat = *matout; 1838 1839 PetscCall(MatDenseGetArrayRead(A,(const PetscScalar**)&v)); 1840 PetscCall(MatDenseGetArray(tmat,&v2)); 1841 tmatd = (Mat_SeqDense*)tmat->data; 1842 M2 = tmatd->lda; 1843 for (j=0; j<n; j++) { 1844 for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M]; 1845 } 1846 PetscCall(MatDenseRestoreArray(tmat,&v2)); 1847 PetscCall(MatDenseRestoreArrayRead(A,(const PetscScalar**)&v)); 1848 PetscCall(MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY)); 1849 PetscCall(MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY)); 1850 *matout = tmat; 1851 } 1852 PetscFunctionReturn(0); 1853 } 1854 1855 static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg) 1856 { 1857 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1858 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1859 PetscInt i; 1860 const PetscScalar *v1,*v2; 1861 1862 PetscFunctionBegin; 1863 if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1864 if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1865 PetscCall(MatDenseGetArrayRead(A1,&v1)); 1866 PetscCall(MatDenseGetArrayRead(A2,&v2)); 1867 for (i=0; i<A1->cmap->n; i++) { 1868 PetscCall(PetscArraycmp(v1,v2,A1->rmap->n,flg)); 1869 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 1870 v1 += mat1->lda; 1871 v2 += mat2->lda; 1872 } 1873 PetscCall(MatDenseRestoreArrayRead(A1,&v1)); 1874 PetscCall(MatDenseRestoreArrayRead(A2,&v2)); 1875 *flg = PETSC_TRUE; 1876 PetscFunctionReturn(0); 1877 } 1878 1879 static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1880 { 1881 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1882 PetscInt i,n,len; 1883 PetscScalar *x; 1884 const PetscScalar *vv; 1885 1886 PetscFunctionBegin; 1887 PetscCall(VecGetSize(v,&n)); 1888 PetscCall(VecGetArray(v,&x)); 1889 len = PetscMin(A->rmap->n,A->cmap->n); 1890 PetscCall(MatDenseGetArrayRead(A,&vv)); 1891 PetscCheck(n == A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1892 for (i=0; i<len; i++) { 1893 x[i] = vv[i*mat->lda + i]; 1894 } 1895 PetscCall(MatDenseRestoreArrayRead(A,&vv)); 1896 PetscCall(VecRestoreArray(v,&x)); 1897 PetscFunctionReturn(0); 1898 } 1899 1900 static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1901 { 1902 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1903 const PetscScalar *l,*r; 1904 PetscScalar x,*v,*vv; 1905 PetscInt i,j,m = A->rmap->n,n = A->cmap->n; 1906 1907 PetscFunctionBegin; 1908 PetscCall(MatDenseGetArray(A,&vv)); 1909 if (ll) { 1910 PetscCall(VecGetSize(ll,&m)); 1911 PetscCall(VecGetArrayRead(ll,&l)); 1912 PetscCheck(m == A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1913 for (i=0; i<m; i++) { 1914 x = l[i]; 1915 v = vv + i; 1916 for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;} 1917 } 1918 PetscCall(VecRestoreArrayRead(ll,&l)); 1919 PetscCall(PetscLogFlops(1.0*n*m)); 1920 } 1921 if (rr) { 1922 PetscCall(VecGetSize(rr,&n)); 1923 PetscCall(VecGetArrayRead(rr,&r)); 1924 PetscCheck(n == A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1925 for (i=0; i<n; i++) { 1926 x = r[i]; 1927 v = vv + i*mat->lda; 1928 for (j=0; j<m; j++) (*v++) *= x; 1929 } 1930 PetscCall(VecRestoreArrayRead(rr,&r)); 1931 PetscCall(PetscLogFlops(1.0*n*m)); 1932 } 1933 PetscCall(MatDenseRestoreArray(A,&vv)); 1934 PetscFunctionReturn(0); 1935 } 1936 1937 PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1938 { 1939 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1940 PetscScalar *v,*vv; 1941 PetscReal sum = 0.0; 1942 PetscInt lda, m=A->rmap->n,i,j; 1943 1944 PetscFunctionBegin; 1945 PetscCall(MatDenseGetArrayRead(A,(const PetscScalar**)&vv)); 1946 PetscCall(MatDenseGetLDA(A,&lda)); 1947 v = vv; 1948 if (type == NORM_FROBENIUS) { 1949 if (lda>m) { 1950 for (j=0; j<A->cmap->n; j++) { 1951 v = vv+j*lda; 1952 for (i=0; i<m; i++) { 1953 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1954 } 1955 } 1956 } else { 1957 #if defined(PETSC_USE_REAL___FP16) 1958 PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n; 1959 PetscStackCallBLAS("BLASnrm2",*nrm = BLASnrm2_(&cnt,v,&one)); 1960 } 1961 #else 1962 for (i=0; i<A->cmap->n*A->rmap->n; i++) { 1963 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1964 } 1965 } 1966 *nrm = PetscSqrtReal(sum); 1967 #endif 1968 PetscCall(PetscLogFlops(2.0*A->cmap->n*A->rmap->n)); 1969 } else if (type == NORM_1) { 1970 *nrm = 0.0; 1971 for (j=0; j<A->cmap->n; j++) { 1972 v = vv + j*mat->lda; 1973 sum = 0.0; 1974 for (i=0; i<A->rmap->n; i++) { 1975 sum += PetscAbsScalar(*v); v++; 1976 } 1977 if (sum > *nrm) *nrm = sum; 1978 } 1979 PetscCall(PetscLogFlops(1.0*A->cmap->n*A->rmap->n)); 1980 } else if (type == NORM_INFINITY) { 1981 *nrm = 0.0; 1982 for (j=0; j<A->rmap->n; j++) { 1983 v = vv + j; 1984 sum = 0.0; 1985 for (i=0; i<A->cmap->n; i++) { 1986 sum += PetscAbsScalar(*v); v += mat->lda; 1987 } 1988 if (sum > *nrm) *nrm = sum; 1989 } 1990 PetscCall(PetscLogFlops(1.0*A->cmap->n*A->rmap->n)); 1991 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 1992 PetscCall(MatDenseRestoreArrayRead(A,(const PetscScalar**)&vv)); 1993 PetscFunctionReturn(0); 1994 } 1995 1996 static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg) 1997 { 1998 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1999 2000 PetscFunctionBegin; 2001 switch (op) { 2002 case MAT_ROW_ORIENTED: 2003 aij->roworiented = flg; 2004 break; 2005 case MAT_NEW_NONZERO_LOCATIONS: 2006 case MAT_NEW_NONZERO_LOCATION_ERR: 2007 case MAT_NEW_NONZERO_ALLOCATION_ERR: 2008 case MAT_FORCE_DIAGONAL_ENTRIES: 2009 case MAT_KEEP_NONZERO_PATTERN: 2010 case MAT_IGNORE_OFF_PROC_ENTRIES: 2011 case MAT_USE_HASH_TABLE: 2012 case MAT_IGNORE_ZERO_ENTRIES: 2013 case MAT_IGNORE_LOWER_TRIANGULAR: 2014 case MAT_SORTED_FULL: 2015 PetscCall(PetscInfo(A,"Option %s ignored\n",MatOptions[op])); 2016 break; 2017 case MAT_SPD: 2018 case MAT_SYMMETRIC: 2019 case MAT_STRUCTURALLY_SYMMETRIC: 2020 case MAT_HERMITIAN: 2021 case MAT_SYMMETRY_ETERNAL: 2022 /* These options are handled directly by MatSetOption() */ 2023 break; 2024 default: 2025 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 2026 } 2027 PetscFunctionReturn(0); 2028 } 2029 2030 PetscErrorCode MatZeroEntries_SeqDense(Mat A) 2031 { 2032 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 2033 PetscInt lda=l->lda,m=A->rmap->n,n=A->cmap->n,j; 2034 PetscScalar *v; 2035 2036 PetscFunctionBegin; 2037 PetscCall(MatDenseGetArrayWrite(A,&v)); 2038 if (lda>m) { 2039 for (j=0; j<n; j++) { 2040 PetscCall(PetscArrayzero(v+j*lda,m)); 2041 } 2042 } else { 2043 PetscCall(PetscArrayzero(v,PetscInt64Mult(m,n))); 2044 } 2045 PetscCall(MatDenseRestoreArrayWrite(A,&v)); 2046 PetscFunctionReturn(0); 2047 } 2048 2049 static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2050 { 2051 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 2052 PetscInt m = l->lda, n = A->cmap->n, i,j; 2053 PetscScalar *slot,*bb,*v; 2054 const PetscScalar *xx; 2055 2056 PetscFunctionBegin; 2057 if (PetscDefined(USE_DEBUG)) { 2058 for (i=0; i<N; i++) { 2059 PetscCheck(rows[i] >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 2060 PetscCheck(rows[i] < A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %" PetscInt_FMT " requested to be zeroed greater than or equal number of rows %" PetscInt_FMT,rows[i],A->rmap->n); 2061 } 2062 } 2063 if (!N) PetscFunctionReturn(0); 2064 2065 /* fix right hand side if needed */ 2066 if (x && b) { 2067 PetscCall(VecGetArrayRead(x,&xx)); 2068 PetscCall(VecGetArray(b,&bb)); 2069 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]]; 2070 PetscCall(VecRestoreArrayRead(x,&xx)); 2071 PetscCall(VecRestoreArray(b,&bb)); 2072 } 2073 2074 PetscCall(MatDenseGetArray(A,&v)); 2075 for (i=0; i<N; i++) { 2076 slot = v + rows[i]; 2077 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 2078 } 2079 if (diag != 0.0) { 2080 PetscCheck(A->rmap->n == A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 2081 for (i=0; i<N; i++) { 2082 slot = v + (m+1)*rows[i]; 2083 *slot = diag; 2084 } 2085 } 2086 PetscCall(MatDenseRestoreArray(A,&v)); 2087 PetscFunctionReturn(0); 2088 } 2089 2090 static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A,PetscInt *lda) 2091 { 2092 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 2093 2094 PetscFunctionBegin; 2095 *lda = mat->lda; 2096 PetscFunctionReturn(0); 2097 } 2098 2099 PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar **array) 2100 { 2101 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 2102 2103 PetscFunctionBegin; 2104 PetscCheck(!mat->matinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 2105 *array = mat->v; 2106 PetscFunctionReturn(0); 2107 } 2108 2109 PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar **array) 2110 { 2111 PetscFunctionBegin; 2112 if (array) *array = NULL; 2113 PetscFunctionReturn(0); 2114 } 2115 2116 /*@ 2117 MatDenseGetLDA - gets the leading dimension of the array returned from MatDenseGetArray() 2118 2119 Not collective 2120 2121 Input Parameter: 2122 . mat - a MATSEQDENSE or MATMPIDENSE matrix 2123 2124 Output Parameter: 2125 . lda - the leading dimension 2126 2127 Level: intermediate 2128 2129 .seealso: `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseSetLDA()` 2130 @*/ 2131 PetscErrorCode MatDenseGetLDA(Mat A,PetscInt *lda) 2132 { 2133 PetscFunctionBegin; 2134 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 2135 PetscValidIntPointer(lda,2); 2136 MatCheckPreallocated(A,1); 2137 PetscUseMethod(A,"MatDenseGetLDA_C",(Mat,PetscInt*),(A,lda)); 2138 PetscFunctionReturn(0); 2139 } 2140 2141 /*@ 2142 MatDenseSetLDA - Sets the leading dimension of the array used by the dense matrix 2143 2144 Not collective 2145 2146 Input Parameters: 2147 + mat - a MATSEQDENSE or MATMPIDENSE matrix 2148 - lda - the leading dimension 2149 2150 Level: intermediate 2151 2152 .seealso: `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetLDA()` 2153 @*/ 2154 PetscErrorCode MatDenseSetLDA(Mat A,PetscInt lda) 2155 { 2156 PetscFunctionBegin; 2157 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 2158 PetscTryMethod(A,"MatDenseSetLDA_C",(Mat,PetscInt),(A,lda)); 2159 PetscFunctionReturn(0); 2160 } 2161 2162 /*@C 2163 MatDenseGetArray - gives read-write access to the array where the data for a dense matrix is stored 2164 2165 Logically Collective on Mat 2166 2167 Input Parameter: 2168 . mat - a dense matrix 2169 2170 Output Parameter: 2171 . array - pointer to the data 2172 2173 Level: intermediate 2174 2175 .seealso: `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()` 2176 @*/ 2177 PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array) 2178 { 2179 PetscFunctionBegin; 2180 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 2181 PetscValidPointer(array,2); 2182 PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array)); 2183 PetscFunctionReturn(0); 2184 } 2185 2186 /*@C 2187 MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray() 2188 2189 Logically Collective on Mat 2190 2191 Input Parameters: 2192 + mat - a dense matrix 2193 - array - pointer to the data 2194 2195 Level: intermediate 2196 2197 .seealso: `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()` 2198 @*/ 2199 PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array) 2200 { 2201 PetscFunctionBegin; 2202 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 2203 PetscValidPointer(array,2); 2204 PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array)); 2205 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 2206 #if defined(PETSC_HAVE_CUDA) 2207 A->offloadmask = PETSC_OFFLOAD_CPU; 2208 #endif 2209 PetscFunctionReturn(0); 2210 } 2211 2212 /*@C 2213 MatDenseGetArrayRead - gives read-only access to the array where the data for a dense matrix is stored 2214 2215 Not Collective 2216 2217 Input Parameter: 2218 . mat - a dense matrix 2219 2220 Output Parameter: 2221 . array - pointer to the data 2222 2223 Level: intermediate 2224 2225 .seealso: `MatDenseRestoreArrayRead()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()` 2226 @*/ 2227 PetscErrorCode MatDenseGetArrayRead(Mat A,const PetscScalar **array) 2228 { 2229 PetscFunctionBegin; 2230 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 2231 PetscValidPointer(array,2); 2232 PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array)); 2233 PetscFunctionReturn(0); 2234 } 2235 2236 /*@C 2237 MatDenseRestoreArrayRead - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArrayRead() 2238 2239 Not Collective 2240 2241 Input Parameters: 2242 + mat - a dense matrix 2243 - array - pointer to the data 2244 2245 Level: intermediate 2246 2247 .seealso: `MatDenseGetArrayRead()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()` 2248 @*/ 2249 PetscErrorCode MatDenseRestoreArrayRead(Mat A,const PetscScalar **array) 2250 { 2251 PetscFunctionBegin; 2252 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 2253 PetscValidPointer(array,2); 2254 PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array)); 2255 PetscFunctionReturn(0); 2256 } 2257 2258 /*@C 2259 MatDenseGetArrayWrite - gives write-only access to the array where the data for a dense matrix is stored 2260 2261 Not Collective 2262 2263 Input Parameter: 2264 . mat - a dense matrix 2265 2266 Output Parameter: 2267 . array - pointer to the data 2268 2269 Level: intermediate 2270 2271 .seealso: `MatDenseRestoreArrayWrite()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()` 2272 @*/ 2273 PetscErrorCode MatDenseGetArrayWrite(Mat A,PetscScalar **array) 2274 { 2275 PetscFunctionBegin; 2276 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 2277 PetscValidPointer(array,2); 2278 PetscUseMethod(A,"MatDenseGetArrayWrite_C",(Mat,PetscScalar**),(A,array)); 2279 PetscFunctionReturn(0); 2280 } 2281 2282 /*@C 2283 MatDenseRestoreArrayWrite - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArrayWrite() 2284 2285 Not Collective 2286 2287 Input Parameters: 2288 + mat - a dense matrix 2289 - array - pointer to the data 2290 2291 Level: intermediate 2292 2293 .seealso: `MatDenseGetArrayWrite()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()` 2294 @*/ 2295 PetscErrorCode MatDenseRestoreArrayWrite(Mat A,PetscScalar **array) 2296 { 2297 PetscFunctionBegin; 2298 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 2299 PetscValidPointer(array,2); 2300 PetscUseMethod(A,"MatDenseRestoreArrayWrite_C",(Mat,PetscScalar**),(A,array)); 2301 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 2302 #if defined(PETSC_HAVE_CUDA) 2303 A->offloadmask = PETSC_OFFLOAD_CPU; 2304 #endif 2305 PetscFunctionReturn(0); 2306 } 2307 2308 static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,MatReuse scall,Mat *B) 2309 { 2310 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 2311 PetscInt i,j,nrows,ncols,ldb; 2312 const PetscInt *irow,*icol; 2313 PetscScalar *av,*bv,*v = mat->v; 2314 Mat newmat; 2315 2316 PetscFunctionBegin; 2317 PetscCall(ISGetIndices(isrow,&irow)); 2318 PetscCall(ISGetIndices(iscol,&icol)); 2319 PetscCall(ISGetLocalSize(isrow,&nrows)); 2320 PetscCall(ISGetLocalSize(iscol,&ncols)); 2321 2322 /* Check submatrixcall */ 2323 if (scall == MAT_REUSE_MATRIX) { 2324 PetscInt n_cols,n_rows; 2325 PetscCall(MatGetSize(*B,&n_rows,&n_cols)); 2326 if (n_rows != nrows || n_cols != ncols) { 2327 /* resize the result matrix to match number of requested rows/columns */ 2328 PetscCall(MatSetSizes(*B,nrows,ncols,nrows,ncols)); 2329 } 2330 newmat = *B; 2331 } else { 2332 /* Create and fill new matrix */ 2333 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&newmat)); 2334 PetscCall(MatSetSizes(newmat,nrows,ncols,nrows,ncols)); 2335 PetscCall(MatSetType(newmat,((PetscObject)A)->type_name)); 2336 PetscCall(MatSeqDenseSetPreallocation(newmat,NULL)); 2337 } 2338 2339 /* Now extract the data pointers and do the copy,column at a time */ 2340 PetscCall(MatDenseGetArray(newmat,&bv)); 2341 PetscCall(MatDenseGetLDA(newmat,&ldb)); 2342 for (i=0; i<ncols; i++) { 2343 av = v + mat->lda*icol[i]; 2344 for (j=0; j<nrows; j++) bv[j] = av[irow[j]]; 2345 bv += ldb; 2346 } 2347 PetscCall(MatDenseRestoreArray(newmat,&bv)); 2348 2349 /* Assemble the matrices so that the correct flags are set */ 2350 PetscCall(MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY)); 2351 PetscCall(MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY)); 2352 2353 /* Free work space */ 2354 PetscCall(ISRestoreIndices(isrow,&irow)); 2355 PetscCall(ISRestoreIndices(iscol,&icol)); 2356 *B = newmat; 2357 PetscFunctionReturn(0); 2358 } 2359 2360 static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2361 { 2362 PetscInt i; 2363 2364 PetscFunctionBegin; 2365 if (scall == MAT_INITIAL_MATRIX) { 2366 PetscCall(PetscCalloc1(n,B)); 2367 } 2368 2369 for (i=0; i<n; i++) { 2370 PetscCall(MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],scall,&(*B)[i])); 2371 } 2372 PetscFunctionReturn(0); 2373 } 2374 2375 static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 2376 { 2377 PetscFunctionBegin; 2378 PetscFunctionReturn(0); 2379 } 2380 2381 static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 2382 { 2383 PetscFunctionBegin; 2384 PetscFunctionReturn(0); 2385 } 2386 2387 PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 2388 { 2389 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data; 2390 const PetscScalar *va; 2391 PetscScalar *vb; 2392 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 2393 2394 PetscFunctionBegin; 2395 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 2396 if (A->ops->copy != B->ops->copy) { 2397 PetscCall(MatCopy_Basic(A,B,str)); 2398 PetscFunctionReturn(0); 2399 } 2400 PetscCheck(m == B->rmap->n && n == B->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 2401 PetscCall(MatDenseGetArrayRead(A,&va)); 2402 PetscCall(MatDenseGetArray(B,&vb)); 2403 if (lda1>m || lda2>m) { 2404 for (j=0; j<n; j++) { 2405 PetscCall(PetscArraycpy(vb+j*lda2,va+j*lda1,m)); 2406 } 2407 } else { 2408 PetscCall(PetscArraycpy(vb,va,A->rmap->n*A->cmap->n)); 2409 } 2410 PetscCall(MatDenseRestoreArray(B,&vb)); 2411 PetscCall(MatDenseRestoreArrayRead(A,&va)); 2412 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 2413 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 2414 PetscFunctionReturn(0); 2415 } 2416 2417 PetscErrorCode MatSetUp_SeqDense(Mat A) 2418 { 2419 PetscFunctionBegin; 2420 PetscCall(PetscLayoutSetUp(A->rmap)); 2421 PetscCall(PetscLayoutSetUp(A->cmap)); 2422 if (!A->preallocated) { 2423 PetscCall(MatSeqDenseSetPreallocation(A,NULL)); 2424 } 2425 PetscFunctionReturn(0); 2426 } 2427 2428 static PetscErrorCode MatConjugate_SeqDense(Mat A) 2429 { 2430 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 2431 PetscInt i,j; 2432 PetscInt min = PetscMin(A->rmap->n,A->cmap->n); 2433 PetscScalar *aa; 2434 2435 PetscFunctionBegin; 2436 PetscCall(MatDenseGetArray(A,&aa)); 2437 for (j=0; j<A->cmap->n; j++) { 2438 for (i=0; i<A->rmap->n; i++) aa[i+j*mat->lda] = PetscConj(aa[i+j*mat->lda]); 2439 } 2440 PetscCall(MatDenseRestoreArray(A,&aa)); 2441 if (mat->tau) for (i = 0; i < min; i++) mat->tau[i] = PetscConj(mat->tau[i]); 2442 PetscFunctionReturn(0); 2443 } 2444 2445 static PetscErrorCode MatRealPart_SeqDense(Mat A) 2446 { 2447 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 2448 PetscInt i,j; 2449 PetscScalar *aa; 2450 2451 PetscFunctionBegin; 2452 PetscCall(MatDenseGetArray(A,&aa)); 2453 for (j=0; j<A->cmap->n; j++) { 2454 for (i=0; i<A->rmap->n; i++) aa[i+j*mat->lda] = PetscRealPart(aa[i+j*mat->lda]); 2455 } 2456 PetscCall(MatDenseRestoreArray(A,&aa)); 2457 PetscFunctionReturn(0); 2458 } 2459 2460 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 2461 { 2462 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 2463 PetscInt i,j; 2464 PetscScalar *aa; 2465 2466 PetscFunctionBegin; 2467 PetscCall(MatDenseGetArray(A,&aa)); 2468 for (j=0; j<A->cmap->n; j++) { 2469 for (i=0; i<A->rmap->n; i++) aa[i+j*mat->lda] = PetscImaginaryPart(aa[i+j*mat->lda]); 2470 } 2471 PetscCall(MatDenseRestoreArray(A,&aa)); 2472 PetscFunctionReturn(0); 2473 } 2474 2475 /* ----------------------------------------------------------------*/ 2476 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C) 2477 { 2478 PetscInt m=A->rmap->n,n=B->cmap->n; 2479 PetscBool cisdense; 2480 2481 PetscFunctionBegin; 2482 PetscCall(MatSetSizes(C,m,n,m,n)); 2483 PetscCall(PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"")); 2484 if (!cisdense) { 2485 PetscBool flg; 2486 2487 PetscCall(PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg)); 2488 PetscCall(MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE)); 2489 } 2490 PetscCall(MatSetUp(C)); 2491 PetscFunctionReturn(0); 2492 } 2493 2494 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2495 { 2496 Mat_SeqDense *a=(Mat_SeqDense*)A->data,*b=(Mat_SeqDense*)B->data,*c=(Mat_SeqDense*)C->data; 2497 PetscBLASInt m,n,k; 2498 const PetscScalar *av,*bv; 2499 PetscScalar *cv; 2500 PetscScalar _DOne=1.0,_DZero=0.0; 2501 2502 PetscFunctionBegin; 2503 PetscCall(PetscBLASIntCast(C->rmap->n,&m)); 2504 PetscCall(PetscBLASIntCast(C->cmap->n,&n)); 2505 PetscCall(PetscBLASIntCast(A->cmap->n,&k)); 2506 if (!m || !n || !k) PetscFunctionReturn(0); 2507 PetscCall(MatDenseGetArrayRead(A,&av)); 2508 PetscCall(MatDenseGetArrayRead(B,&bv)); 2509 PetscCall(MatDenseGetArrayWrite(C,&cv)); 2510 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda)); 2511 PetscCall(PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1))); 2512 PetscCall(MatDenseRestoreArrayRead(A,&av)); 2513 PetscCall(MatDenseRestoreArrayRead(B,&bv)); 2514 PetscCall(MatDenseRestoreArrayWrite(C,&cv)); 2515 PetscFunctionReturn(0); 2516 } 2517 2518 PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C) 2519 { 2520 PetscInt m=A->rmap->n,n=B->rmap->n; 2521 PetscBool cisdense; 2522 2523 PetscFunctionBegin; 2524 PetscCall(MatSetSizes(C,m,n,m,n)); 2525 PetscCall(PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"")); 2526 if (!cisdense) { 2527 PetscBool flg; 2528 2529 PetscCall(PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg)); 2530 PetscCall(MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE)); 2531 } 2532 PetscCall(MatSetUp(C)); 2533 PetscFunctionReturn(0); 2534 } 2535 2536 PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2537 { 2538 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2539 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2540 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2541 const PetscScalar *av,*bv; 2542 PetscScalar *cv; 2543 PetscBLASInt m,n,k; 2544 PetscScalar _DOne=1.0,_DZero=0.0; 2545 2546 PetscFunctionBegin; 2547 PetscCall(PetscBLASIntCast(C->rmap->n,&m)); 2548 PetscCall(PetscBLASIntCast(C->cmap->n,&n)); 2549 PetscCall(PetscBLASIntCast(A->cmap->n,&k)); 2550 if (!m || !n || !k) PetscFunctionReturn(0); 2551 PetscCall(MatDenseGetArrayRead(A,&av)); 2552 PetscCall(MatDenseGetArrayRead(B,&bv)); 2553 PetscCall(MatDenseGetArrayWrite(C,&cv)); 2554 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda)); 2555 PetscCall(MatDenseRestoreArrayRead(A,&av)); 2556 PetscCall(MatDenseRestoreArrayRead(B,&bv)); 2557 PetscCall(MatDenseRestoreArrayWrite(C,&cv)); 2558 PetscCall(PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1))); 2559 PetscFunctionReturn(0); 2560 } 2561 2562 PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C) 2563 { 2564 PetscInt m=A->cmap->n,n=B->cmap->n; 2565 PetscBool cisdense; 2566 2567 PetscFunctionBegin; 2568 PetscCall(MatSetSizes(C,m,n,m,n)); 2569 PetscCall(PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"")); 2570 if (!cisdense) { 2571 PetscBool flg; 2572 2573 PetscCall(PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg)); 2574 PetscCall(MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE)); 2575 } 2576 PetscCall(MatSetUp(C)); 2577 PetscFunctionReturn(0); 2578 } 2579 2580 PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2581 { 2582 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2583 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2584 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2585 const PetscScalar *av,*bv; 2586 PetscScalar *cv; 2587 PetscBLASInt m,n,k; 2588 PetscScalar _DOne=1.0,_DZero=0.0; 2589 2590 PetscFunctionBegin; 2591 PetscCall(PetscBLASIntCast(C->rmap->n,&m)); 2592 PetscCall(PetscBLASIntCast(C->cmap->n,&n)); 2593 PetscCall(PetscBLASIntCast(A->rmap->n,&k)); 2594 if (!m || !n || !k) PetscFunctionReturn(0); 2595 PetscCall(MatDenseGetArrayRead(A,&av)); 2596 PetscCall(MatDenseGetArrayRead(B,&bv)); 2597 PetscCall(MatDenseGetArrayWrite(C,&cv)); 2598 PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda)); 2599 PetscCall(MatDenseRestoreArrayRead(A,&av)); 2600 PetscCall(MatDenseRestoreArrayRead(B,&bv)); 2601 PetscCall(MatDenseRestoreArrayWrite(C,&cv)); 2602 PetscCall(PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1))); 2603 PetscFunctionReturn(0); 2604 } 2605 2606 /* ----------------------------------------------- */ 2607 static PetscErrorCode MatProductSetFromOptions_SeqDense_AB(Mat C) 2608 { 2609 PetscFunctionBegin; 2610 C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqDense; 2611 C->ops->productsymbolic = MatProductSymbolic_AB; 2612 PetscFunctionReturn(0); 2613 } 2614 2615 static PetscErrorCode MatProductSetFromOptions_SeqDense_AtB(Mat C) 2616 { 2617 PetscFunctionBegin; 2618 C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_SeqDense_SeqDense; 2619 C->ops->productsymbolic = MatProductSymbolic_AtB; 2620 PetscFunctionReturn(0); 2621 } 2622 2623 static PetscErrorCode MatProductSetFromOptions_SeqDense_ABt(Mat C) 2624 { 2625 PetscFunctionBegin; 2626 C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqDense_SeqDense; 2627 C->ops->productsymbolic = MatProductSymbolic_ABt; 2628 PetscFunctionReturn(0); 2629 } 2630 2631 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense(Mat C) 2632 { 2633 Mat_Product *product = C->product; 2634 2635 PetscFunctionBegin; 2636 switch (product->type) { 2637 case MATPRODUCT_AB: 2638 PetscCall(MatProductSetFromOptions_SeqDense_AB(C)); 2639 break; 2640 case MATPRODUCT_AtB: 2641 PetscCall(MatProductSetFromOptions_SeqDense_AtB(C)); 2642 break; 2643 case MATPRODUCT_ABt: 2644 PetscCall(MatProductSetFromOptions_SeqDense_ABt(C)); 2645 break; 2646 default: 2647 break; 2648 } 2649 PetscFunctionReturn(0); 2650 } 2651 /* ----------------------------------------------- */ 2652 2653 static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 2654 { 2655 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2656 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2657 PetscScalar *x; 2658 const PetscScalar *aa; 2659 2660 PetscFunctionBegin; 2661 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2662 PetscCall(VecGetArray(v,&x)); 2663 PetscCall(VecGetLocalSize(v,&p)); 2664 PetscCall(MatDenseGetArrayRead(A,&aa)); 2665 PetscCheck(p == A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2666 for (i=0; i<m; i++) { 2667 x[i] = aa[i]; if (idx) idx[i] = 0; 2668 for (j=1; j<n; j++) { 2669 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;} 2670 } 2671 } 2672 PetscCall(MatDenseRestoreArrayRead(A,&aa)); 2673 PetscCall(VecRestoreArray(v,&x)); 2674 PetscFunctionReturn(0); 2675 } 2676 2677 static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 2678 { 2679 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2680 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2681 PetscScalar *x; 2682 PetscReal atmp; 2683 const PetscScalar *aa; 2684 2685 PetscFunctionBegin; 2686 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2687 PetscCall(VecGetArray(v,&x)); 2688 PetscCall(VecGetLocalSize(v,&p)); 2689 PetscCall(MatDenseGetArrayRead(A,&aa)); 2690 PetscCheck(p == A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2691 for (i=0; i<m; i++) { 2692 x[i] = PetscAbsScalar(aa[i]); 2693 for (j=1; j<n; j++) { 2694 atmp = PetscAbsScalar(aa[i+a->lda*j]); 2695 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 2696 } 2697 } 2698 PetscCall(MatDenseRestoreArrayRead(A,&aa)); 2699 PetscCall(VecRestoreArray(v,&x)); 2700 PetscFunctionReturn(0); 2701 } 2702 2703 static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 2704 { 2705 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2706 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2707 PetscScalar *x; 2708 const PetscScalar *aa; 2709 2710 PetscFunctionBegin; 2711 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2712 PetscCall(MatDenseGetArrayRead(A,&aa)); 2713 PetscCall(VecGetArray(v,&x)); 2714 PetscCall(VecGetLocalSize(v,&p)); 2715 PetscCheck(p == A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2716 for (i=0; i<m; i++) { 2717 x[i] = aa[i]; if (idx) idx[i] = 0; 2718 for (j=1; j<n; j++) { 2719 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;} 2720 } 2721 } 2722 PetscCall(VecRestoreArray(v,&x)); 2723 PetscCall(MatDenseRestoreArrayRead(A,&aa)); 2724 PetscFunctionReturn(0); 2725 } 2726 2727 PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 2728 { 2729 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2730 PetscScalar *x; 2731 const PetscScalar *aa; 2732 2733 PetscFunctionBegin; 2734 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2735 PetscCall(MatDenseGetArrayRead(A,&aa)); 2736 PetscCall(VecGetArray(v,&x)); 2737 PetscCall(PetscArraycpy(x,aa+col*a->lda,A->rmap->n)); 2738 PetscCall(VecRestoreArray(v,&x)); 2739 PetscCall(MatDenseRestoreArrayRead(A,&aa)); 2740 PetscFunctionReturn(0); 2741 } 2742 2743 PETSC_INTERN PetscErrorCode MatGetColumnReductions_SeqDense(Mat A,PetscInt type,PetscReal *reductions) 2744 { 2745 PetscInt i,j,m,n; 2746 const PetscScalar *a; 2747 2748 PetscFunctionBegin; 2749 PetscCall(MatGetSize(A,&m,&n)); 2750 PetscCall(PetscArrayzero(reductions,n)); 2751 PetscCall(MatDenseGetArrayRead(A,&a)); 2752 if (type == NORM_2) { 2753 for (i=0; i<n; i++) { 2754 for (j=0; j<m; j++) { 2755 reductions[i] += PetscAbsScalar(a[j]*a[j]); 2756 } 2757 a += m; 2758 } 2759 } else if (type == NORM_1) { 2760 for (i=0; i<n; i++) { 2761 for (j=0; j<m; j++) { 2762 reductions[i] += PetscAbsScalar(a[j]); 2763 } 2764 a += m; 2765 } 2766 } else if (type == NORM_INFINITY) { 2767 for (i=0; i<n; i++) { 2768 for (j=0; j<m; j++) { 2769 reductions[i] = PetscMax(PetscAbsScalar(a[j]),reductions[i]); 2770 } 2771 a += m; 2772 } 2773 } else if (type == REDUCTION_SUM_REALPART || type == REDUCTION_MEAN_REALPART) { 2774 for (i=0; i<n; i++) { 2775 for (j=0; j<m; j++) { 2776 reductions[i] += PetscRealPart(a[j]); 2777 } 2778 a += m; 2779 } 2780 } else if (type == REDUCTION_SUM_IMAGINARYPART || type == REDUCTION_MEAN_IMAGINARYPART) { 2781 for (i=0; i<n; i++) { 2782 for (j=0; j<m; j++) { 2783 reductions[i] += PetscImaginaryPart(a[j]); 2784 } 2785 a += m; 2786 } 2787 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown reduction type"); 2788 PetscCall(MatDenseRestoreArrayRead(A,&a)); 2789 if (type == NORM_2) { 2790 for (i=0; i<n; i++) reductions[i] = PetscSqrtReal(reductions[i]); 2791 } else if (type == REDUCTION_MEAN_REALPART || type == REDUCTION_MEAN_IMAGINARYPART) { 2792 for (i=0; i<n; i++) reductions[i] /= m; 2793 } 2794 PetscFunctionReturn(0); 2795 } 2796 2797 static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx) 2798 { 2799 PetscScalar *a; 2800 PetscInt lda,m,n,i,j; 2801 2802 PetscFunctionBegin; 2803 PetscCall(MatGetSize(x,&m,&n)); 2804 PetscCall(MatDenseGetLDA(x,&lda)); 2805 PetscCall(MatDenseGetArray(x,&a)); 2806 for (j=0; j<n; j++) { 2807 for (i=0; i<m; i++) { 2808 PetscCall(PetscRandomGetValue(rctx,a+j*lda+i)); 2809 } 2810 } 2811 PetscCall(MatDenseRestoreArray(x,&a)); 2812 PetscFunctionReturn(0); 2813 } 2814 2815 static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool *missing,PetscInt *d) 2816 { 2817 PetscFunctionBegin; 2818 *missing = PETSC_FALSE; 2819 PetscFunctionReturn(0); 2820 } 2821 2822 /* vals is not const */ 2823 static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals) 2824 { 2825 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2826 PetscScalar *v; 2827 2828 PetscFunctionBegin; 2829 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2830 PetscCall(MatDenseGetArray(A,&v)); 2831 *vals = v+col*a->lda; 2832 PetscCall(MatDenseRestoreArray(A,&v)); 2833 PetscFunctionReturn(0); 2834 } 2835 2836 static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals) 2837 { 2838 PetscFunctionBegin; 2839 *vals = NULL; /* user cannot accidentally use the array later */ 2840 PetscFunctionReturn(0); 2841 } 2842 2843 /* -------------------------------------------------------------------*/ 2844 static struct _MatOps MatOps_Values = { MatSetValues_SeqDense, 2845 MatGetRow_SeqDense, 2846 MatRestoreRow_SeqDense, 2847 MatMult_SeqDense, 2848 /* 4*/ MatMultAdd_SeqDense, 2849 MatMultTranspose_SeqDense, 2850 MatMultTransposeAdd_SeqDense, 2851 NULL, 2852 NULL, 2853 NULL, 2854 /* 10*/ NULL, 2855 MatLUFactor_SeqDense, 2856 MatCholeskyFactor_SeqDense, 2857 MatSOR_SeqDense, 2858 MatTranspose_SeqDense, 2859 /* 15*/ MatGetInfo_SeqDense, 2860 MatEqual_SeqDense, 2861 MatGetDiagonal_SeqDense, 2862 MatDiagonalScale_SeqDense, 2863 MatNorm_SeqDense, 2864 /* 20*/ MatAssemblyBegin_SeqDense, 2865 MatAssemblyEnd_SeqDense, 2866 MatSetOption_SeqDense, 2867 MatZeroEntries_SeqDense, 2868 /* 24*/ MatZeroRows_SeqDense, 2869 NULL, 2870 NULL, 2871 NULL, 2872 NULL, 2873 /* 29*/ MatSetUp_SeqDense, 2874 NULL, 2875 NULL, 2876 NULL, 2877 NULL, 2878 /* 34*/ MatDuplicate_SeqDense, 2879 NULL, 2880 NULL, 2881 NULL, 2882 NULL, 2883 /* 39*/ MatAXPY_SeqDense, 2884 MatCreateSubMatrices_SeqDense, 2885 NULL, 2886 MatGetValues_SeqDense, 2887 MatCopy_SeqDense, 2888 /* 44*/ MatGetRowMax_SeqDense, 2889 MatScale_SeqDense, 2890 MatShift_SeqDense, 2891 NULL, 2892 MatZeroRowsColumns_SeqDense, 2893 /* 49*/ MatSetRandom_SeqDense, 2894 NULL, 2895 NULL, 2896 NULL, 2897 NULL, 2898 /* 54*/ NULL, 2899 NULL, 2900 NULL, 2901 NULL, 2902 NULL, 2903 /* 59*/ MatCreateSubMatrix_SeqDense, 2904 MatDestroy_SeqDense, 2905 MatView_SeqDense, 2906 NULL, 2907 NULL, 2908 /* 64*/ NULL, 2909 NULL, 2910 NULL, 2911 NULL, 2912 NULL, 2913 /* 69*/ MatGetRowMaxAbs_SeqDense, 2914 NULL, 2915 NULL, 2916 NULL, 2917 NULL, 2918 /* 74*/ NULL, 2919 NULL, 2920 NULL, 2921 NULL, 2922 NULL, 2923 /* 79*/ NULL, 2924 NULL, 2925 NULL, 2926 NULL, 2927 /* 83*/ MatLoad_SeqDense, 2928 MatIsSymmetric_SeqDense, 2929 MatIsHermitian_SeqDense, 2930 NULL, 2931 NULL, 2932 NULL, 2933 /* 89*/ NULL, 2934 NULL, 2935 MatMatMultNumeric_SeqDense_SeqDense, 2936 NULL, 2937 NULL, 2938 /* 94*/ NULL, 2939 NULL, 2940 NULL, 2941 MatMatTransposeMultNumeric_SeqDense_SeqDense, 2942 NULL, 2943 /* 99*/ MatProductSetFromOptions_SeqDense, 2944 NULL, 2945 NULL, 2946 MatConjugate_SeqDense, 2947 NULL, 2948 /*104*/ NULL, 2949 MatRealPart_SeqDense, 2950 MatImaginaryPart_SeqDense, 2951 NULL, 2952 NULL, 2953 /*109*/ NULL, 2954 NULL, 2955 MatGetRowMin_SeqDense, 2956 MatGetColumnVector_SeqDense, 2957 MatMissingDiagonal_SeqDense, 2958 /*114*/ NULL, 2959 NULL, 2960 NULL, 2961 NULL, 2962 NULL, 2963 /*119*/ NULL, 2964 NULL, 2965 NULL, 2966 NULL, 2967 NULL, 2968 /*124*/ NULL, 2969 MatGetColumnReductions_SeqDense, 2970 NULL, 2971 NULL, 2972 NULL, 2973 /*129*/ NULL, 2974 NULL, 2975 NULL, 2976 MatTransposeMatMultNumeric_SeqDense_SeqDense, 2977 NULL, 2978 /*134*/ NULL, 2979 NULL, 2980 NULL, 2981 NULL, 2982 NULL, 2983 /*139*/ NULL, 2984 NULL, 2985 NULL, 2986 NULL, 2987 NULL, 2988 MatCreateMPIMatConcatenateSeqMat_SeqDense, 2989 /*145*/ NULL, 2990 NULL, 2991 NULL 2992 }; 2993 2994 /*@C 2995 MatCreateSeqDense - Creates a sequential dense matrix that 2996 is stored in column major order (the usual Fortran 77 manner). Many 2997 of the matrix operations use the BLAS and LAPACK routines. 2998 2999 Collective 3000 3001 Input Parameters: 3002 + comm - MPI communicator, set to PETSC_COMM_SELF 3003 . m - number of rows 3004 . n - number of columns 3005 - data - optional location of matrix data in column major order. Set data=NULL for PETSc 3006 to control all matrix memory allocation. 3007 3008 Output Parameter: 3009 . A - the matrix 3010 3011 Notes: 3012 The data input variable is intended primarily for Fortran programmers 3013 who wish to allocate their own matrix memory space. Most users should 3014 set data=NULL. 3015 3016 Level: intermediate 3017 3018 .seealso: `MatCreate()`, `MatCreateDense()`, `MatSetValues()` 3019 @*/ 3020 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 3021 { 3022 PetscFunctionBegin; 3023 PetscCall(MatCreate(comm,A)); 3024 PetscCall(MatSetSizes(*A,m,n,m,n)); 3025 PetscCall(MatSetType(*A,MATSEQDENSE)); 3026 PetscCall(MatSeqDenseSetPreallocation(*A,data)); 3027 PetscFunctionReturn(0); 3028 } 3029 3030 /*@C 3031 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 3032 3033 Collective 3034 3035 Input Parameters: 3036 + B - the matrix 3037 - data - the array (or NULL) 3038 3039 Notes: 3040 The data input variable is intended primarily for Fortran programmers 3041 who wish to allocate their own matrix memory space. Most users should 3042 need not call this routine. 3043 3044 Level: intermediate 3045 3046 .seealso: `MatCreate()`, `MatCreateDense()`, `MatSetValues()`, `MatDenseSetLDA()` 3047 3048 @*/ 3049 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 3050 { 3051 PetscFunctionBegin; 3052 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3053 PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data)); 3054 PetscFunctionReturn(0); 3055 } 3056 3057 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 3058 { 3059 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 3060 3061 PetscFunctionBegin; 3062 PetscCheck(!b->matinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 3063 B->preallocated = PETSC_TRUE; 3064 3065 PetscCall(PetscLayoutSetUp(B->rmap)); 3066 PetscCall(PetscLayoutSetUp(B->cmap)); 3067 3068 if (b->lda <= 0) b->lda = B->rmap->n; 3069 3070 if (!data) { /* petsc-allocated storage */ 3071 if (!b->user_alloc) PetscCall(PetscFree(b->v)); 3072 PetscCall(PetscCalloc1((size_t)b->lda*B->cmap->n,&b->v)); 3073 PetscCall(PetscLogObjectMemory((PetscObject)B,b->lda*B->cmap->n*sizeof(PetscScalar))); 3074 3075 b->user_alloc = PETSC_FALSE; 3076 } else { /* user-allocated storage */ 3077 if (!b->user_alloc) PetscCall(PetscFree(b->v)); 3078 b->v = data; 3079 b->user_alloc = PETSC_TRUE; 3080 } 3081 B->assembled = PETSC_TRUE; 3082 PetscFunctionReturn(0); 3083 } 3084 3085 #if defined(PETSC_HAVE_ELEMENTAL) 3086 PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 3087 { 3088 Mat mat_elemental; 3089 const PetscScalar *array; 3090 PetscScalar *v_colwise; 3091 PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols; 3092 3093 PetscFunctionBegin; 3094 PetscCall(PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols)); 3095 PetscCall(MatDenseGetArrayRead(A,&array)); 3096 /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */ 3097 k = 0; 3098 for (j=0; j<N; j++) { 3099 cols[j] = j; 3100 for (i=0; i<M; i++) { 3101 v_colwise[j*M+i] = array[k++]; 3102 } 3103 } 3104 for (i=0; i<M; i++) { 3105 rows[i] = i; 3106 } 3107 PetscCall(MatDenseRestoreArrayRead(A,&array)); 3108 3109 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental)); 3110 PetscCall(MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N)); 3111 PetscCall(MatSetType(mat_elemental,MATELEMENTAL)); 3112 PetscCall(MatSetUp(mat_elemental)); 3113 3114 /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */ 3115 PetscCall(MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES)); 3116 PetscCall(MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY)); 3117 PetscCall(MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY)); 3118 PetscCall(PetscFree3(v_colwise,rows,cols)); 3119 3120 if (reuse == MAT_INPLACE_MATRIX) { 3121 PetscCall(MatHeaderReplace(A,&mat_elemental)); 3122 } else { 3123 *newmat = mat_elemental; 3124 } 3125 PetscFunctionReturn(0); 3126 } 3127 #endif 3128 3129 PetscErrorCode MatDenseSetLDA_SeqDense(Mat B,PetscInt lda) 3130 { 3131 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 3132 PetscBool data; 3133 3134 PetscFunctionBegin; 3135 data = (PetscBool)((B->rmap->n > 0 && B->cmap->n > 0) ? (b->v ? PETSC_TRUE : PETSC_FALSE) : PETSC_FALSE); 3136 PetscCheck(b->user_alloc || !data || b->lda == lda,PETSC_COMM_SELF,PETSC_ERR_ORDER,"LDA cannot be changed after allocation of internal storage"); 3137 PetscCheck(lda >= B->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"LDA %" PetscInt_FMT " must be at least matrix dimension %" PetscInt_FMT,lda,B->rmap->n); 3138 b->lda = lda; 3139 PetscFunctionReturn(0); 3140 } 3141 3142 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 3143 { 3144 PetscMPIInt size; 3145 3146 PetscFunctionBegin; 3147 PetscCallMPI(MPI_Comm_size(comm,&size)); 3148 if (size == 1) { 3149 if (scall == MAT_INITIAL_MATRIX) { 3150 PetscCall(MatDuplicate(inmat,MAT_COPY_VALUES,outmat)); 3151 } else { 3152 PetscCall(MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN)); 3153 } 3154 } else { 3155 PetscCall(MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat)); 3156 } 3157 PetscFunctionReturn(0); 3158 } 3159 3160 PetscErrorCode MatDenseGetColumnVec_SeqDense(Mat A,PetscInt col,Vec *v) 3161 { 3162 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 3163 3164 PetscFunctionBegin; 3165 PetscCheck(!a->vecinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 3166 PetscCheck(!a->matinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 3167 if (!a->cvec) { 3168 PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec)); 3169 PetscCall(PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec)); 3170 } 3171 a->vecinuse = col + 1; 3172 PetscCall(MatDenseGetArray(A,(PetscScalar**)&a->ptrinuse)); 3173 PetscCall(VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda)); 3174 *v = a->cvec; 3175 PetscFunctionReturn(0); 3176 } 3177 3178 PetscErrorCode MatDenseRestoreColumnVec_SeqDense(Mat A,PetscInt col,Vec *v) 3179 { 3180 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 3181 3182 PetscFunctionBegin; 3183 PetscCheck(a->vecinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first"); 3184 PetscCheck(a->cvec,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector"); 3185 a->vecinuse = 0; 3186 PetscCall(MatDenseRestoreArray(A,(PetscScalar**)&a->ptrinuse)); 3187 PetscCall(VecResetArray(a->cvec)); 3188 if (v) *v = NULL; 3189 PetscFunctionReturn(0); 3190 } 3191 3192 PetscErrorCode MatDenseGetColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v) 3193 { 3194 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 3195 3196 PetscFunctionBegin; 3197 PetscCheck(!a->vecinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 3198 PetscCheck(!a->matinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 3199 if (!a->cvec) { 3200 PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec)); 3201 PetscCall(PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec)); 3202 } 3203 a->vecinuse = col + 1; 3204 PetscCall(MatDenseGetArrayRead(A,&a->ptrinuse)); 3205 PetscCall(VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda)); 3206 PetscCall(VecLockReadPush(a->cvec)); 3207 *v = a->cvec; 3208 PetscFunctionReturn(0); 3209 } 3210 3211 PetscErrorCode MatDenseRestoreColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v) 3212 { 3213 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 3214 3215 PetscFunctionBegin; 3216 PetscCheck(a->vecinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first"); 3217 PetscCheck(a->cvec,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector"); 3218 a->vecinuse = 0; 3219 PetscCall(MatDenseRestoreArrayRead(A,&a->ptrinuse)); 3220 PetscCall(VecLockReadPop(a->cvec)); 3221 PetscCall(VecResetArray(a->cvec)); 3222 if (v) *v = NULL; 3223 PetscFunctionReturn(0); 3224 } 3225 3226 PetscErrorCode MatDenseGetColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v) 3227 { 3228 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 3229 3230 PetscFunctionBegin; 3231 PetscCheck(!a->vecinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 3232 PetscCheck(!a->matinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 3233 if (!a->cvec) { 3234 PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec)); 3235 PetscCall(PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec)); 3236 } 3237 a->vecinuse = col + 1; 3238 PetscCall(MatDenseGetArrayWrite(A,(PetscScalar**)&a->ptrinuse)); 3239 PetscCall(VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda)); 3240 *v = a->cvec; 3241 PetscFunctionReturn(0); 3242 } 3243 3244 PetscErrorCode MatDenseRestoreColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v) 3245 { 3246 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 3247 3248 PetscFunctionBegin; 3249 PetscCheck(a->vecinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first"); 3250 PetscCheck(a->cvec,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector"); 3251 a->vecinuse = 0; 3252 PetscCall(MatDenseRestoreArrayWrite(A,(PetscScalar**)&a->ptrinuse)); 3253 PetscCall(VecResetArray(a->cvec)); 3254 if (v) *v = NULL; 3255 PetscFunctionReturn(0); 3256 } 3257 3258 PetscErrorCode MatDenseGetSubMatrix_SeqDense(Mat A,PetscInt rbegin,PetscInt rend,PetscInt cbegin,PetscInt cend,Mat *v) 3259 { 3260 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 3261 3262 PetscFunctionBegin; 3263 PetscCheck(!a->vecinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 3264 PetscCheck(!a->matinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 3265 if (a->cmat && (cend-cbegin != a->cmat->cmap->N || rend-rbegin != a->cmat->rmap->N)) PetscCall(MatDestroy(&a->cmat)); 3266 if (!a->cmat) { 3267 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A),rend-rbegin,PETSC_DECIDE,rend-rbegin,cend-cbegin,a->v+rbegin+(size_t)cbegin*a->lda,&a->cmat)); 3268 PetscCall(PetscLogObjectParent((PetscObject)A,(PetscObject)a->cmat)); 3269 } else { 3270 PetscCall(MatDensePlaceArray(a->cmat,a->v+rbegin+(size_t)cbegin*a->lda)); 3271 } 3272 PetscCall(MatDenseSetLDA(a->cmat,a->lda)); 3273 a->matinuse = cbegin + 1; 3274 *v = a->cmat; 3275 #if defined(PETSC_HAVE_CUDA) 3276 A->offloadmask = PETSC_OFFLOAD_CPU; 3277 #endif 3278 PetscFunctionReturn(0); 3279 } 3280 3281 PetscErrorCode MatDenseRestoreSubMatrix_SeqDense(Mat A,Mat *v) 3282 { 3283 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 3284 3285 PetscFunctionBegin; 3286 PetscCheck(a->matinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetSubMatrix() first"); 3287 PetscCheck(a->cmat,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column matrix"); 3288 PetscCheck(*v == a->cmat,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not the matrix obtained from MatDenseGetSubMatrix()"); 3289 a->matinuse = 0; 3290 PetscCall(MatDenseResetArray(a->cmat)); 3291 *v = NULL; 3292 PetscFunctionReturn(0); 3293 } 3294 3295 /*MC 3296 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 3297 3298 Options Database Keys: 3299 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 3300 3301 Level: beginner 3302 3303 .seealso: `MatCreateSeqDense()` 3304 3305 M*/ 3306 PetscErrorCode MatCreate_SeqDense(Mat B) 3307 { 3308 Mat_SeqDense *b; 3309 PetscMPIInt size; 3310 3311 PetscFunctionBegin; 3312 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B),&size)); 3313 PetscCheck(size <= 1,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 3314 3315 PetscCall(PetscNewLog(B,&b)); 3316 PetscCall(PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps))); 3317 B->data = (void*)b; 3318 3319 b->roworiented = PETSC_TRUE; 3320 3321 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatQRFactor_C",MatQRFactor_SeqDense)); 3322 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseGetLDA_C",MatDenseGetLDA_SeqDense)); 3323 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseSetLDA_C",MatDenseSetLDA_SeqDense)); 3324 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense)); 3325 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense)); 3326 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense)); 3327 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense)); 3328 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseReplaceArray_C",MatDenseReplaceArray_SeqDense)); 3329 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense)); 3330 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense)); 3331 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayWrite_C",MatDenseGetArray_SeqDense)); 3332 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayWrite_C",MatDenseRestoreArray_SeqDense)); 3333 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ)); 3334 #if defined(PETSC_HAVE_ELEMENTAL) 3335 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental)); 3336 #endif 3337 #if defined(PETSC_HAVE_SCALAPACK) 3338 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_scalapack_C",MatConvert_Dense_ScaLAPACK)); 3339 #endif 3340 #if defined(PETSC_HAVE_CUDA) 3341 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqdensecuda_C",MatConvert_SeqDense_SeqDenseCUDA)); 3342 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",MatProductSetFromOptions_SeqDense)); 3343 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdense_C",MatProductSetFromOptions_SeqDense)); 3344 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdensecuda_C",MatProductSetFromOptions_SeqDense)); 3345 #endif 3346 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense)); 3347 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqdense_C",MatProductSetFromOptions_SeqAIJ_SeqDense)); 3348 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdense_C",MatProductSetFromOptions_SeqDense)); 3349 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense)); 3350 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqsbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense)); 3351 3352 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense)); 3353 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense)); 3354 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVec_C",MatDenseGetColumnVec_SeqDense)); 3355 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVec_C",MatDenseRestoreColumnVec_SeqDense)); 3356 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecRead_C",MatDenseGetColumnVecRead_SeqDense)); 3357 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecRead_C",MatDenseRestoreColumnVecRead_SeqDense)); 3358 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecWrite_C",MatDenseGetColumnVecWrite_SeqDense)); 3359 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecWrite_C",MatDenseRestoreColumnVecWrite_SeqDense)); 3360 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseGetSubMatrix_C",MatDenseGetSubMatrix_SeqDense)); 3361 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreSubMatrix_C",MatDenseRestoreSubMatrix_SeqDense)); 3362 PetscCall(PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE)); 3363 PetscFunctionReturn(0); 3364 } 3365 3366 /*@C 3367 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. 3368 3369 Not Collective 3370 3371 Input Parameters: 3372 + mat - a MATSEQDENSE or MATMPIDENSE matrix 3373 - col - column index 3374 3375 Output Parameter: 3376 . vals - pointer to the data 3377 3378 Level: intermediate 3379 3380 .seealso: `MatDenseRestoreColumn()` 3381 @*/ 3382 PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals) 3383 { 3384 PetscFunctionBegin; 3385 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3386 PetscValidLogicalCollectiveInt(A,col,2); 3387 PetscValidPointer(vals,3); 3388 PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals)); 3389 PetscFunctionReturn(0); 3390 } 3391 3392 /*@C 3393 MatDenseRestoreColumn - returns access to a column of a dense matrix which is returned by MatDenseGetColumn(). 3394 3395 Not Collective 3396 3397 Input Parameter: 3398 . mat - a MATSEQDENSE or MATMPIDENSE matrix 3399 3400 Output Parameter: 3401 . vals - pointer to the data 3402 3403 Level: intermediate 3404 3405 .seealso: `MatDenseGetColumn()` 3406 @*/ 3407 PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals) 3408 { 3409 PetscFunctionBegin; 3410 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3411 PetscValidPointer(vals,2); 3412 PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals)); 3413 PetscFunctionReturn(0); 3414 } 3415 3416 /*@ 3417 MatDenseGetColumnVec - Gives read-write access to a column of a dense matrix, represented as a Vec. 3418 3419 Collective 3420 3421 Input Parameters: 3422 + mat - the Mat object 3423 - col - the column index 3424 3425 Output Parameter: 3426 . v - the vector 3427 3428 Notes: 3429 The vector is owned by PETSc. Users need to call MatDenseRestoreColumnVec() when the vector is no longer needed. 3430 Use MatDenseGetColumnVecRead() to obtain read-only access or MatDenseGetColumnVecWrite() for write-only access. 3431 3432 Level: intermediate 3433 3434 .seealso: `MATDENSE`, `MATDENSECUDA`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()` 3435 @*/ 3436 PetscErrorCode MatDenseGetColumnVec(Mat A,PetscInt col,Vec *v) 3437 { 3438 PetscFunctionBegin; 3439 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3440 PetscValidType(A,1); 3441 PetscValidLogicalCollectiveInt(A,col,2); 3442 PetscValidPointer(v,3); 3443 PetscCheck(A->preallocated,PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3444 PetscCheck(col >= 0 && col < A->cmap->N,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %" PetscInt_FMT ", should be in [0,%" PetscInt_FMT ")",col,A->cmap->N); 3445 PetscUseMethod(A,"MatDenseGetColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v)); 3446 PetscFunctionReturn(0); 3447 } 3448 3449 /*@ 3450 MatDenseRestoreColumnVec - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVec(). 3451 3452 Collective 3453 3454 Input Parameters: 3455 + mat - the Mat object 3456 . col - the column index 3457 - v - the Vec object 3458 3459 Level: intermediate 3460 3461 .seealso: `MATDENSE`, `MATDENSECUDA`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()` 3462 @*/ 3463 PetscErrorCode MatDenseRestoreColumnVec(Mat A,PetscInt col,Vec *v) 3464 { 3465 PetscFunctionBegin; 3466 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3467 PetscValidType(A,1); 3468 PetscValidLogicalCollectiveInt(A,col,2); 3469 PetscCheck(A->preallocated,PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3470 PetscCheck(col >= 0 && col < A->cmap->N,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %" PetscInt_FMT ", should be in [0,%" PetscInt_FMT ")",col,A->cmap->N); 3471 PetscUseMethod(A,"MatDenseRestoreColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v)); 3472 PetscFunctionReturn(0); 3473 } 3474 3475 /*@ 3476 MatDenseGetColumnVecRead - Gives read-only access to a column of a dense matrix, represented as a Vec. 3477 3478 Collective 3479 3480 Input Parameters: 3481 + mat - the Mat object 3482 - col - the column index 3483 3484 Output Parameter: 3485 . v - the vector 3486 3487 Notes: 3488 The vector is owned by PETSc and users cannot modify it. 3489 Users need to call MatDenseRestoreColumnVecRead() when the vector is no longer needed. 3490 Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecWrite() for write-only access. 3491 3492 Level: intermediate 3493 3494 .seealso: `MATDENSE`, `MATDENSECUDA`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()` 3495 @*/ 3496 PetscErrorCode MatDenseGetColumnVecRead(Mat A,PetscInt col,Vec *v) 3497 { 3498 PetscFunctionBegin; 3499 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3500 PetscValidType(A,1); 3501 PetscValidLogicalCollectiveInt(A,col,2); 3502 PetscValidPointer(v,3); 3503 PetscCheck(A->preallocated,PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3504 PetscCheck(col >= 0 && col < A->cmap->N,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %" PetscInt_FMT ", should be in [0,%" PetscInt_FMT ")",col,A->cmap->N); 3505 PetscUseMethod(A,"MatDenseGetColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v)); 3506 PetscFunctionReturn(0); 3507 } 3508 3509 /*@ 3510 MatDenseRestoreColumnVecRead - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecRead(). 3511 3512 Collective 3513 3514 Input Parameters: 3515 + mat - the Mat object 3516 . col - the column index 3517 - v - the Vec object 3518 3519 Level: intermediate 3520 3521 .seealso: `MATDENSE`, `MATDENSECUDA`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecWrite()` 3522 @*/ 3523 PetscErrorCode MatDenseRestoreColumnVecRead(Mat A,PetscInt col,Vec *v) 3524 { 3525 PetscFunctionBegin; 3526 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3527 PetscValidType(A,1); 3528 PetscValidLogicalCollectiveInt(A,col,2); 3529 PetscCheck(A->preallocated,PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3530 PetscCheck(col >= 0 && col < A->cmap->N,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %" PetscInt_FMT ", should be in [0,%" PetscInt_FMT ")",col,A->cmap->N); 3531 PetscUseMethod(A,"MatDenseRestoreColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v)); 3532 PetscFunctionReturn(0); 3533 } 3534 3535 /*@ 3536 MatDenseGetColumnVecWrite - Gives write-only access to a column of a dense matrix, represented as a Vec. 3537 3538 Collective 3539 3540 Input Parameters: 3541 + mat - the Mat object 3542 - col - the column index 3543 3544 Output Parameter: 3545 . v - the vector 3546 3547 Notes: 3548 The vector is owned by PETSc. Users need to call MatDenseRestoreColumnVecWrite() when the vector is no longer needed. 3549 Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecRead() for read-only access. 3550 3551 Level: intermediate 3552 3553 .seealso: `MATDENSE`, `MATDENSECUDA`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()` 3554 @*/ 3555 PetscErrorCode MatDenseGetColumnVecWrite(Mat A,PetscInt col,Vec *v) 3556 { 3557 PetscFunctionBegin; 3558 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3559 PetscValidType(A,1); 3560 PetscValidLogicalCollectiveInt(A,col,2); 3561 PetscValidPointer(v,3); 3562 PetscCheck(A->preallocated,PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3563 PetscCheck(col >= 0 && col < A->cmap->N,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %" PetscInt_FMT ", should be in [0,%" PetscInt_FMT ")",col,A->cmap->N); 3564 PetscUseMethod(A,"MatDenseGetColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v)); 3565 PetscFunctionReturn(0); 3566 } 3567 3568 /*@ 3569 MatDenseRestoreColumnVecWrite - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecWrite(). 3570 3571 Collective 3572 3573 Input Parameters: 3574 + mat - the Mat object 3575 . col - the column index 3576 - v - the Vec object 3577 3578 Level: intermediate 3579 3580 .seealso: `MATDENSE`, `MATDENSECUDA`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()` 3581 @*/ 3582 PetscErrorCode MatDenseRestoreColumnVecWrite(Mat A,PetscInt col,Vec *v) 3583 { 3584 PetscFunctionBegin; 3585 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3586 PetscValidType(A,1); 3587 PetscValidLogicalCollectiveInt(A,col,2); 3588 PetscCheck(A->preallocated,PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3589 PetscCheck(col >= 0 && col < A->cmap->N,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %" PetscInt_FMT ", should be in [0,%" PetscInt_FMT ")",col,A->cmap->N); 3590 PetscUseMethod(A,"MatDenseRestoreColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v)); 3591 PetscFunctionReturn(0); 3592 } 3593 3594 /*@ 3595 MatDenseGetSubMatrix - Gives access to a block of rows and columns of a dense matrix, represented as a Mat. 3596 3597 Collective 3598 3599 Input Parameters: 3600 + mat - the Mat object 3601 . rbegin - the first global row index in the block (if PETSC_DECIDE, is 0) 3602 . rend - the global row index past the last one in the block (if PETSC_DECIDE, is M) 3603 . cbegin - the first global column index in the block (if PETSC_DECIDE, is 0) 3604 - cend - the global column index past the last one in the block (if PETSC_DECIDE, is N) 3605 3606 Output Parameter: 3607 . v - the matrix 3608 3609 Notes: 3610 The matrix is owned by PETSc. Users need to call MatDenseRestoreSubMatrix() when the matrix is no longer needed. 3611 The output matrix is not redistributed by PETSc, so depending on the values of rbegin and rend, some processes may have no local rows. 3612 3613 Level: intermediate 3614 3615 .seealso: `MATDENSE`, `MATDENSECUDA`, `MatDenseGetColumnVec()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreSubMatrix()` 3616 @*/ 3617 PetscErrorCode MatDenseGetSubMatrix(Mat A,PetscInt rbegin,PetscInt rend,PetscInt cbegin,PetscInt cend,Mat *v) 3618 { 3619 PetscFunctionBegin; 3620 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3621 PetscValidType(A,1); 3622 PetscValidLogicalCollectiveInt(A,rbegin,2); 3623 PetscValidLogicalCollectiveInt(A,rend,3); 3624 PetscValidLogicalCollectiveInt(A,cbegin,4); 3625 PetscValidLogicalCollectiveInt(A,cend,5); 3626 PetscValidPointer(v,6); 3627 if (rbegin == PETSC_DECIDE) rbegin = 0; 3628 if (rend == PETSC_DECIDE) rend = A->rmap->N; 3629 if (cbegin == PETSC_DECIDE) cbegin = 0; 3630 if (cend == PETSC_DECIDE) cend = A->cmap->N; 3631 PetscCheck(A->preallocated,PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3632 PetscCheck(rbegin >= 0 && rbegin <= A->rmap->N,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid rbegin %" PetscInt_FMT ", should be in [0,%" PetscInt_FMT "]",rbegin,A->rmap->N); 3633 PetscCheck(rend >= rbegin && rend <= A->rmap->N,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid rend %" PetscInt_FMT ", should be in [%" PetscInt_FMT ",%" PetscInt_FMT "]",rend,rbegin,A->rmap->N); 3634 PetscCheck(cbegin >= 0 && cbegin <= A->cmap->N,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid cbegin %" PetscInt_FMT ", should be in [0,%" PetscInt_FMT "]",cbegin,A->cmap->N); 3635 PetscCheck(cend >= cbegin && cend <= A->cmap->N,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid cend %" PetscInt_FMT ", should be in [%" PetscInt_FMT ",%" PetscInt_FMT "]",cend,cbegin,A->cmap->N); 3636 PetscUseMethod(A,"MatDenseGetSubMatrix_C",(Mat,PetscInt,PetscInt,PetscInt,PetscInt,Mat*),(A,rbegin,rend,cbegin,cend,v)); 3637 PetscFunctionReturn(0); 3638 } 3639 3640 /*@ 3641 MatDenseRestoreSubMatrix - Returns access to a block of columns of a dense matrix obtained from MatDenseGetSubMatrix(). 3642 3643 Collective 3644 3645 Input Parameters: 3646 + mat - the Mat object 3647 - v - the Mat object 3648 3649 Level: intermediate 3650 3651 .seealso: `MATDENSE`, `MATDENSECUDA`, `MatDenseGetColumnVec()`, `MatDenseRestoreColumnVec()`, `MatDenseGetSubMatrix()` 3652 @*/ 3653 PetscErrorCode MatDenseRestoreSubMatrix(Mat A,Mat *v) 3654 { 3655 PetscFunctionBegin; 3656 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3657 PetscValidType(A,1); 3658 PetscValidPointer(v,2); 3659 PetscUseMethod(A,"MatDenseRestoreSubMatrix_C",(Mat,Mat*),(A,v)); 3660 PetscFunctionReturn(0); 3661 } 3662 3663 #include <petscblaslapack.h> 3664 #include <petsc/private/kernels/blockinvert.h> 3665 3666 PetscErrorCode MatSeqDenseInvert(Mat A) 3667 { 3668 Mat_SeqDense *a = (Mat_SeqDense*) A->data; 3669 PetscInt bs = A->rmap->n; 3670 MatScalar *values = a->v; 3671 const PetscReal shift = 0.0; 3672 PetscBool allowzeropivot = PetscNot(A->erroriffailure),zeropivotdetected=PETSC_FALSE; 3673 3674 PetscFunctionBegin; 3675 /* factor and invert each block */ 3676 switch (bs) { 3677 case 1: 3678 values[0] = (PetscScalar)1.0 / (values[0] + shift); 3679 break; 3680 case 2: 3681 PetscCall(PetscKernel_A_gets_inverse_A_2(values,shift,allowzeropivot,&zeropivotdetected)); 3682 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3683 break; 3684 case 3: 3685 PetscCall(PetscKernel_A_gets_inverse_A_3(values,shift,allowzeropivot,&zeropivotdetected)); 3686 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3687 break; 3688 case 4: 3689 PetscCall(PetscKernel_A_gets_inverse_A_4(values,shift,allowzeropivot,&zeropivotdetected)); 3690 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3691 break; 3692 case 5: 3693 { 3694 PetscScalar work[25]; 3695 PetscInt ipvt[5]; 3696 3697 PetscCall(PetscKernel_A_gets_inverse_A_5(values,ipvt,work,shift,allowzeropivot,&zeropivotdetected)); 3698 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3699 } 3700 break; 3701 case 6: 3702 PetscCall(PetscKernel_A_gets_inverse_A_6(values,shift,allowzeropivot,&zeropivotdetected)); 3703 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3704 break; 3705 case 7: 3706 PetscCall(PetscKernel_A_gets_inverse_A_7(values,shift,allowzeropivot,&zeropivotdetected)); 3707 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3708 break; 3709 default: 3710 { 3711 PetscInt *v_pivots,*IJ,j; 3712 PetscScalar *v_work; 3713 3714 PetscCall(PetscMalloc3(bs,&v_work,bs,&v_pivots,bs,&IJ)); 3715 for (j=0; j<bs; j++) { 3716 IJ[j] = j; 3717 } 3718 PetscCall(PetscKernel_A_gets_inverse_A(bs,values,v_pivots,v_work,allowzeropivot,&zeropivotdetected)); 3719 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3720 PetscCall(PetscFree3(v_work,v_pivots,IJ)); 3721 } 3722 } 3723 PetscFunctionReturn(0); 3724 } 3725