Actual source code: dense.c


  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      cisdense;

171:   MatSetSizes(C,P->cmap->n,P->cmap->n,P->cmap->N,P->cmap->N);
172:   PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
173:   if (!cisdense) {
174:     PetscBool flg;

176:     PetscObjectTypeCompare((PetscObject)P,((PetscObject)A)->type_name,&flg);
177:     MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
178:   }
179:   MatSetUp(C);
180:   c    = (Mat_SeqDense*)C->data;
181:   MatCreate(PetscObjectComm((PetscObject)A),&c->ptapwork);
182:   MatSetSizes(c->ptapwork,A->rmap->n,P->cmap->n,A->rmap->N,P->cmap->N);
183:   MatSetType(c->ptapwork,((PetscObject)C)->type_name);
184:   MatSetUp(c->ptapwork);
185:   return(0);
186: }

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

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

225:   if (reuse == MAT_INPLACE_MATRIX) {
226:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
227:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
228:     MatHeaderReplace(A,&B);
229:   } else {
230:     if (B) *newmat = B;
231:     MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
232:     MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
233:   }
234:   return(0);
235: }

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

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

269:   if (reuse == MAT_INPLACE_MATRIX) {
270:     MatHeaderReplace(A,&B);
271:   } else if (reuse != MAT_REUSE_MATRIX) *newmat = B;
272:   return(0);
273: }

275: PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str)
276: {
277:   Mat_SeqDense      *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data;
278:   const PetscScalar *xv;
279:   PetscScalar       *yv;
280:   PetscBLASInt      N,m,ldax = 0,lday = 0,one = 1;
281:   PetscErrorCode    ierr;

284:   MatDenseGetArrayRead(X,&xv);
285:   MatDenseGetArray(Y,&yv);
286:   PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);
287:   PetscBLASIntCast(X->rmap->n,&m);
288:   PetscBLASIntCast(x->lda,&ldax);
289:   PetscBLASIntCast(y->lda,&lday);
290:   if (ldax>m || lday>m) {
291:     PetscInt j;

293:     for (j=0; j<X->cmap->n; j++) {
294:       PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&alpha,xv+j*ldax,&one,yv+j*lday,&one));
295:     }
296:   } else {
297:     PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&alpha,xv,&one,yv,&one));
298:   }
299:   MatDenseRestoreArrayRead(X,&xv);
300:   MatDenseRestoreArray(Y,&yv);
301:   PetscLogFlops(PetscMax(2.0*N-1,0));
302:   return(0);
303: }

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

310:   info->block_size        = 1.0;
311:   info->nz_allocated      = N;
312:   info->nz_used           = N;
313:   info->nz_unneeded       = 0;
314:   info->assemblies        = A->num_ass;
315:   info->mallocs           = 0;
316:   info->memory            = ((PetscObject)A)->mem;
317:   info->fill_ratio_given  = 0;
318:   info->fill_ratio_needed = 0;
319:   info->factor_mallocs    = 0;
320:   return(0);
321: }

323: PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha)
324: {
325:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
326:   PetscScalar    *v;
328:   PetscBLASInt   one = 1,j,nz,lda = 0;

331:   MatDenseGetArray(A,&v);
332:   PetscBLASIntCast(a->lda,&lda);
333:   if (lda>A->rmap->n) {
334:     PetscBLASIntCast(A->rmap->n,&nz);
335:     for (j=0; j<A->cmap->n; j++) {
336:       PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&alpha,v+j*lda,&one));
337:     }
338:   } else {
339:     PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);
340:     PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&alpha,v,&one));
341:   }
342:   PetscLogFlops(nz);
343:   MatDenseRestoreArray(A,&v);
344:   return(0);
345: }

347: static PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool  *fl)
348: {
349:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
350:   PetscInt          i,j,m = A->rmap->n,N = a->lda;
351:   const PetscScalar *v;
352:   PetscErrorCode    ierr;

355:   *fl = PETSC_FALSE;
356:   if (A->rmap->n != A->cmap->n) return(0);
357:   MatDenseGetArrayRead(A,&v);
358:   for (i=0; i<m; i++) {
359:     for (j=i; j<m; j++) {
360:       if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) {
361:         goto restore;
362:       }
363:     }
364:   }
365:   *fl  = PETSC_TRUE;
366: restore:
367:   MatDenseRestoreArrayRead(A,&v);
368:   return(0);
369: }

371: static PetscErrorCode MatIsSymmetric_SeqDense(Mat A,PetscReal rtol,PetscBool  *fl)
372: {
373:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
374:   PetscInt          i,j,m = A->rmap->n,N = a->lda;
375:   const PetscScalar *v;
376:   PetscErrorCode    ierr;

379:   *fl = PETSC_FALSE;
380:   if (A->rmap->n != A->cmap->n) return(0);
381:   MatDenseGetArrayRead(A,&v);
382:   for (i=0; i<m; i++) {
383:     for (j=i; j<m; j++) {
384:       if (PetscAbsScalar(v[i+j*N] - v[j+i*N]) > rtol) {
385:         goto restore;
386:       }
387:     }
388:   }
389:   *fl  = PETSC_TRUE;
390: restore:
391:   MatDenseRestoreArrayRead(A,&v);
392:   return(0);
393: }

395: PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues)
396: {
397:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
399:   PetscInt       lda = (PetscInt)mat->lda,j,m,nlda = lda;

402:   PetscLayoutReference(A->rmap,&newi->rmap);
403:   PetscLayoutReference(A->cmap,&newi->cmap);
404:   if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { /* propagate LDA */
405:     MatDenseSetLDA(newi,lda);
406:   }
407:   MatSeqDenseSetPreallocation(newi,NULL);
408:   if (cpvalues == MAT_COPY_VALUES) {
409:     const PetscScalar *av;
410:     PetscScalar       *v;

412:     MatDenseGetArrayRead(A,&av);
413:     MatDenseGetArray(newi,&v);
414:     MatDenseGetLDA(newi,&nlda);
415:     m    = A->rmap->n;
416:     if (lda>m || nlda>m) {
417:       for (j=0; j<A->cmap->n; j++) {
418:         PetscArraycpy(v+j*nlda,av+j*lda,m);
419:       }
420:     } else {
421:       PetscArraycpy(v,av,A->rmap->n*A->cmap->n);
422:     }
423:     MatDenseRestoreArray(newi,&v);
424:     MatDenseRestoreArrayRead(A,&av);
425:   }
426:   return(0);
427: }

429: PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat)
430: {

434:   MatCreate(PetscObjectComm((PetscObject)A),newmat);
435:   MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
436:   MatSetType(*newmat,((PetscObject)A)->type_name);
437:   MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);
438:   return(0);
439: }

441: static PetscErrorCode MatSolve_SeqDense_Internal_LU(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k, PetscBool T)
442: {
443:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
444:   PetscBLASInt    info;
445:   PetscErrorCode  ierr;

448:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
449:   PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_(T ? "T" : "N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
450:   PetscFPTrapPop();
451:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
452:   PetscLogFlops(nrhs*(2.0*m*m - m));
453:   return(0);
454: }

456: static PetscErrorCode MatConjugate_SeqDense(Mat);

458: static PetscErrorCode MatSolve_SeqDense_Internal_Cholesky(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k, PetscBool T)
459: {
460:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
461:   PetscBLASInt    info;
462:   PetscErrorCode  ierr;

465:   if (A->spd) {
466:     if (PetscDefined(USE_COMPLEX) && T) {MatConjugate_SeqDense(A);}
467:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
468:     PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info));
469:     PetscFPTrapPop();
470:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
471:     if (PetscDefined(USE_COMPLEX) && T) {MatConjugate_SeqDense(A);}
472: #if defined(PETSC_USE_COMPLEX)
473:   } else if (A->hermitian) {
474:     if (T) {MatConjugate_SeqDense(A);}
475:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
476:     PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
477:     PetscFPTrapPop();
478:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve");
479:     if (T) {MatConjugate_SeqDense(A);}
480: #endif
481:   } else { /* symmetric case */
482:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
483:     PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
484:     PetscFPTrapPop();
485:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
486:   }
487:   PetscLogFlops(nrhs*(2.0*m*m - m));
488:   return(0);
489: }

491: static PetscErrorCode MatSolve_SeqDense_Internal_QR(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k)
492: {
493:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
494:   PetscBLASInt    info;
495:   char            trans;
496:   PetscErrorCode  ierr;

499:   if (PetscDefined(USE_COMPLEX)) {
500:     trans = 'C';
501:   } else {
502:     trans = 'T';
503:   }
504:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
505:   PetscStackCallBLAS("LAPACKormqr",LAPACKormqr_("L", &trans, &m,&nrhs,&mat->rank,mat->v,&mat->lda,mat->tau,x,&ldx,mat->fwork,&mat->lfwork,&info));
506:   PetscFPTrapPop();
507:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"ORMQR - Bad orthogonal transform");
508:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
509:   PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U", "N", "N", &mat->rank,&nrhs,mat->v,&mat->lda,x,&ldx,&info));
510:   PetscFPTrapPop();
511:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"TRTRS - Bad triangular solve");
512:   for (PetscInt j = 0; j < nrhs; j++) {
513:     for (PetscInt i = mat->rank; i < k; i++) {
514:       x[j*ldx + i] = 0.;
515:     }
516:   }
517:   PetscLogFlops(nrhs*(4.0*m*mat->rank - PetscSqr(mat->rank)));
518:   return(0);
519: }

521: static PetscErrorCode MatSolveTranspose_SeqDense_Internal_QR(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k)
522: {
523:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
524:   PetscBLASInt      info;
525:   PetscErrorCode    ierr;

528:   if (A->rmap->n == A->cmap->n && mat->rank == A->rmap->n) {
529:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
530:     PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U", "T", "N", &m,&nrhs,mat->v,&mat->lda,x,&ldx,&info));
531:     PetscFPTrapPop();
532:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"TRTRS - Bad triangular solve");
533:     if (PetscDefined(USE_COMPLEX)) {MatConjugate_SeqDense(A);}
534:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
535:     PetscStackCallBLAS("LAPACKormqr",LAPACKormqr_("L", "N", &m,&nrhs,&mat->rank,mat->v,&mat->lda,mat->tau,x,&ldx,mat->fwork,&mat->lfwork,&info));
536:     PetscFPTrapPop();
537:     if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"ORMQR - Bad orthogonal transform");
538:     if (PetscDefined(USE_COMPLEX)) {MatConjugate_SeqDense(A);}
539:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"QR factored matrix cannot be used for transpose solve");
540:   PetscLogFlops(nrhs*(4.0*m*mat->rank - PetscSqr(mat->rank)));
541:   return(0);
542: }

544: static PetscErrorCode MatSolve_SeqDense_SetUp(Mat A, Vec xx, Vec yy, PetscScalar **_y, PetscBLASInt *_m, PetscBLASInt *_k)
545: {
546:   Mat_SeqDense      *mat = (Mat_SeqDense *) A->data;
547:   PetscScalar       *y;
548:   PetscBLASInt      m=0, k=0;
549:   PetscErrorCode    ierr;

552:   PetscBLASIntCast(A->rmap->n,&m);
553:   PetscBLASIntCast(A->cmap->n,&k);
554:   if (k < m) {
555:     VecCopy(xx, mat->qrrhs);
556:     VecGetArray(mat->qrrhs,&y);
557:   } else {
558:     VecCopy(xx, yy);
559:     VecGetArray(yy,&y);
560:   }
561:   *_y = y;
562:   *_k = k;
563:   *_m = m;
564:   return(0);
565: }

567: static PetscErrorCode MatSolve_SeqDense_TearDown(Mat A, Vec xx, Vec yy, PetscScalar **_y, PetscBLASInt *_m, PetscBLASInt *_k)
568: {
569:   Mat_SeqDense   *mat = (Mat_SeqDense *) A->data;
570:   PetscScalar    *y = NULL;
571:   PetscBLASInt   m, k;

575:   y   = *_y;
576:   *_y = NULL;
577:   k   = *_k;
578:   m   = *_m;
579:   if (k < m) {
580:     PetscScalar *yv;
581:     VecGetArray(yy,&yv);
582:     PetscArraycpy(yv, y, k);
583:     VecRestoreArray(yy,&yv);
584:     VecRestoreArray(mat->qrrhs, &y);
585:   } else {
586:     VecRestoreArray(yy,&y);
587:   }
588:   return(0);
589: }

591: static PetscErrorCode MatSolve_SeqDense_LU(Mat A, Vec xx, Vec yy)
592: {
593:   PetscScalar    *y = NULL;
594:   PetscBLASInt   m = 0, k = 0;

598:   MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);
599:   MatSolve_SeqDense_Internal_LU(A, y, m, m, 1, k, PETSC_FALSE);
600:   MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);
601:   return(0);
602: }

604: static PetscErrorCode MatSolveTranspose_SeqDense_LU(Mat A, Vec xx, Vec yy)
605: {
606:   PetscScalar    *y = NULL;
607:   PetscBLASInt   m = 0, k = 0;

611:   MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);
612:   MatSolve_SeqDense_Internal_LU(A, y, m, m, 1, k, PETSC_TRUE);
613:   MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);
614:   return(0);
615: }

617: static PetscErrorCode MatSolve_SeqDense_Cholesky(Mat A, Vec xx, Vec yy)
618: {
619:   PetscScalar    *y = NULL;
620:   PetscBLASInt   m = 0, k = 0;

624:   MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);
625:   MatSolve_SeqDense_Internal_Cholesky(A, y, m, m, 1, k, PETSC_FALSE);
626:   MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);
627:   return(0);
628: }

630: static PetscErrorCode MatSolveTranspose_SeqDense_Cholesky(Mat A, Vec xx, Vec yy)
631: {
632:   PetscScalar    *y = NULL;
633:   PetscBLASInt   m = 0, k = 0;

637:   MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);
638:   MatSolve_SeqDense_Internal_Cholesky(A, y, m, m, 1, k, PETSC_TRUE);
639:   MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);
640:   return(0);
641: }

643: static PetscErrorCode MatSolve_SeqDense_QR(Mat A, Vec xx, Vec yy)
644: {
645:   PetscScalar    *y = NULL;
646:   PetscBLASInt   m = 0, k = 0;

650:   MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);
651:   MatSolve_SeqDense_Internal_QR(A, y, PetscMax(m,k), m, 1, k);
652:   MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);
653:   return(0);
654: }

