Actual source code: dense.c

petsc-3.7.3 2016-08-01
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: PetscErrorCode MatZeroRowsColumns_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
 14: {
 15:   PetscErrorCode    ierr;
 16:   Mat_SeqDense      *l = (Mat_SeqDense*)A->data;
 17:   PetscInt          m  = l->lda, n = A->cmap->n,r = A->rmap->n, i,j;
 18:   PetscScalar       *slot,*bb;
 19:   const PetscScalar *xx;

 22: #if defined(PETSC_USE_DEBUG)
 23:   for (i=0; i<N; i++) {
 24:     if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
 25:     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);
 26:     if (rows[i] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Col %D requested to be zeroed greater than or equal number of cols %D",rows[i],A->cmap->n);
 27:   }
 28: #endif

 30:   /* fix right hand side if needed */
 31:   if (x && b) {
 32:     VecGetArrayRead(x,&xx);
 33:     VecGetArray(b,&bb);
 34:     for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
 35:     VecRestoreArrayRead(x,&xx);
 36:     VecRestoreArray(b,&bb);
 37:   }

 39:   for (i=0; i<N; i++) {
 40:     slot = l->v + rows[i]*m;
 41:     PetscMemzero(slot,r*sizeof(PetscScalar));
 42:   }
 43:   for (i=0; i<N; i++) {
 44:     slot = l->v + rows[i];
 45:     for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
 46:   }
 47:   if (diag != 0.0) {
 48:     if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
 49:     for (i=0; i<N; i++) {
 50:       slot  = l->v + (m+1)*rows[i];
 51:       *slot = diag;
 52:     }
 53:   }
 54:   return(0);
 55: }

 59: PetscErrorCode MatPtAPNumeric_SeqDense_SeqDense(Mat A,Mat P,Mat C)
 60: {
 61:   Mat_SeqDense   *c = (Mat_SeqDense*)(C->data);

 65:   MatMatMultNumeric_SeqDense_SeqDense(A,P,c->ptapwork);
 66:   MatTransposeMatMultNumeric_SeqDense_SeqDense(P,c->ptapwork,C);
 67:   return(0);
 68: }

 72: PetscErrorCode MatPtAPSymbolic_SeqDense_SeqDense(Mat A,Mat P,PetscReal fill,Mat *C)
 73: {
 74:   Mat_SeqDense   *c;

 78:   MatCreateSeqDense(PetscObjectComm((PetscObject)A),P->cmap->N,P->cmap->N,NULL,C);
 79:   c = (Mat_SeqDense*)((*C)->data);
 80:   MatCreateSeqDense(PetscObjectComm((PetscObject)A),A->rmap->N,P->cmap->N,NULL,&c->ptapwork);
 81:   return(0);
 82: }

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

 91:   if (reuse == MAT_INITIAL_MATRIX) {
 92:     MatPtAPSymbolic_SeqDense_SeqDense(A,P,fill,C);
 93:   }
 94:   PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);
 95:   (*(*C)->ops->ptapnumeric)(A,P,*C);
 96:   PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);
 97:   return(0);
 98: }

102: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
103: {
104:   Mat            B;
105:   Mat_SeqAIJ     *a = (Mat_SeqAIJ*)A->data;
106:   Mat_SeqDense   *b;
108:   PetscInt       *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->N,i;
109:   MatScalar      *av=a->a;

112:   MatCreate(PetscObjectComm((PetscObject)A),&B);
113:   MatSetSizes(B,m,n,m,n);
114:   MatSetType(B,MATSEQDENSE);
115:   MatSeqDenseSetPreallocation(B,NULL);

117:   b  = (Mat_SeqDense*)(B->data);
118:   for (i=0; i<m; i++) {
119:     PetscInt j;
120:     for (j=0;j<ai[1]-ai[0];j++) {
121:       b->v[*aj*m+i] = *av;
122:       aj++;
123:       av++;
124:     }
125:     ai++;
126:   }
127:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
128:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

130:   if (reuse == MAT_INPLACE_MATRIX) {
131:     MatHeaderReplace(A,&B);
132:   } else {
133:     *newmat = B;
134:   }
135:   return(0);
136: }

140: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
141: {
142:   Mat            B;
143:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
145:   PetscInt       i, j;
146:   PetscInt       *rows, *nnz;
147:   MatScalar      *aa = a->v, *vals;

150:   MatCreate(PetscObjectComm((PetscObject)A),&B);
151:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
152:   MatSetType(B,MATSEQAIJ);
153:   PetscCalloc3(A->rmap->n,&rows,A->rmap->n,&nnz,A->rmap->n,&vals);
154:   for (j=0; j<A->cmap->n; j++) {
155:     for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) ++nnz[i];
156:     aa += a->lda;
157:   }
158:   MatSeqAIJSetPreallocation(B,PETSC_DETERMINE,nnz);
159:   aa = a->v;
160:   for (j=0; j<A->cmap->n; j++) {
161:     PetscInt numRows = 0;
162:     for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) {rows[numRows] = i; vals[numRows++] = aa[i];}
163:     MatSetValues(B,numRows,rows,1,&j,vals,INSERT_VALUES);
164:     aa  += a->lda;
165:   }
166:   PetscFree3(rows,nnz,vals);
167:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
168:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

170:   if (reuse == MAT_INPLACE_MATRIX) {
171:     MatHeaderReplace(A,&B);
172:   } else {
173:     *newmat = B;
174:   }
175:   return(0);
176: }

180: static PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str)
181: {
182:   Mat_SeqDense   *x     = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data;
183:   PetscScalar    oalpha = alpha;
184:   PetscInt       j;
185:   PetscBLASInt   N,m,ldax,lday,one = 1;

189:   PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);
190:   PetscBLASIntCast(X->rmap->n,&m);
191:   PetscBLASIntCast(x->lda,&ldax);
192:   PetscBLASIntCast(y->lda,&lday);
193:   if (ldax>m || lday>m) {
194:     for (j=0; j<X->cmap->n; j++) {
195:       PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&oalpha,x->v+j*ldax,&one,y->v+j*lday,&one));
196:     }
197:   } else {
198:     PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&oalpha,x->v,&one,y->v,&one));
199:   }
200:   PetscObjectStateIncrease((PetscObject)Y);
201:   PetscLogFlops(PetscMax(2*N-1,0));
202:   return(0);
203: }

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

212:   info->block_size        = 1.0;
213:   info->nz_allocated      = (double)N;
214:   info->nz_used           = (double)N;
215:   info->nz_unneeded       = (double)0;
216:   info->assemblies        = (double)A->num_ass;
217:   info->mallocs           = 0;
218:   info->memory            = ((PetscObject)A)->mem;
219:   info->fill_ratio_given  = 0;
220:   info->fill_ratio_needed = 0;
221:   info->factor_mallocs    = 0;
222:   return(0);
223: }

227: static PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha)
228: {
229:   Mat_SeqDense   *a     = (Mat_SeqDense*)A->data;
230:   PetscScalar    oalpha = alpha;
232:   PetscBLASInt   one = 1,j,nz,lda;

235:   PetscBLASIntCast(a->lda,&lda);
236:   if (lda>A->rmap->n) {
237:     PetscBLASIntCast(A->rmap->n,&nz);
238:     for (j=0; j<A->cmap->n; j++) {
239:       PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v+j*lda,&one));
240:     }
241:   } else {
242:     PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);
243:     PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v,&one));
244:   }
245:   PetscLogFlops(nz);
246:   return(0);
247: }

251: static PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool  *fl)
252: {
253:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
254:   PetscInt     i,j,m = A->rmap->n,N;
255:   PetscScalar  *v = a->v;

258:   *fl = PETSC_FALSE;
259:   if (A->rmap->n != A->cmap->n) return(0);
260:   N = a->lda;

262:   for (i=0; i<m; i++) {
263:     for (j=i+1; j<m; j++) {
264:       if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) return(0);
265:     }
266:   }
267:   *fl = PETSC_TRUE;
268:   return(0);
269: }

273: static PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues)
274: {
275:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data,*l;
277:   PetscInt       lda = (PetscInt)mat->lda,j,m;

280:   PetscLayoutReference(A->rmap,&newi->rmap);
281:   PetscLayoutReference(A->cmap,&newi->cmap);
282:   MatSeqDenseSetPreallocation(newi,NULL);
283:   if (cpvalues == MAT_COPY_VALUES) {
284:     l = (Mat_SeqDense*)newi->data;
285:     if (lda>A->rmap->n) {
286:       m = A->rmap->n;
287:       for (j=0; j<A->cmap->n; j++) {
288:         PetscMemcpy(l->v+j*m,mat->v+j*lda,m*sizeof(PetscScalar));
289:       }
290:     } else {
291:       PetscMemcpy(l->v,mat->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
292:     }
293:   }
294:   newi->assembled = PETSC_TRUE;
295:   return(0);
296: }

300: static PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat)
301: {

305:   MatCreate(PetscObjectComm((PetscObject)A),newmat);
306:   MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
307:   MatSetType(*newmat,((PetscObject)A)->type_name);
308:   MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);
309:   return(0);
310: }


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

