Actual source code: dense.c

petsc-dev 2014-04-19
Report Typos and Errors
  2: /*
  3:      Defines the basic matrix operations for sequential dense.
  4: */

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

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

 13: PETSC_EXTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
 14: {
 15:   Mat            B;
 16:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
 18:   PetscInt       i, j;
 19:   PetscInt       *rows, *nnz;
 20:   MatScalar      *aa = a->v, *vals;

 23:   MatCreate(PetscObjectComm((PetscObject)A),&B);
 24:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
 25:   MatSetType(B,MATSEQAIJ);
 26:   PetscCalloc3(A->rmap->n,&rows,A->rmap->n,&nnz,A->rmap->n,&vals);
 27:   for (j=0; j<A->cmap->n; j++) {
 28:     for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) ++nnz[i];
 29:     aa += a->lda;
 30:   }
 31:   MatSeqAIJSetPreallocation(B,PETSC_DETERMINE,nnz);
 32:   aa = a->v;
 33:   for (j=0; j<A->cmap->n; j++) {
 34:     PetscInt numRows = 0;
 35:     for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) {rows[numRows] = i; vals[numRows++] = aa[i];}
 36:     MatSetValues(B,numRows,rows,1,&j,vals,INSERT_VALUES);
 37:     aa  += a->lda;
 38:   }
 39:   PetscFree3(rows,nnz,vals);
 40:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 41:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

 43:   if (reuse == MAT_REUSE_MATRIX) {
 44:     MatHeaderReplace(A,B);
 45:   } else {
 46:     *newmat = B;
 47:   }
 48:   return(0);
 49: }

 53: PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str)
 54: {
 55:   Mat_SeqDense   *x     = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data;
 56:   PetscScalar    oalpha = alpha;
 57:   PetscInt       j;
 58:   PetscBLASInt   N,m,ldax,lday,one = 1;

 62:   PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);
 63:   PetscBLASIntCast(X->rmap->n,&m);
 64:   PetscBLASIntCast(x->lda,&ldax);
 65:   PetscBLASIntCast(y->lda,&lday);
 66:   if (ldax>m || lday>m) {
 67:     for (j=0; j<X->cmap->n; j++) {
 68:       PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&oalpha,x->v+j*ldax,&one,y->v+j*lday,&one));
 69:     }
 70:   } else {
 71:     PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&oalpha,x->v,&one,y->v,&one));
 72:   }
 73:   PetscLogFlops(PetscMax(2*N-1,0));
 74:   return(0);
 75: }

 79: PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info)
 80: {
 81:   PetscInt N = A->rmap->n*A->cmap->n;

 84:   info->block_size        = 1.0;
 85:   info->nz_allocated      = (double)N;
 86:   info->nz_used           = (double)N;
 87:   info->nz_unneeded       = (double)0;
 88:   info->assemblies        = (double)A->num_ass;
 89:   info->mallocs           = 0;
 90:   info->memory            = ((PetscObject)A)->mem;
 91:   info->fill_ratio_given  = 0;
 92:   info->fill_ratio_needed = 0;
 93:   info->factor_mallocs    = 0;
 94:   return(0);
 95: }

 99: PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha)
100: {
101:   Mat_SeqDense   *a     = (Mat_SeqDense*)A->data;
102:   PetscScalar    oalpha = alpha;
104:   PetscBLASInt   one = 1,j,nz,lda;

107:   PetscBLASIntCast(a->lda,&lda);
108:   if (lda>A->rmap->n) {
109:     PetscBLASIntCast(A->rmap->n,&nz);
110:     for (j=0; j<A->cmap->n; j++) {
111:       PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v+j*lda,&one));
112:     }
113:   } else {
114:     PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);
115:     PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v,&one));
116:   }
117:   PetscLogFlops(nz);
118:   return(0);
119: }

123: PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool  *fl)
124: {
125:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
126:   PetscInt     i,j,m = A->rmap->n,N;
127:   PetscScalar  *v = a->v;

130:   *fl = PETSC_FALSE;
131:   if (A->rmap->n != A->cmap->n) return(0);
132:   N = a->lda;

134:   for (i=0; i<m; i++) {
135:     for (j=i+1; j<m; j++) {
136:       if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) return(0);
137:     }
138:   }
139:   *fl = PETSC_TRUE;
140:   return(0);
141: }

145: PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues)
146: {
147:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data,*l;
149:   PetscInt       lda = (PetscInt)mat->lda,j,m;

152:   PetscLayoutReference(A->rmap,&newi->rmap);
153:   PetscLayoutReference(A->cmap,&newi->cmap);
154:   MatSeqDenseSetPreallocation(newi,NULL);
155:   if (cpvalues == MAT_COPY_VALUES) {
156:     l = (Mat_SeqDense*)newi->data;
157:     if (lda>A->rmap->n) {
158:       m = A->rmap->n;
159:       for (j=0; j<A->cmap->n; j++) {
160:         PetscMemcpy(l->v+j*m,mat->v+j*lda,m*sizeof(PetscScalar));
161:       }
162:     } else {
163:       PetscMemcpy(l->v,mat->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
164:     }
165:   }
166:   newi->assembled = PETSC_TRUE;
167:   return(0);
168: }

172: PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat)
173: {

177:   MatCreate(PetscObjectComm((PetscObject)A),newmat);
178:   MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
179:   MatSetType(*newmat,((PetscObject)A)->type_name);
180:   MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);
181:   return(0);
182: }


185: extern PetscErrorCode MatLUFactor_SeqDense(Mat,IS,IS,const MatFactorInfo*);

189: PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
190: {
191:   MatFactorInfo  info;

195:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
196:   MatLUFactor_SeqDense(fact,0,0,&info);
197:   return(0);
198: }

202: PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy)
203: {
204:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
206:   PetscScalar    *x,*y;
207:   PetscBLASInt   one = 1,info,m;

210:   PetscBLASIntCast(A->rmap->n,&m);
211:   VecGetArray(xx,&x);
212:   VecGetArray(yy,&y);
213:   PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
214:   if (A->factortype == MAT_FACTOR_LU) {
215: #if defined(PETSC_MISSING_LAPACK_GETRS)
216:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
217: #else
218:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
219:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
220: #endif
221:   } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
222: #if defined(PETSC_MISSING_LAPACK_POTRS)
223:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
224: #else
225:     PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
226:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
227: #endif
228:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
229:   VecRestoreArray(xx,&x);
230:   VecRestoreArray(yy,&y);
231:   PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
232:   return(0);
233: }

237: PetscErrorCode MatMatSolve_SeqDense(Mat A,Mat B,Mat X)
238: {
239:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
241:   PetscScalar    *b,*x;
242:   PetscInt       n;
243:   PetscBLASInt   nrhs,info,m;
244:   PetscBool      flg;

247:   PetscBLASIntCast(A->rmap->n,&m);
248:   PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
249:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
250:   PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
251:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");

253:   MatGetSize(B,NULL,&n);
254:   PetscBLASIntCast(n,&nrhs);
255:   MatDenseGetArray(B,&b);
256:   MatDenseGetArray(X,&x);

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

260:   if (A->factortype == MAT_FACTOR_LU) {
261: #if defined(PETSC_MISSING_LAPACK_GETRS)
262:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
263: #else
264:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
265:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
266: #endif
267:   } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
268: #if defined(PETSC_MISSING_LAPACK_POTRS)
269:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
270: #else
271:     PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info));
272:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
273: #endif
274:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");

276:   MatDenseRestoreArray(B,&b);
277:   MatDenseRestoreArray(X,&x);
278:   PetscLogFlops(nrhs*(2.0*m*m - m));
279:   return(0);
280: }

