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