317: static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
318: {
319:   MatFactorInfo  info;

323:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
324:   MatLUFactor_SeqDense(fact,0,0,&info);
325:   return(0);
326: }

330: static PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy)
331: {
332:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
333:   PetscErrorCode    ierr;
334:   const PetscScalar *x;
335:   PetscScalar       *y;
336:   PetscBLASInt      one = 1,info,m;

339:   PetscBLASIntCast(A->rmap->n,&m);
340:   VecGetArrayRead(xx,&x);
341:   VecGetArray(yy,&y);
342:   PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
343:   if (A->factortype == MAT_FACTOR_LU) {
344: #if defined(PETSC_MISSING_LAPACK_GETRS)
345:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
346: #else
347:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
348:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
349: #endif
350:   } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
351: #if defined(PETSC_MISSING_LAPACK_POTRS)
352:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
353: #else
354:     PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
355:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
356: #endif
357:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
358:   VecRestoreArrayRead(xx,&x);
359:   VecRestoreArray(yy,&y);
360:   PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
361:   return(0);
362: }

366: static PetscErrorCode MatMatSolve_SeqDense(Mat A,Mat B,Mat X)
367: {
368:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
370:   PetscScalar    *b,*x;
371:   PetscInt       n;
372:   PetscBLASInt   nrhs,info,m;
373:   PetscBool      flg;

376:   PetscBLASIntCast(A->rmap->n,&m);
377:   PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
378:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
379:   PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
380:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");

382:   MatGetSize(B,NULL,&n);
383:   PetscBLASIntCast(n,&nrhs);
384:   MatDenseGetArray(B,&b);
385:   MatDenseGetArray(X,&x);

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

389:   if (A->factortype == MAT_FACTOR_LU) {
390: #if defined(PETSC_MISSING_LAPACK_GETRS)
391:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
392: #else
393:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
394:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
395: #endif
396:   } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
397: #if defined(PETSC_MISSING_LAPACK_POTRS)
398:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
399: #else
400:     PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info));
401:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
402: #endif
403:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");

405:   MatDenseRestoreArray(B,&b);
406:   MatDenseRestoreArray(X,&x);
407:   PetscLogFlops(nrhs*(2.0*m*m - m));
408:   return(0);
409: }

413: static PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy)
414: {
415:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
416:   PetscErrorCode    ierr;
417:   const PetscScalar *x;
418:   PetscScalar       *y;
419:   PetscBLASInt      one = 1,info,m;

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

450: static PetscErrorCode MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
451: {
452:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
453:   PetscErrorCode    ierr;
454:   const PetscScalar *x;
455:   PetscScalar       *y,sone = 1.0;
456:   Vec               tmp = 0;
457:   PetscBLASInt      one = 1,info,m;

460:   PetscBLASIntCast(A->rmap->n,&m);
461:   VecGetArrayRead(xx,&x);
462:   VecGetArray(yy,&y);
463:   if (!A->rmap->n || !A->cmap->n) return(0);
464:   if (yy == zz) {
465:     VecDuplicate(yy,&tmp);
466:     PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);
467:     VecCopy(yy,tmp);
468:   }
469:   PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
470:   /* assume if pivots exist then use LU; else Cholesky */
471:   if (mat->pivots) {
472: #if defined(PETSC_MISSING_LAPACK_GETRS)
473:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
474: #else
475:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
476:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
477: #endif
478:   } else {
479: #if defined(PETSC_MISSING_LAPACK_POTRS)
480:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
481: #else
482:     PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
483:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
484: #endif
485:   }
486:   if (tmp) {
487:     VecAXPY(yy,sone,tmp);
488:     VecDestroy(&tmp);
489:   } else {
490:     VecAXPY(yy,sone,zz);
491:   }
492:   VecRestoreArrayRead(xx,&x);
493:   VecRestoreArray(yy,&y);
494:   PetscLogFlops(2.0*A->cmap->n*A->cmap->n);
495:   return(0);
496: }

500: static PetscErrorCode MatSolveTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
501: {
502:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
503:   PetscErrorCode    ierr;
504:   const PetscScalar *x;
505:   PetscScalar       *y,sone = 1.0;
506:   Vec               tmp = NULL;
507:   PetscBLASInt      one = 1,info,m;

510:   PetscBLASIntCast(A->rmap->n,&m);
511:   if (!A->rmap->n || !A->cmap->n) return(0);
512:   VecGetArrayRead(xx,&x);
513:   VecGetArray(yy,&y);
514:   if (yy == zz) {
515:     VecDuplicate(yy,&tmp);
516:     PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);
517:     VecCopy(yy,tmp);
518:   }
519:   PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
520:   /* assume if pivots exist then use LU; else Cholesky */
521:   if (mat->pivots) {
522: #if defined(PETSC_MISSING_LAPACK_GETRS)
523:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
524: #else
525:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
526:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
527: #endif
528:   } else {
529: #if defined(PETSC_MISSING_LAPACK_POTRS)
530:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
531: #else
532:     PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
533:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
534: #endif
535:   }
536:   if (tmp) {
537:     VecAXPY(yy,sone,tmp);
538:     VecDestroy(&tmp);
539:   } else {
540:     VecAXPY(yy,sone,zz);
541:   }
542:   VecRestoreArrayRead(xx,&x);
543:   VecRestoreArray(yy,&y);
544:   PetscLogFlops(2.0*A->cmap->n*A->cmap->n);
545:   return(0);
546: }

548: /* ---------------------------------------------------------------*/
549: /* COMMENT: I have chosen to hide row permutation in the pivots,
550:    rather than put it in the Mat->row slot.*/
553: static PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo)
554: {
555: #if defined(PETSC_MISSING_LAPACK_GETRF)
557:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable.");
558: #else
559:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
561:   PetscBLASInt   n,m,info;

564:   PetscBLASIntCast(A->cmap->n,&n);
565:   PetscBLASIntCast(A->rmap->n,&m);
566:   if (!mat->pivots) {
567:     PetscMalloc1(A->rmap->n+1,&mat->pivots);
568:     PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
569:   }
570:   if (!A->rmap->n || !A->cmap->n) return(0);
571:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
572:   PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info));
573:   PetscFPTrapPop();

575:   if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization");
576:   if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization");
577:   A->ops->solve             = MatSolve_SeqDense;
578:   A->ops->solvetranspose    = MatSolveTranspose_SeqDense;
579:   A->ops->solveadd          = MatSolveAdd_SeqDense;
580:   A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
581:   A->factortype             = MAT_FACTOR_LU;

583:   PetscFree(A->solvertype);
584:   PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);

586:   PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);
587: #endif
588:   return(0);
589: }

593: static PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo)
594: {
595: #if defined(PETSC_MISSING_LAPACK_POTRF)
597:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable.");
598: #else
599:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
601:   PetscBLASInt   info,n;

604:   PetscBLASIntCast(A->cmap->n,&n);
605:   PetscFree(mat->pivots);

607:   if (!A->rmap->n || !A->cmap->n) return(0);
608:   PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info));
609:   if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1);
610:   A->ops->solve             = MatSolve_SeqDense;
611:   A->ops->solvetranspose    = MatSolveTranspose_SeqDense;
612:   A->ops->solveadd          = MatSolveAdd_SeqDense;
613:   A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
614:   A->factortype             = MAT_FACTOR_CHOLESKY;

616:   PetscFree(A->solvertype);
617:   PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);

619:   PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
620: #endif
621:   return(0);
622: }


627: PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
628: {
630:   MatFactorInfo  info;

633:   info.fill = 1.0;

635:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
636:   MatCholeskyFactor_SeqDense(fact,0,&info);
637:   return(0);
638: }

642: static PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
643: {
645:   fact->assembled                  = PETSC_TRUE;
646:   fact->preallocated               = PETSC_TRUE;
647:   fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
648:   return(0);
649: }

653: static PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info)
654: {
656:   fact->preallocated         = PETSC_TRUE;
657:   fact->assembled            = PETSC_TRUE;
658:   fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense;
659:   return(0);
660: }

