Actual source code: dense.c

petsc-master 2019-07-15
Report Typos and Errors

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

  6:  #include <../src/mat/impls/dense/seq/dense.h>
  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:     Vec xt;

118:     if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
119:     VecDuplicate(x,&xt);
120:     VecCopy(x,xt);
121:     VecScale(xt,-1.0);
122:     MatMultAdd(A,xt,b,b);
123:     VecDestroy(&xt);
124:     VecGetArrayRead(x,&xx);
125:     VecGetArray(b,&bb);
126:     for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
127:     VecRestoreArrayRead(x,&xx);
128:     VecRestoreArray(b,&bb);
129:   }

131:   for (i=0; i<N; i++) {
132:     slot = l->v + rows[i]*m;
133:     PetscArrayzero(slot,r);
134:   }
135:   for (i=0; i<N; i++) {
136:     slot = l->v + rows[i];
137:     for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
138:   }
139:   if (diag != 0.0) {
140:     if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
141:     for (i=0; i<N; i++) {
142:       slot  = l->v + (m+1)*rows[i];
143:       *slot = diag;
144:     }
145:   }
146:   return(0);
147: }

149: PetscErrorCode MatPtAPNumeric_SeqDense_SeqDense(Mat A,Mat P,Mat C)
150: {
151:   Mat_SeqDense   *c = (Mat_SeqDense*)(C->data);

155:   MatMatMultNumeric_SeqDense_SeqDense(A,P,c->ptapwork);
156:   MatTransposeMatMultNumeric_SeqDense_SeqDense(P,c->ptapwork,C);
157:   return(0);
158: }

160: PetscErrorCode MatPtAPSymbolic_SeqDense_SeqDense(Mat A,Mat P,PetscReal fill,Mat *C)
161: {
162:   Mat_SeqDense   *c;

166:   MatCreateSeqDense(PetscObjectComm((PetscObject)A),P->cmap->N,P->cmap->N,NULL,C);
167:   c = (Mat_SeqDense*)((*C)->data);
168:   MatCreateSeqDense(PetscObjectComm((PetscObject)A),A->rmap->N,P->cmap->N,NULL,&c->ptapwork);
169:   return(0);
170: }

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

177:   if (reuse == MAT_INITIAL_MATRIX) {
178:     MatPtAPSymbolic_SeqDense_SeqDense(A,P,fill,C);
179:   }
180:   PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);
181:   (*(*C)->ops->ptapnumeric)(A,P,*C);
182:   PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);
183:   return(0);
184: }

186: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
187: {
188:   Mat            B = NULL;
189:   Mat_SeqAIJ     *a = (Mat_SeqAIJ*)A->data;
190:   Mat_SeqDense   *b;
192:   PetscInt       *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->N,i;
193:   MatScalar      *av=a->a;
194:   PetscBool      isseqdense;

197:   if (reuse == MAT_REUSE_MATRIX) {
198:     PetscObjectTypeCompare((PetscObject)*newmat,MATSEQDENSE,&isseqdense);
199:     if (!isseqdense) SETERRQ1(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_USER,"Cannot reuse matrix of type %s",((PetscObject)(*newmat))->type_name);
200:   }
201:   if (reuse != MAT_REUSE_MATRIX) {
202:     MatCreate(PetscObjectComm((PetscObject)A),&B);
203:     MatSetSizes(B,m,n,m,n);
204:     MatSetType(B,MATSEQDENSE);
205:     MatSeqDenseSetPreallocation(B,NULL);
206:     b    = (Mat_SeqDense*)(B->data);
207:   } else {
208:     b    = (Mat_SeqDense*)((*newmat)->data);
209:     PetscArrayzero(b->v,m*n);
210:   }
211:   for (i=0; i<m; i++) {
212:     PetscInt j;
213:     for (j=0;j<ai[1]-ai[0];j++) {
214:       b->v[*aj*m+i] = *av;
215:       aj++;
216:       av++;
217:     }
218:     ai++;
219:   }

221:   if (reuse == MAT_INPLACE_MATRIX) {
222:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
223:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
224:     MatHeaderReplace(A,&B);
225:   } else {
226:     if (B) *newmat = B;
227:     MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
228:     MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
229:   }
230:   return(0);
231: }

233: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
234: {
235:   Mat            B;
236:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
238:   PetscInt       i, j;
239:   PetscInt       *rows, *nnz;
240:   MatScalar      *aa = a->v, *vals;

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

263:   if (reuse == MAT_INPLACE_MATRIX) {
264:     MatHeaderReplace(A,&B);
265:   } else {
266:     *newmat = B;
267:   }
268:   return(0);
269: }

271: static PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str)
272: {
273:   Mat_SeqDense   *x     = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data;
274:   PetscScalar    oalpha = alpha;
275:   PetscInt       j;
276:   PetscBLASInt   N,m,ldax,lday,one = 1;

280:   PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);
281:   PetscBLASIntCast(X->rmap->n,&m);
282:   PetscBLASIntCast(x->lda,&ldax);
283:   PetscBLASIntCast(y->lda,&lday);
284:   if (ldax>m || lday>m) {
285:     for (j=0; j<X->cmap->n; j++) {
286:       PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&oalpha,x->v+j*ldax,&one,y->v+j*lday,&one));
287:     }
288:   } else {
289:     PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&oalpha,x->v,&one,y->v,&one));
290:   }
291:   PetscObjectStateIncrease((PetscObject)Y);
292:   PetscLogFlops(PetscMax(2*N-1,0));
293:   return(0);
294: }

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

301:   info->block_size        = 1.0;
302:   info->nz_allocated      = (double)N;
303:   info->nz_used           = (double)N;
304:   info->nz_unneeded       = (double)0;
305:   info->assemblies        = (double)A->num_ass;
306:   info->mallocs           = 0;
307:   info->memory            = ((PetscObject)A)->mem;
308:   info->fill_ratio_given  = 0;
309:   info->fill_ratio_needed = 0;
310:   info->factor_mallocs    = 0;
311:   return(0);
312: }

314: static PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha)
315: {
316:   Mat_SeqDense   *a     = (Mat_SeqDense*)A->data;
317:   PetscScalar    oalpha = alpha;
319:   PetscBLASInt   one = 1,j,nz,lda;

322:   PetscBLASIntCast(a->lda,&lda);
323:   if (lda>A->rmap->n) {
324:     PetscBLASIntCast(A->rmap->n,&nz);
325:     for (j=0; j<A->cmap->n; j++) {
326:       PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v+j*lda,&one));
327:     }
328:   } else {
329:     PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);
330:     PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v,&one));
331:   }
332:   PetscLogFlops(nz);
333:   return(0);
334: }

336: static PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool  *fl)
337: {
338:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
339:   PetscInt     i,j,m = A->rmap->n,N;
340:   PetscScalar  *v = a->v;

343:   *fl = PETSC_FALSE;
344:   if (A->rmap->n != A->cmap->n) return(0);
345:   N = a->lda;

347:   for (i=0; i<m; i++) {
348:     for (j=i+1; j<m; j++) {
349:       if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) return(0);
350:     }
351:   }
352:   *fl = PETSC_TRUE;
353:   return(0);
354: }

356: static PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues)
357: {
358:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data,*l;
360:   PetscInt       lda = (PetscInt)mat->lda,j,m;

363:   PetscLayoutReference(A->rmap,&newi->rmap);
364:   PetscLayoutReference(A->cmap,&newi->cmap);
365:   MatSeqDenseSetPreallocation(newi,NULL);
366:   if (cpvalues == MAT_COPY_VALUES) {
367:     l = (Mat_SeqDense*)newi->data;
368:     if (lda>A->rmap->n) {
369:       m = A->rmap->n;
370:       for (j=0; j<A->cmap->n; j++) {
371:         PetscArraycpy(l->v+j*m,mat->v+j*lda,m);
372:       }
373:     } else {
374:       PetscArraycpy(l->v,mat->v,A->rmap->n*A->cmap->n);
375:     }
376:   }
377:   newi->assembled = PETSC_TRUE;
378:   return(0);
379: }

381: static PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat)
382: {

386:   MatCreate(PetscObjectComm((PetscObject)A),newmat);
387:   MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
388:   MatSetType(*newmat,((PetscObject)A)->type_name);
389:   MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);
390:   return(0);
391: }


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

396: static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
397: {
398:   MatFactorInfo  info;

402:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
403:   MatLUFactor_SeqDense(fact,0,0,&info);
404:   return(0);
405: }