284: PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy)
285: {
286:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
288:   PetscScalar    *x,*y;
289:   PetscBLASInt   one = 1,info,m;

292:   PetscBLASIntCast(A->rmap->n,&m);
293:   VecGetArray(xx,&x);
294:   VecGetArray(yy,&y);
295:   PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
296:   /* assume if pivots exist then use LU; else Cholesky */
297:   if (mat->pivots) {
298: #if defined(PETSC_MISSING_LAPACK_GETRS)
299:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
300: #else
301:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
302:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve");
303: #endif
304:   } else {
305: #if defined(PETSC_MISSING_LAPACK_POTRS)
306:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
307: #else
308:     PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
309:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve");
310: #endif
311:   }
312:   VecRestoreArray(xx,&x);
313:   VecRestoreArray(yy,&y);
314:   PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
315:   return(0);
316: }

320: PetscErrorCode MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
321: {
322:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
324:   PetscScalar    *x,*y,sone = 1.0;
325:   Vec            tmp = 0;
326:   PetscBLASInt   one = 1,info,m;

329:   PetscBLASIntCast(A->rmap->n,&m);
330:   VecGetArray(xx,&x);
331:   VecGetArray(yy,&y);
332:   if (!A->rmap->n || !A->cmap->n) return(0);
333:   if (yy == zz) {
334:     VecDuplicate(yy,&tmp);
335:     PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);
336:     VecCopy(yy,tmp);
337:   }
338:   PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
339:   /* assume if pivots exist then use LU; else Cholesky */
340:   if (mat->pivots) {
341: #if defined(PETSC_MISSING_LAPACK_GETRS)
342:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
343: #else
344:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
345:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
346: #endif
347:   } else {
348: #if defined(PETSC_MISSING_LAPACK_POTRS)
349:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
350: #else
351:     PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
352:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
353: #endif
354:   }
355:   if (tmp) {
356:     VecAXPY(yy,sone,tmp);
357:     VecDestroy(&tmp);
358:   } else {
359:     VecAXPY(yy,sone,zz);
360:   }
361:   VecRestoreArray(xx,&x);
362:   VecRestoreArray(yy,&y);
363:   PetscLogFlops(2.0*A->cmap->n*A->cmap->n);
364:   return(0);
365: }

369: PetscErrorCode MatSolveTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
370: {
371:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
373:   PetscScalar    *x,*y,sone = 1.0;
374:   Vec            tmp;
375:   PetscBLASInt   one = 1,info,m;

378:   PetscBLASIntCast(A->rmap->n,&m);
379:   if (!A->rmap->n || !A->cmap->n) return(0);
380:   VecGetArray(xx,&x);
381:   VecGetArray(yy,&y);
382:   if (yy == zz) {
383:     VecDuplicate(yy,&tmp);
384:     PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);
385:     VecCopy(yy,tmp);
386:   }
387:   PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
388:   /* assume if pivots exist then use LU; else Cholesky */
389:   if (mat->pivots) {
390: #if defined(PETSC_MISSING_LAPACK_GETRS)
391:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
392: #else
393:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
394:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
395: #endif
396:   } else {
397: #if defined(PETSC_MISSING_LAPACK_POTRS)
398:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
399: #else
400:     PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
401:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
402: #endif
403:   }
404:   if (tmp) {
405:     VecAXPY(yy,sone,tmp);
406:     VecDestroy(&tmp);
407:   } else {
408:     VecAXPY(yy,sone,zz);
409:   }
410:   VecRestoreArray(xx,&x);
411:   VecRestoreArray(yy,&y);
412:   PetscLogFlops(2.0*A->cmap->n*A->cmap->n);
413:   return(0);
414: }

416: /* ---------------------------------------------------------------*/
417: /* COMMENT: I have chosen to hide row permutation in the pivots,
418:    rather than put it in the Mat->row slot.*/
421: PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo)
422: {
423: #if defined(PETSC_MISSING_LAPACK_GETRF)
425:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable.");
426: #else
427:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
429:   PetscBLASInt   n,m,info;

432:   PetscBLASIntCast(A->cmap->n,&n);
433:   PetscBLASIntCast(A->rmap->n,&m);
434:   if (!mat->pivots) {
435:     PetscMalloc1((A->rmap->n+1),&mat->pivots);
436:     PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
437:   }
438:   if (!A->rmap->n || !A->cmap->n) return(0);
439:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
440:   PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info));
441:   PetscFPTrapPop();

443:   if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization");
444:   if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization");
445:   A->ops->solve             = MatSolve_SeqDense;
446:   A->ops->solvetranspose    = MatSolveTranspose_SeqDense;
447:   A->ops->solveadd          = MatSolveAdd_SeqDense;
448:   A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
449:   A->factortype             = MAT_FACTOR_LU;

451:   PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);
452: #endif
453:   return(0);
454: }

458: PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo)
459: {
460: #if defined(PETSC_MISSING_LAPACK_POTRF)
462:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable.");
463: #else
464:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
466:   PetscBLASInt   info,n;

469:   PetscBLASIntCast(A->cmap->n,&n);
470:   PetscFree(mat->pivots);

472:   if (!A->rmap->n || !A->cmap->n) return(0);
473:   PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info));
474:   if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1);
475:   A->ops->solve             = MatSolve_SeqDense;
476:   A->ops->solvetranspose    = MatSolveTranspose_SeqDense;
477:   A->ops->solveadd          = MatSolveAdd_SeqDense;
478:   A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
479:   A->factortype             = MAT_FACTOR_CHOLESKY;

481:   PetscLogFlops((A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
482: #endif
483:   return(0);
484: }


489: PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
490: {
492:   MatFactorInfo  info;

495:   info.fill = 1.0;

497:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
498:   MatCholeskyFactor_SeqDense(fact,0,&info);
499:   return(0);
500: }

504: PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
505: {
507:   fact->assembled                  = PETSC_TRUE;
508:   fact->preallocated               = PETSC_TRUE;
509:   fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
510:   return(0);
511: }

515: PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info)
516: {
518:   fact->preallocated         = PETSC_TRUE;
519:   fact->assembled            = PETSC_TRUE;
520:   fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense;
521:   return(0);
522: }

526: PETSC_EXTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact)
527: {

531:   MatCreate(PetscObjectComm((PetscObject)A),fact);
532:   MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
533:   MatSetType(*fact,((PetscObject)A)->type_name);
534:   if (ftype == MAT_FACTOR_LU) {
535:     (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
536:   } else {
537:     (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
538:   }
539:   (*fact)->factortype = ftype;
540:   return(0);
541: }

543: /* ------------------------------------------------------------------*/
546: PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx)
547: {
548:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
549:   PetscScalar    *x,*b,*v = mat->v,zero = 0.0,xt;
551:   PetscInt       m = A->rmap->n,i;
552:   PetscBLASInt   o = 1,bm;

555:   PetscBLASIntCast(m,&bm);
556:   if (flag & SOR_ZERO_INITIAL_GUESS) {
557:     /* this is a hack fix, should have another version without the second BLASdot */
558:     VecSet(xx,zero);
559:   }
560:   VecGetArray(xx,&x);
561:   VecGetArray(bb,&b);
562:   its  = its*lits;
563:   if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
564:   while (its--) {
565:     if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
566:       for (i=0; i<m; i++) {
567:         PetscStackCallBLAS("BLASdot",xt   = b[i] - BLASdot_(&bm,v+i,&bm,x,&o));
568:         x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
569:       }
570:     }
571:     if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
572:       for (i=m-1; i>=0; i--) {
573:         PetscStackCallBLAS("BLASdot",xt   = b[i] - BLASdot_(&bm,v+i,&bm,x,&o));
574:         x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
575:       }
576:     }
577:   }
578:   VecRestoreArray(bb,&b);
579:   VecRestoreArray(xx,&x);
580:   return(0);
581: }

