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