Actual source code: dense.c

petsc-3.9.0 2018-04-07
Report Typos and Errors

  2: /*
  3:      Defines the basic matrix operations for sequential dense.
  4: */

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

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

 11: static PetscErrorCode MatSeqDenseSymmetrize_Private(Mat A, PetscBool hermitian)
 12: {
 13:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
 14:   PetscInt      j, k, n = A->rmap->n;

 17:   if (A->rmap->n != A->cmap->n) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot symmetrize a rectangular matrix");
 18:   if (!hermitian) {
 19:     for (k=0;k<n;k++) {
 20:       for (j=k;j<n;j++) {
 21:         mat->v[j*mat->lda + k] = mat->v[k*mat->lda + j];
 22:       }
 23:     }
 24:   } else {
 25:     for (k=0;k<n;k++) {
 26:       for (j=k;j<n;j++) {
 27:         mat->v[j*mat->lda + k] = PetscConj(mat->v[k*mat->lda + j]);
 28:       }
 29:     }
 30:   }
 31:   return(0);
 32: }

 34: PETSC_EXTERN PetscErrorCode MatSeqDenseInvertFactors_Private(Mat A)
 35: {
 36: #if defined(PETSC_MISSING_LAPACK_POTRF)
 38:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable.");
 39: #else
 40:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
 42:   PetscBLASInt   info,n;

 45:   if (!A->rmap->n || !A->cmap->n) return(0);
 46:   PetscBLASIntCast(A->cmap->n,&n);
 47:   if (A->factortype == MAT_FACTOR_LU) {
 48:     if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present");
 49:     if (!mat->fwork) {
 50:       mat->lfwork = n;
 51:       PetscMalloc1(mat->lfwork,&mat->fwork);
 52:       PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
 53:     }
 54:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
 55:     PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
 56:     PetscFPTrapPop();
 57:     PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0); /* TODO CHECK FLOPS */
 58:   } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
 59:     if (A->spd) {
 60:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
 61:       PetscStackCallBLAS("LAPACKpotri",LAPACKpotri_("L",&n,mat->v,&mat->lda,&info));
 62:       PetscFPTrapPop();
 63:       MatSeqDenseSymmetrize_Private(A,PETSC_TRUE);
 64: #if defined(PETSC_USE_COMPLEX)
 65:     } else if (A->hermitian) {
 66:       if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present");
 67:       if (!mat->fwork) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Fwork not present");
 68:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
 69:       PetscStackCallBLAS("LAPACKhetri",LAPACKhetri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info));
 70:       PetscFPTrapPop();
 71:       MatSeqDenseSymmetrize_Private(A,PETSC_TRUE);
 72: #endif
 73:     } else { /* symmetric case */
 74:       if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present");
 75:       if (!mat->fwork) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Fwork not present");
 76:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
 77:       PetscStackCallBLAS("LAPACKsytri",LAPACKsytri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info));
 78:       PetscFPTrapPop();
 79:       MatSeqDenseSymmetrize_Private(A,PETSC_FALSE);
 80:     }
 81:     if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad Inversion: zero pivot in row %D",(PetscInt)info-1);
 82:     PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0); /* TODO CHECK FLOPS */
 83:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
 84: #endif

 86:   A->ops->solve             = NULL;
 87:   A->ops->matsolve          = NULL;
 88:   A->ops->solvetranspose    = NULL;
 89:   A->ops->matsolvetranspose = NULL;
 90:   A->ops->solveadd          = NULL;
 91:   A->ops->solvetransposeadd = NULL;
 92:   A->factortype             = MAT_FACTOR_NONE;
 93:   PetscFree(A->solvertype);
 94:   return(0);
 95: }

 97: PetscErrorCode MatZeroRowsColumns_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
 98: {
 99:   PetscErrorCode    ierr;
100:   Mat_SeqDense      *l = (Mat_SeqDense*)A->data;
101:   PetscInt          m  = l->lda, n = A->cmap->n,r = A->rmap->n, i,j;
102:   PetscScalar       *slot,*bb;
103:   const PetscScalar *xx;

106: #if defined(PETSC_USE_DEBUG)
107:   for (i=0; i<N; i++) {
108:     if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
109:     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);
110:     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);
111:   }
112: #endif

114:   /* fix right hand side if needed */
115:   if (x && b) {
116:     VecGetArrayRead(x,&xx);
117:     VecGetArray(b,&bb);
118:     for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
119:     VecRestoreArrayRead(x,&xx);
120:     VecRestoreArray(b,&bb);
121:   }

123:   for (i=0; i<N; i++) {
124:     slot = l->v + rows[i]*m;
125:     PetscMemzero(slot,r*sizeof(PetscScalar));
126:   }
127:   for (i=0; i<N; i++) {
128:     slot = l->v + rows[i];
129:     for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
130:   }
131:   if (diag != 0.0) {
132:     if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
133:     for (i=0; i<N; i++) {
134:       slot  = l->v + (m+1)*rows[i];
135:       *slot = diag;
136:     }
137:   }
138:   return(0);
139: }

141: PetscErrorCode MatPtAPNumeric_SeqDense_SeqDense(Mat A,Mat P,Mat C)
142: {
143:   Mat_SeqDense   *c = (Mat_SeqDense*)(C->data);

147:   MatMatMultNumeric_SeqDense_SeqDense(A,P,c->ptapwork);
148:   MatTransposeMatMultNumeric_SeqDense_SeqDense(P,c->ptapwork,C);
149:   return(0);
150: }

152: PetscErrorCode MatPtAPSymbolic_SeqDense_SeqDense(Mat A,Mat P,PetscReal fill,Mat *C)
153: {
154:   Mat_SeqDense   *c;

158:   MatCreateSeqDense(PetscObjectComm((PetscObject)A),P->cmap->N,P->cmap->N,NULL,C);
159:   c = (Mat_SeqDense*)((*C)->data);
160:   MatCreateSeqDense(PetscObjectComm((PetscObject)A),A->rmap->N,P->cmap->N,NULL,&c->ptapwork);
161:   return(0);
162: }

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

169:   if (reuse == MAT_INITIAL_MATRIX) {
170:     MatPtAPSymbolic_SeqDense_SeqDense(A,P,fill,C);
171:   }
172:   PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);
173:   (*(*C)->ops->ptapnumeric)(A,P,*C);
174:   PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);
175:   return(0);
176: }

178: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
179: {
180:   Mat            B = NULL;
181:   Mat_SeqAIJ     *a = (Mat_SeqAIJ*)A->data;
182:   Mat_SeqDense   *b;
184:   PetscInt       *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->N,i;
185:   MatScalar      *av=a->a;
186:   PetscBool      isseqdense;

189:   if (reuse == MAT_REUSE_MATRIX) {
190:     PetscObjectTypeCompare((PetscObject)*newmat,MATSEQDENSE,&isseqdense);
191:     if (!isseqdense) SETERRQ1(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_USER,"Cannot reuse matrix of type %s",((PetscObject)(*newmat))->type);
192:   }
193:   if (reuse != MAT_REUSE_MATRIX) {
194:     MatCreate(PetscObjectComm((PetscObject)A),&B);
195:     MatSetSizes(B,m,n,m,n);
196:     MatSetType(B,MATSEQDENSE);
197:     MatSeqDenseSetPreallocation(B,NULL);
198:     b    = (Mat_SeqDense*)(B->data);
199:   } else {
200:     b    = (Mat_SeqDense*)((*newmat)->data);
201:     PetscMemzero(b->v,m*n*sizeof(PetscScalar));
202:   }
203:   for (i=0; i<m; i++) {
204:     PetscInt j;
205:     for (j=0;j<ai[1]-ai[0];j++) {
206:       b->v[*aj*m+i] = *av;
207:       aj++;
208:       av++;
209:     }
210:     ai++;
211:   }

213:   if (reuse == MAT_INPLACE_MATRIX) {
214:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
215:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
216:     MatHeaderReplace(A,&B);
217:   } else {
218:     if (B) *newmat = B;
219:     MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
220:     MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
221:   }
222:   return(0);
223: }

225: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
226: {
227:   Mat            B;
228:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
230:   PetscInt       i, j;
231:   PetscInt       *rows, *nnz;
232:   MatScalar      *aa = a->v, *vals;

235:   MatCreate(PetscObjectComm((PetscObject)A),&B);
236:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
237:   MatSetType(B,MATSEQAIJ);
238:   PetscCalloc3(A->rmap->n,&rows,A->rmap->n,&nnz,A->rmap->n,&vals);
239:   for (j=0; j<A->cmap->n; j++) {
240:     for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) ++nnz[i];
241:     aa += a->lda;
242:   }
243:   MatSeqAIJSetPreallocation(B,PETSC_DETERMINE,nnz);
244:   aa = a->v;
245:   for (j=0; j<A->cmap->n; j++) {
246:     PetscInt numRows = 0;
247:     for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) {rows[numRows] = i; vals[numRows++] = aa[i];}
248:     MatSetValues(B,numRows,rows,1,&j,vals,INSERT_VALUES);
249:     aa  += a->lda;
250:   }
251:   PetscFree3(rows,nnz,vals);
252:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
253:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

255:   if (reuse == MAT_INPLACE_MATRIX) {
256:     MatHeaderReplace(A,&B);
257:   } else {
258:     *newmat = B;
259:   }
260:   return(0);
261: }

263: static PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str)
264: {
265:   Mat_SeqDense   *x     = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data;
266:   PetscScalar    oalpha = alpha;
267:   PetscInt       j;
268:   PetscBLASInt   N,m,ldax,lday,one = 1;

272:   PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);
273:   PetscBLASIntCast(X->rmap->n,&m);
274:   PetscBLASIntCast(x->lda,&ldax);
275:   PetscBLASIntCast(y->lda,&lday);
276:   if (ldax>m || lday>m) {
277:     for (j=0; j<X->cmap->n; j++) {
278:       PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&oalpha,x->v+j*ldax,&one,y->v+j*lday,&one));
279:     }
280:   } else {
281:     PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&oalpha,x->v,&one,y->v,&one));
282:   }
283:   PetscObjectStateIncrease((PetscObject)Y);
284:   PetscLogFlops(PetscMax(2*N-1,0));
285:   return(0);
286: }

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

