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