583: /* -----------------------------------------------------------------*/
586: PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)
587: {
588:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
589:   PetscScalar    *v   = mat->v,*x,*y;
591:   PetscBLASInt   m, n,_One=1;
592:   PetscScalar    _DOne=1.0,_DZero=0.0;

595:   PetscBLASIntCast(A->rmap->n,&m);
596:   PetscBLASIntCast(A->cmap->n,&n);
597:   if (!A->rmap->n || !A->cmap->n) return(0);
598:   VecGetArray(xx,&x);
599:   VecGetArray(yy,&y);
600:   PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One));
601:   VecRestoreArray(xx,&x);
602:   VecRestoreArray(yy,&y);
603:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);
604:   return(0);
605: }

609: PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy)
610: {
611:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
612:   PetscScalar    *v   = mat->v,*x,*y,_DOne=1.0,_DZero=0.0;
614:   PetscBLASInt   m, n, _One=1;

617:   PetscBLASIntCast(A->rmap->n,&m);
618:   PetscBLASIntCast(A->cmap->n,&n);
619:   if (!A->rmap->n || !A->cmap->n) return(0);
620:   VecGetArray(xx,&x);
621:   VecGetArray(yy,&y);
622:   PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One));
623:   VecRestoreArray(xx,&x);
624:   VecRestoreArray(yy,&y);
625:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);
626:   return(0);
627: }

631: PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
632: {
633:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
634:   PetscScalar    *v   = mat->v,*x,*y,_DOne=1.0;
636:   PetscBLASInt   m, n, _One=1;

639:   PetscBLASIntCast(A->rmap->n,&m);
640:   PetscBLASIntCast(A->cmap->n,&n);
641:   if (!A->rmap->n || !A->cmap->n) return(0);
642:   if (zz != yy) {VecCopy(zz,yy);}
643:   VecGetArray(xx,&x);
644:   VecGetArray(yy,&y);
645:   PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
646:   VecRestoreArray(xx,&x);
647:   VecRestoreArray(yy,&y);
648:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
649:   return(0);
650: }

654: PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
655: {
656:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
657:   PetscScalar    *v   = mat->v,*x,*y;
659:   PetscBLASInt   m, n, _One=1;
660:   PetscScalar    _DOne=1.0;

663:   PetscBLASIntCast(A->rmap->n,&m);
664:   PetscBLASIntCast(A->cmap->n,&n);
665:   if (!A->rmap->n || !A->cmap->n) return(0);
666:   if (zz != yy) {VecCopy(zz,yy);}
667:   VecGetArray(xx,&x);
668:   VecGetArray(yy,&y);
669:   PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
670:   VecRestoreArray(xx,&x);
671:   VecRestoreArray(yy,&y);
672:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
673:   return(0);
674: }

676: /* -----------------------------------------------------------------*/
679: PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
680: {
681:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
682:   PetscScalar    *v;
684:   PetscInt       i;

687:   *ncols = A->cmap->n;
688:   if (cols) {
689:     PetscMalloc1((A->cmap->n+1),cols);
690:     for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
691:   }
692:   if (vals) {
693:     PetscMalloc1((A->cmap->n+1),vals);
694:     v    = mat->v + row;
695:     for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
696:   }
697:   return(0);
698: }

702: PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
703: {

707:   if (cols) {PetscFree(*cols);}
708:   if (vals) {PetscFree(*vals); }
709:   return(0);
710: }
711: /* ----------------------------------------------------------------*/
714: PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
715: {
716:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
717:   PetscInt     i,j,idx=0;

720:   if (!mat->roworiented) {
721:     if (addv == INSERT_VALUES) {
722:       for (j=0; j<n; j++) {
723:         if (indexn[j] < 0) {idx += m; continue;}
724: #if defined(PETSC_USE_DEBUG)
725:         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);
726: #endif
727:         for (i=0; i<m; i++) {
728:           if (indexm[i] < 0) {idx++; continue;}
729: #if defined(PETSC_USE_DEBUG)
730:           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);
731: #endif
732:           mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
733:         }
734:       }
735:     } else {
736:       for (j=0; j<n; j++) {
737:         if (indexn[j] < 0) {idx += m; continue;}
738: #if defined(PETSC_USE_DEBUG)
739:         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);
740: #endif
741:         for (i=0; i<m; i++) {
742:           if (indexm[i] < 0) {idx++; continue;}
743: #if defined(PETSC_USE_DEBUG)
744:           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);
745: #endif
746:           mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
747:         }
748:       }
749:     }
750:   } else {
751:     if (addv == INSERT_VALUES) {
752:       for (i=0; i<m; i++) {
753:         if (indexm[i] < 0) { idx += n; continue;}
754: #if defined(PETSC_USE_DEBUG)
755:         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);
756: #endif
757:         for (j=0; j<n; j++) {
758:           if (indexn[j] < 0) { idx++; continue;}
759: #if defined(PETSC_USE_DEBUG)
760:           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);
761: #endif
762:           mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
763:         }
764:       }
765:     } else {
766:       for (i=0; i<m; i++) {
767:         if (indexm[i] < 0) { idx += n; continue;}
768: #if defined(PETSC_USE_DEBUG)
769:         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);
770: #endif
771:         for (j=0; j<n; j++) {
772:           if (indexn[j] < 0) { idx++; continue;}
773: #if defined(PETSC_USE_DEBUG)
774:           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);
775: #endif
776:           mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
777:         }
778:       }
779:     }
780:   }
781:   return(0);
782: }

786: PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
787: {
788:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
789:   PetscInt     i,j;

792:   /* row-oriented output */
793:   for (i=0; i<m; i++) {
794:     if (indexm[i] < 0) {v += n;continue;}
795:     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);
796:     for (j=0; j<n; j++) {
797:       if (indexn[j] < 0) {v++; continue;}
798:       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);
799:       *v++ = mat->v[indexn[j]*mat->lda + indexm[i]];
800:     }
801:   }
802:   return(0);
803: }

805: /* -----------------------------------------------------------------*/

809: PetscErrorCode MatLoad_SeqDense(Mat newmat,PetscViewer viewer)
810: {
811:   Mat_SeqDense   *a;
813:   PetscInt       *scols,i,j,nz,header[4];
814:   int            fd;
815:   PetscMPIInt    size;
816:   PetscInt       *rowlengths = 0,M,N,*cols,grows,gcols;
817:   PetscScalar    *vals,*svals,*v,*w;
818:   MPI_Comm       comm;

821:   PetscObjectGetComm((PetscObject)viewer,&comm);
822:   MPI_Comm_size(comm,&size);
823:   if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor");
824:   PetscViewerBinaryGetDescriptor(viewer,&fd);
825:   PetscBinaryRead(fd,header,4,PETSC_INT);
826:   if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not matrix object");
827:   M = header[1]; N = header[2]; nz = header[3];

829:   /* set global size if not set already*/
830:   if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) {
831:     MatSetSizes(newmat,M,N,M,N);
832:   } else {
833:     /* if sizes and type are already set, check if the vector global sizes are correct */
834:     MatGetSize(newmat,&grows,&gcols);
835:     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);
836:   }
837:   MatSeqDenseSetPreallocation(newmat,NULL);

839:   if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */
840:     a = (Mat_SeqDense*)newmat->data;
841:     v = a->v;
842:     /* Allocate some temp space to read in the values and then flip them
843:        from row major to column major */
844:     PetscMalloc1((M*N > 0 ? M*N : 1),&w);
845:     /* read in nonzero values */
846:     PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);
847:     /* now flip the values and store them in the matrix*/
848:     for (j=0; j<N; j++) {
849:       for (i=0; i<M; i++) {
850:         *v++ =w[i*N+j];
851:       }
852:     }
853:     PetscFree(w);
854:   } else {
855:     /* read row lengths */
856:     PetscMalloc1((M+1),&rowlengths);
857:     PetscBinaryRead(fd,rowlengths,M,PETSC_INT);

859:     a = (Mat_SeqDense*)newmat->data;
860:     v = a->v;

862:     /* read column indices and nonzeros */
863:     PetscMalloc1((nz+1),&scols);
864:     cols = scols;
865:     PetscBinaryRead(fd,cols,nz,PETSC_INT);
866:     PetscMalloc1((nz+1),&svals);
867:     vals = svals;
868:     PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);

