Actual source code: dense.c

petsc-master 2017-08-22
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: static PetscErrorCode MatSeqDenseSymmetrize_Private(Mat A, PetscBool hermitian)
 12: {
 13:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
 14:   PetscInt      j, k, n = A->rmap->n;

 17:   if (A->rmap->n != A->cmap->n) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot symmetrize a rectangular matrix");
 18:   if (!hermitian) {
 19:     for (k=0;k<n;k++) {
 20:       for (j=k;j<n;j++) {
 21:         mat->v[j*mat->lda + k] = mat->v[k*mat->lda + j];
 22:       }
 23:     }
 24:   } else {
 25:     for (k=0;k<n;k++) {
 26:       for (j=k;j<n;j++) {
 27:         mat->v[j*mat->lda + k] = PetscConj(mat->v[k*mat->lda + j]);
 28:       }
 29:     }
 30:   }
 31:   return(0);
 32: }

 34: PETSC_EXTERN PetscErrorCode MatSeqDenseInvertFactors_Private(Mat A)
 35: {
 36: #if defined(PETSC_MISSING_LAPACK_POTRF)
 38:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable.");
 39: #else
 40:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
 42:   PetscBLASInt   info,n;

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

 78:   A->ops->solve             = NULL;
 79:   A->ops->matsolve          = NULL;
 80:   A->ops->solvetranspose    = NULL;
 81:   A->ops->matsolvetranspose = NULL;
 82:   A->ops->solveadd          = NULL;
 83:   A->ops->solvetransposeadd = NULL;
 84:   A->factortype             = MAT_FACTOR_NONE;
 85:   PetscFree(A->solvertype);
 86:   return(0);
 87: }

 89: PetscErrorCode MatZeroRowsColumns_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
 90: {
 91:   PetscErrorCode    ierr;
 92:   Mat_SeqDense      *l = (Mat_SeqDense*)A->data;
 93:   PetscInt          m  = l->lda, n = A->cmap->n,r = A->rmap->n, i,j;
 94:   PetscScalar       *slot,*bb;
 95:   const PetscScalar *xx;

 98: #if defined(PETSC_USE_DEBUG)
 99:   for (i=0; i<N; i++) {
100:     if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
101:     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);
102:     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);
103:   }
104: #endif

106:   /* fix right hand side if needed */
107:   if (x && b) {
108:     VecGetArrayRead(x,&xx);
109:     VecGetArray(b,&bb);
110:     for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
111:     VecRestoreArrayRead(x,&xx);
112:     VecRestoreArray(b,&bb);
113:   }

115:   for (i=0; i<N; i++) {
116:     slot = l->v + rows[i]*m;
117:     PetscMemzero(slot,r*sizeof(PetscScalar));
118:   }
119:   for (i=0; i<N; i++) {
120:     slot = l->v + rows[i];
121:     for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
122:   }
123:   if (diag != 0.0) {
124:     if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
125:     for (i=0; i<N; i++) {
126:       slot  = l->v + (m+1)*rows[i];
127:       *slot = diag;
128:     }
129:   }
130:   return(0);
131: }

133: PetscErrorCode MatPtAPNumeric_SeqDense_SeqDense(Mat A,Mat P,Mat C)
134: {
135:   Mat_SeqDense   *c = (Mat_SeqDense*)(C->data);

139:   MatMatMultNumeric_SeqDense_SeqDense(A,P,c->ptapwork);
140:   MatTransposeMatMultNumeric_SeqDense_SeqDense(P,c->ptapwork,C);
141:   return(0);
142: }

144: PetscErrorCode MatPtAPSymbolic_SeqDense_SeqDense(Mat A,Mat P,PetscReal fill,Mat *C)
145: {
146:   Mat_SeqDense   *c;

150:   MatCreateSeqDense(PetscObjectComm((PetscObject)A),P->cmap->N,P->cmap->N,NULL,C);
151:   c = (Mat_SeqDense*)((*C)->data);
152:   MatCreateSeqDense(PetscObjectComm((PetscObject)A),A->rmap->N,P->cmap->N,NULL,&c->ptapwork);
153:   return(0);
154: }

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

161:   if (reuse == MAT_INITIAL_MATRIX) {
162:     MatPtAPSymbolic_SeqDense_SeqDense(A,P,fill,C);
163:   }
164:   PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);
165:   (*(*C)->ops->ptapnumeric)(A,P,*C);
166:   PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);
167:   return(0);
168: }

170: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
171: {
172:   Mat            B = NULL;
173:   Mat_SeqAIJ     *a = (Mat_SeqAIJ*)A->data;
174:   Mat_SeqDense   *b;
176:   PetscInt       *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->N,i;
177:   MatScalar      *av=a->a;
178:   PetscBool      isseqdense;

181:   if (reuse == MAT_REUSE_MATRIX) {
182:     PetscObjectTypeCompare((PetscObject)*newmat,MATSEQDENSE,&isseqdense);
183:     if (!isseqdense) SETERRQ1(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_USER,"Cannot reuse matrix of type %s",((PetscObject)(*newmat))->type);
184:   }
185:   if (reuse != MAT_REUSE_MATRIX) {
186:     MatCreate(PetscObjectComm((PetscObject)A),&B);
187:     MatSetSizes(B,m,n,m,n);
188:     MatSetType(B,MATSEQDENSE);
189:     MatSeqDenseSetPreallocation(B,NULL);
190:     b    = (Mat_SeqDense*)(B->data);
191:   } else {
192:     b    = (Mat_SeqDense*)((*newmat)->data);
193:     PetscMemzero(b->v,m*n*sizeof(PetscScalar));
194:   }
195:   for (i=0; i<m; i++) {
196:     PetscInt j;
197:     for (j=0;j<ai[1]-ai[0];j++) {
198:       b->v[*aj*m+i] = *av;
199:       aj++;
200:       av++;
201:     }
202:     ai++;
203:   }

205:   if (reuse == MAT_INPLACE_MATRIX) {
206:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
207:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
208:     MatHeaderReplace(A,&B);
209:   } else {
210:     if (B) *newmat = B;
211:     MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
212:     MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
213:   }
214:   return(0);
215: }

217: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
218: {
219:   Mat            B;
220:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
222:   PetscInt       i, j;
223:   PetscInt       *rows, *nnz;
224:   MatScalar      *aa = a->v, *vals;

227:   MatCreate(PetscObjectComm((PetscObject)A),&B);
228:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
229:   MatSetType(B,MATSEQAIJ);
230:   PetscCalloc3(A->rmap->n,&rows,A->rmap->n,&nnz,A->rmap->n,&vals);
231:   for (j=0; j<A->cmap->n; j++) {
232:     for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) ++nnz[i];
233:     aa += a->lda;
234:   }
235:   MatSeqAIJSetPreallocation(B,PETSC_DETERMINE,nnz);
236:   aa = a->v;
237:   for (j=0; j<A->cmap->n; j++) {
238:     PetscInt numRows = 0;
239:     for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) {rows[numRows] = i; vals[numRows++] = aa[i];}
240:     MatSetValues(B,numRows,rows,1,&j,vals,INSERT_VALUES);
241:     aa  += a->lda;
242:   }
243:   PetscFree3(rows,nnz,vals);
244:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
245:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

247:   if (reuse == MAT_INPLACE_MATRIX) {
248:     MatHeaderReplace(A,&B);
249:   } else {
250:     *newmat = B;
251:   }
252:   return(0);
253: }

255: static PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str)
256: {
257:   Mat_SeqDense   *x     = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data;
258:   PetscScalar    oalpha = alpha;
259:   PetscInt       j;
260:   PetscBLASInt   N,m,ldax,lday,one = 1;

264:   PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);
265:   PetscBLASIntCast(X->rmap->n,&m);
266:   PetscBLASIntCast(x->lda,&ldax);
267:   PetscBLASIntCast(y->lda,&lday);
268:   if (ldax>m || lday>m) {
269:     for (j=0; j<X->cmap->n; j++) {
270:       PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&oalpha,x->v+j*ldax,&one,y->v+j*lday,&one));
271:     }
272:   } else {
273:     PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&oalpha,x->v,&one,y->v,&one));
274:   }
275:   PetscObjectStateIncrease((PetscObject)Y);
276:   PetscLogFlops(PetscMax(2*N-1,0));
277:   return(0);
278: }

280: static PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info)
281: {
282:   PetscInt N = A->rmap->n*A->cmap->n;

285:   info->block_size        = 1.0;
286:   info->nz_allocated      = (double)N;
287:   info->nz_used           = (double)N;
288:   info->nz_unneeded       = (double)0;
289:   info->assemblies        = (double)A->num_ass;
290:   info->mallocs           = 0;
291:   info->memory            = ((PetscObject)A)->mem;
292:   info->fill_ratio_given  = 0;
293:   info->fill_ratio_needed = 0;
294:   info->factor_mallocs    = 0;
295:   return(0);
296: }