293:   info->block_size        = 1.0;
294:   info->nz_allocated      = (double)N;
295:   info->nz_used           = (double)N;
296:   info->nz_unneeded       = (double)0;
297:   info->assemblies        = (double)A->num_ass;
298:   info->mallocs           = 0;
299:   info->memory            = ((PetscObject)A)->mem;
300:   info->fill_ratio_given  = 0;
301:   info->fill_ratio_needed = 0;
302:   info->factor_mallocs    = 0;
303:   return(0);
304: }

306: static PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha)
307: {
308:   Mat_SeqDense   *a     = (Mat_SeqDense*)A->data;
309:   PetscScalar    oalpha = alpha;
311:   PetscBLASInt   one = 1,j,nz,lda;

314:   PetscBLASIntCast(a->lda,&lda);
315:   if (lda>A->rmap->n) {
316:     PetscBLASIntCast(A->rmap->n,&nz);
317:     for (j=0; j<A->cmap->n; j++) {
318:       PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v+j*lda,&one));
319:     }
320:   } else {
321:     PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);
322:     PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v,&one));
323:   }
324:   PetscLogFlops(nz);
325:   return(0);
326: }

328: static PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool  *fl)
329: {
330:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
331:   PetscInt     i,j,m = A->rmap->n,N;
332:   PetscScalar  *v = a->v;

335:   *fl = PETSC_FALSE;
336:   if (A->rmap->n != A->cmap->n) return(0);
337:   N = a->lda;

339:   for (i=0; i<m; i++) {
340:     for (j=i+1; j<m; j++) {
341:       if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) return(0);
342:     }
343:   }
344:   *fl = PETSC_TRUE;
345:   return(0);
346: }

348: static PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues)
349: {
350:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data,*l;
352:   PetscInt       lda = (PetscInt)mat->lda,j,m;

355:   PetscLayoutReference(A->rmap,&newi->rmap);
356:   PetscLayoutReference(A->cmap,&newi->cmap);
357:   MatSeqDenseSetPreallocation(newi,NULL);
358:   if (cpvalues == MAT_COPY_VALUES) {
359:     l = (Mat_SeqDense*)newi->data;
360:     if (lda>A->rmap->n) {
361:       m = A->rmap->n;
362:       for (j=0; j<A->cmap->n; j++) {
363:         PetscMemcpy(l->v+j*m,mat->v+j*lda,m*sizeof(PetscScalar));
364:       }
365:     } else {
366:       PetscMemcpy(l->v,mat->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
367:     }
368:   }
369:   newi->assembled = PETSC_TRUE;
370:   return(0);
371: }

373: static PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat)
374: {

378:   MatCreate(PetscObjectComm((PetscObject)A),newmat);
379:   MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
380:   MatSetType(*newmat,((PetscObject)A)->type_name);
381:   MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);
382:   return(0);
383: }


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

388: static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
389: {
390:   MatFactorInfo  info;

394:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
395:   MatLUFactor_SeqDense(fact,0,0,&info);
396:   return(0);
397: }

399: static PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy)
400: {
401:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
402:   PetscErrorCode    ierr;
403:   const PetscScalar *x;
404:   PetscScalar       *y;
405:   PetscBLASInt      one = 1,info,m;

408:   PetscBLASIntCast(A->rmap->n,&m);
409:   VecGetArrayRead(xx,&x);
410:   VecGetArray(yy,&y);
411:   PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
412:   if (A->factortype == MAT_FACTOR_LU) {
413: #if defined(PETSC_MISSING_LAPACK_GETRS)
414:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
415: #else
416:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
417:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
418:     PetscFPTrapPop();
419:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
420: #endif
421:   } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
422: #if defined(PETSC_MISSING_LAPACK_POTRS)
423:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
424: #else
425:     if (A->spd) {
426:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
427:       PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
428:       PetscFPTrapPop();
429:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
430: #if defined(PETSC_USE_COMPLEX)
431:     } else if (A->hermitian) {
432:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
433:       PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
434:       PetscFPTrapPop();
435:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve");
436: #endif
437:     } else { /* symmetric case */
438:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
439:       PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
440:       PetscFPTrapPop();
441:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
442:     }
443: #endif
444:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
445:   VecRestoreArrayRead(xx,&x);
446:   VecRestoreArray(yy,&y);
447:   PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
448:   return(0);
449: }

451: static PetscErrorCode MatMatSolve_SeqDense(Mat A,Mat B,Mat X)
452: {
453:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
455:   PetscScalar    *b,*x;
456:   PetscInt       n;
457:   PetscBLASInt   nrhs,info,m;
458:   PetscBool      flg;

461:   PetscBLASIntCast(A->rmap->n,&m);
462:   PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
463:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
464:   PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
465:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");

467:   MatGetSize(B,NULL,&n);
468:   PetscBLASIntCast(n,&nrhs);
469:   MatDenseGetArray(B,&b);
470:   MatDenseGetArray(X,&x);

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

474:   if (A->factortype == MAT_FACTOR_LU) {
475: #if defined(PETSC_MISSING_LAPACK_GETRS)
476:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
477: #else
478:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
479:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
480:     PetscFPTrapPop();
481:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
482: #endif
483:   } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
484: #if defined(PETSC_MISSING_LAPACK_POTRS)
485:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
486: #else
487:     if (A->spd) {
488:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
489:       PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info));
490:       PetscFPTrapPop();
491:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
492: #if defined(PETSC_USE_COMPLEX)
493:     } else if (A->hermitian) {
494:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
495:       PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
496:       PetscFPTrapPop();
497:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve");
498: #endif
499:     } else { /* symmetric case */
500:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
501:       PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
502:       PetscFPTrapPop();
503:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
504:     }
505: #endif
506:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");

508:   MatDenseRestoreArray(B,&b);
509:   MatDenseRestoreArray(X,&x);
510:   PetscLogFlops(nrhs*(2.0*m*m - m));
511:   return(0);
512: }

514: static PetscErrorCode MatConjugate_SeqDense(Mat);

516: static PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy)
517: {
518:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
519:   PetscErrorCode    ierr;
520:   const PetscScalar *x;
521:   PetscScalar       *y;
522:   PetscBLASInt      one = 1,info,m;

525:   PetscBLASIntCast(A->rmap->n,&m);
526:   VecGetArrayRead(xx,&x);
527:   VecGetArray(yy,&y);
528:   PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
529:   if (A->factortype == MAT_FACTOR_LU) {
530: #if defined(PETSC_MISSING_LAPACK_GETRS)
531:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
532: #else
533:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
534:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
535:     PetscFPTrapPop();
536:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve");
537: #endif
538:   } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
539: #if defined(PETSC_MISSING_LAPACK_POTRS)
540:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
541: #else
542:     if (A->spd) {
543: #if defined(PETSC_USE_COMPLEX)
544:       MatConjugate_SeqDense(A);
545: #endif
546:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
547:       PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
548:       PetscFPTrapPop();
549: #if defined(PETSC_USE_COMPLEX)
550:       MatConjugate_SeqDense(A);
551: #endif
552:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
553: #if defined(PETSC_USE_COMPLEX)
554:     } else if (A->hermitian) {
555:       MatConjugate_SeqDense(A);
556:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
557:       PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
558:       PetscFPTrapPop();
559:       MatConjugate_SeqDense(A);
560: #endif
561:     } else { /* symmetric case */
562:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
563:       PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
564:       PetscFPTrapPop();
565:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
566:     }
567: #endif
568:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
569:   VecRestoreArrayRead(xx,&x);
570:   VecRestoreArray(yy,&y);
571:   PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
572:   return(0);
573: }

575: static PetscErrorCode MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
576: {
577:   PetscErrorCode    ierr;
578:   const PetscScalar *x;
579:   PetscScalar       *y,sone = 1.0;
580:   Vec               tmp = 0;

583:   VecGetArrayRead(xx,&x);
584:   VecGetArray(yy,&y);
585:   if (!A->rmap->n || !A->cmap->n) return(0);
586:   if (yy == zz) {
587:     VecDuplicate(yy,&tmp);
588:     PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);
589:     VecCopy(yy,tmp);
590:   }
591:   MatSolve_SeqDense(A,xx,yy);
592:   if (tmp) {
593:     VecAXPY(yy,sone,tmp);
594:     VecDestroy(&tmp);
595:   } else {
596:     VecAXPY(yy,sone,zz);
597:   }
598:   VecRestoreArrayRead(xx,&x);
599:   VecRestoreArray(yy,&y);
600:   PetscLogFlops(A->cmap->n);
601:   return(0);
602: }

604: static PetscErrorCode MatSolveTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
605: {
606:   PetscErrorCode    ierr;
607:   const PetscScalar *x;
608:   PetscScalar       *y,sone = 1.0;
609:   Vec               tmp = NULL;

612:   if (!A->rmap->n || !A->cmap->n) return(0);
613:   VecGetArrayRead(xx,&x);
614:   VecGetArray(yy,&y);
615:   if (yy == zz) {
616:     VecDuplicate(yy,&tmp);
617:     PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);
618:     VecCopy(yy,tmp);
619:   }
620:   MatSolveTranspose_SeqDense(A,xx,yy);
621:   if (tmp) {
622:     VecAXPY(yy,sone,tmp);
623:     VecDestroy(&tmp);
624:   } else {
625:     VecAXPY(yy,sone,zz);
626:   }
627:   VecRestoreArrayRead(xx,&x);
628:   VecRestoreArray(yy,&y);
629:   PetscLogFlops(A->cmap->n);
630:   return(0);
631: }

633: /* ---------------------------------------------------------------*/
634: /* COMMENT: I have chosen to hide row permutation in the pivots,
635:    rather than put it in the Mat->row slot.*/
636: static PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo)
637: {
638: #if defined(PETSC_MISSING_LAPACK_GETRF)
640:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable.");
641: #else
642:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
644:   PetscBLASInt   n,m,info;

647:   PetscBLASIntCast(A->cmap->n,&n);
648:   PetscBLASIntCast(A->rmap->n,&m);
649:   if (!mat->pivots) {
650:     PetscMalloc1(A->rmap->n,&mat->pivots);
651:     PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
652:   }
653:   if (!A->rmap->n || !A->cmap->n) return(0);
654:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
655:   PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info));
656:   PetscFPTrapPop();

658:   if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization");
659:   if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization");

661:   A->ops->solve             = MatSolve_SeqDense;
662:   A->ops->matsolve          = MatMatSolve_SeqDense;
663:   A->ops->solvetranspose    = MatSolveTranspose_SeqDense;
664:   A->ops->solveadd          = MatSolveAdd_SeqDense;
665:   A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
666:   A->factortype             = MAT_FACTOR_LU;

