Actual source code: dense.c

petsc-master 2020-05-26
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 MatSeqDenseSymmetrize_Private(Mat A, PetscBool hermitian)
 12: {
 13:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
 14:   PetscInt       j, k, n = A->rmap->n;
 15:   PetscScalar    *v;

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

 38: PETSC_EXTERN PetscErrorCode MatSeqDenseInvertFactors_Private(Mat A)
 39: {
 40:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
 42:   PetscBLASInt   info,n;

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

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

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

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

114:   /* fix right hand side if needed */
115:   if (x && b) {
116:     Vec xt;

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

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

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

157:   if (c->ptapwork) {
158:     (*C->ops->matmultnumeric)(A,P,c->ptapwork);
159:     (*C->ops->transposematmultnumeric)(P,c->ptapwork,C);
160:   } else SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"Must call MatPtAPSymbolic_SeqDense_SeqDense() first");
161:   return(0);
162: }

164: PetscErrorCode MatPtAPSymbolic_SeqDense_SeqDense(Mat A,Mat P,PetscReal fill,Mat C)
165: {
166:   Mat_SeqDense   *c;
167:   PetscBool      flg;

171:   PetscObjectTypeCompare((PetscObject)P,((PetscObject)A)->type_name,&flg);
172:   MatSetSizes(C,P->cmap->n,P->cmap->n,P->cmap->N,P->cmap->N);
173:   MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
174:   MatSeqDenseSetPreallocation(C,NULL);
175:   c    = (Mat_SeqDense*)C->data;
176:   MatCreate(PetscObjectComm((PetscObject)A),&c->ptapwork);
177:   MatSetSizes(c->ptapwork,A->rmap->n,P->cmap->n,A->rmap->N,P->cmap->N);
178:   MatSetType(c->ptapwork,flg ? ((PetscObject)A)->type_name : MATDENSE);
179:   MatSeqDenseSetPreallocation(c->ptapwork,NULL);
180:   return(0);
181: }

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

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

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

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

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

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

268: PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str)
269: {
270:   Mat_SeqDense      *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data;
271:   const PetscScalar *xv;
272:   PetscScalar       *yv;
273:   PetscBLASInt      N,m,ldax,lday,one = 1;
274:   PetscErrorCode    ierr;

277:   MatDenseGetArrayRead(X,&xv);
278:   MatDenseGetArray(Y,&yv);
279:   PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);
280:   PetscBLASIntCast(X->rmap->n,&m);
281:   PetscBLASIntCast(x->lda,&ldax);
282:   PetscBLASIntCast(y->lda,&lday);
283:   if (ldax>m || lday>m) {
284:     PetscInt j;

286:     for (j=0; j<X->cmap->n; j++) {
287:       PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&alpha,xv+j*ldax,&one,yv+j*lday,&one));
288:     }
289:   } else {
290:     PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&alpha,xv,&one,yv,&one));
291:   }
292:   MatDenseRestoreArrayRead(X,&xv);
293:   MatDenseRestoreArray(Y,&yv);
294:   PetscLogFlops(PetscMax(2*N-1,0));
295:   return(0);
296: }

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

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

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

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

340: static PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool  *fl)
341: {
342:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
343:   PetscInt          i,j,m = A->rmap->n,N = a->lda;
344:   const PetscScalar *v;
345:   PetscErrorCode    ierr;

348:   *fl = PETSC_FALSE;
349:   if (A->rmap->n != A->cmap->n) return(0);
350:   MatDenseGetArrayRead(A,&v);
351:   for (i=0; i<m; i++) {
352:     for (j=i; j<m; j++) {
353:       if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) return(0);
354:     }
355:   }
356:   MatDenseRestoreArrayRead(A,&v);
357:   *fl  = PETSC_TRUE;
358:   return(0);
359: }

361: PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues)
362: {
363:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
365:   PetscInt       lda = (PetscInt)mat->lda,j,m;

368:   PetscLayoutReference(A->rmap,&newi->rmap);
369:   PetscLayoutReference(A->cmap,&newi->cmap);
370:   MatSeqDenseSetPreallocation(newi,NULL);
371:   if (cpvalues == MAT_COPY_VALUES) {
372:     const PetscScalar *av;
373:     PetscScalar       *v;

375:     MatDenseGetArrayRead(A,&av);
376:     MatDenseGetArray(newi,&v);
377:     if (lda>A->rmap->n) {
378:       m = A->rmap->n;
379:       for (j=0; j<A->cmap->n; j++) {
380:         PetscArraycpy(v+j*m,av+j*lda,m);
381:       }
382:     } else {
383:       PetscArraycpy(v,av,A->rmap->n*A->cmap->n);
384:     }
385:     MatDenseRestoreArray(newi,&v);
386:     MatDenseRestoreArrayRead(A,&av);
387:   }
388:   return(0);
389: }

391: PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat)
392: {

396:   MatCreate(PetscObjectComm((PetscObject)A),newmat);
397:   MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
398:   MatSetType(*newmat,((PetscObject)A)->type_name);
399:   MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);
400:   return(0);
401: }

403: static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
404: {
405:   MatFactorInfo  info;

409:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
410:   (*fact->ops->lufactor)(fact,0,0,&info);
411:   return(0);
412: }

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

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

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

468:   PetscBLASIntCast(A->rmap->n,&m);
469:   MatGetSize(B,NULL,&n);
470:   PetscBLASIntCast(n,&nrhs);
471:   MatDenseGetArrayRead(B,&b);
472:   MatDenseGetArray(X,&x);

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

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

502:   MatDenseRestoreArrayRead(B,&b);
503:   MatDenseRestoreArray(X,&x);
504:   PetscLogFlops(nrhs*(2.0*m*m - m));
505:   return(0);
506: }

508: static PetscErrorCode MatConjugate_SeqDense(Mat);

510: static PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy)
511: {
512:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
513:   PetscErrorCode    ierr;
514:   const PetscScalar *x;
515:   PetscScalar       *y;
516:   PetscBLASInt      one = 1,info,m;

519:   PetscBLASIntCast(A->rmap->n,&m);
520:   VecGetArrayRead(xx,&x);
521:   VecGetArray(yy,&y);
522:   PetscArraycpy(y,x,A->rmap->n);
523:   if (A->factortype == MAT_FACTOR_LU) {
524:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
525:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
526:     PetscFPTrapPop();
527:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve");
528:   } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
529:     if (A->spd) {
530: #if defined(PETSC_USE_COMPLEX)
531:       MatConjugate_SeqDense(A);
532: #endif
533:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
534:       PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
535:       PetscFPTrapPop();
536: #if defined(PETSC_USE_COMPLEX)
537:       MatConjugate_SeqDense(A);
538: #endif
539:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
540: #if defined(PETSC_USE_COMPLEX)
541:     } else if (A->hermitian) {
542:       MatConjugate_SeqDense(A);
543:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
544:       PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
545:       PetscFPTrapPop();
546:       MatConjugate_SeqDense(A);
547: #endif
548:     } else { /* symmetric case */
549:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
550:       PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
551:       PetscFPTrapPop();
552:       if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
553:     }
554:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
555:   VecRestoreArrayRead(xx,&x);
556:   VecRestoreArray(yy,&y);
557:   PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
558:   return(0);
559: }

561: /* ---------------------------------------------------------------*/
562: /* COMMENT: I have chosen to hide row permutation in the pivots,
563:    rather than put it in the Mat->row slot.*/
564: PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo)
565: {
566:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
568:   PetscBLASInt   n,m,info;

571:   PetscBLASIntCast(A->cmap->n,&n);
572:   PetscBLASIntCast(A->rmap->n,&m);
573:   if (!mat->pivots) {
574:     PetscMalloc1(A->rmap->n,&mat->pivots);
575:     PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
576:   }
577:   if (!A->rmap->n || !A->cmap->n) return(0);
578:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
579:   PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info));
580:   PetscFPTrapPop();

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

585:   A->ops->solve             = MatSolve_SeqDense;
586:   A->ops->matsolve          = MatMatSolve_SeqDense;
587:   A->ops->solvetranspose    = MatSolveTranspose_SeqDense;
588:   A->factortype             = MAT_FACTOR_LU;

590:   PetscFree(A->solvertype);
591:   PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);

593:   PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);
594:   return(0);
595: }

