Actual source code: dense.c

petsc-master 2017-11-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:   VecGetArrayRead(xx,&x);
816:   VecGetArray(yy,&y);
817:   if (!A->rmap->n || !A->cmap->n) {
818:     PetscBLASInt i;
819:     for (i=0; i<n; i++) y[i] = 0.0;
820:   } else {
821:     PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One));
822:     PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);
823:   }
824:   VecRestoreArrayRead(xx,&x);
825:   VecRestoreArray(yy,&y);
826:   return(0);
827: }

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

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

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

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

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

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

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

908:   *ncols = A->cmap->n;
909:   if (cols) {
910:     PetscMalloc1(A->cmap->n+1,cols);
911:     for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
912:   }
913:   if (vals) {
914:     PetscMalloc1(A->cmap->n+1,vals);
915:     v    = mat->v + row;
916:     for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
917:   }
918:   return(0);
919: }

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

926:   if (cols) {PetscFree(*cols);}
927:   if (vals) {PetscFree(*vals); }
928:   return(0);
929: }
930: /* ----------------------------------------------------------------*/
931: static PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
932: {
933:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
934:   PetscInt     i,j,idx=0;

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

1001: static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
1002: {
1003:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1004:   PetscInt     i,j;

1007:   /* row-oriented output */
1008:   for (i=0; i<m; i++) {
1009:     if (indexm[i] < 0) {v += n;continue;}
1010:     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);
1011:     for (j=0; j<n; j++) {
1012:       if (indexn[j] < 0) {v++; continue;}
1013:       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);
1014:       *v++ = mat->v[indexn[j]*mat->lda + indexm[i]];
1015:     }
1016:   }
1017:   return(0);
1018: }

1020: /* -----------------------------------------------------------------*/

1022: static PetscErrorCode MatLoad_SeqDense(Mat newmat,PetscViewer viewer)
1023: {
1024:   Mat_SeqDense   *a;
1026:   PetscInt       *scols,i,j,nz,header[4];
1027:   int            fd;
1028:   PetscMPIInt    size;
1029:   PetscInt       *rowlengths = 0,M,N,*cols,grows,gcols;
1030:   PetscScalar    *vals,*svals,*v,*w;
1031:   MPI_Comm       comm;

1034:   /* force binary viewer to load .info file if it has not yet done so */
1035:   PetscViewerSetUp(viewer);
1036:   PetscObjectGetComm((PetscObject)viewer,&comm);
1037:   MPI_Comm_size(comm,&size);
1038:   if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor");
1039:   PetscViewerBinaryGetDescriptor(viewer,&fd);
1040:   PetscBinaryRead(fd,header,4,PETSC_INT);
1041:   if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not matrix object");
1042:   M = header[1]; N = header[2]; nz = header[3];

1044:   /* set global size if not set already*/
1045:   if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) {
1046:     MatSetSizes(newmat,M,N,M,N);
1047:   } else {
1048:     /* if sizes and type are already set, check if the vector global sizes are correct */
1049:     MatGetSize(newmat,&grows,&gcols);
1050:     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);
1051:   }
1052:   a = (Mat_SeqDense*)newmat->data;
1053:   if (!a->user_alloc) {
1054:     MatSeqDenseSetPreallocation(newmat,NULL);
1055:   }

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

1077:     a = (Mat_SeqDense*)newmat->data;
1078:     v = a->v;

1080:     /* read column indices and nonzeros */
1081:     PetscMalloc1(nz+1,&scols);
1082:     cols = scols;
1083:     PetscBinaryRead(fd,cols,nz,PETSC_INT);
1084:     PetscMalloc1(nz+1,&svals);
1085:     vals = svals;
1086:     PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);

1088:     /* insert into matrix */
1089:     for (i=0; i<M; i++) {
1090:       for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j];
1091:       svals += rowlengths[i]; scols += rowlengths[i];
1092:     }
1093:     PetscFree(vals);
1094:     PetscFree(cols);
1095:     PetscFree(rowlengths);
1096:   }
1097:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1098:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
1099:   return(0);
1100: }