664: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact)
665: {

669:   MatCreate(PetscObjectComm((PetscObject)A),fact);
670:   MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
671:   MatSetType(*fact,((PetscObject)A)->type_name);
672:   if (ftype == MAT_FACTOR_LU) {
673:     (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
674:   } else {
675:     (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
676:   }
677:   (*fact)->factortype = ftype;

679:   PetscFree((*fact)->solvertype);
680:   PetscStrallocpy(MATSOLVERPETSC,&(*fact)->solvertype);
681:   return(0);
682: }

684: /* ------------------------------------------------------------------*/
687: static PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx)
688: {
689:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
690:   PetscScalar       *x,*v = mat->v,zero = 0.0,xt;
691:   const PetscScalar *b;
692:   PetscErrorCode    ierr;
693:   PetscInt          m = A->rmap->n,i;
694:   PetscBLASInt      o = 1,bm;

697:   if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */
698:   PetscBLASIntCast(m,&bm);
699:   if (flag & SOR_ZERO_INITIAL_GUESS) {
700:     /* this is a hack fix, should have another version without the second BLASdot */
701:     VecSet(xx,zero);
702:   }
703:   VecGetArray(xx,&x);
704:   VecGetArrayRead(bb,&b);
705:   its  = its*lits;
706:   if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
707:   while (its--) {
708:     if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
709:       for (i=0; i<m; i++) {
710:         PetscStackCallBLAS("BLASdot",xt   = b[i] - BLASdot_(&bm,v+i,&bm,x,&o));
711:         x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
712:       }
713:     }
714:     if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
715:       for (i=m-1; i>=0; i--) {
716:         PetscStackCallBLAS("BLASdot",xt   = b[i] - BLASdot_(&bm,v+i,&bm,x,&o));
717:         x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
718:       }
719:     }
720:   }
721:   VecRestoreArrayRead(bb,&b);
722:   VecRestoreArray(xx,&x);
723:   return(0);
724: }

726: /* -----------------------------------------------------------------*/
729: PETSC_INTERN PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)
730: {
731:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
732:   const PetscScalar *v   = mat->v,*x;
733:   PetscScalar       *y;
734:   PetscErrorCode    ierr;
735:   PetscBLASInt      m, n,_One=1;
736:   PetscScalar       _DOne=1.0,_DZero=0.0;

739:   PetscBLASIntCast(A->rmap->n,&m);
740:   PetscBLASIntCast(A->cmap->n,&n);
741:   if (!A->rmap->n || !A->cmap->n) return(0);
742:   VecGetArrayRead(xx,&x);
743:   VecGetArray(yy,&y);
744:   PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One));
745:   VecRestoreArrayRead(xx,&x);
746:   VecRestoreArray(yy,&y);
747:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);
748:   return(0);
749: }

753: PETSC_INTERN PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy)
754: {
755:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
756:   PetscScalar       *y,_DOne=1.0,_DZero=0.0;
757:   PetscErrorCode    ierr;
758:   PetscBLASInt      m, n, _One=1;
759:   const PetscScalar *v = mat->v,*x;

762:   PetscBLASIntCast(A->rmap->n,&m);
763:   PetscBLASIntCast(A->cmap->n,&n);
764:   if (!A->rmap->n || !A->cmap->n) return(0);
765:   VecGetArrayRead(xx,&x);
766:   VecGetArray(yy,&y);
767:   PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One));
768:   VecRestoreArrayRead(xx,&x);
769:   VecRestoreArray(yy,&y);
770:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);
771:   return(0);
772: }

776: PETSC_INTERN PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
777: {
778:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
779:   const PetscScalar *v = mat->v,*x;
780:   PetscScalar       *y,_DOne=1.0;
781:   PetscErrorCode    ierr;
782:   PetscBLASInt      m, n, _One=1;

785:   PetscBLASIntCast(A->rmap->n,&m);
786:   PetscBLASIntCast(A->cmap->n,&n);
787:   if (!A->rmap->n || !A->cmap->n) return(0);
788:   if (zz != yy) {VecCopy(zz,yy);}
789:   VecGetArrayRead(xx,&x);
790:   VecGetArray(yy,&y);
791:   PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
792:   VecRestoreArrayRead(xx,&x);
793:   VecRestoreArray(yy,&y);
794:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
795:   return(0);
796: }

800: static PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
801: {
802:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
803:   const PetscScalar *v = mat->v,*x;
804:   PetscScalar       *y;
805:   PetscErrorCode    ierr;
806:   PetscBLASInt      m, n, _One=1;
807:   PetscScalar       _DOne=1.0;

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

823: /* -----------------------------------------------------------------*/
826: static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
827: {
828:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
829:   PetscScalar    *v;
831:   PetscInt       i;

834:   *ncols = A->cmap->n;
835:   if (cols) {
836:     PetscMalloc1(A->cmap->n+1,cols);
837:     for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
838:   }
839:   if (vals) {
840:     PetscMalloc1(A->cmap->n+1,vals);
841:     v    = mat->v + row;
842:     for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
843:   }
844:   return(0);
845: }

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

854:   if (cols) {PetscFree(*cols);}
855:   if (vals) {PetscFree(*vals); }
856:   return(0);
857: }
858: /* ----------------------------------------------------------------*/
861: static PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
862: {
863:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
864:   PetscInt     i,j,idx=0;

867:   if (!mat->roworiented) {
868:     if (addv == INSERT_VALUES) {
869:       for (j=0; j<n; j++) {
870:         if (indexn[j] < 0) {idx += m; continue;}
871: #if defined(PETSC_USE_DEBUG)
872:         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);
873: #endif
874:         for (i=0; i<m; i++) {
875:           if (indexm[i] < 0) {idx++; continue;}
876: #if defined(PETSC_USE_DEBUG)
877:           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);
878: #endif
879:           mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
880:         }
881:       }
882:     } else {
883:       for (j=0; j<n; j++) {
884:         if (indexn[j] < 0) {idx += m; continue;}
885: #if defined(PETSC_USE_DEBUG)
886:         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);
887: #endif
888:         for (i=0; i<m; i++) {
889:           if (indexm[i] < 0) {idx++; continue;}
890: #if defined(PETSC_USE_DEBUG)
891:           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);
892: #endif
893:           mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
894:         }
895:       }
896:     }
897:   } else {
898:     if (addv == INSERT_VALUES) {
899:       for (i=0; i<m; i++) {
900:         if (indexm[i] < 0) { idx += n; continue;}
901: #if defined(PETSC_USE_DEBUG)
902:         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);
903: #endif
904:         for (j=0; j<n; j++) {
905:           if (indexn[j] < 0) { idx++; continue;}
906: #if defined(PETSC_USE_DEBUG)
907:           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);
908: #endif
909:           mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
910:         }
911:       }
912:     } else {
913:       for (i=0; i<m; i++) {
914:         if (indexm[i] < 0) { idx += n; continue;}
915: #if defined(PETSC_USE_DEBUG)
916:         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);
917: #endif
918:         for (j=0; j<n; j++) {
919:           if (indexn[j] < 0) { idx++; continue;}
920: #if defined(PETSC_USE_DEBUG)
921:           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);
922: #endif
923:           mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
924:         }
925:       }
926:     }
927:   }
928:   return(0);
929: }

933: static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
934: {
935:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
936:   PetscInt     i,j;

939:   /* row-oriented output */
940:   for (i=0; i<m; i++) {
941:     if (indexm[i] < 0) {v += n;continue;}
942:     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);
943:     for (j=0; j<n; j++) {
944:       if (indexn[j] < 0) {v++; continue;}
945:       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);
946:       *v++ = mat->v[indexn[j]*mat->lda + indexm[i]];
947:     }
948:   }
949:   return(0);
950: }

952: /* -----------------------------------------------------------------*/

956: static PetscErrorCode MatLoad_SeqDense(Mat newmat,PetscViewer viewer)
957: {
958:   Mat_SeqDense   *a;
960:   PetscInt       *scols,i,j,nz,header[4];
961:   int            fd;
962:   PetscMPIInt    size;
963:   PetscInt       *rowlengths = 0,M,N,*cols,grows,gcols;
964:   PetscScalar    *vals,*svals,*v,*w;
965:   MPI_Comm       comm;

968:   /* force binary viewer to load .info file if it has not yet done so */
969:   PetscViewerSetUp(viewer);
970:   PetscObjectGetComm((PetscObject)viewer,&comm);
971:   MPI_Comm_size(comm,&size);
972:   if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor");
973:   PetscViewerBinaryGetDescriptor(viewer,&fd);
974:   PetscBinaryRead(fd,header,4,PETSC_INT);
975:   if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not matrix object");
976:   M = header[1]; N = header[2]; nz = header[3];

978:   /* set global size if not set already*/
979:   if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) {
980:     MatSetSizes(newmat,M,N,M,N);
981:   } else {
982:     /* if sizes and type are already set, check if the vector global sizes are correct */
983:     MatGetSize(newmat,&grows,&gcols);
984:     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);
985:   }
986:   a = (Mat_SeqDense*)newmat->data;
987:   if (!a->user_alloc) {
988:     MatSeqDenseSetPreallocation(newmat,NULL);
989:   }

