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

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

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

 63:   info->block_size = 1.0;

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

 80:   info->nz_unneeded       = 0.0;
 81:   info->assemblies        = A->num_ass;
 82:   info->mallocs           = 0;
 83:   info->memory            = 0; /* REVIEW ME */
 84:   info->fill_ratio_given  = 0; /* determined by Elemental */
 85:   info->fill_ratio_needed = 0;
 86:   info->factor_mallocs    = 0;
 87:   return 0;
 88: }

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

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

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

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

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

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

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

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

218: static PetscErrorCode MatMult_Elemental(Mat A, Vec X, Vec Y)
219: {
220:   Mat_Elemental         *a = (Mat_Elemental *)A->data;
221:   const PetscElemScalar *x;
222:   PetscElemScalar       *y;
223:   PetscElemScalar        one = 1, zero = 0;

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

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

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

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

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

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

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

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

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

314: PetscErrorCode MatMatMultSymbolic_Elemental(Mat A, Mat B, PetscReal fill, Mat Ce)
315: {
316:   MatSetSizes(Ce, A->rmap->n, B->cmap->n, PETSC_DECIDE, PETSC_DECIDE);
317:   MatSetType(Ce, MATELEMENTAL);
318:   MatSetUp(Ce);
319:   Ce->ops->matmultnumeric = MatMatMultNumeric_Elemental;
320:   return 0;
321: }

323: static PetscErrorCode MatMatTransposeMultNumeric_Elemental(Mat A, Mat B, Mat C)
324: {
325:   Mat_Elemental  *a   = (Mat_Elemental *)A->data;
326:   Mat_Elemental  *b   = (Mat_Elemental *)B->data;
327:   Mat_Elemental  *c   = (Mat_Elemental *)C->data;
328:   PetscElemScalar one = 1, zero = 0;

330:   { /* Scoping so that constructor is called before pointer is returned */
331:     El::Gemm(El::NORMAL, El::TRANSPOSE, one, *a->emat, *b->emat, zero, *c->emat);
332:   }
333:   C->assembled = PETSC_TRUE;
334:   return 0;
335: }

337: static PetscErrorCode MatMatTransposeMultSymbolic_Elemental(Mat A, Mat B, PetscReal fill, Mat C)
338: {
339:   MatSetSizes(C, A->rmap->n, B->rmap->n, PETSC_DECIDE, PETSC_DECIDE);
340:   MatSetType(C, MATELEMENTAL);
341:   MatSetUp(C);
342:   return 0;
343: }

345: /* --------------------------------------- */
346: static PetscErrorCode MatProductSetFromOptions_Elemental_AB(Mat C)
347: {
348:   C->ops->matmultsymbolic = MatMatMultSymbolic_Elemental;
349:   C->ops->productsymbolic = MatProductSymbolic_AB;
350:   return 0;
351: }

353: static PetscErrorCode MatProductSetFromOptions_Elemental_ABt(Mat C)
354: {
355:   C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_Elemental;
356:   C->ops->productsymbolic          = MatProductSymbolic_ABt;
357:   return 0;
358: }

360: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Elemental(Mat C)
361: {
362:   Mat_Product *product = C->product;

364:   switch (product->type) {
365:   case MATPRODUCT_AB:
366:     MatProductSetFromOptions_Elemental_AB(C);
367:     break;
368:   case MATPRODUCT_ABt:
369:     MatProductSetFromOptions_Elemental_ABt(C);
370:     break;
371:   default:
372:     break;
373:   }
374:   return 0;
375: }

377: PetscErrorCode MatMatMultNumeric_Elemental_MPIDense(Mat A, Mat B, Mat C)
378: {
379:   Mat Be, Ce;

381:   MatConvert(B, MATELEMENTAL, MAT_INITIAL_MATRIX, &Be);
382:   MatMatMult(A, Be, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Ce);
383:   MatConvert(Ce, MATMPIDENSE, MAT_REUSE_MATRIX, &C);
384:   MatDestroy(&Be);
385:   MatDestroy(&Ce);
386:   return 0;
387: }

389: PetscErrorCode MatMatMultSymbolic_Elemental_MPIDense(Mat A, Mat B, PetscReal fill, Mat C)
390: {
391:   MatSetSizes(C, A->rmap->n, B->cmap->n, PETSC_DECIDE, PETSC_DECIDE);
392:   MatSetType(C, MATMPIDENSE);
393:   MatSetUp(C);
394:   C->ops->matmultnumeric = MatMatMultNumeric_Elemental_MPIDense;
395:   return 0;
396: }

398: PetscErrorCode MatProductSetFromOptions_Elemental_MPIDense_AB(Mat C)
399: {
400:   C->ops->matmultsymbolic = MatMatMultSymbolic_Elemental_MPIDense;
401:   C->ops->productsymbolic = MatProductSymbolic_AB;
402:   return 0;
403: }

405: PetscErrorCode MatProductSetFromOptions_Elemental_MPIDense(Mat C)
406: {
407:   Mat_Product *product = C->product;

409:   if (product->type == MATPRODUCT_AB) MatProductSetFromOptions_Elemental_MPIDense_AB(C);
410:   return 0;
411: }
412: /* --------------------------------------- */

414: static PetscErrorCode MatGetDiagonal_Elemental(Mat A, Vec D)
415: {
416:   PetscInt        i, nrows, ncols, nD, rrank, ridx, crank, cidx;
417:   Mat_Elemental  *a = (Mat_Elemental *)A->data;
418:   PetscElemScalar v;
419:   MPI_Comm        comm;

421:   PetscObjectGetComm((PetscObject)A, &comm);
422:   MatGetSize(A, &nrows, &ncols);
423:   nD = nrows > ncols ? ncols : nrows;
424:   for (i = 0; i < nD; i++) {
425:     PetscInt erow, ecol;
426:     P2RO(A, 0, i, &rrank, &ridx);
427:     RO2E(A, 0, rrank, ridx, &erow);
429:     P2RO(A, 1, i, &crank, &cidx);
430:     RO2E(A, 1, crank, cidx, &ecol);
432:     v = a->emat->Get(erow, ecol);
433:     VecSetValues(D, 1, &i, (PetscScalar *)&v, INSERT_VALUES);
434:   }
435:   VecAssemblyBegin(D);
436:   VecAssemblyEnd(D);
437:   return 0;
438: }

440: static PetscErrorCode MatDiagonalScale_Elemental(Mat X, Vec L, Vec R)
441: {
442:   Mat_Elemental         *x = (Mat_Elemental *)X->data;
443:   const PetscElemScalar *d;

445:   if (R) {
446:     VecGetArrayRead(R, (const PetscScalar **)&d);
447:     El::DistMatrix<PetscElemScalar, El::VC, El::STAR> de;
448:     de.LockedAttach(X->cmap->N, 1, *x->grid, 0, 0, d, X->cmap->n);
449:     El::DiagonalScale(El::RIGHT, El::NORMAL, de, *x->emat);
450:     VecRestoreArrayRead(R, (const PetscScalar **)&d);
451:   }
452:   if (L) {
453:     VecGetArrayRead(L, (const PetscScalar **)&d);
454:     El::DistMatrix<PetscElemScalar, El::VC, El::STAR> de;
455:     de.LockedAttach(X->rmap->N, 1, *x->grid, 0, 0, d, X->rmap->n);
456:     El::DiagonalScale(El::LEFT, El::NORMAL, de, *x->emat);
457:     VecRestoreArrayRead(L, (const PetscScalar **)&d);
458:   }
459:   return 0;
460: }

462: static PetscErrorCode MatMissingDiagonal_Elemental(Mat A, PetscBool *missing, PetscInt *d)
463: {
464:   *missing = PETSC_FALSE;
465:   return 0;
466: }

468: static PetscErrorCode MatScale_Elemental(Mat X, PetscScalar a)
469: {
470:   Mat_Elemental *x = (Mat_Elemental *)X->data;

472:   El::Scale((PetscElemScalar)a, *x->emat);
473:   return 0;
474: }

476: /*
477:   MatAXPY - Computes Y = a*X + Y.
478: */
479: static PetscErrorCode MatAXPY_Elemental(Mat Y, PetscScalar a, Mat X, MatStructure str)
480: {
481:   Mat_Elemental *x = (Mat_Elemental *)X->data;
482:   Mat_Elemental *y = (Mat_Elemental *)Y->data;

484:   El::Axpy((PetscElemScalar)a, *x->emat, *y->emat);
485:   PetscObjectStateIncrease((PetscObject)Y);
486:   return 0;
487: }

489: static PetscErrorCode MatCopy_Elemental(Mat A, Mat B, MatStructure str)
490: {
491:   Mat_Elemental *a = (Mat_Elemental *)A->data;
492:   Mat_Elemental *b = (Mat_Elemental *)B->data;

494:   El::Copy(*a->emat, *b->emat);
495:   PetscObjectStateIncrease((PetscObject)B);
496:   return 0;
497: }

499: static PetscErrorCode MatDuplicate_Elemental(Mat A, MatDuplicateOption op, Mat *B)
500: {
501:   Mat            Be;
502:   MPI_Comm       comm;
503:   Mat_Elemental *a = (Mat_Elemental *)A->data;

505:   PetscObjectGetComm((PetscObject)A, &comm);
506:   MatCreate(comm, &Be);
507:   MatSetSizes(Be, A->rmap->n, A->cmap->n, PETSC_DECIDE, PETSC_DECIDE);
508:   MatSetType(Be, MATELEMENTAL);
509:   MatSetUp(Be);
510:   *B = Be;
511:   if (op == MAT_COPY_VALUES) {
512:     Mat_Elemental *b = (Mat_Elemental *)Be->data;
513:     El::Copy(*a->emat, *b->emat);
514:   }
515:   Be->assembled = PETSC_TRUE;
516:   return 0;
517: }

519: static PetscErrorCode MatTranspose_Elemental(Mat A, MatReuse reuse, Mat *B)
520: {
521:   Mat            Be = *B;
522:   MPI_Comm       comm;
523:   Mat_Elemental *a = (Mat_Elemental *)A->data, *b;

525:   if (reuse == MAT_REUSE_MATRIX) MatTransposeCheckNonzeroState_Private(A, *B);
526:   PetscObjectGetComm((PetscObject)A, &comm);
527:   /* Only out-of-place supported */
529:   if (reuse == MAT_INITIAL_MATRIX) {
530:     MatCreate(comm, &Be);
531:     MatSetSizes(Be, A->cmap->n, A->rmap->n, PETSC_DECIDE, PETSC_DECIDE);
532:     MatSetType(Be, MATELEMENTAL);
533:     MatSetUp(Be);
534:     *B = Be;
535:   }
536:   b = (Mat_Elemental *)Be->data;
537:   El::Transpose(*a->emat, *b->emat);
538:   Be->assembled = PETSC_TRUE;
539:   return 0;
540: }

542: static PetscErrorCode MatConjugate_Elemental(Mat A)
543: {
544:   Mat_Elemental *a = (Mat_Elemental *)A->data;

546:   El::Conjugate(*a->emat);
547:   return 0;
548: }

550: static PetscErrorCode MatHermitianTranspose_Elemental(Mat A, MatReuse reuse, Mat *B)
551: {
552:   Mat            Be = *B;
553:   MPI_Comm       comm;
554:   Mat_Elemental *a = (Mat_Elemental *)A->data, *b;

556:   PetscObjectGetComm((PetscObject)A, &comm);
557:   /* Only out-of-place supported */
558:   if (reuse == MAT_INITIAL_MATRIX) {
559:     MatCreate(comm, &Be);
560:     MatSetSizes(Be, A->cmap->n, A->rmap->n, PETSC_DECIDE, PETSC_DECIDE);
561:     MatSetType(Be, MATELEMENTAL);
562:     MatSetUp(Be);
563:     *B = Be;
564:   }
565:   b = (Mat_Elemental *)Be->data;
566:   El::Adjoint(*a->emat, *b->emat);
567:   Be->assembled = PETSC_TRUE;
568:   return 0;
569: }

571: static PetscErrorCode MatSolve_Elemental(Mat A, Vec B, Vec X)
572: {
573:   Mat_Elemental   *a = (Mat_Elemental *)A->data;
574:   PetscElemScalar *x;
575:   PetscInt         pivoting = a->pivoting;

577:   VecCopy(B, X);
578:   VecGetArray(X, (PetscScalar **)&x);

580:   El::DistMatrix<PetscElemScalar, El::VC, El::STAR> xe;
581:   xe.Attach(A->rmap->N, 1, *a->grid, 0, 0, x, A->rmap->n);
582:   El::DistMatrix<PetscElemScalar, El::MC, El::MR> xer(xe);
583:   switch (A->factortype) {
584:   case MAT_FACTOR_LU:
585:     if (pivoting == 0) {
586:       El::lu::SolveAfter(El::NORMAL, *a->emat, xer);
587:     } else if (pivoting == 1) {
588:       El::lu::SolveAfter(El::NORMAL, *a->emat, *a->P, xer);
589:     } else { /* pivoting == 2 */
590:       El::lu::SolveAfter(El::NORMAL, *a->emat, *a->P, *a->Q, xer);
591:     }
592:     break;
593:   case MAT_FACTOR_CHOLESKY:
594:     El::cholesky::SolveAfter(El::UPPER, El::NORMAL, *a->emat, xer);
595:     break;
596:   default:
597:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unfactored Matrix or Unsupported MatFactorType");
598:     break;
599:   }
600:   El::Copy(xer, xe);

602:   VecRestoreArray(X, (PetscScalar **)&x);
603:   return 0;
604: }

606: static PetscErrorCode MatSolveAdd_Elemental(Mat A, Vec B, Vec Y, Vec X)
607: {
608:   MatSolve_Elemental(A, B, X);
609:   VecAXPY(X, 1, Y);
610:   return 0;
611: }

613: static PetscErrorCode MatMatSolve_Elemental(Mat A, Mat B, Mat X)
614: {
615:   Mat_Elemental *a = (Mat_Elemental *)A->data;
616:   Mat_Elemental *x;
617:   Mat            C;
618:   PetscInt       pivoting = a->pivoting;
619:   PetscBool      flg;
620:   MatType        type;

622:   MatGetType(X, &type);
623:   PetscStrcmp(type, MATELEMENTAL, &flg);
624:   if (!flg) {
625:     MatConvert(B, MATELEMENTAL, MAT_INITIAL_MATRIX, &C);
626:     x = (Mat_Elemental *)C->data;
627:   } else {
628:     x = (Mat_Elemental *)X->data;
629:     El::Copy(*((Mat_Elemental *)B->data)->emat, *x->emat);
630:   }
631:   switch (A->factortype) {
632:   case MAT_FACTOR_LU:
633:     if (pivoting == 0) {
634:       El::lu::SolveAfter(El::NORMAL, *a->emat, *x->emat);
635:     } else if (pivoting == 1) {
636:       El::lu::SolveAfter(El::NORMAL, *a->emat, *a->P, *x->emat);
637:     } else {
638:       El::lu::SolveAfter(El::NORMAL, *a->emat, *a->P, *a->Q, *x->emat);
639:     }
640:     break;
641:   case MAT_FACTOR_CHOLESKY:
642:     El::cholesky::SolveAfter(El::UPPER, El::NORMAL, *a->emat, *x->emat);
643:     break;
644:   default:
645:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unfactored Matrix or Unsupported MatFactorType");
646:     break;
647:   }
648:   if (!flg) {
649:     MatConvert(C, type, MAT_REUSE_MATRIX, &X);
650:     MatDestroy(&C);
651:   }
652:   return 0;
653: }

655: static PetscErrorCode MatLUFactor_Elemental(Mat A, IS row, IS col, const MatFactorInfo *info)
656: {
657:   Mat_Elemental *a        = (Mat_Elemental *)A->data;
658:   PetscInt       pivoting = a->pivoting;

660:   if (pivoting == 0) {
661:     El::LU(*a->emat);
662:   } else if (pivoting == 1) {
663:     El::LU(*a->emat, *a->P);
664:   } else {
665:     El::LU(*a->emat, *a->P, *a->Q);
666:   }
667:   A->factortype = MAT_FACTOR_LU;
668:   A->assembled  = PETSC_TRUE;

670:   PetscFree(A->solvertype);
671:   PetscStrallocpy(MATSOLVERELEMENTAL, &A->solvertype);
672:   return 0;
673: }

675: static PetscErrorCode MatLUFactorNumeric_Elemental(Mat F, Mat A, const MatFactorInfo *info)
676: {
677:   MatCopy(A, F, SAME_NONZERO_PATTERN);
678:   MatLUFactor_Elemental(F, 0, 0, info);
679:   return 0;
680: }

682: static PetscErrorCode MatLUFactorSymbolic_Elemental(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info)
683: {
684:   /* F is created and allocated by MatGetFactor_elemental_petsc(), skip this routine. */
685:   return 0;
686: }

688: static PetscErrorCode MatCholeskyFactor_Elemental(Mat A, IS perm, const MatFactorInfo *info)
689: {
690:   Mat_Elemental                                    *a = (Mat_Elemental *)A->data;
691:   El::DistMatrix<PetscElemScalar, El::MC, El::STAR> d;

693:   El::Cholesky(El::UPPER, *a->emat);
694:   A->factortype = MAT_FACTOR_CHOLESKY;
695:   A->assembled  = PETSC_TRUE;

697:   PetscFree(A->solvertype);
698:   PetscStrallocpy(MATSOLVERELEMENTAL, &A->solvertype);
699:   return 0;
700: }

702: static PetscErrorCode MatCholeskyFactorNumeric_Elemental(Mat F, Mat A, const MatFactorInfo *info)
703: {
704:   MatCopy(A, F, SAME_NONZERO_PATTERN);
705:   MatCholeskyFactor_Elemental(F, 0, info);
706:   return 0;
707: }

709: static PetscErrorCode MatCholeskyFactorSymbolic_Elemental(Mat F, Mat A, IS perm, const MatFactorInfo *info)
710: {
711:   /* F is created and allocated by MatGetFactor_elemental_petsc(), skip this routine. */
712:   return 0;
713: }

715: PetscErrorCode MatFactorGetSolverType_elemental_elemental(Mat A, MatSolverType *type)
716: {
717:   *type = MATSOLVERELEMENTAL;
718:   return 0;
719: }

721: static PetscErrorCode MatGetFactor_elemental_elemental(Mat A, MatFactorType ftype, Mat *F)
722: {
723:   Mat B;

725:   /* Create the factorization matrix */
726:   MatCreate(PetscObjectComm((PetscObject)A), &B);
727:   MatSetSizes(B, A->rmap->n, A->cmap->n, PETSC_DECIDE, PETSC_DECIDE);
728:   MatSetType(B, MATELEMENTAL);
729:   MatSetUp(B);
730:   B->factortype      = ftype;
731:   B->trivialsymbolic = PETSC_TRUE;
732:   PetscFree(B->solvertype);
733:   PetscStrallocpy(MATSOLVERELEMENTAL, &B->solvertype);

735:   PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_elemental_elemental);
736:   *F = B;
737:   return 0;
738: }