668:   PetscFree(A->solvertype);
669:   PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);

671:   PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);
672: #endif
673:   return(0);
674: }

676: /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */
677: static PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo)
678: {
679: #if defined(PETSC_MISSING_LAPACK_POTRF)
681:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable.");
682: #else
683:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
685:   PetscBLASInt   info,n;

688:   PetscBLASIntCast(A->cmap->n,&n);
689:   if (!A->rmap->n || !A->cmap->n) return(0);
690:   if (A->spd) {
691:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
692:     PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info));
693:     PetscFPTrapPop();
694: #if defined(PETSC_USE_COMPLEX)
695:   } else if (A->hermitian) {
696:     if (!mat->pivots) {
697:       PetscMalloc1(A->rmap->n,&mat->pivots);
698:       PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
699:     }
700:     if (!mat->fwork) {
701:       PetscScalar dummy;

703:       mat->lfwork = -1;
704:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
705:       PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
706:       PetscFPTrapPop();
707:       mat->lfwork = (PetscInt)PetscRealPart(dummy);
708:       PetscMalloc1(mat->lfwork,&mat->fwork);
709:       PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
710:     }
711:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
712:     PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
713:     PetscFPTrapPop();
714: #endif
715:   } else { /* symmetric case */
716:     if (!mat->pivots) {
717:       PetscMalloc1(A->rmap->n,&mat->pivots);
718:       PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
719:     }
720:     if (!mat->fwork) {
721:       PetscScalar dummy;

723:       mat->lfwork = -1;
724:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
725:       PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
726:       PetscFPTrapPop();
727:       mat->lfwork = (PetscInt)PetscRealPart(dummy);
728:       PetscMalloc1(mat->lfwork,&mat->fwork);
729:       PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
730:     }
731:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
732:     PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
733:     PetscFPTrapPop();
734:   }
735:   if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1);

737:   A->ops->solve             = MatSolve_SeqDense;
738:   A->ops->matsolve          = MatMatSolve_SeqDense;
739:   A->ops->solvetranspose    = MatSolveTranspose_SeqDense;
740:   A->ops->solveadd          = MatSolveAdd_SeqDense;
741:   A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
742:   A->factortype             = MAT_FACTOR_CHOLESKY;

744:   PetscFree(A->solvertype);
745:   PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);

747:   PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
748: #endif
749:   return(0);
750: }


753: PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
754: {
756:   MatFactorInfo  info;

759:   info.fill = 1.0;

761:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
762:   MatCholeskyFactor_SeqDense(fact,0,&info);
763:   return(0);
764: }

766: static PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
767: {
769:   fact->assembled                  = PETSC_TRUE;
770:   fact->preallocated               = PETSC_TRUE;
771:   fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
772:   fact->ops->solve                 = MatSolve_SeqDense;
773:   fact->ops->matsolve              = MatMatSolve_SeqDense;
774:   fact->ops->solvetranspose        = MatSolveTranspose_SeqDense;
775:   fact->ops->solveadd              = MatSolveAdd_SeqDense;
776:   fact->ops->solvetransposeadd     = MatSolveTransposeAdd_SeqDense;
777:   return(0);
778: }

780: static PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info)
781: {
783:   fact->preallocated           = PETSC_TRUE;
784:   fact->assembled              = PETSC_TRUE;
785:   fact->ops->lufactornumeric   = MatLUFactorNumeric_SeqDense;
786:   fact->ops->solve             = MatSolve_SeqDense;
787:   fact->ops->matsolve          = MatMatSolve_SeqDense;
788:   fact->ops->solvetranspose    = MatSolveTranspose_SeqDense;
789:   fact->ops->solveadd          = MatSolveAdd_SeqDense;
790:   fact->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
791:   return(0);
792: }

794: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact)
795: {

799:   MatCreate(PetscObjectComm((PetscObject)A),fact);
800:   MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
801:   MatSetType(*fact,((PetscObject)A)->type_name);
802:   if (ftype == MAT_FACTOR_LU) {
803:     (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
804:   } else {
805:     (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
806:   }
807:   (*fact)->factortype = ftype;

809:   PetscFree((*fact)->solvertype);
810:   PetscStrallocpy(MATSOLVERPETSC,&(*fact)->solvertype);
811:   return(0);
812: }

814: /* ------------------------------------------------------------------*/
815: static PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx)
816: {
817:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
818:   PetscScalar       *x,*v = mat->v,zero = 0.0,xt;
819:   const PetscScalar *b;
820:   PetscErrorCode    ierr;
821:   PetscInt          m = A->rmap->n,i;
822:   PetscBLASInt      o = 1,bm;

825:   if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */
826:   PetscBLASIntCast(m,&bm);
827:   if (flag & SOR_ZERO_INITIAL_GUESS) {
828:     /* this is a hack fix, should have another version without the second BLASdotu */
829:     VecSet(xx,zero);
830:   }
831:   VecGetArray(xx,&x);
832:   VecGetArrayRead(bb,&b);
833:   its  = its*lits;
834:   if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
835:   while (its--) {
836:     if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
837:       for (i=0; i<m; i++) {
838:         PetscStackCallBLAS("BLASdotu",xt   = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o));
839:         x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
840:       }
841:     }
842:     if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
843:       for (i=m-1; i>=0; i--) {
844:         PetscStackCallBLAS("BLASdotu",xt   = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o));
845:         x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
846:       }
847:     }
848:   }
849:   VecRestoreArrayRead(bb,&b);
850:   VecRestoreArray(xx,&x);
851:   return(0);
852: }

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

865:   PetscBLASIntCast(A->rmap->n,&m);
866:   PetscBLASIntCast(A->cmap->n,&n);
867:   VecGetArrayRead(xx,&x);
868:   VecGetArray(yy,&y);
869:   if (!A->rmap->n || !A->cmap->n) {
870:     PetscBLASInt i;
871:     for (i=0; i<n; i++) y[i] = 0.0;
872:   } else {
873:     PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One));
874:     PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);
875:   }
876:   VecRestoreArrayRead(xx,&x);
877:   VecRestoreArray(yy,&y);
878:   return(0);
879: }

881: PETSC_INTERN PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy)
882: {
883:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
884:   PetscScalar       *y,_DOne=1.0,_DZero=0.0;
885:   PetscErrorCode    ierr;
886:   PetscBLASInt      m, n, _One=1;
887:   const PetscScalar *v = mat->v,*x;

890:   PetscBLASIntCast(A->rmap->n,&m);
891:   PetscBLASIntCast(A->cmap->n,&n);
892:   VecGetArrayRead(xx,&x);
893:   VecGetArray(yy,&y);
894:   if (!A->rmap->n || !A->cmap->n) {
895:     PetscBLASInt i;
896:     for (i=0; i<m; i++) y[i] = 0.0;
897:   } else {
898:     PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One));
899:     PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);
900:   }
901:   VecRestoreArrayRead(xx,&x);
902:   VecRestoreArray(yy,&y);
903:   return(0);
904: }

906: PETSC_INTERN PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
907: {
908:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
909:   const PetscScalar *v = mat->v,*x;
910:   PetscScalar       *y,_DOne=1.0;
911:   PetscErrorCode    ierr;
912:   PetscBLASInt      m, n, _One=1;

915:   PetscBLASIntCast(A->rmap->n,&m);
916:   PetscBLASIntCast(A->cmap->n,&n);
917:   if (!A->rmap->n || !A->cmap->n) return(0);
918:   if (zz != yy) {VecCopy(zz,yy);}
919:   VecGetArrayRead(xx,&x);
920:   VecGetArray(yy,&y);
921:   PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
922:   VecRestoreArrayRead(xx,&x);
923:   VecRestoreArray(yy,&y);
924:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
925:   return(0);
926: }

928: static PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
929: {
930:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
931:   const PetscScalar *v = mat->v,*x;
932:   PetscScalar       *y;
933:   PetscErrorCode    ierr;
934:   PetscBLASInt      m, n, _One=1;
935:   PetscScalar       _DOne=1.0;

938:   PetscBLASIntCast(A->rmap->n,&m);
939:   PetscBLASIntCast(A->cmap->n,&n);
940:   if (!A->rmap->n || !A->cmap->n) return(0);
941:   if (zz != yy) {VecCopy(zz,yy);}
942:   VecGetArrayRead(xx,&x);
943:   VecGetArray(yy,&y);
944:   PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
945:   VecRestoreArrayRead(xx,&x);
946:   VecRestoreArray(yy,&y);
947:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
948:   return(0);
949: }

951: /* -----------------------------------------------------------------*/
952: static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
953: {
954:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
955:   PetscScalar    *v;
957:   PetscInt       i;

960:   *ncols = A->cmap->n;
961:   if (cols) {
962:     PetscMalloc1(A->cmap->n+1,cols);
963:     for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
964:   }
965:   if (vals) {
966:     PetscMalloc1(A->cmap->n+1,vals);
967:     v    = mat->v + row;
968:     for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
969:   }
970:   return(0);
971: }

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

978:   if (cols) {PetscFree(*cols);}
979:   if (vals) {PetscFree(*vals); }
980:   return(0);
981: }
982: /* ----------------------------------------------------------------*/
983: static PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
984: {
985:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
986:   PetscInt     i,j,idx=0;

989:   if (!mat->roworiented) {
990:     if (addv == INSERT_VALUES) {
991:       for (j=0; j<n; j++) {
992:         if (indexn[j] < 0) {idx += m; continue;}
993: #if defined(PETSC_USE_DEBUG)
994:         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);
995: #endif
996:         for (i=0; i<m; i++) {
997:           if (indexm[i] < 0) {idx++; continue;}
998: #if defined(PETSC_USE_DEBUG)
999:           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);
1000: #endif
1001:           mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
1002:         }
1003:       }
1004:     } else {
1005:       for (j=0; j<n; j++) {
1006:         if (indexn[j] < 0) {idx += m; continue;}
1007: #if defined(PETSC_USE_DEBUG)
1008:         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);
1009: #endif
1010:         for (i=0; i<m; i++) {
1011:           if (indexm[i] < 0) {idx++; continue;}
1012: #if defined(PETSC_USE_DEBUG)
1013:           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);
1014: #endif
1015:           mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
1016:         }
1017:       }
1018:     }
1019:   } else {
1020:     if (addv == INSERT_VALUES) {
1021:       for (i=0; i<m; i++) {
1022:         if (indexm[i] < 0) { idx += n; continue;}
1023: #if defined(PETSC_USE_DEBUG)
1024:         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);
1025: #endif
1026:         for (j=0; j<n; j++) {
1027:           if (indexn[j] < 0) { idx++; continue;}
1028: #if defined(PETSC_USE_DEBUG)
1029:           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);
1030: #endif
1031:           mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
1032:         }
1033:       }
1034:     } else {
1035:       for (i=0; i<m; i++) {
1036:         if (indexm[i] < 0) { idx += n; continue;}
1037: #if defined(PETSC_USE_DEBUG)
1038:         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);
1039: #endif
1040:         for (j=0; j<n; j++) {
1041:           if (indexn[j] < 0) { idx++; continue;}
1042: #if defined(PETSC_USE_DEBUG)
1043:           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);
1044: #endif
1045:           mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
1046:         }
1047:       }
1048:     }
1049:   }
1050:   return(0);
1051: }

