Actual source code: dense.c

petsc-master 2015-07-03
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_SeqAIJ_SeqDense(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
 14: {
 15:   Mat            B;
 16:   Mat_SeqAIJ     *a = (Mat_SeqAIJ*)A->data;
 17:   Mat_SeqDense   *b;
 19:   PetscInt       *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->N,i;
 20:   MatScalar      *av=a->a;

 23:   MatCreate(PetscObjectComm((PetscObject)A),&B);
 24:   MatSetSizes(B,m,n,m,n);
 25:   MatSetType(B,MATSEQDENSE);
 26:   MatSeqDenseSetPreallocation(B,NULL);

 28:   b  = (Mat_SeqDense*)(B->data);
 29:   for (i=0; i<m; i++) {
 30:     PetscInt j;
 31:     for (j=0;j<ai[1]-ai[0];j++) {
 32:       b->v[*aj*m+i] = *av;
 33:       aj++;
 34:       av++;
 35:     }
 36:     ai++;
 37:   }
 38:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 39:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

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

 51: PETSC_EXTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
 52: {
 53:   Mat            B;
 54:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
 56:   PetscInt       i, j;
 57:   PetscInt       *rows, *nnz;
 58:   MatScalar      *aa = a->v, *vals;

 61:   MatCreate(PetscObjectComm((PetscObject)A),&B);
 62:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
 63:   MatSetType(B,MATSEQAIJ);
 64:   PetscCalloc3(A->rmap->n,&rows,A->rmap->n,&nnz,A->rmap->n,&vals);
 65:   for (j=0; j<A->cmap->n; j++) {
 66:     for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) ++nnz[i];
 67:     aa += a->lda;
 68:   }
 69:   MatSeqAIJSetPreallocation(B,PETSC_DETERMINE,nnz);
 70:   aa = a->v;
 71:   for (j=0; j<A->cmap->n; j++) {
 72:     PetscInt numRows = 0;
 73:     for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) {rows[numRows] = i; vals[numRows++] = aa[i];}
 74:     MatSetValues(B,numRows,rows,1,&j,vals,INSERT_VALUES);
 75:     aa  += a->lda;
 76:   }
 77:   PetscFree3(rows,nnz,vals);
 78:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 79:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

 81:   if (reuse == MAT_REUSE_MATRIX) {
 82:     MatHeaderReplace(A,B);
 83:   } else {
 84:     *newmat = B;
 85:   }
 86:   return(0);
 87: }

 91: PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str)
 92: {
 93:   Mat_SeqDense   *x     = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data;
 94:   PetscScalar    oalpha = alpha;
 95:   PetscInt       j;
 96:   PetscBLASInt   N,m,ldax,lday,one = 1;

100:   PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);
101:   PetscBLASIntCast(X->rmap->n,&m);
102:   PetscBLASIntCast(x->lda,&ldax);
103:   PetscBLASIntCast(y->lda,&lday);
104:   if (ldax>m || lday>m) {
105:     for (j=0; j<X->cmap->n; j++) {
106:       PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&oalpha,x->v+j*ldax,&one,y->v+j*lday,&one));
107:     }
108:   } else {
109:     PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&oalpha,x->v,&one,y->v,&one));
110:   }
111:   PetscObjectStateIncrease((PetscObject)Y);
112:   PetscLogFlops(PetscMax(2*N-1,0));
113:   return(0);
114: }

118: PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info)
119: {
120:   PetscInt N = A->rmap->n*A->cmap->n;

123:   info->block_size        = 1.0;
124:   info->nz_allocated      = (double)N;
125:   info->nz_used           = (double)N;
126:   info->nz_unneeded       = (double)0;
127:   info->assemblies        = (double)A->num_ass;
128:   info->mallocs           = 0;
129:   info->memory            = ((PetscObject)A)->mem;
130:   info->fill_ratio_given  = 0;
131:   info->fill_ratio_needed = 0;
132:   info->factor_mallocs    = 0;
133:   return(0);
134: }

138: PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha)
139: {
140:   Mat_SeqDense   *a     = (Mat_SeqDense*)A->data;
141:   PetscScalar    oalpha = alpha;
143:   PetscBLASInt   one = 1,j,nz,lda;

146:   PetscBLASIntCast(a->lda,&lda);
147:   if (lda>A->rmap->n) {
148:     PetscBLASIntCast(A->rmap->n,&nz);
149:     for (j=0; j<A->cmap->n; j++) {
150:       PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v+j*lda,&one));
151:     }
152:   } else {
153:     PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);
154:     PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v,&one));
155:   }
156:   PetscLogFlops(nz);
157:   return(0);
158: }

162: PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool  *fl)
163: {
164:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
165:   PetscInt     i,j,m = A->rmap->n,N;
166:   PetscScalar  *v = a->v;

169:   *fl = PETSC_FALSE;
170:   if (A->rmap->n != A->cmap->n) return(0);
171:   N = a->lda;

173:   for (i=0; i<m; i++) {
174:     for (j=i+1; j<m; j++) {
175:       if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) return(0);
176:     }
177:   }
178:   *fl = PETSC_TRUE;
179:   return(0);
180: }

184: PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues)
185: {
186:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data,*l;
188:   PetscInt       lda = (PetscInt)mat->lda,j,m;

191:   PetscLayoutReference(A->rmap,&newi->rmap);
192:   PetscLayoutReference(A->cmap,&newi->cmap);
193:   MatSeqDenseSetPreallocation(newi,NULL);
194:   if (cpvalues == MAT_COPY_VALUES) {
195:     l = (Mat_SeqDense*)newi->data;
196:     if (lda>A->rmap->n) {
197:       m = A->rmap->n;
198:       for (j=0; j<A->cmap->n; j++) {
199:         PetscMemcpy(l->v+j*m,mat->v+j*lda,m*sizeof(PetscScalar));
200:       }
201:     } else {
202:       PetscMemcpy(l->v,mat->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
203:     }
204:   }
205:   newi->assembled = PETSC_TRUE;
206:   return(0);
207: }

211: PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat)
212: {

216:   MatCreate(PetscObjectComm((PetscObject)A),newmat);
217:   MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
218:   MatSetType(*newmat,((PetscObject)A)->type_name);
219:   MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);
220:   return(0);
221: }


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

228: PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
229: {
230:   MatFactorInfo  info;

234:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
235:   MatLUFactor_SeqDense(fact,0,0,&info);
236:   return(0);
237: }

241: PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy)
242: {
243:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
244:   PetscErrorCode    ierr;
245:   const PetscScalar *x;
246:   PetscScalar       *y;
247:   PetscBLASInt      one = 1,info,m;

250:   PetscBLASIntCast(A->rmap->n,&m);
251:   VecGetArrayRead(xx,&x);
252:   VecGetArray(yy,&y);
253:   PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
254:   if (A->factortype == MAT_FACTOR_LU) {
255: #if defined(PETSC_MISSING_LAPACK_GETRS)
256:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
257: #else
258:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
259:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
260: #endif
261:   } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
262: #if defined(PETSC_MISSING_LAPACK_POTRS)
263:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
264: #else
265:     PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
266:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
267: #endif
268:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
269:   VecRestoreArrayRead(xx,&x);
270:   VecRestoreArray(yy,&y);
271:   PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
272:   return(0);
273: }

277: PetscErrorCode MatMatSolve_SeqDense(Mat A,Mat B,Mat X)
278: {
279:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
281:   PetscScalar    *b,*x;
282:   PetscInt       n;
283:   PetscBLASInt   nrhs,info,m;
284:   PetscBool      flg;

287:   PetscBLASIntCast(A->rmap->n,&m);
288:   PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
289:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
290:   PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
291:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");

293:   MatGetSize(B,NULL,&n);
294:   PetscBLASIntCast(n,&nrhs);
295:   MatDenseGetArray(B,&b);
296:   MatDenseGetArray(X,&x);

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

300:   if (A->factortype == MAT_FACTOR_LU) {
301: #if defined(PETSC_MISSING_LAPACK_GETRS)
302:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
303: #else
304:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
305:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
306: #endif
307:   } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
308: #if defined(PETSC_MISSING_LAPACK_POTRS)
309:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
310: #else
311:     PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info));
312:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
313: #endif
314:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");