991:   if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */
992:     a = (Mat_SeqDense*)newmat->data;
993:     v = a->v;
994:     /* Allocate some temp space to read in the values and then flip them
995:        from row major to column major */
996:     PetscMalloc1(M*N > 0 ? M*N : 1,&w);
997:     /* read in nonzero values */
998:     PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);
999:     /* now flip the values and store them in the matrix*/
1000:     for (j=0; j<N; j++) {
1001:       for (i=0; i<M; i++) {
1002:         *v++ =w[i*N+j];
1003:       }
1004:     }
1005:     PetscFree(w);
1006:   } else {
1007:     /* read row lengths */
1008:     PetscMalloc1(M+1,&rowlengths);
1009:     PetscBinaryRead(fd,rowlengths,M,PETSC_INT);

1011:     a = (Mat_SeqDense*)newmat->data;
1012:     v = a->v;

1014:     /* read column indices and nonzeros */
1015:     PetscMalloc1(nz+1,&scols);
1016:     cols = scols;
1017:     PetscBinaryRead(fd,cols,nz,PETSC_INT);
1018:     PetscMalloc1(nz+1,&svals);
1019:     vals = svals;
1020:     PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);

1022:     /* insert into matrix */
1023:     for (i=0; i<M; i++) {
1024:       for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j];
1025:       svals += rowlengths[i]; scols += rowlengths[i];
1026:     }
1027:     PetscFree(vals);
1028:     PetscFree(cols);
1029:     PetscFree(rowlengths);
1030:   }
1031:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1032:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
1033:   return(0);
1034: }

1038: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
1039: {
1040:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
1041:   PetscErrorCode    ierr;
1042:   PetscInt          i,j;
1043:   const char        *name;
1044:   PetscScalar       *v;
1045:   PetscViewerFormat format;
1046: #if defined(PETSC_USE_COMPLEX)
1047:   PetscBool         allreal = PETSC_TRUE;
1048: #endif

1051:   PetscViewerGetFormat(viewer,&format);
1052:   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1053:     return(0);  /* do nothing for now */
1054:   } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1055:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1056:     for (i=0; i<A->rmap->n; i++) {
1057:       v    = a->v + i;
1058:       PetscViewerASCIIPrintf(viewer,"row %D:",i);
1059:       for (j=0; j<A->cmap->n; j++) {
1060: #if defined(PETSC_USE_COMPLEX)
1061:         if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
1062:           PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1063:         } else if (PetscRealPart(*v)) {
1064:           PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));
1065:         }
1066: #else
1067:         if (*v) {
1068:           PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);
1069:         }
1070: #endif
1071:         v += a->lda;
1072:       }
1073:       PetscViewerASCIIPrintf(viewer,"\n");
1074:     }
1075:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1076:   } else {
1077:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1078: #if defined(PETSC_USE_COMPLEX)
1079:     /* determine if matrix has all real values */
1080:     v = a->v;
1081:     for (i=0; i<A->rmap->n*A->cmap->n; i++) {
1082:       if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;}
1083:     }
1084: #endif
1085:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
1086:       PetscObjectGetName((PetscObject)A,&name);
1087:       PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);
1088:       PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);
1089:       PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
1090:     }

1092:     for (i=0; i<A->rmap->n; i++) {
1093:       v = a->v + i;
1094:       for (j=0; j<A->cmap->n; j++) {
1095: #if defined(PETSC_USE_COMPLEX)
1096:         if (allreal) {
1097:           PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));
1098:         } else {
1099:           PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1100:         }
1101: #else
1102:         PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);
1103: #endif
1104:         v += a->lda;
1105:       }
1106:       PetscViewerASCIIPrintf(viewer,"\n");
1107:     }
1108:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
1109:       PetscViewerASCIIPrintf(viewer,"];\n");
1110:     }
1111:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1112:   }
1113:   PetscViewerFlush(viewer);
1114:   return(0);
1115: }

1119: static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer)
1120: {
1121:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
1122:   PetscErrorCode    ierr;
1123:   int               fd;
1124:   PetscInt          ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n;
1125:   PetscScalar       *v,*anonz,*vals;
1126:   PetscViewerFormat format;

1129:   PetscViewerBinaryGetDescriptor(viewer,&fd);

1131:   PetscViewerGetFormat(viewer,&format);
1132:   if (format == PETSC_VIEWER_NATIVE) {
1133:     /* store the matrix as a dense matrix */
1134:     PetscMalloc1(4,&col_lens);

1136:     col_lens[0] = MAT_FILE_CLASSID;
1137:     col_lens[1] = m;
1138:     col_lens[2] = n;
1139:     col_lens[3] = MATRIX_BINARY_FORMAT_DENSE;

1141:     PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);
1142:     PetscFree(col_lens);

1144:     /* write out matrix, by rows */
1145:     PetscMalloc1(m*n+1,&vals);
1146:     v    = a->v;
1147:     for (j=0; j<n; j++) {
1148:       for (i=0; i<m; i++) {
1149:         vals[j + i*n] = *v++;
1150:       }
1151:     }
1152:     PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);
1153:     PetscFree(vals);
1154:   } else {
1155:     PetscMalloc1(4+nz,&col_lens);

1157:     col_lens[0] = MAT_FILE_CLASSID;
1158:     col_lens[1] = m;
1159:     col_lens[2] = n;
1160:     col_lens[3] = nz;

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

1166:     /* Possibly should write in smaller increments, not whole matrix at once? */
1167:     /* store column indices (zero start index) */
1168:     ict = 0;
1169:     for (i=0; i<m; i++) {
1170:       for (j=0; j<n; j++) col_lens[ict++] = j;
1171:     }
1172:     PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);
1173:     PetscFree(col_lens);

1175:     /* store nonzero values */
1176:     PetscMalloc1(nz+1,&anonz);
1177:     ict  = 0;
1178:     for (i=0; i<m; i++) {
1179:       v = a->v + i;
1180:       for (j=0; j<n; j++) {
1181:         anonz[ict++] = *v; v += a->lda;
1182:       }
1183:     }
1184:     PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);
1185:     PetscFree(anonz);
1186:   }
1187:   return(0);
1188: }

1190: #include <petscdraw.h>
1193: static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa)
1194: {
1195:   Mat               A  = (Mat) Aa;
1196:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
1197:   PetscErrorCode    ierr;
1198:   PetscInt          m  = A->rmap->n,n = A->cmap->n,i,j;
1199:   int               color = PETSC_DRAW_WHITE;
1200:   PetscScalar       *v = a->v;
1201:   PetscViewer       viewer;
1202:   PetscReal         xl,yl,xr,yr,x_l,x_r,y_l,y_r;
1203:   PetscViewerFormat format;

1206:   PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
1207:   PetscViewerGetFormat(viewer,&format);
1208:   PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);

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

1212:   if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1213:     PetscDrawCollectiveBegin(draw);
1214:     /* Blue for negative and Red for positive */
1215:     for (j = 0; j < n; j++) {
1216:       x_l = j; x_r = x_l + 1.0;
1217:       for (i = 0; i < m; i++) {
1218:         y_l = m - i - 1.0;
1219:         y_r = y_l + 1.0;
1220:         if (PetscRealPart(v[j*m+i]) >  0.) {
1221:           color = PETSC_DRAW_RED;
1222:         } else if (PetscRealPart(v[j*m+i]) <  0.) {
1223:           color = PETSC_DRAW_BLUE;
1224:         } else {
1225:           continue;
1226:         }
1227:         PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1228:       }
1229:     }
1230:     PetscDrawCollectiveEnd(draw);
1231:   } else {
1232:     /* use contour shading to indicate magnitude of values */
1233:     /* first determine max of all nonzero values */
1234:     PetscReal minv = 0.0, maxv = 0.0;
1235:     PetscDraw popup;

1237:     for (i=0; i < m*n; i++) {
1238:       if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1239:     }
1240:     if (minv >= maxv) maxv = minv + PETSC_SMALL;
1241:     PetscDrawGetPopup(draw,&popup);
1242:     PetscDrawScalePopup(popup,minv,maxv);

1244:     PetscDrawCollectiveBegin(draw);
1245:     for (j=0; j<n; j++) {
1246:       x_l = j;
1247:       x_r = x_l + 1.0;
1248:       for (i=0; i<m; i++) {
1249:         y_l = m - i - 1.0;
1250:         y_r = y_l + 1.0;
1251:         color = PetscDrawRealToColor(PetscAbsScalar(v[j*m+i]),minv,maxv);
1252:         PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1253:       }
1254:     }
1255:     PetscDrawCollectiveEnd(draw);
1256:   }
1257:   return(0);
1258: }