1102: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
1103: {
1104:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
1105:   PetscErrorCode    ierr;
1106:   PetscInt          i,j;
1107:   const char        *name;
1108:   PetscScalar       *v;
1109:   PetscViewerFormat format;
1110: #if defined(PETSC_USE_COMPLEX)
1111:   PetscBool         allreal = PETSC_TRUE;
1112: #endif

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

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

1181: static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer)
1182: {
1183:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
1184:   PetscErrorCode    ierr;
1185:   int               fd;
1186:   PetscInt          ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n;
1187:   PetscScalar       *v,*anonz,*vals;
1188:   PetscViewerFormat format;

1191:   PetscViewerBinaryGetDescriptor(viewer,&fd);

1193:   PetscViewerGetFormat(viewer,&format);
1194:   if (format == PETSC_VIEWER_NATIVE) {
1195:     /* store the matrix as a dense matrix */
1196:     PetscMalloc1(4,&col_lens);

1198:     col_lens[0] = MAT_FILE_CLASSID;
1199:     col_lens[1] = m;
1200:     col_lens[2] = n;
1201:     col_lens[3] = MATRIX_BINARY_FORMAT_DENSE;

1203:     PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);
1204:     PetscFree(col_lens);

1206:     /* write out matrix, by rows */
1207:     PetscMalloc1(m*n+1,&vals);
1208:     v    = a->v;
1209:     for (j=0; j<n; j++) {
1210:       for (i=0; i<m; i++) {
1211:         vals[j + i*n] = *v++;
1212:       }
1213:     }
1214:     PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);
1215:     PetscFree(vals);
1216:   } else {
1217:     PetscMalloc1(4+nz,&col_lens);

1219:     col_lens[0] = MAT_FILE_CLASSID;
1220:     col_lens[1] = m;
1221:     col_lens[2] = n;
1222:     col_lens[3] = nz;

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

1228:     /* Possibly should write in smaller increments, not whole matrix at once? */
1229:     /* store column indices (zero start index) */
1230:     ict = 0;
1231:     for (i=0; i<m; i++) {
1232:       for (j=0; j<n; j++) col_lens[ict++] = j;
1233:     }
1234:     PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);
1235:     PetscFree(col_lens);

1237:     /* store nonzero values */
1238:     PetscMalloc1(nz+1,&anonz);
1239:     ict  = 0;
1240:     for (i=0; i<m; i++) {
1241:       v = a->v + i;
1242:       for (j=0; j<n; j++) {
1243:         anonz[ict++] = *v; v += a->lda;
1244:       }
1245:     }
1246:     PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);
1247:     PetscFree(anonz);
1248:   }
1249:   return(0);
1250: }

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

1266:   PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
1267:   PetscViewerGetFormat(viewer,&format);
1268:   PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);

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

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

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

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

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

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

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

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

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

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

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

1367:   a->unplacedarray       = a->v;
1368:   a->unplaced_user_alloc = a->user_alloc;
1369:   a->v                   = (PetscScalar*) array;
1370:   return(0);
1371: }

1373: static PetscErrorCode MatDenseResetArray_SeqDense(Mat A)
1374: {
1375:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

1378:   a->v             = a->unplacedarray;
1379:   a->user_alloc    = a->unplaced_user_alloc;
1380:   a->unplacedarray = NULL;
1381:   return(0);
1382: }

1384: static PetscErrorCode MatDestroy_SeqDense(Mat mat)
1385: {
1386:   Mat_SeqDense   *l = (Mat_SeqDense*)mat->data;

1390: #if defined(PETSC_USE_LOG)
1391:   PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n);
1392: #endif
1393:   PetscFree(l->pivots);
1394:   PetscFree(l->fwork);
1395:   MatDestroy(&l->ptapwork);
1396:   if (!l->user_alloc) {PetscFree(l->v);}
1397:   PetscFree(mat->data);

