Actual source code: dense.c

petsc-master 2017-06-23
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:     PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
 55:     PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0); /* TODO CHECK FLOPS */
 56:   } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
 57:     if (A->spd) {
 58:       PetscStackCallBLAS("LAPACKpotri",LAPACKpotri_("L",&n,mat->v,&mat->lda,&info));
 59:       MatSeqDenseSymmetrize_Private(A,PETSC_TRUE);
 60: #if defined (PETSC_USE_COMPLEX)
 61:     } else if (A->hermitian) {
 62:       if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present");
 63:       if (!mat->fwork) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Fwork not present");
 64:       PetscStackCallBLAS("LAPACKhetri",LAPACKhetri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info));
 65:       MatSeqDenseSymmetrize_Private(A,PETSC_TRUE);
 66: #endif
 67:     } else { /* symmetric case */
 68:       if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present");
 69:       if (!mat->fwork) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Fwork not present");
 70:       PetscStackCallBLAS("LAPACKsytri",LAPACKsytri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info));
 71:       MatSeqDenseSymmetrize_Private(A,PETSC_FALSE);
 72:     }
 73:     if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad Inversion: zero pivot in row %D",(PetscInt)info-1);
 74:     PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0); /* TODO CHECK FLOPS */
 75:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
 76: #endif

 78:   A->ops->solve             = NULL;
 79:   A->ops->matsolve          = NULL;
 80:   A->ops->solvetranspose    = NULL;
 81:   A->ops->matsolvetranspose = NULL;
 82:   A->ops->solveadd          = NULL;
 83:   A->ops->solvetransposeadd = NULL;
 84:   A->factortype             = MAT_FACTOR_NONE;
 85:   PetscFree(A->solvertype);
 86:   return(0);
 87: }

 89: PetscErrorCode MatZeroRowsColumns_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
 90: {
 91:   PetscErrorCode    ierr;
 92:   Mat_SeqDense      *l = (Mat_SeqDense*)A->data;
 93:   PetscInt          m  = l->lda, n = A->cmap->n,r = A->rmap->n, i,j;
 94:   PetscScalar       *slot,*bb;
 95:   const PetscScalar *xx;

 98: #if defined(PETSC_USE_DEBUG)
 99:   for (i=0; i<N; i++) {
100:     if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
101:     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);
102:     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);
103:   }
104: #endif

106:   /* fix right hand side if needed */
107:   if (x && b) {
108:     VecGetArrayRead(x,&xx);
109:     VecGetArray(b,&bb);
110:     for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
111:     VecRestoreArrayRead(x,&xx);
112:     VecRestoreArray(b,&bb);
113:   }

115:   for (i=0; i<N; i++) {
116:     slot = l->v + rows[i]*m;
117:     PetscMemzero(slot,r*sizeof(PetscScalar));
118:   }
119:   for (i=0; i<N; i++) {
120:     slot = l->v + rows[i];
121:     for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
122:   }
123:   if (diag != 0.0) {
124:     if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
125:     for (i=0; i<N; i++) {
126:       slot  = l->v + (m+1)*rows[i];
127:       *slot = diag;
128:     }
129:   }
130:   return(0);
131: }

133: PetscErrorCode MatPtAPNumeric_SeqDense_SeqDense(Mat A,Mat P,Mat C)
134: {
135:   Mat_SeqDense   *c = (Mat_SeqDense*)(C->data);

139:   MatMatMultNumeric_SeqDense_SeqDense(A,P,c->ptapwork);
140:   MatTransposeMatMultNumeric_SeqDense_SeqDense(P,c->ptapwork,C);
141:   return(0);
142: }

144: PetscErrorCode MatPtAPSymbolic_SeqDense_SeqDense(Mat A,Mat P,PetscReal fill,Mat *C)
145: {
146:   Mat_SeqDense   *c;

150:   MatCreateSeqDense(PetscObjectComm((PetscObject)A),P->cmap->N,P->cmap->N,NULL,C);
151:   c = (Mat_SeqDense*)((*C)->data);
152:   MatCreateSeqDense(PetscObjectComm((PetscObject)A),A->rmap->N,P->cmap->N,NULL,&c->ptapwork);
153:   return(0);
154: }

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

161:   if (reuse == MAT_INITIAL_MATRIX) {
162:     MatPtAPSymbolic_SeqDense_SeqDense(A,P,fill,C);
163:   }
164:   PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);
165:   (*(*C)->ops->ptapnumeric)(A,P,*C);
166:   PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);
167:   return(0);
168: }

170: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
171: {
172:   Mat            B = NULL;
173:   Mat_SeqAIJ     *a = (Mat_SeqAIJ*)A->data;
174:   Mat_SeqDense   *b;
176:   PetscInt       *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->N,i;
177:   MatScalar      *av=a->a;
178:   PetscBool      isseqdense;

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

205:   if (reuse == MAT_INPLACE_MATRIX) {
206:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
207:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
208:     MatHeaderReplace(A,&B);
209:   } else {
210:     if (B) *newmat = B;
211:     MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
212:     MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
213:   }
214:   return(0);
215: }

217: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
218: {
219:   Mat            B;
220:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
222:   PetscInt       i, j;
223:   PetscInt       *rows, *nnz;
224:   MatScalar      *aa = a->v, *vals;

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

247:   if (reuse == MAT_INPLACE_MATRIX) {
248:     MatHeaderReplace(A,&B);
249:   } else {
250:     *newmat = B;
251:   }
252:   return(0);
253: }

255: static PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str)
256: {
257:   Mat_SeqDense   *x     = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data;
258:   PetscScalar    oalpha = alpha;
259:   PetscInt       j;
260:   PetscBLASInt   N,m,ldax,lday,one = 1;

264:   PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);
265:   PetscBLASIntCast(X->rmap->n,&m);
266:   PetscBLASIntCast(x->lda,&ldax);
267:   PetscBLASIntCast(y->lda,&lday);
268:   if (ldax>m || lday>m) {
269:     for (j=0; j<X->cmap->n; j++) {
270:       PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&oalpha,x->v+j*ldax,&one,y->v+j*lday,&one));
271:     }
272:   } else {
273:     PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&oalpha,x->v,&one,y->v,&one));
274:   }
275:   PetscObjectStateIncrease((PetscObject)Y);
276:   PetscLogFlops(PetscMax(2*N-1,0));
277:   return(0);
278: }

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

285:   info->block_size        = 1.0;
286:   info->nz_allocated      = (double)N;
287:   info->nz_used           = (double)N;
288:   info->nz_unneeded       = (double)0;
289:   info->assemblies        = (double)A->num_ass;
290:   info->mallocs           = 0;
291:   info->memory            = ((PetscObject)A)->mem;
292:   info->fill_ratio_given  = 0;
293:   info->fill_ratio_needed = 0;
294:   info->factor_mallocs    = 0;
295:   return(0);
296: }

298: static PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha)
299: {
300:   Mat_SeqDense   *a     = (Mat_SeqDense*)A->data;
301:   PetscScalar    oalpha = alpha;
303:   PetscBLASInt   one = 1,j,nz,lda;

306:   PetscBLASIntCast(a->lda,&lda);
307:   if (lda>A->rmap->n) {
308:     PetscBLASIntCast(A->rmap->n,&nz);
309:     for (j=0; j<A->cmap->n; j++) {
310:       PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v+j*lda,&one));
311:     }
312:   } else {
313:     PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);
314:     PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v,&one));
315:   }
316:   PetscLogFlops(nz);
317:   return(0);
318: }

320: static PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool  *fl)
321: {
322:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
323:   PetscInt     i,j,m = A->rmap->n,N;
324:   PetscScalar  *v = a->v;

327:   *fl = PETSC_FALSE;
328:   if (A->rmap->n != A->cmap->n) return(0);
329:   N = a->lda;

331:   for (i=0; i<m; i++) {
332:     for (j=i+1; j<m; j++) {
333:       if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) return(0);
334:     }
335:   }
336:   *fl = PETSC_TRUE;
337:   return(0);
338: }

