Actual source code: dense.c

petsc-master 2019-11-13
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: #if defined(PETSC_MISSING_LAPACK_POTRF)
 42:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable.");
 43: #else
 44:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
 46:   PetscBLASInt   info,n;

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

 90:   A->ops->solve             = NULL;
 91:   A->ops->matsolve          = NULL;
 92:   A->ops->solvetranspose    = NULL;
 93:   A->ops->matsolvetranspose = NULL;
 94:   A->ops->solveadd          = NULL;
 95:   A->ops->solvetransposeadd = NULL;
 96:   A->factortype             = MAT_FACTOR_NONE;
 97:   PetscFree(A->solvertype);
 98:   return(0);
 99: }

101: PetscErrorCode MatZeroRowsColumns_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
102: {
103:   PetscErrorCode    ierr;
104:   Mat_SeqDense      *l = (Mat_SeqDense*)A->data;
105:   PetscInt          m  = l->lda, n = A->cmap->n,r = A->rmap->n, i,j;
106:   PetscScalar       *slot,*bb,*v;
107:   const PetscScalar *xx;

110: #if defined(PETSC_USE_DEBUG)
111:   for (i=0; i<N; i++) {
112:     if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
113:     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);
114:     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);
115:   }
116: #endif
117:   if (!N) return(0);

119:   /* fix right hand side if needed */
120:   if (x && b) {
121:     Vec xt;

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

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

156: PetscErrorCode MatPtAPNumeric_SeqDense_SeqDense(Mat A,Mat P,Mat C)
157: {
158:   Mat_SeqDense   *c = (Mat_SeqDense*)(C->data);

162:   if (c->ptapwork) {
163:     (*C->ops->matmultnumeric)(A,P,c->ptapwork);
164:     (*C->ops->transposematmultnumeric)(P,c->ptapwork,C);
165:   } else { /* first time went trough the basic. Should we add better dispatching for subclasses? */
166:     MatPtAP_Basic(A,P,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);
167:   }
168:   return(0);
169: }

171: PetscErrorCode MatPtAPSymbolic_SeqDense_SeqDense(Mat A,Mat P,PetscReal fill,Mat *C)
172: {
173:   Mat_SeqDense   *c;
174:   PetscBool      flg;

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

191: PETSC_INTERN PetscErrorCode MatPtAP_SeqDense_SeqDense(Mat A,Mat P,MatReuse reuse,PetscReal fill,Mat *C)
192: {

196:   if (reuse == MAT_INITIAL_MATRIX) {
197:     MatPtAPSymbolic_SeqDense_SeqDense(A,P,fill,C);
198:   }
199:   PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);
200:   (*(*C)->ops->ptapnumeric)(A,P,*C);
201:   PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);
202:   return(0);
203: }

205: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
206: {
207:   Mat            B = NULL;
208:   Mat_SeqAIJ     *a = (Mat_SeqAIJ*)A->data;
209:   Mat_SeqDense   *b;
211:   PetscInt       *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->N,i;
212:   MatScalar      *av=a->a;
213:   PetscBool      isseqdense;

216:   if (reuse == MAT_REUSE_MATRIX) {
217:     PetscObjectTypeCompare((PetscObject)*newmat,MATSEQDENSE,&isseqdense);
218:     if (!isseqdense) SETERRQ1(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_USER,"Cannot reuse matrix of type %s",((PetscObject)(*newmat))->type_name);
219:   }
220:   if (reuse != MAT_REUSE_MATRIX) {
221:     MatCreate(PetscObjectComm((PetscObject)A),&B);
222:     MatSetSizes(B,m,n,m,n);
223:     MatSetType(B,MATSEQDENSE);
224:     MatSeqDenseSetPreallocation(B,NULL);
225:     b    = (Mat_SeqDense*)(B->data);
226:   } else {
227:     b    = (Mat_SeqDense*)((*newmat)->data);
228:     PetscArrayzero(b->v,m*n);
229:   }
230:   for (i=0; i<m; i++) {
231:     PetscInt j;
232:     for (j=0;j<ai[1]-ai[0];j++) {
233:       b->v[*aj*m+i] = *av;
234:       aj++;
235:       av++;
236:     }
237:     ai++;
238:   }

240:   if (reuse == MAT_INPLACE_MATRIX) {
241:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
242:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
243:     MatHeaderReplace(A,&B);
244:   } else {
245:     if (B) *newmat = B;
246:     MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
247:     MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
248:   }
249:   return(0);
250: }

252: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
253: {
254:   Mat            B;
255:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
257:   PetscInt       i, j;
258:   PetscInt       *rows, *nnz;
259:   MatScalar      *aa = a->v, *vals;

262:   MatCreate(PetscObjectComm((PetscObject)A),&B);
263:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
264:   MatSetType(B,MATSEQAIJ);
265:   PetscCalloc3(A->rmap->n,&rows,A->rmap->n,&nnz,A->rmap->n,&vals);
266:   for (j=0; j<A->cmap->n; j++) {
267:     for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) ++nnz[i];
268:     aa += a->lda;
269:   }
270:   MatSeqAIJSetPreallocation(B,PETSC_DETERMINE,nnz);
271:   aa = a->v;
272:   for (j=0; j<A->cmap->n; j++) {
273:     PetscInt numRows = 0;
274:     for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) {rows[numRows] = i; vals[numRows++] = aa[i];}
275:     MatSetValues(B,numRows,rows,1,&j,vals,INSERT_VALUES);
276:     aa  += a->lda;
277:   }
278:   PetscFree3(rows,nnz,vals);
279:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
280:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

282:   if (reuse == MAT_INPLACE_MATRIX) {
283:     MatHeaderReplace(A,&B);
284:   } else {
285:     *newmat = B;
286:   }
287:   return(0);
288: }

290: PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str)
291: {
292:   Mat_SeqDense      *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data;
293:   const PetscScalar *xv;
294:   PetscScalar       *yv;
295:   PetscBLASInt      N,m,ldax,lday,one = 1;
296:   PetscErrorCode    ierr;

299:   MatDenseGetArrayRead(X,&xv);
300:   MatDenseGetArray(Y,&yv);
301:   PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);
302:   PetscBLASIntCast(X->rmap->n,&m);
303:   PetscBLASIntCast(x->lda,&ldax);
304:   PetscBLASIntCast(y->lda,&lday);
305:   if (ldax>m || lday>m) {
306:     PetscInt j;

308:     for (j=0; j<X->cmap->n; j++) {
309:       PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&alpha,xv+j*ldax,&one,yv+j*lday,&one));
310:     }
311:   } else {
312:     PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&alpha,xv,&one,yv,&one));
313:   }
314:   MatDenseRestoreArrayRead(X,&xv);
315:   MatDenseRestoreArray(Y,&yv);
316:   PetscLogFlops(PetscMax(2*N-1,0));
317:   return(0);
318: }

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

325:   info->block_size        = 1.0;
326:   info->nz_allocated      = N;
327:   info->nz_used           = N;
328:   info->nz_unneeded       = 0;
329:   info->assemblies        = A->num_ass;
330:   info->mallocs           = 0;
331:   info->memory            = ((PetscObject)A)->mem;
332:   info->fill_ratio_given  = 0;
333:   info->fill_ratio_needed = 0;
334:   info->factor_mallocs    = 0;
335:   return(0);
336: }

338: static PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha)
339: {
340:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
341:   PetscScalar    *v;
343:   PetscBLASInt   one = 1,j,nz,lda;

346:   MatDenseGetArray(A,&v);
347:   PetscBLASIntCast(a->lda,&lda);
348:   if (lda>A->rmap->n) {
349:     PetscBLASIntCast(A->rmap->n,&nz);
350:     for (j=0; j<A->cmap->n; j++) {
351:       PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&alpha,v+j*lda,&one));
352:     }
353:   } else {
354:     PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);
355:     PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&alpha,v,&one));
356:   }
357:   PetscLogFlops(nz);
358:   MatDenseRestoreArray(A,&v);
359:   return(0);
360: }

362: static PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool  *fl)
363: {
364:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
365:   PetscInt          i,j,m = A->rmap->n,N = a->lda;
366:   const PetscScalar *v;
367:   PetscErrorCode    ierr;

370:   *fl = PETSC_FALSE;
371:   if (A->rmap->n != A->cmap->n) return(0);
372:   MatDenseGetArrayRead(A,&v);
373:   for (i=0; i<m; i++) {
374:     for (j=i; j<m; j++) {
375:       if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) return(0);
376:     }
377:   }
378:   MatDenseRestoreArrayRead(A,&v);
379:   *fl  = PETSC_TRUE;
380:   return(0);
381: }

383: PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues)
384: {
385:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
387:   PetscInt       lda = (PetscInt)mat->lda,j,m;

390:   PetscLayoutReference(A->rmap,&newi->rmap);
391:   PetscLayoutReference(A->cmap,&newi->cmap);
392:   MatSeqDenseSetPreallocation(newi,NULL);
393:   if (cpvalues == MAT_COPY_VALUES) {
394:     const PetscScalar *av;
395:     PetscScalar       *v;

397:     MatDenseGetArrayRead(A,&av);
398:     MatDenseGetArray(newi,&v);
399:     if (lda>A->rmap->n) {
400:       m = A->rmap->n;
401:       for (j=0; j<A->cmap->n; j++) {
402:         PetscArraycpy(v+j*m,av+j*lda,m);
403:       }
404:     } else {
405:       PetscArraycpy(v,av,A->rmap->n*A->cmap->n);
406:     }
407:     MatDenseRestoreArray(newi,&v);
408:     MatDenseRestoreArrayRead(A,&av);
409:   }
410:   return(0);
411: }

413: PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat)
414: {

418:   MatCreate(PetscObjectComm((PetscObject)A),newmat);
419:   MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
420:   MatSetType(*newmat,((PetscObject)A)->type_name);
421:   MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);
422:   return(0);
423: }

425: static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
426: {
427:   MatFactorInfo  info;

431:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
432:   (*fact->ops->lufactor)(fact,0,0,&info);
433:   return(0);
434: }