1053: static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
1054: {
1055:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1056:   PetscInt     i,j;

1059:   /* row-oriented output */
1060:   for (i=0; i<m; i++) {
1061:     if (indexm[i] < 0) {v += n;continue;}
1062:     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);
1063:     for (j=0; j<n; j++) {
1064:       if (indexn[j] < 0) {v++; continue;}
1065:       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);
1066:       *v++ = mat->v[indexn[j]*mat->lda + indexm[i]];
1067:     }
1068:   }
1069:   return(0);
1070: }

1072: /* -----------------------------------------------------------------*/

1074: static PetscErrorCode MatLoad_SeqDense(Mat newmat,PetscViewer viewer)
1075: {
1076:   Mat_SeqDense   *a;
1078:   PetscInt       *scols,i,j,nz,header[4];
1079:   int            fd;
1080:   PetscMPIInt    size;
1081:   PetscInt       *rowlengths = 0,M,N,*cols,grows,gcols;
1082:   PetscScalar    *vals,*svals,*v,*w;
1083:   MPI_Comm       comm;

1086:   /* force binary viewer to load .info file if it has not yet done so */
1087:   PetscViewerSetUp(viewer);
1088:   PetscObjectGetComm((PetscObject)viewer,&comm);
1089:   MPI_Comm_size(comm,&size);
1090:   if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor");
1091:   PetscViewerBinaryGetDescriptor(viewer,&fd);
1092:   PetscBinaryRead(fd,header,4,PETSC_INT);
1093:   if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not matrix object");
1094:   M = header[1]; N = header[2]; nz = header[3];

1096:   /* set global size if not set already*/
1097:   if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) {
1098:     MatSetSizes(newmat,M,N,M,N);
1099:   } else {
1100:     /* if sizes and type are already set, check if the vector global sizes are correct */
1101:     MatGetSize(newmat,&grows,&gcols);
1102:     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);
1103:   }
1104:   a = (Mat_SeqDense*)newmat->data;
1105:   if (!a->user_alloc) {
1106:     MatSeqDenseSetPreallocation(newmat,NULL);
1107:   }

1109:   if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */
1110:     a = (Mat_SeqDense*)newmat->data;
1111:     v = a->v;
1112:     /* Allocate some temp space to read in the values and then flip them
1113:        from row major to column major */
1114:     PetscMalloc1(M*N > 0 ? M*N : 1,&w);
1115:     /* read in nonzero values */
1116:     PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);
1117:     /* now flip the values and store them in the matrix*/
1118:     for (j=0; j<N; j++) {
1119:       for (i=0; i<M; i++) {
1120:         *v++ =w[i*N+j];
1121:       }
1122:     }
1123:     PetscFree(w);
1124:   } else {
1125:     /* read row lengths */
1126:     PetscMalloc1(M+1,&rowlengths);
1127:     PetscBinaryRead(fd,rowlengths,M,PETSC_INT);

1129:     a = (Mat_SeqDense*)newmat->data;
1130:     v = a->v;

1132:     /* read column indices and nonzeros */
1133:     PetscMalloc1(nz+1,&scols);
1134:     cols = scols;
1135:     PetscBinaryRead(fd,cols,nz,PETSC_INT);
1136:     PetscMalloc1(nz+1,&svals);
1137:     vals = svals;
1138:     PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);

1140:     /* insert into matrix */
1141:     for (i=0; i<M; i++) {
1142:       for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j];
1143:       svals += rowlengths[i]; scols += rowlengths[i];
1144:     }
1145:     PetscFree(vals);
1146:     PetscFree(cols);
1147:     PetscFree(rowlengths);
1148:   }
1149:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1150:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
1151:   return(0);
1152: }

1154: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
1155: {
1156:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
1157:   PetscErrorCode    ierr;
1158:   PetscInt          i,j;
1159:   const char        *name;
1160:   PetscScalar       *v;
1161:   PetscViewerFormat format;
1162: #if defined(PETSC_USE_COMPLEX)
1163:   PetscBool         allreal = PETSC_TRUE;
1164: #endif

1167:   PetscViewerGetFormat(viewer,&format);
1168:   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1169:     return(0);  /* do nothing for now */
1170:   } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1171:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1172:     for (i=0; i<A->rmap->n; i++) {
1173:       v    = a->v + i;
1174:       PetscViewerASCIIPrintf(viewer,"row %D:",i);
1175:       for (j=0; j<A->cmap->n; j++) {
1176: #if defined(PETSC_USE_COMPLEX)
1177:         if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
1178:           PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1179:         } else if (PetscRealPart(*v)) {
1180:           PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));
1181:         }
1182: #else
1183:         if (*v) {
1184:           PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);
1185:         }
1186: #endif
1187:         v += a->lda;
1188:       }
1189:       PetscViewerASCIIPrintf(viewer,"\n");
1190:     }
1191:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1192:   } else {
1193:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1194: #if defined(PETSC_USE_COMPLEX)
1195:     /* determine if matrix has all real values */
1196:     v = a->v;
1197:     for (i=0; i<A->rmap->n*A->cmap->n; i++) {
1198:       if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;}
1199:     }
1200: #endif
1201:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
1202:       PetscObjectGetName((PetscObject)A,&name);
1203:       PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);
1204:       PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);
1205:       PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
1206:     }

1208:     for (i=0; i<A->rmap->n; i++) {
1209:       v = a->v + i;
1210:       for (j=0; j<A->cmap->n; j++) {
1211: #if defined(PETSC_USE_COMPLEX)
1212:         if (allreal) {
1213:           PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));
1214:         } else {
1215:           PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1216:         }
1217: #else
1218:         PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);
1219: #endif
1220:         v += a->lda;
1221:       }
1222:       PetscViewerASCIIPrintf(viewer,"\n");
1223:     }
1224:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
1225:       PetscViewerASCIIPrintf(viewer,"];\n");
1226:     }
1227:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1228:   }
1229:   PetscViewerFlush(viewer);
1230:   return(0);
1231: }

1233: static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer)
1234: {
1235:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
1236:   PetscErrorCode    ierr;
1237:   int               fd;
1238:   PetscInt          ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n;
1239:   PetscScalar       *v,*anonz,*vals;
1240:   PetscViewerFormat format;

1243:   PetscViewerBinaryGetDescriptor(viewer,&fd);

1245:   PetscViewerGetFormat(viewer,&format);
1246:   if (format == PETSC_VIEWER_NATIVE) {
1247:     /* store the matrix as a dense matrix */
1248:     PetscMalloc1(4,&col_lens);

1250:     col_lens[0] = MAT_FILE_CLASSID;
1251:     col_lens[1] = m;
1252:     col_lens[2] = n;
1253:     col_lens[3] = MATRIX_BINARY_FORMAT_DENSE;

1255:     PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);
1256:     PetscFree(col_lens);

1258:     /* write out matrix, by rows */
1259:     PetscMalloc1(m*n+1,&vals);
1260:     v    = a->v;
1261:     for (j=0; j<n; j++) {
1262:       for (i=0; i<m; i++) {
1263:         vals[j + i*n] = *v++;
1264:       }
1265:     }
1266:     PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);
1267:     PetscFree(vals);
1268:   } else {
1269:     PetscMalloc1(4+nz,&col_lens);

1271:     col_lens[0] = MAT_FILE_CLASSID;
1272:     col_lens[1] = m;
1273:     col_lens[2] = n;
1274:     col_lens[3] = nz;

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

1280:     /* Possibly should write in smaller increments, not whole matrix at once? */
1281:     /* store column indices (zero start index) */
1282:     ict = 0;
1283:     for (i=0; i<m; i++) {
1284:       for (j=0; j<n; j++) col_lens[ict++] = j;
1285:     }
1286:     PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);
1287:     PetscFree(col_lens);

1289:     /* store nonzero values */
1290:     PetscMalloc1(nz+1,&anonz);
1291:     ict  = 0;
1292:     for (i=0; i<m; i++) {
1293:       v = a->v + i;
1294:       for (j=0; j<n; j++) {
1295:         anonz[ict++] = *v; v += a->lda;
1296:       }
1297:     }
1298:     PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);
1299:     PetscFree(anonz);
1300:   }
1301:   return(0);
1302: }

1304:  #include <petscdraw.h>
1305: static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa)
1306: {
1307:   Mat               A  = (Mat) Aa;
1308:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
1309:   PetscErrorCode    ierr;
1310:   PetscInt          m  = A->rmap->n,n = A->cmap->n,i,j;
1311:   int               color = PETSC_DRAW_WHITE;
1312:   PetscScalar       *v = a->v;
1313:   PetscViewer       viewer;
1314:   PetscReal         xl,yl,xr,yr,x_l,x_r,y_l,y_r;
1315:   PetscViewerFormat format;

1318:   PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
1319:   PetscViewerGetFormat(viewer,&format);
1320:   PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);

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

