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