436: static PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy)
437: {
438:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
439:   PetscErrorCode    ierr;
440:   const PetscScalar *x;
441:   PetscScalar       *y;
442:   PetscBLASInt      one = 1,info,m;

445:   PetscBLASIntCast(A->rmap->n,&m);
446:   VecGetArrayRead(xx,&x);
447:   VecGetArray(yy,&y);
448:   PetscArraycpy(y,x,A->rmap->n);
449:   if (A->factortype == MAT_FACTOR_LU) {
450: #if defined(PETSC_MISSING_LAPACK_GETRS)
451:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
452: #else
453:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
454:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
455:     PetscFPTrapPop();
456:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
457: #endif
458:   } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
459: #if defined(PETSC_MISSING_LAPACK_POTRS)
460:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
461: #else
462:     if (A->spd) {
463:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
464:       PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
465:       PetscFPTrapPop();
466:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
467: #if defined(PETSC_USE_COMPLEX)
468:     } else if (A->hermitian) {
469:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
470:       PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
471:       PetscFPTrapPop();
472:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve");
473: #endif
474:     } else { /* symmetric case */
475:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
476:       PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
477:       PetscFPTrapPop();
478:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
479:     }
480: #endif
481:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
482:   VecRestoreArrayRead(xx,&x);
483:   VecRestoreArray(yy,&y);
484:   PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
485:   return(0);
486: }

488: static PetscErrorCode MatMatSolve_SeqDense(Mat A,Mat B,Mat X)
489: {
490:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
491:   PetscErrorCode    ierr;
492:   const PetscScalar *b;
493:   PetscScalar       *x;
494:   PetscInt          n;
495:   PetscBLASInt      nrhs,info,m;

498:   PetscBLASIntCast(A->rmap->n,&m);
499:   MatGetSize(B,NULL,&n);
500:   PetscBLASIntCast(n,&nrhs);
501:   MatDenseGetArrayRead(B,&b);
502:   MatDenseGetArray(X,&x);

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

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

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

546: static PetscErrorCode MatConjugate_SeqDense(Mat);

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

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

607: /* ---------------------------------------------------------------*/
608: /* COMMENT: I have chosen to hide row permutation in the pivots,
609:    rather than put it in the Mat->row slot.*/
610: PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo)
611: {
612: #if defined(PETSC_MISSING_LAPACK_GETRF)
614:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable.");
615: #else
616:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
618:   PetscBLASInt   n,m,info;

621:   PetscBLASIntCast(A->cmap->n,&n);
622:   PetscBLASIntCast(A->rmap->n,&m);
623:   if (!mat->pivots) {
624:     PetscMalloc1(A->rmap->n,&mat->pivots);
625:     PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
626:   }
627:   if (!A->rmap->n || !A->cmap->n) return(0);
628:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
629:   PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info));
630:   PetscFPTrapPop();

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

635:   A->ops->solve             = MatSolve_SeqDense;
636:   A->ops->matsolve          = MatMatSolve_SeqDense;
637:   A->ops->solvetranspose    = MatSolveTranspose_SeqDense;
638:   A->factortype             = MAT_FACTOR_LU;

640:   PetscFree(A->solvertype);
641:   PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);

643:   PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);
644: #endif
645:   return(0);
646: }

648: /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */
649: PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo)
650: {
651: #if defined(PETSC_MISSING_LAPACK_POTRF)
653:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable.");
654: #else
655:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
657:   PetscBLASInt   info,n;

660:   PetscBLASIntCast(A->cmap->n,&n);
661:   if (!A->rmap->n || !A->cmap->n) return(0);
662:   if (A->spd) {
663:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
664:     PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info));
665:     PetscFPTrapPop();
666: #if defined(PETSC_USE_COMPLEX)
667:   } else if (A->hermitian) {
668:     if (!mat->pivots) {
669:       PetscMalloc1(A->rmap->n,&mat->pivots);
670:       PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
671:     }
672:     if (!mat->fwork) {
673:       PetscScalar dummy;

675:       mat->lfwork = -1;
676:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
677:       PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
678:       PetscFPTrapPop();
679:       mat->lfwork = (PetscInt)PetscRealPart(dummy);
680:       PetscMalloc1(mat->lfwork,&mat->fwork);
681:       PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
682:     }
683:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
684:     PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
685:     PetscFPTrapPop();
686: #endif
687:   } else { /* symmetric case */
688:     if (!mat->pivots) {
689:       PetscMalloc1(A->rmap->n,&mat->pivots);
690:       PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
691:     }
692:     if (!mat->fwork) {
693:       PetscScalar dummy;

695:       mat->lfwork = -1;
696:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
697:       PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
698:       PetscFPTrapPop();
699:       mat->lfwork = (PetscInt)PetscRealPart(dummy);
700:       PetscMalloc1(mat->lfwork,&mat->fwork);
701:       PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
702:     }
703:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
704:     PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
705:     PetscFPTrapPop();
706:   }
707:   if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1);

709:   A->ops->solve             = MatSolve_SeqDense;
710:   A->ops->matsolve          = MatMatSolve_SeqDense;
711:   A->ops->solvetranspose    = MatSolveTranspose_SeqDense;
712:   A->factortype             = MAT_FACTOR_CHOLESKY;

714:   PetscFree(A->solvertype);
715:   PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);

717:   PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
718: #endif
719:   return(0);
720: }


723: PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
724: {
726:   MatFactorInfo  info;

729:   info.fill = 1.0;

731:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
732:   (*fact->ops->choleskyfactor)(fact,0,&info);
733:   return(0);
734: }

736: PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
737: {
739:   fact->assembled                  = PETSC_TRUE;
740:   fact->preallocated               = PETSC_TRUE;
741:   fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
742:   fact->ops->solve                 = MatSolve_SeqDense;
743:   fact->ops->matsolve              = MatMatSolve_SeqDense;
744:   fact->ops->solvetranspose        = MatSolveTranspose_SeqDense;
745:   return(0);
746: }

748: PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info)
749: {
751:   fact->preallocated           = PETSC_TRUE;
752:   fact->assembled              = PETSC_TRUE;
753:   fact->ops->lufactornumeric   = MatLUFactorNumeric_SeqDense;
754:   fact->ops->solve             = MatSolve_SeqDense;
755:   fact->ops->matsolve          = MatMatSolve_SeqDense;
756:   fact->ops->solvetranspose    = MatSolveTranspose_SeqDense;
757:   return(0);
758: }

760: /* uses LAPACK */
761: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact)
762: {

766:   MatCreate(PetscObjectComm((PetscObject)A),fact);
767:   MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
768:   MatSetType(*fact,MATDENSE);
769:   if (ftype == MAT_FACTOR_LU) {
770:     (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
771:   } else {
772:     (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
773:   }
774:   (*fact)->factortype = ftype;

776:   PetscFree((*fact)->solvertype);
777:   PetscStrallocpy(MATSOLVERPETSC,&(*fact)->solvertype);
778:   return(0);
779: }

781: /* ------------------------------------------------------------------*/
782: static PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx)
783: {
784:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
785:   PetscScalar       *x,*v = mat->v,zero = 0.0,xt;
786:   const PetscScalar *b;
787:   PetscErrorCode    ierr;
788:   PetscInt          m = A->rmap->n,i;
789:   PetscBLASInt      o = 1,bm;

792: #if defined(PETSC_HAVE_CUDA)
793:   if (A->offloadmask == PETSC_OFFLOAD_GPU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented");
794: #endif
795:   if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */
796:   PetscBLASIntCast(m,&bm);
797:   if (flag & SOR_ZERO_INITIAL_GUESS) {
798:     /* this is a hack fix, should have another version without the second BLASdotu */
799:     VecSet(xx,zero);
800:   }
801:   VecGetArray(xx,&x);
802:   VecGetArrayRead(bb,&b);
803:   its  = its*lits;
804:   if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
805:   while (its--) {
806:     if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
807:       for (i=0; i<m; i++) {
808:         PetscStackCallBLAS("BLASdotu",xt   = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o));
809:         x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
810:       }
811:     }
812:     if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
813:       for (i=m-1; i>=0; i--) {
814:         PetscStackCallBLAS("BLASdotu",xt   = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o));
815:         x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
816:       }
817:     }
818:   }
819:   VecRestoreArrayRead(bb,&b);
820:   VecRestoreArray(xx,&x);
821:   return(0);
822: }

824: /* -----------------------------------------------------------------*/
825: PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)
826: {
827:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
828:   const PetscScalar *v   = mat->v,*x;
829:   PetscScalar       *y;
830:   PetscErrorCode    ierr;
831:   PetscBLASInt      m, n,_One=1;
832:   PetscScalar       _DOne=1.0,_DZero=0.0;

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

851: PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy)
852: {
853:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
854:   PetscScalar       *y,_DOne=1.0,_DZero=0.0;
855:   PetscErrorCode    ierr;
856:   PetscBLASInt      m, n, _One=1;
857:   const PetscScalar *v = mat->v,*x;

860:   PetscBLASIntCast(A->rmap->n,&m);
861:   PetscBLASIntCast(A->cmap->n,&n);
862:   VecGetArrayRead(xx,&x);
863:   VecGetArrayWrite(yy,&y);
864:   if (!A->rmap->n || !A->cmap->n) {
865:     PetscBLASInt i;
866:     for (i=0; i<m; i++) y[i] = 0.0;
867:   } else {
868:     PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One));
869:     PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);
870:   }
871:   VecRestoreArrayRead(xx,&x);
872:   VecRestoreArrayWrite(yy,&y);
873:   return(0);
874: }

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

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

898: PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
899: {
900:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
901:   const PetscScalar *v = mat->v,*x;
902:   PetscScalar       *y;
903:   PetscErrorCode    ierr;
904:   PetscBLASInt      m, n, _One=1;
905:   PetscScalar       _DOne=1.0;

908:   PetscBLASIntCast(A->rmap->n,&m);
909:   PetscBLASIntCast(A->cmap->n,&n);
910:   if (!A->rmap->n || !A->cmap->n) return(0);
911:   if (zz != yy) {VecCopy(zz,yy);}
912:   VecGetArrayRead(xx,&x);
913:   VecGetArray(yy,&y);
914:   PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
915:   VecRestoreArrayRead(xx,&x);
916:   VecRestoreArray(yy,&y);
917:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
918:   return(0);
919: }

921: /* -----------------------------------------------------------------*/
922: static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
923: {
924:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
926:   PetscInt       i;

929:   *ncols = A->cmap->n;
930:   if (cols) {
931:     PetscMalloc1(A->cmap->n+1,cols);
932:     for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
933:   }
934:   if (vals) {
935:     const PetscScalar *v;

937:     MatDenseGetArrayRead(A,&v);
938:     PetscMalloc1(A->cmap->n+1,vals);
939:     v   += row;
940:     for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
941:     MatDenseRestoreArrayRead(A,&v);
942:   }
943:   return(0);
944: }

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

