Actual source code: dense.c

petsc-master 2020-10-23
Report Typos and Errors

  2: /*
  3:      Defines the basic matrix operations for sequential dense.
  4: */

  6: #include <../src/mat/impls/dense/seq/dense.h>
  7: #include <petscblaslapack.h>

  9: #include <../src/mat/impls/aij/seq/aij.h>

 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;

 19:   if (A->rmap->n != A->cmap->n) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot symmetrize a rectangular matrix");
 20:   MatDenseGetArray(A,&v);
 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:   MatDenseRestoreArray(A,&v);
 35:   return(0);
 36: }

 38: PETSC_EXTERN PetscErrorCode MatSeqDenseInvertFactors_Private(Mat A)
 39: {
 40:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
 42:   PetscBLASInt   info,n;

 45:   if (!A->rmap->n || !A->cmap->n) return(0);
 46:   PetscBLASIntCast(A->cmap->n,&n);
 47:   if (A->factortype == MAT_FACTOR_LU) {
 48:     if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present");
 49:     if (!mat->fwork) {
 50:       mat->lfwork = n;
 51:       PetscMalloc1(mat->lfwork,&mat->fwork);
 52:       PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
 53:     }
 54:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
 55:     PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
 56:     PetscFPTrapPop();
 57:     PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
 58:   } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
 59:     if (A->spd) {
 60:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
 61:       PetscStackCallBLAS("LAPACKpotri",LAPACKpotri_("L",&n,mat->v,&mat->lda,&info));
 62:       PetscFPTrapPop();
 63:       MatSeqDenseSymmetrize_Private(A,PETSC_TRUE);
 64: #if defined(PETSC_USE_COMPLEX)
 65:     } else if (A->hermitian) {
 66:       if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present");
 67:       if (!mat->fwork) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Fwork not present");
 68:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
 69:       PetscStackCallBLAS("LAPACKhetri",LAPACKhetri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info));
 70:       PetscFPTrapPop();
 71:       MatSeqDenseSymmetrize_Private(A,PETSC_TRUE);
 72: #endif
 73:     } else { /* symmetric case */
 74:       if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present");
 75:       if (!mat->fwork) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Fwork not present");
 76:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
 77:       PetscStackCallBLAS("LAPACKsytri",LAPACKsytri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info));
 78:       PetscFPTrapPop();
 79:       MatSeqDenseSymmetrize_Private(A,PETSC_FALSE);
 80:     }
 81:     if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad Inversion: zero pivot in row %D",(PetscInt)info-1);
 82:     PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
 83:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");

 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:   PetscFree(A->solvertype);
 93:   return(0);
 94: }

 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;

105:   if (PetscDefined(USE_DEBUG)) {
106:     for (i=0; i<N; i++) {
107:       if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
108:       if (rows[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D requested to be zeroed greater than or equal number of rows %D",rows[i],A->rmap->n);
109:       if (rows[i] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Col %D requested to be zeroed greater than or equal number of cols %D",rows[i],A->cmap->n);
110:     }
111:   }
112:   if (!N) return(0);

114:   /* fix right hand side if needed */
115:   if (x && b) {
116:     Vec xt;

118:     if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
119:     VecDuplicate(x,&xt);
120:     VecCopy(x,xt);
121:     VecScale(xt,-1.0);
122:     MatMultAdd(A,xt,b,b);
123:     VecDestroy(&xt);
124:     VecGetArrayRead(x,&xx);
125:     VecGetArray(b,&bb);
126:     for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
127:     VecRestoreArrayRead(x,&xx);
128:     VecRestoreArray(b,&bb);
129:   }

131:   MatDenseGetArray(A,&v);
132:   for (i=0; i<N; i++) {
133:     slot = v + rows[i]*m;
134:     PetscArrayzero(slot,r);
135:   }
136:   for (i=0; i<N; i++) {
137:     slot = v + rows[i];
138:     for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
139:   }
140:   if (diag != 0.0) {
141:     if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
142:     for (i=0; i<N; i++) {
143:       slot  = v + (m+1)*rows[i];
144:       *slot = diag;
145:     }
146:   }
147:   MatDenseRestoreArray(A,&v);
148:   return(0);
149: }

151: PetscErrorCode MatPtAPNumeric_SeqDense_SeqDense(Mat A,Mat P,Mat C)
152: {
153:   Mat_SeqDense   *c = (Mat_SeqDense*)(C->data);

157:   if (c->ptapwork) {
158:     (*C->ops->matmultnumeric)(A,P,c->ptapwork);
159:     (*C->ops->transposematmultnumeric)(P,c->ptapwork,C);
160:   } else SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"Must call MatPtAPSymbolic_SeqDense_SeqDense() first");
161:   return(0);
162: }

164: PetscErrorCode MatPtAPSymbolic_SeqDense_SeqDense(Mat A,Mat P,PetscReal fill,Mat C)
165: {
166:   Mat_SeqDense   *c;
167:   PetscBool      cisdense;

171:   MatSetSizes(C,P->cmap->n,P->cmap->n,P->cmap->N,P->cmap->N);
172:   PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
173:   if (!cisdense) {
174:     PetscBool flg;

176:     PetscObjectTypeCompare((PetscObject)P,((PetscObject)A)->type_name,&flg);
177:     MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
178:   }
179:   MatSetUp(C);
180:   c    = (Mat_SeqDense*)C->data;
181:   MatCreate(PetscObjectComm((PetscObject)A),&c->ptapwork);
182:   MatSetSizes(c->ptapwork,A->rmap->n,P->cmap->n,A->rmap->N,P->cmap->N);
183:   MatSetType(c->ptapwork,((PetscObject)C)->type_name);
184:   MatSetUp(c->ptapwork);
185:   return(0);
186: }

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;
194:   PetscInt       *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->N,i;
195:   MatScalar      *av=a->a;
196:   PetscBool      isseqdense;

199:   if (reuse == MAT_REUSE_MATRIX) {
200:     PetscObjectTypeCompare((PetscObject)*newmat,MATSEQDENSE,&isseqdense);
201:     if (!isseqdense) SETERRQ1(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_USER,"Cannot reuse matrix of type %s",((PetscObject)(*newmat))->type_name);
202:   }
203:   if (reuse != MAT_REUSE_MATRIX) {
204:     MatCreate(PetscObjectComm((PetscObject)A),&B);
205:     MatSetSizes(B,m,n,m,n);
206:     MatSetType(B,MATSEQDENSE);
207:     MatSeqDenseSetPreallocation(B,NULL);
208:     b    = (Mat_SeqDense*)(B->data);
209:   } else {
210:     b    = (Mat_SeqDense*)((*newmat)->data);
211:     PetscArrayzero(b->v,m*n);
212:   }
213:   for (i=0; i<m; i++) {
214:     PetscInt j;
215:     for (j=0;j<ai[1]-ai[0];j++) {
216:       b->v[*aj*m+i] = *av;
217:       aj++;
218:       av++;
219:     }
220:     ai++;
221:   }

223:   if (reuse == MAT_INPLACE_MATRIX) {
224:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
225:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
226:     MatHeaderReplace(A,&B);
227:   } else {
228:     if (B) *newmat = B;
229:     MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
230:     MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
231:   }
232:   return(0);
233: }

235: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
236: {
237:   Mat            B;
238:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
240:   PetscInt       i, j;
241:   PetscInt       *rows, *nnz;
242:   MatScalar      *aa = a->v, *vals;

245:   MatCreate(PetscObjectComm((PetscObject)A),&B);
246:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
247:   MatSetType(B,MATSEQAIJ);
248:   PetscCalloc3(A->rmap->n,&rows,A->rmap->n,&nnz,A->rmap->n,&vals);
249:   for (j=0; j<A->cmap->n; j++) {
250:     for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) ++nnz[i];
251:     aa += a->lda;
252:   }
253:   MatSeqAIJSetPreallocation(B,PETSC_DETERMINE,nnz);
254:   aa = a->v;
255:   for (j=0; j<A->cmap->n; j++) {
256:     PetscInt numRows = 0;
257:     for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) {rows[numRows] = i; vals[numRows++] = aa[i];}
258:     MatSetValues(B,numRows,rows,1,&j,vals,INSERT_VALUES);
259:     aa  += a->lda;
260:   }
261:   PetscFree3(rows,nnz,vals);
262:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
263:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

265:   if (reuse == MAT_INPLACE_MATRIX) {
266:     MatHeaderReplace(A,&B);
267:   } else {
268:     *newmat = B;
269:   }
270:   return(0);
271: }

273: PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str)
274: {
275:   Mat_SeqDense      *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data;
276:   const PetscScalar *xv;
277:   PetscScalar       *yv;
278:   PetscBLASInt      N,m,ldax = 0,lday = 0,one = 1;
279:   PetscErrorCode    ierr;

282:   MatDenseGetArrayRead(X,&xv);
283:   MatDenseGetArray(Y,&yv);
284:   PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);
285:   PetscBLASIntCast(X->rmap->n,&m);
286:   PetscBLASIntCast(x->lda,&ldax);
287:   PetscBLASIntCast(y->lda,&lday);
288:   if (ldax>m || lday>m) {
289:     PetscInt j;

291:     for (j=0; j<X->cmap->n; j++) {
292:       PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&alpha,xv+j*ldax,&one,yv+j*lday,&one));
293:     }
294:   } else {
295:     PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&alpha,xv,&one,yv,&one));
296:   }
297:   MatDenseRestoreArrayRead(X,&xv);
298:   MatDenseRestoreArray(Y,&yv);
299:   PetscLogFlops(PetscMax(2.0*N-1,0));
300:   return(0);
301: }

303: static PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info)
304: {
305:   PetscLogDouble N = A->rmap->n*A->cmap->n;

308:   info->block_size        = 1.0;
309:   info->nz_allocated      = N;
310:   info->nz_used           = N;
311:   info->nz_unneeded       = 0;
312:   info->assemblies        = A->num_ass;
313:   info->mallocs           = 0;
314:   info->memory            = ((PetscObject)A)->mem;
315:   info->fill_ratio_given  = 0;
316:   info->fill_ratio_needed = 0;
317:   info->factor_mallocs    = 0;
318:   return(0);
319: }

321: PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha)
322: {
323:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
324:   PetscScalar    *v;
326:   PetscBLASInt   one = 1,j,nz,lda = 0;

329:   MatDenseGetArray(A,&v);
330:   PetscBLASIntCast(a->lda,&lda);
331:   if (lda>A->rmap->n) {
332:     PetscBLASIntCast(A->rmap->n,&nz);
333:     for (j=0; j<A->cmap->n; j++) {
334:       PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&alpha,v+j*lda,&one));
335:     }
336:   } else {
337:     PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);
338:     PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&alpha,v,&one));
339:   }
340:   PetscLogFlops(nz);
341:   MatDenseRestoreArray(A,&v);
342:   return(0);
343: }

345: static PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool  *fl)
346: {
347:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
348:   PetscInt          i,j,m = A->rmap->n,N = a->lda;
349:   const PetscScalar *v;
350:   PetscErrorCode    ierr;

353:   *fl = PETSC_FALSE;
354:   if (A->rmap->n != A->cmap->n) return(0);
355:   MatDenseGetArrayRead(A,&v);
356:   for (i=0; i<m; i++) {
357:     for (j=i; j<m; j++) {
358:       if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) {
359:         goto restore;
360:       }
361:     }
362:   }
363:   *fl  = PETSC_TRUE;
364: restore:
365:   MatDenseRestoreArrayRead(A,&v);
366:   return(0);
367: }

369: static PetscErrorCode MatIsSymmetric_SeqDense(Mat A,PetscReal rtol,PetscBool  *fl)
370: {
371:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
372:   PetscInt          i,j,m = A->rmap->n,N = a->lda;
373:   const PetscScalar *v;
374:   PetscErrorCode    ierr;

377:   *fl = PETSC_FALSE;
378:   if (A->rmap->n != A->cmap->n) return(0);
379:   MatDenseGetArrayRead(A,&v);
380:   for (i=0; i<m; i++) {
381:     for (j=i; j<m; j++) {
382:       if (PetscAbsScalar(v[i+j*N] - v[j+i*N]) > rtol) {
383:         goto restore;
384:       }
385:     }
386:   }
387:   *fl  = PETSC_TRUE;
388: restore:
389:   MatDenseRestoreArrayRead(A,&v);
390:   return(0);
391: }

393: PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues)
394: {
395:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
397:   PetscInt       lda = (PetscInt)mat->lda,j,m,nlda = lda;

400:   PetscLayoutReference(A->rmap,&newi->rmap);
401:   PetscLayoutReference(A->cmap,&newi->cmap);
402:   if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { /* propagate LDA */
403:     MatDenseSetLDA(newi,lda);
404:   }
405:   MatSeqDenseSetPreallocation(newi,NULL);
406:   if (cpvalues == MAT_COPY_VALUES) {
407:     const PetscScalar *av;
408:     PetscScalar       *v;

410:     MatDenseGetArrayRead(A,&av);
411:     MatDenseGetArray(newi,&v);
412:     MatDenseGetLDA(newi,&nlda);
413:     m    = A->rmap->n;
414:     if (lda>m || nlda>m) {
415:       for (j=0; j<A->cmap->n; j++) {
416:         PetscArraycpy(v+j*nlda,av+j*lda,m);
417:       }
418:     } else {
419:       PetscArraycpy(v,av,A->rmap->n*A->cmap->n);
420:     }
421:     MatDenseRestoreArray(newi,&v);
422:     MatDenseRestoreArrayRead(A,&av);
423:   }
424:   return(0);
425: }