656: static PetscErrorCode MatSolveTranspose_SeqDense_QR(Mat A, Vec xx, Vec yy)
657: {
658:   PetscScalar    *y = NULL;
659:   PetscBLASInt   m = 0, k = 0;

663:   MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);
664:   MatSolveTranspose_SeqDense_Internal_QR(A, y, PetscMax(m,k), m, 1, k);
665:   MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);
666:   return(0);
667: }

669: static PetscErrorCode MatMatSolve_SeqDense_SetUp(Mat A, Mat B, Mat X, PetscScalar **_y, PetscBLASInt *_ldy, PetscBLASInt *_m, PetscBLASInt *_nrhs, PetscBLASInt *_k)
670: {
671:   PetscErrorCode    ierr;
672:   const PetscScalar *b;
673:   PetscScalar       *y;
674:   PetscInt          n, _ldb, _ldx;
675:   PetscBLASInt      nrhs=0,m=0,k=0,ldb=0,ldx=0,ldy=0;

678:   *_ldy=0; *_m=0; *_nrhs=0; *_k=0; *_y = NULL;
679:   PetscBLASIntCast(A->rmap->n,&m);
680:   PetscBLASIntCast(A->cmap->n,&k);
681:   MatGetSize(B,NULL,&n);
682:   PetscBLASIntCast(n,&nrhs);
683:   MatDenseGetLDA(B,&_ldb);
684:   PetscBLASIntCast(_ldb, &ldb);
685:   MatDenseGetLDA(X,&_ldx);
686:   PetscBLASIntCast(_ldx, &ldx);
687:   if (ldx < m) {
688:     MatDenseGetArrayRead(B,&b);
689:     PetscMalloc1(nrhs * m, &y);
690:     if (ldb == m) {
691:       PetscArraycpy(y,b,ldb*nrhs);
692:     } else {
693:       for (PetscInt j = 0; j < nrhs; j++) {
694:         PetscArraycpy(&y[j*m],&b[j*ldb],m);
695:       }
696:     }
697:     ldy = m;
698:     MatDenseRestoreArrayRead(B,&b);
699:   } else {
700:     if (ldb == ldx) {
701:       MatCopy(B, X, SAME_NONZERO_PATTERN);
702:       MatDenseGetArray(X,&y);
703:     } else {
704:       MatDenseGetArray(X,&y);
705:       MatDenseGetArrayRead(B,&b);
706:       for (PetscInt j = 0; j < nrhs; j++) {
707:         PetscArraycpy(&y[j*ldx],&b[j*ldb],m);
708:       }
709:       MatDenseRestoreArrayRead(B,&b);
710:     }
711:     ldy = ldx;
712:   }
713:   *_y    = y;
714:   *_ldy = ldy;
715:   *_k    = k;
716:   *_m    = m;
717:   *_nrhs = nrhs;
718:   return(0);
719: }

721: static PetscErrorCode MatMatSolve_SeqDense_TearDown(Mat A, Mat B, Mat X, PetscScalar **_y, PetscBLASInt *_ldy, PetscBLASInt *_m, PetscBLASInt *_nrhs, PetscBLASInt *_k)
722: {
723:   PetscScalar       *y;
724:   PetscInt          _ldx;
725:   PetscBLASInt      k,ldy,nrhs,ldx=0;
726:   PetscErrorCode    ierr;

729:   y    = *_y;
730:   *_y  = NULL;
731:   k    = *_k;
732:   ldy = *_ldy;
733:   nrhs = *_nrhs;
734:   MatDenseGetLDA(X,&_ldx);
735:   PetscBLASIntCast(_ldx, &ldx);
736:   if (ldx != ldy) {
737:     PetscScalar *xv;
738:     MatDenseGetArray(X,&xv);
739:     for (PetscInt j = 0; j < nrhs; j++) {
740:       PetscArraycpy(&xv[j*ldx],&y[j*ldy],k);
741:     }
742:     MatDenseRestoreArray(X,&xv);
743:     PetscFree(y);
744:   } else {
745:     MatDenseRestoreArray(X,&y);
746:   }
747:   return(0);
748: }

750: static PetscErrorCode MatMatSolve_SeqDense_LU(Mat A, Mat B, Mat X)
751: {
752:   PetscScalar    *y;
753:   PetscBLASInt   m, k, ldy, nrhs;

757:   MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k);
758:   MatSolve_SeqDense_Internal_LU(A, y, ldy, m, nrhs, k, PETSC_FALSE);
759:   MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k);
760:   return(0);
761: }

763: static PetscErrorCode MatMatSolveTranspose_SeqDense_LU(Mat A, Mat B, Mat X)
764: {
765:   PetscScalar    *y;
766:   PetscBLASInt   m, k, ldy, nrhs;

770:   MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k);
771:   MatSolve_SeqDense_Internal_LU(A, y, ldy, m, nrhs, k, PETSC_TRUE);
772:   MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k);
773:   return(0);
774: }

776: static PetscErrorCode MatMatSolve_SeqDense_Cholesky(Mat A, Mat B, Mat X)
777: {
778:   PetscScalar    *y;
779:   PetscBLASInt   m, k, ldy, nrhs;

783:   MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k);
784:   MatSolve_SeqDense_Internal_Cholesky(A, y, ldy, m, nrhs, k, PETSC_FALSE);
785:   MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k);
786:   return(0);
787: }

789: static PetscErrorCode MatMatSolveTranspose_SeqDense_Cholesky(Mat A, Mat B, Mat X)
790: {
791:   PetscScalar    *y;
792:   PetscBLASInt   m, k, ldy, nrhs;

796:   MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k);
797:   MatSolve_SeqDense_Internal_Cholesky(A, y, ldy, m, nrhs, k, PETSC_TRUE);
798:   MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k);
799:   return(0);
800: }

802: static PetscErrorCode MatMatSolve_SeqDense_QR(Mat A, Mat B, Mat X)
803: {
804:   PetscScalar    *y;
805:   PetscBLASInt   m, k, ldy, nrhs;

809:   MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k);
810:   MatSolve_SeqDense_Internal_QR(A, y, ldy, m, nrhs, k);
811:   MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k);
812:   return(0);
813: }

815: static PetscErrorCode MatMatSolveTranspose_SeqDense_QR(Mat A, Mat B, Mat X)
816: {
817:   PetscScalar    *y;
818:   PetscBLASInt   m, k, ldy, nrhs;

822:   MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k);
823:   MatSolveTranspose_SeqDense_Internal_QR(A, y, ldy, m, nrhs, k);
824:   MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k);
825:   return(0);
826: }

828: static PetscErrorCode MatConjugate_SeqDense(Mat);

830: /* ---------------------------------------------------------------*/
831: /* COMMENT: I have chosen to hide row permutation in the pivots,
832:    rather than put it in the Mat->row slot.*/
833: PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo)
834: {
835:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
837:   PetscBLASInt   n,m,info;

840:   PetscBLASIntCast(A->cmap->n,&n);
841:   PetscBLASIntCast(A->rmap->n,&m);
842:   if (!mat->pivots) {
843:     PetscMalloc1(A->rmap->n,&mat->pivots);
844:     PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
845:   }
846:   if (!A->rmap->n || !A->cmap->n) return(0);
847:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
848:   PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info));
849:   PetscFPTrapPop();

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

854:   A->ops->solve             = MatSolve_SeqDense_LU;
855:   A->ops->matsolve          = MatMatSolve_SeqDense_LU;
856:   A->ops->solvetranspose    = MatSolveTranspose_SeqDense_LU;
857:   A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_LU;
858:   A->factortype             = MAT_FACTOR_LU;

860:   PetscFree(A->solvertype);
861:   PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);

863:   PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);
864:   return(0);
865: }

867: static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
868: {
869:   MatFactorInfo  info;

873:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
874:   (*fact->ops->lufactor)(fact,NULL,NULL,&info);
875:   return(0);
876: }

878: PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info)
879: {
881:   fact->preallocated           = PETSC_TRUE;
882:   fact->assembled              = PETSC_TRUE;
883:   fact->ops->lufactornumeric   = MatLUFactorNumeric_SeqDense;
884:   return(0);
885: }

887: /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */
888: PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo)
889: {
890:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
892:   PetscBLASInt   info,n;

895:   PetscBLASIntCast(A->cmap->n,&n);
896:   if (!A->rmap->n || !A->cmap->n) return(0);
897:   if (A->spd) {
898:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
899:     PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info));
900:     PetscFPTrapPop();
901: #if defined(PETSC_USE_COMPLEX)
902:   } else if (A->hermitian) {
903:     if (!mat->pivots) {
904:       PetscMalloc1(A->rmap->n,&mat->pivots);
905:       PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
906:     }
907:     if (!mat->fwork) {
908:       PetscScalar dummy;

910:       mat->lfwork = -1;
911:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
912:       PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
913:       PetscFPTrapPop();
914:       mat->lfwork = (PetscInt)PetscRealPart(dummy);
915:       PetscMalloc1(mat->lfwork,&mat->fwork);
916:       PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
917:     }
918:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
919:     PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
920:     PetscFPTrapPop();
921: #endif
922:   } else { /* symmetric case */
923:     if (!mat->pivots) {
924:       PetscMalloc1(A->rmap->n,&mat->pivots);
925:       PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
926:     }
927:     if (!mat->fwork) {
928:       PetscScalar dummy;

930:       mat->lfwork = -1;
931:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
932:       PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
933:       PetscFPTrapPop();
934:       mat->lfwork = (PetscInt)PetscRealPart(dummy);
935:       PetscMalloc1(mat->lfwork,&mat->fwork);
936:       PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
937:     }
938:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
939:     PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
940:     PetscFPTrapPop();
941:   }
942:   if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1);

944:   A->ops->solve             = MatSolve_SeqDense_Cholesky;
945:   A->ops->matsolve          = MatMatSolve_SeqDense_Cholesky;
946:   A->ops->solvetranspose    = MatSolveTranspose_SeqDense_Cholesky;
947:   A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_Cholesky;
948:   A->factortype             = MAT_FACTOR_CHOLESKY;

950:   PetscFree(A->solvertype);
951:   PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);

953:   PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
954:   return(0);
955: }

957: static PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
958: {
960:   MatFactorInfo  info;

963:   info.fill = 1.0;

965:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
966:   (*fact->ops->choleskyfactor)(fact,NULL,&info);
967:   return(0);
968: }

970: PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
971: {
973:   fact->assembled                  = PETSC_TRUE;
974:   fact->preallocated               = PETSC_TRUE;
975:   fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
976:   return(0);
977: }

979: PetscErrorCode MatQRFactor_SeqDense(Mat A,IS col,const MatFactorInfo *minfo)
980: {
981:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
983:   PetscBLASInt   n,m,info, min, max;

986:   PetscBLASIntCast(A->cmap->n,&n);
987:   PetscBLASIntCast(A->rmap->n,&m);
988:   max = PetscMax(m, n);
989:   min = PetscMin(m, n);
990:   if (!mat->tau) {
991:     PetscMalloc1(min,&mat->tau);
992:     PetscLogObjectMemory((PetscObject)A,min*sizeof(PetscScalar));
993:   }
994:   if (!mat->pivots) {
995:     PetscMalloc1(n,&mat->pivots);
996:     PetscLogObjectMemory((PetscObject)A,n*sizeof(PetscScalar));
997:   }
998:   if (!mat->qrrhs) {
999:     MatCreateVecs(A, NULL, &(mat->qrrhs));
1000:   }
1001:   if (!A->rmap->n || !A->cmap->n) return(0);
1002:   if (!mat->fwork) {
1003:     PetscScalar dummy;

1005:     mat->lfwork = -1;
1006:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1007:     PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&m,&n,mat->v,&mat->lda,mat->tau,&dummy,&mat->lfwork,&info));
1008:     PetscFPTrapPop();
1009:     mat->lfwork = (PetscInt)PetscRealPart(dummy);
1010:     PetscMalloc1(mat->lfwork,&mat->fwork);
1011:     PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
1012:   }
1013:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1014:   PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&m,&n,mat->v,&mat->lda,mat->tau,mat->fwork,&mat->lfwork,&info));
1015:   PetscFPTrapPop();
1016:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to QR factorization");
1017:   // TODO: try to estimate rank or test for and use geqp3 for rank revealing QR.  For now just say rank is min of m and n
1018:   mat->rank = min;

1020:   A->ops->solve             = MatSolve_SeqDense_QR;
1021:   A->ops->matsolve          = MatMatSolve_SeqDense_QR;
1022:   A->factortype             = MAT_FACTOR_QR;
1023:   if (m == n) {
1024:     A->ops->solvetranspose    = MatSolveTranspose_SeqDense_QR;
1025:     A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_QR;
1026:   }

1028:   PetscFree(A->solvertype);
1029:   PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);

1031:   PetscLogFlops(2.0*min*min*(max-min/3.0));
1032:   return(0);
1033: }

1035: static PetscErrorCode MatQRFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
1036: {
1038:   MatFactorInfo  info;

1041:   info.fill = 1.0;

1043:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
1044:   PetscUseMethod(fact,"MatQRFactor_C",(Mat,IS,const MatFactorInfo *),(fact,NULL,&info));
1045:   return(0);
1046: }

1048: PetscErrorCode MatQRFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
1049: {

1053:   fact->assembled                  = PETSC_TRUE;
1054:   fact->preallocated               = PETSC_TRUE;
1055:   PetscObjectComposeFunction((PetscObject)fact,"MatQRFactorNumeric_C",MatQRFactorNumeric_SeqDense);
1056:   return(0);
1057: }