597: /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */
598: PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo)
599: {
600:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
602:   PetscBLASInt   info,n;

605:   PetscBLASIntCast(A->cmap->n,&n);
606:   if (!A->rmap->n || !A->cmap->n) return(0);
607:   if (A->spd) {
608:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
609:     PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info));
610:     PetscFPTrapPop();
611: #if defined(PETSC_USE_COMPLEX)
612:   } else if (A->hermitian) {
613:     if (!mat->pivots) {
614:       PetscMalloc1(A->rmap->n,&mat->pivots);
615:       PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
616:     }
617:     if (!mat->fwork) {
618:       PetscScalar dummy;

620:       mat->lfwork = -1;
621:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
622:       PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
623:       PetscFPTrapPop();
624:       mat->lfwork = (PetscInt)PetscRealPart(dummy);
625:       PetscMalloc1(mat->lfwork,&mat->fwork);
626:       PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
627:     }
628:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
629:     PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
630:     PetscFPTrapPop();
631: #endif
632:   } else { /* symmetric case */
633:     if (!mat->pivots) {
634:       PetscMalloc1(A->rmap->n,&mat->pivots);
635:       PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
636:     }
637:     if (!mat->fwork) {
638:       PetscScalar dummy;

640:       mat->lfwork = -1;
641:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
642:       PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
643:       PetscFPTrapPop();
644:       mat->lfwork = (PetscInt)PetscRealPart(dummy);
645:       PetscMalloc1(mat->lfwork,&mat->fwork);
646:       PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
647:     }
648:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
649:     PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
650:     PetscFPTrapPop();
651:   }
652:   if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1);

654:   A->ops->solve             = MatSolve_SeqDense;
655:   A->ops->matsolve          = MatMatSolve_SeqDense;
656:   A->ops->solvetranspose    = MatSolveTranspose_SeqDense;
657:   A->factortype             = MAT_FACTOR_CHOLESKY;

659:   PetscFree(A->solvertype);
660:   PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);

662:   PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
663:   return(0);
664: }


667: PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
668: {
670:   MatFactorInfo  info;

673:   info.fill = 1.0;

675:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
676:   (*fact->ops->choleskyfactor)(fact,0,&info);
677:   return(0);
678: }

680: PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
681: {
683:   fact->assembled                  = PETSC_TRUE;
684:   fact->preallocated               = PETSC_TRUE;
685:   fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
686:   fact->ops->solve                 = MatSolve_SeqDense;
687:   fact->ops->matsolve              = MatMatSolve_SeqDense;
688:   fact->ops->solvetranspose        = MatSolveTranspose_SeqDense;
689:   return(0);
690: }

692: PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info)
693: {
695:   fact->preallocated           = PETSC_TRUE;
696:   fact->assembled              = PETSC_TRUE;
697:   fact->ops->lufactornumeric   = MatLUFactorNumeric_SeqDense;
698:   fact->ops->solve             = MatSolve_SeqDense;
699:   fact->ops->matsolve          = MatMatSolve_SeqDense;
700:   fact->ops->solvetranspose    = MatSolveTranspose_SeqDense;
701:   return(0);
702: }

704: /* uses LAPACK */
705: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact)
706: {

710:   MatCreate(PetscObjectComm((PetscObject)A),fact);
711:   MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
712:   MatSetType(*fact,MATDENSE);
713:   if (ftype == MAT_FACTOR_LU) {
714:     (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
715:   } else {
716:     (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
717:   }
718:   (*fact)->factortype = ftype;

720:   PetscFree((*fact)->solvertype);
721:   PetscStrallocpy(MATSOLVERPETSC,&(*fact)->solvertype);
722:   return(0);
723: }

725: /* ------------------------------------------------------------------*/
726: static PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx)
727: {
728:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
729:   PetscScalar       *x,*v = mat->v,zero = 0.0,xt;
730:   const PetscScalar *b;
731:   PetscErrorCode    ierr;
732:   PetscInt          m = A->rmap->n,i;
733:   PetscBLASInt      o = 1,bm;

736: #if defined(PETSC_HAVE_CUDA)
737:   if (A->offloadmask == PETSC_OFFLOAD_GPU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented");
738: #endif
739:   if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */
740:   PetscBLASIntCast(m,&bm);
741:   if (flag & SOR_ZERO_INITIAL_GUESS) {
742:     /* this is a hack fix, should have another version without the second BLASdotu */
743:     VecSet(xx,zero);
744:   }
745:   VecGetArray(xx,&x);
746:   VecGetArrayRead(bb,&b);
747:   its  = its*lits;
748:   if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
749:   while (its--) {
750:     if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
751:       for (i=0; i<m; i++) {
752:         PetscStackCallBLAS("BLASdotu",xt   = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o));
753:         x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
754:       }
755:     }
756:     if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
757:       for (i=m-1; i>=0; i--) {
758:         PetscStackCallBLAS("BLASdotu",xt   = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o));
759:         x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
760:       }
761:     }
762:   }
763:   VecRestoreArrayRead(bb,&b);
764:   VecRestoreArray(xx,&x);
765:   return(0);
766: }

768: /* -----------------------------------------------------------------*/
769: PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)
770: {
771:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
772:   const PetscScalar *v   = mat->v,*x;
773:   PetscScalar       *y;
774:   PetscErrorCode    ierr;
775:   PetscBLASInt      m, n,_One=1;
776:   PetscScalar       _DOne=1.0,_DZero=0.0;

779:   PetscBLASIntCast(A->rmap->n,&m);
780:   PetscBLASIntCast(A->cmap->n,&n);
781:   VecGetArrayRead(xx,&x);
782:   VecGetArrayWrite(yy,&y);
783:   if (!A->rmap->n || !A->cmap->n) {
784:     PetscBLASInt i;
785:     for (i=0; i<n; i++) y[i] = 0.0;
786:   } else {
787:     PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One));
788:     PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);
789:   }
790:   VecRestoreArrayRead(xx,&x);
791:   VecRestoreArrayWrite(yy,&y);
792:   return(0);
793: }

795: PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy)
796: {
797:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
798:   PetscScalar       *y,_DOne=1.0,_DZero=0.0;
799:   PetscErrorCode    ierr;
800:   PetscBLASInt      m, n, _One=1;
801:   const PetscScalar *v = mat->v,*x;

804:   PetscBLASIntCast(A->rmap->n,&m);
805:   PetscBLASIntCast(A->cmap->n,&n);
806:   VecGetArrayRead(xx,&x);
807:   VecGetArrayWrite(yy,&y);
808:   if (!A->rmap->n || !A->cmap->n) {
809:     PetscBLASInt i;
810:     for (i=0; i<m; i++) y[i] = 0.0;
811:   } else {
812:     PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One));
813:     PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);
814:   }
815:   VecRestoreArrayRead(xx,&x);
816:   VecRestoreArrayWrite(yy,&y);
817:   return(0);
818: }

820: PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
821: {
822:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
823:   const PetscScalar *v = mat->v,*x;
824:   PetscScalar       *y,_DOne=1.0;
825:   PetscErrorCode    ierr;
826:   PetscBLASInt      m, n, _One=1;

829:   PetscBLASIntCast(A->rmap->n,&m);
830:   PetscBLASIntCast(A->cmap->n,&n);
831:   if (!A->rmap->n || !A->cmap->n) return(0);
832:   if (zz != yy) {VecCopy(zz,yy);}
833:   VecGetArrayRead(xx,&x);
834:   VecGetArray(yy,&y);
835:   PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
836:   VecRestoreArrayRead(xx,&x);
837:   VecRestoreArray(yy,&y);
838:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
839:   return(0);
840: }