427: PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat)
428: {

432:   MatCreate(PetscObjectComm((PetscObject)A),newmat);
433:   MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
434:   MatSetType(*newmat,((PetscObject)A)->type_name);
435:   MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);
436:   return(0);
437: }

439: static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
440: {
441:   MatFactorInfo  info;

445:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
446:   (*fact->ops->lufactor)(fact,NULL,NULL,&info);
447:   return(0);
448: }

450: static PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy)
451: {
452:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
453:   PetscErrorCode    ierr;
454:   const PetscScalar *x;
455:   PetscScalar       *y;
456:   PetscBLASInt      one = 1,info,m;

459:   PetscBLASIntCast(A->rmap->n,&m);
460:   VecGetArrayRead(xx,&x);
461:   VecGetArray(yy,&y);
462:   PetscArraycpy(y,x,A->rmap->n);
463:   if (A->factortype == MAT_FACTOR_LU) {
464:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
465:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
466:     PetscFPTrapPop();
467:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
468:   } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
469:     if (A->spd) {
470:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
471:       PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
472:       PetscFPTrapPop();
473:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
474: #if defined(PETSC_USE_COMPLEX)
475:     } else if (A->hermitian) {
476:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
477:       PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
478:       PetscFPTrapPop();
479:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve");
480: #endif
481:     } else { /* symmetric case */
482:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
483:       PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
484:       PetscFPTrapPop();
485:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
486:     }
487:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
488:   VecRestoreArrayRead(xx,&x);
489:   VecRestoreArray(yy,&y);
490:   PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
491:   return(0);
492: }

494: static PetscErrorCode MatMatSolve_SeqDense(Mat A,Mat B,Mat X)
495: {
496:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
497:   PetscErrorCode    ierr;
498:   const PetscScalar *b;
499:   PetscScalar       *x;
500:   PetscInt          n;
501:   PetscBLASInt      nrhs,info,m;

504:   PetscBLASIntCast(A->rmap->n,&m);
505:   MatGetSize(B,NULL,&n);
506:   PetscBLASIntCast(n,&nrhs);
507:   MatDenseGetArrayRead(B,&b);
508:   MatDenseGetArray(X,&x);

510:   PetscArraycpy(x,b,m*nrhs);

512:   if (A->factortype == MAT_FACTOR_LU) {
513:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
514:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
515:     PetscFPTrapPop();
516:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
517:   } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
518:     if (A->spd) {
519:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
520:       PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info));
521:       PetscFPTrapPop();
522:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
523: #if defined(PETSC_USE_COMPLEX)
524:     } else if (A->hermitian) {
525:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
526:       PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
527:       PetscFPTrapPop();
528:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve");
529: #endif
530:     } else { /* symmetric case */
531:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
532:       PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
533:       PetscFPTrapPop();
534:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
535:     }
536:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");

538:   MatDenseRestoreArrayRead(B,&b);
539:   MatDenseRestoreArray(X,&x);
540:   PetscLogFlops(nrhs*(2.0*m*m - m));
541:   return(0);
542: }

544: static PetscErrorCode MatConjugate_SeqDense(Mat);

546: static PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy)
547: {
548:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
549:   PetscErrorCode    ierr;
550:   const PetscScalar *x;
551:   PetscScalar       *y;
552:   PetscBLASInt      one = 1,info,m;

555:   PetscBLASIntCast(A->rmap->n,&m);
556:   VecGetArrayRead(xx,&x);
557:   VecGetArray(yy,&y);
558:   PetscArraycpy(y,x,A->rmap->n);
559:   if (A->factortype == MAT_FACTOR_LU) {
560:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
561:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
562:     PetscFPTrapPop();
563:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve");
564:   } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
565:     if (A->spd) {
566: #if defined(PETSC_USE_COMPLEX)
567:       MatConjugate_SeqDense(A);
568: #endif
569:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
570:       PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
571:       PetscFPTrapPop();
572: #if defined(PETSC_USE_COMPLEX)
573:       MatConjugate_SeqDense(A);
574: #endif
575:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
576: #if defined(PETSC_USE_COMPLEX)
577:     } else if (A->hermitian) {
578:       MatConjugate_SeqDense(A);
579:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
580:       PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
581:       PetscFPTrapPop();
582:       MatConjugate_SeqDense(A);
583: #endif
584:     } else { /* symmetric case */
585:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
586:       PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
587:       PetscFPTrapPop();
588:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
589:     }
590:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
591:   VecRestoreArrayRead(xx,&x);
592:   VecRestoreArray(yy,&y);
593:   PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
594:   return(0);
595: }

597: /* ---------------------------------------------------------------*/
598: /* COMMENT: I have chosen to hide row permutation in the pivots,
599:    rather than put it in the Mat->row slot.*/
600: PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo)
601: {
602:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
604:   PetscBLASInt   n,m,info;

607:   PetscBLASIntCast(A->cmap->n,&n);
608:   PetscBLASIntCast(A->rmap->n,&m);
609:   if (!mat->pivots) {
610:     PetscMalloc1(A->rmap->n,&mat->pivots);
611:     PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
612:   }
613:   if (!A->rmap->n || !A->cmap->n) return(0);
614:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
615:   PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info));
616:   PetscFPTrapPop();

618:   if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization");
619:   if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization");

621:   A->ops->solve             = MatSolve_SeqDense;
622:   A->ops->matsolve          = MatMatSolve_SeqDense;
623:   A->ops->solvetranspose    = MatSolveTranspose_SeqDense;
624:   A->factortype             = MAT_FACTOR_LU;

626:   PetscFree(A->solvertype);
627:   PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);

629:   PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);
630:   return(0);
631: }

633: /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */
634: PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo)
635: {
636:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
638:   PetscBLASInt   info,n;

641:   PetscBLASIntCast(A->cmap->n,&n);
642:   if (!A->rmap->n || !A->cmap->n) return(0);
643:   if (A->spd) {
644:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
645:     PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info));
646:     PetscFPTrapPop();
647: #if defined(PETSC_USE_COMPLEX)
648:   } else if (A->hermitian) {
649:     if (!mat->pivots) {
650:       PetscMalloc1(A->rmap->n,&mat->pivots);
651:       PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
652:     }
653:     if (!mat->fwork) {
654:       PetscScalar dummy;

656:       mat->lfwork = -1;
657:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
658:       PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
659:       PetscFPTrapPop();
660:       mat->lfwork = (PetscInt)PetscRealPart(dummy);
661:       PetscMalloc1(mat->lfwork,&mat->fwork);
662:       PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
663:     }
664:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
665:     PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
666:     PetscFPTrapPop();
667: #endif
668:   } else { /* symmetric case */
669:     if (!mat->pivots) {
670:       PetscMalloc1(A->rmap->n,&mat->pivots);
671:       PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
672:     }
673:     if (!mat->fwork) {
674:       PetscScalar dummy;

676:       mat->lfwork = -1;
677:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
678:       PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
679:       PetscFPTrapPop();
680:       mat->lfwork = (PetscInt)PetscRealPart(dummy);
681:       PetscMalloc1(mat->lfwork,&mat->fwork);
682:       PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
683:     }
684:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
685:     PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
686:     PetscFPTrapPop();
687:   }
688:   if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1);

690:   A->ops->solve             = MatSolve_SeqDense;
691:   A->ops->matsolve          = MatMatSolve_SeqDense;
692:   A->ops->solvetranspose    = MatSolveTranspose_SeqDense;
693:   A->factortype             = MAT_FACTOR_CHOLESKY;

695:   PetscFree(A->solvertype);
696:   PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);

698:   PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
699:   return(0);
700: }

702: PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
703: {
705:   MatFactorInfo  info;

708:   info.fill = 1.0;

710:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
711:   (*fact->ops->choleskyfactor)(fact,NULL,&info);
712:   return(0);
713: }

715: PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
716: {
718:   fact->assembled                  = PETSC_TRUE;
719:   fact->preallocated               = PETSC_TRUE;
720:   fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
721:   fact->ops->solve                 = MatSolve_SeqDense;
722:   fact->ops->matsolve              = MatMatSolve_SeqDense;
723:   fact->ops->solvetranspose        = MatSolveTranspose_SeqDense;
724:   return(0);
725: }

727: PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info)
728: {
730:   fact->preallocated           = PETSC_TRUE;
731:   fact->assembled              = PETSC_TRUE;
732:   fact->ops->lufactornumeric   = MatLUFactorNumeric_SeqDense;
733:   fact->ops->solve             = MatSolve_SeqDense;
734:   fact->ops->matsolve          = MatMatSolve_SeqDense;
735:   fact->ops->solvetranspose    = MatSolveTranspose_SeqDense;
736:   return(0);
737: }

739: /* uses LAPACK */
740: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact)
741: {

745:   MatCreate(PetscObjectComm((PetscObject)A),fact);
746:   MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
747:   MatSetType(*fact,MATDENSE);
748:   if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
749:     (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
750:     (*fact)->ops->ilufactorsymbolic = MatLUFactorSymbolic_SeqDense;
751:   } else {
752:     (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
753:   }
754:   (*fact)->factortype = ftype;

756:   PetscFree((*fact)->solvertype);
757:   PetscStrallocpy(MATSOLVERPETSC,&(*fact)->solvertype);
758:   return(0);
759: }

761: /* ------------------------------------------------------------------*/
762: static PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx)
763: {
764:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
765:   PetscScalar       *x,*v = mat->v,zero = 0.0,xt;
766:   const PetscScalar *b;
767:   PetscErrorCode    ierr;
768:   PetscInt          m = A->rmap->n,i;
769:   PetscBLASInt      o = 1,bm = 0;

772: #if defined(PETSC_HAVE_CUDA)
773:   if (A->offloadmask == PETSC_OFFLOAD_GPU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented");
774: #endif
775:   if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */
776:   PetscBLASIntCast(m,&bm);
777:   if (flag & SOR_ZERO_INITIAL_GUESS) {
778:     /* this is a hack fix, should have another version without the second BLASdotu */
779:     VecSet(xx,zero);
780:   }
781:   VecGetArray(xx,&x);
782:   VecGetArrayRead(bb,&b);
783:   its  = its*lits;
784:   if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
785:   while (its--) {
786:     if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
787:       for (i=0; i<m; i++) {
788:         PetscStackCallBLAS("BLASdotu",xt   = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o));
789:         x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
790:       }
791:     }
792:     if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
793:       for (i=m-1; i>=0; i--) {
794:         PetscStackCallBLAS("BLASdotu",xt   = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o));
795:         x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
796:       }
797:     }
798:   }
799:   VecRestoreArrayRead(bb,&b);
800:   VecRestoreArray(xx,&x);
801:   return(0);
802: }

804: /* -----------------------------------------------------------------*/
805: PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)
806: {
807:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
808:   const PetscScalar *v   = mat->v,*x;
809:   PetscScalar       *y;
810:   PetscErrorCode    ierr;
811:   PetscBLASInt      m, n,_One=1;
812:   PetscScalar       _DOne=1.0,_DZero=0.0;

815:   PetscBLASIntCast(A->rmap->n,&m);
816:   PetscBLASIntCast(A->cmap->n,&n);
817:   VecGetArrayRead(xx,&x);
818:   VecGetArrayWrite(yy,&y);
819:   if (!A->rmap->n || !A->cmap->n) {
820:     PetscBLASInt i;
821:     for (i=0; i<n; i++) y[i] = 0.0;
822:   } else {
823:     PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One));
824:     PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);
825:   }
826:   VecRestoreArrayRead(xx,&x);
827:   VecRestoreArrayWrite(yy,&y);
828:   return(0);
829: }

831: PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy)
832: {
833:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
834:   PetscScalar       *y,_DOne=1.0,_DZero=0.0;
835:   PetscErrorCode    ierr;
836:   PetscBLASInt      m, n, _One=1;
837:   const PetscScalar *v = mat->v,*x;

840:   PetscBLASIntCast(A->rmap->n,&m);
841:   PetscBLASIntCast(A->cmap->n,&n);
842:   VecGetArrayRead(xx,&x);
843:   VecGetArrayWrite(yy,&y);
844:   if (!A->rmap->n || !A->cmap->n) {
845:     PetscBLASInt i;
846:     for (i=0; i<m; i++) y[i] = 0.0;
847:   } else {
848:     PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One));
849:     PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);
850:   }
851:   VecRestoreArrayRead(xx,&x);
852:   VecRestoreArrayWrite(yy,&y);
853:   return(0);
854: }

856: PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
857: {
858:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
859:   const PetscScalar *v = mat->v,*x;
860:   PetscScalar       *y,_DOne=1.0;
861:   PetscErrorCode    ierr;
862:   PetscBLASInt      m, n, _One=1;

865:   PetscBLASIntCast(A->rmap->n,&m);
866:   PetscBLASIntCast(A->cmap->n,&n);
867:   VecCopy(zz,yy);
868:   if (!A->rmap->n || !A->cmap->n) return(0);
869:   VecGetArrayRead(xx,&x);
870:   VecGetArray(yy,&y);
871:   PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
872:   VecRestoreArrayRead(xx,&x);
873:   VecRestoreArray(yy,&y);
874:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
875:   return(0);
876: }