1399:   PetscObjectChangeTypeName((PetscObject)mat,0);
1400:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);
1401:   PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL);
1402:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL);
1403:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);
1404:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);
1405: #if defined(PETSC_HAVE_ELEMENTAL)
1406:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);
1407: #endif
1408:   PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);
1409:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);
1410:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);
1411:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);
1412:   PetscObjectComposeFunction((PetscObject)mat,"MatPtAP_seqaij_seqdense_C",NULL);
1413:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);
1414:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);
1415:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);
1416:   return(0);
1417: }

1419: static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
1420: {
1421:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1423:   PetscInt       k,j,m,n,M;
1424:   PetscScalar    *v,tmp;

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

1445:     if (reuse == MAT_INITIAL_MATRIX) {
1446:       MatCreate(PetscObjectComm((PetscObject)A),&tmat);
1447:       MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);
1448:       MatSetType(tmat,((PetscObject)A)->type_name);
1449:       MatSeqDenseSetPreallocation(tmat,NULL);
1450:     } else {
1451:       tmat = *matout;
1452:     }
1453:     tmatd = (Mat_SeqDense*)tmat->data;
1454:     v = mat->v; v2 = tmatd->v; M2 = tmatd->lda;
1455:     for (j=0; j<n; j++) {
1456:       for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M];
1457:     }
1458:     MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
1459:     MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
1460:     *matout = tmat;
1461:   }
1462:   return(0);
1463: }

1465: static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool  *flg)
1466: {
1467:   Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data;
1468:   Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data;
1469:   PetscInt     i,j;
1470:   PetscScalar  *v1,*v2;

1473:   if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; return(0);}
1474:   if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; return(0);}
1475:   for (i=0; i<A1->rmap->n; i++) {
1476:     v1 = mat1->v+i; v2 = mat2->v+i;
1477:     for (j=0; j<A1->cmap->n; j++) {
1478:       if (*v1 != *v2) {*flg = PETSC_FALSE; return(0);}
1479:       v1 += mat1->lda; v2 += mat2->lda;
1480:     }
1481:   }
1482:   *flg = PETSC_TRUE;
1483:   return(0);
1484: }

1486: static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
1487: {
1488:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1490:   PetscInt       i,n,len;
1491:   PetscScalar    *x,zero = 0.0;

1494:   VecSet(v,zero);
1495:   VecGetSize(v,&n);
1496:   VecGetArray(v,&x);
1497:   len  = PetscMin(A->rmap->n,A->cmap->n);
1498:   if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
1499:   for (i=0; i<len; i++) {
1500:     x[i] = mat->v[i*mat->lda + i];
1501:   }
1502:   VecRestoreArray(v,&x);
1503:   return(0);
1504: }