842: PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
843: {
844:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
845:   const PetscScalar *v = mat->v,*x;
846:   PetscScalar       *y;
847:   PetscErrorCode    ierr;
848:   PetscBLASInt      m, n, _One=1;
849:   PetscScalar       _DOne=1.0;

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

865: /* -----------------------------------------------------------------*/
866: static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
867: {
868:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
870:   PetscInt       i;

873:   *ncols = A->cmap->n;
874:   if (cols) {
875:     PetscMalloc1(A->cmap->n+1,cols);
876:     for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
877:   }
878:   if (vals) {
879:     const PetscScalar *v;

881:     MatDenseGetArrayRead(A,&v);
882:     PetscMalloc1(A->cmap->n+1,vals);
883:     v   += row;
884:     for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
885:     MatDenseRestoreArrayRead(A,&v);
886:   }
887:   return(0);
888: }

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

895:   if (cols) {PetscFree(*cols);}
896:   if (vals) {PetscFree(*vals); }
897:   return(0);
898: }
899: /* ----------------------------------------------------------------*/
900: static PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
901: {
902:   Mat_SeqDense     *mat = (Mat_SeqDense*)A->data;
903:   PetscScalar      *av;
904:   PetscInt         i,j,idx=0;
905: #if defined(PETSC_HAVE_CUDA)
906:   PetscOffloadMask oldf;
907: #endif
908:   PetscErrorCode   ierr;

911:   MatDenseGetArray(A,&av);
912:   if (!mat->roworiented) {
913:     if (addv == INSERT_VALUES) {
914:       for (j=0; j<n; j++) {
915:         if (indexn[j] < 0) {idx += m; continue;}
916:         if (PetscUnlikelyDebug(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);
917:         for (i=0; i<m; i++) {
918:           if (indexm[i] < 0) {idx++; continue;}
919:           if (PetscUnlikelyDebug(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);
920:           av[indexn[j]*mat->lda + indexm[i]] = v[idx++];
921:         }
922:       }
923:     } else {
924:       for (j=0; j<n; j++) {
925:         if (indexn[j] < 0) {idx += m; continue;}
926:         if (PetscUnlikelyDebug(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);
927:         for (i=0; i<m; i++) {
928:           if (indexm[i] < 0) {idx++; continue;}
929:           if (PetscUnlikelyDebug(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);
930:           av[indexn[j]*mat->lda + indexm[i]] += v[idx++];
931:         }
932:       }
933:     }
934:   } else {
935:     if (addv == INSERT_VALUES) {
936:       for (i=0; i<m; i++) {
937:         if (indexm[i] < 0) { idx += n; continue;}
938:         if (PetscUnlikelyDebug(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);
939:         for (j=0; j<n; j++) {
940:           if (indexn[j] < 0) { idx++; continue;}
941:           if (PetscUnlikelyDebug(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);
942:           av[indexn[j]*mat->lda + indexm[i]] = v[idx++];
943:         }
944:       }
945:     } else {
946:       for (i=0; i<m; i++) {
947:         if (indexm[i] < 0) { idx += n; continue;}
948:         if (PetscUnlikelyDebug(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);
949:         for (j=0; j<n; j++) {
950:           if (indexn[j] < 0) { idx++; continue;}
951:           if (PetscUnlikelyDebug(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);
952:           av[indexn[j]*mat->lda + indexm[i]] += v[idx++];
953:         }
954:       }
955:     }
956:   }
957:   /* hack to prevent unneeded copy to the GPU while returning the array */
958: #if defined(PETSC_HAVE_CUDA)
959:   oldf = A->offloadmask;
960:   A->offloadmask = PETSC_OFFLOAD_GPU;
961: #endif
962:   MatDenseRestoreArray(A,&av);
963: #if defined(PETSC_HAVE_CUDA)
964:   A->offloadmask = (oldf == PETSC_OFFLOAD_UNALLOCATED ? PETSC_OFFLOAD_UNALLOCATED : PETSC_OFFLOAD_CPU);
965: #endif
966:   return(0);
967: }

969: static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
970: {
971:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
972:   const PetscScalar *vv;
973:   PetscInt          i,j;
974:   PetscErrorCode    ierr;

977:   MatDenseGetArrayRead(A,&vv);
978:   /* row-oriented output */
979:   for (i=0; i<m; i++) {
980:     if (indexm[i] < 0) {v += n;continue;}
981:     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);
982:     for (j=0; j<n; j++) {
983:       if (indexn[j] < 0) {v++; continue;}
984:       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);
985:       *v++ = vv[indexn[j]*mat->lda + indexm[i]];
986:     }
987:   }
988:   MatDenseRestoreArrayRead(A,&vv);
989:   return(0);
990: }

992: /* -----------------------------------------------------------------*/

994: PetscErrorCode MatView_Dense_Binary(Mat mat,PetscViewer viewer)
995: {
996:   PetscErrorCode    ierr;
997:   PetscBool         skipHeader;
998:   PetscViewerFormat format;
999:   PetscInt          header[4],M,N,m,lda,i,j,k;
1000:   const PetscScalar *v;
1001:   PetscScalar       *vwork;

1004:   PetscViewerSetUp(viewer);
1005:   PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
1006:   PetscViewerGetFormat(viewer,&format);
1007:   if (skipHeader) format = PETSC_VIEWER_NATIVE;

1009:   MatGetSize(mat,&M,&N);

1011:   /* write matrix header */
1012:   header[0] = MAT_FILE_CLASSID; header[1] = M; header[2] = N;
1013:   header[3] = (format == PETSC_VIEWER_NATIVE) ? MATRIX_BINARY_FORMAT_DENSE : M*N;
1014:   if (!skipHeader) {PetscViewerBinaryWrite(viewer,header,4,PETSC_INT);}

1016:   MatGetLocalSize(mat,&m,NULL);
1017:   if (format != PETSC_VIEWER_NATIVE) {
1018:     PetscInt nnz = m*N, *iwork;
1019:     /* store row lengths for each row */
1020:     PetscMalloc1(nnz,&iwork);
1021:     for (i=0; i<m; i++) iwork[i] = N;
1022:     PetscViewerBinaryWriteAll(viewer,iwork,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1023:     /* store column indices (zero start index) */
1024:     for (k=0, i=0; i<m; i++)
1025:       for (j=0; j<N; j++, k++)
1026:         iwork[k] = j;
1027:     PetscViewerBinaryWriteAll(viewer,iwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1028:     PetscFree(iwork);
1029:   }
1030:   /* store matrix values as a dense matrix in row major order */
1031:   PetscMalloc1(m*N,&vwork);
1032:   MatDenseGetArrayRead(mat,&v);
1033:   MatDenseGetLDA(mat,&lda);
1034:   for (k=0, i=0; i<m; i++)
1035:     for (j=0; j<N; j++, k++)
1036:       vwork[k] = v[i+lda*j];
1037:   MatDenseRestoreArrayRead(mat,&v);
1038:   PetscViewerBinaryWriteAll(viewer,vwork,m*N,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);
1039:   PetscFree(vwork);
1040:   return(0);
1041: }

1043: PetscErrorCode MatLoad_Dense_Binary(Mat mat,PetscViewer viewer)
1044: {
1046:   PetscBool      skipHeader;
1047:   PetscInt       header[4],M,N,m,nz,lda,i,j,k;
1048:   PetscInt       rows,cols;
1049:   PetscScalar    *v,*vwork;

1052:   PetscViewerSetUp(viewer);
1053:   PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);

1055:   if (!skipHeader) {
1056:     PetscViewerBinaryRead(viewer,header,4,NULL,PETSC_INT);
1057:     if (header[0] != MAT_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Not a matrix object in file");
1058:     M = header[1]; N = header[2];
1059:     if (M < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix row size (%D) in file is negative",M);
1060:     if (N < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix column size (%D) in file is negative",N);
1061:     nz = header[3];
1062:     if (nz != MATRIX_BINARY_FORMAT_DENSE && nz < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Unknown matrix format %D in file",nz);
1063:   } else {
1064:     MatGetSize(mat,&M,&N);
1065:     if (M < 0 || N < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Matrix binary file header was skipped, thus the user must specify the global sizes of input matrix");
1066:     nz = MATRIX_BINARY_FORMAT_DENSE;
1067:   }

1069:   /* setup global sizes if not set */
1070:   if (mat->rmap->N < 0) mat->rmap->N = M;
1071:   if (mat->cmap->N < 0) mat->cmap->N = N;
1072:   MatSetUp(mat);
1073:   /* check if global sizes are correct */
1074:   MatGetSize(mat,&rows,&cols);
1075:   if (M != rows || N != cols) SETERRQ4(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different sizes (%d, %d) than the input matrix (%d, %d)",M,N,rows,cols);

1077:   MatGetSize(mat,NULL,&N);
1078:   MatGetLocalSize(mat,&m,NULL);
1079:   MatDenseGetArray(mat,&v);
1080:   MatDenseGetLDA(mat,&lda);
1081:   if (nz == MATRIX_BINARY_FORMAT_DENSE) {  /* matrix in file is dense format */
1082:     PetscInt nnz = m*N;
1083:     /* read in matrix values */
1084:     PetscMalloc1(nnz,&vwork);
1085:     PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);
1086:     /* store values in column major order */
1087:     for (j=0; j<N; j++)
1088:       for (i=0; i<m; i++)
1089:         v[i+lda*j] = vwork[i*N+j];
1090:     PetscFree(vwork);
1091:   } else { /* matrix in file is sparse format */
1092:     PetscInt nnz = 0, *rlens, *icols;
1093:     /* read in row lengths */
1094:     PetscMalloc1(m,&rlens);
1095:     PetscViewerBinaryReadAll(viewer,rlens,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1096:     for (i=0; i<m; i++) nnz += rlens[i];
1097:     /* read in column indices and values */
1098:     PetscMalloc2(nnz,&icols,nnz,&vwork);
1099:     PetscViewerBinaryReadAll(viewer,icols,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1100:     PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);
1101:     /* store values in column major order */
1102:     for (k=0, i=0; i<m; i++)
1103:       for (j=0; j<rlens[i]; j++, k++)
1104:         v[i+lda*icols[k]] = vwork[k];
1105:     PetscFree(rlens);
1106:     PetscFree2(icols,vwork);
1107:   }
1108:   MatDenseRestoreArray(mat,&v);
1109:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
1110:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
1111:   return(0);
1112: }

1114: PetscErrorCode MatLoad_SeqDense(Mat newMat, PetscViewer viewer)
1115: {
1116:   PetscBool      isbinary, ishdf5;

1122:   /* force binary viewer to load .info file if it has not yet done so */
1123:   PetscViewerSetUp(viewer);
1124:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1125:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,  &ishdf5);
1126:   if (isbinary) {
1127:     MatLoad_Dense_Binary(newMat,viewer);
1128:   } else if (ishdf5) {
1129: #if defined(PETSC_HAVE_HDF5)
1130:     MatLoad_Dense_HDF5(newMat,viewer);
1131: #else
1132:     SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1133: #endif
1134:   } else {
1135:     SETERRQ2(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"Viewer type %s not yet supported for reading %s matrices",((PetscObject)viewer)->type_name,((PetscObject)newMat)->type_name);
1136:   }
1137:   return(0);
1138: }

1140: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
1141: {
1142:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
1143:   PetscErrorCode    ierr;
1144:   PetscInt          i,j;
1145:   const char        *name;
1146:   PetscScalar       *v,*av;
1147:   PetscViewerFormat format;
1148: #if defined(PETSC_USE_COMPLEX)
1149:   PetscBool         allreal = PETSC_TRUE;
1150: #endif

1153:   MatDenseGetArrayRead(A,(const PetscScalar**)&av);
1154:   PetscViewerGetFormat(viewer,&format);
1155:   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1156:     return(0);  /* do nothing for now */
1157:   } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1158:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1159:     for (i=0; i<A->rmap->n; i++) {
1160:       v    = av + i;
1161:       PetscViewerASCIIPrintf(viewer,"row %D:",i);
1162:       for (j=0; j<A->cmap->n; j++) {
1163: #if defined(PETSC_USE_COMPLEX)
1164:         if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
1165:           PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1166:         } else if (PetscRealPart(*v)) {
1167:           PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));
1168:         }
1169: #else
1170:         if (*v) {
1171:           PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);
1172:         }
1173: #endif
1174:         v += a->lda;
1175:       }
1176:       PetscViewerASCIIPrintf(viewer,"\n");
1177:     }
1178:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1179:   } else {
1180:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1181: #if defined(PETSC_USE_COMPLEX)
1182:     /* determine if matrix has all real values */
1183:     v = av;
1184:     for (i=0; i<A->rmap->n*A->cmap->n; i++) {
1185:       if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;}
1186:     }
1187: #endif
1188:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
1189:       PetscObjectGetName((PetscObject)A,&name);
1190:       PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);
1191:       PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);
1192:       PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
1193:     }

1195:     for (i=0; i<A->rmap->n; i++) {
1196:       v = av + i;
1197:       for (j=0; j<A->cmap->n; j++) {
1198: #if defined(PETSC_USE_COMPLEX)
1199:         if (allreal) {
1200:           PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));
1201:         } else {
1202:           PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1203:         }
1204: #else
1205:         PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);
1206: #endif
1207:         v += a->lda;
1208:       }
1209:       PetscViewerASCIIPrintf(viewer,"\n");
1210:     }
1211:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
1212:       PetscViewerASCIIPrintf(viewer,"];\n");
1213:     }
1214:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1215:   }
1216:   MatDenseRestoreArrayRead(A,(const PetscScalar**)&av);
1217:   PetscViewerFlush(viewer);
1218:   return(0);
1219: }