298: static PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha)
299: {
300:   Mat_SeqDense   *a     = (Mat_SeqDense*)A->data;
301:   PetscScalar    oalpha = alpha;
303:   PetscBLASInt   one = 1,j,nz,lda;

306:   PetscBLASIntCast(a->lda,&lda);
307:   if (lda>A->rmap->n) {
308:     PetscBLASIntCast(A->rmap->n,&nz);
309:     for (j=0; j<A->cmap->n; j++) {
310:       PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v+j*lda,&one));
311:     }
312:   } else {
313:     PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);
314:     PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v,&one));
315:   }
316:   PetscLogFlops(nz);
317:   return(0);
318: }

320: static PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool  *fl)
321: {
322:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
323:   PetscInt     i,j,m = A->rmap->n,N;
324:   PetscScalar  *v = a->v;

327:   *fl = PETSC_FALSE;
328:   if (A->rmap->n != A->cmap->n) return(0);
329:   N = a->lda;

331:   for (i=0; i<m; i++) {
332:     for (j=i+1; j<m; j++) {
333:       if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) return(0);
334:     }
335:   }
336:   *fl = PETSC_TRUE;
337:   return(0);
338: }

340: static PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues)
341: {
342:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data,*l;
344:   PetscInt       lda = (PetscInt)mat->lda,j,m;

347:   PetscLayoutReference(A->rmap,&newi->rmap);
348:   PetscLayoutReference(A->cmap,&newi->cmap);
349:   MatSeqDenseSetPreallocation(newi,NULL);
350:   if (cpvalues == MAT_COPY_VALUES) {
351:     l = (Mat_SeqDense*)newi->data;
352:     if (lda>A->rmap->n) {
353:       m = A->rmap->n;
354:       for (j=0; j<A->cmap->n; j++) {
355:         PetscMemcpy(l->v+j*m,mat->v+j*lda,m*sizeof(PetscScalar));
356:       }
357:     } else {
358:       PetscMemcpy(l->v,mat->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
359:     }
360:   }
361:   newi->assembled = PETSC_TRUE;
362:   return(0);
363: }

365: static PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat)
366: {

370:   MatCreate(PetscObjectComm((PetscObject)A),newmat);
371:   MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
372:   MatSetType(*newmat,((PetscObject)A)->type_name);
373:   MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);
374:   return(0);
375: }


378: static PetscErrorCode MatLUFactor_SeqDense(Mat,IS,IS,const MatFactorInfo*);

380: static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
381: {
382:   MatFactorInfo  info;

386:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
387:   MatLUFactor_SeqDense(fact,0,0,&info);
388:   return(0);
389: }

391: static PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy)
392: {
393:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
394:   PetscErrorCode    ierr;
395:   const PetscScalar *x;
396:   PetscScalar       *y;
397:   PetscBLASInt      one = 1,info,m;

400:   PetscBLASIntCast(A->rmap->n,&m);
401:   VecGetArrayRead(xx,&x);
402:   VecGetArray(yy,&y);
403:   PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
404:   if (A->factortype == MAT_FACTOR_LU) {
405: #if defined(PETSC_MISSING_LAPACK_GETRS)
406:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
407: #else
408:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
409:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
410: #endif
411:   } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
412: #if defined(PETSC_MISSING_LAPACK_POTRS)
413:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
414: #else
415:     if (A->spd) {
416:       PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
417:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
418: #if defined (PETSC_USE_COMPLEX)
419:     } else if (A->hermitian) {
420:       PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
421:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve");
422: #endif
423:     } else { /* symmetric case */
424:       PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
425:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
426:     }
427: #endif
428:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
429:   VecRestoreArrayRead(xx,&x);
430:   VecRestoreArray(yy,&y);
431:   PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
432:   return(0);
433: }

435: static PetscErrorCode MatMatSolve_SeqDense(Mat A,Mat B,Mat X)
436: {
437:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
439:   PetscScalar    *b,*x;
440:   PetscInt       n;
441:   PetscBLASInt   nrhs,info,m;
442:   PetscBool      flg;

445:   PetscBLASIntCast(A->rmap->n,&m);
446:   PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
447:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
448:   PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
449:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");

451:   MatGetSize(B,NULL,&n);
452:   PetscBLASIntCast(n,&nrhs);
453:   MatDenseGetArray(B,&b);
454:   MatDenseGetArray(X,&x);

456:   PetscMemcpy(x,b,m*nrhs*sizeof(PetscScalar));

458:   if (A->factortype == MAT_FACTOR_LU) {
459: #if defined(PETSC_MISSING_LAPACK_GETRS)
460:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
461: #else
462:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
463:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
464: #endif
465:   } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
466: #if defined(PETSC_MISSING_LAPACK_POTRS)
467:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
468: #else
469:     if (A->spd) {
470:       PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info));
471:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
472: #if defined (PETSC_USE_COMPLEX)
473:     } else if (A->hermitian) {
474:       PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
475:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve");
476: #endif
477:     } else { /* symmetric case */
478:       PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
479:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
480:     }
481: #endif
482:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");

484:   MatDenseRestoreArray(B,&b);
485:   MatDenseRestoreArray(X,&x);
486:   PetscLogFlops(nrhs*(2.0*m*m - m));
487:   return(0);
488: }

490: static PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy)
491: {
492:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
493:   PetscErrorCode    ierr;
494:   const PetscScalar *x;
495:   PetscScalar       *y;
496:   PetscBLASInt      one = 1,info,m;

499:   PetscBLASIntCast(A->rmap->n,&m);
500:   VecGetArrayRead(xx,&x);
501:   VecGetArray(yy,&y);
502:   PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
503:   if (A->factortype == MAT_FACTOR_LU) {
504: #if defined(PETSC_MISSING_LAPACK_GETRS)
505:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
506: #else
507:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
508:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve");
509: #endif
510:   } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
511: #if defined(PETSC_MISSING_LAPACK_POTRS)
512:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
513: #else
514:     if (A->spd) {
515:       PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
516:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
517: #if defined (PETSC_USE_COMPLEX)
518:     } else if (A->hermitian) {
519:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatSolveTranspose with Cholesky/Hermitian not available");
520: #endif
521:     } else { /* symmetric case */
522:       PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
523:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
524:     }
525: #endif
526:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
527:   VecRestoreArrayRead(xx,&x);
528:   VecRestoreArray(yy,&y);
529:   PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
530:   return(0);
531: }

533: static PetscErrorCode MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
534: {
535:   PetscErrorCode    ierr;
536:   const PetscScalar *x;
537:   PetscScalar       *y,sone = 1.0;
538:   Vec               tmp = 0;

541:   VecGetArrayRead(xx,&x);
542:   VecGetArray(yy,&y);
543:   if (!A->rmap->n || !A->cmap->n) return(0);
544:   if (yy == zz) {
545:     VecDuplicate(yy,&tmp);
546:     PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);
547:     VecCopy(yy,tmp);
548:   }
549:   MatSolve_SeqDense(A,xx,yy);
550:   if (tmp) {
551:     VecAXPY(yy,sone,tmp);
552:     VecDestroy(&tmp);
553:   } else {
554:     VecAXPY(yy,sone,zz);
555:   }
556:   VecRestoreArrayRead(xx,&x);
557:   VecRestoreArray(yy,&y);
558:   PetscLogFlops(A->cmap->n);
559:   return(0);
560: }

562: static PetscErrorCode MatSolveTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
563: {
564:   PetscErrorCode    ierr;
565:   const PetscScalar *x;
566:   PetscScalar       *y,sone = 1.0;
567:   Vec               tmp = NULL;

570:   if (!A->rmap->n || !A->cmap->n) return(0);
571:   VecGetArrayRead(xx,&x);
572:   VecGetArray(yy,&y);
573:   if (yy == zz) {
574:     VecDuplicate(yy,&tmp);
575:     PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);
576:     VecCopy(yy,tmp);
577:   }
578:   MatSolveTranspose_SeqDense(A,xx,yy);
579:   if (tmp) {
580:     VecAXPY(yy,sone,tmp);
581:     VecDestroy(&tmp);
582:   } else {
583:     VecAXPY(yy,sone,zz);
584:   }
585:   VecRestoreArrayRead(xx,&x);
586:   VecRestoreArray(yy,&y);
587:   PetscLogFlops(A->cmap->n);
588:   return(0);
589: }

591: /* ---------------------------------------------------------------*/
592: /* COMMENT: I have chosen to hide row permutation in the pivots,
593:    rather than put it in the Mat->row slot.*/
594: static PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo)
595: {
596: #if defined(PETSC_MISSING_LAPACK_GETRF)
598:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable.");
599: #else
600:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
602:   PetscBLASInt   n,m,info;

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

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

619:   A->ops->solve             = MatSolve_SeqDense;
620:   A->ops->matsolve          = MatMatSolve_SeqDense;
621:   A->ops->solvetranspose    = MatSolveTranspose_SeqDense;
622:   A->ops->solveadd          = MatSolveAdd_SeqDense;
623:   A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
624:   A->factortype             = MAT_FACTOR_LU;

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

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

634: /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */
635: static PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo)
636: {
637: #if defined(PETSC_MISSING_LAPACK_POTRF)
639:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable.");
640: #else
641:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
643:   PetscBLASInt   info,n;

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

659:       mat->lfwork = -1;
660:       PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
661:       mat->lfwork = (PetscInt)PetscRealPart(dummy);
662:       PetscMalloc1(mat->lfwork,&mat->fwork);
663:       PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
664:     }
665:     PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
666: #endif
667:   } else { /* symmetric case */
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:       PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
677:       mat->lfwork = (PetscInt)PetscRealPart(dummy);
678:       PetscMalloc1(mat->lfwork,&mat->fwork);
679:       PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
680:     }
681:     PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
682:   }
683:   if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1);

