Actual source code: dense.c

petsc-master 2014-12-15
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:   PetscObjectStateIncrease((PetscObject)Y);
 74:   PetscLogFlops(PetscMax(2*N-1,0));
 75:   return(0);
 76: }

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

496:   info.fill = 1.0;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

806: /* -----------------------------------------------------------------*/

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

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

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

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

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

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

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

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

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

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

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

978:   PetscViewerBinaryGetDescriptor(viewer,&fd);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1162:   PetscObjectChangeTypeName((PetscObject)mat,0);
1163:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);
1164:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);
1165:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);
1166: #if defined(PETSC_HAVE_ELEMENTAL)
1167:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);
1168: #endif
1169:   PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);
1170:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);
1171:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);
1172:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);
1173:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);
1174:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);
1175:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);
1176:   return(0);
1177: }

1181: PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
1182: {
1183:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1185:   PetscInt       k,j,m,n,M;
1186:   PetscScalar    *v,tmp;

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

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

1228: PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool  *flg)
1229: {
1230:   Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data;
1231:   Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data;
1232:   PetscInt     i,j;
1233:   PetscScalar  *v1,*v2;

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

1251: PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
1252: {
1253:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1255:   PetscInt       i,n,len;
1256:   PetscScalar    *x,zero = 0.0;

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

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

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

1310: PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1311: {
1312:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1313:   PetscScalar    *v   = mat->v;
1314:   PetscReal      sum  = 0.0;
1315:   PetscInt       lda  =mat->lda,m=A->rmap->n,i,j;

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

1362: PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
1363: {
1364:   Mat_SeqDense   *aij = (Mat_SeqDense*)A->data;

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

1397: PetscErrorCode MatZeroEntries_SeqDense(Mat A)
1398: {
1399:   Mat_SeqDense   *l = (Mat_SeqDense*)A->data;
1401:   PetscInt       lda=l->lda,m=A->rmap->n,j;

1404:   if (lda>m) {
1405:     for (j=0; j<A->cmap->n; j++) {
1406:       PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));
1407:     }
1408:   } else {
1409:     PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1410:   }
1411:   return(0);
1412: }

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

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

1432:   /* fix right hand side if needed */
1433:   if (x && b) {
1434:     VecGetArrayRead(x,&xx);
1435:     VecGetArray(b,&bb);
1436:     for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
1437:     VecRestoreArrayRead(x,&xx);
1438:     VecRestoreArray(b,&bb);
1439:   }

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

1457: PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[])
1458: {
1459:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;

1462:   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");
1463:   *array = mat->v;
1464:   return(0);
1465: }

1469: PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[])
1470: {
1472:   *array = 0; /* user cannot accidently use the array later */
1473:   return(0);
1474: }

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

1481:    Not Collective

1483:    Input Parameter:
1484: .  mat - a MATSEQDENSE or MATMPIDENSE matrix

1486:    Output Parameter:
1487: .   array - pointer to the data

1489:    Level: intermediate

1491: .seealso: MatDenseRestoreArray()
1492: @*/
1493: PetscErrorCode  MatDenseGetArray(Mat A,PetscScalar **array)
1494: {

1498:   PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));
1499:   return(0);
1500: }

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

1507:    Not Collective

1509:    Input Parameters:
1510: .  mat - a MATSEQDENSE or MATMPIDENSE matrix
1511: .  array - pointer to the data

1513:    Level: intermediate

1515: .seealso: MatDenseGetArray()
1516: @*/
1517: PetscErrorCode  MatDenseRestoreArray(Mat A,PetscScalar **array)
1518: {

1522:   PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));
1523:   return(0);
1524: }

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

1538:   ISGetIndices(isrow,&irow);
1539:   ISGetIndices(iscol,&icol);
1540:   ISGetLocalSize(isrow,&nrows);
1541:   ISGetLocalSize(iscol,&ncols);

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

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

1563:   for (i=0; i<ncols; i++) {
1564:     av = v + mat->lda*icol[i];
1565:     for (j=0; j<nrows; j++) *bv++ = av[irow[j]];
1566:   }

1568:   /* Assemble the matrices so that the correct flags are set */
1569:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1570:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);

1572:   /* Free work space */
1573:   ISRestoreIndices(isrow,&irow);
1574:   ISRestoreIndices(iscol,&icol);
1575:   *B   = newmat;
1576:   return(0);
1577: }

1581: PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[])
1582: {
1584:   PetscInt       i;

1587:   if (scall == MAT_INITIAL_MATRIX) {
1588:     PetscMalloc1(n+1,B);
1589:   }

1591:   for (i=0; i<n; i++) {
1592:     MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
1593:   }
1594:   return(0);
1595: }