878: PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
879: {
880:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
881:   const PetscScalar *v = mat->v,*x;
882:   PetscScalar       *y;
883:   PetscErrorCode    ierr;
884:   PetscBLASInt      m, n, _One=1;
885:   PetscScalar       _DOne=1.0;

888:   PetscBLASIntCast(A->rmap->n,&m);
889:   PetscBLASIntCast(A->cmap->n,&n);
890:   VecCopy(zz,yy);
891:   if (!A->rmap->n || !A->cmap->n) return(0);
892:   VecGetArrayRead(xx,&x);
893:   VecGetArray(yy,&y);
894:   PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
895:   VecRestoreArrayRead(xx,&x);
896:   VecRestoreArray(yy,&y);
897:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
898:   return(0);
899: }

901: /* -----------------------------------------------------------------*/
902: static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
903: {
904:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
906:   PetscInt       i;

909:   *ncols = A->cmap->n;
910:   if (cols) {
911:     PetscMalloc1(A->cmap->n+1,cols);
912:     for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
913:   }
914:   if (vals) {
915:     const PetscScalar *v;

917:     MatDenseGetArrayRead(A,&v);
918:     PetscMalloc1(A->cmap->n+1,vals);
919:     v   += row;
920:     for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
921:     MatDenseRestoreArrayRead(A,&v);
922:   }
923:   return(0);
924: }

926: static PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
927: {

931:   if (cols) {PetscFree(*cols);}
932:   if (vals) {PetscFree(*vals); }
933:   return(0);
934: }
935: /* ----------------------------------------------------------------*/
936: static PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
937: {
938:   Mat_SeqDense     *mat = (Mat_SeqDense*)A->data;
939:   PetscScalar      *av;
940:   PetscInt         i,j,idx=0;
941: #if defined(PETSC_HAVE_CUDA)
942:   PetscOffloadMask oldf;
943: #endif
944:   PetscErrorCode   ierr;

947:   MatDenseGetArray(A,&av);
948:   if (!mat->roworiented) {
949:     if (addv == INSERT_VALUES) {
950:       for (j=0; j<n; j++) {
951:         if (indexn[j] < 0) {idx += m; continue;}
952:         if (PetscUnlikelyDebug(indexn[j] >= A->cmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1);
953:         for (i=0; i<m; i++) {
954:           if (indexm[i] < 0) {idx++; continue;}
955:           if (PetscUnlikelyDebug(indexm[i] >= A->rmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1);
956:           av[indexn[j]*mat->lda + indexm[i]] = v[idx++];
957:         }
958:       }
959:     } else {
960:       for (j=0; j<n; j++) {
961:         if (indexn[j] < 0) {idx += m; continue;}
962:         if (PetscUnlikelyDebug(indexn[j] >= A->cmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1);
963:         for (i=0; i<m; i++) {
964:           if (indexm[i] < 0) {idx++; continue;}
965:           if (PetscUnlikelyDebug(indexm[i] >= A->rmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1);
966:           av[indexn[j]*mat->lda + indexm[i]] += v[idx++];
967:         }
968:       }
969:     }
970:   } else {
971:     if (addv == INSERT_VALUES) {
972:       for (i=0; i<m; i++) {
973:         if (indexm[i] < 0) { idx += n; continue;}
974:         if (PetscUnlikelyDebug(indexm[i] >= A->rmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1);
975:         for (j=0; j<n; j++) {
976:           if (indexn[j] < 0) { idx++; continue;}
977:           if (PetscUnlikelyDebug(indexn[j] >= A->cmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1);
978:           av[indexn[j]*mat->lda + indexm[i]] = v[idx++];
979:         }
980:       }
981:     } else {
982:       for (i=0; i<m; i++) {
983:         if (indexm[i] < 0) { idx += n; continue;}
984:         if (PetscUnlikelyDebug(indexm[i] >= A->rmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1);
985:         for (j=0; j<n; j++) {
986:           if (indexn[j] < 0) { idx++; continue;}
987:           if (PetscUnlikelyDebug(indexn[j] >= A->cmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1);
988:           av[indexn[j]*mat->lda + indexm[i]] += v[idx++];
989:         }
990:       }
991:     }
992:   }
993:   /* hack to prevent unneeded copy to the GPU while returning the array */
994: #if defined(PETSC_HAVE_CUDA)
995:   oldf = A->offloadmask;
996:   A->offloadmask = PETSC_OFFLOAD_GPU;
997: #endif
998:   MatDenseRestoreArray(A,&av);
999: #if defined(PETSC_HAVE_CUDA)
1000:   A->offloadmask = (oldf == PETSC_OFFLOAD_UNALLOCATED ? PETSC_OFFLOAD_UNALLOCATED : PETSC_OFFLOAD_CPU);
1001: #endif
1002:   return(0);
1003: }

1005: static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
1006: {
1007:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1008:   const PetscScalar *vv;
1009:   PetscInt          i,j;
1010:   PetscErrorCode    ierr;

1013:   MatDenseGetArrayRead(A,&vv);
1014:   /* row-oriented output */
1015:   for (i=0; i<m; i++) {
1016:     if (indexm[i] < 0) {v += n;continue;}
1017:     if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D requested larger than number rows %D",indexm[i],A->rmap->n);
1018:     for (j=0; j<n; j++) {
1019:       if (indexn[j] < 0) {v++; continue;}
1020:       if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column %D requested larger than number columns %D",indexn[j],A->cmap->n);
1021:       *v++ = vv[indexn[j]*mat->lda + indexm[i]];
1022:     }
1023:   }
1024:   MatDenseRestoreArrayRead(A,&vv);
1025:   return(0);
1026: }

1028: /* -----------------------------------------------------------------*/

1030: PetscErrorCode MatView_Dense_Binary(Mat mat,PetscViewer viewer)
1031: {
1032:   PetscErrorCode    ierr;
1033:   PetscBool         skipHeader;
1034:   PetscViewerFormat format;
1035:   PetscInt          header[4],M,N,m,lda,i,j,k;
1036:   const PetscScalar *v;
1037:   PetscScalar       *vwork;

1040:   PetscViewerSetUp(viewer);
1041:   PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
1042:   PetscViewerGetFormat(viewer,&format);
1043:   if (skipHeader) format = PETSC_VIEWER_NATIVE;

1045:   MatGetSize(mat,&M,&N);

1047:   /* write matrix header */
1048:   header[0] = MAT_FILE_CLASSID; header[1] = M; header[2] = N;
1049:   header[3] = (format == PETSC_VIEWER_NATIVE) ? MATRIX_BINARY_FORMAT_DENSE : M*N;
1050:   if (!skipHeader) {PetscViewerBinaryWrite(viewer,header,4,PETSC_INT);}

1052:   MatGetLocalSize(mat,&m,NULL);
1053:   if (format != PETSC_VIEWER_NATIVE) {
1054:     PetscInt nnz = m*N, *iwork;
1055:     /* store row lengths for each row */
1056:     PetscMalloc1(nnz,&iwork);
1057:     for (i=0; i<m; i++) iwork[i] = N;
1058:     PetscViewerBinaryWriteAll(viewer,iwork,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1059:     /* store column indices (zero start index) */
1060:     for (k=0, i=0; i<m; i++)
1061:       for (j=0; j<N; j++, k++)
1062:         iwork[k] = j;
1063:     PetscViewerBinaryWriteAll(viewer,iwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1064:     PetscFree(iwork);
1065:   }
1066:   /* store matrix values as a dense matrix in row major order */
1067:   PetscMalloc1(m*N,&vwork);
1068:   MatDenseGetArrayRead(mat,&v);
1069:   MatDenseGetLDA(mat,&lda);
1070:   for (k=0, i=0; i<m; i++)
1071:     for (j=0; j<N; j++, k++)
1072:       vwork[k] = v[i+lda*j];
1073:   MatDenseRestoreArrayRead(mat,&v);
1074:   PetscViewerBinaryWriteAll(viewer,vwork,m*N,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);
1075:   PetscFree(vwork);
1076:   return(0);
1077: }

1079: PetscErrorCode MatLoad_Dense_Binary(Mat mat,PetscViewer viewer)
1080: {
1082:   PetscBool      skipHeader;
1083:   PetscInt       header[4],M,N,m,nz,lda,i,j,k;
1084:   PetscInt       rows,cols;
1085:   PetscScalar    *v,*vwork;

1088:   PetscViewerSetUp(viewer);
1089:   PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);

1091:   if (!skipHeader) {
1092:     PetscViewerBinaryRead(viewer,header,4,NULL,PETSC_INT);
1093:     if (header[0] != MAT_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Not a matrix object in file");
1094:     M = header[1]; N = header[2];
1095:     if (M < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix row size (%D) in file is negative",M);
1096:     if (N < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix column size (%D) in file is negative",N);
1097:     nz = header[3];
1098:     if (nz != MATRIX_BINARY_FORMAT_DENSE && nz < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Unknown matrix format %D in file",nz);
1099:   } else {
1100:     MatGetSize(mat,&M,&N);
1101:     if (M < 0 || N < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Matrix binary file header was skipped, thus the user must specify the global sizes of input matrix");
1102:     nz = MATRIX_BINARY_FORMAT_DENSE;
1103:   }

1105:   /* setup global sizes if not set */
1106:   if (mat->rmap->N < 0) mat->rmap->N = M;
1107:   if (mat->cmap->N < 0) mat->cmap->N = N;
1108:   MatSetUp(mat);
1109:   /* check if global sizes are correct */
1110:   MatGetSize(mat,&rows,&cols);
1111:   if (M != rows || N != cols) SETERRQ4(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different sizes (%d, %d) than the input matrix (%d, %d)",M,N,rows,cols);

1113:   MatGetSize(mat,NULL,&N);
1114:   MatGetLocalSize(mat,&m,NULL);
1115:   MatDenseGetArray(mat,&v);
1116:   MatDenseGetLDA(mat,&lda);
1117:   if (nz == MATRIX_BINARY_FORMAT_DENSE) {  /* matrix in file is dense format */
1118:     PetscInt nnz = m*N;
1119:     /* read in matrix values */
1120:     PetscMalloc1(nnz,&vwork);
1121:     PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);
1122:     /* store values in column major order */
1123:     for (j=0; j<N; j++)
1124:       for (i=0; i<m; i++)
1125:         v[i+lda*j] = vwork[i*N+j];
1126:     PetscFree(vwork);
1127:   } else { /* matrix in file is sparse format */
1128:     PetscInt nnz = 0, *rlens, *icols;
1129:     /* read in row lengths */
1130:     PetscMalloc1(m,&rlens);
1131:     PetscViewerBinaryReadAll(viewer,rlens,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1132:     for (i=0; i<m; i++) nnz += rlens[i];
1133:     /* read in column indices and values */
1134:     PetscMalloc2(nnz,&icols,nnz,&vwork);
1135:     PetscViewerBinaryReadAll(viewer,icols,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1136:     PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);
1137:     /* store values in column major order */
1138:     for (k=0, i=0; i<m; i++)
1139:       for (j=0; j<rlens[i]; j++, k++)
1140:         v[i+lda*icols[k]] = vwork[k];
1141:     PetscFree(rlens);
1142:     PetscFree2(icols,vwork);
1143:   }
1144:   MatDenseRestoreArray(mat,&v);
1145:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
1146:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
1147:   return(0);
1148: }

1150: PetscErrorCode MatLoad_SeqDense(Mat newMat, PetscViewer viewer)
1151: {
1152:   PetscBool      isbinary, ishdf5;

1158:   /* force binary viewer to load .info file if it has not yet done so */
1159:   PetscViewerSetUp(viewer);
1160:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1161:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,  &ishdf5);
1162:   if (isbinary) {
1163:     MatLoad_Dense_Binary(newMat,viewer);
1164:   } else if (ishdf5) {
1165: #if defined(PETSC_HAVE_HDF5)
1166:     MatLoad_Dense_HDF5(newMat,viewer);
1167: #else
1168:     SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1169: #endif
1170:   } else {
1171:     SETERRQ2(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"Viewer type %s not yet supported for reading %s matrices",((PetscObject)viewer)->type_name,((PetscObject)newMat)->type_name);
1172:   }
1173:   return(0);
1174: }

1176: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
1177: {
1178:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
1179:   PetscErrorCode    ierr;
1180:   PetscInt          i,j;
1181:   const char        *name;
1182:   PetscScalar       *v,*av;
1183:   PetscViewerFormat format;
1184: #if defined(PETSC_USE_COMPLEX)
1185:   PetscBool         allreal = PETSC_TRUE;
1186: #endif

1189:   MatDenseGetArrayRead(A,(const PetscScalar**)&av);
1190:   PetscViewerGetFormat(viewer,&format);
1191:   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1192:     return(0);  /* do nothing for now */
1193:   } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1194:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1195:     for (i=0; i<A->rmap->n; i++) {
1196:       v    = av + i;
1197:       PetscViewerASCIIPrintf(viewer,"row %D:",i);
1198:       for (j=0; j<A->cmap->n; j++) {
1199: #if defined(PETSC_USE_COMPLEX)
1200:         if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
1201:           PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1202:         } else if (PetscRealPart(*v)) {
1203:           PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));
1204:         }
1205: #else
1206:         if (*v) {
1207:           PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);
1208:         }
1209: #endif
1210:         v += a->lda;
1211:       }
1212:       PetscViewerASCIIPrintf(viewer,"\n");
1213:     }
1214:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1215:   } else {
1216:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1217: #if defined(PETSC_USE_COMPLEX)
1218:     /* determine if matrix has all real values */
1219:     v = av;
1220:     for (i=0; i<A->rmap->n*A->cmap->n; i++) {
1221:       if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;}
1222:     }
1223: #endif
1224:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
1225:       PetscObjectGetName((PetscObject)A,&name);
1226:       PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);
1227:       PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);
1228:       PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
1229:     }

1231:     for (i=0; i<A->rmap->n; i++) {
1232:       v = av + i;
1233:       for (j=0; j<A->cmap->n; j++) {
1234: #if defined(PETSC_USE_COMPLEX)
1235:         if (allreal) {
1236:           PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));
1237:         } else {
1238:           PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1239:         }
1240: #else
1241:         PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);
1242: #endif
1243:         v += a->lda;
1244:       }
1245:       PetscViewerASCIIPrintf(viewer,"\n");
1246:     }
1247:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
1248:       PetscViewerASCIIPrintf(viewer,"];\n");
1249:     }
1250:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1251:   }
1252:   MatDenseRestoreArrayRead(A,(const PetscScalar**)&av);
1253:   PetscViewerFlush(viewer);
1254:   return(0);
1255: }