316:   MatDenseRestoreArray(B,&b);
317:   MatDenseRestoreArray(X,&x);
318:   PetscLogFlops(nrhs*(2.0*m*m - m));
319:   return(0);
320: }

324: PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy)
325: {
326:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
327:   PetscErrorCode    ierr;
328:   const PetscScalar *x;
329:   PetscScalar       *y;
330:   PetscBLASInt      one = 1,info,m;

333:   PetscBLASIntCast(A->rmap->n,&m);
334:   VecGetArrayRead(xx,&x);
335:   VecGetArray(yy,&y);
336:   PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
337:   /* assume if pivots exist then use LU; else Cholesky */
338:   if (mat->pivots) {
339: #if defined(PETSC_MISSING_LAPACK_GETRS)
340:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
341: #else
342:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
343:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve");
344: #endif
345:   } else {
346: #if defined(PETSC_MISSING_LAPACK_POTRS)
347:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
348: #else
349:     PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
350:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve");
351: #endif
352:   }
353:   VecRestoreArrayRead(xx,&x);
354:   VecRestoreArray(yy,&y);
355:   PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
356:   return(0);
357: }

361: PetscErrorCode MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
362: {
363:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
364:   PetscErrorCode    ierr;
365:   const PetscScalar *x;
366:   PetscScalar       *y,sone = 1.0;
367:   Vec               tmp = 0;
368:   PetscBLASInt      one = 1,info,m;

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

411: PetscErrorCode MatSolveTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
412: {
413:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
414:   PetscErrorCode    ierr;
415:   const PetscScalar *x;
416:   PetscScalar       *y,sone = 1.0;
417:   Vec               tmp;
418:   PetscBLASInt      one = 1,info,m;

421:   PetscBLASIntCast(A->rmap->n,&m);
422:   if (!A->rmap->n || !A->cmap->n) return(0);
423:   VecGetArrayRead(xx,&x);
424:   VecGetArray(yy,&y);
425:   if (yy == zz) {
426:     VecDuplicate(yy,&tmp);
427:     PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);
428:     VecCopy(yy,tmp);
429:   }
430:   PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
431:   /* assume if pivots exist then use LU; else Cholesky */
432:   if (mat->pivots) {
433: #if defined(PETSC_MISSING_LAPACK_GETRS)
434:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
435: #else
436:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
437:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
438: #endif
439:   } else {
440: #if defined(PETSC_MISSING_LAPACK_POTRS)
441:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
442: #else
443:     PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
444:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
445: #endif
446:   }
447:   if (tmp) {
448:     VecAXPY(yy,sone,tmp);
449:     VecDestroy(&tmp);
450:   } else {
451:     VecAXPY(yy,sone,zz);
452:   }
453:   VecRestoreArrayRead(xx,&x);
454:   VecRestoreArray(yy,&y);
455:   PetscLogFlops(2.0*A->cmap->n*A->cmap->n);
456:   return(0);
457: }

459: /* ---------------------------------------------------------------*/
460: /* COMMENT: I have chosen to hide row permutation in the pivots,
461:    rather than put it in the Mat->row slot.*/
464: PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo)
465: {
466: #if defined(PETSC_MISSING_LAPACK_GETRF)
468:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable.");
469: #else
470:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
472:   PetscBLASInt   n,m,info;

475:   PetscBLASIntCast(A->cmap->n,&n);
476:   PetscBLASIntCast(A->rmap->n,&m);
477:   if (!mat->pivots) {
478:     PetscMalloc1(A->rmap->n+1,&mat->pivots);
479:     PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
480:   }
481:   if (!A->rmap->n || !A->cmap->n) return(0);
482:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
483:   PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info));
484:   PetscFPTrapPop();

486:   if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization");
487:   if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization");
488:   A->ops->solve             = MatSolve_SeqDense;
489:   A->ops->solvetranspose    = MatSolveTranspose_SeqDense;
490:   A->ops->solveadd          = MatSolveAdd_SeqDense;
491:   A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
492:   A->factortype             = MAT_FACTOR_LU;

494:   PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);
495: #endif
496:   return(0);
497: }

501: PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo)
502: {
503: #if defined(PETSC_MISSING_LAPACK_POTRF)
505:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable.");
506: #else
507:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
509:   PetscBLASInt   info,n;

512:   PetscBLASIntCast(A->cmap->n,&n);
513:   PetscFree(mat->pivots);

515:   if (!A->rmap->n || !A->cmap->n) return(0);
516:   PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info));
517:   if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1);
518:   A->ops->solve             = MatSolve_SeqDense;
519:   A->ops->solvetranspose    = MatSolveTranspose_SeqDense;
520:   A->ops->solveadd          = MatSolveAdd_SeqDense;
521:   A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
522:   A->factortype             = MAT_FACTOR_CHOLESKY;

524:   PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
525: #endif
526:   return(0);
527: }


532: PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
533: {
535:   MatFactorInfo  info;

538:   info.fill = 1.0;

540:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
541:   MatCholeskyFactor_SeqDense(fact,0,&info);
542:   return(0);
543: }

547: PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
548: {
550:   fact->assembled                  = PETSC_TRUE;
551:   fact->preallocated               = PETSC_TRUE;
552:   fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
553:   return(0);
554: }

558: PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info)
559: {
561:   fact->preallocated         = PETSC_TRUE;
562:   fact->assembled            = PETSC_TRUE;
563:   fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense;
564:   return(0);
565: }

569: PETSC_EXTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact)
570: {

574:   MatCreate(PetscObjectComm((PetscObject)A),fact);
575:   MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
576:   MatSetType(*fact,((PetscObject)A)->type_name);
577:   if (ftype == MAT_FACTOR_LU) {
578:     (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
579:   } else {
580:     (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
581:   }
582:   (*fact)->factortype = ftype;
583:   return(0);
584: }

586: /* ------------------------------------------------------------------*/
589: PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx)
590: {
591:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
592:   PetscScalar       *x,*v = mat->v,zero = 0.0,xt;
593:   const PetscScalar *b;
594:   PetscErrorCode    ierr;
595:   PetscInt          m = A->rmap->n,i;
596:   PetscBLASInt      o = 1,bm;

599:   if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */
600:   PetscBLASIntCast(m,&bm);
601:   if (flag & SOR_ZERO_INITIAL_GUESS) {
602:     /* this is a hack fix, should have another version without the second BLASdot */
603:     VecSet(xx,zero);
604:   }
605:   VecGetArray(xx,&x);
606:   VecGetArrayRead(bb,&b);
607:   its  = its*lits;
608:   if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
609:   while (its--) {
610:     if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
611:       for (i=0; i<m; i++) {
612:         PetscStackCallBLAS("BLASdot",xt   = b[i] - BLASdot_(&bm,v+i,&bm,x,&o));
613:         x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
614:       }
615:     }
616:     if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
617:       for (i=m-1; i>=0; i--) {
618:         PetscStackCallBLAS("BLASdot",xt   = b[i] - BLASdot_(&bm,v+i,&bm,x,&o));
619:         x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
620:       }
621:     }
622:   }
623:   VecRestoreArrayRead(bb,&b);
624:   VecRestoreArray(xx,&x);
625:   return(0);
626: }

