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