340: static PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues)
341: {
342:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data,*l;
344:   PetscInt       lda = (PetscInt)mat->lda,j,m;

347:   PetscLayoutReference(A->rmap,&newi->rmap);
348:   PetscLayoutReference(A->cmap,&newi->cmap);
349:   MatSeqDenseSetPreallocation(newi,NULL);
350:   if (cpvalues == MAT_COPY_VALUES) {
351:     l = (Mat_SeqDense*)newi->data;
352:     if (lda>A->rmap->n) {
353:       m = A->rmap->n;
354:       for (j=0; j<A->cmap->n; j++) {
355:         PetscMemcpy(l->v+j*m,mat->v+j*lda,m*sizeof(PetscScalar));
356:       }
357:     } else {
358:       PetscMemcpy(l->v,mat->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
359:     }
360:   }
361:   newi->assembled = PETSC_TRUE;
362:   return(0);
363: }

365: static PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat)
366: {

370:   MatCreate(PetscObjectComm((PetscObject)A),newmat);
371:   MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
372:   MatSetType(*newmat,((PetscObject)A)->type_name);
373:   MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);
374:   return(0);
375: }


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

380: static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
381: {
382:   MatFactorInfo  info;

386:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
387:   MatLUFactor_SeqDense(fact,0,0,&info);
388:   return(0);
389: }

391: static PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy)
392: {
393:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
394:   PetscErrorCode    ierr;
395:   const PetscScalar *x;
396:   PetscScalar       *y;
397:   PetscBLASInt      one = 1,info,m;

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

435: static PetscErrorCode MatMatSolve_SeqDense(Mat A,Mat B,Mat X)
436: {
437:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
439:   PetscScalar    *b,*x;
440:   PetscInt       n;
441:   PetscBLASInt   nrhs,info,m;
442:   PetscBool      flg;

445:   PetscBLASIntCast(A->rmap->n,&m);
446:   PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
447:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
448:   PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
449:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");

451:   MatGetSize(B,NULL,&n);
452:   PetscBLASIntCast(n,&nrhs);
453:   MatDenseGetArray(B,&b);
454:   MatDenseGetArray(X,&x);

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

458:   if (A->factortype == MAT_FACTOR_LU) {
459: #if defined(PETSC_MISSING_LAPACK_GETRS)
460:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
461: #else
462:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
463:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
464: #endif
465:   } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
466: #if defined(PETSC_MISSING_LAPACK_POTRS)
467:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
468: #else
469:     if (A->spd) {
470:       PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info));
471:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
472: #if defined (PETSC_USE_COMPLEX)
473:     } else if (A->hermitian) {
474:       PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
475:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve");
476: #endif
477:     } else { /* symmetric case */
478:       PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
479:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
480:     }
481: #endif
482:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");

484:   MatDenseRestoreArray(B,&b);
485:   MatDenseRestoreArray(X,&x);
486:   PetscLogFlops(nrhs*(2.0*m*m - m));
487:   return(0);
488: }

490: static PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy)
491: {
492:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
493:   PetscErrorCode    ierr;
494:   const PetscScalar *x;
495:   PetscScalar       *y;
496:   PetscBLASInt      one = 1,info,m;

499:   PetscBLASIntCast(A->rmap->n,&m);
500:   VecGetArrayRead(xx,&x);
501:   VecGetArray(yy,&y);
502:   PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
503:   if (A->factortype == MAT_FACTOR_LU) {
504: #if defined(PETSC_MISSING_LAPACK_GETRS)
505:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
506: #else
507:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
508:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve");
509: #endif
510:   } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
511: #if defined(PETSC_MISSING_LAPACK_POTRS)
512:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
513: #else
514:     if (A->spd) {
515:       PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
516:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
517: #if defined (PETSC_USE_COMPLEX)
518:     } else if (A->hermitian) {
519:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatSolveTranspose with Cholesky/Hermitian not available");
520: #endif
521:     } else { /* symmetric case */
522:       PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
523:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
524:     }
525: #endif
526:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
527:   VecRestoreArrayRead(xx,&x);
528:   VecRestoreArray(yy,&y);
529:   PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
530:   return(0);
531: }

533: static PetscErrorCode MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
534: {
535:   PetscErrorCode    ierr;
536:   const PetscScalar *x;
537:   PetscScalar       *y,sone = 1.0;
538:   Vec               tmp = 0;

541:   VecGetArrayRead(xx,&x);
542:   VecGetArray(yy,&y);
543:   if (!A->rmap->n || !A->cmap->n) return(0);
544:   if (yy == zz) {
545:     VecDuplicate(yy,&tmp);
546:     PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);
547:     VecCopy(yy,tmp);
548:   }
549:   MatSolve_SeqDense(A,xx,yy);
550:   if (tmp) {
551:     VecAXPY(yy,sone,tmp);
552:     VecDestroy(&tmp);
553:   } else {
554:     VecAXPY(yy,sone,zz);
555:   }
556:   VecRestoreArrayRead(xx,&x);
557:   VecRestoreArray(yy,&y);
558:   PetscLogFlops(A->cmap->n);
559:   return(0);
560: }

562: static PetscErrorCode MatSolveTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
563: {
564:   PetscErrorCode    ierr;
565:   const PetscScalar *x;
566:   PetscScalar       *y,sone = 1.0;
567:   Vec               tmp = NULL;

570:   if (!A->rmap->n || !A->cmap->n) return(0);
571:   VecGetArrayRead(xx,&x);
572:   VecGetArray(yy,&y);
573:   if (yy == zz) {
574:     VecDuplicate(yy,&tmp);
575:     PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);
576:     VecCopy(yy,tmp);
577:   }
578:   MatSolveTranspose_SeqDense(A,xx,yy);
579:   if (tmp) {
580:     VecAXPY(yy,sone,tmp);
581:     VecDestroy(&tmp);
582:   } else {
583:     VecAXPY(yy,sone,zz);
584:   }
585:   VecRestoreArrayRead(xx,&x);
586:   VecRestoreArray(yy,&y);
587:   PetscLogFlops(A->cmap->n);
588:   return(0);
589: }

591: /* ---------------------------------------------------------------*/
592: /* COMMENT: I have chosen to hide row permutation in the pivots,
593:    rather than put it in the Mat->row slot.*/
594: static PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo)
595: {
596: #if defined(PETSC_MISSING_LAPACK_GETRF)
598:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable.");
599: #else
600:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
602:   PetscBLASInt   n,m,info;

605:   PetscBLASIntCast(A->cmap->n,&n);
606:   PetscBLASIntCast(A->rmap->n,&m);
607:   if (!mat->pivots) {
608:     PetscMalloc1(A->rmap->n,&mat->pivots);
609:     PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
610:   }
611:   if (!A->rmap->n || !A->cmap->n) return(0);
612:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
613:   PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info));
614:   PetscFPTrapPop();

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

619:   A->ops->solve             = MatSolve_SeqDense;
620:   A->ops->matsolve          = MatMatSolve_SeqDense;
621:   A->ops->solvetranspose    = MatSolveTranspose_SeqDense;
622:   A->ops->solveadd          = MatSolveAdd_SeqDense;
623:   A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
624:   A->factortype             = MAT_FACTOR_LU;

626:   PetscFree(A->solvertype);
627:   PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);

629:   PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);
630: #endif
631:   return(0);
632: }