870:     /* insert into matrix */
871:     for (i=0; i<M; i++) {
872:       for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j];
873:       svals += rowlengths[i]; scols += rowlengths[i];
874:     }
875:     PetscFree(vals);
876:     PetscFree(cols);
877:     PetscFree(rowlengths);
878:   }
879:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
880:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
881:   return(0);
882: }

886: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
887: {
888:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
889:   PetscErrorCode    ierr;
890:   PetscInt          i,j;
891:   const char        *name;
892:   PetscScalar       *v;
893:   PetscViewerFormat format;
894: #if defined(PETSC_USE_COMPLEX)
895:   PetscBool         allreal = PETSC_TRUE;
896: #endif

899:   PetscViewerGetFormat(viewer,&format);
900:   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
901:     return(0);  /* do nothing for now */
902:   } else if (format == PETSC_VIEWER_ASCII_COMMON) {
903:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
904:     for (i=0; i<A->rmap->n; i++) {
905:       v    = a->v + i;
906:       PetscViewerASCIIPrintf(viewer,"row %D:",i);
907:       for (j=0; j<A->cmap->n; j++) {
908: #if defined(PETSC_USE_COMPLEX)
909:         if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
910:           PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
911:         } else if (PetscRealPart(*v)) {
912:           PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));
913:         }
914: #else
915:         if (*v) {
916:           PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);
917:         }
918: #endif
919:         v += a->lda;
920:       }
921:       PetscViewerASCIIPrintf(viewer,"\n");
922:     }
923:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
924:   } else {
925:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
926: #if defined(PETSC_USE_COMPLEX)
927:     /* determine if matrix has all real values */
928:     v = a->v;
929:     for (i=0; i<A->rmap->n*A->cmap->n; i++) {
930:       if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;}
931:     }
932: #endif
933:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
934:       PetscObjectGetName((PetscObject)A,&name);
935:       PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);
936:       PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);
937:       PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
938:     }

940:     for (i=0; i<A->rmap->n; i++) {
941:       v = a->v + i;
942:       for (j=0; j<A->cmap->n; j++) {
943: #if defined(PETSC_USE_COMPLEX)
944:         if (allreal) {
945:           PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));
946:         } else {
947:           PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
948:         }
949: #else
950:         PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);
951: #endif
952:         v += a->lda;
953:       }
954:       PetscViewerASCIIPrintf(viewer,"\n");
955:     }
956:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
957:       PetscViewerASCIIPrintf(viewer,"];\n");
958:     }
959:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
960:   }
961:   PetscViewerFlush(viewer);
962:   return(0);
963: }

967: static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer)
968: {
969:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
970:   PetscErrorCode    ierr;
971:   int               fd;
972:   PetscInt          ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n;
973:   PetscScalar       *v,*anonz,*vals;
974:   PetscViewerFormat format;

977:   PetscViewerBinaryGetDescriptor(viewer,&fd);

979:   PetscViewerGetFormat(viewer,&format);
980:   if (format == PETSC_VIEWER_NATIVE) {
981:     /* store the matrix as a dense matrix */
982:     PetscMalloc1(4,&col_lens);

984:     col_lens[0] = MAT_FILE_CLASSID;
985:     col_lens[1] = m;
986:     col_lens[2] = n;
987:     col_lens[3] = MATRIX_BINARY_FORMAT_DENSE;

989:     PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);
990:     PetscFree(col_lens);

992:     /* write out matrix, by rows */
993:     PetscMalloc1((m*n+1),&vals);
994:     v    = a->v;
995:     for (j=0; j<n; j++) {
996:       for (i=0; i<m; i++) {
997:         vals[j + i*n] = *v++;
998:       }
999:     }
1000:     PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);
1001:     PetscFree(vals);
1002:   } else {
1003:     PetscMalloc1((4+nz),&col_lens);

1005:     col_lens[0] = MAT_FILE_CLASSID;
1006:     col_lens[1] = m;
1007:     col_lens[2] = n;
1008:     col_lens[3] = nz;

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

1014:     /* Possibly should write in smaller increments, not whole matrix at once? */
1015:     /* store column indices (zero start index) */
1016:     ict = 0;
1017:     for (i=0; i<m; i++) {
1018:       for (j=0; j<n; j++) col_lens[ict++] = j;
1019:     }
1020:     PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);
1021:     PetscFree(col_lens);

1023:     /* store nonzero values */
1024:     PetscMalloc1((nz+1),&anonz);
1025:     ict  = 0;
1026:     for (i=0; i<m; i++) {
1027:       v = a->v + i;
1028:       for (j=0; j<n; j++) {
1029:         anonz[ict++] = *v; v += a->lda;
1030:       }
1031:     }
1032:     PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);
1033:     PetscFree(anonz);
1034:   }
1035:   return(0);
1036: }

1038: #include <petscdraw.h>
1041: PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa)
1042: {
1043:   Mat               A  = (Mat) Aa;
1044:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
1045:   PetscErrorCode    ierr;
1046:   PetscInt          m  = A->rmap->n,n = A->cmap->n,color,i,j;
1047:   PetscScalar       *v = a->v;
1048:   PetscViewer       viewer;
1049:   PetscDraw         popup;
1050:   PetscReal         xl,yl,xr,yr,x_l,x_r,y_l,y_r,scale,maxv = 0.0;
1051:   PetscViewerFormat format;

1054:   PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
1055:   PetscViewerGetFormat(viewer,&format);
1056:   PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);

1058:   /* Loop over matrix elements drawing boxes */
1059:   if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1060:     /* Blue for negative and Red for positive */
1061:     color = PETSC_DRAW_BLUE;
1062:     for (j = 0; j < n; j++) {
1063:       x_l = j;
1064:       x_r = x_l + 1.0;
1065:       for (i = 0; i < m; i++) {
1066:         y_l = m - i - 1.0;
1067:         y_r = y_l + 1.0;
1068:         if (PetscRealPart(v[j*m+i]) >  0.) {
1069:           color = PETSC_DRAW_RED;
1070:         } else if (PetscRealPart(v[j*m+i]) <  0.) {
1071:           color = PETSC_DRAW_BLUE;
1072:         } else {
1073:           continue;
1074:         }
1075:         PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1076:       }
1077:     }
1078:   } else {
1079:     /* use contour shading to indicate magnitude of values */
1080:     /* first determine max of all nonzero values */
1081:     for (i = 0; i < m*n; i++) {
1082:       if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1083:     }
1084:     scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv;
1085:     PetscDrawGetPopup(draw,&popup);
1086:     if (popup) {PetscDrawScalePopup(popup,0.0,maxv);}
1087:     for (j = 0; j < n; j++) {
1088:       x_l = j;
1089:       x_r = x_l + 1.0;
1090:       for (i = 0; i < m; i++) {
1091:         y_l   = m - i - 1.0;
1092:         y_r   = y_l + 1.0;
1093:         color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(v[j*m+i]));
1094:         PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1095:       }
1096:     }
1097:   }
1098:   return(0);
1099: }

1103: PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
1104: {
1105:   PetscDraw      draw;
1106:   PetscBool      isnull;
1107:   PetscReal      xr,yr,xl,yl,h,w;

1111:   PetscViewerDrawGetDraw(viewer,0,&draw);
1112:   PetscDrawIsNull(draw,&isnull);
1113:   if (isnull) return(0);

1115:   PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
1116:   xr   = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
1117:   xr  += w;    yr += h;  xl = -w;     yl = -h;
1118:   PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
1119:   PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
1120:   PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
1121:   return(0);
1122: }