740: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_Elemental(void)
741: {
742:   MatSolverTypeRegister(MATSOLVERELEMENTAL, MATELEMENTAL, MAT_FACTOR_LU, MatGetFactor_elemental_elemental);
743:   MatSolverTypeRegister(MATSOLVERELEMENTAL, MATELEMENTAL, MAT_FACTOR_CHOLESKY, MatGetFactor_elemental_elemental);
744:   return 0;
745: }

747: static PetscErrorCode MatNorm_Elemental(Mat A, NormType type, PetscReal *nrm)
748: {
749:   Mat_Elemental *a = (Mat_Elemental *)A->data;

751:   switch (type) {
752:   case NORM_1:
753:     *nrm = El::OneNorm(*a->emat);
754:     break;
755:   case NORM_FROBENIUS:
756:     *nrm = El::FrobeniusNorm(*a->emat);
757:     break;
758:   case NORM_INFINITY:
759:     *nrm = El::InfinityNorm(*a->emat);
760:     break;
761:   default:
762:     SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Unsupported norm type");
763:   }
764:   return 0;
765: }

767: static PetscErrorCode MatZeroEntries_Elemental(Mat A)
768: {
769:   Mat_Elemental *a = (Mat_Elemental *)A->data;

771:   El::Zero(*a->emat);
772:   return 0;
773: }