634: /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */
635: static PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo)
636: {
637: #if defined(PETSC_MISSING_LAPACK_POTRF)
639:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable.");
640: #else
641:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
643:   PetscBLASInt   info,n;

646:   PetscBLASIntCast(A->cmap->n,&n);
647:   if (!A->rmap->n || !A->cmap->n) return(0);
648:   if (A->spd) {
649:     PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info));
650: #if defined (PETSC_USE_COMPLEX)
651:   } else if (A->hermitian) {
652:     if (!mat->pivots) {
653:       PetscMalloc1(A->rmap->n,&mat->pivots);
654:       PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
655:     }
656:     if (!mat->fwork) {
657:       PetscScalar dummy;

659:       mat->lfwork = -1;
660:       PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
661:       mat->lfwork = (PetscInt)PetscRealPart(dummy);
662:       PetscMalloc1(mat->lfwork,&mat->fwork);
663:       PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
664:     }
665:     PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
666: #endif
667:   } else { /* symmetric case */
668:     if (!mat->pivots) {
669:       PetscMalloc1(A->rmap->n,&mat->pivots);
670:       PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
671:     }
672:     if (!mat->fwork) {
673:       PetscScalar dummy;

675:       mat->lfwork = -1;
676:       PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
677:       mat->lfwork = (PetscInt)PetscRealPart(dummy);
678:       PetscMalloc1(mat->lfwork,&mat->fwork);
679:       PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
680:     }
681:     PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
682:   }
683:   if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1);

685:   A->ops->solve             = MatSolve_SeqDense;
686:   A->ops->matsolve          = MatMatSolve_SeqDense;
687:   A->ops->solvetranspose    = MatSolveTranspose_SeqDense;
688:   A->ops->solveadd          = MatSolveAdd_SeqDense;
689:   A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
690:   A->factortype             = MAT_FACTOR_CHOLESKY;

692:   PetscFree(A->solvertype);
693:   PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);

695:   PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
696: #endif
697:   return(0);
698: }


701: PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
702: {
704:   MatFactorInfo  info;

707:   info.fill = 1.0;

709:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
710:   MatCholeskyFactor_SeqDense(fact,0,&info);
711:   return(0);
712: }

714: static PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
715: {
717:   fact->assembled                  = PETSC_TRUE;
718:   fact->preallocated               = PETSC_TRUE;
719:   fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
720:   fact->ops->solve                 = MatSolve_SeqDense;
721:   fact->ops->matsolve              = MatMatSolve_SeqDense;
722:   fact->ops->solvetranspose        = MatSolveTranspose_SeqDense;
723:   fact->ops->solveadd              = MatSolveAdd_SeqDense;
724:   fact->ops->solvetransposeadd     = MatSolveTransposeAdd_SeqDense;
725:   return(0);
726: }

728: static PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info)
729: {
731:   fact->preallocated           = PETSC_TRUE;
732:   fact->assembled              = PETSC_TRUE;
733:   fact->ops->lufactornumeric   = MatLUFactorNumeric_SeqDense;
734:   fact->ops->solve             = MatSolve_SeqDense;
735:   fact->ops->matsolve          = MatMatSolve_SeqDense;
736:   fact->ops->solvetranspose    = MatSolveTranspose_SeqDense;
737:   fact->ops->solveadd          = MatSolveAdd_SeqDense;
738:   fact->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
739:   return(0);
740: }

742: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact)
743: {

747:   MatCreate(PetscObjectComm((PetscObject)A),fact);
748:   MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
749:   MatSetType(*fact,((PetscObject)A)->type_name);
750:   if (ftype == MAT_FACTOR_LU) {
751:     (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
752:   } else {
753:     (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
754:   }
755:   (*fact)->factortype = ftype;

757:   PetscFree((*fact)->solvertype);
758:   PetscStrallocpy(MATSOLVERPETSC,&(*fact)->solvertype);
759:   return(0);
760: }

762: /* ------------------------------------------------------------------*/
763: static PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx)
764: {
765:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
766:   PetscScalar       *x,*v = mat->v,zero = 0.0,xt;
767:   const PetscScalar *b;
768:   PetscErrorCode    ierr;
769:   PetscInt          m = A->rmap->n,i;
770:   PetscBLASInt      o = 1,bm;

773:   if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */
774:   PetscBLASIntCast(m,&bm);
775:   if (flag & SOR_ZERO_INITIAL_GUESS) {
776:     /* this is a hack fix, should have another version without the second BLASdot */
777:     VecSet(xx,zero);
778:   }
779:   VecGetArray(xx,&x);
780:   VecGetArrayRead(bb,&b);
781:   its  = its*lits;
782:   if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
783:   while (its--) {
784:     if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
785:       for (i=0; i<m; i++) {
786:         PetscStackCallBLAS("BLASdot",xt   = b[i] - BLASdot_(&bm,v+i,&bm,x,&o));
787:         x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
788:       }
789:     }
790:     if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
791:       for (i=m-1; i>=0; i--) {
792:         PetscStackCallBLAS("BLASdot",xt   = b[i] - BLASdot_(&bm,v+i,&bm,x,&o));
793:         x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
794:       }
795:     }
796:   }
797:   VecRestoreArrayRead(bb,&b);
798:   VecRestoreArray(xx,&x);
799:   return(0);
800: }

802: /* -----------------------------------------------------------------*/
803: PETSC_INTERN PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)
804: {
805:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
806:   const PetscScalar *v   = mat->v,*x;
807:   PetscScalar       *y;
808:   PetscErrorCode    ierr;
809:   PetscBLASInt      m, n,_One=1;
810:   PetscScalar       _DOne=1.0,_DZero=0.0;

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

825: PETSC_INTERN PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy)
826: {
827:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
828:   PetscScalar       *y,_DOne=1.0,_DZero=0.0;
829:   PetscErrorCode    ierr;
830:   PetscBLASInt      m, n, _One=1;
831:   const PetscScalar *v = mat->v,*x;

834:   PetscBLASIntCast(A->rmap->n,&m);
835:   PetscBLASIntCast(A->cmap->n,&n);
836:   if (!A->rmap->n || !A->cmap->n) return(0);
837:   VecGetArrayRead(xx,&x);
838:   VecGetArray(yy,&y);
839:   PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One));
840:   VecRestoreArrayRead(xx,&x);
841:   VecRestoreArray(yy,&y);
842:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);
843:   return(0);
844: }

846: PETSC_INTERN PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
847: {
848:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
849:   const PetscScalar *v = mat->v,*x;
850:   PetscScalar       *y,_DOne=1.0;
851:   PetscErrorCode    ierr;
852:   PetscBLASInt      m, n, _One=1;

855:   PetscBLASIntCast(A->rmap->n,&m);
856:   PetscBLASIntCast(A->cmap->n,&n);
857:   if (!A->rmap->n || !A->cmap->n) return(0);
858:   if (zz != yy) {VecCopy(zz,yy);}
859:   VecGetArrayRead(xx,&x);
860:   VecGetArray(yy,&y);
861:   PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
862:   VecRestoreArrayRead(xx,&x);
863:   VecRestoreArray(yy,&y);
864:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
865:   return(0);
866: }

868: static PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
869: {
870:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
871:   const PetscScalar *v = mat->v,*x;
872:   PetscScalar       *y;
873:   PetscErrorCode    ierr;
874:   PetscBLASInt      m, n, _One=1;
875:   PetscScalar       _DOne=1.0;

878:   PetscBLASIntCast(A->rmap->n,&m);
879:   PetscBLASIntCast(A->cmap->n,&n);
880:   if (!A->rmap->n || !A->cmap->n) return(0);
881:   if (zz != yy) {VecCopy(zz,yy);}
882:   VecGetArrayRead(xx,&x);
883:   VecGetArray(yy,&y);
884:   PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
885:   VecRestoreArrayRead(xx,&x);
886:   VecRestoreArray(yy,&y);
887:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
888:   return(0);
889: }