1257: #include <petscdraw.h>
1258: static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa)
1259: {
1260:   Mat               A  = (Mat) Aa;
1261:   PetscErrorCode    ierr;
1262:   PetscInt          m  = A->rmap->n,n = A->cmap->n,i,j;
1263:   int               color = PETSC_DRAW_WHITE;
1264:   const PetscScalar *v;
1265:   PetscViewer       viewer;
1266:   PetscReal         xl,yl,xr,yr,x_l,x_r,y_l,y_r;
1267:   PetscViewerFormat format;

1270:   PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
1271:   PetscViewerGetFormat(viewer,&format);
1272:   PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);

1274:   /* Loop over matrix elements drawing boxes */
1275:   MatDenseGetArrayRead(A,&v);
1276:   if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1277:     PetscDrawCollectiveBegin(draw);
1278:     /* Blue for negative and Red for positive */
1279:     for (j = 0; j < n; j++) {
1280:       x_l = j; x_r = x_l + 1.0;
1281:       for (i = 0; i < m; i++) {
1282:         y_l = m - i - 1.0;
1283:         y_r = y_l + 1.0;
1284:         if (PetscRealPart(v[j*m+i]) >  0.) color = PETSC_DRAW_RED;
1285:         else if (PetscRealPart(v[j*m+i]) <  0.) color = PETSC_DRAW_BLUE;
1286:         else continue;
1287:         PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1288:       }
1289:     }
1290:     PetscDrawCollectiveEnd(draw);
1291:   } else {
1292:     /* use contour shading to indicate magnitude of values */
1293:     /* first determine max of all nonzero values */
1294:     PetscReal minv = 0.0, maxv = 0.0;
1295:     PetscDraw popup;

1297:     for (i=0; i < m*n; i++) {
1298:       if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1299:     }
1300:     if (minv >= maxv) maxv = minv + PETSC_SMALL;
1301:     PetscDrawGetPopup(draw,&popup);
1302:     PetscDrawScalePopup(popup,minv,maxv);

1304:     PetscDrawCollectiveBegin(draw);
1305:     for (j=0; j<n; j++) {
1306:       x_l = j;
1307:       x_r = x_l + 1.0;
1308:       for (i=0; i<m; i++) {
1309:         y_l   = m - i - 1.0;
1310:         y_r   = y_l + 1.0;
1311:         color = PetscDrawRealToColor(PetscAbsScalar(v[j*m+i]),minv,maxv);
1312:         PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1313:       }
1314:     }
1315:     PetscDrawCollectiveEnd(draw);
1316:   }
1317:   MatDenseRestoreArrayRead(A,&v);
1318:   return(0);
1319: }

1321: static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
1322: {
1323:   PetscDraw      draw;
1324:   PetscBool      isnull;
1325:   PetscReal      xr,yr,xl,yl,h,w;

1329:   PetscViewerDrawGetDraw(viewer,0,&draw);
1330:   PetscDrawIsNull(draw,&isnull);
1331:   if (isnull) return(0);

1333:   xr   = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
1334:   xr  += w;          yr += h;        xl = -w;     yl = -h;
1335:   PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
1336:   PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
1337:   PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
1338:   PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
1339:   PetscDrawSave(draw);
1340:   return(0);
1341: }

1343: PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
1344: {
1346:   PetscBool      iascii,isbinary,isdraw;

1349:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1350:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1351:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);

1353:   if (iascii) {
1354:     MatView_SeqDense_ASCII(A,viewer);
1355:   } else if (isbinary) {
1356:     MatView_Dense_Binary(A,viewer);
1357:   } else if (isdraw) {
1358:     MatView_SeqDense_Draw(A,viewer);
1359:   }
1360:   return(0);
1361: }

1363: static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A,const PetscScalar *array)
1364: {
1365:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;

1368:   if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
1369:   if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
1370:   if (a->unplacedarray) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreArray() first");
1371:   a->unplacedarray       = a->v;
1372:   a->unplaced_user_alloc = a->user_alloc;
1373:   a->v                   = (PetscScalar*) array;
1374:   a->user_alloc          = PETSC_TRUE;
1375: #if defined(PETSC_HAVE_CUDA)
1376:   A->offloadmask = PETSC_OFFLOAD_CPU;
1377: #endif
1378:   return(0);
1379: }

1381: static PetscErrorCode MatDenseResetArray_SeqDense(Mat A)
1382: {
1383:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;

1386:   if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
1387:   if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
1388:   a->v             = a->unplacedarray;
1389:   a->user_alloc    = a->unplaced_user_alloc;
1390:   a->unplacedarray = NULL;
1391: #if defined(PETSC_HAVE_CUDA)
1392:   A->offloadmask = PETSC_OFFLOAD_CPU;
1393: #endif
1394:   return(0);
1395: }

1397: static PetscErrorCode MatDenseReplaceArray_SeqDense(Mat A,const PetscScalar *array)
1398: {
1399:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

1403:   if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
1404:   if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
1405:   if (!a->user_alloc) { PetscFree(a->v); }
1406:   a->v           = (PetscScalar*) array;
1407:   a->user_alloc  = PETSC_FALSE;
1408: #if defined(PETSC_HAVE_CUDA)
1409:   A->offloadmask = PETSC_OFFLOAD_CPU;
1410: #endif
1411:   return(0);
1412: }

1414: PetscErrorCode MatDestroy_SeqDense(Mat mat)
1415: {
1416:   Mat_SeqDense   *l = (Mat_SeqDense*)mat->data;

1420: #if defined(PETSC_USE_LOG)
1421:   PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n);
1422: #endif
1423:   PetscFree(l->pivots);
1424:   PetscFree(l->fwork);
1425:   MatDestroy(&l->ptapwork);
1426:   if (!l->user_alloc) {PetscFree(l->v);}
1427:   if (!l->unplaced_user_alloc) {PetscFree(l->unplacedarray);}
1428:   if (l->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
1429:   if (l->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
1430:   VecDestroy(&l->cvec);
1431:   MatDestroy(&l->cmat);
1432:   PetscFree(mat->data);

1434:   PetscObjectChangeTypeName((PetscObject)mat,NULL);
1435:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetLDA_C",NULL);
1436:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseSetLDA_C",NULL);
1437:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);
1438:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);
1439:   PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL);
1440:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL);
1441:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseReplaceArray_C",NULL);
1442:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayRead_C",NULL);
1443:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayRead_C",NULL);
1444:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayWrite_C",NULL);
1445:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayWrite_C",NULL);
1446:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);
1447: #if defined(PETSC_HAVE_ELEMENTAL)
1448:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);
1449: #endif
1450: #if defined(PETSC_HAVE_SCALAPACK)
1451:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_scalapack_C",NULL);
1452: #endif
1453: #if defined(PETSC_HAVE_CUDA)
1454:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqdensecuda_C",NULL);
1455:   PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",NULL);
1456:   PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdense_C",NULL);
1457: #endif
1458:   PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);
1459:   PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqaij_seqdense_C",NULL);
1460:   PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdense_seqdense_C",NULL);
1461:   PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqbaij_seqdense_C",NULL);
1462:   PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqsbaij_seqdense_C",NULL);

1464:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumn_C",NULL);
1465:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumn_C",NULL);
1466:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVec_C",NULL);
1467:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVec_C",NULL);
1468:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVecRead_C",NULL);
1469:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVecRead_C",NULL);
1470:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVecWrite_C",NULL);
1471:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVecWrite_C",NULL);
1472:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetSubMatrix_C",NULL);
1473:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreSubMatrix_C",NULL);
1474:   return(0);
1475: }

1477: static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
1478: {
1479:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1481:   PetscInt       k,j,m = A->rmap->n, M = mat->lda, n = A->cmap->n;
1482:   PetscScalar    *v,tmp;

1485:   if (reuse == MAT_INPLACE_MATRIX) {
1486:     if (m == n) { /* in place transpose */
1487:       MatDenseGetArray(A,&v);
1488:       for (j=0; j<m; j++) {
1489:         for (k=0; k<j; k++) {
1490:           tmp        = v[j + k*M];
1491:           v[j + k*M] = v[k + j*M];
1492:           v[k + j*M] = tmp;
1493:         }
1494:       }
1495:       MatDenseRestoreArray(A,&v);
1496:     } else { /* reuse memory, temporary allocates new memory */
1497:       PetscScalar *v2;
1498:       PetscLayout tmplayout;

1500:       PetscMalloc1((size_t)m*n,&v2);
1501:       MatDenseGetArray(A,&v);
1502:       for (j=0; j<n; j++) {
1503:         for (k=0; k<m; k++) v2[j + (size_t)k*n] = v[k + (size_t)j*M];
1504:       }
1505:       PetscArraycpy(v,v2,(size_t)m*n);
1506:       PetscFree(v2);
1507:       MatDenseRestoreArray(A,&v);
1508:       /* cleanup size dependent quantities */
1509:       VecDestroy(&mat->cvec);
1510:       MatDestroy(&mat->cmat);
1511:       PetscFree(mat->pivots);
1512:       PetscFree(mat->fwork);
1513:       MatDestroy(&mat->ptapwork);
1514:       /* swap row/col layouts */
1515:       mat->lda  = n;
1516:       tmplayout = A->rmap;
1517:       A->rmap   = A->cmap;
1518:       A->cmap   = tmplayout;
1519:     }
1520:   } else { /* out-of-place transpose */
1521:     Mat          tmat;
1522:     Mat_SeqDense *tmatd;
1523:     PetscScalar  *v2;
1524:     PetscInt     M2;

1526:     if (reuse == MAT_INITIAL_MATRIX) {
1527:       MatCreate(PetscObjectComm((PetscObject)A),&tmat);
1528:       MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);
1529:       MatSetType(tmat,((PetscObject)A)->type_name);
1530:       MatSeqDenseSetPreallocation(tmat,NULL);
1531:     } else tmat = *matout;

1533:     MatDenseGetArrayRead(A,(const PetscScalar**)&v);
1534:     MatDenseGetArray(tmat,&v2);
1535:     tmatd = (Mat_SeqDense*)tmat->data;
1536:     M2    = tmatd->lda;
1537:     for (j=0; j<n; j++) {
1538:       for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M];
1539:     }
1540:     MatDenseRestoreArray(tmat,&v2);
1541:     MatDenseRestoreArrayRead(A,(const PetscScalar**)&v);
1542:     MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
1543:     MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
1544:     *matout = tmat;
1545:   }
1546:   return(0);
1547: }

1549: static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool  *flg)
1550: {
1551:   Mat_SeqDense      *mat1 = (Mat_SeqDense*)A1->data;
1552:   Mat_SeqDense      *mat2 = (Mat_SeqDense*)A2->data;
1553:   PetscInt          i;
1554:   const PetscScalar *v1,*v2;
1555:   PetscErrorCode    ierr;

1558:   if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; return(0);}
1559:   if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; return(0);}
1560:   MatDenseGetArrayRead(A1,&v1);
1561:   MatDenseGetArrayRead(A2,&v2);
1562:   for (i=0; i<A1->cmap->n; i++) {
1563:     PetscArraycmp(v1,v2,A1->rmap->n,flg);
1564:     if (*flg == PETSC_FALSE) return(0);
1565:     v1 += mat1->lda;
1566:     v2 += mat2->lda;
1567:   }
1568:   MatDenseRestoreArrayRead(A1,&v1);
1569:   MatDenseRestoreArrayRead(A2,&v2);
1570:   *flg = PETSC_TRUE;
1571:   return(0);
1572: }

1574: static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
1575: {
1576:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1577:   PetscInt          i,n,len;
1578:   PetscScalar       *x;
1579:   const PetscScalar *vv;
1580:   PetscErrorCode    ierr;

1583:   VecGetSize(v,&n);
1584:   VecGetArray(v,&x);
1585:   len  = PetscMin(A->rmap->n,A->cmap->n);
1586:   MatDenseGetArrayRead(A,&vv);
1587:   if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
1588:   for (i=0; i<len; i++) {
1589:     x[i] = vv[i*mat->lda + i];
1590:   }
1591:   MatDenseRestoreArrayRead(A,&vv);
1592:   VecRestoreArray(v,&x);
1593:   return(0);
1594: }

