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(MatSetUp(C));
403: C->ops->matmultnumeric = MatMatMultNumeric_Elemental_MPIDense;
404: PetscFunctionReturn(PETSC_SUCCESS);
405: }
407: static PetscErrorCode MatProductSetFromOptions_Elemental_MPIDense_AB(Mat C)
408: {
409: PetscFunctionBegin;
410: C->ops->matmultsymbolic = MatMatMultSymbolic_Elemental_MPIDense;
411: C->ops->productsymbolic = MatProductSymbolic_AB;
412: PetscFunctionReturn(PETSC_SUCCESS);
413: }
415: static PetscErrorCode MatProductSetFromOptions_Elemental_MPIDense(Mat C)
416: {
417: Mat_Product *product = C->product;
419: PetscFunctionBegin;
420: if (product->type == MATPRODUCT_AB) PetscCall(MatProductSetFromOptions_Elemental_MPIDense_AB(C));
421: PetscFunctionReturn(PETSC_SUCCESS);
422: }
424: static PetscErrorCode MatGetDiagonal_Elemental(Mat A, Vec D)
425: {
426: PetscInt i, nrows, ncols, nD, rrank, ridx, crank, cidx;
427: Mat_Elemental *a = (Mat_Elemental *)A->data;
428: PetscElemScalar v;
429: MPI_Comm comm;
431: PetscFunctionBegin;
432: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
433: PetscCall(MatGetSize(A, &nrows, &ncols));
434: nD = nrows > ncols ? ncols : nrows;
435: for (i = 0; i < nD; i++) {
436: PetscInt erow, ecol;
437: P2RO(A, 0, i, &rrank, &ridx);
438: RO2E(A, 0, rrank, ridx, &erow);
439: PetscCheck(rrank >= 0 && ridx >= 0 && erow >= 0, comm, PETSC_ERR_PLIB, "Incorrect row translation");
440: P2RO(A, 1, i, &crank, &cidx);
441: RO2E(A, 1, crank, cidx, &ecol);
442: PetscCheck(crank >= 0 && cidx >= 0 && ecol >= 0, comm, PETSC_ERR_PLIB, "Incorrect col translation");
443: v = a->emat->Get(erow, ecol);
444: PetscCall(VecSetValues(D, 1, &i, (PetscScalar *)&v, INSERT_VALUES));
445: }
446: PetscCall(VecAssemblyBegin(D));
447: PetscCall(VecAssemblyEnd(D));
448: PetscFunctionReturn(PETSC_SUCCESS);
449: }
451: static PetscErrorCode MatDiagonalScale_Elemental(Mat X, Vec L, Vec R)
452: {
453: Mat_Elemental *x = (Mat_Elemental *)X->data;
454: const PetscElemScalar *d;
456: PetscFunctionBegin;
457: if (R) {
458: PetscCall(VecGetArrayRead(R, (const PetscScalar **)&d));
459: El::DistMatrix<PetscElemScalar, El::VC, El::STAR> de;
460: de.LockedAttach(X->cmap->N, 1, *x->grid, 0, 0, d, X->cmap->n);
461: El::DiagonalScale(El::RIGHT, El::NORMAL, de, *x->emat);
462: PetscCall(VecRestoreArrayRead(R, (const PetscScalar **)&d));
463: }
464: if (L) {
465: PetscCall(VecGetArrayRead(L, (const PetscScalar **)&d));
466: El::DistMatrix<PetscElemScalar, El::VC, El::STAR> de;
467: de.LockedAttach(X->rmap->N, 1, *x->grid, 0, 0, d, X->rmap->n);
468: El::DiagonalScale(El::LEFT, El::NORMAL, de, *x->emat);
469: PetscCall(VecRestoreArrayRead(L, (const PetscScalar **)&d));
470: }
471: PetscFunctionReturn(PETSC_SUCCESS);
472: }
474: static PetscErrorCode MatScale_Elemental(Mat X, PetscScalar a)
475: {
476: Mat_Elemental *x = (Mat_Elemental *)X->data;
478: PetscFunctionBegin;
479: El::Scale((PetscElemScalar)a, *x->emat);
480: PetscFunctionReturn(PETSC_SUCCESS);
481: }
483: /*
484: MatAXPY - Computes Y = a*X + Y.
485: */
486: static PetscErrorCode MatAXPY_Elemental(Mat Y, PetscScalar a, Mat X, MatStructure)
487: {
488: Mat_Elemental *x = (Mat_Elemental *)X->data;
489: Mat_Elemental *y = (Mat_Elemental *)Y->data;
491: PetscFunctionBegin;
492: El::Axpy((PetscElemScalar)a, *x->emat, *y->emat);
493: PetscCall(PetscObjectStateIncrease((PetscObject)Y));
494: PetscFunctionReturn(PETSC_SUCCESS);
495: }
497: static PetscErrorCode MatCopy_Elemental(Mat A, Mat B, MatStructure)
498: {
499: Mat_Elemental *a = (Mat_Elemental *)A->data;
500: Mat_Elemental *b = (Mat_Elemental *)B->data;
502: PetscFunctionBegin;
503: El::Copy(*a->emat, *b->emat);
504: PetscCall(PetscObjectStateIncrease((PetscObject)B));
505: PetscFunctionReturn(PETSC_SUCCESS);
506: }
508: static PetscErrorCode MatDuplicate_Elemental(Mat A, MatDuplicateOption op, Mat *B)
509: {
510: Mat Be;
511: MPI_Comm comm;
512: Mat_Elemental *a = (Mat_Elemental *)A->data;
514: PetscFunctionBegin;
515: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
516: PetscCall(MatCreate(comm, &Be));
517: PetscCall(MatSetSizes(Be, A->rmap->n, A->cmap->n, PETSC_DECIDE, PETSC_DECIDE));
518: PetscCall(MatSetType(Be, MATELEMENTAL));
519: PetscCall(MatSetUp(Be));
520: *B = Be;
521: if (op == MAT_COPY_VALUES) {
522: Mat_Elemental *b = (Mat_Elemental *)Be->data;
523: El::Copy(*a->emat, *b->emat);
524: }
525: Be->assembled = PETSC_TRUE;
526: PetscFunctionReturn(PETSC_SUCCESS);
527: }
529: static PetscErrorCode MatTranspose_Elemental(Mat A, MatReuse reuse, Mat *B)
530: {
531: Mat Be = *B;
532: MPI_Comm comm;
533: Mat_Elemental *a = (Mat_Elemental *)A->data, *b;
535: PetscFunctionBegin;
536: if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *B));
537: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
538: /* Only out-of-place supported */
539: PetscCheck(reuse != MAT_INPLACE_MATRIX, comm, PETSC_ERR_SUP, "Only out-of-place supported");
540: if (reuse == MAT_INITIAL_MATRIX) {
541: PetscCall(MatCreate(comm, &Be));
542: PetscCall(MatSetSizes(Be, A->cmap->n, A->rmap->n, PETSC_DECIDE, PETSC_DECIDE));
543: PetscCall(MatSetType(Be, MATELEMENTAL));
544: PetscCall(MatSetUp(Be));
545: *B = Be;
546: }
547: b = (Mat_Elemental *)Be->data;
548: El::Transpose(*a->emat, *b->emat);
549: Be->assembled = PETSC_TRUE;
550: PetscFunctionReturn(PETSC_SUCCESS);
551: }
553: static PetscErrorCode MatConjugate_Elemental(Mat A)
554: {
555: Mat_Elemental *a = (Mat_Elemental *)A->data;
557: PetscFunctionBegin;
558: El::Conjugate(*a->emat);
559: PetscFunctionReturn(PETSC_SUCCESS);
560: }
562: static PetscErrorCode MatHermitianTranspose_Elemental(Mat A, MatReuse reuse, Mat *B)
563: {
564: Mat Be = *B;
565: MPI_Comm comm;
566: Mat_Elemental *a = (Mat_Elemental *)A->data, *b;
568: PetscFunctionBegin;
569: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
570: /* Only out-of-place supported */
571: if (reuse == MAT_INITIAL_MATRIX) {
572: PetscCall(MatCreate(comm, &Be));
573: PetscCall(MatSetSizes(Be, A->cmap->n, A->rmap->n, PETSC_DECIDE, PETSC_DECIDE));
574: PetscCall(MatSetType(Be, MATELEMENTAL));
575: PetscCall(MatSetUp(Be));
576: *B = Be;
577: }
578: b = (Mat_Elemental *)Be->data;
579: El::Adjoint(*a->emat, *b->emat);
580: Be->assembled = PETSC_TRUE;
581: PetscFunctionReturn(PETSC_SUCCESS);
582: }
584: static PetscErrorCode MatSolve_Elemental(Mat A, Vec B, Vec X)
585: {
586: Mat_Elemental *a = (Mat_Elemental *)A->data;
587: PetscElemScalar *x;
588: PetscInt pivoting = a->pivoting;
590: PetscFunctionBegin;
591: PetscCall(VecCopy(B, X));
592: PetscCall(VecGetArray(X, (PetscScalar **)&x));
594: El::DistMatrix<PetscElemScalar, El::VC, El::STAR> xe;
595: xe.Attach(A->rmap->N, 1, *a->grid, 0, 0, x, A->rmap->n);
596: El::DistMatrix<PetscElemScalar, El::MC, El::MR> xer(xe);
597: switch (A->factortype) {
598: case MAT_FACTOR_LU:
599: if (pivoting == 0) {
600: El::lu::SolveAfter(El::NORMAL, *a->emat, xer);
601: } else if (pivoting == 1) {
602: El::lu::SolveAfter(El::NORMAL, *a->emat, *a->P, xer);
603: } else { /* pivoting == 2 */
604: El::lu::SolveAfter(El::NORMAL, *a->emat, *a->P, *a->Q, xer);
605: }
606: break;
607: case MAT_FACTOR_CHOLESKY:
608: El::cholesky::SolveAfter(El::UPPER, El::NORMAL, *a->emat, xer);
609: break;
610: default:
611: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unfactored Matrix or Unsupported MatFactorType");
612: break;
613: }
614: El::Copy(xer, xe);
616: PetscCall(VecRestoreArray(X, (PetscScalar **)&x));
617: PetscFunctionReturn(PETSC_SUCCESS);
618: }
620: static PetscErrorCode MatSolveAdd_Elemental(Mat A, Vec B, Vec Y, Vec X)
621: {
622: PetscFunctionBegin;
623: PetscCall(MatSolve_Elemental(A, B, X));
624: PetscCall(VecAXPY(X, 1, Y));
625: PetscFunctionReturn(PETSC_SUCCESS);
626: }
628: static PetscErrorCode MatMatSolve_Elemental(Mat A, Mat B, Mat X)
629: {
630: Mat_Elemental *a = (Mat_Elemental *)A->data;
631: Mat_Elemental *x;
632: Mat C;
633: PetscInt pivoting = a->pivoting;
634: PetscBool flg;
635: MatType type;
637: PetscFunctionBegin;
638: PetscCall(MatGetType(X, &type));
639: PetscCall(PetscStrcmp(type, MATELEMENTAL, &flg));
640: if (!flg) {
641: PetscCall(MatConvert(B, MATELEMENTAL, MAT_INITIAL_MATRIX, &C));
642: x = (Mat_Elemental *)C->data;
643: } else {
644: x = (Mat_Elemental *)X->data;
645: El::Copy(*((Mat_Elemental *)B->data)->emat, *x->emat);
646: }
647: switch (A->factortype) {
648: case MAT_FACTOR_LU:
649: if (pivoting == 0) {
650: El::lu::SolveAfter(El::NORMAL, *a->emat, *x->emat);
651: } else if (pivoting == 1) {
652: El::lu::SolveAfter(El::NORMAL, *a->emat, *a->P, *x->emat);
653: } else {
654: El::lu::SolveAfter(El::NORMAL, *a->emat, *a->P, *a->Q, *x->emat);
655: }
656: break;
657: case MAT_FACTOR_CHOLESKY:
658: El::cholesky::SolveAfter(El::UPPER, El::NORMAL, *a->emat, *x->emat);
659: break;
660: default:
661: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unfactored Matrix or Unsupported MatFactorType");
662: break;
663: }
664: if (!flg) {
665: PetscCall(MatConvert(C, type, MAT_REUSE_MATRIX, &X));
666: PetscCall(MatDestroy(&C));
667: }
668: PetscFunctionReturn(PETSC_SUCCESS);
669: }
671: static PetscErrorCode MatLUFactor_Elemental(Mat A, IS, IS, const MatFactorInfo *)
672: {
673: Mat_Elemental *a = (Mat_Elemental *)A->data;
674: PetscInt pivoting = a->pivoting;
676: PetscFunctionBegin;
677: if (pivoting == 0) {
678: El::LU(*a->emat);
679: } else if (pivoting == 1) {
680: El::LU(*a->emat, *a->P);
681: } else {
682: El::LU(*a->emat, *a->P, *a->Q);
683: }
684: A->factortype = MAT_FACTOR_LU;
685: A->assembled = PETSC_TRUE;
687: PetscCall(PetscFree(A->solvertype));
688: PetscCall(PetscStrallocpy(MATSOLVERELEMENTAL, &A->solvertype));
689: PetscFunctionReturn(PETSC_SUCCESS);
690: }
692: static PetscErrorCode MatLUFactorNumeric_Elemental(Mat F, Mat A, const MatFactorInfo *info)
693: {
694: PetscFunctionBegin;
695: PetscCall(MatCopy(A, F, SAME_NONZERO_PATTERN));
696: PetscCall(MatLUFactor_Elemental(F, nullptr, nullptr, info));
697: PetscFunctionReturn(PETSC_SUCCESS);
698: }
700: static PetscErrorCode MatLUFactorSymbolic_Elemental(Mat, Mat, IS, IS, const MatFactorInfo *)
701: {
702: PetscFunctionBegin;
703: /* F is created and allocated by MatGetFactor_elemental_petsc(), skip this routine. */
704: PetscFunctionReturn(PETSC_SUCCESS);
705: }
707: static PetscErrorCode MatCholeskyFactor_Elemental(Mat A, IS, const MatFactorInfo *)
708: {
709: Mat_Elemental *a = (Mat_Elemental *)A->data;
710: El::DistMatrix<PetscElemScalar, El::MC, El::STAR> d;
712: PetscFunctionBegin;
713: El::Cholesky(El::UPPER, *a->emat);
714: A->factortype = MAT_FACTOR_CHOLESKY;
715: A->assembled = PETSC_TRUE;
717: PetscCall(PetscFree(A->solvertype));
718: PetscCall(PetscStrallocpy(MATSOLVERELEMENTAL, &A->solvertype));
719: PetscFunctionReturn(PETSC_SUCCESS);
720: }
722: static PetscErrorCode MatCholeskyFactorNumeric_Elemental(Mat F, Mat A, const MatFactorInfo *info)
723: {
724: PetscFunctionBegin;
725: PetscCall(MatCopy(A, F, SAME_NONZERO_PATTERN));
726: PetscCall(MatCholeskyFactor_Elemental(F, nullptr, info));
727: PetscFunctionReturn(PETSC_SUCCESS);
728: }
730: static PetscErrorCode MatCholeskyFactorSymbolic_Elemental(Mat, Mat, IS, const MatFactorInfo *)
731: {
732: PetscFunctionBegin;
733: /* F is created and allocated by MatGetFactor_elemental_petsc(), skip this routine. */
734: PetscFunctionReturn(PETSC_SUCCESS);
735: }
737: static PetscErrorCode MatFactorGetSolverType_elemental_elemental(Mat, MatSolverType *type)
738: {
739: PetscFunctionBegin;
740: *type = MATSOLVERELEMENTAL;
741: PetscFunctionReturn(PETSC_SUCCESS);
742: }
744: static PetscErrorCode MatGetFactor_elemental_elemental(Mat A, MatFactorType ftype, Mat *F)
745: {
746: Mat B;
748: PetscFunctionBegin;
749: /* Create the factorization matrix */
750: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
751: PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, PETSC_DECIDE, PETSC_DECIDE));
752: PetscCall(MatSetType(B, MATELEMENTAL));
753: PetscCall(MatSetUp(B));
754: B->factortype = ftype;
755: B->trivialsymbolic = PETSC_TRUE;
756: PetscCall(PetscFree(B->solvertype));
757: PetscCall(PetscStrallocpy(MATSOLVERELEMENTAL, &B->solvertype));
759: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_elemental_elemental));
760: *F = B;
761: PetscFunctionReturn(PETSC_SUCCESS);
762: }
764: PETSC_INTERN PetscErrorCode MatSolverTypeRegister_Elemental(void)
765: {
766: PetscFunctionBegin;
767: PetscCall(MatSolverTypeRegister(MATSOLVERELEMENTAL, MATELEMENTAL, MAT_FACTOR_LU, MatGetFactor_elemental_elemental));
768: PetscCall(MatSolverTypeRegister(MATSOLVERELEMENTAL, MATELEMENTAL, MAT_FACTOR_CHOLESKY, MatGetFactor_elemental_elemental));
769: PetscFunctionReturn(PETSC_SUCCESS);
770: }
772: static PetscErrorCode MatNorm_Elemental(Mat A, NormType type, PetscReal *nrm)
773: {
774: Mat_Elemental *a = (Mat_Elemental *)A->data;
776: PetscFunctionBegin;
777: switch (type) {
778: case NORM_1:
779: *nrm = El::OneNorm(*a->emat);
780: break;
781: case NORM_FROBENIUS:
782: *nrm = El::FrobeniusNorm(*a->emat);
783: break;
784: case NORM_INFINITY:
785: *nrm = El::InfinityNorm(*a->emat);
786: break;
787: default:
788: SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Unsupported norm type");
789: }
790: PetscFunctionReturn(PETSC_SUCCESS);
791: }
793: static PetscErrorCode MatZeroEntries_Elemental(Mat A)
794: {
795: Mat_Elemental *a = (Mat_Elemental *)A->data;
797: PetscFunctionBegin;
798: El::Zero(*a->emat);
799: PetscFunctionReturn(PETSC_SUCCESS);
800: }
802: static PetscErrorCode MatGetOwnershipIS_Elemental(Mat A, IS *rows, IS *cols)
803: {
804: Mat_Elemental *a = (Mat_Elemental *)A->data;
805: PetscInt i, m, shift, stride, *idx;
807: PetscFunctionBegin;
808: if (rows) {
809: m = a->emat->LocalHeight();
810: shift = a->emat->ColShift();
811: stride = a->emat->ColStride();
812: PetscCall(PetscMalloc1(m, &idx));
813: for (i = 0; i < m; i++) {
814: PetscInt rank, offset;
815: E2RO(A, 0, shift + i * stride, &rank, &offset);
816: RO2P(A, 0, rank, offset, &idx[i]);
817: }
818: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, m, idx, PETSC_OWN_POINTER, rows));
819: }
820: if (cols) {
821: m = a->emat->LocalWidth();
822: shift = a->emat->RowShift();
823: stride = a->emat->RowStride();
824: PetscCall(PetscMalloc1(m, &idx));
825: for (i = 0; i < m; i++) {
826: PetscInt rank, offset;
827: E2RO(A, 1, shift + i * stride, &rank, &offset);
828: RO2P(A, 1, rank, offset, &idx[i]);
829: }
830: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, m, idx, PETSC_OWN_POINTER, cols));
831: }
832: PetscFunctionReturn(PETSC_SUCCESS);
833: }
835: static PetscErrorCode MatConvert_Elemental_Dense(Mat A, MatType, MatReuse reuse, Mat *B)
836: {
837: Mat Bmpi;
838: Mat_Elemental *a = (Mat_Elemental *)A->data;
839: MPI_Comm comm;
840: IS isrows, iscols;
841: PetscInt rrank, ridx, crank, cidx, nrows, ncols, i, j, erow, ecol, elrow, elcol;
842: const PetscInt *rows, *cols;
843: PetscElemScalar v;
844: const El::Grid &grid = a->emat->Grid();
846: PetscFunctionBegin;
847: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
849: if (reuse == MAT_REUSE_MATRIX) {
850: Bmpi = *B;
851: } else {
852: PetscCall(MatCreate(comm, &Bmpi));
853: PetscCall(MatSetSizes(Bmpi, A->rmap->n, A->cmap->n, PETSC_DECIDE, PETSC_DECIDE));
854: PetscCall(MatSetType(Bmpi, MATDENSE));
855: PetscCall(MatSetUp(Bmpi));
856: }
858: /* Get local entries of A */
859: PetscCall(MatGetOwnershipIS(A, &isrows, &iscols));
860: PetscCall(ISGetLocalSize(isrows, &nrows));
861: PetscCall(ISGetIndices(isrows, &rows));
862: PetscCall(ISGetLocalSize(iscols, &ncols));
863: PetscCall(ISGetIndices(iscols, &cols));
865: if (a->roworiented) {
866: for (i = 0; i < nrows; i++) {
867: P2RO(A, 0, rows[i], &rrank, &ridx); /* convert indices between PETSc <-> (Rank,Offset) <-> Elemental */
868: RO2E(A, 0, rrank, ridx, &erow);
869: PetscCheck(rrank >= 0 && ridx >= 0 && erow >= 0, comm, PETSC_ERR_PLIB, "Incorrect row translation");
870: for (j = 0; j < ncols; j++) {
871: P2RO(A, 1, cols[j], &crank, &cidx);
872: RO2E(A, 1, crank, cidx, &ecol);
873: PetscCheck(crank >= 0 && cidx >= 0 && ecol >= 0, comm, PETSC_ERR_PLIB, "Incorrect col translation");
875: elrow = erow / grid.MCSize(); /* Elemental local row index */
876: elcol = ecol / grid.MRSize(); /* Elemental local column index */
877: v = a->emat->GetLocal(elrow, elcol);
878: PetscCall(MatSetValues(Bmpi, 1, &rows[i], 1, &cols[j], (PetscScalar *)&v, INSERT_VALUES));
879: }
880: }
881: } else { /* column-oriented */
882: for (j = 0; j < ncols; j++) {
883: P2RO(A, 1, cols[j], &crank, &cidx);
884: RO2E(A, 1, crank, cidx, &ecol);
885: PetscCheck(crank >= 0 && cidx >= 0 && ecol >= 0, comm, PETSC_ERR_PLIB, "Incorrect col translation");
886: for (i = 0; i < nrows; i++) {
887: P2RO(A, 0, rows[i], &rrank, &ridx); /* convert indices between PETSc <-> (Rank,Offset) <-> Elemental */
888: RO2E(A, 0, rrank, ridx, &erow);
889: PetscCheck(rrank >= 0 && ridx >= 0 && erow >= 0, comm, PETSC_ERR_PLIB, "Incorrect row translation");
891: elrow = erow / grid.MCSize(); /* Elemental local row index */
892: elcol = ecol / grid.MRSize(); /* Elemental local column index */
893: v = a->emat->GetLocal(elrow, elcol);
894: PetscCall(MatSetValues(Bmpi, 1, &rows[i], 1, &cols[j], (PetscScalar *)&v, INSERT_VALUES));
895: }
896: }
897: }
898: PetscCall(MatAssemblyBegin(Bmpi, MAT_FINAL_ASSEMBLY));
899: PetscCall(MatAssemblyEnd(Bmpi, MAT_FINAL_ASSEMBLY));
900: if (reuse == MAT_INPLACE_MATRIX) {
901: PetscCall(MatHeaderReplace(A, &Bmpi));
902: } else {
903: *B = Bmpi;
904: }
905: PetscCall(ISDestroy(&isrows));
906: PetscCall(ISDestroy(&iscols));
907: PetscFunctionReturn(PETSC_SUCCESS);
908: }
910: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat A, MatType, MatReuse reuse, Mat *newmat)
911: {
912: Mat mat_elemental;
913: PetscInt M = A->rmap->N, N = A->cmap->N, row, ncols;
914: const PetscInt *cols;
915: const PetscScalar *vals;
917: PetscFunctionBegin;
918: if (reuse == MAT_REUSE_MATRIX) {
919: mat_elemental = *newmat;
920: PetscCall(MatZeroEntries(mat_elemental));
921: } else {
922: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental));
923: PetscCall(MatSetSizes(mat_elemental, PETSC_DECIDE, PETSC_DECIDE, M, N));
924: PetscCall(MatSetType(mat_elemental, MATELEMENTAL));
925: PetscCall(MatSetUp(mat_elemental));
926: }
927: for (row = 0; row < M; row++) {
928: PetscCall(MatGetRow(A, row, &ncols, &cols, &vals));
929: /* PETSc-Elemental interface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
930: PetscCall(MatSetValues(mat_elemental, 1, &row, ncols, cols, vals, ADD_VALUES));
931: PetscCall(MatRestoreRow(A, row, &ncols, &cols, &vals));
932: }
933: PetscCall(MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY));
934: PetscCall(MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY));
936: if (reuse == MAT_INPLACE_MATRIX) {
937: PetscCall(MatHeaderReplace(A, &mat_elemental));
938: } else {
939: *newmat = mat_elemental;
940: }
941: PetscFunctionReturn(PETSC_SUCCESS);
942: }
944: PETSC_INTERN PetscErrorCode MatConvert_MPIAIJ_Elemental(Mat A, MatType, MatReuse reuse, Mat *newmat)
945: {
946: Mat mat_elemental;
947: PetscInt row, ncols, rstart = A->rmap->rstart, rend = A->rmap->rend, j;
948: const PetscInt *cols;
949: const PetscScalar *vals;
951: PetscFunctionBegin;
952: if (reuse == MAT_REUSE_MATRIX) {
953: mat_elemental = *newmat;
954: PetscCall(MatZeroEntries(mat_elemental));
955: } else {
956: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental));
957: PetscCall(MatSetSizes(mat_elemental, PETSC_DECIDE, PETSC_DECIDE, A->rmap->N, A->cmap->N));
958: PetscCall(MatSetType(mat_elemental, MATELEMENTAL));
959: PetscCall(MatSetUp(mat_elemental));
960: }
961: for (row = rstart; row < rend; row++) {
962: PetscCall(MatGetRow(A, row, &ncols, &cols, &vals));
963: for (j = 0; j < ncols; j++) {
964: /* PETSc-Elemental interface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
965: PetscCall(MatSetValues(mat_elemental, 1, &row, 1, &cols[j], &vals[j], ADD_VALUES));
966: }
967: PetscCall(MatRestoreRow(A, row, &ncols, &cols, &vals));
968: }
969: PetscCall(MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY));
970: PetscCall(MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY));
972: if (reuse == MAT_INPLACE_MATRIX) {
973: PetscCall(MatHeaderReplace(A, &mat_elemental));
974: } else {
975: *newmat = mat_elemental;
976: }
977: PetscFunctionReturn(PETSC_SUCCESS);
978: }
980: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_Elemental(Mat A, MatType, MatReuse reuse, Mat *newmat)
981: {
982: Mat mat_elemental;
983: PetscInt M = A->rmap->N, N = A->cmap->N, row, ncols, j;
984: const PetscInt *cols;
985: const PetscScalar *vals;
987: PetscFunctionBegin;
988: if (reuse == MAT_REUSE_MATRIX) {
989: mat_elemental = *newmat;
990: PetscCall(MatZeroEntries(mat_elemental));
991: } else {
992: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental));
993: PetscCall(MatSetSizes(mat_elemental, PETSC_DECIDE, PETSC_DECIDE, M, N));
994: PetscCall(MatSetType(mat_elemental, MATELEMENTAL));
995: PetscCall(MatSetUp(mat_elemental));
996: }
997: PetscCall(MatGetRowUpperTriangular(A));
998: for (row = 0; row < M; row++) {
999: PetscCall(MatGetRow(A, row, &ncols, &cols, &vals));
1000: /* PETSc-Elemental interface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
1001: PetscCall(MatSetValues(mat_elemental, 1, &row, ncols, cols, vals, ADD_VALUES));
1002: for (j = 0; j < ncols; j++) { /* lower triangular part */
1003: PetscScalar v;
1004: if (cols[j] == row) continue;
1005: v = A->hermitian == PETSC_BOOL3_TRUE ? PetscConj(vals[j]) : vals[j];
1006: PetscCall(MatSetValues(mat_elemental, 1, &cols[j], 1, &row, &v, ADD_VALUES));
1007: }
1008: PetscCall(MatRestoreRow(A, row, &ncols, &cols, &vals));
1009: }
1010: PetscCall(MatRestoreRowUpperTriangular(A));
1011: PetscCall(MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY));
1012: PetscCall(MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY));
1014: if (reuse == MAT_INPLACE_MATRIX) {
1015: PetscCall(MatHeaderReplace(A, &mat_elemental));
1016: } else {
1017: *newmat = mat_elemental;
1018: }
1019: PetscFunctionReturn(PETSC_SUCCESS);
1020: }
1022: PETSC_INTERN PetscErrorCode MatConvert_MPISBAIJ_Elemental(Mat A, MatType, MatReuse reuse, Mat *newmat)
1023: {
1024: Mat mat_elemental;
1025: PetscInt M = A->rmap->N, N = A->cmap->N, row, ncols, j, rstart = A->rmap->rstart, rend = A->rmap->rend;
1026: const PetscInt *cols;
1027: const PetscScalar *vals;
1029: PetscFunctionBegin;
1030: if (reuse == MAT_REUSE_MATRIX) {
1031: mat_elemental = *newmat;
1032: PetscCall(MatZeroEntries(mat_elemental));
1033: } else {
1034: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental));
1035: PetscCall(MatSetSizes(mat_elemental, PETSC_DECIDE, PETSC_DECIDE, M, N));
1036: PetscCall(MatSetType(mat_elemental, MATELEMENTAL));
1037: PetscCall(MatSetUp(mat_elemental));
1038: }
1039: PetscCall(MatGetRowUpperTriangular(A));
1040: for (row = rstart; row < rend; row++) {
1041: PetscCall(MatGetRow(A, row, &ncols, &cols, &vals));
1042: /* PETSc-Elemental interface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
1043: PetscCall(MatSetValues(mat_elemental, 1, &row, ncols, cols, vals, ADD_VALUES));
1044: for (j = 0; j < ncols; j++) { /* lower triangular part */
1045: PetscScalar v;
1046: if (cols[j] == row) continue;
1047: v = A->hermitian == PETSC_BOOL3_TRUE ? PetscConj(vals[j]) : vals[j];
1048: PetscCall(MatSetValues(mat_elemental, 1, &cols[j], 1, &row, &v, ADD_VALUES));
1049: }
1050: PetscCall(MatRestoreRow(A, row, &ncols, &cols, &vals));
1051: }
1052: PetscCall(MatRestoreRowUpperTriangular(A));
1053: PetscCall(MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY));
1054: PetscCall(MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY));
1056: if (reuse == MAT_INPLACE_MATRIX) {
1057: PetscCall(MatHeaderReplace(A, &mat_elemental));
1058: } else {
1059: *newmat = mat_elemental;
1060: }
1061: PetscFunctionReturn(PETSC_SUCCESS);
1062: }
1064: static PetscErrorCode MatDestroy_Elemental(Mat A)
1065: {
1066: Mat_Elemental *a = (Mat_Elemental *)A->data;
1067: Mat_Elemental_Grid *commgrid;
1068: PetscMPIInt iflg;
1069: MPI_Comm icomm;
1071: PetscFunctionBegin;
1072: delete a->emat;
1073: delete a->P;
1074: delete a->Q;
1076: El::mpi::Comm cxxcomm(PetscObjectComm((PetscObject)A));
1077: PetscCall(PetscCommDuplicate(cxxcomm.comm, &icomm, nullptr));
1078: PetscCallMPI(MPI_Comm_get_attr(icomm, Petsc_Elemental_keyval, (void **)&commgrid, &iflg));
1079: if (--commgrid->grid_refct == 0) {
1080: delete commgrid->grid;
1081: PetscCall(PetscFree(commgrid));
1082: PetscCallMPI(MPI_Comm_free_keyval(&Petsc_Elemental_keyval));
1083: }
1084: PetscCall(PetscCommDestroy(&icomm));
1085: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatGetOwnershipIS_C", nullptr));
1086: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", nullptr));
1087: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_elemental_mpidense_C", nullptr));
1088: PetscCall(PetscFree(A->data));
1089: PetscFunctionReturn(PETSC_SUCCESS);
1090: }
1092: static PetscErrorCode MatSetUp_Elemental(Mat A)
1093: {
1094: Mat_Elemental *a = (Mat_Elemental *)A->data;
1095: MPI_Comm comm;
1096: PetscMPIInt rsize, csize;
1097: PetscInt n;
1099: PetscFunctionBegin;
1100: PetscCall(PetscLayoutSetUp(A->rmap));
1101: PetscCall(PetscLayoutSetUp(A->cmap));
1103: /* Check if local row and column sizes are equally distributed.
1104: Jed: Elemental uses "element" cyclic ordering so the sizes need to match that
1105: exactly. The strategy in MatElemental is for PETSc to implicitly permute to block ordering (like would be returned by
1106: PetscSplitOwnership(comm,&n,&N), at which point Elemental matrices can act on PETSc vectors without redistributing the vectors. */
1107: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
1108: n = PETSC_DECIDE;
1109: PetscCall(PetscSplitOwnership(comm, &n, &A->rmap->N));
1110: 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);
1112: n = PETSC_DECIDE;
1113: PetscCall(PetscSplitOwnership(comm, &n, &A->cmap->N));
1114: 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);
1116: a->emat->Resize(A->rmap->N, A->cmap->N);
1117: El::Zero(*a->emat);
1119: PetscCallMPI(MPI_Comm_size(A->rmap->comm, &rsize));
1120: PetscCallMPI(MPI_Comm_size(A->cmap->comm, &csize));
1121: PetscCheck(csize == rsize, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Cannot use row and column communicators of different sizes");
1122: a->commsize = rsize;
1123: a->mr[0] = A->rmap->N % rsize;
1124: if (!a->mr[0]) a->mr[0] = rsize;
1125: a->mr[1] = A->cmap->N % csize;
1126: if (!a->mr[1]) a->mr[1] = csize;
1127: a->m[0] = A->rmap->N / rsize + (a->mr[0] != rsize);
1128: a->m[1] = A->cmap->N / csize + (a->mr[1] != csize);
1129: PetscFunctionReturn(PETSC_SUCCESS);
1130: }
1132: static PetscErrorCode MatAssemblyBegin_Elemental(Mat A, MatAssemblyType)
1133: {
1134: Mat_Elemental *a = (Mat_Elemental *)A->data;
1136: PetscFunctionBegin;
1137: /* printf("Calling ProcessQueues\n"); */
1138: a->emat->ProcessQueues();
1139: /* printf("Finished ProcessQueues\n"); */
1140: PetscFunctionReturn(PETSC_SUCCESS);
1141: }
1143: static PetscErrorCode MatAssemblyEnd_Elemental(Mat, MatAssemblyType)
1144: {
1145: PetscFunctionBegin;
1146: /* Currently does nothing */
1147: PetscFunctionReturn(PETSC_SUCCESS);
1148: }
1150: static PetscErrorCode MatLoad_Elemental(Mat newMat, PetscViewer viewer)
1151: {
1152: Mat Adense, Ae;
1153: MPI_Comm comm;
1155: PetscFunctionBegin;
1156: PetscCall(PetscObjectGetComm((PetscObject)newMat, &comm));
1157: PetscCall(MatCreate(comm, &Adense));
1158: PetscCall(MatSetType(Adense, MATDENSE));
1159: PetscCall(MatLoad(Adense, viewer));
1160: PetscCall(MatConvert(Adense, MATELEMENTAL, MAT_INITIAL_MATRIX, &Ae));
1161: PetscCall(MatDestroy(&Adense));
1162: PetscCall(MatHeaderReplace(newMat, &Ae));
1163: PetscFunctionReturn(PETSC_SUCCESS);
1164: }
1166: static struct _MatOps MatOps_Values = {MatSetValues_Elemental,
1167: nullptr,
1168: nullptr,
1169: MatMult_Elemental,
1170: /* 4*/ MatMultAdd_Elemental,
1171: MatMultTranspose_Elemental,
1172: MatMultTransposeAdd_Elemental,
1173: MatSolve_Elemental,
1174: MatSolveAdd_Elemental,
1175: nullptr,
1176: /*10*/ nullptr,
1177: MatLUFactor_Elemental,
1178: MatCholeskyFactor_Elemental,
1179: nullptr,
1180: MatTranspose_Elemental,
1181: /*15*/ MatGetInfo_Elemental,
1182: nullptr,
1183: MatGetDiagonal_Elemental,
1184: MatDiagonalScale_Elemental,
1185: MatNorm_Elemental,
1186: /*20*/ MatAssemblyBegin_Elemental,
1187: MatAssemblyEnd_Elemental,
1188: MatSetOption_Elemental,
1189: MatZeroEntries_Elemental,
1190: /*24*/ nullptr,
1191: MatLUFactorSymbolic_Elemental,
1192: MatLUFactorNumeric_Elemental,
1193: MatCholeskyFactorSymbolic_Elemental,
1194: MatCholeskyFactorNumeric_Elemental,
1195: /*29*/ MatSetUp_Elemental,
1196: nullptr,
1197: nullptr,
1198: nullptr,
1199: nullptr,
1200: /*34*/ MatDuplicate_Elemental,
1201: nullptr,
1202: nullptr,
1203: nullptr,
1204: nullptr,
1205: /*39*/ MatAXPY_Elemental,
1206: nullptr,
1207: nullptr,
1208: nullptr,
1209: MatCopy_Elemental,
1210: /*44*/ nullptr,
1211: MatScale_Elemental,
1212: MatShift_Basic,
1213: nullptr,
1214: nullptr,
1215: /*49*/ nullptr,
1216: nullptr,
1217: nullptr,
1218: nullptr,
1219: nullptr,
1220: /*54*/ nullptr,
1221: nullptr,
1222: nullptr,
1223: nullptr,
1224: nullptr,
1225: /*59*/ nullptr,
1226: MatDestroy_Elemental,
1227: MatView_Elemental,
1228: nullptr,
1229: nullptr,
1230: /*64*/ nullptr,
1231: nullptr,
1232: nullptr,
1233: nullptr,
1234: nullptr,
1235: /*69*/ nullptr,
1236: MatConvert_Elemental_Dense,
1237: nullptr,
1238: nullptr,
1239: nullptr,
1240: /*74*/ nullptr,
1241: nullptr,
1242: nullptr,
1243: nullptr,
1244: MatLoad_Elemental,
1245: /*79*/ nullptr,
1246: nullptr,
1247: nullptr,
1248: nullptr,
1249: nullptr,
1250: /*84*/ nullptr,
1251: MatMatMultNumeric_Elemental,
1252: nullptr,
1253: nullptr,
1254: MatMatTransposeMultNumeric_Elemental,
1255: /*89*/ nullptr,
1256: MatProductSetFromOptions_Elemental,
1257: nullptr,
1258: nullptr,
1259: MatConjugate_Elemental,
1260: /*94*/ nullptr,
1261: nullptr,
1262: nullptr,
1263: nullptr,
1264: nullptr,
1265: /*99*/ nullptr,
1266: MatMatSolve_Elemental,
1267: nullptr,
1268: nullptr,
1269: nullptr,
1270: /*104*/ nullptr,
1271: nullptr,
1272: nullptr,
1273: nullptr,
1274: nullptr,
1275: /*109*/ nullptr,
1276: MatHermitianTranspose_Elemental,
1277: nullptr,
1278: nullptr,
1279: /*114*/ nullptr,
1280: nullptr,
1281: nullptr,
1282: nullptr,
1283: nullptr,
1284: /*119*/ nullptr,
1285: nullptr,
1286: nullptr,
1287: nullptr,
1288: nullptr,
1289: /*124*/ nullptr,
1290: nullptr,
1291: nullptr,
1292: nullptr,
1293: nullptr,
1294: /*129*/ nullptr,
1295: nullptr,
1296: nullptr,
1297: nullptr,
1298: nullptr,
1299: /*134*/ nullptr,
1300: nullptr,
1301: nullptr,
1302: nullptr,
1303: nullptr,
1304: nullptr,
1305: /*140*/ nullptr,
1306: nullptr,
1307: nullptr,
1308: nullptr,
1309: nullptr};
1311: /*MC
1312: MATELEMENTAL = "elemental" - A matrix type for dense matrices using the Elemental package
1314: Use ./configure --download-elemental to install PETSc to use Elemental
1316: Options Database Keys:
1317: + -mat_type elemental - sets the matrix type to "elemental" during a call to MatSetFromOptions()
1318: . -pc_factor_mat_solver_type elemental - to use this direct solver with the option -pc_type lu
1319: - -mat_elemental_grid_height - sets Grid Height for 2D cyclic ordering of internal matrix
1321: Level: beginner
1323: Note:
1324: Note unlike most matrix formats, this format does not store all the matrix entries for a contiguous
1325: range of rows on an MPI rank. Use `MatGetOwnershipIS()` to determine what values are stored on
1326: the given rank.
1328: .seealso: `MATDENSE`, `MATSCALAPACK`, `MatGetOwnershipIS()`
1329: M*/
1330: #if defined(__clang__)
1331: #pragma clang diagnostic push
1332: #pragma clang diagnostic ignored "-Wzero-as-null-pointer-constant"
1333: #endif
1334: PETSC_EXTERN PetscErrorCode MatCreate_Elemental(Mat A)
1335: {
1336: Mat_Elemental *a;
1337: PetscBool flg;
1338: PetscMPIInt iflg;
1339: Mat_Elemental_Grid *commgrid;
1340: MPI_Comm icomm;
1341: PetscInt optv1;
1343: PetscFunctionBegin;
1344: A->ops[0] = MatOps_Values;
1345: A->insertmode = NOT_SET_VALUES;
1347: PetscCall(PetscNew(&a));
1348: A->data = (void *)a;
1350: /* Set up the elemental matrix */
1351: El::mpi::Comm cxxcomm(PetscObjectComm((PetscObject)A));
1353: /* Grid needs to be shared between multiple Mats on the same communicator, implement by attribute caching on the MPI_Comm */
1354: if (Petsc_Elemental_keyval == MPI_KEYVAL_INVALID) {
1355: PetscCallMPI(MPI_Comm_create_keyval(MPI_COMM_NULL_COPY_FN, MPI_COMM_NULL_DELETE_FN, &Petsc_Elemental_keyval, nullptr));
1356: PetscCall(PetscCitationsRegister(ElementalCitation, &ElementalCite));
1357: }
1358: PetscCall(PetscCommDuplicate(cxxcomm.comm, &icomm, NULL));
1359: PetscCallMPI(MPI_Comm_get_attr(icomm, Petsc_Elemental_keyval, (void **)&commgrid, &iflg));
1360: if (!iflg) {
1361: PetscCall(PetscNew(&commgrid));
1363: PetscOptionsBegin(PetscObjectComm((PetscObject)A), ((PetscObject)A)->prefix, "Elemental Options", "Mat");
1364: /* displayed default grid sizes (CommSize,1) are set by us arbitrarily until El::Grid() is called */
1365: PetscCall(PetscOptionsInt("-mat_elemental_grid_height", "Grid Height", "None", El::mpi::Size(cxxcomm), &optv1, &flg));
1366: if (flg) {
1367: 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));
1368: commgrid->grid = new El::Grid(cxxcomm, optv1); /* use user-provided grid height */
1369: } else {
1370: commgrid->grid = new El::Grid(cxxcomm); /* use Elemental default grid sizes */
1371: /* printf("new commgrid->grid = %p\n",commgrid->grid); -- memory leak revealed by valgrind? */
1372: }
1373: commgrid->grid_refct = 1;
1374: PetscCallMPI(MPI_Comm_set_attr(icomm, Petsc_Elemental_keyval, (void *)commgrid));
1376: a->pivoting = 1;
1377: PetscCall(PetscOptionsInt("-mat_elemental_pivoting", "Pivoting", "None", a->pivoting, &a->pivoting, NULL));
1379: PetscOptionsEnd();
1380: } else {
1381: commgrid->grid_refct++;
1382: }
1383: PetscCall(PetscCommDestroy(&icomm));
1384: a->grid = commgrid->grid;
1385: a->emat = new El::DistMatrix<PetscElemScalar>(*a->grid);
1386: a->roworiented = PETSC_TRUE;
1388: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatGetOwnershipIS_C", MatGetOwnershipIS_Elemental));
1389: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_elemental_mpidense_C", MatProductSetFromOptions_Elemental_MPIDense));
1390: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATELEMENTAL));
1391: PetscFunctionReturn(PETSC_SUCCESS);
1392: }
1393: #if defined(__clang__)
1394: #pragma clang diagnostic pop
1395: #endif