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