1059: /* uses LAPACK */
1060: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact)
1061: {

1065:   MatCreate(PetscObjectComm((PetscObject)A),fact);
1066:   MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
1067:   MatSetType(*fact,MATDENSE);
1068:   (*fact)->trivialsymbolic = PETSC_TRUE;
1069:   if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
1070:     (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
1071:     (*fact)->ops->ilufactorsymbolic = MatLUFactorSymbolic_SeqDense;
1072:   } else if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) {
1073:     (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
1074:   } else if (ftype == MAT_FACTOR_QR) {
1075:     PetscObjectComposeFunction((PetscObject)(*fact),"MatQRFactorSymbolic_C",MatQRFactorSymbolic_SeqDense);
1076:   }
1077:   (*fact)->factortype = ftype;

1079:   PetscFree((*fact)->solvertype);
1080:   PetscStrallocpy(MATSOLVERPETSC,&(*fact)->solvertype);
1081:   PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&(*fact)->preferredordering[MAT_FACTOR_LU]);
1082:   PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&(*fact)->preferredordering[MAT_FACTOR_ILU]);
1083:   PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&(*fact)->preferredordering[MAT_FACTOR_CHOLESKY]);
1084:   PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&(*fact)->preferredordering[MAT_FACTOR_ICC]);
1085:   return(0);
1086: }

1088: /* ------------------------------------------------------------------*/
1089: static PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx)
1090: {
1091:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1092:   PetscScalar       *x,*v = mat->v,zero = 0.0,xt;
1093:   const PetscScalar *b;
1094:   PetscErrorCode    ierr;
1095:   PetscInt          m = A->rmap->n,i;
1096:   PetscBLASInt      o = 1,bm = 0;

1099: #if defined(PETSC_HAVE_CUDA)
1100:   if (A->offloadmask == PETSC_OFFLOAD_GPU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented");
1101: #endif
1102:   if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */
1103:   PetscBLASIntCast(m,&bm);
1104:   if (flag & SOR_ZERO_INITIAL_GUESS) {
1105:     /* this is a hack fix, should have another version without the second BLASdotu */
1106:     VecSet(xx,zero);
1107:   }
1108:   VecGetArray(xx,&x);
1109:   VecGetArrayRead(bb,&b);
1110:   its  = its*lits;
1111:   if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
1112:   while (its--) {
1113:     if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
1114:       for (i=0; i<m; i++) {
1115:         PetscStackCallBLAS("BLASdotu",xt   = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o));
1116:         x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
1117:       }
1118:     }
1119:     if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
1120:       for (i=m-1; i>=0; i--) {
1121:         PetscStackCallBLAS("BLASdotu",xt   = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o));
1122:         x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
1123:       }
1124:     }
1125:   }
1126:   VecRestoreArrayRead(bb,&b);
1127:   VecRestoreArray(xx,&x);
1128:   return(0);
1129: }

1131: /* -----------------------------------------------------------------*/
1132: PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)
1133: {
1134:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1135:   const PetscScalar *v   = mat->v,*x;
1136:   PetscScalar       *y;
1137:   PetscErrorCode    ierr;
1138:   PetscBLASInt      m, n,_One=1;
1139:   PetscScalar       _DOne=1.0,_DZero=0.0;

1142:   PetscBLASIntCast(A->rmap->n,&m);
1143:   PetscBLASIntCast(A->cmap->n,&n);
1144:   VecGetArrayRead(xx,&x);
1145:   VecGetArrayWrite(yy,&y);
1146:   if (!A->rmap->n || !A->cmap->n) {
1147:     PetscBLASInt i;
1148:     for (i=0; i<n; i++) y[i] = 0.0;
1149:   } else {
1150:     PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One));
1151:     PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);
1152:   }
1153:   VecRestoreArrayRead(xx,&x);
1154:   VecRestoreArrayWrite(yy,&y);
1155:   return(0);
1156: }

1158: PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy)
1159: {
1160:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1161:   PetscScalar       *y,_DOne=1.0,_DZero=0.0;
1162:   PetscErrorCode    ierr;
1163:   PetscBLASInt      m, n, _One=1;
1164:   const PetscScalar *v = mat->v,*x;

1167:   PetscBLASIntCast(A->rmap->n,&m);
1168:   PetscBLASIntCast(A->cmap->n,&n);
1169:   VecGetArrayRead(xx,&x);
1170:   VecGetArrayWrite(yy,&y);
1171:   if (!A->rmap->n || !A->cmap->n) {
1172:     PetscBLASInt i;
1173:     for (i=0; i<m; i++) y[i] = 0.0;
1174:   } else {
1175:     PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One));
1176:     PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);
1177:   }
1178:   VecRestoreArrayRead(xx,&x);
1179:   VecRestoreArrayWrite(yy,&y);
1180:   return(0);
1181: }

1183: PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
1184: {
1185:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1186:   const PetscScalar *v = mat->v,*x;
1187:   PetscScalar       *y,_DOne=1.0;
1188:   PetscErrorCode    ierr;
1189:   PetscBLASInt      m, n, _One=1;

1192:   PetscBLASIntCast(A->rmap->n,&m);
1193:   PetscBLASIntCast(A->cmap->n,&n);
1194:   VecCopy(zz,yy);
1195:   if (!A->rmap->n || !A->cmap->n) return(0);
1196:   VecGetArrayRead(xx,&x);
1197:   VecGetArray(yy,&y);
1198:   PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
1199:   VecRestoreArrayRead(xx,&x);
1200:   VecRestoreArray(yy,&y);
1201:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
1202:   return(0);
1203: }

1205: PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
1206: {
1207:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1208:   const PetscScalar *v = mat->v,*x;
1209:   PetscScalar       *y;
1210:   PetscErrorCode    ierr;
1211:   PetscBLASInt      m, n, _One=1;
1212:   PetscScalar       _DOne=1.0;

1215:   PetscBLASIntCast(A->rmap->n,&m);
1216:   PetscBLASIntCast(A->cmap->n,&n);
1217:   VecCopy(zz,yy);
1218:   if (!A->rmap->n || !A->cmap->n) return(0);
1219:   VecGetArrayRead(xx,&x);
1220:   VecGetArray(yy,&y);
1221:   PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
1222:   VecRestoreArrayRead(xx,&x);
1223:   VecRestoreArray(yy,&y);
1224:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
1225:   return(0);
1226: }

1228: /* -----------------------------------------------------------------*/
1229: static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
1230: {
1231:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1233:   PetscInt       i;

1236:   *ncols = A->cmap->n;
1237:   if (cols) {
1238:     PetscMalloc1(A->cmap->n,cols);
1239:     for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
1240:   }
1241:   if (vals) {
1242:     const PetscScalar *v;

1244:     MatDenseGetArrayRead(A,&v);
1245:     PetscMalloc1(A->cmap->n,vals);
1246:     v   += row;
1247:     for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
1248:     MatDenseRestoreArrayRead(A,&v);
1249:   }
1250:   return(0);
1251: }

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

1258:   if (ncols) *ncols = 0;
1259:   if (cols) {PetscFree(*cols);}
1260:   if (vals) {PetscFree(*vals); }
1261:   return(0);
1262: }
1263: /* ----------------------------------------------------------------*/
1264: static PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
1265: {
1266:   Mat_SeqDense     *mat = (Mat_SeqDense*)A->data;
1267:   PetscScalar      *av;
1268:   PetscInt         i,j,idx=0;
1269: #if defined(PETSC_HAVE_CUDA)
1270:   PetscOffloadMask oldf;
1271: #endif
1272:   PetscErrorCode   ierr;

1275:   MatDenseGetArray(A,&av);
1276:   if (!mat->roworiented) {
1277:     if (addv == INSERT_VALUES) {
1278:       for (j=0; j<n; j++) {
1279:         if (indexn[j] < 0) {idx += m; continue;}
1280:         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);
1281:         for (i=0; i<m; i++) {
1282:           if (indexm[i] < 0) {idx++; continue;}
1283:           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);
1284:           av[indexn[j]*mat->lda + indexm[i]] = v[idx++];
1285:         }
1286:       }
1287:     } else {
1288:       for (j=0; j<n; j++) {
1289:         if (indexn[j] < 0) {idx += m; continue;}
1290:         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);
1291:         for (i=0; i<m; i++) {
1292:           if (indexm[i] < 0) {idx++; continue;}
1293:           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);
1294:           av[indexn[j]*mat->lda + indexm[i]] += v[idx++];
1295:         }
1296:       }
1297:     }
1298:   } else {
1299:     if (addv == INSERT_VALUES) {
1300:       for (i=0; i<m; i++) {
1301:         if (indexm[i] < 0) { idx += n; continue;}
1302:         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);
1303:         for (j=0; j<n; j++) {
1304:           if (indexn[j] < 0) { idx++; continue;}
1305:           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);
1306:           av[indexn[j]*mat->lda + indexm[i]] = v[idx++];
1307:         }
1308:       }
1309:     } else {
1310:       for (i=0; i<m; i++) {
1311:         if (indexm[i] < 0) { idx += n; continue;}
1312:         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);
1313:         for (j=0; j<n; j++) {
1314:           if (indexn[j] < 0) { idx++; continue;}
1315:           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);
1316:           av[indexn[j]*mat->lda + indexm[i]] += v[idx++];
1317:         }
1318:       }
1319:     }
1320:   }
1321:   /* hack to prevent unneeded copy to the GPU while returning the array */
1322: #if defined(PETSC_HAVE_CUDA)
1323:   oldf = A->offloadmask;
1324:   A->offloadmask = PETSC_OFFLOAD_GPU;
1325: #endif
1326:   MatDenseRestoreArray(A,&av);
1327: #if defined(PETSC_HAVE_CUDA)
1328:   A->offloadmask = (oldf == PETSC_OFFLOAD_UNALLOCATED ? PETSC_OFFLOAD_UNALLOCATED : PETSC_OFFLOAD_CPU);
1329: #endif
1330:   return(0);
1331: }

1333: static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
1334: {
1335:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1336:   const PetscScalar *vv;
1337:   PetscInt          i,j;
1338:   PetscErrorCode    ierr;

1341:   MatDenseGetArrayRead(A,&vv);
1342:   /* row-oriented output */
1343:   for (i=0; i<m; i++) {
1344:     if (indexm[i] < 0) {v += n;continue;}
1345:     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);
1346:     for (j=0; j<n; j++) {
1347:       if (indexn[j] < 0) {v++; continue;}
1348:       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);
1349:       *v++ = vv[indexn[j]*mat->lda + indexm[i]];
1350:     }
1351:   }
1352:   MatDenseRestoreArrayRead(A,&vv);
1353:   return(0);
1354: }

1356: /* -----------------------------------------------------------------*/

1358: PetscErrorCode MatView_Dense_Binary(Mat mat,PetscViewer viewer)
1359: {
1360:   PetscErrorCode    ierr;
1361:   PetscBool         skipHeader;
1362:   PetscViewerFormat format;
1363:   PetscInt          header[4],M,N,m,lda,i,j,k;
1364:   const PetscScalar *v;
1365:   PetscScalar       *vwork;

1368:   PetscViewerSetUp(viewer);
1369:   PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
1370:   PetscViewerGetFormat(viewer,&format);
1371:   if (skipHeader) format = PETSC_VIEWER_NATIVE;

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

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

1380:   MatGetLocalSize(mat,&m,NULL);
1381:   if (format != PETSC_VIEWER_NATIVE) {
1382:     PetscInt nnz = m*N, *iwork;
1383:     /* store row lengths for each row */
1384:     PetscMalloc1(nnz,&iwork);
1385:     for (i=0; i<m; i++) iwork[i] = N;
1386:     PetscViewerBinaryWriteAll(viewer,iwork,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1387:     /* store column indices (zero start index) */
1388:     for (k=0, i=0; i<m; i++)
1389:       for (j=0; j<N; j++, k++)
1390:         iwork[k] = j;
1391:     PetscViewerBinaryWriteAll(viewer,iwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1392:     PetscFree(iwork);
1393:   }
1394:   /* store matrix values as a dense matrix in row major order */
1395:   PetscMalloc1(m*N,&vwork);
1396:   MatDenseGetArrayRead(mat,&v);
1397:   MatDenseGetLDA(mat,&lda);
1398:   for (k=0, i=0; i<m; i++)
1399:     for (j=0; j<N; j++, k++)
1400:       vwork[k] = v[i+lda*j];
1401:   MatDenseRestoreArrayRead(mat,&v);
1402:   PetscViewerBinaryWriteAll(viewer,vwork,m*N,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);
1403:   PetscFree(vwork);
1404:   return(0);
1405: }

1407: PetscErrorCode MatLoad_Dense_Binary(Mat mat,PetscViewer viewer)
1408: {
1410:   PetscBool      skipHeader;
1411:   PetscInt       header[4],M,N,m,nz,lda,i,j,k;
1412:   PetscInt       rows,cols;
1413:   PetscScalar    *v,*vwork;

1416:   PetscViewerSetUp(viewer);
1417:   PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);

1419:   if (!skipHeader) {
1420:     PetscViewerBinaryRead(viewer,header,4,NULL,PETSC_INT);
1421:     if (header[0] != MAT_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Not a matrix object in file");
1422:     M = header[1]; N = header[2];
1423:     if (M < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix row size (%D) in file is negative",M);
1424:     if (N < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix column size (%D) in file is negative",N);
1425:     nz = header[3];
1426:     if (nz != MATRIX_BINARY_FORMAT_DENSE && nz < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Unknown matrix format %D in file",nz);
1427:   } else {
1428:     MatGetSize(mat,&M,&N);
1429:     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");
1430:     nz = MATRIX_BINARY_FORMAT_DENSE;
1431:   }

1433:   /* setup global sizes if not set */
1434:   if (mat->rmap->N < 0) mat->rmap->N = M;
1435:   if (mat->cmap->N < 0) mat->cmap->N = N;
1436:   MatSetUp(mat);
1437:   /* check if global sizes are correct */
1438:   MatGetSize(mat,&rows,&cols);
1439:   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);

1441:   MatGetSize(mat,NULL,&N);
1442:   MatGetLocalSize(mat,&m,NULL);
1443:   MatDenseGetArray(mat,&v);
1444:   MatDenseGetLDA(mat,&lda);
1445:   if (nz == MATRIX_BINARY_FORMAT_DENSE) {  /* matrix in file is dense format */
1446:     PetscInt nnz = m*N;
1447:     /* read in matrix values */
1448:     PetscMalloc1(nnz,&vwork);
1449:     PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);
1450:     /* store values in column major order */
1451:     for (j=0; j<N; j++)
1452:       for (i=0; i<m; i++)
1453:         v[i+lda*j] = vwork[i*N+j];
1454:     PetscFree(vwork);
1455:   } else { /* matrix in file is sparse format */
1456:     PetscInt nnz = 0, *rlens, *icols;
1457:     /* read in row lengths */
1458:     PetscMalloc1(m,&rlens);
1459:     PetscViewerBinaryReadAll(viewer,rlens,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1460:     for (i=0; i<m; i++) nnz += rlens[i];
1461:     /* read in column indices and values */
1462:     PetscMalloc2(nnz,&icols,nnz,&vwork);
1463:     PetscViewerBinaryReadAll(viewer,icols,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1464:     PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);
1465:     /* store values in column major order */
1466:     for (k=0, i=0; i<m; i++)
1467:       for (j=0; j<rlens[i]; j++, k++)
1468:         v[i+lda*icols[k]] = vwork[k];
1469:     PetscFree(rlens);
1470:     PetscFree2(icols,vwork);
1471:   }
1472:   MatDenseRestoreArray(mat,&v);
1473:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
1474:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
1475:   return(0);
1476: }

1478: PetscErrorCode MatLoad_SeqDense(Mat newMat, PetscViewer viewer)
1479: {
1480:   PetscBool      isbinary, ishdf5;

1486:   /* force binary viewer to load .info file if it has not yet done so */
1487:   PetscViewerSetUp(viewer);
1488:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1489:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,  &ishdf5);
1490:   if (isbinary) {
1491:     MatLoad_Dense_Binary(newMat,viewer);
1492:   } else if (ishdf5) {
1493: #if defined(PETSC_HAVE_HDF5)
1494:     MatLoad_Dense_HDF5(newMat,viewer);
1495: #else
1496:     SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1497: #endif
1498:   } else {
1499:     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);
1500:   }
1501:   return(0);
1502: }

1504: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
1505: {
1506:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
1507:   PetscErrorCode    ierr;
1508:   PetscInt          i,j;
1509:   const char        *name;
1510:   PetscScalar       *v,*av;
1511:   PetscViewerFormat format;
1512: #if defined(PETSC_USE_COMPLEX)
1513:   PetscBool         allreal = PETSC_TRUE;
1514: #endif

1517:   MatDenseGetArrayRead(A,(const PetscScalar**)&av);
1518:   PetscViewerGetFormat(viewer,&format);
1519:   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1520:     return(0);  /* do nothing for now */
1521:   } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1522:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1523:     for (i=0; i<A->rmap->n; i++) {
1524:       v    = av + i;
1525:       PetscViewerASCIIPrintf(viewer,"row %D:",i);
1526:       for (j=0; j<A->cmap->n; j++) {
1527: #if defined(PETSC_USE_COMPLEX)
1528:         if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
1529:           PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1530:         } else if (PetscRealPart(*v)) {
1531:           PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));
1532:         }
1533: #else
1534:         if (*v) {
1535:           PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);
1536:         }
1537: #endif
1538:         v += a->lda;
1539:       }
1540:       PetscViewerASCIIPrintf(viewer,"\n");
1541:     }
1542:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1543:   } else {
1544:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1545: #if defined(PETSC_USE_COMPLEX)
1546:     /* determine if matrix has all real values */
1547:     v = av;
1548:     for (i=0; i<A->rmap->n*A->cmap->n; i++) {
1549:       if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;}
1550:     }
1551: #endif
1552:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
1553:       PetscObjectGetName((PetscObject)A,&name);
1554:       PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);
1555:       PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);
1556:       PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
1557:     }

