Actual source code: dense.c

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

 20: #if defined(PETSC_USE_DEBUG)
 21:   for (i=0; i<N; i++) {
 22:     if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
 23:     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);
 24:     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);
 25:   }
 26: #endif

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

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

 55: PetscErrorCode MatPtAPNumeric_SeqDense_SeqDense(Mat A,Mat P,Mat C)
 56: {
 57:   Mat_SeqDense   *c = (Mat_SeqDense*)(C->data);

 61:   MatMatMultNumeric_SeqDense_SeqDense(A,P,c->ptapwork);
 62:   MatTransposeMatMultNumeric_SeqDense_SeqDense(P,c->ptapwork,C);
 63:   return(0);
 64: }

 66: PetscErrorCode MatPtAPSymbolic_SeqDense_SeqDense(Mat A,Mat P,PetscReal fill,Mat *C)
 67: {
 68:   Mat_SeqDense   *c;

 72:   MatCreateSeqDense(PetscObjectComm((PetscObject)A),P->cmap->N,P->cmap->N,NULL,C);
 73:   c = (Mat_SeqDense*)((*C)->data);
 74:   MatCreateSeqDense(PetscObjectComm((PetscObject)A),A->rmap->N,P->cmap->N,NULL,&c->ptapwork);
 75:   return(0);
 76: }

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

 83:   if (reuse == MAT_INITIAL_MATRIX) {
 84:     MatPtAPSymbolic_SeqDense_SeqDense(A,P,fill,C);
 85:   }
 86:   PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);
 87:   (*(*C)->ops->ptapnumeric)(A,P,*C);
 88:   PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);
 89:   return(0);
 90: }

 92: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
 93: {
 94:   Mat            B = NULL;
 95:   Mat_SeqAIJ     *a = (Mat_SeqAIJ*)A->data;
 96:   Mat_SeqDense   *b;
 98:   PetscInt       *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->N,i;
 99:   MatScalar      *av=a->a;
100:   PetscBool      isseqdense;

103:   if (reuse == MAT_REUSE_MATRIX) {
104:     PetscObjectTypeCompare((PetscObject)*newmat,MATSEQDENSE,&isseqdense);
105:     if (!isseqdense) SETERRQ1(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_USER,"Cannot reuse matrix of type %s",((PetscObject)(*newmat))->type);
106:   }
107:   if (reuse != MAT_REUSE_MATRIX) {
108:     MatCreate(PetscObjectComm((PetscObject)A),&B);
109:     MatSetSizes(B,m,n,m,n);
110:     MatSetType(B,MATSEQDENSE);
111:     MatSeqDenseSetPreallocation(B,NULL);
112:     b    = (Mat_SeqDense*)(B->data);
113:   } else {
114:     b    = (Mat_SeqDense*)((*newmat)->data);
115:     PetscMemzero(b->v,m*n*sizeof(PetscScalar));
116:   }
117:   for (i=0; i<m; i++) {
118:     PetscInt j;
119:     for (j=0;j<ai[1]-ai[0];j++) {
120:       b->v[*aj*m+i] = *av;
121:       aj++;
122:       av++;
123:     }
124:     ai++;
125:   }

127:   if (reuse == MAT_INPLACE_MATRIX) {
128:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
129:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
130:     MatHeaderReplace(A,&B);
131:   } else {
132:     if (B) *newmat = B;
133:     MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
134:     MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
135:   }
136:   return(0);
137: }

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

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

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

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

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

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

207:   info->block_size        = 1.0;
208:   info->nz_allocated      = (double)N;
209:   info->nz_used           = (double)N;
210:   info->nz_unneeded       = (double)0;
211:   info->assemblies        = (double)A->num_ass;
212:   info->mallocs           = 0;
213:   info->memory            = ((PetscObject)A)->mem;
214:   info->fill_ratio_given  = 0;
215:   info->fill_ratio_needed = 0;
216:   info->factor_mallocs    = 0;
217:   return(0);
218: }

220: static PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha)
221: {
222:   Mat_SeqDense   *a     = (Mat_SeqDense*)A->data;
223:   PetscScalar    oalpha = alpha;
225:   PetscBLASInt   one = 1,j,nz,lda;

228:   PetscBLASIntCast(a->lda,&lda);
229:   if (lda>A->rmap->n) {
230:     PetscBLASIntCast(A->rmap->n,&nz);
231:     for (j=0; j<A->cmap->n; j++) {
232:       PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v+j*lda,&one));
233:     }
234:   } else {
235:     PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);
236:     PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v,&one));
237:   }
238:   PetscLogFlops(nz);
239:   return(0);
240: }

242: static PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool  *fl)
243: {
244:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
245:   PetscInt     i,j,m = A->rmap->n,N;
246:   PetscScalar  *v = a->v;

249:   *fl = PETSC_FALSE;
250:   if (A->rmap->n != A->cmap->n) return(0);
251:   N = a->lda;

253:   for (i=0; i<m; i++) {
254:     for (j=i+1; j<m; j++) {
255:       if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) return(0);
256:     }
257:   }
258:   *fl = PETSC_TRUE;
259:   return(0);
260: }

262: static PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues)
263: {
264:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data,*l;
266:   PetscInt       lda = (PetscInt)mat->lda,j,m;

269:   PetscLayoutReference(A->rmap,&newi->rmap);
270:   PetscLayoutReference(A->cmap,&newi->cmap);
271:   MatSeqDenseSetPreallocation(newi,NULL);
272:   if (cpvalues == MAT_COPY_VALUES) {
273:     l = (Mat_SeqDense*)newi->data;
274:     if (lda>A->rmap->n) {
275:       m = A->rmap->n;
276:       for (j=0; j<A->cmap->n; j++) {
277:         PetscMemcpy(l->v+j*m,mat->v+j*lda,m*sizeof(PetscScalar));
278:       }
279:     } else {
280:       PetscMemcpy(l->v,mat->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
281:     }
282:   }
283:   newi->assembled = PETSC_TRUE;
284:   return(0);
285: }

287: static PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat)
288: {

292:   MatCreate(PetscObjectComm((PetscObject)A),newmat);
293:   MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
294:   MatSetType(*newmat,((PetscObject)A)->type_name);
295:   MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);
296:   return(0);
297: }


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

302: static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
303: {
304:   MatFactorInfo  info;

308:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
309:   MatLUFactor_SeqDense(fact,0,0,&info);
310:   return(0);
311: }

313: static PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy)
314: {
315:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
316:   PetscErrorCode    ierr;
317:   const PetscScalar *x;
318:   PetscScalar       *y;
319:   PetscBLASInt      one = 1,info,m;

322:   PetscBLASIntCast(A->rmap->n,&m);
323:   VecGetArrayRead(xx,&x);
324:   VecGetArray(yy,&y);
325:   PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
326:   if (A->factortype == MAT_FACTOR_LU) {
327: #if defined(PETSC_MISSING_LAPACK_GETRS)
328:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
329: #else
330:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
331:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
332: #endif
333:   } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
334: #if defined(PETSC_MISSING_LAPACK_POTRS)
335:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
336: #else
337:     PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
338:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
339: #endif
340:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
341:   VecRestoreArrayRead(xx,&x);
342:   VecRestoreArray(yy,&y);
343:   PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
344:   return(0);
345: }

347: static PetscErrorCode MatMatSolve_SeqDense(Mat A,Mat B,Mat X)
348: {
349:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
351:   PetscScalar    *b,*x;
352:   PetscInt       n;
353:   PetscBLASInt   nrhs,info,m;
354:   PetscBool      flg;

357:   PetscBLASIntCast(A->rmap->n,&m);
358:   PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
359:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
360:   PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
361:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");

363:   MatGetSize(B,NULL,&n);
364:   PetscBLASIntCast(n,&nrhs);
365:   MatDenseGetArray(B,&b);
366:   MatDenseGetArray(X,&x);

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

370:   if (A->factortype == MAT_FACTOR_LU) {
371: #if defined(PETSC_MISSING_LAPACK_GETRS)
372:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
373: #else
374:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
375:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
376: #endif
377:   } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
378: #if defined(PETSC_MISSING_LAPACK_POTRS)
379:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
380: #else
381:     PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info));
382:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
383: #endif
384:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");