1126: PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
1127: {
1129:   PetscBool      iascii,isbinary,isdraw;

1132:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1133:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1134:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);

1136:   if (iascii) {
1137:     MatView_SeqDense_ASCII(A,viewer);
1138:   } else if (isbinary) {
1139:     MatView_SeqDense_Binary(A,viewer);
1140:   } else if (isdraw) {
1141:     MatView_SeqDense_Draw(A,viewer);
1142:   }
1143:   return(0);
1144: }

1148: PetscErrorCode MatDestroy_SeqDense(Mat mat)
1149: {
1150:   Mat_SeqDense   *l = (Mat_SeqDense*)mat->data;

1154: #if defined(PETSC_USE_LOG)
1155:   PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n);
1156: #endif
1157:   PetscFree(l->pivots);
1158:   if (!l->user_alloc) {PetscFree(l->v);}
1159:   PetscFree(mat->data);

1161:   PetscObjectChangeTypeName((PetscObject)mat,0);
1162:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);
1163:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);
1164:   PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);
1165:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);
1166:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);
1167:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);
1168:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);
1169:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);
1170:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);
1171:   return(0);
1172: }

1176: PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
1177: {
1178:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1180:   PetscInt       k,j,m,n,M;
1181:   PetscScalar    *v,tmp;

1184:   v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n;
1185:   if (reuse == MAT_REUSE_MATRIX && *matout == A) { /* in place transpose */
1186:     if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can not transpose non-square matrix in place");
1187:     else {
1188:       for (j=0; j<m; j++) {
1189:         for (k=0; k<j; k++) {
1190:           tmp        = v[j + k*M];
1191:           v[j + k*M] = v[k + j*M];
1192:           v[k + j*M] = tmp;
1193:         }
1194:       }
1195:     }
1196:   } else { /* out-of-place transpose */
1197:     Mat          tmat;
1198:     Mat_SeqDense *tmatd;
1199:     PetscScalar  *v2;

1201:     if (reuse == MAT_INITIAL_MATRIX) {
1202:       MatCreate(PetscObjectComm((PetscObject)A),&tmat);
1203:       MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);
1204:       MatSetType(tmat,((PetscObject)A)->type_name);
1205:       MatSeqDenseSetPreallocation(tmat,NULL);
1206:     } else {
1207:       tmat = *matout;
1208:     }
1209:     tmatd = (Mat_SeqDense*)tmat->data;
1210:     v = mat->v; v2 = tmatd->v;
1211:     for (j=0; j<n; j++) {
1212:       for (k=0; k<m; k++) v2[j + k*n] = v[k + j*M];
1213:     }
1214:     MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
1215:     MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
1216:     *matout = tmat;
1217:   }
1218:   return(0);
1219: }

1223: PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool  *flg)
1224: {
1225:   Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data;
1226:   Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data;
1227:   PetscInt     i,j;
1228:   PetscScalar  *v1,*v2;

1231:   if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; return(0);}
1232:   if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; return(0);}
1233:   for (i=0; i<A1->rmap->n; i++) {
1234:     v1 = mat1->v+i; v2 = mat2->v+i;
1235:     for (j=0; j<A1->cmap->n; j++) {
1236:       if (*v1 != *v2) {*flg = PETSC_FALSE; return(0);}
1237:       v1 += mat1->lda; v2 += mat2->lda;
1238:     }
1239:   }
1240:   *flg = PETSC_TRUE;
1241:   return(0);
1242: }

1246: PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
1247: {
1248:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1250:   PetscInt       i,n,len;
1251:   PetscScalar    *x,zero = 0.0;

1254:   VecSet(v,zero);
1255:   VecGetSize(v,&n);
1256:   VecGetArray(v,&x);
1257:   len  = PetscMin(A->rmap->n,A->cmap->n);
1258:   if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
1259:   for (i=0; i<len; i++) {
1260:     x[i] = mat->v[i*mat->lda + i];
1261:   }
1262:   VecRestoreArray(v,&x);
1263:   return(0);
1264: }

1268: PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
1269: {
1270:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1271:   PetscScalar    *l,*r,x,*v;
1273:   PetscInt       i,j,m = A->rmap->n,n = A->cmap->n;

1276:   if (ll) {
1277:     VecGetSize(ll,&m);
1278:     VecGetArray(ll,&l);
1279:     if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size");
1280:     for (i=0; i<m; i++) {
1281:       x = l[i];
1282:       v = mat->v + i;
1283:       for (j=0; j<n; j++) { (*v) *= x; v+= m;}
1284:     }
1285:     VecRestoreArray(ll,&l);
1286:     PetscLogFlops(n*m);
1287:   }
1288:   if (rr) {
1289:     VecGetSize(rr,&n);
1290:     VecGetArray(rr,&r);
1291:     if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size");
1292:     for (i=0; i<n; i++) {
1293:       x = r[i];
1294:       v = mat->v + i*m;
1295:       for (j=0; j<m; j++) (*v++) *= x;
1296:     }
1297:     VecRestoreArray(rr,&r);
1298:     PetscLogFlops(n*m);
1299:   }
1300:   return(0);
1301: }

1305: PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1306: {
1307:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1308:   PetscScalar    *v   = mat->v;
1309:   PetscReal      sum  = 0.0;
1310:   PetscInt       lda  =mat->lda,m=A->rmap->n,i,j;

1314:   if (type == NORM_FROBENIUS) {
1315:     if (lda>m) {
1316:       for (j=0; j<A->cmap->n; j++) {
1317:         v = mat->v+j*lda;
1318:         for (i=0; i<m; i++) {
1319:           sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1320:         }
1321:       }
1322:     } else {
1323:       for (i=0; i<A->cmap->n*A->rmap->n; i++) {
1324:         sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1325:       }
1326:     }
1327:     *nrm = PetscSqrtReal(sum);
1328:     PetscLogFlops(2.0*A->cmap->n*A->rmap->n);
1329:   } else if (type == NORM_1) {
1330:     *nrm = 0.0;
1331:     for (j=0; j<A->cmap->n; j++) {
1332:       v   = mat->v + j*mat->lda;
1333:       sum = 0.0;
1334:       for (i=0; i<A->rmap->n; i++) {
1335:         sum += PetscAbsScalar(*v);  v++;
1336:       }
1337:       if (sum > *nrm) *nrm = sum;
1338:     }
1339:     PetscLogFlops(A->cmap->n*A->rmap->n);
1340:   } else if (type == NORM_INFINITY) {
1341:     *nrm = 0.0;
1342:     for (j=0; j<A->rmap->n; j++) {
1343:       v   = mat->v + j;
1344:       sum = 0.0;
1345:       for (i=0; i<A->cmap->n; i++) {
1346:         sum += PetscAbsScalar(*v); v += mat->lda;
1347:       }
1348:       if (sum > *nrm) *nrm = sum;
1349:     }
1350:     PetscLogFlops(A->cmap->n*A->rmap->n);
1351:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
1352:   return(0);
1353: }

1357: PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
1358: {
1359:   Mat_SeqDense   *aij = (Mat_SeqDense*)A->data;

1363:   switch (op) {
1364:   case MAT_ROW_ORIENTED:
1365:     aij->roworiented = flg;
1366:     break;
1367:   case MAT_NEW_NONZERO_LOCATIONS:
1368:   case MAT_NEW_NONZERO_LOCATION_ERR:
1369:   case MAT_NEW_NONZERO_ALLOCATION_ERR:
1370:   case MAT_NEW_DIAGONALS:
1371:   case MAT_KEEP_NONZERO_PATTERN:
1372:   case MAT_IGNORE_OFF_PROC_ENTRIES:
1373:   case MAT_USE_HASH_TABLE:
1374:   case MAT_IGNORE_LOWER_TRIANGULAR:
1375:     PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
1376:     break;
1377:   case MAT_SPD:
1378:   case MAT_SYMMETRIC:
1379:   case MAT_STRUCTURALLY_SYMMETRIC:
1380:   case MAT_HERMITIAN:
1381:   case MAT_SYMMETRY_ETERNAL:
1382:     /* These options are handled directly by MatSetOption() */
1383:     break;
1384:   default:
1385:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]);
1386:   }
1387:   return(0);
1388: }