1324:   if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1325:     PetscDrawCollectiveBegin(draw);
1326:     /* Blue for negative and Red for positive */
1327:     for (j = 0; j < n; j++) {
1328:       x_l = j; x_r = x_l + 1.0;
1329:       for (i = 0; i < m; i++) {
1330:         y_l = m - i - 1.0;
1331:         y_r = y_l + 1.0;
1332:         if (PetscRealPart(v[j*m+i]) >  0.) {
1333:           color = PETSC_DRAW_RED;
1334:         } else if (PetscRealPart(v[j*m+i]) <  0.) {
1335:           color = PETSC_DRAW_BLUE;
1336:         } else {
1337:           continue;
1338:         }
1339:         PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1340:       }
1341:     }
1342:     PetscDrawCollectiveEnd(draw);
1343:   } else {
1344:     /* use contour shading to indicate magnitude of values */
1345:     /* first determine max of all nonzero values */
1346:     PetscReal minv = 0.0, maxv = 0.0;
1347:     PetscDraw popup;

1349:     for (i=0; i < m*n; i++) {
1350:       if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1351:     }
1352:     if (minv >= maxv) maxv = minv + PETSC_SMALL;
1353:     PetscDrawGetPopup(draw,&popup);
1354:     PetscDrawScalePopup(popup,minv,maxv);

1356:     PetscDrawCollectiveBegin(draw);
1357:     for (j=0; j<n; j++) {
1358:       x_l = j;
1359:       x_r = x_l + 1.0;
1360:       for (i=0; i<m; i++) {
1361:         y_l = m - i - 1.0;
1362:         y_r = y_l + 1.0;
1363:         color = PetscDrawRealToColor(PetscAbsScalar(v[j*m+i]),minv,maxv);
1364:         PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1365:       }
1366:     }
1367:     PetscDrawCollectiveEnd(draw);
1368:   }
1369:   return(0);
1370: }

1372: static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
1373: {
1374:   PetscDraw      draw;
1375:   PetscBool      isnull;
1376:   PetscReal      xr,yr,xl,yl,h,w;

1380:   PetscViewerDrawGetDraw(viewer,0,&draw);
1381:   PetscDrawIsNull(draw,&isnull);
1382:   if (isnull) return(0);

1384:   xr   = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
1385:   xr  += w;          yr += h;        xl = -w;     yl = -h;
1386:   PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
1387:   PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
1388:   PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
1389:   PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
1390:   PetscDrawSave(draw);
1391:   return(0);
1392: }

1394: PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
1395: {
1397:   PetscBool      iascii,isbinary,isdraw;

1400:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1401:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1402:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);

1404:   if (iascii) {
1405:     MatView_SeqDense_ASCII(A,viewer);
1406:   } else if (isbinary) {
1407:     MatView_SeqDense_Binary(A,viewer);
1408:   } else if (isdraw) {
1409:     MatView_SeqDense_Draw(A,viewer);
1410:   }
1411:   return(0);
1412: }

1414: static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A,const PetscScalar array[])
1415: {
1416:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

1419:   a->unplacedarray       = a->v;
1420:   a->unplaced_user_alloc = a->user_alloc;
1421:   a->v                   = (PetscScalar*) array;
1422:   return(0);
1423: }

1425: static PetscErrorCode MatDenseResetArray_SeqDense(Mat A)
1426: {
1427:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

1430:   a->v             = a->unplacedarray;
1431:   a->user_alloc    = a->unplaced_user_alloc;
1432:   a->unplacedarray = NULL;
1433:   return(0);
1434: }

1436: static PetscErrorCode MatDestroy_SeqDense(Mat mat)
1437: {
1438:   Mat_SeqDense   *l = (Mat_SeqDense*)mat->data;

1442: #if defined(PETSC_USE_LOG)
1443:   PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n);
1444: #endif
1445:   PetscFree(l->pivots);
1446:   PetscFree(l->fwork);
1447:   MatDestroy(&l->ptapwork);
1448:   if (!l->user_alloc) {PetscFree(l->v);}
1449:   PetscFree(mat->data);

1451:   PetscObjectChangeTypeName((PetscObject)mat,0);
1452:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);
1453:   PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL);
1454:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL);
1455:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);
1456:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);
1457: #if defined(PETSC_HAVE_ELEMENTAL)
1458:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);
1459: #endif
1460:   PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);
1461:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);
1462:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);
1463:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);
1464:   PetscObjectComposeFunction((PetscObject)mat,"MatPtAP_seqaij_seqdense_C",NULL);
1465:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);
1466:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);
1467:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);
1468:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumn_C",NULL);
1469:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumn_C",NULL);
1470:   return(0);
1471: }

1473: static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
1474: {
1475:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1477:   PetscInt       k,j,m,n,M;
1478:   PetscScalar    *v,tmp;

1481:   v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n;
1482:   if (reuse == MAT_INPLACE_MATRIX) { /* in place transpose */
1483:     if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can not transpose non-square matrix in place");
1484:     else {
1485:       for (j=0; j<m; j++) {
1486:         for (k=0; k<j; k++) {
1487:           tmp        = v[j + k*M];
1488:           v[j + k*M] = v[k + j*M];
1489:           v[k + j*M] = tmp;
1490:         }
1491:       }
1492:     }
1493:   } else { /* out-of-place transpose */
1494:     Mat          tmat;
1495:     Mat_SeqDense *tmatd;
1496:     PetscScalar  *v2;
1497:     PetscInt     M2;

1499:     if (reuse == MAT_INITIAL_MATRIX) {
1500:       MatCreate(PetscObjectComm((PetscObject)A),&tmat);
1501:       MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);
1502:       MatSetType(tmat,((PetscObject)A)->type_name);
1503:       MatSeqDenseSetPreallocation(tmat,NULL);
1504:     } else {
1505:       tmat = *matout;
1506:     }
1507:     tmatd = (Mat_SeqDense*)tmat->data;
1508:     v = mat->v; v2 = tmatd->v; M2 = tmatd->lda;
1509:     for (j=0; j<n; j++) {
1510:       for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M];
1511:     }
1512:     MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
1513:     MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
1514:     *matout = tmat;
1515:   }
1516:   return(0);
1517: }

1519: static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool  *flg)
1520: {
1521:   Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data;
1522:   Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data;
1523:   PetscInt     i,j;
1524:   PetscScalar  *v1,*v2;

1527:   if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; return(0);}
1528:   if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; return(0);}
1529:   for (i=0; i<A1->rmap->n; i++) {
1530:     v1 = mat1->v+i; v2 = mat2->v+i;
1531:     for (j=0; j<A1->cmap->n; j++) {
1532:       if (*v1 != *v2) {*flg = PETSC_FALSE; return(0);}
1533:       v1 += mat1->lda; v2 += mat2->lda;
1534:     }
1535:   }
1536:   *flg = PETSC_TRUE;
1537:   return(0);
1538: }

1540: static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
1541: {
1542:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1544:   PetscInt       i,n,len;
1545:   PetscScalar    *x,zero = 0.0;

1548:   VecSet(v,zero);
1549:   VecGetSize(v,&n);
1550:   VecGetArray(v,&x);
1551:   len  = PetscMin(A->rmap->n,A->cmap->n);
1552:   if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
1553:   for (i=0; i<len; i++) {
1554:     x[i] = mat->v[i*mat->lda + i];
1555:   }
1556:   VecRestoreArray(v,&x);
1557:   return(0);
1558: }

1560: static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
1561: {
1562:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1563:   const PetscScalar *l,*r;
1564:   PetscScalar       x,*v;
1565:   PetscErrorCode    ierr;
1566:   PetscInt          i,j,m = A->rmap->n,n = A->cmap->n;

1569:   if (ll) {
1570:     VecGetSize(ll,&m);
1571:     VecGetArrayRead(ll,&l);
1572:     if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size");
1573:     for (i=0; i<m; i++) {
1574:       x = l[i];
1575:       v = mat->v + i;
1576:       for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;}
1577:     }
1578:     VecRestoreArrayRead(ll,&l);
1579:     PetscLogFlops(1.0*n*m);
1580:   }
1581:   if (rr) {
1582:     VecGetSize(rr,&n);
1583:     VecGetArrayRead(rr,&r);
1584:     if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size");
1585:     for (i=0; i<n; i++) {
1586:       x = r[i];
1587:       v = mat->v + i*mat->lda;
1588:       for (j=0; j<m; j++) (*v++) *= x;
1589:     }
1590:     VecRestoreArrayRead(rr,&r);
1591:     PetscLogFlops(1.0*n*m);
1592:   }
1593:   return(0);
1594: }

1596: static PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1597: {
1598:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1599:   PetscScalar    *v   = mat->v;
1600:   PetscReal      sum  = 0.0;
1601:   PetscInt       lda  =mat->lda,m=A->rmap->n,i,j;

1605:   if (type == NORM_FROBENIUS) {
1606:     if (lda>m) {
1607:       for (j=0; j<A->cmap->n; j++) {
1608:         v = mat->v+j*lda;
1609:         for (i=0; i<m; i++) {
1610:           sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1611:         }
1612:       }
1613:     } else {
1614: #if defined(PETSC_USE_REAL___FP16)
1615:       PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n;
1616:       *nrm = BLASnrm2_(&cnt,v,&one);
1617:     }
1618: #else
1619:       for (i=0; i<A->cmap->n*A->rmap->n; i++) {
1620:         sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1621:       }
1622:     }
1623:     *nrm = PetscSqrtReal(sum);
1624: #endif
1625:     PetscLogFlops(2.0*A->cmap->n*A->rmap->n);
1626:   } else if (type == NORM_1) {
1627:     *nrm = 0.0;
1628:     for (j=0; j<A->cmap->n; j++) {
1629:       v   = mat->v + j*mat->lda;
1630:       sum = 0.0;
1631:       for (i=0; i<A->rmap->n; i++) {
1632:         sum += PetscAbsScalar(*v);  v++;
1633:       }
1634:       if (sum > *nrm) *nrm = sum;
1635:     }
1636:     PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1637:   } else if (type == NORM_INFINITY) {
1638:     *nrm = 0.0;
1639:     for (j=0; j<A->rmap->n; j++) {
1640:       v   = mat->v + j;
1641:       sum = 0.0;
1642:       for (i=0; i<A->cmap->n; i++) {
1643:         sum += PetscAbsScalar(*v); v += mat->lda;
1644:       }
1645:       if (sum > *nrm) *nrm = sum;
1646:     }
1647:     PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1648:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
1649:   return(0);
1650: }