775: static PetscErrorCode MatGetOwnershipIS_Elemental(Mat A, IS *rows, IS *cols)
776: {
777:   Mat_Elemental *a = (Mat_Elemental *)A->data;
778:   PetscInt       i, m, shift, stride, *idx;

780:   if (rows) {
781:     m      = a->emat->LocalHeight();
782:     shift  = a->emat->ColShift();
783:     stride = a->emat->ColStride();
784:     PetscMalloc1(m, &idx);
785:     for (i = 0; i < m; i++) {
786:       PetscInt rank, offset;
787:       E2RO(A, 0, shift + i * stride, &rank, &offset);
788:       RO2P(A, 0, rank, offset, &idx[i]);
789:     }
790:     ISCreateGeneral(PETSC_COMM_SELF, m, idx, PETSC_OWN_POINTER, rows);
791:   }
792:   if (cols) {
793:     m      = a->emat->LocalWidth();
794:     shift  = a->emat->RowShift();
795:     stride = a->emat->RowStride();
796:     PetscMalloc1(m, &idx);
797:     for (i = 0; i < m; i++) {
798:       PetscInt rank, offset;
799:       E2RO(A, 1, shift + i * stride, &rank, &offset);
800:       RO2P(A, 1, rank, offset, &idx[i]);
801:     }
802:     ISCreateGeneral(PETSC_COMM_SELF, m, idx, PETSC_OWN_POINTER, cols);
803:   }
804:   return 0;
805: }