628: /* -----------------------------------------------------------------*/
631: PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)
632: {
633:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
634:   const PetscScalar *v   = mat->v,*x;
635:   PetscScalar       *y;
636:   PetscErrorCode    ierr;
637:   PetscBLASInt      m, n,_One=1;
638:   PetscScalar       _DOne=1.0,_DZero=0.0;

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

655: PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy)
656: {
657:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
658:   PetscScalar       *y,_DOne=1.0,_DZero=0.0;
659:   PetscErrorCode    ierr;
660:   PetscBLASInt      m, n, _One=1;
661:   const PetscScalar *v = mat->v,*x;

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

678: PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
679: {
680:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
681:   const PetscScalar *v = mat->v,*x;
682:   PetscScalar       *y,_DOne=1.0;
683:   PetscErrorCode    ierr;
684:   PetscBLASInt      m, n, _One=1;

687:   PetscBLASIntCast(A->rmap->n,&m);
688:   PetscBLASIntCast(A->cmap->n,&n);
689:   if (!A->rmap->n || !A->cmap->n) return(0);
690:   if (zz != yy) {VecCopy(zz,yy);}
691:   VecGetArrayRead(xx,&x);
692:   VecGetArray(yy,&y);
693:   PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
694:   VecRestoreArrayRead(xx,&x);
695:   VecRestoreArray(yy,&y);
696:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
697:   return(0);
698: }

702: PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
703: {
704:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
705:   const PetscScalar *v = mat->v,*x;
706:   PetscScalar       *y;
707:   PetscErrorCode    ierr;
708:   PetscBLASInt      m, n, _One=1;
709:   PetscScalar       _DOne=1.0;

712:   PetscBLASIntCast(A->rmap->n,&m);
713:   PetscBLASIntCast(A->cmap->n,&n);
714:   if (!A->rmap->n || !A->cmap->n) return(0);
715:   if (zz != yy) {VecCopy(zz,yy);}
716:   VecGetArrayRead(xx,&x);
717:   VecGetArray(yy,&y);
718:   PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
719:   VecRestoreArrayRead(xx,&x);
720:   VecRestoreArray(yy,&y);
721:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
722:   return(0);
723: }

725: /* -----------------------------------------------------------------*/
728: PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
729: {
730:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
731:   PetscScalar    *v;
733:   PetscInt       i;

736:   *ncols = A->cmap->n;
737:   if (cols) {
738:     PetscMalloc1(A->cmap->n+1,cols);
739:     for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
740:   }
741:   if (vals) {
742:     PetscMalloc1(A->cmap->n+1,vals);
743:     v    = mat->v + row;
744:     for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
745:   }
746:   return(0);
747: }

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

756:   if (cols) {PetscFree(*cols);}
757:   if (vals) {PetscFree(*vals); }
758:   return(0);
759: }
760: /* ----------------------------------------------------------------*/
763: PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
764: {
765:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
766:   PetscInt     i,j,idx=0;

769:   if (!mat->roworiented) {
770:     if (addv == INSERT_VALUES) {
771:       for (j=0; j<n; j++) {
772:         if (indexn[j] < 0) {idx += m; 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:         for (i=0; i<m; i++) {
777:           if (indexm[i] < 0) {idx++; continue;}
778: #if defined(PETSC_USE_DEBUG)
779:           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);
780: #endif
781:           mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
782:         }
783:       }
784:     } else {
785:       for (j=0; j<n; j++) {
786:         if (indexn[j] < 0) {idx += m; continue;}
787: #if defined(PETSC_USE_DEBUG)
788:         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);
789: #endif
790:         for (i=0; i<m; i++) {
791:           if (indexm[i] < 0) {idx++; continue;}
792: #if defined(PETSC_USE_DEBUG)
793:           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);
794: #endif
795:           mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
796:         }
797:       }
798:     }
799:   } else {
800:     if (addv == INSERT_VALUES) {
801:       for (i=0; i<m; i++) {
802:         if (indexm[i] < 0) { idx += n; continue;}
803: #if defined(PETSC_USE_DEBUG)
804:         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);
805: #endif
806:         for (j=0; j<n; j++) {
807:           if (indexn[j] < 0) { idx++; continue;}
808: #if defined(PETSC_USE_DEBUG)
809:           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);
810: #endif
811:           mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
812:         }
813:       }
814:     } else {
815:       for (i=0; i<m; i++) {
816:         if (indexm[i] < 0) { idx += n; continue;}
817: #if defined(PETSC_USE_DEBUG)
818:         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);
819: #endif
820:         for (j=0; j<n; j++) {
821:           if (indexn[j] < 0) { idx++; continue;}
822: #if defined(PETSC_USE_DEBUG)
823:           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);
824: #endif
825:           mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
826:         }
827:       }
828:     }
829:   }
830:   return(0);
831: }

835: PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
836: {
837:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
838:   PetscInt     i,j;

841:   /* row-oriented output */
842:   for (i=0; i<m; i++) {
843:     if (indexm[i] < 0) {v += n;continue;}
844:     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);
845:     for (j=0; j<n; j++) {
846:       if (indexn[j] < 0) {v++; continue;}
847:       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);
848:       *v++ = mat->v[indexn[j]*mat->lda + indexm[i]];
849:     }
850:   }
851:   return(0);
852: }

854: /* -----------------------------------------------------------------*/

858: PetscErrorCode MatLoad_SeqDense(Mat newmat,PetscViewer viewer)
859: {
860:   Mat_SeqDense   *a;
862:   PetscInt       *scols,i,j,nz,header[4];
863:   int            fd;
864:   PetscMPIInt    size;
865:   PetscInt       *rowlengths = 0,M,N,*cols,grows,gcols;
866:   PetscScalar    *vals,*svals,*v,*w;
867:   MPI_Comm       comm;

870:   /* force binary viewer to load .info file if it has not yet done so */
871:   PetscViewerSetUp(viewer);
872:   PetscObjectGetComm((PetscObject)viewer,&comm);
873:   MPI_Comm_size(comm,&size);
874:   if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor");
875:   PetscViewerBinaryGetDescriptor(viewer,&fd);
876:   PetscBinaryRead(fd,header,4,PETSC_INT);
877:   if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not matrix object");
878:   M = header[1]; N = header[2]; nz = header[3];

880:   /* set global size if not set already*/
881:   if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) {
882:     MatSetSizes(newmat,M,N,M,N);
883:   } else {
884:     /* if sizes and type are already set, check if the vector global sizes are correct */
885:     MatGetSize(newmat,&grows,&gcols);
886:     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);
887:   }
888:   a = (Mat_SeqDense*)newmat->data;
889:   if (!a->user_alloc) {
890:     MatSeqDenseSetPreallocation(newmat,NULL);
891:   }

893:   if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */
894:     a = (Mat_SeqDense*)newmat->data;
895:     v = a->v;
896:     /* Allocate some temp space to read in the values and then flip them
897:        from row major to column major */
898:     PetscMalloc1(M*N > 0 ? M*N : 1,&w);
899:     /* read in nonzero values */
900:     PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);
901:     /* now flip the values and store them in the matrix*/
902:     for (j=0; j<N; j++) {
903:       for (i=0; i<M; i++) {
904:         *v++ =w[i*N+j];
905:       }
906:     }
907:     PetscFree(w);
908:   } else {
909:     /* read row lengths */
910:     PetscMalloc1(M+1,&rowlengths);
911:     PetscBinaryRead(fd,rowlengths,M,PETSC_INT);

913:     a = (Mat_SeqDense*)newmat->data;
914:     v = a->v;

916:     /* read column indices and nonzeros */
917:     PetscMalloc1(nz+1,&scols);
918:     cols = scols;
919:     PetscBinaryRead(fd,cols,nz,PETSC_INT);
920:     PetscMalloc1(nz+1,&svals);
921:     vals = svals;
922:     PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);