1596: static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
1597: {
1598:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1599:   const PetscScalar *l,*r;
1600:   PetscScalar       x,*v,*vv;
1601:   PetscErrorCode    ierr;
1602:   PetscInt          i,j,m = A->rmap->n,n = A->cmap->n;

1605:   MatDenseGetArray(A,&vv);
1606:   if (ll) {
1607:     VecGetSize(ll,&m);
1608:     VecGetArrayRead(ll,&l);
1609:     if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size");
1610:     for (i=0; i<m; i++) {
1611:       x = l[i];
1612:       v = vv + i;
1613:       for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;}
1614:     }
1615:     VecRestoreArrayRead(ll,&l);
1616:     PetscLogFlops(1.0*n*m);
1617:   }
1618:   if (rr) {
1619:     VecGetSize(rr,&n);
1620:     VecGetArrayRead(rr,&r);
1621:     if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size");
1622:     for (i=0; i<n; i++) {
1623:       x = r[i];
1624:       v = vv + i*mat->lda;
1625:       for (j=0; j<m; j++) (*v++) *= x;
1626:     }
1627:     VecRestoreArrayRead(rr,&r);
1628:     PetscLogFlops(1.0*n*m);
1629:   }
1630:   MatDenseRestoreArray(A,&vv);
1631:   return(0);
1632: }

1634: PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1635: {
1636:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1637:   PetscScalar       *v,*vv;
1638:   PetscReal         sum  = 0.0;
1639:   PetscInt          lda  =mat->lda,m=A->rmap->n,i,j;
1640:   PetscErrorCode    ierr;

1643:   MatDenseGetArrayRead(A,(const PetscScalar**)&vv);
1644:   v    = vv;
1645:   if (type == NORM_FROBENIUS) {
1646:     if (lda>m) {
1647:       for (j=0; j<A->cmap->n; j++) {
1648:         v = vv+j*lda;
1649:         for (i=0; i<m; i++) {
1650:           sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1651:         }
1652:       }
1653:     } else {
1654: #if defined(PETSC_USE_REAL___FP16)
1655:       PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n;
1656:       *nrm = BLASnrm2_(&cnt,v,&one);
1657:     }
1658: #else
1659:       for (i=0; i<A->cmap->n*A->rmap->n; i++) {
1660:         sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1661:       }
1662:     }
1663:     *nrm = PetscSqrtReal(sum);
1664: #endif
1665:     PetscLogFlops(2.0*A->cmap->n*A->rmap->n);
1666:   } else if (type == NORM_1) {
1667:     *nrm = 0.0;
1668:     for (j=0; j<A->cmap->n; j++) {
1669:       v   = vv + j*mat->lda;
1670:       sum = 0.0;
1671:       for (i=0; i<A->rmap->n; i++) {
1672:         sum += PetscAbsScalar(*v);  v++;
1673:       }
1674:       if (sum > *nrm) *nrm = sum;
1675:     }
1676:     PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1677:   } else if (type == NORM_INFINITY) {
1678:     *nrm = 0.0;
1679:     for (j=0; j<A->rmap->n; j++) {
1680:       v   = vv + j;
1681:       sum = 0.0;
1682:       for (i=0; i<A->cmap->n; i++) {
1683:         sum += PetscAbsScalar(*v); v += mat->lda;
1684:       }
1685:       if (sum > *nrm) *nrm = sum;
1686:     }
1687:     PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1688:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
1689:   MatDenseRestoreArrayRead(A,(const PetscScalar**)&vv);
1690:   return(0);
1691: }

1693: static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
1694: {
1695:   Mat_SeqDense   *aij = (Mat_SeqDense*)A->data;

1699:   switch (op) {
1700:   case MAT_ROW_ORIENTED:
1701:     aij->roworiented = flg;
1702:     break;
1703:   case MAT_NEW_NONZERO_LOCATIONS:
1704:   case MAT_NEW_NONZERO_LOCATION_ERR:
1705:   case MAT_NEW_NONZERO_ALLOCATION_ERR:
1706:   case MAT_NEW_DIAGONALS:
1707:   case MAT_KEEP_NONZERO_PATTERN:
1708:   case MAT_IGNORE_OFF_PROC_ENTRIES:
1709:   case MAT_USE_HASH_TABLE:
1710:   case MAT_IGNORE_ZERO_ENTRIES:
1711:   case MAT_IGNORE_LOWER_TRIANGULAR:
1712:   case MAT_SORTED_FULL:
1713:     PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
1714:     break;
1715:   case MAT_SPD:
1716:   case MAT_SYMMETRIC:
1717:   case MAT_STRUCTURALLY_SYMMETRIC:
1718:   case MAT_HERMITIAN:
1719:   case MAT_SYMMETRY_ETERNAL:
1720:     /* These options are handled directly by MatSetOption() */
1721:     break;
1722:   default:
1723:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]);
1724:   }
1725:   return(0);
1726: }

1728: static PetscErrorCode MatZeroEntries_SeqDense(Mat A)
1729: {
1730:   Mat_SeqDense   *l = (Mat_SeqDense*)A->data;
1732:   PetscInt       lda=l->lda,m=A->rmap->n,j;
1733:   PetscScalar    *v;

1736:   MatDenseGetArray(A,&v);
1737:   if (lda>m) {
1738:     for (j=0; j<A->cmap->n; j++) {
1739:       PetscArrayzero(v+j*lda,m);
1740:     }
1741:   } else {
1742:     PetscArrayzero(v,A->rmap->n*A->cmap->n);
1743:   }
1744:   MatDenseRestoreArray(A,&v);
1745:   return(0);
1746: }

1748: static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1749: {
1750:   PetscErrorCode    ierr;
1751:   Mat_SeqDense      *l = (Mat_SeqDense*)A->data;
1752:   PetscInt          m  = l->lda, n = A->cmap->n, i,j;
1753:   PetscScalar       *slot,*bb,*v;
1754:   const PetscScalar *xx;

1757:   if (PetscDefined(USE_DEBUG)) {
1758:     for (i=0; i<N; i++) {
1759:       if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
1760:       if (rows[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D requested to be zeroed greater than or equal number of rows %D",rows[i],A->rmap->n);
1761:     }
1762:   }
1763:   if (!N) return(0);

1765:   /* fix right hand side if needed */
1766:   if (x && b) {
1767:     VecGetArrayRead(x,&xx);
1768:     VecGetArray(b,&bb);
1769:     for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
1770:     VecRestoreArrayRead(x,&xx);
1771:     VecRestoreArray(b,&bb);
1772:   }

1774:   MatDenseGetArray(A,&v);
1775:   for (i=0; i<N; i++) {
1776:     slot = v + rows[i];
1777:     for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
1778:   }
1779:   if (diag != 0.0) {
1780:     if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
1781:     for (i=0; i<N; i++) {
1782:       slot  = v + (m+1)*rows[i];
1783:       *slot = diag;
1784:     }
1785:   }
1786:   MatDenseRestoreArray(A,&v);
1787:   return(0);
1788: }

1790: static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A,PetscInt *lda)
1791: {
1792:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;

1795:   *lda = mat->lda;
1796:   return(0);
1797: }

1799: PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar **array)
1800: {
1801:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;

1804:   if (mat->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
1805:   *array = mat->v;
1806:   return(0);
1807: }

1809: PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar **array)
1810: {
1812:   *array = NULL;
1813:   return(0);
1814: }

1816: /*@C
1817:    MatDenseGetLDA - gets the leading dimension of the array returned from MatDenseGetArray()

1819:    Not collective

1821:    Input Parameter:
1822: .  mat - a MATSEQDENSE or MATMPIDENSE matrix

1824:    Output Parameter:
1825: .   lda - the leading dimension

1827:    Level: intermediate

1829: .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseSetLDA()
1830: @*/
1831: PetscErrorCode  MatDenseGetLDA(Mat A,PetscInt *lda)
1832: {

1838:   PetscUseMethod(A,"MatDenseGetLDA_C",(Mat,PetscInt*),(A,lda));
1839:   return(0);
1840: }

1842: /*@C
1843:    MatDenseSetLDA - Sets the leading dimension of the array used by the dense matrix

1845:    Not collective

1847:    Input Parameter:
1848: +  mat - a MATSEQDENSE or MATMPIDENSE matrix
1849: -  lda - the leading dimension

1851:    Level: intermediate

1853: .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetLDA()
1854: @*/
1855: PetscErrorCode  MatDenseSetLDA(Mat A,PetscInt lda)
1856: {

1861:   PetscTryMethod(A,"MatDenseSetLDA_C",(Mat,PetscInt),(A,lda));
1862:   return(0);
1863: }

1865: /*@C
1866:    MatDenseGetArray - gives read-write access to the array where the data for a dense matrix is stored

1868:    Logically Collective on Mat

1870:    Input Parameter:
1871: .  mat - a dense matrix

1873:    Output Parameter:
1874: .   array - pointer to the data

1876:    Level: intermediate

1878: .seealso: MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
1879: @*/
1880: PetscErrorCode  MatDenseGetArray(Mat A,PetscScalar **array)
1881: {

1887:   PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));
1888:   return(0);
1889: }

1891: /*@C
1892:    MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray()

1894:    Logically Collective on Mat

1896:    Input Parameters:
1897: +  mat - a dense matrix
1898: -  array - pointer to the data

1900:    Level: intermediate

1902: .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
1903: @*/
1904: PetscErrorCode  MatDenseRestoreArray(Mat A,PetscScalar **array)
1905: {

1911:   PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));
1912:   PetscObjectStateIncrease((PetscObject)A);
1913: #if defined(PETSC_HAVE_CUDA)
1914:   A->offloadmask = PETSC_OFFLOAD_CPU;
1915: #endif
1916:   return(0);
1917: }

1919: /*@C
1920:    MatDenseGetArrayRead - gives read-only access to the array where the data for a dense matrix is stored

1922:    Not Collective

1924:    Input Parameter:
1925: .  mat - a dense matrix

1927:    Output Parameter:
1928: .   array - pointer to the data

1930:    Level: intermediate

1932: .seealso: MatDenseRestoreArrayRead(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
1933: @*/
1934: PetscErrorCode  MatDenseGetArrayRead(Mat A,const PetscScalar **array)
1935: {

1941:   PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array));
1942:   return(0);
1943: }

1945: /*@C
1946:    MatDenseRestoreArrayRead - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArrayRead()

1948:    Not Collective

1950:    Input Parameters:
1951: +  mat - a dense matrix
1952: -  array - pointer to the data

1954:    Level: intermediate

1956: .seealso: MatDenseGetArrayRead(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
1957: @*/
1958: PetscErrorCode  MatDenseRestoreArrayRead(Mat A,const PetscScalar **array)
1959: {

1965:   PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array));
1966:   return(0);
1967: }

1969: /*@C
1970:    MatDenseGetArrayWrite - gives write-only access to the array where the data for a dense matrix is stored

1972:    Not Collective

1974:    Input Parameter:
1975: .  mat - a dense matrix

1977:    Output Parameter:
1978: .   array - pointer to the data

1980:    Level: intermediate

1982: .seealso: MatDenseRestoreArrayWrite(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead()
1983: @*/
1984: PetscErrorCode  MatDenseGetArrayWrite(Mat A,PetscScalar **array)
1985: {

1991:   PetscUseMethod(A,"MatDenseGetArrayWrite_C",(Mat,PetscScalar**),(A,array));
1992:   return(0);
1993: }

1995: /*@C
1996:    MatDenseRestoreArrayWrite - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArrayWrite()

1998:    Not Collective

2000:    Input Parameters:
2001: +  mat - a dense matrix
2002: -  array - pointer to the data

2004:    Level: intermediate

2006: .seealso: MatDenseGetArrayWrite(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead()
2007: @*/
2008: PetscErrorCode  MatDenseRestoreArrayWrite(Mat A,PetscScalar **array)
2009: {

2015:   PetscUseMethod(A,"MatDenseRestoreArrayWrite_C",(Mat,PetscScalar**),(A,array));
2016:   PetscObjectStateIncrease((PetscObject)A);
2017: #if defined(PETSC_HAVE_CUDA)
2018:   A->offloadmask = PETSC_OFFLOAD_CPU;
2019: #endif
2020:   return(0);
2021: }

2023: static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B)
2024: {
2025:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
2027:   PetscInt       i,j,nrows,ncols,blda;
2028:   const PetscInt *irow,*icol;
2029:   PetscScalar    *av,*bv,*v = mat->v;
2030:   Mat            newmat;

2033:   ISGetIndices(isrow,&irow);
2034:   ISGetIndices(iscol,&icol);
2035:   ISGetLocalSize(isrow,&nrows);
2036:   ISGetLocalSize(iscol,&ncols);

2038:   /* Check submatrixcall */
2039:   if (scall == MAT_REUSE_MATRIX) {
2040:     PetscInt n_cols,n_rows;
2041:     MatGetSize(*B,&n_rows,&n_cols);
2042:     if (n_rows != nrows || n_cols != ncols) {
2043:       /* resize the result matrix to match number of requested rows/columns */
2044:       MatSetSizes(*B,nrows,ncols,nrows,ncols);
2045:     }
2046:     newmat = *B;
2047:   } else {
2048:     /* Create and fill new matrix */
2049:     MatCreate(PetscObjectComm((PetscObject)A),&newmat);
2050:     MatSetSizes(newmat,nrows,ncols,nrows,ncols);
2051:     MatSetType(newmat,((PetscObject)A)->type_name);
2052:     MatSeqDenseSetPreallocation(newmat,NULL);
2053:   }

2055:   /* Now extract the data pointers and do the copy,column at a time */
2056:   MatDenseGetArray(newmat,&bv);
2057:   MatDenseGetLDA(newmat,&blda);
2058:   for (i=0; i<ncols; i++) {
2059:     av = v + mat->lda*icol[i];
2060:     for (j=0; j<nrows; j++) bv[j] = av[irow[j]];
2061:     bv += blda;
2062:   }
2063:   MatDenseRestoreArray(newmat,&bv);

2065:   /* Assemble the matrices so that the correct flags are set */
2066:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
2067:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);