1221: static PetscErrorCode MatView_SeqDense_Binary(Mat mat,PetscViewer viewer)
1222: {
1223:   PetscErrorCode    ierr;
1224:   PetscViewerFormat format;
1225:   PetscInt          header[4],M,N,m,lda,i,j,k;
1226:   const PetscScalar *v;
1227:   PetscScalar       *vwork;

1230:   PetscViewerSetUp(viewer);

1232:   PetscViewerGetFormat(viewer,&format);
1233:   MatGetSize(mat,&M,&N);

1235:   /* write matrix header */
1236:   header[0] = MAT_FILE_CLASSID; header[1] = M; header[2] = N;
1237:   header[3] = (format == PETSC_VIEWER_NATIVE) ? MATRIX_BINARY_FORMAT_DENSE : M*N;
1238:   PetscViewerBinaryWrite(viewer,header,4,PETSC_INT);

1240:   MatGetLocalSize(mat,&m,NULL);
1241:   if (format != PETSC_VIEWER_NATIVE) {
1242:     PetscInt nnz = m*N, *iwork;
1243:     /* store row lengths for each row */
1244:     PetscMalloc1(nnz,&iwork);
1245:     for (i=0; i<m; i++) iwork[i] = N;
1246:     PetscViewerBinaryWrite(viewer,iwork,m,PETSC_INT);
1247:     /* store column indices (zero start index) */
1248:     for (k=0, i=0; i<m; i++)
1249:       for (j=0; j<N; j++, k++)
1250:         iwork[k] = j;
1251:     PetscViewerBinaryWrite(viewer,iwork,nnz,PETSC_INT);
1252:     PetscFree(iwork);
1253:   }
1254:   /* store the matrix values as a dense matrix in row major order */
1255:   PetscMalloc1(m*N,&vwork);
1256:   MatDenseGetArrayRead(mat,&v);
1257:   MatDenseGetLDA(mat,&lda);
1258:   for (k=0, i=0; i<m; i++)
1259:     for (j=0; j<N; j++, k++)
1260:       vwork[k] = v[i+lda*j];
1261:   MatDenseRestoreArrayRead(mat,&v);
1262:   PetscViewerBinaryWrite(viewer,vwork,m*N,PETSC_SCALAR);
1263:   PetscFree(vwork);
1264:   return(0);
1265: }

1267:  #include <petscdraw.h>
1268: static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa)
1269: {
1270:   Mat               A  = (Mat) Aa;
1271:   PetscErrorCode    ierr;
1272:   PetscInt          m  = A->rmap->n,n = A->cmap->n,i,j;
1273:   int               color = PETSC_DRAW_WHITE;
1274:   const PetscScalar *v;
1275:   PetscViewer       viewer;
1276:   PetscReal         xl,yl,xr,yr,x_l,x_r,y_l,y_r;
1277:   PetscViewerFormat format;

1280:   PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
1281:   PetscViewerGetFormat(viewer,&format);
1282:   PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);

1284:   /* Loop over matrix elements drawing boxes */
1285:   MatDenseGetArrayRead(A,&v);
1286:   if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1287:     PetscDrawCollectiveBegin(draw);
1288:     /* Blue for negative and Red for positive */
1289:     for (j = 0; j < n; j++) {
1290:       x_l = j; x_r = x_l + 1.0;
1291:       for (i = 0; i < m; i++) {
1292:         y_l = m - i - 1.0;
1293:         y_r = y_l + 1.0;
1294:         if (PetscRealPart(v[j*m+i]) >  0.) color = PETSC_DRAW_RED;
1295:         else if (PetscRealPart(v[j*m+i]) <  0.) color = PETSC_DRAW_BLUE;
1296:         else continue;
1297:         PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1298:       }
1299:     }
1300:     PetscDrawCollectiveEnd(draw);
1301:   } else {
1302:     /* use contour shading to indicate magnitude of values */
1303:     /* first determine max of all nonzero values */
1304:     PetscReal minv = 0.0, maxv = 0.0;
1305:     PetscDraw popup;

1307:     for (i=0; i < m*n; i++) {
1308:       if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1309:     }
1310:     if (minv >= maxv) maxv = minv + PETSC_SMALL;
1311:     PetscDrawGetPopup(draw,&popup);
1312:     PetscDrawScalePopup(popup,minv,maxv);

1314:     PetscDrawCollectiveBegin(draw);
1315:     for (j=0; j<n; j++) {
1316:       x_l = j;
1317:       x_r = x_l + 1.0;
1318:       for (i=0; i<m; i++) {
1319:         y_l   = m - i - 1.0;
1320:         y_r   = y_l + 1.0;
1321:         color = PetscDrawRealToColor(PetscAbsScalar(v[j*m+i]),minv,maxv);
1322:         PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1323:       }
1324:     }
1325:     PetscDrawCollectiveEnd(draw);
1326:   }
1327:   MatDenseRestoreArrayRead(A,&v);
1328:   return(0);
1329: }

1331: static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
1332: {
1333:   PetscDraw      draw;
1334:   PetscBool      isnull;
1335:   PetscReal      xr,yr,xl,yl,h,w;

1339:   PetscViewerDrawGetDraw(viewer,0,&draw);
1340:   PetscDrawIsNull(draw,&isnull);
1341:   if (isnull) return(0);

1343:   xr   = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
1344:   xr  += w;          yr += h;        xl = -w;     yl = -h;
1345:   PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
1346:   PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
1347:   PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
1348:   PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
1349:   PetscDrawSave(draw);
1350:   return(0);
1351: }

1353: PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
1354: {
1356:   PetscBool      iascii,isbinary,isdraw;

1359:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1360:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1361:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);

1363:   if (iascii) {
1364:     MatView_SeqDense_ASCII(A,viewer);
1365:   } else if (isbinary) {
1366:     MatView_SeqDense_Binary(A,viewer);
1367:   } else if (isdraw) {
1368:     MatView_SeqDense_Draw(A,viewer);
1369:   }
1370:   return(0);
1371: }

1373: static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A,const PetscScalar array[])
1374: {
1375:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

1378:   a->unplacedarray       = a->v;
1379:   a->unplaced_user_alloc = a->user_alloc;
1380:   a->v                   = (PetscScalar*) array;
1381: #if defined(PETSC_HAVE_CUDA)
1382:   A->offloadmask = PETSC_OFFLOAD_CPU;
1383: #endif
1384:   return(0);
1385: }