1559:     for (i=0; i<A->rmap->n; i++) {
1560:       v = av + i;
1561:       for (j=0; j<A->cmap->n; j++) {
1562: #if defined(PETSC_USE_COMPLEX)
1563:         if (allreal) {
1564:           PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));
1565:         } else {
1566:           PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1567:         }
1568: #else
1569:         PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);
1570: #endif
1571:         v += a->lda;
1572:       }
1573:       PetscViewerASCIIPrintf(viewer,"\n");
1574:     }
1575:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
1576:       PetscViewerASCIIPrintf(viewer,"];\n");
1577:     }
1578:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1579:   }
1580:   MatDenseRestoreArrayRead(A,(const PetscScalar**)&av);
1581:   PetscViewerFlush(viewer);
1582:   return(0);
1583: }

1585: #include <petscdraw.h>
1586: static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa)
1587: {
1588:   Mat               A  = (Mat) Aa;
1589:   PetscErrorCode    ierr;
1590:   PetscInt          m  = A->rmap->n,n = A->cmap->n,i,j;
1591:   int               color = PETSC_DRAW_WHITE;
1592:   const PetscScalar *v;
1593:   PetscViewer       viewer;
1594:   PetscReal         xl,yl,xr,yr,x_l,x_r,y_l,y_r;
1595:   PetscViewerFormat format;

1598:   PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
1599:   PetscViewerGetFormat(viewer,&format);
1600:   PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);

1602:   /* Loop over matrix elements drawing boxes */
1603:   MatDenseGetArrayRead(A,&v);
1604:   if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1605:     PetscDrawCollectiveBegin(draw);
1606:     /* Blue for negative and Red for positive */
1607:     for (j = 0; j < n; j++) {
1608:       x_l = j; x_r = x_l + 1.0;
1609:       for (i = 0; i < m; i++) {
1610:         y_l = m - i - 1.0;
1611:         y_r = y_l + 1.0;
1612:         if (PetscRealPart(v[j*m+i]) >  0.) color = PETSC_DRAW_RED;
1613:         else if (PetscRealPart(v[j*m+i]) <  0.) color = PETSC_DRAW_BLUE;
1614:         else continue;
1615:         PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1616:       }
1617:     }
1618:     PetscDrawCollectiveEnd(draw);
1619:   } else {
1620:     /* use contour shading to indicate magnitude of values */
1621:     /* first determine max of all nonzero values */
1622:     PetscReal minv = 0.0, maxv = 0.0;
1623:     PetscDraw popup;

1625:     for (i=0; i < m*n; i++) {
1626:       if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1627:     }
1628:     if (minv >= maxv) maxv = minv + PETSC_SMALL;
1629:     PetscDrawGetPopup(draw,&popup);
1630:     PetscDrawScalePopup(popup,minv,maxv);

1632:     PetscDrawCollectiveBegin(draw);
1633:     for (j=0; j<n; j++) {
1634:       x_l = j;
1635:       x_r = x_l + 1.0;
1636:       for (i=0; i<m; i++) {
1637:         y_l   = m - i - 1.0;
1638:         y_r   = y_l + 1.0;
1639:         color = PetscDrawRealToColor(PetscAbsScalar(v[j*m+i]),minv,maxv);
1640:         PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1641:       }
1642:     }
1643:     PetscDrawCollectiveEnd(draw);
1644:   }
1645:   MatDenseRestoreArrayRead(A,&v);
1646:   return(0);
1647: }

1649: static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
1650: {
1651:   PetscDraw      draw;
1652:   PetscBool      isnull;
1653:   PetscReal      xr,yr,xl,yl,h,w;

1657:   PetscViewerDrawGetDraw(viewer,0,&draw);
1658:   PetscDrawIsNull(draw,&isnull);
1659:   if (isnull) return(0);

1661:   xr   = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
1662:   xr  += w;          yr += h;        xl = -w;     yl = -h;
1663:   PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
1664:   PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
1665:   PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
1666:   PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
1667:   PetscDrawSave(draw);
1668:   return(0);
1669: }

1671: PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
1672: {
1674:   PetscBool      iascii,isbinary,isdraw;

1677:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1678:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1679:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);

1681:   if (iascii) {
1682:     MatView_SeqDense_ASCII(A,viewer);
1683:   } else if (isbinary) {
1684:     MatView_Dense_Binary(A,viewer);
1685:   } else if (isdraw) {
1686:     MatView_SeqDense_Draw(A,viewer);
1687:   }
1688:   return(0);
1689: }

1691: static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A,const PetscScalar *array)
1692: {
1693:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;

1696:   if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
1697:   if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
1698:   if (a->unplacedarray) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreArray() first");
1699:   a->unplacedarray       = a->v;
1700:   a->unplaced_user_alloc = a->user_alloc;
1701:   a->v                   = (PetscScalar*) array;
1702:   a->user_alloc          = PETSC_TRUE;
1703: #if defined(PETSC_HAVE_CUDA)
1704:   A->offloadmask = PETSC_OFFLOAD_CPU;
1705: #endif
1706:   return(0);
1707: }

1709: static PetscErrorCode MatDenseResetArray_SeqDense(Mat A)
1710: {
1711:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;

1714:   if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
1715:   if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
1716:   a->v             = a->unplacedarray;
1717:   a->user_alloc    = a->unplaced_user_alloc;
1718:   a->unplacedarray = NULL;
1719: #if defined(PETSC_HAVE_CUDA)
1720:   A->offloadmask = PETSC_OFFLOAD_CPU;
1721: #endif
1722:   return(0);
1723: }

1725: static PetscErrorCode MatDenseReplaceArray_SeqDense(Mat A,const PetscScalar *array)
1726: {
1727:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

1731:   if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
1732:   if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
1733:   if (!a->user_alloc) { PetscFree(a->v); }
1734:   a->v           = (PetscScalar*) array;
1735:   a->user_alloc  = PETSC_FALSE;
1736: #if defined(PETSC_HAVE_CUDA)
1737:   A->offloadmask = PETSC_OFFLOAD_CPU;
1738: #endif
1739:   return(0);
1740: }

1742: PetscErrorCode MatDestroy_SeqDense(Mat mat)
1743: {
1744:   Mat_SeqDense   *l = (Mat_SeqDense*)mat->data;

1748: #if defined(PETSC_USE_LOG)
1749:   PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n);
1750: #endif
1751:   VecDestroy(&(l->qrrhs));
1752:   PetscFree(l->tau);
1753:   PetscFree(l->pivots);
1754:   PetscFree(l->fwork);
1755:   MatDestroy(&l->ptapwork);
1756:   if (!l->user_alloc) {PetscFree(l->v);}
1757:   if (!l->unplaced_user_alloc) {PetscFree(l->unplacedarray);}
1758:   if (l->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
1759:   if (l->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
1760:   VecDestroy(&l->cvec);
1761:   MatDestroy(&l->cmat);
1762:   PetscFree(mat->data);

1764:   PetscObjectChangeTypeName((PetscObject)mat,NULL);
1765:   PetscObjectComposeFunction((PetscObject)mat,"MatQRFactor_C",NULL);
1766:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetLDA_C",NULL);
1767:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseSetLDA_C",NULL);
1768:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);
1769:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);
1770:   PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL);
1771:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL);
1772:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseReplaceArray_C",NULL);
1773:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayRead_C",NULL);
1774:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayRead_C",NULL);
1775:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayWrite_C",NULL);
1776:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayWrite_C",NULL);
1777:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);
1778: #if defined(PETSC_HAVE_ELEMENTAL)
1779:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);
1780: #endif
1781: #if defined(PETSC_HAVE_SCALAPACK)
1782:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_scalapack_C",NULL);
1783: #endif
1784: #if defined(PETSC_HAVE_CUDA)
1785:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqdensecuda_C",NULL);
1786:   PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",NULL);
1787:   PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdense_C",NULL);
1788: #endif
1789:   PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);
1790:   PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqaij_seqdense_C",NULL);
1791:   PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdense_seqdense_C",NULL);
1792:   PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqbaij_seqdense_C",NULL);
1793:   PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqsbaij_seqdense_C",NULL);

1795:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumn_C",NULL);
1796:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumn_C",NULL);
1797:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVec_C",NULL);
1798:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVec_C",NULL);
1799:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVecRead_C",NULL);
1800:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVecRead_C",NULL);
1801:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVecWrite_C",NULL);
1802:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVecWrite_C",NULL);
1803:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetSubMatrix_C",NULL);
1804:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreSubMatrix_C",NULL);
1805:   return(0);
1806: }

1808: static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
1809: {
1810:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1812:   PetscInt       k,j,m = A->rmap->n, M = mat->lda, n = A->cmap->n;
1813:   PetscScalar    *v,tmp;

1816:   if (reuse == MAT_INPLACE_MATRIX) {
1817:     if (m == n) { /* in place transpose */
1818:       MatDenseGetArray(A,&v);
1819:       for (j=0; j<m; j++) {
1820:         for (k=0; k<j; k++) {
1821:           tmp        = v[j + k*M];
1822:           v[j + k*M] = v[k + j*M];
1823:           v[k + j*M] = tmp;
1824:         }
1825:       }
1826:       MatDenseRestoreArray(A,&v);
1827:     } else { /* reuse memory, temporary allocates new memory */
1828:       PetscScalar *v2;
1829:       PetscLayout tmplayout;

1831:       PetscMalloc1((size_t)m*n,&v2);
1832:       MatDenseGetArray(A,&v);
1833:       for (j=0; j<n; j++) {
1834:         for (k=0; k<m; k++) v2[j + (size_t)k*n] = v[k + (size_t)j*M];
1835:       }
1836:       PetscArraycpy(v,v2,(size_t)m*n);
1837:       PetscFree(v2);
1838:       MatDenseRestoreArray(A,&v);
1839:       /* cleanup size dependent quantities */
1840:       VecDestroy(&mat->cvec);
1841:       MatDestroy(&mat->cmat);
1842:       PetscFree(mat->pivots);
1843:       PetscFree(mat->fwork);
1844:       MatDestroy(&mat->ptapwork);
1845:       /* swap row/col layouts */
1846:       mat->lda  = n;
1847:       tmplayout = A->rmap;
1848:       A->rmap   = A->cmap;
1849:       A->cmap   = tmplayout;
1850:     }
1851:   } else { /* out-of-place transpose */
1852:     Mat          tmat;
1853:     Mat_SeqDense *tmatd;
1854:     PetscScalar  *v2;
1855:     PetscInt     M2;

1857:     if (reuse == MAT_INITIAL_MATRIX) {
1858:       MatCreate(PetscObjectComm((PetscObject)A),&tmat);
1859:       MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);
1860:       MatSetType(tmat,((PetscObject)A)->type_name);
1861:       MatSeqDenseSetPreallocation(tmat,NULL);
1862:     } else tmat = *matout;

1864:     MatDenseGetArrayRead(A,(const PetscScalar**)&v);
1865:     MatDenseGetArray(tmat,&v2);
1866:     tmatd = (Mat_SeqDense*)tmat->data;
1867:     M2    = tmatd->lda;
1868:     for (j=0; j<n; j++) {
1869:       for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M];
1870:     }
1871:     MatDenseRestoreArray(tmat,&v2);
1872:     MatDenseRestoreArrayRead(A,(const PetscScalar**)&v);
1873:     MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
1874:     MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
1875:     *matout = tmat;
1876:   }
1877:   return(0);
1878: }