891: /* -----------------------------------------------------------------*/
892: static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
893: {
894:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
895:   PetscScalar    *v;
897:   PetscInt       i;

900:   *ncols = A->cmap->n;
901:   if (cols) {
902:     PetscMalloc1(A->cmap->n+1,cols);
903:     for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
904:   }
905:   if (vals) {
906:     PetscMalloc1(A->cmap->n+1,vals);
907:     v    = mat->v + row;
908:     for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
909:   }
910:   return(0);
911: }

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

918:   if (cols) {PetscFree(*cols);}
919:   if (vals) {PetscFree(*vals); }
920:   return(0);
921: }
922: /* ----------------------------------------------------------------*/
923: static PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
924: {
925:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
926:   PetscInt     i,j,idx=0;

929:   if (!mat->roworiented) {
930:     if (addv == INSERT_VALUES) {
931:       for (j=0; j<n; j++) {
932:         if (indexn[j] < 0) {idx += m; continue;}
933: #if defined(PETSC_USE_DEBUG)
934:         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);
935: #endif
936:         for (i=0; i<m; i++) {
937:           if (indexm[i] < 0) {idx++; continue;}
938: #if defined(PETSC_USE_DEBUG)
939:           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);
940: #endif
941:           mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
942:         }
943:       }
944:     } else {
945:       for (j=0; j<n; j++) {
946:         if (indexn[j] < 0) {idx += m; continue;}
947: #if defined(PETSC_USE_DEBUG)
948:         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);
949: #endif
950:         for (i=0; i<m; i++) {
951:           if (indexm[i] < 0) {idx++; continue;}
952: #if defined(PETSC_USE_DEBUG)
953:           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);
954: #endif
955:           mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
956:         }
957:       }
958:     }
959:   } else {
960:     if (addv == INSERT_VALUES) {
961:       for (i=0; i<m; i++) {
962:         if (indexm[i] < 0) { idx += n; continue;}
963: #if defined(PETSC_USE_DEBUG)
964:         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);
965: #endif
966:         for (j=0; j<n; j++) {
967:           if (indexn[j] < 0) { idx++; continue;}
968: #if defined(PETSC_USE_DEBUG)
969:           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);
970: #endif
971:           mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
972:         }
973:       }
974:     } else {
975:       for (i=0; i<m; i++) {
976:         if (indexm[i] < 0) { idx += n; continue;}
977: #if defined(PETSC_USE_DEBUG)
978:         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);
979: #endif
980:         for (j=0; j<n; j++) {
981:           if (indexn[j] < 0) { idx++; continue;}
982: #if defined(PETSC_USE_DEBUG)
983:           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);
984: #endif
985:           mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
986:         }
987:       }
988:     }
989:   }
990:   return(0);
991: }

993: static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
994: {
995:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
996:   PetscInt     i,j;

999:   /* row-oriented output */
1000:   for (i=0; i<m; i++) {
1001:     if (indexm[i] < 0) {v += n;continue;}
1002:     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);
1003:     for (j=0; j<n; j++) {
1004:       if (indexn[j] < 0) {v++; continue;}
1005:       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);
1006:       *v++ = mat->v[indexn[j]*mat->lda + indexm[i]];
1007:     }
1008:   }
1009:   return(0);
1010: }

1012: /* -----------------------------------------------------------------*/

1014: static PetscErrorCode MatLoad_SeqDense(Mat newmat,PetscViewer viewer)
1015: {
1016:   Mat_SeqDense   *a;
1018:   PetscInt       *scols,i,j,nz,header[4];
1019:   int            fd;
1020:   PetscMPIInt    size;
1021:   PetscInt       *rowlengths = 0,M,N,*cols,grows,gcols;
1022:   PetscScalar    *vals,*svals,*v,*w;
1023:   MPI_Comm       comm;

1026:   /* force binary viewer to load .info file if it has not yet done so */
1027:   PetscViewerSetUp(viewer);
1028:   PetscObjectGetComm((PetscObject)viewer,&comm);
1029:   MPI_Comm_size(comm,&size);
1030:   if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor");
1031:   PetscViewerBinaryGetDescriptor(viewer,&fd);
1032:   PetscBinaryRead(fd,header,4,PETSC_INT);
1033:   if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not matrix object");
1034:   M = header[1]; N = header[2]; nz = header[3];

1036:   /* set global size if not set already*/
1037:   if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) {
1038:     MatSetSizes(newmat,M,N,M,N);
1039:   } else {
1040:     /* if sizes and type are already set, check if the vector global sizes are correct */
1041:     MatGetSize(newmat,&grows,&gcols);
1042:     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);
1043:   }
1044:   a = (Mat_SeqDense*)newmat->data;
1045:   if (!a->user_alloc) {
1046:     MatSeqDenseSetPreallocation(newmat,NULL);
1047:   }

1049:   if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */
1050:     a = (Mat_SeqDense*)newmat->data;
1051:     v = a->v;
1052:     /* Allocate some temp space to read in the values and then flip them
1053:        from row major to column major */
1054:     PetscMalloc1(M*N > 0 ? M*N : 1,&w);
1055:     /* read in nonzero values */
1056:     PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);
1057:     /* now flip the values and store them in the matrix*/
1058:     for (j=0; j<N; j++) {
1059:       for (i=0; i<M; i++) {
1060:         *v++ =w[i*N+j];
1061:       }
1062:     }
1063:     PetscFree(w);
1064:   } else {
1065:     /* read row lengths */
1066:     PetscMalloc1(M+1,&rowlengths);
1067:     PetscBinaryRead(fd,rowlengths,M,PETSC_INT);

1069:     a = (Mat_SeqDense*)newmat->data;
1070:     v = a->v;

1072:     /* read column indices and nonzeros */
1073:     PetscMalloc1(nz+1,&scols);
1074:     cols = scols;
1075:     PetscBinaryRead(fd,cols,nz,PETSC_INT);
1076:     PetscMalloc1(nz+1,&svals);
1077:     vals = svals;
1078:     PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);

1080:     /* insert into matrix */
1081:     for (i=0; i<M; i++) {
1082:       for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j];
1083:       svals += rowlengths[i]; scols += rowlengths[i];
1084:     }
1085:     PetscFree(vals);
1086:     PetscFree(cols);
1087:     PetscFree(rowlengths);
1088:   }
1089:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1090:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
1091:   return(0);
1092: }

1094: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
1095: {
1096:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
1097:   PetscErrorCode    ierr;
1098:   PetscInt          i,j;
1099:   const char        *name;
1100:   PetscScalar       *v;
1101:   PetscViewerFormat format;
1102: #if defined(PETSC_USE_COMPLEX)
1103:   PetscBool         allreal = PETSC_TRUE;
1104: #endif

1107:   PetscViewerGetFormat(viewer,&format);
1108:   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1109:     return(0);  /* do nothing for now */
1110:   } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1111:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1112:     for (i=0; i<A->rmap->n; i++) {
1113:       v    = a->v + i;
1114:       PetscViewerASCIIPrintf(viewer,"row %D:",i);
1115:       for (j=0; j<A->cmap->n; j++) {
1116: #if defined(PETSC_USE_COMPLEX)
1117:         if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
1118:           PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1119:         } else if (PetscRealPart(*v)) {
1120:           PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));
1121:         }
1122: #else
1123:         if (*v) {
1124:           PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);
1125:         }
1126: #endif
1127:         v += a->lda;
1128:       }
1129:       PetscViewerASCIIPrintf(viewer,"\n");
1130:     }
1131:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1132:   } else {
1133:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1134: #if defined(PETSC_USE_COMPLEX)
1135:     /* determine if matrix has all real values */
1136:     v = a->v;
1137:     for (i=0; i<A->rmap->n*A->cmap->n; i++) {
1138:       if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;}
1139:     }
1140: #endif
1141:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
1142:       PetscObjectGetName((PetscObject)A,&name);
1143:       PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);
1144:       PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);
1145:       PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
1146:     }