1392: PetscErrorCode MatZeroEntries_SeqDense(Mat A)
1393: {
1394:   Mat_SeqDense   *l = (Mat_SeqDense*)A->data;
1396:   PetscInt       lda=l->lda,m=A->rmap->n,j;

1399:   if (lda>m) {
1400:     for (j=0; j<A->cmap->n; j++) {
1401:       PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));
1402:     }
1403:   } else {
1404:     PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1405:   }
1406:   return(0);
1407: }

1411: PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1412: {
1413:   PetscErrorCode    ierr;
1414:   Mat_SeqDense      *l = (Mat_SeqDense*)A->data;
1415:   PetscInt          m  = l->lda, n = A->cmap->n, i,j;
1416:   PetscScalar       *slot,*bb;
1417:   const PetscScalar *xx;

1420: #if defined(PETSC_USE_DEBUG)
1421:   for (i=0; i<N; i++) {
1422:     if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
1423:     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);
1424:   }
1425: #endif

1427:   /* fix right hand side if needed */
1428:   if (x && b) {
1429:     VecGetArrayRead(x,&xx);
1430:     VecGetArray(b,&bb);
1431:     for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
1432:     VecRestoreArrayRead(x,&xx);
1433:     VecRestoreArray(b,&bb);
1434:   }

1436:   for (i=0; i<N; i++) {
1437:     slot = l->v + rows[i];
1438:     for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
1439:   }
1440:   if (diag != 0.0) {
1441:     if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
1442:     for (i=0; i<N; i++) {
1443:       slot  = l->v + (m+1)*rows[i];
1444:       *slot = diag;
1445:     }
1446:   }
1447:   return(0);
1448: }

1452: PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[])
1453: {
1454:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;

1457:   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");
1458:   *array = mat->v;
1459:   return(0);
1460: }

1464: PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[])
1465: {
1467:   *array = 0; /* user cannot accidently use the array later */
1468:   return(0);
1469: }

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

1476:    Not Collective

1478:    Input Parameter:
1479: .  mat - a MATSEQDENSE matrix

1481:    Output Parameter:
1482: .   array - pointer to the data

1484:    Level: intermediate

1486: .seealso: MatDenseRestoreArray()
1487: @*/
1488: PetscErrorCode  MatDenseGetArray(Mat A,PetscScalar **array)
1489: {

1493:   PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));
1494:   return(0);
1495: }

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

1502:    Not Collective

1504:    Input Parameters:
1505: .  mat - a MATSEQDENSE matrix
1506: .  array - pointer to the data

1508:    Level: intermediate

1510: .seealso: MatDenseGetArray()
1511: @*/
1512: PetscErrorCode  MatDenseRestoreArray(Mat A,PetscScalar **array)
1513: {

1517:   PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));
1518:   return(0);
1519: }

1523: static PetscErrorCode MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B)
1524: {
1525:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1527:   PetscInt       i,j,nrows,ncols;
1528:   const PetscInt *irow,*icol;
1529:   PetscScalar    *av,*bv,*v = mat->v;
1530:   Mat            newmat;

1533:   ISGetIndices(isrow,&irow);
1534:   ISGetIndices(iscol,&icol);
1535:   ISGetLocalSize(isrow,&nrows);
1536:   ISGetLocalSize(iscol,&ncols);

1538:   /* Check submatrixcall */
1539:   if (scall == MAT_REUSE_MATRIX) {
1540:     PetscInt n_cols,n_rows;
1541:     MatGetSize(*B,&n_rows,&n_cols);
1542:     if (n_rows != nrows || n_cols != ncols) {
1543:       /* resize the result matrix to match number of requested rows/columns */
1544:       MatSetSizes(*B,nrows,ncols,nrows,ncols);
1545:     }
1546:     newmat = *B;
1547:   } else {
1548:     /* Create and fill new matrix */
1549:     MatCreate(PetscObjectComm((PetscObject)A),&newmat);
1550:     MatSetSizes(newmat,nrows,ncols,nrows,ncols);
1551:     MatSetType(newmat,((PetscObject)A)->type_name);
1552:     MatSeqDenseSetPreallocation(newmat,NULL);
1553:   }

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

1558:   for (i=0; i<ncols; i++) {
1559:     av = v + mat->lda*icol[i];
1560:     for (j=0; j<nrows; j++) *bv++ = av[irow[j]];
1561:   }

1563:   /* Assemble the matrices so that the correct flags are set */
1564:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1565:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);

1567:   /* Free work space */
1568:   ISRestoreIndices(isrow,&irow);
1569:   ISRestoreIndices(iscol,&icol);
1570:   *B   = newmat;
1571:   return(0);
1572: }

1576: PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[])
1577: {
1579:   PetscInt       i;

1582:   if (scall == MAT_INITIAL_MATRIX) {
1583:     PetscMalloc1((n+1),B);
1584:   }

1586:   for (i=0; i<n; i++) {
1587:     MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
1588:   }
1589:   return(0);
1590: }

1594: PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
1595: {
1597:   return(0);
1598: }

1602: PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
1603: {
1605:   return(0);
1606: }

1610: PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
1611: {
1612:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data;
1614:   PetscInt       lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j;

1617:   /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
1618:   if (A->ops->copy != B->ops->copy) {
1619:     MatCopy_Basic(A,B,str);
1620:     return(0);
1621:   }
1622:   if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)");
1623:   if (lda1>m || lda2>m) {
1624:     for (j=0; j<n; j++) {
1625:       PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));
1626:     }
1627:   } else {
1628:     PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1629:   }
1630:   return(0);
1631: }

1635: PetscErrorCode MatSetUp_SeqDense(Mat A)
1636: {

1640:    MatSeqDenseSetPreallocation(A,0);
1641:   return(0);
1642: }

1646: static PetscErrorCode MatConjugate_SeqDense(Mat A)
1647: {
1648:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1649:   PetscInt     i,nz = A->rmap->n*A->cmap->n;
1650:   PetscScalar  *aa = a->v;

1653:   for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
1654:   return(0);
1655: }

1659: static PetscErrorCode MatRealPart_SeqDense(Mat A)
1660: {
1661:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1662:   PetscInt     i,nz = A->rmap->n*A->cmap->n;
1663:   PetscScalar  *aa = a->v;

1666:   for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
1667:   return(0);
1668: }

1672: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
1673: {
1674:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1675:   PetscInt     i,nz = A->rmap->n*A->cmap->n;
1676:   PetscScalar  *aa = a->v;

1679:   for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
1680:   return(0);
1681: }

1683: /* ----------------------------------------------------------------*/
1686: PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1687: {

1691:   if (scall == MAT_INITIAL_MATRIX) {
1692:     PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
1693:     MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
1694:     PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
1695:   }
1696:   PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
1697:   MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);
1698:   PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
1699:   return(0);
1700: }

1704: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1705: {
1707:   PetscInt       m=A->rmap->n,n=B->cmap->n;
1708:   Mat            Cmat;

1711:   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);
1712:   MatCreate(PETSC_COMM_SELF,&Cmat);
1713:   MatSetSizes(Cmat,m,n,m,n);
1714:   MatSetType(Cmat,MATSEQDENSE);
1715:   MatSeqDenseSetPreallocation(Cmat,NULL);