1880: static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool  *flg)
1881: {
1882:   Mat_SeqDense      *mat1 = (Mat_SeqDense*)A1->data;
1883:   Mat_SeqDense      *mat2 = (Mat_SeqDense*)A2->data;
1884:   PetscInt          i;
1885:   const PetscScalar *v1,*v2;
1886:   PetscErrorCode    ierr;

1889:   if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; return(0);}
1890:   if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; return(0);}
1891:   MatDenseGetArrayRead(A1,&v1);
1892:   MatDenseGetArrayRead(A2,&v2);
1893:   for (i=0; i<A1->cmap->n; i++) {
1894:     PetscArraycmp(v1,v2,A1->rmap->n,flg);
1895:     if (*flg == PETSC_FALSE) return(0);
1896:     v1 += mat1->lda;
1897:     v2 += mat2->lda;
1898:   }
1899:   MatDenseRestoreArrayRead(A1,&v1);
1900:   MatDenseRestoreArrayRead(A2,&v2);
1901:   *flg = PETSC_TRUE;
1902:   return(0);
1903: }

1905: static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
1906: {
1907:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1908:   PetscInt          i,n,len;
1909:   PetscScalar       *x;
1910:   const PetscScalar *vv;
1911:   PetscErrorCode    ierr;

1914:   VecGetSize(v,&n);
1915:   VecGetArray(v,&x);
1916:   len  = PetscMin(A->rmap->n,A->cmap->n);
1917:   MatDenseGetArrayRead(A,&vv);
1918:   if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
1919:   for (i=0; i<len; i++) {
1920:     x[i] = vv[i*mat->lda + i];
1921:   }
1922:   MatDenseRestoreArrayRead(A,&vv);
1923:   VecRestoreArray(v,&x);
1924:   return(0);
1925: }

1927: static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
1928: {
1929:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1930:   const PetscScalar *l,*r;
1931:   PetscScalar       x,*v,*vv;
1932:   PetscErrorCode    ierr;
1933:   PetscInt          i,j,m = A->rmap->n,n = A->cmap->n;

1936:   MatDenseGetArray(A,&vv);
1937:   if (ll) {
1938:     VecGetSize(ll,&m);
1939:     VecGetArrayRead(ll,&l);
1940:     if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size");
1941:     for (i=0; i<m; i++) {
1942:       x = l[i];
1943:       v = vv + i;
1944:       for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;}
1945:     }
1946:     VecRestoreArrayRead(ll,&l);
1947:     PetscLogFlops(1.0*n*m);
1948:   }
1949:   if (rr) {
1950:     VecGetSize(rr,&n);
1951:     VecGetArrayRead(rr,&r);
1952:     if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size");
1953:     for (i=0; i<n; i++) {
1954:       x = r[i];
1955:       v = vv + i*mat->lda;
1956:       for (j=0; j<m; j++) (*v++) *= x;
1957:     }
1958:     VecRestoreArrayRead(rr,&r);
1959:     PetscLogFlops(1.0*n*m);
1960:   }
1961:   MatDenseRestoreArray(A,&vv);
1962:   return(0);
1963: }

1965: PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1966: {
1967:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1968:   PetscScalar       *v,*vv;
1969:   PetscReal         sum  = 0.0;
1970:   PetscInt          lda  =mat->lda,m=A->rmap->n,i,j;
1971:   PetscErrorCode    ierr;

1974:   MatDenseGetArrayRead(A,(const PetscScalar**)&vv);
1975:   v    = vv;
1976:   if (type == NORM_FROBENIUS) {
1977:     if (lda>m) {
1978:       for (j=0; j<A->cmap->n; j++) {
1979:         v = vv+j*lda;
1980:         for (i=0; i<m; i++) {
1981:           sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1982:         }
1983:       }
1984:     } else {
1985: #if defined(PETSC_USE_REAL___FP16)
1986:       PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n;
1987:       PetscStackCallBLAS("BLASnrm2",*nrm = BLASnrm2_(&cnt,v,&one));
1988:     }
1989: #else
1990:       for (i=0; i<A->cmap->n*A->rmap->n; i++) {
1991:         sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1992:       }
1993:     }
1994:     *nrm = PetscSqrtReal(sum);
1995: #endif
1996:     PetscLogFlops(2.0*A->cmap->n*A->rmap->n);
1997:   } else if (type == NORM_1) {
1998:     *nrm = 0.0;
1999:     for (j=0; j<A->cmap->n; j++) {
2000:       v   = vv + j*mat->lda;
2001:       sum = 0.0;
2002:       for (i=0; i<A->rmap->n; i++) {
2003:         sum += PetscAbsScalar(*v);  v++;
2004:       }
2005:       if (sum > *nrm) *nrm = sum;
2006:     }
2007:     PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
2008:   } else if (type == NORM_INFINITY) {
2009:     *nrm = 0.0;
2010:     for (j=0; j<A->rmap->n; j++) {
2011:       v   = vv + j;
2012:       sum = 0.0;
2013:       for (i=0; i<A->cmap->n; i++) {
2014:         sum += PetscAbsScalar(*v); v += mat->lda;
2015:       }
2016:       if (sum > *nrm) *nrm = sum;
2017:     }
2018:     PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
2019:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
2020:   MatDenseRestoreArrayRead(A,(const PetscScalar**)&vv);
2021:   return(0);
2022: }

2024: static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
2025: {
2026:   Mat_SeqDense   *aij = (Mat_SeqDense*)A->data;

2030:   switch (op) {
2031:   case MAT_ROW_ORIENTED:
2032:     aij->roworiented = flg;
2033:     break;
2034:   case MAT_NEW_NONZERO_LOCATIONS:
2035:   case MAT_NEW_NONZERO_LOCATION_ERR:
2036:   case MAT_NEW_NONZERO_ALLOCATION_ERR:
2037:   case MAT_FORCE_DIAGONAL_ENTRIES:
2038:   case MAT_KEEP_NONZERO_PATTERN:
2039:   case MAT_IGNORE_OFF_PROC_ENTRIES:
2040:   case MAT_USE_HASH_TABLE:
2041:   case MAT_IGNORE_ZERO_ENTRIES:
2042:   case MAT_IGNORE_LOWER_TRIANGULAR:
2043:   case MAT_SORTED_FULL:
2044:     PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
2045:     break;
2046:   case MAT_SPD:
2047:   case MAT_SYMMETRIC:
2048:   case MAT_STRUCTURALLY_SYMMETRIC:
2049:   case MAT_HERMITIAN:
2050:   case MAT_SYMMETRY_ETERNAL:
2051:     /* These options are handled directly by MatSetOption() */
2052:     break;
2053:   default:
2054:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]);
2055:   }
2056:   return(0);
2057: }

2059: PetscErrorCode MatZeroEntries_SeqDense(Mat A)
2060: {
2061:   Mat_SeqDense   *l = (Mat_SeqDense*)A->data;
2063:   PetscInt       lda=l->lda,m=A->rmap->n,n=A->cmap->n,j;
2064:   PetscScalar    *v;

2067:   MatDenseGetArrayWrite(A,&v);
2068:   if (lda>m) {
2069:     for (j=0; j<n; j++) {
2070:       PetscArrayzero(v+j*lda,m);
2071:     }
2072:   } else {
2073:     PetscArrayzero(v,PetscInt64Mult(m,n));
2074:   }
2075:   MatDenseRestoreArrayWrite(A,&v);
2076:   return(0);
2077: }

2079: static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
2080: {
2081:   PetscErrorCode    ierr;
2082:   Mat_SeqDense      *l = (Mat_SeqDense*)A->data;
2083:   PetscInt          m  = l->lda, n = A->cmap->n, i,j;
2084:   PetscScalar       *slot,*bb,*v;
2085:   const PetscScalar *xx;

2088:   if (PetscDefined(USE_DEBUG)) {
2089:     for (i=0; i<N; i++) {
2090:       if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
2091:       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);
2092:     }
2093:   }
2094:   if (!N) return(0);

2096:   /* fix right hand side if needed */
2097:   if (x && b) {
2098:     VecGetArrayRead(x,&xx);
2099:     VecGetArray(b,&bb);
2100:     for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
2101:     VecRestoreArrayRead(x,&xx);
2102:     VecRestoreArray(b,&bb);
2103:   }

2105:   MatDenseGetArray(A,&v);
2106:   for (i=0; i<N; i++) {
2107:     slot = v + rows[i];
2108:     for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
2109:   }
2110:   if (diag != 0.0) {
2111:     if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
2112:     for (i=0; i<N; i++) {
2113:       slot  = v + (m+1)*rows[i];
2114:       *slot = diag;
2115:     }
2116:   }
2117:   MatDenseRestoreArray(A,&v);
2118:   return(0);
2119: }

2121: static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A,PetscInt *lda)
2122: {
2123:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;

2126:   *lda = mat->lda;
2127:   return(0);
2128: }

2130: PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar **array)
2131: {
2132:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;

2135:   if (mat->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
2136:   *array = mat->v;
2137:   return(0);
2138: }

2140: PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar **array)
2141: {
2143:   *array = NULL;
2144:   return(0);
2145: }

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

2150:    Not collective

2152:    Input Parameter:
2153: .  mat - a MATSEQDENSE or MATMPIDENSE matrix

2155:    Output Parameter:
2156: .   lda - the leading dimension

2158:    Level: intermediate

2160: .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseSetLDA()
2161: @*/
2162: PetscErrorCode  MatDenseGetLDA(Mat A,PetscInt *lda)
2163: {

2169:   PetscUseMethod(A,"MatDenseGetLDA_C",(Mat,PetscInt*),(A,lda));
2170:   return(0);
2171: }

2173: /*@C
2174:    MatDenseSetLDA - Sets the leading dimension of the array used by the dense matrix

2176:    Not collective

2178:    Input Parameters:
2179: +  mat - a MATSEQDENSE or MATMPIDENSE matrix
2180: -  lda - the leading dimension

2182:    Level: intermediate

2184: .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetLDA()
2185: @*/
2186: PetscErrorCode  MatDenseSetLDA(Mat A,PetscInt lda)
2187: {

2192:   PetscTryMethod(A,"MatDenseSetLDA_C",(Mat,PetscInt),(A,lda));
2193:   return(0);
2194: }

2196: /*@C
2197:    MatDenseGetArray - gives read-write access to the array where the data for a dense matrix is stored

2199:    Logically Collective on Mat

2201:    Input Parameter:
2202: .  mat - a dense matrix

2204:    Output Parameter:
2205: .   array - pointer to the data

2207:    Level: intermediate

2209: .seealso: MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
2210: @*/
2211: PetscErrorCode  MatDenseGetArray(Mat A,PetscScalar **array)
2212: {

2218:   PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));
2219:   return(0);
2220: }

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

2225:    Logically Collective on Mat

2227:    Input Parameters:
2228: +  mat - a dense matrix
2229: -  array - pointer to the data

2231:    Level: intermediate

2233: .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
2234: @*/
2235: PetscErrorCode  MatDenseRestoreArray(Mat A,PetscScalar **array)
2236: {

2242:   PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));
2243:   PetscObjectStateIncrease((PetscObject)A);
2244: #if defined(PETSC_HAVE_CUDA)
2245:   A->offloadmask = PETSC_OFFLOAD_CPU;
2246: #endif
2247:   return(0);
2248: }

2250: /*@C
2251:    MatDenseGetArrayRead - gives read-only access to the array where the data for a dense matrix is stored

2253:    Not Collective

2255:    Input Parameter:
2256: .  mat - a dense matrix

2258:    Output Parameter:
2259: .   array - pointer to the data

2261:    Level: intermediate

2263: .seealso: MatDenseRestoreArrayRead(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
2264: @*/
2265: PetscErrorCode  MatDenseGetArrayRead(Mat A,const PetscScalar **array)
2266: {

2272:   PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array));
2273:   return(0);
2274: }

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

2279:    Not Collective

2281:    Input Parameters:
2282: +  mat - a dense matrix
2283: -  array - pointer to the data

2285:    Level: intermediate

2287: .seealso: MatDenseGetArrayRead(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
2288: @*/
2289: PetscErrorCode  MatDenseRestoreArrayRead(Mat A,const PetscScalar **array)
2290: {

2296:   PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array));
2297:   return(0);
2298: }

2300: /*@C
2301:    MatDenseGetArrayWrite - gives write-only access to the array where the data for a dense matrix is stored

2303:    Not Collective

2305:    Input Parameter:
2306: .  mat - a dense matrix

2308:    Output Parameter:
2309: .   array - pointer to the data

2311:    Level: intermediate

2313: .seealso: MatDenseRestoreArrayWrite(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead()
2314: @*/
2315: PetscErrorCode  MatDenseGetArrayWrite(Mat A,PetscScalar **array)
2316: {

2322:   PetscUseMethod(A,"MatDenseGetArrayWrite_C",(Mat,PetscScalar**),(A,array));
2323:   return(0);
2324: }

2326: /*@C
2327:    MatDenseRestoreArrayWrite - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArrayWrite()

2329:    Not Collective

2331:    Input Parameters:
2332: +  mat - a dense matrix
2333: -  array - pointer to the data

2335:    Level: intermediate

2337: .seealso: MatDenseGetArrayWrite(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead()
2338: @*/
2339: PetscErrorCode  MatDenseRestoreArrayWrite(Mat A,PetscScalar **array)
2340: {

2346:   PetscUseMethod(A,"MatDenseRestoreArrayWrite_C",(Mat,PetscScalar**),(A,array));
2347:   PetscObjectStateIncrease((PetscObject)A);
2348: #if defined(PETSC_HAVE_CUDA)
2349:   A->offloadmask = PETSC_OFFLOAD_CPU;
2350: #endif
2351:   return(0);
2352: }

