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