407: static PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy)
408: {
409:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
410:   PetscErrorCode    ierr;
411:   const PetscScalar *x;
412:   PetscScalar       *y;
413:   PetscBLASInt      one = 1,info,m;

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

459: static PetscErrorCode MatMatSolve_SeqDense(Mat A,Mat B,Mat X)
460: {
461:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
462:   PetscErrorCode    ierr;
463:   const PetscScalar *b;
464:   PetscScalar       *x;
465:   PetscInt          n;
466:   PetscBLASInt      nrhs,info,m;
467:   PetscBool         flg;

470:   PetscBLASIntCast(A->rmap->n,&m);
471:   PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
472:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
473:   PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
474:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");

476:   MatGetSize(B,NULL,&n);
477:   PetscBLASIntCast(n,&nrhs);
478:   MatDenseGetArrayRead(B,&b);
479:   MatDenseGetArray(X,&x);

481:   PetscArraycpy(x,b,m*nrhs);

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

517:   MatDenseRestoreArrayRead(B,&b);
518:   MatDenseRestoreArray(X,&x);
519:   PetscLogFlops(nrhs*(2.0*m*m - m));
520:   return(0);
521: }

523: static PetscErrorCode MatConjugate_SeqDense(Mat);

525: static PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy)
526: {
527:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
528:   PetscErrorCode    ierr;
529:   const PetscScalar *x;
530:   PetscScalar       *y;
531:   PetscBLASInt      one = 1,info,m;

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

584: static PetscErrorCode MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
585: {
586:   PetscErrorCode    ierr;
587:   const PetscScalar *x;
588:   PetscScalar       *y,sone = 1.0;
589:   Vec               tmp = 0;

592:   VecGetArrayRead(xx,&x);
593:   VecGetArray(yy,&y);
594:   if (!A->rmap->n || !A->cmap->n) return(0);
595:   if (yy == zz) {
596:     VecDuplicate(yy,&tmp);
597:     PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);
598:     VecCopy(yy,tmp);
599:   }
600:   MatSolve_SeqDense(A,xx,yy);
601:   if (tmp) {
602:     VecAXPY(yy,sone,tmp);
603:     VecDestroy(&tmp);
604:   } else {
605:     VecAXPY(yy,sone,zz);
606:   }
607:   VecRestoreArrayRead(xx,&x);
608:   VecRestoreArray(yy,&y);
609:   PetscLogFlops(A->cmap->n);
610:   return(0);
611: }

613: static PetscErrorCode MatSolveTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
614: {
615:   PetscErrorCode    ierr;
616:   const PetscScalar *x;
617:   PetscScalar       *y,sone = 1.0;
618:   Vec               tmp = NULL;

621:   if (!A->rmap->n || !A->cmap->n) return(0);
622:   VecGetArrayRead(xx,&x);
623:   VecGetArray(yy,&y);
624:   if (yy == zz) {
625:     VecDuplicate(yy,&tmp);
626:     PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);
627:     VecCopy(yy,tmp);
628:   }
629:   MatSolveTranspose_SeqDense(A,xx,yy);
630:   if (tmp) {
631:     VecAXPY(yy,sone,tmp);
632:     VecDestroy(&tmp);
633:   } else {
634:     VecAXPY(yy,sone,zz);
635:   }
636:   VecRestoreArrayRead(xx,&x);
637:   VecRestoreArray(yy,&y);
638:   PetscLogFlops(A->cmap->n);
639:   return(0);
640: }

642: /* ---------------------------------------------------------------*/
643: /* COMMENT: I have chosen to hide row permutation in the pivots,
644:    rather than put it in the Mat->row slot.*/
645: static PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo)
646: {
647: #if defined(PETSC_MISSING_LAPACK_GETRF)
649:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable.");
650: #else
651:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
653:   PetscBLASInt   n,m,info;

656:   PetscBLASIntCast(A->cmap->n,&n);
657:   PetscBLASIntCast(A->rmap->n,&m);
658:   if (!mat->pivots) {
659:     PetscMalloc1(A->rmap->n,&mat->pivots);
660:     PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
661:   }
662:   if (!A->rmap->n || !A->cmap->n) return(0);
663:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
664:   PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info));
665:   PetscFPTrapPop();

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

670:   A->ops->solve             = MatSolve_SeqDense;
671:   A->ops->matsolve          = MatMatSolve_SeqDense;
672:   A->ops->solvetranspose    = MatSolveTranspose_SeqDense;
673:   A->ops->solveadd          = MatSolveAdd_SeqDense;
674:   A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
675:   A->factortype             = MAT_FACTOR_LU;

677:   PetscFree(A->solvertype);
678:   PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);

680:   PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);
681: #endif
682:   return(0);
683: }

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

697:   PetscBLASIntCast(A->cmap->n,&n);
698:   if (!A->rmap->n || !A->cmap->n) return(0);
699:   if (A->spd) {
700:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
701:     PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info));
702:     PetscFPTrapPop();
703: #if defined(PETSC_USE_COMPLEX)
704:   } else if (A->hermitian) {
705:     if (!mat->pivots) {
706:       PetscMalloc1(A->rmap->n,&mat->pivots);
707:       PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
708:     }
709:     if (!mat->fwork) {
710:       PetscScalar dummy;

712:       mat->lfwork = -1;
713:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
714:       PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
715:       PetscFPTrapPop();
716:       mat->lfwork = (PetscInt)PetscRealPart(dummy);
717:       PetscMalloc1(mat->lfwork,&mat->fwork);
718:       PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
719:     }
720:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
721:     PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
722:     PetscFPTrapPop();
723: #endif
724:   } else { /* symmetric case */
725:     if (!mat->pivots) {
726:       PetscMalloc1(A->rmap->n,&mat->pivots);
727:       PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
728:     }
729:     if (!mat->fwork) {
730:       PetscScalar dummy;

732:       mat->lfwork = -1;
733:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
734:       PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
735:       PetscFPTrapPop();
736:       mat->lfwork = (PetscInt)PetscRealPart(dummy);
737:       PetscMalloc1(mat->lfwork,&mat->fwork);
738:       PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
739:     }
740:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
741:     PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
742:     PetscFPTrapPop();
743:   }
744:   if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1);

746:   A->ops->solve             = MatSolve_SeqDense;
747:   A->ops->matsolve          = MatMatSolve_SeqDense;
748:   A->ops->solvetranspose    = MatSolveTranspose_SeqDense;
749:   A->ops->solveadd          = MatSolveAdd_SeqDense;
750:   A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
751:   A->factortype             = MAT_FACTOR_CHOLESKY;

753:   PetscFree(A->solvertype);
754:   PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);

756:   PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
757: #endif
758:   return(0);
759: }


762: PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
763: {
765:   MatFactorInfo  info;

768:   info.fill = 1.0;

770:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
771:   MatCholeskyFactor_SeqDense(fact,0,&info);
772:   return(0);
773: }

775: static PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
776: {
778:   fact->assembled                  = PETSC_TRUE;
779:   fact->preallocated               = PETSC_TRUE;
780:   fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
781:   fact->ops->solve                 = MatSolve_SeqDense;
782:   fact->ops->matsolve              = MatMatSolve_SeqDense;
783:   fact->ops->solvetranspose        = MatSolveTranspose_SeqDense;
784:   fact->ops->solveadd              = MatSolveAdd_SeqDense;
785:   fact->ops->solvetransposeadd     = MatSolveTransposeAdd_SeqDense;
786:   return(0);
787: }

789: static PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info)
790: {
792:   fact->preallocated           = PETSC_TRUE;
793:   fact->assembled              = PETSC_TRUE;
794:   fact->ops->lufactornumeric   = MatLUFactorNumeric_SeqDense;
795:   fact->ops->solve             = MatSolve_SeqDense;
796:   fact->ops->matsolve          = MatMatSolve_SeqDense;
797:   fact->ops->solvetranspose    = MatSolveTranspose_SeqDense;
798:   fact->ops->solveadd          = MatSolveAdd_SeqDense;
799:   fact->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
800:   return(0);
801: }

803: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact)
804: {

808:   MatCreate(PetscObjectComm((PetscObject)A),fact);
809:   MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
810:   MatSetType(*fact,((PetscObject)A)->type_name);
811:   if (ftype == MAT_FACTOR_LU) {
812:     (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
813:   } else {
814:     (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
815:   }
816:   (*fact)->factortype = ftype;

818:   PetscFree((*fact)->solvertype);
819:   PetscStrallocpy(MATSOLVERPETSC,&(*fact)->solvertype);
820:   return(0);
821: }

823: /* ------------------------------------------------------------------*/
824: static PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx)
825: {
826:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
827:   PetscScalar       *x,*v = mat->v,zero = 0.0,xt;
828:   const PetscScalar *b;
829:   PetscErrorCode    ierr;
830:   PetscInt          m = A->rmap->n,i;
831:   PetscBLASInt      o = 1,bm;

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

863: /* -----------------------------------------------------------------*/
864: PETSC_INTERN PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)
865: {
866:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
867:   const PetscScalar *v   = mat->v,*x;
868:   PetscScalar       *y;
869:   PetscErrorCode    ierr;
870:   PetscBLASInt      m, n,_One=1;
871:   PetscScalar       _DOne=1.0,_DZero=0.0;

874:   PetscBLASIntCast(A->rmap->n,&m);
875:   PetscBLASIntCast(A->cmap->n,&n);
876:   VecGetArrayRead(xx,&x);
877:   VecGetArray(yy,&y);
878:   if (!A->rmap->n || !A->cmap->n) {
879:     PetscBLASInt i;
880:     for (i=0; i<n; i++) y[i] = 0.0;
881:   } else {
882:     PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One));
883:     PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);
884:   }
885:   VecRestoreArrayRead(xx,&x);
886:   VecRestoreArray(yy,&y);
887:   return(0);
888: }