685:   A->ops->solve             = MatSolve_SeqDense;
686:   A->ops->matsolve          = MatMatSolve_SeqDense;
687:   A->ops->solvetranspose    = MatSolveTranspose_SeqDense;
688:   A->ops->solveadd          = MatSolveAdd_SeqDense;
689:   A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
690:   A->factortype             = MAT_FACTOR_CHOLESKY;

692:   PetscFree(A->solvertype);
693:   PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);

695:   PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
696: #endif
697:   return(0);
698: }


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

707:   info.fill = 1.0;

709:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
710:   MatCholeskyFactor_SeqDense(fact,0,&info);
711:   return(0);
712: }

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

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

742: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact)
743: {

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

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

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

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

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

813:   PetscBLASIntCast(A->rmap->n,&m);
814:   PetscBLASIntCast(A->cmap->n,&n);
815:   if (!A->rmap->n || !A->cmap->n) return(0);
816:   VecGetArrayRead(xx,&x);
817:   VecGetArray(yy,&y);
818:   PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One));
819:   VecRestoreArrayRead(xx,&x);
820:   VecRestoreArray(yy,&y);
821:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);
822:   return(0);
823: }

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

834:   PetscBLASIntCast(A->rmap->n,&m);
835:   PetscBLASIntCast(A->cmap->n,&n);
836:   if (!A->rmap->n || !A->cmap->n) return(0);
837:   VecGetArrayRead(xx,&x);
838:   VecGetArray(yy,&y);
839:   PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One));
840:   VecRestoreArrayRead(xx,&x);
841:   VecRestoreArray(yy,&y);
842:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);
843:   return(0);
844: }

846: PETSC_INTERN PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
847: {
848:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
849:   const PetscScalar *v = mat->v,*x;
850:   PetscScalar       *y,_DOne=1.0;
851:   PetscErrorCode    ierr;
852:   PetscBLASInt      m, n, _One=1;

855:   PetscBLASIntCast(A->rmap->n,&m);
856:   PetscBLASIntCast(A->cmap->n,&n);
857:   if (!A->rmap->n || !A->cmap->n) return(0);
858:   if (zz != yy) {VecCopy(zz,yy);}
859:   VecGetArrayRead(xx,&x);
860:   VecGetArray(yy,&y);
861:   PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
862:   VecRestoreArrayRead(xx,&x);
863:   VecRestoreArray(yy,&y);
864:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
865:   return(0);
866: }

868: static PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
869: {
870:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
871:   const PetscScalar *v = mat->v,*x;
872:   PetscScalar       *y;
873:   PetscErrorCode    ierr;
874:   PetscBLASInt      m, n, _One=1;
875:   PetscScalar       _DOne=1.0;

878:   PetscBLASIntCast(A->rmap->n,&m);
879:   PetscBLASIntCast(A->cmap->n,&n);
880:   if (!A->rmap->n || !A->cmap->n) return(0);
881:   if (zz != yy) {VecCopy(zz,yy);}
882:   VecGetArrayRead(xx,&x);
883:   VecGetArray(yy,&y);
884:   PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
885:   VecRestoreArrayRead(xx,&x);
886:   VecRestoreArray(yy,&y);
887:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
888:   return(0);
889: }

891: /* -----------------------------------------------------------------*/
892: static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
893: {
894:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
895:   PetscScalar    *v;
897:   PetscInt       i;

900:   *ncols = A->cmap->n;
901:   if (cols) {
902:     PetscMalloc1(A->cmap->n+1,cols);
903:     for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
904:   }
905:   if (vals) {
906:     PetscMalloc1(A->cmap->n+1,vals);
907:     v    = mat->v + row;
908:     for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
909:   }
910:   return(0);
911: }

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

918:   if (cols) {PetscFree(*cols);}
919:   if (vals) {PetscFree(*vals); }
920:   return(0);
921: }
922: /* ----------------------------------------------------------------*/
923: static PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
924: {
925:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
926:   PetscInt     i,j,idx=0;

929:   if (!mat->roworiented) {
930:     if (addv == INSERT_VALUES) {
931:       for (j=0; j<n; j++) {
932:         if (indexn[j] < 0) {idx += m; continue;}
933: #if defined(PETSC_USE_DEBUG)
934:         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);
935: #endif
936:         for (i=0; i<m; i++) {
937:           if (indexm[i] < 0) {idx++; continue;}
938: #if defined(PETSC_USE_DEBUG)
939:           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);
940: #endif
941:           mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
942:         }
943:       }
944:     } else {
945:       for (j=0; j<n; j++) {
946:         if (indexn[j] < 0) {idx += m; continue;}
947: #if defined(PETSC_USE_DEBUG)
948:         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);
949: #endif
950:         for (i=0; i<m; i++) {
951:           if (indexm[i] < 0) {idx++; continue;}
952: #if defined(PETSC_USE_DEBUG)
953:           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);
954: #endif
955:           mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
956:         }
957:       }
958:     }
959:   } else {
960:     if (addv == INSERT_VALUES) {
961:       for (i=0; i<m; i++) {
962:         if (indexm[i] < 0) { idx += n; continue;}
963: #if defined(PETSC_USE_DEBUG)
964:         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);
965: #endif
966:         for (j=0; j<n; j++) {
967:           if (indexn[j] < 0) { idx++; continue;}
968: #if defined(PETSC_USE_DEBUG)
969:           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);
970: #endif
971:           mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
972:         }
973:       }
974:     } else {
975:       for (i=0; i<m; i++) {
976:         if (indexm[i] < 0) { idx += n; 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:         for (j=0; j<n; j++) {
981:           if (indexn[j] < 0) { idx++; continue;}
982: #if defined(PETSC_USE_DEBUG)
983:           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);
984: #endif
985:           mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
986:         }
987:       }
988:     }
989:   }
990:   return(0);
991: }

993: static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
994: {
995:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
996:   PetscInt     i,j;

999:   /* row-oriented output */
1000:   for (i=0; i<m; i++) {
1001:     if (indexm[i] < 0) {v += n;continue;}
1002:     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);
1003:     for (j=0; j<n; j++) {
1004:       if (indexn[j] < 0) {v++; continue;}
1005:       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);
1006:       *v++ = mat->v[indexn[j]*mat->lda + indexm[i]];
1007:     }
1008:   }
1009:   return(0);
1010: }

1012: /* -----------------------------------------------------------------*/

1014: static PetscErrorCode MatLoad_SeqDense(Mat newmat,PetscViewer viewer)
1015: {
1016:   Mat_SeqDense   *a;
1018:   PetscInt       *scols,i,j,nz,header[4];
1019:   int            fd;
1020:   PetscMPIInt    size;
1021:   PetscInt       *rowlengths = 0,M,N,*cols,grows,gcols;
1022:   PetscScalar    *vals,*svals,*v,*w;
1023:   MPI_Comm       comm;

1026:   /* force binary viewer to load .info file if it has not yet done so */
1027:   PetscViewerSetUp(viewer);
1028:   PetscObjectGetComm((PetscObject)viewer,&comm);
1029:   MPI_Comm_size(comm,&size);
1030:   if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor");
1031:   PetscViewerBinaryGetDescriptor(viewer,&fd);
1032:   PetscBinaryRead(fd,header,4,PETSC_INT);
1033:   if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not matrix object");
1034:   M = header[1]; N = header[2]; nz = header[3];

1036:   /* set global size if not set already*/
1037:   if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) {
1038:     MatSetSizes(newmat,M,N,M,N);
1039:   } else {
1040:     /* if sizes and type are already set, check if the vector global sizes are correct */
1041:     MatGetSize(newmat,&grows,&gcols);
1042:     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);
1043:   }
1044:   a = (Mat_SeqDense*)newmat->data;
1045:   if (!a->user_alloc) {
1046:     MatSeqDenseSetPreallocation(newmat,NULL);
1047:   }

1049:   if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */
1050:     a = (Mat_SeqDense*)newmat->data;
1051:     v = a->v;
1052:     /* Allocate some temp space to read in the values and then flip them
1053:        from row major to column major */
1054:     PetscMalloc1(M*N > 0 ? M*N : 1,&w);
1055:     /* read in nonzero values */
1056:     PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);
1057:     /* now flip the values and store them in the matrix*/
1058:     for (j=0; j<N; j++) {
1059:       for (i=0; i<M; i++) {
1060:         *v++ =w[i*N+j];
1061:       }
1062:     }
1063:     PetscFree(w);
1064:   } else {
1065:     /* read row lengths */
1066:     PetscMalloc1(M+1,&rowlengths);
1067:     PetscBinaryRead(fd,rowlengths,M,PETSC_INT);

1069:     a = (Mat_SeqDense*)newmat->data;
1070:     v = a->v;

1072:     /* read column indices and nonzeros */
1073:     PetscMalloc1(nz+1,&scols);
1074:     cols = scols;
1075:     PetscBinaryRead(fd,cols,nz,PETSC_INT);
1076:     PetscMalloc1(nz+1,&svals);
1077:     vals = svals;
1078:     PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);