951:   if (cols) {PetscFree(*cols);}
952:   if (vals) {PetscFree(*vals); }
953:   return(0);
954: }
955: /* ----------------------------------------------------------------*/
956: static PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
957: {
958:   Mat_SeqDense     *mat = (Mat_SeqDense*)A->data;
959:   PetscScalar      *av;
960:   PetscInt         i,j,idx=0;
961: #if defined(PETSC_HAVE_CUDA)
962:   PetscOffloadMask oldf;
963: #endif
964:   PetscErrorCode   ierr;

967:   MatDenseGetArray(A,&av);
968:   if (!mat->roworiented) {
969:     if (addv == INSERT_VALUES) {
970:       for (j=0; j<n; j++) {
971:         if (indexn[j] < 0) {idx += m; continue;}
972: #if defined(PETSC_USE_DEBUG)
973:         if (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);
974: #endif
975:         for (i=0; i<m; i++) {
976:           if (indexm[i] < 0) {idx++; continue;}
977: #if defined(PETSC_USE_DEBUG)
978:           if (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);
979: #endif
980:           av[indexn[j]*mat->lda + indexm[i]] = v[idx++];
981:         }
982:       }
983:     } else {
984:       for (j=0; j<n; j++) {
985:         if (indexn[j] < 0) {idx += m; continue;}
986: #if defined(PETSC_USE_DEBUG)
987:         if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1);
988: #endif
989:         for (i=0; i<m; i++) {
990:           if (indexm[i] < 0) {idx++; continue;}
991: #if defined(PETSC_USE_DEBUG)
992:           if (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);
993: #endif
994:           av[indexn[j]*mat->lda + indexm[i]] += v[idx++];
995:         }
996:       }
997:     }
998:   } else {
999:     if (addv == INSERT_VALUES) {
1000:       for (i=0; i<m; i++) {
1001:         if (indexm[i] < 0) { idx += n; continue;}
1002: #if defined(PETSC_USE_DEBUG)
1003:         if (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);
1004: #endif
1005:         for (j=0; j<n; j++) {
1006:           if (indexn[j] < 0) { idx++; continue;}
1007: #if defined(PETSC_USE_DEBUG)
1008:           if (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);
1009: #endif
1010:           av[indexn[j]*mat->lda + indexm[i]] = v[idx++];
1011:         }
1012:       }
1013:     } else {
1014:       for (i=0; i<m; i++) {
1015:         if (indexm[i] < 0) { idx += n; continue;}
1016: #if defined(PETSC_USE_DEBUG)
1017:         if (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);
1018: #endif
1019:         for (j=0; j<n; j++) {
1020:           if (indexn[j] < 0) { idx++; continue;}
1021: #if defined(PETSC_USE_DEBUG)
1022:           if (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);
1023: #endif
1024:           av[indexn[j]*mat->lda + indexm[i]] += v[idx++];
1025:         }
1026:       }
1027:     }
1028:   }
1029:   /* hack to prevent unneeded copy to the GPU while returning the array */
1030: #if defined(PETSC_HAVE_CUDA)
1031:   oldf = A->offloadmask;
1032:   A->offloadmask = PETSC_OFFLOAD_GPU;
1033: #endif
1034:   MatDenseRestoreArray(A,&av);
1035: #if defined(PETSC_HAVE_CUDA)
1036:   A->offloadmask = (oldf == PETSC_OFFLOAD_UNALLOCATED ? PETSC_OFFLOAD_UNALLOCATED : PETSC_OFFLOAD_CPU);
1037: #endif
1038:   return(0);
1039: }

1041: static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
1042: {
1043:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1044:   const PetscScalar *vv;
1045:   PetscInt          i,j;
1046:   PetscErrorCode    ierr;

1049:   MatDenseGetArrayRead(A,&vv);
1050:   /* row-oriented output */
1051:   for (i=0; i<m; i++) {
1052:     if (indexm[i] < 0) {v += n;continue;}
1053:     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);
1054:     for (j=0; j<n; j++) {
1055:       if (indexn[j] < 0) {v++; continue;}
1056:       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);
1057:       *v++ = vv[indexn[j]*mat->lda + indexm[i]];
1058:     }
1059:   }
1060:   MatDenseRestoreArrayRead(A,&vv);
1061:   return(0);
1062: }

1064: /* -----------------------------------------------------------------*/

1066: static PetscErrorCode MatLoad_SeqDense_Binary(Mat newmat,PetscViewer viewer)
1067: {
1068:   Mat_SeqDense   *a;
1070:   PetscInt       *scols,i,j,nz,header[4];
1071:   int            fd;
1072:   PetscMPIInt    size;
1073:   PetscInt       *rowlengths = 0,M,N,*cols,grows,gcols;
1074:   PetscScalar    *vals,*svals,*v,*w;
1075:   MPI_Comm       comm;

1078:   PetscObjectGetComm((PetscObject)viewer,&comm);
1079:   MPI_Comm_size(comm,&size);
1080:   if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor");
1081:   PetscViewerBinaryGetDescriptor(viewer,&fd);
1082:   PetscBinaryRead(fd,header,4,NULL,PETSC_INT);
1083:   if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not matrix object");
1084:   M = header[1]; N = header[2]; nz = header[3];

1086:   /* set global size if not set already*/
1087:   if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) {
1088:     MatSetSizes(newmat,M,N,M,N);
1089:   } else {
1090:     /* if sizes and type are already set, check if the vector global sizes are correct */
1091:     MatGetSize(newmat,&grows,&gcols);
1092:     if (M != grows ||  N != gcols) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different length (%d, %d) than the input matrix (%d, %d)",M,N,grows,gcols);
1093:   }
1094:   a = (Mat_SeqDense*)newmat->data;
1095:   if (!a->user_alloc) {
1096:     MatSeqDenseSetPreallocation(newmat,NULL);
1097:   }

1099:   if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */
1100:     a = (Mat_SeqDense*)newmat->data;
1101:     v = a->v;
1102:     /* Allocate some temp space to read in the values and then flip them
1103:        from row major to column major */
1104:     PetscMalloc1(M*N > 0 ? M*N : 1,&w);
1105:     /* read in nonzero values */
1106:     PetscBinaryRead(fd,w,M*N,NULL,PETSC_SCALAR);
1107:     /* now flip the values and store them in the matrix*/
1108:     for (j=0; j<N; j++) {
1109:       for (i=0; i<M; i++) {
1110:         *v++ =w[i*N+j];
1111:       }
1112:     }
1113:     PetscFree(w);
1114:   } else {
1115:     /* read row lengths */
1116:     PetscMalloc1(M+1,&rowlengths);
1117:     PetscBinaryRead(fd,rowlengths,M,NULL,PETSC_INT);

1119:     a = (Mat_SeqDense*)newmat->data;
1120:     v = a->v;

1122:     /* read column indices and nonzeros */
1123:     PetscMalloc1(nz+1,&scols);
1124:     cols = scols;
1125:     PetscBinaryRead(fd,cols,nz,NULL,PETSC_INT);
1126:     PetscMalloc1(nz+1,&svals);
1127:     vals = svals;
1128:     PetscBinaryRead(fd,vals,nz,NULL,PETSC_SCALAR);

1130:     /* insert into matrix */
1131:     for (i=0; i<M; i++) {
1132:       for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j];
1133:       svals += rowlengths[i]; scols += rowlengths[i];
1134:     }
1135:     PetscFree(vals);
1136:     PetscFree(cols);
1137:     PetscFree(rowlengths);
1138:   }
1139:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1140:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
1141:   return(0);
1142: }

1144: PetscErrorCode MatLoad_SeqDense(Mat newMat, PetscViewer viewer)
1145: {
1146:   PetscBool      isbinary, ishdf5;

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

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

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

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

1251: static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer)
1252: {
1253:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
1254:   PetscErrorCode    ierr;
1255:   int               fd;
1256:   PetscInt          ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n;
1257:   PetscScalar       *av,*v,*anonz,*vals;
1258:   PetscViewerFormat format;

1261:   PetscViewerBinaryGetDescriptor(viewer,&fd);
1262:   MatDenseGetArrayRead(A,(const PetscScalar**)&av);
1263:   PetscViewerGetFormat(viewer,&format);
1264:   if (format == PETSC_VIEWER_NATIVE) {
1265:     /* store the matrix as a dense matrix */
1266:     PetscMalloc1(4,&col_lens);

1268:     col_lens[0] = MAT_FILE_CLASSID;
1269:     col_lens[1] = m;
1270:     col_lens[2] = n;
1271:     col_lens[3] = MATRIX_BINARY_FORMAT_DENSE;

1273:     PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);
1274:     PetscFree(col_lens);

1276:     /* write out matrix, by rows */
1277:     PetscMalloc1(m*n+1,&vals);
1278:     v    = av;
1279:     for (j=0; j<n; j++) {
1280:       for (i=0; i<m; i++) {
1281:         vals[j + i*n] = *v++;
1282:       }
1283:     }
1284:     PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);
1285:     PetscFree(vals);
1286:   } else {
1287:     PetscMalloc1(4+nz,&col_lens);

1289:     col_lens[0] = MAT_FILE_CLASSID;
1290:     col_lens[1] = m;
1291:     col_lens[2] = n;
1292:     col_lens[3] = nz;

1294:     /* store lengths of each row and write (including header) to file */
1295:     for (i=0; i<m; i++) col_lens[4+i] = n;
1296:     PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);

1298:     /* Possibly should write in smaller increments, not whole matrix at once? */
1299:     /* store column indices (zero start index) */
1300:     ict = 0;
1301:     for (i=0; i<m; i++) {
1302:       for (j=0; j<n; j++) col_lens[ict++] = j;
1303:     }
1304:     PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);
1305:     PetscFree(col_lens);

1307:     /* store nonzero values */
1308:     PetscMalloc1(nz+1,&anonz);
1309:     ict  = 0;
1310:     for (i=0; i<m; i++) {
1311:       v = av + i;
1312:       for (j=0; j<n; j++) {
1313:         anonz[ict++] = *v; v += a->lda;
1314:       }
1315:     }
1316:     PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);
1317:     PetscFree(anonz);
1318:   }
1319:   MatDenseRestoreArrayRead(A,(const PetscScalar**)&av);
1320:   return(0);
1321: }

