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;

 18:   MatDenseGetArray(A,&v);
 19:   if (!hermitian) {
 20:     for (k=0;k<n;k++) {
 21:       for (j=k;j<n;j++) {
 22:         v[j*mat->lda + k] = v[k*mat->lda + j];
 23:       }
 24:     }
 25:   } else {
 26:     for (k=0;k<n;k++) {
 27:       for (j=k;j<n;j++) {
 28:         v[j*mat->lda + k] = PetscConj(v[k*mat->lda + j]);
 29:       }
 30:     }
 31:   }
 32:   MatDenseRestoreArray(A,&v);
 33:   return 0;
 34: }

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

 41:   if (!A->rmap->n || !A->cmap->n) return 0;
 42:   PetscBLASIntCast(A->cmap->n,&n);
 43:   if (A->factortype == MAT_FACTOR_LU) {
 45:     if (!mat->fwork) {
 46:       mat->lfwork = n;
 47:       PetscMalloc1(mat->lfwork,&mat->fwork);
 48:       PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
 49:     }
 50:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
 51:     PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
 52:     PetscFPTrapPop();
 53:     PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
 54:   } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
 55:     if (A->spd) {
 56:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
 57:       PetscStackCallBLAS("LAPACKpotri",LAPACKpotri_("L",&n,mat->v,&mat->lda,&info));
 58:       PetscFPTrapPop();
 59:       MatSeqDenseSymmetrize_Private(A,PETSC_TRUE);
 60: #if defined(PETSC_USE_COMPLEX)
 61:     } else if (A->hermitian) {
 64:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
 65:       PetscStackCallBLAS("LAPACKhetri",LAPACKhetri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info));
 66:       PetscFPTrapPop();
 67:       MatSeqDenseSymmetrize_Private(A,PETSC_TRUE);
 68: #endif
 69:     } else { /* symmetric case */
 72:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
 73:       PetscStackCallBLAS("LAPACKsytri",LAPACKsytri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info));
 74:       PetscFPTrapPop();
 75:       MatSeqDenseSymmetrize_Private(A,PETSC_FALSE);
 76:     }
 78:     PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
 79:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");

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

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

 99:   if (PetscDefined(USE_DEBUG)) {
100:     for (i=0; i<N; i++) {
104:     }
105:   }
106:   if (!N) return 0;

108:   /* fix right hand side if needed */
109:   if (x && b) {
110:     Vec xt;

113:     VecDuplicate(x,&xt);
114:     VecCopy(x,xt);
115:     VecScale(xt,-1.0);
116:     MatMultAdd(A,xt,b,b);
117:     VecDestroy(&xt);
118:     VecGetArrayRead(x,&xx);
119:     VecGetArray(b,&bb);
120:     for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
121:     VecRestoreArrayRead(x,&xx);
122:     VecRestoreArray(b,&bb);
123:   }

125:   MatDenseGetArray(A,&v);
126:   for (i=0; i<N; i++) {
127:     slot = v + rows[i]*m;
128:     PetscArrayzero(slot,r);
129:   }
130:   for (i=0; i<N; i++) {
131:     slot = v + rows[i];
132:     for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
133:   }
134:   if (diag != 0.0) {
136:     for (i=0; i<N; i++) {
137:       slot  = v + (m+1)*rows[i];
138:       *slot = diag;
139:     }
140:   }
141:   MatDenseRestoreArray(A,&v);
142:   return 0;
143: }

145: PetscErrorCode MatPtAPNumeric_SeqDense_SeqDense(Mat A,Mat P,Mat C)
146: {
147:   Mat_SeqDense   *c = (Mat_SeqDense*)(C->data);

149:   if (c->ptapwork) {
150:     (*C->ops->matmultnumeric)(A,P,c->ptapwork);
151:     (*C->ops->transposematmultnumeric)(P,c->ptapwork,C);
152:   } else SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"Must call MatPtAPSymbolic_SeqDense_SeqDense() first");
153:   return 0;
154: }

156: PetscErrorCode MatPtAPSymbolic_SeqDense_SeqDense(Mat A,Mat P,PetscReal fill,Mat C)
157: {
158:   Mat_SeqDense   *c;
159:   PetscBool      cisdense;

161:   MatSetSizes(C,P->cmap->n,P->cmap->n,P->cmap->N,P->cmap->N);
162:   PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
163:   if (!cisdense) {
164:     PetscBool flg;

166:     PetscObjectTypeCompare((PetscObject)P,((PetscObject)A)->type_name,&flg);
167:     MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
168:   }
169:   MatSetUp(C);
170:   c    = (Mat_SeqDense*)C->data;
171:   MatCreate(PetscObjectComm((PetscObject)A),&c->ptapwork);
172:   MatSetSizes(c->ptapwork,A->rmap->n,P->cmap->n,A->rmap->N,P->cmap->N);
173:   MatSetType(c->ptapwork,((PetscObject)C)->type_name);
174:   MatSetUp(c->ptapwork);
175:   return 0;
176: }

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

187:   if (reuse == MAT_REUSE_MATRIX) {
188:     PetscObjectTypeCompare((PetscObject)*newmat,MATSEQDENSE,&isseqdense);
190:   }
191:   if (reuse != MAT_REUSE_MATRIX) {
192:     MatCreate(PetscObjectComm((PetscObject)A),&B);
193:     MatSetSizes(B,m,n,m,n);
194:     MatSetType(B,MATSEQDENSE);
195:     MatSeqDenseSetPreallocation(B,NULL);
196:     b    = (Mat_SeqDense*)(B->data);
197:   } else {
198:     b    = (Mat_SeqDense*)((*newmat)->data);
199:     PetscArrayzero(b->v,m*n);
200:   }
201:   MatSeqAIJGetArrayRead(A,&av);
202:   for (i=0; i<m; i++) {
203:     PetscInt j;
204:     for (j=0;j<ai[1]-ai[0];j++) {
205:       b->v[*aj*m+i] = *av;
206:       aj++;
207:       av++;
208:     }
209:     ai++;
210:   }
211:   MatSeqAIJRestoreArrayRead(A,&av);

213:   if (reuse == MAT_INPLACE_MATRIX) {
214:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
215:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
216:     MatHeaderReplace(A,&B);
217:   } else {
218:     if (B) *newmat = B;
219:     MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
220:     MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
221:   }
222:   return 0;
223: }

225: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
226: {
227:   Mat            B = NULL;
228:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
229:   PetscInt       i, j;
230:   PetscInt       *rows, *nnz;
231:   MatScalar      *aa = a->v, *vals;

233:   PetscCalloc3(A->rmap->n,&rows,A->rmap->n,&nnz,A->rmap->n,&vals);
234:   if (reuse != MAT_REUSE_MATRIX) {
235:     MatCreate(PetscObjectComm((PetscObject)A),&B);
236:     MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
237:     MatSetType(B,MATSEQAIJ);
238:     for (j=0; j<A->cmap->n; j++) {
239:       for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || (i == j && A->cmap->n == A->rmap->n)) ++nnz[i];
240:       aa += a->lda;
241:     }
242:     MatSeqAIJSetPreallocation(B,PETSC_DETERMINE,nnz);
243:   } else B = *newmat;
244:   aa = a->v;
245:   for (j=0; j<A->cmap->n; j++) {
246:     PetscInt numRows = 0;
247:     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];}
248:     MatSetValues(B,numRows,rows,1,&j,vals,INSERT_VALUES);
249:     aa  += a->lda;
250:   }
251:   PetscFree3(rows,nnz,vals);
252:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
253:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

255:   if (reuse == MAT_INPLACE_MATRIX) {
256:     MatHeaderReplace(A,&B);
257:   } else if (reuse != MAT_REUSE_MATRIX) *newmat = B;
258:   return 0;
259: }

261: PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str)
262: {
263:   Mat_SeqDense      *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data;
264:   const PetscScalar *xv;
265:   PetscScalar       *yv;
266:   PetscBLASInt      N,m,ldax = 0,lday = 0,one = 1;

268:   MatDenseGetArrayRead(X,&xv);
269:   MatDenseGetArray(Y,&yv);
270:   PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);
271:   PetscBLASIntCast(X->rmap->n,&m);
272:   PetscBLASIntCast(x->lda,&ldax);
273:   PetscBLASIntCast(y->lda,&lday);
274:   if (ldax>m || lday>m) {
275:     PetscInt j;

277:     for (j=0; j<X->cmap->n; j++) {
278:       PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&alpha,xv+j*ldax,&one,yv+j*lday,&one));
279:     }
280:   } else {
281:     PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&alpha,xv,&one,yv,&one));
282:   }
283:   MatDenseRestoreArrayRead(X,&xv);
284:   MatDenseRestoreArray(Y,&yv);
285:   PetscLogFlops(PetscMax(2.0*N-1,0));
286:   return 0;
287: }

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

293:   info->block_size        = 1.0;
294:   info->nz_allocated      = N;
295:   info->nz_used           = N;
296:   info->nz_unneeded       = 0;
297:   info->assemblies        = A->num_ass;
298:   info->mallocs           = 0;
299:   info->memory            = ((PetscObject)A)->mem;
300:   info->fill_ratio_given  = 0;
301:   info->fill_ratio_needed = 0;
302:   info->factor_mallocs    = 0;
303:   return 0;
304: }

306: PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha)
307: {
308:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
309:   PetscScalar    *v;
310:   PetscBLASInt   one = 1,j,nz,lda = 0;

312:   MatDenseGetArray(A,&v);
313:   PetscBLASIntCast(a->lda,&lda);
314:   if (lda>A->rmap->n) {
315:     PetscBLASIntCast(A->rmap->n,&nz);
316:     for (j=0; j<A->cmap->n; j++) {
317:       PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&alpha,v+j*lda,&one));
318:     }
319:   } else {
320:     PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);
321:     PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&alpha,v,&one));
322:   }
323:   PetscLogFlops(nz);
324:   MatDenseRestoreArray(A,&v);
325:   return 0;
326: }

328: static PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool  *fl)
329: {
330:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
331:   PetscInt          i,j,m = A->rmap->n,N = a->lda;
332:   const PetscScalar *v;

334:   *fl = PETSC_FALSE;
335:   if (A->rmap->n != A->cmap->n) return 0;
336:   MatDenseGetArrayRead(A,&v);
337:   for (i=0; i<m; i++) {
338:     for (j=i; j<m; j++) {
339:       if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) {
340:         goto restore;
341:       }
342:     }
343:   }
344:   *fl  = PETSC_TRUE;
345: restore:
346:   MatDenseRestoreArrayRead(A,&v);
347:   return 0;
348: }

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

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

372: PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues)
373: {
374:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
375:   PetscInt       lda = (PetscInt)mat->lda,j,m,nlda = lda;
376:   PetscBool      isdensecpu;

378:   PetscLayoutReference(A->rmap,&newi->rmap);
379:   PetscLayoutReference(A->cmap,&newi->cmap);
380:   if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { /* propagate LDA */
381:     MatDenseSetLDA(newi,lda);
382:   }
383:   PetscObjectTypeCompare((PetscObject)newi,MATSEQDENSE,&isdensecpu);
384:   if (isdensecpu) MatSeqDenseSetPreallocation(newi,NULL);
385:   if (cpvalues == MAT_COPY_VALUES) {
386:     const PetscScalar *av;
387:     PetscScalar       *v;

389:     MatDenseGetArrayRead(A,&av);
390:     MatDenseGetArrayWrite(newi,&v);
391:     MatDenseGetLDA(newi,&nlda);
392:     m    = A->rmap->n;
393:     if (lda>m || nlda>m) {
394:       for (j=0; j<A->cmap->n; j++) {
395:         PetscArraycpy(v+j*nlda,av+j*lda,m);
396:       }
397:     } else {
398:       PetscArraycpy(v,av,A->rmap->n*A->cmap->n);
399:     }
400:     MatDenseRestoreArrayWrite(newi,&v);
401:     MatDenseRestoreArrayRead(A,&av);
402:   }
403:   return 0;
404: }

406: PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat)
407: {
408:   MatCreate(PetscObjectComm((PetscObject)A),newmat);
409:   MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
410:   MatSetType(*newmat,((PetscObject)A)->type_name);
411:   MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);
412:   return 0;
413: }

415: static PetscErrorCode MatSolve_SeqDense_Internal_LU(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k, PetscBool T)
416: {
417:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
418:   PetscBLASInt    info;

420:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
421:   PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_(T ? "T" : "N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
422:   PetscFPTrapPop();
424:   PetscLogFlops(nrhs*(2.0*m*m - m));
425:   return 0;
426: }

428: static PetscErrorCode MatConjugate_SeqDense(Mat);

430: static PetscErrorCode MatSolve_SeqDense_Internal_Cholesky(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k, PetscBool T)
431: {
432:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
433:   PetscBLASInt    info;

435:   if (A->spd) {
436:     if (PetscDefined(USE_COMPLEX) && T) MatConjugate_SeqDense(A);
437:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
438:     PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info));
439:     PetscFPTrapPop();
441:     if (PetscDefined(USE_COMPLEX) && T) MatConjugate_SeqDense(A);
442: #if defined(PETSC_USE_COMPLEX)
443:   } else if (A->hermitian) {
444:     if (T) MatConjugate_SeqDense(A);
445:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
446:     PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
447:     PetscFPTrapPop();
449:     if (T) MatConjugate_SeqDense(A);
450: #endif
451:   } else { /* symmetric case */
452:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
453:     PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
454:     PetscFPTrapPop();
456:   }
457:   PetscLogFlops(nrhs*(2.0*m*m - m));
458:   return 0;
459: }

461: static PetscErrorCode MatSolve_SeqDense_Internal_QR(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k)
462: {
463:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
464:   PetscBLASInt    info;
465:   char            trans;

467:   if (PetscDefined(USE_COMPLEX)) {
468:     trans = 'C';
469:   } else {
470:     trans = 'T';
471:   }
472:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
473:   PetscStackCallBLAS("LAPACKormqr",LAPACKormqr_("L", &trans, &m,&nrhs,&mat->rank,mat->v,&mat->lda,mat->tau,x,&ldx,mat->fwork,&mat->lfwork,&info));
474:   PetscFPTrapPop();
476:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
477:   PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U", "N", "N", &mat->rank,&nrhs,mat->v,&mat->lda,x,&ldx,&info));
478:   PetscFPTrapPop();
480:   for (PetscInt j = 0; j < nrhs; j++) {
481:     for (PetscInt i = mat->rank; i < k; i++) {
482:       x[j*ldx + i] = 0.;
483:     }
484:   }
485:   PetscLogFlops(nrhs*(4.0*m*mat->rank - PetscSqr(mat->rank)));
486:   return 0;
487: }