1652: static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
1653: {
1654:   Mat_SeqDense   *aij = (Mat_SeqDense*)A->data;

1658:   switch (op) {
1659:   case MAT_ROW_ORIENTED:
1660:     aij->roworiented = flg;
1661:     break;
1662:   case MAT_NEW_NONZERO_LOCATIONS:
1663:   case MAT_NEW_NONZERO_LOCATION_ERR:
1664:   case MAT_NEW_NONZERO_ALLOCATION_ERR:
1665:   case MAT_NEW_DIAGONALS:
1666:   case MAT_KEEP_NONZERO_PATTERN:
1667:   case MAT_IGNORE_OFF_PROC_ENTRIES:
1668:   case MAT_USE_HASH_TABLE:
1669:   case MAT_IGNORE_ZERO_ENTRIES:
1670:   case MAT_IGNORE_LOWER_TRIANGULAR:
1671:     PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
1672:     break;
1673:   case MAT_SPD:
1674:   case MAT_SYMMETRIC:
1675:   case MAT_STRUCTURALLY_SYMMETRIC:
1676:   case MAT_HERMITIAN:
1677:   case MAT_SYMMETRY_ETERNAL:
1678:     /* These options are handled directly by MatSetOption() */
1679:     break;
1680:   default:
1681:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]);
1682:   }
1683:   return(0);
1684: }

1686: static PetscErrorCode MatZeroEntries_SeqDense(Mat A)
1687: {
1688:   Mat_SeqDense   *l = (Mat_SeqDense*)A->data;
1690:   PetscInt       lda=l->lda,m=A->rmap->n,j;

1693:   if (lda>m) {
1694:     for (j=0; j<A->cmap->n; j++) {
1695:       PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));
1696:     }
1697:   } else {
1698:     PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1699:   }
1700:   return(0);
1701: }

1703: static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1704: {
1705:   PetscErrorCode    ierr;
1706:   Mat_SeqDense      *l = (Mat_SeqDense*)A->data;
1707:   PetscInt          m  = l->lda, n = A->cmap->n, i,j;
1708:   PetscScalar       *slot,*bb;
1709:   const PetscScalar *xx;

1712: #if defined(PETSC_USE_DEBUG)
1713:   for (i=0; i<N; i++) {
1714:     if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
1715:     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);
1716:   }
1717: #endif

1719:   /* fix right hand side if needed */
1720:   if (x && b) {
1721:     VecGetArrayRead(x,&xx);
1722:     VecGetArray(b,&bb);
1723:     for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
1724:     VecRestoreArrayRead(x,&xx);
1725:     VecRestoreArray(b,&bb);
1726:   }

1728:   for (i=0; i<N; i++) {
1729:     slot = l->v + rows[i];
1730:     for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
1731:   }
1732:   if (diag != 0.0) {
1733:     if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
1734:     for (i=0; i<N; i++) {
1735:       slot  = l->v + (m+1)*rows[i];
1736:       *slot = diag;
1737:     }
1738:   }
1739:   return(0);
1740: }

1742: static PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[])
1743: {
1744:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;

1747:   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");
1748:   *array = mat->v;
1749:   return(0);
1750: }

1752: static PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[])
1753: {
1755:   return(0);
1756: }

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

1761:    Logically Collective on Mat

1763:    Input Parameter:
1764: .  mat - a MATSEQDENSE or MATMPIDENSE matrix

1766:    Output Parameter:
1767: .   array - pointer to the data

1769:    Level: intermediate

1771: .seealso: MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead()
1772: @*/
1773: PetscErrorCode  MatDenseGetArray(Mat A,PetscScalar **array)
1774: {

1778:   PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));
1779:   return(0);
1780: }

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

1785:    Logically Collective on Mat

1787:    Input Parameters:
1788: .  mat - a MATSEQDENSE or MATMPIDENSE matrix
1789: .  array - pointer to the data

1791:    Level: intermediate

1793: .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead()
1794: @*/
1795: PetscErrorCode  MatDenseRestoreArray(Mat A,PetscScalar **array)
1796: {

1800:   PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));
1801:   if (array) *array = NULL;
1802:   PetscObjectStateIncrease((PetscObject)A);
1803:   return(0);
1804: }

1806: /*@C
1807:    MatDenseGetArrayRead - gives access to the array where the data for a SeqDense matrix is stored

1809:    Not Collective

1811:    Input Parameter:
1812: .  mat - a MATSEQDENSE or MATMPIDENSE matrix

1814:    Output Parameter:
1815: .   array - pointer to the data

1817:    Level: intermediate

1819: .seealso: MatDenseRestoreArray(), MatDenseGetArray(), MatDenseRestoreArrayRead()
1820: @*/
1821: PetscErrorCode  MatDenseGetArrayRead(Mat A,const PetscScalar **array)
1822: {

1826:   PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array));
1827:   return(0);
1828: }

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

1833:    Not Collective

1835:    Input Parameters:
1836: .  mat - a MATSEQDENSE or MATMPIDENSE matrix
1837: .  array - pointer to the data

1839:    Level: intermediate

1841: .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArray()
1842: @*/
1843: PetscErrorCode  MatDenseRestoreArrayRead(Mat A,const PetscScalar **array)
1844: {

1848:   PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array));
1849:   if (array) *array = NULL;
1850:   return(0);
1851: }

1853: static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B)
1854: {
1855:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1857:   PetscInt       i,j,nrows,ncols;
1858:   const PetscInt *irow,*icol;
1859:   PetscScalar    *av,*bv,*v = mat->v;
1860:   Mat            newmat;

1863:   ISGetIndices(isrow,&irow);
1864:   ISGetIndices(iscol,&icol);
1865:   ISGetLocalSize(isrow,&nrows);
1866:   ISGetLocalSize(iscol,&ncols);

1868:   /* Check submatrixcall */
1869:   if (scall == MAT_REUSE_MATRIX) {
1870:     PetscInt n_cols,n_rows;
1871:     MatGetSize(*B,&n_rows,&n_cols);
1872:     if (n_rows != nrows || n_cols != ncols) {
1873:       /* resize the result matrix to match number of requested rows/columns */
1874:       MatSetSizes(*B,nrows,ncols,nrows,ncols);
1875:     }
1876:     newmat = *B;
1877:   } else {
1878:     /* Create and fill new matrix */
1879:     MatCreate(PetscObjectComm((PetscObject)A),&newmat);
1880:     MatSetSizes(newmat,nrows,ncols,nrows,ncols);
1881:     MatSetType(newmat,((PetscObject)A)->type_name);
1882:     MatSeqDenseSetPreallocation(newmat,NULL);
1883:   }

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

1888:   for (i=0; i<ncols; i++) {
1889:     av = v + mat->lda*icol[i];
1890:     for (j=0; j<nrows; j++) *bv++ = av[irow[j]];
1891:   }

1893:   /* Assemble the matrices so that the correct flags are set */
1894:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1895:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);

1897:   /* Free work space */
1898:   ISRestoreIndices(isrow,&irow);
1899:   ISRestoreIndices(iscol,&icol);
1900:   *B   = newmat;
1901:   return(0);
1902: }

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

1910:   if (scall == MAT_INITIAL_MATRIX) {
1911:     PetscCalloc1(n+1,B);
1912:   }

1914:   for (i=0; i<n; i++) {
1915:     MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
1916:   }
1917:   return(0);
1918: }

1920: static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
1921: {
1923:   return(0);
1924: }

1926: static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
1927: {
1929:   return(0);
1930: }

1932: static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
1933: {
1934:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data;
1936:   PetscInt       lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j;

1939:   /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
1940:   if (A->ops->copy != B->ops->copy) {
1941:     MatCopy_Basic(A,B,str);
1942:     return(0);
1943:   }
1944:   if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)");
1945:   if (lda1>m || lda2>m) {
1946:     for (j=0; j<n; j++) {
1947:       PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));
1948:     }
1949:   } else {
1950:     PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1951:   }
1952:   PetscObjectStateIncrease((PetscObject)B);
1953:   return(0);
1954: }

1956: static PetscErrorCode MatSetUp_SeqDense(Mat A)
1957: {

1961:    MatSeqDenseSetPreallocation(A,0);
1962:   return(0);
1963: }

1965: static PetscErrorCode MatConjugate_SeqDense(Mat A)
1966: {
1967:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1968:   PetscInt     i,nz = A->rmap->n*A->cmap->n;
1969:   PetscScalar  *aa = a->v;

1972:   for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
1973:   return(0);
1974: }

1976: static PetscErrorCode MatRealPart_SeqDense(Mat A)
1977: {
1978:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1979:   PetscInt     i,nz = A->rmap->n*A->cmap->n;
1980:   PetscScalar  *aa = a->v;

1983:   for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
1984:   return(0);
1985: }

1987: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
1988: {
1989:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1990:   PetscInt     i,nz = A->rmap->n*A->cmap->n;
1991:   PetscScalar  *aa = a->v;

1994:   for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
1995:   return(0);
1996: }

1998: /* ----------------------------------------------------------------*/
1999: PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
2000: {

2004:   if (scall == MAT_INITIAL_MATRIX) {
2005:     PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
2006:     MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
2007:     PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
2008:   }
2009:   PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
2010:   MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);
2011:   PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
2012:   return(0);
2013: }

2015: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
2016: {
2018:   PetscInt       m=A->rmap->n,n=B->cmap->n;
2019:   Mat            Cmat;

2022:   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);
2023:   MatCreate(PETSC_COMM_SELF,&Cmat);
2024:   MatSetSizes(Cmat,m,n,m,n);
2025:   MatSetType(Cmat,MATSEQDENSE);
2026:   MatSeqDenseSetPreallocation(Cmat,NULL);

2028:   *C = Cmat;
2029:   return(0);
2030: }

2032: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2033: {
2034:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2035:   Mat_SeqDense   *b = (Mat_SeqDense*)B->data;
2036:   Mat_SeqDense   *c = (Mat_SeqDense*)C->data;
2037:   PetscBLASInt   m,n,k;
2038:   PetscScalar    _DOne=1.0,_DZero=0.0;
2040:   PetscBool      flg;

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

2046:   /* Handle case where where user provided the final C matrix rather than calling MatMatMult() with MAT_INITIAL_MATRIX*/
2047:   PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
2048:   if (flg) {
2049:     C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqDense;
2050:     (*C->ops->matmultnumeric)(A,B,C);
2051:     return(0);
2052:   }

2054:   PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&flg);
2055:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"First matrix must be dense");
2056:   PetscBLASIntCast(C->rmap->n,&m);
2057:   PetscBLASIntCast(C->cmap->n,&n);
2058:   PetscBLASIntCast(A->cmap->n,&k);
2059:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
2060:   return(0);
2061: }

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

