Actual source code: dense.c

      2: /*
      3:      Defines the basic matrix operations for sequential dense.
      4: */

      6: #include <../src/mat/impls/dense/seq/dense.h>
      7: #include <petscblaslapack.h>

      9: #include <../src/mat/impls/aij/seq/aij.h>

     11: PetscErrorCode MatSeqDenseSymmetrize_Private(Mat A, PetscBool hermitian)
     12: {
     13:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
     14:   PetscInt       j, k, n = A->rmap->n;
     15:   PetscScalar    *v;

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

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

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

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

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

    105:   if (PetscDefined(USE_DEBUG)) {
    106:     for (i=0; i<N; i++) {
    107:       if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
    108:       if (rows[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D requested to be zeroed greater than or equal number of rows %D",rows[i],A->rmap->n);
    109:       if (rows[i] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Col %D requested to be zeroed greater than or equal number of cols %D",rows[i],A->cmap->n);
    110:     }
    111:   }
    112:   if (!N) return(0);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    456: static PetscErrorCode MatConjugate_SeqDense(Mat);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    669: static PetscErrorCode MatMatSolve_SeqDense_SetUp(Mat A, Mat B, Mat X, PetscScalar **_y, PetscBLASInt *_ylda, PetscBLASInt *_m, PetscBLASInt *_nrhs, PetscBLASInt *_k)
    670: {
    671:   PetscErrorCode    ierr;
    672:   const PetscScalar *b;
    673:   PetscScalar       *y;
    674:   PetscInt          n, _blda, _xlda;
    675:   PetscBLASInt      nrhs=0,m=0,k=0,blda=0,xlda=0,ylda=0;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    828: static PetscErrorCode MatConjugate_SeqDense(Mat);

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

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

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

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

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

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

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

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

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

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

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

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

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

    947:   A->ops->solve             = MatSolve_SeqDense_Cholesky;
    948:   A->ops->matsolve          = MatMatSolve_SeqDense_Cholesky;
    949:   A->ops->solvetranspose    = MatSolveTranspose_SeqDense_Cholesky;
    950:   A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_Cholesky;
    951:   A->factortype             = MAT_FACTOR_CHOLESKY;

    953:   PetscFree(A->solvertype);
    954:   PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);

    956:   PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
    957:   return(0);
    958: }

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

    966:   info.fill = 1.0;

    968:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
    969:   (*fact->ops->choleskyfactor)(fact,NULL,&info);
    970:   return(0);
    971: }

    973: PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
    974: {
    976:   fact->assembled                  = PETSC_TRUE;
    977:   fact->preallocated               = PETSC_TRUE;
    978:   fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
    979:   fact->ops->solve                 = MatSolve_SeqDense_Cholesky;
    980:   fact->ops->matsolve              = MatMatSolve_SeqDense_Cholesky;
    981:   fact->ops->solvetranspose        = MatSolve_SeqDense_Cholesky;
    982:   return(0);
    983: }

    985: static PetscErrorCode MatQRFactor_SeqDense(Mat A,IS col,const MatFactorInfo *minfo)
    986: {
    987:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
    989:   PetscBLASInt   n,m,info, min, max;

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

    1011:     mat->lfwork = -1;
    1012:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
    1013:     PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&m,&n,mat->v,&mat->lda,mat->tau,&dummy,&mat->lfwork,&info));
    1014:     PetscFPTrapPop();
    1015:     mat->lfwork = (PetscInt)PetscRealPart(dummy);
    1016:     PetscMalloc1(mat->lfwork,&mat->fwork);
    1017:     PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
    1018:   }
    1019:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
    1020:   PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&m,&n,mat->v,&mat->lda,mat->tau,mat->fwork,&mat->lfwork,&info));
    1021:   PetscFPTrapPop();
    1022:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to QR factorization");
    1023:   // 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
    1024:   mat->rank = min;

    1026:   A->ops->solve             = MatSolve_SeqDense_QR;
    1027:   A->ops->matsolve          = MatMatSolve_SeqDense_QR;
    1028:   A->factortype             = MAT_FACTOR_QR;
    1029:   if (m == n) {
    1030:     A->ops->solvetranspose    = MatSolveTranspose_SeqDense_QR;
    1031:     A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_QR;
    1032:   }

    1034:   PetscFree(A->solvertype);
    1035:   PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);

    1037:   PetscLogFlops(2.0*min*min*(max-min/3.0));
    1038:   return(0);
    1039: }

    1041: static PetscErrorCode MatQRFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
    1042: {
    1044:   MatFactorInfo  info;

    1047:   info.fill = 1.0;

    1049:   MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
    1050:   MatQRFactor_SeqDense(fact,NULL,&info);
    1051:   return(0);
    1052: }

    1054: static PetscErrorCode MatQRFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
    1055: {

    1059:   fact->assembled                  = PETSC_TRUE;
    1060:   fact->preallocated               = PETSC_TRUE;
    1061:   fact->ops->solve                 = MatSolve_SeqDense_QR;
    1062:   fact->ops->matsolve              = MatMatSolve_SeqDense_QR;
    1063:   if (A->cmap->n == A->rmap->n) {
    1064:     fact->ops->solvetranspose    = MatSolveTranspose_SeqDense_QR;
    1065:     fact->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_QR;
    1066:   }
    1067:   PetscObjectComposeFunction((PetscObject)fact,"MatQRFactorNumeric_C",MatQRFactorNumeric_SeqDense);
    1068:   return(0);
    1069: }

    1072: /* uses LAPACK */
    1073: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact)
    1074: {

    1078:   MatCreate(PetscObjectComm((PetscObject)A),fact);
    1079:   MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
    1080:   MatSetType(*fact,MATDENSE);
    1081:   if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
    1082:     (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
    1083:     (*fact)->ops->ilufactorsymbolic = MatLUFactorSymbolic_SeqDense;
    1084:   } else {
    1085:     (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
    1086:   }
    1087:   (*fact)->factortype = ftype;

    1089:   PetscFree((*fact)->solvertype);
    1090:   PetscStrallocpy(MATSOLVERPETSC,&(*fact)->solvertype);
    1091:   return(0);
    1092: }

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

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

    1137: /* -----------------------------------------------------------------*/
    1138: PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)
    1139: {
    1140:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
    1141:   const PetscScalar *v   = mat->v,*x;
    1142:   PetscScalar       *y;
    1143:   PetscErrorCode    ierr;
    1144:   PetscBLASInt      m, n,_One=1;
    1145:   PetscScalar       _DOne=1.0,_DZero=0.0;

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

    1164: PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy)
    1165: {
    1166:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
    1167:   PetscScalar       *y,_DOne=1.0,_DZero=0.0;
    1168:   PetscErrorCode    ierr;
    1169:   PetscBLASInt      m, n, _One=1;
    1170:   const PetscScalar *v = mat->v,*x;

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

    1189: PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
    1190: {
    1191:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
    1192:   const PetscScalar *v = mat->v,*x;
    1193:   PetscScalar       *y,_DOne=1.0;
    1194:   PetscErrorCode    ierr;
    1195:   PetscBLASInt      m, n, _One=1;

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

    1211: PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
    1212: {
    1213:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
    1214:   const PetscScalar *v = mat->v,*x;
    1215:   PetscScalar       *y;
    1216:   PetscErrorCode    ierr;
    1217:   PetscBLASInt      m, n, _One=1;
    1218:   PetscScalar       _DOne=1.0;

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

    1234: /* -----------------------------------------------------------------*/
    1235: static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
    1236: {
    1237:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
    1239:   PetscInt       i;

    1242:   *ncols = A->cmap->n;
    1243:   if (cols) {
    1244:     PetscMalloc1(A->cmap->n,cols);
    1245:     for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
    1246:   }
    1247:   if (vals) {
    1248:     const PetscScalar *v;

    1250:     MatDenseGetArrayRead(A,&v);
    1251:     PetscMalloc1(A->cmap->n,vals);
    1252:     v   += row;
    1253:     for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
    1254:     MatDenseRestoreArrayRead(A,&v);
    1255:   }
    1256:   return(0);
    1257: }

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

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

    1281:   MatDenseGetArray(A,&av);
    1282:   if (!mat->roworiented) {
    1283:     if (addv == INSERT_VALUES) {
    1284:       for (j=0; j<n; j++) {
    1285:         if (indexn[j] < 0) {idx += m; continue;}
    1286:         if (PetscUnlikelyDebug(indexn[j] >= A->cmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1);
    1287:         for (i=0; i<m; i++) {
    1288:           if (indexm[i] < 0) {idx++; continue;}
    1289:           if (PetscUnlikelyDebug(indexm[i] >= A->rmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1);
    1290:           av[indexn[j]*mat->lda + indexm[i]] = v[idx++];
    1291:         }
    1292:       }
    1293:     } else {
    1294:       for (j=0; j<n; j++) {
    1295:         if (indexn[j] < 0) {idx += m; continue;}
    1296:         if (PetscUnlikelyDebug(indexn[j] >= A->cmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1);
    1297:         for (i=0; i<m; i++) {
    1298:           if (indexm[i] < 0) {idx++; continue;}
    1299:           if (PetscUnlikelyDebug(indexm[i] >= A->rmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1);
    1300:           av[indexn[j]*mat->lda + indexm[i]] += v[idx++];
    1301:         }
    1302:       }
    1303:     }
    1304:   } else {
    1305:     if (addv == INSERT_VALUES) {
    1306:       for (i=0; i<m; i++) {
    1307:         if (indexm[i] < 0) { idx += n; continue;}
    1308:         if (PetscUnlikelyDebug(indexm[i] >= A->rmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1);
    1309:         for (j=0; j<n; j++) {
    1310:           if (indexn[j] < 0) { idx++; continue;}
    1311:           if (PetscUnlikelyDebug(indexn[j] >= A->cmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1);
    1312:           av[indexn[j]*mat->lda + indexm[i]] = v[idx++];
    1313:         }
    1314:       }
    1315:     } else {
    1316:       for (i=0; i<m; i++) {
    1317:         if (indexm[i] < 0) { idx += n; continue;}
    1318:         if (PetscUnlikelyDebug(indexm[i] >= A->rmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1);
    1319:         for (j=0; j<n; j++) {
    1320:           if (indexn[j] < 0) { idx++; continue;}
    1321:           if (PetscUnlikelyDebug(indexn[j] >= A->cmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1);
    1322:           av[indexn[j]*mat->lda + indexm[i]] += v[idx++];
    1323:         }
    1324:       }
    1325:     }
    1326:   }
    1327:   /* hack to prevent unneeded copy to the GPU while returning the array */
    1328: #if defined(PETSC_HAVE_CUDA)
    1329:   oldf = A->offloadmask;
    1330:   A->offloadmask = PETSC_OFFLOAD_GPU;
    1331: #endif
    1332:   MatDenseRestoreArray(A,&av);
    1333: #if defined(PETSC_HAVE_CUDA)
    1334:   A->offloadmask = (oldf == PETSC_OFFLOAD_UNALLOCATED ? PETSC_OFFLOAD_UNALLOCATED : PETSC_OFFLOAD_CPU);
    1335: #endif
    1336:   return(0);
    1337: }

    1339: static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
    1340: {
    1341:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
    1342:   const PetscScalar *vv;
    1343:   PetscInt          i,j;
    1344:   PetscErrorCode    ierr;

    1347:   MatDenseGetArrayRead(A,&vv);
    1348:   /* row-oriented output */
    1349:   for (i=0; i<m; i++) {
    1350:     if (indexm[i] < 0) {v += n;continue;}
    1351:     if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D requested larger than number rows %D",indexm[i],A->rmap->n);
    1352:     for (j=0; j<n; j++) {
    1353:       if (indexn[j] < 0) {v++; continue;}
    1354:       if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column %D requested larger than number columns %D",indexn[j],A->cmap->n);
    1355:       *v++ = vv[indexn[j]*mat->lda + indexm[i]];
    1356:     }
    1357:   }
    1358:   MatDenseRestoreArrayRead(A,&vv);
    1359:   return(0);
    1360: }

    1362: /* -----------------------------------------------------------------*/

    1364: PetscErrorCode MatView_Dense_Binary(Mat mat,PetscViewer viewer)
    1365: {
    1366:   PetscErrorCode    ierr;
    1367:   PetscBool         skipHeader;
    1368:   PetscViewerFormat format;
    1369:   PetscInt          header[4],M,N,m,lda,i,j,k;
    1370:   const PetscScalar *v;
    1371:   PetscScalar       *vwork;

    1374:   PetscViewerSetUp(viewer);
    1375:   PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
    1376:   PetscViewerGetFormat(viewer,&format);
    1377:   if (skipHeader) format = PETSC_VIEWER_NATIVE;

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

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

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

    1413: PetscErrorCode MatLoad_Dense_Binary(Mat mat,PetscViewer viewer)
    1414: {
    1416:   PetscBool      skipHeader;
    1417:   PetscInt       header[4],M,N,m,nz,lda,i,j,k;
    1418:   PetscInt       rows,cols;
    1419:   PetscScalar    *v,*vwork;

    1422:   PetscViewerSetUp(viewer);
    1423:   PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);

    1425:   if (!skipHeader) {
    1426:     PetscViewerBinaryRead(viewer,header,4,NULL,PETSC_INT);
    1427:     if (header[0] != MAT_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Not a matrix object in file");
    1428:     M = header[1]; N = header[2];
    1429:     if (M < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix row size (%D) in file is negative",M);
    1430:     if (N < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix column size (%D) in file is negative",N);
    1431:     nz = header[3];
    1432:     if (nz != MATRIX_BINARY_FORMAT_DENSE && nz < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Unknown matrix format %D in file",nz);
    1433:   } else {
    1434:     MatGetSize(mat,&M,&N);
    1435:     if (M < 0 || N < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Matrix binary file header was skipped, thus the user must specify the global sizes of input matrix");
    1436:     nz = MATRIX_BINARY_FORMAT_DENSE;
    1437:   }

    1439:   /* setup global sizes if not set */
    1440:   if (mat->rmap->N < 0) mat->rmap->N = M;
    1441:   if (mat->cmap->N < 0) mat->cmap->N = N;
    1442:   MatSetUp(mat);
    1443:   /* check if global sizes are correct */
    1444:   MatGetSize(mat,&rows,&cols);
    1445:   if (M != rows || N != cols) SETERRQ4(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different sizes (%d, %d) than the input matrix (%d, %d)",M,N,rows,cols);

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

    1484: PetscErrorCode MatLoad_SeqDense(Mat newMat, PetscViewer viewer)
    1485: {
    1486:   PetscBool      isbinary, ishdf5;

    1492:   /* force binary viewer to load .info file if it has not yet done so */
    1493:   PetscViewerSetUp(viewer);
    1494:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
    1495:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,  &ishdf5);
    1496:   if (isbinary) {
    1497:     MatLoad_Dense_Binary(newMat,viewer);
    1498:   } else if (ishdf5) {
    1499: #if defined(PETSC_HAVE_HDF5)
    1500:     MatLoad_Dense_HDF5(newMat,viewer);
    1501: #else
    1502:     SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
    1503: #endif
    1504:   } else {
    1505:     SETERRQ2(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"Viewer type %s not yet supported for reading %s matrices",((PetscObject)viewer)->type_name,((PetscObject)newMat)->type_name);
    1506:   }
    1507:   return(0);
    1508: }

    1510: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
    1511: {
    1512:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
    1513:   PetscErrorCode    ierr;
    1514:   PetscInt          i,j;
    1515:   const char        *name;
    1516:   PetscScalar       *v,*av;
    1517:   PetscViewerFormat format;
    1518: #if defined(PETSC_USE_COMPLEX)
    1519:   PetscBool         allreal = PETSC_TRUE;
    1520: #endif

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

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

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

    1604:   PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
    1605:   PetscViewerGetFormat(viewer,&format);
    1606:   PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);

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

    1631:     for (i=0; i < m*n; i++) {
    1632:       if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
    1633:     }
    1634:     if (minv >= maxv) maxv = minv + PETSC_SMALL;
    1635:     PetscDrawGetPopup(draw,&popup);
    1636:     PetscDrawScalePopup(popup,minv,maxv);

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

    1655: static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
    1656: {
    1657:   PetscDraw      draw;
    1658:   PetscBool      isnull;
    1659:   PetscReal      xr,yr,xl,yl,h,w;

    1663:   PetscViewerDrawGetDraw(viewer,0,&draw);
    1664:   PetscDrawIsNull(draw,&isnull);
    1665:   if (isnull) return(0);

    1667:   xr   = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
    1668:   xr  += w;          yr += h;        xl = -w;     yl = -h;
    1669:   PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
    1670:   PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
    1671:   PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
    1672:   PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
    1673:   PetscDrawSave(draw);
    1674:   return(0);
    1675: }

    1677: PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
    1678: {
    1680:   PetscBool      iascii,isbinary,isdraw;

    1683:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
    1684:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
    1685:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);

    1687:   if (iascii) {
    1688:     MatView_SeqDense_ASCII(A,viewer);
    1689:   } else if (isbinary) {
    1690:     MatView_Dense_Binary(A,viewer);
    1691:   } else if (isdraw) {
    1692:     MatView_SeqDense_Draw(A,viewer);
    1693:   }
    1694:   return(0);
    1695: }

    1697: static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A,const PetscScalar *array)
    1698: {
    1699:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;

    1702:   if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
    1703:   if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
    1704:   if (a->unplacedarray) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreArray() first");
    1705:   a->unplacedarray       = a->v;
    1706:   a->unplaced_user_alloc = a->user_alloc;
    1707:   a->v                   = (PetscScalar*) array;
    1708:   a->user_alloc          = PETSC_TRUE;
    1709: #if defined(PETSC_HAVE_CUDA)
    1710:   A->offloadmask = PETSC_OFFLOAD_CPU;
    1711: #endif
    1712:   return(0);
    1713: }

    1715: static PetscErrorCode MatDenseResetArray_SeqDense(Mat A)
    1716: {
    1717:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;

    1720:   if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
    1721:   if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
    1722:   a->v             = a->unplacedarray;
    1723:   a->user_alloc    = a->unplaced_user_alloc;
    1724:   a->unplacedarray = NULL;
    1725: #if defined(PETSC_HAVE_CUDA)
    1726:   A->offloadmask = PETSC_OFFLOAD_CPU;
    1727: #endif
    1728:   return(0);
    1729: }

    1731: static PetscErrorCode MatDenseReplaceArray_SeqDense(Mat A,const PetscScalar *array)
    1732: {
    1733:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

    1737:   if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
    1738:   if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
    1739:   if (!a->user_alloc) { PetscFree(a->v); }
    1740:   a->v           = (PetscScalar*) array;
    1741:   a->user_alloc  = PETSC_FALSE;
    1742: #if defined(PETSC_HAVE_CUDA)
    1743:   A->offloadmask = PETSC_OFFLOAD_CPU;
    1744: #endif
    1745:   return(0);
    1746: }

    1748: PetscErrorCode MatDestroy_SeqDense(Mat mat)
    1749: {
    1750:   Mat_SeqDense   *l = (Mat_SeqDense*)mat->data;

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

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

    1803:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumn_C",NULL);
    1804:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumn_C",NULL);
    1805:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVec_C",NULL);
    1806:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVec_C",NULL);
    1807:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVecRead_C",NULL);
    1808:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVecRead_C",NULL);
    1809:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVecWrite_C",NULL);
    1810:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVecWrite_C",NULL);
    1811:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetSubMatrix_C",NULL);
    1812:   PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreSubMatrix_C",NULL);
    1813:   return(0);
    1814: }

    1816: static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
    1817: {
    1818:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
    1820:   PetscInt       k,j,m = A->rmap->n, M = mat->lda, n = A->cmap->n;
    1821:   PetscScalar    *v,tmp;

    1824:   if (reuse == MAT_INPLACE_MATRIX) {
    1825:     if (m == n) { /* in place transpose */
    1826:       MatDenseGetArray(A,&v);
    1827:       for (j=0; j<m; j++) {
    1828:         for (k=0; k<j; k++) {
    1829:           tmp        = v[j + k*M];
    1830:           v[j + k*M] = v[k + j*M];
    1831:           v[k + j*M] = tmp;
    1832:         }
    1833:       }
    1834:       MatDenseRestoreArray(A,&v);
    1835:     } else { /* reuse memory, temporary allocates new memory */
    1836:       PetscScalar *v2;
    1837:       PetscLayout tmplayout;

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

    1865:     if (reuse == MAT_INITIAL_MATRIX) {
    1866:       MatCreate(PetscObjectComm((PetscObject)A),&tmat);
    1867:       MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);
    1868:       MatSetType(tmat,((PetscObject)A)->type_name);
    1869:       MatSeqDenseSetPreallocation(tmat,NULL);
    1870:     } else tmat = *matout;

    1872:     MatDenseGetArrayRead(A,(const PetscScalar**)&v);
    1873:     MatDenseGetArray(tmat,&v2);
    1874:     tmatd = (Mat_SeqDense*)tmat->data;
    1875:     M2    = tmatd->lda;
    1876:     for (j=0; j<n; j++) {
    1877:       for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M];
    1878:     }
    1879:     MatDenseRestoreArray(tmat,&v2);
    1880:     MatDenseRestoreArrayRead(A,(const PetscScalar**)&v);
    1881:     MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
    1882:     MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
    1883:     *matout = tmat;
    1884:   }
    1885:   return(0);
    1886: }

    1888: static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool  *flg)
    1889: {
    1890:   Mat_SeqDense      *mat1 = (Mat_SeqDense*)A1->data;
    1891:   Mat_SeqDense      *mat2 = (Mat_SeqDense*)A2->data;
    1892:   PetscInt          i;
    1893:   const PetscScalar *v1,*v2;
    1894:   PetscErrorCode    ierr;

    1897:   if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; return(0);}
    1898:   if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; return(0);}
    1899:   MatDenseGetArrayRead(A1,&v1);
    1900:   MatDenseGetArrayRead(A2,&v2);
    1901:   for (i=0; i<A1->cmap->n; i++) {
    1902:     PetscArraycmp(v1,v2,A1->rmap->n,flg);
    1903:     if (*flg == PETSC_FALSE) return(0);
    1904:     v1 += mat1->lda;
    1905:     v2 += mat2->lda;
    1906:   }
    1907:   MatDenseRestoreArrayRead(A1,&v1);
    1908:   MatDenseRestoreArrayRead(A2,&v2);
    1909:   *flg = PETSC_TRUE;
    1910:   return(0);
    1911: }

    1913: static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
    1914: {
    1915:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
    1916:   PetscInt          i,n,len;
    1917:   PetscScalar       *x;
    1918:   const PetscScalar *vv;
    1919:   PetscErrorCode    ierr;

    1922:   VecGetSize(v,&n);
    1923:   VecGetArray(v,&x);
    1924:   len  = PetscMin(A->rmap->n,A->cmap->n);
    1925:   MatDenseGetArrayRead(A,&vv);
    1926:   if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
    1927:   for (i=0; i<len; i++) {
    1928:     x[i] = vv[i*mat->lda + i];
    1929:   }
    1930:   MatDenseRestoreArrayRead(A,&vv);
    1931:   VecRestoreArray(v,&x);
    1932:   return(0);
    1933: }

    1935: static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
    1936: {
    1937:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
    1938:   const PetscScalar *l,*r;
    1939:   PetscScalar       x,*v,*vv;
    1940:   PetscErrorCode    ierr;
    1941:   PetscInt          i,j,m = A->rmap->n,n = A->cmap->n;

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

    1973: PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
    1974: {
    1975:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
    1976:   PetscScalar       *v,*vv;
    1977:   PetscReal         sum  = 0.0;
    1978:   PetscInt          lda  =mat->lda,m=A->rmap->n,i,j;
    1979:   PetscErrorCode    ierr;

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

    2032: static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
    2033: {
    2034:   Mat_SeqDense   *aij = (Mat_SeqDense*)A->data;

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

    2067: static PetscErrorCode MatZeroEntries_SeqDense(Mat A)
    2068: {
    2069:   Mat_SeqDense   *l = (Mat_SeqDense*)A->data;
    2071:   PetscInt       lda=l->lda,m=A->rmap->n,j;
    2072:   PetscScalar    *v;

    2075:   MatDenseGetArray(A,&v);
    2076:   if (lda>m) {
    2077:     for (j=0; j<A->cmap->n; j++) {
    2078:       PetscArrayzero(v+j*lda,m);
    2079:     }
    2080:   } else {
    2081:     PetscArrayzero(v,A->rmap->n*A->cmap->n);
    2082:   }
    2083:   MatDenseRestoreArray(A,&v);
    2084:   return(0);
    2085: }

    2087: static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
    2088: {
    2089:   PetscErrorCode    ierr;
    2090:   Mat_SeqDense      *l = (Mat_SeqDense*)A->data;
    2091:   PetscInt          m  = l->lda, n = A->cmap->n, i,j;
    2092:   PetscScalar       *slot,*bb,*v;
    2093:   const PetscScalar *xx;

    2096:   if (PetscDefined(USE_DEBUG)) {
    2097:     for (i=0; i<N; i++) {
    2098:       if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
    2099:       if (rows[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D requested to be zeroed greater than or equal number of rows %D",rows[i],A->rmap->n);
    2100:     }
    2101:   }
    2102:   if (!N) return(0);

    2104:   /* fix right hand side if needed */
    2105:   if (x && b) {
    2106:     VecGetArrayRead(x,&xx);
    2107:     VecGetArray(b,&bb);
    2108:     for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
    2109:     VecRestoreArrayRead(x,&xx);
    2110:     VecRestoreArray(b,&bb);
    2111:   }

    2113:   MatDenseGetArray(A,&v);
    2114:   for (i=0; i<N; i++) {
    2115:     slot = v + rows[i];
    2116:     for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
    2117:   }
    2118:   if (diag != 0.0) {
    2119:     if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
    2120:     for (i=0; i<N; i++) {
    2121:       slot  = v + (m+1)*rows[i];
    2122:       *slot = diag;
    2123:     }
    2124:   }
    2125:   MatDenseRestoreArray(A,&v);
    2126:   return(0);
    2127: }

    2129: static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A,PetscInt *lda)
    2130: {
    2131:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;

    2134:   *lda = mat->lda;
    2135:   return(0);
    2136: }

    2138: PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar **array)
    2139: {
    2140:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;

    2143:   if (mat->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
    2144:   *array = mat->v;
    2145:   return(0);
    2146: }

    2148: PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar **array)
    2149: {
    2151:   *array = NULL;
    2152:   return(0);
    2153: }

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

    2158:    Not collective

    2160:    Input Parameter:
    2161: .  mat - a MATSEQDENSE or MATMPIDENSE matrix

    2163:    Output Parameter:
    2164: .   lda - the leading dimension

    2166:    Level: intermediate

    2168: .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseSetLDA()
    2169: @*/
    2170: PetscErrorCode  MatDenseGetLDA(Mat A,PetscInt *lda)
    2171: {

    2177:   PetscUseMethod(A,"MatDenseGetLDA_C",(Mat,PetscInt*),(A,lda));
    2178:   return(0);
    2179: }

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

    2184:    Not collective

    2186:    Input Parameter:
    2187: +  mat - a MATSEQDENSE or MATMPIDENSE matrix
    2188: -  lda - the leading dimension

    2190:    Level: intermediate

    2192: .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetLDA()
    2193: @*/
    2194: PetscErrorCode  MatDenseSetLDA(Mat A,PetscInt lda)
    2195: {

    2200:   PetscTryMethod(A,"MatDenseSetLDA_C",(Mat,PetscInt),(A,lda));
    2201:   return(0);
    2202: }

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

    2207:    Logically Collective on Mat

    2209:    Input Parameter:
    2210: .  mat - a dense matrix

    2212:    Output Parameter:
    2213: .   array - pointer to the data

    2215:    Level: intermediate

    2217: .seealso: MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
    2218: @*/
    2219: PetscErrorCode  MatDenseGetArray(Mat A,PetscScalar **array)
    2220: {

    2226:   PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));
    2227:   return(0);
    2228: }

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

    2233:    Logically Collective on Mat

    2235:    Input Parameters:
    2236: +  mat - a dense matrix
    2237: -  array - pointer to the data

    2239:    Level: intermediate

    2241: .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
    2242: @*/
    2243: PetscErrorCode  MatDenseRestoreArray(Mat A,PetscScalar **array)
    2244: {

    2250:   PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));
    2251:   PetscObjectStateIncrease((PetscObject)A);
    2252: #if defined(PETSC_HAVE_CUDA)
    2253:   A->offloadmask = PETSC_OFFLOAD_CPU;
    2254: #endif
    2255:   return(0);
    2256: }

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

    2261:    Not Collective

    2263:    Input Parameter:
    2264: .  mat - a dense matrix

    2266:    Output Parameter:
    2267: .   array - pointer to the data

    2269:    Level: intermediate

    2271: .seealso: MatDenseRestoreArrayRead(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
    2272: @*/
    2273: PetscErrorCode  MatDenseGetArrayRead(Mat A,const PetscScalar **array)
    2274: {

    2280:   PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array));
    2281:   return(0);
    2282: }

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

    2287:    Not Collective

    2289:    Input Parameters:
    2290: +  mat - a dense matrix
    2291: -  array - pointer to the data

    2293:    Level: intermediate

    2295: .seealso: MatDenseGetArrayRead(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
    2296: @*/
    2297: PetscErrorCode  MatDenseRestoreArrayRead(Mat A,const PetscScalar **array)
    2298: {

    2304:   PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array));
    2305:   return(0);
    2306: }

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

    2311:    Not Collective

    2313:    Input Parameter:
    2314: .  mat - a dense matrix

    2316:    Output Parameter:
    2317: .   array - pointer to the data

    2319:    Level: intermediate

    2321: .seealso: MatDenseRestoreArrayWrite(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead()
    2322: @*/
    2323: PetscErrorCode  MatDenseGetArrayWrite(Mat A,PetscScalar **array)
    2324: {

    2330:   PetscUseMethod(A,"MatDenseGetArrayWrite_C",(Mat,PetscScalar**),(A,array));
    2331:   return(0);
    2332: }

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

    2337:    Not Collective

    2339:    Input Parameters:
    2340: +  mat - a dense matrix
    2341: -  array - pointer to the data

    2343:    Level: intermediate

    2345: .seealso: MatDenseGetArrayWrite(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead()
    2346: @*/
    2347: PetscErrorCode  MatDenseRestoreArrayWrite(Mat A,PetscScalar **array)
    2348: {

    2354:   PetscUseMethod(A,"MatDenseRestoreArrayWrite_C",(Mat,PetscScalar**),(A,array));
    2355:   PetscObjectStateIncrease((PetscObject)A);
    2356: #if defined(PETSC_HAVE_CUDA)
    2357:   A->offloadmask = PETSC_OFFLOAD_CPU;
    2358: #endif
    2359:   return(0);
    2360: }

    2362: static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,MatReuse scall,Mat *B)
    2363: {
    2364:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
    2366:   PetscInt       i,j,nrows,ncols,blda;
    2367:   const PetscInt *irow,*icol;
    2368:   PetscScalar    *av,*bv,*v = mat->v;
    2369:   Mat            newmat;

    2372:   ISGetIndices(isrow,&irow);
    2373:   ISGetIndices(iscol,&icol);
    2374:   ISGetLocalSize(isrow,&nrows);
    2375:   ISGetLocalSize(iscol,&ncols);

    2377:   /* Check submatrixcall */
    2378:   if (scall == MAT_REUSE_MATRIX) {
    2379:     PetscInt n_cols,n_rows;
    2380:     MatGetSize(*B,&n_rows,&n_cols);
    2381:     if (n_rows != nrows || n_cols != ncols) {
    2382:       /* resize the result matrix to match number of requested rows/columns */
    2383:       MatSetSizes(*B,nrows,ncols,nrows,ncols);
    2384:     }
    2385:     newmat = *B;
    2386:   } else {
    2387:     /* Create and fill new matrix */
    2388:     MatCreate(PetscObjectComm((PetscObject)A),&newmat);
    2389:     MatSetSizes(newmat,nrows,ncols,nrows,ncols);
    2390:     MatSetType(newmat,((PetscObject)A)->type_name);
    2391:     MatSeqDenseSetPreallocation(newmat,NULL);
    2392:   }

    2394:   /* Now extract the data pointers and do the copy,column at a time */
    2395:   MatDenseGetArray(newmat,&bv);
    2396:   MatDenseGetLDA(newmat,&blda);
    2397:   for (i=0; i<ncols; i++) {
    2398:     av = v + mat->lda*icol[i];
    2399:     for (j=0; j<nrows; j++) bv[j] = av[irow[j]];
    2400:     bv += blda;
    2401:   }
    2402:   MatDenseRestoreArray(newmat,&bv);

    2404:   /* Assemble the matrices so that the correct flags are set */
    2405:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
    2406:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);

    2408:   /* Free work space */
    2409:   ISRestoreIndices(isrow,&irow);
    2410:   ISRestoreIndices(iscol,&icol);
    2411:   *B   = newmat;
    2412:   return(0);
    2413: }

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

    2421:   if (scall == MAT_INITIAL_MATRIX) {
    2422:     PetscCalloc1(n,B);
    2423:   }

    2425:   for (i=0; i<n; i++) {
    2426:     MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],scall,&(*B)[i]);
    2427:   }
    2428:   return(0);
    2429: }

    2431: static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
    2432: {
    2434:   return(0);
    2435: }

    2437: static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
    2438: {
    2440:   return(0);
    2441: }

    2443: PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
    2444: {
    2445:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data;
    2446:   PetscErrorCode    ierr;
    2447:   const PetscScalar *va;
    2448:   PetscScalar       *vb;
    2449:   PetscInt          lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j;

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

    2474: static PetscErrorCode MatSetUp_SeqDense(Mat A)
    2475: {

    2479:   PetscLayoutSetUp(A->rmap);
    2480:   PetscLayoutSetUp(A->cmap);
    2481:   if (!A->preallocated) {
    2482:     MatSeqDenseSetPreallocation(A,NULL);
    2483:   }
    2484:   return(0);
    2485: }

    2487: static PetscErrorCode MatConjugate_SeqDense(Mat A)
    2488: {
    2489:   Mat_SeqDense   *mat = (Mat_SeqDense *) A->data;
    2490:   PetscInt       i,nz = A->rmap->n*A->cmap->n;
    2491:   PetscInt       min = PetscMin(A->rmap->n,A->cmap->n);
    2492:   PetscScalar    *aa;

    2496:   MatDenseGetArray(A,&aa);
    2497:   for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
    2498:   MatDenseRestoreArray(A,&aa);
    2499:   if (mat->tau) for (i = 0; i < min; i++) mat->tau[i] = PetscConj(mat->tau[i]);
    2500:   return(0);
    2501: }

    2503: static PetscErrorCode MatRealPart_SeqDense(Mat A)
    2504: {
    2505:   PetscInt       i,nz = A->rmap->n*A->cmap->n;
    2506:   PetscScalar    *aa;

    2510:   MatDenseGetArray(A,&aa);
    2511:   for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
    2512:   MatDenseRestoreArray(A,&aa);
    2513:   return(0);
    2514: }

    2516: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
    2517: {
    2518:   PetscInt       i,nz = A->rmap->n*A->cmap->n;
    2519:   PetscScalar    *aa;

    2523:   MatDenseGetArray(A,&aa);
    2524:   for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
    2525:   MatDenseRestoreArray(A,&aa);
    2526:   return(0);
    2527: }

    2529: /* ----------------------------------------------------------------*/
    2530: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
    2531: {
    2533:   PetscInt       m=A->rmap->n,n=B->cmap->n;
    2534:   PetscBool      cisdense;

    2537:   MatSetSizes(C,m,n,m,n);
    2538:   PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
    2539:   if (!cisdense) {
    2540:     PetscBool flg;

    2542:     PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
    2543:     MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
    2544:   }
    2545:   MatSetUp(C);
    2546:   return(0);
    2547: }

    2549: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
    2550: {
    2551:   Mat_SeqDense       *a=(Mat_SeqDense*)A->data,*b=(Mat_SeqDense*)B->data,*c=(Mat_SeqDense*)C->data;
    2552:   PetscBLASInt       m,n,k;
    2553:   const PetscScalar *av,*bv;
    2554:   PetscScalar       *cv;
    2555:   PetscScalar       _DOne=1.0,_DZero=0.0;
    2556:   PetscErrorCode    ierr;

    2559:   PetscBLASIntCast(C->rmap->n,&m);
    2560:   PetscBLASIntCast(C->cmap->n,&n);
    2561:   PetscBLASIntCast(A->cmap->n,&k);
    2562:   if (!m || !n || !k) return(0);
    2563:   MatDenseGetArrayRead(A,&av);
    2564:   MatDenseGetArrayRead(B,&bv);
    2565:   MatDenseGetArrayWrite(C,&cv);
    2566:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
    2567:   PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
    2568:   MatDenseRestoreArrayRead(A,&av);
    2569:   MatDenseRestoreArrayRead(B,&bv);
    2570:   MatDenseRestoreArrayWrite(C,&cv);
    2571:   return(0);
    2572: }

    2574: PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
    2575: {
    2577:   PetscInt       m=A->rmap->n,n=B->rmap->n;
    2578:   PetscBool      cisdense;

    2581:   MatSetSizes(C,m,n,m,n);
    2582:   PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
    2583:   if (!cisdense) {
    2584:     PetscBool flg;

    2586:     PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
    2587:     MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
    2588:   }
    2589:   MatSetUp(C);
    2590:   return(0);
    2591: }

    2593: PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
    2594: {
    2595:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
    2596:   Mat_SeqDense      *b = (Mat_SeqDense*)B->data;
    2597:   Mat_SeqDense      *c = (Mat_SeqDense*)C->data;
    2598:   const PetscScalar *av,*bv;
    2599:   PetscScalar       *cv;
    2600:   PetscBLASInt      m,n,k;
    2601:   PetscScalar       _DOne=1.0,_DZero=0.0;
    2602:   PetscErrorCode    ierr;

    2605:   PetscBLASIntCast(C->rmap->n,&m);
    2606:   PetscBLASIntCast(C->cmap->n,&n);
    2607:   PetscBLASIntCast(A->cmap->n,&k);
    2608:   if (!m || !n || !k) return(0);
    2609:   MatDenseGetArrayRead(A,&av);
    2610:   MatDenseGetArrayRead(B,&bv);
    2611:   MatDenseGetArrayWrite(C,&cv);
    2612:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
    2613:   MatDenseRestoreArrayRead(A,&av);
    2614:   MatDenseRestoreArrayRead(B,&bv);
    2615:   MatDenseRestoreArrayWrite(C,&cv);
    2616:   PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
    2617:   return(0);
    2618: }

    2620: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
    2621: {
    2623:   PetscInt       m=A->cmap->n,n=B->cmap->n;
    2624:   PetscBool      cisdense;

    2627:   MatSetSizes(C,m,n,m,n);
    2628:   PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
    2629:   if (!cisdense) {
    2630:     PetscBool flg;

    2632:     PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
    2633:     MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
    2634:   }
    2635:   MatSetUp(C);
    2636:   return(0);
    2637: }

    2639: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
    2640: {
    2641:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
    2642:   Mat_SeqDense      *b = (Mat_SeqDense*)B->data;
    2643:   Mat_SeqDense      *c = (Mat_SeqDense*)C->data;
    2644:   const PetscScalar *av,*bv;
    2645:   PetscScalar       *cv;
    2646:   PetscBLASInt      m,n,k;
    2647:   PetscScalar       _DOne=1.0,_DZero=0.0;
    2648:   PetscErrorCode    ierr;

    2651:   PetscBLASIntCast(C->rmap->n,&m);
    2652:   PetscBLASIntCast(C->cmap->n,&n);
    2653:   PetscBLASIntCast(A->rmap->n,&k);
    2654:   if (!m || !n || !k) return(0);
    2655:   MatDenseGetArrayRead(A,&av);
    2656:   MatDenseGetArrayRead(B,&bv);
    2657:   MatDenseGetArrayWrite(C,&cv);
    2658:   PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
    2659:   MatDenseRestoreArrayRead(A,&av);
    2660:   MatDenseRestoreArrayRead(B,&bv);
    2661:   MatDenseRestoreArrayWrite(C,&cv);
    2662:   PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
    2663:   return(0);
    2664: }

    2666: /* ----------------------------------------------- */
    2667: static PetscErrorCode MatProductSetFromOptions_SeqDense_AB(Mat C)
    2668: {
    2670:   C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqDense;
    2671:   C->ops->productsymbolic = MatProductSymbolic_AB;
    2672:   return(0);
    2673: }

    2675: static PetscErrorCode MatProductSetFromOptions_SeqDense_AtB(Mat C)
    2676: {
    2678:   C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_SeqDense_SeqDense;
    2679:   C->ops->productsymbolic          = MatProductSymbolic_AtB;
    2680:   return(0);
    2681: }

    2683: static PetscErrorCode MatProductSetFromOptions_SeqDense_ABt(Mat C)
    2684: {
    2686:   C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqDense_SeqDense;
    2687:   C->ops->productsymbolic          = MatProductSymbolic_ABt;
    2688:   return(0);
    2689: }

    2691: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense(Mat C)
    2692: {
    2694:   Mat_Product    *product = C->product;

    2697:   switch (product->type) {
    2698:   case MATPRODUCT_AB:
    2699:     MatProductSetFromOptions_SeqDense_AB(C);
    2700:     break;
    2701:   case MATPRODUCT_AtB:
    2702:     MatProductSetFromOptions_SeqDense_AtB(C);
    2703:     break;
    2704:   case MATPRODUCT_ABt:
    2705:     MatProductSetFromOptions_SeqDense_ABt(C);
    2706:     break;
    2707:   default:
    2708:     break;
    2709:   }
    2710:   return(0);
    2711: }
    2712: /* ----------------------------------------------- */

    2714: static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[])
    2715: {
    2716:   Mat_SeqDense       *a = (Mat_SeqDense*)A->data;
    2717:   PetscErrorCode     ierr;
    2718:   PetscInt           i,j,m = A->rmap->n,n = A->cmap->n,p;
    2719:   PetscScalar        *x;
    2720:   const PetscScalar *aa;

    2723:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    2724:   VecGetArray(v,&x);
    2725:   VecGetLocalSize(v,&p);
    2726:   MatDenseGetArrayRead(A,&aa);
    2727:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
    2728:   for (i=0; i<m; i++) {
    2729:     x[i] = aa[i]; if (idx) idx[i] = 0;
    2730:     for (j=1; j<n; j++) {
    2731:       if (PetscRealPart(x[i]) < PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;}
    2732:     }
    2733:   }
    2734:   MatDenseRestoreArrayRead(A,&aa);
    2735:   VecRestoreArray(v,&x);
    2736:   return(0);
    2737: }

    2739: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])
    2740: {
    2741:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
    2742:   PetscErrorCode    ierr;
    2743:   PetscInt          i,j,m = A->rmap->n,n = A->cmap->n,p;
    2744:   PetscScalar       *x;
    2745:   PetscReal         atmp;
    2746:   const PetscScalar *aa;

    2749:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    2750:   VecGetArray(v,&x);
    2751:   VecGetLocalSize(v,&p);
    2752:   MatDenseGetArrayRead(A,&aa);
    2753:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
    2754:   for (i=0; i<m; i++) {
    2755:     x[i] = PetscAbsScalar(aa[i]);
    2756:     for (j=1; j<n; j++) {
    2757:       atmp = PetscAbsScalar(aa[i+a->lda*j]);
    2758:       if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;}
    2759:     }
    2760:   }
    2761:   MatDenseRestoreArrayRead(A,&aa);
    2762:   VecRestoreArray(v,&x);
    2763:   return(0);
    2764: }

    2766: static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])
    2767: {
    2768:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
    2769:   PetscErrorCode    ierr;
    2770:   PetscInt          i,j,m = A->rmap->n,n = A->cmap->n,p;
    2771:   PetscScalar       *x;
    2772:   const PetscScalar *aa;

    2775:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    2776:   MatDenseGetArrayRead(A,&aa);
    2777:   VecGetArray(v,&x);
    2778:   VecGetLocalSize(v,&p);
    2779:   if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
    2780:   for (i=0; i<m; i++) {
    2781:     x[i] = aa[i]; if (idx) idx[i] = 0;
    2782:     for (j=1; j<n; j++) {
    2783:       if (PetscRealPart(x[i]) > PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;}
    2784:     }
    2785:   }
    2786:   VecRestoreArray(v,&x);
    2787:   MatDenseRestoreArrayRead(A,&aa);
    2788:   return(0);
    2789: }

    2791: PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
    2792: {
    2793:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
    2794:   PetscErrorCode    ierr;
    2795:   PetscScalar       *x;
    2796:   const PetscScalar *aa;

    2799:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    2800:   MatDenseGetArrayRead(A,&aa);
    2801:   VecGetArray(v,&x);
    2802:   PetscArraycpy(x,aa+col*a->lda,A->rmap->n);
    2803:   VecRestoreArray(v,&x);
    2804:   MatDenseRestoreArrayRead(A,&aa);
    2805:   return(0);
    2806: }

    2808: PETSC_INTERN PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms)
    2809: {
    2810:   PetscErrorCode    ierr;
    2811:   PetscInt          i,j,m,n;
    2812:   const PetscScalar *a;

    2815:   MatGetSize(A,&m,&n);
    2816:   PetscArrayzero(norms,n);
    2817:   MatDenseGetArrayRead(A,&a);
    2818:   if (type == NORM_2) {
    2819:     for (i=0; i<n; i++) {
    2820:       for (j=0; j<m; j++) {
    2821:         norms[i] += PetscAbsScalar(a[j]*a[j]);
    2822:       }
    2823:       a += m;
    2824:     }
    2825:   } else if (type == NORM_1) {
    2826:     for (i=0; i<n; i++) {
    2827:       for (j=0; j<m; j++) {
    2828:         norms[i] += PetscAbsScalar(a[j]);
    2829:       }
    2830:       a += m;
    2831:     }
    2832:   } else if (type == NORM_INFINITY) {
    2833:     for (i=0; i<n; i++) {
    2834:       for (j=0; j<m; j++) {
    2835:         norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]);
    2836:       }
    2837:       a += m;
    2838:     }
    2839:   } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType");
    2840:   MatDenseRestoreArrayRead(A,&a);
    2841:   if (type == NORM_2) {
    2842:     for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]);
    2843:   }
    2844:   return(0);
    2845: }

    2847: static PetscErrorCode  MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
    2848: {
    2850:   PetscScalar    *a;
    2851:   PetscInt       lda,m,n,i,j;

    2854:   MatGetSize(x,&m,&n);
    2855:   MatDenseGetLDA(x,&lda);
    2856:   MatDenseGetArray(x,&a);
    2857:   for (j=0; j<n; j++) {
    2858:     for (i=0; i<m; i++) {
    2859:       PetscRandomGetValue(rctx,a+j*lda+i);
    2860:     }
    2861:   }
    2862:   MatDenseRestoreArray(x,&a);
    2863:   return(0);
    2864: }

    2866: static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool  *missing,PetscInt *d)
    2867: {
    2869:   *missing = PETSC_FALSE;
    2870:   return(0);
    2871: }

    2873: /* vals is not const */
    2874: static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals)
    2875: {
    2877:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
    2878:   PetscScalar    *v;

    2881:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    2882:   MatDenseGetArray(A,&v);
    2883:   *vals = v+col*a->lda;
    2884:   MatDenseRestoreArray(A,&v);
    2885:   return(0);
    2886: }

    2888: static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals)
    2889: {
    2891:   *vals = NULL; /* user cannot accidently use the array later */
    2892:   return(0);
    2893: }

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

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

    3051:    Collective

    3053:    Input Parameters:
    3054: +  comm - MPI communicator, set to PETSC_COMM_SELF
    3055: .  m - number of rows
    3056: .  n - number of columns
    3057: -  data - optional location of matrix data in column major order.  Set data=NULL for PETSc
    3058:    to control all matrix memory allocation.

    3060:    Output Parameter:
    3061: .  A - the matrix

    3063:    Notes:
    3064:    The data input variable is intended primarily for Fortran programmers
    3065:    who wish to allocate their own matrix memory space.  Most users should
    3066:    set data=NULL.

    3068:    Level: intermediate

    3070: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
    3071: @*/
    3072: PetscErrorCode  MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
    3073: {

    3077:   MatCreate(comm,A);
    3078:   MatSetSizes(*A,m,n,m,n);
    3079:   MatSetType(*A,MATSEQDENSE);
    3080:   MatSeqDenseSetPreallocation(*A,data);
    3081:   return(0);
    3082: }

    3084: /*@C
    3085:    MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements

    3087:    Collective

    3089:    Input Parameters:
    3090: +  B - the matrix
    3091: -  data - the array (or NULL)

    3093:    Notes:
    3094:    The data input variable is intended primarily for Fortran programmers
    3095:    who wish to allocate their own matrix memory space.  Most users should
    3096:    need not call this routine.

    3098:    Level: intermediate

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

    3102: @*/
    3103: PetscErrorCode  MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
    3104: {

    3109:   PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
    3110:   return(0);
    3111: }

    3113: PetscErrorCode  MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
    3114: {
    3115:   Mat_SeqDense   *b = (Mat_SeqDense*)B->data;

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

    3122:   PetscLayoutSetUp(B->rmap);
    3123:   PetscLayoutSetUp(B->cmap);

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

    3127:   PetscIntMultError(b->lda,B->cmap->n,NULL);
    3128:   if (!data) { /* petsc-allocated storage */
    3129:     if (!b->user_alloc) { PetscFree(b->v); }
    3130:     PetscCalloc1((size_t)b->lda*B->cmap->n,&b->v);
    3131:     PetscLogObjectMemory((PetscObject)B,b->lda*B->cmap->n*sizeof(PetscScalar));

    3133:     b->user_alloc = PETSC_FALSE;
    3134:   } else { /* user-allocated storage */
    3135:     if (!b->user_alloc) { PetscFree(b->v); }
    3136:     b->v          = data;
    3137:     b->user_alloc = PETSC_TRUE;
    3138:   }
    3139:   B->assembled = PETSC_TRUE;
    3140:   return(0);
    3141: }

    3143: #if defined(PETSC_HAVE_ELEMENTAL)
    3144: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
    3145: {
    3146:   Mat               mat_elemental;
    3147:   PetscErrorCode    ierr;
    3148:   const PetscScalar *array;
    3149:   PetscScalar       *v_colwise;
    3150:   PetscInt          M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols;

    3153:   PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);
    3154:   MatDenseGetArrayRead(A,&array);
    3155:   /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
    3156:   k = 0;
    3157:   for (j=0; j<N; j++) {
    3158:     cols[j] = j;
    3159:     for (i=0; i<M; i++) {
    3160:       v_colwise[j*M+i] = array[k++];
    3161:     }
    3162:   }
    3163:   for (i=0; i<M; i++) {
    3164:     rows[i] = i;
    3165:   }
    3166:   MatDenseRestoreArrayRead(A,&array);

    3168:   MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
    3169:   MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);
    3170:   MatSetType(mat_elemental,MATELEMENTAL);
    3171:   MatSetUp(mat_elemental);

    3173:   /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
    3174:   MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);
    3175:   MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
    3176:   MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);
    3177:   PetscFree3(v_colwise,rows,cols);

    3179:   if (reuse == MAT_INPLACE_MATRIX) {
    3180:     MatHeaderReplace(A,&mat_elemental);
    3181:   } else {
    3182:     *newmat = mat_elemental;
    3183:   }
    3184:   return(0);
    3185: }
    3186: #endif

    3188: PetscErrorCode  MatDenseSetLDA_SeqDense(Mat B,PetscInt lda)
    3189: {
    3190:   Mat_SeqDense *b = (Mat_SeqDense*)B->data;
    3191:   PetscBool    data;

    3194:   data = (PetscBool)((B->rmap->n > 0 && B->cmap->n > 0) ? (b->v ? PETSC_TRUE : PETSC_FALSE) : PETSC_FALSE);
    3195:   if (!b->user_alloc && data && b->lda!=lda) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"LDA cannot be changed after allocation of internal storage");
    3196:   if (lda < B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"LDA %D must be at least matrix dimension %D",lda,B->rmap->n);
    3197:   b->lda = lda;
    3198:   return(0);
    3199: }

    3201: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat)
    3202: {
    3204:   PetscMPIInt    size;

    3207:   MPI_Comm_size(comm,&size);
    3208:   if (size == 1) {
    3209:     if (scall == MAT_INITIAL_MATRIX) {
    3210:       MatDuplicate(inmat,MAT_COPY_VALUES,outmat);
    3211:     } else {
    3212:       MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);
    3213:     }
    3214:   } else {
    3215:     MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);
    3216:   }
    3217:   return(0);
    3218: }

    3220: PetscErrorCode MatDenseGetColumnVec_SeqDense(Mat A,PetscInt col,Vec *v)
    3221: {
    3222:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

    3226:   if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
    3227:   if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
    3228:   if (!a->cvec) {
    3229:     VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);
    3230:     PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);
    3231:   }
    3232:   a->vecinuse = col + 1;
    3233:   MatDenseGetArray(A,(PetscScalar**)&a->ptrinuse);
    3234:   VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);
    3235:   *v   = a->cvec;
    3236:   return(0);
    3237: }

    3239: PetscErrorCode MatDenseRestoreColumnVec_SeqDense(Mat A,PetscInt col,Vec *v)
    3240: {
    3241:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

    3245:   if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first");
    3246:   if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector");
    3247:   a->vecinuse = 0;
    3248:   MatDenseRestoreArray(A,(PetscScalar**)&a->ptrinuse);
    3249:   VecResetArray(a->cvec);
    3250:   *v   = NULL;
    3251:   return(0);
    3252: }

    3254: PetscErrorCode MatDenseGetColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v)
    3255: {
    3256:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

    3260:   if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
    3261:   if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
    3262:   if (!a->cvec) {
    3263:     VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);
    3264:     PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);
    3265:   }
    3266:   a->vecinuse = col + 1;
    3267:   MatDenseGetArrayRead(A,&a->ptrinuse);
    3268:   VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);
    3269:   VecLockReadPush(a->cvec);
    3270:   *v   = a->cvec;
    3271:   return(0);
    3272: }

    3274: PetscErrorCode MatDenseRestoreColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v)
    3275: {
    3276:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

    3280:   if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first");
    3281:   if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector");
    3282:   a->vecinuse = 0;
    3283:   MatDenseRestoreArrayRead(A,&a->ptrinuse);
    3284:   VecLockReadPop(a->cvec);
    3285:   VecResetArray(a->cvec);
    3286:   *v   = NULL;
    3287:   return(0);
    3288: }

    3290: PetscErrorCode MatDenseGetColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v)
    3291: {
    3292:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

    3296:   if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
    3297:   if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
    3298:   if (!a->cvec) {
    3299:     VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);
    3300:     PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);
    3301:   }
    3302:   a->vecinuse = col + 1;
    3303:   MatDenseGetArrayWrite(A,(PetscScalar**)&a->ptrinuse);
    3304:   VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);
    3305:   *v   = a->cvec;
    3306:   return(0);
    3307: }

    3309: PetscErrorCode MatDenseRestoreColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v)
    3310: {
    3311:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

    3315:   if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first");
    3316:   if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector");
    3317:   a->vecinuse = 0;
    3318:   MatDenseRestoreArrayWrite(A,(PetscScalar**)&a->ptrinuse);
    3319:   VecResetArray(a->cvec);
    3320:   *v   = NULL;
    3321:   return(0);
    3322: }

    3324: PetscErrorCode MatDenseGetSubMatrix_SeqDense(Mat A,PetscInt cbegin,PetscInt cend,Mat *v)
    3325: {
    3326:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

    3330:   if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
    3331:   if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
    3332:   if (a->cmat && cend-cbegin != a->cmat->cmap->N) {
    3333:     MatDestroy(&a->cmat);
    3334:   }
    3335:   if (!a->cmat) {
    3336:     MatCreateDense(PetscObjectComm((PetscObject)A),A->rmap->n,PETSC_DECIDE,A->rmap->N,cend-cbegin,(PetscScalar*)a->v + (size_t)cbegin * (size_t)a->lda,&a->cmat);
    3337:     PetscLogObjectParent((PetscObject)A,(PetscObject)a->cmat);
    3338:   } else {
    3339:     MatDensePlaceArray(a->cmat,a->v + (size_t)cbegin * (size_t)a->lda);
    3340:   }
    3341:   MatDenseSetLDA(a->cmat,a->lda);
    3342:   a->matinuse = cbegin + 1;
    3343:   *v = a->cmat;
    3344:   return(0);
    3345: }

    3347: PetscErrorCode MatDenseRestoreSubMatrix_SeqDense(Mat A,Mat *v)
    3348: {
    3349:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

    3353:   if (!a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetSubMatrix() first");
    3354:   if (!a->cmat) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column matrix");
    3355:   if (*v != a->cmat) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not the matrix obtained from MatDenseGetSubMatrix()");
    3356:   a->matinuse = 0;
    3357:   MatDenseResetArray(a->cmat);
    3358:   *v   = NULL;
    3359:   return(0);
    3360: }

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

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

    3368:   Level: beginner

    3370: .seealso: MatCreateSeqDense()

    3372: M*/
    3373: PetscErrorCode MatCreate_SeqDense(Mat B)
    3374: {
    3375:   Mat_SeqDense   *b;
    3377:   PetscMPIInt    size;

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

    3383:   PetscNewLog(B,&b);
    3384:   PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
    3385:   B->data = (void*)b;

    3387:   b->roworiented = PETSC_TRUE;

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

    3422:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);
    3423:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);
    3424:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVec_C",MatDenseGetColumnVec_SeqDense);
    3425:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVec_C",MatDenseRestoreColumnVec_SeqDense);
    3426:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecRead_C",MatDenseGetColumnVecRead_SeqDense);
    3427:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecRead_C",MatDenseRestoreColumnVecRead_SeqDense);
    3428:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecWrite_C",MatDenseGetColumnVecWrite_SeqDense);
    3429:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecWrite_C",MatDenseRestoreColumnVecWrite_SeqDense);
    3430:   PetscObjectComposeFunction((PetscObject)B,"MatDenseGetSubMatrix_C",MatDenseGetSubMatrix_SeqDense);
    3431:   PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreSubMatrix_C",MatDenseRestoreSubMatrix_SeqDense);
    3432:   PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
    3433:   return(0);
    3434: }

    3436: /*@C
    3437:    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.

    3439:    Not Collective

    3441:    Input Parameters:
    3442: +  mat - a MATSEQDENSE or MATMPIDENSE matrix
    3443: -  col - column index

    3445:    Output Parameter:
    3446: .  vals - pointer to the data

    3448:    Level: intermediate

    3450: .seealso: MatDenseRestoreColumn()
    3451: @*/
    3452: PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals)
    3453: {

    3460:   PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));
    3461:   return(0);
    3462: }

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

    3467:    Not Collective

    3469:    Input Parameter:
    3470: .  mat - a MATSEQDENSE or MATMPIDENSE matrix

    3472:    Output Parameter:
    3473: .  vals - pointer to the data

    3475:    Level: intermediate

    3477: .seealso: MatDenseGetColumn()
    3478: @*/
    3479: PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals)
    3480: {

    3486:   PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));
    3487:   return(0);
    3488: }

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

    3493:    Collective

    3495:    Input Parameters:
    3496: +  mat - the Mat object
    3497: -  col - the column index

    3499:    Output Parameter:
    3500: .  v - the vector

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

    3506:    Level: intermediate

    3508: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
    3509: @*/
    3510: PetscErrorCode MatDenseGetColumnVec(Mat A,PetscInt col,Vec *v)
    3511: {

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

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

    3528:    Collective

    3530:    Input Parameters:
    3531: +  mat - the Mat object
    3532: .  col - the column index
    3533: -  v - the Vec object

    3535:    Level: intermediate

    3537: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
    3538: @*/
    3539: PetscErrorCode MatDenseRestoreColumnVec(Mat A,PetscInt col,Vec *v)
    3540: {

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

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

    3557:    Collective

    3559:    Input Parameters:
    3560: +  mat - the Mat object
    3561: -  col - the column index

    3563:    Output Parameter:
    3564: .  v - the vector

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

    3571:    Level: intermediate

    3573: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
    3574: @*/
    3575: PetscErrorCode MatDenseGetColumnVecRead(Mat A,PetscInt col,Vec *v)
    3576: {

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

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

    3593:    Collective

    3595:    Input Parameters:
    3596: +  mat - the Mat object
    3597: .  col - the column index
    3598: -  v - the Vec object

    3600:    Level: intermediate

    3602: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecWrite()
    3603: @*/
    3604: PetscErrorCode MatDenseRestoreColumnVecRead(Mat A,PetscInt col,Vec *v)
    3605: {

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

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

    3622:    Collective

    3624:    Input Parameters:
    3625: +  mat - the Mat object
    3626: -  col - the column index

    3628:    Output Parameter:
    3629: .  v - the vector

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

    3635:    Level: intermediate

    3637: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
    3638: @*/
    3639: PetscErrorCode MatDenseGetColumnVecWrite(Mat A,PetscInt col,Vec *v)
    3640: {

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

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

    3657:    Collective

    3659:    Input Parameters:
    3660: +  mat - the Mat object
    3661: .  col - the column index
    3662: -  v - the Vec object

    3664:    Level: intermediate

    3666: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead()
    3667: @*/
    3668: PetscErrorCode MatDenseRestoreColumnVecWrite(Mat A,PetscInt col,Vec *v)
    3669: {

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

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

    3686:    Collective

    3688:    Input Parameters:
    3689: +  mat - the Mat object
    3690: .  cbegin - the first index in the block
    3691: -  cend - the last index in the block

    3693:    Output Parameter:
    3694: .  v - the matrix

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

    3699:    Level: intermediate

    3701: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseRestoreColumnVec(), MatDenseRestoreSubMatrix()
    3702: @*/
    3703: PetscErrorCode MatDenseGetSubMatrix(Mat A,PetscInt cbegin,PetscInt cend,Mat *v)
    3704: {

    3713:   if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
    3714:   if (cbegin < 0 || cbegin > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid cbegin %D, should be in [0,%D)",cbegin,A->cmap->N);
    3715:   if (cend < cbegin || cend > A->cmap->N) SETERRQ3(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid cend %D, should be in [%D,%D)",cend,cbegin,A->cmap->N);
    3716:   PetscUseMethod(A,"MatDenseGetSubMatrix_C",(Mat,PetscInt,PetscInt,Mat*),(A,cbegin,cend,v));
    3717:   return(0);
    3718: }

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

    3723:    Collective

    3725:    Input Parameters:
    3726: +  mat - the Mat object
    3727: -  v - the Mat object

    3729:    Level: intermediate

    3731: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseRestoreColumnVec(), MatDenseGetSubMatrix()
    3732: @*/
    3733: PetscErrorCode MatDenseRestoreSubMatrix(Mat A,Mat *v)
    3734: {

    3741:   PetscUseMethod(A,"MatDenseRestoreSubMatrix_C",(Mat,Mat*),(A,v));
    3742:   return(0);
    3743: }