489: static PetscErrorCode MatSolveTranspose_SeqDense_Internal_QR(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k)
490: {
491:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
492:   PetscBLASInt      info;

494:   if (A->rmap->n == A->cmap->n && mat->rank == A->rmap->n) {
495:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
496:     PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U", "T", "N", &m,&nrhs,mat->v,&mat->lda,x,&ldx,&info));
497:     PetscFPTrapPop();
499:     if (PetscDefined(USE_COMPLEX)) MatConjugate_SeqDense(A);
500:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
501:     PetscStackCallBLAS("LAPACKormqr",LAPACKormqr_("L", "N", &m,&nrhs,&mat->rank,mat->v,&mat->lda,mat->tau,x,&ldx,mat->fwork,&mat->lfwork,&info));
502:     PetscFPTrapPop();
504:     if (PetscDefined(USE_COMPLEX)) MatConjugate_SeqDense(A);
505:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"QR factored matrix cannot be used for transpose solve");
506:   PetscLogFlops(nrhs*(4.0*m*mat->rank - PetscSqr(mat->rank)));
507:   return 0;
508: }

510: static PetscErrorCode MatSolve_SeqDense_SetUp(Mat A, Vec xx, Vec yy, PetscScalar **_y, PetscBLASInt *_m, PetscBLASInt *_k)
511: {
512:   Mat_SeqDense      *mat = (Mat_SeqDense *) A->data;
513:   PetscScalar       *y;
514:   PetscBLASInt      m=0, k=0;

516:   PetscBLASIntCast(A->rmap->n,&m);
517:   PetscBLASIntCast(A->cmap->n,&k);
518:   if (k < m) {
519:     VecCopy(xx, mat->qrrhs);
520:     VecGetArray(mat->qrrhs,&y);
521:   } else {
522:     VecCopy(xx, yy);
523:     VecGetArray(yy,&y);
524:   }
525:   *_y = y;
526:   *_k = k;
527:   *_m = m;
528:   return 0;
529: }

531: static PetscErrorCode MatSolve_SeqDense_TearDown(Mat A, Vec xx, Vec yy, PetscScalar **_y, PetscBLASInt *_m, PetscBLASInt *_k)
532: {
533:   Mat_SeqDense   *mat = (Mat_SeqDense *) A->data;
534:   PetscScalar    *y = NULL;
535:   PetscBLASInt   m, k;

537:   y   = *_y;
538:   *_y = NULL;
539:   k   = *_k;
540:   m   = *_m;
541:   if (k < m) {
542:     PetscScalar *yv;
543:     VecGetArray(yy,&yv);
544:     PetscArraycpy(yv, y, k);
545:     VecRestoreArray(yy,&yv);
546:     VecRestoreArray(mat->qrrhs, &y);
547:   } else {
548:     VecRestoreArray(yy,&y);
549:   }
550:   return 0;
551: }

553: static PetscErrorCode MatSolve_SeqDense_LU(Mat A, Vec xx, Vec yy)
554: {
555:   PetscScalar    *y = NULL;
556:   PetscBLASInt   m = 0, k = 0;

558:   MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);
559:   MatSolve_SeqDense_Internal_LU(A, y, m, m, 1, k, PETSC_FALSE);
560:   MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);
561:   return 0;
562: }

564: static PetscErrorCode MatSolveTranspose_SeqDense_LU(Mat A, Vec xx, Vec yy)
565: {
566:   PetscScalar    *y = NULL;
567:   PetscBLASInt   m = 0, k = 0;

569:   MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);
570:   MatSolve_SeqDense_Internal_LU(A, y, m, m, 1, k, PETSC_TRUE);
571:   MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);
572:   return 0;
573: }

575: static PetscErrorCode MatSolve_SeqDense_Cholesky(Mat A, Vec xx, Vec yy)
576: {
577:   PetscScalar    *y = NULL;
578:   PetscBLASInt   m = 0, k = 0;

580:   MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);
581:   MatSolve_SeqDense_Internal_Cholesky(A, y, m, m, 1, k, PETSC_FALSE);
582:   MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);
583:   return 0;
584: }

586: static PetscErrorCode MatSolveTranspose_SeqDense_Cholesky(Mat A, Vec xx, Vec yy)
587: {
588:   PetscScalar    *y = NULL;
589:   PetscBLASInt   m = 0, k = 0;

591:   MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);
592:   MatSolve_SeqDense_Internal_Cholesky(A, y, m, m, 1, k, PETSC_TRUE);
593:   MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);
594:   return 0;
595: }

597: static PetscErrorCode MatSolve_SeqDense_QR(Mat A, Vec xx, Vec yy)
598: {
599:   PetscScalar    *y = NULL;
600:   PetscBLASInt   m = 0, k = 0;

602:   MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);
603:   MatSolve_SeqDense_Internal_QR(A, y, PetscMax(m,k), m, 1, k);
604:   MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);
605:   return 0;
606: }

608: static PetscErrorCode MatSolveTranspose_SeqDense_QR(Mat A, Vec xx, Vec yy)
609: {
610:   PetscScalar    *y = NULL;
611:   PetscBLASInt   m = 0, k = 0;

613:   MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);
614:   MatSolveTranspose_SeqDense_Internal_QR(A, y, PetscMax(m,k), m, 1, k);
615:   MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);
616:   return 0;
617: }

619: static PetscErrorCode MatMatSolve_SeqDense_SetUp(Mat A, Mat B, Mat X, PetscScalar **_y, PetscBLASInt *_ldy, PetscBLASInt *_m, PetscBLASInt *_nrhs, PetscBLASInt *_k)
620: {
621:   const PetscScalar *b;
622:   PetscScalar       *y;
623:   PetscInt          n, _ldb, _ldx;
624:   PetscBLASInt      nrhs=0,m=0,k=0,ldb=0,ldx=0,ldy=0;

626:   *_ldy=0; *_m=0; *_nrhs=0; *_k=0; *_y = NULL;
627:   PetscBLASIntCast(A->rmap->n,&m);
628:   PetscBLASIntCast(A->cmap->n,&k);
629:   MatGetSize(B,NULL,&n);
630:   PetscBLASIntCast(n,&nrhs);
631:   MatDenseGetLDA(B,&_ldb);
632:   PetscBLASIntCast(_ldb, &ldb);
633:   MatDenseGetLDA(X,&_ldx);
634:   PetscBLASIntCast(_ldx, &ldx);
635:   if (ldx < m) {
636:     MatDenseGetArrayRead(B,&b);
637:     PetscMalloc1(nrhs * m, &y);
638:     if (ldb == m) {
639:       PetscArraycpy(y,b,ldb*nrhs);
640:     } else {
641:       for (PetscInt j = 0; j < nrhs; j++) {
642:         PetscArraycpy(&y[j*m],&b[j*ldb],m);
643:       }
644:     }
645:     ldy = m;
646:     MatDenseRestoreArrayRead(B,&b);
647:   } else {
648:     if (ldb == ldx) {
649:       MatCopy(B, X, SAME_NONZERO_PATTERN);
650:       MatDenseGetArray(X,&y);
651:     } else {
652:       MatDenseGetArray(X,&y);
653:       MatDenseGetArrayRead(B,&b);
654:       for (PetscInt j = 0; j < nrhs; j++) {
655:         PetscArraycpy(&y[j*ldx],&b[j*ldb],m);
656:       }
657:       MatDenseRestoreArrayRead(B,&b);
658:     }
659:     ldy = ldx;
660:   }
661:   *_y    = y;
662:   *_ldy = ldy;
663:   *_k    = k;
664:   *_m    = m;
665:   *_nrhs = nrhs;
666:   return 0;
667: }

669: static PetscErrorCode MatMatSolve_SeqDense_TearDown(Mat A, Mat B, Mat X, PetscScalar **_y, PetscBLASInt *_ldy, PetscBLASInt *_m, PetscBLASInt *_nrhs, PetscBLASInt *_k)
670: {
671:   PetscScalar       *y;
672:   PetscInt          _ldx;
673:   PetscBLASInt      k,ldy,nrhs,ldx=0;

675:   y    = *_y;
676:   *_y  = NULL;
677:   k    = *_k;
678:   ldy = *_ldy;
679:   nrhs = *_nrhs;
680:   MatDenseGetLDA(X,&_ldx);
681:   PetscBLASIntCast(_ldx, &ldx);
682:   if (ldx != ldy) {
683:     PetscScalar *xv;
684:     MatDenseGetArray(X,&xv);
685:     for (PetscInt j = 0; j < nrhs; j++) {
686:       PetscArraycpy(&xv[j*ldx],&y[j*ldy],k);
687:     }
688:     MatDenseRestoreArray(X,&xv);
689:     PetscFree(y);
690:   } else {
691:     MatDenseRestoreArray(X,&y);
692:   }
693:   return 0;
694: }

696: static PetscErrorCode MatMatSolve_SeqDense_LU(Mat A, Mat B, Mat X)
697: {
698:   PetscScalar    *y;
699:   PetscBLASInt   m, k, ldy, nrhs;

701:   MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k);
702:   MatSolve_SeqDense_Internal_LU(A, y, ldy, m, nrhs, k, PETSC_FALSE);
703:   MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k);
704:   return 0;
705: }

707: static PetscErrorCode MatMatSolveTranspose_SeqDense_LU(Mat A, Mat B, Mat X)
708: {
709:   PetscScalar    *y;
710:   PetscBLASInt   m, k, ldy, nrhs;

712:   MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k);
713:   MatSolve_SeqDense_Internal_LU(A, y, ldy, m, nrhs, k, PETSC_TRUE);
714:   MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k);
715:   return 0;
716: }

718: static PetscErrorCode MatMatSolve_SeqDense_Cholesky(Mat A, Mat B, Mat X)
719: {
720:   PetscScalar    *y;
721:   PetscBLASInt   m, k, ldy, nrhs;

723:   MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k);
724:   MatSolve_SeqDense_Internal_Cholesky(A, y, ldy, m, nrhs, k, PETSC_FALSE);
725:   MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k);
726:   return 0;
727: }

729: static PetscErrorCode MatMatSolveTranspose_SeqDense_Cholesky(Mat A, Mat B, Mat X)
730: {
731:   PetscScalar    *y;
732:   PetscBLASInt   m, k, ldy, nrhs;

734:   MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k);
735:   MatSolve_SeqDense_Internal_Cholesky(A, y, ldy, m, nrhs, k, PETSC_TRUE);
736:   MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k);
737:   return 0;
738: }

740: static PetscErrorCode MatMatSolve_SeqDense_QR(Mat A, Mat B, Mat X)
741: {
742:   PetscScalar    *y;
743:   PetscBLASInt   m, k, ldy, nrhs;

745:   MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k);
746:   MatSolve_SeqDense_Internal_QR(A, y, ldy, m, nrhs, k);
747:   MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k);
748:   return 0;
749: }

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

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

762: static PetscErrorCode MatConjugate_SeqDense(Mat);

764: /* ---------------------------------------------------------------*/
765: /* COMMENT: I have chosen to hide row permutation in the pivots,
766:    rather than put it in the Mat->row slot.*/
767: PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo)
768: {
769:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
770:   PetscBLASInt   n,m,info;

772:   PetscBLASIntCast(A->cmap->n,&n);
773:   PetscBLASIntCast(A->rmap->n,&m);
774:   if (!mat->pivots) {
775:     PetscMalloc1(A->rmap->n,&mat->pivots);
776:     PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
777:   }
778:   if (!A->rmap->n || !A->cmap->n) return 0;
779:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
780:   PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info));
781:   PetscFPTrapPop();


786:   A->ops->solve             = MatSolve_SeqDense_LU;
787:   A->ops->matsolve          = MatMatSolve_SeqDense_LU;
788:   A->ops->solvetranspose    = MatSolveTranspose_SeqDense_LU;
789:   A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_LU;
790:   A->factortype             = MAT_FACTOR_LU;

792:   PetscFree(A->solvertype);
793:   PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);

795:   PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);
796:   return 0;
797: }

799: static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
800: {
801:   MatFactorInfo  info;

803:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
804:   (*fact->ops->lufactor)(fact,NULL,NULL,&info);
805:   return 0;
806: }

808: PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info)
809: {
810:   fact->preallocated           = PETSC_TRUE;
811:   fact->assembled              = PETSC_TRUE;
812:   fact->ops->lufactornumeric   = MatLUFactorNumeric_SeqDense;
813:   return 0;
814: }

816: /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */
817: PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo)
818: {
819:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
820:   PetscBLASInt   info,n;

822:   PetscBLASIntCast(A->cmap->n,&n);
823:   if (!A->rmap->n || !A->cmap->n) return 0;
824:   if (A->spd) {
825:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
826:     PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info));
827:     PetscFPTrapPop();
828: #if defined(PETSC_USE_COMPLEX)
829:   } else if (A->hermitian) {
830:     if (!mat->pivots) {
831:       PetscMalloc1(A->rmap->n,&mat->pivots);
832:       PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
833:     }
834:     if (!mat->fwork) {
835:       PetscScalar dummy;

837:       mat->lfwork = -1;
838:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
839:       PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
840:       PetscFPTrapPop();
841:       mat->lfwork = (PetscInt)PetscRealPart(dummy);
842:       PetscMalloc1(mat->lfwork,&mat->fwork);
843:       PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
844:     }
845:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
846:     PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
847:     PetscFPTrapPop();
848: #endif
849:   } else { /* symmetric case */
850:     if (!mat->pivots) {
851:       PetscMalloc1(A->rmap->n,&mat->pivots);
852:       PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
853:     }
854:     if (!mat->fwork) {
855:       PetscScalar dummy;

857:       mat->lfwork = -1;
858:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
859:       PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
860:       PetscFPTrapPop();
861:       mat->lfwork = (PetscInt)PetscRealPart(dummy);
862:       PetscMalloc1(mat->lfwork,&mat->fwork);
863:       PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
864:     }
865:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
866:     PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
867:     PetscFPTrapPop();
868:   }