386:   MatDenseRestoreArray(B,&b);
387:   MatDenseRestoreArray(X,&x);
388:   PetscLogFlops(nrhs*(2.0*m*m - m));
389:   return(0);
390: }

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

401:   PetscBLASIntCast(A->rmap->n,&m);
402:   VecGetArrayRead(xx,&x);
403:   VecGetArray(yy,&y);
404:   PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
405:   /* assume if pivots exist then use LU; else Cholesky */
406:   if (mat->pivots) {
407: #if defined(PETSC_MISSING_LAPACK_GETRS)
408:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
409: #else
410:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
411:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve");
412: #endif
413:   } else {
414: #if defined(PETSC_MISSING_LAPACK_POTRS)
415:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
416: #else
417:     PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
418:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve");
419: #endif
420:   }
421:   VecRestoreArrayRead(xx,&x);
422:   VecRestoreArray(yy,&y);
423:   PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
424:   return(0);
425: }

427: static PetscErrorCode MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
428: {
429:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
430:   PetscErrorCode    ierr;
431:   const PetscScalar *x;
432:   PetscScalar       *y,sone = 1.0;
433:   Vec               tmp = 0;
434:   PetscBLASInt      one = 1,info,m;

437:   PetscBLASIntCast(A->rmap->n,&m);
438:   VecGetArrayRead(xx,&x);
439:   VecGetArray(yy,&y);
440:   if (!A->rmap->n || !A->cmap->n) return(0);
441:   if (yy == zz) {
442:     VecDuplicate(yy,&tmp);
443:     PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);
444:     VecCopy(yy,tmp);
445:   }
446:   PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
447:   /* assume if pivots exist then use LU; else Cholesky */
448:   if (mat->pivots) {
449: #if defined(PETSC_MISSING_LAPACK_GETRS)
450:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
451: #else
452:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
453:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
454: #endif
455:   } else {
456: #if defined(PETSC_MISSING_LAPACK_POTRS)
457:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
458: #else
459:     PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
460:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
461: #endif
462:   }
463:   if (tmp) {
464:     VecAXPY(yy,sone,tmp);
465:     VecDestroy(&tmp);
466:   } else {
467:     VecAXPY(yy,sone,zz);
468:   }
469:   VecRestoreArrayRead(xx,&x);
470:   VecRestoreArray(yy,&y);
471:   PetscLogFlops(2.0*A->cmap->n*A->cmap->n);
472:   return(0);
473: }

475: static PetscErrorCode MatSolveTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
476: {
477:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
478:   PetscErrorCode    ierr;
479:   const PetscScalar *x;
480:   PetscScalar       *y,sone = 1.0;
481:   Vec               tmp = NULL;
482:   PetscBLASInt      one = 1,info,m;

485:   PetscBLASIntCast(A->rmap->n,&m);
486:   if (!A->rmap->n || !A->cmap->n) return(0);
487:   VecGetArrayRead(xx,&x);
488:   VecGetArray(yy,&y);
489:   if (yy == zz) {
490:     VecDuplicate(yy,&tmp);
491:     PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);
492:     VecCopy(yy,tmp);
493:   }
494:   PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
495:   /* assume if pivots exist then use LU; else Cholesky */
496:   if (mat->pivots) {
497: #if defined(PETSC_MISSING_LAPACK_GETRS)
498:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
499: #else
500:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
501:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
502: #endif
503:   } else {
504: #if defined(PETSC_MISSING_LAPACK_POTRS)
505:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
506: #else
507:     PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
508:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
509: #endif
510:   }
511:   if (tmp) {
512:     VecAXPY(yy,sone,tmp);
513:     VecDestroy(&tmp);
514:   } else {
515:     VecAXPY(yy,sone,zz);
516:   }
517:   VecRestoreArrayRead(xx,&x);
518:   VecRestoreArray(yy,&y);
519:   PetscLogFlops(2.0*A->cmap->n*A->cmap->n);
520:   return(0);
521: }

523: /* ---------------------------------------------------------------*/
524: /* COMMENT: I have chosen to hide row permutation in the pivots,
525:    rather than put it in the Mat->row slot.*/
526: static PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo)
527: {
528: #if defined(PETSC_MISSING_LAPACK_GETRF)
530:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable.");
531: #else
532:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
534:   PetscBLASInt   n,m,info;

537:   PetscBLASIntCast(A->cmap->n,&n);
538:   PetscBLASIntCast(A->rmap->n,&m);
539:   if (!mat->pivots) {
540:     PetscMalloc1(A->rmap->n+1,&mat->pivots);
541:     PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
542:   }
543:   if (!A->rmap->n || !A->cmap->n) return(0);
544:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
545:   PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info));
546:   PetscFPTrapPop();

548:   if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization");
549:   if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization");
550:   A->ops->solve             = MatSolve_SeqDense;
551:   A->ops->solvetranspose    = MatSolveTranspose_SeqDense;
552:   A->ops->solveadd          = MatSolveAdd_SeqDense;
553:   A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
554:   A->factortype             = MAT_FACTOR_LU;

556:   PetscFree(A->solvertype);
557:   PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);

559:   PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);
560: #endif
561:   return(0);
562: }

564: static PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo)
565: {
566: #if defined(PETSC_MISSING_LAPACK_POTRF)
568:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable.");
569: #else
570:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
572:   PetscBLASInt   info,n;

575:   PetscBLASIntCast(A->cmap->n,&n);
576:   PetscFree(mat->pivots);

578:   if (!A->rmap->n || !A->cmap->n) return(0);
579:   PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info));
580:   if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1);
581:   A->ops->solve             = MatSolve_SeqDense;
582:   A->ops->solvetranspose    = MatSolveTranspose_SeqDense;
583:   A->ops->solveadd          = MatSolveAdd_SeqDense;
584:   A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
585:   A->factortype             = MAT_FACTOR_CHOLESKY;

587:   PetscFree(A->solvertype);
588:   PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);

590:   PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
591: #endif
592:   return(0);
593: }


596: PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
597: {
599:   MatFactorInfo  info;

602:   info.fill = 1.0;

604:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
605:   MatCholeskyFactor_SeqDense(fact,0,&info);
606:   return(0);
607: }

609: static PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
610: {
612:   fact->assembled                  = PETSC_TRUE;
613:   fact->preallocated               = PETSC_TRUE;
614:   fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
615:   return(0);
616: }

618: static PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info)
619: {
621:   fact->preallocated         = PETSC_TRUE;
622:   fact->assembled            = PETSC_TRUE;
623:   fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense;
624:   return(0);
625: }

627: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact)
628: {

632:   MatCreate(PetscObjectComm((PetscObject)A),fact);
633:   MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
634:   MatSetType(*fact,((PetscObject)A)->type_name);
635:   if (ftype == MAT_FACTOR_LU) {
636:     (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
637:   } else {
638:     (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
639:   }
640:   (*fact)->factortype = ftype;

642:   PetscFree((*fact)->solvertype);
643:   PetscStrallocpy(MATSOLVERPETSC,&(*fact)->solvertype);
644:   return(0);
645: }

647: /* ------------------------------------------------------------------*/
648: static PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx)
649: {
650:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
651:   PetscScalar       *x,*v = mat->v,zero = 0.0,xt;
652:   const PetscScalar *b;
653:   PetscErrorCode    ierr;
654:   PetscInt          m = A->rmap->n,i;
655:   PetscBLASInt      o = 1,bm;

658:   if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */
659:   PetscBLASIntCast(m,&bm);
660:   if (flag & SOR_ZERO_INITIAL_GUESS) {
661:     /* this is a hack fix, should have another version without the second BLASdot */
662:     VecSet(xx,zero);
663:   }
664:   VecGetArray(xx,&x);
665:   VecGetArrayRead(bb,&b);
666:   its  = its*lits;
667:   if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
668:   while (its--) {
669:     if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
670:       for (i=0; i<m; i++) {
671:         PetscStackCallBLAS("BLASdot",xt   = b[i] - BLASdot_(&bm,v+i,&bm,x,&o));
672:         x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
673:       }
674:     }
675:     if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
676:       for (i=m-1; i>=0; i--) {
677:         PetscStackCallBLAS("BLASdot",xt   = b[i] - BLASdot_(&bm,v+i,&bm,x,&o));
678:         x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
679:       }
680:     }
681:   }
682:   VecRestoreArrayRead(bb,&b);
683:   VecRestoreArray(xx,&x);
684:   return(0);
685: }