1262: static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
1263: {
1264:   PetscDraw      draw;
1265:   PetscBool      isnull;
1266:   PetscReal      xr,yr,xl,yl,h,w;

1270:   PetscViewerDrawGetDraw(viewer,0,&draw);
1271:   PetscDrawIsNull(draw,&isnull);
1272:   if (isnull) return(0);

1274:   xr   = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
1275:   xr  += w;          yr += h;        xl = -w;     yl = -h;
1276:   PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
1277:   PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
1278:   PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
1279:   PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
1280:   PetscDrawSave(draw);
1281:   return(0);
1282: }

1286: PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
1287: {
1289:   PetscBool      iascii,isbinary,isdraw;

1292:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1293:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1294:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);

1296:   if (iascii) {
1297:     MatView_SeqDense_ASCII(A,viewer);
1298:   } else if (isbinary) {
1299:     MatView_SeqDense_Binary(A,viewer);
1300:   } else if (isdraw) {
1301:     MatView_SeqDense_Draw(A,viewer);
1302:   }
1303:   return(0);
1304: }

1308: static PetscErrorCode MatDestroy_SeqDense(Mat mat)
1309: {
1310:   Mat_SeqDense   *l = (Mat_SeqDense*)mat->data;

1314: #if defined(PETSC_USE_LOG)
1315:   PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n);
1316: #endif
1317:   PetscFree(l->pivots);
1318:   MatDestroy(&l->ptapwork);
1319:   if (!l->user_alloc) {PetscFree(l->v);}
1320:   PetscFree(mat->data);

1322:   PetscObjectChangeTypeName((PetscObject)mat,0);
1323:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);
1324:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);
1325:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);
1326: #if defined(PETSC_HAVE_ELEMENTAL)
1327:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);
1328: #endif
1329:   PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);
1330:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);
1331:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);
1332:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);
1333:   PetscObjectComposeFunction((PetscObject)mat,"MatPtAP_seqaij_seqdense_C",NULL);
1334:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);
1335:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);
1336:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);
1337:   return(0);
1338: }

1342: static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
1343: {
1344:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1346:   PetscInt       k,j,m,n,M;
1347:   PetscScalar    *v,tmp;

1350:   v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n;
1351:   if (reuse == MAT_REUSE_MATRIX && *matout == A) { /* in place transpose */
1352:     if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can not transpose non-square matrix in place");
1353:     else {
1354:       for (j=0; j<m; j++) {
1355:         for (k=0; k<j; k++) {
1356:           tmp        = v[j + k*M];
1357:           v[j + k*M] = v[k + j*M];
1358:           v[k + j*M] = tmp;
1359:         }
1360:       }
1361:     }
1362:   } else { /* out-of-place transpose */
1363:     Mat          tmat;
1364:     Mat_SeqDense *tmatd;
1365:     PetscScalar  *v2;
1366:     PetscInt     M2;

1368:     if (reuse == MAT_INITIAL_MATRIX) {
1369:       MatCreate(PetscObjectComm((PetscObject)A),&tmat);
1370:       MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);
1371:       MatSetType(tmat,((PetscObject)A)->type_name);
1372:       MatSeqDenseSetPreallocation(tmat,NULL);
1373:     } else {
1374:       tmat = *matout;
1375:     }
1376:     tmatd = (Mat_SeqDense*)tmat->data;
1377:     v = mat->v; v2 = tmatd->v; M2 = tmatd->lda;
1378:     for (j=0; j<n; j++) {
1379:       for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M];
1380:     }
1381:     MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
1382:     MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
1383:     *matout = tmat;
1384:   }
1385:   return(0);
1386: }

1390: static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool  *flg)
1391: {
1392:   Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data;
1393:   Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data;
1394:   PetscInt     i,j;
1395:   PetscScalar  *v1,*v2;

1398:   if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; return(0);}
1399:   if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; return(0);}
1400:   for (i=0; i<A1->rmap->n; i++) {
1401:     v1 = mat1->v+i; v2 = mat2->v+i;
1402:     for (j=0; j<A1->cmap->n; j++) {
1403:       if (*v1 != *v2) {*flg = PETSC_FALSE; return(0);}
1404:       v1 += mat1->lda; v2 += mat2->lda;
1405:     }
1406:   }
1407:   *flg = PETSC_TRUE;
1408:   return(0);
1409: }

1413: static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
1414: {
1415:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1417:   PetscInt       i,n,len;
1418:   PetscScalar    *x,zero = 0.0;

1421:   VecSet(v,zero);
1422:   VecGetSize(v,&n);
1423:   VecGetArray(v,&x);
1424:   len  = PetscMin(A->rmap->n,A->cmap->n);
1425:   if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
1426:   for (i=0; i<len; i++) {
1427:     x[i] = mat->v[i*mat->lda + i];
1428:   }
1429:   VecRestoreArray(v,&x);
1430:   return(0);
1431: }

1435: static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
1436: {
1437:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1438:   const PetscScalar *l,*r;
1439:   PetscScalar       x,*v;
1440:   PetscErrorCode    ierr;
1441:   PetscInt          i,j,m = A->rmap->n,n = A->cmap->n;

1444:   if (ll) {
1445:     VecGetSize(ll,&m);
1446:     VecGetArrayRead(ll,&l);
1447:     if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size");
1448:     for (i=0; i<m; i++) {
1449:       x = l[i];
1450:       v = mat->v + i;
1451:       for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;}
1452:     }
1453:     VecRestoreArrayRead(ll,&l);
1454:     PetscLogFlops(1.0*n*m);
1455:   }
1456:   if (rr) {
1457:     VecGetSize(rr,&n);
1458:     VecGetArrayRead(rr,&r);
1459:     if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size");
1460:     for (i=0; i<n; i++) {
1461:       x = r[i];
1462:       v = mat->v + i*mat->lda;
1463:       for (j=0; j<m; j++) (*v++) *= x;
1464:     }
1465:     VecRestoreArrayRead(rr,&r);
1466:     PetscLogFlops(1.0*n*m);
1467:   }
1468:   return(0);
1469: }

1473: static PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1474: {
1475:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1476:   PetscScalar    *v   = mat->v;
1477:   PetscReal      sum  = 0.0;
1478:   PetscInt       lda  =mat->lda,m=A->rmap->n,i,j;

1482:   if (type == NORM_FROBENIUS) {
1483:     if (lda>m) {
1484:       for (j=0; j<A->cmap->n; j++) {
1485:         v = mat->v+j*lda;
1486:         for (i=0; i<m; i++) {
1487:           sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1488:         }
1489:       }
1490:     } else {
1491:       for (i=0; i<A->cmap->n*A->rmap->n; i++) {
1492:         sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1493:       }
1494:     }
1495:     *nrm = PetscSqrtReal(sum);
1496:     PetscLogFlops(2.0*A->cmap->n*A->rmap->n);
1497:   } else if (type == NORM_1) {
1498:     *nrm = 0.0;
1499:     for (j=0; j<A->cmap->n; j++) {
1500:       v   = mat->v + j*mat->lda;
1501:       sum = 0.0;
1502:       for (i=0; i<A->rmap->n; i++) {
1503:         sum += PetscAbsScalar(*v);  v++;
1504:       }
1505:       if (sum > *nrm) *nrm = sum;
1506:     }
1507:     PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1508:   } else if (type == NORM_INFINITY) {
1509:     *nrm = 0.0;
1510:     for (j=0; j<A->rmap->n; j++) {
1511:       v   = mat->v + j;
1512:       sum = 0.0;
1513:       for (i=0; i<A->cmap->n; i++) {
1514:         sum += PetscAbsScalar(*v); v += mat->lda;
1515:       }
1516:       if (sum > *nrm) *nrm = sum;
1517:     }
1518:     PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1519:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
1520:   return(0);
1521: }

1525: static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
1526: {
1527:   Mat_SeqDense   *aij = (Mat_SeqDense*)A->data;

1531:   switch (op) {
1532:   case MAT_ROW_ORIENTED:
1533:     aij->roworiented = flg;
1534:     break;
1535:   case MAT_NEW_NONZERO_LOCATIONS:
1536:   case MAT_NEW_NONZERO_LOCATION_ERR:
1537:   case MAT_NEW_NONZERO_ALLOCATION_ERR:
1538:   case MAT_NEW_DIAGONALS:
1539:   case MAT_KEEP_NONZERO_PATTERN:
1540:   case MAT_IGNORE_OFF_PROC_ENTRIES:
1541:   case MAT_USE_HASH_TABLE:
1542:   case MAT_IGNORE_LOWER_TRIANGULAR:
1543:     PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
1544:     break;
1545:   case MAT_SPD:
1546:   case MAT_SYMMETRIC:
1547:   case MAT_STRUCTURALLY_SYMMETRIC:
1548:   case MAT_HERMITIAN:
1549:   case MAT_SYMMETRY_ETERNAL:
1550:     /* These options are handled directly by MatSetOption() */
1551:     break;
1552:   default:
1553:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]);
1554:   }
1555:   return(0);
1556: }

