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,"MatQRFactorSymbolic_C",NULL)); 1744 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatQRFactorNumeric_C",NULL)); 1745 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetLDA_C",NULL)); 1746 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseSetLDA_C",NULL)); 1747 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL)); 1748 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL)); 1749 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL)); 1750 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL)); 1751 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseReplaceArray_C",NULL)); 1752 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayRead_C",NULL)); 1753 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayRead_C",NULL)); 1754 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayWrite_C",NULL)); 1755 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayWrite_C",NULL)); 1756 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL)); 1757 #if defined(PETSC_HAVE_ELEMENTAL) 1758 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL)); 1759 #endif 1760 #if defined(PETSC_HAVE_SCALAPACK) 1761 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_scalapack_C",NULL)); 1762 #endif 1763 #if defined(PETSC_HAVE_CUDA) 1764 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqdensecuda_C",NULL)); 1765 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",NULL)); 1766 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdense_C",NULL)); 1767 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdense_seqdensecuda_C",NULL)); 1768 #endif 1769 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL)); 1770 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqaij_seqdense_C",NULL)); 1771 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdense_seqdense_C",NULL)); 1772 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqbaij_seqdense_C",NULL)); 1773 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqsbaij_seqdense_C",NULL)); 1774 1775 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumn_C",NULL)); 1776 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumn_C",NULL)); 1777 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVec_C",NULL)); 1778 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVec_C",NULL)); 1779 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVecRead_C",NULL)); 1780 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVecRead_C",NULL)); 1781 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVecWrite_C",NULL)); 1782 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVecWrite_C",NULL)); 1783 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetSubMatrix_C",NULL)); 1784 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreSubMatrix_C",NULL)); 1785 PetscFunctionReturn(0); 1786 } 1787 1788 static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout) 1789 { 1790 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1791 PetscInt k,j,m = A->rmap->n, M = mat->lda, n = A->cmap->n; 1792 PetscScalar *v,tmp; 1793 1794 PetscFunctionBegin; 1795 if (reuse == MAT_INPLACE_MATRIX) { 1796 if (m == n) { /* in place transpose */ 1797 PetscCall(MatDenseGetArray(A,&v)); 1798 for (j=0; j<m; j++) { 1799 for (k=0; k<j; k++) { 1800 tmp = v[j + k*M]; 1801 v[j + k*M] = v[k + j*M]; 1802 v[k + j*M] = tmp; 1803 } 1804 } 1805 PetscCall(MatDenseRestoreArray(A,&v)); 1806 } else { /* reuse memory, temporary allocates new memory */ 1807 PetscScalar *v2; 1808 PetscLayout tmplayout; 1809 1810 PetscCall(PetscMalloc1((size_t)m*n,&v2)); 1811 PetscCall(MatDenseGetArray(A,&v)); 1812 for (j=0; j<n; j++) { 1813 for (k=0; k<m; k++) v2[j + (size_t)k*n] = v[k + (size_t)j*M]; 1814 } 1815 PetscCall(PetscArraycpy(v,v2,(size_t)m*n)); 1816 PetscCall(PetscFree(v2)); 1817 PetscCall(MatDenseRestoreArray(A,&v)); 1818 /* cleanup size dependent quantities */ 1819 PetscCall(VecDestroy(&mat->cvec)); 1820 PetscCall(MatDestroy(&mat->cmat)); 1821 PetscCall(PetscFree(mat->pivots)); 1822 PetscCall(PetscFree(mat->fwork)); 1823 PetscCall(MatDestroy(&mat->ptapwork)); 1824 /* swap row/col layouts */ 1825 mat->lda = n; 1826 tmplayout = A->rmap; 1827 A->rmap = A->cmap; 1828 A->cmap = tmplayout; 1829 } 1830 } else { /* out-of-place transpose */ 1831 Mat tmat; 1832 Mat_SeqDense *tmatd; 1833 PetscScalar *v2; 1834 PetscInt M2; 1835 1836 if (reuse == MAT_INITIAL_MATRIX) { 1837 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&tmat)); 1838 PetscCall(MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n)); 1839 PetscCall(MatSetType(tmat,((PetscObject)A)->type_name)); 1840 PetscCall(MatSeqDenseSetPreallocation(tmat,NULL)); 1841 } else tmat = *matout; 1842 1843 PetscCall(MatDenseGetArrayRead(A,(const PetscScalar**)&v)); 1844 PetscCall(MatDenseGetArray(tmat,&v2)); 1845 tmatd = (Mat_SeqDense*)tmat->data; 1846 M2 = tmatd->lda; 1847 for (j=0; j<n; j++) { 1848 for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M]; 1849 } 1850 PetscCall(MatDenseRestoreArray(tmat,&v2)); 1851 PetscCall(MatDenseRestoreArrayRead(A,(const PetscScalar**)&v)); 1852 PetscCall(MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY)); 1853 PetscCall(MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY)); 1854 *matout = tmat; 1855 } 1856 PetscFunctionReturn(0); 1857 } 1858 1859 static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg) 1860 { 1861 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1862 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1863 PetscInt i; 1864 const PetscScalar *v1,*v2; 1865 1866 PetscFunctionBegin; 1867 if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1868 if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1869 PetscCall(MatDenseGetArrayRead(A1,&v1)); 1870 PetscCall(MatDenseGetArrayRead(A2,&v2)); 1871 for (i=0; i<A1->cmap->n; i++) { 1872 PetscCall(PetscArraycmp(v1,v2,A1->rmap->n,flg)); 1873 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 1874 v1 += mat1->lda; 1875 v2 += mat2->lda; 1876 } 1877 PetscCall(MatDenseRestoreArrayRead(A1,&v1)); 1878 PetscCall(MatDenseRestoreArrayRead(A2,&v2)); 1879 *flg = PETSC_TRUE; 1880 PetscFunctionReturn(0); 1881 } 1882 1883 static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1884 { 1885 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1886 PetscInt i,n,len; 1887 PetscScalar *x; 1888 const PetscScalar *vv; 1889 1890 PetscFunctionBegin; 1891 PetscCall(VecGetSize(v,&n)); 1892 PetscCall(VecGetArray(v,&x)); 1893 len = PetscMin(A->rmap->n,A->cmap->n); 1894 PetscCall(MatDenseGetArrayRead(A,&vv)); 1895 PetscCheck(n == A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1896 for (i=0; i<len; i++) { 1897 x[i] = vv[i*mat->lda + i]; 1898 } 1899 PetscCall(MatDenseRestoreArrayRead(A,&vv)); 1900 PetscCall(VecRestoreArray(v,&x)); 1901 PetscFunctionReturn(0); 1902 } 1903 1904 static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1905 { 1906 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1907 const PetscScalar *l,*r; 1908 PetscScalar x,*v,*vv; 1909 PetscInt i,j,m = A->rmap->n,n = A->cmap->n; 1910 1911 PetscFunctionBegin; 1912 PetscCall(MatDenseGetArray(A,&vv)); 1913 if (ll) { 1914 PetscCall(VecGetSize(ll,&m)); 1915 PetscCall(VecGetArrayRead(ll,&l)); 1916 PetscCheck(m == A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1917 for (i=0; i<m; i++) { 1918 x = l[i]; 1919 v = vv + i; 1920 for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;} 1921 } 1922 PetscCall(VecRestoreArrayRead(ll,&l)); 1923 PetscCall(PetscLogFlops(1.0*n*m)); 1924 } 1925 if (rr) { 1926 PetscCall(VecGetSize(rr,&n)); 1927 PetscCall(VecGetArrayRead(rr,&r)); 1928 PetscCheck(n == A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1929 for (i=0; i<n; i++) { 1930 x = r[i]; 1931 v = vv + i*mat->lda; 1932 for (j=0; j<m; j++) (*v++) *= x; 1933 } 1934 PetscCall(VecRestoreArrayRead(rr,&r)); 1935 PetscCall(PetscLogFlops(1.0*n*m)); 1936 } 1937 PetscCall(MatDenseRestoreArray(A,&vv)); 1938 PetscFunctionReturn(0); 1939 } 1940 1941 PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1942 { 1943 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1944 PetscScalar *v,*vv; 1945 PetscReal sum = 0.0; 1946 PetscInt lda, m=A->rmap->n,i,j; 1947 1948 PetscFunctionBegin; 1949 PetscCall(MatDenseGetArrayRead(A,(const PetscScalar**)&vv)); 1950 PetscCall(MatDenseGetLDA(A,&lda)); 1951 v = vv; 1952 if (type == NORM_FROBENIUS) { 1953 if (lda>m) { 1954 for (j=0; j<A->cmap->n; j++) { 1955 v = vv+j*lda; 1956 for (i=0; i<m; i++) { 1957 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1958 } 1959 } 1960 } else { 1961 #if defined(PETSC_USE_REAL___FP16) 1962 PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n; 1963 PetscStackCallBLAS("BLASnrm2",*nrm = BLASnrm2_(&cnt,v,&one)); 1964 } 1965 #else 1966 for (i=0; i<A->cmap->n*A->rmap->n; i++) { 1967 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1968 } 1969 } 1970 *nrm = PetscSqrtReal(sum); 1971 #endif 1972 PetscCall(PetscLogFlops(2.0*A->cmap->n*A->rmap->n)); 1973 } else if (type == NORM_1) { 1974 *nrm = 0.0; 1975 for (j=0; j<A->cmap->n; j++) { 1976 v = vv + j*mat->lda; 1977 sum = 0.0; 1978 for (i=0; i<A->rmap->n; i++) { 1979 sum += PetscAbsScalar(*v); v++; 1980 } 1981 if (sum > *nrm) *nrm = sum; 1982 } 1983 PetscCall(PetscLogFlops(1.0*A->cmap->n*A->rmap->n)); 1984 } else if (type == NORM_INFINITY) { 1985 *nrm = 0.0; 1986 for (j=0; j<A->rmap->n; j++) { 1987 v = vv + j; 1988 sum = 0.0; 1989 for (i=0; i<A->cmap->n; i++) { 1990 sum += PetscAbsScalar(*v); v += mat->lda; 1991 } 1992 if (sum > *nrm) *nrm = sum; 1993 } 1994 PetscCall(PetscLogFlops(1.0*A->cmap->n*A->rmap->n)); 1995 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 1996 PetscCall(MatDenseRestoreArrayRead(A,(const PetscScalar**)&vv)); 1997 PetscFunctionReturn(0); 1998 } 1999 2000 static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg) 2001 { 2002 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 2003 2004 PetscFunctionBegin; 2005 switch (op) { 2006 case MAT_ROW_ORIENTED: 2007 aij->roworiented = flg; 2008 break; 2009 case MAT_NEW_NONZERO_LOCATIONS: 2010 case MAT_NEW_NONZERO_LOCATION_ERR: 2011 case MAT_NEW_NONZERO_ALLOCATION_ERR: 2012 case MAT_FORCE_DIAGONAL_ENTRIES: 2013 case MAT_KEEP_NONZERO_PATTERN: 2014 case MAT_IGNORE_OFF_PROC_ENTRIES: 2015 case MAT_USE_HASH_TABLE: 2016 case MAT_IGNORE_ZERO_ENTRIES: 2017 case MAT_IGNORE_LOWER_TRIANGULAR: 2018 case MAT_SORTED_FULL: 2019 PetscCall(PetscInfo(A,"Option %s ignored\n",MatOptions[op])); 2020 break; 2021 case MAT_SPD: 2022 case MAT_SYMMETRIC: 2023 case MAT_STRUCTURALLY_SYMMETRIC: 2024 case MAT_HERMITIAN: 2025 case MAT_SYMMETRY_ETERNAL: 2026 /* These options are handled directly by MatSetOption() */ 2027 break; 2028 default: 2029 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 2030 } 2031 PetscFunctionReturn(0); 2032 } 2033 2034 PetscErrorCode MatZeroEntries_SeqDense(Mat A) 2035 { 2036 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 2037 PetscInt lda=l->lda,m=A->rmap->n,n=A->cmap->n,j; 2038 PetscScalar *v; 2039 2040 PetscFunctionBegin; 2041 PetscCall(MatDenseGetArrayWrite(A,&v)); 2042 if (lda>m) { 2043 for (j=0; j<n; j++) { 2044 PetscCall(PetscArrayzero(v+j*lda,m)); 2045 } 2046 } else { 2047 PetscCall(PetscArrayzero(v,PetscInt64Mult(m,n))); 2048 } 2049 PetscCall(MatDenseRestoreArrayWrite(A,&v)); 2050 PetscFunctionReturn(0); 2051 } 2052 2053 static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2054 { 2055 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 2056 PetscInt m = l->lda, n = A->cmap->n, i,j; 2057 PetscScalar *slot,*bb,*v; 2058 const PetscScalar *xx; 2059 2060 PetscFunctionBegin; 2061 if (PetscDefined(USE_DEBUG)) { 2062 for (i=0; i<N; i++) { 2063 PetscCheck(rows[i] >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 2064 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); 2065 } 2066 } 2067 if (!N) PetscFunctionReturn(0); 2068 2069 /* fix right hand side if needed */ 2070 if (x && b) { 2071 PetscCall(VecGetArrayRead(x,&xx)); 2072 PetscCall(VecGetArray(b,&bb)); 2073 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]]; 2074 PetscCall(VecRestoreArrayRead(x,&xx)); 2075 PetscCall(VecRestoreArray(b,&bb)); 2076 } 2077 2078 PetscCall(MatDenseGetArray(A,&v)); 2079 for (i=0; i<N; i++) { 2080 slot = v + rows[i]; 2081 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 2082 } 2083 if (diag != 0.0) { 2084 PetscCheck(A->rmap->n == A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 2085 for (i=0; i<N; i++) { 2086 slot = v + (m+1)*rows[i]; 2087 *slot = diag; 2088 } 2089 } 2090 PetscCall(MatDenseRestoreArray(A,&v)); 2091 PetscFunctionReturn(0); 2092 } 2093 2094 static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A,PetscInt *lda) 2095 { 2096 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 2097 2098 PetscFunctionBegin; 2099 *lda = mat->lda; 2100 PetscFunctionReturn(0); 2101 } 2102 2103 PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar **array) 2104 { 2105 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 2106 2107 PetscFunctionBegin; 2108 PetscCheck(!mat->matinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 2109 *array = mat->v; 2110 PetscFunctionReturn(0); 2111 } 2112 2113 PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar **array) 2114 { 2115 PetscFunctionBegin; 2116 if (array) *array = NULL; 2117 PetscFunctionReturn(0); 2118 } 2119 2120 /*@ 2121 MatDenseGetLDA - gets the leading dimension of the array returned from MatDenseGetArray() 2122 2123 Not collective 2124 2125 Input Parameter: 2126 . mat - a MATSEQDENSE or MATMPIDENSE matrix 2127 2128 Output Parameter: 2129 . lda - the leading dimension 2130 2131 Level: intermediate 2132 2133 .seealso: `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseSetLDA()` 2134 @*/ 2135 PetscErrorCode MatDenseGetLDA(Mat A,PetscInt *lda) 2136 { 2137 PetscFunctionBegin; 2138 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 2139 PetscValidIntPointer(lda,2); 2140 MatCheckPreallocated(A,1); 2141 PetscUseMethod(A,"MatDenseGetLDA_C",(Mat,PetscInt*),(A,lda)); 2142 PetscFunctionReturn(0); 2143 } 2144 2145 /*@ 2146 MatDenseSetLDA - Sets the leading dimension of the array used by the dense matrix 2147 2148 Not collective 2149 2150 Input Parameters: 2151 + mat - a MATSEQDENSE or MATMPIDENSE matrix 2152 - lda - the leading dimension 2153 2154 Level: intermediate 2155 2156 .seealso: `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetLDA()` 2157 @*/ 2158 PetscErrorCode MatDenseSetLDA(Mat A,PetscInt lda) 2159 { 2160 PetscFunctionBegin; 2161 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 2162 PetscTryMethod(A,"MatDenseSetLDA_C",(Mat,PetscInt),(A,lda)); 2163 PetscFunctionReturn(0); 2164 } 2165 2166 /*@C 2167 MatDenseGetArray - gives read-write access to the array where the data for a dense matrix is stored 2168 2169 Logically Collective on Mat 2170 2171 Input Parameter: 2172 . mat - a dense matrix 2173 2174 Output Parameter: 2175 . array - pointer to the data 2176 2177 Level: intermediate 2178 2179 .seealso: `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()` 2180 @*/ 2181 PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array) 2182 { 2183 PetscFunctionBegin; 2184 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 2185 PetscValidPointer(array,2); 2186 PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array)); 2187 PetscFunctionReturn(0); 2188 } 2189 2190 /*@C 2191 MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray() 2192 2193 Logically Collective on Mat 2194 2195 Input Parameters: 2196 + mat - a dense matrix 2197 - array - pointer to the data 2198 2199 Level: intermediate 2200 2201 .seealso: `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()` 2202 @*/ 2203 PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array) 2204 { 2205 PetscFunctionBegin; 2206 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 2207 PetscValidPointer(array,2); 2208 PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array)); 2209 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 2210 #if defined(PETSC_HAVE_CUDA) 2211 A->offloadmask = PETSC_OFFLOAD_CPU; 2212 #endif 2213 PetscFunctionReturn(0); 2214 } 2215 2216 /*@C 2217 MatDenseGetArrayRead - gives read-only access to the array where the data for a dense matrix is stored 2218 2219 Not Collective 2220 2221 Input Parameter: 2222 . mat - a dense matrix 2223 2224 Output Parameter: 2225 . array - pointer to the data 2226 2227 Level: intermediate 2228 2229 .seealso: `MatDenseRestoreArrayRead()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()` 2230 @*/ 2231 PetscErrorCode MatDenseGetArrayRead(Mat A,const PetscScalar **array) 2232 { 2233 PetscFunctionBegin; 2234 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 2235 PetscValidPointer(array,2); 2236 PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array)); 2237 PetscFunctionReturn(0); 2238 } 2239 2240 /*@C 2241 MatDenseRestoreArrayRead - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArrayRead() 2242 2243 Not Collective 2244 2245 Input Parameters: 2246 + mat - a dense matrix 2247 - array - pointer to the data 2248 2249 Level: intermediate 2250 2251 .seealso: `MatDenseGetArrayRead()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()` 2252 @*/ 2253 PetscErrorCode MatDenseRestoreArrayRead(Mat A,const PetscScalar **array) 2254 { 2255 PetscFunctionBegin; 2256 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 2257 PetscValidPointer(array,2); 2258 PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array)); 2259 PetscFunctionReturn(0); 2260 } 2261 2262 /*@C 2263 MatDenseGetArrayWrite - gives write-only access to the array where the data for a dense matrix is stored 2264 2265 Not Collective 2266 2267 Input Parameter: 2268 . mat - a dense matrix 2269 2270 Output Parameter: 2271 . array - pointer to the data 2272 2273 Level: intermediate 2274 2275 .seealso: `MatDenseRestoreArrayWrite()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()` 2276 @*/ 2277 PetscErrorCode MatDenseGetArrayWrite(Mat A,PetscScalar **array) 2278 { 2279 PetscFunctionBegin; 2280 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 2281 PetscValidPointer(array,2); 2282 PetscUseMethod(A,"MatDenseGetArrayWrite_C",(Mat,PetscScalar**),(A,array)); 2283 PetscFunctionReturn(0); 2284 } 2285 2286 /*@C 2287 MatDenseRestoreArrayWrite - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArrayWrite() 2288 2289 Not Collective 2290 2291 Input Parameters: 2292 + mat - a dense matrix 2293 - array - pointer to the data 2294 2295 Level: intermediate 2296 2297 .seealso: `MatDenseGetArrayWrite()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()` 2298 @*/ 2299 PetscErrorCode MatDenseRestoreArrayWrite(Mat A,PetscScalar **array) 2300 { 2301 PetscFunctionBegin; 2302 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 2303 PetscValidPointer(array,2); 2304 PetscUseMethod(A,"MatDenseRestoreArrayWrite_C",(Mat,PetscScalar**),(A,array)); 2305 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 2306 #if defined(PETSC_HAVE_CUDA) 2307 A->offloadmask = PETSC_OFFLOAD_CPU; 2308 #endif 2309 PetscFunctionReturn(0); 2310 } 2311 2312 static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,MatReuse scall,Mat *B) 2313 { 2314 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 2315 PetscInt i,j,nrows,ncols,ldb; 2316 const PetscInt *irow,*icol; 2317 PetscScalar *av,*bv,*v = mat->v; 2318 Mat newmat; 2319 2320 PetscFunctionBegin; 2321 PetscCall(ISGetIndices(isrow,&irow)); 2322 PetscCall(ISGetIndices(iscol,&icol)); 2323 PetscCall(ISGetLocalSize(isrow,&nrows)); 2324 PetscCall(ISGetLocalSize(iscol,&ncols)); 2325 2326 /* Check submatrixcall */ 2327 if (scall == MAT_REUSE_MATRIX) { 2328 PetscInt n_cols,n_rows; 2329 PetscCall(MatGetSize(*B,&n_rows,&n_cols)); 2330 if (n_rows != nrows || n_cols != ncols) { 2331 /* resize the result matrix to match number of requested rows/columns */ 2332 PetscCall(MatSetSizes(*B,nrows,ncols,nrows,ncols)); 2333 } 2334 newmat = *B; 2335 } else { 2336 /* Create and fill new matrix */ 2337 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&newmat)); 2338 PetscCall(MatSetSizes(newmat,nrows,ncols,nrows,ncols)); 2339 PetscCall(MatSetType(newmat,((PetscObject)A)->type_name)); 2340 PetscCall(MatSeqDenseSetPreallocation(newmat,NULL)); 2341 } 2342 2343 /* Now extract the data pointers and do the copy,column at a time */ 2344 PetscCall(MatDenseGetArray(newmat,&bv)); 2345 PetscCall(MatDenseGetLDA(newmat,&ldb)); 2346 for (i=0; i<ncols; i++) { 2347 av = v + mat->lda*icol[i]; 2348 for (j=0; j<nrows; j++) bv[j] = av[irow[j]]; 2349 bv += ldb; 2350 } 2351 PetscCall(MatDenseRestoreArray(newmat,&bv)); 2352 2353 /* Assemble the matrices so that the correct flags are set */ 2354 PetscCall(MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY)); 2355 PetscCall(MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY)); 2356 2357 /* Free work space */ 2358 PetscCall(ISRestoreIndices(isrow,&irow)); 2359 PetscCall(ISRestoreIndices(iscol,&icol)); 2360 *B = newmat; 2361 PetscFunctionReturn(0); 2362 } 2363 2364 static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2365 { 2366 PetscInt i; 2367 2368 PetscFunctionBegin; 2369 if (scall == MAT_INITIAL_MATRIX) { 2370 PetscCall(PetscCalloc1(n,B)); 2371 } 2372 2373 for (i=0; i<n; i++) { 2374 PetscCall(MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],scall,&(*B)[i])); 2375 } 2376 PetscFunctionReturn(0); 2377 } 2378 2379 static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 2380 { 2381 PetscFunctionBegin; 2382 PetscFunctionReturn(0); 2383 } 2384 2385 static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 2386 { 2387 PetscFunctionBegin; 2388 PetscFunctionReturn(0); 2389 } 2390 2391 PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 2392 { 2393 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data; 2394 const PetscScalar *va; 2395 PetscScalar *vb; 2396 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 2397 2398 PetscFunctionBegin; 2399 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 2400 if (A->ops->copy != B->ops->copy) { 2401 PetscCall(MatCopy_Basic(A,B,str)); 2402 PetscFunctionReturn(0); 2403 } 2404 PetscCheck(m == B->rmap->n && n == B->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 2405 PetscCall(MatDenseGetArrayRead(A,&va)); 2406 PetscCall(MatDenseGetArray(B,&vb)); 2407 if (lda1>m || lda2>m) { 2408 for (j=0; j<n; j++) { 2409 PetscCall(PetscArraycpy(vb+j*lda2,va+j*lda1,m)); 2410 } 2411 } else { 2412 PetscCall(PetscArraycpy(vb,va,A->rmap->n*A->cmap->n)); 2413 } 2414 PetscCall(MatDenseRestoreArray(B,&vb)); 2415 PetscCall(MatDenseRestoreArrayRead(A,&va)); 2416 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 2417 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 2418 PetscFunctionReturn(0); 2419 } 2420 2421 PetscErrorCode MatSetUp_SeqDense(Mat A) 2422 { 2423 PetscFunctionBegin; 2424 PetscCall(PetscLayoutSetUp(A->rmap)); 2425 PetscCall(PetscLayoutSetUp(A->cmap)); 2426 if (!A->preallocated) { 2427 PetscCall(MatSeqDenseSetPreallocation(A,NULL)); 2428 } 2429 PetscFunctionReturn(0); 2430 } 2431 2432 static PetscErrorCode MatConjugate_SeqDense(Mat A) 2433 { 2434 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 2435 PetscInt i,j; 2436 PetscInt min = PetscMin(A->rmap->n,A->cmap->n); 2437 PetscScalar *aa; 2438 2439 PetscFunctionBegin; 2440 PetscCall(MatDenseGetArray(A,&aa)); 2441 for (j=0; j<A->cmap->n; j++) { 2442 for (i=0; i<A->rmap->n; i++) aa[i+j*mat->lda] = PetscConj(aa[i+j*mat->lda]); 2443 } 2444 PetscCall(MatDenseRestoreArray(A,&aa)); 2445 if (mat->tau) for (i = 0; i < min; i++) mat->tau[i] = PetscConj(mat->tau[i]); 2446 PetscFunctionReturn(0); 2447 } 2448 2449 static PetscErrorCode MatRealPart_SeqDense(Mat A) 2450 { 2451 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 2452 PetscInt i,j; 2453 PetscScalar *aa; 2454 2455 PetscFunctionBegin; 2456 PetscCall(MatDenseGetArray(A,&aa)); 2457 for (j=0; j<A->cmap->n; j++) { 2458 for (i=0; i<A->rmap->n; i++) aa[i+j*mat->lda] = PetscRealPart(aa[i+j*mat->lda]); 2459 } 2460 PetscCall(MatDenseRestoreArray(A,&aa)); 2461 PetscFunctionReturn(0); 2462 } 2463 2464 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 2465 { 2466 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 2467 PetscInt i,j; 2468 PetscScalar *aa; 2469 2470 PetscFunctionBegin; 2471 PetscCall(MatDenseGetArray(A,&aa)); 2472 for (j=0; j<A->cmap->n; j++) { 2473 for (i=0; i<A->rmap->n; i++) aa[i+j*mat->lda] = PetscImaginaryPart(aa[i+j*mat->lda]); 2474 } 2475 PetscCall(MatDenseRestoreArray(A,&aa)); 2476 PetscFunctionReturn(0); 2477 } 2478 2479 /* ----------------------------------------------------------------*/ 2480 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C) 2481 { 2482 PetscInt m=A->rmap->n,n=B->cmap->n; 2483 PetscBool cisdense; 2484 2485 PetscFunctionBegin; 2486 PetscCall(MatSetSizes(C,m,n,m,n)); 2487 PetscCall(PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"")); 2488 if (!cisdense) { 2489 PetscBool flg; 2490 2491 PetscCall(PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg)); 2492 PetscCall(MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE)); 2493 } 2494 PetscCall(MatSetUp(C)); 2495 PetscFunctionReturn(0); 2496 } 2497 2498 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2499 { 2500 Mat_SeqDense *a=(Mat_SeqDense*)A->data,*b=(Mat_SeqDense*)B->data,*c=(Mat_SeqDense*)C->data; 2501 PetscBLASInt m,n,k; 2502 const PetscScalar *av,*bv; 2503 PetscScalar *cv; 2504 PetscScalar _DOne=1.0,_DZero=0.0; 2505 2506 PetscFunctionBegin; 2507 PetscCall(PetscBLASIntCast(C->rmap->n,&m)); 2508 PetscCall(PetscBLASIntCast(C->cmap->n,&n)); 2509 PetscCall(PetscBLASIntCast(A->cmap->n,&k)); 2510 if (!m || !n || !k) PetscFunctionReturn(0); 2511 PetscCall(MatDenseGetArrayRead(A,&av)); 2512 PetscCall(MatDenseGetArrayRead(B,&bv)); 2513 PetscCall(MatDenseGetArrayWrite(C,&cv)); 2514 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda)); 2515 PetscCall(PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1))); 2516 PetscCall(MatDenseRestoreArrayRead(A,&av)); 2517 PetscCall(MatDenseRestoreArrayRead(B,&bv)); 2518 PetscCall(MatDenseRestoreArrayWrite(C,&cv)); 2519 PetscFunctionReturn(0); 2520 } 2521 2522 PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C) 2523 { 2524 PetscInt m=A->rmap->n,n=B->rmap->n; 2525 PetscBool cisdense; 2526 2527 PetscFunctionBegin; 2528 PetscCall(MatSetSizes(C,m,n,m,n)); 2529 PetscCall(PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"")); 2530 if (!cisdense) { 2531 PetscBool flg; 2532 2533 PetscCall(PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg)); 2534 PetscCall(MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE)); 2535 } 2536 PetscCall(MatSetUp(C)); 2537 PetscFunctionReturn(0); 2538 } 2539 2540 PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2541 { 2542 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2543 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2544 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2545 const PetscScalar *av,*bv; 2546 PetscScalar *cv; 2547 PetscBLASInt m,n,k; 2548 PetscScalar _DOne=1.0,_DZero=0.0; 2549 2550 PetscFunctionBegin; 2551 PetscCall(PetscBLASIntCast(C->rmap->n,&m)); 2552 PetscCall(PetscBLASIntCast(C->cmap->n,&n)); 2553 PetscCall(PetscBLASIntCast(A->cmap->n,&k)); 2554 if (!m || !n || !k) PetscFunctionReturn(0); 2555 PetscCall(MatDenseGetArrayRead(A,&av)); 2556 PetscCall(MatDenseGetArrayRead(B,&bv)); 2557 PetscCall(MatDenseGetArrayWrite(C,&cv)); 2558 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda)); 2559 PetscCall(MatDenseRestoreArrayRead(A,&av)); 2560 PetscCall(MatDenseRestoreArrayRead(B,&bv)); 2561 PetscCall(MatDenseRestoreArrayWrite(C,&cv)); 2562 PetscCall(PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1))); 2563 PetscFunctionReturn(0); 2564 } 2565 2566 PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C) 2567 { 2568 PetscInt m=A->cmap->n,n=B->cmap->n; 2569 PetscBool cisdense; 2570 2571 PetscFunctionBegin; 2572 PetscCall(MatSetSizes(C,m,n,m,n)); 2573 PetscCall(PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"")); 2574 if (!cisdense) { 2575 PetscBool flg; 2576 2577 PetscCall(PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg)); 2578 PetscCall(MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE)); 2579 } 2580 PetscCall(MatSetUp(C)); 2581 PetscFunctionReturn(0); 2582 } 2583 2584 PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2585 { 2586 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2587 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2588 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2589 const PetscScalar *av,*bv; 2590 PetscScalar *cv; 2591 PetscBLASInt m,n,k; 2592 PetscScalar _DOne=1.0,_DZero=0.0; 2593 2594 PetscFunctionBegin; 2595 PetscCall(PetscBLASIntCast(C->rmap->n,&m)); 2596 PetscCall(PetscBLASIntCast(C->cmap->n,&n)); 2597 PetscCall(PetscBLASIntCast(A->rmap->n,&k)); 2598 if (!m || !n || !k) PetscFunctionReturn(0); 2599 PetscCall(MatDenseGetArrayRead(A,&av)); 2600 PetscCall(MatDenseGetArrayRead(B,&bv)); 2601 PetscCall(MatDenseGetArrayWrite(C,&cv)); 2602 PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda)); 2603 PetscCall(MatDenseRestoreArrayRead(A,&av)); 2604 PetscCall(MatDenseRestoreArrayRead(B,&bv)); 2605 PetscCall(MatDenseRestoreArrayWrite(C,&cv)); 2606 PetscCall(PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1))); 2607 PetscFunctionReturn(0); 2608 } 2609 2610 /* ----------------------------------------------- */ 2611 static PetscErrorCode MatProductSetFromOptions_SeqDense_AB(Mat C) 2612 { 2613 PetscFunctionBegin; 2614 C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqDense; 2615 C->ops->productsymbolic = MatProductSymbolic_AB; 2616 PetscFunctionReturn(0); 2617 } 2618 2619 static PetscErrorCode MatProductSetFromOptions_SeqDense_AtB(Mat C) 2620 { 2621 PetscFunctionBegin; 2622 C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_SeqDense_SeqDense; 2623 C->ops->productsymbolic = MatProductSymbolic_AtB; 2624 PetscFunctionReturn(0); 2625 } 2626 2627 static PetscErrorCode MatProductSetFromOptions_SeqDense_ABt(Mat C) 2628 { 2629 PetscFunctionBegin; 2630 C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqDense_SeqDense; 2631 C->ops->productsymbolic = MatProductSymbolic_ABt; 2632 PetscFunctionReturn(0); 2633 } 2634 2635 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense(Mat C) 2636 { 2637 Mat_Product *product = C->product; 2638 2639 PetscFunctionBegin; 2640 switch (product->type) { 2641 case MATPRODUCT_AB: 2642 PetscCall(MatProductSetFromOptions_SeqDense_AB(C)); 2643 break; 2644 case MATPRODUCT_AtB: 2645 PetscCall(MatProductSetFromOptions_SeqDense_AtB(C)); 2646 break; 2647 case MATPRODUCT_ABt: 2648 PetscCall(MatProductSetFromOptions_SeqDense_ABt(C)); 2649 break; 2650 default: 2651 break; 2652 } 2653 PetscFunctionReturn(0); 2654 } 2655 /* ----------------------------------------------- */ 2656 2657 static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 2658 { 2659 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2660 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2661 PetscScalar *x; 2662 const PetscScalar *aa; 2663 2664 PetscFunctionBegin; 2665 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2666 PetscCall(VecGetArray(v,&x)); 2667 PetscCall(VecGetLocalSize(v,&p)); 2668 PetscCall(MatDenseGetArrayRead(A,&aa)); 2669 PetscCheck(p == A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2670 for (i=0; i<m; i++) { 2671 x[i] = aa[i]; if (idx) idx[i] = 0; 2672 for (j=1; j<n; j++) { 2673 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;} 2674 } 2675 } 2676 PetscCall(MatDenseRestoreArrayRead(A,&aa)); 2677 PetscCall(VecRestoreArray(v,&x)); 2678 PetscFunctionReturn(0); 2679 } 2680 2681 static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 2682 { 2683 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2684 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2685 PetscScalar *x; 2686 PetscReal atmp; 2687 const PetscScalar *aa; 2688 2689 PetscFunctionBegin; 2690 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2691 PetscCall(VecGetArray(v,&x)); 2692 PetscCall(VecGetLocalSize(v,&p)); 2693 PetscCall(MatDenseGetArrayRead(A,&aa)); 2694 PetscCheck(p == A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2695 for (i=0; i<m; i++) { 2696 x[i] = PetscAbsScalar(aa[i]); 2697 for (j=1; j<n; j++) { 2698 atmp = PetscAbsScalar(aa[i+a->lda*j]); 2699 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 2700 } 2701 } 2702 PetscCall(MatDenseRestoreArrayRead(A,&aa)); 2703 PetscCall(VecRestoreArray(v,&x)); 2704 PetscFunctionReturn(0); 2705 } 2706 2707 static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 2708 { 2709 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2710 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2711 PetscScalar *x; 2712 const PetscScalar *aa; 2713 2714 PetscFunctionBegin; 2715 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2716 PetscCall(MatDenseGetArrayRead(A,&aa)); 2717 PetscCall(VecGetArray(v,&x)); 2718 PetscCall(VecGetLocalSize(v,&p)); 2719 PetscCheck(p == A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2720 for (i=0; i<m; i++) { 2721 x[i] = aa[i]; if (idx) idx[i] = 0; 2722 for (j=1; j<n; j++) { 2723 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;} 2724 } 2725 } 2726 PetscCall(VecRestoreArray(v,&x)); 2727 PetscCall(MatDenseRestoreArrayRead(A,&aa)); 2728 PetscFunctionReturn(0); 2729 } 2730 2731 PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 2732 { 2733 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2734 PetscScalar *x; 2735 const PetscScalar *aa; 2736 2737 PetscFunctionBegin; 2738 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2739 PetscCall(MatDenseGetArrayRead(A,&aa)); 2740 PetscCall(VecGetArray(v,&x)); 2741 PetscCall(PetscArraycpy(x,aa+col*a->lda,A->rmap->n)); 2742 PetscCall(VecRestoreArray(v,&x)); 2743 PetscCall(MatDenseRestoreArrayRead(A,&aa)); 2744 PetscFunctionReturn(0); 2745 } 2746 2747 PETSC_INTERN PetscErrorCode MatGetColumnReductions_SeqDense(Mat A,PetscInt type,PetscReal *reductions) 2748 { 2749 PetscInt i,j,m,n; 2750 const PetscScalar *a; 2751 2752 PetscFunctionBegin; 2753 PetscCall(MatGetSize(A,&m,&n)); 2754 PetscCall(PetscArrayzero(reductions,n)); 2755 PetscCall(MatDenseGetArrayRead(A,&a)); 2756 if (type == NORM_2) { 2757 for (i=0; i<n; i++) { 2758 for (j=0; j<m; j++) { 2759 reductions[i] += PetscAbsScalar(a[j]*a[j]); 2760 } 2761 a += m; 2762 } 2763 } else if (type == NORM_1) { 2764 for (i=0; i<n; i++) { 2765 for (j=0; j<m; j++) { 2766 reductions[i] += PetscAbsScalar(a[j]); 2767 } 2768 a += m; 2769 } 2770 } else if (type == NORM_INFINITY) { 2771 for (i=0; i<n; i++) { 2772 for (j=0; j<m; j++) { 2773 reductions[i] = PetscMax(PetscAbsScalar(a[j]),reductions[i]); 2774 } 2775 a += m; 2776 } 2777 } else if (type == REDUCTION_SUM_REALPART || type == REDUCTION_MEAN_REALPART) { 2778 for (i=0; i<n; i++) { 2779 for (j=0; j<m; j++) { 2780 reductions[i] += PetscRealPart(a[j]); 2781 } 2782 a += m; 2783 } 2784 } else if (type == REDUCTION_SUM_IMAGINARYPART || type == REDUCTION_MEAN_IMAGINARYPART) { 2785 for (i=0; i<n; i++) { 2786 for (j=0; j<m; j++) { 2787 reductions[i] += PetscImaginaryPart(a[j]); 2788 } 2789 a += m; 2790 } 2791 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown reduction type"); 2792 PetscCall(MatDenseRestoreArrayRead(A,&a)); 2793 if (type == NORM_2) { 2794 for (i=0; i<n; i++) reductions[i] = PetscSqrtReal(reductions[i]); 2795 } else if (type == REDUCTION_MEAN_REALPART || type == REDUCTION_MEAN_IMAGINARYPART) { 2796 for (i=0; i<n; i++) reductions[i] /= m; 2797 } 2798 PetscFunctionReturn(0); 2799 } 2800 2801 static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx) 2802 { 2803 PetscScalar *a; 2804 PetscInt lda,m,n,i,j; 2805 2806 PetscFunctionBegin; 2807 PetscCall(MatGetSize(x,&m,&n)); 2808 PetscCall(MatDenseGetLDA(x,&lda)); 2809 PetscCall(MatDenseGetArray(x,&a)); 2810 for (j=0; j<n; j++) { 2811 for (i=0; i<m; i++) { 2812 PetscCall(PetscRandomGetValue(rctx,a+j*lda+i)); 2813 } 2814 } 2815 PetscCall(MatDenseRestoreArray(x,&a)); 2816 PetscFunctionReturn(0); 2817 } 2818 2819 static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool *missing,PetscInt *d) 2820 { 2821 PetscFunctionBegin; 2822 *missing = PETSC_FALSE; 2823 PetscFunctionReturn(0); 2824 } 2825 2826 /* vals is not const */ 2827 static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals) 2828 { 2829 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2830 PetscScalar *v; 2831 2832 PetscFunctionBegin; 2833 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2834 PetscCall(MatDenseGetArray(A,&v)); 2835 *vals = v+col*a->lda; 2836 PetscCall(MatDenseRestoreArray(A,&v)); 2837 PetscFunctionReturn(0); 2838 } 2839 2840 static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals) 2841 { 2842 PetscFunctionBegin; 2843 *vals = NULL; /* user cannot accidentally use the array later */ 2844 PetscFunctionReturn(0); 2845 } 2846 2847 /* -------------------------------------------------------------------*/ 2848 static struct _MatOps MatOps_Values = { MatSetValues_SeqDense, 2849 MatGetRow_SeqDense, 2850 MatRestoreRow_SeqDense, 2851 MatMult_SeqDense, 2852 /* 4*/ MatMultAdd_SeqDense, 2853 MatMultTranspose_SeqDense, 2854 MatMultTransposeAdd_SeqDense, 2855 NULL, 2856 NULL, 2857 NULL, 2858 /* 10*/ NULL, 2859 MatLUFactor_SeqDense, 2860 MatCholeskyFactor_SeqDense, 2861 MatSOR_SeqDense, 2862 MatTranspose_SeqDense, 2863 /* 15*/ MatGetInfo_SeqDense, 2864 MatEqual_SeqDense, 2865 MatGetDiagonal_SeqDense, 2866 MatDiagonalScale_SeqDense, 2867 MatNorm_SeqDense, 2868 /* 20*/ MatAssemblyBegin_SeqDense, 2869 MatAssemblyEnd_SeqDense, 2870 MatSetOption_SeqDense, 2871 MatZeroEntries_SeqDense, 2872 /* 24*/ MatZeroRows_SeqDense, 2873 NULL, 2874 NULL, 2875 NULL, 2876 NULL, 2877 /* 29*/ MatSetUp_SeqDense, 2878 NULL, 2879 NULL, 2880 NULL, 2881 NULL, 2882 /* 34*/ MatDuplicate_SeqDense, 2883 NULL, 2884 NULL, 2885 NULL, 2886 NULL, 2887 /* 39*/ MatAXPY_SeqDense, 2888 MatCreateSubMatrices_SeqDense, 2889 NULL, 2890 MatGetValues_SeqDense, 2891 MatCopy_SeqDense, 2892 /* 44*/ MatGetRowMax_SeqDense, 2893 MatScale_SeqDense, 2894 MatShift_SeqDense, 2895 NULL, 2896 MatZeroRowsColumns_SeqDense, 2897 /* 49*/ MatSetRandom_SeqDense, 2898 NULL, 2899 NULL, 2900 NULL, 2901 NULL, 2902 /* 54*/ NULL, 2903 NULL, 2904 NULL, 2905 NULL, 2906 NULL, 2907 /* 59*/ MatCreateSubMatrix_SeqDense, 2908 MatDestroy_SeqDense, 2909 MatView_SeqDense, 2910 NULL, 2911 NULL, 2912 /* 64*/ NULL, 2913 NULL, 2914 NULL, 2915 NULL, 2916 NULL, 2917 /* 69*/ MatGetRowMaxAbs_SeqDense, 2918 NULL, 2919 NULL, 2920 NULL, 2921 NULL, 2922 /* 74*/ NULL, 2923 NULL, 2924 NULL, 2925 NULL, 2926 NULL, 2927 /* 79*/ NULL, 2928 NULL, 2929 NULL, 2930 NULL, 2931 /* 83*/ MatLoad_SeqDense, 2932 MatIsSymmetric_SeqDense, 2933 MatIsHermitian_SeqDense, 2934 NULL, 2935 NULL, 2936 NULL, 2937 /* 89*/ NULL, 2938 NULL, 2939 MatMatMultNumeric_SeqDense_SeqDense, 2940 NULL, 2941 NULL, 2942 /* 94*/ NULL, 2943 NULL, 2944 NULL, 2945 MatMatTransposeMultNumeric_SeqDense_SeqDense, 2946 NULL, 2947 /* 99*/ MatProductSetFromOptions_SeqDense, 2948 NULL, 2949 NULL, 2950 MatConjugate_SeqDense, 2951 NULL, 2952 /*104*/ NULL, 2953 MatRealPart_SeqDense, 2954 MatImaginaryPart_SeqDense, 2955 NULL, 2956 NULL, 2957 /*109*/ NULL, 2958 NULL, 2959 MatGetRowMin_SeqDense, 2960 MatGetColumnVector_SeqDense, 2961 MatMissingDiagonal_SeqDense, 2962 /*114*/ NULL, 2963 NULL, 2964 NULL, 2965 NULL, 2966 NULL, 2967 /*119*/ NULL, 2968 NULL, 2969 NULL, 2970 NULL, 2971 NULL, 2972 /*124*/ NULL, 2973 MatGetColumnReductions_SeqDense, 2974 NULL, 2975 NULL, 2976 NULL, 2977 /*129*/ NULL, 2978 NULL, 2979 NULL, 2980 MatTransposeMatMultNumeric_SeqDense_SeqDense, 2981 NULL, 2982 /*134*/ NULL, 2983 NULL, 2984 NULL, 2985 NULL, 2986 NULL, 2987 /*139*/ NULL, 2988 NULL, 2989 NULL, 2990 NULL, 2991 NULL, 2992 MatCreateMPIMatConcatenateSeqMat_SeqDense, 2993 /*145*/ NULL, 2994 NULL, 2995 NULL 2996 }; 2997 2998 /*@C 2999 MatCreateSeqDense - Creates a sequential dense matrix that 3000 is stored in column major order (the usual Fortran 77 manner). Many 3001 of the matrix operations use the BLAS and LAPACK routines. 3002 3003 Collective 3004 3005 Input Parameters: 3006 + comm - MPI communicator, set to PETSC_COMM_SELF 3007 . m - number of rows 3008 . n - number of columns 3009 - data - optional location of matrix data in column major order. Set data=NULL for PETSc 3010 to control all matrix memory allocation. 3011 3012 Output Parameter: 3013 . A - the matrix 3014 3015 Notes: 3016 The data input variable is intended primarily for Fortran programmers 3017 who wish to allocate their own matrix memory space. Most users should 3018 set data=NULL. 3019 3020 Level: intermediate 3021 3022 .seealso: `MatCreate()`, `MatCreateDense()`, `MatSetValues()` 3023 @*/ 3024 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 3025 { 3026 PetscFunctionBegin; 3027 PetscCall(MatCreate(comm,A)); 3028 PetscCall(MatSetSizes(*A,m,n,m,n)); 3029 PetscCall(MatSetType(*A,MATSEQDENSE)); 3030 PetscCall(MatSeqDenseSetPreallocation(*A,data)); 3031 PetscFunctionReturn(0); 3032 } 3033 3034 /*@C 3035 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 3036 3037 Collective 3038 3039 Input Parameters: 3040 + B - the matrix 3041 - data - the array (or NULL) 3042 3043 Notes: 3044 The data input variable is intended primarily for Fortran programmers 3045 who wish to allocate their own matrix memory space. Most users should 3046 need not call this routine. 3047 3048 Level: intermediate 3049 3050 .seealso: `MatCreate()`, `MatCreateDense()`, `MatSetValues()`, `MatDenseSetLDA()` 3051 3052 @*/ 3053 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 3054 { 3055 PetscFunctionBegin; 3056 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3057 PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data)); 3058 PetscFunctionReturn(0); 3059 } 3060 3061 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 3062 { 3063 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 3064 3065 PetscFunctionBegin; 3066 PetscCheck(!b->matinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 3067 B->preallocated = PETSC_TRUE; 3068 3069 PetscCall(PetscLayoutSetUp(B->rmap)); 3070 PetscCall(PetscLayoutSetUp(B->cmap)); 3071 3072 if (b->lda <= 0) b->lda = B->rmap->n; 3073 3074 if (!data) { /* petsc-allocated storage */ 3075 if (!b->user_alloc) PetscCall(PetscFree(b->v)); 3076 PetscCall(PetscCalloc1((size_t)b->lda*B->cmap->n,&b->v)); 3077 PetscCall(PetscLogObjectMemory((PetscObject)B,b->lda*B->cmap->n*sizeof(PetscScalar))); 3078 3079 b->user_alloc = PETSC_FALSE; 3080 } else { /* user-allocated storage */ 3081 if (!b->user_alloc) PetscCall(PetscFree(b->v)); 3082 b->v = data; 3083 b->user_alloc = PETSC_TRUE; 3084 } 3085 B->assembled = PETSC_TRUE; 3086 PetscFunctionReturn(0); 3087 } 3088 3089 #if defined(PETSC_HAVE_ELEMENTAL) 3090 PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 3091 { 3092 Mat mat_elemental; 3093 const PetscScalar *array; 3094 PetscScalar *v_colwise; 3095 PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols; 3096 3097 PetscFunctionBegin; 3098 PetscCall(PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols)); 3099 PetscCall(MatDenseGetArrayRead(A,&array)); 3100 /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */ 3101 k = 0; 3102 for (j=0; j<N; j++) { 3103 cols[j] = j; 3104 for (i=0; i<M; i++) { 3105 v_colwise[j*M+i] = array[k++]; 3106 } 3107 } 3108 for (i=0; i<M; i++) { 3109 rows[i] = i; 3110 } 3111 PetscCall(MatDenseRestoreArrayRead(A,&array)); 3112 3113 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental)); 3114 PetscCall(MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N)); 3115 PetscCall(MatSetType(mat_elemental,MATELEMENTAL)); 3116 PetscCall(MatSetUp(mat_elemental)); 3117 3118 /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */ 3119 PetscCall(MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES)); 3120 PetscCall(MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY)); 3121 PetscCall(MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY)); 3122 PetscCall(PetscFree3(v_colwise,rows,cols)); 3123 3124 if (reuse == MAT_INPLACE_MATRIX) { 3125 PetscCall(MatHeaderReplace(A,&mat_elemental)); 3126 } else { 3127 *newmat = mat_elemental; 3128 } 3129 PetscFunctionReturn(0); 3130 } 3131 #endif 3132 3133 PetscErrorCode MatDenseSetLDA_SeqDense(Mat B,PetscInt lda) 3134 { 3135 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 3136 PetscBool data; 3137 3138 PetscFunctionBegin; 3139 data = (PetscBool)((B->rmap->n > 0 && B->cmap->n > 0) ? (b->v ? PETSC_TRUE : PETSC_FALSE) : PETSC_FALSE); 3140 PetscCheck(b->user_alloc || !data || b->lda == lda,PETSC_COMM_SELF,PETSC_ERR_ORDER,"LDA cannot be changed after allocation of internal storage"); 3141 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); 3142 b->lda = lda; 3143 PetscFunctionReturn(0); 3144 } 3145 3146 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 3147 { 3148 PetscMPIInt size; 3149 3150 PetscFunctionBegin; 3151 PetscCallMPI(MPI_Comm_size(comm,&size)); 3152 if (size == 1) { 3153 if (scall == MAT_INITIAL_MATRIX) { 3154 PetscCall(MatDuplicate(inmat,MAT_COPY_VALUES,outmat)); 3155 } else { 3156 PetscCall(MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN)); 3157 } 3158 } else { 3159 PetscCall(MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat)); 3160 } 3161 PetscFunctionReturn(0); 3162 } 3163 3164 PetscErrorCode MatDenseGetColumnVec_SeqDense(Mat A,PetscInt col,Vec *v) 3165 { 3166 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 3167 3168 PetscFunctionBegin; 3169 PetscCheck(!a->vecinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 3170 PetscCheck(!a->matinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 3171 if (!a->cvec) { 3172 PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec)); 3173 PetscCall(PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec)); 3174 } 3175 a->vecinuse = col + 1; 3176 PetscCall(MatDenseGetArray(A,(PetscScalar**)&a->ptrinuse)); 3177 PetscCall(VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda)); 3178 *v = a->cvec; 3179 PetscFunctionReturn(0); 3180 } 3181 3182 PetscErrorCode MatDenseRestoreColumnVec_SeqDense(Mat A,PetscInt col,Vec *v) 3183 { 3184 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 3185 3186 PetscFunctionBegin; 3187 PetscCheck(a->vecinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first"); 3188 PetscCheck(a->cvec,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector"); 3189 a->vecinuse = 0; 3190 PetscCall(MatDenseRestoreArray(A,(PetscScalar**)&a->ptrinuse)); 3191 PetscCall(VecResetArray(a->cvec)); 3192 if (v) *v = NULL; 3193 PetscFunctionReturn(0); 3194 } 3195 3196 PetscErrorCode MatDenseGetColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v) 3197 { 3198 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 3199 3200 PetscFunctionBegin; 3201 PetscCheck(!a->vecinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 3202 PetscCheck(!a->matinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 3203 if (!a->cvec) { 3204 PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec)); 3205 PetscCall(PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec)); 3206 } 3207 a->vecinuse = col + 1; 3208 PetscCall(MatDenseGetArrayRead(A,&a->ptrinuse)); 3209 PetscCall(VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda)); 3210 PetscCall(VecLockReadPush(a->cvec)); 3211 *v = a->cvec; 3212 PetscFunctionReturn(0); 3213 } 3214 3215 PetscErrorCode MatDenseRestoreColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v) 3216 { 3217 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 3218 3219 PetscFunctionBegin; 3220 PetscCheck(a->vecinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first"); 3221 PetscCheck(a->cvec,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector"); 3222 a->vecinuse = 0; 3223 PetscCall(MatDenseRestoreArrayRead(A,&a->ptrinuse)); 3224 PetscCall(VecLockReadPop(a->cvec)); 3225 PetscCall(VecResetArray(a->cvec)); 3226 if (v) *v = NULL; 3227 PetscFunctionReturn(0); 3228 } 3229 3230 PetscErrorCode MatDenseGetColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v) 3231 { 3232 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 3233 3234 PetscFunctionBegin; 3235 PetscCheck(!a->vecinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 3236 PetscCheck(!a->matinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 3237 if (!a->cvec) { 3238 PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec)); 3239 PetscCall(PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec)); 3240 } 3241 a->vecinuse = col + 1; 3242 PetscCall(MatDenseGetArrayWrite(A,(PetscScalar**)&a->ptrinuse)); 3243 PetscCall(VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda)); 3244 *v = a->cvec; 3245 PetscFunctionReturn(0); 3246 } 3247 3248 PetscErrorCode MatDenseRestoreColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v) 3249 { 3250 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 3251 3252 PetscFunctionBegin; 3253 PetscCheck(a->vecinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first"); 3254 PetscCheck(a->cvec,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector"); 3255 a->vecinuse = 0; 3256 PetscCall(MatDenseRestoreArrayWrite(A,(PetscScalar**)&a->ptrinuse)); 3257 PetscCall(VecResetArray(a->cvec)); 3258 if (v) *v = NULL; 3259 PetscFunctionReturn(0); 3260 } 3261 3262 PetscErrorCode MatDenseGetSubMatrix_SeqDense(Mat A,PetscInt rbegin,PetscInt rend,PetscInt cbegin,PetscInt cend,Mat *v) 3263 { 3264 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 3265 3266 PetscFunctionBegin; 3267 PetscCheck(!a->vecinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 3268 PetscCheck(!a->matinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 3269 if (a->cmat && (cend-cbegin != a->cmat->cmap->N || rend-rbegin != a->cmat->rmap->N)) PetscCall(MatDestroy(&a->cmat)); 3270 if (!a->cmat) { 3271 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A),rend-rbegin,PETSC_DECIDE,rend-rbegin,cend-cbegin,a->v+rbegin+(size_t)cbegin*a->lda,&a->cmat)); 3272 PetscCall(PetscLogObjectParent((PetscObject)A,(PetscObject)a->cmat)); 3273 } else { 3274 PetscCall(MatDensePlaceArray(a->cmat,a->v+rbegin+(size_t)cbegin*a->lda)); 3275 } 3276 PetscCall(MatDenseSetLDA(a->cmat,a->lda)); 3277 a->matinuse = cbegin + 1; 3278 *v = a->cmat; 3279 #if defined(PETSC_HAVE_CUDA) 3280 A->offloadmask = PETSC_OFFLOAD_CPU; 3281 #endif 3282 PetscFunctionReturn(0); 3283 } 3284 3285 PetscErrorCode MatDenseRestoreSubMatrix_SeqDense(Mat A,Mat *v) 3286 { 3287 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 3288 3289 PetscFunctionBegin; 3290 PetscCheck(a->matinuse,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetSubMatrix() first"); 3291 PetscCheck(a->cmat,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column matrix"); 3292 PetscCheck(*v == a->cmat,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not the matrix obtained from MatDenseGetSubMatrix()"); 3293 a->matinuse = 0; 3294 PetscCall(MatDenseResetArray(a->cmat)); 3295 *v = NULL; 3296 PetscFunctionReturn(0); 3297 } 3298 3299 /*MC 3300 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 3301 3302 Options Database Keys: 3303 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 3304 3305 Level: beginner 3306 3307 .seealso: `MatCreateSeqDense()` 3308 3309 M*/ 3310 PetscErrorCode MatCreate_SeqDense(Mat B) 3311 { 3312 Mat_SeqDense *b; 3313 PetscMPIInt size; 3314 3315 PetscFunctionBegin; 3316 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B),&size)); 3317 PetscCheck(size <= 1,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 3318 3319 PetscCall(PetscNewLog(B,&b)); 3320 PetscCall(PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps))); 3321 B->data = (void*)b; 3322 3323 b->roworiented = PETSC_TRUE; 3324 3325 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatQRFactor_C",MatQRFactor_SeqDense)); 3326 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseGetLDA_C",MatDenseGetLDA_SeqDense)); 3327 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseSetLDA_C",MatDenseSetLDA_SeqDense)); 3328 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense)); 3329 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense)); 3330 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense)); 3331 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense)); 3332 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseReplaceArray_C",MatDenseReplaceArray_SeqDense)); 3333 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense)); 3334 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense)); 3335 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayWrite_C",MatDenseGetArray_SeqDense)); 3336 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayWrite_C",MatDenseRestoreArray_SeqDense)); 3337 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ)); 3338 #if defined(PETSC_HAVE_ELEMENTAL) 3339 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental)); 3340 #endif 3341 #if defined(PETSC_HAVE_SCALAPACK) 3342 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_scalapack_C",MatConvert_Dense_ScaLAPACK)); 3343 #endif 3344 #if defined(PETSC_HAVE_CUDA) 3345 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqdensecuda_C",MatConvert_SeqDense_SeqDenseCUDA)); 3346 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",MatProductSetFromOptions_SeqDense)); 3347 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdense_C",MatProductSetFromOptions_SeqDense)); 3348 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdensecuda_C",MatProductSetFromOptions_SeqDense)); 3349 #endif 3350 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense)); 3351 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqdense_C",MatProductSetFromOptions_SeqAIJ_SeqDense)); 3352 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdense_C",MatProductSetFromOptions_SeqDense)); 3353 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense)); 3354 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqsbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense)); 3355 3356 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense)); 3357 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense)); 3358 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVec_C",MatDenseGetColumnVec_SeqDense)); 3359 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVec_C",MatDenseRestoreColumnVec_SeqDense)); 3360 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecRead_C",MatDenseGetColumnVecRead_SeqDense)); 3361 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecRead_C",MatDenseRestoreColumnVecRead_SeqDense)); 3362 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecWrite_C",MatDenseGetColumnVecWrite_SeqDense)); 3363 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecWrite_C",MatDenseRestoreColumnVecWrite_SeqDense)); 3364 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseGetSubMatrix_C",MatDenseGetSubMatrix_SeqDense)); 3365 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreSubMatrix_C",MatDenseRestoreSubMatrix_SeqDense)); 3366 PetscCall(PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE)); 3367 PetscFunctionReturn(0); 3368 } 3369 3370 /*@C 3371 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. 3372 3373 Not Collective 3374 3375 Input Parameters: 3376 + mat - a MATSEQDENSE or MATMPIDENSE matrix 3377 - col - column index 3378 3379 Output Parameter: 3380 . vals - pointer to the data 3381 3382 Level: intermediate 3383 3384 .seealso: `MatDenseRestoreColumn()` 3385 @*/ 3386 PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals) 3387 { 3388 PetscFunctionBegin; 3389 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3390 PetscValidLogicalCollectiveInt(A,col,2); 3391 PetscValidPointer(vals,3); 3392 PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals)); 3393 PetscFunctionReturn(0); 3394 } 3395 3396 /*@C 3397 MatDenseRestoreColumn - returns access to a column of a dense matrix which is returned by MatDenseGetColumn(). 3398 3399 Not Collective 3400 3401 Input Parameter: 3402 . mat - a MATSEQDENSE or MATMPIDENSE matrix 3403 3404 Output Parameter: 3405 . vals - pointer to the data 3406 3407 Level: intermediate 3408 3409 .seealso: `MatDenseGetColumn()` 3410 @*/ 3411 PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals) 3412 { 3413 PetscFunctionBegin; 3414 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3415 PetscValidPointer(vals,2); 3416 PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals)); 3417 PetscFunctionReturn(0); 3418 } 3419 3420 /*@ 3421 MatDenseGetColumnVec - Gives read-write access to a column of a dense matrix, represented as a Vec. 3422 3423 Collective 3424 3425 Input Parameters: 3426 + mat - the Mat object 3427 - col - the column index 3428 3429 Output Parameter: 3430 . v - the vector 3431 3432 Notes: 3433 The vector is owned by PETSc. Users need to call MatDenseRestoreColumnVec() when the vector is no longer needed. 3434 Use MatDenseGetColumnVecRead() to obtain read-only access or MatDenseGetColumnVecWrite() for write-only access. 3435 3436 Level: intermediate 3437 3438 .seealso: `MATDENSE`, `MATDENSECUDA`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()` 3439 @*/ 3440 PetscErrorCode MatDenseGetColumnVec(Mat A,PetscInt col,Vec *v) 3441 { 3442 PetscFunctionBegin; 3443 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3444 PetscValidType(A,1); 3445 PetscValidLogicalCollectiveInt(A,col,2); 3446 PetscValidPointer(v,3); 3447 PetscCheck(A->preallocated,PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3448 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); 3449 PetscUseMethod(A,"MatDenseGetColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v)); 3450 PetscFunctionReturn(0); 3451 } 3452 3453 /*@ 3454 MatDenseRestoreColumnVec - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVec(). 3455 3456 Collective 3457 3458 Input Parameters: 3459 + mat - the Mat object 3460 . col - the column index 3461 - v - the Vec object 3462 3463 Level: intermediate 3464 3465 .seealso: `MATDENSE`, `MATDENSECUDA`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()` 3466 @*/ 3467 PetscErrorCode MatDenseRestoreColumnVec(Mat A,PetscInt col,Vec *v) 3468 { 3469 PetscFunctionBegin; 3470 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3471 PetscValidType(A,1); 3472 PetscValidLogicalCollectiveInt(A,col,2); 3473 PetscCheck(A->preallocated,PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3474 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); 3475 PetscUseMethod(A,"MatDenseRestoreColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v)); 3476 PetscFunctionReturn(0); 3477 } 3478 3479 /*@ 3480 MatDenseGetColumnVecRead - Gives read-only access to a column of a dense matrix, represented as a Vec. 3481 3482 Collective 3483 3484 Input Parameters: 3485 + mat - the Mat object 3486 - col - the column index 3487 3488 Output Parameter: 3489 . v - the vector 3490 3491 Notes: 3492 The vector is owned by PETSc and users cannot modify it. 3493 Users need to call MatDenseRestoreColumnVecRead() when the vector is no longer needed. 3494 Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecWrite() for write-only access. 3495 3496 Level: intermediate 3497 3498 .seealso: `MATDENSE`, `MATDENSECUDA`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()` 3499 @*/ 3500 PetscErrorCode MatDenseGetColumnVecRead(Mat A,PetscInt col,Vec *v) 3501 { 3502 PetscFunctionBegin; 3503 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3504 PetscValidType(A,1); 3505 PetscValidLogicalCollectiveInt(A,col,2); 3506 PetscValidPointer(v,3); 3507 PetscCheck(A->preallocated,PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3508 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); 3509 PetscUseMethod(A,"MatDenseGetColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v)); 3510 PetscFunctionReturn(0); 3511 } 3512 3513 /*@ 3514 MatDenseRestoreColumnVecRead - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecRead(). 3515 3516 Collective 3517 3518 Input Parameters: 3519 + mat - the Mat object 3520 . col - the column index 3521 - v - the Vec object 3522 3523 Level: intermediate 3524 3525 .seealso: `MATDENSE`, `MATDENSECUDA`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecWrite()` 3526 @*/ 3527 PetscErrorCode MatDenseRestoreColumnVecRead(Mat A,PetscInt col,Vec *v) 3528 { 3529 PetscFunctionBegin; 3530 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3531 PetscValidType(A,1); 3532 PetscValidLogicalCollectiveInt(A,col,2); 3533 PetscCheck(A->preallocated,PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3534 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); 3535 PetscUseMethod(A,"MatDenseRestoreColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v)); 3536 PetscFunctionReturn(0); 3537 } 3538 3539 /*@ 3540 MatDenseGetColumnVecWrite - Gives write-only access to a column of a dense matrix, represented as a Vec. 3541 3542 Collective 3543 3544 Input Parameters: 3545 + mat - the Mat object 3546 - col - the column index 3547 3548 Output Parameter: 3549 . v - the vector 3550 3551 Notes: 3552 The vector is owned by PETSc. Users need to call MatDenseRestoreColumnVecWrite() when the vector is no longer needed. 3553 Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecRead() for read-only access. 3554 3555 Level: intermediate 3556 3557 .seealso: `MATDENSE`, `MATDENSECUDA`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()` 3558 @*/ 3559 PetscErrorCode MatDenseGetColumnVecWrite(Mat A,PetscInt col,Vec *v) 3560 { 3561 PetscFunctionBegin; 3562 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3563 PetscValidType(A,1); 3564 PetscValidLogicalCollectiveInt(A,col,2); 3565 PetscValidPointer(v,3); 3566 PetscCheck(A->preallocated,PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3567 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); 3568 PetscUseMethod(A,"MatDenseGetColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v)); 3569 PetscFunctionReturn(0); 3570 } 3571 3572 /*@ 3573 MatDenseRestoreColumnVecWrite - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecWrite(). 3574 3575 Collective 3576 3577 Input Parameters: 3578 + mat - the Mat object 3579 . col - the column index 3580 - v - the Vec object 3581 3582 Level: intermediate 3583 3584 .seealso: `MATDENSE`, `MATDENSECUDA`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()` 3585 @*/ 3586 PetscErrorCode MatDenseRestoreColumnVecWrite(Mat A,PetscInt col,Vec *v) 3587 { 3588 PetscFunctionBegin; 3589 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3590 PetscValidType(A,1); 3591 PetscValidLogicalCollectiveInt(A,col,2); 3592 PetscCheck(A->preallocated,PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3593 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); 3594 PetscUseMethod(A,"MatDenseRestoreColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v)); 3595 PetscFunctionReturn(0); 3596 } 3597 3598 /*@ 3599 MatDenseGetSubMatrix - Gives access to a block of rows and columns of a dense matrix, represented as a Mat. 3600 3601 Collective 3602 3603 Input Parameters: 3604 + mat - the Mat object 3605 . rbegin - the first global row index in the block (if PETSC_DECIDE, is 0) 3606 . rend - the global row index past the last one in the block (if PETSC_DECIDE, is M) 3607 . cbegin - the first global column index in the block (if PETSC_DECIDE, is 0) 3608 - cend - the global column index past the last one in the block (if PETSC_DECIDE, is N) 3609 3610 Output Parameter: 3611 . v - the matrix 3612 3613 Notes: 3614 The matrix is owned by PETSc. Users need to call MatDenseRestoreSubMatrix() when the matrix is no longer needed. 3615 The output matrix is not redistributed by PETSc, so depending on the values of rbegin and rend, some processes may have no local rows. 3616 3617 Level: intermediate 3618 3619 .seealso: `MATDENSE`, `MATDENSECUDA`, `MatDenseGetColumnVec()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreSubMatrix()` 3620 @*/ 3621 PetscErrorCode MatDenseGetSubMatrix(Mat A,PetscInt rbegin,PetscInt rend,PetscInt cbegin,PetscInt cend,Mat *v) 3622 { 3623 PetscFunctionBegin; 3624 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3625 PetscValidType(A,1); 3626 PetscValidLogicalCollectiveInt(A,rbegin,2); 3627 PetscValidLogicalCollectiveInt(A,rend,3); 3628 PetscValidLogicalCollectiveInt(A,cbegin,4); 3629 PetscValidLogicalCollectiveInt(A,cend,5); 3630 PetscValidPointer(v,6); 3631 if (rbegin == PETSC_DECIDE) rbegin = 0; 3632 if (rend == PETSC_DECIDE) rend = A->rmap->N; 3633 if (cbegin == PETSC_DECIDE) cbegin = 0; 3634 if (cend == PETSC_DECIDE) cend = A->cmap->N; 3635 PetscCheck(A->preallocated,PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3636 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); 3637 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); 3638 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); 3639 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); 3640 PetscUseMethod(A,"MatDenseGetSubMatrix_C",(Mat,PetscInt,PetscInt,PetscInt,PetscInt,Mat*),(A,rbegin,rend,cbegin,cend,v)); 3641 PetscFunctionReturn(0); 3642 } 3643 3644 /*@ 3645 MatDenseRestoreSubMatrix - Returns access to a block of columns of a dense matrix obtained from MatDenseGetSubMatrix(). 3646 3647 Collective 3648 3649 Input Parameters: 3650 + mat - the Mat object 3651 - v - the Mat object 3652 3653 Level: intermediate 3654 3655 .seealso: `MATDENSE`, `MATDENSECUDA`, `MatDenseGetColumnVec()`, `MatDenseRestoreColumnVec()`, `MatDenseGetSubMatrix()` 3656 @*/ 3657 PetscErrorCode MatDenseRestoreSubMatrix(Mat A,Mat *v) 3658 { 3659 PetscFunctionBegin; 3660 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3661 PetscValidType(A,1); 3662 PetscValidPointer(v,2); 3663 PetscUseMethod(A,"MatDenseRestoreSubMatrix_C",(Mat,Mat*),(A,v)); 3664 PetscFunctionReturn(0); 3665 } 3666 3667 #include <petscblaslapack.h> 3668 #include <petsc/private/kernels/blockinvert.h> 3669 3670 PetscErrorCode MatSeqDenseInvert(Mat A) 3671 { 3672 Mat_SeqDense *a = (Mat_SeqDense*) A->data; 3673 PetscInt bs = A->rmap->n; 3674 MatScalar *values = a->v; 3675 const PetscReal shift = 0.0; 3676 PetscBool allowzeropivot = PetscNot(A->erroriffailure),zeropivotdetected=PETSC_FALSE; 3677 3678 PetscFunctionBegin; 3679 /* factor and invert each block */ 3680 switch (bs) { 3681 case 1: 3682 values[0] = (PetscScalar)1.0 / (values[0] + shift); 3683 break; 3684 case 2: 3685 PetscCall(PetscKernel_A_gets_inverse_A_2(values,shift,allowzeropivot,&zeropivotdetected)); 3686 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3687 break; 3688 case 3: 3689 PetscCall(PetscKernel_A_gets_inverse_A_3(values,shift,allowzeropivot,&zeropivotdetected)); 3690 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3691 break; 3692 case 4: 3693 PetscCall(PetscKernel_A_gets_inverse_A_4(values,shift,allowzeropivot,&zeropivotdetected)); 3694 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3695 break; 3696 case 5: 3697 { 3698 PetscScalar work[25]; 3699 PetscInt ipvt[5]; 3700 3701 PetscCall(PetscKernel_A_gets_inverse_A_5(values,ipvt,work,shift,allowzeropivot,&zeropivotdetected)); 3702 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3703 } 3704 break; 3705 case 6: 3706 PetscCall(PetscKernel_A_gets_inverse_A_6(values,shift,allowzeropivot,&zeropivotdetected)); 3707 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3708 break; 3709 case 7: 3710 PetscCall(PetscKernel_A_gets_inverse_A_7(values,shift,allowzeropivot,&zeropivotdetected)); 3711 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3712 break; 3713 default: 3714 { 3715 PetscInt *v_pivots,*IJ,j; 3716 PetscScalar *v_work; 3717 3718 PetscCall(PetscMalloc3(bs,&v_work,bs,&v_pivots,bs,&IJ)); 3719 for (j=0; j<bs; j++) { 3720 IJ[j] = j; 3721 } 3722 PetscCall(PetscKernel_A_gets_inverse_A(bs,values,v_pivots,v_work,allowzeropivot,&zeropivotdetected)); 3723 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3724 PetscCall(PetscFree3(v_work,v_pivots,IJ)); 3725 } 3726 } 3727 PetscFunctionReturn(0); 3728 } 3729