687: /* -----------------------------------------------------------------*/
688: PETSC_INTERN PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)
689: {
690:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
691:   const PetscScalar *v   = mat->v,*x;
692:   PetscScalar       *y;
693:   PetscErrorCode    ierr;
694:   PetscBLASInt      m, n,_One=1;
695:   PetscScalar       _DOne=1.0,_DZero=0.0;

698:   PetscBLASIntCast(A->rmap->n,&m);
699:   PetscBLASIntCast(A->cmap->n,&n);
700:   if (!A->rmap->n || !A->cmap->n) return(0);
701:   VecGetArrayRead(xx,&x);
702:   VecGetArray(yy,&y);
703:   PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One));
704:   VecRestoreArrayRead(xx,&x);
705:   VecRestoreArray(yy,&y);
706:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);
707:   return(0);
708: }

710: PETSC_INTERN PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy)
711: {
712:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
713:   PetscScalar       *y,_DOne=1.0,_DZero=0.0;
714:   PetscErrorCode    ierr;
715:   PetscBLASInt      m, n, _One=1;
716:   const PetscScalar *v = mat->v,*x;

719:   PetscBLASIntCast(A->rmap->n,&m);
720:   PetscBLASIntCast(A->cmap->n,&n);
721:   if (!A->rmap->n || !A->cmap->n) return(0);
722:   VecGetArrayRead(xx,&x);
723:   VecGetArray(yy,&y);
724:   PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One));
725:   VecRestoreArrayRead(xx,&x);
726:   VecRestoreArray(yy,&y);
727:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);
728:   return(0);
729: }

731: PETSC_INTERN PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
732: {
733:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
734:   const PetscScalar *v = mat->v,*x;
735:   PetscScalar       *y,_DOne=1.0;
736:   PetscErrorCode    ierr;
737:   PetscBLASInt      m, n, _One=1;

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

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

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

776: /* -----------------------------------------------------------------*/
777: static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
778: {
779:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
780:   PetscScalar    *v;
782:   PetscInt       i;

785:   *ncols = A->cmap->n;
786:   if (cols) {
787:     PetscMalloc1(A->cmap->n+1,cols);
788:     for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
789:   }
790:   if (vals) {
791:     PetscMalloc1(A->cmap->n+1,vals);
792:     v    = mat->v + row;
793:     for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
794:   }
795:   return(0);
796: }

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

803:   if (cols) {PetscFree(*cols);}
804:   if (vals) {PetscFree(*vals); }
805:   return(0);
806: }
807: /* ----------------------------------------------------------------*/
808: static PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
809: {
810:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
811:   PetscInt     i,j,idx=0;

814:   if (!mat->roworiented) {
815:     if (addv == INSERT_VALUES) {
816:       for (j=0; j<n; j++) {
817:         if (indexn[j] < 0) {idx += m; continue;}
818: #if defined(PETSC_USE_DEBUG)
819:         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);
820: #endif
821:         for (i=0; i<m; i++) {
822:           if (indexm[i] < 0) {idx++; continue;}
823: #if defined(PETSC_USE_DEBUG)
824:           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);
825: #endif
826:           mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
827:         }
828:       }
829:     } else {
830:       for (j=0; j<n; j++) {
831:         if (indexn[j] < 0) {idx += m; continue;}
832: #if defined(PETSC_USE_DEBUG)
833:         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);
834: #endif
835:         for (i=0; i<m; i++) {
836:           if (indexm[i] < 0) {idx++; continue;}
837: #if defined(PETSC_USE_DEBUG)
838:           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);
839: #endif
840:           mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
841:         }
842:       }
843:     }
844:   } else {
845:     if (addv == INSERT_VALUES) {
846:       for (i=0; i<m; i++) {
847:         if (indexm[i] < 0) { idx += n; continue;}
848: #if defined(PETSC_USE_DEBUG)
849:         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);
850: #endif
851:         for (j=0; j<n; j++) {
852:           if (indexn[j] < 0) { idx++; continue;}
853: #if defined(PETSC_USE_DEBUG)
854:           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);
855: #endif
856:           mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
857:         }
858:       }
859:     } else {
860:       for (i=0; i<m; i++) {
861:         if (indexm[i] < 0) { idx += n; continue;}
862: #if defined(PETSC_USE_DEBUG)
863:         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);
864: #endif
865:         for (j=0; j<n; j++) {
866:           if (indexn[j] < 0) { idx++; continue;}
867: #if defined(PETSC_USE_DEBUG)
868:           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);
869: #endif
870:           mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
871:         }
872:       }
873:     }
874:   }
875:   return(0);
876: }

878: static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
879: {
880:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
881:   PetscInt     i,j;

884:   /* row-oriented output */
885:   for (i=0; i<m; i++) {
886:     if (indexm[i] < 0) {v += n;continue;}
887:     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);
888:     for (j=0; j<n; j++) {
889:       if (indexn[j] < 0) {v++; continue;}
890:       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);
891:       *v++ = mat->v[indexn[j]*mat->lda + indexm[i]];
892:     }
893:   }
894:   return(0);
895: }

897: /* -----------------------------------------------------------------*/

899: static PetscErrorCode MatLoad_SeqDense(Mat newmat,PetscViewer viewer)
900: {
901:   Mat_SeqDense   *a;
903:   PetscInt       *scols,i,j,nz,header[4];
904:   int            fd;
905:   PetscMPIInt    size;
906:   PetscInt       *rowlengths = 0,M,N,*cols,grows,gcols;
907:   PetscScalar    *vals,*svals,*v,*w;
908:   MPI_Comm       comm;

911:   /* force binary viewer to load .info file if it has not yet done so */
912:   PetscViewerSetUp(viewer);
913:   PetscObjectGetComm((PetscObject)viewer,&comm);
914:   MPI_Comm_size(comm,&size);
915:   if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor");
916:   PetscViewerBinaryGetDescriptor(viewer,&fd);
917:   PetscBinaryRead(fd,header,4,PETSC_INT);
918:   if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not matrix object");
919:   M = header[1]; N = header[2]; nz = header[3];

921:   /* set global size if not set already*/
922:   if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) {
923:     MatSetSizes(newmat,M,N,M,N);
924:   } else {
925:     /* if sizes and type are already set, check if the vector global sizes are correct */
926:     MatGetSize(newmat,&grows,&gcols);
927:     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);
928:   }
929:   a = (Mat_SeqDense*)newmat->data;
930:   if (!a->user_alloc) {
931:     MatSeqDenseSetPreallocation(newmat,NULL);
932:   }

934:   if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */
935:     a = (Mat_SeqDense*)newmat->data;
936:     v = a->v;
937:     /* Allocate some temp space to read in the values and then flip them
938:        from row major to column major */
939:     PetscMalloc1(M*N > 0 ? M*N : 1,&w);
940:     /* read in nonzero values */
941:     PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);
942:     /* now flip the values and store them in the matrix*/
943:     for (j=0; j<N; j++) {
944:       for (i=0; i<M; i++) {
945:         *v++ =w[i*N+j];
946:       }
947:     }
948:     PetscFree(w);
949:   } else {
950:     /* read row lengths */
951:     PetscMalloc1(M+1,&rowlengths);
952:     PetscBinaryRead(fd,rowlengths,M,PETSC_INT);

954:     a = (Mat_SeqDense*)newmat->data;
955:     v = a->v;

957:     /* read column indices and nonzeros */
958:     PetscMalloc1(nz+1,&scols);
959:     cols = scols;
960:     PetscBinaryRead(fd,cols,nz,PETSC_INT);
961:     PetscMalloc1(nz+1,&svals);
962:     vals = svals;
963:     PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);

