Actual source code: matelem.cxx
1: #include <petsc/private/petscelemental.h>
3: const char ElementalCitation[] = "@Article{Elemental2012,\n"
4: " author = {Jack Poulson and Bryan Marker and Jeff R. Hammond and Nichols A. Romero and Robert {v}an~{d}e~{G}eijn},\n"
5: " title = {Elemental: A New Framework for Distributed Memory Dense Matrix Computations},\n"
6: " journal = {{ACM} Transactions on Mathematical Software},\n"
7: " volume = {39},\n"
8: " number = {2},\n"
9: " year = {2013}\n"
10: "}\n";
11: static PetscBool ElementalCite = PETSC_FALSE;
13: /*
14: The variable Petsc_Elemental_keyval is used to indicate an MPI attribute that
15: is attached to a communicator, in this case the attribute is a Mat_Elemental_Grid
16: */
17: static PetscMPIInt Petsc_Elemental_keyval = MPI_KEYVAL_INVALID;
19: static PetscErrorCode MatView_Elemental(Mat A, PetscViewer viewer)
20: {
21: Mat_Elemental *a = (Mat_Elemental *)A->data;
22: PetscBool iascii;
24: PetscFunctionBegin;
25: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
26: if (iascii) {
27: PetscViewerFormat format;
28: PetscCall(PetscViewerGetFormat(viewer, &format));
29: if (format == PETSC_VIEWER_ASCII_INFO) {
30: /* call elemental viewing function */
31: PetscCall(PetscViewerASCIIPrintf(viewer, "Elemental run parameters:\n"));
32: PetscCall(PetscViewerASCIIPrintf(viewer, " allocated entries=%zu\n", (*a->emat).AllocatedMemory()));
33: PetscCall(PetscViewerASCIIPrintf(viewer, " grid height=%d, grid width=%d\n", (*a->emat).Grid().Height(), (*a->emat).Grid().Width()));
34: if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) {
35: /* call elemental viewing function */
36: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)viewer), "test matview_elemental 2\n"));
37: }
39: } else if (format == PETSC_VIEWER_DEFAULT) {
40: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
41: El::Print(*a->emat, "Elemental matrix (cyclic ordering)");
42: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
43: if (A->factortype == MAT_FACTOR_NONE) {
44: Mat Adense;
45: PetscCall(MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &Adense));
46: PetscCall(MatView(Adense, viewer));
47: PetscCall(MatDestroy(&Adense));
48: }
49: } else SETERRQ(PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Format");
50: } else {
51: /* convert to dense format and call MatView() */
52: Mat Adense;
53: PetscCall(MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &Adense));
54: PetscCall(MatView(Adense, viewer));
55: PetscCall(MatDestroy(&Adense));
56: }
57: PetscFunctionReturn(PETSC_SUCCESS);
58: }
60: static PetscErrorCode MatGetInfo_Elemental(Mat A, MatInfoType flag, MatInfo *info)
61: {
62: Mat_Elemental *a = (Mat_Elemental *)A->data;
64: PetscFunctionBegin;
65: info->block_size = 1.0;
67: if (flag == MAT_LOCAL) {
68: info->nz_allocated = (*a->emat).AllocatedMemory(); /* locally allocated */
69: info->nz_used = info->nz_allocated;
70: } else if (flag == MAT_GLOBAL_MAX) {
71: //PetscCall(MPIU_Allreduce(isend,irecv,5,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)matin)));
72: /* see MatGetInfo_MPIAIJ() for getting global info->nz_allocated! */
73: //SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP," MAT_GLOBAL_MAX not written yet");
74: } else if (flag == MAT_GLOBAL_SUM) {
75: //SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP," MAT_GLOBAL_SUM not written yet");
76: info->nz_allocated = (*a->emat).AllocatedMemory(); /* locally allocated */
77: info->nz_used = info->nz_allocated; /* assume Elemental does accurate allocation */
78: //PetscCall(MPIU_Allreduce(isend,irecv,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)A)));
79: //PetscPrintf(PETSC_COMM_SELF," ... [%d] locally allocated %g\n",rank,info->nz_allocated);
80: }
82: info->nz_unneeded = 0.0;
83: info->assemblies = A->num_ass;
84: info->mallocs = 0;
85: info->memory = 0; /* REVIEW ME */
86: info->fill_ratio_given = 0; /* determined by Elemental */
87: info->fill_ratio_needed = 0;
88: info->factor_mallocs = 0;
89: PetscFunctionReturn(PETSC_SUCCESS);
90: }
92: static PetscErrorCode MatSetOption_Elemental(Mat A, MatOption op, PetscBool flg)
93: {
94: Mat_Elemental *a = (Mat_Elemental *)A->data;
96: PetscFunctionBegin;
97: switch (op) {
98: case MAT_NEW_NONZERO_LOCATIONS:
99: case MAT_NEW_NONZERO_LOCATION_ERR:
100: case MAT_NEW_NONZERO_ALLOCATION_ERR:
101: case MAT_SYMMETRIC:
102: case MAT_SORTED_FULL:
103: case MAT_HERMITIAN:
104: break;
105: case MAT_ROW_ORIENTED:
106: a->roworiented = flg;
107: break;
108: default:
109: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "unknown option %s", MatOptions[op]);
110: }
111: PetscFunctionReturn(PETSC_SUCCESS);
112: }
114: static PetscErrorCode MatSetValues_Elemental(Mat A, PetscInt nr, const PetscInt *rows, PetscInt nc, const PetscInt *cols, const PetscScalar *vals, InsertMode imode)
115: {
116: Mat_Elemental *a = (Mat_Elemental *)A->data;
117: PetscInt i, j, rrank, ridx, crank, cidx, erow, ecol, numQueues = 0;
119: PetscFunctionBegin;
120: // TODO: Initialize matrix to all zeros?
122: // Count the number of queues from this process
123: if (a->roworiented) {
124: for (i = 0; i < nr; i++) {
125: if (rows[i] < 0) continue;
126: P2RO(A, 0, rows[i], &rrank, &ridx);
127: RO2E(A, 0, rrank, ridx, &erow);
128: PetscCheck(rrank >= 0 && ridx >= 0 && erow >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Incorrect row translation");
129: for (j = 0; j < nc; j++) {
130: if (cols[j] < 0) continue;
131: P2RO(A, 1, cols[j], &crank, &cidx);
132: RO2E(A, 1, crank, cidx, &ecol);
133: PetscCheck(crank >= 0 && cidx >= 0 && ecol >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Incorrect col translation");
134: if (!a->emat->IsLocal(erow, ecol)) { /* off-proc entry */
135: /* printf("Will later remotely update (%d,%d)\n",erow,ecol); */
136: PetscCheck(imode == ADD_VALUES, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only ADD_VALUES to off-processor entry is supported");
137: ++numQueues;
138: continue;
139: }
140: /* printf("Locally updating (%d,%d)\n",erow,ecol); */
141: switch (imode) {
142: case INSERT_VALUES:
143: a->emat->Set(erow, ecol, (PetscElemScalar)vals[i * nc + j]);
144: break;
145: case ADD_VALUES:
146: a->emat->Update(erow, ecol, (PetscElemScalar)vals[i * nc + j]);
147: break;
148: default:
149: SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "No support for InsertMode %d", (int)imode);
150: }
151: }
152: }
154: /* printf("numQueues=%d\n",numQueues); */
155: a->emat->Reserve(numQueues);
156: for (i = 0; i < nr; i++) {
157: if (rows[i] < 0) continue;
158: P2RO(A, 0, rows[i], &rrank, &ridx);
159: RO2E(A, 0, rrank, ridx, &erow);
160: for (j = 0; j < nc; j++) {
161: if (cols[j] < 0) continue;
162: P2RO(A, 1, cols[j], &crank, &cidx);
163: RO2E(A, 1, crank, cidx, &ecol);
164: if (!a->emat->IsLocal(erow, ecol)) { /*off-proc entry*/
165: /* printf("Queueing remotely update of (%d,%d)\n",erow,ecol); */
166: a->emat->QueueUpdate(erow, ecol, vals[i * nc + j]);
167: }
168: }
169: }
170: } else { /* column-oriented */
171: for (j = 0; j < nc; j++) {
172: if (cols[j] < 0) continue;
173: P2RO(A, 1, cols[j], &crank, &cidx);
174: RO2E(A, 1, crank, cidx, &ecol);
175: PetscCheck(crank >= 0 && cidx >= 0 && ecol >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Incorrect col translation");
176: for (i = 0; i < nr; i++) {
177: if (rows[i] < 0) continue;
178: P2RO(A, 0, rows[i], &rrank, &ridx);
179: RO2E(A, 0, rrank, ridx, &erow);
180: PetscCheck(rrank >= 0 && ridx >= 0 && erow >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Incorrect row translation");
181: if (!a->emat->IsLocal(erow, ecol)) { /* off-proc entry */
182: /* printf("Will later remotely update (%d,%d)\n",erow,ecol); */
183: PetscCheck(imode == ADD_VALUES, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only ADD_VALUES to off-processor entry is supported");
184: ++numQueues;
185: continue;
186: }
187: /* printf("Locally updating (%d,%d)\n",erow,ecol); */
188: switch (imode) {
189: case INSERT_VALUES:
190: a->emat->Set(erow, ecol, (PetscElemScalar)vals[i + j * nr]);
191: break;
192: case ADD_VALUES:
193: a->emat->Update(erow, ecol, (PetscElemScalar)vals[i + j * nr]);
194: break;
195: default:
196: SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "No support for InsertMode %d", (int)imode);
197: }
198: }
199: }
201: /* printf("numQueues=%d\n",numQueues); */
202: a->emat->Reserve(numQueues);
203: for (j = 0; j < nc; j++) {
204: if (cols[j] < 0) continue;
205: P2RO(A, 1, cols[j], &crank, &cidx);
206: RO2E(A, 1, crank, cidx, &ecol);
208: for (i = 0; i < nr; i++) {
209: if (rows[i] < 0) continue;
210: P2RO(A, 0, rows[i], &rrank, &ridx);
211: RO2E(A, 0, rrank, ridx, &erow);
212: if (!a->emat->IsLocal(erow, ecol)) { /*off-proc entry*/
213: /* printf("Queueing remotely update of (%d,%d)\n",erow,ecol); */
214: a->emat->QueueUpdate(erow, ecol, vals[i + j * nr]);
215: }
216: }
217: }
218: }
219: PetscFunctionReturn(PETSC_SUCCESS);
220: }
222: static PetscErrorCode MatMult_Elemental(Mat A, Vec X, Vec Y)
223: {
224: Mat_Elemental *a = (Mat_Elemental *)A->data;
225: const PetscElemScalar *x;
226: PetscElemScalar *y;
227: PetscElemScalar one = 1, zero = 0;
229: PetscFunctionBegin;
230: PetscCall(VecGetArrayRead(X, (const PetscScalar **)&x));
231: PetscCall(VecGetArray(Y, (PetscScalar **)&y));
232: { /* Scoping so that constructor is called before pointer is returned */
233: El::DistMatrix<PetscElemScalar, El::VC, El::STAR> xe, ye;
234: xe.LockedAttach(A->cmap->N, 1, *a->grid, 0, 0, x, A->cmap->n);
235: ye.Attach(A->rmap->N, 1, *a->grid, 0, 0, y, A->rmap->n);
236: El::Gemv(El::NORMAL, one, *a->emat, xe, zero, ye);
237: }
238: PetscCall(VecRestoreArrayRead(X, (const PetscScalar **)&x));
239: PetscCall(VecRestoreArray(Y, (PetscScalar **)&y));
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: static PetscErrorCode MatMultTranspose_Elemental(Mat A, Vec X, Vec Y)
244: {
245: Mat_Elemental *a = (Mat_Elemental *)A->data;
246: const PetscElemScalar *x;
247: PetscElemScalar *y;
248: PetscElemScalar one = 1, zero = 0;
250: PetscFunctionBegin;
251: PetscCall(VecGetArrayRead(X, (const PetscScalar **)&x));
252: PetscCall(VecGetArray(Y, (PetscScalar **)&y));
253: { /* Scoping so that constructor is called before pointer is returned */
254: El::DistMatrix<PetscElemScalar, El::VC, El::STAR> xe, ye;
255: xe.LockedAttach(A->rmap->N, 1, *a->grid, 0, 0, x, A->rmap->n);
256: ye.Attach(A->cmap->N, 1, *a->grid, 0, 0, y, A->cmap->n);
257: El::Gemv(El::TRANSPOSE, one, *a->emat, xe, zero, ye);
258: }
259: PetscCall(VecRestoreArrayRead(X, (const PetscScalar **)&x));
260: PetscCall(VecRestoreArray(Y, (PetscScalar **)&y));
261: PetscFunctionReturn(PETSC_SUCCESS);
262: }
264: static PetscErrorCode MatMultAdd_Elemental(Mat A, Vec X, Vec Y, Vec Z)
265: {
266: Mat_Elemental *a = (Mat_Elemental *)A->data;
267: const PetscElemScalar *x;
268: PetscElemScalar *z;
269: PetscElemScalar one = 1;
271: PetscFunctionBegin;
272: if (Y != Z) PetscCall(VecCopy(Y, Z));
273: PetscCall(VecGetArrayRead(X, (const PetscScalar **)&x));
274: PetscCall(VecGetArray(Z, (PetscScalar **)&z));
275: { /* Scoping so that constructor is called before pointer is returned */
276: El::DistMatrix<PetscElemScalar, El::VC, El::STAR> xe, ze;
277: xe.LockedAttach(A->cmap->N, 1, *a->grid, 0, 0, x, A->cmap->n);
278: ze.Attach(A->rmap->N, 1, *a->grid, 0, 0, z, A->rmap->n);
279: El::Gemv(El::NORMAL, one, *a->emat, xe, one, ze);
280: }
281: PetscCall(VecRestoreArrayRead(X, (const PetscScalar **)&x));
282: PetscCall(VecRestoreArray(Z, (PetscScalar **)&z));
283: PetscFunctionReturn(PETSC_SUCCESS);
284: }
286: static PetscErrorCode MatMultTransposeAdd_Elemental(Mat A, Vec X, Vec Y, Vec Z)
287: {
288: Mat_Elemental *a = (Mat_Elemental *)A->data;
289: const PetscElemScalar *x;
290: PetscElemScalar *z;
291: PetscElemScalar one = 1;
293: PetscFunctionBegin;
294: if (Y != Z) PetscCall(VecCopy(Y, Z));
295: PetscCall(VecGetArrayRead(X, (const PetscScalar **)&x));
296: PetscCall(VecGetArray(Z, (PetscScalar **)&z));
297: { /* Scoping so that constructor is called before pointer is returned */
298: El::DistMatrix<PetscElemScalar, El::VC, El::STAR> xe, ze;
299: xe.LockedAttach(A->rmap->N, 1, *a->grid, 0, 0, x, A->rmap->n);
300: ze.Attach(A->cmap->N, 1, *a->grid, 0, 0, z, A->cmap->n);
301: El::Gemv(El::TRANSPOSE, one, *a->emat, xe, one, ze);
302: }
303: PetscCall(VecRestoreArrayRead(X, (const PetscScalar **)&x));
304: PetscCall(VecRestoreArray(Z, (PetscScalar **)&z));
305: PetscFunctionReturn(PETSC_SUCCESS);
306: }
308: PetscErrorCode MatMatMultNumeric_Elemental(Mat A, Mat B, Mat C)
309: {
310: Mat_Elemental *a = (Mat_Elemental *)A->data;
311: Mat_Elemental *b = (Mat_Elemental *)B->data;
312: Mat_Elemental *c = (Mat_Elemental *)C->data;
313: PetscElemScalar one = 1, zero = 0;
315: PetscFunctionBegin;
316: { /* Scoping so that constructor is called before pointer is returned */
317: El::Gemm(El::NORMAL, El::NORMAL, one, *a->emat, *b->emat, zero, *c->emat);
318: }
319: C->assembled = PETSC_TRUE;
320: PetscFunctionReturn(PETSC_SUCCESS);
321: }
323: PetscErrorCode MatMatMultSymbolic_Elemental(Mat A, Mat B, PetscReal, Mat Ce)
324: {
325: PetscFunctionBegin;
326: PetscCall(MatSetSizes(Ce, A->rmap->n, B->cmap->n, PETSC_DECIDE, PETSC_DECIDE));
327: PetscCall(MatSetType(Ce, MATELEMENTAL));
328: PetscCall(MatSetUp(Ce));
329: Ce->ops->matmultnumeric = MatMatMultNumeric_Elemental;
330: PetscFunctionReturn(PETSC_SUCCESS);
331: }
333: static PetscErrorCode MatMatTransposeMultNumeric_Elemental(Mat A, Mat B, Mat C)
334: {
335: Mat_Elemental *a = (Mat_Elemental *)A->data;
336: Mat_Elemental *b = (Mat_Elemental *)B->data;
337: Mat_Elemental *c = (Mat_Elemental *)C->data;
338: PetscElemScalar one = 1, zero = 0;
340: PetscFunctionBegin;
341: { /* Scoping so that constructor is called before pointer is returned */
342: El::Gemm(El::NORMAL, El::TRANSPOSE, one, *a->emat, *b->emat, zero, *c->emat);
343: }
344: C->assembled = PETSC_TRUE;
345: PetscFunctionReturn(PETSC_SUCCESS);
346: }
348: static PetscErrorCode MatMatTransposeMultSymbolic_Elemental(Mat A, Mat B, PetscReal, Mat C)
349: {
350: PetscFunctionBegin;
351: PetscCall(MatSetSizes(C, A->rmap->n, B->rmap->n, PETSC_DECIDE, PETSC_DECIDE));
352: PetscCall(MatSetType(C, MATELEMENTAL));
353: PetscCall(MatSetUp(C));
354: PetscFunctionReturn(PETSC_SUCCESS);
355: }
357: static PetscErrorCode MatProductSetFromOptions_Elemental_AB(Mat C)
358: {
359: PetscFunctionBegin;
360: C->ops->matmultsymbolic = MatMatMultSymbolic_Elemental;
361: C->ops->productsymbolic = MatProductSymbolic_AB;
362: PetscFunctionReturn(PETSC_SUCCESS);
363: }
365: static PetscErrorCode MatProductSetFromOptions_Elemental_ABt(Mat C)
366: {
367: PetscFunctionBegin;
368: C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_Elemental;
369: C->ops->productsymbolic = MatProductSymbolic_ABt;
370: PetscFunctionReturn(PETSC_SUCCESS);
371: }
373: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Elemental(Mat C)
374: {
375: Mat_Product *product = C->product;
377: PetscFunctionBegin;
378: switch (product->type) {
379: case MATPRODUCT_AB:
380: PetscCall(MatProductSetFromOptions_Elemental_AB(C));
381: break;
382: case MATPRODUCT_ABt:
383: PetscCall(MatProductSetFromOptions_Elemental_ABt(C));
384: break;
385: default:
386: break;
387: }
388: PetscFunctionReturn(PETSC_SUCCESS);
389: }
391: static PetscErrorCode MatMatMultNumeric_Elemental_MPIDense(Mat A, Mat B, Mat C)
392: {
393: Mat Be, Ce;
395: PetscFunctionBegin;
396: PetscCall(MatConvert(B, MATELEMENTAL, MAT_INITIAL_MATRIX, &Be));
397: PetscCall(MatMatMult(A, Be, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Ce));
398: PetscCall(MatConvert(Ce, MATMPIDENSE, MAT_REUSE_MATRIX, &C));
399: PetscCall(MatDestroy(&Be));
400: PetscCall(MatDestroy(&Ce));
401: PetscFunctionReturn(PETSC_SUCCESS);
402: }
404: static PetscErrorCode MatMatMultSymbolic_Elemental_MPIDense(Mat A, Mat B, PetscReal, Mat C)
405: {
406: PetscFunctionBegin;
407: PetscCall(MatSetSizes(C, A->rmap->n, B->cmap->n, PETSC_DECIDE, PETSC_DECIDE));
408: PetscCall(MatSetType(C, MATMPIDENSE));
409: PetscCall(MatSetUp(C));
410: C->ops->matmultnumeric = MatMatMultNumeric_Elemental_MPIDense;
411: PetscFunctionReturn(PETSC_SUCCESS);
412: }
414: static PetscErrorCode MatProductSetFromOptions_Elemental_MPIDense_AB(Mat C)
415: {
416: PetscFunctionBegin;
417: C->ops->matmultsymbolic = MatMatMultSymbolic_Elemental_MPIDense;
418: C->ops->productsymbolic = MatProductSymbolic_AB;
419: PetscFunctionReturn(PETSC_SUCCESS);
420: }
422: static PetscErrorCode MatProductSetFromOptions_Elemental_MPIDense(Mat C)
423: {
424: Mat_Product *product = C->product;
426: PetscFunctionBegin;
427: if (product->type == MATPRODUCT_AB) PetscCall(MatProductSetFromOptions_Elemental_MPIDense_AB(C));
428: PetscFunctionReturn(PETSC_SUCCESS);
429: }
431: static PetscErrorCode MatGetDiagonal_Elemental(Mat A, Vec D)
432: {
433: PetscInt i, nrows, ncols, nD, rrank, ridx, crank, cidx;
434: Mat_Elemental *a = (Mat_Elemental *)A->data;
435: PetscElemScalar v;
436: MPI_Comm comm;
438: PetscFunctionBegin;
439: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
440: PetscCall(MatGetSize(A, &nrows, &ncols));
441: nD = nrows > ncols ? ncols : nrows;
442: for (i = 0; i < nD; i++) {
443: PetscInt erow, ecol;
444: P2RO(A, 0, i, &rrank, &ridx);
445: RO2E(A, 0, rrank, ridx, &erow);
446: PetscCheck(rrank >= 0 && ridx >= 0 && erow >= 0, comm, PETSC_ERR_PLIB, "Incorrect row translation");
447: P2RO(A, 1, i, &crank, &cidx);
448: RO2E(A, 1, crank, cidx, &ecol);
449: PetscCheck(crank >= 0 && cidx >= 0 && ecol >= 0, comm, PETSC_ERR_PLIB, "Incorrect col translation");
450: v = a->emat->Get(erow, ecol);
451: PetscCall(VecSetValues(D, 1, &i, (PetscScalar *)&v, INSERT_VALUES));
452: }
453: PetscCall(VecAssemblyBegin(D));
454: PetscCall(VecAssemblyEnd(D));
455: PetscFunctionReturn(PETSC_SUCCESS);
456: }
458: static PetscErrorCode MatDiagonalScale_Elemental(Mat X, Vec L, Vec R)
459: {
460: Mat_Elemental *x = (Mat_Elemental *)X->data;
461: const PetscElemScalar *d;
463: PetscFunctionBegin;
464: if (R) {
465: PetscCall(VecGetArrayRead(R, (const PetscScalar **)&d));
466: El::DistMatrix<PetscElemScalar, El::VC, El::STAR> de;
467: de.LockedAttach(X->cmap->N, 1, *x->grid, 0, 0, d, X->cmap->n);
468: El::DiagonalScale(El::RIGHT, El::NORMAL, de, *x->emat);
469: PetscCall(VecRestoreArrayRead(R, (const PetscScalar **)&d));
470: }
471: if (L) {
472: PetscCall(VecGetArrayRead(L, (const PetscScalar **)&d));
473: El::DistMatrix<PetscElemScalar, El::VC, El::STAR> de;
474: de.LockedAttach(X->rmap->N, 1, *x->grid, 0, 0, d, X->rmap->n);
475: El::DiagonalScale(El::LEFT, El::NORMAL, de, *x->emat);
476: PetscCall(VecRestoreArrayRead(L, (const PetscScalar **)&d));
477: }
478: PetscFunctionReturn(PETSC_SUCCESS);
479: }
481: static PetscErrorCode MatMissingDiagonal_Elemental(Mat, PetscBool *missing, PetscInt *)
482: {
483: PetscFunctionBegin;
484: *missing = PETSC_FALSE;
485: PetscFunctionReturn(PETSC_SUCCESS);
486: }
488: static PetscErrorCode MatScale_Elemental(Mat X, PetscScalar a)
489: {
490: Mat_Elemental *x = (Mat_Elemental *)X->data;
492: PetscFunctionBegin;
493: El::Scale((PetscElemScalar)a, *x->emat);
494: PetscFunctionReturn(PETSC_SUCCESS);
495: }
497: /*
498: MatAXPY - Computes Y = a*X + Y.
499: */
500: static PetscErrorCode MatAXPY_Elemental(Mat Y, PetscScalar a, Mat X, MatStructure)
501: {
502: Mat_Elemental *x = (Mat_Elemental *)X->data;
503: Mat_Elemental *y = (Mat_Elemental *)Y->data;
505: PetscFunctionBegin;
506: El::Axpy((PetscElemScalar)a, *x->emat, *y->emat);
507: PetscCall(PetscObjectStateIncrease((PetscObject)Y));
508: PetscFunctionReturn(PETSC_SUCCESS);
509: }
511: static PetscErrorCode MatCopy_Elemental(Mat A, Mat B, MatStructure)
512: {
513: Mat_Elemental *a = (Mat_Elemental *)A->data;
514: Mat_Elemental *b = (Mat_Elemental *)B->data;
516: PetscFunctionBegin;
517: El::Copy(*a->emat, *b->emat);
518: PetscCall(PetscObjectStateIncrease((PetscObject)B));
519: PetscFunctionReturn(PETSC_SUCCESS);
520: }
522: static PetscErrorCode MatDuplicate_Elemental(Mat A, MatDuplicateOption op, Mat *B)
523: {
524: Mat Be;
525: MPI_Comm comm;
526: Mat_Elemental *a = (Mat_Elemental *)A->data;
528: PetscFunctionBegin;
529: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
530: PetscCall(MatCreate(comm, &Be));
531: PetscCall(MatSetSizes(Be, A->rmap->n, A->cmap->n, PETSC_DECIDE, PETSC_DECIDE));
532: PetscCall(MatSetType(Be, MATELEMENTAL));
533: PetscCall(MatSetUp(Be));
534: *B = Be;
535: if (op == MAT_COPY_VALUES) {
536: Mat_Elemental *b = (Mat_Elemental *)Be->data;
537: El::Copy(*a->emat, *b->emat);
538: }
539: Be->assembled = PETSC_TRUE;
540: PetscFunctionReturn(PETSC_SUCCESS);
541: }
543: static PetscErrorCode MatTranspose_Elemental(Mat A, MatReuse reuse, Mat *B)
544: {
545: Mat Be = *B;
546: MPI_Comm comm;
547: Mat_Elemental *a = (Mat_Elemental *)A->data, *b;
549: PetscFunctionBegin;
550: if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *B));
551: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
552: /* Only out-of-place supported */
553: PetscCheck(reuse != MAT_INPLACE_MATRIX, comm, PETSC_ERR_SUP, "Only out-of-place supported");
554: if (reuse == MAT_INITIAL_MATRIX) {
555: PetscCall(MatCreate(comm, &Be));
556: PetscCall(MatSetSizes(Be, A->cmap->n, A->rmap->n, PETSC_DECIDE, PETSC_DECIDE));
557: PetscCall(MatSetType(Be, MATELEMENTAL));
558: PetscCall(MatSetUp(Be));
559: *B = Be;
560: }
561: b = (Mat_Elemental *)Be->data;
562: El::Transpose(*a->emat, *b->emat);
563: Be->assembled = PETSC_TRUE;
564: PetscFunctionReturn(PETSC_SUCCESS);
565: }
567: static PetscErrorCode MatConjugate_Elemental(Mat A)
568: {
569: Mat_Elemental *a = (Mat_Elemental *)A->data;
571: PetscFunctionBegin;
572: El::Conjugate(*a->emat);
573: PetscFunctionReturn(PETSC_SUCCESS);
574: }
576: static PetscErrorCode MatHermitianTranspose_Elemental(Mat A, MatReuse reuse, Mat *B)
577: {
578: Mat Be = *B;
579: MPI_Comm comm;
580: Mat_Elemental *a = (Mat_Elemental *)A->data, *b;
582: PetscFunctionBegin;
583: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
584: /* Only out-of-place supported */
585: if (reuse == MAT_INITIAL_MATRIX) {
586: PetscCall(MatCreate(comm, &Be));
587: PetscCall(MatSetSizes(Be, A->cmap->n, A->rmap->n, PETSC_DECIDE, PETSC_DECIDE));
588: PetscCall(MatSetType(Be, MATELEMENTAL));
589: PetscCall(MatSetUp(Be));
590: *B = Be;
591: }
592: b = (Mat_Elemental *)Be->data;
593: El::Adjoint(*a->emat, *b->emat);
594: Be->assembled = PETSC_TRUE;
595: PetscFunctionReturn(PETSC_SUCCESS);
596: }
598: static PetscErrorCode MatSolve_Elemental(Mat A, Vec B, Vec X)
599: {
600: Mat_Elemental *a = (Mat_Elemental *)A->data;
601: PetscElemScalar *x;
602: PetscInt pivoting = a->pivoting;
604: PetscFunctionBegin;
605: PetscCall(VecCopy(B, X));
606: PetscCall(VecGetArray(X, (PetscScalar **)&x));
608: El::DistMatrix<PetscElemScalar, El::VC, El::STAR> xe;
609: xe.Attach(A->rmap->N, 1, *a->grid, 0, 0, x, A->rmap->n);
610: El::DistMatrix<PetscElemScalar, El::MC, El::MR> xer(xe);
611: switch (A->factortype) {
612: case MAT_FACTOR_LU:
613: if (pivoting == 0) {
614: El::lu::SolveAfter(El::NORMAL, *a->emat, xer);
615: } else if (pivoting == 1) {
616: El::lu::SolveAfter(El::NORMAL, *a->emat, *a->P, xer);
617: } else { /* pivoting == 2 */
618: El::lu::SolveAfter(El::NORMAL, *a->emat, *a->P, *a->Q, xer);
619: }
620: break;
621: case MAT_FACTOR_CHOLESKY:
622: El::cholesky::SolveAfter(El::UPPER, El::NORMAL, *a->emat, xer);
623: break;
624: default:
625: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unfactored Matrix or Unsupported MatFactorType");
626: break;
627: }
628: El::Copy(xer, xe);
630: PetscCall(VecRestoreArray(X, (PetscScalar **)&x));
631: PetscFunctionReturn(PETSC_SUCCESS);
632: }
634: static PetscErrorCode MatSolveAdd_Elemental(Mat A, Vec B, Vec Y, Vec X)
635: {
636: PetscFunctionBegin;
637: PetscCall(MatSolve_Elemental(A, B, X));
638: PetscCall(VecAXPY(X, 1, Y));
639: PetscFunctionReturn(PETSC_SUCCESS);
640: }
642: static PetscErrorCode MatMatSolve_Elemental(Mat A, Mat B, Mat X)
643: {
644: Mat_Elemental *a = (Mat_Elemental *)A->data;
645: Mat_Elemental *x;
646: Mat C;
647: PetscInt pivoting = a->pivoting;
648: PetscBool flg;
649: MatType type;
651: PetscFunctionBegin;
652: PetscCall(MatGetType(X, &type));
653: PetscCall(PetscStrcmp(type, MATELEMENTAL, &flg));
654: if (!flg) {
655: PetscCall(MatConvert(B, MATELEMENTAL, MAT_INITIAL_MATRIX, &C));
656: x = (Mat_Elemental *)C->data;
657: } else {
658: x = (Mat_Elemental *)X->data;
659: El::Copy(*((Mat_Elemental *)B->data)->emat, *x->emat);
660: }
661: switch (A->factortype) {
662: case MAT_FACTOR_LU:
663: if (pivoting == 0) {
664: El::lu::SolveAfter(El::NORMAL, *a->emat, *x->emat);
665: } else if (pivoting == 1) {
666: El::lu::SolveAfter(El::NORMAL, *a->emat, *a->P, *x->emat);
667: } else {
668: El::lu::SolveAfter(El::NORMAL, *a->emat, *a->P, *a->Q, *x->emat);
669: }
670: break;
671: case MAT_FACTOR_CHOLESKY:
672: El::cholesky::SolveAfter(El::UPPER, El::NORMAL, *a->emat, *x->emat);
673: break;
674: default:
675: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unfactored Matrix or Unsupported MatFactorType");
676: break;
677: }
678: if (!flg) {
679: PetscCall(MatConvert(C, type, MAT_REUSE_MATRIX, &X));
680: PetscCall(MatDestroy(&C));
681: }
682: PetscFunctionReturn(PETSC_SUCCESS);
683: }
685: static PetscErrorCode MatLUFactor_Elemental(Mat A, IS, IS, const MatFactorInfo *)
686: {
687: Mat_Elemental *a = (Mat_Elemental *)A->data;
688: PetscInt pivoting = a->pivoting;
690: PetscFunctionBegin;
691: if (pivoting == 0) {
692: El::LU(*a->emat);
693: } else if (pivoting == 1) {
694: El::LU(*a->emat, *a->P);
695: } else {
696: El::LU(*a->emat, *a->P, *a->Q);
697: }
698: A->factortype = MAT_FACTOR_LU;
699: A->assembled = PETSC_TRUE;
701: PetscCall(PetscFree(A->solvertype));
702: PetscCall(PetscStrallocpy(MATSOLVERELEMENTAL, &A->solvertype));
703: PetscFunctionReturn(PETSC_SUCCESS);
704: }
706: static PetscErrorCode MatLUFactorNumeric_Elemental(Mat F, Mat A, const MatFactorInfo *info)
707: {
708: PetscFunctionBegin;
709: PetscCall(MatCopy(A, F, SAME_NONZERO_PATTERN));
710: PetscCall(MatLUFactor_Elemental(F, nullptr, nullptr, info));
711: PetscFunctionReturn(PETSC_SUCCESS);
712: }
714: static PetscErrorCode MatLUFactorSymbolic_Elemental(Mat, Mat, IS, IS, const MatFactorInfo *)
715: {
716: PetscFunctionBegin;
717: /* F is created and allocated by MatGetFactor_elemental_petsc(), skip this routine. */
718: PetscFunctionReturn(PETSC_SUCCESS);
719: }
721: static PetscErrorCode MatCholeskyFactor_Elemental(Mat A, IS, const MatFactorInfo *)
722: {
723: Mat_Elemental *a = (Mat_Elemental *)A->data;
724: El::DistMatrix<PetscElemScalar, El::MC, El::STAR> d;
726: PetscFunctionBegin;
727: El::Cholesky(El::UPPER, *a->emat);
728: A->factortype = MAT_FACTOR_CHOLESKY;
729: A->assembled = PETSC_TRUE;
731: PetscCall(PetscFree(A->solvertype));
732: PetscCall(PetscStrallocpy(MATSOLVERELEMENTAL, &A->solvertype));
733: PetscFunctionReturn(PETSC_SUCCESS);
734: }
736: static PetscErrorCode MatCholeskyFactorNumeric_Elemental(Mat F, Mat A, const MatFactorInfo *info)
737: {
738: PetscFunctionBegin;
739: PetscCall(MatCopy(A, F, SAME_NONZERO_PATTERN));
740: PetscCall(MatCholeskyFactor_Elemental(F, nullptr, info));
741: PetscFunctionReturn(PETSC_SUCCESS);
742: }
744: static PetscErrorCode MatCholeskyFactorSymbolic_Elemental(Mat, Mat, IS, const MatFactorInfo *)
745: {
746: PetscFunctionBegin;
747: /* F is created and allocated by MatGetFactor_elemental_petsc(), skip this routine. */
748: PetscFunctionReturn(PETSC_SUCCESS);
749: }
751: static PetscErrorCode MatFactorGetSolverType_elemental_elemental(Mat, MatSolverType *type)
752: {
753: PetscFunctionBegin;
754: *type = MATSOLVERELEMENTAL;
755: PetscFunctionReturn(PETSC_SUCCESS);
756: }
758: static PetscErrorCode MatGetFactor_elemental_elemental(Mat A, MatFactorType ftype, Mat *F)
759: {
760: Mat B;
762: PetscFunctionBegin;
763: /* Create the factorization matrix */
764: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
765: PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, PETSC_DECIDE, PETSC_DECIDE));
766: PetscCall(MatSetType(B, MATELEMENTAL));
767: PetscCall(MatSetUp(B));
768: B->factortype = ftype;
769: B->trivialsymbolic = PETSC_TRUE;
770: PetscCall(PetscFree(B->solvertype));
771: PetscCall(PetscStrallocpy(MATSOLVERELEMENTAL, &B->solvertype));
773: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_elemental_elemental));
774: *F = B;
775: PetscFunctionReturn(PETSC_SUCCESS);
776: }
778: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_Elemental(void)
779: {
780: PetscFunctionBegin;
781: PetscCall(MatSolverTypeRegister(MATSOLVERELEMENTAL, MATELEMENTAL, MAT_FACTOR_LU, MatGetFactor_elemental_elemental));
782: PetscCall(MatSolverTypeRegister(MATSOLVERELEMENTAL, MATELEMENTAL, MAT_FACTOR_CHOLESKY, MatGetFactor_elemental_elemental));
783: PetscFunctionReturn(PETSC_SUCCESS);
784: }
786: static PetscErrorCode MatNorm_Elemental(Mat A, NormType type, PetscReal *nrm)
787: {
788: Mat_Elemental *a = (Mat_Elemental *)A->data;
790: PetscFunctionBegin;
791: switch (type) {
792: case NORM_1:
793: *nrm = El::OneNorm(*a->emat);
794: break;
795: case NORM_FROBENIUS:
796: *nrm = El::FrobeniusNorm(*a->emat);
797: break;
798: case NORM_INFINITY:
799: *nrm = El::InfinityNorm(*a->emat);
800: break;
801: default:
802: SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Unsupported norm type");
803: }
804: PetscFunctionReturn(PETSC_SUCCESS);
805: }
807: static PetscErrorCode MatZeroEntries_Elemental(Mat A)
808: {
809: Mat_Elemental *a = (Mat_Elemental *)A->data;
811: PetscFunctionBegin;
812: El::Zero(*a->emat);
813: PetscFunctionReturn(PETSC_SUCCESS);
814: }
816: static PetscErrorCode MatGetOwnershipIS_Elemental(Mat A, IS *rows, IS *cols)
817: {
818: Mat_Elemental *a = (Mat_Elemental *)A->data;
819: PetscInt i, m, shift, stride, *idx;
821: PetscFunctionBegin;
822: if (rows) {
823: m = a->emat->LocalHeight();
824: shift = a->emat->ColShift();
825: stride = a->emat->ColStride();
826: PetscCall(PetscMalloc1(m, &idx));
827: for (i = 0; i < m; i++) {
828: PetscInt rank, offset;
829: E2RO(A, 0, shift + i * stride, &rank, &offset);
830: RO2P(A, 0, rank, offset, &idx[i]);
831: }
832: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, m, idx, PETSC_OWN_POINTER, rows));
833: }
834: if (cols) {
835: m = a->emat->LocalWidth();
836: shift = a->emat->RowShift();
837: stride = a->emat->RowStride();
838: PetscCall(PetscMalloc1(m, &idx));
839: for (i = 0; i < m; i++) {
840: PetscInt rank, offset;
841: E2RO(A, 1, shift + i * stride, &rank, &offset);
842: RO2P(A, 1, rank, offset, &idx[i]);
843: }
844: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, m, idx, PETSC_OWN_POINTER, cols));
845: }
846: PetscFunctionReturn(PETSC_SUCCESS);
847: }
849: static PetscErrorCode MatConvert_Elemental_Dense(Mat A, MatType, MatReuse reuse, Mat *B)
850: {
851: Mat Bmpi;
852: Mat_Elemental *a = (Mat_Elemental *)A->data;
853: MPI_Comm comm;
854: IS isrows, iscols;
855: PetscInt rrank, ridx, crank, cidx, nrows, ncols, i, j, erow, ecol, elrow, elcol;
856: const PetscInt *rows, *cols;
857: PetscElemScalar v;
858: const El::Grid &grid = a->emat->Grid();
860: PetscFunctionBegin;
861: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
863: if (reuse == MAT_REUSE_MATRIX) {
864: Bmpi = *B;
865: } else {
866: PetscCall(MatCreate(comm, &Bmpi));
867: PetscCall(MatSetSizes(Bmpi, A->rmap->n, A->cmap->n, PETSC_DECIDE, PETSC_DECIDE));
868: PetscCall(MatSetType(Bmpi, MATDENSE));
869: PetscCall(MatSetUp(Bmpi));
870: }
872: /* Get local entries of A */
873: PetscCall(MatGetOwnershipIS(A, &isrows, &iscols));
874: PetscCall(ISGetLocalSize(isrows, &nrows));
875: PetscCall(ISGetIndices(isrows, &rows));
876: PetscCall(ISGetLocalSize(iscols, &ncols));
877: PetscCall(ISGetIndices(iscols, &cols));
879: if (a->roworiented) {
880: for (i = 0; i < nrows; i++) {
881: P2RO(A, 0, rows[i], &rrank, &ridx); /* convert indices between PETSc <-> (Rank,Offset) <-> Elemental */
882: RO2E(A, 0, rrank, ridx, &erow);
883: PetscCheck(rrank >= 0 && ridx >= 0 && erow >= 0, comm, PETSC_ERR_PLIB, "Incorrect row translation");
884: for (j = 0; j < ncols; j++) {
885: P2RO(A, 1, cols[j], &crank, &cidx);
886: RO2E(A, 1, crank, cidx, &ecol);
887: PetscCheck(crank >= 0 && cidx >= 0 && ecol >= 0, comm, PETSC_ERR_PLIB, "Incorrect col translation");
889: elrow = erow / grid.MCSize(); /* Elemental local row index */
890: elcol = ecol / grid.MRSize(); /* Elemental local column index */
891: v = a->emat->GetLocal(elrow, elcol);
892: PetscCall(MatSetValues(Bmpi, 1, &rows[i], 1, &cols[j], (PetscScalar *)&v, INSERT_VALUES));
893: }
894: }
895: } else { /* column-oriented */
896: for (j = 0; j < ncols; j++) {
897: P2RO(A, 1, cols[j], &crank, &cidx);
898: RO2E(A, 1, crank, cidx, &ecol);
899: PetscCheck(crank >= 0 && cidx >= 0 && ecol >= 0, comm, PETSC_ERR_PLIB, "Incorrect col translation");
900: for (i = 0; i < nrows; i++) {
901: P2RO(A, 0, rows[i], &rrank, &ridx); /* convert indices between PETSc <-> (Rank,Offset) <-> Elemental */
902: RO2E(A, 0, rrank, ridx, &erow);
903: PetscCheck(rrank >= 0 && ridx >= 0 && erow >= 0, comm, PETSC_ERR_PLIB, "Incorrect row translation");
905: elrow = erow / grid.MCSize(); /* Elemental local row index */
906: elcol = ecol / grid.MRSize(); /* Elemental local column index */
907: v = a->emat->GetLocal(elrow, elcol);
908: PetscCall(MatSetValues(Bmpi, 1, &rows[i], 1, &cols[j], (PetscScalar *)&v, INSERT_VALUES));
909: }
910: }
911: }
912: PetscCall(MatAssemblyBegin(Bmpi, MAT_FINAL_ASSEMBLY));
913: PetscCall(MatAssemblyEnd(Bmpi, MAT_FINAL_ASSEMBLY));
914: if (reuse == MAT_INPLACE_MATRIX) {
915: PetscCall(MatHeaderReplace(A, &Bmpi));
916: } else {
917: *B = Bmpi;
918: }
919: PetscCall(ISDestroy(&isrows));
920: PetscCall(ISDestroy(&iscols));
921: PetscFunctionReturn(PETSC_SUCCESS);
922: }
924: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat A, MatType, MatReuse reuse, Mat *newmat)
925: {
926: Mat mat_elemental;
927: PetscInt M = A->rmap->N, N = A->cmap->N, row, ncols;
928: const PetscInt *cols;
929: const PetscScalar *vals;
931: PetscFunctionBegin;
932: if (reuse == MAT_REUSE_MATRIX) {
933: mat_elemental = *newmat;
934: PetscCall(MatZeroEntries(mat_elemental));
935: } else {
936: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental));
937: PetscCall(MatSetSizes(mat_elemental, PETSC_DECIDE, PETSC_DECIDE, M, N));
938: PetscCall(MatSetType(mat_elemental, MATELEMENTAL));
939: PetscCall(MatSetUp(mat_elemental));
940: }
941: for (row = 0; row < M; row++) {
942: PetscCall(MatGetRow(A, row, &ncols, &cols, &vals));
943: /* PETSc-Elemental interface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
944: PetscCall(MatSetValues(mat_elemental, 1, &row, ncols, cols, vals, ADD_VALUES));
945: PetscCall(MatRestoreRow(A, row, &ncols, &cols, &vals));
946: }
947: PetscCall(MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY));
948: PetscCall(MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY));
950: if (reuse == MAT_INPLACE_MATRIX) {
951: PetscCall(MatHeaderReplace(A, &mat_elemental));
952: } else {
953: *newmat = mat_elemental;
954: }
955: PetscFunctionReturn(PETSC_SUCCESS);
956: }
958: PETSC_INTERN PetscErrorCode MatConvert_MPIAIJ_Elemental(Mat A, MatType, MatReuse reuse, Mat *newmat)
959: {
960: Mat mat_elemental;
961: PetscInt row, ncols, rstart = A->rmap->rstart, rend = A->rmap->rend, j;
962: const PetscInt *cols;
963: const PetscScalar *vals;
965: PetscFunctionBegin;
966: if (reuse == MAT_REUSE_MATRIX) {
967: mat_elemental = *newmat;
968: PetscCall(MatZeroEntries(mat_elemental));
969: } else {
970: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental));
971: PetscCall(MatSetSizes(mat_elemental, PETSC_DECIDE, PETSC_DECIDE, A->rmap->N, A->cmap->N));
972: PetscCall(MatSetType(mat_elemental, MATELEMENTAL));
973: PetscCall(MatSetUp(mat_elemental));
974: }
975: for (row = rstart; row < rend; row++) {
976: PetscCall(MatGetRow(A, row, &ncols, &cols, &vals));
977: for (j = 0; j < ncols; j++) {
978: /* PETSc-Elemental interface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
979: PetscCall(MatSetValues(mat_elemental, 1, &row, 1, &cols[j], &vals[j], ADD_VALUES));
980: }
981: PetscCall(MatRestoreRow(A, row, &ncols, &cols, &vals));
982: }
983: PetscCall(MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY));
984: PetscCall(MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY));
986: if (reuse == MAT_INPLACE_MATRIX) {
987: PetscCall(MatHeaderReplace(A, &mat_elemental));
988: } else {
989: *newmat = mat_elemental;
990: }
991: PetscFunctionReturn(PETSC_SUCCESS);
992: }
994: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_Elemental(Mat A, MatType, MatReuse reuse, Mat *newmat)
995: {
996: Mat mat_elemental;
997: PetscInt M = A->rmap->N, N = A->cmap->N, row, ncols, j;
998: const PetscInt *cols;
999: const PetscScalar *vals;
1001: PetscFunctionBegin;
1002: if (reuse == MAT_REUSE_MATRIX) {
1003: mat_elemental = *newmat;
1004: PetscCall(MatZeroEntries(mat_elemental));
1005: } else {
1006: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental));
1007: PetscCall(MatSetSizes(mat_elemental, PETSC_DECIDE, PETSC_DECIDE, M, N));
1008: PetscCall(MatSetType(mat_elemental, MATELEMENTAL));
1009: PetscCall(MatSetUp(mat_elemental));
1010: }
1011: PetscCall(MatGetRowUpperTriangular(A));
1012: for (row = 0; row < M; row++) {
1013: PetscCall(MatGetRow(A, row, &ncols, &cols, &vals));
1014: /* PETSc-Elemental interface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
1015: PetscCall(MatSetValues(mat_elemental, 1, &row, ncols, cols, vals, ADD_VALUES));
1016: for (j = 0; j < ncols; j++) { /* lower triangular part */
1017: PetscScalar v;
1018: if (cols[j] == row) continue;
1019: v = A->hermitian == PETSC_BOOL3_TRUE ? PetscConj(vals[j]) : vals[j];
1020: PetscCall(MatSetValues(mat_elemental, 1, &cols[j], 1, &row, &v, ADD_VALUES));
1021: }
1022: PetscCall(MatRestoreRow(A, row, &ncols, &cols, &vals));
1023: }
1024: PetscCall(MatRestoreRowUpperTriangular(A));
1025: PetscCall(MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY));
1026: PetscCall(MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY));
1028: if (reuse == MAT_INPLACE_MATRIX) {
1029: PetscCall(MatHeaderReplace(A, &mat_elemental));
1030: } else {
1031: *newmat = mat_elemental;
1032: }
1033: PetscFunctionReturn(PETSC_SUCCESS);
1034: }
1036: PETSC_INTERN PetscErrorCode MatConvert_MPISBAIJ_Elemental(Mat A, MatType, MatReuse reuse, Mat *newmat)
1037: {
1038: Mat mat_elemental;
1039: PetscInt M = A->rmap->N, N = A->cmap->N, row, ncols, j, rstart = A->rmap->rstart, rend = A->rmap->rend;
1040: const PetscInt *cols;
1041: const PetscScalar *vals;
1043: PetscFunctionBegin;
1044: if (reuse == MAT_REUSE_MATRIX) {
1045: mat_elemental = *newmat;
1046: PetscCall(MatZeroEntries(mat_elemental));
1047: } else {
1048: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental));
1049: PetscCall(MatSetSizes(mat_elemental, PETSC_DECIDE, PETSC_DECIDE, M, N));
1050: PetscCall(MatSetType(mat_elemental, MATELEMENTAL));
1051: PetscCall(MatSetUp(mat_elemental));
1052: }
1053: PetscCall(MatGetRowUpperTriangular(A));
1054: for (row = rstart; row < rend; row++) {
1055: PetscCall(MatGetRow(A, row, &ncols, &cols, &vals));
1056: /* PETSc-Elemental interface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
1057: PetscCall(MatSetValues(mat_elemental, 1, &row, ncols, cols, vals, ADD_VALUES));
1058: for (j = 0; j < ncols; j++) { /* lower triangular part */
1059: PetscScalar v;
1060: if (cols[j] == row) continue;
1061: v = A->hermitian == PETSC_BOOL3_TRUE ? PetscConj(vals[j]) : vals[j];
1062: PetscCall(MatSetValues(mat_elemental, 1, &cols[j], 1, &row, &v, ADD_VALUES));
1063: }
1064: PetscCall(MatRestoreRow(A, row, &ncols, &cols, &vals));
1065: }
1066: PetscCall(MatRestoreRowUpperTriangular(A));
1067: PetscCall(MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY));
1068: PetscCall(MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY));
1070: if (reuse == MAT_INPLACE_MATRIX) {
1071: PetscCall(MatHeaderReplace(A, &mat_elemental));
1072: } else {
1073: *newmat = mat_elemental;
1074: }
1075: PetscFunctionReturn(PETSC_SUCCESS);
1076: }
1078: static PetscErrorCode MatDestroy_Elemental(Mat A)
1079: {
1080: Mat_Elemental *a = (Mat_Elemental *)A->data;
1081: Mat_Elemental_Grid *commgrid;
1082: PetscBool flg;
1083: MPI_Comm icomm;
1085: PetscFunctionBegin;
1086: delete a->emat;
1087: delete a->P;
1088: delete a->Q;
1090: El::mpi::Comm cxxcomm(PetscObjectComm((PetscObject)A));
1091: PetscCall(PetscCommDuplicate(cxxcomm.comm, &icomm, nullptr));
1092: PetscCallMPI(MPI_Comm_get_attr(icomm, Petsc_Elemental_keyval, (void **)&commgrid, (int *)&flg));
1093: if (--commgrid->grid_refct == 0) {
1094: delete commgrid->grid;
1095: PetscCall(PetscFree(commgrid));
1096: PetscCallMPI(MPI_Comm_free_keyval(&Petsc_Elemental_keyval));
1097: }
1098: PetscCall(PetscCommDestroy(&icomm));
1099: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatGetOwnershipIS_C", nullptr));
1100: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", nullptr));
1101: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_elemental_mpidense_C", nullptr));
1102: PetscCall(PetscFree(A->data));
1103: PetscFunctionReturn(PETSC_SUCCESS);
1104: }
1106: static PetscErrorCode MatSetUp_Elemental(Mat A)
1107: {
1108: Mat_Elemental *a = (Mat_Elemental *)A->data;
1109: MPI_Comm comm;
1110: PetscMPIInt rsize, csize;
1111: PetscInt n;
1113: PetscFunctionBegin;
1114: PetscCall(PetscLayoutSetUp(A->rmap));
1115: PetscCall(PetscLayoutSetUp(A->cmap));
1117: /* Check if local row and column sizes are equally distributed.
1118: Jed: Elemental uses "element" cyclic ordering so the sizes need to match that
1119: exactly. The strategy in MatElemental is for PETSc to implicitly permute to block ordering (like would be returned by
1120: PetscSplitOwnership(comm,&n,&N), at which point Elemental matrices can act on PETSc vectors without redistributing the vectors. */
1121: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
1122: n = PETSC_DECIDE;
1123: PetscCall(PetscSplitOwnership(comm, &n, &A->rmap->N));
1124: PetscCheck(n == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local row size %" PetscInt_FMT " of ELEMENTAL matrix must be equally distributed", A->rmap->n);
1126: n = PETSC_DECIDE;
1127: PetscCall(PetscSplitOwnership(comm, &n, &A->cmap->N));
1128: PetscCheck(n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local column size %" PetscInt_FMT " of ELEMENTAL matrix must be equally distributed", A->cmap->n);
1130: a->emat->Resize(A->rmap->N, A->cmap->N);
1131: El::Zero(*a->emat);
1133: PetscCallMPI(MPI_Comm_size(A->rmap->comm, &rsize));
1134: PetscCallMPI(MPI_Comm_size(A->cmap->comm, &csize));
1135: PetscCheck(csize == rsize, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Cannot use row and column communicators of different sizes");
1136: a->commsize = rsize;
1137: a->mr[0] = A->rmap->N % rsize;
1138: if (!a->mr[0]) a->mr[0] = rsize;
1139: a->mr[1] = A->cmap->N % csize;
1140: if (!a->mr[1]) a->mr[1] = csize;
1141: a->m[0] = A->rmap->N / rsize + (a->mr[0] != rsize);
1142: a->m[1] = A->cmap->N / csize + (a->mr[1] != csize);
1143: PetscFunctionReturn(PETSC_SUCCESS);
1144: }
1146: static PetscErrorCode MatAssemblyBegin_Elemental(Mat A, MatAssemblyType)
1147: {
1148: Mat_Elemental *a = (Mat_Elemental *)A->data;
1150: PetscFunctionBegin;
1151: /* printf("Calling ProcessQueues\n"); */
1152: a->emat->ProcessQueues();
1153: /* printf("Finished ProcessQueues\n"); */
1154: PetscFunctionReturn(PETSC_SUCCESS);
1155: }
1157: static PetscErrorCode MatAssemblyEnd_Elemental(Mat, MatAssemblyType)
1158: {
1159: PetscFunctionBegin;
1160: /* Currently does nothing */
1161: PetscFunctionReturn(PETSC_SUCCESS);
1162: }
1164: static PetscErrorCode MatLoad_Elemental(Mat newMat, PetscViewer viewer)
1165: {
1166: Mat Adense, Ae;
1167: MPI_Comm comm;
1169: PetscFunctionBegin;
1170: PetscCall(PetscObjectGetComm((PetscObject)newMat, &comm));
1171: PetscCall(MatCreate(comm, &Adense));
1172: PetscCall(MatSetType(Adense, MATDENSE));
1173: PetscCall(MatLoad(Adense, viewer));
1174: PetscCall(MatConvert(Adense, MATELEMENTAL, MAT_INITIAL_MATRIX, &Ae));
1175: PetscCall(MatDestroy(&Adense));
1176: PetscCall(MatHeaderReplace(newMat, &Ae));
1177: PetscFunctionReturn(PETSC_SUCCESS);
1178: }
1180: static struct _MatOps MatOps_Values = {MatSetValues_Elemental,
1181: nullptr,
1182: nullptr,
1183: MatMult_Elemental,
1184: /* 4*/ MatMultAdd_Elemental,
1185: MatMultTranspose_Elemental,
1186: MatMultTransposeAdd_Elemental,
1187: MatSolve_Elemental,
1188: MatSolveAdd_Elemental,
1189: nullptr,
1190: /*10*/ nullptr,
1191: MatLUFactor_Elemental,
1192: MatCholeskyFactor_Elemental,
1193: nullptr,
1194: MatTranspose_Elemental,
1195: /*15*/ MatGetInfo_Elemental,
1196: nullptr,
1197: MatGetDiagonal_Elemental,
1198: MatDiagonalScale_Elemental,
1199: MatNorm_Elemental,
1200: /*20*/ MatAssemblyBegin_Elemental,
1201: MatAssemblyEnd_Elemental,
1202: MatSetOption_Elemental,
1203: MatZeroEntries_Elemental,
1204: /*24*/ nullptr,
1205: MatLUFactorSymbolic_Elemental,
1206: MatLUFactorNumeric_Elemental,
1207: MatCholeskyFactorSymbolic_Elemental,
1208: MatCholeskyFactorNumeric_Elemental,
1209: /*29*/ MatSetUp_Elemental,
1210: nullptr,
1211: nullptr,
1212: nullptr,
1213: nullptr,
1214: /*34*/ MatDuplicate_Elemental,
1215: nullptr,
1216: nullptr,
1217: nullptr,
1218: nullptr,
1219: /*39*/ MatAXPY_Elemental,
1220: nullptr,
1221: nullptr,
1222: nullptr,
1223: MatCopy_Elemental,
1224: /*44*/ nullptr,
1225: MatScale_Elemental,
1226: MatShift_Basic,
1227: nullptr,
1228: nullptr,
1229: /*49*/ nullptr,
1230: nullptr,
1231: nullptr,
1232: nullptr,
1233: nullptr,
1234: /*54*/ nullptr,
1235: nullptr,
1236: nullptr,
1237: nullptr,
1238: nullptr,
1239: /*59*/ nullptr,
1240: MatDestroy_Elemental,
1241: MatView_Elemental,
1242: nullptr,
1243: nullptr,
1244: /*64*/ nullptr,
1245: nullptr,
1246: nullptr,
1247: nullptr,
1248: nullptr,
1249: /*69*/ nullptr,
1250: nullptr,
1251: MatConvert_Elemental_Dense,
1252: nullptr,
1253: nullptr,
1254: /*74*/ nullptr,
1255: nullptr,
1256: nullptr,
1257: nullptr,
1258: nullptr,
1259: /*79*/ nullptr,
1260: nullptr,
1261: nullptr,
1262: nullptr,
1263: MatLoad_Elemental,
1264: /*84*/ nullptr,
1265: nullptr,
1266: nullptr,
1267: nullptr,
1268: nullptr,
1269: /*89*/ nullptr,
1270: nullptr,
1271: MatMatMultNumeric_Elemental,
1272: nullptr,
1273: nullptr,
1274: /*94*/ nullptr,
1275: nullptr,
1276: nullptr,
1277: MatMatTransposeMultNumeric_Elemental,
1278: nullptr,
1279: /*99*/ MatProductSetFromOptions_Elemental,
1280: nullptr,
1281: nullptr,
1282: MatConjugate_Elemental,
1283: nullptr,
1284: /*104*/ nullptr,
1285: nullptr,
1286: nullptr,
1287: nullptr,
1288: nullptr,
1289: /*109*/ MatMatSolve_Elemental,
1290: nullptr,
1291: nullptr,
1292: nullptr,
1293: MatMissingDiagonal_Elemental,
1294: /*114*/ nullptr,
1295: nullptr,
1296: nullptr,
1297: nullptr,
1298: nullptr,
1299: /*119*/ nullptr,
1300: MatHermitianTranspose_Elemental,
1301: nullptr,
1302: nullptr,
1303: nullptr,
1304: /*124*/ nullptr,
1305: nullptr,
1306: nullptr,
1307: nullptr,
1308: nullptr,
1309: /*129*/ nullptr,
1310: nullptr,
1311: nullptr,
1312: nullptr,
1313: nullptr,
1314: /*134*/ nullptr,
1315: nullptr,
1316: nullptr,
1317: nullptr,
1318: nullptr,
1319: nullptr,
1320: /*140*/ nullptr,
1321: nullptr,
1322: nullptr,
1323: nullptr,
1324: nullptr,
1325: /*145*/ nullptr,
1326: nullptr,
1327: nullptr,
1328: nullptr,
1329: nullptr,
1330: /*150*/ nullptr,
1331: nullptr,
1332: nullptr};
1334: /*MC
1335: MATELEMENTAL = "elemental" - A matrix type for dense matrices using the Elemental package
1337: Use ./configure --download-elemental to install PETSc to use Elemental
1339: Options Database Keys:
1340: + -mat_type elemental - sets the matrix type to "elemental" during a call to MatSetFromOptions()
1341: . -pc_factor_mat_solver_type elemental - to use this direct solver with the option -pc_type lu
1342: - -mat_elemental_grid_height - sets Grid Height for 2D cyclic ordering of internal matrix
1344: Level: beginner
1346: Note:
1347: Note unlike most matrix formats, this format does not store all the matrix entries for a contiguous
1348: range of rows on an MPI rank. Use `MatGetOwnershipIS()` to determine what values are stored on
1349: the given rank.
1351: .seealso: `MATDENSE`, `MATSCALAPACK`, `MatGetOwnershipIS()`
1352: M*/
1353: #if defined(__clang__)
1354: #pragma clang diagnostic push
1355: #pragma clang diagnostic ignored "-Wzero-as-null-pointer-constant"
1356: #endif
1357: PETSC_EXTERN PetscErrorCode MatCreate_Elemental(Mat A)
1358: {
1359: Mat_Elemental *a;
1360: PetscBool flg, flg1;
1361: Mat_Elemental_Grid *commgrid;
1362: MPI_Comm icomm;
1363: PetscInt optv1;
1365: PetscFunctionBegin;
1366: A->ops[0] = MatOps_Values;
1367: A->insertmode = NOT_SET_VALUES;
1369: PetscCall(PetscNew(&a));
1370: A->data = (void *)a;
1372: /* Set up the elemental matrix */
1373: El::mpi::Comm cxxcomm(PetscObjectComm((PetscObject)A));
1375: /* Grid needs to be shared between multiple Mats on the same communicator, implement by attribute caching on the MPI_Comm */
1376: if (Petsc_Elemental_keyval == MPI_KEYVAL_INVALID) {
1377: PetscCallMPI(MPI_Comm_create_keyval(MPI_COMM_NULL_COPY_FN, MPI_COMM_NULL_DELETE_FN, &Petsc_Elemental_keyval, (void *)0));
1378: PetscCall(PetscCitationsRegister(ElementalCitation, &ElementalCite));
1379: }
1380: PetscCall(PetscCommDuplicate(cxxcomm.comm, &icomm, NULL));
1381: PetscCallMPI(MPI_Comm_get_attr(icomm, Petsc_Elemental_keyval, (void **)&commgrid, (int *)&flg));
1382: if (!flg) {
1383: PetscCall(PetscNew(&commgrid));
1385: PetscOptionsBegin(PetscObjectComm((PetscObject)A), ((PetscObject)A)->prefix, "Elemental Options", "Mat");
1386: /* displayed default grid sizes (CommSize,1) are set by us arbitrarily until El::Grid() is called */
1387: PetscCall(PetscOptionsInt("-mat_elemental_grid_height", "Grid Height", "None", El::mpi::Size(cxxcomm), &optv1, &flg1));
1388: if (flg1) {
1389: PetscCheck((El::mpi::Size(cxxcomm) % optv1) == 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Grid Height %" PetscInt_FMT " must evenly divide CommSize %" PetscInt_FMT, optv1, (PetscInt)El::mpi::Size(cxxcomm));
1390: commgrid->grid = new El::Grid(cxxcomm, optv1); /* use user-provided grid height */
1391: } else {
1392: commgrid->grid = new El::Grid(cxxcomm); /* use Elemental default grid sizes */
1393: /* printf("new commgrid->grid = %p\n",commgrid->grid); -- memory leak revealed by valgrind? */
1394: }
1395: commgrid->grid_refct = 1;
1396: PetscCallMPI(MPI_Comm_set_attr(icomm, Petsc_Elemental_keyval, (void *)commgrid));
1398: a->pivoting = 1;
1399: PetscCall(PetscOptionsInt("-mat_elemental_pivoting", "Pivoting", "None", a->pivoting, &a->pivoting, NULL));
1401: PetscOptionsEnd();
1402: } else {
1403: commgrid->grid_refct++;
1404: }
1405: PetscCall(PetscCommDestroy(&icomm));
1406: a->grid = commgrid->grid;
1407: a->emat = new El::DistMatrix<PetscElemScalar>(*a->grid);
1408: a->roworiented = PETSC_TRUE;
1410: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatGetOwnershipIS_C", MatGetOwnershipIS_Elemental));
1411: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_elemental_mpidense_C", MatProductSetFromOptions_Elemental_MPIDense));
1412: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATELEMENTAL));
1413: PetscFunctionReturn(PETSC_SUCCESS);
1414: }
1415: #if defined(__clang__)
1416: #pragma clang diagnostic pop
1417: #endif