807: static PetscErrorCode MatConvert_Elemental_Dense(Mat A, MatType newtype, MatReuse reuse, Mat *B)
808: {
809:   Mat             Bmpi;
810:   Mat_Elemental  *a = (Mat_Elemental *)A->data;
811:   MPI_Comm        comm;
812:   IS              isrows, iscols;
813:   PetscInt        rrank, ridx, crank, cidx, nrows, ncols, i, j, erow, ecol, elrow, elcol;
814:   const PetscInt *rows, *cols;
815:   PetscElemScalar v;
816:   const El::Grid &grid = a->emat->Grid();

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

820:   if (reuse == MAT_REUSE_MATRIX) {
821:     Bmpi = *B;
822:   } else {
823:     MatCreate(comm, &Bmpi);
824:     MatSetSizes(Bmpi, A->rmap->n, A->cmap->n, PETSC_DECIDE, PETSC_DECIDE);
825:     MatSetType(Bmpi, MATDENSE);
826:     MatSetUp(Bmpi);
827:   }

829:   /* Get local entries of A */
830:   MatGetOwnershipIS(A, &isrows, &iscols);
831:   ISGetLocalSize(isrows, &nrows);
832:   ISGetIndices(isrows, &rows);
833:   ISGetLocalSize(iscols, &ncols);
834:   ISGetIndices(iscols, &cols);

836:   if (a->roworiented) {
837:     for (i = 0; i < nrows; i++) {
838:       P2RO(A, 0, rows[i], &rrank, &ridx); /* convert indices between PETSc <-> (Rank,Offset) <-> Elemental */
839:       RO2E(A, 0, rrank, ridx, &erow);
841:       for (j = 0; j < ncols; j++) {
842:         P2RO(A, 1, cols[j], &crank, &cidx);
843:         RO2E(A, 1, crank, cidx, &ecol);

846:         elrow = erow / grid.MCSize(); /* Elemental local row index */
847:         elcol = ecol / grid.MRSize(); /* Elemental local column index */
848:         v     = a->emat->GetLocal(elrow, elcol);
849:         MatSetValues(Bmpi, 1, &rows[i], 1, &cols[j], (PetscScalar *)&v, INSERT_VALUES);
850:       }
851:     }
852:   } else { /* column-oriented */
853:     for (j = 0; j < ncols; j++) {
854:       P2RO(A, 1, cols[j], &crank, &cidx);
855:       RO2E(A, 1, crank, cidx, &ecol);
857:       for (i = 0; i < nrows; i++) {
858:         P2RO(A, 0, rows[i], &rrank, &ridx); /* convert indices between PETSc <-> (Rank,Offset) <-> Elemental */
859:         RO2E(A, 0, rrank, ridx, &erow);

862:         elrow = erow / grid.MCSize(); /* Elemental local row index */
863:         elcol = ecol / grid.MRSize(); /* Elemental local column index */
864:         v     = a->emat->GetLocal(elrow, elcol);
865:         MatSetValues(Bmpi, 1, &rows[i], 1, &cols[j], (PetscScalar *)&v, INSERT_VALUES);
866:       }
867:     }
868:   }
869:   MatAssemblyBegin(Bmpi, MAT_FINAL_ASSEMBLY);
870:   MatAssemblyEnd(Bmpi, MAT_FINAL_ASSEMBLY);
871:   if (reuse == MAT_INPLACE_MATRIX) {
872:     MatHeaderReplace(A, &Bmpi);
873:   } else {
874:     *B = Bmpi;
875:   }
876:   ISDestroy(&isrows);
877:   ISDestroy(&iscols);
878:   return 0;
879: }