965:     /* insert into matrix */
966:     for (i=0; i<M; i++) {
967:       for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j];
968:       svals += rowlengths[i]; scols += rowlengths[i];
969:     }
970:     PetscFree(vals);
971:     PetscFree(cols);
972:     PetscFree(rowlengths);
973:   }
974:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
975:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
976:   return(0);
977: }

979: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
980: {
981:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
982:   PetscErrorCode    ierr;
983:   PetscInt          i,j;
984:   const char        *name;
985:   PetscScalar       *v;
986:   PetscViewerFormat format;
987: #if defined(PETSC_USE_COMPLEX)
988:   PetscBool         allreal = PETSC_TRUE;
989: #endif

992:   PetscViewerGetFormat(viewer,&format);
993:   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
994:     return(0);  /* do nothing for now */
995:   } else if (format == PETSC_VIEWER_ASCII_COMMON) {
996:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
997:     for (i=0; i<A->rmap->n; i++) {
998:       v    = a->v + i;
999:       PetscViewerASCIIPrintf(viewer,"row %D:",i);
1000:       for (j=0; j<A->cmap->n; j++) {
1001: #if defined(PETSC_USE_COMPLEX)
1002:         if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
1003:           PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1004:         } else if (PetscRealPart(*v)) {
1005:           PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));
1006:         }
1007: #else
1008:         if (*v) {
1009:           PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);
1010:         }
1011: #endif
1012:         v += a->lda;
1013:       }
1014:       PetscViewerASCIIPrintf(viewer,"\n");
1015:     }
1016:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1017:   } else {
1018:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1019: #if defined(PETSC_USE_COMPLEX)
1020:     /* determine if matrix has all real values */
1021:     v = a->v;
1022:     for (i=0; i<A->rmap->n*A->cmap->n; i++) {
1023:       if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;}
1024:     }
1025: #endif
1026:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
1027:       PetscObjectGetName((PetscObject)A,&name);
1028:       PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);
1029:       PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);
1030:       PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
1031:     }

1033:     for (i=0; i<A->rmap->n; i++) {
1034:       v = a->v + i;
1035:       for (j=0; j<A->cmap->n; j++) {
1036: #if defined(PETSC_USE_COMPLEX)
1037:         if (allreal) {
1038:           PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));
1039:         } else {
1040:           PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1041:         }
1042: #else
1043:         PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);
1044: #endif
1045:         v += a->lda;
1046:       }
1047:       PetscViewerASCIIPrintf(viewer,"\n");
1048:     }
1049:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
1050:       PetscViewerASCIIPrintf(viewer,"];\n");
1051:     }
1052:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1053:   }
1054:   PetscViewerFlush(viewer);
1055:   return(0);
1056: }

1058: static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer)
1059: {
1060:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
1061:   PetscErrorCode    ierr;
1062:   int               fd;
1063:   PetscInt          ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n;
1064:   PetscScalar       *v,*anonz,*vals;
1065:   PetscViewerFormat format;

1068:   PetscViewerBinaryGetDescriptor(viewer,&fd);

1070:   PetscViewerGetFormat(viewer,&format);
1071:   if (format == PETSC_VIEWER_NATIVE) {
1072:     /* store the matrix as a dense matrix */
1073:     PetscMalloc1(4,&col_lens);

1075:     col_lens[0] = MAT_FILE_CLASSID;
1076:     col_lens[1] = m;
1077:     col_lens[2] = n;
1078:     col_lens[3] = MATRIX_BINARY_FORMAT_DENSE;

1080:     PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);
1081:     PetscFree(col_lens);

1083:     /* write out matrix, by rows */
1084:     PetscMalloc1(m*n+1,&vals);
1085:     v    = a->v;
1086:     for (j=0; j<n; j++) {
1087:       for (i=0; i<m; i++) {
1088:         vals[j + i*n] = *v++;
1089:       }
1090:     }
1091:     PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);
1092:     PetscFree(vals);
1093:   } else {
1094:     PetscMalloc1(4+nz,&col_lens);

1096:     col_lens[0] = MAT_FILE_CLASSID;
1097:     col_lens[1] = m;
1098:     col_lens[2] = n;
1099:     col_lens[3] = nz;

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

1105:     /* Possibly should write in smaller increments, not whole matrix at once? */
1106:     /* store column indices (zero start index) */
1107:     ict = 0;
1108:     for (i=0; i<m; i++) {
1109:       for (j=0; j<n; j++) col_lens[ict++] = j;
1110:     }
1111:     PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);
1112:     PetscFree(col_lens);

1114:     /* store nonzero values */
1115:     PetscMalloc1(nz+1,&anonz);
1116:     ict  = 0;
1117:     for (i=0; i<m; i++) {
1118:       v = a->v + i;
1119:       for (j=0; j<n; j++) {
1120:         anonz[ict++] = *v; v += a->lda;
1121:       }
1122:     }
1123:     PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);
1124:     PetscFree(anonz);
1125:   }
1126:   return(0);
1127: }

1129:  #include <petscdraw.h>
1130: static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa)
1131: {
1132:   Mat               A  = (Mat) Aa;
1133:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
1134:   PetscErrorCode    ierr;
1135:   PetscInt          m  = A->rmap->n,n = A->cmap->n,i,j;
1136:   int               color = PETSC_DRAW_WHITE;
1137:   PetscScalar       *v = a->v;
1138:   PetscViewer       viewer;
1139:   PetscReal         xl,yl,xr,yr,x_l,x_r,y_l,y_r;
1140:   PetscViewerFormat format;

1143:   PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
1144:   PetscViewerGetFormat(viewer,&format);
1145:   PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);

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

1149:   if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1150:     PetscDrawCollectiveBegin(draw);
1151:     /* Blue for negative and Red for positive */
1152:     for (j = 0; j < n; j++) {
1153:       x_l = j; x_r = x_l + 1.0;
1154:       for (i = 0; i < m; i++) {
1155:         y_l = m - i - 1.0;
1156:         y_r = y_l + 1.0;
1157:         if (PetscRealPart(v[j*m+i]) >  0.) {
1158:           color = PETSC_DRAW_RED;
1159:         } else if (PetscRealPart(v[j*m+i]) <  0.) {
1160:           color = PETSC_DRAW_BLUE;
1161:         } else {
1162:           continue;
1163:         }
1164:         PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1165:       }
1166:     }
1167:     PetscDrawCollectiveEnd(draw);
1168:   } else {
1169:     /* use contour shading to indicate magnitude of values */
1170:     /* first determine max of all nonzero values */
1171:     PetscReal minv = 0.0, maxv = 0.0;
1172:     PetscDraw popup;

1174:     for (i=0; i < m*n; i++) {
1175:       if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1176:     }
1177:     if (minv >= maxv) maxv = minv + PETSC_SMALL;
1178:     PetscDrawGetPopup(draw,&popup);
1179:     PetscDrawScalePopup(popup,minv,maxv);

1181:     PetscDrawCollectiveBegin(draw);
1182:     for (j=0; j<n; j++) {
1183:       x_l = j;
1184:       x_r = x_l + 1.0;
1185:       for (i=0; i<m; i++) {
1186:         y_l = m - i - 1.0;
1187:         y_r = y_l + 1.0;
1188:         color = PetscDrawRealToColor(PetscAbsScalar(v[j*m+i]),minv,maxv);
1189:         PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1190:       }
1191:     }
1192:     PetscDrawCollectiveEnd(draw);
1193:   }
1194:   return(0);
1195: }

1197: static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
1198: {
1199:   PetscDraw      draw;
1200:   PetscBool      isnull;
1201:   PetscReal      xr,yr,xl,yl,h,w;

1205:   PetscViewerDrawGetDraw(viewer,0,&draw);
1206:   PetscDrawIsNull(draw,&isnull);
1207:   if (isnull) return(0);

1209:   xr   = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
1210:   xr  += w;          yr += h;        xl = -w;     yl = -h;
1211:   PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
1212:   PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
1213:   PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
1214:   PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
1215:   PetscDrawSave(draw);
1216:   return(0);
1217: }