1148:     for (i=0; i<A->rmap->n; i++) {
1149:       v = a->v + i;
1150:       for (j=0; j<A->cmap->n; j++) {
1151: #if defined(PETSC_USE_COMPLEX)
1152:         if (allreal) {
1153:           PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));
1154:         } else {
1155:           PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1156:         }
1157: #else
1158:         PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);
1159: #endif
1160:         v += a->lda;
1161:       }
1162:       PetscViewerASCIIPrintf(viewer,"\n");
1163:     }
1164:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
1165:       PetscViewerASCIIPrintf(viewer,"];\n");
1166:     }
1167:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1168:   }
1169:   PetscViewerFlush(viewer);
1170:   return(0);
1171: }

1173: static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer)
1174: {
1175:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
1176:   PetscErrorCode    ierr;
1177:   int               fd;
1178:   PetscInt          ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n;
1179:   PetscScalar       *v,*anonz,*vals;
1180:   PetscViewerFormat format;

1183:   PetscViewerBinaryGetDescriptor(viewer,&fd);

1185:   PetscViewerGetFormat(viewer,&format);
1186:   if (format == PETSC_VIEWER_NATIVE) {
1187:     /* store the matrix as a dense matrix */
1188:     PetscMalloc1(4,&col_lens);

1190:     col_lens[0] = MAT_FILE_CLASSID;
1191:     col_lens[1] = m;
1192:     col_lens[2] = n;
1193:     col_lens[3] = MATRIX_BINARY_FORMAT_DENSE;

1195:     PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);
1196:     PetscFree(col_lens);

1198:     /* write out matrix, by rows */
1199:     PetscMalloc1(m*n+1,&vals);
1200:     v    = a->v;
1201:     for (j=0; j<n; j++) {
1202:       for (i=0; i<m; i++) {
1203:         vals[j + i*n] = *v++;
1204:       }
1205:     }
1206:     PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);
1207:     PetscFree(vals);
1208:   } else {
1209:     PetscMalloc1(4+nz,&col_lens);

1211:     col_lens[0] = MAT_FILE_CLASSID;
1212:     col_lens[1] = m;
1213:     col_lens[2] = n;
1214:     col_lens[3] = nz;

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

1220:     /* Possibly should write in smaller increments, not whole matrix at once? */
1221:     /* store column indices (zero start index) */
1222:     ict = 0;
1223:     for (i=0; i<m; i++) {
1224:       for (j=0; j<n; j++) col_lens[ict++] = j;
1225:     }
1226:     PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);
1227:     PetscFree(col_lens);

1229:     /* store nonzero values */
1230:     PetscMalloc1(nz+1,&anonz);
1231:     ict  = 0;
1232:     for (i=0; i<m; i++) {
1233:       v = a->v + i;
1234:       for (j=0; j<n; j++) {
1235:         anonz[ict++] = *v; v += a->lda;
1236:       }
1237:     }
1238:     PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);
1239:     PetscFree(anonz);
1240:   }
1241:   return(0);
1242: }

1244:  #include <petscdraw.h>
1245: static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa)
1246: {
1247:   Mat               A  = (Mat) Aa;
1248:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
1249:   PetscErrorCode    ierr;
1250:   PetscInt          m  = A->rmap->n,n = A->cmap->n,i,j;
1251:   int               color = PETSC_DRAW_WHITE;
1252:   PetscScalar       *v = a->v;
1253:   PetscViewer       viewer;
1254:   PetscReal         xl,yl,xr,yr,x_l,x_r,y_l,y_r;
1255:   PetscViewerFormat format;

1258:   PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
1259:   PetscViewerGetFormat(viewer,&format);
1260:   PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);

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

1264:   if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1265:     PetscDrawCollectiveBegin(draw);
1266:     /* Blue for negative and Red for positive */
1267:     for (j = 0; j < n; j++) {
1268:       x_l = j; x_r = x_l + 1.0;
1269:       for (i = 0; i < m; i++) {
1270:         y_l = m - i - 1.0;
1271:         y_r = y_l + 1.0;
1272:         if (PetscRealPart(v[j*m+i]) >  0.) {
1273:           color = PETSC_DRAW_RED;
1274:         } else if (PetscRealPart(v[j*m+i]) <  0.) {
1275:           color = PETSC_DRAW_BLUE;
1276:         } else {
1277:           continue;
1278:         }
1279:         PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1280:       }
1281:     }
1282:     PetscDrawCollectiveEnd(draw);
1283:   } else {
1284:     /* use contour shading to indicate magnitude of values */
1285:     /* first determine max of all nonzero values */
1286:     PetscReal minv = 0.0, maxv = 0.0;
1287:     PetscDraw popup;

1289:     for (i=0; i < m*n; i++) {
1290:       if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1291:     }
1292:     if (minv >= maxv) maxv = minv + PETSC_SMALL;
1293:     PetscDrawGetPopup(draw,&popup);
1294:     PetscDrawScalePopup(popup,minv,maxv);

1296:     PetscDrawCollectiveBegin(draw);
1297:     for (j=0; j<n; j++) {
1298:       x_l = j;
1299:       x_r = x_l + 1.0;
1300:       for (i=0; i<m; i++) {
1301:         y_l = m - i - 1.0;
1302:         y_r = y_l + 1.0;
1303:         color = PetscDrawRealToColor(PetscAbsScalar(v[j*m+i]),minv,maxv);
1304:         PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1305:       }
1306:     }
1307:     PetscDrawCollectiveEnd(draw);
1308:   }
1309:   return(0);
1310: }

1312: static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
1313: {
1314:   PetscDraw      draw;
1315:   PetscBool      isnull;
1316:   PetscReal      xr,yr,xl,yl,h,w;

1320:   PetscViewerDrawGetDraw(viewer,0,&draw);
1321:   PetscDrawIsNull(draw,&isnull);
1322:   if (isnull) return(0);

1324:   xr   = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
1325:   xr  += w;          yr += h;        xl = -w;     yl = -h;
1326:   PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
1327:   PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
1328:   PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
1329:   PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
1330:   PetscDrawSave(draw);
1331:   return(0);
1332: }

1334: PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
1335: {
1337:   PetscBool      iascii,isbinary,isdraw;

1340:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1341:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1342:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);

1344:   if (iascii) {
1345:     MatView_SeqDense_ASCII(A,viewer);
1346:   } else if (isbinary) {
1347:     MatView_SeqDense_Binary(A,viewer);
1348:   } else if (isdraw) {
1349:     MatView_SeqDense_Draw(A,viewer);
1350:   }
1351:   return(0);
1352: }

1354: static PetscErrorCode MatDestroy_SeqDense(Mat mat)
1355: {
1356:   Mat_SeqDense   *l = (Mat_SeqDense*)mat->data;

1360: #if defined(PETSC_USE_LOG)
1361:   PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n);
1362: #endif
1363:   PetscFree(l->pivots);
1364:   PetscFree(l->fwork);
1365:   MatDestroy(&l->ptapwork);
1366:   if (!l->user_alloc) {PetscFree(l->v);}
1367:   PetscFree(mat->data);

1369:   PetscObjectChangeTypeName((PetscObject)mat,0);
1370:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);
1371:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);
1372:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);
1373: #if defined(PETSC_HAVE_ELEMENTAL)
1374:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);
1375: #endif
1376:   PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);
1377:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);
1378:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);
1379:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);
1380:   PetscObjectComposeFunction((PetscObject)mat,"MatPtAP_seqaij_seqdense_C",NULL);
1381:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);
1382:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);
1383:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);
1384:   return(0);
1385: }