881: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
882: {
883:   Mat                mat_elemental;
884:   PetscInt           M = A->rmap->N, N = A->cmap->N, row, ncols;
885:   const PetscInt    *cols;
886:   const PetscScalar *vals;

888:   if (reuse == MAT_REUSE_MATRIX) {
889:     mat_elemental = *newmat;
890:     MatZeroEntries(mat_elemental);
891:   } else {
892:     MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
893:     MatSetSizes(mat_elemental, PETSC_DECIDE, PETSC_DECIDE, M, N);
894:     MatSetType(mat_elemental, MATELEMENTAL);
895:     MatSetUp(mat_elemental);
896:   }
897:   for (row = 0; row < M; row++) {
898:     MatGetRow(A, row, &ncols, &cols, &vals);
899:     /* PETSc-Elemental interface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
900:     MatSetValues(mat_elemental, 1, &row, ncols, cols, vals, ADD_VALUES);
901:     MatRestoreRow(A, row, &ncols, &cols, &vals);
902:   }
903:   MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
904:   MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);

906:   if (reuse == MAT_INPLACE_MATRIX) {
907:     MatHeaderReplace(A, &mat_elemental);
908:   } else {
909:     *newmat = mat_elemental;
910:   }
911:   return 0;
912: }

914: PETSC_INTERN PetscErrorCode MatConvert_MPIAIJ_Elemental(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
915: {
916:   Mat                mat_elemental;
917:   PetscInt           row, ncols, rstart = A->rmap->rstart, rend = A->rmap->rend, j;
918:   const PetscInt    *cols;
919:   const PetscScalar *vals;

921:   if (reuse == MAT_REUSE_MATRIX) {
922:     mat_elemental = *newmat;
923:     MatZeroEntries(mat_elemental);
924:   } else {
925:     MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
926:     MatSetSizes(mat_elemental, PETSC_DECIDE, PETSC_DECIDE, A->rmap->N, A->cmap->N);
927:     MatSetType(mat_elemental, MATELEMENTAL);
928:     MatSetUp(mat_elemental);
929:   }
930:   for (row = rstart; row < rend; row++) {
931:     MatGetRow(A, row, &ncols, &cols, &vals);
932:     for (j = 0; j < ncols; j++) {
933:       /* PETSc-Elemental interface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
934:       MatSetValues(mat_elemental, 1, &row, 1, &cols[j], &vals[j], ADD_VALUES);
935:     }
936:     MatRestoreRow(A, row, &ncols, &cols, &vals);
937:   }
938:   MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
939:   MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);

941:   if (reuse == MAT_INPLACE_MATRIX) {
942:     MatHeaderReplace(A, &mat_elemental);
943:   } else {
944:     *newmat = mat_elemental;
945:   }
946:   return 0;
947: }

949: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_Elemental(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
950: {
951:   Mat                mat_elemental;
952:   PetscInt           M = A->rmap->N, N = A->cmap->N, row, ncols, j;
953:   const PetscInt    *cols;
954:   const PetscScalar *vals;

956:   if (reuse == MAT_REUSE_MATRIX) {
957:     mat_elemental = *newmat;
958:     MatZeroEntries(mat_elemental);
959:   } else {
960:     MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
961:     MatSetSizes(mat_elemental, PETSC_DECIDE, PETSC_DECIDE, M, N);
962:     MatSetType(mat_elemental, MATELEMENTAL);
963:     MatSetUp(mat_elemental);
964:   }
965:   MatGetRowUpperTriangular(A);
966:   for (row = 0; row < M; row++) {
967:     MatGetRow(A, row, &ncols, &cols, &vals);
968:     /* PETSc-Elemental interface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
969:     MatSetValues(mat_elemental, 1, &row, ncols, cols, vals, ADD_VALUES);
970:     for (j = 0; j < ncols; j++) { /* lower triangular part */
971:       PetscScalar v;
972:       if (cols[j] == row) continue;
973:       v = A->hermitian == PETSC_BOOL3_TRUE ? PetscConj(vals[j]) : vals[j];
974:       MatSetValues(mat_elemental, 1, &cols[j], 1, &row, &v, ADD_VALUES);
975:     }
976:     MatRestoreRow(A, row, &ncols, &cols, &vals);
977:   }
978:   MatRestoreRowUpperTriangular(A);
979:   MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
980:   MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);

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