924:     /* insert into matrix */
925:     for (i=0; i<M; i++) {
926:       for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j];
927:       svals += rowlengths[i]; scols += rowlengths[i];
928:     }
929:     PetscFree(vals);
930:     PetscFree(cols);
931:     PetscFree(rowlengths);
932:   }
933:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
934:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
935:   return(0);
936: }

940: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
941: {
942:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
943:   PetscErrorCode    ierr;
944:   PetscInt          i,j;
945:   const char        *name;
946:   PetscScalar       *v;
947:   PetscViewerFormat format;
948: #if defined(PETSC_USE_COMPLEX)
949:   PetscBool         allreal = PETSC_TRUE;
950: #endif

953:   PetscViewerGetFormat(viewer,&format);
954:   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
955:     return(0);  /* do nothing for now */
956:   } else if (format == PETSC_VIEWER_ASCII_COMMON) {
957:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
958:     for (i=0; i<A->rmap->n; i++) {
959:       v    = a->v + i;
960:       PetscViewerASCIIPrintf(viewer,"row %D:",i);
961:       for (j=0; j<A->cmap->n; j++) {
962: #if defined(PETSC_USE_COMPLEX)
963:         if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
964:           PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
965:         } else if (PetscRealPart(*v)) {
966:           PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));
967:         }
968: #else
969:         if (*v) {
970:           PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);
971:         }
972: #endif
973:         v += a->lda;
974:       }
975:       PetscViewerASCIIPrintf(viewer,"\n");
976:     }
977:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
978:   } else {
979:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
980: #if defined(PETSC_USE_COMPLEX)
981:     /* determine if matrix has all real values */
982:     v = a->v;
983:     for (i=0; i<A->rmap->n*A->cmap->n; i++) {
984:       if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;}
985:     }
986: #endif
987:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
988:       PetscObjectGetName((PetscObject)A,&name);
989:       PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);
990:       PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);
991:       PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
992:     }

994:     for (i=0; i<A->rmap->n; i++) {
995:       v = a->v + i;
996:       for (j=0; j<A->cmap->n; j++) {
997: #if defined(PETSC_USE_COMPLEX)
998:         if (allreal) {
999:           PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));
1000:         } else {
1001:           PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1002:         }
1003: #else
1004:         PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);
1005: #endif
1006:         v += a->lda;
1007:       }
1008:       PetscViewerASCIIPrintf(viewer,"\n");
1009:     }
1010:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
1011:       PetscViewerASCIIPrintf(viewer,"];\n");
1012:     }
1013:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1014:   }
1015:   PetscViewerFlush(viewer);
1016:   return(0);
1017: }

1021: static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer)
1022: {
1023:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
1024:   PetscErrorCode    ierr;
1025:   int               fd;
1026:   PetscInt          ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n;
1027:   PetscScalar       *v,*anonz,*vals;
1028:   PetscViewerFormat format;

1031:   PetscViewerBinaryGetDescriptor(viewer,&fd);

1033:   PetscViewerGetFormat(viewer,&format);
1034:   if (format == PETSC_VIEWER_NATIVE) {
1035:     /* store the matrix as a dense matrix */
1036:     PetscMalloc1(4,&col_lens);

1038:     col_lens[0] = MAT_FILE_CLASSID;
1039:     col_lens[1] = m;
1040:     col_lens[2] = n;
1041:     col_lens[3] = MATRIX_BINARY_FORMAT_DENSE;

1043:     PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);
1044:     PetscFree(col_lens);

1046:     /* write out matrix, by rows */
1047:     PetscMalloc1(m*n+1,&vals);
1048:     v    = a->v;
1049:     for (j=0; j<n; j++) {
1050:       for (i=0; i<m; i++) {
1051:         vals[j + i*n] = *v++;
1052:       }
1053:     }
1054:     PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);
1055:     PetscFree(vals);
1056:   } else {
1057:     PetscMalloc1(4+nz,&col_lens);

1059:     col_lens[0] = MAT_FILE_CLASSID;
1060:     col_lens[1] = m;
1061:     col_lens[2] = n;
1062:     col_lens[3] = nz;

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

1068:     /* Possibly should write in smaller increments, not whole matrix at once? */
1069:     /* store column indices (zero start index) */
1070:     ict = 0;
1071:     for (i=0; i<m; i++) {
1072:       for (j=0; j<n; j++) col_lens[ict++] = j;
1073:     }
1074:     PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);
1075:     PetscFree(col_lens);

1077:     /* store nonzero values */
1078:     PetscMalloc1(nz+1,&anonz);
1079:     ict  = 0;
1080:     for (i=0; i<m; i++) {
1081:       v = a->v + i;
1082:       for (j=0; j<n; j++) {
1083:         anonz[ict++] = *v; v += a->lda;
1084:       }
1085:     }
1086:     PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);
1087:     PetscFree(anonz);
1088:   }
1089:   return(0);
1090: }

1092: #include <petscdraw.h>
1095: PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa)
1096: {
1097:   Mat               A  = (Mat) Aa;
1098:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
1099:   PetscErrorCode    ierr;
1100:   PetscInt          m  = A->rmap->n,n = A->cmap->n,color,i,j;
1101:   PetscScalar       *v = a->v;
1102:   PetscViewer       viewer;
1103:   PetscDraw         popup;
1104:   PetscReal         xl,yl,xr,yr,x_l,x_r,y_l,y_r,scale,maxv = 0.0;
1105:   PetscViewerFormat format;

1108:   PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
1109:   PetscViewerGetFormat(viewer,&format);
1110:   PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);

1112:   /* Loop over matrix elements drawing boxes */
1113:   if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1114:     /* Blue for negative and Red for positive */
1115:     color = PETSC_DRAW_BLUE;
1116:     for (j = 0; j < n; j++) {
1117:       x_l = j;
1118:       x_r = x_l + 1.0;
1119:       for (i = 0; i < m; i++) {
1120:         y_l = m - i - 1.0;
1121:         y_r = y_l + 1.0;
1122:         if (PetscRealPart(v[j*m+i]) >  0.) {
1123:           color = PETSC_DRAW_RED;
1124:         } else if (PetscRealPart(v[j*m+i]) <  0.) {
1125:           color = PETSC_DRAW_BLUE;
1126:         } else {
1127:           continue;
1128:         }
1129:         PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1130:       }
1131:     }
1132:   } else {
1133:     /* use contour shading to indicate magnitude of values */
1134:     /* first determine max of all nonzero values */
1135:     for (i = 0; i < m*n; i++) {
1136:       if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1137:     }
1138:     scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv;
1139:     PetscDrawGetPopup(draw,&popup);
1140:     if (popup) {PetscDrawScalePopup(popup,0.0,maxv);}
1141:     for (j = 0; j < n; j++) {
1142:       x_l = j;
1143:       x_r = x_l + 1.0;
1144:       for (i = 0; i < m; i++) {
1145:         y_l   = m - i - 1.0;
1146:         y_r   = y_l + 1.0;
1147:         color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(v[j*m+i]));
1148:         PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1149:       }
1150:     }
1151:   }
1152:   return(0);
1153: }

1157: PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
1158: {
1159:   PetscDraw      draw;
1160:   PetscBool      isnull;
1161:   PetscReal      xr,yr,xl,yl,h,w;

1165:   PetscViewerDrawGetDraw(viewer,0,&draw);
1166:   PetscDrawIsNull(draw,&isnull);
1167:   if (isnull) return(0);

1169:   PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
1170:   xr   = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
1171:   xr  += w;    yr += h;  xl = -w;     yl = -h;
1172:   PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
1173:   PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
1174:   PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
1175:   return(0);
1176: }

1180: PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
1181: {
1183:   PetscBool      iascii,isbinary,isdraw;

1186:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1187:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1188:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);

