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