871:   A->ops->solve             = MatSolve_SeqDense_Cholesky;
872:   A->ops->matsolve          = MatMatSolve_SeqDense_Cholesky;
873:   A->ops->solvetranspose    = MatSolveTranspose_SeqDense_Cholesky;
874:   A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_Cholesky;
875:   A->factortype             = MAT_FACTOR_CHOLESKY;

877:   PetscFree(A->solvertype);
878:   PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);

880:   PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
881:   return 0;
882: }

884: static PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
885: {
886:   MatFactorInfo  info;

888:   info.fill = 1.0;

890:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
891:   (*fact->ops->choleskyfactor)(fact,NULL,&info);
892:   return 0;
893: }

895: PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
896: {
897:   fact->assembled                  = PETSC_TRUE;
898:   fact->preallocated               = PETSC_TRUE;
899:   fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
900:   return 0;
901: }

903: PetscErrorCode MatQRFactor_SeqDense(Mat A,IS col,const MatFactorInfo *minfo)
904: {
905:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
906:   PetscBLASInt   n,m,info, min, max;

908:   PetscBLASIntCast(A->cmap->n,&n);
909:   PetscBLASIntCast(A->rmap->n,&m);
910:   max = PetscMax(m, n);
911:   min = PetscMin(m, n);
912:   if (!mat->tau) {
913:     PetscMalloc1(min,&mat->tau);
914:     PetscLogObjectMemory((PetscObject)A,min*sizeof(PetscScalar));
915:   }
916:   if (!mat->pivots) {
917:     PetscMalloc1(n,&mat->pivots);
918:     PetscLogObjectMemory((PetscObject)A,n*sizeof(PetscScalar));
919:   }
920:   if (!mat->qrrhs) {
921:     MatCreateVecs(A, NULL, &(mat->qrrhs));
922:   }
923:   if (!A->rmap->n || !A->cmap->n) return 0;
924:   if (!mat->fwork) {
925:     PetscScalar dummy;

927:     mat->lfwork = -1;
928:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
929:     PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&m,&n,mat->v,&mat->lda,mat->tau,&dummy,&mat->lfwork,&info));
930:     PetscFPTrapPop();
931:     mat->lfwork = (PetscInt)PetscRealPart(dummy);
932:     PetscMalloc1(mat->lfwork,&mat->fwork);
933:     PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
934:   }
935:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
936:   PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&m,&n,mat->v,&mat->lda,mat->tau,mat->fwork,&mat->lfwork,&info));
937:   PetscFPTrapPop();
939:   // 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
940:   mat->rank = min;

942:   A->ops->solve             = MatSolve_SeqDense_QR;
943:   A->ops->matsolve          = MatMatSolve_SeqDense_QR;
944:   A->factortype             = MAT_FACTOR_QR;
945:   if (m == n) {
946:     A->ops->solvetranspose    = MatSolveTranspose_SeqDense_QR;
947:     A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_QR;
948:   }

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

953:   PetscLogFlops(2.0*min*min*(max-min/3.0));
954:   return 0;
955: }

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

961:   info.fill = 1.0;

963:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
964:   PetscUseMethod(fact,"MatQRFactor_C",(Mat,IS,const MatFactorInfo *),(fact,NULL,&info));
965:   return 0;
966: }

968: PetscErrorCode MatQRFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
969: {
970:   fact->assembled                  = PETSC_TRUE;
971:   fact->preallocated               = PETSC_TRUE;
972:   PetscObjectComposeFunction((PetscObject)fact,"MatQRFactorNumeric_C",MatQRFactorNumeric_SeqDense);
973:   return 0;
974: }

976: /* uses LAPACK */
977: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact)
978: {
979:   MatCreate(PetscObjectComm((PetscObject)A),fact);
980:   MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
981:   MatSetType(*fact,MATDENSE);
982:   (*fact)->trivialsymbolic = PETSC_TRUE;
983:   if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
984:     (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
985:     (*fact)->ops->ilufactorsymbolic = MatLUFactorSymbolic_SeqDense;
986:   } else if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) {
987:     (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
988:   } else if (ftype == MAT_FACTOR_QR) {
989:     PetscObjectComposeFunction((PetscObject)(*fact),"MatQRFactorSymbolic_C",MatQRFactorSymbolic_SeqDense);
990:   }
991:   (*fact)->factortype = ftype;

993:   PetscFree((*fact)->solvertype);
994:   PetscStrallocpy(MATSOLVERPETSC,&(*fact)->solvertype);
995:   PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&(*fact)->preferredordering[MAT_FACTOR_LU]);
996:   PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&(*fact)->preferredordering[MAT_FACTOR_ILU]);
997:   PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&(*fact)->preferredordering[MAT_FACTOR_CHOLESKY]);
998:   PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&(*fact)->preferredordering[MAT_FACTOR_ICC]);
999:   return 0;
1000: }

1002: /* ------------------------------------------------------------------*/
1003: static PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx)
1004: {
1005:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1006:   PetscScalar       *x,*v = mat->v,zero = 0.0,xt;
1007:   const PetscScalar *b;
1008:   PetscInt          m = A->rmap->n,i;
1009:   PetscBLASInt      o = 1,bm = 0;

1011: #if defined(PETSC_HAVE_CUDA)
1013: #endif
1014:   if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */
1015:   PetscBLASIntCast(m,&bm);
1016:   if (flag & SOR_ZERO_INITIAL_GUESS) {
1017:     /* this is a hack fix, should have another version without the second BLASdotu */
1018:     VecSet(xx,zero);
1019:   }
1020:   VecGetArray(xx,&x);
1021:   VecGetArrayRead(bb,&b);
1022:   its  = its*lits;
1024:   while (its--) {
1025:     if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
1026:       for (i=0; i<m; i++) {
1027:         PetscStackCallBLAS("BLASdotu",xt   = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o));
1028:         x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
1029:       }
1030:     }
1031:     if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
1032:       for (i=m-1; i>=0; i--) {
1033:         PetscStackCallBLAS("BLASdotu",xt   = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o));
1034:         x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
1035:       }
1036:     }
1037:   }
1038:   VecRestoreArrayRead(bb,&b);
1039:   VecRestoreArray(xx,&x);
1040:   return 0;
1041: }

1043: /* -----------------------------------------------------------------*/
1044: PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)
1045: {
1046:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1047:   const PetscScalar *v   = mat->v,*x;
1048:   PetscScalar       *y;
1049:   PetscBLASInt      m, n,_One=1;
1050:   PetscScalar       _DOne=1.0,_DZero=0.0;

1052:   PetscBLASIntCast(A->rmap->n,&m);
1053:   PetscBLASIntCast(A->cmap->n,&n);
1054:   VecGetArrayRead(xx,&x);
1055:   VecGetArrayWrite(yy,&y);
1056:   if (!A->rmap->n || !A->cmap->n) {
1057:     PetscBLASInt i;
1058:     for (i=0; i<n; i++) y[i] = 0.0;
1059:   } else {
1060:     PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One));
1061:     PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);
1062:   }
1063:   VecRestoreArrayRead(xx,&x);
1064:   VecRestoreArrayWrite(yy,&y);
1065:   return 0;
1066: }

1068: PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy)
1069: {
1070:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1071:   PetscScalar       *y,_DOne=1.0,_DZero=0.0;
1072:   PetscBLASInt      m, n, _One=1;
1073:   const PetscScalar *v = mat->v,*x;

1075:   PetscBLASIntCast(A->rmap->n,&m);
1076:   PetscBLASIntCast(A->cmap->n,&n);
1077:   VecGetArrayRead(xx,&x);
1078:   VecGetArrayWrite(yy,&y);
1079:   if (!A->rmap->n || !A->cmap->n) {
1080:     PetscBLASInt i;
1081:     for (i=0; i<m; i++) y[i] = 0.0;
1082:   } else {
1083:     PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One));
1084:     PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);
1085:   }
1086:   VecRestoreArrayRead(xx,&x);
1087:   VecRestoreArrayWrite(yy,&y);
1088:   return 0;
1089: }

1091: PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
1092: {
1093:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1094:   const PetscScalar *v = mat->v,*x;
1095:   PetscScalar       *y,_DOne=1.0;
1096:   PetscBLASInt      m, n, _One=1;

1098:   PetscBLASIntCast(A->rmap->n,&m);
1099:   PetscBLASIntCast(A->cmap->n,&n);
1100:   VecCopy(zz,yy);
1101:   if (!A->rmap->n || !A->cmap->n) return 0;
1102:   VecGetArrayRead(xx,&x);
1103:   VecGetArray(yy,&y);
1104:   PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
1105:   VecRestoreArrayRead(xx,&x);
1106:   VecRestoreArray(yy,&y);
1107:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
1108:   return 0;
1109: }

1111: PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
1112: {
1113:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1114:   const PetscScalar *v = mat->v,*x;
1115:   PetscScalar       *y;
1116:   PetscBLASInt      m, n, _One=1;
1117:   PetscScalar       _DOne=1.0;

1119:   PetscBLASIntCast(A->rmap->n,&m);
1120:   PetscBLASIntCast(A->cmap->n,&n);
1121:   VecCopy(zz,yy);
1122:   if (!A->rmap->n || !A->cmap->n) return 0;
1123:   VecGetArrayRead(xx,&x);
1124:   VecGetArray(yy,&y);
1125:   PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
1126:   VecRestoreArrayRead(xx,&x);
1127:   VecRestoreArray(yy,&y);
1128:   PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
1129:   return 0;
1130: }

1132: /* -----------------------------------------------------------------*/
1133: static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
1134: {
1135:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1136:   PetscInt       i;

1138:   *ncols = A->cmap->n;
1139:   if (cols) {
1140:     PetscMalloc1(A->cmap->n,cols);
1141:     for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
1142:   }
1143:   if (vals) {
1144:     const PetscScalar *v;

1146:     MatDenseGetArrayRead(A,&v);
1147:     PetscMalloc1(A->cmap->n,vals);
1148:     v   += row;
1149:     for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
1150:     MatDenseRestoreArrayRead(A,&v);
1151:   }
1152:   return 0;
1153: }

1155: static PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
1156: {
1157:   if (ncols) *ncols = 0;
1158:   if (cols) PetscFree(*cols);
1159:   if (vals) PetscFree(*vals);
1160:   return 0;
1161: }
1162: /* ----------------------------------------------------------------*/
1163: static PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
1164: {
1165:   Mat_SeqDense     *mat = (Mat_SeqDense*)A->data;
1166:   PetscScalar      *av;
1167:   PetscInt         i,j,idx=0;
1168: #if defined(PETSC_HAVE_CUDA)
1169:   PetscOffloadMask oldf;
1170: #endif

1172:   MatDenseGetArray(A,&av);
1173:   if (!mat->roworiented) {
1174:     if (addv == INSERT_VALUES) {
1175:       for (j=0; j<n; j++) {
1176:         if (indexn[j] < 0) {idx += m; continue;}
1178:         for (i=0; i<m; i++) {
1179:           if (indexm[i] < 0) {idx++; continue;}
1181:           av[indexn[j]*mat->lda + indexm[i]] = v[idx++];
1182:         }
1183:       }
1184:     } else {
1185:       for (j=0; j<n; j++) {
1186:         if (indexn[j] < 0) {idx += m; continue;}
1188:         for (i=0; i<m; i++) {
1189:           if (indexm[i] < 0) {idx++; continue;}
1191:           av[indexn[j]*mat->lda + indexm[i]] += v[idx++];
1192:         }
1193:       }
1194:     }
1195:   } else {
1196:     if (addv == INSERT_VALUES) {
1197:       for (i=0; i<m; i++) {
1198:         if (indexm[i] < 0) { idx += n; continue;}
1200:         for (j=0; j<n; j++) {
1201:           if (indexn[j] < 0) { idx++; continue;}
1203:           av[indexn[j]*mat->lda + indexm[i]] = v[idx++];
1204:         }
1205:       }
1206:     } else {
1207:       for (i=0; i<m; i++) {
1208:         if (indexm[i] < 0) { idx += n; continue;}
1210:         for (j=0; j<n; j++) {
1211:           if (indexn[j] < 0) { idx++; continue;}
1213:           av[indexn[j]*mat->lda + indexm[i]] += v[idx++];
1214:         }
1215:       }
1216:     }
1217:   }
1218:   /* hack to prevent unneeded copy to the GPU while returning the array */
1219: #if defined(PETSC_HAVE_CUDA)
1220:   oldf = A->offloadmask;
1221:   A->offloadmask = PETSC_OFFLOAD_GPU;
1222: #endif
1223:   MatDenseRestoreArray(A,&av);
1224: #if defined(PETSC_HAVE_CUDA)
1225:   A->offloadmask = (oldf == PETSC_OFFLOAD_UNALLOCATED ? PETSC_OFFLOAD_UNALLOCATED : PETSC_OFFLOAD_CPU);
1226: #endif
1227:   return 0;
1228: }

1230: static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
1231: {
1232:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1233:   const PetscScalar *vv;
1234:   PetscInt          i,j;

1236:   MatDenseGetArrayRead(A,&vv);
1237:   /* row-oriented output */
1238:   for (i=0; i<m; i++) {
1239:     if (indexm[i] < 0) {v += n;continue;}
1241:     for (j=0; j<n; j++) {
1242:       if (indexn[j] < 0) {v++; continue;}
1244:       *v++ = vv[indexn[j]*mat->lda + indexm[i]];
1245:     }
1246:   }
1247:   MatDenseRestoreArrayRead(A,&vv);
1248:   return 0;
1249: }

1251: /* -----------------------------------------------------------------*/