1506: static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
1507: {
1508:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1509:   const PetscScalar *l,*r;
1510:   PetscScalar       x,*v;
1511:   PetscErrorCode    ierr;
1512:   PetscInt          i,j,m = A->rmap->n,n = A->cmap->n;

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

1542: static PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1543: {
1544:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1545:   PetscScalar    *v   = mat->v;
1546:   PetscReal      sum  = 0.0;
1547:   PetscInt       lda  =mat->lda,m=A->rmap->n,i,j;

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

1598: static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
1599: {
1600:   Mat_SeqDense   *aij = (Mat_SeqDense*)A->data;

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

1632: static PetscErrorCode MatZeroEntries_SeqDense(Mat A)
1633: {
1634:   Mat_SeqDense   *l = (Mat_SeqDense*)A->data;
1636:   PetscInt       lda=l->lda,m=A->rmap->n,j;

1639:   if (lda>m) {
1640:     for (j=0; j<A->cmap->n; j++) {
1641:       PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));
1642:     }
1643:   } else {
1644:     PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1645:   }
1646:   return(0);
1647: }

1649: static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1650: {
1651:   PetscErrorCode    ierr;
1652:   Mat_SeqDense      *l = (Mat_SeqDense*)A->data;
1653:   PetscInt          m  = l->lda, n = A->cmap->n, i,j;
1654:   PetscScalar       *slot,*bb;
1655:   const PetscScalar *xx;

1658: #if defined(PETSC_USE_DEBUG)
1659:   for (i=0; i<N; i++) {
1660:     if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
1661:     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);
1662:   }
1663: #endif

1665:   /* fix right hand side if needed */
1666:   if (x && b) {
1667:     VecGetArrayRead(x,&xx);
1668:     VecGetArray(b,&bb);
1669:     for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
1670:     VecRestoreArrayRead(x,&xx);
1671:     VecRestoreArray(b,&bb);
1672:   }

1674:   for (i=0; i<N; i++) {
1675:     slot = l->v + rows[i];
1676:     for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
1677:   }
1678:   if (diag != 0.0) {
1679:     if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
1680:     for (i=0; i<N; i++) {
1681:       slot  = l->v + (m+1)*rows[i];
1682:       *slot = diag;
1683:     }
1684:   }
1685:   return(0);
1686: }

1688: static PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[])
1689: {
1690:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;

1693:   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");
1694:   *array = mat->v;
1695:   return(0);
1696: }

1698: static PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[])
1699: {
1701:   *array = 0; /* user cannot accidently use the array later */
1702:   return(0);
1703: }

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

1708:    Not Collective

1710:    Input Parameter:
1711: .  mat - a MATSEQDENSE or MATMPIDENSE matrix

1713:    Output Parameter:
1714: .   array - pointer to the data

1716:    Level: intermediate

1718: .seealso: MatDenseRestoreArray()
1719: @*/
1720: PetscErrorCode  MatDenseGetArray(Mat A,PetscScalar **array)
1721: {

1725:   PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));
1726:   return(0);
1727: }

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

1732:    Not Collective

1734:    Input Parameters:
1735: .  mat - a MATSEQDENSE or MATMPIDENSE matrix
1736: .  array - pointer to the data

1738:    Level: intermediate

1740: .seealso: MatDenseGetArray()
1741: @*/
1742: PetscErrorCode  MatDenseRestoreArray(Mat A,PetscScalar **array)
1743: {

1747:   PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));
1748:   return(0);
1749: }

1751: static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B)
1752: {
1753:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1755:   PetscInt       i,j,nrows,ncols;
1756:   const PetscInt *irow,*icol;
1757:   PetscScalar    *av,*bv,*v = mat->v;
1758:   Mat            newmat;

1761:   ISGetIndices(isrow,&irow);
1762:   ISGetIndices(iscol,&icol);
1763:   ISGetLocalSize(isrow,&nrows);
1764:   ISGetLocalSize(iscol,&ncols);

1766:   /* Check submatrixcall */
1767:   if (scall == MAT_REUSE_MATRIX) {
1768:     PetscInt n_cols,n_rows;
1769:     MatGetSize(*B,&n_rows,&n_cols);
1770:     if (n_rows != nrows || n_cols != ncols) {
1771:       /* resize the result matrix to match number of requested rows/columns */
1772:       MatSetSizes(*B,nrows,ncols,nrows,ncols);
1773:     }
1774:     newmat = *B;
1775:   } else {
1776:     /* Create and fill new matrix */
1777:     MatCreate(PetscObjectComm((PetscObject)A),&newmat);
1778:     MatSetSizes(newmat,nrows,ncols,nrows,ncols);
1779:     MatSetType(newmat,((PetscObject)A)->type_name);
1780:     MatSeqDenseSetPreallocation(newmat,NULL);
1781:   }

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

1786:   for (i=0; i<ncols; i++) {
1787:     av = v + mat->lda*icol[i];
1788:     for (j=0; j<nrows; j++) *bv++ = av[irow[j]];
1789:   }

1791:   /* Assemble the matrices so that the correct flags are set */
1792:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1793:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);

