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