1387: static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
1388: {
1389:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1391:   PetscInt       k,j,m,n,M;
1392:   PetscScalar    *v,tmp;

1395:   v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n;
1396:   if (reuse == MAT_INPLACE_MATRIX) { /* in place transpose */
1397:     if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can not transpose non-square matrix in place");
1398:     else {
1399:       for (j=0; j<m; j++) {
1400:         for (k=0; k<j; k++) {
1401:           tmp        = v[j + k*M];
1402:           v[j + k*M] = v[k + j*M];
1403:           v[k + j*M] = tmp;
1404:         }
1405:       }
1406:     }
1407:   } else { /* out-of-place transpose */
1408:     Mat          tmat;
1409:     Mat_SeqDense *tmatd;
1410:     PetscScalar  *v2;
1411:     PetscInt     M2;

1413:     if (reuse == MAT_INITIAL_MATRIX) {
1414:       MatCreate(PetscObjectComm((PetscObject)A),&tmat);
1415:       MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);
1416:       MatSetType(tmat,((PetscObject)A)->type_name);
1417:       MatSeqDenseSetPreallocation(tmat,NULL);
1418:     } else {
1419:       tmat = *matout;
1420:     }
1421:     tmatd = (Mat_SeqDense*)tmat->data;
1422:     v = mat->v; v2 = tmatd->v; M2 = tmatd->lda;
1423:     for (j=0; j<n; j++) {
1424:       for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M];
1425:     }
1426:     MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
1427:     MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
1428:     *matout = tmat;
1429:   }
1430:   return(0);
1431: }

1433: static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool  *flg)
1434: {
1435:   Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data;
1436:   Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data;
1437:   PetscInt     i,j;
1438:   PetscScalar  *v1,*v2;

1441:   if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; return(0);}
1442:   if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; return(0);}
1443:   for (i=0; i<A1->rmap->n; i++) {
1444:     v1 = mat1->v+i; v2 = mat2->v+i;
1445:     for (j=0; j<A1->cmap->n; j++) {
1446:       if (*v1 != *v2) {*flg = PETSC_FALSE; return(0);}
1447:       v1 += mat1->lda; v2 += mat2->lda;
1448:     }
1449:   }
1450:   *flg = PETSC_TRUE;
1451:   return(0);
1452: }

1454: static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
1455: {
1456:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1458:   PetscInt       i,n,len;
1459:   PetscScalar    *x,zero = 0.0;

1462:   VecSet(v,zero);
1463:   VecGetSize(v,&n);
1464:   VecGetArray(v,&x);
1465:   len  = PetscMin(A->rmap->n,A->cmap->n);
1466:   if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
1467:   for (i=0; i<len; i++) {
1468:     x[i] = mat->v[i*mat->lda + i];
1469:   }
1470:   VecRestoreArray(v,&x);
1471:   return(0);
1472: }

1474: static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
1475: {
1476:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1477:   const PetscScalar *l,*r;
1478:   PetscScalar       x,*v;
1479:   PetscErrorCode    ierr;
1480:   PetscInt          i,j,m = A->rmap->n,n = A->cmap->n;

1483:   if (ll) {
1484:     VecGetSize(ll,&m);
1485:     VecGetArrayRead(ll,&l);
1486:     if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size");
1487:     for (i=0; i<m; i++) {
1488:       x = l[i];
1489:       v = mat->v + i;
1490:       for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;}
1491:     }
1492:     VecRestoreArrayRead(ll,&l);
1493:     PetscLogFlops(1.0*n*m);
1494:   }
1495:   if (rr) {
1496:     VecGetSize(rr,&n);
1497:     VecGetArrayRead(rr,&r);
1498:     if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size");
1499:     for (i=0; i<n; i++) {
1500:       x = r[i];
1501:       v = mat->v + i*mat->lda;
1502:       for (j=0; j<m; j++) (*v++) *= x;
1503:     }
1504:     VecRestoreArrayRead(rr,&r);
1505:     PetscLogFlops(1.0*n*m);
1506:   }
1507:   return(0);
1508: }

1510: static PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1511: {
1512:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1513:   PetscScalar    *v   = mat->v;
1514:   PetscReal      sum  = 0.0;
1515:   PetscInt       lda  =mat->lda,m=A->rmap->n,i,j;

1519:   if (type == NORM_FROBENIUS) {
1520:     if (lda>m) {
1521:       for (j=0; j<A->cmap->n; j++) {
1522:         v = mat->v+j*lda;
1523:         for (i=0; i<m; i++) {
1524:           sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1525:         }
1526:       }
1527:     } else {
1528: #if defined(PETSC_USE_REAL___FP16)
1529:       PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n;
1530:       *nrm = BLASnrm2_(&cnt,v,&one);
1531:     }
1532: #else
1533:       for (i=0; i<A->cmap->n*A->rmap->n; i++) {
1534:         sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1535:       }
1536:     }
1537:     *nrm = PetscSqrtReal(sum);
1538: #endif
1539:     PetscLogFlops(2.0*A->cmap->n*A->rmap->n);
1540:   } else if (type == NORM_1) {
1541:     *nrm = 0.0;
1542:     for (j=0; j<A->cmap->n; j++) {
1543:       v   = mat->v + j*mat->lda;
1544:       sum = 0.0;
1545:       for (i=0; i<A->rmap->n; i++) {
1546:         sum += PetscAbsScalar(*v);  v++;
1547:       }
1548:       if (sum > *nrm) *nrm = sum;
1549:     }
1550:     PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1551:   } else if (type == NORM_INFINITY) {
1552:     *nrm = 0.0;
1553:     for (j=0; j<A->rmap->n; j++) {
1554:       v   = mat->v + j;
1555:       sum = 0.0;
1556:       for (i=0; i<A->cmap->n; i++) {
1557:         sum += PetscAbsScalar(*v); v += mat->lda;
1558:       }
1559:       if (sum > *nrm) *nrm = sum;
1560:     }
1561:     PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1562:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
1563:   return(0);
1564: }

1566: static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
1567: {
1568:   Mat_SeqDense   *aij = (Mat_SeqDense*)A->data;

1572:   switch (op) {
1573:   case MAT_ROW_ORIENTED:
1574:     aij->roworiented = flg;
1575:     break;
1576:   case MAT_NEW_NONZERO_LOCATIONS:
1577:   case MAT_NEW_NONZERO_LOCATION_ERR:
1578:   case MAT_NEW_NONZERO_ALLOCATION_ERR:
1579:   case MAT_NEW_DIAGONALS:
1580:   case MAT_KEEP_NONZERO_PATTERN:
1581:   case MAT_IGNORE_OFF_PROC_ENTRIES:
1582:   case MAT_USE_HASH_TABLE:
1583:   case MAT_IGNORE_LOWER_TRIANGULAR:
1584:     PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
1585:     break;
1586:   case MAT_SPD:
1587:   case MAT_SYMMETRIC:
1588:   case MAT_STRUCTURALLY_SYMMETRIC:
1589:   case MAT_HERMITIAN:
1590:   case MAT_SYMMETRY_ETERNAL:
1591:     /* These options are handled directly by MatSetOption() */
1592:     break;
1593:   default:
1594:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]);
1595:   }
1596:   return(0);
1597: }

1599: static PetscErrorCode MatZeroEntries_SeqDense(Mat A)
1600: {
1601:   Mat_SeqDense   *l = (Mat_SeqDense*)A->data;
1603:   PetscInt       lda=l->lda,m=A->rmap->n,j;

1606:   if (lda>m) {
1607:     for (j=0; j<A->cmap->n; j++) {
1608:       PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));
1609:     }
1610:   } else {
1611:     PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1612:   }
1613:   return(0);
1614: }

1616: static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1617: {
1618:   PetscErrorCode    ierr;
1619:   Mat_SeqDense      *l = (Mat_SeqDense*)A->data;
1620:   PetscInt          m  = l->lda, n = A->cmap->n, i,j;
1621:   PetscScalar       *slot,*bb;
1622:   const PetscScalar *xx;

1625: #if defined(PETSC_USE_DEBUG)
1626:   for (i=0; i<N; i++) {
1627:     if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
1628:     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);
1629:   }
1630: #endif