1795:   /* Free work space */
1796:   ISRestoreIndices(isrow,&irow);
1797:   ISRestoreIndices(iscol,&icol);
1798:   *B   = newmat;
1799:   return(0);
1800: }

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

1808:   if (scall == MAT_INITIAL_MATRIX) {
1809:     PetscCalloc1(n+1,B);
1810:   }

1812:   for (i=0; i<n; i++) {
1813:     MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
1814:   }
1815:   return(0);
1816: }

1818: static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
1819: {
1821:   return(0);
1822: }

1824: static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
1825: {
1827:   return(0);
1828: }

1830: static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
1831: {
1832:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data;
1834:   PetscInt       lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j;

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

1854: static PetscErrorCode MatSetUp_SeqDense(Mat A)
1855: {

1859:    MatSeqDenseSetPreallocation(A,0);
1860:   return(0);
1861: }

1863: static PetscErrorCode MatConjugate_SeqDense(Mat A)
1864: {
1865:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1866:   PetscInt     i,nz = A->rmap->n*A->cmap->n;
1867:   PetscScalar  *aa = a->v;

1870:   for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
1871:   return(0);
1872: }

1874: static PetscErrorCode MatRealPart_SeqDense(Mat A)
1875: {
1876:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1877:   PetscInt     i,nz = A->rmap->n*A->cmap->n;
1878:   PetscScalar  *aa = a->v;

1881:   for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
1882:   return(0);
1883: }

1885: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
1886: {
1887:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1888:   PetscInt     i,nz = A->rmap->n*A->cmap->n;
1889:   PetscScalar  *aa = a->v;

1892:   for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
1893:   return(0);
1894: }

1896: /* ----------------------------------------------------------------*/
1897: PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1898: {

1902:   if (scall == MAT_INITIAL_MATRIX) {
1903:     PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
1904:     MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
1905:     PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
1906:   }
1907:   PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
1908:   MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);
1909:   PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
1910:   return(0);
1911: }

1913: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1914: {
1916:   PetscInt       m=A->rmap->n,n=B->cmap->n;
1917:   Mat            Cmat;

1920:   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);
1921:   MatCreate(PETSC_COMM_SELF,&Cmat);
1922:   MatSetSizes(Cmat,m,n,m,n);
1923:   MatSetType(Cmat,MATSEQDENSE);
1924:   MatSeqDenseSetPreallocation(Cmat,NULL);

1926:   *C = Cmat;
1927:   return(0);
1928: }

1930: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
1931: {
1932:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
1933:   Mat_SeqDense   *b = (Mat_SeqDense*)B->data;
1934:   Mat_SeqDense   *c = (Mat_SeqDense*)C->data;
1935:   PetscBLASInt   m,n,k;
1936:   PetscScalar    _DOne=1.0,_DZero=0.0;
1938:   PetscBool      flg;

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

1944:   /* Handle case where where user provided the final C matrix rather than calling MatMatMult() with MAT_INITIAL_MATRIX*/
1945:   PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
1946:   if (flg) {
1947:     C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqDense;
1948:     (*C->ops->matmultnumeric)(A,B,C);
1949:     return(0);
1950:   }

1952:   PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&flg);
1953:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"First matrix must be dense");
1954:   PetscBLASIntCast(C->rmap->n,&m);
1955:   PetscBLASIntCast(C->cmap->n,&n);
1956:   PetscBLASIntCast(A->cmap->n,&k);
1957:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
1958:   return(0);
1959: }

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

1966:   if (scall == MAT_INITIAL_MATRIX) {
1967:     PetscLogEventBegin(MAT_MatTransposeMultSymbolic,A,B,0,0);
1968:     MatMatTransposeMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
1969:     PetscLogEventEnd(MAT_MatTransposeMultSymbolic,A,B,0,0);
1970:   }
1971:   PetscLogEventBegin(MAT_MatTransposeMultNumeric,A,B,0,0);
1972:   MatMatTransposeMultNumeric_SeqDense_SeqDense(A,B,*C);
1973:   PetscLogEventEnd(MAT_MatTransposeMultNumeric,A,B,0,0);
1974:   return(0);
1975: }