1323:  #include <petscdraw.h>
1324: static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa)
1325: {
1326:   Mat               A  = (Mat) Aa;
1327:   PetscErrorCode    ierr;
1328:   PetscInt          m  = A->rmap->n,n = A->cmap->n,i,j;
1329:   int               color = PETSC_DRAW_WHITE;
1330:   const PetscScalar *v;
1331:   PetscViewer       viewer;
1332:   PetscReal         xl,yl,xr,yr,x_l,x_r,y_l,y_r;
1333:   PetscViewerFormat format;

1336:   PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
1337:   PetscViewerGetFormat(viewer,&format);
1338:   PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);

1340:   /* Loop over matrix elements drawing boxes */
1341:   MatDenseGetArrayRead(A,&v);
1342:   if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1343:     PetscDrawCollectiveBegin(draw);
1344:     /* Blue for negative and Red for positive */
1345:     for (j = 0; j < n; j++) {
1346:       x_l = j; x_r = x_l + 1.0;
1347:       for (i = 0; i < m; i++) {
1348:         y_l = m - i - 1.0;
1349:         y_r = y_l + 1.0;
1350:         if (PetscRealPart(v[j*m+i]) >  0.) color = PETSC_DRAW_RED;
1351:         else if (PetscRealPart(v[j*m+i]) <  0.) color = PETSC_DRAW_BLUE;
1352:         else continue;
1353:         PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1354:       }
1355:     }
1356:     PetscDrawCollectiveEnd(draw);
1357:   } else {
1358:     /* use contour shading to indicate magnitude of values */
1359:     /* first determine max of all nonzero values */
1360:     PetscReal minv = 0.0, maxv = 0.0;
1361:     PetscDraw popup;

1363:     for (i=0; i < m*n; i++) {
1364:       if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1365:     }
1366:     if (minv >= maxv) maxv = minv + PETSC_SMALL;
1367:     PetscDrawGetPopup(draw,&popup);
1368:     PetscDrawScalePopup(popup,minv,maxv);

1370:     PetscDrawCollectiveBegin(draw);
1371:     for (j=0; j<n; j++) {
1372:       x_l = j;
1373:       x_r = x_l + 1.0;
1374:       for (i=0; i<m; i++) {
1375:         y_l   = m - i - 1.0;
1376:         y_r   = y_l + 1.0;
1377:         color = PetscDrawRealToColor(PetscAbsScalar(v[j*m+i]),minv,maxv);
1378:         PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1379:       }
1380:     }
1381:     PetscDrawCollectiveEnd(draw);
1382:   }
1383:   MatDenseRestoreArrayRead(A,&v);
1384:   return(0);
1385: }

1387: static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
1388: {
1389:   PetscDraw      draw;
1390:   PetscBool      isnull;
1391:   PetscReal      xr,yr,xl,yl,h,w;

1395:   PetscViewerDrawGetDraw(viewer,0,&draw);
1396:   PetscDrawIsNull(draw,&isnull);
1397:   if (isnull) return(0);

1399:   xr   = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
1400:   xr  += w;          yr += h;        xl = -w;     yl = -h;
1401:   PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
1402:   PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
1403:   PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
1404:   PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
1405:   PetscDrawSave(draw);
1406:   return(0);
1407: }

1409: PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
1410: {
1412:   PetscBool      iascii,isbinary,isdraw;

1415:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1416:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1417:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);

1419:   if (iascii) {
1420:     MatView_SeqDense_ASCII(A,viewer);
1421:   } else if (isbinary) {
1422:     MatView_SeqDense_Binary(A,viewer);
1423:   } else if (isdraw) {
1424:     MatView_SeqDense_Draw(A,viewer);
1425:   }
1426:   return(0);
1427: }

1429: static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A,const PetscScalar array[])
1430: {
1431:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

1434:   a->unplacedarray       = a->v;
1435:   a->unplaced_user_alloc = a->user_alloc;
1436:   a->v                   = (PetscScalar*) array;
1437: #if defined(PETSC_HAVE_CUDA)
1438:   A->offloadmask = PETSC_OFFLOAD_CPU;
1439: #endif
1440:   return(0);
1441: }

1443: static PetscErrorCode MatDenseResetArray_SeqDense(Mat A)
1444: {
1445:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

1448:   a->v             = a->unplacedarray;
1449:   a->user_alloc    = a->unplaced_user_alloc;
1450:   a->unplacedarray = NULL;
1451: #if defined(PETSC_HAVE_CUDA)
1452:   A->offloadmask = PETSC_OFFLOAD_CPU;
1453: #endif
1454:   return(0);
1455: }

1457: PetscErrorCode MatDestroy_SeqDense(Mat mat)
1458: {
1459:   Mat_SeqDense   *l = (Mat_SeqDense*)mat->data;

1463: #if defined(PETSC_USE_LOG)
1464:   PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n);
1465: #endif
1466:   PetscFree(l->pivots);
1467:   PetscFree(l->fwork);
1468:   MatDestroy(&l->ptapwork);
1469:   if (!l->user_alloc) {PetscFree(l->v);}
1470:   PetscFree(mat->data);

1472:   PetscObjectChangeTypeName((PetscObject)mat,0);
1473:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetLDA_C",NULL);
1474:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);
1475:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);
1476:   PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL);
1477:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL);
1478:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayRead_C",NULL);
1479:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayRead_C",NULL);
1480:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);
1481: #if defined(PETSC_HAVE_ELEMENTAL)
1482:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);
1483: #endif
1484: #if defined(PETSC_HAVE_CUDA)
1485:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqdensecuda_C",NULL);
1486:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaijcusparse_seqdense_C",NULL);
1487: #endif
1488: #if defined(PETSC_HAVE_VIENNACL)
1489:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaijviennacl_seqdense_C",NULL);
1490: #endif
1491:   PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);
1492:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);
1493:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);
1494:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);
1495:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqbaij_seqdense_C",NULL);
1496:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqbaij_seqdense_C",NULL);
1497:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqbaij_seqdense_C",NULL);
1498:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqsbaij_seqdense_C",NULL);
1499:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqsbaij_seqdense_C",NULL);
1500:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqsbaij_seqdense_C",NULL);
1501:   PetscObjectComposeFunction((PetscObject)mat,"MatPtAP_seqaij_seqdense_C",NULL);
1502:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaijperm_seqdense_C",NULL);
1503:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaijperm_seqdense_C",NULL);
1504:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaijperm_seqdense_C",NULL);
1505:   PetscObjectComposeFunction((PetscObject)mat,"MatPtAP_seqaijperm_seqdense_C",NULL);
1506:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaijsell_seqdense_C",NULL);
1507:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaijsell_seqdense_C",NULL);
1508:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaijsell_seqdense_C",NULL);
1509:   PetscObjectComposeFunction((PetscObject)mat,"MatPtAP_seqaijsell_seqdense_C",NULL);
1510:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaijmkl_seqdense_C",NULL);
1511:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaijmkl_seqdense_C",NULL);
1512:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaijmkl_seqdense_C",NULL);
1513:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_nest_seqdense_C",NULL);
1514:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_nest_seqdense_C",NULL);
1515:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_nest_seqdense_C",NULL);
1516:   PetscObjectComposeFunction((PetscObject)mat,"MatPtAP_seqaijmkl_seqdense_C",NULL);

1518:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);
1519:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);
1520:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);
1521:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaijperm_seqdense_C",NULL);
1522:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaijperm_seqdense_C",NULL);
1523:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaijperm_seqdense_C",NULL);
1524:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaijsell_seqdense_C",NULL);
1525:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaijsell_seqdense_C",NULL);
1526:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaijsell_seqdense_C",NULL);

1528:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaijmkl_seqdense_C",NULL);
1529:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaijmkl_seqdense_C",NULL);
1530:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaijmkl_seqdense_C",NULL);
1531:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumn_C",NULL);
1532:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumn_C",NULL);
1533:   return(0);
1534: }

1536: static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
1537: {
1538:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1540:   PetscInt       k,j,m,n,M;
1541:   PetscScalar    *v,tmp;

1544:   m = A->rmap->n; M = mat->lda; n = A->cmap->n;
1545:   if (reuse == MAT_INPLACE_MATRIX && m == n) { /* in place transpose */
1546:     MatDenseGetArray(A,&v);
1547:     for (j=0; j<m; j++) {
1548:       for (k=0; k<j; k++) {
1549:         tmp        = v[j + k*M];
1550:         v[j + k*M] = v[k + j*M];
1551:         v[k + j*M] = tmp;
1552:       }
1553:     }
1554:     MatDenseRestoreArray(A,&v);
1555:   } else { /* out-of-place transpose */
1556:     Mat          tmat;
1557:     Mat_SeqDense *tmatd;
1558:     PetscScalar  *v2;
1559:     PetscInt     M2;

1561:     if (reuse != MAT_REUSE_MATRIX) {
1562:       MatCreate(PetscObjectComm((PetscObject)A),&tmat);
1563:       MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);
1564:       MatSetType(tmat,((PetscObject)A)->type_name);
1565:       MatSeqDenseSetPreallocation(tmat,NULL);
1566:     } else tmat = *matout;

1568:     MatDenseGetArrayRead(A,(const PetscScalar**)&v);
1569:     MatDenseGetArray(tmat,&v2);
1570:     tmatd = (Mat_SeqDense*)tmat->data;
1571:     M2    = tmatd->lda;
1572:     for (j=0; j<n; j++) {
1573:       for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M];
1574:     }
1575:     MatDenseRestoreArray(tmat,&v2);
1576:     MatDenseRestoreArrayRead(A,(const PetscScalar**)&v);
1577:     MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
1578:     MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
1579:     if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) *matout = tmat;
1580:     else {
1581:       MatHeaderMerge(A,&tmat);
1582:     }
1583:   }
1584:   return(0);
1585: }

1587: static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool  *flg)
1588: {
1589:   Mat_SeqDense      *mat1 = (Mat_SeqDense*)A1->data;
1590:   Mat_SeqDense      *mat2 = (Mat_SeqDense*)A2->data;
1591:   PetscInt          i;
1592:   const PetscScalar *v1,*v2;
1593:   PetscErrorCode    ierr;

1596:   if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; return(0);}
1597:   if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; return(0);}
1598:   MatDenseGetArrayRead(A1,&v1);
1599:   MatDenseGetArrayRead(A2,&v2);
1600:   for (i=0; i<A1->cmap->n; i++) {
1601:     PetscArraycmp(v1,v2,A1->rmap->n,flg);
1602:     if (*flg == PETSC_FALSE) return(0);
1603:     v1 += mat1->lda;
1604:     v2 += mat2->lda;
1605:   }
1606:   MatDenseRestoreArrayRead(A1,&v1);
1607:   MatDenseRestoreArrayRead(A2,&v2);
1608:   *flg = PETSC_TRUE;
1609:   return(0);
1610: }

