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