1253: PetscErrorCode MatView_Dense_Binary(Mat mat,PetscViewer viewer)
1254: {
1255:   PetscBool         skipHeader;
1256:   PetscViewerFormat format;
1257:   PetscInt          header[4],M,N,m,lda,i,j,k;
1258:   const PetscScalar *v;
1259:   PetscScalar       *vwork;

1261:   PetscViewerSetUp(viewer);
1262:   PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
1263:   PetscViewerGetFormat(viewer,&format);
1264:   if (skipHeader) format = PETSC_VIEWER_NATIVE;

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

1268:   /* write matrix header */
1269:   header[0] = MAT_FILE_CLASSID; header[1] = M; header[2] = N;
1270:   header[3] = (format == PETSC_VIEWER_NATIVE) ? MATRIX_BINARY_FORMAT_DENSE : M*N;
1271:   if (!skipHeader) PetscViewerBinaryWrite(viewer,header,4,PETSC_INT);

1273:   MatGetLocalSize(mat,&m,NULL);
1274:   if (format != PETSC_VIEWER_NATIVE) {
1275:     PetscInt nnz = m*N, *iwork;
1276:     /* store row lengths for each row */
1277:     PetscMalloc1(nnz,&iwork);
1278:     for (i=0; i<m; i++) iwork[i] = N;
1279:     PetscViewerBinaryWriteAll(viewer,iwork,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1280:     /* store column indices (zero start index) */
1281:     for (k=0, i=0; i<m; i++)
1282:       for (j=0; j<N; j++, k++)
1283:         iwork[k] = j;
1284:     PetscViewerBinaryWriteAll(viewer,iwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1285:     PetscFree(iwork);
1286:   }
1287:   /* store matrix values as a dense matrix in row major order */
1288:   PetscMalloc1(m*N,&vwork);
1289:   MatDenseGetArrayRead(mat,&v);
1290:   MatDenseGetLDA(mat,&lda);
1291:   for (k=0, i=0; i<m; i++)
1292:     for (j=0; j<N; j++, k++)
1293:       vwork[k] = v[i+lda*j];
1294:   MatDenseRestoreArrayRead(mat,&v);
1295:   PetscViewerBinaryWriteAll(viewer,vwork,m*N,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);
1296:   PetscFree(vwork);
1297:   return 0;
1298: }

1300: PetscErrorCode MatLoad_Dense_Binary(Mat mat,PetscViewer viewer)
1301: {
1302:   PetscBool      skipHeader;
1303:   PetscInt       header[4],M,N,m,nz,lda,i,j,k;
1304:   PetscInt       rows,cols;
1305:   PetscScalar    *v,*vwork;

1307:   PetscViewerSetUp(viewer);
1308:   PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);

1310:   if (!skipHeader) {
1311:     PetscViewerBinaryRead(viewer,header,4,NULL,PETSC_INT);
1313:     M = header[1]; N = header[2];
1316:     nz = header[3];
1318:   } else {
1319:     MatGetSize(mat,&M,&N);
1321:     nz = MATRIX_BINARY_FORMAT_DENSE;
1322:   }

1324:   /* setup global sizes if not set */
1325:   if (mat->rmap->N < 0) mat->rmap->N = M;
1326:   if (mat->cmap->N < 0) mat->cmap->N = N;
1327:   MatSetUp(mat);
1328:   /* check if global sizes are correct */
1329:   MatGetSize(mat,&rows,&cols);

1332:   MatGetSize(mat,NULL,&N);
1333:   MatGetLocalSize(mat,&m,NULL);
1334:   MatDenseGetArray(mat,&v);
1335:   MatDenseGetLDA(mat,&lda);
1336:   if (nz == MATRIX_BINARY_FORMAT_DENSE) {  /* matrix in file is dense format */
1337:     PetscInt nnz = m*N;
1338:     /* read in matrix values */
1339:     PetscMalloc1(nnz,&vwork);
1340:     PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);
1341:     /* store values in column major order */
1342:     for (j=0; j<N; j++)
1343:       for (i=0; i<m; i++)
1344:         v[i+lda*j] = vwork[i*N+j];
1345:     PetscFree(vwork);
1346:   } else { /* matrix in file is sparse format */
1347:     PetscInt nnz = 0, *rlens, *icols;
1348:     /* read in row lengths */
1349:     PetscMalloc1(m,&rlens);
1350:     PetscViewerBinaryReadAll(viewer,rlens,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1351:     for (i=0; i<m; i++) nnz += rlens[i];
1352:     /* read in column indices and values */
1353:     PetscMalloc2(nnz,&icols,nnz,&vwork);
1354:     PetscViewerBinaryReadAll(viewer,icols,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1355:     PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);
1356:     /* store values in column major order */
1357:     for (k=0, i=0; i<m; i++)
1358:       for (j=0; j<rlens[i]; j++, k++)
1359:         v[i+lda*icols[k]] = vwork[k];
1360:     PetscFree(rlens);
1361:     PetscFree2(icols,vwork);
1362:   }
1363:   MatDenseRestoreArray(mat,&v);
1364:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
1365:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
1366:   return 0;
1367: }

1369: PetscErrorCode MatLoad_SeqDense(Mat newMat, PetscViewer viewer)
1370: {
1371:   PetscBool      isbinary, ishdf5;

1375:   /* force binary viewer to load .info file if it has not yet done so */
1376:   PetscViewerSetUp(viewer);
1377:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1378:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,  &ishdf5);
1379:   if (isbinary) {
1380:     MatLoad_Dense_Binary(newMat,viewer);
1381:   } else if (ishdf5) {
1382: #if defined(PETSC_HAVE_HDF5)
1383:     MatLoad_Dense_HDF5(newMat,viewer);
1384: #else
1385:     SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1386: #endif
1387:   } else {
1388:     SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"Viewer type %s not yet supported for reading %s matrices",((PetscObject)viewer)->type_name,((PetscObject)newMat)->type_name);
1389:   }
1390:   return 0;
1391: }

1393: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
1394: {
1395:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
1396:   PetscInt          i,j;
1397:   const char        *name;
1398:   PetscScalar       *v,*av;
1399:   PetscViewerFormat format;
1400: #if defined(PETSC_USE_COMPLEX)
1401:   PetscBool         allreal = PETSC_TRUE;
1402: #endif

1404:   MatDenseGetArrayRead(A,(const PetscScalar**)&av);
1405:   PetscViewerGetFormat(viewer,&format);
1406:   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1407:     return 0;  /* do nothing for now */
1408:   } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1409:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1410:     for (i=0; i<A->rmap->n; i++) {
1411:       v    = av + i;
1412:       PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i);
1413:       for (j=0; j<A->cmap->n; j++) {
1414: #if defined(PETSC_USE_COMPLEX)
1415:         if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
1416:           PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1417:         } else if (PetscRealPart(*v)) {
1418:           PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",j,(double)PetscRealPart(*v));
1419:         }
1420: #else
1421:         if (*v) {
1422:           PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",j,(double)*v);
1423:         }
1424: #endif
1425:         v += a->lda;
1426:       }
1427:       PetscViewerASCIIPrintf(viewer,"\n");
1428:     }
1429:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1430:   } else {
1431:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1432: #if defined(PETSC_USE_COMPLEX)
1433:     /* determine if matrix has all real values */
1434:     v = av;
1435:     for (i=0; i<A->rmap->n*A->cmap->n; i++) {
1436:       if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;}
1437:     }
1438: #endif
1439:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
1440:       PetscObjectGetName((PetscObject)A,&name);
1441:       PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT " \n",A->rmap->n,A->cmap->n);
1442:       PetscViewerASCIIPrintf(viewer,"%s = zeros(%" PetscInt_FMT ",%" PetscInt_FMT ");\n",name,A->rmap->n,A->cmap->n);
1443:       PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
1444:     }

1446:     for (i=0; i<A->rmap->n; i++) {
1447:       v = av + i;
1448:       for (j=0; j<A->cmap->n; j++) {
1449: #if defined(PETSC_USE_COMPLEX)
1450:         if (allreal) {
1451:           PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));
1452:         } else {
1453:           PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1454:         }
1455: #else
1456:         PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);
1457: #endif
1458:         v += a->lda;
1459:       }
1460:       PetscViewerASCIIPrintf(viewer,"\n");
1461:     }
1462:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
1463:       PetscViewerASCIIPrintf(viewer,"];\n");
1464:     }
1465:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1466:   }
1467:   MatDenseRestoreArrayRead(A,(const PetscScalar**)&av);
1468:   PetscViewerFlush(viewer);
1469:   return 0;
1470: }

1472: #include <petscdraw.h>
1473: static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa)
1474: {
1475:   Mat               A  = (Mat) Aa;
1476:   PetscInt          m  = A->rmap->n,n = A->cmap->n,i,j;
1477:   int               color = PETSC_DRAW_WHITE;
1478:   const PetscScalar *v;
1479:   PetscViewer       viewer;
1480:   PetscReal         xl,yl,xr,yr,x_l,x_r,y_l,y_r;
1481:   PetscViewerFormat format;
1482:   PetscErrorCode    ierr;

1484:   PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
1485:   PetscViewerGetFormat(viewer,&format);
1486:   PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);

1488:   /* Loop over matrix elements drawing boxes */
1489:   MatDenseGetArrayRead(A,&v);
1490:   if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1491:     PetscDrawCollectiveBegin(draw);
1492:     /* Blue for negative and Red for positive */
1493:     for (j = 0; j < n; j++) {
1494:       x_l = j; x_r = x_l + 1.0;
1495:       for (i = 0; i < m; i++) {
1496:         y_l = m - i - 1.0;
1497:         y_r = y_l + 1.0;
1498:         if (PetscRealPart(v[j*m+i]) >  0.) color = PETSC_DRAW_RED;
1499:         else if (PetscRealPart(v[j*m+i]) <  0.) color = PETSC_DRAW_BLUE;
1500:         else continue;
1501:         PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1502:       }
1503:     }
1504:     PetscDrawCollectiveEnd(draw);
1505:   } else {
1506:     /* use contour shading to indicate magnitude of values */
1507:     /* first determine max of all nonzero values */
1508:     PetscReal minv = 0.0, maxv = 0.0;
1509:     PetscDraw popup;

1511:     for (i=0; i < m*n; i++) {
1512:       if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1513:     }
1514:     if (minv >= maxv) maxv = minv + PETSC_SMALL;
1515:     PetscDrawGetPopup(draw,&popup);
1516:     PetscDrawScalePopup(popup,minv,maxv);

1518:     PetscDrawCollectiveBegin(draw);
1519:     for (j=0; j<n; j++) {
1520:       x_l = j;
1521:       x_r = x_l + 1.0;
1522:       for (i=0; i<m; i++) {
1523:         y_l   = m - i - 1.0;
1524:         y_r   = y_l + 1.0;
1525:         color = PetscDrawRealToColor(PetscAbsScalar(v[j*m+i]),minv,maxv);
1526:         PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1527:       }
1528:     }
1529:     PetscDrawCollectiveEnd(draw);
1530:   }
1531:   MatDenseRestoreArrayRead(A,&v);
1532:   return 0;
1533: }

1535: static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
1536: {
1537:   PetscDraw      draw;
1538:   PetscBool      isnull;
1539:   PetscReal      xr,yr,xl,yl,h,w;

1541:   PetscViewerDrawGetDraw(viewer,0,&draw);
1542:   PetscDrawIsNull(draw,&isnull);
1543:   if (isnull) return 0;

1545:   xr   = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
1546:   xr  += w;          yr += h;        xl = -w;     yl = -h;
1547:   PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
1548:   PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
1549:   PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
1550:   PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
1551:   PetscDrawSave(draw);
1552:   return 0;
1553: }

1555: PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
1556: {
1557:   PetscBool      iascii,isbinary,isdraw;

1559:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1560:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1561:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1562:   if (iascii) {
1563:     MatView_SeqDense_ASCII(A,viewer);
1564:   } else if (isbinary) {
1565:     MatView_Dense_Binary(A,viewer);
1566:   } else if (isdraw) {
1567:     MatView_SeqDense_Draw(A,viewer);
1568:   }
1569:   return 0;
1570: }

1572: static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A,const PetscScalar *array)
1573: {
1574:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;

1579:   a->unplacedarray       = a->v;
1580:   a->unplaced_user_alloc = a->user_alloc;
1581:   a->v                   = (PetscScalar*) array;
1582:   a->user_alloc          = PETSC_TRUE;
1583: #if defined(PETSC_HAVE_CUDA)
1584:   A->offloadmask = PETSC_OFFLOAD_CPU;
1585: #endif
1586:   return 0;
1587: }

1589: static PetscErrorCode MatDenseResetArray_SeqDense(Mat A)
1590: {
1591:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;

1595:   a->v             = a->unplacedarray;
1596:   a->user_alloc    = a->unplaced_user_alloc;
1597:   a->unplacedarray = NULL;
1598: #if defined(PETSC_HAVE_CUDA)
1599:   A->offloadmask = PETSC_OFFLOAD_CPU;
1600: #endif
1601:   return 0;
1602: }

1604: static PetscErrorCode MatDenseReplaceArray_SeqDense(Mat A,const PetscScalar *array)
1605: {
1606:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

1610:   if (!a->user_alloc) PetscFree(a->v);
1611:   a->v           = (PetscScalar*) array;
1612:   a->user_alloc  = PETSC_FALSE;
1613: #if defined(PETSC_HAVE_CUDA)
1614:   A->offloadmask = PETSC_OFFLOAD_CPU;
1615: #endif
1616:   return 0;
1617: }

1619: PetscErrorCode MatDestroy_SeqDense(Mat mat)
1620: {
1621:   Mat_SeqDense   *l = (Mat_SeqDense*)mat->data;

1623: #if defined(PETSC_USE_LOG)
1624:   PetscLogObjectState((PetscObject)mat,"Rows %" PetscInt_FMT " Cols %" PetscInt_FMT,mat->rmap->n,mat->cmap->n);
1625: #endif
1626:   VecDestroy(&(l->qrrhs));
1627:   PetscFree(l->tau);
1628:   PetscFree(l->pivots);
1629:   PetscFree(l->fwork);
1630:   MatDestroy(&l->ptapwork);
1631:   if (!l->user_alloc) PetscFree(l->v);
1632:   if (!l->unplaced_user_alloc) PetscFree(l->unplacedarray);
1635:   VecDestroy(&l->cvec);
1636:   MatDestroy(&l->cmat);
1637:   PetscFree(mat->data);

1639:   PetscObjectChangeTypeName((PetscObject)mat,NULL);
1640:   PetscObjectComposeFunction((PetscObject)mat,"MatQRFactor_C",NULL);
1641:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetLDA_C",NULL);
1642:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseSetLDA_C",NULL);
1643:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);
1644:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);
1645:   PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL);
1646:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL);
1647:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseReplaceArray_C",NULL);
1648:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayRead_C",NULL);
1649:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayRead_C",NULL);
1650:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayWrite_C",NULL);
1651:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayWrite_C",NULL);
1652:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);
1653: #if defined(PETSC_HAVE_ELEMENTAL)
1654:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);
1655: #endif
1656: #if defined(PETSC_HAVE_SCALAPACK)
1657:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_scalapack_C",NULL);
1658: #endif
1659: #if defined(PETSC_HAVE_CUDA)
1660:   PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqdensecuda_C",NULL);
1661:   PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",NULL);
1662:   PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdense_C",NULL);
1663: #endif
1664:   PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);
1665:   PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqaij_seqdense_C",NULL);
1666:   PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdense_seqdense_C",NULL);
1667:   PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqbaij_seqdense_C",NULL);
1668:   PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqsbaij_seqdense_C",NULL);