1080:     /* insert into matrix */
1081:     for (i=0; i<M; i++) {
1082:       for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j];
1083:       svals += rowlengths[i]; scols += rowlengths[i];
1084:     }
1085:     PetscFree(vals);
1086:     PetscFree(cols);
1087:     PetscFree(rowlengths);
1088:   }
1089:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1090:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
1091:   return(0);
1092: }

1094: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
1095: {
1096:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
1097:   PetscErrorCode    ierr;
1098:   PetscInt          i,j;
1099:   const char        *name;
1100:   PetscScalar       *v;
1101:   PetscViewerFormat format;
1102: #if defined(PETSC_USE_COMPLEX)
1103:   PetscBool         allreal = PETSC_TRUE;
1104: #endif

1107:   PetscViewerGetFormat(viewer,&format);
1108:   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1109:     return(0);  /* do nothing for now */
1110:   } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1111:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1112:     for (i=0; i<A->rmap->n; i++) {
1113:       v    = a->v + i;
1114:       PetscViewerASCIIPrintf(viewer,"row %D:",i);
1115:       for (j=0; j<A->cmap->n; j++) {
1116: #if defined(PETSC_USE_COMPLEX)
1117:         if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
1118:           PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1119:         } else if (PetscRealPart(*v)) {
1120:           PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));
1121:         }
1122: #else
1123:         if (*v) {
1124:           PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);
1125:         }
1126: #endif
1127:         v += a->lda;
1128:       }
1129:       PetscViewerASCIIPrintf(viewer,"\n");
1130:     }
1131:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1132:   } else {
1133:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1134: #if defined(PETSC_USE_COMPLEX)
1135:     /* determine if matrix has all real values */
1136:     v = a->v;
1137:     for (i=0; i<A->rmap->n*A->cmap->n; i++) {
1138:       if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;}
1139:     }
1140: #endif
1141:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
1142:       PetscObjectGetName((PetscObject)A,&name);
1143:       PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);
1144:       PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);
1145:       PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
1146:     }

1148:     for (i=0; i<A->rmap->n; i++) {
1149:       v = a->v + i;
1150:       for (j=0; j<A->cmap->n; j++) {
1151: #if defined(PETSC_USE_COMPLEX)
1152:         if (allreal) {
1153:           PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));
1154:         } else {
1155:           PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1156:         }
1157: #else
1158:         PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);
1159: #endif
1160:         v += a->lda;
1161:       }
1162:       PetscViewerASCIIPrintf(viewer,"\n");
1163:     }
1164:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
1165:       PetscViewerASCIIPrintf(viewer,"];\n");
1166:     }
1167:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1168:   }
1169:   PetscViewerFlush(viewer);
1170:   return(0);
1171: }

1173: static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer)
1174: {
1175:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
1176:   PetscErrorCode    ierr;
1177:   int               fd;
1178:   PetscInt          ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n;
1179:   PetscScalar       *v,*anonz,*vals;
1180:   PetscViewerFormat format;

1183:   PetscViewerBinaryGetDescriptor(viewer,&fd);

1185:   PetscViewerGetFormat(viewer,&format);
1186:   if (format == PETSC_VIEWER_NATIVE) {
1187:     /* store the matrix as a dense matrix */
1188:     PetscMalloc1(4,&col_lens);

1190:     col_lens[0] = MAT_FILE_CLASSID;
1191:     col_lens[1] = m;
1192:     col_lens[2] = n;
1193:     col_lens[3] = MATRIX_BINARY_FORMAT_DENSE;

1195:     PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);
1196:     PetscFree(col_lens);

1198:     /* write out matrix, by rows */
1199:     PetscMalloc1(m*n+1,&vals);
1200:     v    = a->v;
1201:     for (j=0; j<n; j++) {
1202:       for (i=0; i<m; i++) {
1203:         vals[j + i*n] = *v++;
1204:       }
1205:     }
1206:     PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);
1207:     PetscFree(vals);
1208:   } else {
1209:     PetscMalloc1(4+nz,&col_lens);

1211:     col_lens[0] = MAT_FILE_CLASSID;
1212:     col_lens[1] = m;
1213:     col_lens[2] = n;
1214:     col_lens[3] = nz;

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

1220:     /* Possibly should write in smaller increments, not whole matrix at once? */
1221:     /* store column indices (zero start index) */
1222:     ict = 0;
1223:     for (i=0; i<m; i++) {
1224:       for (j=0; j<n; j++) col_lens[ict++] = j;
1225:     }
1226:     PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);
1227:     PetscFree(col_lens);

1229:     /* store nonzero values */
1230:     PetscMalloc1(nz+1,&anonz);
1231:     ict  = 0;
1232:     for (i=0; i<m; i++) {
1233:       v = a->v + i;
1234:       for (j=0; j<n; j++) {
1235:         anonz[ict++] = *v; v += a->lda;
1236:       }
1237:     }
1238:     PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);
1239:     PetscFree(anonz);
1240:   }
1241:   return(0);
1242: }

1244:  #include <petscdraw.h>
1245: static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa)
1246: {
1247:   Mat               A  = (Mat) Aa;
1248:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
1249:   PetscErrorCode    ierr;
1250:   PetscInt          m  = A->rmap->n,n = A->cmap->n,i,j;
1251:   int               color = PETSC_DRAW_WHITE;
1252:   PetscScalar       *v = a->v;
1253:   PetscViewer       viewer;
1254:   PetscReal         xl,yl,xr,yr,x_l,x_r,y_l,y_r;
1255:   PetscViewerFormat format;

1258:   PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
1259:   PetscViewerGetFormat(viewer,&format);
1260:   PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);

1262:   /* Loop over matrix elements drawing boxes */

1264:   if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1265:     PetscDrawCollectiveBegin(draw);
1266:     /* Blue for negative and Red for positive */
1267:     for (j = 0; j < n; j++) {
1268:       x_l = j; x_r = x_l + 1.0;
1269:       for (i = 0; i < m; i++) {
1270:         y_l = m - i - 1.0;
1271:         y_r = y_l + 1.0;
1272:         if (PetscRealPart(v[j*m+i]) >  0.) {
1273:           color = PETSC_DRAW_RED;
1274:         } else if (PetscRealPart(v[j*m+i]) <  0.) {
1275:           color = PETSC_DRAW_BLUE;
1276:         } else {
1277:           continue;
1278:         }
1279:         PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1280:       }
1281:     }
1282:     PetscDrawCollectiveEnd(draw);
1283:   } else {
1284:     /* use contour shading to indicate magnitude of values */
1285:     /* first determine max of all nonzero values */
1286:     PetscReal minv = 0.0, maxv = 0.0;
1287:     PetscDraw popup;

1289:     for (i=0; i < m*n; i++) {
1290:       if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1291:     }
1292:     if (minv >= maxv) maxv = minv + PETSC_SMALL;
1293:     PetscDrawGetPopup(draw,&popup);
1294:     PetscDrawScalePopup(popup,minv,maxv);

1296:     PetscDrawCollectiveBegin(draw);
1297:     for (j=0; j<n; j++) {
1298:       x_l = j;
1299:       x_r = x_l + 1.0;
1300:       for (i=0; i<m; i++) {
1301:         y_l = m - i - 1.0;
1302:         y_r = y_l + 1.0;
1303:         color = PetscDrawRealToColor(PetscAbsScalar(v[j*m+i]),minv,maxv);
1304:         PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1305:       }
1306:     }
1307:     PetscDrawCollectiveEnd(draw);
1308:   }
1309:   return(0);
1310: }

1312: static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
1313: {
1314:   PetscDraw      draw;
1315:   PetscBool      isnull;
1316:   PetscReal      xr,yr,xl,yl,h,w;

1320:   PetscViewerDrawGetDraw(viewer,0,&draw);
1321:   PetscDrawIsNull(draw,&isnull);
1322:   if (isnull) return(0);

1324:   xr   = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
1325:   xr  += w;          yr += h;        xl = -w;     yl = -h;
1326:   PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
1327:   PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
1328:   PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
1329:   PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
1330:   PetscDrawSave(draw);
1331:   return(0);
1332: }

1334: PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
1335: {
1337:   PetscBool      iascii,isbinary,isdraw;

1340:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1341:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1342:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);

1344:   if (iascii) {
1345:     MatView_SeqDense_ASCII(A,viewer);
1346:   } else if (isbinary) {
1347:     MatView_SeqDense_Binary(A,viewer);
1348:   } else if (isdraw) {
1349:     MatView_SeqDense_Draw(A,viewer);
1350:   }
1351:   return(0);
1352: }

1354: static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A,const PetscScalar array[])
1355: {
1356:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

1359:   a->unplacedarray       = a->v;
1360:   a->unplaced_user_alloc = a->user_alloc;
1361:   a->v                   = (PetscScalar*) array;
1362:   return(0);
1363: }