2354: static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,MatReuse scall,Mat *B)
2355: {
2356:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
2358:   PetscInt       i,j,nrows,ncols,ldb;
2359:   const PetscInt *irow,*icol;
2360:   PetscScalar    *av,*bv,*v = mat->v;
2361:   Mat            newmat;

2364:   ISGetIndices(isrow,&irow);
2365:   ISGetIndices(iscol,&icol);
2366:   ISGetLocalSize(isrow,&nrows);
2367:   ISGetLocalSize(iscol,&ncols);

2369:   /* Check submatrixcall */
2370:   if (scall == MAT_REUSE_MATRIX) {
2371:     PetscInt n_cols,n_rows;
2372:     MatGetSize(*B,&n_rows,&n_cols);
2373:     if (n_rows != nrows || n_cols != ncols) {
2374:       /* resize the result matrix to match number of requested rows/columns */
2375:       MatSetSizes(*B,nrows,ncols,nrows,ncols);
2376:     }
2377:     newmat = *B;
2378:   } else {
2379:     /* Create and fill new matrix */
2380:     MatCreate(PetscObjectComm((PetscObject)A),&newmat);
2381:     MatSetSizes(newmat,nrows,ncols,nrows,ncols);
2382:     MatSetType(newmat,((PetscObject)A)->type_name);
2383:     MatSeqDenseSetPreallocation(newmat,NULL);
2384:   }

2386:   /* Now extract the data pointers and do the copy,column at a time */
2387:   MatDenseGetArray(newmat,&bv);
2388:   MatDenseGetLDA(newmat,&ldb);
2389:   for (i=0; i<ncols; i++) {
2390:     av = v + mat->lda*icol[i];
2391:     for (j=0; j<nrows; j++) bv[j] = av[irow[j]];
2392:     bv += ldb;
2393:   }
2394:   MatDenseRestoreArray(newmat,&bv);

2396:   /* Assemble the matrices so that the correct flags are set */
2397:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
2398:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);

2400:   /* Free work space */
2401:   ISRestoreIndices(isrow,&irow);
2402:   ISRestoreIndices(iscol,&icol);
2403:   *B   = newmat;
2404:   return(0);
2405: }

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

2413:   if (scall == MAT_INITIAL_MATRIX) {
2414:     PetscCalloc1(n,B);
2415:   }

2417:   for (i=0; i<n; i++) {
2418:     MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],scall,&(*B)[i]);
2419:   }
2420:   return(0);
2421: }

2423: static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
2424: {
2426:   return(0);
2427: }

2429: static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
2430: {
2432:   return(0);
2433: }

2435: PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
2436: {
2437:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data;
2438:   PetscErrorCode    ierr;
2439:   const PetscScalar *va;
2440:   PetscScalar       *vb;
2441:   PetscInt          lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j;

2444:   /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
2445:   if (A->ops->copy != B->ops->copy) {
2446:     MatCopy_Basic(A,B,str);
2447:     return(0);
2448:   }
2449:   if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)");
2450:   MatDenseGetArrayRead(A,&va);
2451:   MatDenseGetArray(B,&vb);
2452:   if (lda1>m || lda2>m) {
2453:     for (j=0; j<n; j++) {
2454:       PetscArraycpy(vb+j*lda2,va+j*lda1,m);
2455:     }
2456:   } else {
2457:     PetscArraycpy(vb,va,A->rmap->n*A->cmap->n);
2458:   }
2459:   MatDenseRestoreArray(B,&vb);
2460:   MatDenseRestoreArrayRead(A,&va);
2461:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
2462:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
2463:   return(0);
2464: }

2466: static PetscErrorCode MatSetUp_SeqDense(Mat A)
2467: {

2471:   PetscLayoutSetUp(A->rmap);
2472:   PetscLayoutSetUp(A->cmap);
2473:   if (!A->preallocated) {
2474:     MatSeqDenseSetPreallocation(A,NULL);
2475:   }
2476:   return(0);
2477: }

2479: static PetscErrorCode MatConjugate_SeqDense(Mat A)
2480: {
2481:   Mat_SeqDense   *mat = (Mat_SeqDense *) A->data;
2482:   PetscInt       i,nz = A->rmap->n*A->cmap->n;
2483:   PetscInt       min = PetscMin(A->rmap->n,A->cmap->n);
2484:   PetscScalar    *aa;

2488:   MatDenseGetArray(A,&aa);
2489:   for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
2490:   MatDenseRestoreArray(A,&aa);
2491:   if (mat->tau) for (i = 0; i < min; i++) mat->tau[i] = PetscConj(mat->tau[i]);
2492:   return(0);
2493: }

2495: static PetscErrorCode MatRealPart_SeqDense(Mat A)
2496: {
2497:   PetscInt       i,nz = A->rmap->n*A->cmap->n;
2498:   PetscScalar    *aa;

2502:   MatDenseGetArray(A,&aa);
2503:   for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
2504:   MatDenseRestoreArray(A,&aa);
2505:   return(0);
2506: }

2508: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
2509: {
2510:   PetscInt       i,nz = A->rmap->n*A->cmap->n;
2511:   PetscScalar    *aa;

2515:   MatDenseGetArray(A,&aa);
2516:   for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
2517:   MatDenseRestoreArray(A,&aa);
2518:   return(0);
2519: }

2521: /* ----------------------------------------------------------------*/
2522: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2523: {
2525:   PetscInt       m=A->rmap->n,n=B->cmap->n;
2526:   PetscBool      cisdense;

2529:   MatSetSizes(C,m,n,m,n);
2530:   PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
2531:   if (!cisdense) {
2532:     PetscBool flg;

2534:     PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2535:     MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
2536:   }
2537:   MatSetUp(C);
2538:   return(0);
2539: }

2541: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2542: {
2543:   Mat_SeqDense       *a=(Mat_SeqDense*)A->data,*b=(Mat_SeqDense*)B->data,*c=(Mat_SeqDense*)C->data;
2544:   PetscBLASInt       m,n,k;
2545:   const PetscScalar *av,*bv;
2546:   PetscScalar       *cv;
2547:   PetscScalar       _DOne=1.0,_DZero=0.0;
2548:   PetscErrorCode    ierr;

2551:   PetscBLASIntCast(C->rmap->n,&m);
2552:   PetscBLASIntCast(C->cmap->n,&n);
2553:   PetscBLASIntCast(A->cmap->n,&k);
2554:   if (!m || !n || !k) return(0);
2555:   MatDenseGetArrayRead(A,&av);
2556:   MatDenseGetArrayRead(B,&bv);
2557:   MatDenseGetArrayWrite(C,&cv);
2558:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2559:   PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2560:   MatDenseRestoreArrayRead(A,&av);
2561:   MatDenseRestoreArrayRead(B,&bv);
2562:   MatDenseRestoreArrayWrite(C,&cv);
2563:   return(0);
2564: }

2566: PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2567: {
2569:   PetscInt       m=A->rmap->n,n=B->rmap->n;
2570:   PetscBool      cisdense;

2573:   MatSetSizes(C,m,n,m,n);
2574:   PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
2575:   if (!cisdense) {
2576:     PetscBool flg;

2578:     PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2579:     MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
2580:   }
2581:   MatSetUp(C);
2582:   return(0);
2583: }

2585: PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2586: {
2587:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
2588:   Mat_SeqDense      *b = (Mat_SeqDense*)B->data;
2589:   Mat_SeqDense      *c = (Mat_SeqDense*)C->data;
2590:   const PetscScalar *av,*bv;
2591:   PetscScalar       *cv;
2592:   PetscBLASInt      m,n,k;
2593:   PetscScalar       _DOne=1.0,_DZero=0.0;
2594:   PetscErrorCode    ierr;

2597:   PetscBLASIntCast(C->rmap->n,&m);
2598:   PetscBLASIntCast(C->cmap->n,&n);
2599:   PetscBLASIntCast(A->cmap->n,&k);
2600:   if (!m || !n || !k) return(0);
2601:   MatDenseGetArrayRead(A,&av);
2602:   MatDenseGetArrayRead(B,&bv);
2603:   MatDenseGetArrayWrite(C,&cv);
2604:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2605:   MatDenseRestoreArrayRead(A,&av);
2606:   MatDenseRestoreArrayRead(B,&bv);
2607:   MatDenseRestoreArrayWrite(C,&cv);
2608:   PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2609:   return(0);
2610: }

2612: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2613: {
2615:   PetscInt       m=A->cmap->n,n=B->cmap->n;
2616:   PetscBool      cisdense;

2619:   MatSetSizes(C,m,n,m,n);
2620:   PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
2621:   if (!cisdense) {
2622:     PetscBool flg;

2624:     PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2625:     MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
2626:   }
2627:   MatSetUp(C);
2628:   return(0);
2629: }

2631: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2632: {
2633:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
2634:   Mat_SeqDense      *b = (Mat_SeqDense*)B->data;
2635:   Mat_SeqDense      *c = (Mat_SeqDense*)C->data;
2636:   const PetscScalar *av,*bv;
2637:   PetscScalar       *cv;
2638:   PetscBLASInt      m,n,k;
2639:   PetscScalar       _DOne=1.0,_DZero=0.0;
2640:   PetscErrorCode    ierr;

2643:   PetscBLASIntCast(C->rmap->n,&m);
2644:   PetscBLASIntCast(C->cmap->n,&n);
2645:   PetscBLASIntCast(A->rmap->n,&k);
2646:   if (!m || !n || !k) return(0);
2647:   MatDenseGetArrayRead(A,&av);
2648:   MatDenseGetArrayRead(B,&bv);
2649:   MatDenseGetArrayWrite(C,&cv);
2650:   PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2651:   MatDenseRestoreArrayRead(A,&av);
2652:   MatDenseRestoreArrayRead(B,&bv);
2653:   MatDenseRestoreArrayWrite(C,&cv);
2654:   PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2655:   return(0);
2656: }

2658: /* ----------------------------------------------- */
2659: static PetscErrorCode MatProductSetFromOptions_SeqDense_AB(Mat C)
2660: {
2662:   C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqDense;
2663:   C->ops->productsymbolic = MatProductSymbolic_AB;
2664:   return(0);
2665: }

2667: static PetscErrorCode MatProductSetFromOptions_SeqDense_AtB(Mat C)
2668: {
2670:   C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_SeqDense_SeqDense;
2671:   C->ops->productsymbolic          = MatProductSymbolic_AtB;
2672:   return(0);
2673: }

2675: static PetscErrorCode MatProductSetFromOptions_SeqDense_ABt(Mat C)
2676: {
2678:   C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqDense_SeqDense;
2679:   C->ops->productsymbolic          = MatProductSymbolic_ABt;
2680:   return(0);
2681: }

2683: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense(Mat C)
2684: {
2686:   Mat_Product    *product = C->product;

2689:   switch (product->type) {
2690:   case MATPRODUCT_AB:
2691:     MatProductSetFromOptions_SeqDense_AB(C);
2692:     break;
2693:   case MATPRODUCT_AtB:
2694:     MatProductSetFromOptions_SeqDense_AtB(C);
2695:     break;
2696:   case MATPRODUCT_ABt:
2697:     MatProductSetFromOptions_SeqDense_ABt(C);
2698:     break;
2699:   default:
2700:     break;
2701:   }
2702:   return(0);
2703: }
2704: /* ----------------------------------------------- */

2706: static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[])
2707: {
2708:   Mat_SeqDense       *a = (Mat_SeqDense*)A->data;
2709:   PetscErrorCode     ierr;
2710:   PetscInt           i,j,m = A->rmap->n,n = A->cmap->n,p;
2711:   PetscScalar        *x;
2712:   const PetscScalar *aa;

2715:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2716:   VecGetArray(v,&x);
2717:   VecGetLocalSize(v,&p);
2718:   MatDenseGetArrayRead(A,&aa);
2719:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2720:   for (i=0; i<m; i++) {
2721:     x[i] = aa[i]; if (idx) idx[i] = 0;
2722:     for (j=1; j<n; j++) {
2723:       if (PetscRealPart(x[i]) < PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;}
2724:     }
2725:   }
2726:   MatDenseRestoreArrayRead(A,&aa);
2727:   VecRestoreArray(v,&x);
2728:   return(0);
2729: }

2731: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])
2732: {
2733:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
2734:   PetscErrorCode    ierr;
2735:   PetscInt          i,j,m = A->rmap->n,n = A->cmap->n,p;
2736:   PetscScalar       *x;
2737:   PetscReal         atmp;
2738:   const PetscScalar *aa;

2741:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2742:   VecGetArray(v,&x);
2743:   VecGetLocalSize(v,&p);
2744:   MatDenseGetArrayRead(A,&aa);
2745:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2746:   for (i=0; i<m; i++) {
2747:     x[i] = PetscAbsScalar(aa[i]);
2748:     for (j=1; j<n; j++) {
2749:       atmp = PetscAbsScalar(aa[i+a->lda*j]);
2750:       if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;}
2751:     }
2752:   }
2753:   MatDenseRestoreArrayRead(A,&aa);
2754:   VecRestoreArray(v,&x);
2755:   return(0);
2756: }

2758: static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])
2759: {
2760:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
2761:   PetscErrorCode    ierr;
2762:   PetscInt          i,j,m = A->rmap->n,n = A->cmap->n,p;
2763:   PetscScalar       *x;
2764:   const PetscScalar *aa;

2767:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2768:   MatDenseGetArrayRead(A,&aa);
2769:   VecGetArray(v,&x);
2770:   VecGetLocalSize(v,&p);
2771:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2772:   for (i=0; i<m; i++) {
2773:     x[i] = aa[i]; if (idx) idx[i] = 0;
2774:     for (j=1; j<n; j++) {
2775:       if (PetscRealPart(x[i]) > PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;}
2776:     }
2777:   }
2778:   VecRestoreArray(v,&x);
2779:   MatDenseRestoreArrayRead(A,&aa);
2780:   return(0);
2781: }

2783: PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
2784: {
2785:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
2786:   PetscErrorCode    ierr;
2787:   PetscScalar       *x;
2788:   const PetscScalar *aa;

2791:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2792:   MatDenseGetArrayRead(A,&aa);
2793:   VecGetArray(v,&x);
2794:   PetscArraycpy(x,aa+col*a->lda,A->rmap->n);
2795:   VecRestoreArray(v,&x);
2796:   MatDenseRestoreArrayRead(A,&aa);
2797:   return(0);
2798: }