1670:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumn_C",NULL);
1671:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumn_C",NULL);
1672:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVec_C",NULL);
1673:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVec_C",NULL);
1674:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVecRead_C",NULL);
1675:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVecRead_C",NULL);
1676:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVecWrite_C",NULL);
1677:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVecWrite_C",NULL);
1678:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetSubMatrix_C",NULL);
1679:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreSubMatrix_C",NULL);
1680:   return 0;
1681: }

1683: static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
1684: {
1685:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
1686:   PetscInt       k,j,m = A->rmap->n, M = mat->lda, n = A->cmap->n;
1687:   PetscScalar    *v,tmp;

1689:   if (reuse == MAT_INPLACE_MATRIX) {
1690:     if (m == n) { /* in place transpose */
1691:       MatDenseGetArray(A,&v);
1692:       for (j=0; j<m; j++) {
1693:         for (k=0; k<j; k++) {
1694:           tmp        = v[j + k*M];
1695:           v[j + k*M] = v[k + j*M];
1696:           v[k + j*M] = tmp;
1697:         }
1698:       }
1699:       MatDenseRestoreArray(A,&v);
1700:     } else { /* reuse memory, temporary allocates new memory */
1701:       PetscScalar *v2;
1702:       PetscLayout tmplayout;

1704:       PetscMalloc1((size_t)m*n,&v2);
1705:       MatDenseGetArray(A,&v);
1706:       for (j=0; j<n; j++) {
1707:         for (k=0; k<m; k++) v2[j + (size_t)k*n] = v[k + (size_t)j*M];
1708:       }
1709:       PetscArraycpy(v,v2,(size_t)m*n);
1710:       PetscFree(v2);
1711:       MatDenseRestoreArray(A,&v);
1712:       /* cleanup size dependent quantities */
1713:       VecDestroy(&mat->cvec);
1714:       MatDestroy(&mat->cmat);
1715:       PetscFree(mat->pivots);
1716:       PetscFree(mat->fwork);
1717:       MatDestroy(&mat->ptapwork);
1718:       /* swap row/col layouts */
1719:       mat->lda  = n;
1720:       tmplayout = A->rmap;
1721:       A->rmap   = A->cmap;
1722:       A->cmap   = tmplayout;
1723:     }
1724:   } else { /* out-of-place transpose */
1725:     Mat          tmat;
1726:     Mat_SeqDense *tmatd;
1727:     PetscScalar  *v2;
1728:     PetscInt     M2;

1730:     if (reuse == MAT_INITIAL_MATRIX) {
1731:       MatCreate(PetscObjectComm((PetscObject)A),&tmat);
1732:       MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);
1733:       MatSetType(tmat,((PetscObject)A)->type_name);
1734:       MatSeqDenseSetPreallocation(tmat,NULL);
1735:     } else tmat = *matout;

1737:     MatDenseGetArrayRead(A,(const PetscScalar**)&v);
1738:     MatDenseGetArray(tmat,&v2);
1739:     tmatd = (Mat_SeqDense*)tmat->data;
1740:     M2    = tmatd->lda;
1741:     for (j=0; j<n; j++) {
1742:       for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M];
1743:     }
1744:     MatDenseRestoreArray(tmat,&v2);
1745:     MatDenseRestoreArrayRead(A,(const PetscScalar**)&v);
1746:     MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
1747:     MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
1748:     *matout = tmat;
1749:   }
1750:   return 0;
1751: }

1753: static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool  *flg)
1754: {
1755:   Mat_SeqDense      *mat1 = (Mat_SeqDense*)A1->data;
1756:   Mat_SeqDense      *mat2 = (Mat_SeqDense*)A2->data;
1757:   PetscInt          i;
1758:   const PetscScalar *v1,*v2;

1760:   if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; return 0;}
1761:   if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; return 0;}
1762:   MatDenseGetArrayRead(A1,&v1);
1763:   MatDenseGetArrayRead(A2,&v2);
1764:   for (i=0; i<A1->cmap->n; i++) {
1765:     PetscArraycmp(v1,v2,A1->rmap->n,flg);
1766:     if (*flg == PETSC_FALSE) return 0;
1767:     v1 += mat1->lda;
1768:     v2 += mat2->lda;
1769:   }
1770:   MatDenseRestoreArrayRead(A1,&v1);
1771:   MatDenseRestoreArrayRead(A2,&v2);
1772:   *flg = PETSC_TRUE;
1773:   return 0;
1774: }

1776: static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
1777: {
1778:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1779:   PetscInt          i,n,len;
1780:   PetscScalar       *x;
1781:   const PetscScalar *vv;

1783:   VecGetSize(v,&n);
1784:   VecGetArray(v,&x);
1785:   len  = PetscMin(A->rmap->n,A->cmap->n);
1786:   MatDenseGetArrayRead(A,&vv);
1788:   for (i=0; i<len; i++) {
1789:     x[i] = vv[i*mat->lda + i];
1790:   }
1791:   MatDenseRestoreArrayRead(A,&vv);
1792:   VecRestoreArray(v,&x);
1793:   return 0;
1794: }

1796: static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
1797: {
1798:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1799:   const PetscScalar *l,*r;
1800:   PetscScalar       x,*v,*vv;
1801:   PetscInt          i,j,m = A->rmap->n,n = A->cmap->n;

1803:   MatDenseGetArray(A,&vv);
1804:   if (ll) {
1805:     VecGetSize(ll,&m);
1806:     VecGetArrayRead(ll,&l);
1808:     for (i=0; i<m; i++) {
1809:       x = l[i];
1810:       v = vv + i;
1811:       for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;}
1812:     }
1813:     VecRestoreArrayRead(ll,&l);
1814:     PetscLogFlops(1.0*n*m);
1815:   }
1816:   if (rr) {
1817:     VecGetSize(rr,&n);
1818:     VecGetArrayRead(rr,&r);
1820:     for (i=0; i<n; i++) {
1821:       x = r[i];
1822:       v = vv + i*mat->lda;
1823:       for (j=0; j<m; j++) (*v++) *= x;
1824:     }
1825:     VecRestoreArrayRead(rr,&r);
1826:     PetscLogFlops(1.0*n*m);
1827:   }
1828:   MatDenseRestoreArray(A,&vv);
1829:   return 0;
1830: }

1832: PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1833: {
1834:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
1835:   PetscScalar       *v,*vv;
1836:   PetscReal         sum = 0.0;
1837:   PetscInt          lda, m=A->rmap->n,i,j;

1839:   MatDenseGetArrayRead(A,(const PetscScalar**)&vv);
1840:   MatDenseGetLDA(A,&lda);
1841:   v    = vv;
1842:   if (type == NORM_FROBENIUS) {
1843:     if (lda>m) {
1844:       for (j=0; j<A->cmap->n; j++) {
1845:         v = vv+j*lda;
1846:         for (i=0; i<m; i++) {
1847:           sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1848:         }
1849:       }
1850:     } else {
1851: #if defined(PETSC_USE_REAL___FP16)
1852:       PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n;
1853:       PetscStackCallBLAS("BLASnrm2",*nrm = BLASnrm2_(&cnt,v,&one));
1854:     }
1855: #else
1856:       for (i=0; i<A->cmap->n*A->rmap->n; i++) {
1857:         sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1858:       }
1859:     }
1860:     *nrm = PetscSqrtReal(sum);
1861: #endif
1862:     PetscLogFlops(2.0*A->cmap->n*A->rmap->n);
1863:   } else if (type == NORM_1) {
1864:     *nrm = 0.0;
1865:     for (j=0; j<A->cmap->n; j++) {
1866:       v   = vv + j*mat->lda;
1867:       sum = 0.0;
1868:       for (i=0; i<A->rmap->n; i++) {
1869:         sum += PetscAbsScalar(*v);  v++;
1870:       }
1871:       if (sum > *nrm) *nrm = sum;
1872:     }
1873:     PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1874:   } else if (type == NORM_INFINITY) {
1875:     *nrm = 0.0;
1876:     for (j=0; j<A->rmap->n; j++) {
1877:       v   = vv + j;
1878:       sum = 0.0;
1879:       for (i=0; i<A->cmap->n; i++) {
1880:         sum += PetscAbsScalar(*v); v += mat->lda;
1881:       }
1882:       if (sum > *nrm) *nrm = sum;
1883:     }
1884:     PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1885:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
1886:   MatDenseRestoreArrayRead(A,(const PetscScalar**)&vv);
1887:   return 0;
1888: }

1890: static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
1891: {
1892:   Mat_SeqDense   *aij = (Mat_SeqDense*)A->data;

1894:   switch (op) {
1895:   case MAT_ROW_ORIENTED:
1896:     aij->roworiented = flg;
1897:     break;
1898:   case MAT_NEW_NONZERO_LOCATIONS:
1899:   case MAT_NEW_NONZERO_LOCATION_ERR:
1900:   case MAT_NEW_NONZERO_ALLOCATION_ERR:
1901:   case MAT_FORCE_DIAGONAL_ENTRIES:
1902:   case MAT_KEEP_NONZERO_PATTERN:
1903:   case MAT_IGNORE_OFF_PROC_ENTRIES:
1904:   case MAT_USE_HASH_TABLE:
1905:   case MAT_IGNORE_ZERO_ENTRIES:
1906:   case MAT_IGNORE_LOWER_TRIANGULAR:
1907:   case MAT_SORTED_FULL:
1908:     PetscInfo(A,"Option %s ignored\n",MatOptions[op]);
1909:     break;
1910:   case MAT_SPD:
1911:   case MAT_SYMMETRIC:
1912:   case MAT_STRUCTURALLY_SYMMETRIC:
1913:   case MAT_HERMITIAN:
1914:   case MAT_SYMMETRY_ETERNAL:
1915:     /* These options are handled directly by MatSetOption() */
1916:     break;
1917:   default:
1918:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]);
1919:   }
1920:   return 0;
1921: }

1923: PetscErrorCode MatZeroEntries_SeqDense(Mat A)
1924: {
1925:   Mat_SeqDense   *l = (Mat_SeqDense*)A->data;
1926:   PetscInt       lda=l->lda,m=A->rmap->n,n=A->cmap->n,j;
1927:   PetscScalar    *v;

1929:   MatDenseGetArrayWrite(A,&v);
1930:   if (lda>m) {
1931:     for (j=0; j<n; j++) {
1932:       PetscArrayzero(v+j*lda,m);
1933:     }
1934:   } else {
1935:     PetscArrayzero(v,PetscInt64Mult(m,n));
1936:   }
1937:   MatDenseRestoreArrayWrite(A,&v);
1938:   return 0;
1939: }

1941: static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1942: {
1943:   Mat_SeqDense      *l = (Mat_SeqDense*)A->data;
1944:   PetscInt          m  = l->lda, n = A->cmap->n, i,j;
1945:   PetscScalar       *slot,*bb,*v;
1946:   const PetscScalar *xx;

1948:   if (PetscDefined(USE_DEBUG)) {
1949:     for (i=0; i<N; i++) {
1952:     }
1953:   }
1954:   if (!N) return 0;

1956:   /* fix right hand side if needed */
1957:   if (x && b) {
1958:     VecGetArrayRead(x,&xx);
1959:     VecGetArray(b,&bb);
1960:     for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
1961:     VecRestoreArrayRead(x,&xx);
1962:     VecRestoreArray(b,&bb);
1963:   }

1965:   MatDenseGetArray(A,&v);
1966:   for (i=0; i<N; i++) {
1967:     slot = v + rows[i];
1968:     for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
1969:   }
1970:   if (diag != 0.0) {
1972:     for (i=0; i<N; i++) {
1973:       slot  = v + (m+1)*rows[i];
1974:       *slot = diag;
1975:     }
1976:   }
1977:   MatDenseRestoreArray(A,&v);
1978:   return 0;
1979: }

1981: static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A,PetscInt *lda)
1982: {
1983:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;

1985:   *lda = mat->lda;
1986:   return 0;
1987: }

1989: PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar **array)
1990: {
1991:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;

1994:   *array = mat->v;
1995:   return 0;
1996: }

1998: PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar **array)
1999: {
2000:   if (array) *array = NULL;
2001:   return 0;
2002: }

2004: /*@
2005:    MatDenseGetLDA - gets the leading dimension of the array returned from MatDenseGetArray()

2007:    Not collective

2009:    Input Parameter:
2010: .  mat - a MATSEQDENSE or MATMPIDENSE matrix

2012:    Output Parameter:
2013: .   lda - the leading dimension

2015:    Level: intermediate

2017: .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseSetLDA()
2018: @*/
2019: PetscErrorCode  MatDenseGetLDA(Mat A,PetscInt *lda)
2020: {
2023:   MatCheckPreallocated(A,1);
2024:   PetscUseMethod(A,"MatDenseGetLDA_C",(Mat,PetscInt*),(A,lda));
2025:   return 0;
2026: }

2028: /*@
2029:    MatDenseSetLDA - Sets the leading dimension of the array used by the dense matrix

2031:    Not collective

2033:    Input Parameters:
2034: +  mat - a MATSEQDENSE or MATMPIDENSE matrix
2035: -  lda - the leading dimension

2037:    Level: intermediate

2039: .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetLDA()
2040: @*/
2041: PetscErrorCode  MatDenseSetLDA(Mat A,PetscInt lda)
2042: {
2044:   PetscTryMethod(A,"MatDenseSetLDA_C",(Mat,PetscInt),(A,lda));
2045:   return 0;
2046: }

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