1190:   if (iascii) {
1191:     MatView_SeqDense_ASCII(A,viewer);
1192:   } else if (isbinary) {
1193:     MatView_SeqDense_Binary(A,viewer);
1194:   } else if (isdraw) {
1195:     MatView_SeqDense_Draw(A,viewer);
1196:   }
1197:   return(0);
1198: }

1202: PetscErrorCode MatDestroy_SeqDense(Mat mat)
1203: {
1204:   Mat_SeqDense   *l = (Mat_SeqDense*)mat->data;

1208: #if defined(PETSC_USE_LOG)
1209:   PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n);
1210: #endif
1211:   PetscFree(l->pivots);
1212:   if (!l->user_alloc) {PetscFree(l->v);}
1213:   PetscFree(mat->data);

1215:   PetscObjectChangeTypeName((PetscObject)mat,0);
1216:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);
1217:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);
1218:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);
1219: #if defined(PETSC_HAVE_ELEMENTAL)
1220:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);
1221: #endif
1222:   PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);
1223:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);
1224:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);
1225:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);
1226:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);
1227:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);
1228:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);
1229:   return(0);
1230: }

1234: PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
1235: {
1236:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1238:   PetscInt       k,j,m,n,M;
1239:   PetscScalar    *v,tmp;

1242:   v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n;
1243:   if (reuse == MAT_REUSE_MATRIX && *matout == A) { /* in place transpose */
1244:     if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can not transpose non-square matrix in place");
1245:     else {
1246:       for (j=0; j<m; j++) {
1247:         for (k=0; k<j; k++) {
1248:           tmp        = v[j + k*M];
1249:           v[j + k*M] = v[k + j*M];
1250:           v[k + j*M] = tmp;
1251:         }
1252:       }
1253:     }
1254:   } else { /* out-of-place transpose */
1255:     Mat          tmat;
1256:     Mat_SeqDense *tmatd;
1257:     PetscScalar  *v2;
1258:     PetscInt     M2;

1260:     if (reuse == MAT_INITIAL_MATRIX) {
1261:       MatCreate(PetscObjectComm((PetscObject)A),&tmat);
1262:       MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);
1263:       MatSetType(tmat,((PetscObject)A)->type_name);
1264:       MatSeqDenseSetPreallocation(tmat,NULL);
1265:     } else {
1266:       tmat = *matout;
1267:     }
1268:     tmatd = (Mat_SeqDense*)tmat->data;
1269:     v = mat->v; v2 = tmatd->v; M2 = tmatd->lda;
1270:     for (j=0; j<n; j++) {
1271:       for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M];
1272:     }
1273:     MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
1274:     MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
1275:     *matout = tmat;
1276:   }
1277:   return(0);
1278: }

1282: PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool  *flg)
1283: {
1284:   Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data;
1285:   Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data;
1286:   PetscInt     i,j;
1287:   PetscScalar  *v1,*v2;

1290:   if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; return(0);}
1291:   if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; return(0);}
1292:   for (i=0; i<A1->rmap->n; i++) {
1293:     v1 = mat1->v+i; v2 = mat2->v+i;
1294:     for (j=0; j<A1->cmap->n; j++) {
1295:       if (*v1 != *v2) {*flg = PETSC_FALSE; return(0);}
1296:       v1 += mat1->lda; v2 += mat2->lda;
1297:     }
1298:   }
1299:   *flg = PETSC_TRUE;
1300:   return(0);
1301: }

1305: PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
1306: {
1307:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1309:   PetscInt       i,n,len;
1310:   PetscScalar    *x,zero = 0.0;

1313:   VecSet(v,zero);
1314:   VecGetSize(v,&n);
1315:   VecGetArray(v,&x);
1316:   len  = PetscMin(A->rmap->n,A->cmap->n);
1317:   if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
1318:   for (i=0; i<len; i++) {
1319:     x[i] = mat->v[i*mat->lda + i];
1320:   }
1321:   VecRestoreArray(v,&x);
1322:   return(0);
1323: }

1327: PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
1328: {
1329:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1330:   const PetscScalar *l,*r;
1331:   PetscScalar       x,*v;
1332:   PetscErrorCode    ierr;
1333:   PetscInt          i,j,m = A->rmap->n,n = A->cmap->n;

1336:   if (ll) {
1337:     VecGetSize(ll,&m);
1338:     VecGetArrayRead(ll,&l);
1339:     if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size");
1340:     for (i=0; i<m; i++) {
1341:       x = l[i];
1342:       v = mat->v + i;
1343:       for (j=0; j<n; j++) { (*v) *= x; v+= m;}
1344:     }
1345:     VecRestoreArrayRead(ll,&l);
1346:     PetscLogFlops(1.0*n*m);
1347:   }
1348:   if (rr) {
1349:     VecGetSize(rr,&n);
1350:     VecGetArrayRead(rr,&r);
1351:     if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size");
1352:     for (i=0; i<n; i++) {
1353:       x = r[i];
1354:       v = mat->v + i*m;
1355:       for (j=0; j<m; j++) (*v++) *= x;
1356:     }
1357:     VecRestoreArrayRead(rr,&r);
1358:     PetscLogFlops(1.0*n*m);
1359:   }
1360:   return(0);
1361: }

1365: PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1366: {
1367:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1368:   PetscScalar    *v   = mat->v;
1369:   PetscReal      sum  = 0.0;
1370:   PetscInt       lda  =mat->lda,m=A->rmap->n,i,j;

1374:   if (type == NORM_FROBENIUS) {
1375:     if (lda>m) {
1376:       for (j=0; j<A->cmap->n; j++) {
1377:         v = mat->v+j*lda;
1378:         for (i=0; i<m; i++) {
1379:           sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1380:         }
1381:       }
1382:     } else {
1383:       for (i=0; i<A->cmap->n*A->rmap->n; i++) {
1384:         sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1385:       }
1386:     }
1387:     *nrm = PetscSqrtReal(sum);
1388:     PetscLogFlops(2.0*A->cmap->n*A->rmap->n);
1389:   } else if (type == NORM_1) {
1390:     *nrm = 0.0;
1391:     for (j=0; j<A->cmap->n; j++) {
1392:       v   = mat->v + j*mat->lda;
1393:       sum = 0.0;
1394:       for (i=0; i<A->rmap->n; i++) {
1395:         sum += PetscAbsScalar(*v);  v++;
1396:       }
1397:       if (sum > *nrm) *nrm = sum;
1398:     }
1399:     PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1400:   } else if (type == NORM_INFINITY) {
1401:     *nrm = 0.0;
1402:     for (j=0; j<A->rmap->n; j++) {
1403:       v   = mat->v + j;
1404:       sum = 0.0;
1405:       for (i=0; i<A->cmap->n; i++) {
1406:         sum += PetscAbsScalar(*v); v += mat->lda;
1407:       }
1408:       if (sum > *nrm) *nrm = sum;
1409:     }
1410:     PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1411:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
1412:   return(0);
1413: }

1417: PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
1418: {
1419:   Mat_SeqDense   *aij = (Mat_SeqDense*)A->data;

1423:   switch (op) {
1424:   case MAT_ROW_ORIENTED:
1425:     aij->roworiented = flg;
1426:     break;
1427:   case MAT_NEW_NONZERO_LOCATIONS:
1428:   case MAT_NEW_NONZERO_LOCATION_ERR:
1429:   case MAT_NEW_NONZERO_ALLOCATION_ERR:
1430:   case MAT_NEW_DIAGONALS:
1431:   case MAT_KEEP_NONZERO_PATTERN:
1432:   case MAT_IGNORE_OFF_PROC_ENTRIES:
1433:   case MAT_USE_HASH_TABLE:
1434:   case MAT_IGNORE_LOWER_TRIANGULAR:
1435:     PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
1436:     break;
1437:   case MAT_SPD:
1438:   case MAT_SYMMETRIC:
1439:   case MAT_STRUCTURALLY_SYMMETRIC:
1440:   case MAT_HERMITIAN:
1441:   case MAT_SYMMETRY_ETERNAL:
1442:     /* These options are handled directly by MatSetOption() */
1443:     break;
1444:   default:
1445:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]);
1446:   }
1447:   return(0);
1448: }