890: PETSC_INTERN PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy)
891: {
892:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
893:   PetscScalar       *y,_DOne=1.0,_DZero=0.0;
894:   PetscErrorCode    ierr;
895:   PetscBLASInt      m, n, _One=1;
896:   const PetscScalar *v = mat->v,*x;

899:   PetscBLASIntCast(A->rmap->n,&m);
900:   PetscBLASIntCast(A->cmap->n,&n);
901:   VecGetArrayRead(xx,&x);
902:   VecGetArray(yy,&y);
903:   if (!A->rmap->n || !A->cmap->n) {
904:     PetscBLASInt i;
905:     for (i=0; i<m; i++) y[i] = 0.0;
906:   } else {
907:     PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One));
908:     PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);
909:   }
910:   VecRestoreArrayRead(xx,&x);
911:   VecRestoreArray(yy,&y);
912:   return(0);
913: }

915: PETSC_INTERN PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
916: {
917:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
918:   const PetscScalar *v = mat->v,*x;
919:   PetscScalar       *y,_DOne=1.0;
920:   PetscErrorCode    ierr;
921:   PetscBLASInt      m, n, _One=1;

924:   PetscBLASIntCast(A->rmap->n,&m);
925:   PetscBLASIntCast(A->cmap->n,&n);
926:   if (!A->rmap->n || !A->cmap->n) return(0);
927:   if (zz != yy) {VecCopy(zz,yy);}
928:   VecGetArrayRead(xx,&x);
929:   VecGetArray(yy,&y);
930:   PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
931:   VecRestoreArrayRead(xx,&x);
932:   VecRestoreArray(yy,&y);
933:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
934:   return(0);
935: }

937: static PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
938: {
939:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
940:   const PetscScalar *v = mat->v,*x;
941:   PetscScalar       *y;
942:   PetscErrorCode    ierr;
943:   PetscBLASInt      m, n, _One=1;
944:   PetscScalar       _DOne=1.0;

947:   PetscBLASIntCast(A->rmap->n,&m);
948:   PetscBLASIntCast(A->cmap->n,&n);
949:   if (!A->rmap->n || !A->cmap->n) return(0);
950:   if (zz != yy) {VecCopy(zz,yy);}
951:   VecGetArrayRead(xx,&x);
952:   VecGetArray(yy,&y);
953:   PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
954:   VecRestoreArrayRead(xx,&x);
955:   VecRestoreArray(yy,&y);
956:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
957:   return(0);
958: }

960: /* -----------------------------------------------------------------*/
961: static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
962: {
963:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
964:   PetscScalar    *v;
966:   PetscInt       i;

969:   *ncols = A->cmap->n;
970:   if (cols) {
971:     PetscMalloc1(A->cmap->n+1,cols);
972:     for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
973:   }
974:   if (vals) {
975:     PetscMalloc1(A->cmap->n+1,vals);
976:     v    = mat->v + row;
977:     for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
978:   }
979:   return(0);
980: }

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

987:   if (cols) {PetscFree(*cols);}
988:   if (vals) {PetscFree(*vals); }
989:   return(0);
990: }
991: /* ----------------------------------------------------------------*/
992: static PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
993: {
994:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
995:   PetscInt     i,j,idx=0;

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

1062: static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
1063: {
1064:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1065:   PetscInt     i,j;

1068:   /* row-oriented output */
1069:   for (i=0; i<m; i++) {
1070:     if (indexm[i] < 0) {v += n;continue;}
1071:     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);
1072:     for (j=0; j<n; j++) {
1073:       if (indexn[j] < 0) {v++; continue;}
1074:       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);
1075:       *v++ = mat->v[indexn[j]*mat->lda + indexm[i]];
1076:     }
1077:   }
1078:   return(0);
1079: }

1081: /* -----------------------------------------------------------------*/

1083: static PetscErrorCode MatLoad_SeqDense_Binary(Mat newmat,PetscViewer viewer)
1084: {
1085:   Mat_SeqDense   *a;
1087:   PetscInt       *scols,i,j,nz,header[4];
1088:   int            fd;
1089:   PetscMPIInt    size;
1090:   PetscInt       *rowlengths = 0,M,N,*cols,grows,gcols;
1091:   PetscScalar    *vals,*svals,*v,*w;
1092:   MPI_Comm       comm;

1095:   PetscObjectGetComm((PetscObject)viewer,&comm);
1096:   MPI_Comm_size(comm,&size);
1097:   if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor");
1098:   PetscViewerBinaryGetDescriptor(viewer,&fd);
1099:   PetscBinaryRead(fd,header,4,NULL,PETSC_INT);
1100:   if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not matrix object");
1101:   M = header[1]; N = header[2]; nz = header[3];

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

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

1136:     a = (Mat_SeqDense*)newmat->data;
1137:     v = a->v;

1139:     /* read column indices and nonzeros */
1140:     PetscMalloc1(nz+1,&scols);
1141:     cols = scols;
1142:     PetscBinaryRead(fd,cols,nz,NULL,PETSC_INT);
1143:     PetscMalloc1(nz+1,&svals);
1144:     vals = svals;
1145:     PetscBinaryRead(fd,vals,nz,NULL,PETSC_SCALAR);

1147:     /* insert into matrix */
1148:     for (i=0; i<M; i++) {
1149:       for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j];
1150:       svals += rowlengths[i]; scols += rowlengths[i];
1151:     }
1152:     PetscFree(vals);
1153:     PetscFree(cols);
1154:     PetscFree(rowlengths);
1155:   }
1156:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1157:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
1158:   return(0);
1159: }

1161: PetscErrorCode MatLoad_SeqDense(Mat newMat, PetscViewer viewer)
1162: {
1163:   PetscBool      isbinary, ishdf5;

1169:   /* force binary viewer to load .info file if it has not yet done so */
1170:   PetscViewerSetUp(viewer);
1171:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1172:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,  &ishdf5);
1173:   if (isbinary) {
1174:     MatLoad_SeqDense_Binary(newMat,viewer);
1175:   } else if (ishdf5) {
1176: #if defined(PETSC_HAVE_HDF5)
1177:     MatLoad_Dense_HDF5(newMat,viewer);
1178: #else
1179:     SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1180: #endif
1181:   } else {
1182:     SETERRQ2(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"Viewer type %s not yet supported for reading %s matrices",((PetscObject)viewer)->type_name,((PetscObject)newMat)->type_name);
1183:   }
1184:   return(0);
1185: }

1187: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
1188: {
1189:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
1190:   PetscErrorCode    ierr;
1191:   PetscInt          i,j;
1192:   const char        *name;
1193:   PetscScalar       *v;
1194:   PetscViewerFormat format;
1195: #if defined(PETSC_USE_COMPLEX)
1196:   PetscBool         allreal = PETSC_TRUE;
1197: #endif

1200:   PetscViewerGetFormat(viewer,&format);
1201:   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1202:     return(0);  /* do nothing for now */
1203:   } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1204:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1205:     for (i=0; i<A->rmap->n; i++) {
1206:       v    = a->v + i;
1207:       PetscViewerASCIIPrintf(viewer,"row %D:",i);
1208:       for (j=0; j<A->cmap->n; j++) {
1209: #if defined(PETSC_USE_COMPLEX)
1210:         if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
1211:           PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1212:         } else if (PetscRealPart(*v)) {
1213:           PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));
1214:         }
1215: #else
1216:         if (*v) {
1217:           PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);
1218:         }
1219: #endif
1220:         v += a->lda;
1221:       }
1222:       PetscViewerASCIIPrintf(viewer,"\n");
1223:     }
1224:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1225:   } else {
1226:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1227: #if defined(PETSC_USE_COMPLEX)
1228:     /* determine if matrix has all real values */
1229:     v = a->v;
1230:     for (i=0; i<A->rmap->n*A->cmap->n; i++) {
1231:       if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;}
1232:     }
1233: #endif
1234:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
1235:       PetscObjectGetName((PetscObject)A,&name);
1236:       PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);
1237:       PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);
1238:       PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
1239:     }

1241:     for (i=0; i<A->rmap->n; i++) {
1242:       v = a->v + i;
1243:       for (j=0; j<A->cmap->n; j++) {
1244: #if defined(PETSC_USE_COMPLEX)
1245:         if (allreal) {
1246:           PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));
1247:         } else {
1248:           PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1249:         }
1250: #else
1251:         PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);
1252: #endif
1253:         v += a->lda;
1254:       }
1255:       PetscViewerASCIIPrintf(viewer,"\n");
1256:     }
1257:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
1258:       PetscViewerASCIIPrintf(viewer,"];\n");
1259:     }
1260:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1261:   }
1262:   PetscViewerFlush(viewer);
1263:   return(0);
1264: }