1387: static PetscErrorCode MatDenseResetArray_SeqDense(Mat A)
1388: {
1389:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

1392:   a->v             = a->unplacedarray;
1393:   a->user_alloc    = a->unplaced_user_alloc;
1394:   a->unplacedarray = NULL;
1395: #if defined(PETSC_HAVE_CUDA)
1396:   A->offloadmask = PETSC_OFFLOAD_CPU;
1397: #endif
1398:   return(0);
1399: }

1401: PetscErrorCode MatDestroy_SeqDense(Mat mat)
1402: {
1403:   Mat_SeqDense   *l = (Mat_SeqDense*)mat->data;

1407: #if defined(PETSC_USE_LOG)
1408:   PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n);
1409: #endif
1410:   PetscFree(l->pivots);
1411:   PetscFree(l->fwork);
1412:   MatDestroy(&l->ptapwork);
1413:   if (!l->user_alloc) {PetscFree(l->v);}
1414:   PetscFree(mat->data);

1416:   PetscObjectChangeTypeName((PetscObject)mat,0);
1417:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetLDA_C",NULL);
1418:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);
1419:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);
1420:   PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL);
1421:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL);
1422:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayRead_C",NULL);
1423:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayRead_C",NULL);
1424:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);
1425: #if defined(PETSC_HAVE_ELEMENTAL)
1426:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);
1427: #endif
1428: #if defined(PETSC_HAVE_CUDA)
1429:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqdensecuda_C",NULL);
1430:   PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",NULL);
1431:   PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdense_C",NULL);
1432:   PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqaijcusparse_seqdense_C",NULL);
1433: #endif
1434:   PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);
1435:   PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqaij_seqdense_C",NULL);
1436:   PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdense_seqdense_C",NULL);
1437:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);
1438:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);
1439:   PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqbaij_seqdense_C",NULL);
1440:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqbaij_seqdense_C",NULL);
1441:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqbaij_seqdense_C",NULL);
1442:   PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqsbaij_seqdense_C",NULL);
1443:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqsbaij_seqdense_C",NULL);
1444:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqsbaij_seqdense_C",NULL);
1445:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaijperm_seqdense_C",NULL);
1446:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaijperm_seqdense_C",NULL);
1447:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaijsell_seqdense_C",NULL);
1448:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaijsell_seqdense_C",NULL);
1449:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaijmkl_seqdense_C",NULL);
1450:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaijmkl_seqdense_C",NULL);
1451:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_nest_seqdense_C",NULL);
1452:   PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_nest_seqdense_C",NULL);

1454:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);
1455:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);
1456:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaijperm_seqdense_C",NULL);
1457:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaijperm_seqdense_C",NULL);
1458:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaijsell_seqdense_C",NULL);
1459:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaijsell_seqdense_C",NULL);
1460:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaijmkl_seqdense_C",NULL);
1461:   PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaijmkl_seqdense_C",NULL);
1462:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumn_C",NULL);
1463:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumn_C",NULL);
1464:   return(0);
1465: }

1467: static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
1468: {
1469:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1471:   PetscInt       k,j,m,n,M;
1472:   PetscScalar    *v,tmp;

1475:   m = A->rmap->n; M = mat->lda; n = A->cmap->n;
1476:   if (reuse == MAT_INPLACE_MATRIX && m == n) { /* in place transpose */
1477:     MatDenseGetArray(A,&v);
1478:     for (j=0; j<m; j++) {
1479:       for (k=0; k<j; k++) {
1480:         tmp        = v[j + k*M];
1481:         v[j + k*M] = v[k + j*M];
1482:         v[k + j*M] = tmp;
1483:       }
1484:     }
1485:     MatDenseRestoreArray(A,&v);
1486:   } else { /* out-of-place transpose */
1487:     Mat          tmat;
1488:     Mat_SeqDense *tmatd;
1489:     PetscScalar  *v2;
1490:     PetscInt     M2;

1492:     if (reuse != MAT_REUSE_MATRIX) {
1493:       MatCreate(PetscObjectComm((PetscObject)A),&tmat);
1494:       MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);
1495:       MatSetType(tmat,((PetscObject)A)->type_name);
1496:       MatSeqDenseSetPreallocation(tmat,NULL);
1497:     } else tmat = *matout;

1499:     MatDenseGetArrayRead(A,(const PetscScalar**)&v);
1500:     MatDenseGetArray(tmat,&v2);
1501:     tmatd = (Mat_SeqDense*)tmat->data;
1502:     M2    = tmatd->lda;
1503:     for (j=0; j<n; j++) {
1504:       for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M];
1505:     }
1506:     MatDenseRestoreArray(tmat,&v2);
1507:     MatDenseRestoreArrayRead(A,(const PetscScalar**)&v);
1508:     MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
1509:     MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
1510:     if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) *matout = tmat;
1511:     else {
1512:       MatHeaderMerge(A,&tmat);
1513:     }
1514:   }
1515:   return(0);
1516: }

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

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

1543: static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
1544: {
1545:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1546:   PetscInt          i,n,len;
1547:   PetscScalar       *x;
1548:   const PetscScalar *vv;
1549:   PetscErrorCode    ierr;

1552:   VecGetSize(v,&n);
1553:   VecGetArray(v,&x);
1554:   len  = PetscMin(A->rmap->n,A->cmap->n);
1555:   MatDenseGetArrayRead(A,&vv);
1556:   if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
1557:   for (i=0; i<len; i++) {
1558:     x[i] = vv[i*mat->lda + i];
1559:   }
1560:   MatDenseRestoreArrayRead(A,&vv);
1561:   VecRestoreArray(v,&x);
1562:   return(0);
1563: }

1565: static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
1566: {
1567:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1568:   const PetscScalar *l,*r;
1569:   PetscScalar       x,*v,*vv;
1570:   PetscErrorCode    ierr;
1571:   PetscInt          i,j,m = A->rmap->n,n = A->cmap->n;

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

1603: PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1604: {
1605:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1606:   PetscScalar       *v,*vv;
1607:   PetscReal         sum  = 0.0;
1608:   PetscInt          lda  =mat->lda,m=A->rmap->n,i,j;
1609:   PetscErrorCode    ierr;

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

1662: static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
1663: {
1664:   Mat_SeqDense   *aij = (Mat_SeqDense*)A->data;

1668:   switch (op) {
1669:   case MAT_ROW_ORIENTED:
1670:     aij->roworiented = flg;
1671:     break;
1672:   case MAT_NEW_NONZERO_LOCATIONS:
1673:   case MAT_NEW_NONZERO_LOCATION_ERR:
1674:   case MAT_NEW_NONZERO_ALLOCATION_ERR:
1675:   case MAT_NEW_DIAGONALS:
1676:   case MAT_KEEP_NONZERO_PATTERN:
1677:   case MAT_IGNORE_OFF_PROC_ENTRIES:
1678:   case MAT_USE_HASH_TABLE:
1679:   case MAT_IGNORE_ZERO_ENTRIES:
1680:   case MAT_IGNORE_LOWER_TRIANGULAR:
1681:   case MAT_SORTED_FULL:
1682:     PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
1683:     break;
1684:   case MAT_SPD:
1685:   case MAT_SYMMETRIC:
1686:   case MAT_STRUCTURALLY_SYMMETRIC:
1687:   case MAT_HERMITIAN:
1688:   case MAT_SYMMETRY_ETERNAL:
1689:     /* These options are handled directly by MatSetOption() */
1690:     break;
1691:   default:
1692:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]);
1693:   }
1694:   return(0);
1695: }

1697: static PetscErrorCode MatZeroEntries_SeqDense(Mat A)
1698: {
1699:   Mat_SeqDense   *l = (Mat_SeqDense*)A->data;
1701:   PetscInt       lda=l->lda,m=A->rmap->n,j;
1702:   PetscScalar    *v;

1705:   MatDenseGetArray(A,&v);
1706:   if (lda>m) {
1707:     for (j=0; j<A->cmap->n; j++) {
1708:       PetscArrayzero(v+j*lda,m);
1709:     }
1710:   } else {
1711:     PetscArrayzero(v,A->rmap->n*A->cmap->n);
1712:   }
1713:   MatDenseRestoreArray(A,&v);
1714:   return(0);
1715: }

1717: static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1718: {
1719:   PetscErrorCode    ierr;
1720:   Mat_SeqDense      *l = (Mat_SeqDense*)A->data;
1721:   PetscInt          m  = l->lda, n = A->cmap->n, i,j;
1722:   PetscScalar       *slot,*bb,*v;
1723:   const PetscScalar *xx;

1726:   if (PetscDefined(USE_DEBUG)) {
1727:     for (i=0; i<N; i++) {
1728:       if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
1729:       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);
1730:     }
1731:   }
1732:   if (!N) return(0);

1734:   /* fix right hand side if needed */
1735:   if (x && b) {
1736:     VecGetArrayRead(x,&xx);
1737:     VecGetArray(b,&bb);
1738:     for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
1739:     VecRestoreArrayRead(x,&xx);
1740:     VecRestoreArray(b,&bb);
1741:   }