1452: PetscErrorCode MatZeroEntries_SeqDense(Mat A)
1453: {
1454:   Mat_SeqDense   *l = (Mat_SeqDense*)A->data;
1456:   PetscInt       lda=l->lda,m=A->rmap->n,j;

1459:   if (lda>m) {
1460:     for (j=0; j<A->cmap->n; j++) {
1461:       PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));
1462:     }
1463:   } else {
1464:     PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1465:   }
1466:   return(0);
1467: }

1471: PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1472: {
1473:   PetscErrorCode    ierr;
1474:   Mat_SeqDense      *l = (Mat_SeqDense*)A->data;
1475:   PetscInt          m  = l->lda, n = A->cmap->n, i,j;
1476:   PetscScalar       *slot,*bb;
1477:   const PetscScalar *xx;

1480: #if defined(PETSC_USE_DEBUG)
1481:   for (i=0; i<N; i++) {
1482:     if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
1483:     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);
1484:   }
1485: #endif

1487:   /* fix right hand side if needed */
1488:   if (x && b) {
1489:     VecGetArrayRead(x,&xx);
1490:     VecGetArray(b,&bb);
1491:     for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
1492:     VecRestoreArrayRead(x,&xx);
1493:     VecRestoreArray(b,&bb);
1494:   }

1496:   for (i=0; i<N; i++) {
1497:     slot = l->v + rows[i];
1498:     for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
1499:   }
1500:   if (diag != 0.0) {
1501:     if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
1502:     for (i=0; i<N; i++) {
1503:       slot  = l->v + (m+1)*rows[i];
1504:       *slot = diag;
1505:     }
1506:   }
1507:   return(0);
1508: }

1512: PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[])
1513: {
1514:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;

1517:   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");
1518:   *array = mat->v;
1519:   return(0);
1520: }

1524: PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[])
1525: {
1527:   *array = 0; /* user cannot accidently use the array later */
1528:   return(0);
1529: }

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

1536:    Not Collective

1538:    Input Parameter:
1539: .  mat - a MATSEQDENSE or MATMPIDENSE matrix

1541:    Output Parameter:
1542: .   array - pointer to the data

1544:    Level: intermediate

1546: .seealso: MatDenseRestoreArray()
1547: @*/
1548: PetscErrorCode  MatDenseGetArray(Mat A,PetscScalar **array)
1549: {

1553:   PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));
1554:   return(0);
1555: }

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

1562:    Not Collective

1564:    Input Parameters:
1565: .  mat - a MATSEQDENSE or MATMPIDENSE matrix
1566: .  array - pointer to the data

1568:    Level: intermediate

1570: .seealso: MatDenseGetArray()
1571: @*/
1572: PetscErrorCode  MatDenseRestoreArray(Mat A,PetscScalar **array)
1573: {

1577:   PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));
1578:   return(0);
1579: }

1583: static PetscErrorCode MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B)
1584: {
1585:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1587:   PetscInt       i,j,nrows,ncols;
1588:   const PetscInt *irow,*icol;
1589:   PetscScalar    *av,*bv,*v = mat->v;
1590:   Mat            newmat;

1593:   ISGetIndices(isrow,&irow);
1594:   ISGetIndices(iscol,&icol);
1595:   ISGetLocalSize(isrow,&nrows);
1596:   ISGetLocalSize(iscol,&ncols);

1598:   /* Check submatrixcall */
1599:   if (scall == MAT_REUSE_MATRIX) {
1600:     PetscInt n_cols,n_rows;
1601:     MatGetSize(*B,&n_rows,&n_cols);
1602:     if (n_rows != nrows || n_cols != ncols) {
1603:       /* resize the result matrix to match number of requested rows/columns */
1604:       MatSetSizes(*B,nrows,ncols,nrows,ncols);
1605:     }
1606:     newmat = *B;
1607:   } else {
1608:     /* Create and fill new matrix */
1609:     MatCreate(PetscObjectComm((PetscObject)A),&newmat);
1610:     MatSetSizes(newmat,nrows,ncols,nrows,ncols);
1611:     MatSetType(newmat,((PetscObject)A)->type_name);
1612:     MatSeqDenseSetPreallocation(newmat,NULL);
1613:   }

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

1618:   for (i=0; i<ncols; i++) {
1619:     av = v + mat->lda*icol[i];
1620:     for (j=0; j<nrows; j++) *bv++ = av[irow[j]];
1621:   }

1623:   /* Assemble the matrices so that the correct flags are set */
1624:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1625:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);

1627:   /* Free work space */
1628:   ISRestoreIndices(isrow,&irow);
1629:   ISRestoreIndices(iscol,&icol);
1630:   *B   = newmat;
1631:   return(0);
1632: }

1636: PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[])
1637: {
1639:   PetscInt       i;

1642:   if (scall == MAT_INITIAL_MATRIX) {
1643:     PetscMalloc1(n+1,B);
1644:   }

1646:   for (i=0; i<n; i++) {
1647:     MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
1648:   }
1649:   return(0);
1650: }

1654: PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
1655: {
1657:   return(0);
1658: }

1662: PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
1663: {
1665:   return(0);
1666: }

1670: PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
1671: {
1672:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data;
1674:   PetscInt       lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j;

1677:   /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
1678:   if (A->ops->copy != B->ops->copy) {
1679:     MatCopy_Basic(A,B,str);
1680:     return(0);
1681:   }
1682:   if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)");
1683:   if (lda1>m || lda2>m) {
1684:     for (j=0; j<n; j++) {
1685:       PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));
1686:     }
1687:   } else {
1688:     PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1689:   }
1690:   return(0);
1691: }

1695: PetscErrorCode MatSetUp_SeqDense(Mat A)
1696: {

1700:    MatSeqDenseSetPreallocation(A,0);
1701:   return(0);
1702: }

1706: static PetscErrorCode MatConjugate_SeqDense(Mat A)
1707: {
1708:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1709:   PetscInt     i,nz = A->rmap->n*A->cmap->n;
1710:   PetscScalar  *aa = a->v;

1713:   for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
1714:   return(0);
1715: }

1719: static PetscErrorCode MatRealPart_SeqDense(Mat A)
1720: {
1721:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1722:   PetscInt     i,nz = A->rmap->n*A->cmap->n;
1723:   PetscScalar  *aa = a->v;

1726:   for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
1727:   return(0);
1728: }

1732: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
1733: {
1734:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1735:   PetscInt     i,nz = A->rmap->n*A->cmap->n;
1736:   PetscScalar  *aa = a->v;

1739:   for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
1740:   return(0);
1741: }

1743: /* ----------------------------------------------------------------*/
1746: PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1747: {

1751:   if (scall == MAT_INITIAL_MATRIX) {
1752:     PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
1753:     MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
1754:     PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
1755:   }
1756:   PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
1757:   MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);
1758:   PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
1759:   return(0);
1760: }

1764: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1765: {
1767:   PetscInt       m=A->rmap->n,n=B->cmap->n;
1768:   Mat            Cmat;

1771:   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);
1772:   MatCreate(PETSC_COMM_SELF,&Cmat);
1773:   MatSetSizes(Cmat,m,n,m,n);
1774:   MatSetType(Cmat,MATSEQDENSE);
1775:   MatSeqDenseSetPreallocation(Cmat,NULL);

1777:   *C = Cmat;
1778:   return(0);
1779: }

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