1266: static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer)
1267: {
1268:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
1269:   PetscErrorCode    ierr;
1270:   int               fd;
1271:   PetscInt          ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n;
1272:   PetscScalar       *v,*anonz,*vals;
1273:   PetscViewerFormat format;

1276:   PetscViewerBinaryGetDescriptor(viewer,&fd);

1278:   PetscViewerGetFormat(viewer,&format);
1279:   if (format == PETSC_VIEWER_NATIVE) {
1280:     /* store the matrix as a dense matrix */
1281:     PetscMalloc1(4,&col_lens);

1283:     col_lens[0] = MAT_FILE_CLASSID;
1284:     col_lens[1] = m;
1285:     col_lens[2] = n;
1286:     col_lens[3] = MATRIX_BINARY_FORMAT_DENSE;

1288:     PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);
1289:     PetscFree(col_lens);

1291:     /* write out matrix, by rows */
1292:     PetscMalloc1(m*n+1,&vals);
1293:     v    = a->v;
1294:     for (j=0; j<n; j++) {
1295:       for (i=0; i<m; i++) {
1296:         vals[j + i*n] = *v++;
1297:       }
1298:     }
1299:     PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);
1300:     PetscFree(vals);
1301:   } else {
1302:     PetscMalloc1(4+nz,&col_lens);

1304:     col_lens[0] = MAT_FILE_CLASSID;
1305:     col_lens[1] = m;
1306:     col_lens[2] = n;
1307:     col_lens[3] = nz;

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

1313:     /* Possibly should write in smaller increments, not whole matrix at once? */
1314:     /* store column indices (zero start index) */
1315:     ict = 0;
1316:     for (i=0; i<m; i++) {
1317:       for (j=0; j<n; j++) col_lens[ict++] = j;
1318:     }
1319:     PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);
1320:     PetscFree(col_lens);

1322:     /* store nonzero values */
1323:     PetscMalloc1(nz+1,&anonz);
1324:     ict  = 0;
1325:     for (i=0; i<m; i++) {
1326:       v = a->v + i;
1327:       for (j=0; j<n; j++) {
1328:         anonz[ict++] = *v; v += a->lda;
1329:       }
1330:     }
1331:     PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);
1332:     PetscFree(anonz);
1333:   }
1334:   return(0);
1335: }

1337:  #include <petscdraw.h>
1338: static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa)
1339: {
1340:   Mat               A  = (Mat) Aa;
1341:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
1342:   PetscErrorCode    ierr;
1343:   PetscInt          m  = A->rmap->n,n = A->cmap->n,i,j;
1344:   int               color = PETSC_DRAW_WHITE;
1345:   PetscScalar       *v = a->v;
1346:   PetscViewer       viewer;
1347:   PetscReal         xl,yl,xr,yr,x_l,x_r,y_l,y_r;
1348:   PetscViewerFormat format;

1351:   PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
1352:   PetscViewerGetFormat(viewer,&format);
1353:   PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);

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

1357:   if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1358:     PetscDrawCollectiveBegin(draw);
1359:     /* Blue for negative and Red for positive */
1360:     for (j = 0; j < n; j++) {
1361:       x_l = j; x_r = x_l + 1.0;
1362:       for (i = 0; i < m; i++) {
1363:         y_l = m - i - 1.0;
1364:         y_r = y_l + 1.0;
1365:         if (PetscRealPart(v[j*m+i]) >  0.) {
1366:           color = PETSC_DRAW_RED;
1367:         } else if (PetscRealPart(v[j*m+i]) <  0.) {
1368:           color = PETSC_DRAW_BLUE;
1369:         } else {
1370:           continue;
1371:         }
1372:         PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1373:       }
1374:     }
1375:     PetscDrawCollectiveEnd(draw);
1376:   } else {
1377:     /* use contour shading to indicate magnitude of values */
1378:     /* first determine max of all nonzero values */
1379:     PetscReal minv = 0.0, maxv = 0.0;
1380:     PetscDraw popup;

1382:     for (i=0; i < m*n; i++) {
1383:       if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1384:     }
1385:     if (minv >= maxv) maxv = minv + PETSC_SMALL;
1386:     PetscDrawGetPopup(draw,&popup);
1387:     PetscDrawScalePopup(popup,minv,maxv);

1389:     PetscDrawCollectiveBegin(draw);
1390:     for (j=0; j<n; j++) {
1391:       x_l = j;
1392:       x_r = x_l + 1.0;
1393:       for (i=0; i<m; i++) {
1394:         y_l = m - i - 1.0;
1395:         y_r = y_l + 1.0;
1396:         color = PetscDrawRealToColor(PetscAbsScalar(v[j*m+i]),minv,maxv);
1397:         PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1398:       }
1399:     }
1400:     PetscDrawCollectiveEnd(draw);
1401:   }
1402:   return(0);
1403: }

1405: static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
1406: {
1407:   PetscDraw      draw;
1408:   PetscBool      isnull;
1409:   PetscReal      xr,yr,xl,yl,h,w;

1413:   PetscViewerDrawGetDraw(viewer,0,&draw);
1414:   PetscDrawIsNull(draw,&isnull);
1415:   if (isnull) return(0);

1417:   xr   = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
1418:   xr  += w;          yr += h;        xl = -w;     yl = -h;
1419:   PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
1420:   PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
1421:   PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
1422:   PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
1423:   PetscDrawSave(draw);
1424:   return(0);
1425: }

1427: PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
1428: {
1430:   PetscBool      iascii,isbinary,isdraw;

1433:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1434:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1435:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);

1437:   if (iascii) {
1438:     MatView_SeqDense_ASCII(A,viewer);
1439:   } else if (isbinary) {
1440:     MatView_SeqDense_Binary(A,viewer);
1441:   } else if (isdraw) {
1442:     MatView_SeqDense_Draw(A,viewer);
1443:   }
1444:   return(0);
1445: }

1447: static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A,const PetscScalar array[])
1448: {
1449:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

1452:   a->unplacedarray       = a->v;
1453:   a->unplaced_user_alloc = a->user_alloc;
1454:   a->v                   = (PetscScalar*) array;
1455:   return(0);
1456: }

1458: static PetscErrorCode MatDenseResetArray_SeqDense(Mat A)
1459: {
1460:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

1463:   a->v             = a->unplacedarray;
1464:   a->user_alloc    = a->unplaced_user_alloc;
1465:   a->unplacedarray = NULL;
1466:   return(0);
1467: }

1469: static PetscErrorCode MatDestroy_SeqDense(Mat mat)
1470: {
1471:   Mat_SeqDense   *l = (Mat_SeqDense*)mat->data;

1475: #if defined(PETSC_USE_LOG)
1476:   PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n);
1477: #endif
1478:   PetscFree(l->pivots);
1479:   PetscFree(l->fwork);
1480:   MatDestroy(&l->ptapwork);
1481:   if (!l->user_alloc) {PetscFree(l->v);}
1482:   PetscFree(mat->data);

1484:   PetscObjectChangeTypeName((PetscObject)mat,0);
1485:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetLDA_C",NULL);
1486:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);
1487:   PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL);
1488:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL);
1489:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);
1490:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);
1491: #if defined(PETSC_HAVE_ELEMENTAL)
1492:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);
1493: #endif
1494:   PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);
1495:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);
1496:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);
1497:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);
1498:   PetscObjectComposeFunction((PetscObject)mat,"MatPtAP_seqaij_seqdense_C",NULL);
1499:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);
1500:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);
1501:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);
1502:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumn_C",NULL);
1503:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumn_C",NULL);
1504:   return(0);
1505: }

1507: static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
1508: {
1509:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1511:   PetscInt       k,j,m,n,M;
1512:   PetscScalar    *v,tmp;

1515:   v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n;
1516:   if (reuse == MAT_INPLACE_MATRIX && m == n) { /* in place transpose */
1517:     for (j=0; j<m; j++) {
1518:       for (k=0; k<j; k++) {
1519:         tmp        = v[j + k*M];
1520:         v[j + k*M] = v[k + j*M];
1521:         v[k + j*M] = tmp;
1522:       }
1523:     }
1524:   } else { /* out-of-place transpose */
1525:     Mat          tmat;
1526:     Mat_SeqDense *tmatd;
1527:     PetscScalar  *v2;
1528:     PetscInt     M2;

1530:     if (reuse != MAT_REUSE_MATRIX) {
1531:       MatCreate(PetscObjectComm((PetscObject)A),&tmat);
1532:       MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);
1533:       MatSetType(tmat,((PetscObject)A)->type_name);
1534:       MatSeqDenseSetPreallocation(tmat,NULL);
1535:     } else {
1536:       tmat = *matout;
1537:     }
1538:     tmatd = (Mat_SeqDense*)tmat->data;
1539:     v = mat->v; v2 = tmatd->v; M2 = tmatd->lda;
1540:     for (j=0; j<n; j++) {
1541:       for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M];
1542:     }
1543:     MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
1544:     MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
1545:     if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) *matout = tmat;
1546:     else {
1547:       MatHeaderMerge(A,&tmat);
1548:     }
1549:   }
1550:   return(0);
1551: }

1553: static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool  *flg)
1554: {
1555:   Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data;
1556:   Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data;
1557:   PetscInt     i,j;
1558:   PetscScalar  *v1,*v2;

1561:   if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; return(0);}
1562:   if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; return(0);}
1563:   for (i=0; i<A1->rmap->n; i++) {
1564:     v1 = mat1->v+i; v2 = mat2->v+i;
1565:     for (j=0; j<A1->cmap->n; j++) {
1566:       if (*v1 != *v2) {*flg = PETSC_FALSE; return(0);}
1567:       v1 += mat1->lda; v2 += mat2->lda;
1568:     }
1569:   }
1570:   *flg = PETSC_TRUE;
1571:   return(0);
1572: }