1717:   *C = Cmat;
1718:   return(0);
1719: }

1723: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
1724: {
1725:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
1726:   Mat_SeqDense   *b = (Mat_SeqDense*)B->data;
1727:   Mat_SeqDense   *c = (Mat_SeqDense*)C->data;
1728:   PetscBLASInt   m,n,k;
1729:   PetscScalar    _DOne=1.0,_DZero=0.0;

1733:   PetscBLASIntCast(A->rmap->n,&m);
1734:   PetscBLASIntCast(B->cmap->n,&n);
1735:   PetscBLASIntCast(A->cmap->n,&k);
1736:   PetscStackCallBLAS("BLASgemv",BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
1737:   return(0);
1738: }

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

1747:   if (scall == MAT_INITIAL_MATRIX) {
1748:     PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);
1749:     MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
1750:     PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);
1751:   }
1752:   PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);
1753:   MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);
1754:   PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);
1755:   return(0);
1756: }

1760: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1761: {
1763:   PetscInt       m=A->cmap->n,n=B->cmap->n;
1764:   Mat            Cmat;

1767:   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);
1768:   MatCreate(PETSC_COMM_SELF,&Cmat);
1769:   MatSetSizes(Cmat,m,n,m,n);
1770:   MatSetType(Cmat,MATSEQDENSE);
1771:   MatSeqDenseSetPreallocation(Cmat,NULL);

1773:   Cmat->assembled = PETSC_TRUE;

1775:   *C = Cmat;
1776:   return(0);
1777: }

1781: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
1782: {
1783:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
1784:   Mat_SeqDense   *b = (Mat_SeqDense*)B->data;
1785:   Mat_SeqDense   *c = (Mat_SeqDense*)C->data;
1786:   PetscBLASInt   m,n,k;
1787:   PetscScalar    _DOne=1.0,_DZero=0.0;

1791:   PetscBLASIntCast(A->cmap->n,&m);
1792:   PetscBLASIntCast(B->cmap->n,&n);
1793:   PetscBLASIntCast(A->rmap->n,&k);
1794:   /*
1795:      Note the m and n arguments below are the number rows and columns of A', not A!
1796:   */
1797:   PetscStackCallBLAS("BLASgemv",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
1798:   return(0);
1799: }

1803: PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[])
1804: {
1805:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
1807:   PetscInt       i,j,m = A->rmap->n,n = A->cmap->n,p;
1808:   PetscScalar    *x;
1809:   MatScalar      *aa = a->v;

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

1814:   VecSet(v,0.0);
1815:   VecGetArray(v,&x);
1816:   VecGetLocalSize(v,&p);
1817:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
1818:   for (i=0; i<m; i++) {
1819:     x[i] = aa[i]; if (idx) idx[i] = 0;
1820:     for (j=1; j<n; j++) {
1821:       if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
1822:     }
1823:   }
1824:   VecRestoreArray(v,&x);
1825:   return(0);
1826: }

1830: PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])
1831: {
1832:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
1834:   PetscInt       i,j,m = A->rmap->n,n = A->cmap->n,p;
1835:   PetscScalar    *x;
1836:   PetscReal      atmp;
1837:   MatScalar      *aa = a->v;

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

1842:   VecSet(v,0.0);
1843:   VecGetArray(v,&x);
1844:   VecGetLocalSize(v,&p);
1845:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
1846:   for (i=0; i<m; i++) {
1847:     x[i] = PetscAbsScalar(aa[i]);
1848:     for (j=1; j<n; j++) {
1849:       atmp = PetscAbsScalar(aa[i+m*j]);
1850:       if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;}
1851:     }
1852:   }
1853:   VecRestoreArray(v,&x);
1854:   return(0);
1855: }

1859: PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])
1860: {
1861:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
1863:   PetscInt       i,j,m = A->rmap->n,n = A->cmap->n,p;
1864:   PetscScalar    *x;
1865:   MatScalar      *aa = a->v;

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

1870:   VecSet(v,0.0);
1871:   VecGetArray(v,&x);
1872:   VecGetLocalSize(v,&p);
1873:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
1874:   for (i=0; i<m; i++) {
1875:     x[i] = aa[i]; if (idx) idx[i] = 0;
1876:     for (j=1; j<n; j++) {
1877:       if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
1878:     }
1879:   }
1880:   VecRestoreArray(v,&x);
1881:   return(0);
1882: }

1886: PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
1887: {
1888:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
1890:   PetscScalar    *x;

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

1895:   VecGetArray(v,&x);
1896:   PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));
1897:   VecRestoreArray(v,&x);
1898:   return(0);
1899: }


1904: PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms)
1905: {
1907:   PetscInt       i,j,m,n;
1908:   PetscScalar    *a;

1911:   MatGetSize(A,&m,&n);
1912:   PetscMemzero(norms,n*sizeof(PetscReal));
1913:   MatDenseGetArray(A,&a);
1914:   if (type == NORM_2) {
1915:     for (i=0; i<n; i++) {
1916:       for (j=0; j<m; j++) {
1917:         norms[i] += PetscAbsScalar(a[j]*a[j]);
1918:       }
1919:       a += m;
1920:     }
1921:   } else if (type == NORM_1) {
1922:     for (i=0; i<n; i++) {
1923:       for (j=0; j<m; j++) {
1924:         norms[i] += PetscAbsScalar(a[j]);
1925:       }
1926:       a += m;
1927:     }
1928:   } else if (type == NORM_INFINITY) {
1929:     for (i=0; i<n; i++) {
1930:       for (j=0; j<m; j++) {
1931:         norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]);
1932:       }
1933:       a += m;
1934:     }
1935:   } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType");
1936:   MatDenseRestoreArray(A,&a);
1937:   if (type == NORM_2) {
1938:     for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]);
1939:   }
1940:   return(0);
1941: }

1945: static PetscErrorCode  MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
1946: {
1948:   PetscScalar    *a;
1949:   PetscInt       m,n,i;

1952:   MatGetSize(x,&m,&n);
1953:   MatDenseGetArray(x,&a);
1954:   for (i=0; i<m*n; i++) {
1955:     PetscRandomGetValue(rctx,a+i);
1956:   }
1957:   MatDenseRestoreArray(x,&a);
1958:   return(0);
1959: }