1977: PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1978: {
1980:   PetscInt       m=A->rmap->n,n=B->rmap->n;
1981:   Mat            Cmat;

1984:   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);
1985:   MatCreate(PETSC_COMM_SELF,&Cmat);
1986:   MatSetSizes(Cmat,m,n,m,n);
1987:   MatSetType(Cmat,MATSEQDENSE);
1988:   MatSeqDenseSetPreallocation(Cmat,NULL);

1990:   Cmat->assembled = PETSC_TRUE;

1992:   *C = Cmat;
1993:   return(0);
1994: }

1996: PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
1997: {
1998:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
1999:   Mat_SeqDense   *b = (Mat_SeqDense*)B->data;
2000:   Mat_SeqDense   *c = (Mat_SeqDense*)C->data;
2001:   PetscBLASInt   m,n,k;
2002:   PetscScalar    _DOne=1.0,_DZero=0.0;

2006:   PetscBLASIntCast(A->rmap->n,&m);
2007:   PetscBLASIntCast(B->rmap->n,&n);
2008:   PetscBLASIntCast(A->cmap->n,&k);
2009:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
2010:   return(0);
2011: }

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

2018:   if (scall == MAT_INITIAL_MATRIX) {
2019:     PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);
2020:     MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
2021:     PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);
2022:   }
2023:   PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);
2024:   MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);
2025:   PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);
2026:   return(0);
2027: }

2029: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
2030: {
2032:   PetscInt       m=A->cmap->n,n=B->cmap->n;
2033:   Mat            Cmat;

2036:   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);
2037:   MatCreate(PETSC_COMM_SELF,&Cmat);
2038:   MatSetSizes(Cmat,m,n,m,n);
2039:   MatSetType(Cmat,MATSEQDENSE);
2040:   MatSeqDenseSetPreallocation(Cmat,NULL);

2042:   Cmat->assembled = PETSC_TRUE;

2044:   *C = Cmat;
2045:   return(0);
2046: }

2048: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2049: {
2050:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2051:   Mat_SeqDense   *b = (Mat_SeqDense*)B->data;
2052:   Mat_SeqDense   *c = (Mat_SeqDense*)C->data;
2053:   PetscBLASInt   m,n,k;
2054:   PetscScalar    _DOne=1.0,_DZero=0.0;

2058:   PetscBLASIntCast(C->rmap->n,&m);
2059:   PetscBLASIntCast(C->cmap->n,&n);
2060:   PetscBLASIntCast(A->rmap->n,&k);
2061:   PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
2062:   return(0);
2063: }

2065: static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[])
2066: {
2067:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2069:   PetscInt       i,j,m = A->rmap->n,n = A->cmap->n,p;
2070:   PetscScalar    *x;
2071:   MatScalar      *aa = a->v;

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

2076:   VecSet(v,0.0);
2077:   VecGetArray(v,&x);
2078:   VecGetLocalSize(v,&p);
2079:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2080:   for (i=0; i<m; i++) {
2081:     x[i] = aa[i]; if (idx) idx[i] = 0;
2082:     for (j=1; j<n; j++) {
2083:       if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
2084:     }
2085:   }
2086:   VecRestoreArray(v,&x);
2087:   return(0);
2088: }

2090: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])
2091: {
2092:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2094:   PetscInt       i,j,m = A->rmap->n,n = A->cmap->n,p;
2095:   PetscScalar    *x;
2096:   PetscReal      atmp;
2097:   MatScalar      *aa = a->v;

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

2102:   VecSet(v,0.0);
2103:   VecGetArray(v,&x);
2104:   VecGetLocalSize(v,&p);
2105:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2106:   for (i=0; i<m; i++) {
2107:     x[i] = PetscAbsScalar(aa[i]);
2108:     for (j=1; j<n; j++) {
2109:       atmp = PetscAbsScalar(aa[i+m*j]);
2110:       if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;}
2111:     }
2112:   }
2113:   VecRestoreArray(v,&x);
2114:   return(0);
2115: }

