Actual source code: dense.c

petsc-master 2020-08-09
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,lday,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;

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,0,0,&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,0,&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) {
749:     (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
750:   } else {
751:     (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
752:   }
753:   (*fact)->factortype = ftype;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1012:   MatDenseGetArrayRead(A,&vv);
1013:   /* row-oriented output */
1014:   for (i=0; i<m; i++) {
1015:     if (indexm[i] < 0) {v += n;continue;}
1016:     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);
1017:     for (j=0; j<n; j++) {
1018:       if (indexn[j] < 0) {v++; continue;}
1019:       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);
1020:       *v++ = vv[indexn[j]*mat->lda + indexm[i]];
1021:     }
1022:   }
1023:   MatDenseRestoreArrayRead(A,&vv);
1024:   return(0);
1025: }

1027: /* -----------------------------------------------------------------*/

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

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

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

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

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

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

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

1090:   if (!skipHeader) {
1091:     PetscViewerBinaryRead(viewer,header,4,NULL,PETSC_INT);
1092:     if (header[0] != MAT_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Not a matrix object in file");
1093:     M = header[1]; N = header[2];
1094:     if (M < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix row size (%D) in file is negative",M);
1095:     if (N < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix column size (%D) in file is negative",N);
1096:     nz = header[3];
1097:     if (nz != MATRIX_BINARY_FORMAT_DENSE && nz < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Unknown matrix format %D in file",nz);
1098:   } else {
1099:     MatGetSize(mat,&M,&N);
1100:     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");
1101:     nz = MATRIX_BINARY_FORMAT_DENSE;
1102:   }

1104:   /* setup global sizes if not set */
1105:   if (mat->rmap->N < 0) mat->rmap->N = M;
1106:   if (mat->cmap->N < 0) mat->cmap->N = N;
1107:   MatSetUp(mat);
1108:   /* check if global sizes are correct */
1109:   MatGetSize(mat,&rows,&cols);
1110:   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);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1756:   if (PetscDefined(USE_DEBUG)) {
1757:     for (i=0; i<N; i++) {
1758:       if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
1759:       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);
1760:     }
1761:   }
1762:   if (!N) return(0);

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

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

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

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

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

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

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

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

1818:    Not collective

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

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

1826:    Level: intermediate

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

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

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

1844:    Not collective

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

1850:    Level: intermediate

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

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

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

1867:    Logically Collective on Mat

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

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

1875:    Level: intermediate

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

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

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

1893:    Logically Collective on Mat

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

1899:    Level: intermediate

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

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

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

1921:    Not Collective

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

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

1929:    Level: intermediate

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

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

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

1947:    Not Collective

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

1953:    Level: intermediate

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

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

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

1971:    Not Collective

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

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

1979:    Level: intermediate

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

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

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

1997:    Not Collective

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

2003:    Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2708:    Collective

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

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

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

2725:    Level: intermediate

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

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

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

2744:    Collective

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

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

2755:    Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2851:   data = (PetscBool)((B->rmap->n > 0 && B->cmap->n > 0) ? (b->v ? PETSC_TRUE : PETSC_FALSE) : PETSC_FALSE);
2852:   if (!b->user_alloc && data && b->lda!=lda) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"LDA cannot be changed after allocation of internal storage");
2853:   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);
2854:   b->lda = lda;
2855:   return(0);
2856: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3025:   Level: beginner

3027: .seealso: MatCreateSeqDense()

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

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

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

3044:   b->roworiented = PETSC_TRUE;

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

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

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

3093:    Not Collective

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

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

3102:    Level: intermediate

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

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

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

3121:    Not Collective

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

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

3129:    Level: intermediate

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

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

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

3147:    Collective

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

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

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

3160:    Level: intermediate

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

3173:   if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3174:   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);
3175:   PetscUseMethod(A,"MatDenseGetColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));
3176:   return(0);
3177: }

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

3182:    Collective

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

3189:    Level: intermediate

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

3202:   if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3203:   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);
3204:   PetscUseMethod(A,"MatDenseRestoreColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));
3205:   return(0);
3206: }

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

3211:    Collective

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

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

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

3225:    Level: intermediate

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

3238:   if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3239:   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);
3240:   PetscUseMethod(A,"MatDenseGetColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));
3241:   return(0);
3242: }

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

3247:    Collective

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

3254:    Level: intermediate

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

3267:   if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3268:   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);
3269:   PetscUseMethod(A,"MatDenseRestoreColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));
3270:   return(0);
3271: }

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

3276:    Collective

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

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

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

3289:    Level: intermediate

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

3302:   if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3303:   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);
3304:   PetscUseMethod(A,"MatDenseGetColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));
3305:   return(0);
3306: }

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

3311:    Collective

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

3318:    Level: intermediate

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

3331:   if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3332:   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);
3333:   PetscUseMethod(A,"MatDenseRestoreColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));
3334:   return(0);
3335: }

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

3340:    Collective

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

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

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

3353:    Level: intermediate

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

3367:   if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3368:   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);
3369:   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);
3370:   PetscUseMethod(A,"MatDenseGetSubMatrix_C",(Mat,PetscInt,PetscInt,Mat*),(A,cbegin,cend,v));
3371:   return(0);
3372: }

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

3377:    Collective

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

3383:    Level: intermediate

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

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