1219: PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
1220: {
1222:   PetscBool      iascii,isbinary,isdraw;

1225:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1226:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1227:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);

1229:   if (iascii) {
1230:     MatView_SeqDense_ASCII(A,viewer);
1231:   } else if (isbinary) {
1232:     MatView_SeqDense_Binary(A,viewer);
1233:   } else if (isdraw) {
1234:     MatView_SeqDense_Draw(A,viewer);
1235:   }
1236:   return(0);
1237: }

1239: static PetscErrorCode MatDestroy_SeqDense(Mat mat)
1240: {
1241:   Mat_SeqDense   *l = (Mat_SeqDense*)mat->data;

1245: #if defined(PETSC_USE_LOG)
1246:   PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n);
1247: #endif
1248:   PetscFree(l->pivots);
1249:   MatDestroy(&l->ptapwork);
1250:   if (!l->user_alloc) {PetscFree(l->v);}
1251:   PetscFree(mat->data);

1253:   PetscObjectChangeTypeName((PetscObject)mat,0);
1254:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);
1255:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);
1256:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);
1257: #if defined(PETSC_HAVE_ELEMENTAL)
1258:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);
1259: #endif
1260:   PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);
1261:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);
1262:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);
1263:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);
1264:   PetscObjectComposeFunction((PetscObject)mat,"MatPtAP_seqaij_seqdense_C",NULL);
1265:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);
1266:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);
1267:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);
1268:   return(0);
1269: }

1271: static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
1272: {
1273:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1275:   PetscInt       k,j,m,n,M;
1276:   PetscScalar    *v,tmp;

1279:   v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n;
1280:   if (reuse == MAT_INPLACE_MATRIX) { /* in place transpose */
1281:     if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can not transpose non-square matrix in place");
1282:     else {
1283:       for (j=0; j<m; j++) {
1284:         for (k=0; k<j; k++) {
1285:           tmp        = v[j + k*M];
1286:           v[j + k*M] = v[k + j*M];
1287:           v[k + j*M] = tmp;
1288:         }
1289:       }
1290:     }
1291:   } else { /* out-of-place transpose */
1292:     Mat          tmat;
1293:     Mat_SeqDense *tmatd;
1294:     PetscScalar  *v2;
1295:     PetscInt     M2;

1297:     if (reuse == MAT_INITIAL_MATRIX) {
1298:       MatCreate(PetscObjectComm((PetscObject)A),&tmat);
1299:       MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);
1300:       MatSetType(tmat,((PetscObject)A)->type_name);
1301:       MatSeqDenseSetPreallocation(tmat,NULL);
1302:     } else {
1303:       tmat = *matout;
1304:     }
1305:     tmatd = (Mat_SeqDense*)tmat->data;
1306:     v = mat->v; v2 = tmatd->v; M2 = tmatd->lda;
1307:     for (j=0; j<n; j++) {
1308:       for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M];
1309:     }
1310:     MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
1311:     MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
1312:     *matout = tmat;
1313:   }
1314:   return(0);
1315: }

1317: static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool  *flg)
1318: {
1319:   Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data;
1320:   Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data;
1321:   PetscInt     i,j;
1322:   PetscScalar  *v1,*v2;

1325:   if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; return(0);}
1326:   if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; return(0);}
1327:   for (i=0; i<A1->rmap->n; i++) {
1328:     v1 = mat1->v+i; v2 = mat2->v+i;
1329:     for (j=0; j<A1->cmap->n; j++) {
1330:       if (*v1 != *v2) {*flg = PETSC_FALSE; return(0);}
1331:       v1 += mat1->lda; v2 += mat2->lda;
1332:     }
1333:   }
1334:   *flg = PETSC_TRUE;
1335:   return(0);
1336: }

1338: static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
1339: {
1340:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1342:   PetscInt       i,n,len;
1343:   PetscScalar    *x,zero = 0.0;

1346:   VecSet(v,zero);
1347:   VecGetSize(v,&n);
1348:   VecGetArray(v,&x);
1349:   len  = PetscMin(A->rmap->n,A->cmap->n);
1350:   if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
1351:   for (i=0; i<len; i++) {
1352:     x[i] = mat->v[i*mat->lda + i];
1353:   }
1354:   VecRestoreArray(v,&x);
1355:   return(0);
1356: }

1358: static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
1359: {
1360:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1361:   const PetscScalar *l,*r;
1362:   PetscScalar       x,*v;
1363:   PetscErrorCode    ierr;
1364:   PetscInt          i,j,m = A->rmap->n,n = A->cmap->n;

1367:   if (ll) {
1368:     VecGetSize(ll,&m);
1369:     VecGetArrayRead(ll,&l);
1370:     if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size");
1371:     for (i=0; i<m; i++) {
1372:       x = l[i];
1373:       v = mat->v + i;
1374:       for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;}
1375:     }
1376:     VecRestoreArrayRead(ll,&l);
1377:     PetscLogFlops(1.0*n*m);
1378:   }
1379:   if (rr) {
1380:     VecGetSize(rr,&n);
1381:     VecGetArrayRead(rr,&r);
1382:     if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size");
1383:     for (i=0; i<n; i++) {
1384:       x = r[i];
1385:       v = mat->v + i*mat->lda;
1386:       for (j=0; j<m; j++) (*v++) *= x;
1387:     }
1388:     VecRestoreArrayRead(rr,&r);
1389:     PetscLogFlops(1.0*n*m);
1390:   }
1391:   return(0);
1392: }

1394: static PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1395: {
1396:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1397:   PetscScalar    *v   = mat->v;
1398:   PetscReal      sum  = 0.0;
1399:   PetscInt       lda  =mat->lda,m=A->rmap->n,i,j;

1403:   if (type == NORM_FROBENIUS) {
1404:     if (lda>m) {
1405:       for (j=0; j<A->cmap->n; j++) {
1406:         v = mat->v+j*lda;
1407:         for (i=0; i<m; i++) {
1408:           sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1409:         }
1410:       }
1411:     } else {
1412: #if defined(PETSC_USE_REAL___FP16)
1413:       PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n;
1414:       *nrm = BLASnrm2_(&cnt,v,&one);
1415:     }
1416: #else
1417:       for (i=0; i<A->cmap->n*A->rmap->n; i++) {
1418:         sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1419:       }
1420:     }
1421:     *nrm = PetscSqrtReal(sum);
1422: #endif
1423:     PetscLogFlops(2.0*A->cmap->n*A->rmap->n);
1424:   } else if (type == NORM_1) {
1425:     *nrm = 0.0;
1426:     for (j=0; j<A->cmap->n; j++) {
1427:       v   = mat->v + j*mat->lda;
1428:       sum = 0.0;
1429:       for (i=0; i<A->rmap->n; i++) {
1430:         sum += PetscAbsScalar(*v);  v++;
1431:       }
1432:       if (sum > *nrm) *nrm = sum;
1433:     }
1434:     PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1435:   } else if (type == NORM_INFINITY) {
1436:     *nrm = 0.0;
1437:     for (j=0; j<A->rmap->n; j++) {
1438:       v   = mat->v + j;
1439:       sum = 0.0;
1440:       for (i=0; i<A->cmap->n; i++) {
1441:         sum += PetscAbsScalar(*v); v += mat->lda;
1442:       }
1443:       if (sum > *nrm) *nrm = sum;
1444:     }
1445:     PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1446:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
1447:   return(0);
1448: }

1450: static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
1451: {
1452:   Mat_SeqDense   *aij = (Mat_SeqDense*)A->data;

1456:   switch (op) {
1457:   case MAT_ROW_ORIENTED:
1458:     aij->roworiented = flg;
1459:     break;
1460:   case MAT_NEW_NONZERO_LOCATIONS:
1461:   case MAT_NEW_NONZERO_LOCATION_ERR:
1462:   case MAT_NEW_NONZERO_ALLOCATION_ERR:
1463:   case MAT_NEW_DIAGONALS:
1464:   case MAT_KEEP_NONZERO_PATTERN:
1465:   case MAT_IGNORE_OFF_PROC_ENTRIES:
1466:   case MAT_USE_HASH_TABLE:
1467:   case MAT_IGNORE_LOWER_TRIANGULAR:
1468:     PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
1469:     break;
1470:   case MAT_SPD:
1471:   case MAT_SYMMETRIC:
1472:   case MAT_STRUCTURALLY_SYMMETRIC:
1473:   case MAT_HERMITIAN:
1474:   case MAT_SYMMETRY_ETERNAL:
1475:     /* These options are handled directly by MatSetOption() */
1476:     break;
1477:   default:
1478:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]);
1479:   }
1480:   return(0);
1481: }