2051:    Logically Collective on Mat

2053:    Input Parameter:
2054: .  mat - a dense matrix

2056:    Output Parameter:
2057: .   array - pointer to the data

2059:    Level: intermediate

2061: .seealso: MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
2062: @*/
2063: PetscErrorCode  MatDenseGetArray(Mat A,PetscScalar **array)
2064: {
2067:   PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));
2068:   return 0;
2069: }

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

2074:    Logically Collective on Mat

2076:    Input Parameters:
2077: +  mat - a dense matrix
2078: -  array - pointer to the data

2080:    Level: intermediate

2082: .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
2083: @*/
2084: PetscErrorCode  MatDenseRestoreArray(Mat A,PetscScalar **array)
2085: {
2088:   PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));
2089:   PetscObjectStateIncrease((PetscObject)A);
2090: #if defined(PETSC_HAVE_CUDA)
2091:   A->offloadmask = PETSC_OFFLOAD_CPU;
2092: #endif
2093:   return 0;
2094: }

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

2099:    Not Collective

2101:    Input Parameter:
2102: .  mat - a dense matrix

2104:    Output Parameter:
2105: .   array - pointer to the data

2107:    Level: intermediate

2109: .seealso: MatDenseRestoreArrayRead(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
2110: @*/
2111: PetscErrorCode  MatDenseGetArrayRead(Mat A,const PetscScalar **array)
2112: {
2115:   PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array));
2116:   return 0;
2117: }

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

2122:    Not Collective

2124:    Input Parameters:
2125: +  mat - a dense matrix
2126: -  array - pointer to the data

2128:    Level: intermediate

2130: .seealso: MatDenseGetArrayRead(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
2131: @*/
2132: PetscErrorCode  MatDenseRestoreArrayRead(Mat A,const PetscScalar **array)
2133: {
2136:   PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array));
2137:   return 0;
2138: }

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

2143:    Not Collective

2145:    Input Parameter:
2146: .  mat - a dense matrix

2148:    Output Parameter:
2149: .   array - pointer to the data

2151:    Level: intermediate

2153: .seealso: MatDenseRestoreArrayWrite(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead()
2154: @*/
2155: PetscErrorCode  MatDenseGetArrayWrite(Mat A,PetscScalar **array)
2156: {
2159:   PetscUseMethod(A,"MatDenseGetArrayWrite_C",(Mat,PetscScalar**),(A,array));
2160:   return 0;
2161: }

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

2166:    Not Collective

2168:    Input Parameters:
2169: +  mat - a dense matrix
2170: -  array - pointer to the data

2172:    Level: intermediate

2174: .seealso: MatDenseGetArrayWrite(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead()
2175: @*/
2176: PetscErrorCode  MatDenseRestoreArrayWrite(Mat A,PetscScalar **array)
2177: {
2180:   PetscUseMethod(A,"MatDenseRestoreArrayWrite_C",(Mat,PetscScalar**),(A,array));
2181:   PetscObjectStateIncrease((PetscObject)A);
2182: #if defined(PETSC_HAVE_CUDA)
2183:   A->offloadmask = PETSC_OFFLOAD_CPU;
2184: #endif
2185:   return 0;
2186: }

2188: static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,MatReuse scall,Mat *B)
2189: {
2190:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
2191:   PetscInt       i,j,nrows,ncols,ldb;
2192:   const PetscInt *irow,*icol;
2193:   PetscScalar    *av,*bv,*v = mat->v;
2194:   Mat            newmat;

2196:   ISGetIndices(isrow,&irow);
2197:   ISGetIndices(iscol,&icol);
2198:   ISGetLocalSize(isrow,&nrows);
2199:   ISGetLocalSize(iscol,&ncols);

2201:   /* Check submatrixcall */
2202:   if (scall == MAT_REUSE_MATRIX) {
2203:     PetscInt n_cols,n_rows;
2204:     MatGetSize(*B,&n_rows,&n_cols);
2205:     if (n_rows != nrows || n_cols != ncols) {
2206:       /* resize the result matrix to match number of requested rows/columns */
2207:       MatSetSizes(*B,nrows,ncols,nrows,ncols);
2208:     }
2209:     newmat = *B;
2210:   } else {
2211:     /* Create and fill new matrix */
2212:     MatCreate(PetscObjectComm((PetscObject)A),&newmat);
2213:     MatSetSizes(newmat,nrows,ncols,nrows,ncols);
2214:     MatSetType(newmat,((PetscObject)A)->type_name);
2215:     MatSeqDenseSetPreallocation(newmat,NULL);
2216:   }

2218:   /* Now extract the data pointers and do the copy,column at a time */
2219:   MatDenseGetArray(newmat,&bv);
2220:   MatDenseGetLDA(newmat,&ldb);
2221:   for (i=0; i<ncols; i++) {
2222:     av = v + mat->lda*icol[i];
2223:     for (j=0; j<nrows; j++) bv[j] = av[irow[j]];
2224:     bv += ldb;
2225:   }
2226:   MatDenseRestoreArray(newmat,&bv);

2228:   /* Assemble the matrices so that the correct flags are set */
2229:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
2230:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);

2232:   /* Free work space */
2233:   ISRestoreIndices(isrow,&irow);
2234:   ISRestoreIndices(iscol,&icol);
2235:   *B   = newmat;
2236:   return 0;
2237: }

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

2243:   if (scall == MAT_INITIAL_MATRIX) {
2244:     PetscCalloc1(n,B);
2245:   }

2247:   for (i=0; i<n; i++) {
2248:     MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],scall,&(*B)[i]);
2249:   }
2250:   return 0;
2251: }

2253: static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
2254: {
2255:   return 0;
2256: }

2258: static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
2259: {
2260:   return 0;
2261: }

2263: PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
2264: {
2265:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data;
2266:   const PetscScalar *va;
2267:   PetscScalar       *vb;
2268:   PetscInt          lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j;

2270:   /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
2271:   if (A->ops->copy != B->ops->copy) {
2272:     MatCopy_Basic(A,B,str);
2273:     return 0;
2274:   }
2276:   MatDenseGetArrayRead(A,&va);
2277:   MatDenseGetArray(B,&vb);
2278:   if (lda1>m || lda2>m) {
2279:     for (j=0; j<n; j++) {
2280:       PetscArraycpy(vb+j*lda2,va+j*lda1,m);
2281:     }
2282:   } else {
2283:     PetscArraycpy(vb,va,A->rmap->n*A->cmap->n);
2284:   }
2285:   MatDenseRestoreArray(B,&vb);
2286:   MatDenseRestoreArrayRead(A,&va);
2287:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
2288:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
2289:   return 0;
2290: }

2292: PetscErrorCode MatSetUp_SeqDense(Mat A)
2293: {
2294:   PetscLayoutSetUp(A->rmap);
2295:   PetscLayoutSetUp(A->cmap);
2296:   if (!A->preallocated) {
2297:     MatSeqDenseSetPreallocation(A,NULL);
2298:   }
2299:   return 0;
2300: }

2302: static PetscErrorCode MatConjugate_SeqDense(Mat A)
2303: {
2304:   Mat_SeqDense   *mat = (Mat_SeqDense *) A->data;
2305:   PetscInt       i,nz = A->rmap->n*A->cmap->n;
2306:   PetscInt       min = PetscMin(A->rmap->n,A->cmap->n);
2307:   PetscScalar    *aa;

2309:   MatDenseGetArray(A,&aa);
2310:   for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
2311:   MatDenseRestoreArray(A,&aa);
2312:   if (mat->tau) for (i = 0; i < min; i++) mat->tau[i] = PetscConj(mat->tau[i]);
2313:   return 0;
2314: }

2316: static PetscErrorCode MatRealPart_SeqDense(Mat A)
2317: {
2318:   PetscInt       i,nz = A->rmap->n*A->cmap->n;
2319:   PetscScalar    *aa;

2321:   MatDenseGetArray(A,&aa);
2322:   for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
2323:   MatDenseRestoreArray(A,&aa);
2324:   return 0;
2325: }

2327: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
2328: {
2329:   PetscInt       i,nz = A->rmap->n*A->cmap->n;
2330:   PetscScalar    *aa;

2332:   MatDenseGetArray(A,&aa);
2333:   for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
2334:   MatDenseRestoreArray(A,&aa);
2335:   return 0;
2336: }

2338: /* ----------------------------------------------------------------*/
2339: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2340: {
2341:   PetscInt       m=A->rmap->n,n=B->cmap->n;
2342:   PetscBool      cisdense;

2344:   MatSetSizes(C,m,n,m,n);
2345:   PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
2346:   if (!cisdense) {
2347:     PetscBool flg;

2349:     PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2350:     MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
2351:   }
2352:   MatSetUp(C);
2353:   return 0;
2354: }

2356: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2357: {
2358:   Mat_SeqDense       *a=(Mat_SeqDense*)A->data,*b=(Mat_SeqDense*)B->data,*c=(Mat_SeqDense*)C->data;
2359:   PetscBLASInt       m,n,k;
2360:   const PetscScalar *av,*bv;
2361:   PetscScalar       *cv;
2362:   PetscScalar       _DOne=1.0,_DZero=0.0;

2364:   PetscBLASIntCast(C->rmap->n,&m);
2365:   PetscBLASIntCast(C->cmap->n,&n);
2366:   PetscBLASIntCast(A->cmap->n,&k);
2367:   if (!m || !n || !k) return 0;
2368:   MatDenseGetArrayRead(A,&av);
2369:   MatDenseGetArrayRead(B,&bv);
2370:   MatDenseGetArrayWrite(C,&cv);
2371:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2372:   PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2373:   MatDenseRestoreArrayRead(A,&av);
2374:   MatDenseRestoreArrayRead(B,&bv);
2375:   MatDenseRestoreArrayWrite(C,&cv);
2376:   return 0;
2377: }

2379: PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2380: {
2381:   PetscInt       m=A->rmap->n,n=B->rmap->n;
2382:   PetscBool      cisdense;

2384:   MatSetSizes(C,m,n,m,n);
2385:   PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
2386:   if (!cisdense) {
2387:     PetscBool flg;

2389:     PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2390:     MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
2391:   }
2392:   MatSetUp(C);
2393:   return 0;
2394: }

2396: PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2397: {
2398:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
2399:   Mat_SeqDense      *b = (Mat_SeqDense*)B->data;
2400:   Mat_SeqDense      *c = (Mat_SeqDense*)C->data;
2401:   const PetscScalar *av,*bv;
2402:   PetscScalar       *cv;
2403:   PetscBLASInt      m,n,k;
2404:   PetscScalar       _DOne=1.0,_DZero=0.0;

2406:   PetscBLASIntCast(C->rmap->n,&m);
2407:   PetscBLASIntCast(C->cmap->n,&n);
2408:   PetscBLASIntCast(A->cmap->n,&k);
2409:   if (!m || !n || !k) return 0;
2410:   MatDenseGetArrayRead(A,&av);
2411:   MatDenseGetArrayRead(B,&bv);
2412:   MatDenseGetArrayWrite(C,&cv);
2413:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2414:   MatDenseRestoreArrayRead(A,&av);
2415:   MatDenseRestoreArrayRead(B,&bv);
2416:   MatDenseRestoreArrayWrite(C,&cv);
2417:   PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2418:   return 0;
2419: }

2421: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2422: {
2423:   PetscInt       m=A->cmap->n,n=B->cmap->n;
2424:   PetscBool      cisdense;

2426:   MatSetSizes(C,m,n,m,n);
2427:   PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
2428:   if (!cisdense) {
2429:     PetscBool flg;

2431:     PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2432:     MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
2433:   }
2434:   MatSetUp(C);
2435:   return 0;
2436: }

2438: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2439: {
2440:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
2441:   Mat_SeqDense      *b = (Mat_SeqDense*)B->data;
2442:   Mat_SeqDense      *c = (Mat_SeqDense*)C->data;
2443:   const PetscScalar *av,*bv;
2444:   PetscScalar       *cv;
2445:   PetscBLASInt      m,n,k;
2446:   PetscScalar       _DOne=1.0,_DZero=0.0;

2448:   PetscBLASIntCast(C->rmap->n,&m);
2449:   PetscBLASIntCast(C->cmap->n,&n);
2450:   PetscBLASIntCast(A->rmap->n,&k);
2451:   if (!m || !n || !k) return 0;
2452:   MatDenseGetArrayRead(A,&av);
2453:   MatDenseGetArrayRead(B,&bv);
2454:   MatDenseGetArrayWrite(C,&cv);
2455:   PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2456:   MatDenseRestoreArrayRead(A,&av);
2457:   MatDenseRestoreArrayRead(B,&bv);
2458:   MatDenseRestoreArrayWrite(C,&cv);
2459:   PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2460:   return 0;
2461: }

2463: /* ----------------------------------------------- */
2464: static PetscErrorCode MatProductSetFromOptions_SeqDense_AB(Mat C)
2465: {
2466:   C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqDense;
2467:   C->ops->productsymbolic = MatProductSymbolic_AB;
2468:   return 0;
2469: }

2471: static PetscErrorCode MatProductSetFromOptions_SeqDense_AtB(Mat C)
2472: {
2473:   C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_SeqDense_SeqDense;
2474:   C->ops->productsymbolic          = MatProductSymbolic_AtB;
2475:   return 0;
2476: }

2478: static PetscErrorCode MatProductSetFromOptions_SeqDense_ABt(Mat C)
2479: {
2480:   C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqDense_SeqDense;
2481:   C->ops->productsymbolic          = MatProductSymbolic_ABt;
2482:   return 0;
2483: }

2485: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense(Mat C)
2486: {
2487:   Mat_Product    *product = C->product;

2489:   switch (product->type) {
2490:   case MATPRODUCT_AB:
2491:     MatProductSetFromOptions_SeqDense_AB(C);
2492:     break;
2493:   case MATPRODUCT_AtB:
2494:     MatProductSetFromOptions_SeqDense_AtB(C);
2495:     break;
2496:   case MATPRODUCT_ABt:
2497:     MatProductSetFromOptions_SeqDense_ABt(C);
2498:     break;
2499:   default:
2500:     break;
2501:   }
2502:   return 0;
2503: }
2504: /* ----------------------------------------------- */