1632:   /* fix right hand side if needed */
1633:   if (x && b) {
1634:     VecGetArrayRead(x,&xx);
1635:     VecGetArray(b,&bb);
1636:     for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
1637:     VecRestoreArrayRead(x,&xx);
1638:     VecRestoreArray(b,&bb);
1639:   }

1641:   for (i=0; i<N; i++) {
1642:     slot = l->v + rows[i];
1643:     for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
1644:   }
1645:   if (diag != 0.0) {
1646:     if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
1647:     for (i=0; i<N; i++) {
1648:       slot  = l->v + (m+1)*rows[i];
1649:       *slot = diag;
1650:     }
1651:   }
1652:   return(0);
1653: }

1655: static PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[])
1656: {
1657:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;

1660:   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");
1661:   *array = mat->v;
1662:   return(0);
1663: }

1665: static PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[])
1666: {
1668:   *array = 0; /* user cannot accidently use the array later */
1669:   return(0);
1670: }

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

1675:    Not Collective

1677:    Input Parameter:
1678: .  mat - a MATSEQDENSE or MATMPIDENSE matrix

1680:    Output Parameter:
1681: .   array - pointer to the data

1683:    Level: intermediate

1685: .seealso: MatDenseRestoreArray()
1686: @*/
1687: PetscErrorCode  MatDenseGetArray(Mat A,PetscScalar **array)
1688: {

1692:   PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));
1693:   return(0);
1694: }

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

1699:    Not Collective

1701:    Input Parameters:
1702: .  mat - a MATSEQDENSE or MATMPIDENSE matrix
1703: .  array - pointer to the data

1705:    Level: intermediate

1707: .seealso: MatDenseGetArray()
1708: @*/
1709: PetscErrorCode  MatDenseRestoreArray(Mat A,PetscScalar **array)
1710: {

1714:   PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));
1715:   return(0);
1716: }

1718: static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B)
1719: {
1720:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1722:   PetscInt       i,j,nrows,ncols;
1723:   const PetscInt *irow,*icol;
1724:   PetscScalar    *av,*bv,*v = mat->v;
1725:   Mat            newmat;

1728:   ISGetIndices(isrow,&irow);
1729:   ISGetIndices(iscol,&icol);
1730:   ISGetLocalSize(isrow,&nrows);
1731:   ISGetLocalSize(iscol,&ncols);

1733:   /* Check submatrixcall */
1734:   if (scall == MAT_REUSE_MATRIX) {
1735:     PetscInt n_cols,n_rows;
1736:     MatGetSize(*B,&n_rows,&n_cols);
1737:     if (n_rows != nrows || n_cols != ncols) {
1738:       /* resize the result matrix to match number of requested rows/columns */
1739:       MatSetSizes(*B,nrows,ncols,nrows,ncols);
1740:     }
1741:     newmat = *B;
1742:   } else {
1743:     /* Create and fill new matrix */
1744:     MatCreate(PetscObjectComm((PetscObject)A),&newmat);
1745:     MatSetSizes(newmat,nrows,ncols,nrows,ncols);
1746:     MatSetType(newmat,((PetscObject)A)->type_name);
1747:     MatSeqDenseSetPreallocation(newmat,NULL);
1748:   }

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

1753:   for (i=0; i<ncols; i++) {
1754:     av = v + mat->lda*icol[i];
1755:     for (j=0; j<nrows; j++) *bv++ = av[irow[j]];
1756:   }

1758:   /* Assemble the matrices so that the correct flags are set */
1759:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1760:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);

1762:   /* Free work space */
1763:   ISRestoreIndices(isrow,&irow);
1764:   ISRestoreIndices(iscol,&icol);
1765:   *B   = newmat;
1766:   return(0);
1767: }

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

1775:   if (scall == MAT_INITIAL_MATRIX) {
1776:     PetscCalloc1(n+1,B);
1777:   }

1779:   for (i=0; i<n; i++) {
1780:     MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
1781:   }
1782:   return(0);
1783: }

1785: static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
1786: {
1788:   return(0);
1789: }

1791: static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
1792: {
1794:   return(0);
1795: }

1797: static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
1798: {
1799:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data;
1801:   PetscInt       lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j;

1804:   /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
1805:   if (A->ops->copy != B->ops->copy) {
1806:     MatCopy_Basic(A,B,str);
1807:     return(0);
1808:   }
1809:   if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)");
1810:   if (lda1>m || lda2>m) {
1811:     for (j=0; j<n; j++) {
1812:       PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));
1813:     }
1814:   } else {
1815:     PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1816:   }
1817:   return(0);
1818: }

1820: static PetscErrorCode MatSetUp_SeqDense(Mat A)
1821: {

1825:    MatSeqDenseSetPreallocation(A,0);
1826:   return(0);
1827: }

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

1836:   for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
1837:   return(0);
1838: }

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

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

1851: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
1852: {
1853:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1854:   PetscInt     i,nz = A->rmap->n*A->cmap->n;
1855:   PetscScalar  *aa = a->v;

1858:   for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
1859:   return(0);
1860: }

1862: /* ----------------------------------------------------------------*/
1863: PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1864: {

1868:   if (scall == MAT_INITIAL_MATRIX) {
1869:     PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
1870:     MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
1871:     PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
1872:   }
1873:   PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
1874:   MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);
1875:   PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
1876:   return(0);
1877: }

1879: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1880: {
1882:   PetscInt       m=A->rmap->n,n=B->cmap->n;
1883:   Mat            Cmat;

1886:   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);
1887:   MatCreate(PETSC_COMM_SELF,&Cmat);
1888:   MatSetSizes(Cmat,m,n,m,n);
1889:   MatSetType(Cmat,MATSEQDENSE);
1890:   MatSeqDenseSetPreallocation(Cmat,NULL);

1892:   *C = Cmat;
1893:   return(0);
1894: }

1896: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
1897: {
1898:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
1899:   Mat_SeqDense   *b = (Mat_SeqDense*)B->data;
1900:   Mat_SeqDense   *c = (Mat_SeqDense*)C->data;
1901:   PetscBLASInt   m,n,k;
1902:   PetscScalar    _DOne=1.0,_DZero=0.0;
1904:   PetscBool      flg;

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

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

1918:   PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&flg);
1919:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"First matrix must be dense");
1920:   PetscBLASIntCast(C->rmap->n,&m);
1921:   PetscBLASIntCast(C->cmap->n,&n);
1922:   PetscBLASIntCast(A->cmap->n,&k);
1923:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
1924:   return(0);
1925: }

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

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

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

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

1956:   Cmat->assembled = PETSC_TRUE;

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

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

1972:   PetscBLASIntCast(C->rmap->n,&m);
1973:   PetscBLASIntCast(C->cmap->n,&n);
1974:   PetscBLASIntCast(A->rmap->n,&k);
1975:   PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
1976:   return(0);
1977: }

1979: static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[])
1980: {
1981:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
1983:   PetscInt       i,j,m = A->rmap->n,n = A->cmap->n,p;
1984:   PetscScalar    *x;
1985:   MatScalar      *aa = a->v;

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

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

2004: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])
2005: {
2006:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2008:   PetscInt       i,j,m = A->rmap->n,n = A->cmap->n,p;
2009:   PetscScalar    *x;
2010:   PetscReal      atmp;
2011:   MatScalar      *aa = a->v;

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

2016:   VecSet(v,0.0);
2017:   VecGetArray(v,&x);
2018:   VecGetLocalSize(v,&p);
2019:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2020:   for (i=0; i<m; i++) {
2021:     x[i] = PetscAbsScalar(aa[i]);
2022:     for (j=1; j<n; j++) {
2023:       atmp = PetscAbsScalar(aa[i+m*j]);
2024:       if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;}
2025:     }
2026:   }
2027:   VecRestoreArray(v,&x);
2028:   return(0);
2029: }