1612: static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
1613: {
1614:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1615:   PetscInt          i,n,len;
1616:   PetscScalar       *x;
1617:   const PetscScalar *vv;
1618:   PetscErrorCode    ierr;

1621:   VecGetSize(v,&n);
1622:   VecGetArray(v,&x);
1623:   len  = PetscMin(A->rmap->n,A->cmap->n);
1624:   MatDenseGetArrayRead(A,&vv);
1625:   if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
1626:   for (i=0; i<len; i++) {
1627:     x[i] = vv[i*mat->lda + i];
1628:   }
1629:   MatDenseRestoreArrayRead(A,&vv);
1630:   VecRestoreArray(v,&x);
1631:   return(0);
1632: }

1634: static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
1635: {
1636:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1637:   const PetscScalar *l,*r;
1638:   PetscScalar       x,*v,*vv;
1639:   PetscErrorCode    ierr;
1640:   PetscInt          i,j,m = A->rmap->n,n = A->cmap->n;

1643:   MatDenseGetArray(A,&vv);
1644:   if (ll) {
1645:     VecGetSize(ll,&m);
1646:     VecGetArrayRead(ll,&l);
1647:     if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size");
1648:     for (i=0; i<m; i++) {
1649:       x = l[i];
1650:       v = vv + i;
1651:       for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;}
1652:     }
1653:     VecRestoreArrayRead(ll,&l);
1654:     PetscLogFlops(1.0*n*m);
1655:   }
1656:   if (rr) {
1657:     VecGetSize(rr,&n);
1658:     VecGetArrayRead(rr,&r);
1659:     if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size");
1660:     for (i=0; i<n; i++) {
1661:       x = r[i];
1662:       v = vv + i*mat->lda;
1663:       for (j=0; j<m; j++) (*v++) *= x;
1664:     }
1665:     VecRestoreArrayRead(rr,&r);
1666:     PetscLogFlops(1.0*n*m);
1667:   }
1668:   MatDenseRestoreArray(A,&vv);
1669:   return(0);
1670: }

1672: PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1673: {
1674:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1675:   PetscScalar       *v,*vv;
1676:   PetscReal         sum  = 0.0;
1677:   PetscInt          lda  =mat->lda,m=A->rmap->n,i,j;
1678:   PetscErrorCode    ierr;

1681:   MatDenseGetArrayRead(A,(const PetscScalar**)&vv);
1682:   v    = vv;
1683:   if (type == NORM_FROBENIUS) {
1684:     if (lda>m) {
1685:       for (j=0; j<A->cmap->n; j++) {
1686:         v = vv+j*lda;
1687:         for (i=0; i<m; i++) {
1688:           sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1689:         }
1690:       }
1691:     } else {
1692: #if defined(PETSC_USE_REAL___FP16)
1693:       PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n;
1694:       *nrm = BLASnrm2_(&cnt,v,&one);
1695:     }
1696: #else
1697:       for (i=0; i<A->cmap->n*A->rmap->n; i++) {
1698:         sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1699:       }
1700:     }
1701:     *nrm = PetscSqrtReal(sum);
1702: #endif
1703:     PetscLogFlops(2.0*A->cmap->n*A->rmap->n);
1704:   } else if (type == NORM_1) {
1705:     *nrm = 0.0;
1706:     for (j=0; j<A->cmap->n; j++) {
1707:       v   = vv + j*mat->lda;
1708:       sum = 0.0;
1709:       for (i=0; i<A->rmap->n; i++) {
1710:         sum += PetscAbsScalar(*v);  v++;
1711:       }
1712:       if (sum > *nrm) *nrm = sum;
1713:     }
1714:     PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1715:   } else if (type == NORM_INFINITY) {
1716:     *nrm = 0.0;
1717:     for (j=0; j<A->rmap->n; j++) {
1718:       v   = vv + j;
1719:       sum = 0.0;
1720:       for (i=0; i<A->cmap->n; i++) {
1721:         sum += PetscAbsScalar(*v); v += mat->lda;
1722:       }
1723:       if (sum > *nrm) *nrm = sum;
1724:     }
1725:     PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1726:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
1727:   MatDenseRestoreArrayRead(A,(const PetscScalar**)&vv);
1728:   return(0);
1729: }

1731: static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
1732: {
1733:   Mat_SeqDense   *aij = (Mat_SeqDense*)A->data;

1737:   switch (op) {
1738:   case MAT_ROW_ORIENTED:
1739:     aij->roworiented = flg;
1740:     break;
1741:   case MAT_NEW_NONZERO_LOCATIONS:
1742:   case MAT_NEW_NONZERO_LOCATION_ERR:
1743:   case MAT_NEW_NONZERO_ALLOCATION_ERR:
1744:   case MAT_NEW_DIAGONALS:
1745:   case MAT_KEEP_NONZERO_PATTERN:
1746:   case MAT_IGNORE_OFF_PROC_ENTRIES:
1747:   case MAT_USE_HASH_TABLE:
1748:   case MAT_IGNORE_ZERO_ENTRIES:
1749:   case MAT_IGNORE_LOWER_TRIANGULAR:
1750:   case MAT_SORTED_FULL:
1751:     PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
1752:     break;
1753:   case MAT_SPD:
1754:   case MAT_SYMMETRIC:
1755:   case MAT_STRUCTURALLY_SYMMETRIC:
1756:   case MAT_HERMITIAN:
1757:   case MAT_SYMMETRY_ETERNAL:
1758:     /* These options are handled directly by MatSetOption() */
1759:     break;
1760:   default:
1761:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]);
1762:   }
1763:   return(0);
1764: }

1766: static PetscErrorCode MatZeroEntries_SeqDense(Mat A)
1767: {
1768:   Mat_SeqDense   *l = (Mat_SeqDense*)A->data;
1770:   PetscInt       lda=l->lda,m=A->rmap->n,j;
1771:   PetscScalar    *v;

1774:   MatDenseGetArray(A,&v);
1775:   if (lda>m) {
1776:     for (j=0; j<A->cmap->n; j++) {
1777:       PetscArrayzero(v+j*lda,m);
1778:     }
1779:   } else {
1780:     PetscArrayzero(v,A->rmap->n*A->cmap->n);
1781:   }
1782:   MatDenseRestoreArray(A,&v);
1783:   return(0);
1784: }

1786: static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1787: {
1788:   PetscErrorCode    ierr;
1789:   Mat_SeqDense      *l = (Mat_SeqDense*)A->data;
1790:   PetscInt          m  = l->lda, n = A->cmap->n, i,j;
1791:   PetscScalar       *slot,*bb,*v;
1792:   const PetscScalar *xx;

1795: #if defined(PETSC_USE_DEBUG)
1796:   for (i=0; i<N; i++) {
1797:     if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
1798:     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);
1799:   }
1800: #endif
1801:   if (!N) return(0);

1803:   /* fix right hand side if needed */
1804:   if (x && b) {
1805:     VecGetArrayRead(x,&xx);
1806:     VecGetArray(b,&bb);
1807:     for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
1808:     VecRestoreArrayRead(x,&xx);
1809:     VecRestoreArray(b,&bb);
1810:   }

1812:   MatDenseGetArray(A,&v);
1813:   for (i=0; i<N; i++) {
1814:     slot = v + rows[i];
1815:     for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
1816:   }
1817:   if (diag != 0.0) {
1818:     if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
1819:     for (i=0; i<N; i++) {
1820:       slot  = v + (m+1)*rows[i];
1821:       *slot = diag;
1822:     }
1823:   }
1824:   MatDenseRestoreArray(A,&v);
1825:   return(0);
1826: }

1828: static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A,PetscInt *lda)
1829: {
1830:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;

1833:   *lda = mat->lda;
1834:   return(0);
1835: }

1837: PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[])
1838: {
1839:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;

1842:   *array = mat->v;
1843:   return(0);
1844: }

1846: PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[])
1847: {
1849:   return(0);
1850: }

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

1855:    Logically Collective on Mat

1857:    Input Parameter:
1858: .  mat - a MATSEQDENSE or MATMPIDENSE matrix

1860:    Output Parameter:
1861: .   lda - the leading dimension

1863:    Level: intermediate

1865: .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatSeqDenseSetLDA()
1866: @*/
1867: PetscErrorCode  MatDenseGetLDA(Mat A,PetscInt *lda)
1868: {

1872:   PetscUseMethod(A,"MatDenseGetLDA_C",(Mat,PetscInt*),(A,lda));
1873:   return(0);
1874: }

1876: /*@C
1877:    MatDenseGetArray - gives access to the array where the data for a SeqDense matrix is stored

1879:    Logically Collective on Mat

1881:    Input Parameter:
1882: .  mat - a MATSEQDENSE or MATMPIDENSE matrix

1884:    Output Parameter:
1885: .   array - pointer to the data

1887:    Level: intermediate

1889: .seealso: MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead()
1890: @*/
1891: PetscErrorCode  MatDenseGetArray(Mat A,PetscScalar **array)
1892: {

1896:   PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));
1897:   return(0);
1898: }

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

1903:    Logically Collective on Mat

1905:    Input Parameters:
1906: +  mat - a MATSEQDENSE or MATMPIDENSE matrix
1907: -  array - pointer to the data

1909:    Level: intermediate

1911: .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead()
1912: @*/
1913: PetscErrorCode  MatDenseRestoreArray(Mat A,PetscScalar **array)
1914: {

1918:   PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));
1919:   if (array) *array = NULL;
1920:   PetscObjectStateIncrease((PetscObject)A);
1921:   return(0);
1922: }

1924: /*@C
1925:    MatDenseGetArrayRead - gives access to the array where the data for a SeqDense matrix is stored

1927:    Not Collective

1929:    Input Parameter:
1930: .  mat - a MATSEQDENSE or MATMPIDENSE matrix

1932:    Output Parameter:
1933: .   array - pointer to the data

1935:    Level: intermediate

1937: .seealso: MatDenseRestoreArray(), MatDenseGetArray(), MatDenseRestoreArrayRead()
1938: @*/
1939: PetscErrorCode  MatDenseGetArrayRead(Mat A,const PetscScalar **array)
1940: {

1944:   PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array));
1945:   return(0);
1946: }

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