1599: PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
1600: {
1602:   return(0);
1603: }

1607: PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
1608: {
1610:   return(0);
1611: }

1615: PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
1616: {
1617:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data;
1619:   PetscInt       lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j;

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

1640: PetscErrorCode MatSetUp_SeqDense(Mat A)
1641: {

1645:    MatSeqDenseSetPreallocation(A,0);
1646:   return(0);
1647: }

1651: static PetscErrorCode MatConjugate_SeqDense(Mat A)
1652: {
1653:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1654:   PetscInt     i,nz = A->rmap->n*A->cmap->n;
1655:   PetscScalar  *aa = a->v;

1658:   for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
1659:   return(0);
1660: }

1664: static PetscErrorCode MatRealPart_SeqDense(Mat A)
1665: {
1666:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1667:   PetscInt     i,nz = A->rmap->n*A->cmap->n;
1668:   PetscScalar  *aa = a->v;

1671:   for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
1672:   return(0);
1673: }

1677: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
1678: {
1679:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1680:   PetscInt     i,nz = A->rmap->n*A->cmap->n;
1681:   PetscScalar  *aa = a->v;

1684:   for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
1685:   return(0);
1686: }

1688: /* ----------------------------------------------------------------*/
1691: PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1692: {

1696:   if (scall == MAT_INITIAL_MATRIX) {
1697:     PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
1698:     MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
1699:     PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
1700:   }
1701:   PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
1702:   MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);
1703:   PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
1704:   return(0);
1705: }

1709: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1710: {
1712:   PetscInt       m=A->rmap->n,n=B->cmap->n;
1713:   Mat            Cmat;

1716:   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);
1717:   MatCreate(PETSC_COMM_SELF,&Cmat);
1718:   MatSetSizes(Cmat,m,n,m,n);
1719:   MatSetType(Cmat,MATSEQDENSE);
1720:   MatSeqDenseSetPreallocation(Cmat,NULL);

1722:   *C = Cmat;
1723:   return(0);
1724: }

1728: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
1729: {
1730:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
1731:   Mat_SeqDense   *b = (Mat_SeqDense*)B->data;
1732:   Mat_SeqDense   *c = (Mat_SeqDense*)C->data;
1733:   PetscBLASInt   m,n,k;
1734:   PetscScalar    _DOne=1.0,_DZero=0.0;

1738:   PetscBLASIntCast(A->rmap->n,&m);
1739:   PetscBLASIntCast(B->cmap->n,&n);
1740:   PetscBLASIntCast(A->cmap->n,&k);
1741:   PetscStackCallBLAS("BLASgemv",BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
1742:   return(0);
1743: }

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

1752:   if (scall == MAT_INITIAL_MATRIX) {
1753:     PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);
1754:     MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
1755:     PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);
1756:   }
1757:   PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);
1758:   MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);
1759:   PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);
1760:   return(0);
1761: }

1765: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1766: {
1768:   PetscInt       m=A->cmap->n,n=B->cmap->n;
1769:   Mat            Cmat;

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

1778:   Cmat->assembled = PETSC_TRUE;

1780:   *C = Cmat;
1781:   return(0);
1782: }

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

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

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

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

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

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

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

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

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

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

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

1891: PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
1892: {
1893:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
1895:   PetscScalar    *x;

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

1900:   VecGetArray(v,&x);
1901:   PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));
1902:   VecRestoreArray(v,&x);
1903:   return(0);
1904: }