1483: static PetscErrorCode MatZeroEntries_SeqDense(Mat A)
1484: {
1485:   Mat_SeqDense   *l = (Mat_SeqDense*)A->data;
1487:   PetscInt       lda=l->lda,m=A->rmap->n,j;

1490:   if (lda>m) {
1491:     for (j=0; j<A->cmap->n; j++) {
1492:       PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));
1493:     }
1494:   } else {
1495:     PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1496:   }
1497:   return(0);
1498: }

1500: static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1501: {
1502:   PetscErrorCode    ierr;
1503:   Mat_SeqDense      *l = (Mat_SeqDense*)A->data;
1504:   PetscInt          m  = l->lda, n = A->cmap->n, i,j;
1505:   PetscScalar       *slot,*bb;
1506:   const PetscScalar *xx;

1509: #if defined(PETSC_USE_DEBUG)
1510:   for (i=0; i<N; i++) {
1511:     if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
1512:     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);
1513:   }
1514: #endif

1516:   /* fix right hand side if needed */
1517:   if (x && b) {
1518:     VecGetArrayRead(x,&xx);
1519:     VecGetArray(b,&bb);
1520:     for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
1521:     VecRestoreArrayRead(x,&xx);
1522:     VecRestoreArray(b,&bb);
1523:   }

1525:   for (i=0; i<N; i++) {
1526:     slot = l->v + rows[i];
1527:     for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
1528:   }
1529:   if (diag != 0.0) {
1530:     if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
1531:     for (i=0; i<N; i++) {
1532:       slot  = l->v + (m+1)*rows[i];
1533:       *slot = diag;
1534:     }
1535:   }
1536:   return(0);
1537: }

1539: static PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[])
1540: {
1541:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;

1544:   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");
1545:   *array = mat->v;
1546:   return(0);
1547: }

1549: static PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[])
1550: {
1552:   *array = 0; /* user cannot accidently use the array later */
1553:   return(0);
1554: }

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

1559:    Not Collective

1561:    Input Parameter:
1562: .  mat - a MATSEQDENSE or MATMPIDENSE matrix

1564:    Output Parameter:
1565: .   array - pointer to the data

1567:    Level: intermediate

1569: .seealso: MatDenseRestoreArray()
1570: @*/
1571: PetscErrorCode  MatDenseGetArray(Mat A,PetscScalar **array)
1572: {

1576:   PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));
1577:   return(0);
1578: }

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

1583:    Not Collective

1585:    Input Parameters:
1586: .  mat - a MATSEQDENSE or MATMPIDENSE matrix
1587: .  array - pointer to the data

1589:    Level: intermediate

1591: .seealso: MatDenseGetArray()
1592: @*/
1593: PetscErrorCode  MatDenseRestoreArray(Mat A,PetscScalar **array)
1594: {

1598:   PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));
1599:   return(0);
1600: }

1602: static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B)
1603: {
1604:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1606:   PetscInt       i,j,nrows,ncols;
1607:   const PetscInt *irow,*icol;
1608:   PetscScalar    *av,*bv,*v = mat->v;
1609:   Mat            newmat;

1612:   ISGetIndices(isrow,&irow);
1613:   ISGetIndices(iscol,&icol);
1614:   ISGetLocalSize(isrow,&nrows);
1615:   ISGetLocalSize(iscol,&ncols);

1617:   /* Check submatrixcall */
1618:   if (scall == MAT_REUSE_MATRIX) {
1619:     PetscInt n_cols,n_rows;
1620:     MatGetSize(*B,&n_rows,&n_cols);
1621:     if (n_rows != nrows || n_cols != ncols) {
1622:       /* resize the result matrix to match number of requested rows/columns */
1623:       MatSetSizes(*B,nrows,ncols,nrows,ncols);
1624:     }
1625:     newmat = *B;
1626:   } else {
1627:     /* Create and fill new matrix */
1628:     MatCreate(PetscObjectComm((PetscObject)A),&newmat);
1629:     MatSetSizes(newmat,nrows,ncols,nrows,ncols);
1630:     MatSetType(newmat,((PetscObject)A)->type_name);
1631:     MatSeqDenseSetPreallocation(newmat,NULL);
1632:   }

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

1637:   for (i=0; i<ncols; i++) {
1638:     av = v + mat->lda*icol[i];
1639:     for (j=0; j<nrows; j++) *bv++ = av[irow[j]];
1640:   }

1642:   /* Assemble the matrices so that the correct flags are set */
1643:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1644:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);

1646:   /* Free work space */
1647:   ISRestoreIndices(isrow,&irow);
1648:   ISRestoreIndices(iscol,&icol);
1649:   *B   = newmat;
1650:   return(0);
1651: }

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

1659:   if (scall == MAT_INITIAL_MATRIX) {
1660:     PetscCalloc1(n+1,B);
1661:   }

1663:   for (i=0; i<n; i++) {
1664:     MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
1665:   }
1666:   return(0);
1667: }

1669: static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
1670: {
1672:   return(0);
1673: }

1675: static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
1676: {
1678:   return(0);
1679: }

1681: static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
1682: {
1683:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data;
1685:   PetscInt       lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j;

1688:   /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
1689:   if (A->ops->copy != B->ops->copy) {
1690:     MatCopy_Basic(A,B,str);
1691:     return(0);
1692:   }
1693:   if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)");
1694:   if (lda1>m || lda2>m) {
1695:     for (j=0; j<n; j++) {
1696:       PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));
1697:     }
1698:   } else {
1699:     PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1700:   }
1701:   return(0);
1702: }

1704: static PetscErrorCode MatSetUp_SeqDense(Mat A)
1705: {

1709:    MatSeqDenseSetPreallocation(A,0);
1710:   return(0);
1711: }

1713: static PetscErrorCode MatConjugate_SeqDense(Mat A)
1714: {
1715:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1716:   PetscInt     i,nz = A->rmap->n*A->cmap->n;
1717:   PetscScalar  *aa = a->v;

1720:   for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
1721:   return(0);
1722: }

1724: static PetscErrorCode MatRealPart_SeqDense(Mat A)
1725: {
1726:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1727:   PetscInt     i,nz = A->rmap->n*A->cmap->n;
1728:   PetscScalar  *aa = a->v;

1731:   for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
1732:   return(0);
1733: }

1735: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
1736: {
1737:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1738:   PetscInt     i,nz = A->rmap->n*A->cmap->n;
1739:   PetscScalar  *aa = a->v;

1742:   for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
1743:   return(0);
1744: }

1746: /* ----------------------------------------------------------------*/
1747: PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1748: {

1752:   if (scall == MAT_INITIAL_MATRIX) {
1753:     PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
1754:     MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
1755:     PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
1756:   }
1757:   PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
1758:   MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);
1759:   PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
1760:   return(0);
1761: }

1763: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1764: {
1766:   PetscInt       m=A->rmap->n,n=B->cmap->n;
1767:   Mat            Cmat;

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

1776:   *C = Cmat;
1777:   return(0);
1778: }

1780: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
1781: {
1782:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
1783:   Mat_SeqDense   *b = (Mat_SeqDense*)B->data;
1784:   Mat_SeqDense   *c = (Mat_SeqDense*)C->data;
1785:   PetscBLASInt   m,n,k;
1786:   PetscScalar    _DOne=1.0,_DZero=0.0;
1788:   PetscBool      flg;

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

1794:   /* Handle case where where user provided the final C matrix rather than calling MatMatMult() with MAT_INITIAL_MATRIX*/
1795:   PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
1796:   if (flg) {
1797:     C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqDense;
1798:     (*C->ops->matmultnumeric)(A,B,C);
1799:     return(0);
1800:   }

1802:   PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&flg);
1803:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"First matrix must be dense");
1804:   PetscBLASIntCast(A->rmap->n,&m);
1805:   PetscBLASIntCast(B->cmap->n,&n);
1806:   PetscBLASIntCast(A->cmap->n,&k);
1807:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
1808:   return(0);
1809: }

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