1951:    Not Collective

1953:    Input Parameters:
1954: +  mat - a MATSEQDENSE or MATMPIDENSE matrix
1955: -  array - pointer to the data

1957:    Level: intermediate

1959: .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArray()
1960: @*/
1961: PetscErrorCode  MatDenseRestoreArrayRead(Mat A,const PetscScalar **array)
1962: {

1966:   PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array));
1967:   if (array) *array = NULL;
1968:   return(0);
1969: }

1971: static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B)
1972: {
1973:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1975:   PetscInt       i,j,nrows,ncols,blda;
1976:   const PetscInt *irow,*icol;
1977:   PetscScalar    *av,*bv,*v = mat->v;
1978:   Mat            newmat;

1981:   ISGetIndices(isrow,&irow);
1982:   ISGetIndices(iscol,&icol);
1983:   ISGetLocalSize(isrow,&nrows);
1984:   ISGetLocalSize(iscol,&ncols);

1986:   /* Check submatrixcall */
1987:   if (scall == MAT_REUSE_MATRIX) {
1988:     PetscInt n_cols,n_rows;
1989:     MatGetSize(*B,&n_rows,&n_cols);
1990:     if (n_rows != nrows || n_cols != ncols) {
1991:       /* resize the result matrix to match number of requested rows/columns */
1992:       MatSetSizes(*B,nrows,ncols,nrows,ncols);
1993:     }
1994:     newmat = *B;
1995:   } else {
1996:     /* Create and fill new matrix */
1997:     MatCreate(PetscObjectComm((PetscObject)A),&newmat);
1998:     MatSetSizes(newmat,nrows,ncols,nrows,ncols);
1999:     MatSetType(newmat,((PetscObject)A)->type_name);
2000:     MatSeqDenseSetPreallocation(newmat,NULL);
2001:   }

2003:   /* Now extract the data pointers and do the copy,column at a time */
2004:   MatDenseGetArray(newmat,&bv);
2005:   MatDenseGetLDA(newmat,&blda);
2006:   for (i=0; i<ncols; i++) {
2007:     av = v + mat->lda*icol[i];
2008:     for (j=0; j<nrows; j++) bv[j] = av[irow[j]];
2009:     bv += blda;
2010:   }
2011:   MatDenseRestoreArray(newmat,&bv);

2013:   /* Assemble the matrices so that the correct flags are set */
2014:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
2015:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);

2017:   /* Free work space */
2018:   ISRestoreIndices(isrow,&irow);
2019:   ISRestoreIndices(iscol,&icol);
2020:   *B   = newmat;
2021:   return(0);
2022: }

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

2030:   if (scall == MAT_INITIAL_MATRIX) {
2031:     PetscCalloc1(n+1,B);
2032:   }

2034:   for (i=0; i<n; i++) {
2035:     MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
2036:   }
2037:   return(0);
2038: }

2040: static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
2041: {
2043:   return(0);
2044: }

2046: static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
2047: {
2049:   return(0);
2050: }

2052: static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
2053: {
2054:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data;
2055:   PetscErrorCode    ierr;
2056:   const PetscScalar *va;
2057:   PetscScalar       *vb;
2058:   PetscInt          lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j;

2061:   /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
2062:   if (A->ops->copy != B->ops->copy) {
2063:     MatCopy_Basic(A,B,str);
2064:     return(0);
2065:   }
2066:   if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)");
2067:   MatDenseGetArrayRead(A,&va);
2068:   MatDenseGetArray(B,&vb);
2069:   if (lda1>m || lda2>m) {
2070:     for (j=0; j<n; j++) {
2071:       PetscArraycpy(vb+j*lda2,va+j*lda1,m);
2072:     }
2073:   } else {
2074:     PetscArraycpy(vb,va,A->rmap->n*A->cmap->n);
2075:   }
2076:   MatDenseRestoreArray(B,&vb);
2077:   MatDenseRestoreArrayRead(A,&va);
2078:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
2079:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
2080:   return(0);
2081: }

2083: static PetscErrorCode MatSetUp_SeqDense(Mat A)
2084: {

2088:   MatSeqDenseSetPreallocation(A,0);
2089:   return(0);
2090: }

2092: static PetscErrorCode MatConjugate_SeqDense(Mat A)
2093: {
2094:   PetscInt       i,nz = A->rmap->n*A->cmap->n;
2095:   PetscScalar    *aa;

2099:   MatDenseGetArray(A,&aa);
2100:   for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
2101:   MatDenseRestoreArray(A,&aa);
2102:   return(0);
2103: }

2105: static PetscErrorCode MatRealPart_SeqDense(Mat A)
2106: {
2107:   PetscInt       i,nz = A->rmap->n*A->cmap->n;
2108:   PetscScalar    *aa;

2112:   MatDenseGetArray(A,&aa);
2113:   for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
2114:   MatDenseRestoreArray(A,&aa);
2115:   return(0);
2116: }

2118: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
2119: {
2120:   PetscInt       i,nz = A->rmap->n*A->cmap->n;
2121:   PetscScalar    *aa;

2125:   MatDenseGetArray(A,&aa);
2126:   for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
2127:   MatDenseRestoreArray(A,&aa);
2128:   return(0);
2129: }

2131: /* ----------------------------------------------------------------*/
2132: PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
2133: {

2137:   if (scall == MAT_INITIAL_MATRIX) {
2138:     PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
2139:     MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
2140:     PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
2141:   }
2142:   PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
2143:   if ((*C)->ops->matmultnumeric) {
2144:     (*(*C)->ops->matmultnumeric)(A,B,*C);
2145:   } else {
2146:     MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);
2147:   }
2148:   PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
2149:   return(0);
2150: }

2152: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
2153: {
2155:   PetscInt       m=A->rmap->n,n=B->cmap->n;
2156:   Mat            Cmat;
2157:   PetscBool      flg;

2160:   MatCreate(PETSC_COMM_SELF,&Cmat);
2161:   MatSetSizes(Cmat,m,n,m,n);
2162:   PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2163:   MatSetType(Cmat,flg ? ((PetscObject)A)->type_name : MATDENSE);
2164:   MatSeqDenseSetPreallocation(Cmat,NULL);
2165:   *C   = Cmat;
2166:   return(0);
2167: }

2169: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2170: {
2171:   Mat_SeqDense       *a,*b=(Mat_SeqDense*)B->data,*c=(Mat_SeqDense*)C->data;
2172:   PetscBLASInt       m,n,k;
2173:   const PetscScalar *av,*bv;
2174:   PetscScalar       *cv;
2175:   PetscScalar       _DOne=1.0,_DZero=0.0;
2176:   PetscBool         flg;
2177:   PetscErrorCode    (*numeric)(Mat,Mat,Mat)=NULL;
2178:   PetscErrorCode    ierr;

2181:   /* Handle case where where user provided the final C matrix rather than calling MatMatMult() with MAT_INITIAL_MATRIX*/
2182:   PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
2183:   if (flg) numeric = MatMatMultNumeric_SeqAIJ_SeqDense;
2184:   PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&flg);
2185:   if (flg) numeric = MatMatMultNumeric_SeqBAIJ_SeqDense;
2186:   PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&flg);
2187:   if (flg) numeric = MatMatMultNumeric_SeqSBAIJ_SeqDense;
2188:   PetscObjectTypeCompare((PetscObject)A,MATNEST,&flg);
2189:   if (flg) numeric = MatMatMultNumeric_Nest_Dense;
2190:   if (numeric) {
2191:     C->ops->matmultnumeric = numeric;
2192:     (*numeric)(A,B,C);
2193:     return(0);
2194:   }
2195:   a = (Mat_SeqDense*)A->data;
2196:   PetscBLASIntCast(C->rmap->n,&m);
2197:   PetscBLASIntCast(C->cmap->n,&n);
2198:   PetscBLASIntCast(A->cmap->n,&k);
2199:   if (!m || !n || !k) return(0);
2200:   MatDenseGetArrayRead(A,&av);
2201:   MatDenseGetArrayRead(B,&bv);
2202:   MatDenseGetArray(C,&cv);
2203:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2204:   PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2205:   MatDenseRestoreArrayRead(A,&av);
2206:   MatDenseRestoreArrayRead(B,&bv);
2207:   MatDenseRestoreArray(C,&cv);
2208:   return(0);
2209: }

2211: PetscErrorCode MatMatTransposeMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
2212: {

2216:   if (scall == MAT_INITIAL_MATRIX) {
2217:     PetscLogEventBegin(MAT_MatTransposeMultSymbolic,A,B,0,0);
2218:     MatMatTransposeMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
2219:     PetscLogEventEnd(MAT_MatTransposeMultSymbolic,A,B,0,0);
2220:   }
2221:   PetscLogEventBegin(MAT_MatTransposeMultNumeric,A,B,0,0);
2222:   MatMatTransposeMultNumeric_SeqDense_SeqDense(A,B,*C);
2223:   PetscLogEventEnd(MAT_MatTransposeMultNumeric,A,B,0,0);
2224:   return(0);
2225: }

2227: PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
2228: {
2230:   PetscInt       m=A->rmap->n,n=B->rmap->n;
2231:   Mat            Cmat;
2232:   PetscBool      flg;

2235:   MatCreate(PETSC_COMM_SELF,&Cmat);
2236:   MatSetSizes(Cmat,m,n,m,n);
2237:   PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2238:   MatSetType(Cmat,flg ? ((PetscObject)A)->type_name : MATDENSE);
2239:   MatSeqDenseSetPreallocation(Cmat,NULL);
2240:   *C   = Cmat;
2241:   return(0);
2242: }

2244: PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2245: {
2246:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2247:   Mat_SeqDense   *b = (Mat_SeqDense*)B->data;
2248:   Mat_SeqDense   *c = (Mat_SeqDense*)C->data;
2249:   PetscBLASInt   m,n,k;
2250:   PetscScalar    _DOne=1.0,_DZero=0.0;

2254:   PetscBLASIntCast(C->rmap->n,&m);
2255:   PetscBLASIntCast(C->cmap->n,&n);
2256:   PetscBLASIntCast(A->cmap->n,&k);
2257:   if (!m || !n || !k) return(0);
2258:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
2259:   PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2260:   return(0);
2261: }

