Actual source code: matelem.cxx

      1: #include <petsc/private/petscelemental.h>

      3: const char ElementalCitation[] = "@Article{Elemental2012,\n"
      4: "  author  = {Jack Poulson and Bryan Marker and Jeff R. Hammond and Nichols A. Romero and Robert {v}an~{d}e~{G}eijn},\n"
      5: "  title   = {Elemental: A New Framework for Distributed Memory Dense Matrix Computations},\n"
      6: "  journal = {{ACM} Transactions on Mathematical Software},\n"
      7: "  volume  = {39},\n"
      8: "  number  = {2},\n"
      9: "  year    = {2013}\n"
     10: "}\n";
     11: static PetscBool ElementalCite = PETSC_FALSE;

     13: /*
     14:     The variable Petsc_Elemental_keyval is used to indicate an MPI attribute that
     15:   is attached to a communicator, in this case the attribute is a Mat_Elemental_Grid
     16: */
     17: static PetscMPIInt Petsc_Elemental_keyval = MPI_KEYVAL_INVALID;

     19: static PetscErrorCode MatView_Elemental(Mat A,PetscViewer viewer)
     20: {
     22:   Mat_Elemental  *a = (Mat_Elemental*)A->data;
     23:   PetscBool      iascii;

     26:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
     27:   if (iascii) {
     28:     PetscViewerFormat format;
     29:     PetscViewerGetFormat(viewer,&format);
     30:     if (format == PETSC_VIEWER_ASCII_INFO) {
     31:       /* call elemental viewing function */
     32:       PetscViewerASCIIPrintf(viewer,"Elemental run parameters:\n");
     33:       PetscViewerASCIIPrintf(viewer,"  allocated entries=%d\n",(*a->emat).AllocatedMemory());
     34:       PetscViewerASCIIPrintf(viewer,"  grid height=%d, grid width=%d\n",(*a->emat).Grid().Height(),(*a->emat).Grid().Width());
     35:       if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) {
     36:         /* call elemental viewing function */
     37:         PetscPrintf(PetscObjectComm((PetscObject)viewer),"test matview_elemental 2\n");
     38:       }

     40:     } else if (format == PETSC_VIEWER_DEFAULT) {
     41:       PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
     42:       El::Print( *a->emat, "Elemental matrix (cyclic ordering)");
     43:       PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
     44:       if (A->factortype == MAT_FACTOR_NONE){
     45:         Mat Adense;
     46:         MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&Adense);
     47:         MatView(Adense,viewer);
     48:         MatDestroy(&Adense);
     49:       }
     50:     } else SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Format");
     51:   } else {
     52:     /* convert to dense format and call MatView() */
     53:     Mat Adense;
     54:     MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&Adense);
     55:     MatView(Adense,viewer);
     56:     MatDestroy(&Adense);
     57:   }
     58:   return(0);
     59: }

     61: static PetscErrorCode MatGetInfo_Elemental(Mat A,MatInfoType flag,MatInfo *info)
     62: {
     63:   Mat_Elemental  *a = (Mat_Elemental*)A->data;

     66:   info->block_size = 1.0;

     68:   if (flag == MAT_LOCAL) {
     69:     info->nz_allocated   = (*a->emat).AllocatedMemory(); /* locally allocated */
     70:     info->nz_used        = info->nz_allocated;
     71:   } else if (flag == MAT_GLOBAL_MAX) {
     72:     //MPIU_Allreduce(isend,irecv,5,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)matin));
     73:     /* see MatGetInfo_MPIAIJ() for getting global info->nz_allocated! */
     74:     //SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP," MAT_GLOBAL_MAX not written yet");
     75:   } else if (flag == MAT_GLOBAL_SUM) {
     76:     //SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP," MAT_GLOBAL_SUM not written yet");
     77:     info->nz_allocated   = (*a->emat).AllocatedMemory(); /* locally allocated */
     78:     info->nz_used        = info->nz_allocated; /* assume Elemental does accurate allocation */
     79:     //MPIU_Allreduce(isend,irecv,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)A));
     80:     //PetscPrintf(PETSC_COMM_SELF,"    ... [%d] locally allocated %g\n",rank,info->nz_allocated);
     81:   }

     83:   info->nz_unneeded       = 0.0;
     84:   info->assemblies        = A->num_ass;
     85:   info->mallocs           = 0;
     86:   info->memory            = ((PetscObject)A)->mem;
     87:   info->fill_ratio_given  = 0; /* determined by Elemental */
     88:   info->fill_ratio_needed = 0;
     89:   info->factor_mallocs    = 0;
     90:   return(0);
     91: }

     93: PetscErrorCode MatSetOption_Elemental(Mat A,MatOption op,PetscBool flg)
     94: {
     95:   Mat_Elemental  *a = (Mat_Elemental*)A->data;

     98:   switch (op) {
     99:   case MAT_NEW_NONZERO_LOCATIONS:
    100:   case MAT_NEW_NONZERO_LOCATION_ERR:
    101:   case MAT_NEW_NONZERO_ALLOCATION_ERR:
    102:   case MAT_SYMMETRIC:
    103:   case MAT_SORTED_FULL:
    104:   case MAT_HERMITIAN:
    105:     break;
    106:   case MAT_ROW_ORIENTED:
    107:     a->roworiented = flg;
    108:     break;
    109:   default:
    110:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]);
    111:   }
    112:   return(0);
    113: }

    115: static PetscErrorCode MatSetValues_Elemental(Mat A,PetscInt nr,const PetscInt *rows,PetscInt nc,const PetscInt *cols,const PetscScalar *vals,InsertMode imode)
    116: {
    117:   Mat_Elemental  *a = (Mat_Elemental*)A->data;
    118:   PetscInt       i,j,rrank,ridx,crank,cidx,erow,ecol,numQueues=0;

    121:   // TODO: Initialize matrix to all zeros?

    123:   // Count the number of queues from this process
    124:   if (a->roworiented) {
    125:     for (i=0; i<nr; i++) {
    126:       if (rows[i] < 0) continue;
    127:       P2RO(A,0,rows[i],&rrank,&ridx);
    128:       RO2E(A,0,rrank,ridx,&erow);
    129:       if (rrank < 0 || ridx < 0 || erow < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Incorrect row translation");
    130:       for (j=0; j<nc; j++) {
    131:         if (cols[j] < 0) continue;
    132:         P2RO(A,1,cols[j],&crank,&cidx);
    133:         RO2E(A,1,crank,cidx,&ecol);
    134:         if (crank < 0 || cidx < 0 || ecol < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Incorrect col translation");
    135:         if (!a->emat->IsLocal(erow,ecol)){ /* off-proc entry */
    136:           /* printf("Will later remotely update (%d,%d)\n",erow,ecol); */
    137:           if (imode != ADD_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only ADD_VALUES to off-processor entry is supported");
    138:           ++numQueues;
    139:           continue;
    140:         }
    141:         /* printf("Locally updating (%d,%d)\n",erow,ecol); */
    142:         switch (imode) {
    143:         case INSERT_VALUES: a->emat->Set(erow,ecol,(PetscElemScalar)vals[i*nc+j]); break;
    144:         case ADD_VALUES: a->emat->Update(erow,ecol,(PetscElemScalar)vals[i*nc+j]); break;
    145:         default: SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for InsertMode %d",(int)imode);
    146:         }
    147:       }
    148:     }

    150:     /* printf("numQueues=%d\n",numQueues); */
    151:     a->emat->Reserve( numQueues);
    152:     for (i=0; i<nr; i++) {
    153:       if (rows[i] < 0) continue;
    154:       P2RO(A,0,rows[i],&rrank,&ridx);
    155:       RO2E(A,0,rrank,ridx,&erow);
    156:       for (j=0; j<nc; j++) {
    157:         if (cols[j] < 0) continue;
    158:         P2RO(A,1,cols[j],&crank,&cidx);
    159:         RO2E(A,1,crank,cidx,&ecol);
    160:         if (!a->emat->IsLocal(erow,ecol)) { /*off-proc entry*/
    161:           /* printf("Queueing remotely update of (%d,%d)\n",erow,ecol); */
    162:           a->emat->QueueUpdate( erow, ecol, vals[i*nc+j]);
    163:         }
    164:       }
    165:     }
    166:   } else { /* columnoriented */
    167:     for (j=0; j<nc; j++) {
    168:       if (cols[j] < 0) continue;
    169:       P2RO(A,1,cols[j],&crank,&cidx);
    170:       RO2E(A,1,crank,cidx,&ecol);
    171:       if (crank < 0 || cidx < 0 || ecol < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Incorrect col translation");
    172:       for (i=0; i<nr; i++) {
    173:         if (rows[i] < 0) continue;
    174:         P2RO(A,0,rows[i],&rrank,&ridx);
    175:         RO2E(A,0,rrank,ridx,&erow);
    176:         if (rrank < 0 || ridx < 0 || erow < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Incorrect row translation");
    177:         if (!a->emat->IsLocal(erow,ecol)){ /* off-proc entry */
    178:           /* printf("Will later remotely update (%d,%d)\n",erow,ecol); */
    179:           if (imode != ADD_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only ADD_VALUES to off-processor entry is supported");
    180:           ++numQueues;
    181:           continue;
    182:         }
    183:         /* printf("Locally updating (%d,%d)\n",erow,ecol); */
    184:         switch (imode) {
    185:         case INSERT_VALUES: a->emat->Set(erow,ecol,(PetscElemScalar)vals[i+j*nr]); break;
    186:         case ADD_VALUES: a->emat->Update(erow,ecol,(PetscElemScalar)vals[i+j*nr]); break;
    187:         default: SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for InsertMode %d",(int)imode);
    188:         }
    189:       }
    190:     }

    192:     /* printf("numQueues=%d\n",numQueues); */
    193:     a->emat->Reserve( numQueues);
    194:     for (j=0; j<nc; j++) {
    195:       if (cols[j] < 0) continue;
    196:       P2RO(A,1,cols[j],&crank,&cidx);
    197:       RO2E(A,1,crank,cidx,&ecol);

    199:       for (i=0; i<nr; i++) {
    200:         if (rows[i] < 0) continue;
    201:         P2RO(A,0,rows[i],&rrank,&ridx);
    202:         RO2E(A,0,rrank,ridx,&erow);
    203:         if (!a->emat->IsLocal(erow,ecol)) { /*off-proc entry*/
    204:           /* printf("Queueing remotely update of (%d,%d)\n",erow,ecol); */
    205:           a->emat->QueueUpdate( erow, ecol, vals[i+j*nr]);
    206:         }
    207:       }
    208:     }
    209:   }
    210:   return(0);
    211: }

    213: static PetscErrorCode MatMult_Elemental(Mat A,Vec X,Vec Y)
    214: {
    215:   Mat_Elemental         *a = (Mat_Elemental*)A->data;
    216:   PetscErrorCode        ierr;
    217:   const PetscElemScalar *x;
    218:   PetscElemScalar       *y;
    219:   PetscElemScalar       one = 1,zero = 0;

    222:   VecGetArrayRead(X,(const PetscScalar **)&x);
    223:   VecGetArray(Y,(PetscScalar **)&y);
    224:   { /* Scoping so that constructor is called before pointer is returned */
    225:     El::DistMatrix<PetscElemScalar,El::VC,El::STAR> xe, ye;
    226:     xe.LockedAttach(A->cmap->N,1,*a->grid,0,0,x,A->cmap->n);
    227:     ye.Attach(A->rmap->N,1,*a->grid,0,0,y,A->rmap->n);
    228:     El::Gemv(El::NORMAL,one,*a->emat,xe,zero,ye);
    229:   }
    230:   VecRestoreArrayRead(X,(const PetscScalar **)&x);
    231:   VecRestoreArray(Y,(PetscScalar **)&y);
    232:   return(0);
    233: }

    235: static PetscErrorCode MatMultTranspose_Elemental(Mat A,Vec X,Vec Y)
    236: {
    237:   Mat_Elemental         *a = (Mat_Elemental*)A->data;
    238:   PetscErrorCode        ierr;
    239:   const PetscElemScalar *x;
    240:   PetscElemScalar       *y;
    241:   PetscElemScalar       one = 1,zero = 0;

    244:   VecGetArrayRead(X,(const PetscScalar **)&x);
    245:   VecGetArray(Y,(PetscScalar **)&y);
    246:   { /* Scoping so that constructor is called before pointer is returned */
    247:     El::DistMatrix<PetscElemScalar,El::VC,El::STAR> xe, ye;
    248:     xe.LockedAttach(A->rmap->N,1,*a->grid,0,0,x,A->rmap->n);
    249:     ye.Attach(A->cmap->N,1,*a->grid,0,0,y,A->cmap->n);
    250:     El::Gemv(El::TRANSPOSE,one,*a->emat,xe,zero,ye);
    251:   }
    252:   VecRestoreArrayRead(X,(const PetscScalar **)&x);
    253:   VecRestoreArray(Y,(PetscScalar **)&y);
    254:   return(0);
    255: }

    257: static PetscErrorCode MatMultAdd_Elemental(Mat A,Vec X,Vec Y,Vec Z)
    258: {
    259:   Mat_Elemental         *a = (Mat_Elemental*)A->data;
    260:   PetscErrorCode        ierr;
    261:   const PetscElemScalar *x;
    262:   PetscElemScalar       *z;
    263:   PetscElemScalar       one = 1;

    266:   if (Y != Z) {VecCopy(Y,Z);}
    267:   VecGetArrayRead(X,(const PetscScalar **)&x);
    268:   VecGetArray(Z,(PetscScalar **)&z);
    269:   { /* Scoping so that constructor is called before pointer is returned */
    270:     El::DistMatrix<PetscElemScalar,El::VC,El::STAR> xe, ze;
    271:     xe.LockedAttach(A->cmap->N,1,*a->grid,0,0,x,A->cmap->n);
    272:     ze.Attach(A->rmap->N,1,*a->grid,0,0,z,A->rmap->n);
    273:     El::Gemv(El::NORMAL,one,*a->emat,xe,one,ze);
    274:   }
    275:   VecRestoreArrayRead(X,(const PetscScalar **)&x);
    276:   VecRestoreArray(Z,(PetscScalar **)&z);
    277:   return(0);
    278: }

    280: static PetscErrorCode MatMultTransposeAdd_Elemental(Mat A,Vec X,Vec Y,Vec Z)
    281: {
    282:   Mat_Elemental         *a = (Mat_Elemental*)A->data;
    283:   PetscErrorCode        ierr;
    284:   const PetscElemScalar *x;
    285:   PetscElemScalar       *z;
    286:   PetscElemScalar       one = 1;

    289:   if (Y != Z) {VecCopy(Y,Z);}
    290:   VecGetArrayRead(X,(const PetscScalar **)&x);
    291:   VecGetArray(Z,(PetscScalar **)&z);
    292:   { /* Scoping so that constructor is called before pointer is returned */
    293:     El::DistMatrix<PetscElemScalar,El::VC,El::STAR> xe, ze;
    294:     xe.LockedAttach(A->rmap->N,1,*a->grid,0,0,x,A->rmap->n);
    295:     ze.Attach(A->cmap->N,1,*a->grid,0,0,z,A->cmap->n);
    296:     El::Gemv(El::TRANSPOSE,one,*a->emat,xe,one,ze);
    297:   }
    298:   VecRestoreArrayRead(X,(const PetscScalar **)&x);
    299:   VecRestoreArray(Z,(PetscScalar **)&z);
    300:   return(0);
    301: }

    303: PetscErrorCode MatMatMultNumeric_Elemental(Mat A,Mat B,Mat C)
    304: {
    305:   Mat_Elemental    *a = (Mat_Elemental*)A->data;
    306:   Mat_Elemental    *b = (Mat_Elemental*)B->data;
    307:   Mat_Elemental    *c = (Mat_Elemental*)C->data;
    308:   PetscElemScalar  one = 1,zero = 0;

    311:   { /* Scoping so that constructor is called before pointer is returned */
    312:     El::Gemm(El::NORMAL,El::NORMAL,one,*a->emat,*b->emat,zero,*c->emat);
    313:   }
    314:   C->assembled = PETSC_TRUE;
    315:   return(0);
    316: }

    318: PetscErrorCode MatMatMultSymbolic_Elemental(Mat A,Mat B,PetscReal fill,Mat Ce)
    319: {

    323:   MatSetSizes(Ce,A->rmap->n,B->cmap->n,PETSC_DECIDE,PETSC_DECIDE);
    324:   MatSetType(Ce,MATELEMENTAL);
    325:   MatSetUp(Ce);
    326:   Ce->ops->matmultnumeric = MatMatMultNumeric_Elemental;
    327:   return(0);
    328: }

    330: static PetscErrorCode MatMatTransposeMultNumeric_Elemental(Mat A,Mat B,Mat C)
    331: {
    332:   Mat_Elemental      *a = (Mat_Elemental*)A->data;
    333:   Mat_Elemental      *b = (Mat_Elemental*)B->data;
    334:   Mat_Elemental      *c = (Mat_Elemental*)C->data;
    335:   PetscElemScalar    one = 1,zero = 0;

    338:   { /* Scoping so that constructor is called before pointer is returned */
    339:     El::Gemm(El::NORMAL,El::TRANSPOSE,one,*a->emat,*b->emat,zero,*c->emat);
    340:   }
    341:   C->assembled = PETSC_TRUE;
    342:   return(0);
    343: }

    345: static PetscErrorCode MatMatTransposeMultSymbolic_Elemental(Mat A,Mat B,PetscReal fill,Mat C)
    346: {

    350:   MatSetSizes(C,A->rmap->n,B->rmap->n,PETSC_DECIDE,PETSC_DECIDE);
    351:   MatSetType(C,MATELEMENTAL);
    352:   MatSetUp(C);
    353:   return(0);
    354: }

    356: /* --------------------------------------- */
    357: static PetscErrorCode MatProductSetFromOptions_Elemental_AB(Mat C)
    358: {
    360:   C->ops->matmultsymbolic = MatMatMultSymbolic_Elemental;
    361:   C->ops->productsymbolic = MatProductSymbolic_AB;
    362:   return(0);
    363: }

    365: static PetscErrorCode MatProductSetFromOptions_Elemental_ABt(Mat C)
    366: {
    368:   C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_Elemental;
    369:   C->ops->productsymbolic          = MatProductSymbolic_ABt;
    370:   return(0);
    371: }

    373: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Elemental(Mat C)
    374: {
    376:   Mat_Product    *product = C->product;

    379:   switch (product->type) {
    380:   case MATPRODUCT_AB:
    381:     MatProductSetFromOptions_Elemental_AB(C);
    382:     break;
    383:   case MATPRODUCT_ABt:
    384:     MatProductSetFromOptions_Elemental_ABt(C);
    385:     break;
    386:   default:
    387:     break;
    388:   }
    389:   return(0);
    390: }

    392: PetscErrorCode MatMatMultNumeric_Elemental_MPIDense(Mat A,Mat B,Mat C)
    393: {
    394:   Mat            Be,Ce;

    398:   MatConvert(B,MATELEMENTAL,MAT_INITIAL_MATRIX,&Be);
    399:   MatMatMult(A,Be,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&Ce);
    400:   MatConvert(Ce,MATMPIDENSE,MAT_REUSE_MATRIX,&C);
    401:   MatDestroy(&Be);
    402:   MatDestroy(&Ce);
    403:   return(0);
    404: }

    406: PetscErrorCode MatMatMultSymbolic_Elemental_MPIDense(Mat A,Mat B,PetscReal fill,Mat C)
    407: {

    411:   MatSetSizes(C,A->rmap->n,B->cmap->n,PETSC_DECIDE,PETSC_DECIDE);
    412:   MatSetType(C,MATMPIDENSE);
    413:   MatSetUp(C);
    414:   C->ops->matmultnumeric = MatMatMultNumeric_Elemental_MPIDense;
    415:   return(0);
    416: }

    418: PetscErrorCode MatProductSetFromOptions_Elemental_MPIDense_AB(Mat C)
    419: {
    421:   C->ops->matmultsymbolic = MatMatMultSymbolic_Elemental_MPIDense;
    422:   C->ops->productsymbolic = MatProductSymbolic_AB;
    423:   return(0);
    424: }

    426: PetscErrorCode MatProductSetFromOptions_Elemental_MPIDense(Mat C)
    427: {
    429:   Mat_Product    *product = C->product;

    432:   if (product->type == MATPRODUCT_AB) {
    433:     MatProductSetFromOptions_Elemental_MPIDense_AB(C);
    434:   }
    435:   return(0);
    436: }
    437: /* --------------------------------------- */

    439: static PetscErrorCode MatGetDiagonal_Elemental(Mat A,Vec D)
    440: {
    441:   PetscInt        i,nrows,ncols,nD,rrank,ridx,crank,cidx;
    442:   Mat_Elemental   *a = (Mat_Elemental*)A->data;
    443:   PetscErrorCode  ierr;
    444:   PetscElemScalar v;
    445:   MPI_Comm        comm;

    448:   PetscObjectGetComm((PetscObject)A,&comm);
    449:   MatGetSize(A,&nrows,&ncols);
    450:   nD = nrows>ncols ? ncols : nrows;
    451:   for (i=0; i<nD; i++) {
    452:     PetscInt erow,ecol;
    453:     P2RO(A,0,i,&rrank,&ridx);
    454:     RO2E(A,0,rrank,ridx,&erow);
    455:     if (rrank < 0 || ridx < 0 || erow < 0) SETERRQ(comm,PETSC_ERR_PLIB,"Incorrect row translation");
    456:     P2RO(A,1,i,&crank,&cidx);
    457:     RO2E(A,1,crank,cidx,&ecol);
    458:     if (crank < 0 || cidx < 0 || ecol < 0) SETERRQ(comm,PETSC_ERR_PLIB,"Incorrect col translation");
    459:     v = a->emat->Get(erow,ecol);
    460:     VecSetValues(D,1,&i,(PetscScalar*)&v,INSERT_VALUES);
    461:   }
    462:   VecAssemblyBegin(D);
    463:   VecAssemblyEnd(D);
    464:   return(0);
    465: }

    467: static PetscErrorCode MatDiagonalScale_Elemental(Mat X,Vec L,Vec R)
    468: {
    469:   Mat_Elemental         *x = (Mat_Elemental*)X->data;
    470:   const PetscElemScalar *d;
    471:   PetscErrorCode        ierr;

    474:   if (R) {
    475:     VecGetArrayRead(R,(const PetscScalar **)&d);
    476:     El::DistMatrix<PetscElemScalar,El::VC,El::STAR> de;
    477:     de.LockedAttach(X->cmap->N,1,*x->grid,0,0,d,X->cmap->n);
    478:     El::DiagonalScale(El::RIGHT,El::NORMAL,de,*x->emat);
    479:     VecRestoreArrayRead(R,(const PetscScalar **)&d);
    480:   }
    481:   if (L) {
    482:     VecGetArrayRead(L,(const PetscScalar **)&d);
    483:     El::DistMatrix<PetscElemScalar,El::VC,El::STAR> de;
    484:     de.LockedAttach(X->rmap->N,1,*x->grid,0,0,d,X->rmap->n);
    485:     El::DiagonalScale(El::LEFT,El::NORMAL,de,*x->emat);
    486:     VecRestoreArrayRead(L,(const PetscScalar **)&d);
    487:   }
    488:   return(0);
    489: }

    491: static PetscErrorCode MatMissingDiagonal_Elemental(Mat A,PetscBool *missing,PetscInt *d)
    492: {
    494:   *missing = PETSC_FALSE;
    495:   return(0);
    496: }

    498: static PetscErrorCode MatScale_Elemental(Mat X,PetscScalar a)
    499: {
    500:   Mat_Elemental  *x = (Mat_Elemental*)X->data;

    503:   El::Scale((PetscElemScalar)a,*x->emat);
    504:   return(0);
    505: }

    507: /*
    508:   MatAXPY - Computes Y = a*X + Y.
    509: */
    510: static PetscErrorCode MatAXPY_Elemental(Mat Y,PetscScalar a,Mat X,MatStructure str)
    511: {
    512:   Mat_Elemental  *x = (Mat_Elemental*)X->data;
    513:   Mat_Elemental  *y = (Mat_Elemental*)Y->data;

    517:   El::Axpy((PetscElemScalar)a,*x->emat,*y->emat);
    518:   PetscObjectStateIncrease((PetscObject)Y);
    519:   return(0);
    520: }

    522: static PetscErrorCode MatCopy_Elemental(Mat A,Mat B,MatStructure str)
    523: {
    524:   Mat_Elemental *a=(Mat_Elemental*)A->data;
    525:   Mat_Elemental *b=(Mat_Elemental*)B->data;

    529:   El::Copy(*a->emat,*b->emat);
    530:   PetscObjectStateIncrease((PetscObject)B);
    531:   return(0);
    532: }

    534: static PetscErrorCode MatDuplicate_Elemental(Mat A,MatDuplicateOption op,Mat *B)
    535: {
    536:   Mat            Be;
    537:   MPI_Comm       comm;
    538:   Mat_Elemental  *a=(Mat_Elemental*)A->data;

    542:   PetscObjectGetComm((PetscObject)A,&comm);
    543:   MatCreate(comm,&Be);
    544:   MatSetSizes(Be,A->rmap->n,A->cmap->n,PETSC_DECIDE,PETSC_DECIDE);
    545:   MatSetType(Be,MATELEMENTAL);
    546:   MatSetUp(Be);
    547:   *B = Be;
    548:   if (op == MAT_COPY_VALUES) {
    549:     Mat_Elemental *b=(Mat_Elemental*)Be->data;
    550:     El::Copy(*a->emat,*b->emat);
    551:   }
    552:   Be->assembled = PETSC_TRUE;
    553:   return(0);
    554: }

    556: static PetscErrorCode MatTranspose_Elemental(Mat A,MatReuse reuse,Mat *B)
    557: {
    558:   Mat            Be = *B;
    560:   MPI_Comm       comm;
    561:   Mat_Elemental  *a = (Mat_Elemental*)A->data, *b;

    564:   PetscObjectGetComm((PetscObject)A,&comm);
    565:   /* Only out-of-place supported */
    566:   if (reuse == MAT_INPLACE_MATRIX) SETERRQ(comm,PETSC_ERR_SUP,"Only out-of-place supported");
    567:   if (reuse == MAT_INITIAL_MATRIX) {
    568:     MatCreate(comm,&Be);
    569:     MatSetSizes(Be,A->cmap->n,A->rmap->n,PETSC_DECIDE,PETSC_DECIDE);
    570:     MatSetType(Be,MATELEMENTAL);
    571:     MatSetUp(Be);
    572:     *B = Be;
    573:   }
    574:   b = (Mat_Elemental*)Be->data;
    575:   El::Transpose(*a->emat,*b->emat);
    576:   Be->assembled = PETSC_TRUE;
    577:   return(0);
    578: }

    580: static PetscErrorCode MatConjugate_Elemental(Mat A)
    581: {
    582:   Mat_Elemental  *a = (Mat_Elemental*)A->data;

    585:   El::Conjugate(*a->emat);
    586:   return(0);
    587: }

    589: static PetscErrorCode MatHermitianTranspose_Elemental(Mat A,MatReuse reuse,Mat *B)
    590: {
    591:   Mat            Be = *B;
    593:   MPI_Comm       comm;
    594:   Mat_Elemental  *a = (Mat_Elemental*)A->data, *b;

    597:   PetscObjectGetComm((PetscObject)A,&comm);
    598:   /* Only out-of-place supported */
    599:   if (reuse == MAT_INITIAL_MATRIX){
    600:     MatCreate(comm,&Be);
    601:     MatSetSizes(Be,A->cmap->n,A->rmap->n,PETSC_DECIDE,PETSC_DECIDE);
    602:     MatSetType(Be,MATELEMENTAL);
    603:     MatSetUp(Be);
    604:     *B = Be;
    605:   }
    606:   b = (Mat_Elemental*)Be->data;
    607:   El::Adjoint(*a->emat,*b->emat);
    608:   Be->assembled = PETSC_TRUE;
    609:   return(0);
    610: }

    612: static PetscErrorCode MatSolve_Elemental(Mat A,Vec B,Vec X)
    613: {
    614:   Mat_Elemental     *a = (Mat_Elemental*)A->data;
    615:   PetscErrorCode    ierr;
    616:   PetscElemScalar   *x;
    617:   PetscInt          pivoting = a->pivoting;

    620:   VecCopy(B,X);
    621:   VecGetArray(X,(PetscScalar **)&x);

    623:   El::DistMatrix<PetscElemScalar,El::VC,El::STAR> xe;
    624:   xe.Attach(A->rmap->N,1,*a->grid,0,0,x,A->rmap->n);
    625:   El::DistMatrix<PetscElemScalar,El::MC,El::MR> xer(xe);
    626:   switch (A->factortype) {
    627:   case MAT_FACTOR_LU:
    628:     if (pivoting == 0) {
    629:       El::lu::SolveAfter(El::NORMAL,*a->emat,xer);
    630:     } else if (pivoting == 1) {
    631:       El::lu::SolveAfter(El::NORMAL,*a->emat,*a->P,xer);
    632:     } else { /* pivoting == 2 */
    633:       El::lu::SolveAfter(El::NORMAL,*a->emat,*a->P,*a->Q,xer);
    634:     }
    635:     break;
    636:   case MAT_FACTOR_CHOLESKY:
    637:     El::cholesky::SolveAfter(El::UPPER,El::NORMAL,*a->emat,xer);
    638:     break;
    639:   default:
    640:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unfactored Matrix or Unsupported MatFactorType");
    641:     break;
    642:   }
    643:   El::Copy(xer,xe);

    645:   VecRestoreArray(X,(PetscScalar **)&x);
    646:   return(0);
    647: }

    649: static PetscErrorCode MatSolveAdd_Elemental(Mat A,Vec B,Vec Y,Vec X)
    650: {
    651:   PetscErrorCode    ierr;

    654:   MatSolve_Elemental(A,B,X);
    655:   VecAXPY(X,1,Y);
    656:   return(0);
    657: }

    659: static PetscErrorCode MatMatSolve_Elemental(Mat A,Mat B,Mat X)
    660: {
    661:   Mat_Elemental  *a = (Mat_Elemental*)A->data;
    662:   Mat_Elemental  *x;
    663:   Mat            C;
    664:   PetscInt       pivoting = a->pivoting;
    665:   PetscBool      flg;
    666:   MatType        type;

    670:   MatGetType(X,&type);
    671:   PetscStrcmp(type,MATELEMENTAL,&flg);
    672:   if (!flg) {
    673:     MatConvert(B,MATELEMENTAL,MAT_INITIAL_MATRIX,&C);
    674:     x = (Mat_Elemental*)C->data;
    675:   } else {
    676:     x = (Mat_Elemental*)X->data;
    677:     El::Copy(*((Mat_Elemental*)B->data)->emat,*x->emat);
    678:   }
    679:   switch (A->factortype) {
    680:   case MAT_FACTOR_LU:
    681:     if (pivoting == 0) {
    682:       El::lu::SolveAfter(El::NORMAL,*a->emat,*x->emat);
    683:     } else if (pivoting == 1) {
    684:       El::lu::SolveAfter(El::NORMAL,*a->emat,*a->P,*x->emat);
    685:     } else {
    686:       El::lu::SolveAfter(El::NORMAL,*a->emat,*a->P,*a->Q,*x->emat);
    687:     }
    688:     break;
    689:   case MAT_FACTOR_CHOLESKY:
    690:     El::cholesky::SolveAfter(El::UPPER,El::NORMAL,*a->emat,*x->emat);
    691:     break;
    692:   default:
    693:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unfactored Matrix or Unsupported MatFactorType");
    694:     break;
    695:   }
    696:   if (!flg) {
    697:     MatConvert(C,type,MAT_REUSE_MATRIX,&X);
    698:     MatDestroy(&C);
    699:   }
    700:   return(0);
    701: }

    703: static PetscErrorCode MatLUFactor_Elemental(Mat A,IS row,IS col,const MatFactorInfo *info)
    704: {
    705:   Mat_Elemental  *a = (Mat_Elemental*)A->data;
    707:   PetscInt       pivoting = a->pivoting;

    710:   if (pivoting == 0) {
    711:     El::LU(*a->emat);
    712:   } else if (pivoting == 1) {
    713:     El::LU(*a->emat,*a->P);
    714:   } else {
    715:     El::LU(*a->emat,*a->P,*a->Q);
    716:   }
    717:   A->factortype = MAT_FACTOR_LU;
    718:   A->assembled  = PETSC_TRUE;

    720:   PetscFree(A->solvertype);
    721:   PetscStrallocpy(MATSOLVERELEMENTAL,&A->solvertype);
    722:   return(0);
    723: }

    725: static PetscErrorCode  MatLUFactorNumeric_Elemental(Mat F,Mat A,const MatFactorInfo *info)
    726: {

    730:   MatCopy(A,F,SAME_NONZERO_PATTERN);
    731:   MatLUFactor_Elemental(F,0,0,info);
    732:   return(0);
    733: }

    735: static PetscErrorCode  MatLUFactorSymbolic_Elemental(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
    736: {
    738:   /* F is created and allocated by MatGetFactor_elemental_petsc(), skip this routine. */
    739:   return(0);
    740: }

    742: static PetscErrorCode MatCholeskyFactor_Elemental(Mat A,IS perm,const MatFactorInfo *info)
    743: {
    744:   Mat_Elemental  *a = (Mat_Elemental*)A->data;
    745:   El::DistMatrix<PetscElemScalar,El::MC,El::STAR> d;

    749:   El::Cholesky(El::UPPER,*a->emat);
    750:   A->factortype = MAT_FACTOR_CHOLESKY;
    751:   A->assembled  = PETSC_TRUE;

    753:   PetscFree(A->solvertype);
    754:   PetscStrallocpy(MATSOLVERELEMENTAL,&A->solvertype);
    755:   return(0);
    756: }

    758: static PetscErrorCode MatCholeskyFactorNumeric_Elemental(Mat F,Mat A,const MatFactorInfo *info)
    759: {

    763:   MatCopy(A,F,SAME_NONZERO_PATTERN);
    764:   MatCholeskyFactor_Elemental(F,0,info);
    765:   return(0);
    766: }

    768: static PetscErrorCode MatCholeskyFactorSymbolic_Elemental(Mat F,Mat A,IS perm,const MatFactorInfo *info)
    769: {
    771:   /* F is created and allocated by MatGetFactor_elemental_petsc(), skip this routine. */
    772:   return(0);
    773: }

    775: PetscErrorCode MatFactorGetSolverType_elemental_elemental(Mat A,MatSolverType *type)
    776: {
    778:   *type = MATSOLVERELEMENTAL;
    779:   return(0);
    780: }

    782: static PetscErrorCode MatGetFactor_elemental_elemental(Mat A,MatFactorType ftype,Mat *F)
    783: {
    784:   Mat            B;

    788:   /* Create the factorization matrix */
    789:   MatCreate(PetscObjectComm((PetscObject)A),&B);
    790:   MatSetSizes(B,A->rmap->n,A->cmap->n,PETSC_DECIDE,PETSC_DECIDE);
    791:   MatSetType(B,MATELEMENTAL);
    792:   MatSetUp(B);
    793:   B->factortype = ftype;
    794:   PetscFree(B->solvertype);
    795:   PetscStrallocpy(MATSOLVERELEMENTAL,&B->solvertype);

    797:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_elemental_elemental);
    798:   *F            = B;
    799:   return(0);
    800: }

    802: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_Elemental(void)
    803: {

    807:   MatSolverTypeRegister(MATSOLVERELEMENTAL,MATELEMENTAL,        MAT_FACTOR_LU,MatGetFactor_elemental_elemental);
    808:   MatSolverTypeRegister(MATSOLVERELEMENTAL,MATELEMENTAL,        MAT_FACTOR_CHOLESKY,MatGetFactor_elemental_elemental);
    809:   return(0);
    810: }

    812: static PetscErrorCode MatNorm_Elemental(Mat A,NormType type,PetscReal *nrm)
    813: {
    814:   Mat_Elemental *a=(Mat_Elemental*)A->data;

    817:   switch (type){
    818:   case NORM_1:
    819:     *nrm = El::OneNorm(*a->emat);
    820:     break;
    821:   case NORM_FROBENIUS:
    822:     *nrm = El::FrobeniusNorm(*a->emat);
    823:     break;
    824:   case NORM_INFINITY:
    825:     *nrm = El::InfinityNorm(*a->emat);
    826:     break;
    827:   default:
    828:     SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Unsupported norm type");
    829:   }
    830:   return(0);
    831: }

    833: static PetscErrorCode MatZeroEntries_Elemental(Mat A)
    834: {
    835:   Mat_Elemental *a=(Mat_Elemental*)A->data;

    838:   El::Zero(*a->emat);
    839:   return(0);
    840: }

    842: static PetscErrorCode MatGetOwnershipIS_Elemental(Mat A,IS *rows,IS *cols)
    843: {
    844:   Mat_Elemental  *a = (Mat_Elemental*)A->data;
    846:   PetscInt       i,m,shift,stride,*idx;

    849:   if (rows) {
    850:     m = a->emat->LocalHeight();
    851:     shift = a->emat->ColShift();
    852:     stride = a->emat->ColStride();
    853:     PetscMalloc1(m,&idx);
    854:     for (i=0; i<m; i++) {
    855:       PetscInt rank,offset;
    856:       E2RO(A,0,shift+i*stride,&rank,&offset);
    857:       RO2P(A,0,rank,offset,&idx[i]);
    858:     }
    859:     ISCreateGeneral(PETSC_COMM_SELF,m,idx,PETSC_OWN_POINTER,rows);
    860:   }
    861:   if (cols) {
    862:     m = a->emat->LocalWidth();
    863:     shift = a->emat->RowShift();
    864:     stride = a->emat->RowStride();
    865:     PetscMalloc1(m,&idx);
    866:     for (i=0; i<m; i++) {
    867:       PetscInt rank,offset;
    868:       E2RO(A,1,shift+i*stride,&rank,&offset);
    869:       RO2P(A,1,rank,offset,&idx[i]);
    870:     }
    871:     ISCreateGeneral(PETSC_COMM_SELF,m,idx,PETSC_OWN_POINTER,cols);
    872:   }
    873:   return(0);
    874: }

    876: static PetscErrorCode MatConvert_Elemental_Dense(Mat A,MatType newtype,MatReuse reuse,Mat *B)
    877: {
    878:   Mat                Bmpi;
    879:   Mat_Elemental      *a = (Mat_Elemental*)A->data;
    880:   MPI_Comm           comm;
    881:   PetscErrorCode     ierr;
    882:   IS                 isrows,iscols;
    883:   PetscInt           rrank,ridx,crank,cidx,nrows,ncols,i,j,erow,ecol,elrow,elcol;
    884:   const PetscInt     *rows,*cols;
    885:   PetscElemScalar    v;
    886:   const El::Grid     &grid = a->emat->Grid();

    889:   PetscObjectGetComm((PetscObject)A,&comm);

    891:   if (reuse == MAT_REUSE_MATRIX) {
    892:     Bmpi = *B;
    893:   } else {
    894:     MatCreate(comm,&Bmpi);
    895:     MatSetSizes(Bmpi,A->rmap->n,A->cmap->n,PETSC_DECIDE,PETSC_DECIDE);
    896:     MatSetType(Bmpi,MATDENSE);
    897:     MatSetUp(Bmpi);
    898:   }

    900:   /* Get local entries of A */
    901:   MatGetOwnershipIS(A,&isrows,&iscols);
    902:   ISGetLocalSize(isrows,&nrows);
    903:   ISGetIndices(isrows,&rows);
    904:   ISGetLocalSize(iscols,&ncols);
    905:   ISGetIndices(iscols,&cols);

    907:   if (a->roworiented) {
    908:     for (i=0; i<nrows; i++) {
    909:       P2RO(A,0,rows[i],&rrank,&ridx); /* convert indices between PETSc <-> (Rank,Offset) <-> Elemental */
    910:       RO2E(A,0,rrank,ridx,&erow);
    911:       if (rrank < 0 || ridx < 0 || erow < 0) SETERRQ(comm,PETSC_ERR_PLIB,"Incorrect row translation");
    912:       for (j=0; j<ncols; j++) {
    913:         P2RO(A,1,cols[j],&crank,&cidx);
    914:         RO2E(A,1,crank,cidx,&ecol);
    915:         if (crank < 0 || cidx < 0 || ecol < 0) SETERRQ(comm,PETSC_ERR_PLIB,"Incorrect col translation");

    917:         elrow = erow / grid.MCSize(); /* Elemental local row index */
    918:         elcol = ecol / grid.MRSize(); /* Elemental local column index */
    919:         v = a->emat->GetLocal(elrow,elcol);
    920:         MatSetValues(Bmpi,1,&rows[i],1,&cols[j],(PetscScalar *)&v,INSERT_VALUES);
    921:       }
    922:     }
    923:   } else { /* column-oriented */
    924:     for (j=0; j<ncols; j++) {
    925:       P2RO(A,1,cols[j],&crank,&cidx);
    926:       RO2E(A,1,crank,cidx,&ecol);
    927:       if (crank < 0 || cidx < 0 || ecol < 0) SETERRQ(comm,PETSC_ERR_PLIB,"Incorrect col translation");
    928:       for (i=0; i<nrows; i++) {
    929:         P2RO(A,0,rows[i],&rrank,&ridx); /* convert indices between PETSc <-> (Rank,Offset) <-> Elemental */
    930:         RO2E(A,0,rrank,ridx,&erow);
    931:         if (rrank < 0 || ridx < 0 || erow < 0) SETERRQ(comm,PETSC_ERR_PLIB,"Incorrect row translation");

    933:         elrow = erow / grid.MCSize(); /* Elemental local row index */
    934:         elcol = ecol / grid.MRSize(); /* Elemental local column index */
    935:         v = a->emat->GetLocal(elrow,elcol);
    936:         MatSetValues(Bmpi,1,&rows[i],1,&cols[j],(PetscScalar *)&v,INSERT_VALUES);
    937:       }
    938:     }
    939:   }
    940:   MatAssemblyBegin(Bmpi,MAT_FINAL_ASSEMBLY);
    941:   MatAssemblyEnd(Bmpi,MAT_FINAL_ASSEMBLY);
    942:   if (reuse == MAT_INPLACE_MATRIX) {
    943:     MatHeaderReplace(A,&Bmpi);
    944:   } else {
    945:     *B = Bmpi;
    946:   }
    947:   ISDestroy(&isrows);
    948:   ISDestroy(&iscols);
    949:   return(0);
    950: }

    952: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
    953: {
    954:   Mat               mat_elemental;
    955:   PetscErrorCode    ierr;
    956:   PetscInt          M=A->rmap->N,N=A->cmap->N,row,ncols;
    957:   const PetscInt    *cols;
    958:   const PetscScalar *vals;

    961:   if (reuse == MAT_REUSE_MATRIX) {
    962:     mat_elemental = *newmat;
    963:     MatZeroEntries(mat_elemental);
    964:   } else {
    965:     MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
    966:     MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);
    967:     MatSetType(mat_elemental,MATELEMENTAL);
    968:     MatSetUp(mat_elemental);
    969:   }
    970:   for (row=0; row<M; row++) {
    971:     MatGetRow(A,row,&ncols,&cols,&vals);
    972:     /* PETSc-Elemental interface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
    973:     MatSetValues(mat_elemental,1,&row,ncols,cols,vals,ADD_VALUES);
    974:     MatRestoreRow(A,row,&ncols,&cols,&vals);
    975:   }
    976:   MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
    977:   MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);

    979:   if (reuse == MAT_INPLACE_MATRIX) {
    980:     MatHeaderReplace(A,&mat_elemental);
    981:   } else {
    982:     *newmat = mat_elemental;
    983:   }
    984:   return(0);
    985: }

    987: PETSC_INTERN PetscErrorCode MatConvert_MPIAIJ_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
    988: {
    989:   Mat               mat_elemental;
    990:   PetscErrorCode    ierr;
    991:   PetscInt          row,ncols,rstart=A->rmap->rstart,rend=A->rmap->rend,j;
    992:   const PetscInt    *cols;
    993:   const PetscScalar *vals;

    996:   if (reuse == MAT_REUSE_MATRIX) {
    997:     mat_elemental = *newmat;
    998:     MatZeroEntries(mat_elemental);
    999:   } else {
    1000:     MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
    1001:     MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,A->rmap->N,A->cmap->N);
    1002:     MatSetType(mat_elemental,MATELEMENTAL);
    1003:     MatSetUp(mat_elemental);
    1004:   }
    1005:   for (row=rstart; row<rend; row++) {
    1006:     MatGetRow(A,row,&ncols,&cols,&vals);
    1007:     for (j=0; j<ncols; j++) {
    1008:       /* PETSc-Elemental interface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
    1009:       MatSetValues(mat_elemental,1,&row,1,&cols[j],&vals[j],ADD_VALUES);
    1010:     }
    1011:     MatRestoreRow(A,row,&ncols,&cols,&vals);
    1012:   }
    1013:   MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
    1014:   MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);

    1016:   if (reuse == MAT_INPLACE_MATRIX) {
    1017:     MatHeaderReplace(A,&mat_elemental);
    1018:   } else {
    1019:     *newmat = mat_elemental;
    1020:   }
    1021:   return(0);
    1022: }

    1024: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
    1025: {
    1026:   Mat               mat_elemental;
    1027:   PetscErrorCode    ierr;
    1028:   PetscInt          M=A->rmap->N,N=A->cmap->N,row,ncols,j;
    1029:   const PetscInt    *cols;
    1030:   const PetscScalar *vals;

    1033:   if (reuse == MAT_REUSE_MATRIX) {
    1034:     mat_elemental = *newmat;
    1035:     MatZeroEntries(mat_elemental);
    1036:   } else {
    1037:     MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
    1038:     MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);
    1039:     MatSetType(mat_elemental,MATELEMENTAL);
    1040:     MatSetUp(mat_elemental);
    1041:   }
    1042:   MatGetRowUpperTriangular(A);
    1043:   for (row=0; row<M; row++) {
    1044:     MatGetRow(A,row,&ncols,&cols,&vals);
    1045:     /* PETSc-Elemental interface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
    1046:     MatSetValues(mat_elemental,1,&row,ncols,cols,vals,ADD_VALUES);
    1047:     for (j=0; j<ncols; j++) { /* lower triangular part */
    1048:       PetscScalar v;
    1049:       if (cols[j] == row) continue;
    1050:       v    = A->hermitian ? PetscConj(vals[j]) : vals[j];
    1051:       MatSetValues(mat_elemental,1,&cols[j],1,&row,&v,ADD_VALUES);
    1052:     }
    1053:     MatRestoreRow(A,row,&ncols,&cols,&vals);
    1054:   }
    1055:   MatRestoreRowUpperTriangular(A);
    1056:   MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
    1057:   MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);

    1059:   if (reuse == MAT_INPLACE_MATRIX) {
    1060:     MatHeaderReplace(A,&mat_elemental);
    1061:   } else {
    1062:     *newmat = mat_elemental;
    1063:   }
    1064:   return(0);
    1065: }

    1067: PETSC_INTERN PetscErrorCode MatConvert_MPISBAIJ_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
    1068: {
    1069:   Mat               mat_elemental;
    1070:   PetscErrorCode    ierr;
    1071:   PetscInt          M=A->rmap->N,N=A->cmap->N,row,ncols,j,rstart=A->rmap->rstart,rend=A->rmap->rend;
    1072:   const PetscInt    *cols;
    1073:   const PetscScalar *vals;

    1076:   if (reuse == MAT_REUSE_MATRIX) {
    1077:     mat_elemental = *newmat;
    1078:     MatZeroEntries(mat_elemental);
    1079:   } else {
    1080:     MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
    1081:     MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);
    1082:     MatSetType(mat_elemental,MATELEMENTAL);
    1083:     MatSetUp(mat_elemental);
    1084:   }
    1085:   MatGetRowUpperTriangular(A);
    1086:   for (row=rstart; row<rend; row++) {
    1087:     MatGetRow(A,row,&ncols,&cols,&vals);
    1088:     /* PETSc-Elemental interface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
    1089:     MatSetValues(mat_elemental,1,&row,ncols,cols,vals,ADD_VALUES);
    1090:     for (j=0; j<ncols; j++) { /* lower triangular part */
    1091:       PetscScalar v;
    1092:       if (cols[j] == row) continue;
    1093:       v    = A->hermitian ? PetscConj(vals[j]) : vals[j];
    1094:       MatSetValues(mat_elemental,1,&cols[j],1,&row,&v,ADD_VALUES);
    1095:     }
    1096:     MatRestoreRow(A,row,&ncols,&cols,&vals);
    1097:   }
    1098:   MatRestoreRowUpperTriangular(A);
    1099:   MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
    1100:   MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);

    1102:   if (reuse == MAT_INPLACE_MATRIX) {
    1103:     MatHeaderReplace(A,&mat_elemental);
    1104:   } else {
    1105:     *newmat = mat_elemental;
    1106:   }
    1107:   return(0);
    1108: }

    1110: static PetscErrorCode MatDestroy_Elemental(Mat A)
    1111: {
    1112:   Mat_Elemental      *a = (Mat_Elemental*)A->data;
    1113:   PetscErrorCode     ierr;
    1114:   Mat_Elemental_Grid *commgrid;
    1115:   PetscBool          flg;
    1116:   MPI_Comm           icomm;

    1119:   delete a->emat;
    1120:   delete a->P;
    1121:   delete a->Q;

    1123:   El::mpi::Comm cxxcomm(PetscObjectComm((PetscObject)A));
    1124:   PetscCommDuplicate(cxxcomm.comm,&icomm,NULL);
    1125:   MPI_Comm_get_attr(icomm,Petsc_Elemental_keyval,(void**)&commgrid,(int*)&flg);
    1126:   if (--commgrid->grid_refct == 0) {
    1127:     delete commgrid->grid;
    1128:     PetscFree(commgrid);
    1129:     MPI_Comm_free_keyval(&Petsc_Elemental_keyval);
    1130:   }
    1131:   PetscCommDestroy(&icomm);
    1132:   PetscObjectComposeFunction((PetscObject)A,"MatGetOwnershipIS_C",NULL);
    1133:   PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverType_C",NULL);
    1134:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_elemental_mpidense_C",NULL);
    1135:   PetscFree(A->data);
    1136:   return(0);
    1137: }

    1139: PetscErrorCode MatSetUp_Elemental(Mat A)
    1140: {
    1141:   Mat_Elemental  *a = (Mat_Elemental*)A->data;
    1143:   MPI_Comm       comm;
    1144:   PetscMPIInt    rsize,csize;
    1145:   PetscInt       n;

    1148:   PetscLayoutSetUp(A->rmap);
    1149:   PetscLayoutSetUp(A->cmap);

    1151:   /* Check if local row and column sizes are equally distributed.
    1152:      Jed: Elemental uses "element" cyclic ordering so the sizes need to match that
    1153:      exactly.  The strategy in MatElemental is for PETSc to implicitly permute to block ordering (like would be returned by
    1154:      PetscSplitOwnership(comm,&n,&N), at which point Elemental matrices can act on PETSc vectors without redistributing the vectors. */
    1155:   PetscObjectGetComm((PetscObject)A,&comm);
    1156:   n = PETSC_DECIDE;
    1157:   PetscSplitOwnership(comm,&n,&A->rmap->N);
    1158:   if (n != A->rmap->n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local row size %D of ELEMENTAL matrix must be equally distributed",A->rmap->n);

    1160:   n = PETSC_DECIDE;
    1161:   PetscSplitOwnership(comm,&n,&A->cmap->N);
    1162:   if (n != A->cmap->n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local column size %D of ELEMENTAL matrix must be equally distributed",A->cmap->n);

    1164:   a->emat->Resize(A->rmap->N,A->cmap->N);
    1165:   El::Zero(*a->emat);

    1167:   MPI_Comm_size(A->rmap->comm,&rsize);
    1168:   MPI_Comm_size(A->cmap->comm,&csize);
    1169:   if (csize != rsize) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Cannot use row and column communicators of different sizes");
    1170:   a->commsize = rsize;
    1171:   a->mr[0] = A->rmap->N % rsize; if (!a->mr[0]) a->mr[0] = rsize;
    1172:   a->mr[1] = A->cmap->N % csize; if (!a->mr[1]) a->mr[1] = csize;
    1173:   a->m[0]  = A->rmap->N / rsize + (a->mr[0] != rsize);
    1174:   a->m[1]  = A->cmap->N / csize + (a->mr[1] != csize);
    1175:   return(0);
    1176: }

    1178: PetscErrorCode MatAssemblyBegin_Elemental(Mat A, MatAssemblyType type)
    1179: {
    1180:   Mat_Elemental  *a = (Mat_Elemental*)A->data;

    1183:   /* printf("Calling ProcessQueues\n"); */
    1184:   a->emat->ProcessQueues();
    1185:   /* printf("Finished ProcessQueues\n"); */
    1186:   return(0);
    1187: }

    1189: PetscErrorCode MatAssemblyEnd_Elemental(Mat A, MatAssemblyType type)
    1190: {
    1192:   /* Currently does nothing */
    1193:   return(0);
    1194: }

    1196: PetscErrorCode MatLoad_Elemental(Mat newMat, PetscViewer viewer)
    1197: {
    1199:   Mat            Adense,Ae;
    1200:   MPI_Comm       comm;

    1203:   PetscObjectGetComm((PetscObject)newMat,&comm);
    1204:   MatCreate(comm,&Adense);
    1205:   MatSetType(Adense,MATDENSE);
    1206:   MatLoad(Adense,viewer);
    1207:   MatConvert(Adense, MATELEMENTAL, MAT_INITIAL_MATRIX,&Ae);
    1208:   MatDestroy(&Adense);
    1209:   MatHeaderReplace(newMat,&Ae);
    1210:   return(0);
    1211: }

    1213: /* -------------------------------------------------------------------*/
    1214: static struct _MatOps MatOps_Values = {
    1215:        MatSetValues_Elemental,
    1216:        0,
    1217:        0,
    1218:        MatMult_Elemental,
    1219: /* 4*/ MatMultAdd_Elemental,
    1220:        MatMultTranspose_Elemental,
    1221:        MatMultTransposeAdd_Elemental,
    1222:        MatSolve_Elemental,
    1223:        MatSolveAdd_Elemental,
    1224:        0,
    1225: /*10*/ 0,
    1226:        MatLUFactor_Elemental,
    1227:        MatCholeskyFactor_Elemental,
    1228:        0,
    1229:        MatTranspose_Elemental,
    1230: /*15*/ MatGetInfo_Elemental,
    1231:        0,
    1232:        MatGetDiagonal_Elemental,
    1233:        MatDiagonalScale_Elemental,
    1234:        MatNorm_Elemental,
    1235: /*20*/ MatAssemblyBegin_Elemental,
    1236:        MatAssemblyEnd_Elemental,
    1237:        MatSetOption_Elemental,
    1238:        MatZeroEntries_Elemental,
    1239: /*24*/ 0,
    1240:        MatLUFactorSymbolic_Elemental,
    1241:        MatLUFactorNumeric_Elemental,
    1242:        MatCholeskyFactorSymbolic_Elemental,
    1243:        MatCholeskyFactorNumeric_Elemental,
    1244: /*29*/ MatSetUp_Elemental,
    1245:        0,
    1246:        0,
    1247:        0,
    1248:        0,
    1249: /*34*/ MatDuplicate_Elemental,
    1250:        0,
    1251:        0,
    1252:        0,
    1253:        0,
    1254: /*39*/ MatAXPY_Elemental,
    1255:        0,
    1256:        0,
    1257:        0,
    1258:        MatCopy_Elemental,
    1259: /*44*/ 0,
    1260:        MatScale_Elemental,
    1261:        MatShift_Basic,
    1262:        0,
    1263:        0,
    1264: /*49*/ 0,
    1265:        0,
    1266:        0,
    1267:        0,
    1268:        0,
    1269: /*54*/ 0,
    1270:        0,
    1271:        0,
    1272:        0,
    1273:        0,
    1274: /*59*/ 0,
    1275:        MatDestroy_Elemental,
    1276:        MatView_Elemental,
    1277:        0,
    1278:        0,
    1279: /*64*/ 0,
    1280:        0,
    1281:        0,
    1282:        0,
    1283:        0,
    1284: /*69*/ 0,
    1285:        0,
    1286:        MatConvert_Elemental_Dense,
    1287:        0,
    1288:        0,
    1289: /*74*/ 0,
    1290:        0,
    1291:        0,
    1292:        0,
    1293:        0,
    1294: /*79*/ 0,
    1295:        0,
    1296:        0,
    1297:        0,
    1298:        MatLoad_Elemental,
    1299: /*84*/ 0,
    1300:        0,
    1301:        0,
    1302:        0,
    1303:        0,
    1304: /*89*/ 0,
    1305:        0,
    1306:        MatMatMultNumeric_Elemental,
    1307:        0,
    1308:        0,
    1309: /*94*/ 0,
    1310:        0,
    1311:        0,
    1312:        MatMatTransposeMultNumeric_Elemental,
    1313:        0,
    1314: /*99*/ MatProductSetFromOptions_Elemental,
    1315:        0,
    1316:        0,
    1317:        MatConjugate_Elemental,
    1318:        0,
    1319: /*104*/0,
    1320:        0,
    1321:        0,
    1322:        0,
    1323:        0,
    1324: /*109*/MatMatSolve_Elemental,
    1325:        0,
    1326:        0,
    1327:        0,
    1328:        MatMissingDiagonal_Elemental,
    1329: /*114*/0,
    1330:        0,
    1331:        0,
    1332:        0,
    1333:        0,
    1334: /*119*/0,
    1335:        MatHermitianTranspose_Elemental,
    1336:        0,
    1337:        0,
    1338:        0,
    1339: /*124*/0,
    1340:        0,
    1341:        0,
    1342:        0,
    1343:        0,
    1344: /*129*/0,
    1345:        0,
    1346:        0,
    1347:        0,
    1348:        0,
    1349: /*134*/0,
    1350:        0,
    1351:        0,
    1352:        0,
    1353:        0,
    1354:        0,
    1355: /*140*/0,
    1356:        0,
    1357:        0,
    1358:        0,
    1359:        0,
    1360: /*145*/0,
    1361:        0,
    1362:        0
    1363: };

    1365: /*MC
    1366:    MATELEMENTAL = "elemental" - A matrix type for dense matrices using the Elemental package

    1368:   Use ./configure --download-elemental to install PETSc to use Elemental

    1370:   Use -pc_type lu -pc_factor_mat_solver_type elemental to use this direct solver

    1372:    Options Database Keys:
    1373: + -mat_type elemental - sets the matrix type to "elemental" during a call to MatSetFromOptions()
    1374: - -mat_elemental_grid_height - sets Grid Height for 2D cyclic ordering of internal matrix

    1376:   Level: beginner

    1378: .seealso: MATDENSE
    1379: M*/

    1381: PETSC_EXTERN PetscErrorCode MatCreate_Elemental(Mat A)
    1382: {
    1383:   Mat_Elemental      *a;
    1384:   PetscErrorCode     ierr;
    1385:   PetscBool          flg,flg1;
    1386:   Mat_Elemental_Grid *commgrid;
    1387:   MPI_Comm           icomm;
    1388:   PetscInt           optv1;

    1391:   PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));
    1392:   A->insertmode = NOT_SET_VALUES;

    1394:   PetscNewLog(A,&a);
    1395:   A->data = (void*)a;

    1397:   /* Set up the elemental matrix */
    1398:   El::mpi::Comm cxxcomm(PetscObjectComm((PetscObject)A));

    1400:   /* Grid needs to be shared between multiple Mats on the same communicator, implement by attribute caching on the MPI_Comm */
    1401:   if (Petsc_Elemental_keyval == MPI_KEYVAL_INVALID) {
    1402:     MPI_Comm_create_keyval(MPI_COMM_NULL_COPY_FN,MPI_COMM_NULL_DELETE_FN,&Petsc_Elemental_keyval,(void*)0);
    1403:     PetscCitationsRegister(ElementalCitation,&ElementalCite);
    1404:   }
    1405:   PetscCommDuplicate(cxxcomm.comm,&icomm,NULL);
    1406:   MPI_Comm_get_attr(icomm,Petsc_Elemental_keyval,(void**)&commgrid,(int*)&flg);
    1407:   if (!flg) {
    1408:     PetscNewLog(A,&commgrid);

    1410:     PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"Elemental Options","Mat");
    1411:     /* displayed default grid sizes (CommSize,1) are set by us arbitrarily until El::Grid() is called */
    1412:     PetscOptionsInt("-mat_elemental_grid_height","Grid Height","None",El::mpi::Size(cxxcomm),&optv1,&flg1);
    1413:     if (flg1) {
    1414:       if (El::mpi::Size(cxxcomm) % optv1) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Grid Height %D must evenly divide CommSize %D",optv1,(PetscInt)El::mpi::Size(cxxcomm));
    1415:       commgrid->grid = new El::Grid(cxxcomm,optv1); /* use user-provided grid height */
    1416:     } else {
    1417:       commgrid->grid = new El::Grid(cxxcomm); /* use Elemental default grid sizes */
    1418:       /* printf("new commgrid->grid = %p\n",commgrid->grid);  -- memory leak revealed by valgrind? */
    1419:     }
    1420:     commgrid->grid_refct = 1;
    1421:     MPI_Comm_set_attr(icomm,Petsc_Elemental_keyval,(void*)commgrid);

    1423:     a->pivoting    = 1;
    1424:     PetscOptionsInt("-mat_elemental_pivoting","Pivoting","None",a->pivoting,&a->pivoting,NULL);

    1426:     PetscOptionsEnd();
    1427:   } else {
    1428:     commgrid->grid_refct++;
    1429:   }
    1430:   PetscCommDestroy(&icomm);
    1431:   a->grid        = commgrid->grid;
    1432:   a->emat        = new El::DistMatrix<PetscElemScalar>(*a->grid);
    1433:   a->roworiented = PETSC_TRUE;

    1435:   PetscObjectComposeFunction((PetscObject)A,"MatGetOwnershipIS_C",MatGetOwnershipIS_Elemental);
    1436:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_elemental_mpidense_C",MatProductSetFromOptions_Elemental_MPIDense);
    1437:   PetscObjectChangeTypeName((PetscObject)A,MATELEMENTAL);
    1438:   return(0);
    1439: }