1365: static PetscErrorCode MatDenseResetArray_SeqDense(Mat A)
1366: {
1367:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

1370:   a->v             = a->unplacedarray;
1371:   a->user_alloc    = a->unplaced_user_alloc;
1372:   a->unplacedarray = NULL;
1373:   return(0);
1374: }

1376: static PetscErrorCode MatDestroy_SeqDense(Mat mat)
1377: {
1378:   Mat_SeqDense   *l = (Mat_SeqDense*)mat->data;

1382: #if defined(PETSC_USE_LOG)
1383:   PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n);
1384: #endif
1385:   PetscFree(l->pivots);
1386:   PetscFree(l->fwork);
1387:   MatDestroy(&l->ptapwork);
1388:   if (!l->user_alloc) {PetscFree(l->v);}
1389:   PetscFree(mat->data);

1391:   PetscObjectChangeTypeName((PetscObject)mat,0);
1392:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);
1393:   PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL);
1394:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL);
1395:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);
1396:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);
1397: #if defined(PETSC_HAVE_ELEMENTAL)
1398:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);
1399: #endif
1400:   PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);
1401:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);
1402:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);
1403:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);
1404:   PetscObjectComposeFunction((PetscObject)mat,"MatPtAP_seqaij_seqdense_C",NULL);
1405:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);
1406:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);
1407:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);
1408:   return(0);
1409: }

1411: static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
1412: {
1413:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1415:   PetscInt       k,j,m,n,M;
1416:   PetscScalar    *v,tmp;

1419:   v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n;
1420:   if (reuse == MAT_INPLACE_MATRIX) { /* in place transpose */
1421:     if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can not transpose non-square matrix in place");
1422:     else {
1423:       for (j=0; j<m; j++) {
1424:         for (k=0; k<j; k++) {
1425:           tmp        = v[j + k*M];
1426:           v[j + k*M] = v[k + j*M];
1427:           v[k + j*M] = tmp;
1428:         }
1429:       }
1430:     }
1431:   } else { /* out-of-place transpose */
1432:     Mat          tmat;
1433:     Mat_SeqDense *tmatd;
1434:     PetscScalar  *v2;
1435:     PetscInt     M2;

1437:     if (reuse == MAT_INITIAL_MATRIX) {
1438:       MatCreate(PetscObjectComm((PetscObject)A),&tmat);
1439:       MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);
1440:       MatSetType(tmat,((PetscObject)A)->type_name);
1441:       MatSeqDenseSetPreallocation(tmat,NULL);
1442:     } else {
1443:       tmat = *matout;
1444:     }
1445:     tmatd = (Mat_SeqDense*)tmat->data;
1446:     v = mat->v; v2 = tmatd->v; M2 = tmatd->lda;
1447:     for (j=0; j<n; j++) {
1448:       for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M];
1449:     }
1450:     MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
1451:     MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
1452:     *matout = tmat;
1453:   }
1454:   return(0);
1455: }

1457: static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool  *flg)
1458: {
1459:   Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data;
1460:   Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data;
1461:   PetscInt     i,j;
1462:   PetscScalar  *v1,*v2;

1465:   if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; return(0);}
1466:   if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; return(0);}
1467:   for (i=0; i<A1->rmap->n; i++) {
1468:     v1 = mat1->v+i; v2 = mat2->v+i;
1469:     for (j=0; j<A1->cmap->n; j++) {
1470:       if (*v1 != *v2) {*flg = PETSC_FALSE; return(0);}
1471:       v1 += mat1->lda; v2 += mat2->lda;
1472:     }
1473:   }
1474:   *flg = PETSC_TRUE;
1475:   return(0);
1476: }

1478: static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
1479: {
1480:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1482:   PetscInt       i,n,len;
1483:   PetscScalar    *x,zero = 0.0;

1486:   VecSet(v,zero);
1487:   VecGetSize(v,&n);
1488:   VecGetArray(v,&x);
1489:   len  = PetscMin(A->rmap->n,A->cmap->n);
1490:   if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
1491:   for (i=0; i<len; i++) {
1492:     x[i] = mat->v[i*mat->lda + i];
1493:   }
1494:   VecRestoreArray(v,&x);
1495:   return(0);
1496: }

1498: static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
1499: {
1500:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1501:   const PetscScalar *l,*r;
1502:   PetscScalar       x,*v;
1503:   PetscErrorCode    ierr;
1504:   PetscInt          i,j,m = A->rmap->n,n = A->cmap->n;

1507:   if (ll) {
1508:     VecGetSize(ll,&m);
1509:     VecGetArrayRead(ll,&l);
1510:     if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size");
1511:     for (i=0; i<m; i++) {
1512:       x = l[i];
1513:       v = mat->v + i;
1514:       for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;}
1515:     }
1516:     VecRestoreArrayRead(ll,&l);
1517:     PetscLogFlops(1.0*n*m);
1518:   }
1519:   if (rr) {
1520:     VecGetSize(rr,&n);
1521:     VecGetArrayRead(rr,&r);
1522:     if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size");
1523:     for (i=0; i<n; i++) {
1524:       x = r[i];
1525:       v = mat->v + i*mat->lda;
1526:       for (j=0; j<m; j++) (*v++) *= x;
1527:     }
1528:     VecRestoreArrayRead(rr,&r);
1529:     PetscLogFlops(1.0*n*m);
1530:   }
1531:   return(0);
1532: }

1534: static PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1535: {
1536:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1537:   PetscScalar    *v   = mat->v;
1538:   PetscReal      sum  = 0.0;
1539:   PetscInt       lda  =mat->lda,m=A->rmap->n,i,j;

1543:   if (type == NORM_FROBENIUS) {
1544:     if (lda>m) {
1545:       for (j=0; j<A->cmap->n; j++) {
1546:         v = mat->v+j*lda;
1547:         for (i=0; i<m; i++) {
1548:           sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1549:         }
1550:       }
1551:     } else {
1552: #if defined(PETSC_USE_REAL___FP16)
1553:       PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n;
1554:       *nrm = BLASnrm2_(&cnt,v,&one);
1555:     }
1556: #else
1557:       for (i=0; i<A->cmap->n*A->rmap->n; i++) {
1558:         sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1559:       }
1560:     }
1561:     *nrm = PetscSqrtReal(sum);
1562: #endif
1563:     PetscLogFlops(2.0*A->cmap->n*A->rmap->n);
1564:   } else if (type == NORM_1) {
1565:     *nrm = 0.0;
1566:     for (j=0; j<A->cmap->n; j++) {
1567:       v   = mat->v + j*mat->lda;
1568:       sum = 0.0;
1569:       for (i=0; i<A->rmap->n; i++) {
1570:         sum += PetscAbsScalar(*v);  v++;
1571:       }
1572:       if (sum > *nrm) *nrm = sum;
1573:     }
1574:     PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1575:   } else if (type == NORM_INFINITY) {
1576:     *nrm = 0.0;
1577:     for (j=0; j<A->rmap->n; j++) {
1578:       v   = mat->v + j;
1579:       sum = 0.0;
1580:       for (i=0; i<A->cmap->n; i++) {
1581:         sum += PetscAbsScalar(*v); v += mat->lda;
1582:       }
1583:       if (sum > *nrm) *nrm = sum;
1584:     }
1585:     PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1586:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
1587:   return(0);
1588: }

1590: static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
1591: {
1592:   Mat_SeqDense   *aij = (Mat_SeqDense*)A->data;

1596:   switch (op) {
1597:   case MAT_ROW_ORIENTED:
1598:     aij->roworiented = flg;
1599:     break;
1600:   case MAT_NEW_NONZERO_LOCATIONS:
1601:   case MAT_NEW_NONZERO_LOCATION_ERR:
1602:   case MAT_NEW_NONZERO_ALLOCATION_ERR:
1603:   case MAT_NEW_DIAGONALS:
1604:   case MAT_KEEP_NONZERO_PATTERN:
1605:   case MAT_IGNORE_OFF_PROC_ENTRIES:
1606:   case MAT_USE_HASH_TABLE:
1607:   case MAT_IGNORE_ZERO_ENTRIES:
1608:   case MAT_IGNORE_LOWER_TRIANGULAR:
1609:     PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
1610:     break;
1611:   case MAT_SPD:
1612:   case MAT_SYMMETRIC:
1613:   case MAT_STRUCTURALLY_SYMMETRIC:
1614:   case MAT_HERMITIAN:
1615:   case MAT_SYMMETRY_ETERNAL:
1616:     /* These options are handled directly by MatSetOption() */
1617:     break;
1618:   default:
1619:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]);
1620:   }
1621:   return(0);
1622: }

1624: static PetscErrorCode MatZeroEntries_SeqDense(Mat A)
1625: {
1626:   Mat_SeqDense   *l = (Mat_SeqDense*)A->data;
1628:   PetscInt       lda=l->lda,m=A->rmap->n,j;

1631:   if (lda>m) {
1632:     for (j=0; j<A->cmap->n; j++) {
1633:       PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));
1634:     }
1635:   } else {
1636:     PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1637:   }
1638:   return(0);
1639: }