2117: static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])
2118: {
2119:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2121:   PetscInt       i,j,m = A->rmap->n,n = A->cmap->n,p;
2122:   PetscScalar    *x;
2123:   MatScalar      *aa = a->v;

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

2128:   VecSet(v,0.0);
2129:   VecGetArray(v,&x);
2130:   VecGetLocalSize(v,&p);
2131:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2132:   for (i=0; i<m; i++) {
2133:     x[i] = aa[i]; if (idx) idx[i] = 0;
2134:     for (j=1; j<n; j++) {
2135:       if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
2136:     }
2137:   }
2138:   VecRestoreArray(v,&x);
2139:   return(0);
2140: }

2142: static PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
2143: {
2144:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2146:   PetscScalar    *x;

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

2151:   VecGetArray(v,&x);
2152:   PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));
2153:   VecRestoreArray(v,&x);
2154:   return(0);
2155: }


2158: PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms)
2159: {
2161:   PetscInt       i,j,m,n;
2162:   PetscScalar    *a;

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

2197: static PetscErrorCode  MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
2198: {
2200:   PetscScalar    *a;
2201:   PetscInt       m,n,i;

2204:   MatGetSize(x,&m,&n);
2205:   MatDenseGetArray(x,&a);
2206:   for (i=0; i<m*n; i++) {
2207:     PetscRandomGetValue(rctx,a+i);
2208:   }
2209:   MatDenseRestoreArray(x,&a);
2210:   return(0);
2211: }

2213: static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool  *missing,PetscInt *d)
2214: {
2216:   *missing = PETSC_FALSE;
2217:   return(0);
2218: }


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

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

2371:    Collective on MPI_Comm

2373:    Input Parameters:
2374: +  comm - MPI communicator, set to PETSC_COMM_SELF
2375: .  m - number of rows
2376: .  n - number of columns
2377: -  data - optional location of matrix data in column major order.  Set data=NULL for PETSc
2378:    to control all matrix memory allocation.

2380:    Output Parameter:
2381: .  A - the matrix

2383:    Notes:
2384:    The data input variable is intended primarily for Fortran programmers
2385:    who wish to allocate their own matrix memory space.  Most users should
2386:    set data=NULL.

2388:    Level: intermediate

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

2392: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
2393: @*/
2394: PetscErrorCode  MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
2395: {

2399:   MatCreate(comm,A);
2400:   MatSetSizes(*A,m,n,m,n);
2401:   MatSetType(*A,MATSEQDENSE);
2402:   MatSeqDenseSetPreallocation(*A,data);
2403:   return(0);
2404: }

2406: /*@C
2407:    MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements

2409:    Collective on MPI_Comm

2411:    Input Parameters:
2412: +  B - the matrix
2413: -  data - the array (or NULL)

2415:    Notes:
2416:    The data input variable is intended primarily for Fortran programmers
2417:    who wish to allocate their own matrix memory space.  Most users should
2418:    need not call this routine.

2420:    Level: intermediate

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

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

2426: @*/
2427: PetscErrorCode  MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
2428: {

2432:   PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
2433:   return(0);
2434: }

2436: PetscErrorCode  MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
2437: {
2438:   Mat_SeqDense   *b;

2442:   B->preallocated = PETSC_TRUE;

2444:   PetscLayoutSetUp(B->rmap);
2445:   PetscLayoutSetUp(B->cmap);

2447:   b       = (Mat_SeqDense*)B->data;
2448:   b->Mmax = B->rmap->n;
2449:   b->Nmax = B->cmap->n;
2450:   if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n;

2452:   PetscIntMultError(b->lda,b->Nmax,NULL);
2453:   if (!data) { /* petsc-allocated storage */
2454:     if (!b->user_alloc) { PetscFree(b->v); }
2455:     PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);
2456:     PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));