1560: static PetscErrorCode MatZeroEntries_SeqDense(Mat A)
1561: {
1562:   Mat_SeqDense   *l = (Mat_SeqDense*)A->data;
1564:   PetscInt       lda=l->lda,m=A->rmap->n,j;

1567:   if (lda>m) {
1568:     for (j=0; j<A->cmap->n; j++) {
1569:       PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));
1570:     }
1571:   } else {
1572:     PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1573:   }
1574:   return(0);
1575: }

1579: static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1580: {
1581:   PetscErrorCode    ierr;
1582:   Mat_SeqDense      *l = (Mat_SeqDense*)A->data;
1583:   PetscInt          m  = l->lda, n = A->cmap->n, i,j;
1584:   PetscScalar       *slot,*bb;
1585:   const PetscScalar *xx;

1588: #if defined(PETSC_USE_DEBUG)
1589:   for (i=0; i<N; i++) {
1590:     if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
1591:     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);
1592:   }
1593: #endif

1595:   /* fix right hand side if needed */
1596:   if (x && b) {
1597:     VecGetArrayRead(x,&xx);
1598:     VecGetArray(b,&bb);
1599:     for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
1600:     VecRestoreArrayRead(x,&xx);
1601:     VecRestoreArray(b,&bb);
1602:   }

1604:   for (i=0; i<N; i++) {
1605:     slot = l->v + rows[i];
1606:     for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
1607:   }
1608:   if (diag != 0.0) {
1609:     if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
1610:     for (i=0; i<N; i++) {
1611:       slot  = l->v + (m+1)*rows[i];
1612:       *slot = diag;
1613:     }
1614:   }
1615:   return(0);
1616: }

1620: static PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[])
1621: {
1622:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;

1625:   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");
1626:   *array = mat->v;
1627:   return(0);
1628: }

1632: static PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[])
1633: {
1635:   *array = 0; /* user cannot accidently use the array later */
1636:   return(0);
1637: }

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

1644:    Not Collective

1646:    Input Parameter:
1647: .  mat - a MATSEQDENSE or MATMPIDENSE matrix

1649:    Output Parameter:
1650: .   array - pointer to the data

1652:    Level: intermediate

1654: .seealso: MatDenseRestoreArray()
1655: @*/
1656: PetscErrorCode  MatDenseGetArray(Mat A,PetscScalar **array)
1657: {

1661:   PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));
1662:   return(0);
1663: }

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

1670:    Not Collective

1672:    Input Parameters:
1673: .  mat - a MATSEQDENSE or MATMPIDENSE matrix
1674: .  array - pointer to the data

1676:    Level: intermediate

1678: .seealso: MatDenseGetArray()
1679: @*/
1680: PetscErrorCode  MatDenseRestoreArray(Mat A,PetscScalar **array)
1681: {

1685:   PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));
1686:   return(0);
1687: }

1691: static PetscErrorCode MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B)
1692: {
1693:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1695:   PetscInt       i,j,nrows,ncols;
1696:   const PetscInt *irow,*icol;
1697:   PetscScalar    *av,*bv,*v = mat->v;
1698:   Mat            newmat;

1701:   ISGetIndices(isrow,&irow);
1702:   ISGetIndices(iscol,&icol);
1703:   ISGetLocalSize(isrow,&nrows);
1704:   ISGetLocalSize(iscol,&ncols);

1706:   /* Check submatrixcall */
1707:   if (scall == MAT_REUSE_MATRIX) {
1708:     PetscInt n_cols,n_rows;
1709:     MatGetSize(*B,&n_rows,&n_cols);
1710:     if (n_rows != nrows || n_cols != ncols) {
1711:       /* resize the result matrix to match number of requested rows/columns */
1712:       MatSetSizes(*B,nrows,ncols,nrows,ncols);
1713:     }
1714:     newmat = *B;
1715:   } else {
1716:     /* Create and fill new matrix */
1717:     MatCreate(PetscObjectComm((PetscObject)A),&newmat);
1718:     MatSetSizes(newmat,nrows,ncols,nrows,ncols);
1719:     MatSetType(newmat,((PetscObject)A)->type_name);
1720:     MatSeqDenseSetPreallocation(newmat,NULL);
1721:   }

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

1726:   for (i=0; i<ncols; i++) {
1727:     av = v + mat->lda*icol[i];
1728:     for (j=0; j<nrows; j++) *bv++ = av[irow[j]];
1729:   }

1731:   /* Assemble the matrices so that the correct flags are set */
1732:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1733:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);

1735:   /* Free work space */
1736:   ISRestoreIndices(isrow,&irow);
1737:   ISRestoreIndices(iscol,&icol);
1738:   *B   = newmat;
1739:   return(0);
1740: }

1744: static PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[])
1745: {
1747:   PetscInt       i;

1750:   if (scall == MAT_INITIAL_MATRIX) {
1751:     PetscMalloc1(n+1,B);
1752:   }

1754:   for (i=0; i<n; i++) {
1755:     MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
1756:   }
1757:   return(0);
1758: }

1762: static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
1763: {
1765:   return(0);
1766: }

1770: static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
1771: {
1773:   return(0);
1774: }

1778: static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
1779: {
1780:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data;
1782:   PetscInt       lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j;

1785:   /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
1786:   if (A->ops->copy != B->ops->copy) {
1787:     MatCopy_Basic(A,B,str);
1788:     return(0);
1789:   }
1790:   if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)");
1791:   if (lda1>m || lda2>m) {
1792:     for (j=0; j<n; j++) {
1793:       PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));
1794:     }
1795:   } else {
1796:     PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1797:   }
1798:   return(0);
1799: }

1803: static PetscErrorCode MatSetUp_SeqDense(Mat A)
1804: {

1808:    MatSeqDenseSetPreallocation(A,0);
1809:   return(0);
1810: }

1814: static PetscErrorCode MatConjugate_SeqDense(Mat A)
1815: {
1816:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1817:   PetscInt     i,nz = A->rmap->n*A->cmap->n;
1818:   PetscScalar  *aa = a->v;

1821:   for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
1822:   return(0);
1823: }

1827: static PetscErrorCode MatRealPart_SeqDense(Mat A)
1828: {
1829:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1830:   PetscInt     i,nz = A->rmap->n*A->cmap->n;
1831:   PetscScalar  *aa = a->v;

1834:   for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
1835:   return(0);
1836: }

1840: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
1841: {
1842:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1843:   PetscInt     i,nz = A->rmap->n*A->cmap->n;
1844:   PetscScalar  *aa = a->v;

1847:   for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
1848:   return(0);
1849: }

1851: /* ----------------------------------------------------------------*/
1854: PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1855: {

1859:   if (scall == MAT_INITIAL_MATRIX) {
1860:     PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
1861:     MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
1862:     PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
1863:   }
1864:   PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
1865:   MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);
1866:   PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
1867:   return(0);
1868: }

1872: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1873: {
1875:   PetscInt       m=A->rmap->n,n=B->cmap->n;
1876:   Mat            Cmat;

1879:   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);
1880:   MatCreate(PETSC_COMM_SELF,&Cmat);
1881:   MatSetSizes(Cmat,m,n,m,n);
1882:   MatSetType(Cmat,MATSEQDENSE);
1883:   MatSeqDenseSetPreallocation(Cmat,NULL);

1885:   *C = Cmat;
1886:   return(0);
1887: }