1574: static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
1575: {
1576:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1578:   PetscInt       i,n,len;
1579:   PetscScalar    *x,zero = 0.0;

1582:   VecSet(v,zero);
1583:   VecGetSize(v,&n);
1584:   VecGetArray(v,&x);
1585:   len  = PetscMin(A->rmap->n,A->cmap->n);
1586:   if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
1587:   for (i=0; i<len; i++) {
1588:     x[i] = mat->v[i*mat->lda + i];
1589:   }
1590:   VecRestoreArray(v,&x);
1591:   return(0);
1592: }

1594: static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
1595: {
1596:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1597:   const PetscScalar *l,*r;
1598:   PetscScalar       x,*v;
1599:   PetscErrorCode    ierr;
1600:   PetscInt          i,j,m = A->rmap->n,n = A->cmap->n;

1603:   if (ll) {
1604:     VecGetSize(ll,&m);
1605:     VecGetArrayRead(ll,&l);
1606:     if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size");
1607:     for (i=0; i<m; i++) {
1608:       x = l[i];
1609:       v = mat->v + i;
1610:       for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;}
1611:     }
1612:     VecRestoreArrayRead(ll,&l);
1613:     PetscLogFlops(1.0*n*m);
1614:   }
1615:   if (rr) {
1616:     VecGetSize(rr,&n);
1617:     VecGetArrayRead(rr,&r);
1618:     if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size");
1619:     for (i=0; i<n; i++) {
1620:       x = r[i];
1621:       v = mat->v + i*mat->lda;
1622:       for (j=0; j<m; j++) (*v++) *= x;
1623:     }
1624:     VecRestoreArrayRead(rr,&r);
1625:     PetscLogFlops(1.0*n*m);
1626:   }
1627:   return(0);
1628: }

1630: static PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1631: {
1632:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1633:   PetscScalar    *v   = mat->v;
1634:   PetscReal      sum  = 0.0;
1635:   PetscInt       lda  =mat->lda,m=A->rmap->n,i,j;

1639:   if (type == NORM_FROBENIUS) {
1640:     if (lda>m) {
1641:       for (j=0; j<A->cmap->n; j++) {
1642:         v = mat->v+j*lda;
1643:         for (i=0; i<m; i++) {
1644:           sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1645:         }
1646:       }
1647:     } else {
1648: #if defined(PETSC_USE_REAL___FP16)
1649:       PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n;
1650:       *nrm = BLASnrm2_(&cnt,v,&one);
1651:     }
1652: #else
1653:       for (i=0; i<A->cmap->n*A->rmap->n; i++) {
1654:         sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1655:       }
1656:     }
1657:     *nrm = PetscSqrtReal(sum);
1658: #endif
1659:     PetscLogFlops(2.0*A->cmap->n*A->rmap->n);
1660:   } else if (type == NORM_1) {
1661:     *nrm = 0.0;
1662:     for (j=0; j<A->cmap->n; j++) {
1663:       v   = mat->v + j*mat->lda;
1664:       sum = 0.0;
1665:       for (i=0; i<A->rmap->n; i++) {
1666:         sum += PetscAbsScalar(*v);  v++;
1667:       }
1668:       if (sum > *nrm) *nrm = sum;
1669:     }
1670:     PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1671:   } else if (type == NORM_INFINITY) {
1672:     *nrm = 0.0;
1673:     for (j=0; j<A->rmap->n; j++) {
1674:       v   = mat->v + j;
1675:       sum = 0.0;
1676:       for (i=0; i<A->cmap->n; i++) {
1677:         sum += PetscAbsScalar(*v); v += mat->lda;
1678:       }
1679:       if (sum > *nrm) *nrm = sum;
1680:     }
1681:     PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1682:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
1683:   return(0);
1684: }

1686: static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
1687: {
1688:   Mat_SeqDense   *aij = (Mat_SeqDense*)A->data;

1692:   switch (op) {
1693:   case MAT_ROW_ORIENTED:
1694:     aij->roworiented = flg;
1695:     break;
1696:   case MAT_NEW_NONZERO_LOCATIONS:
1697:   case MAT_NEW_NONZERO_LOCATION_ERR:
1698:   case MAT_NEW_NONZERO_ALLOCATION_ERR:
1699:   case MAT_NEW_DIAGONALS:
1700:   case MAT_KEEP_NONZERO_PATTERN:
1701:   case MAT_IGNORE_OFF_PROC_ENTRIES:
1702:   case MAT_USE_HASH_TABLE:
1703:   case MAT_IGNORE_ZERO_ENTRIES:
1704:   case MAT_IGNORE_LOWER_TRIANGULAR:
1705:   case MAT_SORTED_FULL:
1706:     PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
1707:     break;
1708:   case MAT_SPD:
1709:   case MAT_SYMMETRIC:
1710:   case MAT_STRUCTURALLY_SYMMETRIC:
1711:   case MAT_HERMITIAN:
1712:   case MAT_SYMMETRY_ETERNAL:
1713:     /* These options are handled directly by MatSetOption() */
1714:     break;
1715:   default:
1716:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]);
1717:   }
1718:   return(0);
1719: }

1721: static PetscErrorCode MatZeroEntries_SeqDense(Mat A)
1722: {
1723:   Mat_SeqDense   *l = (Mat_SeqDense*)A->data;
1725:   PetscInt       lda=l->lda,m=A->rmap->n,j;

1728:   if (lda>m) {
1729:     for (j=0; j<A->cmap->n; j++) {
1730:       PetscArrayzero(l->v+j*lda,m);
1731:     }
1732:   } else {
1733:     PetscArrayzero(l->v,A->rmap->n*A->cmap->n);
1734:   }
1735:   return(0);
1736: }

1738: static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1739: {
1740:   PetscErrorCode    ierr;
1741:   Mat_SeqDense      *l = (Mat_SeqDense*)A->data;
1742:   PetscInt          m  = l->lda, n = A->cmap->n, i,j;
1743:   PetscScalar       *slot,*bb;
1744:   const PetscScalar *xx;

1747: #if defined(PETSC_USE_DEBUG)
1748:   for (i=0; i<N; i++) {
1749:     if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
1750:     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);
1751:   }
1752: #endif

1754:   /* fix right hand side if needed */
1755:   if (x && b) {
1756:     VecGetArrayRead(x,&xx);
1757:     VecGetArray(b,&bb);
1758:     for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
1759:     VecRestoreArrayRead(x,&xx);
1760:     VecRestoreArray(b,&bb);
1761:   }

1763:   for (i=0; i<N; i++) {
1764:     slot = l->v + rows[i];
1765:     for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
1766:   }
1767:   if (diag != 0.0) {
1768:     if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
1769:     for (i=0; i<N; i++) {
1770:       slot  = l->v + (m+1)*rows[i];
1771:       *slot = diag;
1772:     }
1773:   }
1774:   return(0);
1775: }

1777: static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A,PetscInt *lda)
1778: {
1779:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;

1782:   *lda = mat->lda;
1783:   return(0);
1784: }

1786: static PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[])
1787: {
1788:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;

1791:   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");
1792:   *array = mat->v;
1793:   return(0);
1794: }

1796: static PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[])
1797: {
1799:   return(0);
1800: }

1802: /*@C
1803:    MatDenseGetLDA - gets the leading dimension of the array returned from MatDenseGetArray()

1805:    Logically Collective on Mat

1807:    Input Parameter:
1808: .  mat - a MATSEQDENSE or MATMPIDENSE matrix

1810:    Output Parameter:
1811: .   lda - the leading dimension

1813:    Level: intermediate

1815: .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatSeqDenseSetLDA()
1816: @*/
1817: PetscErrorCode  MatDenseGetLDA(Mat A,PetscInt *lda)
1818: {

1822:   PetscUseMethod(A,"MatDenseGetLDA_C",(Mat,PetscInt*),(A,lda));
1823:   return(0);
1824: }

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

1829:    Logically Collective on Mat

1831:    Input Parameter:
1832: .  mat - a MATSEQDENSE or MATMPIDENSE matrix

1834:    Output Parameter:
1835: .   array - pointer to the data

1837:    Level: intermediate

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

1846:   PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));
1847:   return(0);
1848: }

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

1853:    Logically Collective on Mat

1855:    Input Parameters:
1856: +  mat - a MATSEQDENSE or MATMPIDENSE matrix
1857: -  array - pointer to the data

1859:    Level: intermediate

1861: .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead()
1862: @*/
1863: PetscErrorCode  MatDenseRestoreArray(Mat A,PetscScalar **array)
1864: {

1868:   PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));
1869:   if (array) *array = NULL;
1870:   PetscObjectStateIncrease((PetscObject)A);
1871:   return(0);
1872: }

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

1877:    Not Collective

1879:    Input Parameter:
1880: .  mat - a MATSEQDENSE or MATMPIDENSE matrix

1882:    Output Parameter:
1883: .   array - pointer to the data

1885:    Level: intermediate