2031: static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])
2032: {
2033:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2035:   PetscInt       i,j,m = A->rmap->n,n = A->cmap->n,p;
2036:   PetscScalar    *x;
2037:   MatScalar      *aa = a->v;

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

2042:   VecSet(v,0.0);
2043:   VecGetArray(v,&x);
2044:   VecGetLocalSize(v,&p);
2045:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2046:   for (i=0; i<m; i++) {
2047:     x[i] = aa[i]; if (idx) idx[i] = 0;
2048:     for (j=1; j<n; j++) {
2049:       if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
2050:     }
2051:   }
2052:   VecRestoreArray(v,&x);
2053:   return(0);
2054: }

2056: static PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
2057: {
2058:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2060:   PetscScalar    *x;

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

2065:   VecGetArray(v,&x);
2066:   PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));
2067:   VecRestoreArray(v,&x);
2068:   return(0);
2069: }


2072: PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms)
2073: {
2075:   PetscInt       i,j,m,n;
2076:   PetscScalar    *a;

2079:   MatGetSize(A,&m,&n);
2080:   PetscMemzero(norms,n*sizeof(PetscReal));
2081:   MatDenseGetArray(A,&a);
2082:   if (type == NORM_2) {
2083:     for (i=0; i<n; i++) {
2084:       for (j=0; j<m; j++) {
2085:         norms[i] += PetscAbsScalar(a[j]*a[j]);
2086:       }
2087:       a += m;
2088:     }
2089:   } else if (type == NORM_1) {
2090:     for (i=0; i<n; i++) {
2091:       for (j=0; j<m; j++) {
2092:         norms[i] += PetscAbsScalar(a[j]);
2093:       }
2094:       a += m;
2095:     }
2096:   } else if (type == NORM_INFINITY) {
2097:     for (i=0; i<n; i++) {
2098:       for (j=0; j<m; j++) {
2099:         norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]);
2100:       }
2101:       a += m;
2102:     }
2103:   } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType");
2104:   MatDenseRestoreArray(A,&a);
2105:   if (type == NORM_2) {
2106:     for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]);
2107:   }
2108:   return(0);
2109: }

2111: static PetscErrorCode  MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
2112: {
2114:   PetscScalar    *a;
2115:   PetscInt       m,n,i;

2118:   MatGetSize(x,&m,&n);
2119:   MatDenseGetArray(x,&a);
2120:   for (i=0; i<m*n; i++) {
2121:     PetscRandomGetValue(rctx,a+i);
2122:   }
2123:   MatDenseRestoreArray(x,&a);
2124:   return(0);
2125: }

2127: static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool  *missing,PetscInt *d)
2128: {
2130:   *missing = PETSC_FALSE;
2131:   return(0);
2132: }


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

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

2285:    Collective on MPI_Comm

2287:    Input Parameters:
2288: +  comm - MPI communicator, set to PETSC_COMM_SELF
2289: .  m - number of rows
2290: .  n - number of columns
2291: -  data - optional location of matrix data in column major order.  Set data=NULL for PETSc
2292:    to control all matrix memory allocation.

2294:    Output Parameter:
2295: .  A - the matrix

2297:    Notes:
2298:    The data input variable is intended primarily for Fortran programmers
2299:    who wish to allocate their own matrix memory space.  Most users should
2300:    set data=NULL.

2302:    Level: intermediate

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

2306: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
2307: @*/
2308: PetscErrorCode  MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
2309: {

2313:   MatCreate(comm,A);
2314:   MatSetSizes(*A,m,n,m,n);
2315:   MatSetType(*A,MATSEQDENSE);
2316:   MatSeqDenseSetPreallocation(*A,data);
2317:   return(0);
2318: }

2320: /*@C
2321:    MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements

2323:    Collective on MPI_Comm

2325:    Input Parameters:
2326: +  B - the matrix
2327: -  data - the array (or NULL)

2329:    Notes:
2330:    The data input variable is intended primarily for Fortran programmers
2331:    who wish to allocate their own matrix memory space.  Most users should
2332:    need not call this routine.

2334:    Level: intermediate

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

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

2340: @*/
2341: PetscErrorCode  MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
2342: {

2346:   PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
2347:   return(0);
2348: }

2350: PetscErrorCode  MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
2351: {
2352:   Mat_SeqDense   *b;

2356:   B->preallocated = PETSC_TRUE;

2358:   PetscLayoutSetUp(B->rmap);
2359:   PetscLayoutSetUp(B->cmap);

2361:   b       = (Mat_SeqDense*)B->data;
2362:   b->Mmax = B->rmap->n;
2363:   b->Nmax = B->cmap->n;
2364:   if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n;

2366:   PetscIntMultError(b->lda,b->Nmax,NULL);
2367:   if (!data) { /* petsc-allocated storage */
2368:     if (!b->user_alloc) { PetscFree(b->v); }
2369:     PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);
2370:     PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));

2372:     b->user_alloc = PETSC_FALSE;
2373:   } else { /* user-allocated storage */
2374:     if (!b->user_alloc) { PetscFree(b->v); }
2375:     b->v          = data;
2376:     b->user_alloc = PETSC_TRUE;
2377:   }
2378:   B->assembled = PETSC_TRUE;
2379:   return(0);
2380: }

2382: #if defined(PETSC_HAVE_ELEMENTAL)
2383: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
2384: {
2385:   Mat            mat_elemental;
2387:   PetscScalar    *array,*v_colwise;
2388:   PetscInt       M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols;

2391:   PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);
2392:   MatDenseGetArray(A,&array);
2393:   /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
2394:   k = 0;
2395:   for (j=0; j<N; j++) {
2396:     cols[j] = j;
2397:     for (i=0; i<M; i++) {
2398:       v_colwise[j*M+i] = array[k++];
2399:     }
2400:   }
2401:   for (i=0; i<M; i++) {
2402:     rows[i] = i;
2403:   }
2404:   MatDenseRestoreArray(A,&array);

2406:   MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
2407:   MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);
2408:   MatSetType(mat_elemental,MATELEMENTAL);
2409:   MatSetUp(mat_elemental);

2411:   /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
2412:   MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);
2413:   MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
2414:   MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);
2415:   PetscFree3(v_colwise,rows,cols);

2417:   if (reuse == MAT_INPLACE_MATRIX) {
2418:     MatHeaderReplace(A,&mat_elemental);
2419:   } else {
2420:     *newmat = mat_elemental;
2421:   }
2422:   return(0);
2423: }
2424: #endif

2426: /*@C
2427:   MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array

2429:   Input parameter:
2430: + A - the matrix
2431: - lda - the leading dimension

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

2438:   Level: intermediate

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

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

2444: @*/
2445: PetscErrorCode  MatSeqDenseSetLDA(Mat B,PetscInt lda)
2446: {
2447:   Mat_SeqDense *b = (Mat_SeqDense*)B->data;

2450:   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);
2451:   b->lda       = lda;
2452:   b->changelda = PETSC_FALSE;
2453:   b->Mmax      = PetscMax(b->Mmax,lda);
2454:   return(0);
2455: }

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

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

2463:   Level: beginner

2465: .seealso: MatCreateSeqDense()

2467: M*/

2469: PETSC_EXTERN PetscErrorCode MatCreate_SeqDense(Mat B)
2470: {
2471:   Mat_SeqDense   *b;
2473:   PetscMPIInt    size;

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

2479:   PetscNewLog(B,&b);
2480:   PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
2481:   B->data = (void*)b;

2483:   b->roworiented = PETSC_TRUE;

2485:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);
2486:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);
2487:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);
2488: #if defined(PETSC_HAVE_ELEMENTAL)
2489:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);
2490: #endif
2491:   PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);
2492:   PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2493:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2494:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2495:   PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaij_seqdense_C",MatPtAP_SeqDense_SeqDense);
2496:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2497:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2498:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2499:   PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
2500:   return(0);
2501: }