1793:   PetscBLASIntCast(A->rmap->n,&m);
1794:   PetscBLASIntCast(B->cmap->n,&n);
1795:   PetscBLASIntCast(A->cmap->n,&k);
1796:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
1797:   return(0);
1798: }

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

1807:   if (scall == MAT_INITIAL_MATRIX) {
1808:     PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);
1809:     MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
1810:     PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);
1811:   }
1812:   PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);
1813:   MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);
1814:   PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);
1815:   return(0);
1816: }

1820: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1821: {
1823:   PetscInt       m=A->cmap->n,n=B->cmap->n;
1824:   Mat            Cmat;

1827:   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);
1828:   MatCreate(PETSC_COMM_SELF,&Cmat);
1829:   MatSetSizes(Cmat,m,n,m,n);
1830:   MatSetType(Cmat,MATSEQDENSE);
1831:   MatSeqDenseSetPreallocation(Cmat,NULL);

1833:   Cmat->assembled = PETSC_TRUE;

1835:   *C = Cmat;
1836:   return(0);
1837: }

1841: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
1842: {
1843:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
1844:   Mat_SeqDense   *b = (Mat_SeqDense*)B->data;
1845:   Mat_SeqDense   *c = (Mat_SeqDense*)C->data;
1846:   PetscBLASInt   m,n,k;
1847:   PetscScalar    _DOne=1.0,_DZero=0.0;

1851:   PetscBLASIntCast(A->cmap->n,&m);
1852:   PetscBLASIntCast(B->cmap->n,&n);
1853:   PetscBLASIntCast(A->rmap->n,&k);
1854:   /*
1855:      Note the m and n arguments below are the number rows and columns of A', not A!
1856:   */
1857:   PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
1858:   return(0);
1859: }

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

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

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

1890: PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])
1891: {
1892:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
1894:   PetscInt       i,j,m = A->rmap->n,n = A->cmap->n,p;
1895:   PetscScalar    *x;
1896:   PetscReal      atmp;
1897:   MatScalar      *aa = a->v;

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

1902:   VecSet(v,0.0);
1903:   VecGetArray(v,&x);
1904:   VecGetLocalSize(v,&p);
1905:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
1906:   for (i=0; i<m; i++) {
1907:     x[i] = PetscAbsScalar(aa[i]);
1908:     for (j=1; j<n; j++) {
1909:       atmp = PetscAbsScalar(aa[i+m*j]);
1910:       if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;}
1911:     }
1912:   }
1913:   VecRestoreArray(v,&x);
1914:   return(0);
1915: }

1919: PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])
1920: {
1921:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
1923:   PetscInt       i,j,m = A->rmap->n,n = A->cmap->n,p;
1924:   PetscScalar    *x;
1925:   MatScalar      *aa = a->v;

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

1930:   VecSet(v,0.0);
1931:   VecGetArray(v,&x);
1932:   VecGetLocalSize(v,&p);
1933:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
1934:   for (i=0; i<m; i++) {
1935:     x[i] = aa[i]; if (idx) idx[i] = 0;
1936:     for (j=1; j<n; j++) {
1937:       if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
1938:     }
1939:   }
1940:   VecRestoreArray(v,&x);
1941:   return(0);
1942: }

1946: PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
1947: {
1948:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
1950:   PetscScalar    *x;

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

1955:   VecGetArray(v,&x);
1956:   PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));
1957:   VecRestoreArray(v,&x);
1958:   return(0);
1959: }


1964: PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms)
1965: {
1967:   PetscInt       i,j,m,n;
1968:   PetscScalar    *a;

1971:   MatGetSize(A,&m,&n);
1972:   PetscMemzero(norms,n*sizeof(PetscReal));
1973:   MatDenseGetArray(A,&a);
1974:   if (type == NORM_2) {
1975:     for (i=0; i<n; i++) {
1976:       for (j=0; j<m; j++) {
1977:         norms[i] += PetscAbsScalar(a[j]*a[j]);
1978:       }
1979:       a += m;
1980:     }
1981:   } else if (type == NORM_1) {
1982:     for (i=0; i<n; i++) {
1983:       for (j=0; j<m; j++) {
1984:         norms[i] += PetscAbsScalar(a[j]);
1985:       }
1986:       a += m;
1987:     }
1988:   } else if (type == NORM_INFINITY) {
1989:     for (i=0; i<n; i++) {
1990:       for (j=0; j<m; j++) {
1991:         norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]);
1992:       }
1993:       a += m;
1994:     }
1995:   } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType");
1996:   MatDenseRestoreArray(A,&a);
1997:   if (type == NORM_2) {
1998:     for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]);
1999:   }
2000:   return(0);
2001: }

2005: static PetscErrorCode  MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
2006: {
2008:   PetscScalar    *a;
2009:   PetscInt       m,n,i;

2012:   MatGetSize(x,&m,&n);
2013:   MatDenseGetArray(x,&a);
2014:   for (i=0; i<m*n; i++) {
2015:     PetscRandomGetValue(rctx,a+i);
2016:   }
2017:   MatDenseRestoreArray(x,&a);
2018:   return(0);
2019: }


2022: /* -------------------------------------------------------------------*/
2023: static struct _MatOps MatOps_Values = { MatSetValues_SeqDense,
2024:                                         MatGetRow_SeqDense,
2025:                                         MatRestoreRow_SeqDense,
2026:                                         MatMult_SeqDense,
2027:                                 /*  4*/ MatMultAdd_SeqDense,
2028:                                         MatMultTranspose_SeqDense,
2029:                                         MatMultTransposeAdd_SeqDense,
2030:                                         0,
2031:                                         0,
2032:                                         0,
2033:                                 /* 10*/ 0,
2034:                                         MatLUFactor_SeqDense,
2035:                                         MatCholeskyFactor_SeqDense,
2036:                                         MatSOR_SeqDense,
2037:                                         MatTranspose_SeqDense,
2038:                                 /* 15*/ MatGetInfo_SeqDense,
2039:                                         MatEqual_SeqDense,
2040:                                         MatGetDiagonal_SeqDense,
2041:                                         MatDiagonalScale_SeqDense,
2042:                                         MatNorm_SeqDense,
2043:                                 /* 20*/ MatAssemblyBegin_SeqDense,
2044:                                         MatAssemblyEnd_SeqDense,
2045:                                         MatSetOption_SeqDense,
2046:                                         MatZeroEntries_SeqDense,
2047:                                 /* 24*/ MatZeroRows_SeqDense,
2048:                                         0,
2049:                                         0,
2050:                                         0,
2051:                                         0,
2052:                                 /* 29*/ MatSetUp_SeqDense,
2053:                                         0,
2054:                                         0,
2055:                                         0,
2056:                                         0,
2057:                                 /* 34*/ MatDuplicate_SeqDense,
2058:                                         0,
2059:                                         0,
2060:                                         0,
2061:                                         0,
2062:                                 /* 39*/ MatAXPY_SeqDense,
2063:                                         MatGetSubMatrices_SeqDense,
2064:                                         0,
2065:                                         MatGetValues_SeqDense,
2066:                                         MatCopy_SeqDense,
2067:                                 /* 44*/ MatGetRowMax_SeqDense,
2068:                                         MatScale_SeqDense,
2069:                                         MatShift_Basic,
2070:                                         0,
2071:                                         0,
2072:                                 /* 49*/ MatSetRandom_SeqDense,
2073:                                         0,
2074:                                         0,
2075:                                         0,
2076:                                         0,
2077:                                 /* 54*/ 0,
2078:                                         0,
2079:                                         0,
2080:                                         0,
2081:                                         0,
2082:                                 /* 59*/ 0,
2083:                                         MatDestroy_SeqDense,
2084:                                         MatView_SeqDense,
2085:                                         0,
2086:                                         0,
2087:                                 /* 64*/ 0,
2088:                                         0,
2089:                                         0,
2090:                                         0,
2091:                                         0,
2092:                                 /* 69*/ MatGetRowMaxAbs_SeqDense,
2093:                                         0,
2094:                                         0,
2095:                                         0,
2096:                                         0,
2097:                                 /* 74*/ 0,
2098:                                         0,
2099:                                         0,
2100:                                         0,
2101:                                         0,
2102:                                 /* 79*/ 0,
2103:                                         0,
2104:                                         0,
2105:                                         0,
2106:                                 /* 83*/ MatLoad_SeqDense,
2107:                                         0,
2108:                                         MatIsHermitian_SeqDense,
2109:                                         0,
2110:                                         0,
2111:                                         0,
2112:                                 /* 89*/ MatMatMult_SeqDense_SeqDense,
2113:                                         MatMatMultSymbolic_SeqDense_SeqDense,
2114:                                         MatMatMultNumeric_SeqDense_SeqDense,
2115:                                         0,
2116:                                         0,
2117:                                 /* 94*/ 0,
2118:                                         0,
2119:                                         0,
2120:                                         0,
2121:                                         0,
2122:                                 /* 99*/ 0,
2123:                                         0,
2124:                                         0,
2125:                                         MatConjugate_SeqDense,
2126:                                         0,
2127:                                 /*104*/ 0,
2128:                                         MatRealPart_SeqDense,
2129:                                         MatImaginaryPart_SeqDense,
2130:                                         0,
2131:                                         0,
2132:                                 /*109*/ MatMatSolve_SeqDense,
2133:                                         0,
2134:                                         MatGetRowMin_SeqDense,
2135:                                         MatGetColumnVector_SeqDense,
2136:                                         0,
2137:                                 /*114*/ 0,
2138:                                         0,
2139:                                         0,
2140:                                         0,
2141:                                         0,
2142:                                 /*119*/ 0,
2143:                                         0,
2144:                                         0,
2145:                                         0,
2146:                                         0,
2147:                                 /*124*/ 0,
2148:                                         MatGetColumnNorms_SeqDense,
2149:                                         0,
2150:                                         0,
2151:                                         0,
2152:                                 /*129*/ 0,
2153:                                         MatTransposeMatMult_SeqDense_SeqDense,
2154:                                         MatTransposeMatMultSymbolic_SeqDense_SeqDense,
2155:                                         MatTransposeMatMultNumeric_SeqDense_SeqDense,
2156:                                         0,
2157:                                 /*134*/ 0,
2158:                                         0,
2159:                                         0,
2160:                                         0,
2161:                                         0,
2162:                                 /*139*/ 0,
2163:                                         0,
2164:                                         0
2165: };

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