2506: static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[])
2507: {
2508:   Mat_SeqDense       *a = (Mat_SeqDense*)A->data;
2509:   PetscInt           i,j,m = A->rmap->n,n = A->cmap->n,p;
2510:   PetscScalar        *x;
2511:   const PetscScalar *aa;

2514:   VecGetArray(v,&x);
2515:   VecGetLocalSize(v,&p);
2516:   MatDenseGetArrayRead(A,&aa);
2518:   for (i=0; i<m; i++) {
2519:     x[i] = aa[i]; if (idx) idx[i] = 0;
2520:     for (j=1; j<n; j++) {
2521:       if (PetscRealPart(x[i]) < PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;}
2522:     }
2523:   }
2524:   MatDenseRestoreArrayRead(A,&aa);
2525:   VecRestoreArray(v,&x);
2526:   return 0;
2527: }

2529: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])
2530: {
2531:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
2532:   PetscInt          i,j,m = A->rmap->n,n = A->cmap->n,p;
2533:   PetscScalar       *x;
2534:   PetscReal         atmp;
2535:   const PetscScalar *aa;

2538:   VecGetArray(v,&x);
2539:   VecGetLocalSize(v,&p);
2540:   MatDenseGetArrayRead(A,&aa);
2542:   for (i=0; i<m; i++) {
2543:     x[i] = PetscAbsScalar(aa[i]);
2544:     for (j=1; j<n; j++) {
2545:       atmp = PetscAbsScalar(aa[i+a->lda*j]);
2546:       if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;}
2547:     }
2548:   }
2549:   MatDenseRestoreArrayRead(A,&aa);
2550:   VecRestoreArray(v,&x);
2551:   return 0;
2552: }

2554: static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])
2555: {
2556:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
2557:   PetscInt          i,j,m = A->rmap->n,n = A->cmap->n,p;
2558:   PetscScalar       *x;
2559:   const PetscScalar *aa;

2562:   MatDenseGetArrayRead(A,&aa);
2563:   VecGetArray(v,&x);
2564:   VecGetLocalSize(v,&p);
2566:   for (i=0; i<m; i++) {
2567:     x[i] = aa[i]; if (idx) idx[i] = 0;
2568:     for (j=1; j<n; j++) {
2569:       if (PetscRealPart(x[i]) > PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;}
2570:     }
2571:   }
2572:   VecRestoreArray(v,&x);
2573:   MatDenseRestoreArrayRead(A,&aa);
2574:   return 0;
2575: }

2577: PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
2578: {
2579:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
2580:   PetscScalar       *x;
2581:   const PetscScalar *aa;

2584:   MatDenseGetArrayRead(A,&aa);
2585:   VecGetArray(v,&x);
2586:   PetscArraycpy(x,aa+col*a->lda,A->rmap->n);
2587:   VecRestoreArray(v,&x);
2588:   MatDenseRestoreArrayRead(A,&aa);
2589:   return 0;
2590: }

2592: PETSC_INTERN PetscErrorCode MatGetColumnReductions_SeqDense(Mat A,PetscInt type,PetscReal *reductions)
2593: {
2594:   PetscInt          i,j,m,n;
2595:   const PetscScalar *a;

2597:   MatGetSize(A,&m,&n);
2598:   PetscArrayzero(reductions,n);
2599:   MatDenseGetArrayRead(A,&a);
2600:   if (type == NORM_2) {
2601:     for (i=0; i<n; i++) {
2602:       for (j=0; j<m; j++) {
2603:         reductions[i] += PetscAbsScalar(a[j]*a[j]);
2604:       }
2605:       a += m;
2606:     }
2607:   } else if (type == NORM_1) {
2608:     for (i=0; i<n; i++) {
2609:       for (j=0; j<m; j++) {
2610:         reductions[i] += PetscAbsScalar(a[j]);
2611:       }
2612:       a += m;
2613:     }
2614:   } else if (type == NORM_INFINITY) {
2615:     for (i=0; i<n; i++) {
2616:       for (j=0; j<m; j++) {
2617:         reductions[i] = PetscMax(PetscAbsScalar(a[j]),reductions[i]);
2618:       }
2619:       a += m;
2620:     }
2621:   } else if (type == REDUCTION_SUM_REALPART || type == REDUCTION_MEAN_REALPART) {
2622:     for (i=0; i<n; i++) {
2623:       for (j=0; j<m; j++) {
2624:         reductions[i] += PetscRealPart(a[j]);
2625:       }
2626:       a += m;
2627:     }
2628:   } else if (type == REDUCTION_SUM_IMAGINARYPART || type == REDUCTION_MEAN_IMAGINARYPART) {
2629:     for (i=0; i<n; i++) {
2630:       for (j=0; j<m; j++) {
2631:         reductions[i] += PetscImaginaryPart(a[j]);
2632:       }
2633:       a += m;
2634:     }
2635:   } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown reduction type");
2636:   MatDenseRestoreArrayRead(A,&a);
2637:   if (type == NORM_2) {
2638:     for (i=0; i<n; i++) reductions[i] = PetscSqrtReal(reductions[i]);
2639:   } else if (type == REDUCTION_MEAN_REALPART || type == REDUCTION_MEAN_IMAGINARYPART) {
2640:     for (i=0; i<n; i++) reductions[i] /= m;
2641:   }
2642:   return 0;
2643: }

2645: static PetscErrorCode  MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
2646: {
2647:   PetscScalar    *a;
2648:   PetscInt       lda,m,n,i,j;

2650:   MatGetSize(x,&m,&n);
2651:   MatDenseGetLDA(x,&lda);
2652:   MatDenseGetArray(x,&a);
2653:   for (j=0; j<n; j++) {
2654:     for (i=0; i<m; i++) {
2655:       PetscRandomGetValue(rctx,a+j*lda+i);
2656:     }
2657:   }
2658:   MatDenseRestoreArray(x,&a);
2659:   return 0;
2660: }

2662: static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool  *missing,PetscInt *d)
2663: {
2664:   *missing = PETSC_FALSE;
2665:   return 0;
2666: }

2668: /* vals is not const */
2669: static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals)
2670: {
2671:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
2672:   PetscScalar    *v;

2675:   MatDenseGetArray(A,&v);
2676:   *vals = v+col*a->lda;
2677:   MatDenseRestoreArray(A,&v);
2678:   return 0;
2679: }

2681: static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals)
2682: {
2683:   *vals = NULL; /* user cannot accidentally use the array later */
2684:   return 0;
2685: }

2687: /* -------------------------------------------------------------------*/
2688: static struct _MatOps MatOps_Values = { MatSetValues_SeqDense,
2689:                                         MatGetRow_SeqDense,
2690:                                         MatRestoreRow_SeqDense,
2691:                                         MatMult_SeqDense,
2692:                                 /*  4*/ MatMultAdd_SeqDense,
2693:                                         MatMultTranspose_SeqDense,
2694:                                         MatMultTransposeAdd_SeqDense,
2695:                                         NULL,
2696:                                         NULL,
2697:                                         NULL,
2698:                                 /* 10*/ NULL,
2699:                                         MatLUFactor_SeqDense,
2700:                                         MatCholeskyFactor_SeqDense,
2701:                                         MatSOR_SeqDense,
2702:                                         MatTranspose_SeqDense,
2703:                                 /* 15*/ MatGetInfo_SeqDense,
2704:                                         MatEqual_SeqDense,
2705:                                         MatGetDiagonal_SeqDense,
2706:                                         MatDiagonalScale_SeqDense,
2707:                                         MatNorm_SeqDense,
2708:                                 /* 20*/ MatAssemblyBegin_SeqDense,
2709:                                         MatAssemblyEnd_SeqDense,
2710:                                         MatSetOption_SeqDense,
2711:                                         MatZeroEntries_SeqDense,
2712:                                 /* 24*/ MatZeroRows_SeqDense,
2713:                                         NULL,
2714:                                         NULL,
2715:                                         NULL,
2716:                                         NULL,
2717:                                 /* 29*/ MatSetUp_SeqDense,
2718:                                         NULL,
2719:                                         NULL,
2720:                                         NULL,
2721:                                         NULL,
2722:                                 /* 34*/ MatDuplicate_SeqDense,
2723:                                         NULL,
2724:                                         NULL,
2725:                                         NULL,
2726:                                         NULL,
2727:                                 /* 39*/ MatAXPY_SeqDense,
2728:                                         MatCreateSubMatrices_SeqDense,
2729:                                         NULL,
2730:                                         MatGetValues_SeqDense,
2731:                                         MatCopy_SeqDense,
2732:                                 /* 44*/ MatGetRowMax_SeqDense,
2733:                                         MatScale_SeqDense,
2734:                                         MatShift_Basic,
2735:                                         NULL,
2736:                                         MatZeroRowsColumns_SeqDense,
2737:                                 /* 49*/ MatSetRandom_SeqDense,
2738:                                         NULL,
2739:                                         NULL,
2740:                                         NULL,
2741:                                         NULL,
2742:                                 /* 54*/ NULL,
2743:                                         NULL,
2744:                                         NULL,
2745:                                         NULL,
2746:                                         NULL,
2747:                                 /* 59*/ MatCreateSubMatrix_SeqDense,
2748:                                         MatDestroy_SeqDense,
2749:                                         MatView_SeqDense,
2750:                                         NULL,
2751:                                         NULL,
2752:                                 /* 64*/ NULL,
2753:                                         NULL,
2754:                                         NULL,
2755:                                         NULL,
2756:                                         NULL,
2757:                                 /* 69*/ MatGetRowMaxAbs_SeqDense,
2758:                                         NULL,
2759:                                         NULL,
2760:                                         NULL,
2761:                                         NULL,
2762:                                 /* 74*/ NULL,
2763:                                         NULL,
2764:                                         NULL,
2765:                                         NULL,
2766:                                         NULL,
2767:                                 /* 79*/ NULL,
2768:                                         NULL,
2769:                                         NULL,
2770:                                         NULL,
2771:                                 /* 83*/ MatLoad_SeqDense,
2772:                                         MatIsSymmetric_SeqDense,
2773:                                         MatIsHermitian_SeqDense,
2774:                                         NULL,
2775:                                         NULL,
2776:                                         NULL,
2777:                                 /* 89*/ NULL,
2778:                                         NULL,
2779:                                         MatMatMultNumeric_SeqDense_SeqDense,
2780:                                         NULL,
2781:                                         NULL,
2782:                                 /* 94*/ NULL,
2783:                                         NULL,
2784:                                         NULL,
2785:                                         MatMatTransposeMultNumeric_SeqDense_SeqDense,
2786:                                         NULL,
2787:                                 /* 99*/ MatProductSetFromOptions_SeqDense,
2788:                                         NULL,
2789:                                         NULL,
2790:                                         MatConjugate_SeqDense,
2791:                                         NULL,
2792:                                 /*104*/ NULL,
2793:                                         MatRealPart_SeqDense,
2794:                                         MatImaginaryPart_SeqDense,
2795:                                         NULL,
2796:                                         NULL,
2797:                                 /*109*/ NULL,
2798:                                         NULL,
2799:                                         MatGetRowMin_SeqDense,
2800:                                         MatGetColumnVector_SeqDense,
2801:                                         MatMissingDiagonal_SeqDense,
2802:                                 /*114*/ NULL,
2803:                                         NULL,
2804:                                         NULL,
2805:                                         NULL,
2806:                                         NULL,
2807:                                 /*119*/ NULL,
2808:                                         NULL,
2809:                                         NULL,
2810:                                         NULL,
2811:                                         NULL,
2812:                                 /*124*/ NULL,
2813:                                         MatGetColumnReductions_SeqDense,
2814:                                         NULL,
2815:                                         NULL,
2816:                                         NULL,
2817:                                 /*129*/ NULL,
2818:                                         NULL,
2819:                                         NULL,
2820:                                         MatTransposeMatMultNumeric_SeqDense_SeqDense,
2821:                                         NULL,
2822:                                 /*134*/ NULL,
2823:                                         NULL,
2824:                                         NULL,
2825:                                         NULL,
2826:                                         NULL,
2827:                                 /*139*/ NULL,
2828:                                         NULL,
2829:                                         NULL,
2830:                                         NULL,
2831:                                         NULL,
2832:                                         MatCreateMPIMatConcatenateSeqMat_SeqDense,
2833:                                 /*145*/ NULL,
2834:                                         NULL,
2835:                                         NULL
2836: };

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

2843:    Collective

2845:    Input Parameters:
2846: +  comm - MPI communicator, set to PETSC_COMM_SELF
2847: .  m - number of rows
2848: .  n - number of columns
2849: -  data - optional location of matrix data in column major order.  Set data=NULL for PETSc
2850:    to control all matrix memory allocation.

2852:    Output Parameter:
2853: .  A - the matrix

2855:    Notes:
2856:    The data input variable is intended primarily for Fortran programmers
2857:    who wish to allocate their own matrix memory space.  Most users should
2858:    set data=NULL.

2860:    Level: intermediate

2862: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
2863: @*/
2864: PetscErrorCode  MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
2865: {
2866:   MatCreate(comm,A);
2867:   MatSetSizes(*A,m,n,m,n);
2868:   MatSetType(*A,MATSEQDENSE);
2869:   MatSeqDenseSetPreallocation(*A,data);
2870:   return 0;
2871: }

2873: /*@C
2874:    MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements

2876:    Collective

2878:    Input Parameters:
2879: +  B - the matrix
2880: -  data - the array (or NULL)

2882:    Notes:
2883:    The data input variable is intended primarily for Fortran programmers
2884:    who wish to allocate their own matrix memory space.  Most users should
2885:    need not call this routine.

2887:    Level: intermediate

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

2891: @*/
2892: PetscErrorCode  MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
2893: {
2895:   PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
2896:   return 0;
2897: }

