1 #include <petsc/private/petscelemental.h>
2
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;
12
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;
18
MatView_Elemental(Mat A,PetscViewer viewer)19 static PetscErrorCode MatView_Elemental(Mat A, PetscViewer viewer)
20 {
21 Mat_Elemental *a = (Mat_Elemental *)A->data;
22 PetscBool isascii;
23
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 }
38
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 }
59
MatGetInfo_Elemental(Mat A,MatInfoType flag,MatInfo * info)60 static PetscErrorCode MatGetInfo_Elemental(Mat A, MatInfoType flag, MatInfo *info)
61 {
62 Mat_Elemental *a = (Mat_Elemental *)A->data;
63
64 PetscFunctionBegin;
65 info->block_size = 1.0;
66
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 }
81
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 }
91
MatSetOption_Elemental(Mat A,MatOption op,PetscBool flg)92 static PetscErrorCode MatSetOption_Elemental(Mat A, MatOption op, PetscBool flg)
93 {
94 Mat_Elemental *a = (Mat_Elemental *)A->data;
95
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 }
106
MatSetValues_Elemental(Mat A,PetscInt nr,const PetscInt * rows,PetscInt nc,const PetscInt * cols,const PetscScalar * vals,InsertMode imode)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;
111
112 PetscFunctionBegin;
113 // TODO: Initialize matrix to all zeros?
114
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 }
146
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 }
193
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);
200
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 }
214
MatMult_Elemental(Mat A,Vec X,Vec Y)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;
221
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 }
235
MatMultTranspose_Elemental(Mat A,Vec X,Vec Y)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;
242
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 }
256
MatMultAdd_Elemental(Mat A,Vec X,Vec Y,Vec Z)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;
263
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 }
278
MatMultTransposeAdd_Elemental(Mat A,Vec X,Vec Y,Vec Z)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;
285
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 }
300
MatMatMultNumeric_Elemental(Mat A,Mat B,Mat C)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;
307
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 }
315
MatMatMultSymbolic_Elemental(Mat A,Mat B,PetscReal,Mat Ce)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 }
325
MatMatTransposeMultNumeric_Elemental(Mat A,Mat B,Mat C)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;
332
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 }
340
MatMatTransposeMultSymbolic_Elemental(Mat A,Mat B,PetscReal,Mat C)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 }
349
MatProductSetFromOptions_Elemental_AB(Mat C)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 }
357
MatProductSetFromOptions_Elemental_ABt(Mat C)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 }
365
MatProductSetFromOptions_Elemental(Mat C)366 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Elemental(Mat C)
367 {
368 Mat_Product *product = C->product;
369
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 }
383
MatMatMultNumeric_Elemental_MPIDense(Mat A,Mat B,Mat C)384 static PetscErrorCode MatMatMultNumeric_Elemental_MPIDense(Mat A, Mat B, Mat C)
385 {
386 Mat Be, Ce;
387
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 }
396
MatMatMultSymbolic_Elemental_MPIDense(Mat A,Mat B,PetscReal,Mat C)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 }
406
MatProductSetFromOptions_Elemental_MPIDense_AB(Mat C)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 }
414
MatProductSetFromOptions_Elemental_MPIDense(Mat C)415 static PetscErrorCode MatProductSetFromOptions_Elemental_MPIDense(Mat C)
416 {
417 Mat_Product *product = C->product;
418
419 PetscFunctionBegin;
420 if (product->type == MATPRODUCT_AB) PetscCall(MatProductSetFromOptions_Elemental_MPIDense_AB(C));
421 PetscFunctionReturn(PETSC_SUCCESS);
422 }
423
MatGetDiagonal_Elemental(Mat A,Vec D)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;
430
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 }
450
MatDiagonalScale_Elemental(Mat X,Vec L,Vec R)451 static PetscErrorCode MatDiagonalScale_Elemental(Mat X, Vec L, Vec R)
452 {
453 Mat_Elemental *x = (Mat_Elemental *)X->data;
454 const PetscElemScalar *d;
455
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 }
473
MatScale_Elemental(Mat X,PetscScalar a)474 static PetscErrorCode MatScale_Elemental(Mat X, PetscScalar a)
475 {
476 Mat_Elemental *x = (Mat_Elemental *)X->data;
477
478 PetscFunctionBegin;
479 El::Scale((PetscElemScalar)a, *x->emat);
480 PetscFunctionReturn(PETSC_SUCCESS);
481 }
482
483 /*
484 MatAXPY - Computes Y = a*X + Y.
485 */
MatAXPY_Elemental(Mat Y,PetscScalar a,Mat X,MatStructure)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;
490
491 PetscFunctionBegin;
492 El::Axpy((PetscElemScalar)a, *x->emat, *y->emat);
493 PetscCall(PetscObjectStateIncrease((PetscObject)Y));
494 PetscFunctionReturn(PETSC_SUCCESS);
495 }
496
MatCopy_Elemental(Mat A,Mat B,MatStructure)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;
501
502 PetscFunctionBegin;
503 El::Copy(*a->emat, *b->emat);
504 PetscCall(PetscObjectStateIncrease((PetscObject)B));
505 PetscFunctionReturn(PETSC_SUCCESS);
506 }
507
MatDuplicate_Elemental(Mat A,MatDuplicateOption op,Mat * B)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;
513
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 }
528
MatTranspose_Elemental(Mat A,MatReuse reuse,Mat * B)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;
534
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 }
552
MatConjugate_Elemental(Mat A)553 static PetscErrorCode MatConjugate_Elemental(Mat A)
554 {
555 Mat_Elemental *a = (Mat_Elemental *)A->data;
556
557 PetscFunctionBegin;
558 El::Conjugate(*a->emat);
559 PetscFunctionReturn(PETSC_SUCCESS);
560 }
561
MatHermitianTranspose_Elemental(Mat A,MatReuse reuse,Mat * B)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;
567
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 }
583
MatSolve_Elemental(Mat A,Vec B,Vec X)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;
589
590 PetscFunctionBegin;
591 PetscCall(VecCopy(B, X));
592 PetscCall(VecGetArray(X, (PetscScalar **)&x));
593
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);
615
616 PetscCall(VecRestoreArray(X, (PetscScalar **)&x));
617 PetscFunctionReturn(PETSC_SUCCESS);
618 }
619
MatSolveAdd_Elemental(Mat A,Vec B,Vec Y,Vec X)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 }
627
MatMatSolve_Elemental(Mat A,Mat B,Mat X)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;
636
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 }
670
MatLUFactor_Elemental(Mat A,IS,IS,const MatFactorInfo *)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;
675
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;
686
687 PetscCall(PetscFree(A->solvertype));
688 PetscCall(PetscStrallocpy(MATSOLVERELEMENTAL, &A->solvertype));
689 PetscFunctionReturn(PETSC_SUCCESS);
690 }
691
MatLUFactorNumeric_Elemental(Mat F,Mat A,const MatFactorInfo * info)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 }
699
MatLUFactorSymbolic_Elemental(Mat,Mat,IS,IS,const MatFactorInfo *)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 }
706
MatCholeskyFactor_Elemental(Mat A,IS,const MatFactorInfo *)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;
711
712 PetscFunctionBegin;
713 El::Cholesky(El::UPPER, *a->emat);
714 A->factortype = MAT_FACTOR_CHOLESKY;
715 A->assembled = PETSC_TRUE;
716
717 PetscCall(PetscFree(A->solvertype));
718 PetscCall(PetscStrallocpy(MATSOLVERELEMENTAL, &A->solvertype));
719 PetscFunctionReturn(PETSC_SUCCESS);
720 }
721
MatCholeskyFactorNumeric_Elemental(Mat F,Mat A,const MatFactorInfo * info)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 }
729
MatCholeskyFactorSymbolic_Elemental(Mat,Mat,IS,const MatFactorInfo *)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 }
736
MatFactorGetSolverType_elemental_elemental(Mat,MatSolverType * type)737 static PetscErrorCode MatFactorGetSolverType_elemental_elemental(Mat, MatSolverType *type)
738 {
739 PetscFunctionBegin;
740 *type = MATSOLVERELEMENTAL;
741 PetscFunctionReturn(PETSC_SUCCESS);
742 }
743
MatGetFactor_elemental_elemental(Mat A,MatFactorType ftype,Mat * F)744 static PetscErrorCode MatGetFactor_elemental_elemental(Mat A, MatFactorType ftype, Mat *F)
745 {
746 Mat B;
747
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));
758
759 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_elemental_elemental));
760 *F = B;
761 PetscFunctionReturn(PETSC_SUCCESS);
762 }
763
MatSolverTypeRegister_Elemental(void)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 }
771
MatNorm_Elemental(Mat A,NormType type,PetscReal * nrm)772 static PetscErrorCode MatNorm_Elemental(Mat A, NormType type, PetscReal *nrm)
773 {
774 Mat_Elemental *a = (Mat_Elemental *)A->data;
775
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 }
792
MatZeroEntries_Elemental(Mat A)793 static PetscErrorCode MatZeroEntries_Elemental(Mat A)
794 {
795 Mat_Elemental *a = (Mat_Elemental *)A->data;
796
797 PetscFunctionBegin;
798 El::Zero(*a->emat);
799 PetscFunctionReturn(PETSC_SUCCESS);
800 }
801
MatGetOwnershipIS_Elemental(Mat A,IS * rows,IS * cols)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;
806
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 }
834
MatConvert_Elemental_Dense(Mat A,MatType,MatReuse reuse,Mat * B)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();
845
846 PetscFunctionBegin;
847 PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
848
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 }
857
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));
864
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");
874
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");
890
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 }
909
MatConvert_SeqAIJ_Elemental(Mat A,MatType,MatReuse reuse,Mat * newmat)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;
916
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));
935
936 if (reuse == MAT_INPLACE_MATRIX) {
937 PetscCall(MatHeaderReplace(A, &mat_elemental));
938 } else {
939 *newmat = mat_elemental;
940 }
941 PetscFunctionReturn(PETSC_SUCCESS);
942 }
943
MatConvert_MPIAIJ_Elemental(Mat A,MatType,MatReuse reuse,Mat * newmat)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;
950
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));
971
972 if (reuse == MAT_INPLACE_MATRIX) {
973 PetscCall(MatHeaderReplace(A, &mat_elemental));
974 } else {
975 *newmat = mat_elemental;
976 }
977 PetscFunctionReturn(PETSC_SUCCESS);
978 }
979
MatConvert_SeqSBAIJ_Elemental(Mat A,MatType,MatReuse reuse,Mat * newmat)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;
986
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));
1013
1014 if (reuse == MAT_INPLACE_MATRIX) {
1015 PetscCall(MatHeaderReplace(A, &mat_elemental));
1016 } else {
1017 *newmat = mat_elemental;
1018 }
1019 PetscFunctionReturn(PETSC_SUCCESS);
1020 }
1021
MatConvert_MPISBAIJ_Elemental(Mat A,MatType,MatReuse reuse,Mat * newmat)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;
1028
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));
1055
1056 if (reuse == MAT_INPLACE_MATRIX) {
1057 PetscCall(MatHeaderReplace(A, &mat_elemental));
1058 } else {
1059 *newmat = mat_elemental;
1060 }
1061 PetscFunctionReturn(PETSC_SUCCESS);
1062 }
1063
MatDestroy_Elemental(Mat A)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;
1070
1071 PetscFunctionBegin;
1072 delete a->emat;
1073 delete a->P;
1074 delete a->Q;
1075
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 }
1091
MatSetUp_Elemental(Mat A)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;
1098
1099 PetscFunctionBegin;
1100 PetscCall(PetscLayoutSetUp(A->rmap));
1101 PetscCall(PetscLayoutSetUp(A->cmap));
1102
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);
1111
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);
1115
1116 a->emat->Resize(A->rmap->N, A->cmap->N);
1117 El::Zero(*a->emat);
1118
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 }
1131
MatAssemblyBegin_Elemental(Mat A,MatAssemblyType)1132 static PetscErrorCode MatAssemblyBegin_Elemental(Mat A, MatAssemblyType)
1133 {
1134 Mat_Elemental *a = (Mat_Elemental *)A->data;
1135
1136 PetscFunctionBegin;
1137 /* printf("Calling ProcessQueues\n"); */
1138 a->emat->ProcessQueues();
1139 /* printf("Finished ProcessQueues\n"); */
1140 PetscFunctionReturn(PETSC_SUCCESS);
1141 }
1142
MatAssemblyEnd_Elemental(Mat,MatAssemblyType)1143 static PetscErrorCode MatAssemblyEnd_Elemental(Mat, MatAssemblyType)
1144 {
1145 PetscFunctionBegin;
1146 /* Currently does nothing */
1147 PetscFunctionReturn(PETSC_SUCCESS);
1148 }
1149
MatLoad_Elemental(Mat newMat,PetscViewer viewer)1150 static PetscErrorCode MatLoad_Elemental(Mat newMat, PetscViewer viewer)
1151 {
1152 Mat Adense, Ae;
1153 MPI_Comm comm;
1154
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 }
1165
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};
1310
1311 /*MC
1312 MATELEMENTAL = "elemental" - A matrix type for dense matrices using the Elemental package
1313
1314 Use ./configure --download-elemental to install PETSc to use Elemental
1315
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
1320
1321 Level: beginner
1322
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.
1327
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
MatCreate_Elemental(Mat A)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;
1342
1343 PetscFunctionBegin;
1344 A->ops[0] = MatOps_Values;
1345 A->insertmode = NOT_SET_VALUES;
1346
1347 PetscCall(PetscNew(&a));
1348 A->data = (void *)a;
1349
1350 /* Set up the elemental matrix */
1351 El::mpi::Comm cxxcomm(PetscObjectComm((PetscObject)A));
1352
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));
1362
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));
1375
1376 a->pivoting = 1;
1377 PetscCall(PetscOptionsInt("-mat_elemental_pivoting", "Pivoting", "None", a->pivoting, &a->pivoting, NULL));
1378
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;
1387
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
1396