990: PETSC_INTERN PetscErrorCode MatConvert_MPISBAIJ_Elemental(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
991: {
992:   Mat                mat_elemental;
993:   PetscInt           M = A->rmap->N, N = A->cmap->N, row, ncols, j, rstart = A->rmap->rstart, rend = A->rmap->rend;
994:   const PetscInt    *cols;
995:   const PetscScalar *vals;

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

1023:   if (reuse == MAT_INPLACE_MATRIX) {
1024:     MatHeaderReplace(A, &mat_elemental);
1025:   } else {
1026:     *newmat = mat_elemental;
1027:   }
1028:   return 0;
1029: }

1031: static PetscErrorCode MatDestroy_Elemental(Mat A)
1032: {
1033:   Mat_Elemental      *a = (Mat_Elemental *)A->data;
1034:   Mat_Elemental_Grid *commgrid;
1035:   PetscBool           flg;
1036:   MPI_Comm            icomm;

1038:   delete a->emat;
1039:   delete a->P;
1040:   delete a->Q;

1042:   El::mpi::Comm cxxcomm(PetscObjectComm((PetscObject)A));
1043:   PetscCommDuplicate(cxxcomm.comm, &icomm, NULL);
1044:   MPI_Comm_get_attr(icomm, Petsc_Elemental_keyval, (void **)&commgrid, (int *)&flg);
1045:   if (--commgrid->grid_refct == 0) {
1046:     delete commgrid->grid;
1047:     PetscFree(commgrid);
1048:     MPI_Comm_free_keyval(&Petsc_Elemental_keyval);
1049:   }
1050:   PetscCommDestroy(&icomm);
1051:   PetscObjectComposeFunction((PetscObject)A, "MatGetOwnershipIS_C", NULL);
1052:   PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL);
1053:   PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_elemental_mpidense_C", NULL);
1054:   PetscFree(A->data);
1055:   return 0;
1056: }