2069:   /* Free work space */
2070:   ISRestoreIndices(isrow,&irow);
2071:   ISRestoreIndices(iscol,&icol);
2072:   *B   = newmat;
2073:   return(0);
2074: }

2076: static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[])
2077: {
2079:   PetscInt       i;

2082:   if (scall == MAT_INITIAL_MATRIX) {
2083:     PetscCalloc1(n+1,B);
2084:   }

2086:   for (i=0; i<n; i++) {
2087:     MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
2088:   }
2089:   return(0);
2090: }

2092: static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
2093: {
2095:   return(0);
2096: }

2098: static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
2099: {
2101:   return(0);
2102: }

2104: static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
2105: {
2106:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data;
2107:   PetscErrorCode    ierr;
2108:   const PetscScalar *va;
2109:   PetscScalar       *vb;
2110:   PetscInt          lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j;

2113:   /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
2114:   if (A->ops->copy != B->ops->copy) {
2115:     MatCopy_Basic(A,B,str);
2116:     return(0);
2117:   }
2118:   if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)");
2119:   MatDenseGetArrayRead(A,&va);
2120:   MatDenseGetArray(B,&vb);
2121:   if (lda1>m || lda2>m) {
2122:     for (j=0; j<n; j++) {
2123:       PetscArraycpy(vb+j*lda2,va+j*lda1,m);
2124:     }
2125:   } else {
2126:     PetscArraycpy(vb,va,A->rmap->n*A->cmap->n);
2127:   }
2128:   MatDenseRestoreArray(B,&vb);
2129:   MatDenseRestoreArrayRead(A,&va);
2130:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
2131:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
2132:   return(0);
2133: }

2135: static PetscErrorCode MatSetUp_SeqDense(Mat A)
2136: {

2140:   PetscLayoutSetUp(A->rmap);
2141:   PetscLayoutSetUp(A->cmap);
2142:   if (!A->preallocated) {
2143:     MatSeqDenseSetPreallocation(A,NULL);
2144:   }
2145:   return(0);
2146: }

2148: static PetscErrorCode MatConjugate_SeqDense(Mat A)
2149: {
2150:   PetscInt       i,nz = A->rmap->n*A->cmap->n;
2151:   PetscScalar    *aa;

2155:   MatDenseGetArray(A,&aa);
2156:   for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
2157:   MatDenseRestoreArray(A,&aa);
2158:   return(0);
2159: }

2161: static PetscErrorCode MatRealPart_SeqDense(Mat A)
2162: {
2163:   PetscInt       i,nz = A->rmap->n*A->cmap->n;
2164:   PetscScalar    *aa;

2168:   MatDenseGetArray(A,&aa);
2169:   for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
2170:   MatDenseRestoreArray(A,&aa);
2171:   return(0);
2172: }

2174: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
2175: {
2176:   PetscInt       i,nz = A->rmap->n*A->cmap->n;
2177:   PetscScalar    *aa;

2181:   MatDenseGetArray(A,&aa);
2182:   for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
2183:   MatDenseRestoreArray(A,&aa);
2184:   return(0);
2185: }

2187: /* ----------------------------------------------------------------*/
2188: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2189: {
2191:   PetscInt       m=A->rmap->n,n=B->cmap->n;
2192:   PetscBool      cisdense;

2195:   MatSetSizes(C,m,n,m,n);
2196:   PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
2197:   if (!cisdense) {
2198:     PetscBool flg;

2200:     PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2201:     MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
2202:   }
2203:   MatSetUp(C);
2204:   return(0);
2205: }

2207: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2208: {
2209:   Mat_SeqDense       *a=(Mat_SeqDense*)A->data,*b=(Mat_SeqDense*)B->data,*c=(Mat_SeqDense*)C->data;
2210:   PetscBLASInt       m,n,k;
2211:   const PetscScalar *av,*bv;
2212:   PetscScalar       *cv;
2213:   PetscScalar       _DOne=1.0,_DZero=0.0;
2214:   PetscErrorCode    ierr;

2217:   PetscBLASIntCast(C->rmap->n,&m);
2218:   PetscBLASIntCast(C->cmap->n,&n);
2219:   PetscBLASIntCast(A->cmap->n,&k);
2220:   if (!m || !n || !k) return(0);
2221:   MatDenseGetArrayRead(A,&av);
2222:   MatDenseGetArrayRead(B,&bv);
2223:   MatDenseGetArrayWrite(C,&cv);
2224:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2225:   PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2226:   MatDenseRestoreArrayRead(A,&av);
2227:   MatDenseRestoreArrayRead(B,&bv);
2228:   MatDenseRestoreArrayWrite(C,&cv);
2229:   return(0);
2230: }

2232: PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2233: {
2235:   PetscInt       m=A->rmap->n,n=B->rmap->n;
2236:   PetscBool      cisdense;

2239:   MatSetSizes(C,m,n,m,n);
2240:   PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
2241:   if (!cisdense) {
2242:     PetscBool flg;

2244:     PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2245:     MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
2246:   }
2247:   MatSetUp(C);
2248:   return(0);
2249: }

2251: PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2252: {
2253:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
2254:   Mat_SeqDense      *b = (Mat_SeqDense*)B->data;
2255:   Mat_SeqDense      *c = (Mat_SeqDense*)C->data;
2256:   const PetscScalar *av,*bv;
2257:   PetscScalar       *cv;
2258:   PetscBLASInt      m,n,k;
2259:   PetscScalar       _DOne=1.0,_DZero=0.0;
2260:   PetscErrorCode    ierr;

2263:   PetscBLASIntCast(C->rmap->n,&m);
2264:   PetscBLASIntCast(C->cmap->n,&n);
2265:   PetscBLASIntCast(A->cmap->n,&k);
2266:   if (!m || !n || !k) return(0);
2267:   MatDenseGetArrayRead(A,&av);
2268:   MatDenseGetArrayRead(B,&bv);
2269:   MatDenseGetArrayWrite(C,&cv);
2270:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2271:   MatDenseRestoreArrayRead(A,&av);
2272:   MatDenseRestoreArrayRead(B,&bv);
2273:   MatDenseRestoreArrayWrite(C,&cv);
2274:   PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2275:   return(0);
2276: }

2278: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2279: {
2281:   PetscInt       m=A->cmap->n,n=B->cmap->n;
2282:   PetscBool      cisdense;

2285:   MatSetSizes(C,m,n,m,n);
2286:   PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
2287:   if (!cisdense) {
2288:     PetscBool flg;

2290:     PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2291:     MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
2292:   }
2293:   MatSetUp(C);
2294:   return(0);
2295: }

2297: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2298: {
2299:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
2300:   Mat_SeqDense      *b = (Mat_SeqDense*)B->data;
2301:   Mat_SeqDense      *c = (Mat_SeqDense*)C->data;
2302:   const PetscScalar *av,*bv;
2303:   PetscScalar       *cv;
2304:   PetscBLASInt      m,n,k;
2305:   PetscScalar       _DOne=1.0,_DZero=0.0;
2306:   PetscErrorCode    ierr;

2309:   PetscBLASIntCast(C->rmap->n,&m);
2310:   PetscBLASIntCast(C->cmap->n,&n);
2311:   PetscBLASIntCast(A->rmap->n,&k);
2312:   if (!m || !n || !k) return(0);
2313:   MatDenseGetArrayRead(A,&av);
2314:   MatDenseGetArrayRead(B,&bv);
2315:   MatDenseGetArrayWrite(C,&cv);
2316:   PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2317:   MatDenseRestoreArrayRead(A,&av);
2318:   MatDenseRestoreArrayRead(B,&bv);
2319:   MatDenseRestoreArrayWrite(C,&cv);
2320:   PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2321:   return(0);
2322: }

2324: /* ----------------------------------------------- */
2325: static PetscErrorCode MatProductSetFromOptions_SeqDense_AB(Mat C)
2326: {
2328:   C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqDense;
2329:   C->ops->productsymbolic = MatProductSymbolic_AB;
2330:   return(0);
2331: }

2333: static PetscErrorCode MatProductSetFromOptions_SeqDense_AtB(Mat C)
2334: {
2336:   C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_SeqDense_SeqDense;
2337:   C->ops->productsymbolic          = MatProductSymbolic_AtB;
2338:   return(0);
2339: }

2341: static PetscErrorCode MatProductSetFromOptions_SeqDense_ABt(Mat C)
2342: {
2344:   C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqDense_SeqDense;
2345:   C->ops->productsymbolic          = MatProductSymbolic_ABt;
2346:   return(0);
2347: }

2349: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense(Mat C)
2350: {
2352:   Mat_Product    *product = C->product;

2355:   switch (product->type) {
2356:   case MATPRODUCT_AB:
2357:     MatProductSetFromOptions_SeqDense_AB(C);
2358:     break;
2359:   case MATPRODUCT_AtB:
2360:     MatProductSetFromOptions_SeqDense_AtB(C);
2361:     break;
2362:   case MATPRODUCT_ABt:
2363:     MatProductSetFromOptions_SeqDense_ABt(C);
2364:     break;
2365:   default:
2366:     break;
2367:   }
2368:   return(0);
2369: }
2370: /* ----------------------------------------------- */

2372: static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[])
2373: {
2374:   Mat_SeqDense       *a = (Mat_SeqDense*)A->data;
2375:   PetscErrorCode     ierr;
2376:   PetscInt           i,j,m = A->rmap->n,n = A->cmap->n,p;
2377:   PetscScalar        *x;
2378:   const PetscScalar *aa;

2381:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2382:   VecGetArray(v,&x);
2383:   VecGetLocalSize(v,&p);
2384:   MatDenseGetArrayRead(A,&aa);
2385:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2386:   for (i=0; i<m; i++) {
2387:     x[i] = aa[i]; if (idx) idx[i] = 0;
2388:     for (j=1; j<n; j++) {
2389:       if (PetscRealPart(x[i]) < PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;}
2390:     }
2391:   }
2392:   MatDenseRestoreArrayRead(A,&aa);
2393:   VecRestoreArray(v,&x);
2394:   return(0);
2395: }

2397: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])
2398: {
2399:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
2400:   PetscErrorCode    ierr;
2401:   PetscInt          i,j,m = A->rmap->n,n = A->cmap->n,p;
2402:   PetscScalar       *x;
2403:   PetscReal         atmp;
2404:   const PetscScalar *aa;

2407:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2408:   VecGetArray(v,&x);
2409:   VecGetLocalSize(v,&p);
2410:   MatDenseGetArrayRead(A,&aa);
2411:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2412:   for (i=0; i<m; i++) {
2413:     x[i] = PetscAbsScalar(aa[i]);
2414:     for (j=1; j<n; j++) {
2415:       atmp = PetscAbsScalar(aa[i+a->lda*j]);
2416:       if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;}
2417:     }
2418:   }
2419:   MatDenseRestoreArrayRead(A,&aa);
2420:   VecRestoreArray(v,&x);
2421:   return(0);
2422: }

2424: static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])
2425: {
2426:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
2427:   PetscErrorCode    ierr;
2428:   PetscInt          i,j,m = A->rmap->n,n = A->cmap->n,p;
2429:   PetscScalar       *x;
2430:   const PetscScalar *aa;

2433:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2434:   MatDenseGetArrayRead(A,&aa);
2435:   VecGetArray(v,&x);
2436:   VecGetLocalSize(v,&p);
2437:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2438:   for (i=0; i<m; i++) {
2439:     x[i] = aa[i]; if (idx) idx[i] = 0;
2440:     for (j=1; j<n; j++) {
2441:       if (PetscRealPart(x[i]) > PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;}
2442:     }
2443:   }
2444:   VecRestoreArray(v,&x);
2445:   MatDenseRestoreArrayRead(A,&aa);
2446:   return(0);
2447: }

2449: PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
2450: {
2451:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
2452:   PetscErrorCode    ierr;
2453:   PetscScalar       *x;
2454:   const PetscScalar *aa;

2457:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2458:   MatDenseGetArrayRead(A,&aa);
2459:   VecGetArray(v,&x);
2460:   PetscArraycpy(x,aa+col*a->lda,A->rmap->n);
2461:   VecRestoreArray(v,&x);
2462:   MatDenseRestoreArrayRead(A,&aa);
2463:   return(0);
2464: }

2466: PETSC_INTERN PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms)
2467: {
2468:   PetscErrorCode    ierr;
2469:   PetscInt          i,j,m,n;
2470:   const PetscScalar *a;

2473:   MatGetSize(A,&m,&n);
2474:   PetscArrayzero(norms,n);
2475:   MatDenseGetArrayRead(A,&a);
2476:   if (type == NORM_2) {
2477:     for (i=0; i<n; i++) {
2478:       for (j=0; j<m; j++) {
2479:         norms[i] += PetscAbsScalar(a[j]*a[j]);
2480:       }
2481:       a += m;
2482:     }
2483:   } else if (type == NORM_1) {
2484:     for (i=0; i<n; i++) {
2485:       for (j=0; j<m; j++) {
2486:         norms[i] += PetscAbsScalar(a[j]);
2487:       }
2488:       a += m;
2489:     }
2490:   } else if (type == NORM_INFINITY) {
2491:     for (i=0; i<n; i++) {
2492:       for (j=0; j<m; j++) {
2493:         norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]);
2494:       }
2495:       a += m;
2496:     }
2497:   } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType");
2498:   MatDenseRestoreArrayRead(A,&a);
2499:   if (type == NORM_2) {
2500:     for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]);
2501:   }
2502:   return(0);
2503: }