1816:   if (scall == MAT_INITIAL_MATRIX) {
1817:     PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);
1818:     MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
1819:     PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);
1820:   }
1821:   PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);
1822:   MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);
1823:   PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);
1824:   return(0);
1825: }

1827: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1828: {
1830:   PetscInt       m=A->cmap->n,n=B->cmap->n;
1831:   Mat            Cmat;

1834:   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);
1835:   MatCreate(PETSC_COMM_SELF,&Cmat);
1836:   MatSetSizes(Cmat,m,n,m,n);
1837:   MatSetType(Cmat,MATSEQDENSE);
1838:   MatSeqDenseSetPreallocation(Cmat,NULL);

1840:   Cmat->assembled = PETSC_TRUE;

1842:   *C = Cmat;
1843:   return(0);
1844: }

1846: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
1847: {
1848:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
1849:   Mat_SeqDense   *b = (Mat_SeqDense*)B->data;
1850:   Mat_SeqDense   *c = (Mat_SeqDense*)C->data;
1851:   PetscBLASInt   m,n,k;
1852:   PetscScalar    _DOne=1.0,_DZero=0.0;

1856:   PetscBLASIntCast(A->cmap->n,&m);
1857:   PetscBLASIntCast(B->cmap->n,&n);
1858:   PetscBLASIntCast(A->rmap->n,&k);
1859:   /*
1860:      Note the m and n arguments below are the number rows and columns of A', not A!
1861:   */
1862:   PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
1863:   return(0);
1864: }

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

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

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

1891: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])
1892: {
1893:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
1895:   PetscInt       i,j,m = A->rmap->n,n = A->cmap->n,p;
1896:   PetscScalar    *x;
1897:   PetscReal      atmp;
1898:   MatScalar      *aa = a->v;

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

1903:   VecSet(v,0.0);
1904:   VecGetArray(v,&x);
1905:   VecGetLocalSize(v,&p);
1906:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
1907:   for (i=0; i<m; i++) {
1908:     x[i] = PetscAbsScalar(aa[i]);
1909:     for (j=1; j<n; j++) {
1910:       atmp = PetscAbsScalar(aa[i+m*j]);
1911:       if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;}
1912:     }
1913:   }
1914:   VecRestoreArray(v,&x);
1915:   return(0);
1916: }

1918: static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])
1919: {
1920:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
1922:   PetscInt       i,j,m = A->rmap->n,n = A->cmap->n,p;
1923:   PetscScalar    *x;
1924:   MatScalar      *aa = a->v;

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

1929:   VecSet(v,0.0);
1930:   VecGetArray(v,&x);
1931:   VecGetLocalSize(v,&p);
1932:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
1933:   for (i=0; i<m; i++) {
1934:     x[i] = aa[i]; if (idx) idx[i] = 0;
1935:     for (j=1; j<n; j++) {
1936:       if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
1937:     }
1938:   }
1939:   VecRestoreArray(v,&x);
1940:   return(0);
1941: }

1943: static PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
1944: {
1945:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
1947:   PetscScalar    *x;

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

1952:   VecGetArray(v,&x);
1953:   PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));
1954:   VecRestoreArray(v,&x);
1955:   return(0);
1956: }


1959: PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms)
1960: {
1962:   PetscInt       i,j,m,n;
1963:   PetscScalar    *a;

1966:   MatGetSize(A,&m,&n);
1967:   PetscMemzero(norms,n*sizeof(PetscReal));
1968:   MatDenseGetArray(A,&a);
1969:   if (type == NORM_2) {
1970:     for (i=0; i<n; i++) {
1971:       for (j=0; j<m; j++) {
1972:         norms[i] += PetscAbsScalar(a[j]*a[j]);
1973:       }
1974:       a += m;
1975:     }
1976:   } else if (type == NORM_1) {
1977:     for (i=0; i<n; i++) {
1978:       for (j=0; j<m; j++) {
1979:         norms[i] += PetscAbsScalar(a[j]);
1980:       }
1981:       a += m;
1982:     }
1983:   } else if (type == NORM_INFINITY) {
1984:     for (i=0; i<n; i++) {
1985:       for (j=0; j<m; j++) {
1986:         norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]);
1987:       }
1988:       a += m;
1989:     }
1990:   } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType");
1991:   MatDenseRestoreArray(A,&a);
1992:   if (type == NORM_2) {
1993:     for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]);
1994:   }
1995:   return(0);
1996: }

1998: static PetscErrorCode  MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
1999: {
2001:   PetscScalar    *a;
2002:   PetscInt       m,n,i;

2005:   MatGetSize(x,&m,&n);
2006:   MatDenseGetArray(x,&a);
2007:   for (i=0; i<m*n; i++) {
2008:     PetscRandomGetValue(rctx,a+i);
2009:   }
2010:   MatDenseRestoreArray(x,&a);
2011:   return(0);
2012: }

2014: static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool  *missing,PetscInt *d)
2015: {
2017:   *missing = PETSC_FALSE;
2018:   return(0);
2019: }