2800: PETSC_INTERN PetscErrorCode MatGetColumnReductions_SeqDense(Mat A,PetscInt type,PetscReal *reductions)
2801: {
2802:   PetscErrorCode    ierr;
2803:   PetscInt          i,j,m,n;
2804:   const PetscScalar *a;

2807:   MatGetSize(A,&m,&n);
2808:   PetscArrayzero(reductions,n);
2809:   MatDenseGetArrayRead(A,&a);
2810:   if (type == NORM_2) {
2811:     for (i=0; i<n; i++) {
2812:       for (j=0; j<m; j++) {
2813:         reductions[i] += PetscAbsScalar(a[j]*a[j]);
2814:       }
2815:       a += m;
2816:     }
2817:   } else if (type == NORM_1) {
2818:     for (i=0; i<n; i++) {
2819:       for (j=0; j<m; j++) {
2820:         reductions[i] += PetscAbsScalar(a[j]);
2821:       }
2822:       a += m;
2823:     }
2824:   } else if (type == NORM_INFINITY) {
2825:     for (i=0; i<n; i++) {
2826:       for (j=0; j<m; j++) {
2827:         reductions[i] = PetscMax(PetscAbsScalar(a[j]),reductions[i]);
2828:       }
2829:       a += m;
2830:     }
2831:   } else if (type == REDUCTION_SUM_REALPART || type == REDUCTION_MEAN_REALPART) {
2832:     for (i=0; i<n; i++) {
2833:       for (j=0; j<m; j++) {
2834:         reductions[i] += PetscRealPart(a[j]);
2835:       }
2836:       a += m;
2837:     }
2838:   } else if (type == REDUCTION_SUM_IMAGINARYPART || type == REDUCTION_MEAN_IMAGINARYPART) {
2839:     for (i=0; i<n; i++) {
2840:       for (j=0; j<m; j++) {
2841:         reductions[i] += PetscImaginaryPart(a[j]);
2842:       }
2843:       a += m;
2844:     }
2845:   } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown reduction type");
2846:   MatDenseRestoreArrayRead(A,&a);
2847:   if (type == NORM_2) {
2848:     for (i=0; i<n; i++) reductions[i] = PetscSqrtReal(reductions[i]);
2849:   } else if (type == REDUCTION_MEAN_REALPART || type == REDUCTION_MEAN_IMAGINARYPART) {
2850:     for (i=0; i<n; i++) reductions[i] /= m;
2851:   }
2852:   return(0);
2853: }

2855: static PetscErrorCode  MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
2856: {
2858:   PetscScalar    *a;
2859:   PetscInt       lda,m,n,i,j;

2862:   MatGetSize(x,&m,&n);
2863:   MatDenseGetLDA(x,&lda);
2864:   MatDenseGetArray(x,&a);
2865:   for (j=0; j<n; j++) {
2866:     for (i=0; i<m; i++) {
2867:       PetscRandomGetValue(rctx,a+j*lda+i);
2868:     }
2869:   }
2870:   MatDenseRestoreArray(x,&a);
2871:   return(0);
2872: }

2874: static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool  *missing,PetscInt *d)
2875: {
2877:   *missing = PETSC_FALSE;
2878:   return(0);
2879: }

2881: /* vals is not const */
2882: static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals)
2883: {
2885:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2886:   PetscScalar    *v;

2889:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2890:   MatDenseGetArray(A,&v);
2891:   *vals = v+col*a->lda;
2892:   MatDenseRestoreArray(A,&v);
2893:   return(0);
2894: }

2896: static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals)
2897: {
2899:   *vals = NULL; /* user cannot accidentally use the array later */
2900:   return(0);
2901: }

2903: /* -------------------------------------------------------------------*/
2904: static struct _MatOps MatOps_Values = { MatSetValues_SeqDense,
2905:                                         MatGetRow_SeqDense,
2906:                                         MatRestoreRow_SeqDense,
2907:                                         MatMult_SeqDense,
2908:                                 /*  4*/ MatMultAdd_SeqDense,
2909:                                         MatMultTranspose_SeqDense,
2910:                                         MatMultTransposeAdd_SeqDense,
2911:                                         NULL,
2912:                                         NULL,
2913:                                         NULL,
2914:                                 /* 10*/ NULL,
2915:                                         MatLUFactor_SeqDense,
2916:                                         MatCholeskyFactor_SeqDense,
2917:                                         MatSOR_SeqDense,
2918:                                         MatTranspose_SeqDense,
2919:                                 /* 15*/ MatGetInfo_SeqDense,
2920:                                         MatEqual_SeqDense,
2921:                                         MatGetDiagonal_SeqDense,
2922:                                         MatDiagonalScale_SeqDense,
2923:                                         MatNorm_SeqDense,
2924:                                 /* 20*/ MatAssemblyBegin_SeqDense,
2925:                                         MatAssemblyEnd_SeqDense,
2926:                                         MatSetOption_SeqDense,
2927:                                         MatZeroEntries_SeqDense,
2928:                                 /* 24*/ MatZeroRows_SeqDense,
2929:                                         NULL,
2930:                                         NULL,
2931:                                         NULL,
2932:                                         NULL,
2933:                                 /* 29*/ MatSetUp_SeqDense,
2934:                                         NULL,
2935:                                         NULL,
2936:                                         NULL,
2937:                                         NULL,
2938:                                 /* 34*/ MatDuplicate_SeqDense,
2939:                                         NULL,
2940:                                         NULL,
2941:                                         NULL,
2942:                                         NULL,
2943:                                 /* 39*/ MatAXPY_SeqDense,
2944:                                         MatCreateSubMatrices_SeqDense,
2945:                                         NULL,
2946:                                         MatGetValues_SeqDense,
2947:                                         MatCopy_SeqDense,
2948:                                 /* 44*/ MatGetRowMax_SeqDense,
2949:                                         MatScale_SeqDense,
2950:                                         MatShift_Basic,
2951:                                         NULL,
2952:                                         MatZeroRowsColumns_SeqDense,
2953:                                 /* 49*/ MatSetRandom_SeqDense,
2954:                                         NULL,
2955:                                         NULL,
2956:                                         NULL,
2957:                                         NULL,
2958:                                 /* 54*/ NULL,
2959:                                         NULL,
2960:                                         NULL,
2961:                                         NULL,
2962:                                         NULL,
2963:                                 /* 59*/ MatCreateSubMatrix_SeqDense,
2964:                                         MatDestroy_SeqDense,
2965:                                         MatView_SeqDense,
2966:                                         NULL,
2967:                                         NULL,
2968:                                 /* 64*/ NULL,
2969:                                         NULL,
2970:                                         NULL,
2971:                                         NULL,
2972:                                         NULL,
2973:                                 /* 69*/ MatGetRowMaxAbs_SeqDense,
2974:                                         NULL,
2975:                                         NULL,
2976:                                         NULL,
2977:                                         NULL,
2978:                                 /* 74*/ NULL,
2979:                                         NULL,
2980:                                         NULL,
2981:                                         NULL,
2982:                                         NULL,
2983:                                 /* 79*/ NULL,
2984:                                         NULL,
2985:                                         NULL,
2986:                                         NULL,
2987:                                 /* 83*/ MatLoad_SeqDense,
2988:                                         MatIsSymmetric_SeqDense,
2989:                                         MatIsHermitian_SeqDense,
2990:                                         NULL,
2991:                                         NULL,
2992:                                         NULL,
2993:                                 /* 89*/ NULL,
2994:                                         NULL,
2995:                                         MatMatMultNumeric_SeqDense_SeqDense,
2996:                                         NULL,
2997:                                         NULL,
2998:                                 /* 94*/ NULL,
2999:                                         NULL,
3000:                                         NULL,
3001:                                         MatMatTransposeMultNumeric_SeqDense_SeqDense,
3002:                                         NULL,
3003:                                 /* 99*/ MatProductSetFromOptions_SeqDense,
3004:                                         NULL,
3005:                                         NULL,
3006:                                         MatConjugate_SeqDense,
3007:                                         NULL,
3008:                                 /*104*/ NULL,
3009:                                         MatRealPart_SeqDense,
3010:                                         MatImaginaryPart_SeqDense,
3011:                                         NULL,
3012:                                         NULL,
3013:                                 /*109*/ NULL,
3014:                                         NULL,
3015:                                         MatGetRowMin_SeqDense,
3016:                                         MatGetColumnVector_SeqDense,
3017:                                         MatMissingDiagonal_SeqDense,
3018:                                 /*114*/ NULL,
3019:                                         NULL,
3020:                                         NULL,
3021:                                         NULL,
3022:                                         NULL,
3023:                                 /*119*/ NULL,
3024:                                         NULL,
3025:                                         NULL,
3026:                                         NULL,
3027:                                         NULL,
3028:                                 /*124*/ NULL,
3029:                                         MatGetColumnReductions_SeqDense,
3030:                                         NULL,
3031:                                         NULL,
3032:                                         NULL,
3033:                                 /*129*/ NULL,
3034:                                         NULL,
3035:                                         NULL,
3036:                                         MatTransposeMatMultNumeric_SeqDense_SeqDense,
3037:                                         NULL,
3038:                                 /*134*/ NULL,
3039:                                         NULL,
3040:                                         NULL,
3041:                                         NULL,
3042:                                         NULL,
3043:                                 /*139*/ NULL,
3044:                                         NULL,
3045:                                         NULL,
3046:                                         NULL,
3047:                                         NULL,
3048:                                         MatCreateMPIMatConcatenateSeqMat_SeqDense,
3049:                                 /*145*/ NULL,
3050:                                         NULL,
3051:                                         NULL
3052: };

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

3059:    Collective

3061:    Input Parameters:
3062: +  comm - MPI communicator, set to PETSC_COMM_SELF
3063: .  m - number of rows
3064: .  n - number of columns
3065: -  data - optional location of matrix data in column major order.  Set data=NULL for PETSc
3066:    to control all matrix memory allocation.

3068:    Output Parameter:
3069: .  A - the matrix

3071:    Notes:
3072:    The data input variable is intended primarily for Fortran programmers
3073:    who wish to allocate their own matrix memory space.  Most users should
3074:    set data=NULL.

3076:    Level: intermediate

3078: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
3079: @*/
3080: PetscErrorCode  MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
3081: {

3085:   MatCreate(comm,A);
3086:   MatSetSizes(*A,m,n,m,n);
3087:   MatSetType(*A,MATSEQDENSE);
3088:   MatSeqDenseSetPreallocation(*A,data);
3089:   return(0);
3090: }

3092: /*@C
3093:    MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements

3095:    Collective

3097:    Input Parameters:
3098: +  B - the matrix
3099: -  data - the array (or NULL)

3101:    Notes:
3102:    The data input variable is intended primarily for Fortran programmers
3103:    who wish to allocate their own matrix memory space.  Most users should
3104:    need not call this routine.

3106:    Level: intermediate

3108: .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatDenseSetLDA()

3110: @*/
3111: PetscErrorCode  MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
3112: {

3117:   PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
3118:   return(0);
3119: }

3121: PetscErrorCode  MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
3122: {
3123:   Mat_SeqDense   *b = (Mat_SeqDense*)B->data;

3127:   if (b->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
3128:   B->preallocated = PETSC_TRUE;

3130:   PetscLayoutSetUp(B->rmap);
3131:   PetscLayoutSetUp(B->cmap);

3133:   if (b->lda <= 0) b->lda = B->rmap->n;

3135:   PetscIntMultError(b->lda,B->cmap->n,NULL);
3136:   if (!data) { /* petsc-allocated storage */
3137:     if (!b->user_alloc) { PetscFree(b->v); }
3138:     PetscCalloc1((size_t)b->lda*B->cmap->n,&b->v);
3139:     PetscLogObjectMemory((PetscObject)B,b->lda*B->cmap->n*sizeof(PetscScalar));

3141:     b->user_alloc = PETSC_FALSE;
3142:   } else { /* user-allocated storage */
3143:     if (!b->user_alloc) { PetscFree(b->v); }
3144:     b->v          = data;
3145:     b->user_alloc = PETSC_TRUE;
3146:   }
3147:   B->assembled = PETSC_TRUE;
3148:   return(0);
3149: }

3151: #if defined(PETSC_HAVE_ELEMENTAL)
3152: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
3153: {
3154:   Mat               mat_elemental;
3155:   PetscErrorCode    ierr;
3156:   const PetscScalar *array;
3157:   PetscScalar       *v_colwise;
3158:   PetscInt          M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols;

3161:   PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);
3162:   MatDenseGetArrayRead(A,&array);
3163:   /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
3164:   k = 0;
3165:   for (j=0; j<N; j++) {
3166:     cols[j] = j;
3167:     for (i=0; i<M; i++) {
3168:       v_colwise[j*M+i] = array[k++];
3169:     }
3170:   }
3171:   for (i=0; i<M; i++) {
3172:     rows[i] = i;
3173:   }
3174:   MatDenseRestoreArrayRead(A,&array);

3176:   MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
3177:   MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);
3178:   MatSetType(mat_elemental,MATELEMENTAL);
3179:   MatSetUp(mat_elemental);

3181:   /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
3182:   MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);
3183:   MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
3184:   MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);
3185:   PetscFree3(v_colwise,rows,cols);

3187:   if (reuse == MAT_INPLACE_MATRIX) {
3188:     MatHeaderReplace(A,&mat_elemental);
3189:   } else {
3190:     *newmat = mat_elemental;
3191:   }
3192:   return(0);
3193: }
3194: #endif

3196: PetscErrorCode  MatDenseSetLDA_SeqDense(Mat B,PetscInt lda)
3197: {
3198:   Mat_SeqDense *b = (Mat_SeqDense*)B->data;
3199:   PetscBool    data;

3202:   data = (PetscBool)((B->rmap->n > 0 && B->cmap->n > 0) ? (b->v ? PETSC_TRUE : PETSC_FALSE) : PETSC_FALSE);
3203:   if (!b->user_alloc && data && b->lda!=lda) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"LDA cannot be changed after allocation of internal storage");
3204:   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);
3205:   b->lda = lda;
3206:   return(0);
3207: }