1909: PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms)
1910: {
1912:   PetscInt       i,j,m,n;
1913:   PetscScalar    *a;

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

1950: static PetscErrorCode  MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
1951: {
1953:   PetscScalar    *a;
1954:   PetscInt       m,n,i;

1957:   MatGetSize(x,&m,&n);
1958:   MatDenseGetArray(x,&a);
1959:   for (i=0; i<m*n; i++) {
1960:     PetscRandomGetValue(rctx,a+i);
1961:   }
1962:   MatDenseRestoreArray(x,&a);
1963:   return(0);
1964: }


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

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

2119:    Collective on MPI_Comm

2121:    Input Parameters:
2122: +  comm - MPI communicator, set to PETSC_COMM_SELF
2123: .  m - number of rows
2124: .  n - number of columns
2125: -  data - optional location of matrix data in column major order.  Set data=NULL for PETSc
2126:    to control all matrix memory allocation.

2128:    Output Parameter:
2129: .  A - the matrix

2131:    Notes:
2132:    The data input variable is intended primarily for Fortran programmers
2133:    who wish to allocate their own matrix memory space.  Most users should
2134:    set data=NULL.

2136:    Level: intermediate

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

2140: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
2141: @*/
2142: PetscErrorCode  MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
2143: {

2147:   MatCreate(comm,A);
2148:   MatSetSizes(*A,m,n,m,n);
2149:   MatSetType(*A,MATSEQDENSE);
2150:   MatSeqDenseSetPreallocation(*A,data);
2151:   return(0);
2152: }

2156: /*@C
2157:    MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements

2159:    Collective on MPI_Comm

2161:    Input Parameters:
2162: +  B - the matrix
2163: -  data - the array (or NULL)

2165:    Notes:
2166:    The data input variable is intended primarily for Fortran programmers
2167:    who wish to allocate their own matrix memory space.  Most users should
2168:    need not call this routine.

2170:    Level: intermediate

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

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

2176: @*/
2177: PetscErrorCode  MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
2178: {

2182:   PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
2183:   return(0);
2184: }

2188: PetscErrorCode  MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
2189: {
2190:   Mat_SeqDense   *b;

2194:   B->preallocated = PETSC_TRUE;

2196:   PetscLayoutSetUp(B->rmap);
2197:   PetscLayoutSetUp(B->cmap);

2199:   b       = (Mat_SeqDense*)B->data;
2200:   b->Mmax = B->rmap->n;
2201:   b->Nmax = B->cmap->n;
2202:   if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n;

2204:   if (!data) { /* petsc-allocated storage */
2205:     if (!b->user_alloc) { PetscFree(b->v); }
2206:     PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);
2207:     PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));

2209:     b->user_alloc = PETSC_FALSE;
2210:   } else { /* user-allocated storage */
2211:     if (!b->user_alloc) { PetscFree(b->v); }
2212:     b->v          = data;
2213:     b->user_alloc = PETSC_TRUE;
2214:   }
2215:   B->assembled = PETSC_TRUE;
2216:   return(0);
2217: }

2221: PETSC_EXTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
2222: {
2223:   Mat            mat_elemental;
2225:   PetscScalar    *array,*v_colwise;
2226:   PetscInt       M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols;

2229:   PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);
2230:   MatDenseGetArray(A,&array);
2231:   /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
2232:   k = 0;
2233:   for (j=0; j<N; j++) {
2234:     cols[j] = j;
2235:     for (i=0; i<M; i++) {
2236:       v_colwise[j*M+i] = array[k++];
2237:     }
2238:   }
2239:   for (i=0; i<M; i++) {
2240:     rows[i] = i;
2241:   }
2242:   MatDenseRestoreArray(A,&array);

2244:   MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
2245:   MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);
2246:   MatSetType(mat_elemental,MATELEMENTAL);
2247:   MatSetUp(mat_elemental);

2249:   /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
2250:   MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);
2251:   MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
2252:   MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);
2253:   PetscFree3(v_colwise,rows,cols);

2255:   if (reuse == MAT_REUSE_MATRIX) {
2256:     MatHeaderReplace(A,mat_elemental);
2257:   } else {
2258:     *newmat = mat_elemental;
2259:   }
2260:   return(0);
2261: }

2265: /*@C
2266:   MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array

2268:   Input parameter:
2269: + A - the matrix
2270: - lda - the leading dimension

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

2277:   Level: intermediate

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

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

2283: @*/
2284: PetscErrorCode  MatSeqDenseSetLDA(Mat B,PetscInt lda)
2285: {
2286:   Mat_SeqDense *b = (Mat_SeqDense*)B->data;

2289:   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);
2290:   b->lda       = lda;
2291:   b->changelda = PETSC_FALSE;
2292:   b->Mmax      = PetscMax(b->Mmax,lda);
2293:   return(0);
2294: }

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

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

2302:   Level: beginner

2304: .seealso: MatCreateSeqDense()

2306: M*/

2310: PETSC_EXTERN PetscErrorCode MatCreate_SeqDense(Mat B)
2311: {
2312:   Mat_SeqDense   *b;
2314:   PetscMPIInt    size;

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

2320:   PetscNewLog(B,&b);
2321:   PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
2322:   B->data = (void*)b;

2324:   b->pivots      = 0;
2325:   b->roworiented = PETSC_TRUE;
2326:   b->v           = 0;
2327:   b->changelda   = PETSC_FALSE;

2329:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);
2330:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);
2331:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);
2332: #if defined(PETSC_HAVE_ELEMENTAL)
2333:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);
2334: #endif
2335:   PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);
2336:   PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2337:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2338:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2339:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2340:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2341:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2342:   PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
2343:   return(0);
2344: }