1891: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
1892: {
1893:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
1894:   Mat_SeqDense   *b = (Mat_SeqDense*)B->data;
1895:   Mat_SeqDense   *c = (Mat_SeqDense*)C->data;
1896:   PetscBLASInt   m,n,k;
1897:   PetscScalar    _DOne=1.0,_DZero=0.0;
1899:   PetscBool      flg;

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

1905:   /* Handle case where where user provided the final C matrix rather than calling MatMatMult() with MAT_INITIAL_MATRIX*/
1906:   PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
1907:   if (flg) {
1908:     C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqDense;
1909:     (*C->ops->matmultnumeric)(A,B,C);
1910:     return(0);
1911:   }

1913:   PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&flg);
1914:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"First matrix must be dense");
1915:   PetscBLASIntCast(A->rmap->n,&m);
1916:   PetscBLASIntCast(B->cmap->n,&n);
1917:   PetscBLASIntCast(A->cmap->n,&k);
1918:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
1919:   return(0);
1920: }

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

1929:   if (scall == MAT_INITIAL_MATRIX) {
1930:     PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);
1931:     MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
1932:     PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);
1933:   }
1934:   PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);
1935:   MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);
1936:   PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);
1937:   return(0);
1938: }

1942: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1943: {
1945:   PetscInt       m=A->cmap->n,n=B->cmap->n;
1946:   Mat            Cmat;

1949:   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);
1950:   MatCreate(PETSC_COMM_SELF,&Cmat);
1951:   MatSetSizes(Cmat,m,n,m,n);
1952:   MatSetType(Cmat,MATSEQDENSE);
1953:   MatSeqDenseSetPreallocation(Cmat,NULL);

1955:   Cmat->assembled = PETSC_TRUE;

1957:   *C = Cmat;
1958:   return(0);
1959: }

1963: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
1964: {
1965:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
1966:   Mat_SeqDense   *b = (Mat_SeqDense*)B->data;
1967:   Mat_SeqDense   *c = (Mat_SeqDense*)C->data;
1968:   PetscBLASInt   m,n,k;
1969:   PetscScalar    _DOne=1.0,_DZero=0.0;

1973:   PetscBLASIntCast(A->cmap->n,&m);
1974:   PetscBLASIntCast(B->cmap->n,&n);
1975:   PetscBLASIntCast(A->rmap->n,&k);
1976:   /*
1977:      Note the m and n arguments below are the number rows and columns of A', not A!
1978:   */
1979:   PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
1980:   return(0);
1981: }

1985: static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[])
1986: {
1987:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
1989:   PetscInt       i,j,m = A->rmap->n,n = A->cmap->n,p;
1990:   PetscScalar    *x;
1991:   MatScalar      *aa = a->v;

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

1996:   VecSet(v,0.0);
1997:   VecGetArray(v,&x);
1998:   VecGetLocalSize(v,&p);
1999:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2000:   for (i=0; i<m; i++) {
2001:     x[i] = aa[i]; if (idx) idx[i] = 0;
2002:     for (j=1; j<n; j++) {
2003:       if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
2004:     }
2005:   }
2006:   VecRestoreArray(v,&x);
2007:   return(0);
2008: }

2012: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])
2013: {
2014:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2016:   PetscInt       i,j,m = A->rmap->n,n = A->cmap->n,p;
2017:   PetscScalar    *x;
2018:   PetscReal      atmp;
2019:   MatScalar      *aa = a->v;

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

2024:   VecSet(v,0.0);
2025:   VecGetArray(v,&x);
2026:   VecGetLocalSize(v,&p);
2027:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2028:   for (i=0; i<m; i++) {
2029:     x[i] = PetscAbsScalar(aa[i]);
2030:     for (j=1; j<n; j++) {
2031:       atmp = PetscAbsScalar(aa[i+m*j]);
2032:       if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;}
2033:     }
2034:   }
2035:   VecRestoreArray(v,&x);
2036:   return(0);
2037: }

2041: static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])
2042: {
2043:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2045:   PetscInt       i,j,m = A->rmap->n,n = A->cmap->n,p;
2046:   PetscScalar    *x;
2047:   MatScalar      *aa = a->v;

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

2052:   VecSet(v,0.0);
2053:   VecGetArray(v,&x);
2054:   VecGetLocalSize(v,&p);
2055:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2056:   for (i=0; i<m; i++) {
2057:     x[i] = aa[i]; if (idx) idx[i] = 0;
2058:     for (j=1; j<n; j++) {
2059:       if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
2060:     }
2061:   }
2062:   VecRestoreArray(v,&x);
2063:   return(0);
2064: }

2068: static PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
2069: {
2070:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2072:   PetscScalar    *x;

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

2077:   VecGetArray(v,&x);
2078:   PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));
2079:   VecRestoreArray(v,&x);
2080:   return(0);
2081: }


2086: PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms)
2087: {
2089:   PetscInt       i,j,m,n;
2090:   PetscScalar    *a;

2093:   MatGetSize(A,&m,&n);
2094:   PetscMemzero(norms,n*sizeof(PetscReal));
2095:   MatDenseGetArray(A,&a);
2096:   if (type == NORM_2) {
2097:     for (i=0; i<n; i++) {
2098:       for (j=0; j<m; j++) {
2099:         norms[i] += PetscAbsScalar(a[j]*a[j]);
2100:       }
2101:       a += m;
2102:     }
2103:   } else if (type == NORM_1) {
2104:     for (i=0; i<n; i++) {
2105:       for (j=0; j<m; j++) {
2106:         norms[i] += PetscAbsScalar(a[j]);
2107:       }
2108:       a += m;
2109:     }
2110:   } else if (type == NORM_INFINITY) {
2111:     for (i=0; i<n; i++) {
2112:       for (j=0; j<m; j++) {
2113:         norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]);
2114:       }
2115:       a += m;
2116:     }
2117:   } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType");
2118:   MatDenseRestoreArray(A,&a);
2119:   if (type == NORM_2) {
2120:     for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]);
2121:   }
2122:   return(0);
2123: }

2127: static PetscErrorCode  MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
2128: {
2130:   PetscScalar    *a;
2131:   PetscInt       m,n,i;

2134:   MatGetSize(x,&m,&n);
2135:   MatDenseGetArray(x,&a);
2136:   for (i=0; i<m*n; i++) {
2137:     PetscRandomGetValue(rctx,a+i);
2138:   }
2139:   MatDenseRestoreArray(x,&a);
2140:   return(0);
2141: }

2145: static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool  *missing,PetscInt *d)
2146: {
2148:   *missing = PETSC_FALSE;
2149:   return(0);
2150: }