1058: PetscErrorCode MatSetUp_Elemental(Mat A)
1059: {
1060:   Mat_Elemental *a = (Mat_Elemental *)A->data;
1061:   MPI_Comm       comm;
1062:   PetscMPIInt    rsize, csize;
1063:   PetscInt       n;

1065:   PetscLayoutSetUp(A->rmap);
1066:   PetscLayoutSetUp(A->cmap);

1068:   /* Check if local row and column sizes are equally distributed.
1069:      Jed: Elemental uses "element" cyclic ordering so the sizes need to match that
1070:      exactly.  The strategy in MatElemental is for PETSc to implicitly permute to block ordering (like would be returned by
1071:      PetscSplitOwnership(comm,&n,&N), at which point Elemental matrices can act on PETSc vectors without redistributing the vectors. */
1072:   PetscObjectGetComm((PetscObject)A, &comm);
1073:   n = PETSC_DECIDE;
1074:   PetscSplitOwnership(comm, &n, &A->rmap->N);

1077:   n = PETSC_DECIDE;
1078:   PetscSplitOwnership(comm, &n, &A->cmap->N);

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

1084:   MPI_Comm_size(A->rmap->comm, &rsize);
1085:   MPI_Comm_size(A->cmap->comm, &csize);
1087:   a->commsize = rsize;
1088:   a->mr[0]    = A->rmap->N % rsize;
1089:   if (!a->mr[0]) a->mr[0] = rsize;
1090:   a->mr[1] = A->cmap->N % csize;
1091:   if (!a->mr[1]) a->mr[1] = csize;
1092:   a->m[0] = A->rmap->N / rsize + (a->mr[0] != rsize);
1093:   a->m[1] = A->cmap->N / csize + (a->mr[1] != csize);
1094:   return 0;
1095: }

1097: PetscErrorCode MatAssemblyBegin_Elemental(Mat A, MatAssemblyType type)
1098: {
1099:   Mat_Elemental *a = (Mat_Elemental *)A->data;

1101:   /* printf("Calling ProcessQueues\n"); */
1102:   a->emat->ProcessQueues();
1103:   /* printf("Finished ProcessQueues\n"); */
1104:   return 0;
1105: }

1107: PetscErrorCode MatAssemblyEnd_Elemental(Mat A, MatAssemblyType type)
1108: {
1109:   /* Currently does nothing */
1110:   return 0;
1111: }

1113: PetscErrorCode MatLoad_Elemental(Mat newMat, PetscViewer viewer)
1114: {
1115:   Mat      Adense, Ae;
1116:   MPI_Comm comm;

1118:   PetscObjectGetComm((PetscObject)newMat, &comm);
1119:   MatCreate(comm, &Adense);
1120:   MatSetType(Adense, MATDENSE);
1121:   MatLoad(Adense, viewer);
1122:   MatConvert(Adense, MATELEMENTAL, MAT_INITIAL_MATRIX, &Ae);
1123:   MatDestroy(&Adense);
1124:   MatHeaderReplace(newMat, &Ae);
1125:   return 0;
1126: }

1128: /* -------------------------------------------------------------------*/
1129: static struct _MatOps MatOps_Values = {MatSetValues_Elemental,
1130:                                        0,
1131:                                        0,
1132:                                        MatMult_Elemental,
1133:                                        /* 4*/ MatMultAdd_Elemental,
1134:                                        MatMultTranspose_Elemental,
1135:                                        MatMultTransposeAdd_Elemental,
1136:                                        MatSolve_Elemental,
1137:                                        MatSolveAdd_Elemental,
1138:                                        0,
1139:                                        /*10*/ 0,
1140:                                        MatLUFactor_Elemental,
1141:                                        MatCholeskyFactor_Elemental,
1142:                                        0,
1143:                                        MatTranspose_Elemental,
1144:                                        /*15*/ MatGetInfo_Elemental,
1145:                                        0,
1146:                                        MatGetDiagonal_Elemental,
1147:                                        MatDiagonalScale_Elemental,
1148:                                        MatNorm_Elemental,
1149:                                        /*20*/ MatAssemblyBegin_Elemental,
1150:                                        MatAssemblyEnd_Elemental,
1151:                                        MatSetOption_Elemental,
1152:                                        MatZeroEntries_Elemental,
1153:                                        /*24*/ 0,
1154:                                        MatLUFactorSymbolic_Elemental,
1155:                                        MatLUFactorNumeric_Elemental,
1156:                                        MatCholeskyFactorSymbolic_Elemental,
1157:                                        MatCholeskyFactorNumeric_Elemental,
1158:                                        /*29*/ MatSetUp_Elemental,
1159:                                        0,
1160:                                        0,
1161:                                        0,
1162:                                        0,
1163:                                        /*34*/ MatDuplicate_Elemental,
1164:                                        0,
1165:                                        0,
1166:                                        0,
1167:                                        0,
1168:                                        /*39*/ MatAXPY_Elemental,
1169:                                        0,
1170:                                        0,
1171:                                        0,
1172:                                        MatCopy_Elemental,
1173:                                        /*44*/ 0,
1174:                                        MatScale_Elemental,
1175:                                        MatShift_Basic,
1176:                                        0,
1177:                                        0,
1178:                                        /*49*/ 0,
1179:                                        0,
1180:                                        0,
1181:                                        0,
1182:                                        0,
1183:                                        /*54*/ 0,
1184:                                        0,
1185:                                        0,
1186:                                        0,
1187:                                        0,
1188:                                        /*59*/ 0,
1189:                                        MatDestroy_Elemental,
1190:                                        MatView_Elemental,
1191:                                        0,
1192:                                        0,
1193:                                        /*64*/ 0,
1194:                                        0,
1195:                                        0,
1196:                                        0,
1197:                                        0,
1198:                                        /*69*/ 0,
1199:                                        0,
1200:                                        MatConvert_Elemental_Dense,
1201:                                        0,
1202:                                        0,
1203:                                        /*74*/ 0,
1204:                                        0,
1205:                                        0,
1206:                                        0,
1207:                                        0,
1208:                                        /*79*/ 0,
1209:                                        0,
1210:                                        0,
1211:                                        0,
1212:                                        MatLoad_Elemental,
1213:                                        /*84*/ 0,
1214:                                        0,
1215:                                        0,
1216:                                        0,
1217:                                        0,
1218:                                        /*89*/ 0,
1219:                                        0,
1220:                                        MatMatMultNumeric_Elemental,
1221:                                        0,
1222:                                        0,
1223:                                        /*94*/ 0,
1224:                                        0,
1225:                                        0,
1226:                                        MatMatTransposeMultNumeric_Elemental,
1227:                                        0,
1228:                                        /*99*/ MatProductSetFromOptions_Elemental,
1229:                                        0,
1230:                                        0,
1231:                                        MatConjugate_Elemental,
1232:                                        0,
1233:                                        /*104*/ 0,
1234:                                        0,
1235:                                        0,
1236:                                        0,
1237:                                        0,
1238:                                        /*109*/ MatMatSolve_Elemental,
1239:                                        0,
1240:                                        0,
1241:                                        0,
1242:                                        MatMissingDiagonal_Elemental,
1243:                                        /*114*/ 0,
1244:                                        0,
1245:                                        0,
1246:                                        0,
1247:                                        0,
1248:                                        /*119*/ 0,
1249:                                        MatHermitianTranspose_Elemental,
1250:                                        0,
1251:                                        0,
1252:                                        0,
1253:                                        /*124*/ 0,
1254:                                        0,
1255:                                        0,
1256:                                        0,
1257:                                        0,
1258:                                        /*129*/ 0,
1259:                                        0,
1260:                                        0,
1261:                                        0,
1262:                                        0,
1263:                                        /*134*/ 0,
1264:                                        0,
1265:                                        0,
1266:                                        0,
1267:                                        0,
1268:                                        0,
1269:                                        /*140*/ 0,
1270:                                        0,
1271:                                        0,
1272:                                        0,
1273:                                        0,
1274:                                        /*145*/ 0,
1275:                                        0,
1276:                                        0,
1277:                                        0,
1278:                                        0,
1279:                                        /*150*/ 0};

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

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

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