1887: .seealso: MatDenseRestoreArray(), MatDenseGetArray(), MatDenseRestoreArrayRead()
1888: @*/
1889: PetscErrorCode  MatDenseGetArrayRead(Mat A,const PetscScalar **array)
1890: {

1894:   PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array));
1895:   return(0);
1896: }

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

1901:    Not Collective

1903:    Input Parameters:
1904: +  mat - a MATSEQDENSE or MATMPIDENSE matrix
1905: -  array - pointer to the data

1907:    Level: intermediate

1909: .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArray()
1910: @*/
1911: PetscErrorCode  MatDenseRestoreArrayRead(Mat A,const PetscScalar **array)
1912: {

1916:   PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array));
1917:   if (array) *array = NULL;
1918:   return(0);
1919: }

1921: static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B)
1922: {
1923:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1925:   PetscInt       i,j,nrows,ncols;
1926:   const PetscInt *irow,*icol;
1927:   PetscScalar    *av,*bv,*v = mat->v;
1928:   Mat            newmat;

1931:   ISGetIndices(isrow,&irow);
1932:   ISGetIndices(iscol,&icol);
1933:   ISGetLocalSize(isrow,&nrows);
1934:   ISGetLocalSize(iscol,&ncols);

1936:   /* Check submatrixcall */
1937:   if (scall == MAT_REUSE_MATRIX) {
1938:     PetscInt n_cols,n_rows;
1939:     MatGetSize(*B,&n_rows,&n_cols);
1940:     if (n_rows != nrows || n_cols != ncols) {
1941:       /* resize the result matrix to match number of requested rows/columns */
1942:       MatSetSizes(*B,nrows,ncols,nrows,ncols);
1943:     }
1944:     newmat = *B;
1945:   } else {
1946:     /* Create and fill new matrix */
1947:     MatCreate(PetscObjectComm((PetscObject)A),&newmat);
1948:     MatSetSizes(newmat,nrows,ncols,nrows,ncols);
1949:     MatSetType(newmat,((PetscObject)A)->type_name);
1950:     MatSeqDenseSetPreallocation(newmat,NULL);
1951:   }

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

1956:   for (i=0; i<ncols; i++) {
1957:     av = v + mat->lda*icol[i];
1958:     for (j=0; j<nrows; j++) *bv++ = av[irow[j]];
1959:   }

1961:   /* Assemble the matrices so that the correct flags are set */
1962:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1963:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);

1965:   /* Free work space */
1966:   ISRestoreIndices(isrow,&irow);
1967:   ISRestoreIndices(iscol,&icol);
1968:   *B   = newmat;
1969:   return(0);
1970: }

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

1978:   if (scall == MAT_INITIAL_MATRIX) {
1979:     PetscCalloc1(n+1,B);
1980:   }

1982:   for (i=0; i<n; i++) {
1983:     MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
1984:   }
1985:   return(0);
1986: }

1988: static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
1989: {
1991:   return(0);
1992: }

1994: static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
1995: {
1997:   return(0);
1998: }

2000: static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
2001: {
2002:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data;
2004:   PetscInt       lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j;

2007:   /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
2008:   if (A->ops->copy != B->ops->copy) {
2009:     MatCopy_Basic(A,B,str);
2010:     return(0);
2011:   }
2012:   if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)");
2013:   if (lda1>m || lda2>m) {
2014:     for (j=0; j<n; j++) {
2015:       PetscArraycpy(b->v+j*lda2,a->v+j*lda1,m);
2016:     }
2017:   } else {
2018:     PetscArraycpy(b->v,a->v,A->rmap->n*A->cmap->n);
2019:   }
2020:   PetscObjectStateIncrease((PetscObject)B);
2021:   return(0);
2022: }

2024: static PetscErrorCode MatSetUp_SeqDense(Mat A)
2025: {

2029:    MatSeqDenseSetPreallocation(A,0);
2030:   return(0);
2031: }

2033: static PetscErrorCode MatConjugate_SeqDense(Mat A)
2034: {
2035:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2036:   PetscInt     i,nz = A->rmap->n*A->cmap->n;
2037:   PetscScalar  *aa = a->v;

2040:   for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
2041:   return(0);
2042: }

2044: static PetscErrorCode MatRealPart_SeqDense(Mat A)
2045: {
2046:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2047:   PetscInt     i,nz = A->rmap->n*A->cmap->n;
2048:   PetscScalar  *aa = a->v;

2051:   for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
2052:   return(0);
2053: }

2055: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
2056: {
2057:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2058:   PetscInt     i,nz = A->rmap->n*A->cmap->n;
2059:   PetscScalar  *aa = a->v;

2062:   for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
2063:   return(0);
2064: }

2066: /* ----------------------------------------------------------------*/
2067: PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
2068: {

2072:   if (scall == MAT_INITIAL_MATRIX) {
2073:     PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
2074:     MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
2075:     PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
2076:   }
2077:   PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
2078:   MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);
2079:   PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
2080:   return(0);
2081: }

2083: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
2084: {
2086:   PetscInt       m=A->rmap->n,n=B->cmap->n;
2087:   Mat            Cmat;

2090:   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);
2091:   MatCreate(PETSC_COMM_SELF,&Cmat);
2092:   MatSetSizes(Cmat,m,n,m,n);
2093:   MatSetType(Cmat,MATSEQDENSE);
2094:   MatSeqDenseSetPreallocation(Cmat,NULL);

2096:   *C = Cmat;
2097:   return(0);
2098: }

2100: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2101: {
2102:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2103:   Mat_SeqDense   *b = (Mat_SeqDense*)B->data;
2104:   Mat_SeqDense   *c = (Mat_SeqDense*)C->data;
2105:   PetscBLASInt   m,n,k;
2106:   PetscScalar    _DOne=1.0,_DZero=0.0;
2108:   PetscBool      flg;

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

2114:   /* Handle case where where user provided the final C matrix rather than calling MatMatMult() with MAT_INITIAL_MATRIX*/
2115:   PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
2116:   if (flg) {
2117:     C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqDense;
2118:     (*C->ops->matmultnumeric)(A,B,C);
2119:     return(0);
2120:   }

2122:   PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&flg);
2123:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"First matrix must be dense");
2124:   PetscBLASIntCast(C->rmap->n,&m);
2125:   PetscBLASIntCast(C->cmap->n,&n);
2126:   PetscBLASIntCast(A->cmap->n,&k);
2127:   if (!m || !n || !k) return(0);
2128:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
2129:   return(0);
2130: }

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

2137:   if (scall == MAT_INITIAL_MATRIX) {
2138:     PetscLogEventBegin(MAT_MatTransposeMultSymbolic,A,B,0,0);
2139:     MatMatTransposeMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
2140:     PetscLogEventEnd(MAT_MatTransposeMultSymbolic,A,B,0,0);
2141:   }
2142:   PetscLogEventBegin(MAT_MatTransposeMultNumeric,A,B,0,0);
2143:   MatMatTransposeMultNumeric_SeqDense_SeqDense(A,B,*C);
2144:   PetscLogEventEnd(MAT_MatTransposeMultNumeric,A,B,0,0);
2145:   return(0);
2146: }

2148: PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
2149: {
2151:   PetscInt       m=A->rmap->n,n=B->rmap->n;
2152:   Mat            Cmat;

2155:   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);
2156:   MatCreate(PETSC_COMM_SELF,&Cmat);
2157:   MatSetSizes(Cmat,m,n,m,n);
2158:   MatSetType(Cmat,MATSEQDENSE);
2159:   MatSeqDenseSetPreallocation(Cmat,NULL);

2161:   Cmat->assembled = PETSC_TRUE;

2163:   *C = Cmat;
2164:   return(0);
2165: }

2167: PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2168: {
2169:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2170:   Mat_SeqDense   *b = (Mat_SeqDense*)B->data;
2171:   Mat_SeqDense   *c = (Mat_SeqDense*)C->data;
2172:   PetscBLASInt   m,n,k;
2173:   PetscScalar    _DOne=1.0,_DZero=0.0;

2177:   PetscBLASIntCast(C->rmap->n,&m);
2178:   PetscBLASIntCast(C->cmap->n,&n);
2179:   PetscBLASIntCast(A->cmap->n,&k);
2180:   if (!m || !n || !k) return(0);
2181:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
2182:   return(0);
2183: }

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

2190:   if (scall == MAT_INITIAL_MATRIX) {
2191:     PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);
2192:     MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
2193:     PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);
2194:   }
2195:   PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);
2196:   MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);
2197:   PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);
2198:   return(0);
2199: }

2201: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
2202: {
2204:   PetscInt       m=A->cmap->n,n=B->cmap->n;
2205:   Mat            Cmat;

2208:   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);
2209:   MatCreate(PETSC_COMM_SELF,&Cmat);
2210:   MatSetSizes(Cmat,m,n,m,n);
2211:   MatSetType(Cmat,MATSEQDENSE);
2212:   MatSeqDenseSetPreallocation(Cmat,NULL);

2214:   Cmat->assembled = PETSC_TRUE;

2216:   *C = Cmat;
2217:   return(0);
2218: }