2263: PetscErrorCode MatTransposeMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
2264: {

2268:   if (scall == MAT_INITIAL_MATRIX) {
2269:     PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);
2270:     MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
2271:     PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);
2272:   }
2273:   PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);
2274:   MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);
2275:   PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);
2276:   return(0);
2277: }

2279: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
2280: {
2282:   PetscInt       m=A->cmap->n,n=B->cmap->n;
2283:   Mat            Cmat;
2284:   PetscBool      flg;

2287:   MatCreate(PETSC_COMM_SELF,&Cmat);
2288:   MatSetSizes(Cmat,m,n,m,n);
2289:   PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2290:   MatSetType(Cmat,flg ? ((PetscObject)A)->type_name : MATDENSE);
2291:   MatSeqDenseSetPreallocation(Cmat,NULL);
2292:   *C   = Cmat;
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:   PetscBLASInt   m,n,k;
2302:   PetscScalar    _DOne=1.0,_DZero=0.0;

2306:   PetscBLASIntCast(C->rmap->n,&m);
2307:   PetscBLASIntCast(C->cmap->n,&n);
2308:   PetscBLASIntCast(A->rmap->n,&k);
2309:   if (!m || !n || !k) return(0);
2310:   PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
2311:   PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2312:   return(0);
2313: }

2315: static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[])
2316: {
2317:   Mat_SeqDense       *a = (Mat_SeqDense*)A->data;
2318:   PetscErrorCode     ierr;
2319:   PetscInt           i,j,m = A->rmap->n,n = A->cmap->n,p;
2320:   PetscScalar        *x;
2321:   const PetscScalar *aa;

2324:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2325:   VecGetArray(v,&x);
2326:   VecGetLocalSize(v,&p);
2327:   MatDenseGetArrayRead(A,&aa);
2328:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2329:   for (i=0; i<m; i++) {
2330:     x[i] = aa[i]; if (idx) idx[i] = 0;
2331:     for (j=1; j<n; j++) {
2332:       if (PetscRealPart(x[i]) < PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;}
2333:     }
2334:   }
2335:   MatDenseRestoreArrayRead(A,&aa);
2336:   VecRestoreArray(v,&x);
2337:   return(0);
2338: }

2340: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])
2341: {
2342:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
2343:   PetscErrorCode    ierr;
2344:   PetscInt          i,j,m = A->rmap->n,n = A->cmap->n,p;
2345:   PetscScalar       *x;
2346:   PetscReal         atmp;
2347:   const PetscScalar *aa;

2350:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2351:   VecGetArray(v,&x);
2352:   VecGetLocalSize(v,&p);
2353:   MatDenseGetArrayRead(A,&aa);
2354:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2355:   for (i=0; i<m; i++) {
2356:     x[i] = PetscAbsScalar(aa[i]);
2357:     for (j=1; j<n; j++) {
2358:       atmp = PetscAbsScalar(aa[i+a->lda*j]);
2359:       if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;}
2360:     }
2361:   }
2362:   MatDenseRestoreArrayRead(A,&aa);
2363:   VecRestoreArray(v,&x);
2364:   return(0);
2365: }

2367: static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])
2368: {
2369:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
2370:   PetscErrorCode    ierr;
2371:   PetscInt          i,j,m = A->rmap->n,n = A->cmap->n,p;
2372:   PetscScalar       *x;
2373:   const PetscScalar *aa;

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

2392: static PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
2393: {
2394:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
2395:   PetscErrorCode    ierr;
2396:   PetscScalar       *x;
2397:   const PetscScalar *aa;

2400:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2401:   MatDenseGetArrayRead(A,&aa);
2402:   VecGetArray(v,&x);
2403:   PetscArraycpy(x,aa+col*a->lda,A->rmap->n);
2404:   VecRestoreArray(v,&x);
2405:   MatDenseRestoreArrayRead(A,&aa);
2406:   return(0);
2407: }

2409: PETSC_INTERN PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms)
2410: {
2411:   PetscErrorCode    ierr;
2412:   PetscInt          i,j,m,n;
2413:   const PetscScalar *a;

2416:   MatGetSize(A,&m,&n);
2417:   PetscArrayzero(norms,n);
2418:   MatDenseGetArrayRead(A,&a);
2419:   if (type == NORM_2) {
2420:     for (i=0; i<n; i++) {
2421:       for (j=0; j<m; j++) {
2422:         norms[i] += PetscAbsScalar(a[j]*a[j]);
2423:       }
2424:       a += m;
2425:     }
2426:   } else if (type == NORM_1) {
2427:     for (i=0; i<n; i++) {
2428:       for (j=0; j<m; j++) {
2429:         norms[i] += PetscAbsScalar(a[j]);
2430:       }
2431:       a += m;
2432:     }
2433:   } else if (type == NORM_INFINITY) {
2434:     for (i=0; i<n; i++) {
2435:       for (j=0; j<m; j++) {
2436:         norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]);
2437:       }
2438:       a += m;
2439:     }
2440:   } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType");
2441:   MatDenseRestoreArrayRead(A,&a);
2442:   if (type == NORM_2) {
2443:     for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]);
2444:   }
2445:   return(0);
2446: }

2448: static PetscErrorCode  MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
2449: {
2451:   PetscScalar    *a;
2452:   PetscInt       m,n,i;

2455:   MatGetSize(x,&m,&n);
2456:   MatDenseGetArray(x,&a);
2457:   for (i=0; i<m*n; i++) {
2458:     PetscRandomGetValue(rctx,a+i);
2459:   }
2460:   MatDenseRestoreArray(x,&a);
2461:   return(0);
2462: }

2464: static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool  *missing,PetscInt *d)
2465: {
2467:   *missing = PETSC_FALSE;
2468:   return(0);
2469: }

2471: /* vals is not const */
2472: static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals)
2473: {
2475:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2476:   PetscScalar    *v;

2479:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2480:   MatDenseGetArray(A,&v);
2481:   *vals = v+col*a->lda;
2482:   MatDenseRestoreArray(A,&v);
2483:   return(0);
2484: }

2486: static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals)
2487: {
2489:   *vals = 0; /* user cannot accidently use the array later */
2490:   return(0);
2491: }

2493: /* -------------------------------------------------------------------*/
2494: static struct _MatOps MatOps_Values = { MatSetValues_SeqDense,
2495:                                         MatGetRow_SeqDense,
2496:                                         MatRestoreRow_SeqDense,
2497:                                         MatMult_SeqDense,
2498:                                 /*  4*/ MatMultAdd_SeqDense,
2499:                                         MatMultTranspose_SeqDense,
2500:                                         MatMultTransposeAdd_SeqDense,
2501:                                         0,
2502:                                         0,
2503:                                         0,
2504:                                 /* 10*/ 0,
2505:                                         MatLUFactor_SeqDense,
2506:                                         MatCholeskyFactor_SeqDense,
2507:                                         MatSOR_SeqDense,
2508:                                         MatTranspose_SeqDense,
2509:                                 /* 15*/ MatGetInfo_SeqDense,
2510:                                         MatEqual_SeqDense,
2511:                                         MatGetDiagonal_SeqDense,
2512:                                         MatDiagonalScale_SeqDense,
2513:                                         MatNorm_SeqDense,
2514:                                 /* 20*/ MatAssemblyBegin_SeqDense,
2515:                                         MatAssemblyEnd_SeqDense,
2516:                                         MatSetOption_SeqDense,
2517:                                         MatZeroEntries_SeqDense,
2518:                                 /* 24*/ MatZeroRows_SeqDense,
2519:                                         0,
2520:                                         0,
2521:                                         0,
2522:                                         0,
2523:                                 /* 29*/ MatSetUp_SeqDense,
2524:                                         0,
2525:                                         0,
2526:                                         0,
2527:                                         0,
2528:                                 /* 34*/ MatDuplicate_SeqDense,
2529:                                         0,
2530:                                         0,
2531:                                         0,
2532:                                         0,
2533:                                 /* 39*/ MatAXPY_SeqDense,
2534:                                         MatCreateSubMatrices_SeqDense,
2535:                                         0,
2536:                                         MatGetValues_SeqDense,
2537:                                         MatCopy_SeqDense,
2538:                                 /* 44*/ MatGetRowMax_SeqDense,
2539:                                         MatScale_SeqDense,
2540:                                         MatShift_Basic,
2541:                                         0,
2542:                                         MatZeroRowsColumns_SeqDense,
2543:                                 /* 49*/ MatSetRandom_SeqDense,
2544:                                         0,
2545:                                         0,
2546:                                         0,
2547:                                         0,
2548:                                 /* 54*/ 0,
2549:                                         0,
2550:                                         0,
2551:                                         0,
2552:                                         0,
2553:                                 /* 59*/ 0,
2554:                                         MatDestroy_SeqDense,
2555:                                         MatView_SeqDense,
2556:                                         0,
2557:                                         0,
2558:                                 /* 64*/ 0,
2559:                                         0,
2560:                                         0,
2561:                                         0,
2562:                                         0,
2563:                                 /* 69*/ MatGetRowMaxAbs_SeqDense,
2564:                                         0,
2565:                                         0,
2566:                                         0,
2567:                                         0,
2568:                                 /* 74*/ 0,
2569:                                         0,
2570:                                         0,
2571:                                         0,
2572:                                         0,
2573:                                 /* 79*/ 0,
2574:                                         0,
2575:                                         0,
2576:                                         0,
2577:                                 /* 83*/ MatLoad_SeqDense,
2578:                                         0,
2579:                                         MatIsHermitian_SeqDense,
2580:                                         0,
2581:                                         0,
2582:                                         0,
2583:                                 /* 89*/ MatMatMult_SeqDense_SeqDense,
2584:                                         MatMatMultSymbolic_SeqDense_SeqDense,
2585:                                         MatMatMultNumeric_SeqDense_SeqDense,
2586:                                         MatPtAP_SeqDense_SeqDense,
2587:                                         MatPtAPSymbolic_SeqDense_SeqDense,
2588:                                 /* 94*/ MatPtAPNumeric_SeqDense_SeqDense,
2589:                                         MatMatTransposeMult_SeqDense_SeqDense,
2590:                                         MatMatTransposeMultSymbolic_SeqDense_SeqDense,
2591:                                         MatMatTransposeMultNumeric_SeqDense_SeqDense,
2592:                                         0,
2593:                                 /* 99*/ 0,
2594:                                         0,
2595:                                         0,
2596:                                         MatConjugate_SeqDense,
2597:                                         0,
2598:                                 /*104*/ 0,
2599:                                         MatRealPart_SeqDense,
2600:                                         MatImaginaryPart_SeqDense,
2601:                                         0,
2602:                                         0,
2603:                                 /*109*/ 0,
2604:                                         0,
2605:                                         MatGetRowMin_SeqDense,
2606:                                         MatGetColumnVector_SeqDense,
2607:                                         MatMissingDiagonal_SeqDense,
2608:                                 /*114*/ 0,
2609:                                         0,
2610:                                         0,
2611:                                         0,
2612:                                         0,
2613:                                 /*119*/ 0,
2614:                                         0,
2615:                                         0,
2616:                                         0,
2617:                                         0,
2618:                                 /*124*/ 0,
2619:                                         MatGetColumnNorms_SeqDense,
2620:                                         0,
2621:                                         0,
2622:                                         0,
2623:                                 /*129*/ 0,
2624:                                         MatTransposeMatMult_SeqDense_SeqDense,
2625:                                         MatTransposeMatMultSymbolic_SeqDense_SeqDense,
2626:                                         MatTransposeMatMultNumeric_SeqDense_SeqDense,
2627:                                         0,
2628:                                 /*134*/ 0,
2629:                                         0,
2630:                                         0,
2631:                                         0,
2632:                                         0,
2633:                                 /*139*/ 0,
2634:                                         0,
2635:                                         0,
2636:                                         0,
2637:                                         0,
2638:                                 /*144*/ MatCreateMPIMatConcatenateSeqMat_SeqDense
2639: };

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