2505: static PetscErrorCode  MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
2506: {
2508:   PetscScalar    *a;
2509:   PetscInt       lda,m,n,i,j;

2512:   MatGetSize(x,&m,&n);
2513:   MatDenseGetLDA(x,&lda);
2514:   MatDenseGetArray(x,&a);
2515:   for (j=0; j<n; j++) {
2516:     for (i=0; i<m; i++) {
2517:       PetscRandomGetValue(rctx,a+j*lda+i);
2518:     }
2519:   }
2520:   MatDenseRestoreArray(x,&a);
2521:   return(0);
2522: }

2524: static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool  *missing,PetscInt *d)
2525: {
2527:   *missing = PETSC_FALSE;
2528:   return(0);
2529: }

2531: /* vals is not const */
2532: static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals)
2533: {
2535:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2536:   PetscScalar    *v;

2539:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2540:   MatDenseGetArray(A,&v);
2541:   *vals = v+col*a->lda;
2542:   MatDenseRestoreArray(A,&v);
2543:   return(0);
2544: }

2546: static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals)
2547: {
2549:   *vals = NULL; /* user cannot accidently use the array later */
2550:   return(0);
2551: }

2553: /* -------------------------------------------------------------------*/
2554: static struct _MatOps MatOps_Values = { MatSetValues_SeqDense,
2555:                                         MatGetRow_SeqDense,
2556:                                         MatRestoreRow_SeqDense,
2557:                                         MatMult_SeqDense,
2558:                                 /*  4*/ MatMultAdd_SeqDense,
2559:                                         MatMultTranspose_SeqDense,
2560:                                         MatMultTransposeAdd_SeqDense,
2561:                                         NULL,
2562:                                         NULL,
2563:                                         NULL,
2564:                                 /* 10*/ NULL,
2565:                                         MatLUFactor_SeqDense,
2566:                                         MatCholeskyFactor_SeqDense,
2567:                                         MatSOR_SeqDense,
2568:                                         MatTranspose_SeqDense,
2569:                                 /* 15*/ MatGetInfo_SeqDense,
2570:                                         MatEqual_SeqDense,
2571:                                         MatGetDiagonal_SeqDense,
2572:                                         MatDiagonalScale_SeqDense,
2573:                                         MatNorm_SeqDense,
2574:                                 /* 20*/ MatAssemblyBegin_SeqDense,
2575:                                         MatAssemblyEnd_SeqDense,
2576:                                         MatSetOption_SeqDense,
2577:                                         MatZeroEntries_SeqDense,
2578:                                 /* 24*/ MatZeroRows_SeqDense,
2579:                                         NULL,
2580:                                         NULL,
2581:                                         NULL,
2582:                                         NULL,
2583:                                 /* 29*/ MatSetUp_SeqDense,
2584:                                         NULL,
2585:                                         NULL,
2586:                                         NULL,
2587:                                         NULL,
2588:                                 /* 34*/ MatDuplicate_SeqDense,
2589:                                         NULL,
2590:                                         NULL,
2591:                                         NULL,
2592:                                         NULL,
2593:                                 /* 39*/ MatAXPY_SeqDense,
2594:                                         MatCreateSubMatrices_SeqDense,
2595:                                         NULL,
2596:                                         MatGetValues_SeqDense,
2597:                                         MatCopy_SeqDense,
2598:                                 /* 44*/ MatGetRowMax_SeqDense,
2599:                                         MatScale_SeqDense,
2600:                                         MatShift_Basic,
2601:                                         NULL,
2602:                                         MatZeroRowsColumns_SeqDense,
2603:                                 /* 49*/ MatSetRandom_SeqDense,
2604:                                         NULL,
2605:                                         NULL,
2606:                                         NULL,
2607:                                         NULL,
2608:                                 /* 54*/ NULL,
2609:                                         NULL,
2610:                                         NULL,
2611:                                         NULL,
2612:                                         NULL,
2613:                                 /* 59*/ NULL,
2614:                                         MatDestroy_SeqDense,
2615:                                         MatView_SeqDense,
2616:                                         NULL,
2617:                                         NULL,
2618:                                 /* 64*/ NULL,
2619:                                         NULL,
2620:                                         NULL,
2621:                                         NULL,
2622:                                         NULL,
2623:                                 /* 69*/ MatGetRowMaxAbs_SeqDense,
2624:                                         NULL,
2625:                                         NULL,
2626:                                         NULL,
2627:                                         NULL,
2628:                                 /* 74*/ NULL,
2629:                                         NULL,
2630:                                         NULL,
2631:                                         NULL,
2632:                                         NULL,
2633:                                 /* 79*/ NULL,
2634:                                         NULL,
2635:                                         NULL,
2636:                                         NULL,
2637:                                 /* 83*/ MatLoad_SeqDense,
2638:                                         MatIsSymmetric_SeqDense,
2639:                                         MatIsHermitian_SeqDense,
2640:                                         NULL,
2641:                                         NULL,
2642:                                         NULL,
2643:                                 /* 89*/ NULL,
2644:                                         NULL,
2645:                                         MatMatMultNumeric_SeqDense_SeqDense,
2646:                                         NULL,
2647:                                         NULL,
2648:                                 /* 94*/ NULL,
2649:                                         NULL,
2650:                                         NULL,
2651:                                         MatMatTransposeMultNumeric_SeqDense_SeqDense,
2652:                                         NULL,
2653:                                 /* 99*/ MatProductSetFromOptions_SeqDense,
2654:                                         NULL,
2655:                                         NULL,
2656:                                         MatConjugate_SeqDense,
2657:                                         NULL,
2658:                                 /*104*/ NULL,
2659:                                         MatRealPart_SeqDense,
2660:                                         MatImaginaryPart_SeqDense,
2661:                                         NULL,
2662:                                         NULL,
2663:                                 /*109*/ NULL,
2664:                                         NULL,
2665:                                         MatGetRowMin_SeqDense,
2666:                                         MatGetColumnVector_SeqDense,
2667:                                         MatMissingDiagonal_SeqDense,
2668:                                 /*114*/ NULL,
2669:                                         NULL,
2670:                                         NULL,
2671:                                         NULL,
2672:                                         NULL,
2673:                                 /*119*/ NULL,
2674:                                         NULL,
2675:                                         NULL,
2676:                                         NULL,
2677:                                         NULL,
2678:                                 /*124*/ NULL,
2679:                                         MatGetColumnNorms_SeqDense,
2680:                                         NULL,
2681:                                         NULL,
2682:                                         NULL,
2683:                                 /*129*/ NULL,
2684:                                         NULL,
2685:                                         NULL,
2686:                                         MatTransposeMatMultNumeric_SeqDense_SeqDense,
2687:                                         NULL,
2688:                                 /*134*/ NULL,
2689:                                         NULL,
2690:                                         NULL,
2691:                                         NULL,
2692:                                         NULL,
2693:                                 /*139*/ NULL,
2694:                                         NULL,
2695:                                         NULL,
2696:                                         NULL,
2697:                                         NULL,
2698:                                         MatCreateMPIMatConcatenateSeqMat_SeqDense,
2699:                                 /*145*/ NULL,
2700:                                         NULL,
2701:                                         NULL
2702: };

2704: /*@C
2705:    MatCreateSeqDense - Creates a sequential dense matrix that
2706:    is stored in column major order (the usual Fortran 77 manner). Many
2707:    of the matrix operations use the BLAS and LAPACK routines.

2709:    Collective

2711:    Input Parameters:
2712: +  comm - MPI communicator, set to PETSC_COMM_SELF
2713: .  m - number of rows
2714: .  n - number of columns
2715: -  data - optional location of matrix data in column major order.  Set data=NULL for PETSc
2716:    to control all matrix memory allocation.

2718:    Output Parameter:
2719: .  A - the matrix

2721:    Notes:
2722:    The data input variable is intended primarily for Fortran programmers
2723:    who wish to allocate their own matrix memory space.  Most users should
2724:    set data=NULL.

2726:    Level: intermediate

2728: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
2729: @*/
2730: PetscErrorCode  MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
2731: {

2735:   MatCreate(comm,A);
2736:   MatSetSizes(*A,m,n,m,n);
2737:   MatSetType(*A,MATSEQDENSE);
2738:   MatSeqDenseSetPreallocation(*A,data);
2739:   return(0);
2740: }

2742: /*@C
2743:    MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements

2745:    Collective

2747:    Input Parameters:
2748: +  B - the matrix
2749: -  data - the array (or NULL)

2751:    Notes:
2752:    The data input variable is intended primarily for Fortran programmers
2753:    who wish to allocate their own matrix memory space.  Most users should
2754:    need not call this routine.

2756:    Level: intermediate

2758: .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatDenseSetLDA()

2760: @*/
2761: PetscErrorCode  MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
2762: {

2767:   PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
2768:   return(0);
2769: }

2771: PetscErrorCode  MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
2772: {
2773:   Mat_SeqDense   *b = (Mat_SeqDense*)B->data;

2777:   if (b->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
2778:   B->preallocated = PETSC_TRUE;

2780:   PetscLayoutSetUp(B->rmap);
2781:   PetscLayoutSetUp(B->cmap);

2783:   if (b->lda <= 0) b->lda = B->rmap->n;

2785:   PetscIntMultError(b->lda,B->cmap->n,NULL);
2786:   if (!data) { /* petsc-allocated storage */
2787:     if (!b->user_alloc) { PetscFree(b->v); }
2788:     PetscCalloc1((size_t)b->lda*B->cmap->n,&b->v);
2789:     PetscLogObjectMemory((PetscObject)B,b->lda*B->cmap->n*sizeof(PetscScalar));

2791:     b->user_alloc = PETSC_FALSE;
2792:   } else { /* user-allocated storage */
2793:     if (!b->user_alloc) { PetscFree(b->v); }
2794:     b->v          = data;
2795:     b->user_alloc = PETSC_TRUE;
2796:   }
2797:   B->assembled = PETSC_TRUE;
2798:   return(0);
2799: }

2801: #if defined(PETSC_HAVE_ELEMENTAL)
2802: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
2803: {
2804:   Mat               mat_elemental;
2805:   PetscErrorCode    ierr;
2806:   const PetscScalar *array;
2807:   PetscScalar       *v_colwise;
2808:   PetscInt          M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols;

2811:   PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);
2812:   MatDenseGetArrayRead(A,&array);
2813:   /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
2814:   k = 0;
2815:   for (j=0; j<N; j++) {
2816:     cols[j] = j;
2817:     for (i=0; i<M; i++) {
2818:       v_colwise[j*M+i] = array[k++];
2819:     }
2820:   }
2821:   for (i=0; i<M; i++) {
2822:     rows[i] = i;
2823:   }
2824:   MatDenseRestoreArrayRead(A,&array);

2826:   MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
2827:   MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);
2828:   MatSetType(mat_elemental,MATELEMENTAL);
2829:   MatSetUp(mat_elemental);

2831:   /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
2832:   MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);
2833:   MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
2834:   MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);
2835:   PetscFree3(v_colwise,rows,cols);

2837:   if (reuse == MAT_INPLACE_MATRIX) {
2838:     MatHeaderReplace(A,&mat_elemental);
2839:   } else {
2840:     *newmat = mat_elemental;
2841:   }
2842:   return(0);
2843: }
2844: #endif

2846: PetscErrorCode  MatDenseSetLDA_SeqDense(Mat B,PetscInt lda)
2847: {
2848:   Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2849:   PetscBool    data;

2852:   data = (PetscBool)((B->rmap->n > 0 && B->cmap->n > 0) ? (b->v ? PETSC_TRUE : PETSC_FALSE) : PETSC_FALSE);
2853:   if (!b->user_alloc && data && b->lda!=lda) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"LDA cannot be changed after allocation of internal storage");
2854:   if (lda < B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"LDA %D must be at least matrix dimension %D",lda,B->rmap->n);
2855:   b->lda = lda;
2856:   return(0);
2857: }

2859: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat)
2860: {
2862:   PetscMPIInt    size;

2865:   MPI_Comm_size(comm,&size);
2866:   if (size == 1) {
2867:     if (scall == MAT_INITIAL_MATRIX) {
2868:       MatDuplicate(inmat,MAT_COPY_VALUES,outmat);
2869:     } else {
2870:       MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);
2871:     }
2872:   } else {
2873:     MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);
2874:   }
2875:   return(0);
2876: }

2878: PetscErrorCode MatDenseGetColumnVec_SeqDense(Mat A,PetscInt col,Vec *v)
2879: {
2880:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

2884:   if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
2885:   if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
2886:   if (!a->cvec) {
2887:     VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);
2888:     PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);
2889:   }
2890:   a->vecinuse = col + 1;
2891:   MatDenseGetArray(A,(PetscScalar**)&a->ptrinuse);
2892:   VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);
2893:   *v   = a->cvec;
2894:   return(0);
2895: }