1641: static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1642: {
1643:   PetscErrorCode    ierr;
1644:   Mat_SeqDense      *l = (Mat_SeqDense*)A->data;
1645:   PetscInt          m  = l->lda, n = A->cmap->n, i,j;
1646:   PetscScalar       *slot,*bb;
1647:   const PetscScalar *xx;

1650: #if defined(PETSC_USE_DEBUG)
1651:   for (i=0; i<N; i++) {
1652:     if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
1653:     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);
1654:   }
1655: #endif

1657:   /* fix right hand side if needed */
1658:   if (x && b) {
1659:     VecGetArrayRead(x,&xx);
1660:     VecGetArray(b,&bb);
1661:     for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
1662:     VecRestoreArrayRead(x,&xx);
1663:     VecRestoreArray(b,&bb);
1664:   }

1666:   for (i=0; i<N; i++) {
1667:     slot = l->v + rows[i];
1668:     for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
1669:   }
1670:   if (diag != 0.0) {
1671:     if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
1672:     for (i=0; i<N; i++) {
1673:       slot  = l->v + (m+1)*rows[i];
1674:       *slot = diag;
1675:     }
1676:   }
1677:   return(0);
1678: }

1680: static PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[])
1681: {
1682:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;

1685:   if (mat->lda != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get array for Dense matrices with LDA different from number of rows");
1686:   *array = mat->v;
1687:   return(0);
1688: }

1690: static PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[])
1691: {
1693:   *array = 0; /* user cannot accidently use the array later */
1694:   return(0);
1695: }

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

1700:    Not Collective

1702:    Input Parameter:
1703: .  mat - a MATSEQDENSE or MATMPIDENSE matrix

1705:    Output Parameter:
1706: .   array - pointer to the data

1708:    Level: intermediate

1710: .seealso: MatDenseRestoreArray()
1711: @*/
1712: PetscErrorCode  MatDenseGetArray(Mat A,PetscScalar **array)
1713: {

1717:   PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));
1718:   return(0);
1719: }

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

1724:    Not Collective

1726:    Input Parameters:
1727: .  mat - a MATSEQDENSE or MATMPIDENSE matrix
1728: .  array - pointer to the data

1730:    Level: intermediate

1732: .seealso: MatDenseGetArray()
1733: @*/
1734: PetscErrorCode  MatDenseRestoreArray(Mat A,PetscScalar **array)
1735: {

1739:   PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));
1740:   return(0);
1741: }

1743: static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B)
1744: {
1745:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1747:   PetscInt       i,j,nrows,ncols;
1748:   const PetscInt *irow,*icol;
1749:   PetscScalar    *av,*bv,*v = mat->v;
1750:   Mat            newmat;

1753:   ISGetIndices(isrow,&irow);
1754:   ISGetIndices(iscol,&icol);
1755:   ISGetLocalSize(isrow,&nrows);
1756:   ISGetLocalSize(iscol,&ncols);

1758:   /* Check submatrixcall */
1759:   if (scall == MAT_REUSE_MATRIX) {
1760:     PetscInt n_cols,n_rows;
1761:     MatGetSize(*B,&n_rows,&n_cols);
1762:     if (n_rows != nrows || n_cols != ncols) {
1763:       /* resize the result matrix to match number of requested rows/columns */
1764:       MatSetSizes(*B,nrows,ncols,nrows,ncols);
1765:     }
1766:     newmat = *B;
1767:   } else {
1768:     /* Create and fill new matrix */
1769:     MatCreate(PetscObjectComm((PetscObject)A),&newmat);
1770:     MatSetSizes(newmat,nrows,ncols,nrows,ncols);
1771:     MatSetType(newmat,((PetscObject)A)->type_name);
1772:     MatSeqDenseSetPreallocation(newmat,NULL);
1773:   }

1775:   /* Now extract the data pointers and do the copy,column at a time */
1776:   bv = ((Mat_SeqDense*)newmat->data)->v;

1778:   for (i=0; i<ncols; i++) {
1779:     av = v + mat->lda*icol[i];
1780:     for (j=0; j<nrows; j++) *bv++ = av[irow[j]];
1781:   }

1783:   /* Assemble the matrices so that the correct flags are set */
1784:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1785:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);

1787:   /* Free work space */
1788:   ISRestoreIndices(isrow,&irow);
1789:   ISRestoreIndices(iscol,&icol);
1790:   *B   = newmat;
1791:   return(0);
1792: }

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

1800:   if (scall == MAT_INITIAL_MATRIX) {
1801:     PetscCalloc1(n+1,B);
1802:   }

1804:   for (i=0; i<n; i++) {
1805:     MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
1806:   }
1807:   return(0);
1808: }

1810: static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
1811: {
1813:   return(0);
1814: }

1816: static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
1817: {
1819:   return(0);
1820: }

1822: static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
1823: {
1824:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data;
1826:   PetscInt       lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j;

1829:   /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
1830:   if (A->ops->copy != B->ops->copy) {
1831:     MatCopy_Basic(A,B,str);
1832:     return(0);
1833:   }
1834:   if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)");
1835:   if (lda1>m || lda2>m) {
1836:     for (j=0; j<n; j++) {
1837:       PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));
1838:     }
1839:   } else {
1840:     PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1841:   }
1842:   PetscObjectStateIncrease((PetscObject)B);
1843:   return(0);
1844: }

1846: static PetscErrorCode MatSetUp_SeqDense(Mat A)
1847: {

1851:    MatSeqDenseSetPreallocation(A,0);
1852:   return(0);
1853: }

1855: static PetscErrorCode MatConjugate_SeqDense(Mat A)
1856: {
1857:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1858:   PetscInt     i,nz = A->rmap->n*A->cmap->n;
1859:   PetscScalar  *aa = a->v;

1862:   for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
1863:   return(0);
1864: }

1866: static PetscErrorCode MatRealPart_SeqDense(Mat A)
1867: {
1868:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1869:   PetscInt     i,nz = A->rmap->n*A->cmap->n;
1870:   PetscScalar  *aa = a->v;

1873:   for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
1874:   return(0);
1875: }

1877: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
1878: {
1879:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1880:   PetscInt     i,nz = A->rmap->n*A->cmap->n;
1881:   PetscScalar  *aa = a->v;

1884:   for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
1885:   return(0);
1886: }

1888: /* ----------------------------------------------------------------*/
1889: PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1890: {

1894:   if (scall == MAT_INITIAL_MATRIX) {
1895:     PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
1896:     MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
1897:     PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
1898:   }
1899:   PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
1900:   MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);
1901:   PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
1902:   return(0);
1903: }

1905: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1906: {
1908:   PetscInt       m=A->rmap->n,n=B->cmap->n;
1909:   Mat            Cmat;

1912:   if (A->cmap->n != B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"A->cmap->n %d != B->rmap->n %d\n",A->cmap->n,B->rmap->n);
1913:   MatCreate(PETSC_COMM_SELF,&Cmat);
1914:   MatSetSizes(Cmat,m,n,m,n);
1915:   MatSetType(Cmat,MATSEQDENSE);
1916:   MatSeqDenseSetPreallocation(Cmat,NULL);

1918:   *C = Cmat;
1919:   return(0);
1920: }

1922: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
1923: {
1924:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
1925:   Mat_SeqDense   *b = (Mat_SeqDense*)B->data;
1926:   Mat_SeqDense   *c = (Mat_SeqDense*)C->data;
1927:   PetscBLASInt   m,n,k;
1928:   PetscScalar    _DOne=1.0,_DZero=0.0;
1930:   PetscBool      flg;

1933:   PetscObjectTypeCompare((PetscObject)B,MATSEQDENSE,&flg);
1934:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Second matrix must be dense");

1936:   /* Handle case where where user provided the final C matrix rather than calling MatMatMult() with MAT_INITIAL_MATRIX*/
1937:   PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
1938:   if (flg) {
1939:     C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqDense;
1940:     (*C->ops->matmultnumeric)(A,B,C);
1941:     return(0);
1942:   }

1944:   PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&flg);
1945:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"First matrix must be dense");
1946:   PetscBLASIntCast(C->rmap->n,&m);
1947:   PetscBLASIntCast(C->cmap->n,&n);
1948:   PetscBLASIntCast(A->cmap->n,&k);
1949:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
1950:   return(0);
1951: }

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

1958:   if (scall == MAT_INITIAL_MATRIX) {
1959:     PetscLogEventBegin(MAT_MatTransposeMultSymbolic,A,B,0,0);
1960:     MatMatTransposeMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
1961:     PetscLogEventEnd(MAT_MatTransposeMultSymbolic,A,B,0,0);
1962:   }
1963:   PetscLogEventBegin(MAT_MatTransposeMultNumeric,A,B,0,0);
1964:   MatMatTransposeMultNumeric_SeqDense_SeqDense(A,B,*C);
1965:   PetscLogEventEnd(MAT_MatTransposeMultNumeric,A,B,0,0);
1966:   return(0);
1967: }

