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: {
     21:   Mat_Elemental *a = (Mat_Elemental *)A->data;
     22:   PetscBool      iascii;

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

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

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

     64:   PetscFunctionBegin;
     65:   info->block_size = 1.0;

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

     82:   info->nz_unneeded       = 0.0;
     83:   info->assemblies        = A->num_ass;
     84:   info->mallocs           = 0;
     85:   info->memory            = 0; /* REVIEW ME */
     86:   info->fill_ratio_given  = 0; /* determined by Elemental */
     87:   info->fill_ratio_needed = 0;
     88:   info->factor_mallocs    = 0;
     89:   PetscFunctionReturn(PETSC_SUCCESS);
     90: }

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

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

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

    119:   PetscFunctionBegin;
    120:   // TODO: Initialize matrix to all zeros?

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

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

    201:     /* printf("numQueues=%d\n",numQueues); */
    202:     a->emat->Reserve(numQueues);
    203:     for (j = 0; j < nc; j++) {
    204:       if (cols[j] < 0) continue;
    205:       P2RO(A, 1, cols[j], &crank, &cidx);
    206:       RO2E(A, 1, crank, cidx, &ecol);

    208:       for (i = 0; i < nr; i++) {
    209:         if (rows[i] < 0) continue;
    210:         P2RO(A, 0, rows[i], &rrank, &ridx);
    211:         RO2E(A, 0, rrank, ridx, &erow);
    212:         if (!a->emat->IsLocal(erow, ecol)) { /*off-proc entry*/
    213:           /* printf("Queueing remotely update of (%d,%d)\n",erow,ecol); */
    214:           a->emat->QueueUpdate(erow, ecol, vals[i + j * nr]);
    215:         }
    216:       }
    217:     }
    218:   }
    219:   PetscFunctionReturn(PETSC_SUCCESS);
    220: }

    222: static PetscErrorCode MatMult_Elemental(Mat A, Vec X, Vec Y)
    223: {
    224:   Mat_Elemental         *a = (Mat_Elemental *)A->data;
    225:   const PetscElemScalar *x;
    226:   PetscElemScalar       *y;
    227:   PetscElemScalar        one = 1, zero = 0;

    229:   PetscFunctionBegin;
    230:   PetscCall(VecGetArrayRead(X, (const PetscScalar **)&x));
    231:   PetscCall(VecGetArray(Y, (PetscScalar **)&y));
    232:   { /* Scoping so that constructor is called before pointer is returned */
    233:     El::DistMatrix<PetscElemScalar, El::VC, El::STAR> xe, ye;
    234:     xe.LockedAttach(A->cmap->N, 1, *a->grid, 0, 0, x, A->cmap->n);
    235:     ye.Attach(A->rmap->N, 1, *a->grid, 0, 0, y, A->rmap->n);
    236:     El::Gemv(El::NORMAL, one, *a->emat, xe, zero, ye);
    237:   }
    238:   PetscCall(VecRestoreArrayRead(X, (const PetscScalar **)&x));
    239:   PetscCall(VecRestoreArray(Y, (PetscScalar **)&y));
    240:   PetscFunctionReturn(PETSC_SUCCESS);
    241: }

    243: static PetscErrorCode MatMultTranspose_Elemental(Mat A, Vec X, Vec Y)
    244: {
    245:   Mat_Elemental         *a = (Mat_Elemental *)A->data;
    246:   const PetscElemScalar *x;
    247:   PetscElemScalar       *y;
    248:   PetscElemScalar        one = 1, zero = 0;

    250:   PetscFunctionBegin;
    251:   PetscCall(VecGetArrayRead(X, (const PetscScalar **)&x));
    252:   PetscCall(VecGetArray(Y, (PetscScalar **)&y));
    253:   { /* Scoping so that constructor is called before pointer is returned */
    254:     El::DistMatrix<PetscElemScalar, El::VC, El::STAR> xe, ye;
    255:     xe.LockedAttach(A->rmap->N, 1, *a->grid, 0, 0, x, A->rmap->n);
    256:     ye.Attach(A->cmap->N, 1, *a->grid, 0, 0, y, A->cmap->n);
    257:     El::Gemv(El::TRANSPOSE, one, *a->emat, xe, zero, ye);
    258:   }
    259:   PetscCall(VecRestoreArrayRead(X, (const PetscScalar **)&x));
    260:   PetscCall(VecRestoreArray(Y, (PetscScalar **)&y));
    261:   PetscFunctionReturn(PETSC_SUCCESS);
    262: }

    264: static PetscErrorCode MatMultAdd_Elemental(Mat A, Vec X, Vec Y, Vec Z)
    265: {
    266:   Mat_Elemental         *a = (Mat_Elemental *)A->data;
    267:   const PetscElemScalar *x;
    268:   PetscElemScalar       *z;
    269:   PetscElemScalar        one = 1;

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

    286: static PetscErrorCode MatMultTransposeAdd_Elemental(Mat A, Vec X, Vec Y, Vec Z)
    287: {
    288:   Mat_Elemental         *a = (Mat_Elemental *)A->data;
    289:   const PetscElemScalar *x;
    290:   PetscElemScalar       *z;
    291:   PetscElemScalar        one = 1;

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

    308: PetscErrorCode MatMatMultNumeric_Elemental(Mat A, Mat B, Mat C)
    309: {
    310:   Mat_Elemental  *a   = (Mat_Elemental *)A->data;
    311:   Mat_Elemental  *b   = (Mat_Elemental *)B->data;
    312:   Mat_Elemental  *c   = (Mat_Elemental *)C->data;
    313:   PetscElemScalar one = 1, zero = 0;

    315:   PetscFunctionBegin;
    316:   { /* Scoping so that constructor is called before pointer is returned */
    317:     El::Gemm(El::NORMAL, El::NORMAL, one, *a->emat, *b->emat, zero, *c->emat);
    318:   }
    319:   C->assembled = PETSC_TRUE;
    320:   PetscFunctionReturn(PETSC_SUCCESS);
    321: }

    323: PetscErrorCode MatMatMultSymbolic_Elemental(Mat A, Mat B, PetscReal, Mat Ce)
    324: {
    325:   PetscFunctionBegin;
    326:   PetscCall(MatSetSizes(Ce, A->rmap->n, B->cmap->n, PETSC_DECIDE, PETSC_DECIDE));
    327:   PetscCall(MatSetType(Ce, MATELEMENTAL));
    328:   PetscCall(MatSetUp(Ce));
    329:   Ce->ops->matmultnumeric = MatMatMultNumeric_Elemental;
    330:   PetscFunctionReturn(PETSC_SUCCESS);
    331: }

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

    340:   PetscFunctionBegin;
    341:   { /* Scoping so that constructor is called before pointer is returned */
    342:     El::Gemm(El::NORMAL, El::TRANSPOSE, one, *a->emat, *b->emat, zero, *c->emat);
    343:   }
    344:   C->assembled = PETSC_TRUE;
    345:   PetscFunctionReturn(PETSC_SUCCESS);
    346: }

    348: static PetscErrorCode MatMatTransposeMultSymbolic_Elemental(Mat A, Mat B, PetscReal, Mat C)
    349: {
    350:   PetscFunctionBegin;
    351:   PetscCall(MatSetSizes(C, A->rmap->n, B->rmap->n, PETSC_DECIDE, PETSC_DECIDE));
    352:   PetscCall(MatSetType(C, MATELEMENTAL));
    353:   PetscCall(MatSetUp(C));
    354:   PetscFunctionReturn(PETSC_SUCCESS);
    355: }

    357: static PetscErrorCode MatProductSetFromOptions_Elemental_AB(Mat C)
    358: {
    359:   PetscFunctionBegin;
    360:   C->ops->matmultsymbolic = MatMatMultSymbolic_Elemental;
    361:   C->ops->productsymbolic = MatProductSymbolic_AB;
    362:   PetscFunctionReturn(PETSC_SUCCESS);
    363: }

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

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

    377:   PetscFunctionBegin;
    378:   switch (product->type) {
    379:   case MATPRODUCT_AB:
    380:     PetscCall(MatProductSetFromOptions_Elemental_AB(C));
    381:     break;
    382:   case MATPRODUCT_ABt:
    383:     PetscCall(MatProductSetFromOptions_Elemental_ABt(C));
    384:     break;
    385:   default:
    386:     break;
    387:   }
    388:   PetscFunctionReturn(PETSC_SUCCESS);
    389: }

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

    395:   PetscFunctionBegin;
    396:   PetscCall(MatConvert(B, MATELEMENTAL, MAT_INITIAL_MATRIX, &Be));
    397:   PetscCall(MatMatMult(A, Be, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &Ce));
    398:   PetscCall(MatConvert(Ce, MATMPIDENSE, MAT_REUSE_MATRIX, &C));
    399:   PetscCall(MatDestroy(&Be));
    400:   PetscCall(MatDestroy(&Ce));
    401:   PetscFunctionReturn(PETSC_SUCCESS);
    402: }

    404: static PetscErrorCode MatMatMultSymbolic_Elemental_MPIDense(Mat A, Mat B, PetscReal, Mat C)
    405: {
    406:   PetscFunctionBegin;
    407:   PetscCall(MatSetSizes(C, A->rmap->n, B->cmap->n, PETSC_DECIDE, PETSC_DECIDE));
    408:   PetscCall(MatSetType(C, MATMPIDENSE));
    409:   PetscCall(MatSetUp(C));
    410:   C->ops->matmultnumeric = MatMatMultNumeric_Elemental_MPIDense;
    411:   PetscFunctionReturn(PETSC_SUCCESS);
    412: }

    414: static PetscErrorCode MatProductSetFromOptions_Elemental_MPIDense_AB(Mat C)
    415: {
    416:   PetscFunctionBegin;
    417:   C->ops->matmultsymbolic = MatMatMultSymbolic_Elemental_MPIDense;
    418:   C->ops->productsymbolic = MatProductSymbolic_AB;
    419:   PetscFunctionReturn(PETSC_SUCCESS);
    420: }

    422: static PetscErrorCode MatProductSetFromOptions_Elemental_MPIDense(Mat C)
    423: {
    424:   Mat_Product *product = C->product;

    426:   PetscFunctionBegin;
    427:   if (product->type == MATPRODUCT_AB) PetscCall(MatProductSetFromOptions_Elemental_MPIDense_AB(C));
    428:   PetscFunctionReturn(PETSC_SUCCESS);
    429: }

    431: static PetscErrorCode MatGetDiagonal_Elemental(Mat A, Vec D)
    432: {
    433:   PetscInt        i, nrows, ncols, nD, rrank, ridx, crank, cidx;
    434:   Mat_Elemental  *a = (Mat_Elemental *)A->data;
    435:   PetscElemScalar v;
    436:   MPI_Comm        comm;

    438:   PetscFunctionBegin;
    439:   PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
    440:   PetscCall(MatGetSize(A, &nrows, &ncols));
    441:   nD = nrows > ncols ? ncols : nrows;
    442:   for (i = 0; i < nD; i++) {
    443:     PetscInt erow, ecol;
    444:     P2RO(A, 0, i, &rrank, &ridx);
    445:     RO2E(A, 0, rrank, ridx, &erow);
    446:     PetscCheck(rrank >= 0 && ridx >= 0 && erow >= 0, comm, PETSC_ERR_PLIB, "Incorrect row translation");
    447:     P2RO(A, 1, i, &crank, &cidx);
    448:     RO2E(A, 1, crank, cidx, &ecol);
    449:     PetscCheck(crank >= 0 && cidx >= 0 && ecol >= 0, comm, PETSC_ERR_PLIB, "Incorrect col translation");
    450:     v = a->emat->Get(erow, ecol);
    451:     PetscCall(VecSetValues(D, 1, &i, (PetscScalar *)&v, INSERT_VALUES));
    452:   }
    453:   PetscCall(VecAssemblyBegin(D));
    454:   PetscCall(VecAssemblyEnd(D));
    455:   PetscFunctionReturn(PETSC_SUCCESS);
    456: }

    458: static PetscErrorCode MatDiagonalScale_Elemental(Mat X, Vec L, Vec R)
    459: {
    460:   Mat_Elemental         *x = (Mat_Elemental *)X->data;
    461:   const PetscElemScalar *d;

    463:   PetscFunctionBegin;
    464:   if (R) {
    465:     PetscCall(VecGetArrayRead(R, (const PetscScalar **)&d));
    466:     El::DistMatrix<PetscElemScalar, El::VC, El::STAR> de;
    467:     de.LockedAttach(X->cmap->N, 1, *x->grid, 0, 0, d, X->cmap->n);
    468:     El::DiagonalScale(El::RIGHT, El::NORMAL, de, *x->emat);
    469:     PetscCall(VecRestoreArrayRead(R, (const PetscScalar **)&d));
    470:   }
    471:   if (L) {
    472:     PetscCall(VecGetArrayRead(L, (const PetscScalar **)&d));
    473:     El::DistMatrix<PetscElemScalar, El::VC, El::STAR> de;
    474:     de.LockedAttach(X->rmap->N, 1, *x->grid, 0, 0, d, X->rmap->n);
    475:     El::DiagonalScale(El::LEFT, El::NORMAL, de, *x->emat);
    476:     PetscCall(VecRestoreArrayRead(L, (const PetscScalar **)&d));
    477:   }
    478:   PetscFunctionReturn(PETSC_SUCCESS);
    479: }

    481: static PetscErrorCode MatMissingDiagonal_Elemental(Mat, PetscBool *missing, PetscInt *)
    482: {
    483:   PetscFunctionBegin;
    484:   *missing = PETSC_FALSE;
    485:   PetscFunctionReturn(PETSC_SUCCESS);
    486: }

    488: static PetscErrorCode MatScale_Elemental(Mat X, PetscScalar a)
    489: {
    490:   Mat_Elemental *x = (Mat_Elemental *)X->data;

    492:   PetscFunctionBegin;
    493:   El::Scale((PetscElemScalar)a, *x->emat);
    494:   PetscFunctionReturn(PETSC_SUCCESS);
    495: }

    497: /*
    498:   MatAXPY - Computes Y = a*X + Y.
    499: */
    500: static PetscErrorCode MatAXPY_Elemental(Mat Y, PetscScalar a, Mat X, MatStructure)
    501: {
    502:   Mat_Elemental *x = (Mat_Elemental *)X->data;
    503:   Mat_Elemental *y = (Mat_Elemental *)Y->data;

    505:   PetscFunctionBegin;
    506:   El::Axpy((PetscElemScalar)a, *x->emat, *y->emat);
    507:   PetscCall(PetscObjectStateIncrease((PetscObject)Y));
    508:   PetscFunctionReturn(PETSC_SUCCESS);
    509: }

    511: static PetscErrorCode MatCopy_Elemental(Mat A, Mat B, MatStructure)
    512: {
    513:   Mat_Elemental *a = (Mat_Elemental *)A->data;
    514:   Mat_Elemental *b = (Mat_Elemental *)B->data;

    516:   PetscFunctionBegin;
    517:   El::Copy(*a->emat, *b->emat);
    518:   PetscCall(PetscObjectStateIncrease((PetscObject)B));
    519:   PetscFunctionReturn(PETSC_SUCCESS);
    520: }

    522: static PetscErrorCode MatDuplicate_Elemental(Mat A, MatDuplicateOption op, Mat *B)
    523: {
    524:   Mat            Be;
    525:   MPI_Comm       comm;
    526:   Mat_Elemental *a = (Mat_Elemental *)A->data;

    528:   PetscFunctionBegin;
    529:   PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
    530:   PetscCall(MatCreate(comm, &Be));
    531:   PetscCall(MatSetSizes(Be, A->rmap->n, A->cmap->n, PETSC_DECIDE, PETSC_DECIDE));
    532:   PetscCall(MatSetType(Be, MATELEMENTAL));
    533:   PetscCall(MatSetUp(Be));
    534:   *B = Be;
    535:   if (op == MAT_COPY_VALUES) {
    536:     Mat_Elemental *b = (Mat_Elemental *)Be->data;
    537:     El::Copy(*a->emat, *b->emat);
    538:   }
    539:   Be->assembled = PETSC_TRUE;
    540:   PetscFunctionReturn(PETSC_SUCCESS);
    541: }

    543: static PetscErrorCode MatTranspose_Elemental(Mat A, MatReuse reuse, Mat *B)
    544: {
    545:   Mat            Be = *B;
    546:   MPI_Comm       comm;
    547:   Mat_Elemental *a = (Mat_Elemental *)A->data, *b;

    549:   PetscFunctionBegin;
    550:   if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *B));
    551:   PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
    552:   /* Only out-of-place supported */
    553:   PetscCheck(reuse != MAT_INPLACE_MATRIX, comm, PETSC_ERR_SUP, "Only out-of-place supported");
    554:   if (reuse == MAT_INITIAL_MATRIX) {
    555:     PetscCall(MatCreate(comm, &Be));
    556:     PetscCall(MatSetSizes(Be, A->cmap->n, A->rmap->n, PETSC_DECIDE, PETSC_DECIDE));
    557:     PetscCall(MatSetType(Be, MATELEMENTAL));
    558:     PetscCall(MatSetUp(Be));
    559:     *B = Be;
    560:   }
    561:   b = (Mat_Elemental *)Be->data;
    562:   El::Transpose(*a->emat, *b->emat);
    563:   Be->assembled = PETSC_TRUE;
    564:   PetscFunctionReturn(PETSC_SUCCESS);
    565: }

    567: static PetscErrorCode MatConjugate_Elemental(Mat A)
    568: {
    569:   Mat_Elemental *a = (Mat_Elemental *)A->data;

    571:   PetscFunctionBegin;
    572:   El::Conjugate(*a->emat);
    573:   PetscFunctionReturn(PETSC_SUCCESS);
    574: }

    576: static PetscErrorCode MatHermitianTranspose_Elemental(Mat A, MatReuse reuse, Mat *B)
    577: {
    578:   Mat            Be = *B;
    579:   MPI_Comm       comm;
    580:   Mat_Elemental *a = (Mat_Elemental *)A->data, *b;

    582:   PetscFunctionBegin;
    583:   PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
    584:   /* Only out-of-place supported */
    585:   if (reuse == MAT_INITIAL_MATRIX) {
    586:     PetscCall(MatCreate(comm, &Be));
    587:     PetscCall(MatSetSizes(Be, A->cmap->n, A->rmap->n, PETSC_DECIDE, PETSC_DECIDE));
    588:     PetscCall(MatSetType(Be, MATELEMENTAL));
    589:     PetscCall(MatSetUp(Be));
    590:     *B = Be;
    591:   }
    592:   b = (Mat_Elemental *)Be->data;
    593:   El::Adjoint(*a->emat, *b->emat);
    594:   Be->assembled = PETSC_TRUE;
    595:   PetscFunctionReturn(PETSC_SUCCESS);
    596: }

    598: static PetscErrorCode MatSolve_Elemental(Mat A, Vec B, Vec X)
    599: {
    600:   Mat_Elemental   *a = (Mat_Elemental *)A->data;
    601:   PetscElemScalar *x;
    602:   PetscInt         pivoting = a->pivoting;

    604:   PetscFunctionBegin;
    605:   PetscCall(VecCopy(B, X));
    606:   PetscCall(VecGetArray(X, (PetscScalar **)&x));

    608:   El::DistMatrix<PetscElemScalar, El::VC, El::STAR> xe;
    609:   xe.Attach(A->rmap->N, 1, *a->grid, 0, 0, x, A->rmap->n);
    610:   El::DistMatrix<PetscElemScalar, El::MC, El::MR> xer(xe);
    611:   switch (A->factortype) {
    612:   case MAT_FACTOR_LU:
    613:     if (pivoting == 0) {
    614:       El::lu::SolveAfter(El::NORMAL, *a->emat, xer);
    615:     } else if (pivoting == 1) {
    616:       El::lu::SolveAfter(El::NORMAL, *a->emat, *a->P, xer);
    617:     } else { /* pivoting == 2 */
    618:       El::lu::SolveAfter(El::NORMAL, *a->emat, *a->P, *a->Q, xer);
    619:     }
    620:     break;
    621:   case MAT_FACTOR_CHOLESKY:
    622:     El::cholesky::SolveAfter(El::UPPER, El::NORMAL, *a->emat, xer);
    623:     break;
    624:   default:
    625:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unfactored Matrix or Unsupported MatFactorType");
    626:     break;
    627:   }
    628:   El::Copy(xer, xe);

    630:   PetscCall(VecRestoreArray(X, (PetscScalar **)&x));
    631:   PetscFunctionReturn(PETSC_SUCCESS);
    632: }

    634: static PetscErrorCode MatSolveAdd_Elemental(Mat A, Vec B, Vec Y, Vec X)
    635: {
    636:   PetscFunctionBegin;
    637:   PetscCall(MatSolve_Elemental(A, B, X));
    638:   PetscCall(VecAXPY(X, 1, Y));
    639:   PetscFunctionReturn(PETSC_SUCCESS);
    640: }

    642: static PetscErrorCode MatMatSolve_Elemental(Mat A, Mat B, Mat X)
    643: {
    644:   Mat_Elemental *a = (Mat_Elemental *)A->data;
    645:   Mat_Elemental *x;
    646:   Mat            C;
    647:   PetscInt       pivoting = a->pivoting;
    648:   PetscBool      flg;
    649:   MatType        type;

    651:   PetscFunctionBegin;
    652:   PetscCall(MatGetType(X, &type));
    653:   PetscCall(PetscStrcmp(type, MATELEMENTAL, &flg));
    654:   if (!flg) {
    655:     PetscCall(MatConvert(B, MATELEMENTAL, MAT_INITIAL_MATRIX, &C));
    656:     x = (Mat_Elemental *)C->data;
    657:   } else {
    658:     x = (Mat_Elemental *)X->data;
    659:     El::Copy(*((Mat_Elemental *)B->data)->emat, *x->emat);
    660:   }
    661:   switch (A->factortype) {
    662:   case MAT_FACTOR_LU:
    663:     if (pivoting == 0) {
    664:       El::lu::SolveAfter(El::NORMAL, *a->emat, *x->emat);
    665:     } else if (pivoting == 1) {
    666:       El::lu::SolveAfter(El::NORMAL, *a->emat, *a->P, *x->emat);
    667:     } else {
    668:       El::lu::SolveAfter(El::NORMAL, *a->emat, *a->P, *a->Q, *x->emat);
    669:     }
    670:     break;
    671:   case MAT_FACTOR_CHOLESKY:
    672:     El::cholesky::SolveAfter(El::UPPER, El::NORMAL, *a->emat, *x->emat);
    673:     break;
    674:   default:
    675:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unfactored Matrix or Unsupported MatFactorType");
    676:     break;
    677:   }
    678:   if (!flg) {
    679:     PetscCall(MatConvert(C, type, MAT_REUSE_MATRIX, &X));
    680:     PetscCall(MatDestroy(&C));
    681:   }
    682:   PetscFunctionReturn(PETSC_SUCCESS);
    683: }

    685: static PetscErrorCode MatLUFactor_Elemental(Mat A, IS, IS, const MatFactorInfo *)
    686: {
    687:   Mat_Elemental *a        = (Mat_Elemental *)A->data;
    688:   PetscInt       pivoting = a->pivoting;

    690:   PetscFunctionBegin;
    691:   if (pivoting == 0) {
    692:     El::LU(*a->emat);
    693:   } else if (pivoting == 1) {
    694:     El::LU(*a->emat, *a->P);
    695:   } else {
    696:     El::LU(*a->emat, *a->P, *a->Q);
    697:   }
    698:   A->factortype = MAT_FACTOR_LU;
    699:   A->assembled  = PETSC_TRUE;

    701:   PetscCall(PetscFree(A->solvertype));
    702:   PetscCall(PetscStrallocpy(MATSOLVERELEMENTAL, &A->solvertype));
    703:   PetscFunctionReturn(PETSC_SUCCESS);
    704: }

    706: static PetscErrorCode MatLUFactorNumeric_Elemental(Mat F, Mat A, const MatFactorInfo *info)
    707: {
    708:   PetscFunctionBegin;
    709:   PetscCall(MatCopy(A, F, SAME_NONZERO_PATTERN));
    710:   PetscCall(MatLUFactor_Elemental(F, nullptr, nullptr, info));
    711:   PetscFunctionReturn(PETSC_SUCCESS);
    712: }

    714: static PetscErrorCode MatLUFactorSymbolic_Elemental(Mat, Mat, IS, IS, const MatFactorInfo *)
    715: {
    716:   PetscFunctionBegin;
    717:   /* F is created and allocated by MatGetFactor_elemental_petsc(), skip this routine. */
    718:   PetscFunctionReturn(PETSC_SUCCESS);
    719: }

    721: static PetscErrorCode MatCholeskyFactor_Elemental(Mat A, IS, const MatFactorInfo *)
    722: {
    723:   Mat_Elemental                                    *a = (Mat_Elemental *)A->data;
    724:   El::DistMatrix<PetscElemScalar, El::MC, El::STAR> d;

    726:   PetscFunctionBegin;
    727:   El::Cholesky(El::UPPER, *a->emat);
    728:   A->factortype = MAT_FACTOR_CHOLESKY;
    729:   A->assembled  = PETSC_TRUE;

    731:   PetscCall(PetscFree(A->solvertype));
    732:   PetscCall(PetscStrallocpy(MATSOLVERELEMENTAL, &A->solvertype));
    733:   PetscFunctionReturn(PETSC_SUCCESS);
    734: }

    736: static PetscErrorCode MatCholeskyFactorNumeric_Elemental(Mat F, Mat A, const MatFactorInfo *info)
    737: {
    738:   PetscFunctionBegin;
    739:   PetscCall(MatCopy(A, F, SAME_NONZERO_PATTERN));
    740:   PetscCall(MatCholeskyFactor_Elemental(F, nullptr, info));
    741:   PetscFunctionReturn(PETSC_SUCCESS);
    742: }

    744: static PetscErrorCode MatCholeskyFactorSymbolic_Elemental(Mat, Mat, IS, const MatFactorInfo *)
    745: {
    746:   PetscFunctionBegin;
    747:   /* F is created and allocated by MatGetFactor_elemental_petsc(), skip this routine. */
    748:   PetscFunctionReturn(PETSC_SUCCESS);
    749: }

    751: static PetscErrorCode MatFactorGetSolverType_elemental_elemental(Mat, MatSolverType *type)
    752: {
    753:   PetscFunctionBegin;
    754:   *type = MATSOLVERELEMENTAL;
    755:   PetscFunctionReturn(PETSC_SUCCESS);
    756: }

    758: static PetscErrorCode MatGetFactor_elemental_elemental(Mat A, MatFactorType ftype, Mat *F)
    759: {
    760:   Mat B;

    762:   PetscFunctionBegin;
    763:   /* Create the factorization matrix */
    764:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
    765:   PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, PETSC_DECIDE, PETSC_DECIDE));
    766:   PetscCall(MatSetType(B, MATELEMENTAL));
    767:   PetscCall(MatSetUp(B));
    768:   B->factortype      = ftype;
    769:   B->trivialsymbolic = PETSC_TRUE;
    770:   PetscCall(PetscFree(B->solvertype));
    771:   PetscCall(PetscStrallocpy(MATSOLVERELEMENTAL, &B->solvertype));

    773:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_elemental_elemental));
    774:   *F = B;
    775:   PetscFunctionReturn(PETSC_SUCCESS);
    776: }

    778: PETSC_INTERN PetscErrorCode MatSolverTypeRegister_Elemental(void)
    779: {
    780:   PetscFunctionBegin;
    781:   PetscCall(MatSolverTypeRegister(MATSOLVERELEMENTAL, MATELEMENTAL, MAT_FACTOR_LU, MatGetFactor_elemental_elemental));
    782:   PetscCall(MatSolverTypeRegister(MATSOLVERELEMENTAL, MATELEMENTAL, MAT_FACTOR_CHOLESKY, MatGetFactor_elemental_elemental));
    783:   PetscFunctionReturn(PETSC_SUCCESS);
    784: }

    786: static PetscErrorCode MatNorm_Elemental(Mat A, NormType type, PetscReal *nrm)
    787: {
    788:   Mat_Elemental *a = (Mat_Elemental *)A->data;

    790:   PetscFunctionBegin;
    791:   switch (type) {
    792:   case NORM_1:
    793:     *nrm = El::OneNorm(*a->emat);
    794:     break;
    795:   case NORM_FROBENIUS:
    796:     *nrm = El::FrobeniusNorm(*a->emat);
    797:     break;
    798:   case NORM_INFINITY:
    799:     *nrm = El::InfinityNorm(*a->emat);
    800:     break;
    801:   default:
    802:     SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Unsupported norm type");
    803:   }
    804:   PetscFunctionReturn(PETSC_SUCCESS);
    805: }

    807: static PetscErrorCode MatZeroEntries_Elemental(Mat A)
    808: {
    809:   Mat_Elemental *a = (Mat_Elemental *)A->data;

    811:   PetscFunctionBegin;
    812:   El::Zero(*a->emat);
    813:   PetscFunctionReturn(PETSC_SUCCESS);
    814: }

    816: static PetscErrorCode MatGetOwnershipIS_Elemental(Mat A, IS *rows, IS *cols)
    817: {
    818:   Mat_Elemental *a = (Mat_Elemental *)A->data;
    819:   PetscInt       i, m, shift, stride, *idx;

    821:   PetscFunctionBegin;
    822:   if (rows) {
    823:     m      = a->emat->LocalHeight();
    824:     shift  = a->emat->ColShift();
    825:     stride = a->emat->ColStride();
    826:     PetscCall(PetscMalloc1(m, &idx));
    827:     for (i = 0; i < m; i++) {
    828:       PetscInt rank, offset;
    829:       E2RO(A, 0, shift + i * stride, &rank, &offset);
    830:       RO2P(A, 0, rank, offset, &idx[i]);
    831:     }
    832:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, m, idx, PETSC_OWN_POINTER, rows));
    833:   }
    834:   if (cols) {
    835:     m      = a->emat->LocalWidth();
    836:     shift  = a->emat->RowShift();
    837:     stride = a->emat->RowStride();
    838:     PetscCall(PetscMalloc1(m, &idx));
    839:     for (i = 0; i < m; i++) {
    840:       PetscInt rank, offset;
    841:       E2RO(A, 1, shift + i * stride, &rank, &offset);
    842:       RO2P(A, 1, rank, offset, &idx[i]);
    843:     }
    844:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, m, idx, PETSC_OWN_POINTER, cols));
    845:   }
    846:   PetscFunctionReturn(PETSC_SUCCESS);
    847: }

    849: static PetscErrorCode MatConvert_Elemental_Dense(Mat A, MatType, MatReuse reuse, Mat *B)
    850: {
    851:   Mat             Bmpi;
    852:   Mat_Elemental  *a = (Mat_Elemental *)A->data;
    853:   MPI_Comm        comm;
    854:   IS              isrows, iscols;
    855:   PetscInt        rrank, ridx, crank, cidx, nrows, ncols, i, j, erow, ecol, elrow, elcol;
    856:   const PetscInt *rows, *cols;
    857:   PetscElemScalar v;
    858:   const El::Grid &grid = a->emat->Grid();

    860:   PetscFunctionBegin;
    861:   PetscCall(PetscObjectGetComm((PetscObject)A, &comm));

    863:   if (reuse == MAT_REUSE_MATRIX) {
    864:     Bmpi = *B;
    865:   } else {
    866:     PetscCall(MatCreate(comm, &Bmpi));
    867:     PetscCall(MatSetSizes(Bmpi, A->rmap->n, A->cmap->n, PETSC_DECIDE, PETSC_DECIDE));
    868:     PetscCall(MatSetType(Bmpi, MATDENSE));
    869:     PetscCall(MatSetUp(Bmpi));
    870:   }

    872:   /* Get local entries of A */
    873:   PetscCall(MatGetOwnershipIS(A, &isrows, &iscols));
    874:   PetscCall(ISGetLocalSize(isrows, &nrows));
    875:   PetscCall(ISGetIndices(isrows, &rows));
    876:   PetscCall(ISGetLocalSize(iscols, &ncols));
    877:   PetscCall(ISGetIndices(iscols, &cols));

    879:   if (a->roworiented) {
    880:     for (i = 0; i < nrows; i++) {
    881:       P2RO(A, 0, rows[i], &rrank, &ridx); /* convert indices between PETSc <-> (Rank,Offset) <-> Elemental */
    882:       RO2E(A, 0, rrank, ridx, &erow);
    883:       PetscCheck(rrank >= 0 && ridx >= 0 && erow >= 0, comm, PETSC_ERR_PLIB, "Incorrect row translation");
    884:       for (j = 0; j < ncols; j++) {
    885:         P2RO(A, 1, cols[j], &crank, &cidx);
    886:         RO2E(A, 1, crank, cidx, &ecol);
    887:         PetscCheck(crank >= 0 && cidx >= 0 && ecol >= 0, comm, PETSC_ERR_PLIB, "Incorrect col translation");

    889:         elrow = erow / grid.MCSize(); /* Elemental local row index */
    890:         elcol = ecol / grid.MRSize(); /* Elemental local column index */
    891:         v     = a->emat->GetLocal(elrow, elcol);
    892:         PetscCall(MatSetValues(Bmpi, 1, &rows[i], 1, &cols[j], (PetscScalar *)&v, INSERT_VALUES));
    893:       }
    894:     }
    895:   } else { /* column-oriented */
    896:     for (j = 0; j < ncols; j++) {
    897:       P2RO(A, 1, cols[j], &crank, &cidx);
    898:       RO2E(A, 1, crank, cidx, &ecol);
    899:       PetscCheck(crank >= 0 && cidx >= 0 && ecol >= 0, comm, PETSC_ERR_PLIB, "Incorrect col translation");
    900:       for (i = 0; i < nrows; i++) {
    901:         P2RO(A, 0, rows[i], &rrank, &ridx); /* convert indices between PETSc <-> (Rank,Offset) <-> Elemental */
    902:         RO2E(A, 0, rrank, ridx, &erow);
    903:         PetscCheck(rrank >= 0 && ridx >= 0 && erow >= 0, comm, PETSC_ERR_PLIB, "Incorrect row translation");

    905:         elrow = erow / grid.MCSize(); /* Elemental local row index */
    906:         elcol = ecol / grid.MRSize(); /* Elemental local column index */
    907:         v     = a->emat->GetLocal(elrow, elcol);
    908:         PetscCall(MatSetValues(Bmpi, 1, &rows[i], 1, &cols[j], (PetscScalar *)&v, INSERT_VALUES));
    909:       }
    910:     }
    911:   }
    912:   PetscCall(MatAssemblyBegin(Bmpi, MAT_FINAL_ASSEMBLY));
    913:   PetscCall(MatAssemblyEnd(Bmpi, MAT_FINAL_ASSEMBLY));
    914:   if (reuse == MAT_INPLACE_MATRIX) {
    915:     PetscCall(MatHeaderReplace(A, &Bmpi));
    916:   } else {
    917:     *B = Bmpi;
    918:   }
    919:   PetscCall(ISDestroy(&isrows));
    920:   PetscCall(ISDestroy(&iscols));
    921:   PetscFunctionReturn(PETSC_SUCCESS);
    922: }

    924: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat A, MatType, MatReuse reuse, Mat *newmat)
    925: {
    926:   Mat                mat_elemental;
    927:   PetscInt           M = A->rmap->N, N = A->cmap->N, row, ncols;
    928:   const PetscInt    *cols;
    929:   const PetscScalar *vals;

    931:   PetscFunctionBegin;
    932:   if (reuse == MAT_REUSE_MATRIX) {
    933:     mat_elemental = *newmat;
    934:     PetscCall(MatZeroEntries(mat_elemental));
    935:   } else {
    936:     PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental));
    937:     PetscCall(MatSetSizes(mat_elemental, PETSC_DECIDE, PETSC_DECIDE, M, N));
    938:     PetscCall(MatSetType(mat_elemental, MATELEMENTAL));
    939:     PetscCall(MatSetUp(mat_elemental));
    940:   }
    941:   for (row = 0; row < M; row++) {
    942:     PetscCall(MatGetRow(A, row, &ncols, &cols, &vals));
    943:     /* PETSc-Elemental interface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
    944:     PetscCall(MatSetValues(mat_elemental, 1, &row, ncols, cols, vals, ADD_VALUES));
    945:     PetscCall(MatRestoreRow(A, row, &ncols, &cols, &vals));
    946:   }
    947:   PetscCall(MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY));
    948:   PetscCall(MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY));

    950:   if (reuse == MAT_INPLACE_MATRIX) {
    951:     PetscCall(MatHeaderReplace(A, &mat_elemental));
    952:   } else {
    953:     *newmat = mat_elemental;
    954:   }
    955:   PetscFunctionReturn(PETSC_SUCCESS);
    956: }

    958: PETSC_INTERN PetscErrorCode MatConvert_MPIAIJ_Elemental(Mat A, MatType, MatReuse reuse, Mat *newmat)
    959: {
    960:   Mat                mat_elemental;
    961:   PetscInt           row, ncols, rstart = A->rmap->rstart, rend = A->rmap->rend, j;
    962:   const PetscInt    *cols;
    963:   const PetscScalar *vals;

    965:   PetscFunctionBegin;
    966:   if (reuse == MAT_REUSE_MATRIX) {
    967:     mat_elemental = *newmat;
    968:     PetscCall(MatZeroEntries(mat_elemental));
    969:   } else {
    970:     PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental));
    971:     PetscCall(MatSetSizes(mat_elemental, PETSC_DECIDE, PETSC_DECIDE, A->rmap->N, A->cmap->N));
    972:     PetscCall(MatSetType(mat_elemental, MATELEMENTAL));
    973:     PetscCall(MatSetUp(mat_elemental));
    974:   }
    975:   for (row = rstart; row < rend; row++) {
    976:     PetscCall(MatGetRow(A, row, &ncols, &cols, &vals));
    977:     for (j = 0; j < ncols; j++) {
    978:       /* PETSc-Elemental interface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
    979:       PetscCall(MatSetValues(mat_elemental, 1, &row, 1, &cols[j], &vals[j], ADD_VALUES));
    980:     }
    981:     PetscCall(MatRestoreRow(A, row, &ncols, &cols, &vals));
    982:   }
    983:   PetscCall(MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY));
    984:   PetscCall(MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY));

    986:   if (reuse == MAT_INPLACE_MATRIX) {
    987:     PetscCall(MatHeaderReplace(A, &mat_elemental));
    988:   } else {
    989:     *newmat = mat_elemental;
    990:   }
    991:   PetscFunctionReturn(PETSC_SUCCESS);
    992: }

    994: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_Elemental(Mat A, MatType, MatReuse reuse, Mat *newmat)
    995: {
    996:   Mat                mat_elemental;
    997:   PetscInt           M = A->rmap->N, N = A->cmap->N, row, ncols, j;
    998:   const PetscInt    *cols;
    999:   const PetscScalar *vals;

    1001:   PetscFunctionBegin;
    1002:   if (reuse == MAT_REUSE_MATRIX) {
    1003:     mat_elemental = *newmat;
    1004:     PetscCall(MatZeroEntries(mat_elemental));
    1005:   } else {
    1006:     PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental));
    1007:     PetscCall(MatSetSizes(mat_elemental, PETSC_DECIDE, PETSC_DECIDE, M, N));
    1008:     PetscCall(MatSetType(mat_elemental, MATELEMENTAL));
    1009:     PetscCall(MatSetUp(mat_elemental));
    1010:   }
    1011:   PetscCall(MatGetRowUpperTriangular(A));
    1012:   for (row = 0; row < M; row++) {
    1013:     PetscCall(MatGetRow(A, row, &ncols, &cols, &vals));
    1014:     /* PETSc-Elemental interface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
    1015:     PetscCall(MatSetValues(mat_elemental, 1, &row, ncols, cols, vals, ADD_VALUES));
    1016:     for (j = 0; j < ncols; j++) { /* lower triangular part */
    1017:       PetscScalar v;
    1018:       if (cols[j] == row) continue;
    1019:       v = A->hermitian == PETSC_BOOL3_TRUE ? PetscConj(vals[j]) : vals[j];
    1020:       PetscCall(MatSetValues(mat_elemental, 1, &cols[j], 1, &row, &v, ADD_VALUES));
    1021:     }
    1022:     PetscCall(MatRestoreRow(A, row, &ncols, &cols, &vals));
    1023:   }
    1024:   PetscCall(MatRestoreRowUpperTriangular(A));
    1025:   PetscCall(MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY));
    1026:   PetscCall(MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY));

    1028:   if (reuse == MAT_INPLACE_MATRIX) {
    1029:     PetscCall(MatHeaderReplace(A, &mat_elemental));
    1030:   } else {
    1031:     *newmat = mat_elemental;
    1032:   }
    1033:   PetscFunctionReturn(PETSC_SUCCESS);
    1034: }

    1036: PETSC_INTERN PetscErrorCode MatConvert_MPISBAIJ_Elemental(Mat A, MatType, MatReuse reuse, Mat *newmat)
    1037: {
    1038:   Mat                mat_elemental;
    1039:   PetscInt           M = A->rmap->N, N = A->cmap->N, row, ncols, j, rstart = A->rmap->rstart, rend = A->rmap->rend;
    1040:   const PetscInt    *cols;
    1041:   const PetscScalar *vals;

    1043:   PetscFunctionBegin;
    1044:   if (reuse == MAT_REUSE_MATRIX) {
    1045:     mat_elemental = *newmat;
    1046:     PetscCall(MatZeroEntries(mat_elemental));
    1047:   } else {
    1048:     PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental));
    1049:     PetscCall(MatSetSizes(mat_elemental, PETSC_DECIDE, PETSC_DECIDE, M, N));
    1050:     PetscCall(MatSetType(mat_elemental, MATELEMENTAL));
    1051:     PetscCall(MatSetUp(mat_elemental));
    1052:   }
    1053:   PetscCall(MatGetRowUpperTriangular(A));
    1054:   for (row = rstart; row < rend; row++) {
    1055:     PetscCall(MatGetRow(A, row, &ncols, &cols, &vals));
    1056:     /* PETSc-Elemental interface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
    1057:     PetscCall(MatSetValues(mat_elemental, 1, &row, ncols, cols, vals, ADD_VALUES));
    1058:     for (j = 0; j < ncols; j++) { /* lower triangular part */
    1059:       PetscScalar v;
    1060:       if (cols[j] == row) continue;
    1061:       v = A->hermitian == PETSC_BOOL3_TRUE ? PetscConj(vals[j]) : vals[j];
    1062:       PetscCall(MatSetValues(mat_elemental, 1, &cols[j], 1, &row, &v, ADD_VALUES));
    1063:     }
    1064:     PetscCall(MatRestoreRow(A, row, &ncols, &cols, &vals));
    1065:   }
    1066:   PetscCall(MatRestoreRowUpperTriangular(A));
    1067:   PetscCall(MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY));
    1068:   PetscCall(MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY));

    1070:   if (reuse == MAT_INPLACE_MATRIX) {
    1071:     PetscCall(MatHeaderReplace(A, &mat_elemental));
    1072:   } else {
    1073:     *newmat = mat_elemental;
    1074:   }
    1075:   PetscFunctionReturn(PETSC_SUCCESS);
    1076: }

    1078: static PetscErrorCode MatDestroy_Elemental(Mat A)
    1079: {
    1080:   Mat_Elemental      *a = (Mat_Elemental *)A->data;
    1081:   Mat_Elemental_Grid *commgrid;
    1082:   PetscBool           flg;
    1083:   MPI_Comm            icomm;

    1085:   PetscFunctionBegin;
    1086:   delete a->emat;
    1087:   delete a->P;
    1088:   delete a->Q;

    1090:   El::mpi::Comm cxxcomm(PetscObjectComm((PetscObject)A));
    1091:   PetscCall(PetscCommDuplicate(cxxcomm.comm, &icomm, nullptr));
    1092:   PetscCallMPI(MPI_Comm_get_attr(icomm, Petsc_Elemental_keyval, (void **)&commgrid, (int *)&flg));
    1093:   if (--commgrid->grid_refct == 0) {
    1094:     delete commgrid->grid;
    1095:     PetscCall(PetscFree(commgrid));
    1096:     PetscCallMPI(MPI_Comm_free_keyval(&Petsc_Elemental_keyval));
    1097:   }
    1098:   PetscCall(PetscCommDestroy(&icomm));
    1099:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatGetOwnershipIS_C", nullptr));
    1100:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", nullptr));
    1101:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_elemental_mpidense_C", nullptr));
    1102:   PetscCall(PetscFree(A->data));
    1103:   PetscFunctionReturn(PETSC_SUCCESS);
    1104: }

    1106: static PetscErrorCode MatSetUp_Elemental(Mat A)
    1107: {
    1108:   Mat_Elemental *a = (Mat_Elemental *)A->data;
    1109:   MPI_Comm       comm;
    1110:   PetscMPIInt    rsize, csize;
    1111:   PetscInt       n;

    1113:   PetscFunctionBegin;
    1114:   PetscCall(PetscLayoutSetUp(A->rmap));
    1115:   PetscCall(PetscLayoutSetUp(A->cmap));

    1117:   /* Check if local row and column sizes are equally distributed.
    1118:      Jed: Elemental uses "element" cyclic ordering so the sizes need to match that
    1119:      exactly.  The strategy in MatElemental is for PETSc to implicitly permute to block ordering (like would be returned by
    1120:      PetscSplitOwnership(comm,&n,&N), at which point Elemental matrices can act on PETSc vectors without redistributing the vectors. */
    1121:   PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
    1122:   n = PETSC_DECIDE;
    1123:   PetscCall(PetscSplitOwnership(comm, &n, &A->rmap->N));
    1124:   PetscCheck(n == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local row size %" PetscInt_FMT " of ELEMENTAL matrix must be equally distributed", A->rmap->n);

    1126:   n = PETSC_DECIDE;
    1127:   PetscCall(PetscSplitOwnership(comm, &n, &A->cmap->N));
    1128:   PetscCheck(n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local column size %" PetscInt_FMT " of ELEMENTAL matrix must be equally distributed", A->cmap->n);

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

    1133:   PetscCallMPI(MPI_Comm_size(A->rmap->comm, &rsize));
    1134:   PetscCallMPI(MPI_Comm_size(A->cmap->comm, &csize));
    1135:   PetscCheck(csize == rsize, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Cannot use row and column communicators of different sizes");
    1136:   a->commsize = rsize;
    1137:   a->mr[0]    = A->rmap->N % rsize;
    1138:   if (!a->mr[0]) a->mr[0] = rsize;
    1139:   a->mr[1] = A->cmap->N % csize;
    1140:   if (!a->mr[1]) a->mr[1] = csize;
    1141:   a->m[0] = A->rmap->N / rsize + (a->mr[0] != rsize);
    1142:   a->m[1] = A->cmap->N / csize + (a->mr[1] != csize);
    1143:   PetscFunctionReturn(PETSC_SUCCESS);
    1144: }

    1146: static PetscErrorCode MatAssemblyBegin_Elemental(Mat A, MatAssemblyType)
    1147: {
    1148:   Mat_Elemental *a = (Mat_Elemental *)A->data;

    1150:   PetscFunctionBegin;
    1151:   /* printf("Calling ProcessQueues\n"); */
    1152:   a->emat->ProcessQueues();
    1153:   /* printf("Finished ProcessQueues\n"); */
    1154:   PetscFunctionReturn(PETSC_SUCCESS);
    1155: }

    1157: static PetscErrorCode MatAssemblyEnd_Elemental(Mat, MatAssemblyType)
    1158: {
    1159:   PetscFunctionBegin;
    1160:   /* Currently does nothing */
    1161:   PetscFunctionReturn(PETSC_SUCCESS);
    1162: }

    1164: static PetscErrorCode MatLoad_Elemental(Mat newMat, PetscViewer viewer)
    1165: {
    1166:   Mat      Adense, Ae;
    1167:   MPI_Comm comm;

    1169:   PetscFunctionBegin;
    1170:   PetscCall(PetscObjectGetComm((PetscObject)newMat, &comm));
    1171:   PetscCall(MatCreate(comm, &Adense));
    1172:   PetscCall(MatSetType(Adense, MATDENSE));
    1173:   PetscCall(MatLoad(Adense, viewer));
    1174:   PetscCall(MatConvert(Adense, MATELEMENTAL, MAT_INITIAL_MATRIX, &Ae));
    1175:   PetscCall(MatDestroy(&Adense));
    1176:   PetscCall(MatHeaderReplace(newMat, &Ae));
    1177:   PetscFunctionReturn(PETSC_SUCCESS);
    1178: }

    1180: static struct _MatOps MatOps_Values = {MatSetValues_Elemental,
    1181:                                        nullptr,
    1182:                                        nullptr,
    1183:                                        MatMult_Elemental,
    1184:                                        /* 4*/ MatMultAdd_Elemental,
    1185:                                        MatMultTranspose_Elemental,
    1186:                                        MatMultTransposeAdd_Elemental,
    1187:                                        MatSolve_Elemental,
    1188:                                        MatSolveAdd_Elemental,
    1189:                                        nullptr,
    1190:                                        /*10*/ nullptr,
    1191:                                        MatLUFactor_Elemental,
    1192:                                        MatCholeskyFactor_Elemental,
    1193:                                        nullptr,
    1194:                                        MatTranspose_Elemental,
    1195:                                        /*15*/ MatGetInfo_Elemental,
    1196:                                        nullptr,
    1197:                                        MatGetDiagonal_Elemental,
    1198:                                        MatDiagonalScale_Elemental,
    1199:                                        MatNorm_Elemental,
    1200:                                        /*20*/ MatAssemblyBegin_Elemental,
    1201:                                        MatAssemblyEnd_Elemental,
    1202:                                        MatSetOption_Elemental,
    1203:                                        MatZeroEntries_Elemental,
    1204:                                        /*24*/ nullptr,
    1205:                                        MatLUFactorSymbolic_Elemental,
    1206:                                        MatLUFactorNumeric_Elemental,
    1207:                                        MatCholeskyFactorSymbolic_Elemental,
    1208:                                        MatCholeskyFactorNumeric_Elemental,
    1209:                                        /*29*/ MatSetUp_Elemental,
    1210:                                        nullptr,
    1211:                                        nullptr,
    1212:                                        nullptr,
    1213:                                        nullptr,
    1214:                                        /*34*/ MatDuplicate_Elemental,
    1215:                                        nullptr,
    1216:                                        nullptr,
    1217:                                        nullptr,
    1218:                                        nullptr,
    1219:                                        /*39*/ MatAXPY_Elemental,
    1220:                                        nullptr,
    1221:                                        nullptr,
    1222:                                        nullptr,
    1223:                                        MatCopy_Elemental,
    1224:                                        /*44*/ nullptr,
    1225:                                        MatScale_Elemental,
    1226:                                        MatShift_Basic,
    1227:                                        nullptr,
    1228:                                        nullptr,
    1229:                                        /*49*/ nullptr,
    1230:                                        nullptr,
    1231:                                        nullptr,
    1232:                                        nullptr,
    1233:                                        nullptr,
    1234:                                        /*54*/ nullptr,
    1235:                                        nullptr,
    1236:                                        nullptr,
    1237:                                        nullptr,
    1238:                                        nullptr,
    1239:                                        /*59*/ nullptr,
    1240:                                        MatDestroy_Elemental,
    1241:                                        MatView_Elemental,
    1242:                                        nullptr,
    1243:                                        nullptr,
    1244:                                        /*64*/ nullptr,
    1245:                                        nullptr,
    1246:                                        nullptr,
    1247:                                        nullptr,
    1248:                                        nullptr,
    1249:                                        /*69*/ nullptr,
    1250:                                        nullptr,
    1251:                                        MatConvert_Elemental_Dense,
    1252:                                        nullptr,
    1253:                                        nullptr,
    1254:                                        /*74*/ nullptr,
    1255:                                        nullptr,
    1256:                                        nullptr,
    1257:                                        nullptr,
    1258:                                        nullptr,
    1259:                                        /*79*/ nullptr,
    1260:                                        nullptr,
    1261:                                        nullptr,
    1262:                                        nullptr,
    1263:                                        MatLoad_Elemental,
    1264:                                        /*84*/ nullptr,
    1265:                                        nullptr,
    1266:                                        nullptr,
    1267:                                        nullptr,
    1268:                                        nullptr,
    1269:                                        /*89*/ nullptr,
    1270:                                        nullptr,
    1271:                                        MatMatMultNumeric_Elemental,
    1272:                                        nullptr,
    1273:                                        nullptr,
    1274:                                        /*94*/ nullptr,
    1275:                                        nullptr,
    1276:                                        nullptr,
    1277:                                        MatMatTransposeMultNumeric_Elemental,
    1278:                                        nullptr,
    1279:                                        /*99*/ MatProductSetFromOptions_Elemental,
    1280:                                        nullptr,
    1281:                                        nullptr,
    1282:                                        MatConjugate_Elemental,
    1283:                                        nullptr,
    1284:                                        /*104*/ nullptr,
    1285:                                        nullptr,
    1286:                                        nullptr,
    1287:                                        nullptr,
    1288:                                        nullptr,
    1289:                                        /*109*/ MatMatSolve_Elemental,
    1290:                                        nullptr,
    1291:                                        nullptr,
    1292:                                        nullptr,
    1293:                                        MatMissingDiagonal_Elemental,
    1294:                                        /*114*/ nullptr,
    1295:                                        nullptr,
    1296:                                        nullptr,
    1297:                                        nullptr,
    1298:                                        nullptr,
    1299:                                        /*119*/ nullptr,
    1300:                                        MatHermitianTranspose_Elemental,
    1301:                                        nullptr,
    1302:                                        nullptr,
    1303:                                        nullptr,
    1304:                                        /*124*/ nullptr,
    1305:                                        nullptr,
    1306:                                        nullptr,
    1307:                                        nullptr,
    1308:                                        nullptr,
    1309:                                        /*129*/ nullptr,
    1310:                                        nullptr,
    1311:                                        nullptr,
    1312:                                        nullptr,
    1313:                                        nullptr,
    1314:                                        /*134*/ nullptr,
    1315:                                        nullptr,
    1316:                                        nullptr,
    1317:                                        nullptr,
    1318:                                        nullptr,
    1319:                                        nullptr,
    1320:                                        /*140*/ nullptr,
    1321:                                        nullptr,
    1322:                                        nullptr,
    1323:                                        nullptr,
    1324:                                        nullptr,
    1325:                                        /*145*/ nullptr,
    1326:                                        nullptr,
    1327:                                        nullptr,
    1328:                                        nullptr,
    1329:                                        nullptr,
    1330:                                        /*150*/ nullptr,
    1331:                                        nullptr,
    1332:                                        nullptr,
    1333:                                        nullptr,
    1334:                                        nullptr,
    1335:                                        nullptr};

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

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

    1342:    Options Database Keys:
    1343: + -mat_type elemental - sets the matrix type to "elemental" during a call to MatSetFromOptions()
    1344: . -pc_factor_mat_solver_type elemental - to use this direct solver with the option -pc_type lu
    1345: - -mat_elemental_grid_height - sets Grid Height for 2D cyclic ordering of internal matrix

    1347:   Level: beginner

    1349:   Note:
    1350:    Note unlike most matrix formats, this format does not store all the matrix entries for a contiguous
    1351:    range of rows on an MPI rank. Use `MatGetOwnershipIS()` to determine what values are stored on
    1352:    the given rank.

    1354: .seealso: `MATDENSE`, `MATSCALAPACK`, `MatGetOwnershipIS()`
    1355: M*/
    1356: #if defined(__clang__)
    1357:   #pragma clang diagnostic push
    1358:   #pragma clang diagnostic ignored "-Wzero-as-null-pointer-constant"
    1359: #endif
    1360: PETSC_EXTERN PetscErrorCode MatCreate_Elemental(Mat A)
    1361: {
    1362:   Mat_Elemental      *a;
    1363:   PetscBool           flg, flg1;
    1364:   Mat_Elemental_Grid *commgrid;
    1365:   MPI_Comm            icomm;
    1366:   PetscInt            optv1;

    1368:   PetscFunctionBegin;
    1369:   A->ops[0]     = MatOps_Values;
    1370:   A->insertmode = NOT_SET_VALUES;

    1372:   PetscCall(PetscNew(&a));
    1373:   A->data = (void *)a;

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

    1378:   /* Grid needs to be shared between multiple Mats on the same communicator, implement by attribute caching on the MPI_Comm */
    1379:   if (Petsc_Elemental_keyval == MPI_KEYVAL_INVALID) {
    1380:     PetscCallMPI(MPI_Comm_create_keyval(MPI_COMM_NULL_COPY_FN, MPI_COMM_NULL_DELETE_FN, &Petsc_Elemental_keyval, nullptr));
    1381:     PetscCall(PetscCitationsRegister(ElementalCitation, &ElementalCite));
    1382:   }
    1383:   PetscCall(PetscCommDuplicate(cxxcomm.comm, &icomm, NULL));
    1384:   PetscCallMPI(MPI_Comm_get_attr(icomm, Petsc_Elemental_keyval, (void **)&commgrid, (int *)&flg));
    1385:   if (!flg) {
    1386:     PetscCall(PetscNew(&commgrid));

    1388:     PetscOptionsBegin(PetscObjectComm((PetscObject)A), ((PetscObject)A)->prefix, "Elemental Options", "Mat");
    1389:     /* displayed default grid sizes (CommSize,1) are set by us arbitrarily until El::Grid() is called */
    1390:     PetscCall(PetscOptionsInt("-mat_elemental_grid_height", "Grid Height", "None", El::mpi::Size(cxxcomm), &optv1, &flg1));
    1391:     if (flg1) {
    1392:       PetscCheck((El::mpi::Size(cxxcomm) % optv1) == 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Grid Height %" PetscInt_FMT " must evenly divide CommSize %" PetscInt_FMT, optv1, (PetscInt)El::mpi::Size(cxxcomm));
    1393:       commgrid->grid = new El::Grid(cxxcomm, optv1); /* use user-provided grid height */
    1394:     } else {
    1395:       commgrid->grid = new El::Grid(cxxcomm); /* use Elemental default grid sizes */
    1396:       /* printf("new commgrid->grid = %p\n",commgrid->grid);  -- memory leak revealed by valgrind? */
    1397:     }
    1398:     commgrid->grid_refct = 1;
    1399:     PetscCallMPI(MPI_Comm_set_attr(icomm, Petsc_Elemental_keyval, (void *)commgrid));

    1401:     a->pivoting = 1;
    1402:     PetscCall(PetscOptionsInt("-mat_elemental_pivoting", "Pivoting", "None", a->pivoting, &a->pivoting, NULL));

    1404:     PetscOptionsEnd();
    1405:   } else {
    1406:     commgrid->grid_refct++;
    1407:   }
    1408:   PetscCall(PetscCommDestroy(&icomm));
    1409:   a->grid        = commgrid->grid;
    1410:   a->emat        = new El::DistMatrix<PetscElemScalar>(*a->grid);
    1411:   a->roworiented = PETSC_TRUE;

    1413:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatGetOwnershipIS_C", MatGetOwnershipIS_Elemental));
    1414:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_elemental_mpidense_C", MatProductSetFromOptions_Elemental_MPIDense));
    1415:   PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATELEMENTAL));
    1416:   PetscFunctionReturn(PETSC_SUCCESS);
    1417: }
    1418: #if defined(__clang__)
    1419:   #pragma clang diagnostic pop
    1420: #endif