2068:   if (scall == MAT_INITIAL_MATRIX) {
2069:     PetscLogEventBegin(MAT_MatTransposeMultSymbolic,A,B,0,0);
2070:     MatMatTransposeMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
2071:     PetscLogEventEnd(MAT_MatTransposeMultSymbolic,A,B,0,0);
2072:   }
2073:   PetscLogEventBegin(MAT_MatTransposeMultNumeric,A,B,0,0);
2074:   MatMatTransposeMultNumeric_SeqDense_SeqDense(A,B,*C);
2075:   PetscLogEventEnd(MAT_MatTransposeMultNumeric,A,B,0,0);
2076:   return(0);
2077: }

2079: PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
2080: {
2082:   PetscInt       m=A->rmap->n,n=B->rmap->n;
2083:   Mat            Cmat;

2086:   if (A->cmap->n != B->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"A->cmap->n %d != B->cmap->n %d\n",A->cmap->n,B->cmap->n);
2087:   MatCreate(PETSC_COMM_SELF,&Cmat);
2088:   MatSetSizes(Cmat,m,n,m,n);
2089:   MatSetType(Cmat,MATSEQDENSE);
2090:   MatSeqDenseSetPreallocation(Cmat,NULL);

2092:   Cmat->assembled = PETSC_TRUE;

2094:   *C = Cmat;
2095:   return(0);
2096: }

2098: PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2099: {
2100:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2101:   Mat_SeqDense   *b = (Mat_SeqDense*)B->data;
2102:   Mat_SeqDense   *c = (Mat_SeqDense*)C->data;
2103:   PetscBLASInt   m,n,k;
2104:   PetscScalar    _DOne=1.0,_DZero=0.0;

2108:   PetscBLASIntCast(A->rmap->n,&m);
2109:   PetscBLASIntCast(B->rmap->n,&n);
2110:   PetscBLASIntCast(A->cmap->n,&k);
2111:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
2112:   return(0);
2113: }

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

2120:   if (scall == MAT_INITIAL_MATRIX) {
2121:     PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);
2122:     MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
2123:     PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);
2124:   }
2125:   PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);
2126:   MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);
2127:   PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);
2128:   return(0);
2129: }

2131: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
2132: {
2134:   PetscInt       m=A->cmap->n,n=B->cmap->n;
2135:   Mat            Cmat;

2138:   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);
2139:   MatCreate(PETSC_COMM_SELF,&Cmat);
2140:   MatSetSizes(Cmat,m,n,m,n);
2141:   MatSetType(Cmat,MATSEQDENSE);
2142:   MatSeqDenseSetPreallocation(Cmat,NULL);

2144:   Cmat->assembled = PETSC_TRUE;

2146:   *C = Cmat;
2147:   return(0);
2148: }

2150: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2151: {
2152:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2153:   Mat_SeqDense   *b = (Mat_SeqDense*)B->data;
2154:   Mat_SeqDense   *c = (Mat_SeqDense*)C->data;
2155:   PetscBLASInt   m,n,k;
2156:   PetscScalar    _DOne=1.0,_DZero=0.0;

2160:   PetscBLASIntCast(C->rmap->n,&m);
2161:   PetscBLASIntCast(C->cmap->n,&n);
2162:   PetscBLASIntCast(A->rmap->n,&k);
2163:   PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
2164:   return(0);
2165: }

2167: static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[])
2168: {
2169:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2171:   PetscInt       i,j,m = A->rmap->n,n = A->cmap->n,p;
2172:   PetscScalar    *x;
2173:   MatScalar      *aa = a->v;

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

2178:   VecSet(v,0.0);
2179:   VecGetArray(v,&x);
2180:   VecGetLocalSize(v,&p);
2181:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2182:   for (i=0; i<m; i++) {
2183:     x[i] = aa[i]; if (idx) idx[i] = 0;
2184:     for (j=1; j<n; j++) {
2185:       if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
2186:     }
2187:   }
2188:   VecRestoreArray(v,&x);
2189:   return(0);
2190: }

2192: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])
2193: {
2194:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2196:   PetscInt       i,j,m = A->rmap->n,n = A->cmap->n,p;
2197:   PetscScalar    *x;
2198:   PetscReal      atmp;
2199:   MatScalar      *aa = a->v;

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

2204:   VecSet(v,0.0);
2205:   VecGetArray(v,&x);
2206:   VecGetLocalSize(v,&p);
2207:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2208:   for (i=0; i<m; i++) {
2209:     x[i] = PetscAbsScalar(aa[i]);
2210:     for (j=1; j<n; j++) {
2211:       atmp = PetscAbsScalar(aa[i+m*j]);
2212:       if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;}
2213:     }
2214:   }
2215:   VecRestoreArray(v,&x);
2216:   return(0);
2217: }

2219: static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])
2220: {
2221:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2223:   PetscInt       i,j,m = A->rmap->n,n = A->cmap->n,p;
2224:   PetscScalar    *x;
2225:   MatScalar      *aa = a->v;

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

2230:   VecSet(v,0.0);
2231:   VecGetArray(v,&x);
2232:   VecGetLocalSize(v,&p);
2233:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2234:   for (i=0; i<m; i++) {
2235:     x[i] = aa[i]; if (idx) idx[i] = 0;
2236:     for (j=1; j<n; j++) {
2237:       if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
2238:     }
2239:   }
2240:   VecRestoreArray(v,&x);
2241:   return(0);
2242: }

2244: static PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
2245: {
2246:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2248:   PetscScalar    *x;

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

2253:   VecGetArray(v,&x);
2254:   PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));
2255:   VecRestoreArray(v,&x);
2256:   return(0);
2257: }

2259: PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms)
2260: {
2262:   PetscInt       i,j,m,n;
2263:   PetscScalar    *a;

2266:   MatGetSize(A,&m,&n);
2267:   PetscMemzero(norms,n*sizeof(PetscReal));
2268:   MatDenseGetArray(A,&a);
2269:   if (type == NORM_2) {
2270:     for (i=0; i<n; i++) {
2271:       for (j=0; j<m; j++) {
2272:         norms[i] += PetscAbsScalar(a[j]*a[j]);
2273:       }
2274:       a += m;
2275:     }
2276:   } else if (type == NORM_1) {
2277:     for (i=0; i<n; i++) {
2278:       for (j=0; j<m; j++) {
2279:         norms[i] += PetscAbsScalar(a[j]);
2280:       }
2281:       a += m;
2282:     }
2283:   } else if (type == NORM_INFINITY) {
2284:     for (i=0; i<n; i++) {
2285:       for (j=0; j<m; j++) {
2286:         norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]);
2287:       }
2288:       a += m;
2289:     }
2290:   } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType");
2291:   MatDenseRestoreArray(A,&a);
2292:   if (type == NORM_2) {
2293:     for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]);
2294:   }
2295:   return(0);
2296: }

2298: static PetscErrorCode  MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
2299: {
2301:   PetscScalar    *a;
2302:   PetscInt       m,n,i;

2305:   MatGetSize(x,&m,&n);
2306:   MatDenseGetArray(x,&a);
2307:   for (i=0; i<m*n; i++) {
2308:     PetscRandomGetValue(rctx,a+i);
2309:   }
2310:   MatDenseRestoreArray(x,&a);
2311:   return(0);
2312: }

2314: static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool  *missing,PetscInt *d)
2315: {
2317:   *missing = PETSC_FALSE;
2318:   return(0);
2319: }

2321: static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals)
2322: {
2323:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;

2326:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2327:   *vals = a->v+col*a->lda;
2328:   return(0);
2329: }

2331: static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals)
2332: {
2334:   *vals = 0; /* user cannot accidently use the array later */
2335:   return(0);
2336: }

