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