1962: /* -------------------------------------------------------------------*/
1963: static struct _MatOps MatOps_Values = { MatSetValues_SeqDense,
1964:                                         MatGetRow_SeqDense,
1965:                                         MatRestoreRow_SeqDense,
1966:                                         MatMult_SeqDense,
1967:                                 /*  4*/ MatMultAdd_SeqDense,
1968:                                         MatMultTranspose_SeqDense,
1969:                                         MatMultTransposeAdd_SeqDense,
1970:                                         0,
1971:                                         0,
1972:                                         0,
1973:                                 /* 10*/ 0,
1974:                                         MatLUFactor_SeqDense,
1975:                                         MatCholeskyFactor_SeqDense,
1976:                                         MatSOR_SeqDense,
1977:                                         MatTranspose_SeqDense,
1978:                                 /* 15*/ MatGetInfo_SeqDense,
1979:                                         MatEqual_SeqDense,
1980:                                         MatGetDiagonal_SeqDense,
1981:                                         MatDiagonalScale_SeqDense,
1982:                                         MatNorm_SeqDense,
1983:                                 /* 20*/ MatAssemblyBegin_SeqDense,
1984:                                         MatAssemblyEnd_SeqDense,
1985:                                         MatSetOption_SeqDense,
1986:                                         MatZeroEntries_SeqDense,
1987:                                 /* 24*/ MatZeroRows_SeqDense,
1988:                                         0,
1989:                                         0,
1990:                                         0,
1991:                                         0,
1992:                                 /* 29*/ MatSetUp_SeqDense,
1993:                                         0,
1994:                                         0,
1995:                                         0,
1996:                                         0,
1997:                                 /* 34*/ MatDuplicate_SeqDense,
1998:                                         0,
1999:                                         0,
2000:                                         0,
2001:                                         0,
2002:                                 /* 39*/ MatAXPY_SeqDense,
2003:                                         MatGetSubMatrices_SeqDense,
2004:                                         0,
2005:                                         MatGetValues_SeqDense,
2006:                                         MatCopy_SeqDense,
2007:                                 /* 44*/ MatGetRowMax_SeqDense,
2008:                                         MatScale_SeqDense,
2009:                                         0,
2010:                                         0,
2011:                                         0,
2012:                                 /* 49*/ MatSetRandom_SeqDense,
2013:                                         0,
2014:                                         0,
2015:                                         0,
2016:                                         0,
2017:                                 /* 54*/ 0,
2018:                                         0,
2019:                                         0,
2020:                                         0,
2021:                                         0,
2022:                                 /* 59*/ 0,
2023:                                         MatDestroy_SeqDense,
2024:                                         MatView_SeqDense,
2025:                                         0,
2026:                                         0,
2027:                                 /* 64*/ 0,
2028:                                         0,
2029:                                         0,
2030:                                         0,
2031:                                         0,
2032:                                 /* 69*/ MatGetRowMaxAbs_SeqDense,
2033:                                         0,
2034:                                         0,
2035:                                         0,
2036:                                         0,
2037:                                 /* 74*/ 0,
2038:                                         0,
2039:                                         0,
2040:                                         0,
2041:                                         0,
2042:                                 /* 79*/ 0,
2043:                                         0,
2044:                                         0,
2045:                                         0,
2046:                                 /* 83*/ MatLoad_SeqDense,
2047:                                         0,
2048:                                         MatIsHermitian_SeqDense,
2049:                                         0,
2050:                                         0,
2051:                                         0,
2052:                                 /* 89*/ MatMatMult_SeqDense_SeqDense,
2053:                                         MatMatMultSymbolic_SeqDense_SeqDense,
2054:                                         MatMatMultNumeric_SeqDense_SeqDense,
2055:                                         0,
2056:                                         0,
2057:                                 /* 94*/ 0,
2058:                                         0,
2059:                                         0,
2060:                                         0,
2061:                                         0,
2062:                                 /* 99*/ 0,
2063:                                         0,
2064:                                         0,
2065:                                         MatConjugate_SeqDense,
2066:                                         0,
2067:                                 /*104*/ 0,
2068:                                         MatRealPart_SeqDense,
2069:                                         MatImaginaryPart_SeqDense,
2070:                                         0,
2071:                                         0,
2072:                                 /*109*/ MatMatSolve_SeqDense,
2073:                                         0,
2074:                                         MatGetRowMin_SeqDense,
2075:                                         MatGetColumnVector_SeqDense,
2076:                                         0,
2077:                                 /*114*/ 0,
2078:                                         0,
2079:                                         0,
2080:                                         0,
2081:                                         0,
2082:                                 /*119*/ 0,
2083:                                         0,
2084:                                         0,
2085:                                         0,
2086:                                         0,
2087:                                 /*124*/ 0,
2088:                                         MatGetColumnNorms_SeqDense,
2089:                                         0,
2090:                                         0,
2091:                                         0,
2092:                                 /*129*/ 0,
2093:                                         MatTransposeMatMult_SeqDense_SeqDense,
2094:                                         MatTransposeMatMultSymbolic_SeqDense_SeqDense,
2095:                                         MatTransposeMatMultNumeric_SeqDense_SeqDense,
2096:                                         0,
2097:                                 /*134*/ 0,
2098:                                         0,
2099:                                         0,
2100:                                         0,
2101:                                         0,
2102:                                 /*139*/ 0,
2103:                                         0,
2104:                                         0
2105: };

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

2114:    Collective on MPI_Comm

2116:    Input Parameters:
2117: +  comm - MPI communicator, set to PETSC_COMM_SELF
2118: .  m - number of rows
2119: .  n - number of columns
2120: -  data - optional location of matrix data in column major order.  Set data=NULL for PETSc
2121:    to control all matrix memory allocation.

2123:    Output Parameter:
2124: .  A - the matrix

2126:    Notes:
2127:    The data input variable is intended primarily for Fortran programmers
2128:    who wish to allocate their own matrix memory space.  Most users should
2129:    set data=NULL.

2131:    Level: intermediate

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

2135: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
2136: @*/
2137: PetscErrorCode  MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
2138: {

2142:   MatCreate(comm,A);
2143:   MatSetSizes(*A,m,n,m,n);
2144:   MatSetType(*A,MATSEQDENSE);
2145:   MatSeqDenseSetPreallocation(*A,data);
2146:   return(0);
2147: }

2151: /*@C
2152:    MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements

2154:    Collective on MPI_Comm

2156:    Input Parameters:
2157: +  A - the matrix
2158: -  data - the array (or NULL)

2160:    Notes:
2161:    The data input variable is intended primarily for Fortran programmers
2162:    who wish to allocate their own matrix memory space.  Most users should
2163:    need not call this routine.

2165:    Level: intermediate

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

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

2171: @*/
2172: PetscErrorCode  MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
2173: {

2177:   PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
2178:   return(0);
2179: }

2183: PetscErrorCode  MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
2184: {
2185:   Mat_SeqDense   *b;

2189:   B->preallocated = PETSC_TRUE;

2191:   PetscLayoutSetUp(B->rmap);
2192:   PetscLayoutSetUp(B->cmap);

2194:   b       = (Mat_SeqDense*)B->data;
2195:   b->Mmax = B->rmap->n;
2196:   b->Nmax = B->cmap->n;
2197:   if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n;

2199:   if (!data) { /* petsc-allocated storage */
2200:     if (!b->user_alloc) { PetscFree(b->v); }
2201:     PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);
2202:     PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));

2204:     b->user_alloc = PETSC_FALSE;
2205:   } else { /* user-allocated storage */
2206:     if (!b->user_alloc) { PetscFree(b->v); }
2207:     b->v          = data;
2208:     b->user_alloc = PETSC_TRUE;
2209:   }
2210:   B->assembled = PETSC_TRUE;
2211:   return(0);
2212: }

2216: /*@C
2217:   MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array

2219:   Input parameter:
2220: + A - the matrix
2221: - lda - the leading dimension

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

2228:   Level: intermediate

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

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

2234: @*/
2235: PetscErrorCode  MatSeqDenseSetLDA(Mat B,PetscInt lda)
2236: {
2237:   Mat_SeqDense *b = (Mat_SeqDense*)B->data;

2240:   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);
2241:   b->lda       = lda;
2242:   b->changelda = PETSC_FALSE;
2243:   b->Mmax      = PetscMax(b->Mmax,lda);
2244:   return(0);
2245: }

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

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

2253:   Level: beginner

2255: .seealso: MatCreateSeqDense()

2257: M*/

2261: PETSC_EXTERN PetscErrorCode MatCreate_SeqDense(Mat B)
2262: {
2263:   Mat_SeqDense   *b;
2265:   PetscMPIInt    size;

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

2271:   PetscNewLog(B,&b);
2272:   PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
2273:   B->data = (void*)b;

2275:   b->pivots      = 0;
2276:   b->roworiented = PETSC_TRUE;
2277:   b->v           = 0;
2278:   b->changelda   = PETSC_FALSE;

2280:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);
2281:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);

2283:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);
2284:   PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_petsc_C",MatGetFactor_seqdense_petsc);
2285:   PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);
2286:   PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2287:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2288:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2289:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2290:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2291:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2292:   PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
2293:   return(0);
2294: }