2458:     b->user_alloc = PETSC_FALSE;
2459:   } else { /* user-allocated storage */
2460:     if (!b->user_alloc) { PetscFree(b->v); }
2461:     b->v          = data;
2462:     b->user_alloc = PETSC_TRUE;
2463:   }
2464:   B->assembled = PETSC_TRUE;
2465:   return(0);
2466: }

2468: #if defined(PETSC_HAVE_ELEMENTAL)
2469: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
2470: {
2471:   Mat            mat_elemental;
2473:   PetscScalar    *array,*v_colwise;
2474:   PetscInt       M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols;

2477:   PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);
2478:   MatDenseGetArray(A,&array);
2479:   /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
2480:   k = 0;
2481:   for (j=0; j<N; j++) {
2482:     cols[j] = j;
2483:     for (i=0; i<M; i++) {
2484:       v_colwise[j*M+i] = array[k++];
2485:     }
2486:   }
2487:   for (i=0; i<M; i++) {
2488:     rows[i] = i;
2489:   }
2490:   MatDenseRestoreArray(A,&array);

2492:   MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
2493:   MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);
2494:   MatSetType(mat_elemental,MATELEMENTAL);
2495:   MatSetUp(mat_elemental);

2497:   /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
2498:   MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);
2499:   MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
2500:   MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);
2501:   PetscFree3(v_colwise,rows,cols);

2503:   if (reuse == MAT_INPLACE_MATRIX) {
2504:     MatHeaderReplace(A,&mat_elemental);
2505:   } else {
2506:     *newmat = mat_elemental;
2507:   }
2508:   return(0);
2509: }
2510: #endif

2512: /*@C
2513:   MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array

2515:   Input parameter:
2516: + A - the matrix
2517: - lda - the leading dimension

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

2524:   Level: intermediate

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

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

2530: @*/
2531: PetscErrorCode  MatSeqDenseSetLDA(Mat B,PetscInt lda)
2532: {
2533:   Mat_SeqDense *b = (Mat_SeqDense*)B->data;

2536:   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);
2537:   b->lda       = lda;
2538:   b->changelda = PETSC_FALSE;
2539:   b->Mmax      = PetscMax(b->Mmax,lda);
2540:   return(0);
2541: }

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

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

2549:   Level: beginner

2551: .seealso: MatCreateSeqDense()

2553: M*/

2555: PETSC_EXTERN PetscErrorCode MatCreate_SeqDense(Mat B)
2556: {
2557:   Mat_SeqDense   *b;
2559:   PetscMPIInt    size;

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

2565:   PetscNewLog(B,&b);
2566:   PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
2567:   B->data = (void*)b;

2569:   b->roworiented = PETSC_TRUE;

2571:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);
2572:   PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);
2573:   PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);
2574:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);
2575:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);
2576: #if defined(PETSC_HAVE_ELEMENTAL)
2577:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);
2578: #endif
2579:   PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);
2580:   PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2581:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2582:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2583:   PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaij_seqdense_C",MatPtAP_SeqDense_SeqDense);
2584:   PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijperm_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2585:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijperm_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2586:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijperm_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2587:   PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijperm_seqdense_C",MatPtAP_SeqDense_SeqDense);
2588:   PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijmkl_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2589:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijmkl_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2590:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijmkl_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2591:   PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijmkl_seqdense_C",MatPtAP_SeqDense_SeqDense);

2593:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2594:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2595:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2596:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijperm_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2597:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijperm_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2598:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijperm_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);

2600:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijmkl_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2601:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijmkl_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2602:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijmkl_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2603:   PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
2604:   return(0);
2605: }