2220: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2221: {
2222:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2223:   Mat_SeqDense   *b = (Mat_SeqDense*)B->data;
2224:   Mat_SeqDense   *c = (Mat_SeqDense*)C->data;
2225:   PetscBLASInt   m,n,k;
2226:   PetscScalar    _DOne=1.0,_DZero=0.0;

2230:   PetscBLASIntCast(C->rmap->n,&m);
2231:   PetscBLASIntCast(C->cmap->n,&n);
2232:   PetscBLASIntCast(A->rmap->n,&k);
2233:   if (!m || !n || !k) return(0);
2234:   PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
2235:   return(0);
2236: }

2238: static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[])
2239: {
2240:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2242:   PetscInt       i,j,m = A->rmap->n,n = A->cmap->n,p;
2243:   PetscScalar    *x;
2244:   MatScalar      *aa = a->v;

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

2249:   VecSet(v,0.0);
2250:   VecGetArray(v,&x);
2251:   VecGetLocalSize(v,&p);
2252:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2253:   for (i=0; i<m; i++) {
2254:     x[i] = aa[i]; if (idx) idx[i] = 0;
2255:     for (j=1; j<n; j++) {
2256:       if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
2257:     }
2258:   }
2259:   VecRestoreArray(v,&x);
2260:   return(0);
2261: }

2263: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])
2264: {
2265:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2267:   PetscInt       i,j,m = A->rmap->n,n = A->cmap->n,p;
2268:   PetscScalar    *x;
2269:   PetscReal      atmp;
2270:   MatScalar      *aa = a->v;

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

2275:   VecSet(v,0.0);
2276:   VecGetArray(v,&x);
2277:   VecGetLocalSize(v,&p);
2278:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2279:   for (i=0; i<m; i++) {
2280:     x[i] = PetscAbsScalar(aa[i]);
2281:     for (j=1; j<n; j++) {
2282:       atmp = PetscAbsScalar(aa[i+m*j]);
2283:       if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;}
2284:     }
2285:   }
2286:   VecRestoreArray(v,&x);
2287:   return(0);
2288: }

2290: static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])
2291: {
2292:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2294:   PetscInt       i,j,m = A->rmap->n,n = A->cmap->n,p;
2295:   PetscScalar    *x;
2296:   MatScalar      *aa = a->v;

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

2301:   VecSet(v,0.0);
2302:   VecGetArray(v,&x);
2303:   VecGetLocalSize(v,&p);
2304:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2305:   for (i=0; i<m; i++) {
2306:     x[i] = aa[i]; if (idx) idx[i] = 0;
2307:     for (j=1; j<n; j++) {
2308:       if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
2309:     }
2310:   }
2311:   VecRestoreArray(v,&x);
2312:   return(0);
2313: }

2315: static PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
2316: {
2317:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2319:   PetscScalar    *x;

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

2324:   VecGetArray(v,&x);
2325:   PetscArraycpy(x,a->v+col*a->lda,A->rmap->n);
2326:   VecRestoreArray(v,&x);
2327:   return(0);
2328: }

2330: PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms)
2331: {
2332:   PetscErrorCode    ierr;
2333:   PetscInt          i,j,m,n;
2334:   const PetscScalar *a;

2337:   MatGetSize(A,&m,&n);
2338:   PetscArrayzero(norms,n);
2339:   MatDenseGetArrayRead(A,&a);
2340:   if (type == NORM_2) {
2341:     for (i=0; i<n; i++) {
2342:       for (j=0; j<m; j++) {
2343:         norms[i] += PetscAbsScalar(a[j]*a[j]);
2344:       }
2345:       a += m;
2346:     }
2347:   } else if (type == NORM_1) {
2348:     for (i=0; i<n; i++) {
2349:       for (j=0; j<m; j++) {
2350:         norms[i] += PetscAbsScalar(a[j]);
2351:       }
2352:       a += m;
2353:     }
2354:   } else if (type == NORM_INFINITY) {
2355:     for (i=0; i<n; i++) {
2356:       for (j=0; j<m; j++) {
2357:         norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]);
2358:       }
2359:       a += m;
2360:     }
2361:   } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType");
2362:   MatDenseRestoreArrayRead(A,&a);
2363:   if (type == NORM_2) {
2364:     for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]);
2365:   }
2366:   return(0);
2367: }

2369: static PetscErrorCode  MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
2370: {
2372:   PetscScalar    *a;
2373:   PetscInt       m,n,i;

2376:   MatGetSize(x,&m,&n);
2377:   MatDenseGetArray(x,&a);
2378:   for (i=0; i<m*n; i++) {
2379:     PetscRandomGetValue(rctx,a+i);
2380:   }
2381:   MatDenseRestoreArray(x,&a);
2382:   return(0);
2383: }

2385: static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool  *missing,PetscInt *d)
2386: {
2388:   *missing = PETSC_FALSE;
2389:   return(0);
2390: }

2392: static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals)
2393: {
2394:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;

2397:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2398:   *vals = a->v+col*a->lda;
2399:   return(0);
2400: }

2402: static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals)
2403: {
2405:   *vals = 0; /* user cannot accidently use the array later */
2406:   return(0);
2407: }

2409: /* -------------------------------------------------------------------*/
2410: static struct _MatOps MatOps_Values = { MatSetValues_SeqDense,
2411:                                         MatGetRow_SeqDense,
2412:                                         MatRestoreRow_SeqDense,
2413:                                         MatMult_SeqDense,
2414:                                 /*  4*/ MatMultAdd_SeqDense,
2415:                                         MatMultTranspose_SeqDense,
2416:                                         MatMultTransposeAdd_SeqDense,
2417:                                         0,
2418:                                         0,
2419:                                         0,
2420:                                 /* 10*/ 0,
2421:                                         MatLUFactor_SeqDense,
2422:                                         MatCholeskyFactor_SeqDense,
2423:                                         MatSOR_SeqDense,
2424:                                         MatTranspose_SeqDense,
2425:                                 /* 15*/ MatGetInfo_SeqDense,
2426:                                         MatEqual_SeqDense,
2427:                                         MatGetDiagonal_SeqDense,
2428:                                         MatDiagonalScale_SeqDense,
2429:                                         MatNorm_SeqDense,
2430:                                 /* 20*/ MatAssemblyBegin_SeqDense,
2431:                                         MatAssemblyEnd_SeqDense,
2432:                                         MatSetOption_SeqDense,
2433:                                         MatZeroEntries_SeqDense,
2434:                                 /* 24*/ MatZeroRows_SeqDense,
2435:                                         0,
2436:                                         0,
2437:                                         0,
2438:                                         0,
2439:                                 /* 29*/ MatSetUp_SeqDense,
2440:                                         0,
2441:                                         0,
2442:                                         0,
2443:                                         0,
2444:                                 /* 34*/ MatDuplicate_SeqDense,
2445:                                         0,
2446:                                         0,
2447:                                         0,
2448:                                         0,
2449:                                 /* 39*/ MatAXPY_SeqDense,
2450:                                         MatCreateSubMatrices_SeqDense,
2451:                                         0,
2452:                                         MatGetValues_SeqDense,
2453:                                         MatCopy_SeqDense,
2454:                                 /* 44*/ MatGetRowMax_SeqDense,
2455:                                         MatScale_SeqDense,
2456:                                         MatShift_Basic,
2457:                                         0,
2458:                                         MatZeroRowsColumns_SeqDense,
2459:                                 /* 49*/ MatSetRandom_SeqDense,
2460:                                         0,
2461:                                         0,
2462:                                         0,
2463:                                         0,
2464:                                 /* 54*/ 0,
2465:                                         0,
2466:                                         0,
2467:                                         0,
2468:                                         0,
2469:                                 /* 59*/ 0,
2470:                                         MatDestroy_SeqDense,
2471:                                         MatView_SeqDense,
2472:                                         0,
2473:                                         0,
2474:                                 /* 64*/ 0,
2475:                                         0,
2476:                                         0,
2477:                                         0,
2478:                                         0,
2479:                                 /* 69*/ MatGetRowMaxAbs_SeqDense,
2480:                                         0,
2481:                                         0,
2482:                                         0,
2483:                                         0,
2484:                                 /* 74*/ 0,
2485:                                         0,
2486:                                         0,
2487:                                         0,
2488:                                         0,
2489:                                 /* 79*/ 0,
2490:                                         0,
2491:                                         0,
2492:                                         0,
2493:                                 /* 83*/ MatLoad_SeqDense,
2494:                                         0,
2495:                                         MatIsHermitian_SeqDense,
2496:                                         0,
2497:                                         0,
2498:                                         0,
2499:                                 /* 89*/ MatMatMult_SeqDense_SeqDense,
2500:                                         MatMatMultSymbolic_SeqDense_SeqDense,
2501:                                         MatMatMultNumeric_SeqDense_SeqDense,
2502:                                         MatPtAP_SeqDense_SeqDense,
2503:                                         MatPtAPSymbolic_SeqDense_SeqDense,
2504:                                 /* 94*/ MatPtAPNumeric_SeqDense_SeqDense,
2505:                                         MatMatTransposeMult_SeqDense_SeqDense,
2506:                                         MatMatTransposeMultSymbolic_SeqDense_SeqDense,
2507:                                         MatMatTransposeMultNumeric_SeqDense_SeqDense,
2508:                                         0,
2509:                                 /* 99*/ 0,
2510:                                         0,
2511:                                         0,
2512:                                         MatConjugate_SeqDense,
2513:                                         0,
2514:                                 /*104*/ 0,
2515:                                         MatRealPart_SeqDense,
2516:                                         MatImaginaryPart_SeqDense,
2517:                                         0,
2518:                                         0,
2519:                                 /*109*/ 0,
2520:                                         0,
2521:                                         MatGetRowMin_SeqDense,
2522:                                         MatGetColumnVector_SeqDense,
2523:                                         MatMissingDiagonal_SeqDense,
2524:                                 /*114*/ 0,
2525:                                         0,
2526:                                         0,
2527:                                         0,
2528:                                         0,
2529:                                 /*119*/ 0,
2530:                                         0,
2531:                                         0,
2532:                                         0,
2533:                                         0,
2534:                                 /*124*/ 0,
2535:                                         MatGetColumnNorms_SeqDense,
2536:                                         0,
2537:                                         0,
2538:                                         0,
2539:                                 /*129*/ 0,
2540:                                         MatTransposeMatMult_SeqDense_SeqDense,
2541:                                         MatTransposeMatMultSymbolic_SeqDense_SeqDense,
2542:                                         MatTransposeMatMultNumeric_SeqDense_SeqDense,
2543:                                         0,
2544:                                 /*134*/ 0,
2545:                                         0,
2546:                                         0,
2547:                                         0,
2548:                                         0,
2549:                                 /*139*/ 0,
2550:                                         0,
2551:                                         0,
2552:                                         0,
2553:                                         0,
2554:                                 /*144*/ MatCreateMPIMatConcatenateSeqMat_SeqDense
2555: };

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