1291:   Level: beginner

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

1298: .seealso: `MATDENSE`, `MATSCALAPACK`, `MatGetOwnershipIS()`
1299: M*/

1301: PETSC_EXTERN PetscErrorCode MatCreate_Elemental(Mat A)
1302: {
1303:   Mat_Elemental      *a;
1304:   PetscBool           flg, flg1;
1305:   Mat_Elemental_Grid *commgrid;
1306:   MPI_Comm            icomm;
1307:   PetscInt            optv1;

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

1312:   PetscNew(&a);
1313:   A->data = (void *)a;

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

1318:   /* Grid needs to be shared between multiple Mats on the same communicator, implement by attribute caching on the MPI_Comm */
1319:   if (Petsc_Elemental_keyval == MPI_KEYVAL_INVALID) {
1320:     MPI_Comm_create_keyval(MPI_COMM_NULL_COPY_FN, MPI_COMM_NULL_DELETE_FN, &Petsc_Elemental_keyval, (void *)0);
1321:     PetscCitationsRegister(ElementalCitation, &ElementalCite);
1322:   }
1323:   PetscCommDuplicate(cxxcomm.comm, &icomm, NULL);
1324:   MPI_Comm_get_attr(icomm, Petsc_Elemental_keyval, (void **)&commgrid, (int *)&flg);
1325:   if (!flg) {
1326:     PetscNew(&commgrid);

1328:     PetscOptionsBegin(PetscObjectComm((PetscObject)A), ((PetscObject)A)->prefix, "Elemental Options", "Mat");
1329:     /* displayed default grid sizes (CommSize,1) are set by us arbitrarily until El::Grid() is called */
1330:     PetscOptionsInt("-mat_elemental_grid_height", "Grid Height", "None", El::mpi::Size(cxxcomm), &optv1, &flg1);
1331:     if (flg1) {
1333:       commgrid->grid = new El::Grid(cxxcomm, optv1); /* use user-provided grid height */
1334:     } else {
1335:       commgrid->grid = new El::Grid(cxxcomm); /* use Elemental default grid sizes */
1336:       /* printf("new commgrid->grid = %p\n",commgrid->grid);  -- memory leak revealed by valgrind? */
1337:     }
1338:     commgrid->grid_refct = 1;
1339:     MPI_Comm_set_attr(icomm, Petsc_Elemental_keyval, (void *)commgrid);

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

1344:     PetscOptionsEnd();
1345:   } else {
1346:     commgrid->grid_refct++;
1347:   }
1348:   PetscCommDestroy(&icomm);
1349:   a->grid        = commgrid->grid;
1350:   a->emat        = new El::DistMatrix<PetscElemScalar>(*a->grid);
1351:   a->roworiented = PETSC_TRUE;

1353:   PetscObjectComposeFunction((PetscObject)A, "MatGetOwnershipIS_C", MatGetOwnershipIS_Elemental);
1354:   PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_elemental_mpidense_C", MatProductSetFromOptions_Elemental_MPIDense);
1355:   PetscObjectChangeTypeName((PetscObject)A, MATELEMENTAL);
1356:   return 0;
1357: }