2646:    Collective

2648:    Input Parameters:
2649: +  comm - MPI communicator, set to PETSC_COMM_SELF
2650: .  m - number of rows
2651: .  n - number of columns
2652: -  data - optional location of matrix data in column major order.  Set data=NULL for PETSc
2653:    to control all matrix memory allocation.

2655:    Output Parameter:
2656: .  A - the matrix

2658:    Notes:
2659:    The data input variable is intended primarily for Fortran programmers
2660:    who wish to allocate their own matrix memory space.  Most users should
2661:    set data=NULL.

2663:    Level: intermediate

2665: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
2666: @*/
2667: PetscErrorCode  MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
2668: {

2672:   MatCreate(comm,A);
2673:   MatSetSizes(*A,m,n,m,n);
2674:   MatSetType(*A,MATSEQDENSE);
2675:   MatSeqDenseSetPreallocation(*A,data);
2676:   return(0);
2677: }

2679: /*@C
2680:    MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements

2682:    Collective

2684:    Input Parameters:
2685: +  B - the matrix
2686: -  data - the array (or NULL)

2688:    Notes:
2689:    The data input variable is intended primarily for Fortran programmers
2690:    who wish to allocate their own matrix memory space.  Most users should
2691:    need not call this routine.

2693:    Level: intermediate

2695: .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA()

2697: @*/
2698: PetscErrorCode  MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
2699: {

2703:   PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
2704:   return(0);
2705: }

2707: PetscErrorCode  MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
2708: {
2709:   Mat_SeqDense   *b;

2713:   B->preallocated = PETSC_TRUE;

2715:   PetscLayoutSetUp(B->rmap);
2716:   PetscLayoutSetUp(B->cmap);

2718:   b       = (Mat_SeqDense*)B->data;
2719:   b->Mmax = B->rmap->n;
2720:   b->Nmax = B->cmap->n;
2721:   if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n;

2723:   PetscIntMultError(b->lda,b->Nmax,NULL);
2724:   if (!data) { /* petsc-allocated storage */
2725:     if (!b->user_alloc) { PetscFree(b->v); }
2726:     PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);
2727:     PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));

2729:     b->user_alloc = PETSC_FALSE;
2730:   } else { /* user-allocated storage */
2731:     if (!b->user_alloc) { PetscFree(b->v); }
2732:     b->v          = data;
2733:     b->user_alloc = PETSC_TRUE;
2734:   }
2735:   B->assembled = PETSC_TRUE;
2736:   return(0);
2737: }

2739: #if defined(PETSC_HAVE_ELEMENTAL)
2740: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
2741: {
2742:   Mat               mat_elemental;
2743:   PetscErrorCode    ierr;
2744:   const PetscScalar *array;
2745:   PetscScalar       *v_colwise;
2746:   PetscInt          M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols;

2749:   PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);
2750:   MatDenseGetArrayRead(A,&array);
2751:   /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
2752:   k = 0;
2753:   for (j=0; j<N; j++) {
2754:     cols[j] = j;
2755:     for (i=0; i<M; i++) {
2756:       v_colwise[j*M+i] = array[k++];
2757:     }
2758:   }
2759:   for (i=0; i<M; i++) {
2760:     rows[i] = i;
2761:   }
2762:   MatDenseRestoreArrayRead(A,&array);

2764:   MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
2765:   MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);
2766:   MatSetType(mat_elemental,MATELEMENTAL);
2767:   MatSetUp(mat_elemental);

2769:   /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
2770:   MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);
2771:   MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
2772:   MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);
2773:   PetscFree3(v_colwise,rows,cols);

2775:   if (reuse == MAT_INPLACE_MATRIX) {
2776:     MatHeaderReplace(A,&mat_elemental);
2777:   } else {
2778:     *newmat = mat_elemental;
2779:   }
2780:   return(0);
2781: }
2782: #endif

2784: /*@C
2785:   MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array

2787:   Input parameter:
2788: + A - the matrix
2789: - lda - the leading dimension

2791:   Notes:
2792:   This routine is to be used in conjunction with MatSeqDenseSetPreallocation();
2793:   it asserts that the preallocation has a leading dimension (the LDA parameter
2794:   of Blas and Lapack fame) larger than M, the first dimension of the matrix.

2796:   Level: intermediate

2798: .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize()

2800: @*/
2801: PetscErrorCode  MatSeqDenseSetLDA(Mat B,PetscInt lda)
2802: {
2803:   Mat_SeqDense *b = (Mat_SeqDense*)B->data;

2806:   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);
2807:   b->lda       = lda;
2808:   b->changelda = PETSC_FALSE;
2809:   b->Mmax      = PetscMax(b->Mmax,lda);
2810:   return(0);
2811: }

2813: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat)
2814: {
2816:   PetscMPIInt    size;

2819:   MPI_Comm_size(comm,&size);
2820:   if (size == 1) {
2821:     if (scall == MAT_INITIAL_MATRIX) {
2822:       MatDuplicate(inmat,MAT_COPY_VALUES,outmat);
2823:     } else {
2824:       MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);
2825:     }
2826:   } else {
2827:     MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);
2828:   }
2829:   return(0);
2830: }

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

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

2838:   Level: beginner

2840: .seealso: MatCreateSeqDense()

2842: M*/
2843: PetscErrorCode MatCreate_SeqDense(Mat B)
2844: {
2845:   Mat_SeqDense   *b;
2847:   PetscMPIInt    size;

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

2853:   PetscNewLog(B,&b);
2854:   PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
2855:   B->data = (void*)b;

2857:   b->roworiented = PETSC_TRUE;

2859:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetLDA_C",MatDenseGetLDA_SeqDense);
2860:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);
2861:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);
2862:   PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);
2863:   PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);
2864:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);
2865:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense);
2866:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);
2867: #if defined(PETSC_HAVE_ELEMENTAL)
2868:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);
2869: #endif
2870: #if defined(PETSC_HAVE_CUDA)
2871:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqdensecuda_C",MatConvert_SeqDense_SeqDenseCUDA);
2872:   PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijcusparse_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2873: #endif
2874:   PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);
2875:   PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2876: #if defined(PETSC_HAVE_VIENNACL)
2877:   PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijviennacl_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2878: #endif
2879:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2880:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2881:   PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqbaij_seqdense_C",MatMatMult_SeqBAIJ_SeqDense);
2882:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqbaij_seqdense_C",MatMatMultSymbolic_SeqBAIJ_SeqDense);
2883:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqbaij_seqdense_C",MatMatMultNumeric_SeqBAIJ_SeqDense);
2884:   PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqsbaij_seqdense_C",MatMatMult_SeqSBAIJ_SeqDense);
2885:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqsbaij_seqdense_C",MatMatMultSymbolic_SeqSBAIJ_SeqDense);
2886:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqsbaij_seqdense_C",MatMatMultNumeric_SeqSBAIJ_SeqDense);
2887:   PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaij_seqdense_C",MatPtAP_SeqDense_SeqDense);
2888:   PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijperm_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2889:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijperm_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2890:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijperm_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2891:   PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijperm_seqdense_C",MatPtAP_SeqDense_SeqDense);
2892:   PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijsell_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2893:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijsell_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2894:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijsell_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2895:   PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijsell_seqdense_C",MatPtAP_SeqDense_SeqDense);
2896:   PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijmkl_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2897:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijmkl_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2898:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijmkl_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2899:   PetscObjectComposeFunction((PetscObject)B,"MatMatMult_nest_seqdense_C",MatMatMult_Nest_Dense);
2900:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_nest_seqdense_C",MatMatMultSymbolic_Nest_Dense);
2901:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_nest_seqdense_C",MatMatMultNumeric_Nest_Dense);
2902:   PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijmkl_seqdense_C",MatPtAP_SeqDense_SeqDense);

2904:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2905:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2906:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2907:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijperm_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2908:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijperm_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2909:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijperm_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2910:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijsell_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2911:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijsell_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2912:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijsell_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);

2914:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijmkl_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2915:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijmkl_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2916:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijmkl_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2917:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);
2918:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);
2919:   PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
2920:   return(0);
2921: }

2923: /*@C
2924:    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.

2926:    Not Collective

2928:    Input Parameter:
2929: +  mat - a MATSEQDENSE or MATMPIDENSE matrix
2930: -  col - column index

2932:    Output Parameter:
2933: .  vals - pointer to the data

2935:    Level: intermediate

2937: .seealso: MatDenseRestoreColumn()
2938: @*/
2939: PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals)
2940: {

2944:   PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));
2945:   return(0);
2946: }

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

2951:    Not Collective

2953:    Input Parameter:
2954: .  mat - a MATSEQDENSE or MATMPIDENSE matrix

2956:    Output Parameter:
2957: .  vals - pointer to the data

2959:    Level: intermediate

2961: .seealso: MatDenseGetColumn()
2962: @*/
2963: PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals)
2964: {

2968:   PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));
2969:   return(0);
2970: }