1743:   MatDenseGetArray(A,&v);
1744:   for (i=0; i<N; i++) {
1745:     slot = v + rows[i];
1746:     for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
1747:   }
1748:   if (diag != 0.0) {
1749:     if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
1750:     for (i=0; i<N; i++) {
1751:       slot  = v + (m+1)*rows[i];
1752:       *slot = diag;
1753:     }
1754:   }
1755:   MatDenseRestoreArray(A,&v);
1756:   return(0);
1757: }

1759: static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A,PetscInt *lda)
1760: {
1761:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;

1764:   *lda = mat->lda;
1765:   return(0);
1766: }

1768: PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[])
1769: {
1770:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;

1773:   *array = mat->v;
1774:   return(0);
1775: }

1777: PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[])
1778: {
1780:   return(0);
1781: }

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

1786:    Logically Collective on Mat

1788:    Input Parameter:
1789: .  mat - a MATSEQDENSE or MATMPIDENSE matrix

1791:    Output Parameter:
1792: .   lda - the leading dimension

1794:    Level: intermediate

1796: .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatSeqDenseSetLDA()
1797: @*/
1798: PetscErrorCode  MatDenseGetLDA(Mat A,PetscInt *lda)
1799: {

1803:   PetscUseMethod(A,"MatDenseGetLDA_C",(Mat,PetscInt*),(A,lda));
1804:   return(0);
1805: }

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

1810:    Logically Collective on Mat

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

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

1818:    Level: intermediate

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

1827:   PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));
1828:   return(0);
1829: }

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

1834:    Logically Collective on Mat

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

1840:    Level: intermediate

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

1849:   PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));
1850:   if (array) *array = NULL;
1851:   PetscObjectStateIncrease((PetscObject)A);
1852:   return(0);
1853: }

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

1858:    Not Collective

1860:    Input Parameter:
1861: .  mat - a MATSEQDENSE or MATMPIDENSE matrix

1863:    Output Parameter:
1864: .   array - pointer to the data

1866:    Level: intermediate

1868: .seealso: MatDenseRestoreArray(), MatDenseGetArray(), MatDenseRestoreArrayRead()
1869: @*/
1870: PetscErrorCode  MatDenseGetArrayRead(Mat A,const PetscScalar **array)
1871: {

1875:   PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array));
1876:   return(0);
1877: }

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

1882:    Not Collective

1884:    Input Parameters:
1885: +  mat - a MATSEQDENSE or MATMPIDENSE matrix
1886: -  array - pointer to the data

1888:    Level: intermediate

1890: .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArray()
1891: @*/
1892: PetscErrorCode  MatDenseRestoreArrayRead(Mat A,const PetscScalar **array)
1893: {

1897:   PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array));
1898:   if (array) *array = NULL;
1899:   return(0);
1900: }

1902: static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B)
1903: {
1904:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1906:   PetscInt       i,j,nrows,ncols,blda;
1907:   const PetscInt *irow,*icol;
1908:   PetscScalar    *av,*bv,*v = mat->v;
1909:   Mat            newmat;

1912:   ISGetIndices(isrow,&irow);
1913:   ISGetIndices(iscol,&icol);
1914:   ISGetLocalSize(isrow,&nrows);
1915:   ISGetLocalSize(iscol,&ncols);

1917:   /* Check submatrixcall */
1918:   if (scall == MAT_REUSE_MATRIX) {
1919:     PetscInt n_cols,n_rows;
1920:     MatGetSize(*B,&n_rows,&n_cols);
1921:     if (n_rows != nrows || n_cols != ncols) {
1922:       /* resize the result matrix to match number of requested rows/columns */
1923:       MatSetSizes(*B,nrows,ncols,nrows,ncols);
1924:     }
1925:     newmat = *B;
1926:   } else {
1927:     /* Create and fill new matrix */
1928:     MatCreate(PetscObjectComm((PetscObject)A),&newmat);
1929:     MatSetSizes(newmat,nrows,ncols,nrows,ncols);
1930:     MatSetType(newmat,((PetscObject)A)->type_name);
1931:     MatSeqDenseSetPreallocation(newmat,NULL);
1932:   }

1934:   /* Now extract the data pointers and do the copy,column at a time */
1935:   MatDenseGetArray(newmat,&bv);
1936:   MatDenseGetLDA(newmat,&blda);
1937:   for (i=0; i<ncols; i++) {
1938:     av = v + mat->lda*icol[i];
1939:     for (j=0; j<nrows; j++) bv[j] = av[irow[j]];
1940:     bv += blda;
1941:   }
1942:   MatDenseRestoreArray(newmat,&bv);

1944:   /* Assemble the matrices so that the correct flags are set */
1945:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1946:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);

1948:   /* Free work space */
1949:   ISRestoreIndices(isrow,&irow);
1950:   ISRestoreIndices(iscol,&icol);
1951:   *B   = newmat;
1952:   return(0);
1953: }

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

1961:   if (scall == MAT_INITIAL_MATRIX) {
1962:     PetscCalloc1(n+1,B);
1963:   }

1965:   for (i=0; i<n; i++) {
1966:     MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
1967:   }
1968:   return(0);
1969: }

1971: static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
1972: {
1974:   return(0);
1975: }

1977: static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
1978: {
1980:   return(0);
1981: }

1983: static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
1984: {
1985:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data;
1986:   PetscErrorCode    ierr;
1987:   const PetscScalar *va;
1988:   PetscScalar       *vb;
1989:   PetscInt          lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j;

1992:   /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
1993:   if (A->ops->copy != B->ops->copy) {
1994:     MatCopy_Basic(A,B,str);
1995:     return(0);
1996:   }
1997:   if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)");
1998:   MatDenseGetArrayRead(A,&va);
1999:   MatDenseGetArray(B,&vb);
2000:   if (lda1>m || lda2>m) {
2001:     for (j=0; j<n; j++) {
2002:       PetscArraycpy(vb+j*lda2,va+j*lda1,m);
2003:     }
2004:   } else {
2005:     PetscArraycpy(vb,va,A->rmap->n*A->cmap->n);
2006:   }
2007:   MatDenseRestoreArray(B,&vb);
2008:   MatDenseRestoreArrayRead(A,&va);
2009:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
2010:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
2011:   return(0);
2012: }

2014: static PetscErrorCode MatSetUp_SeqDense(Mat A)
2015: {

2019:   PetscLayoutSetUp(A->rmap);
2020:   PetscLayoutSetUp(A->cmap);
2021:   if (!A->preallocated) {
2022:     MatSeqDenseSetPreallocation(A,0);
2023:   }
2024:   return(0);
2025: }

2027: static PetscErrorCode MatConjugate_SeqDense(Mat A)
2028: {
2029:   PetscInt       i,nz = A->rmap->n*A->cmap->n;
2030:   PetscScalar    *aa;

2034:   MatDenseGetArray(A,&aa);
2035:   for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
2036:   MatDenseRestoreArray(A,&aa);
2037:   return(0);
2038: }

2040: static PetscErrorCode MatRealPart_SeqDense(Mat A)
2041: {
2042:   PetscInt       i,nz = A->rmap->n*A->cmap->n;
2043:   PetscScalar    *aa;

2047:   MatDenseGetArray(A,&aa);
2048:   for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
2049:   MatDenseRestoreArray(A,&aa);
2050:   return(0);
2051: }

2053: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
2054: {
2055:   PetscInt       i,nz = A->rmap->n*A->cmap->n;
2056:   PetscScalar    *aa;

2060:   MatDenseGetArray(A,&aa);
2061:   for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
2062:   MatDenseRestoreArray(A,&aa);
2063:   return(0);
2064: }

2066: /* ----------------------------------------------------------------*/
2067: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2068: {
2070:   PetscInt       m=A->rmap->n,n=B->cmap->n;
2071:   PetscBool      flg;

2074:   MatSetSizes(C,m,n,m,n);
2075:   PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2076:   MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
2077:   MatSetUp(C);
2078:   return(0);
2079: }

2081: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2082: {
2083:   Mat_SeqDense       *a,*b=(Mat_SeqDense*)B->data,*c=(Mat_SeqDense*)C->data;
2084:   PetscBLASInt       m,n,k;
2085:   const PetscScalar *av,*bv;
2086:   PetscScalar       *cv;
2087:   PetscScalar       _DOne=1.0,_DZero=0.0;
2088:   PetscBool         flg;
2089:   PetscErrorCode    (*numeric)(Mat,Mat,Mat)=NULL;
2090:   PetscErrorCode    ierr;

2093:   /* Handle case where where user provided the final C matrix rather than calling MatMatMult() with MAT_INITIAL_MATRIX*/
2094:   PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
2095:   if (flg) numeric = MatMatMultNumeric_SeqAIJ_SeqDense;
2096:   PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&flg);
2097:   if (flg) numeric = MatMatMultNumeric_SeqBAIJ_SeqDense;
2098:   PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&flg);
2099:   if (flg) numeric = MatMatMultNumeric_SeqSBAIJ_SeqDense;
2100:   PetscObjectTypeCompare((PetscObject)A,MATNEST,&flg);
2101:   if (flg) numeric = MatMatMultNumeric_Nest_Dense;
2102:   if (numeric) {
2103:     C->ops->matmultnumeric = numeric;
2104:     (*numeric)(A,B,C);
2105:     return(0);
2106:   }
2107:   a = (Mat_SeqDense*)A->data;
2108:   PetscBLASIntCast(C->rmap->n,&m);
2109:   PetscBLASIntCast(C->cmap->n,&n);
2110:   PetscBLASIntCast(A->cmap->n,&k);
2111:   if (!m || !n || !k) return(0);
2112:   MatDenseGetArrayRead(A,&av);
2113:   MatDenseGetArrayRead(B,&bv);
2114:   MatDenseGetArray(C,&cv);
2115:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2116:   PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2117:   MatDenseRestoreArrayRead(A,&av);
2118:   MatDenseRestoreArrayRead(B,&bv);
2119:   MatDenseRestoreArray(C,&cv);
2120:   return(0);
2121: }