2153: /* -------------------------------------------------------------------*/
2154: static struct _MatOps MatOps_Values = { MatSetValues_SeqDense,
2155:                                         MatGetRow_SeqDense,
2156:                                         MatRestoreRow_SeqDense,
2157:                                         MatMult_SeqDense,
2158:                                 /*  4*/ MatMultAdd_SeqDense,
2159:                                         MatMultTranspose_SeqDense,
2160:                                         MatMultTransposeAdd_SeqDense,
2161:                                         0,
2162:                                         0,
2163:                                         0,
2164:                                 /* 10*/ 0,
2165:                                         MatLUFactor_SeqDense,
2166:                                         MatCholeskyFactor_SeqDense,
2167:                                         MatSOR_SeqDense,
2168:                                         MatTranspose_SeqDense,
2169:                                 /* 15*/ MatGetInfo_SeqDense,
2170:                                         MatEqual_SeqDense,
2171:                                         MatGetDiagonal_SeqDense,
2172:                                         MatDiagonalScale_SeqDense,
2173:                                         MatNorm_SeqDense,
2174:                                 /* 20*/ MatAssemblyBegin_SeqDense,
2175:                                         MatAssemblyEnd_SeqDense,
2176:                                         MatSetOption_SeqDense,
2177:                                         MatZeroEntries_SeqDense,
2178:                                 /* 24*/ MatZeroRows_SeqDense,
2179:                                         0,
2180:                                         0,
2181:                                         0,
2182:                                         0,
2183:                                 /* 29*/ MatSetUp_SeqDense,
2184:                                         0,
2185:                                         0,
2186:                                         0,
2187:                                         0,
2188:                                 /* 34*/ MatDuplicate_SeqDense,
2189:                                         0,
2190:                                         0,
2191:                                         0,
2192:                                         0,
2193:                                 /* 39*/ MatAXPY_SeqDense,
2194:                                         MatGetSubMatrices_SeqDense,
2195:                                         0,
2196:                                         MatGetValues_SeqDense,
2197:                                         MatCopy_SeqDense,
2198:                                 /* 44*/ MatGetRowMax_SeqDense,
2199:                                         MatScale_SeqDense,
2200:                                         MatShift_Basic,
2201:                                         0,
2202:                                         MatZeroRowsColumns_SeqDense,
2203:                                 /* 49*/ MatSetRandom_SeqDense,
2204:                                         0,
2205:                                         0,
2206:                                         0,
2207:                                         0,
2208:                                 /* 54*/ 0,
2209:                                         0,
2210:                                         0,
2211:                                         0,
2212:                                         0,
2213:                                 /* 59*/ 0,
2214:                                         MatDestroy_SeqDense,
2215:                                         MatView_SeqDense,
2216:                                         0,
2217:                                         0,
2218:                                 /* 64*/ 0,
2219:                                         0,
2220:                                         0,
2221:                                         0,
2222:                                         0,
2223:                                 /* 69*/ MatGetRowMaxAbs_SeqDense,
2224:                                         0,
2225:                                         0,
2226:                                         0,
2227:                                         0,
2228:                                 /* 74*/ 0,
2229:                                         0,
2230:                                         0,
2231:                                         0,
2232:                                         0,
2233:                                 /* 79*/ 0,
2234:                                         0,
2235:                                         0,
2236:                                         0,
2237:                                 /* 83*/ MatLoad_SeqDense,
2238:                                         0,
2239:                                         MatIsHermitian_SeqDense,
2240:                                         0,
2241:                                         0,
2242:                                         0,
2243:                                 /* 89*/ MatMatMult_SeqDense_SeqDense,
2244:                                         MatMatMultSymbolic_SeqDense_SeqDense,
2245:                                         MatMatMultNumeric_SeqDense_SeqDense,
2246:                                         MatPtAP_SeqDense_SeqDense,
2247:                                         MatPtAPSymbolic_SeqDense_SeqDense,
2248:                                 /* 94*/ MatPtAPNumeric_SeqDense_SeqDense,
2249:                                         0,
2250:                                         0,
2251:                                         0,
2252:                                         0,
2253:                                 /* 99*/ 0,
2254:                                         0,
2255:                                         0,
2256:                                         MatConjugate_SeqDense,
2257:                                         0,
2258:                                 /*104*/ 0,
2259:                                         MatRealPart_SeqDense,
2260:                                         MatImaginaryPart_SeqDense,
2261:                                         0,
2262:                                         0,
2263:                                 /*109*/ MatMatSolve_SeqDense,
2264:                                         0,
2265:                                         MatGetRowMin_SeqDense,
2266:                                         MatGetColumnVector_SeqDense,
2267:                                         MatMissingDiagonal_SeqDense,
2268:                                 /*114*/ 0,
2269:                                         0,
2270:                                         0,
2271:                                         0,
2272:                                         0,
2273:                                 /*119*/ 0,
2274:                                         0,
2275:                                         0,
2276:                                         0,
2277:                                         0,
2278:                                 /*124*/ 0,
2279:                                         MatGetColumnNorms_SeqDense,
2280:                                         0,
2281:                                         0,
2282:                                         0,
2283:                                 /*129*/ 0,
2284:                                         MatTransposeMatMult_SeqDense_SeqDense,
2285:                                         MatTransposeMatMultSymbolic_SeqDense_SeqDense,
2286:                                         MatTransposeMatMultNumeric_SeqDense_SeqDense,
2287:                                         0,
2288:                                 /*134*/ 0,
2289:                                         0,
2290:                                         0,
2291:                                         0,
2292:                                         0,
2293:                                 /*139*/ 0,
2294:                                         0,
2295:                                         0
2296: };

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

2305:    Collective on MPI_Comm

2307:    Input Parameters:
2308: +  comm - MPI communicator, set to PETSC_COMM_SELF
2309: .  m - number of rows
2310: .  n - number of columns
2311: -  data - optional location of matrix data in column major order.  Set data=NULL for PETSc
2312:    to control all matrix memory allocation.

2314:    Output Parameter:
2315: .  A - the matrix

2317:    Notes:
2318:    The data input variable is intended primarily for Fortran programmers
2319:    who wish to allocate their own matrix memory space.  Most users should
2320:    set data=NULL.

2322:    Level: intermediate

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

2326: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
2327: @*/
2328: PetscErrorCode  MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
2329: {

2333:   MatCreate(comm,A);
2334:   MatSetSizes(*A,m,n,m,n);
2335:   MatSetType(*A,MATSEQDENSE);
2336:   MatSeqDenseSetPreallocation(*A,data);
2337:   return(0);
2338: }

2342: /*@C
2343:    MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements

2345:    Collective on MPI_Comm

2347:    Input Parameters:
2348: +  B - the matrix
2349: -  data - the array (or NULL)

2351:    Notes:
2352:    The data input variable is intended primarily for Fortran programmers
2353:    who wish to allocate their own matrix memory space.  Most users should
2354:    need not call this routine.

2356:    Level: intermediate

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

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

2362: @*/
2363: PetscErrorCode  MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
2364: {

2368:   PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
2369:   return(0);
2370: }

2374: PetscErrorCode  MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
2375: {
2376:   Mat_SeqDense   *b;

2380:   B->preallocated = PETSC_TRUE;

2382:   PetscLayoutSetUp(B->rmap);
2383:   PetscLayoutSetUp(B->cmap);

2385:   b       = (Mat_SeqDense*)B->data;
2386:   b->Mmax = B->rmap->n;
2387:   b->Nmax = B->cmap->n;
2388:   if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n;

2390:   if (!data) { /* petsc-allocated storage */
2391:     if (!b->user_alloc) { PetscFree(b->v); }
2392:     PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);
2393:     PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));

2395:     b->user_alloc = PETSC_FALSE;
2396:   } else { /* user-allocated storage */
2397:     if (!b->user_alloc) { PetscFree(b->v); }
2398:     b->v          = data;
2399:     b->user_alloc = PETSC_TRUE;
2400:   }
2401:   B->assembled = PETSC_TRUE;
2402:   return(0);
2403: }

2405: #if defined(PETSC_HAVE_ELEMENTAL)
2408: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
2409: {
2410:   Mat            mat_elemental;
2412:   PetscScalar    *array,*v_colwise;
2413:   PetscInt       M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols;

2416:   PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);
2417:   MatDenseGetArray(A,&array);
2418:   /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
2419:   k = 0;
2420:   for (j=0; j<N; j++) {
2421:     cols[j] = j;
2422:     for (i=0; i<M; i++) {
2423:       v_colwise[j*M+i] = array[k++];
2424:     }
2425:   }
2426:   for (i=0; i<M; i++) {
2427:     rows[i] = i;
2428:   }
2429:   MatDenseRestoreArray(A,&array);

2431:   MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
2432:   MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);
2433:   MatSetType(mat_elemental,MATELEMENTAL);
2434:   MatSetUp(mat_elemental);

2436:   /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
2437:   MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);
2438:   MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
2439:   MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);
2440:   PetscFree3(v_colwise,rows,cols);

2442:   if (reuse == MAT_INPLACE_MATRIX) {
2443:     MatHeaderReplace(A,&mat_elemental);
2444:   } else {
2445:     *newmat = mat_elemental;
2446:   }
2447:   return(0);
2448: }
2449: #endif

2453: /*@C
2454:   MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array

2456:   Input parameter:
2457: + A - the matrix
2458: - lda - the leading dimension

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

2465:   Level: intermediate

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

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

2471: @*/
2472: PetscErrorCode  MatSeqDenseSetLDA(Mat B,PetscInt lda)
2473: {
2474:   Mat_SeqDense *b = (Mat_SeqDense*)B->data;

2477:   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);
2478:   b->lda       = lda;
2479:   b->changelda = PETSC_FALSE;
2480:   b->Mmax      = PetscMax(b->Mmax,lda);
2481:   return(0);
2482: }

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

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

2490:   Level: beginner

2492: .seealso: MatCreateSeqDense()

2494: M*/

2498: PETSC_EXTERN PetscErrorCode MatCreate_SeqDense(Mat B)
2499: {
2500:   Mat_SeqDense   *b;
2502:   PetscMPIInt    size;

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

2508:   PetscNewLog(B,&b);
2509:   PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
2510:   B->data = (void*)b;

2512:   b->pivots      = 0;
2513:   b->roworiented = PETSC_TRUE;
2514:   b->v           = 0;
2515:   b->changelda   = PETSC_FALSE;

2517:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);
2518:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);
2519:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);
2520: #if defined(PETSC_HAVE_ELEMENTAL)
2521:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);
2522: #endif
2523:   PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);
2524:   PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2525:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2526:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2527:   PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaij_seqdense_C",MatPtAP_SeqDense_SeqDense);
2528:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2529:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2530:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2531:   PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
2532:   return(0);
2533: }