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