2338: /* -------------------------------------------------------------------*/
2339: static struct _MatOps MatOps_Values = { MatSetValues_SeqDense,
2340:                                         MatGetRow_SeqDense,
2341:                                         MatRestoreRow_SeqDense,
2342:                                         MatMult_SeqDense,
2343:                                 /*  4*/ MatMultAdd_SeqDense,
2344:                                         MatMultTranspose_SeqDense,
2345:                                         MatMultTransposeAdd_SeqDense,
2346:                                         0,
2347:                                         0,
2348:                                         0,
2349:                                 /* 10*/ 0,
2350:                                         MatLUFactor_SeqDense,
2351:                                         MatCholeskyFactor_SeqDense,
2352:                                         MatSOR_SeqDense,
2353:                                         MatTranspose_SeqDense,
2354:                                 /* 15*/ MatGetInfo_SeqDense,
2355:                                         MatEqual_SeqDense,
2356:                                         MatGetDiagonal_SeqDense,
2357:                                         MatDiagonalScale_SeqDense,
2358:                                         MatNorm_SeqDense,
2359:                                 /* 20*/ MatAssemblyBegin_SeqDense,
2360:                                         MatAssemblyEnd_SeqDense,
2361:                                         MatSetOption_SeqDense,
2362:                                         MatZeroEntries_SeqDense,
2363:                                 /* 24*/ MatZeroRows_SeqDense,
2364:                                         0,
2365:                                         0,
2366:                                         0,
2367:                                         0,
2368:                                 /* 29*/ MatSetUp_SeqDense,
2369:                                         0,
2370:                                         0,
2371:                                         0,
2372:                                         0,
2373:                                 /* 34*/ MatDuplicate_SeqDense,
2374:                                         0,
2375:                                         0,
2376:                                         0,
2377:                                         0,
2378:                                 /* 39*/ MatAXPY_SeqDense,
2379:                                         MatCreateSubMatrices_SeqDense,
2380:                                         0,
2381:                                         MatGetValues_SeqDense,
2382:                                         MatCopy_SeqDense,
2383:                                 /* 44*/ MatGetRowMax_SeqDense,
2384:                                         MatScale_SeqDense,
2385:                                         MatShift_Basic,
2386:                                         0,
2387:                                         MatZeroRowsColumns_SeqDense,
2388:                                 /* 49*/ MatSetRandom_SeqDense,
2389:                                         0,
2390:                                         0,
2391:                                         0,
2392:                                         0,
2393:                                 /* 54*/ 0,
2394:                                         0,
2395:                                         0,
2396:                                         0,
2397:                                         0,
2398:                                 /* 59*/ 0,
2399:                                         MatDestroy_SeqDense,
2400:                                         MatView_SeqDense,
2401:                                         0,
2402:                                         0,
2403:                                 /* 64*/ 0,
2404:                                         0,
2405:                                         0,
2406:                                         0,
2407:                                         0,
2408:                                 /* 69*/ MatGetRowMaxAbs_SeqDense,
2409:                                         0,
2410:                                         0,
2411:                                         0,
2412:                                         0,
2413:                                 /* 74*/ 0,
2414:                                         0,
2415:                                         0,
2416:                                         0,
2417:                                         0,
2418:                                 /* 79*/ 0,
2419:                                         0,
2420:                                         0,
2421:                                         0,
2422:                                 /* 83*/ MatLoad_SeqDense,
2423:                                         0,
2424:                                         MatIsHermitian_SeqDense,
2425:                                         0,
2426:                                         0,
2427:                                         0,
2428:                                 /* 89*/ MatMatMult_SeqDense_SeqDense,
2429:                                         MatMatMultSymbolic_SeqDense_SeqDense,
2430:                                         MatMatMultNumeric_SeqDense_SeqDense,
2431:                                         MatPtAP_SeqDense_SeqDense,
2432:                                         MatPtAPSymbolic_SeqDense_SeqDense,
2433:                                 /* 94*/ MatPtAPNumeric_SeqDense_SeqDense,
2434:                                         MatMatTransposeMult_SeqDense_SeqDense,
2435:                                         MatMatTransposeMultSymbolic_SeqDense_SeqDense,
2436:                                         MatMatTransposeMultNumeric_SeqDense_SeqDense,
2437:                                         0,
2438:                                 /* 99*/ 0,
2439:                                         0,
2440:                                         0,
2441:                                         MatConjugate_SeqDense,
2442:                                         0,
2443:                                 /*104*/ 0,
2444:                                         MatRealPart_SeqDense,
2445:                                         MatImaginaryPart_SeqDense,
2446:                                         0,
2447:                                         0,
2448:                                 /*109*/ 0,
2449:                                         0,
2450:                                         MatGetRowMin_SeqDense,
2451:                                         MatGetColumnVector_SeqDense,
2452:                                         MatMissingDiagonal_SeqDense,
2453:                                 /*114*/ 0,
2454:                                         0,
2455:                                         0,
2456:                                         0,
2457:                                         0,
2458:                                 /*119*/ 0,
2459:                                         0,
2460:                                         0,
2461:                                         0,
2462:                                         0,
2463:                                 /*124*/ 0,
2464:                                         MatGetColumnNorms_SeqDense,
2465:                                         0,
2466:                                         0,
2467:                                         0,
2468:                                 /*129*/ 0,
2469:                                         MatTransposeMatMult_SeqDense_SeqDense,
2470:                                         MatTransposeMatMultSymbolic_SeqDense_SeqDense,
2471:                                         MatTransposeMatMultNumeric_SeqDense_SeqDense,
2472:                                         0,
2473:                                 /*134*/ 0,
2474:                                         0,
2475:                                         0,
2476:                                         0,
2477:                                         0,
2478:                                 /*139*/ 0,
2479:                                         0,
2480:                                         0,
2481:                                         0,
2482:                                         0,
2483:                                 /*144*/ MatCreateMPIMatConcatenateSeqMat_SeqDense
2484: };

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

2491:    Collective on MPI_Comm

2493:    Input Parameters:
2494: +  comm - MPI communicator, set to PETSC_COMM_SELF
2495: .  m - number of rows
2496: .  n - number of columns
2497: -  data - optional location of matrix data in column major order.  Set data=NULL for PETSc
2498:    to control all matrix memory allocation.

2500:    Output Parameter:
2501: .  A - the matrix

2503:    Notes:
2504:    The data input variable is intended primarily for Fortran programmers
2505:    who wish to allocate their own matrix memory space.  Most users should
2506:    set data=NULL.

2508:    Level: intermediate

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

2512: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
2513: @*/
2514: PetscErrorCode  MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
2515: {

2519:   MatCreate(comm,A);
2520:   MatSetSizes(*A,m,n,m,n);
2521:   MatSetType(*A,MATSEQDENSE);
2522:   MatSeqDenseSetPreallocation(*A,data);
2523:   return(0);
2524: }

2526: /*@C
2527:    MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements

2529:    Collective on MPI_Comm

2531:    Input Parameters:
2532: +  B - the matrix
2533: -  data - the array (or NULL)

2535:    Notes:
2536:    The data input variable is intended primarily for Fortran programmers
2537:    who wish to allocate their own matrix memory space.  Most users should
2538:    need not call this routine.

2540:    Level: intermediate

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

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

2546: @*/
2547: PetscErrorCode  MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
2548: {

2552:   PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
2553:   return(0);
2554: }

2556: PetscErrorCode  MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
2557: {
2558:   Mat_SeqDense   *b;

2562:   B->preallocated = PETSC_TRUE;

2564:   PetscLayoutSetUp(B->rmap);
2565:   PetscLayoutSetUp(B->cmap);

2567:   b       = (Mat_SeqDense*)B->data;
2568:   b->Mmax = B->rmap->n;
2569:   b->Nmax = B->cmap->n;
2570:   if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n;

2572:   PetscIntMultError(b->lda,b->Nmax,NULL);
2573:   if (!data) { /* petsc-allocated storage */
2574:     if (!b->user_alloc) { PetscFree(b->v); }
2575:     PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);
2576:     PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));

2578:     b->user_alloc = PETSC_FALSE;
2579:   } else { /* user-allocated storage */
2580:     if (!b->user_alloc) { PetscFree(b->v); }
2581:     b->v          = data;
2582:     b->user_alloc = PETSC_TRUE;
2583:   }
2584:   B->assembled = PETSC_TRUE;
2585:   return(0);
2586: }

2588: #if defined(PETSC_HAVE_ELEMENTAL)
2589: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
2590: {
2591:   Mat            mat_elemental;
2593:   PetscScalar    *array,*v_colwise;
2594:   PetscInt       M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols;

2597:   PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);
2598:   MatDenseGetArray(A,&array);
2599:   /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
2600:   k = 0;
2601:   for (j=0; j<N; j++) {
2602:     cols[j] = j;
2603:     for (i=0; i<M; i++) {
2604:       v_colwise[j*M+i] = array[k++];
2605:     }
2606:   }
2607:   for (i=0; i<M; i++) {
2608:     rows[i] = i;
2609:   }
2610:   MatDenseRestoreArray(A,&array);

2612:   MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
2613:   MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);
2614:   MatSetType(mat_elemental,MATELEMENTAL);
2615:   MatSetUp(mat_elemental);

2617:   /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
2618:   MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);
2619:   MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
2620:   MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);
2621:   PetscFree3(v_colwise,rows,cols);

2623:   if (reuse == MAT_INPLACE_MATRIX) {
2624:     MatHeaderReplace(A,&mat_elemental);
2625:   } else {
2626:     *newmat = mat_elemental;
2627:   }
2628:   return(0);
2629: }
2630: #endif

2632: /*@C
2633:   MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array

2635:   Input parameter:
2636: + A - the matrix
2637: - lda - the leading dimension

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

2644:   Level: intermediate

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

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

2650: @*/
2651: PetscErrorCode  MatSeqDenseSetLDA(Mat B,PetscInt lda)
2652: {
2653:   Mat_SeqDense *b = (Mat_SeqDense*)B->data;

2656:   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);
2657:   b->lda       = lda;
2658:   b->changelda = PETSC_FALSE;
2659:   b->Mmax      = PetscMax(b->Mmax,lda);
2660:   return(0);
2661: }

2663: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat)
2664: {
2666:   PetscMPIInt    size;

2669:   MPI_Comm_size(comm,&size);
2670:   if (size == 1) {
2671:     if (scall == MAT_INITIAL_MATRIX) {
2672:       MatDuplicate(inmat,MAT_COPY_VALUES,outmat);
2673:     } else {
2674:       MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);
2675:     }
2676:   } else {
2677:     MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);
2678:   }
2679:   return(0);
2680: }

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

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

2688:   Level: beginner

2690: .seealso: MatCreateSeqDense()

2692: M*/

2694: PETSC_EXTERN PetscErrorCode MatCreate_SeqDense(Mat B)
2695: {
2696:   Mat_SeqDense   *b;
2698:   PetscMPIInt    size;

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

2704:   PetscNewLog(B,&b);
2705:   PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
2706:   B->data = (void*)b;

2708:   b->roworiented = PETSC_TRUE;

2710:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);
2711:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);
2712:   PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);
2713:   PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);
2714:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);
2715:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);
2716:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);
2717: #if defined(PETSC_HAVE_ELEMENTAL)
2718:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);
2719: #endif
2720:   PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);
2721:   PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2722:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2723:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2724:   PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaij_seqdense_C",MatPtAP_SeqDense_SeqDense);
2725:   PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijperm_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2726:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijperm_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2727:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijperm_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2728:   PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijperm_seqdense_C",MatPtAP_SeqDense_SeqDense);
2729:   PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijmkl_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2730:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijmkl_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2731:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijmkl_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2732:   PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijmkl_seqdense_C",MatPtAP_SeqDense_SeqDense);

2734:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2735:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2736:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2737:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijperm_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2738:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijperm_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2739:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijperm_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);

2741:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijmkl_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2742:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijmkl_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2743:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijmkl_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2744:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);
2745:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);
2746:   PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
2747:   return(0);
2748: }

2750: /*@C
2751:    MatDenseGetColumn - gives access to a column of a dense matrix. This is only the local part of the column. You MUST call MatDenseRestoreColumn() to avoid memory bleeding.

2753:    Not Collective

2755:    Input Parameter:
2756: +  mat - a MATSEQDENSE or MATMPIDENSE matrix
2757: -  col - column index

2759:    Output Parameter:
2760: .  vals - pointer to the data

2762:    Level: intermediate

2764: .seealso: MatDenseRestoreColumn()
2765: @*/
2766: PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals)
2767: {

2771:   PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));
2772:   return(0);
2773: }

2775: /*@C
2776:    MatDenseRestoreColumn - returns access to a column of a dense matrix which is returned by MatDenseGetColumn().

2778:    Not Collective

2780:    Input Parameter:
2781: .  mat - a MATSEQDENSE or MATMPIDENSE matrix

2783:    Output Parameter:
2784: .  vals - pointer to the data

2786:    Level: intermediate

2788: .seealso: MatDenseGetColumn()
2789: @*/
2790: PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals)
2791: {

2795:   PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));
2796:   return(0);
2797: }