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