2897: PetscErrorCode MatDenseRestoreColumnVec_SeqDense(Mat A,PetscInt col,Vec *v)
2898: {
2899:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

2903:   if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first");
2904:   if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector");
2905:   a->vecinuse = 0;
2906:   MatDenseRestoreArray(A,(PetscScalar**)&a->ptrinuse);
2907:   VecResetArray(a->cvec);
2908:   *v   = NULL;
2909:   return(0);
2910: }

2912: PetscErrorCode MatDenseGetColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v)
2913: {
2914:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

2918:   if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
2919:   if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
2920:   if (!a->cvec) {
2921:     VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);
2922:     PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);
2923:   }
2924:   a->vecinuse = col + 1;
2925:   MatDenseGetArrayRead(A,&a->ptrinuse);
2926:   VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);
2927:   VecLockReadPush(a->cvec);
2928:   *v   = a->cvec;
2929:   return(0);
2930: }

2932: PetscErrorCode MatDenseRestoreColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v)
2933: {
2934:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

2938:   if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first");
2939:   if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector");
2940:   a->vecinuse = 0;
2941:   MatDenseRestoreArrayRead(A,&a->ptrinuse);
2942:   VecLockReadPop(a->cvec);
2943:   VecResetArray(a->cvec);
2944:   *v   = NULL;
2945:   return(0);
2946: }

2948: PetscErrorCode MatDenseGetColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v)
2949: {
2950:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

2954:   if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
2955:   if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
2956:   if (!a->cvec) {
2957:     VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);
2958:     PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);
2959:   }
2960:   a->vecinuse = col + 1;
2961:   MatDenseGetArrayWrite(A,(PetscScalar**)&a->ptrinuse);
2962:   VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);
2963:   *v   = a->cvec;
2964:   return(0);
2965: }

2967: PetscErrorCode MatDenseRestoreColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v)
2968: {
2969:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

2973:   if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first");
2974:   if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector");
2975:   a->vecinuse = 0;
2976:   MatDenseRestoreArrayWrite(A,(PetscScalar**)&a->ptrinuse);
2977:   VecResetArray(a->cvec);
2978:   *v   = NULL;
2979:   return(0);
2980: }

2982: PetscErrorCode MatDenseGetSubMatrix_SeqDense(Mat A,PetscInt cbegin,PetscInt cend,Mat *v)
2983: {
2984:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

2988:   if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
2989:   if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
2990:   if (a->cmat && cend-cbegin != a->cmat->cmap->N) {
2991:     MatDestroy(&a->cmat);
2992:   }
2993:   if (!a->cmat) {
2994:     MatCreateDense(PetscObjectComm((PetscObject)A),A->rmap->n,PETSC_DECIDE,A->rmap->N,cend-cbegin,(PetscScalar*)a->v + (size_t)cbegin * (size_t)a->lda,&a->cmat);
2995:     PetscLogObjectParent((PetscObject)A,(PetscObject)a->cmat);
2996:   } else {
2997:     MatDensePlaceArray(a->cmat,a->v + (size_t)cbegin * (size_t)a->lda);
2998:   }
2999:   MatDenseSetLDA(a->cmat,a->lda);
3000:   a->matinuse = cbegin + 1;
3001:   *v = a->cmat;
3002:   return(0);
3003: }

3005: PetscErrorCode MatDenseRestoreSubMatrix_SeqDense(Mat A,Mat *v)
3006: {
3007:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

3011:   if (!a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetSubMatrix() first");
3012:   if (!a->cmat) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column matrix");
3013:   if (*v != a->cmat) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not the matrix obtained from MatDenseGetSubMatrix()");
3014:   a->matinuse = 0;
3015:   MatDenseResetArray(a->cmat);
3016:   *v   = NULL;
3017:   return(0);
3018: }

3020: /*MC
3021:    MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices.

3023:    Options Database Keys:
3024: . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions()

3026:   Level: beginner

3028: .seealso: MatCreateSeqDense()

3030: M*/
3031: PetscErrorCode MatCreate_SeqDense(Mat B)
3032: {
3033:   Mat_SeqDense   *b;
3035:   PetscMPIInt    size;

3038:   MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
3039:   if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1");

3041:   PetscNewLog(B,&b);
3042:   PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
3043:   B->data = (void*)b;

3045:   b->roworiented = PETSC_TRUE;

3047:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetLDA_C",MatDenseGetLDA_SeqDense);
3048:   PetscObjectComposeFunction((PetscObject)B,"MatDenseSetLDA_C",MatDenseSetLDA_SeqDense);
3049:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);
3050:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);
3051:   PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);
3052:   PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);
3053:   PetscObjectComposeFunction((PetscObject)B,"MatDenseReplaceArray_C",MatDenseReplaceArray_SeqDense);
3054:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);
3055:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense);
3056:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayWrite_C",MatDenseGetArray_SeqDense);
3057:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayWrite_C",MatDenseRestoreArray_SeqDense);
3058:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);
3059: #if defined(PETSC_HAVE_ELEMENTAL)
3060:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);
3061: #endif
3062: #if defined(PETSC_HAVE_SCALAPACK)
3063:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_scalapack_C",MatConvert_Dense_ScaLAPACK);
3064: #endif
3065: #if defined(PETSC_HAVE_CUDA)
3066:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqdensecuda_C",MatConvert_SeqDense_SeqDenseCUDA);
3067:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",MatProductSetFromOptions_SeqDense);
3068:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdense_C",MatProductSetFromOptions_SeqDense);
3069:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdensecuda_C",MatProductSetFromOptions_SeqDense);
3070: #endif
3071:   PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);
3072:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqdense_C",MatProductSetFromOptions_SeqAIJ_SeqDense);
3073:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdense_C",MatProductSetFromOptions_SeqDense);
3074:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);
3075:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqsbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);

3077:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);
3078:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);
3079:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVec_C",MatDenseGetColumnVec_SeqDense);
3080:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVec_C",MatDenseRestoreColumnVec_SeqDense);
3081:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecRead_C",MatDenseGetColumnVecRead_SeqDense);
3082:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecRead_C",MatDenseRestoreColumnVecRead_SeqDense);
3083:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecWrite_C",MatDenseGetColumnVecWrite_SeqDense);
3084:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecWrite_C",MatDenseRestoreColumnVecWrite_SeqDense);
3085:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetSubMatrix_C",MatDenseGetSubMatrix_SeqDense);
3086:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreSubMatrix_C",MatDenseRestoreSubMatrix_SeqDense);
3087:   PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
3088:   return(0);
3089: }

3091: /*@C
3092:    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.

3094:    Not Collective

3096:    Input Parameters:
3097: +  mat - a MATSEQDENSE or MATMPIDENSE matrix
3098: -  col - column index

3100:    Output Parameter:
3101: .  vals - pointer to the data

3103:    Level: intermediate

3105: .seealso: MatDenseRestoreColumn()
3106: @*/
3107: PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals)
3108: {

3115:   PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));
3116:   return(0);
3117: }

3119: /*@C
3120:    MatDenseRestoreColumn - returns access to a column of a dense matrix which is returned by MatDenseGetColumn().

3122:    Not Collective

3124:    Input Parameter:
3125: .  mat - a MATSEQDENSE or MATMPIDENSE matrix

3127:    Output Parameter:
3128: .  vals - pointer to the data

3130:    Level: intermediate

3132: .seealso: MatDenseGetColumn()
3133: @*/
3134: PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals)
3135: {

3141:   PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));
3142:   return(0);
3143: }

3145: /*@C
3146:    MatDenseGetColumnVec - Gives read-write access to a column of a dense matrix, represented as a Vec.

3148:    Collective

3150:    Input Parameters:
3151: +  mat - the Mat object
3152: -  col - the column index

3154:    Output Parameter:
3155: .  v - the vector

3157:    Notes:
3158:      The vector is owned by PETSc. Users need to call MatDenseRestoreColumnVec() when the vector is no longer needed.
3159:      Use MatDenseGetColumnVecRead() to obtain read-only access or MatDenseGetColumnVecWrite() for write-only access.

3161:    Level: intermediate

3163: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3164: @*/
3165: PetscErrorCode MatDenseGetColumnVec(Mat A,PetscInt col,Vec *v)
3166: {

3174:   if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3175:   if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N);
3176:   PetscUseMethod(A,"MatDenseGetColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));
3177:   return(0);
3178: }

3180: /*@C
3181:    MatDenseRestoreColumnVec - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVec().

3183:    Collective

3185:    Input Parameters:
3186: +  mat - the Mat object
3187: .  col - the column index
3188: -  v - the Vec object

3190:    Level: intermediate

3192: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3193: @*/
3194: PetscErrorCode MatDenseRestoreColumnVec(Mat A,PetscInt col,Vec *v)
3195: {

3203:   if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3204:   if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N);
3205:   PetscUseMethod(A,"MatDenseRestoreColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));
3206:   return(0);
3207: }

3209: /*@C
3210:    MatDenseGetColumnVecRead - Gives read-only access to a column of a dense matrix, represented as a Vec.

3212:    Collective

3214:    Input Parameters:
3215: +  mat - the Mat object
3216: -  col - the column index

3218:    Output Parameter:
3219: .  v - the vector

3221:    Notes:
3222:      The vector is owned by PETSc and users cannot modify it.
3223:      Users need to call MatDenseRestoreColumnVecRead() when the vector is no longer needed.
3224:      Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecWrite() for write-only access.

3226:    Level: intermediate

3228: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3229: @*/
3230: PetscErrorCode MatDenseGetColumnVecRead(Mat A,PetscInt col,Vec *v)
3231: {

3239:   if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3240:   if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N);
3241:   PetscUseMethod(A,"MatDenseGetColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));
3242:   return(0);
3243: }

3245: /*@C
3246:    MatDenseRestoreColumnVecRead - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecRead().

3248:    Collective

3250:    Input Parameters:
3251: +  mat - the Mat object
3252: .  col - the column index
3253: -  v - the Vec object

3255:    Level: intermediate

3257: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecWrite()
3258: @*/
3259: PetscErrorCode MatDenseRestoreColumnVecRead(Mat A,PetscInt col,Vec *v)
3260: {

3268:   if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3269:   if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N);
3270:   PetscUseMethod(A,"MatDenseRestoreColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));
3271:   return(0);
3272: }

3274: /*@C
3275:    MatDenseGetColumnVecWrite - Gives write-only access to a column of a dense matrix, represented as a Vec.

3277:    Collective

3279:    Input Parameters:
3280: +  mat - the Mat object
3281: -  col - the column index

3283:    Output Parameter:
3284: .  v - the vector

3286:    Notes:
3287:      The vector is owned by PETSc. Users need to call MatDenseRestoreColumnVecWrite() when the vector is no longer needed.
3288:      Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecRead() for read-only access.

3290:    Level: intermediate

3292: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3293: @*/
3294: PetscErrorCode MatDenseGetColumnVecWrite(Mat A,PetscInt col,Vec *v)
3295: {

3303:   if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3304:   if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N);
3305:   PetscUseMethod(A,"MatDenseGetColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));
3306:   return(0);
3307: }

3309: /*@C
3310:    MatDenseRestoreColumnVecWrite - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecWrite().

3312:    Collective

3314:    Input Parameters:
3315: +  mat - the Mat object
3316: .  col - the column index
3317: -  v - the Vec object

3319:    Level: intermediate

3321: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead()
3322: @*/
3323: PetscErrorCode MatDenseRestoreColumnVecWrite(Mat A,PetscInt col,Vec *v)
3324: {

3332:   if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3333:   if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N);
3334:   PetscUseMethod(A,"MatDenseRestoreColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));
3335:   return(0);
3336: }

3338: /*@C
3339:    MatDenseGetSubMatrix - Gives access to a block of columns of a dense matrix, represented as a Mat.

3341:    Collective

3343:    Input Parameters:
3344: +  mat - the Mat object
3345: .  cbegin - the first index in the block
3346: -  cend - the last index in the block

3348:    Output Parameter:
3349: .  v - the matrix

3351:    Notes:
3352:      The matrix is owned by PETSc. Users need to call MatDenseRestoreSubMatrix() when the matrix is no longer needed.

3354:    Level: intermediate

3356: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseRestoreColumnVec(), MatDenseRestoreSubMatrix()
3357: @*/
3358: PetscErrorCode MatDenseGetSubMatrix(Mat A,PetscInt cbegin,PetscInt cend,Mat *v)
3359: {

3368:   if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3369:   if (cbegin < 0 || cbegin > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid cbegin %D, should be in [0,%D)",cbegin,A->cmap->N);
3370:   if (cend < cbegin || cend > A->cmap->N) SETERRQ3(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid cend %D, should be in [%D,%D)",cend,cbegin,A->cmap->N);
3371:   PetscUseMethod(A,"MatDenseGetSubMatrix_C",(Mat,PetscInt,PetscInt,Mat*),(A,cbegin,cend,v));
3372:   return(0);
3373: }

3375: /*@C
3376:    MatDenseRestoreSubMatrix - Returns access to a block of columns of a dense matrix obtained from MatDenseGetSubMatrix().

3378:    Collective

3380:    Input Parameters:
3381: +  mat - the Mat object
3382: -  v - the Mat object

3384:    Level: intermediate

3386: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseRestoreColumnVec(), MatDenseGetSubMatrix()
3387: @*/
3388: PetscErrorCode MatDenseRestoreSubMatrix(Mat A,Mat *v)
3389: {

3396:   PetscUseMethod(A,"MatDenseRestoreSubMatrix_C",(Mat,Mat*),(A,v));
3397:   return(0);
3398: }