1969: PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1970: {
1972:   PetscInt       m=A->rmap->n,n=B->rmap->n;
1973:   Mat            Cmat;

1976:   if (A->cmap->n != B->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"A->cmap->n %d != B->cmap->n %d\n",A->cmap->n,B->cmap->n);
1977:   MatCreate(PETSC_COMM_SELF,&Cmat);
1978:   MatSetSizes(Cmat,m,n,m,n);
1979:   MatSetType(Cmat,MATSEQDENSE);
1980:   MatSeqDenseSetPreallocation(Cmat,NULL);

1982:   Cmat->assembled = PETSC_TRUE;

1984:   *C = Cmat;
1985:   return(0);
1986: }

1988: PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
1989: {
1990:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
1991:   Mat_SeqDense   *b = (Mat_SeqDense*)B->data;
1992:   Mat_SeqDense   *c = (Mat_SeqDense*)C->data;
1993:   PetscBLASInt   m,n,k;
1994:   PetscScalar    _DOne=1.0,_DZero=0.0;

1998:   PetscBLASIntCast(A->rmap->n,&m);
1999:   PetscBLASIntCast(B->rmap->n,&n);
2000:   PetscBLASIntCast(A->cmap->n,&k);
2001:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
2002:   return(0);
2003: }

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

2010:   if (scall == MAT_INITIAL_MATRIX) {
2011:     PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);
2012:     MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
2013:     PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);
2014:   }
2015:   PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);
2016:   MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);
2017:   PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);
2018:   return(0);
2019: }

2021: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
2022: {
2024:   PetscInt       m=A->cmap->n,n=B->cmap->n;
2025:   Mat            Cmat;

2028:   if (A->rmap->n != B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"A->rmap->n %d != B->rmap->n %d\n",A->rmap->n,B->rmap->n);
2029:   MatCreate(PETSC_COMM_SELF,&Cmat);
2030:   MatSetSizes(Cmat,m,n,m,n);
2031:   MatSetType(Cmat,MATSEQDENSE);
2032:   MatSeqDenseSetPreallocation(Cmat,NULL);

2034:   Cmat->assembled = PETSC_TRUE;

2036:   *C = Cmat;
2037:   return(0);
2038: }

2040: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2041: {
2042:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2043:   Mat_SeqDense   *b = (Mat_SeqDense*)B->data;
2044:   Mat_SeqDense   *c = (Mat_SeqDense*)C->data;
2045:   PetscBLASInt   m,n,k;
2046:   PetscScalar    _DOne=1.0,_DZero=0.0;

2050:   PetscBLASIntCast(C->rmap->n,&m);
2051:   PetscBLASIntCast(C->cmap->n,&n);
2052:   PetscBLASIntCast(A->rmap->n,&k);
2053:   PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
2054:   return(0);
2055: }

2057: static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[])
2058: {
2059:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2061:   PetscInt       i,j,m = A->rmap->n,n = A->cmap->n,p;
2062:   PetscScalar    *x;
2063:   MatScalar      *aa = a->v;

2066:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

2068:   VecSet(v,0.0);
2069:   VecGetArray(v,&x);
2070:   VecGetLocalSize(v,&p);
2071:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2072:   for (i=0; i<m; i++) {
2073:     x[i] = aa[i]; if (idx) idx[i] = 0;
2074:     for (j=1; j<n; j++) {
2075:       if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
2076:     }
2077:   }
2078:   VecRestoreArray(v,&x);
2079:   return(0);
2080: }

2082: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])
2083: {
2084:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2086:   PetscInt       i,j,m = A->rmap->n,n = A->cmap->n,p;
2087:   PetscScalar    *x;
2088:   PetscReal      atmp;
2089:   MatScalar      *aa = a->v;

2092:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

2094:   VecSet(v,0.0);
2095:   VecGetArray(v,&x);
2096:   VecGetLocalSize(v,&p);
2097:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2098:   for (i=0; i<m; i++) {
2099:     x[i] = PetscAbsScalar(aa[i]);
2100:     for (j=1; j<n; j++) {
2101:       atmp = PetscAbsScalar(aa[i+m*j]);
2102:       if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;}
2103:     }
2104:   }
2105:   VecRestoreArray(v,&x);
2106:   return(0);
2107: }

2109: static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])
2110: {
2111:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2113:   PetscInt       i,j,m = A->rmap->n,n = A->cmap->n,p;
2114:   PetscScalar    *x;
2115:   MatScalar      *aa = a->v;

2118:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

2120:   VecSet(v,0.0);
2121:   VecGetArray(v,&x);
2122:   VecGetLocalSize(v,&p);
2123:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2124:   for (i=0; i<m; i++) {
2125:     x[i] = aa[i]; if (idx) idx[i] = 0;
2126:     for (j=1; j<n; j++) {
2127:       if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
2128:     }
2129:   }
2130:   VecRestoreArray(v,&x);
2131:   return(0);
2132: }

2134: static PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
2135: {
2136:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2138:   PetscScalar    *x;

2141:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

2143:   VecGetArray(v,&x);
2144:   PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));
2145:   VecRestoreArray(v,&x);
2146:   return(0);
2147: }


2150: PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms)
2151: {
2153:   PetscInt       i,j,m,n;
2154:   PetscScalar    *a;

2157:   MatGetSize(A,&m,&n);
2158:   PetscMemzero(norms,n*sizeof(PetscReal));
2159:   MatDenseGetArray(A,&a);
2160:   if (type == NORM_2) {
2161:     for (i=0; i<n; i++) {
2162:       for (j=0; j<m; j++) {
2163:         norms[i] += PetscAbsScalar(a[j]*a[j]);
2164:       }
2165:       a += m;
2166:     }
2167:   } else if (type == NORM_1) {
2168:     for (i=0; i<n; i++) {
2169:       for (j=0; j<m; j++) {
2170:         norms[i] += PetscAbsScalar(a[j]);
2171:       }
2172:       a += m;
2173:     }
2174:   } else if (type == NORM_INFINITY) {
2175:     for (i=0; i<n; i++) {
2176:       for (j=0; j<m; j++) {
2177:         norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]);
2178:       }
2179:       a += m;
2180:     }
2181:   } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType");
2182:   MatDenseRestoreArray(A,&a);
2183:   if (type == NORM_2) {
2184:     for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]);
2185:   }
2186:   return(0);
2187: }

2189: static PetscErrorCode  MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
2190: {
2192:   PetscScalar    *a;
2193:   PetscInt       m,n,i;

2196:   MatGetSize(x,&m,&n);
2197:   MatDenseGetArray(x,&a);
2198:   for (i=0; i<m*n; i++) {
2199:     PetscRandomGetValue(rctx,a+i);
2200:   }
2201:   MatDenseRestoreArray(x,&a);
2202:   return(0);
2203: }

2205: static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool  *missing,PetscInt *d)
2206: {
2208:   *missing = PETSC_FALSE;
2209:   return(0);
2210: }