2174:    Collective on MPI_Comm

2176:    Input Parameters:
2177: +  comm - MPI communicator, set to PETSC_COMM_SELF
2178: .  m - number of rows
2179: .  n - number of columns
2180: -  data - optional location of matrix data in column major order.  Set data=NULL for PETSc
2181:    to control all matrix memory allocation.

2183:    Output Parameter:
2184: .  A - the matrix

2186:    Notes:
2187:    The data input variable is intended primarily for Fortran programmers
2188:    who wish to allocate their own matrix memory space.  Most users should
2189:    set data=NULL.

2191:    Level: intermediate

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

2195: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
2196: @*/
2197: PetscErrorCode  MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
2198: {

2202:   MatCreate(comm,A);
2203:   MatSetSizes(*A,m,n,m,n);
2204:   MatSetType(*A,MATSEQDENSE);
2205:   MatSeqDenseSetPreallocation(*A,data);
2206:   return(0);
2207: }

2211: /*@C
2212:    MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements

2214:    Collective on MPI_Comm

2216:    Input Parameters:
2217: +  B - the matrix
2218: -  data - the array (or NULL)

2220:    Notes:
2221:    The data input variable is intended primarily for Fortran programmers
2222:    who wish to allocate their own matrix memory space.  Most users should
2223:    need not call this routine.

2225:    Level: intermediate

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

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

2231: @*/
2232: PetscErrorCode  MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
2233: {

2237:   PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
2238:   return(0);
2239: }

2243: PetscErrorCode  MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
2244: {
2245:   Mat_SeqDense   *b;

2249:   B->preallocated = PETSC_TRUE;

2251:   PetscLayoutSetUp(B->rmap);
2252:   PetscLayoutSetUp(B->cmap);

2254:   b       = (Mat_SeqDense*)B->data;
2255:   b->Mmax = B->rmap->n;
2256:   b->Nmax = B->cmap->n;
2257:   if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n;

2259:   if (!data) { /* petsc-allocated storage */
2260:     if (!b->user_alloc) { PetscFree(b->v); }
2261:     PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);
2262:     PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));

2264:     b->user_alloc = PETSC_FALSE;
2265:   } else { /* user-allocated storage */
2266:     if (!b->user_alloc) { PetscFree(b->v); }
2267:     b->v          = data;
2268:     b->user_alloc = PETSC_TRUE;
2269:   }
2270:   B->assembled = PETSC_TRUE;
2271:   return(0);
2272: }

2274: #if defined(PETSC_HAVE_ELEMENTAL)
2277: PETSC_EXTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
2278: {
2279:   Mat            mat_elemental;
2281:   PetscScalar    *array,*v_colwise;
2282:   PetscInt       M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols;

2285:   PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);
2286:   MatDenseGetArray(A,&array);
2287:   /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
2288:   k = 0;
2289:   for (j=0; j<N; j++) {
2290:     cols[j] = j;
2291:     for (i=0; i<M; i++) {
2292:       v_colwise[j*M+i] = array[k++];
2293:     }
2294:   }
2295:   for (i=0; i<M; i++) {
2296:     rows[i] = i;
2297:   }
2298:   MatDenseRestoreArray(A,&array);

2300:   MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
2301:   MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);
2302:   MatSetType(mat_elemental,MATELEMENTAL);
2303:   MatSetUp(mat_elemental);

2305:   /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
2306:   MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);
2307:   MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
2308:   MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);
2309:   PetscFree3(v_colwise,rows,cols);

2311:   if (reuse == MAT_REUSE_MATRIX) {
2312:     MatHeaderReplace(A,mat_elemental);
2313:   } else {
2314:     *newmat = mat_elemental;
2315:   }
2316:   return(0);
2317: }
2318: #endif

2322: /*@C
2323:   MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array

2325:   Input parameter:
2326: + A - the matrix
2327: - lda - the leading dimension

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

2334:   Level: intermediate

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

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

2340: @*/
2341: PetscErrorCode  MatSeqDenseSetLDA(Mat B,PetscInt lda)
2342: {
2343:   Mat_SeqDense *b = (Mat_SeqDense*)B->data;

2346:   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);
2347:   b->lda       = lda;
2348:   b->changelda = PETSC_FALSE;
2349:   b->Mmax      = PetscMax(b->Mmax,lda);
2350:   return(0);
2351: }

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

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

2359:   Level: beginner

2361: .seealso: MatCreateSeqDense()

2363: M*/

2367: PETSC_EXTERN PetscErrorCode MatCreate_SeqDense(Mat B)
2368: {
2369:   Mat_SeqDense   *b;
2371:   PetscMPIInt    size;

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

2377:   PetscNewLog(B,&b);
2378:   PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
2379:   B->data = (void*)b;

2381:   b->pivots      = 0;
2382:   b->roworiented = PETSC_TRUE;
2383:   b->v           = 0;
2384:   b->changelda   = PETSC_FALSE;

2386:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);
2387:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);
2388:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);
2389: #if defined(PETSC_HAVE_ELEMENTAL)
2390:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);
2391: #endif
2392:   PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);
2393:   PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2394:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2395:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2396:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2397:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2398:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2399:   PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
2400:   return(0);
2401: }