2899: PetscErrorCode  MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
2900: {
2901:   Mat_SeqDense   *b = (Mat_SeqDense*)B->data;

2904:   B->preallocated = PETSC_TRUE;

2906:   PetscLayoutSetUp(B->rmap);
2907:   PetscLayoutSetUp(B->cmap);

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

2911:   if (!data) { /* petsc-allocated storage */
2912:     if (!b->user_alloc) PetscFree(b->v);
2913:     PetscCalloc1((size_t)b->lda*B->cmap->n,&b->v);
2914:     PetscLogObjectMemory((PetscObject)B,b->lda*B->cmap->n*sizeof(PetscScalar));

2916:     b->user_alloc = PETSC_FALSE;
2917:   } else { /* user-allocated storage */
2918:     if (!b->user_alloc) PetscFree(b->v);
2919:     b->v          = data;
2920:     b->user_alloc = PETSC_TRUE;
2921:   }
2922:   B->assembled = PETSC_TRUE;
2923:   return 0;
2924: }

2926: #if defined(PETSC_HAVE_ELEMENTAL)
2927: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
2928: {
2929:   Mat               mat_elemental;
2930:   const PetscScalar *array;
2931:   PetscScalar       *v_colwise;
2932:   PetscInt          M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols;

2934:   PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);
2935:   MatDenseGetArrayRead(A,&array);
2936:   /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
2937:   k = 0;
2938:   for (j=0; j<N; j++) {
2939:     cols[j] = j;
2940:     for (i=0; i<M; i++) {
2941:       v_colwise[j*M+i] = array[k++];
2942:     }
2943:   }
2944:   for (i=0; i<M; i++) {
2945:     rows[i] = i;
2946:   }
2947:   MatDenseRestoreArrayRead(A,&array);

2949:   MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
2950:   MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);
2951:   MatSetType(mat_elemental,MATELEMENTAL);
2952:   MatSetUp(mat_elemental);

2954:   /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
2955:   MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);
2956:   MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
2957:   MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);
2958:   PetscFree3(v_colwise,rows,cols);

2960:   if (reuse == MAT_INPLACE_MATRIX) {
2961:     MatHeaderReplace(A,&mat_elemental);
2962:   } else {
2963:     *newmat = mat_elemental;
2964:   }
2965:   return 0;
2966: }
2967: #endif

2969: PetscErrorCode  MatDenseSetLDA_SeqDense(Mat B,PetscInt lda)
2970: {
2971:   Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2972:   PetscBool    data;

2974:   data = (PetscBool)((B->rmap->n > 0 && B->cmap->n > 0) ? (b->v ? PETSC_TRUE : PETSC_FALSE) : PETSC_FALSE);
2977:   b->lda = lda;
2978:   return 0;
2979: }

2981: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat)
2982: {
2983:   PetscMPIInt    size;

2985:   MPI_Comm_size(comm,&size);
2986:   if (size == 1) {
2987:     if (scall == MAT_INITIAL_MATRIX) {
2988:       MatDuplicate(inmat,MAT_COPY_VALUES,outmat);
2989:     } else {
2990:       MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);
2991:     }
2992:   } else {
2993:     MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);
2994:   }
2995:   return 0;
2996: }

2998: PetscErrorCode MatDenseGetColumnVec_SeqDense(Mat A,PetscInt col,Vec *v)
2999: {
3000:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

3004:   if (!a->cvec) {
3005:     VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);
3006:     PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);
3007:   }
3008:   a->vecinuse = col + 1;
3009:   MatDenseGetArray(A,(PetscScalar**)&a->ptrinuse);
3010:   VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);
3011:   *v   = a->cvec;
3012:   return 0;
3013: }

3015: PetscErrorCode MatDenseRestoreColumnVec_SeqDense(Mat A,PetscInt col,Vec *v)
3016: {
3017:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

3021:   a->vecinuse = 0;
3022:   MatDenseRestoreArray(A,(PetscScalar**)&a->ptrinuse);
3023:   VecResetArray(a->cvec);
3024:   if (v) *v = NULL;
3025:   return 0;
3026: }

3028: PetscErrorCode MatDenseGetColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v)
3029: {
3030:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

3034:   if (!a->cvec) {
3035:     VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);
3036:     PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);
3037:   }
3038:   a->vecinuse = col + 1;
3039:   MatDenseGetArrayRead(A,&a->ptrinuse);
3040:   VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);
3041:   VecLockReadPush(a->cvec);
3042:   *v   = a->cvec;
3043:   return 0;
3044: }

3046: PetscErrorCode MatDenseRestoreColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v)
3047: {
3048:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

3052:   a->vecinuse = 0;
3053:   MatDenseRestoreArrayRead(A,&a->ptrinuse);
3054:   VecLockReadPop(a->cvec);
3055:   VecResetArray(a->cvec);
3056:   if (v) *v = NULL;
3057:   return 0;
3058: }

3060: PetscErrorCode MatDenseGetColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v)
3061: {
3062:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

3066:   if (!a->cvec) {
3067:     VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);
3068:     PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);
3069:   }
3070:   a->vecinuse = col + 1;
3071:   MatDenseGetArrayWrite(A,(PetscScalar**)&a->ptrinuse);
3072:   VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);
3073:   *v   = a->cvec;
3074:   return 0;
3075: }

3077: PetscErrorCode MatDenseRestoreColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v)
3078: {
3079:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

3083:   a->vecinuse = 0;
3084:   MatDenseRestoreArrayWrite(A,(PetscScalar**)&a->ptrinuse);
3085:   VecResetArray(a->cvec);
3086:   if (v) *v = NULL;
3087:   return 0;
3088: }

3090: PetscErrorCode MatDenseGetSubMatrix_SeqDense(Mat A,PetscInt cbegin,PetscInt cend,Mat *v)
3091: {
3092:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

3096:   if (a->cmat && cend-cbegin != a->cmat->cmap->N) {
3097:     MatDestroy(&a->cmat);
3098:   }
3099:   if (!a->cmat) {
3100:     MatCreateDense(PetscObjectComm((PetscObject)A),A->rmap->n,PETSC_DECIDE,A->rmap->N,cend-cbegin,a->v+(size_t)cbegin*a->lda,&a->cmat);
3101:     PetscLogObjectParent((PetscObject)A,(PetscObject)a->cmat);
3102:   } else {
3103:     MatDensePlaceArray(a->cmat,a->v+(size_t)cbegin*a->lda);
3104:   }
3105:   MatDenseSetLDA(a->cmat,a->lda);
3106:   a->matinuse = cbegin + 1;
3107:   *v = a->cmat;
3108: #if defined(PETSC_HAVE_CUDA)
3109:   A->offloadmask = PETSC_OFFLOAD_CPU;
3110: #endif
3111:   return 0;
3112: }

3114: PetscErrorCode MatDenseRestoreSubMatrix_SeqDense(Mat A,Mat *v)
3115: {
3116:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

3121:   a->matinuse = 0;
3122:   MatDenseResetArray(a->cmat);
3123:   *v   = NULL;
3124:   return 0;
3125: }

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

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

3133:   Level: beginner

3135: .seealso: MatCreateSeqDense()

3137: M*/
3138: PetscErrorCode MatCreate_SeqDense(Mat B)
3139: {
3140:   Mat_SeqDense   *b;
3141:   PetscMPIInt    size;

3143:   MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);

3146:   PetscNewLog(B,&b);
3147:   PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
3148:   B->data = (void*)b;

3150:   b->roworiented = PETSC_TRUE;

3152:   PetscObjectComposeFunction((PetscObject)B,"MatQRFactor_C",MatQRFactor_SeqDense);
3153:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetLDA_C",MatDenseGetLDA_SeqDense);
3154:   PetscObjectComposeFunction((PetscObject)B,"MatDenseSetLDA_C",MatDenseSetLDA_SeqDense);
3155:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);
3156:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);
3157:   PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);
3158:   PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);
3159:   PetscObjectComposeFunction((PetscObject)B,"MatDenseReplaceArray_C",MatDenseReplaceArray_SeqDense);
3160:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);
3161:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense);
3162:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayWrite_C",MatDenseGetArray_SeqDense);
3163:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayWrite_C",MatDenseRestoreArray_SeqDense);
3164:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);
3165: #if defined(PETSC_HAVE_ELEMENTAL)
3166:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);
3167: #endif
3168: #if defined(PETSC_HAVE_SCALAPACK)
3169:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_scalapack_C",MatConvert_Dense_ScaLAPACK);
3170: #endif
3171: #if defined(PETSC_HAVE_CUDA)
3172:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqdensecuda_C",MatConvert_SeqDense_SeqDenseCUDA);
3173:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",MatProductSetFromOptions_SeqDense);
3174:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdense_C",MatProductSetFromOptions_SeqDense);
3175:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdensecuda_C",MatProductSetFromOptions_SeqDense);
3176: #endif
3177:   PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);
3178:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqdense_C",MatProductSetFromOptions_SeqAIJ_SeqDense);
3179:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdense_C",MatProductSetFromOptions_SeqDense);
3180:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);
3181:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqsbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);

3183:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);
3184:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);
3185:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVec_C",MatDenseGetColumnVec_SeqDense);
3186:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVec_C",MatDenseRestoreColumnVec_SeqDense);
3187:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecRead_C",MatDenseGetColumnVecRead_SeqDense);
3188:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecRead_C",MatDenseRestoreColumnVecRead_SeqDense);
3189:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecWrite_C",MatDenseGetColumnVecWrite_SeqDense);
3190:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecWrite_C",MatDenseRestoreColumnVecWrite_SeqDense);
3191:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetSubMatrix_C",MatDenseGetSubMatrix_SeqDense);
3192:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreSubMatrix_C",MatDenseRestoreSubMatrix_SeqDense);
3193:   PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
3194:   return 0;
3195: }

3197: /*@C
3198:    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.

3200:    Not Collective

3202:    Input Parameters:
3203: +  mat - a MATSEQDENSE or MATMPIDENSE matrix
3204: -  col - column index

3206:    Output Parameter:
3207: .  vals - pointer to the data

3209:    Level: intermediate

3211: .seealso: MatDenseRestoreColumn()
3212: @*/
3213: PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals)
3214: {
3218:   PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));
3219:   return 0;
3220: }

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

3225:    Not Collective

3227:    Input Parameter:
3228: .  mat - a MATSEQDENSE or MATMPIDENSE matrix

3230:    Output Parameter:
3231: .  vals - pointer to the data

3233:    Level: intermediate

3235: .seealso: MatDenseGetColumn()
3236: @*/
3237: PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals)
3238: {
3241:   PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));
3242:   return 0;
3243: }

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

3248:    Collective

3250:    Input Parameters:
3251: +  mat - the Mat object
3252: -  col - the column index

3254:    Output Parameter:
3255: .  v - the vector

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

3261:    Level: intermediate

3263: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3264: @*/
3265: PetscErrorCode MatDenseGetColumnVec(Mat A,PetscInt col,Vec *v)
3266: {
3273:   PetscUseMethod(A,"MatDenseGetColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));
3274:   return 0;
3275: }

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

3280:    Collective

3282:    Input Parameters:
3283: +  mat - the Mat object
3284: .  col - the column index
3285: -  v - the Vec object

3287:    Level: intermediate

3289: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3290: @*/
3291: PetscErrorCode MatDenseRestoreColumnVec(Mat A,PetscInt col,Vec *v)
3292: {
3298:   PetscUseMethod(A,"MatDenseRestoreColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));
3299:   return 0;
3300: }

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

3305:    Collective

3307:    Input Parameters:
3308: +  mat - the Mat object
3309: -  col - the column index

3311:    Output Parameter:
3312: .  v - the vector

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

3319:    Level: intermediate

3321: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3322: @*/
3323: PetscErrorCode MatDenseGetColumnVecRead(Mat A,PetscInt col,Vec *v)
3324: {
3331:   PetscUseMethod(A,"MatDenseGetColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));
3332:   return 0;
3333: }

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

3338:    Collective

3340:    Input Parameters:
3341: +  mat - the Mat object
3342: .  col - the column index
3343: -  v - the Vec object

3345:    Level: intermediate

3347: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecWrite()
3348: @*/
3349: PetscErrorCode MatDenseRestoreColumnVecRead(Mat A,PetscInt col,Vec *v)
3350: {
3356:   PetscUseMethod(A,"MatDenseRestoreColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));
3357:   return 0;
3358: }

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

3363:    Collective

3365:    Input Parameters:
3366: +  mat - the Mat object
3367: -  col - the column index

3369:    Output Parameter:
3370: .  v - the vector

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

3376:    Level: intermediate

3378: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3379: @*/
3380: PetscErrorCode MatDenseGetColumnVecWrite(Mat A,PetscInt col,Vec *v)
3381: {
3388:   PetscUseMethod(A,"MatDenseGetColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));
3389:   return 0;
3390: }

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

3395:    Collective

3397:    Input Parameters:
3398: +  mat - the Mat object
3399: .  col - the column index
3400: -  v - the Vec object

3402:    Level: intermediate

3404: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead()
3405: @*/
3406: PetscErrorCode MatDenseRestoreColumnVecWrite(Mat A,PetscInt col,Vec *v)
3407: {
3413:   PetscUseMethod(A,"MatDenseRestoreColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));
3414:   return 0;
3415: }

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

3420:    Collective

3422:    Input Parameters:
3423: +  mat - the Mat object
3424: .  cbegin - the first index in the block
3425: -  cend - the last index in the block

3427:    Output Parameter:
3428: .  v - the matrix

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

3433:    Level: intermediate

3435: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseRestoreColumnVec(), MatDenseRestoreSubMatrix()
3436: @*/
3437: PetscErrorCode MatDenseGetSubMatrix(Mat A,PetscInt cbegin,PetscInt cend,Mat *v)
3438: {
3447:   PetscUseMethod(A,"MatDenseGetSubMatrix_C",(Mat,PetscInt,PetscInt,Mat*),(A,cbegin,cend,v));
3448:   return 0;
3449: }

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

3454:    Collective

3456:    Input Parameters:
3457: +  mat - the Mat object
3458: -  v - the Mat object

3460:    Level: intermediate

3462: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseRestoreColumnVec(), MatDenseGetSubMatrix()
3463: @*/
3464: PetscErrorCode MatDenseRestoreSubMatrix(Mat A,Mat *v)
3465: {
3469:   PetscUseMethod(A,"MatDenseRestoreSubMatrix_C",(Mat,Mat*),(A,v));
3470:   return 0;
3471: }