2213: /* -------------------------------------------------------------------*/
2214: static struct _MatOps MatOps_Values = { MatSetValues_SeqDense,
2215:                                         MatGetRow_SeqDense,
2216:                                         MatRestoreRow_SeqDense,
2217:                                         MatMult_SeqDense,
2218:                                 /*  4*/ MatMultAdd_SeqDense,
2219:                                         MatMultTranspose_SeqDense,
2220:                                         MatMultTransposeAdd_SeqDense,
2221:                                         0,
2222:                                         0,
2223:                                         0,
2224:                                 /* 10*/ 0,
2225:                                         MatLUFactor_SeqDense,
2226:                                         MatCholeskyFactor_SeqDense,
2227:                                         MatSOR_SeqDense,
2228:                                         MatTranspose_SeqDense,
2229:                                 /* 15*/ MatGetInfo_SeqDense,
2230:                                         MatEqual_SeqDense,
2231:                                         MatGetDiagonal_SeqDense,
2232:                                         MatDiagonalScale_SeqDense,
2233:                                         MatNorm_SeqDense,
2234:                                 /* 20*/ MatAssemblyBegin_SeqDense,
2235:                                         MatAssemblyEnd_SeqDense,
2236:                                         MatSetOption_SeqDense,
2237:                                         MatZeroEntries_SeqDense,
2238:                                 /* 24*/ MatZeroRows_SeqDense,
2239:                                         0,
2240:                                         0,
2241:                                         0,
2242:                                         0,
2243:                                 /* 29*/ MatSetUp_SeqDense,
2244:                                         0,
2245:                                         0,
2246:                                         0,
2247:                                         0,
2248:                                 /* 34*/ MatDuplicate_SeqDense,
2249:                                         0,
2250:                                         0,
2251:                                         0,
2252:                                         0,
2253:                                 /* 39*/ MatAXPY_SeqDense,
2254:                                         MatCreateSubMatrices_SeqDense,
2255:                                         0,
2256:                                         MatGetValues_SeqDense,
2257:                                         MatCopy_SeqDense,
2258:                                 /* 44*/ MatGetRowMax_SeqDense,
2259:                                         MatScale_SeqDense,
2260:                                         MatShift_Basic,
2261:                                         0,
2262:                                         MatZeroRowsColumns_SeqDense,
2263:                                 /* 49*/ MatSetRandom_SeqDense,
2264:                                         0,
2265:                                         0,
2266:                                         0,
2267:                                         0,
2268:                                 /* 54*/ 0,
2269:                                         0,
2270:                                         0,
2271:                                         0,
2272:                                         0,
2273:                                 /* 59*/ 0,
2274:                                         MatDestroy_SeqDense,
2275:                                         MatView_SeqDense,
2276:                                         0,
2277:                                         0,
2278:                                 /* 64*/ 0,
2279:                                         0,
2280:                                         0,
2281:                                         0,
2282:                                         0,
2283:                                 /* 69*/ MatGetRowMaxAbs_SeqDense,
2284:                                         0,
2285:                                         0,
2286:                                         0,
2287:                                         0,
2288:                                 /* 74*/ 0,
2289:                                         0,
2290:                                         0,
2291:                                         0,
2292:                                         0,
2293:                                 /* 79*/ 0,
2294:                                         0,
2295:                                         0,
2296:                                         0,
2297:                                 /* 83*/ MatLoad_SeqDense,
2298:                                         0,
2299:                                         MatIsHermitian_SeqDense,
2300:                                         0,
2301:                                         0,
2302:                                         0,
2303:                                 /* 89*/ MatMatMult_SeqDense_SeqDense,
2304:                                         MatMatMultSymbolic_SeqDense_SeqDense,
2305:                                         MatMatMultNumeric_SeqDense_SeqDense,
2306:                                         MatPtAP_SeqDense_SeqDense,
2307:                                         MatPtAPSymbolic_SeqDense_SeqDense,
2308:                                 /* 94*/ MatPtAPNumeric_SeqDense_SeqDense,
2309:                                         MatMatTransposeMult_SeqDense_SeqDense,
2310:                                         MatMatTransposeMultSymbolic_SeqDense_SeqDense,
2311:                                         MatMatTransposeMultNumeric_SeqDense_SeqDense,
2312:                                         0,
2313:                                 /* 99*/ 0,
2314:                                         0,
2315:                                         0,
2316:                                         MatConjugate_SeqDense,
2317:                                         0,
2318:                                 /*104*/ 0,
2319:                                         MatRealPart_SeqDense,
2320:                                         MatImaginaryPart_SeqDense,
2321:                                         0,
2322:                                         0,
2323:                                 /*109*/ 0,
2324:                                         0,
2325:                                         MatGetRowMin_SeqDense,
2326:                                         MatGetColumnVector_SeqDense,
2327:                                         MatMissingDiagonal_SeqDense,
2328:                                 /*114*/ 0,
2329:                                         0,
2330:                                         0,
2331:                                         0,
2332:                                         0,
2333:                                 /*119*/ 0,
2334:                                         0,
2335:                                         0,
2336:                                         0,
2337:                                         0,
2338:                                 /*124*/ 0,
2339:                                         MatGetColumnNorms_SeqDense,
2340:                                         0,
2341:                                         0,
2342:                                         0,
2343:                                 /*129*/ 0,
2344:                                         MatTransposeMatMult_SeqDense_SeqDense,
2345:                                         MatTransposeMatMultSymbolic_SeqDense_SeqDense,
2346:                                         MatTransposeMatMultNumeric_SeqDense_SeqDense,
2347:                                         0,
2348:                                 /*134*/ 0,
2349:                                         0,
2350:                                         0,
2351:                                         0,
2352:                                         0,
2353:                                 /*139*/ 0,
2354:                                         0,
2355:                                         0
2356: };

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

2363:    Collective on MPI_Comm

2365:    Input Parameters:
2366: +  comm - MPI communicator, set to PETSC_COMM_SELF
2367: .  m - number of rows
2368: .  n - number of columns
2369: -  data - optional location of matrix data in column major order.  Set data=NULL for PETSc
2370:    to control all matrix memory allocation.

2372:    Output Parameter:
2373: .  A - the matrix

2375:    Notes:
2376:    The data input variable is intended primarily for Fortran programmers
2377:    who wish to allocate their own matrix memory space.  Most users should
2378:    set data=NULL.

2380:    Level: intermediate

2382: .keywords: dense, matrix, LAPACK, BLAS

2384: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
2385: @*/
2386: PetscErrorCode  MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
2387: {

2391:   MatCreate(comm,A);
2392:   MatSetSizes(*A,m,n,m,n);
2393:   MatSetType(*A,MATSEQDENSE);
2394:   MatSeqDenseSetPreallocation(*A,data);
2395:   return(0);
2396: }

2398: /*@C
2399:    MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements

2401:    Collective on MPI_Comm

2403:    Input Parameters:
2404: +  B - the matrix
2405: -  data - the array (or NULL)

2407:    Notes:
2408:    The data input variable is intended primarily for Fortran programmers
2409:    who wish to allocate their own matrix memory space.  Most users should
2410:    need not call this routine.

2412:    Level: intermediate

2414: .keywords: dense, matrix, LAPACK, BLAS

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

2418: @*/
2419: PetscErrorCode  MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
2420: {

2424:   PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
2425:   return(0);
2426: }

2428: PetscErrorCode  MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
2429: {
2430:   Mat_SeqDense   *b;

2434:   B->preallocated = PETSC_TRUE;

2436:   PetscLayoutSetUp(B->rmap);
2437:   PetscLayoutSetUp(B->cmap);

2439:   b       = (Mat_SeqDense*)B->data;
2440:   b->Mmax = B->rmap->n;
2441:   b->Nmax = B->cmap->n;
2442:   if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n;

2444:   PetscIntMultError(b->lda,b->Nmax,NULL);
2445:   if (!data) { /* petsc-allocated storage */
2446:     if (!b->user_alloc) { PetscFree(b->v); }
2447:     PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);
2448:     PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));

2450:     b->user_alloc = PETSC_FALSE;
2451:   } else { /* user-allocated storage */
2452:     if (!b->user_alloc) { PetscFree(b->v); }
2453:     b->v          = data;
2454:     b->user_alloc = PETSC_TRUE;
2455:   }
2456:   B->assembled = PETSC_TRUE;
2457:   return(0);
2458: }

2460: #if defined(PETSC_HAVE_ELEMENTAL)
2461: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
2462: {
2463:   Mat            mat_elemental;
2465:   PetscScalar    *array,*v_colwise;
2466:   PetscInt       M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols;

2469:   PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);
2470:   MatDenseGetArray(A,&array);
2471:   /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
2472:   k = 0;
2473:   for (j=0; j<N; j++) {
2474:     cols[j] = j;
2475:     for (i=0; i<M; i++) {
2476:       v_colwise[j*M+i] = array[k++];
2477:     }
2478:   }
2479:   for (i=0; i<M; i++) {
2480:     rows[i] = i;
2481:   }
2482:   MatDenseRestoreArray(A,&array);

2484:   MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
2485:   MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);
2486:   MatSetType(mat_elemental,MATELEMENTAL);
2487:   MatSetUp(mat_elemental);

2489:   /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
2490:   MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);
2491:   MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
2492:   MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);
2493:   PetscFree3(v_colwise,rows,cols);

2495:   if (reuse == MAT_INPLACE_MATRIX) {
2496:     MatHeaderReplace(A,&mat_elemental);
2497:   } else {
2498:     *newmat = mat_elemental;
2499:   }
2500:   return(0);
2501: }
2502: #endif

2504: /*@C
2505:   MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array

2507:   Input parameter:
2508: + A - the matrix
2509: - lda - the leading dimension

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

2516:   Level: intermediate

2518: .keywords: dense, matrix, LAPACK, BLAS

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

2522: @*/
2523: PetscErrorCode  MatSeqDenseSetLDA(Mat B,PetscInt lda)
2524: {
2525:   Mat_SeqDense *b = (Mat_SeqDense*)B->data;

2528:   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);
2529:   b->lda       = lda;
2530:   b->changelda = PETSC_FALSE;
2531:   b->Mmax      = PetscMax(b->Mmax,lda);
2532:   return(0);
2533: }

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

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

2541:   Level: beginner

2543: .seealso: MatCreateSeqDense()

2545: M*/

2547: PETSC_EXTERN PetscErrorCode MatCreate_SeqDense(Mat B)
2548: {
2549:   Mat_SeqDense   *b;
2551:   PetscMPIInt    size;

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

2557:   PetscNewLog(B,&b);
2558:   PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
2559:   B->data = (void*)b;

2561:   b->roworiented = PETSC_TRUE;

2563:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);
2564:   PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);
2565:   PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);
2566:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);
2567:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);
2568: #if defined(PETSC_HAVE_ELEMENTAL)
2569:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);
2570: #endif
2571:   PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);
2572:   PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2573:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2574:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2575:   PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaij_seqdense_C",MatPtAP_SeqDense_SeqDense);
2576:   PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijperm_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2577:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijperm_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2578:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijperm_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2579:   PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijperm_seqdense_C",MatPtAP_SeqDense_SeqDense);
2580:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2581:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2582:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2583:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijperm_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2584:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijperm_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2585:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijperm_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2586:   PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
2587:   return(0);
2588: }