2123: PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2124: {
2126:   PetscInt       m=A->rmap->n,n=B->rmap->n;
2127:   PetscBool      flg;

2130:   MatSetSizes(C,m,n,m,n);
2131:   PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2132:   MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
2133:   MatSetUp(C);
2134:   return(0);
2135: }

2137: PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2138: {
2139:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2140:   Mat_SeqDense   *b = (Mat_SeqDense*)B->data;
2141:   Mat_SeqDense   *c = (Mat_SeqDense*)C->data;
2142:   PetscBLASInt   m,n,k;
2143:   PetscScalar    _DOne=1.0,_DZero=0.0;

2147:   PetscBLASIntCast(C->rmap->n,&m);
2148:   PetscBLASIntCast(C->cmap->n,&n);
2149:   PetscBLASIntCast(A->cmap->n,&k);
2150:   if (!m || !n || !k) return(0);
2151:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
2152:   PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2153:   return(0);
2154: }

2156: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2157: {
2159:   PetscInt       m=A->cmap->n,n=B->cmap->n;
2160:   PetscBool      flg;

2163:   MatSetSizes(C,m,n,m,n);
2164:   PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2165:   MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
2166:   MatSetUp(C);
2167:   return(0);
2168: }

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

2180:   PetscBLASIntCast(C->rmap->n,&m);
2181:   PetscBLASIntCast(C->cmap->n,&n);
2182:   PetscBLASIntCast(A->rmap->n,&k);
2183:   if (!m || !n || !k) return(0);
2184:   PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
2185:   PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2186:   return(0);
2187: }

2189: /* ----------------------------------------------- */
2190: static PetscErrorCode MatProductSetFromOptions_SeqDense_AB(Mat C)
2191: {
2193:   C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqDense;
2194:   C->ops->productsymbolic = MatProductSymbolic_AB;
2195:   /* dense mat may not call MatProductSymbolic(), thus set C->ops->productnumeric here */
2196:   C->ops->productnumeric  = MatProductNumeric_AB;
2197:   return(0);
2198: }

2200: static PetscErrorCode MatProductSetFromOptions_SeqDense_AtB(Mat C)
2201: {
2203:   C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_SeqDense_SeqDense;
2204:   C->ops->productsymbolic          = MatProductSymbolic_AtB;
2205:   C->ops->productnumeric           = MatProductNumeric_AtB;
2206:   return(0);
2207: }

2209: static PetscErrorCode MatProductSetFromOptions_SeqDense_ABt(Mat C)
2210: {
2212:   C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqDense_SeqDense;
2213:   C->ops->productsymbolic          = MatProductSymbolic_ABt;
2214:   C->ops->productnumeric           = MatProductNumeric_ABt;
2215:   return(0);
2216: }

2218: static PetscErrorCode MatProductSetFromOptions_SeqDense_PtAP(Mat C)
2219: {
2221:   C->ops->productsymbolic = MatProductSymbolic_Basic;
2222:   return(0);
2223: }

2225: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense(Mat C)
2226: {
2228:   Mat_Product    *product = C->product;

2231:   switch (product->type) {
2232:   case MATPRODUCT_AB:
2233:     MatProductSetFromOptions_SeqDense_AB(C);
2234:     break;
2235:   case MATPRODUCT_AtB:
2236:     MatProductSetFromOptions_SeqDense_AtB(C);
2237:     break;
2238:   case MATPRODUCT_ABt:
2239:     MatProductSetFromOptions_SeqDense_ABt(C);
2240:     break;
2241:   case MATPRODUCT_PtAP:
2242:     MatProductSetFromOptions_SeqDense_PtAP(C);
2243:     break;
2244:   default: SETERRQ1(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"MatProduct type %s is not supported for SeqDense and SeqDense matrices",MatProductTypes[product->type]);
2245:   }
2246:   return(0);
2247: }
2248: /* ----------------------------------------------- */

2250: static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[])
2251: {
2252:   Mat_SeqDense       *a = (Mat_SeqDense*)A->data;
2253:   PetscErrorCode     ierr;
2254:   PetscInt           i,j,m = A->rmap->n,n = A->cmap->n,p;
2255:   PetscScalar        *x;
2256:   const PetscScalar *aa;

2259:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2260:   VecGetArray(v,&x);
2261:   VecGetLocalSize(v,&p);
2262:   MatDenseGetArrayRead(A,&aa);
2263:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2264:   for (i=0; i<m; i++) {
2265:     x[i] = aa[i]; if (idx) idx[i] = 0;
2266:     for (j=1; j<n; j++) {
2267:       if (PetscRealPart(x[i]) < PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;}
2268:     }
2269:   }
2270:   MatDenseRestoreArrayRead(A,&aa);
2271:   VecRestoreArray(v,&x);
2272:   return(0);
2273: }

2275: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])
2276: {
2277:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
2278:   PetscErrorCode    ierr;
2279:   PetscInt          i,j,m = A->rmap->n,n = A->cmap->n,p;
2280:   PetscScalar       *x;
2281:   PetscReal         atmp;
2282:   const PetscScalar *aa;

2285:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2286:   VecGetArray(v,&x);
2287:   VecGetLocalSize(v,&p);
2288:   MatDenseGetArrayRead(A,&aa);
2289:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2290:   for (i=0; i<m; i++) {
2291:     x[i] = PetscAbsScalar(aa[i]);
2292:     for (j=1; j<n; j++) {
2293:       atmp = PetscAbsScalar(aa[i+a->lda*j]);
2294:       if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;}
2295:     }
2296:   }
2297:   MatDenseRestoreArrayRead(A,&aa);
2298:   VecRestoreArray(v,&x);
2299:   return(0);
2300: }

2302: static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])
2303: {
2304:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
2305:   PetscErrorCode    ierr;
2306:   PetscInt          i,j,m = A->rmap->n,n = A->cmap->n,p;
2307:   PetscScalar       *x;
2308:   const PetscScalar *aa;

2311:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2312:   MatDenseGetArrayRead(A,&aa);
2313:   VecGetArray(v,&x);
2314:   VecGetLocalSize(v,&p);
2315:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2316:   for (i=0; i<m; i++) {
2317:     x[i] = aa[i]; if (idx) idx[i] = 0;
2318:     for (j=1; j<n; j++) {
2319:       if (PetscRealPart(x[i]) > PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;}
2320:     }
2321:   }
2322:   VecRestoreArray(v,&x);
2323:   MatDenseRestoreArrayRead(A,&aa);
2324:   return(0);
2325: }

2327: static PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
2328: {
2329:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
2330:   PetscErrorCode    ierr;
2331:   PetscScalar       *x;
2332:   const PetscScalar *aa;

2335:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2336:   MatDenseGetArrayRead(A,&aa);
2337:   VecGetArray(v,&x);
2338:   PetscArraycpy(x,aa+col*a->lda,A->rmap->n);
2339:   VecRestoreArray(v,&x);
2340:   MatDenseRestoreArrayRead(A,&aa);
2341:   return(0);
2342: }

2344: PETSC_INTERN PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms)
2345: {
2346:   PetscErrorCode    ierr;
2347:   PetscInt          i,j,m,n;
2348:   const PetscScalar *a;

2351:   MatGetSize(A,&m,&n);
2352:   PetscArrayzero(norms,n);
2353:   MatDenseGetArrayRead(A,&a);
2354:   if (type == NORM_2) {
2355:     for (i=0; i<n; i++) {
2356:       for (j=0; j<m; j++) {
2357:         norms[i] += PetscAbsScalar(a[j]*a[j]);
2358:       }
2359:       a += m;
2360:     }
2361:   } else if (type == NORM_1) {
2362:     for (i=0; i<n; i++) {
2363:       for (j=0; j<m; j++) {
2364:         norms[i] += PetscAbsScalar(a[j]);
2365:       }
2366:       a += m;
2367:     }
2368:   } else if (type == NORM_INFINITY) {
2369:     for (i=0; i<n; i++) {
2370:       for (j=0; j<m; j++) {
2371:         norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]);
2372:       }
2373:       a += m;
2374:     }
2375:   } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType");
2376:   MatDenseRestoreArrayRead(A,&a);
2377:   if (type == NORM_2) {
2378:     for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]);
2379:   }
2380:   return(0);
2381: }