2022: /* -------------------------------------------------------------------*/
2023: static struct _MatOps MatOps_Values = { MatSetValues_SeqDense,
2024:                                         MatGetRow_SeqDense,
2025:                                         MatRestoreRow_SeqDense,
2026:                                         MatMult_SeqDense,
2027:                                 /*  4*/ MatMultAdd_SeqDense,
2028:                                         MatMultTranspose_SeqDense,
2029:                                         MatMultTransposeAdd_SeqDense,
2030:                                         0,
2031:                                         0,
2032:                                         0,
2033:                                 /* 10*/ 0,
2034:                                         MatLUFactor_SeqDense,
2035:                                         MatCholeskyFactor_SeqDense,
2036:                                         MatSOR_SeqDense,
2037:                                         MatTranspose_SeqDense,
2038:                                 /* 15*/ MatGetInfo_SeqDense,
2039:                                         MatEqual_SeqDense,
2040:                                         MatGetDiagonal_SeqDense,
2041:                                         MatDiagonalScale_SeqDense,
2042:                                         MatNorm_SeqDense,
2043:                                 /* 20*/ MatAssemblyBegin_SeqDense,
2044:                                         MatAssemblyEnd_SeqDense,
2045:                                         MatSetOption_SeqDense,
2046:                                         MatZeroEntries_SeqDense,
2047:                                 /* 24*/ MatZeroRows_SeqDense,
2048:                                         0,
2049:                                         0,
2050:                                         0,
2051:                                         0,
2052:                                 /* 29*/ MatSetUp_SeqDense,
2053:                                         0,
2054:                                         0,
2055:                                         0,
2056:                                         0,
2057:                                 /* 34*/ MatDuplicate_SeqDense,
2058:                                         0,
2059:                                         0,
2060:                                         0,
2061:                                         0,
2062:                                 /* 39*/ MatAXPY_SeqDense,
2063:                                         MatCreateSubMatrices_SeqDense,
2064:                                         0,
2065:                                         MatGetValues_SeqDense,
2066:                                         MatCopy_SeqDense,
2067:                                 /* 44*/ MatGetRowMax_SeqDense,
2068:                                         MatScale_SeqDense,
2069:                                         MatShift_Basic,
2070:                                         0,
2071:                                         MatZeroRowsColumns_SeqDense,
2072:                                 /* 49*/ MatSetRandom_SeqDense,
2073:                                         0,
2074:                                         0,
2075:                                         0,
2076:                                         0,
2077:                                 /* 54*/ 0,
2078:                                         0,
2079:                                         0,
2080:                                         0,
2081:                                         0,
2082:                                 /* 59*/ 0,
2083:                                         MatDestroy_SeqDense,
2084:                                         MatView_SeqDense,
2085:                                         0,
2086:                                         0,
2087:                                 /* 64*/ 0,
2088:                                         0,
2089:                                         0,
2090:                                         0,
2091:                                         0,
2092:                                 /* 69*/ MatGetRowMaxAbs_SeqDense,
2093:                                         0,
2094:                                         0,
2095:                                         0,
2096:                                         0,
2097:                                 /* 74*/ 0,
2098:                                         0,
2099:                                         0,
2100:                                         0,
2101:                                         0,
2102:                                 /* 79*/ 0,
2103:                                         0,
2104:                                         0,
2105:                                         0,
2106:                                 /* 83*/ MatLoad_SeqDense,
2107:                                         0,
2108:                                         MatIsHermitian_SeqDense,
2109:                                         0,
2110:                                         0,
2111:                                         0,
2112:                                 /* 89*/ MatMatMult_SeqDense_SeqDense,
2113:                                         MatMatMultSymbolic_SeqDense_SeqDense,
2114:                                         MatMatMultNumeric_SeqDense_SeqDense,
2115:                                         MatPtAP_SeqDense_SeqDense,
2116:                                         MatPtAPSymbolic_SeqDense_SeqDense,
2117:                                 /* 94*/ MatPtAPNumeric_SeqDense_SeqDense,
2118:                                         0,
2119:                                         0,
2120:                                         0,
2121:                                         0,
2122:                                 /* 99*/ 0,
2123:                                         0,
2124:                                         0,
2125:                                         MatConjugate_SeqDense,
2126:                                         0,
2127:                                 /*104*/ 0,
2128:                                         MatRealPart_SeqDense,
2129:                                         MatImaginaryPart_SeqDense,
2130:                                         0,
2131:                                         0,
2132:                                 /*109*/ MatMatSolve_SeqDense,
2133:                                         0,
2134:                                         MatGetRowMin_SeqDense,
2135:                                         MatGetColumnVector_SeqDense,
2136:                                         MatMissingDiagonal_SeqDense,
2137:                                 /*114*/ 0,
2138:                                         0,
2139:                                         0,
2140:                                         0,
2141:                                         0,
2142:                                 /*119*/ 0,
2143:                                         0,
2144:                                         0,
2145:                                         0,
2146:                                         0,
2147:                                 /*124*/ 0,
2148:                                         MatGetColumnNorms_SeqDense,
2149:                                         0,
2150:                                         0,
2151:                                         0,
2152:                                 /*129*/ 0,
2153:                                         MatTransposeMatMult_SeqDense_SeqDense,
2154:                                         MatTransposeMatMultSymbolic_SeqDense_SeqDense,
2155:                                         MatTransposeMatMultNumeric_SeqDense_SeqDense,
2156:                                         0,
2157:                                 /*134*/ 0,
2158:                                         0,
2159:                                         0,
2160:                                         0,
2161:                                         0,
2162:                                 /*139*/ 0,
2163:                                         0,
2164:                                         0
2165: };

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

2172:    Collective on MPI_Comm

2174:    Input Parameters:
2175: +  comm - MPI communicator, set to PETSC_COMM_SELF
2176: .  m - number of rows
2177: .  n - number of columns
2178: -  data - optional location of matrix data in column major order.  Set data=NULL for PETSc
2179:    to control all matrix memory allocation.

2181:    Output Parameter:
2182: .  A - the matrix

2184:    Notes:
2185:    The data input variable is intended primarily for Fortran programmers
2186:    who wish to allocate their own matrix memory space.  Most users should
2187:    set data=NULL.

2189:    Level: intermediate

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

2193: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
2194: @*/
2195: PetscErrorCode  MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
2196: {

2200:   MatCreate(comm,A);
2201:   MatSetSizes(*A,m,n,m,n);
2202:   MatSetType(*A,MATSEQDENSE);
2203:   MatSeqDenseSetPreallocation(*A,data);
2204:   return(0);
2205: }

2207: /*@C
2208:    MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements

2210:    Collective on MPI_Comm

2212:    Input Parameters:
2213: +  B - the matrix
2214: -  data - the array (or NULL)

2216:    Notes:
2217:    The data input variable is intended primarily for Fortran programmers
2218:    who wish to allocate their own matrix memory space.  Most users should
2219:    need not call this routine.

2221:    Level: intermediate

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

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

2227: @*/
2228: PetscErrorCode  MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
2229: {

2233:   PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
2234:   return(0);
2235: }

2237: PetscErrorCode  MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
2238: {
2239:   Mat_SeqDense   *b;

2243:   B->preallocated = PETSC_TRUE;

2245:   PetscLayoutSetUp(B->rmap);
2246:   PetscLayoutSetUp(B->cmap);

2248:   b       = (Mat_SeqDense*)B->data;
2249:   b->Mmax = B->rmap->n;
2250:   b->Nmax = B->cmap->n;
2251:   if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n;

2253:   PetscIntMultError(b->lda,b->Nmax,NULL);
2254:   if (!data) { /* petsc-allocated storage */
2255:     if (!b->user_alloc) { PetscFree(b->v); }
2256:     PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);
2257:     PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));

2259:     b->user_alloc = PETSC_FALSE;
2260:   } else { /* user-allocated storage */
2261:     if (!b->user_alloc) { PetscFree(b->v); }
2262:     b->v          = data;
2263:     b->user_alloc = PETSC_TRUE;
2264:   }
2265:   B->assembled = PETSC_TRUE;
2266:   return(0);
2267: }

2269: #if defined(PETSC_HAVE_ELEMENTAL)
2270: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
2271: {
2272:   Mat            mat_elemental;
2274:   PetscScalar    *array,*v_colwise;
2275:   PetscInt       M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols;

2278:   PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);
2279:   MatDenseGetArray(A,&array);
2280:   /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
2281:   k = 0;
2282:   for (j=0; j<N; j++) {
2283:     cols[j] = j;
2284:     for (i=0; i<M; i++) {
2285:       v_colwise[j*M+i] = array[k++];
2286:     }
2287:   }
2288:   for (i=0; i<M; i++) {
2289:     rows[i] = i;
2290:   }
2291:   MatDenseRestoreArray(A,&array);

2293:   MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
2294:   MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);
2295:   MatSetType(mat_elemental,MATELEMENTAL);
2296:   MatSetUp(mat_elemental);

2298:   /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
2299:   MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);
2300:   MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
2301:   MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);
2302:   PetscFree3(v_colwise,rows,cols);

2304:   if (reuse == MAT_INPLACE_MATRIX) {
2305:     MatHeaderReplace(A,&mat_elemental);
2306:   } else {
2307:     *newmat = mat_elemental;
2308:   }
2309:   return(0);
2310: }
2311: #endif

2313: /*@C
2314:   MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array

2316:   Input parameter:
2317: + A - the matrix
2318: - lda - the leading dimension

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

2325:   Level: intermediate

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

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

2331: @*/
2332: PetscErrorCode  MatSeqDenseSetLDA(Mat B,PetscInt lda)
2333: {
2334:   Mat_SeqDense *b = (Mat_SeqDense*)B->data;

2337:   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);
2338:   b->lda       = lda;
2339:   b->changelda = PETSC_FALSE;
2340:   b->Mmax      = PetscMax(b->Mmax,lda);
2341:   return(0);
2342: }

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

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

2350:   Level: beginner

2352: .seealso: MatCreateSeqDense()

2354: M*/

2356: PETSC_EXTERN PetscErrorCode MatCreate_SeqDense(Mat B)
2357: {
2358:   Mat_SeqDense   *b;
2360:   PetscMPIInt    size;

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

2366:   PetscNewLog(B,&b);
2367:   PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
2368:   B->data = (void*)b;

2370:   b->pivots      = 0;
2371:   b->roworiented = PETSC_TRUE;
2372:   b->v           = 0;
2373:   b->changelda   = PETSC_FALSE;

2375:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);
2376:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);
2377:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);
2378: #if defined(PETSC_HAVE_ELEMENTAL)
2379:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);
2380: #endif
2381:   PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);
2382:   PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2383:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2384:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2385:   PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaij_seqdense_C",MatPtAP_SeqDense_SeqDense);
2386:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2387:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2388:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2389:   PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
2390:   return(0);
2391: }