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: }