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:     //PetscCall(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:     //PetscCall(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 { /* columnoriented */
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_DEFAULT, &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_EXTERN 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};

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

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

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

1343:   Level: beginner

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

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

1364:   PetscFunctionBegin;
1365:   A->ops[0]     = MatOps_Values;
1366:   A->insertmode = NOT_SET_VALUES;

1368:   PetscCall(PetscNew(&a));
1369:   A->data = (void *)a;

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

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

1384:     PetscOptionsBegin(PetscObjectComm((PetscObject)A), ((PetscObject)A)->prefix, "Elemental Options", "Mat");
1385:     /* displayed default grid sizes (CommSize,1) are set by us arbitrarily until El::Grid() is called */
1386:     PetscCall(PetscOptionsInt("-mat_elemental_grid_height", "Grid Height", "None", El::mpi::Size(cxxcomm), &optv1, &flg1));
1387:     if (flg1) {
1388:       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));
1389:       commgrid->grid = new El::Grid(cxxcomm, optv1); /* use user-provided grid height */
1390:     } else {
1391:       commgrid->grid = new El::Grid(cxxcomm); /* use Elemental default grid sizes */
1392:       /* printf("new commgrid->grid = %p\n",commgrid->grid);  -- memory leak revealed by valgrind? */
1393:     }
1394:     commgrid->grid_refct = 1;
1395:     PetscCallMPI(MPI_Comm_set_attr(icomm, Petsc_Elemental_keyval, (void *)commgrid));

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

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

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