2562:    Collective

2564:    Input Parameters:
2565: +  comm - MPI communicator, set to PETSC_COMM_SELF
2566: .  m - number of rows
2567: .  n - number of columns
2568: -  data - optional location of matrix data in column major order.  Set data=NULL for PETSc
2569:    to control all matrix memory allocation.

2571:    Output Parameter:
2572: .  A - the matrix

2574:    Notes:
2575:    The data input variable is intended primarily for Fortran programmers
2576:    who wish to allocate their own matrix memory space.  Most users should
2577:    set data=NULL.

2579:    Level: intermediate

2581: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
2582: @*/
2583: PetscErrorCode  MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
2584: {

2588:   MatCreate(comm,A);
2589:   MatSetSizes(*A,m,n,m,n);
2590:   MatSetType(*A,MATSEQDENSE);
2591:   MatSeqDenseSetPreallocation(*A,data);
2592:   return(0);
2593: }

2595: /*@C
2596:    MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements

2598:    Collective

2600:    Input Parameters:
2601: +  B - the matrix
2602: -  data - the array (or NULL)

2604:    Notes:
2605:    The data input variable is intended primarily for Fortran programmers
2606:    who wish to allocate their own matrix memory space.  Most users should
2607:    need not call this routine.

2609:    Level: intermediate

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

2613: @*/
2614: PetscErrorCode  MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
2615: {

2619:   PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
2620:   return(0);
2621: }

2623: PetscErrorCode  MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
2624: {
2625:   Mat_SeqDense   *b;

2629:   B->preallocated = PETSC_TRUE;

2631:   PetscLayoutSetUp(B->rmap);
2632:   PetscLayoutSetUp(B->cmap);

2634:   b       = (Mat_SeqDense*)B->data;
2635:   b->Mmax = B->rmap->n;
2636:   b->Nmax = B->cmap->n;
2637:   if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n;

2639:   PetscIntMultError(b->lda,b->Nmax,NULL);
2640:   if (!data) { /* petsc-allocated storage */
2641:     if (!b->user_alloc) { PetscFree(b->v); }
2642:     PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);
2643:     PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));

2645:     b->user_alloc = PETSC_FALSE;
2646:   } else { /* user-allocated storage */
2647:     if (!b->user_alloc) { PetscFree(b->v); }
2648:     b->v          = data;
2649:     b->user_alloc = PETSC_TRUE;
2650:   }
2651:   B->assembled = PETSC_TRUE;
2652:   return(0);
2653: }

2655: #if defined(PETSC_HAVE_ELEMENTAL)
2656: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
2657: {
2658:   Mat               mat_elemental;
2659:   PetscErrorCode    ierr;
2660:   const PetscScalar *array;
2661:   PetscScalar       *v_colwise;
2662:   PetscInt          M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols;

2665:   PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);
2666:   MatDenseGetArrayRead(A,&array);
2667:   /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
2668:   k = 0;
2669:   for (j=0; j<N; j++) {
2670:     cols[j] = j;
2671:     for (i=0; i<M; i++) {
2672:       v_colwise[j*M+i] = array[k++];
2673:     }
2674:   }
2675:   for (i=0; i<M; i++) {
2676:     rows[i] = i;
2677:   }
2678:   MatDenseRestoreArrayRead(A,&array);

2680:   MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
2681:   MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);
2682:   MatSetType(mat_elemental,MATELEMENTAL);
2683:   MatSetUp(mat_elemental);

2685:   /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
2686:   MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);
2687:   MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
2688:   MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);
2689:   PetscFree3(v_colwise,rows,cols);

2691:   if (reuse == MAT_INPLACE_MATRIX) {
2692:     MatHeaderReplace(A,&mat_elemental);
2693:   } else {
2694:     *newmat = mat_elemental;
2695:   }
2696:   return(0);
2697: }
2698: #endif

2700: /*@C
2701:   MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array

2703:   Input parameter:
2704: + A - the matrix
2705: - lda - the leading dimension

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

2712:   Level: intermediate

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

2716: @*/
2717: PetscErrorCode  MatSeqDenseSetLDA(Mat B,PetscInt lda)
2718: {
2719:   Mat_SeqDense *b = (Mat_SeqDense*)B->data;

2722:   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);
2723:   b->lda       = lda;
2724:   b->changelda = PETSC_FALSE;
2725:   b->Mmax      = PetscMax(b->Mmax,lda);
2726:   return(0);
2727: }

2729: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat)
2730: {
2732:   PetscMPIInt    size;

2735:   MPI_Comm_size(comm,&size);
2736:   if (size == 1) {
2737:     if (scall == MAT_INITIAL_MATRIX) {
2738:       MatDuplicate(inmat,MAT_COPY_VALUES,outmat);
2739:     } else {
2740:       MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);
2741:     }
2742:   } else {
2743:     MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);
2744:   }
2745:   return(0);
2746: }

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

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

2754:   Level: beginner

2756: .seealso: MatCreateSeqDense()

2758: M*/

2760: PETSC_EXTERN PetscErrorCode MatCreate_SeqDense(Mat B)
2761: {
2762:   Mat_SeqDense   *b;
2764:   PetscMPIInt    size;

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

2770:   PetscNewLog(B,&b);
2771:   PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
2772:   B->data = (void*)b;

2774:   b->roworiented = PETSC_TRUE;

2776:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetLDA_C",MatDenseGetLDA_SeqDense);
2777:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);
2778:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);
2779:   PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);
2780:   PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);
2781:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);
2782:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense);
2783:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);
2784: #if defined(PETSC_HAVE_ELEMENTAL)
2785:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);
2786: #endif
2787:   PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);
2788:   PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2789:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2790:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2791:   PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaij_seqdense_C",MatPtAP_SeqDense_SeqDense);
2792:   PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijperm_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2793:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijperm_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2794:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijperm_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2795:   PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijperm_seqdense_C",MatPtAP_SeqDense_SeqDense);
2796:   PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijsell_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2797:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijsell_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2798:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijsell_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2799:   PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijsell_seqdense_C",MatPtAP_SeqDense_SeqDense);
2800:   PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijmkl_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2801:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijmkl_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2802:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijmkl_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2803:   PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijmkl_seqdense_C",MatPtAP_SeqDense_SeqDense);

2805:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2806:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2807:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2808:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijperm_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2809:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijperm_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2810:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijperm_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2811:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijsell_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2812:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijsell_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2813:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijsell_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);

2815:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijmkl_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2816:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijmkl_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2817:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijmkl_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2818:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);
2819:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);
2820:   PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
2821:   return(0);
2822: }

2824: /*@C
2825:    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.

2827:    Not Collective

2829:    Input Parameter:
2830: +  mat - a MATSEQDENSE or MATMPIDENSE matrix
2831: -  col - column index

2833:    Output Parameter:
2834: .  vals - pointer to the data

2836:    Level: intermediate

2838: .seealso: MatDenseRestoreColumn()
2839: @*/
2840: PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals)
2841: {

2845:   PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));
2846:   return(0);
2847: }

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

2852:    Not Collective

2854:    Input Parameter:
2855: .  mat - a MATSEQDENSE or MATMPIDENSE matrix

2857:    Output Parameter:
2858: .  vals - pointer to the data

2860:    Level: intermediate

2862: .seealso: MatDenseGetColumn()
2863: @*/
2864: PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals)
2865: {

2869:   PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));
2870:   return(0);
2871: }