3209: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat)
3210: {
3212:   PetscMPIInt    size;

3215:   MPI_Comm_size(comm,&size);
3216:   if (size == 1) {
3217:     if (scall == MAT_INITIAL_MATRIX) {
3218:       MatDuplicate(inmat,MAT_COPY_VALUES,outmat);
3219:     } else {
3220:       MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);
3221:     }
3222:   } else {
3223:     MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);
3224:   }
3225:   return(0);
3226: }

3228: PetscErrorCode MatDenseGetColumnVec_SeqDense(Mat A,PetscInt col,Vec *v)
3229: {
3230:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

3234:   if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
3235:   if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
3236:   if (!a->cvec) {
3237:     VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);
3238:     PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);
3239:   }
3240:   a->vecinuse = col + 1;
3241:   MatDenseGetArray(A,(PetscScalar**)&a->ptrinuse);
3242:   VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);
3243:   *v   = a->cvec;
3244:   return(0);
3245: }

3247: PetscErrorCode MatDenseRestoreColumnVec_SeqDense(Mat A,PetscInt col,Vec *v)
3248: {
3249:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

3253:   if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first");
3254:   if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector");
3255:   a->vecinuse = 0;
3256:   MatDenseRestoreArray(A,(PetscScalar**)&a->ptrinuse);
3257:   VecResetArray(a->cvec);
3258:   *v   = NULL;
3259:   return(0);
3260: }

3262: PetscErrorCode MatDenseGetColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v)
3263: {
3264:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

3268:   if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
3269:   if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
3270:   if (!a->cvec) {
3271:     VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);
3272:     PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);
3273:   }
3274:   a->vecinuse = col + 1;
3275:   MatDenseGetArrayRead(A,&a->ptrinuse);
3276:   VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);
3277:   VecLockReadPush(a->cvec);
3278:   *v   = a->cvec;
3279:   return(0);
3280: }

3282: PetscErrorCode MatDenseRestoreColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v)
3283: {
3284:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

3288:   if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first");
3289:   if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector");
3290:   a->vecinuse = 0;
3291:   MatDenseRestoreArrayRead(A,&a->ptrinuse);
3292:   VecLockReadPop(a->cvec);
3293:   VecResetArray(a->cvec);
3294:   *v   = NULL;
3295:   return(0);
3296: }

3298: PetscErrorCode MatDenseGetColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v)
3299: {
3300:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

3304:   if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
3305:   if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
3306:   if (!a->cvec) {
3307:     VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);
3308:     PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);
3309:   }
3310:   a->vecinuse = col + 1;
3311:   MatDenseGetArrayWrite(A,(PetscScalar**)&a->ptrinuse);
3312:   VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);
3313:   *v   = a->cvec;
3314:   return(0);
3315: }

3317: PetscErrorCode MatDenseRestoreColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v)
3318: {
3319:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

3323:   if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first");
3324:   if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector");
3325:   a->vecinuse = 0;
3326:   MatDenseRestoreArrayWrite(A,(PetscScalar**)&a->ptrinuse);
3327:   VecResetArray(a->cvec);
3328:   *v   = NULL;
3329:   return(0);
3330: }

3332: PetscErrorCode MatDenseGetSubMatrix_SeqDense(Mat A,PetscInt cbegin,PetscInt cend,Mat *v)
3333: {
3334:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

3338:   if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
3339:   if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
3340:   if (a->cmat && cend-cbegin != a->cmat->cmap->N) {
3341:     MatDestroy(&a->cmat);
3342:   }
3343:   if (!a->cmat) {
3344:     MatCreateDense(PetscObjectComm((PetscObject)A),A->rmap->n,PETSC_DECIDE,A->rmap->N,cend-cbegin,(PetscScalar*)a->v + (size_t)cbegin * (size_t)a->lda,&a->cmat);
3345:     PetscLogObjectParent((PetscObject)A,(PetscObject)a->cmat);
3346:   } else {
3347:     MatDensePlaceArray(a->cmat,a->v + (size_t)cbegin * (size_t)a->lda);
3348:   }
3349:   MatDenseSetLDA(a->cmat,a->lda);
3350:   a->matinuse = cbegin + 1;
3351:   *v = a->cmat;
3352:   return(0);
3353: }

3355: PetscErrorCode MatDenseRestoreSubMatrix_SeqDense(Mat A,Mat *v)
3356: {
3357:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

3361:   if (!a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetSubMatrix() first");
3362:   if (!a->cmat) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column matrix");
3363:   if (*v != a->cmat) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not the matrix obtained from MatDenseGetSubMatrix()");
3364:   a->matinuse = 0;
3365:   MatDenseResetArray(a->cmat);
3366:   *v   = NULL;
3367:   return(0);
3368: }

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

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

3376:   Level: beginner

3378: .seealso: MatCreateSeqDense()

3380: M*/
3381: PetscErrorCode MatCreate_SeqDense(Mat B)
3382: {
3383:   Mat_SeqDense   *b;
3385:   PetscMPIInt    size;

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

3391:   PetscNewLog(B,&b);
3392:   PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
3393:   B->data = (void*)b;

3395:   b->roworiented = PETSC_TRUE;

3397:   PetscObjectComposeFunction((PetscObject)B,"MatQRFactor_C",MatQRFactor_SeqDense);
3398:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetLDA_C",MatDenseGetLDA_SeqDense);
3399:   PetscObjectComposeFunction((PetscObject)B,"MatDenseSetLDA_C",MatDenseSetLDA_SeqDense);
3400:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);
3401:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);
3402:   PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);
3403:   PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);
3404:   PetscObjectComposeFunction((PetscObject)B,"MatDenseReplaceArray_C",MatDenseReplaceArray_SeqDense);
3405:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);
3406:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense);
3407:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayWrite_C",MatDenseGetArray_SeqDense);
3408:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayWrite_C",MatDenseRestoreArray_SeqDense);
3409:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);
3410: #if defined(PETSC_HAVE_ELEMENTAL)
3411:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);
3412: #endif
3413: #if defined(PETSC_HAVE_SCALAPACK)
3414:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_scalapack_C",MatConvert_Dense_ScaLAPACK);
3415: #endif
3416: #if defined(PETSC_HAVE_CUDA)
3417:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqdensecuda_C",MatConvert_SeqDense_SeqDenseCUDA);
3418:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",MatProductSetFromOptions_SeqDense);
3419:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdense_C",MatProductSetFromOptions_SeqDense);
3420:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdensecuda_C",MatProductSetFromOptions_SeqDense);
3421: #endif
3422:   PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);
3423:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqdense_C",MatProductSetFromOptions_SeqAIJ_SeqDense);
3424:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdense_C",MatProductSetFromOptions_SeqDense);
3425:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);
3426:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqsbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);

3428:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);
3429:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);
3430:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVec_C",MatDenseGetColumnVec_SeqDense);
3431:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVec_C",MatDenseRestoreColumnVec_SeqDense);
3432:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecRead_C",MatDenseGetColumnVecRead_SeqDense);
3433:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecRead_C",MatDenseRestoreColumnVecRead_SeqDense);
3434:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecWrite_C",MatDenseGetColumnVecWrite_SeqDense);
3435:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecWrite_C",MatDenseRestoreColumnVecWrite_SeqDense);
3436:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetSubMatrix_C",MatDenseGetSubMatrix_SeqDense);
3437:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreSubMatrix_C",MatDenseRestoreSubMatrix_SeqDense);
3438:   PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
3439:   return(0);
3440: }

3442: /*@C
3443:    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.

3445:    Not Collective

3447:    Input Parameters:
3448: +  mat - a MATSEQDENSE or MATMPIDENSE matrix
3449: -  col - column index

3451:    Output Parameter:
3452: .  vals - pointer to the data

3454:    Level: intermediate

3456: .seealso: MatDenseRestoreColumn()
3457: @*/
3458: PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals)
3459: {

3466:   PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));
3467:   return(0);
3468: }

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

3473:    Not Collective

3475:    Input Parameter:
3476: .  mat - a MATSEQDENSE or MATMPIDENSE matrix

3478:    Output Parameter:
3479: .  vals - pointer to the data

3481:    Level: intermediate

3483: .seealso: MatDenseGetColumn()
3484: @*/
3485: PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals)
3486: {

3492:   PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));
3493:   return(0);
3494: }

3496: /*@C
3497:    MatDenseGetColumnVec - Gives read-write access to a column of a dense matrix, represented as a Vec.

3499:    Collective

3501:    Input Parameters:
3502: +  mat - the Mat object
3503: -  col - the column index

3505:    Output Parameter:
3506: .  v - the vector

3508:    Notes:
3509:      The vector is owned by PETSc. Users need to call MatDenseRestoreColumnVec() when the vector is no longer needed.
3510:      Use MatDenseGetColumnVecRead() to obtain read-only access or MatDenseGetColumnVecWrite() for write-only access.

3512:    Level: intermediate

3514: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3515: @*/
3516: PetscErrorCode MatDenseGetColumnVec(Mat A,PetscInt col,Vec *v)
3517: {

3525:   if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3526:   if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N);
3527:   PetscUseMethod(A,"MatDenseGetColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));
3528:   return(0);
3529: }

3531: /*@C
3532:    MatDenseRestoreColumnVec - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVec().

3534:    Collective

3536:    Input Parameters:
3537: +  mat - the Mat object
3538: .  col - the column index
3539: -  v - the Vec object

3541:    Level: intermediate

3543: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3544: @*/
3545: PetscErrorCode MatDenseRestoreColumnVec(Mat A,PetscInt col,Vec *v)
3546: {

3554:   if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3555:   if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N);
3556:   PetscUseMethod(A,"MatDenseRestoreColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));
3557:   return(0);
3558: }

3560: /*@C
3561:    MatDenseGetColumnVecRead - Gives read-only access to a column of a dense matrix, represented as a Vec.

3563:    Collective

3565:    Input Parameters:
3566: +  mat - the Mat object
3567: -  col - the column index

3569:    Output Parameter:
3570: .  v - the vector

3572:    Notes:
3573:      The vector is owned by PETSc and users cannot modify it.
3574:      Users need to call MatDenseRestoreColumnVecRead() when the vector is no longer needed.
3575:      Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecWrite() for write-only access.

3577:    Level: intermediate

3579: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3580: @*/
3581: PetscErrorCode MatDenseGetColumnVecRead(Mat A,PetscInt col,Vec *v)
3582: {

3590:   if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3591:   if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N);
3592:   PetscUseMethod(A,"MatDenseGetColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));
3593:   return(0);
3594: }

3596: /*@C
3597:    MatDenseRestoreColumnVecRead - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecRead().

3599:    Collective

3601:    Input Parameters:
3602: +  mat - the Mat object
3603: .  col - the column index
3604: -  v - the Vec object

3606:    Level: intermediate

3608: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecWrite()
3609: @*/
3610: PetscErrorCode MatDenseRestoreColumnVecRead(Mat A,PetscInt col,Vec *v)
3611: {

3619:   if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3620:   if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N);
3621:   PetscUseMethod(A,"MatDenseRestoreColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));
3622:   return(0);
3623: }

3625: /*@C
3626:    MatDenseGetColumnVecWrite - Gives write-only access to a column of a dense matrix, represented as a Vec.

3628:    Collective

3630:    Input Parameters:
3631: +  mat - the Mat object
3632: -  col - the column index

3634:    Output Parameter:
3635: .  v - the vector

3637:    Notes:
3638:      The vector is owned by PETSc. Users need to call MatDenseRestoreColumnVecWrite() when the vector is no longer needed.
3639:      Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecRead() for read-only access.

3641:    Level: intermediate

3643: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3644: @*/
3645: PetscErrorCode MatDenseGetColumnVecWrite(Mat A,PetscInt col,Vec *v)
3646: {

3654:   if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3655:   if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N);
3656:   PetscUseMethod(A,"MatDenseGetColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));
3657:   return(0);
3658: }

3660: /*@C
3661:    MatDenseRestoreColumnVecWrite - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecWrite().

3663:    Collective

3665:    Input Parameters:
3666: +  mat - the Mat object
3667: .  col - the column index
3668: -  v - the Vec object

3670:    Level: intermediate

3672: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead()
3673: @*/
3674: PetscErrorCode MatDenseRestoreColumnVecWrite(Mat A,PetscInt col,Vec *v)
3675: {

3683:   if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3684:   if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N);
3685:   PetscUseMethod(A,"MatDenseRestoreColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));
3686:   return(0);
3687: }

3689: /*@C
3690:    MatDenseGetSubMatrix - Gives access to a block of columns of a dense matrix, represented as a Mat.

3692:    Collective

3694:    Input Parameters:
3695: +  mat - the Mat object
3696: .  cbegin - the first index in the block
3697: -  cend - the last index in the block

3699:    Output Parameter:
3700: .  v - the matrix

3702:    Notes:
3703:      The matrix is owned by PETSc. Users need to call MatDenseRestoreSubMatrix() when the matrix is no longer needed.

3705:    Level: intermediate

3707: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseRestoreColumnVec(), MatDenseRestoreSubMatrix()
3708: @*/
3709: PetscErrorCode MatDenseGetSubMatrix(Mat A,PetscInt cbegin,PetscInt cend,Mat *v)
3710: {

3719:   if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3720:   if (cbegin < 0 || cbegin > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid cbegin %D, should be in [0,%D)",cbegin,A->cmap->N);
3721:   if (cend < cbegin || cend > A->cmap->N) SETERRQ3(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid cend %D, should be in [%D,%D)",cend,cbegin,A->cmap->N);
3722:   PetscUseMethod(A,"MatDenseGetSubMatrix_C",(Mat,PetscInt,PetscInt,Mat*),(A,cbegin,cend,v));
3723:   return(0);
3724: }

3726: /*@C
3727:    MatDenseRestoreSubMatrix - Returns access to a block of columns of a dense matrix obtained from MatDenseGetSubMatrix().

3729:    Collective

3731:    Input Parameters:
3732: +  mat - the Mat object
3733: -  v - the Mat object

3735:    Level: intermediate

3737: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseRestoreColumnVec(), MatDenseGetSubMatrix()
3738: @*/
3739: PetscErrorCode MatDenseRestoreSubMatrix(Mat A,Mat *v)
3740: {

3747:   PetscUseMethod(A,"MatDenseRestoreSubMatrix_C",(Mat,Mat*),(A,v));
3748:   return(0);
3749: }