2383: static PetscErrorCode  MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
2384: {
2386:   PetscScalar    *a;
2387:   PetscInt       m,n,i;

2390:   MatGetSize(x,&m,&n);
2391:   MatDenseGetArray(x,&a);
2392:   for (i=0; i<m*n; i++) {
2393:     PetscRandomGetValue(rctx,a+i);
2394:   }
2395:   MatDenseRestoreArray(x,&a);
2396:   return(0);
2397: }

2399: static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool  *missing,PetscInt *d)
2400: {
2402:   *missing = PETSC_FALSE;
2403:   return(0);
2404: }

2406: /* vals is not const */
2407: static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals)
2408: {
2410:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2411:   PetscScalar    *v;

2414:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2415:   MatDenseGetArray(A,&v);
2416:   *vals = v+col*a->lda;
2417:   MatDenseRestoreArray(A,&v);
2418:   return(0);
2419: }

2421: static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals)
2422: {
2424:   *vals = 0; /* user cannot accidently use the array later */
2425:   return(0);
2426: }

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

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

2584:    Collective

2586:    Input Parameters:
2587: +  comm - MPI communicator, set to PETSC_COMM_SELF
2588: .  m - number of rows
2589: .  n - number of columns
2590: -  data - optional location of matrix data in column major order.  Set data=NULL for PETSc
2591:    to control all matrix memory allocation.

2593:    Output Parameter:
2594: .  A - the matrix

2596:    Notes:
2597:    The data input variable is intended primarily for Fortran programmers
2598:    who wish to allocate their own matrix memory space.  Most users should
2599:    set data=NULL.

2601:    Level: intermediate

2603: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
2604: @*/
2605: PetscErrorCode  MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
2606: {

2610:   MatCreate(comm,A);
2611:   MatSetSizes(*A,m,n,m,n);
2612:   MatSetType(*A,MATSEQDENSE);
2613:   MatSeqDenseSetPreallocation(*A,data);
2614:   return(0);
2615: }

2617: /*@C
2618:    MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements

2620:    Collective

2622:    Input Parameters:
2623: +  B - the matrix
2624: -  data - the array (or NULL)

2626:    Notes:
2627:    The data input variable is intended primarily for Fortran programmers
2628:    who wish to allocate their own matrix memory space.  Most users should
2629:    need not call this routine.

2631:    Level: intermediate

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

2635: @*/
2636: PetscErrorCode  MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
2637: {

2641:   PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
2642:   return(0);
2643: }

2645: PetscErrorCode  MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
2646: {
2647:   Mat_SeqDense   *b;

2651:   B->preallocated = PETSC_TRUE;

2653:   PetscLayoutSetUp(B->rmap);
2654:   PetscLayoutSetUp(B->cmap);

2656:   b       = (Mat_SeqDense*)B->data;
2657:   b->Mmax = B->rmap->n;
2658:   b->Nmax = B->cmap->n;
2659:   if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n;

2661:   PetscIntMultError(b->lda,b->Nmax,NULL);
2662:   if (!data) { /* petsc-allocated storage */
2663:     if (!b->user_alloc) { PetscFree(b->v); }
2664:     PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);
2665:     PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));

2667:     b->user_alloc = PETSC_FALSE;
2668:   } else { /* user-allocated storage */
2669:     if (!b->user_alloc) { PetscFree(b->v); }
2670:     b->v          = data;
2671:     b->user_alloc = PETSC_TRUE;
2672:   }
2673:   B->assembled = PETSC_TRUE;
2674:   return(0);
2675: }

2677: #if defined(PETSC_HAVE_ELEMENTAL)
2678: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
2679: {
2680:   Mat               mat_elemental;
2681:   PetscErrorCode    ierr;
2682:   const PetscScalar *array;
2683:   PetscScalar       *v_colwise;
2684:   PetscInt          M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols;

2687:   PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);
2688:   MatDenseGetArrayRead(A,&array);
2689:   /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
2690:   k = 0;
2691:   for (j=0; j<N; j++) {
2692:     cols[j] = j;
2693:     for (i=0; i<M; i++) {
2694:       v_colwise[j*M+i] = array[k++];
2695:     }
2696:   }
2697:   for (i=0; i<M; i++) {
2698:     rows[i] = i;
2699:   }
2700:   MatDenseRestoreArrayRead(A,&array);

2702:   MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
2703:   MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);
2704:   MatSetType(mat_elemental,MATELEMENTAL);
2705:   MatSetUp(mat_elemental);

2707:   /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
2708:   MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);
2709:   MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
2710:   MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);
2711:   PetscFree3(v_colwise,rows,cols);

2713:   if (reuse == MAT_INPLACE_MATRIX) {
2714:     MatHeaderReplace(A,&mat_elemental);
2715:   } else {
2716:     *newmat = mat_elemental;
2717:   }
2718:   return(0);
2719: }
2720: #endif

2722: /*@C
2723:   MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array

2725:   Input parameter:
2726: + A - the matrix
2727: - lda - the leading dimension

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

2734:   Level: intermediate

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

2738: @*/
2739: PetscErrorCode  MatSeqDenseSetLDA(Mat B,PetscInt lda)
2740: {
2741:   Mat_SeqDense *b = (Mat_SeqDense*)B->data;

2744:   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);
2745:   b->lda       = lda;
2746:   b->changelda = PETSC_FALSE;
2747:   b->Mmax      = PetscMax(b->Mmax,lda);
2748:   return(0);
2749: }

2751: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat)
2752: {
2754:   PetscMPIInt    size;

2757:   MPI_Comm_size(comm,&size);
2758:   if (size == 1) {
2759:     if (scall == MAT_INITIAL_MATRIX) {
2760:       MatDuplicate(inmat,MAT_COPY_VALUES,outmat);
2761:     } else {
2762:       MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);
2763:     }
2764:   } else {
2765:     MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);
2766:   }
2767:   return(0);
2768: }

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

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

2776:   Level: beginner

2778: .seealso: MatCreateSeqDense()

2780: M*/
2781: PetscErrorCode MatCreate_SeqDense(Mat B)
2782: {
2783:   Mat_SeqDense   *b;
2785:   PetscMPIInt    size;

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

2791:   PetscNewLog(B,&b);
2792:   PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
2793:   B->data = (void*)b;

2795:   b->roworiented = PETSC_TRUE;

2797:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetLDA_C",MatDenseGetLDA_SeqDense);
2798:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);
2799:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);
2800:   PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);
2801:   PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);
2802:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);
2803:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense);
2804:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);
2805: #if defined(PETSC_HAVE_ELEMENTAL)
2806:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);
2807: #endif
2808: #if defined(PETSC_HAVE_CUDA)
2809:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqdensecuda_C",MatConvert_SeqDense_SeqDenseCUDA);
2810:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",MatProductSetFromOptions_SeqDense);
2811:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdense_C",MatProductSetFromOptions_SeqDense);
2812:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaijcusparse_seqdense_C",MatProductSetFromOptions_SeqAIJ_SeqDense);
2813: #endif
2814:   PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);
2815:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqdense_C",MatProductSetFromOptions_SeqAIJ_SeqDense);
2816:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdense_C",MatProductSetFromOptions_SeqDense);
2817:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2818:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2819:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);
2820:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqbaij_seqdense_C",MatMatMultSymbolic_SeqBAIJ_SeqDense);
2821:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqbaij_seqdense_C",MatMatMultNumeric_SeqBAIJ_SeqDense);
2822:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqsbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);
2823:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqsbaij_seqdense_C",MatMatMultSymbolic_SeqSBAIJ_SeqDense);
2824:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqsbaij_seqdense_C",MatMatMultNumeric_SeqSBAIJ_SeqDense);
2825:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijperm_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2826:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijperm_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2827:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijsell_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2828:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijsell_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2829:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijmkl_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2830:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijmkl_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2831:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_nest_seqdense_C",MatMatMultSymbolic_Nest_Dense);
2832:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_nest_seqdense_C",MatMatMultNumeric_Nest_Dense);

2834:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2835:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2836:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijperm_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2837:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijperm_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2838:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijsell_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2839:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijsell_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);

2841:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijmkl_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2842:   PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijmkl_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2843:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);
2844:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);
2845:   PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
2846:   return(0);
2847: }

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

2852:    Not Collective

2854:    Input Parameter:
2855: +  mat - a MATSEQDENSE or MATMPIDENSE matrix
2856: -  col - column index

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

2861:    Level: intermediate

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

2870:   PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));
2871:   return(0);
2872: }

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

2877:    Not Collective

2879:    Input Parameter:
2880: .  mat - a MATSEQDENSE or MATMPIDENSE matrix

2882:    Output Parameter:
2883: .  vals - pointer to the data

2885:    Level: intermediate

2887: .seealso: MatDenseGetColumn()
2888: @*/
2889: PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals)
2890: {

2894:   PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));
2895:   return(0);
2896: }