Actual source code: matnest.c
1: #include <../src/mat/impls/nest/matnestimpl.h>
2: #include <../src/mat/impls/aij/seq/aij.h>
3: #include <petscsf.h>
5: static PetscErrorCode MatSetUp_NestIS_Private(Mat, PetscInt, const IS[], PetscInt, const IS[]);
6: static PetscErrorCode MatCreateVecs_Nest(Mat, Vec *, Vec *);
7: static PetscErrorCode MatReset_Nest(Mat);
9: PETSC_INTERN PetscErrorCode MatConvert_Nest_IS(Mat, MatType, MatReuse, Mat *);
11: /* private functions */
12: static PetscErrorCode MatNestGetSizes_Private(Mat A, PetscInt *m, PetscInt *n, PetscInt *M, PetscInt *N)
13: {
14: Mat_Nest *bA = (Mat_Nest *)A->data;
15: PetscInt i, j;
17: PetscFunctionBegin;
18: *m = *n = *M = *N = 0;
19: for (i = 0; i < bA->nr; i++) { /* rows */
20: PetscInt sm, sM;
21: PetscCall(ISGetLocalSize(bA->isglobal.row[i], &sm));
22: PetscCall(ISGetSize(bA->isglobal.row[i], &sM));
23: *m += sm;
24: *M += sM;
25: }
26: for (j = 0; j < bA->nc; j++) { /* cols */
27: PetscInt sn, sN;
28: PetscCall(ISGetLocalSize(bA->isglobal.col[j], &sn));
29: PetscCall(ISGetSize(bA->isglobal.col[j], &sN));
30: *n += sn;
31: *N += sN;
32: }
33: PetscFunctionReturn(PETSC_SUCCESS);
34: }
36: /* operations */
37: static PetscErrorCode MatMult_Nest(Mat A, Vec x, Vec y)
38: {
39: Mat_Nest *bA = (Mat_Nest *)A->data;
40: Vec *bx = bA->right, *by = bA->left;
41: PetscInt i, j, nr = bA->nr, nc = bA->nc;
43: PetscFunctionBegin;
44: for (i = 0; i < nr; i++) PetscCall(VecGetSubVector(y, bA->isglobal.row[i], &by[i]));
45: for (i = 0; i < nc; i++) PetscCall(VecGetSubVector(x, bA->isglobal.col[i], &bx[i]));
46: for (i = 0; i < nr; i++) {
47: PetscCall(VecZeroEntries(by[i]));
48: for (j = 0; j < nc; j++) {
49: if (!bA->m[i][j]) continue;
50: /* y[i] <- y[i] + A[i][j] * x[j] */
51: PetscCall(MatMultAdd(bA->m[i][j], bx[j], by[i], by[i]));
52: }
53: }
54: for (i = 0; i < nr; i++) PetscCall(VecRestoreSubVector(y, bA->isglobal.row[i], &by[i]));
55: for (i = 0; i < nc; i++) PetscCall(VecRestoreSubVector(x, bA->isglobal.col[i], &bx[i]));
56: PetscFunctionReturn(PETSC_SUCCESS);
57: }
59: static PetscErrorCode MatMultAdd_Nest(Mat A, Vec x, Vec y, Vec z)
60: {
61: Mat_Nest *bA = (Mat_Nest *)A->data;
62: Vec *bx = bA->right, *bz = bA->left;
63: PetscInt i, j, nr = bA->nr, nc = bA->nc;
65: PetscFunctionBegin;
66: for (i = 0; i < nr; i++) PetscCall(VecGetSubVector(z, bA->isglobal.row[i], &bz[i]));
67: for (i = 0; i < nc; i++) PetscCall(VecGetSubVector(x, bA->isglobal.col[i], &bx[i]));
68: for (i = 0; i < nr; i++) {
69: if (y != z) {
70: Vec by;
71: PetscCall(VecGetSubVector(y, bA->isglobal.row[i], &by));
72: PetscCall(VecCopy(by, bz[i]));
73: PetscCall(VecRestoreSubVector(y, bA->isglobal.row[i], &by));
74: }
75: for (j = 0; j < nc; j++) {
76: if (!bA->m[i][j]) continue;
77: /* y[i] <- y[i] + A[i][j] * x[j] */
78: PetscCall(MatMultAdd(bA->m[i][j], bx[j], bz[i], bz[i]));
79: }
80: }
81: for (i = 0; i < nr; i++) PetscCall(VecRestoreSubVector(z, bA->isglobal.row[i], &bz[i]));
82: for (i = 0; i < nc; i++) PetscCall(VecRestoreSubVector(x, bA->isglobal.col[i], &bx[i]));
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: typedef struct {
87: Mat *workC; /* array of Mat with specific containers depending on the underlying MatMatMult implementation */
88: PetscScalar *tarray; /* buffer for storing all temporary products A[i][j] B[j] */
89: PetscInt *dm, *dn, k; /* displacements and number of submatrices */
90: } Nest_Dense;
92: static PetscErrorCode MatProductNumeric_Nest_Dense(Mat C)
93: {
94: Mat_Nest *bA;
95: Nest_Dense *contents;
96: Mat viewB, viewC, productB, workC;
97: const PetscScalar *barray;
98: PetscScalar *carray;
99: PetscInt i, j, M, N, nr, nc, ldb, ldc;
100: Mat A, B;
102: PetscFunctionBegin;
103: MatCheckProduct(C, 1);
104: A = C->product->A;
105: B = C->product->B;
106: PetscCall(MatGetSize(B, NULL, &N));
107: if (!N) {
108: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
109: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
110: PetscFunctionReturn(PETSC_SUCCESS);
111: }
112: contents = (Nest_Dense *)C->product->data;
113: PetscCheck(contents, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data empty");
114: bA = (Mat_Nest *)A->data;
115: nr = bA->nr;
116: nc = bA->nc;
117: PetscCall(MatDenseGetLDA(B, &ldb));
118: PetscCall(MatDenseGetLDA(C, &ldc));
119: PetscCall(MatZeroEntries(C));
120: PetscCall(MatDenseGetArrayRead(B, &barray));
121: PetscCall(MatDenseGetArray(C, &carray));
122: for (i = 0; i < nr; i++) {
123: PetscCall(ISGetSize(bA->isglobal.row[i], &M));
124: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A), contents->dm[i + 1] - contents->dm[i], PETSC_DECIDE, M, N, PetscSafePointerPlusOffset(carray, contents->dm[i]), &viewC));
125: PetscCall(MatDenseSetLDA(viewC, ldc));
126: for (j = 0; j < nc; j++) {
127: if (!bA->m[i][j]) continue;
128: PetscCall(ISGetSize(bA->isglobal.col[j], &M));
129: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A), contents->dn[j + 1] - contents->dn[j], PETSC_DECIDE, M, N, PetscSafePointerPlusOffset((PetscScalar *)barray, contents->dn[j]), &viewB));
130: PetscCall(MatDenseSetLDA(viewB, ldb));
132: /* MatMatMultNumeric(bA->m[i][j],viewB,contents->workC[i*nc + j]); */
133: workC = contents->workC[i * nc + j];
134: productB = workC->product->B;
135: workC->product->B = viewB; /* use newly created dense matrix viewB */
136: PetscCall(MatProductNumeric(workC));
137: PetscCall(MatDestroy(&viewB));
138: workC->product->B = productB; /* resume original B */
140: /* C[i] <- workC + C[i] */
141: PetscCall(MatAXPY(viewC, 1.0, contents->workC[i * nc + j], SAME_NONZERO_PATTERN));
142: }
143: PetscCall(MatDestroy(&viewC));
144: }
145: PetscCall(MatDenseRestoreArray(C, &carray));
146: PetscCall(MatDenseRestoreArrayRead(B, &barray));
148: PetscCall(MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
149: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
150: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
151: PetscFunctionReturn(PETSC_SUCCESS);
152: }
154: static PetscErrorCode MatNest_DenseDestroy(void *ctx)
155: {
156: Nest_Dense *contents = (Nest_Dense *)ctx;
157: PetscInt i;
159: PetscFunctionBegin;
160: PetscCall(PetscFree(contents->tarray));
161: for (i = 0; i < contents->k; i++) PetscCall(MatDestroy(contents->workC + i));
162: PetscCall(PetscFree3(contents->dm, contents->dn, contents->workC));
163: PetscCall(PetscFree(contents));
164: PetscFunctionReturn(PETSC_SUCCESS);
165: }
167: static PetscErrorCode MatProductSymbolic_Nest_Dense(Mat C)
168: {
169: Mat_Nest *bA;
170: Mat viewB, workC;
171: const PetscScalar *barray;
172: PetscInt i, j, M, N, m, n, nr, nc, maxm = 0, ldb;
173: Nest_Dense *contents = NULL;
174: PetscBool cisdense;
175: Mat A, B;
176: PetscReal fill;
178: PetscFunctionBegin;
179: MatCheckProduct(C, 1);
180: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
181: A = C->product->A;
182: B = C->product->B;
183: fill = C->product->fill;
184: bA = (Mat_Nest *)A->data;
185: nr = bA->nr;
186: nc = bA->nc;
187: PetscCall(MatGetLocalSize(C, &m, &n));
188: PetscCall(MatGetSize(C, &M, &N));
189: if (m == PETSC_DECIDE || n == PETSC_DECIDE || M == PETSC_DECIDE || N == PETSC_DECIDE) {
190: PetscCall(MatGetLocalSize(B, NULL, &n));
191: PetscCall(MatGetSize(B, NULL, &N));
192: PetscCall(MatGetLocalSize(A, &m, NULL));
193: PetscCall(MatGetSize(A, &M, NULL));
194: PetscCall(MatSetSizes(C, m, n, M, N));
195: }
196: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATMPIDENSE, MATSEQDENSECUDA, MATMPIDENSECUDA, ""));
197: if (!cisdense) PetscCall(MatSetType(C, ((PetscObject)B)->type_name));
198: PetscCall(MatSetUp(C));
199: if (!N) {
200: C->ops->productnumeric = MatProductNumeric_Nest_Dense;
201: PetscFunctionReturn(PETSC_SUCCESS);
202: }
204: PetscCall(PetscNew(&contents));
205: C->product->data = contents;
206: C->product->destroy = MatNest_DenseDestroy;
207: PetscCall(PetscCalloc3(nr + 1, &contents->dm, nc + 1, &contents->dn, nr * nc, &contents->workC));
208: contents->k = nr * nc;
209: for (i = 0; i < nr; i++) {
210: PetscCall(ISGetLocalSize(bA->isglobal.row[i], contents->dm + i + 1));
211: maxm = PetscMax(maxm, contents->dm[i + 1]);
212: contents->dm[i + 1] += contents->dm[i];
213: }
214: for (i = 0; i < nc; i++) {
215: PetscCall(ISGetLocalSize(bA->isglobal.col[i], contents->dn + i + 1));
216: contents->dn[i + 1] += contents->dn[i];
217: }
218: PetscCall(PetscMalloc1(maxm * N, &contents->tarray));
219: PetscCall(MatDenseGetLDA(B, &ldb));
220: PetscCall(MatGetSize(B, NULL, &N));
221: PetscCall(MatDenseGetArrayRead(B, &barray));
222: /* loops are permuted compared to MatMatMultNumeric so that viewB is created only once per column of A */
223: for (j = 0; j < nc; j++) {
224: PetscCall(ISGetSize(bA->isglobal.col[j], &M));
225: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A), contents->dn[j + 1] - contents->dn[j], PETSC_DECIDE, M, N, PetscSafePointerPlusOffset((PetscScalar *)barray, contents->dn[j]), &viewB));
226: PetscCall(MatDenseSetLDA(viewB, ldb));
227: for (i = 0; i < nr; i++) {
228: if (!bA->m[i][j]) continue;
229: /* MatMatMultSymbolic may attach a specific container (depending on MatType of bA->m[i][j]) to workC[i][j] */
231: PetscCall(MatProductCreate(bA->m[i][j], viewB, NULL, &contents->workC[i * nc + j]));
232: workC = contents->workC[i * nc + j];
233: PetscCall(MatProductSetType(workC, MATPRODUCT_AB));
234: PetscCall(MatProductSetAlgorithm(workC, "default"));
235: PetscCall(MatProductSetFill(workC, fill));
236: PetscCall(MatProductSetFromOptions(workC));
237: PetscCall(MatProductSymbolic(workC));
239: /* since tarray will be shared by all Mat */
240: PetscCall(MatSeqDenseSetPreallocation(workC, contents->tarray));
241: PetscCall(MatMPIDenseSetPreallocation(workC, contents->tarray));
242: }
243: PetscCall(MatDestroy(&viewB));
244: }
245: PetscCall(MatDenseRestoreArrayRead(B, &barray));
247: C->ops->productnumeric = MatProductNumeric_Nest_Dense;
248: PetscFunctionReturn(PETSC_SUCCESS);
249: }
251: static PetscErrorCode MatProductSetFromOptions_Nest_Dense(Mat C)
252: {
253: Mat_Product *product = C->product;
255: PetscFunctionBegin;
256: if (product->type == MATPRODUCT_AB) C->ops->productsymbolic = MatProductSymbolic_Nest_Dense;
257: PetscFunctionReturn(PETSC_SUCCESS);
258: }
260: static PetscErrorCode MatMultTransposeKernel_Nest(Mat A, Vec x, Vec y, PetscBool herm)
261: {
262: Mat_Nest *bA = (Mat_Nest *)A->data;
263: Vec *bx = bA->left, *by = bA->right;
264: PetscInt i, j, nr = bA->nr, nc = bA->nc;
266: PetscFunctionBegin;
267: for (i = 0; i < nr; i++) PetscCall(VecGetSubVector(x, bA->isglobal.row[i], &bx[i]));
268: for (i = 0; i < nc; i++) PetscCall(VecGetSubVector(y, bA->isglobal.col[i], &by[i]));
269: for (j = 0; j < nc; j++) {
270: PetscCall(VecZeroEntries(by[j]));
271: for (i = 0; i < nr; i++) {
272: if (!bA->m[i][j]) continue;
273: if (herm) PetscCall(MatMultHermitianTransposeAdd(bA->m[i][j], bx[i], by[j], by[j])); /* y[j] <- y[j] + (A[i][j])^H * x[i] */
274: else PetscCall(MatMultTransposeAdd(bA->m[i][j], bx[i], by[j], by[j])); /* y[j] <- y[j] + (A[i][j])^T * x[i] */
275: }
276: }
277: for (i = 0; i < nr; i++) PetscCall(VecRestoreSubVector(x, bA->isglobal.row[i], &bx[i]));
278: for (i = 0; i < nc; i++) PetscCall(VecRestoreSubVector(y, bA->isglobal.col[i], &by[i]));
279: PetscFunctionReturn(PETSC_SUCCESS);
280: }
282: static PetscErrorCode MatMultTranspose_Nest(Mat A, Vec x, Vec y)
283: {
284: PetscFunctionBegin;
285: PetscCall(MatMultTransposeKernel_Nest(A, x, y, PETSC_FALSE));
286: PetscFunctionReturn(PETSC_SUCCESS);
287: }
289: static PetscErrorCode MatMultHermitianTranspose_Nest(Mat A, Vec x, Vec y)
290: {
291: PetscFunctionBegin;
292: PetscCall(MatMultTransposeKernel_Nest(A, x, y, PETSC_TRUE));
293: PetscFunctionReturn(PETSC_SUCCESS);
294: }
296: static PetscErrorCode MatMultTransposeAddKernel_Nest(Mat A, Vec x, Vec y, Vec z, PetscBool herm)
297: {
298: Mat_Nest *bA = (Mat_Nest *)A->data;
299: Vec *bx = bA->left, *bz = bA->right;
300: PetscInt i, j, nr = bA->nr, nc = bA->nc;
302: PetscFunctionBegin;
303: for (i = 0; i < nr; i++) PetscCall(VecGetSubVector(x, bA->isglobal.row[i], &bx[i]));
304: for (i = 0; i < nc; i++) PetscCall(VecGetSubVector(z, bA->isglobal.col[i], &bz[i]));
305: for (j = 0; j < nc; j++) {
306: if (y != z) {
307: Vec by;
308: PetscCall(VecGetSubVector(y, bA->isglobal.col[j], &by));
309: PetscCall(VecCopy(by, bz[j]));
310: PetscCall(VecRestoreSubVector(y, bA->isglobal.col[j], &by));
311: }
312: for (i = 0; i < nr; i++) {
313: if (!bA->m[i][j]) continue;
314: if (herm) PetscCall(MatMultHermitianTransposeAdd(bA->m[i][j], bx[i], bz[j], bz[j])); /* z[j] <- y[j] + (A[i][j])^H * x[i] */
315: else PetscCall(MatMultTransposeAdd(bA->m[i][j], bx[i], bz[j], bz[j])); /* z[j] <- y[j] + (A[i][j])^T * x[i] */
316: }
317: }
318: for (i = 0; i < nr; i++) PetscCall(VecRestoreSubVector(x, bA->isglobal.row[i], &bx[i]));
319: for (i = 0; i < nc; i++) PetscCall(VecRestoreSubVector(z, bA->isglobal.col[i], &bz[i]));
320: PetscFunctionReturn(PETSC_SUCCESS);
321: }
323: static PetscErrorCode MatMultTransposeAdd_Nest(Mat A, Vec x, Vec y, Vec z)
324: {
325: PetscFunctionBegin;
326: PetscCall(MatMultTransposeAddKernel_Nest(A, x, y, z, PETSC_FALSE));
327: PetscFunctionReturn(PETSC_SUCCESS);
328: }
330: static PetscErrorCode MatMultHermitianTransposeAdd_Nest(Mat A, Vec x, Vec y, Vec z)
331: {
332: PetscFunctionBegin;
333: PetscCall(MatMultTransposeAddKernel_Nest(A, x, y, z, PETSC_TRUE));
334: PetscFunctionReturn(PETSC_SUCCESS);
335: }
337: static PetscErrorCode MatTranspose_Nest(Mat A, MatReuse reuse, Mat *B)
338: {
339: Mat_Nest *bA = (Mat_Nest *)A->data, *bC;
340: Mat C;
341: PetscInt i, j, nr = bA->nr, nc = bA->nc;
343: PetscFunctionBegin;
344: if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *B));
345: PetscCheck(reuse != MAT_INPLACE_MATRIX || nr == nc, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_SIZ, "Square nested matrix only for in-place");
347: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
348: Mat *subs;
349: IS *is_row, *is_col;
351: PetscCall(PetscCalloc1(nr * nc, &subs));
352: PetscCall(PetscMalloc2(nr, &is_row, nc, &is_col));
353: PetscCall(MatNestGetISs(A, is_row, is_col));
354: if (reuse == MAT_INPLACE_MATRIX) {
355: for (i = 0; i < nr; i++) {
356: for (j = 0; j < nc; j++) subs[i + nr * j] = bA->m[i][j];
357: }
358: }
360: PetscCall(MatCreateNest(PetscObjectComm((PetscObject)A), nc, is_col, nr, is_row, subs, &C));
361: PetscCall(PetscFree(subs));
362: PetscCall(PetscFree2(is_row, is_col));
363: } else {
364: C = *B;
365: }
367: bC = (Mat_Nest *)C->data;
368: for (i = 0; i < nr; i++) {
369: for (j = 0; j < nc; j++) {
370: if (bA->m[i][j]) {
371: PetscCall(MatTranspose(bA->m[i][j], reuse, &bC->m[j][i]));
372: } else {
373: bC->m[j][i] = NULL;
374: }
375: }
376: }
378: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
379: *B = C;
380: } else {
381: PetscCall(MatHeaderMerge(A, &C));
382: }
383: PetscFunctionReturn(PETSC_SUCCESS);
384: }
386: static PetscErrorCode MatNestDestroyISList(PetscInt n, IS **list)
387: {
388: IS *lst = *list;
389: PetscInt i;
391: PetscFunctionBegin;
392: if (!lst) PetscFunctionReturn(PETSC_SUCCESS);
393: for (i = 0; i < n; i++)
394: if (lst[i]) PetscCall(ISDestroy(&lst[i]));
395: PetscCall(PetscFree(lst));
396: *list = NULL;
397: PetscFunctionReturn(PETSC_SUCCESS);
398: }
400: static PetscErrorCode MatReset_Nest(Mat A)
401: {
402: Mat_Nest *vs = (Mat_Nest *)A->data;
403: PetscInt i, j;
405: PetscFunctionBegin;
406: /* release the matrices and the place holders */
407: PetscCall(MatNestDestroyISList(vs->nr, &vs->isglobal.row));
408: PetscCall(MatNestDestroyISList(vs->nc, &vs->isglobal.col));
409: PetscCall(MatNestDestroyISList(vs->nr, &vs->islocal.row));
410: PetscCall(MatNestDestroyISList(vs->nc, &vs->islocal.col));
412: PetscCall(PetscFree(vs->row_len));
413: PetscCall(PetscFree(vs->col_len));
414: PetscCall(PetscFree(vs->nnzstate));
416: PetscCall(PetscFree2(vs->left, vs->right));
418: /* release the matrices and the place holders */
419: if (vs->m) {
420: for (i = 0; i < vs->nr; i++) {
421: for (j = 0; j < vs->nc; j++) PetscCall(MatDestroy(&vs->m[i][j]));
422: }
423: PetscCall(PetscFree(vs->m[0]));
424: PetscCall(PetscFree(vs->m));
425: }
427: /* restore defaults */
428: vs->nr = 0;
429: vs->nc = 0;
430: vs->splitassembly = PETSC_FALSE;
431: PetscFunctionReturn(PETSC_SUCCESS);
432: }
434: static PetscErrorCode MatDestroy_Nest(Mat A)
435: {
436: PetscFunctionBegin;
437: PetscCall(MatReset_Nest(A));
438: PetscCall(PetscFree(A->data));
439: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSubMat_C", NULL));
440: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetSubMat_C", NULL));
441: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSubMats_C", NULL));
442: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSize_C", NULL));
443: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetISs_C", NULL));
444: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetLocalISs_C", NULL));
445: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetVecType_C", NULL));
446: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetSubMats_C", NULL));
447: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_mpiaij_C", NULL));
448: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_seqaij_C", NULL));
449: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_aij_C", NULL));
450: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_is_C", NULL));
451: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_mpidense_C", NULL));
452: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_seqdense_C", NULL));
453: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_seqdense_C", NULL));
454: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_mpidense_C", NULL));
455: PetscFunctionReturn(PETSC_SUCCESS);
456: }
458: static PetscErrorCode MatMissingDiagonal_Nest(Mat mat, PetscBool *missing, PetscInt *dd)
459: {
460: Mat_Nest *vs = (Mat_Nest *)mat->data;
461: PetscInt i;
463: PetscFunctionBegin;
464: if (dd) *dd = 0;
465: if (!vs->nr) {
466: *missing = PETSC_TRUE;
467: PetscFunctionReturn(PETSC_SUCCESS);
468: }
469: *missing = PETSC_FALSE;
470: for (i = 0; i < vs->nr && !(*missing); i++) {
471: *missing = PETSC_TRUE;
472: if (vs->m[i][i]) {
473: PetscCall(MatMissingDiagonal(vs->m[i][i], missing, NULL));
474: PetscCheck(!*missing || !dd, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "First missing entry not yet implemented");
475: }
476: }
477: PetscFunctionReturn(PETSC_SUCCESS);
478: }
480: static PetscErrorCode MatAssemblyBegin_Nest(Mat A, MatAssemblyType type)
481: {
482: Mat_Nest *vs = (Mat_Nest *)A->data;
483: PetscInt i, j;
484: PetscBool nnzstate = PETSC_FALSE;
486: PetscFunctionBegin;
487: for (i = 0; i < vs->nr; i++) {
488: for (j = 0; j < vs->nc; j++) {
489: PetscObjectState subnnzstate = 0;
490: if (vs->m[i][j]) {
491: PetscCall(MatAssemblyBegin(vs->m[i][j], type));
492: if (!vs->splitassembly) {
493: /* Note: split assembly will fail if the same block appears more than once (even indirectly through a nested
494: * sub-block). This could be fixed by adding a flag to Mat so that there was a way to check if a Mat was
495: * already performing an assembly, but the result would by more complicated and appears to offer less
496: * potential for diagnostics and correctness checking. Split assembly should be fixed once there is an
497: * interface for libraries to make asynchronous progress in "user-defined non-blocking collectives".
498: */
499: PetscCall(MatAssemblyEnd(vs->m[i][j], type));
500: PetscCall(MatGetNonzeroState(vs->m[i][j], &subnnzstate));
501: }
502: }
503: nnzstate = (PetscBool)(nnzstate || vs->nnzstate[i * vs->nc + j] != subnnzstate);
504: vs->nnzstate[i * vs->nc + j] = subnnzstate;
505: }
506: }
507: if (nnzstate) A->nonzerostate++;
508: PetscFunctionReturn(PETSC_SUCCESS);
509: }
511: static PetscErrorCode MatAssemblyEnd_Nest(Mat A, MatAssemblyType type)
512: {
513: Mat_Nest *vs = (Mat_Nest *)A->data;
514: PetscInt i, j;
516: PetscFunctionBegin;
517: for (i = 0; i < vs->nr; i++) {
518: for (j = 0; j < vs->nc; j++) {
519: if (vs->m[i][j]) {
520: if (vs->splitassembly) PetscCall(MatAssemblyEnd(vs->m[i][j], type));
521: }
522: }
523: }
524: PetscFunctionReturn(PETSC_SUCCESS);
525: }
527: static PetscErrorCode MatNestFindNonzeroSubMatRow(Mat A, PetscInt row, Mat *B)
528: {
529: Mat_Nest *vs = (Mat_Nest *)A->data;
530: PetscInt j;
531: Mat sub;
533: PetscFunctionBegin;
534: sub = (row < vs->nc) ? vs->m[row][row] : (Mat)NULL; /* Prefer to find on the diagonal */
535: for (j = 0; !sub && j < vs->nc; j++) sub = vs->m[row][j];
536: if (sub) PetscCall(MatSetUp(sub)); /* Ensure that the sizes are available */
537: *B = sub;
538: PetscFunctionReturn(PETSC_SUCCESS);
539: }
541: static PetscErrorCode MatNestFindNonzeroSubMatCol(Mat A, PetscInt col, Mat *B)
542: {
543: Mat_Nest *vs = (Mat_Nest *)A->data;
544: PetscInt i;
545: Mat sub;
547: PetscFunctionBegin;
548: sub = (col < vs->nr) ? vs->m[col][col] : (Mat)NULL; /* Prefer to find on the diagonal */
549: for (i = 0; !sub && i < vs->nr; i++) sub = vs->m[i][col];
550: if (sub) PetscCall(MatSetUp(sub)); /* Ensure that the sizes are available */
551: *B = sub;
552: PetscFunctionReturn(PETSC_SUCCESS);
553: }
555: static PetscErrorCode MatNestFindISRange(Mat A, PetscInt n, const IS list[], IS is, PetscInt *begin, PetscInt *end)
556: {
557: PetscInt i, j, size, m;
558: PetscBool flg;
559: IS out, concatenate[2];
561: PetscFunctionBegin;
562: PetscAssertPointer(list, 3);
564: if (begin) {
565: PetscAssertPointer(begin, 5);
566: *begin = -1;
567: }
568: if (end) {
569: PetscAssertPointer(end, 6);
570: *end = -1;
571: }
572: for (i = 0; i < n; i++) {
573: if (!list[i]) continue;
574: PetscCall(ISEqualUnsorted(list[i], is, &flg));
575: if (flg) {
576: if (begin) *begin = i;
577: if (end) *end = i + 1;
578: PetscFunctionReturn(PETSC_SUCCESS);
579: }
580: }
581: PetscCall(ISGetSize(is, &size));
582: for (i = 0; i < n - 1; i++) {
583: if (!list[i]) continue;
584: m = 0;
585: PetscCall(ISConcatenate(PetscObjectComm((PetscObject)A), 2, list + i, &out));
586: PetscCall(ISGetSize(out, &m));
587: for (j = i + 2; j < n && m < size; j++) {
588: if (list[j]) {
589: concatenate[0] = out;
590: concatenate[1] = list[j];
591: PetscCall(ISConcatenate(PetscObjectComm((PetscObject)A), 2, concatenate, &out));
592: PetscCall(ISDestroy(concatenate));
593: PetscCall(ISGetSize(out, &m));
594: }
595: }
596: if (m == size) {
597: PetscCall(ISEqualUnsorted(out, is, &flg));
598: if (flg) {
599: if (begin) *begin = i;
600: if (end) *end = j;
601: PetscCall(ISDestroy(&out));
602: PetscFunctionReturn(PETSC_SUCCESS);
603: }
604: }
605: PetscCall(ISDestroy(&out));
606: }
607: PetscFunctionReturn(PETSC_SUCCESS);
608: }
610: static PetscErrorCode MatNestFillEmptyMat_Private(Mat A, PetscInt i, PetscInt j, Mat *B)
611: {
612: Mat_Nest *vs = (Mat_Nest *)A->data;
613: PetscInt lr, lc;
615: PetscFunctionBegin;
616: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B));
617: PetscCall(ISGetLocalSize(vs->isglobal.row[i], &lr));
618: PetscCall(ISGetLocalSize(vs->isglobal.col[j], &lc));
619: PetscCall(MatSetSizes(*B, lr, lc, PETSC_DECIDE, PETSC_DECIDE));
620: PetscCall(MatSetType(*B, MATAIJ));
621: PetscCall(MatSeqAIJSetPreallocation(*B, 0, NULL));
622: PetscCall(MatMPIAIJSetPreallocation(*B, 0, NULL, 0, NULL));
623: PetscCall(MatSetUp(*B));
624: PetscCall(MatSetOption(*B, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
625: PetscCall(MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY));
626: PetscCall(MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY));
627: PetscFunctionReturn(PETSC_SUCCESS);
628: }
630: static PetscErrorCode MatNestGetBlock_Private(Mat A, PetscInt rbegin, PetscInt rend, PetscInt cbegin, PetscInt cend, Mat *B)
631: {
632: Mat_Nest *vs = (Mat_Nest *)A->data;
633: Mat *a;
634: PetscInt i, j, k, l, nr = rend - rbegin, nc = cend - cbegin;
635: char keyname[256];
636: PetscBool *b;
637: PetscBool flg;
639: PetscFunctionBegin;
640: *B = NULL;
641: PetscCall(PetscSNPrintf(keyname, sizeof(keyname), "NestBlock_%" PetscInt_FMT "-%" PetscInt_FMT "x%" PetscInt_FMT "-%" PetscInt_FMT, rbegin, rend, cbegin, cend));
642: PetscCall(PetscObjectQuery((PetscObject)A, keyname, (PetscObject *)B));
643: if (*B) PetscFunctionReturn(PETSC_SUCCESS);
645: PetscCall(PetscMalloc2(nr * nc, &a, nr * nc, &b));
646: for (i = 0; i < nr; i++) {
647: for (j = 0; j < nc; j++) {
648: a[i * nc + j] = vs->m[rbegin + i][cbegin + j];
649: b[i * nc + j] = PETSC_FALSE;
650: }
651: }
652: if (nc != vs->nc && nr != vs->nr) {
653: for (i = 0; i < nr; i++) {
654: for (j = 0; j < nc; j++) {
655: flg = PETSC_FALSE;
656: for (k = 0; (k < nr && !flg); k++) {
657: if (a[j + k * nc]) flg = PETSC_TRUE;
658: }
659: if (flg) {
660: flg = PETSC_FALSE;
661: for (l = 0; (l < nc && !flg); l++) {
662: if (a[i * nc + l]) flg = PETSC_TRUE;
663: }
664: }
665: if (!flg) {
666: b[i * nc + j] = PETSC_TRUE;
667: PetscCall(MatNestFillEmptyMat_Private(A, rbegin + i, cbegin + j, a + i * nc + j));
668: }
669: }
670: }
671: }
672: PetscCall(MatCreateNest(PetscObjectComm((PetscObject)A), nr, nr != vs->nr ? NULL : vs->isglobal.row, nc, nc != vs->nc ? NULL : vs->isglobal.col, a, B));
673: for (i = 0; i < nr; i++) {
674: for (j = 0; j < nc; j++) {
675: if (b[i * nc + j]) PetscCall(MatDestroy(a + i * nc + j));
676: }
677: }
678: PetscCall(PetscFree2(a, b));
679: (*B)->assembled = A->assembled;
680: PetscCall(PetscObjectCompose((PetscObject)A, keyname, (PetscObject)*B));
681: PetscCall(PetscObjectDereference((PetscObject)*B)); /* Leave the only remaining reference in the composition */
682: PetscFunctionReturn(PETSC_SUCCESS);
683: }
685: static PetscErrorCode MatNestFindSubMat(Mat A, struct MatNestISPair *is, IS isrow, IS iscol, Mat *B)
686: {
687: Mat_Nest *vs = (Mat_Nest *)A->data;
688: PetscInt rbegin, rend, cbegin, cend;
690: PetscFunctionBegin;
691: PetscCall(MatNestFindISRange(A, vs->nr, is->row, isrow, &rbegin, &rend));
692: PetscCall(MatNestFindISRange(A, vs->nc, is->col, iscol, &cbegin, &cend));
693: if (rend == rbegin + 1 && cend == cbegin + 1) {
694: if (!vs->m[rbegin][cbegin]) PetscCall(MatNestFillEmptyMat_Private(A, rbegin, cbegin, vs->m[rbegin] + cbegin));
695: *B = vs->m[rbegin][cbegin];
696: } else if (rbegin != -1 && cbegin != -1) {
697: PetscCall(MatNestGetBlock_Private(A, rbegin, rend, cbegin, cend, B));
698: } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Could not find index set");
699: PetscFunctionReturn(PETSC_SUCCESS);
700: }
702: /*
703: TODO: This does not actually returns a submatrix we can modify
704: */
705: static PetscErrorCode MatCreateSubMatrix_Nest(Mat A, IS isrow, IS iscol, MatReuse reuse, Mat *B)
706: {
707: Mat_Nest *vs = (Mat_Nest *)A->data;
708: Mat sub;
710: PetscFunctionBegin;
711: PetscCall(MatNestFindSubMat(A, &vs->isglobal, isrow, iscol, &sub));
712: switch (reuse) {
713: case MAT_INITIAL_MATRIX:
714: if (sub) PetscCall(PetscObjectReference((PetscObject)sub));
715: *B = sub;
716: break;
717: case MAT_REUSE_MATRIX:
718: PetscCheck(sub == *B, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Submatrix was not used before in this call");
719: break;
720: case MAT_IGNORE_MATRIX: /* Nothing to do */
721: break;
722: case MAT_INPLACE_MATRIX: /* Nothing to do */
723: SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MAT_INPLACE_MATRIX is not supported yet");
724: }
725: PetscFunctionReturn(PETSC_SUCCESS);
726: }
728: static PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A, IS isrow, IS iscol, Mat *B)
729: {
730: Mat_Nest *vs = (Mat_Nest *)A->data;
731: Mat sub;
733: PetscFunctionBegin;
734: PetscCall(MatNestFindSubMat(A, &vs->islocal, isrow, iscol, &sub));
735: /* We allow the submatrix to be NULL, perhaps it would be better for the user to return an empty matrix instead */
736: if (sub) PetscCall(PetscObjectReference((PetscObject)sub));
737: *B = sub;
738: PetscFunctionReturn(PETSC_SUCCESS);
739: }
741: static PetscErrorCode MatRestoreLocalSubMatrix_Nest(Mat A, IS isrow, IS iscol, Mat *B)
742: {
743: Mat_Nest *vs = (Mat_Nest *)A->data;
744: Mat sub;
746: PetscFunctionBegin;
747: PetscCall(MatNestFindSubMat(A, &vs->islocal, isrow, iscol, &sub));
748: PetscCheck(*B == sub, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Local submatrix has not been gotten");
749: if (sub) {
750: PetscCheck(((PetscObject)sub)->refct > 1, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Local submatrix has had reference count decremented too many times");
751: PetscCall(MatDestroy(B));
752: }
753: PetscFunctionReturn(PETSC_SUCCESS);
754: }
756: static PetscErrorCode MatGetDiagonal_Nest(Mat A, Vec v)
757: {
758: Mat_Nest *bA = (Mat_Nest *)A->data;
759: PetscInt i;
761: PetscFunctionBegin;
762: for (i = 0; i < bA->nr; i++) {
763: Vec bv;
764: PetscCall(VecGetSubVector(v, bA->isglobal.row[i], &bv));
765: if (bA->m[i][i]) {
766: PetscCall(MatGetDiagonal(bA->m[i][i], bv));
767: } else {
768: PetscCall(VecSet(bv, 0.0));
769: }
770: PetscCall(VecRestoreSubVector(v, bA->isglobal.row[i], &bv));
771: }
772: PetscFunctionReturn(PETSC_SUCCESS);
773: }
775: static PetscErrorCode MatDiagonalScale_Nest(Mat A, Vec l, Vec r)
776: {
777: Mat_Nest *bA = (Mat_Nest *)A->data;
778: Vec bl, *br;
779: PetscInt i, j;
781: PetscFunctionBegin;
782: PetscCall(PetscCalloc1(bA->nc, &br));
783: if (r) {
784: for (j = 0; j < bA->nc; j++) PetscCall(VecGetSubVector(r, bA->isglobal.col[j], &br[j]));
785: }
786: bl = NULL;
787: for (i = 0; i < bA->nr; i++) {
788: if (l) PetscCall(VecGetSubVector(l, bA->isglobal.row[i], &bl));
789: for (j = 0; j < bA->nc; j++) {
790: if (bA->m[i][j]) PetscCall(MatDiagonalScale(bA->m[i][j], bl, br[j]));
791: }
792: if (l) PetscCall(VecRestoreSubVector(l, bA->isglobal.row[i], &bl));
793: }
794: if (r) {
795: for (j = 0; j < bA->nc; j++) PetscCall(VecRestoreSubVector(r, bA->isglobal.col[j], &br[j]));
796: }
797: PetscCall(PetscFree(br));
798: PetscFunctionReturn(PETSC_SUCCESS);
799: }
801: static PetscErrorCode MatScale_Nest(Mat A, PetscScalar a)
802: {
803: Mat_Nest *bA = (Mat_Nest *)A->data;
804: PetscInt i, j;
806: PetscFunctionBegin;
807: for (i = 0; i < bA->nr; i++) {
808: for (j = 0; j < bA->nc; j++) {
809: if (bA->m[i][j]) PetscCall(MatScale(bA->m[i][j], a));
810: }
811: }
812: PetscFunctionReturn(PETSC_SUCCESS);
813: }
815: static PetscErrorCode MatShift_Nest(Mat A, PetscScalar a)
816: {
817: Mat_Nest *bA = (Mat_Nest *)A->data;
818: PetscInt i;
819: PetscBool nnzstate = PETSC_FALSE;
821: PetscFunctionBegin;
822: for (i = 0; i < bA->nr; i++) {
823: PetscObjectState subnnzstate = 0;
824: PetscCheck(bA->m[i][i], PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "No support for shifting an empty diagonal block, insert a matrix in block (%" PetscInt_FMT ",%" PetscInt_FMT ")", i, i);
825: PetscCall(MatShift(bA->m[i][i], a));
826: PetscCall(MatGetNonzeroState(bA->m[i][i], &subnnzstate));
827: nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i * bA->nc + i] != subnnzstate);
828: bA->nnzstate[i * bA->nc + i] = subnnzstate;
829: }
830: if (nnzstate) A->nonzerostate++;
831: PetscFunctionReturn(PETSC_SUCCESS);
832: }
834: static PetscErrorCode MatDiagonalSet_Nest(Mat A, Vec D, InsertMode is)
835: {
836: Mat_Nest *bA = (Mat_Nest *)A->data;
837: PetscInt i;
838: PetscBool nnzstate = PETSC_FALSE;
840: PetscFunctionBegin;
841: for (i = 0; i < bA->nr; i++) {
842: PetscObjectState subnnzstate = 0;
843: Vec bv;
844: PetscCall(VecGetSubVector(D, bA->isglobal.row[i], &bv));
845: if (bA->m[i][i]) {
846: PetscCall(MatDiagonalSet(bA->m[i][i], bv, is));
847: PetscCall(MatGetNonzeroState(bA->m[i][i], &subnnzstate));
848: }
849: PetscCall(VecRestoreSubVector(D, bA->isglobal.row[i], &bv));
850: nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i * bA->nc + i] != subnnzstate);
851: bA->nnzstate[i * bA->nc + i] = subnnzstate;
852: }
853: if (nnzstate) A->nonzerostate++;
854: PetscFunctionReturn(PETSC_SUCCESS);
855: }
857: static PetscErrorCode MatSetRandom_Nest(Mat A, PetscRandom rctx)
858: {
859: Mat_Nest *bA = (Mat_Nest *)A->data;
860: PetscInt i, j;
862: PetscFunctionBegin;
863: for (i = 0; i < bA->nr; i++) {
864: for (j = 0; j < bA->nc; j++) {
865: if (bA->m[i][j]) PetscCall(MatSetRandom(bA->m[i][j], rctx));
866: }
867: }
868: PetscFunctionReturn(PETSC_SUCCESS);
869: }
871: static PetscErrorCode MatCreateVecs_Nest(Mat A, Vec *right, Vec *left)
872: {
873: Mat_Nest *bA = (Mat_Nest *)A->data;
874: Vec *L, *R;
875: MPI_Comm comm;
876: PetscInt i, j;
878: PetscFunctionBegin;
879: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
880: if (right) {
881: /* allocate R */
882: PetscCall(PetscMalloc1(bA->nc, &R));
883: /* Create the right vectors */
884: for (j = 0; j < bA->nc; j++) {
885: for (i = 0; i < bA->nr; i++) {
886: if (bA->m[i][j]) {
887: PetscCall(MatCreateVecs(bA->m[i][j], &R[j], NULL));
888: break;
889: }
890: }
891: PetscCheck(i != bA->nr, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null column.");
892: }
893: PetscCall(VecCreateNest(comm, bA->nc, bA->isglobal.col, R, right));
894: /* hand back control to the nest vector */
895: for (j = 0; j < bA->nc; j++) PetscCall(VecDestroy(&R[j]));
896: PetscCall(PetscFree(R));
897: }
899: if (left) {
900: /* allocate L */
901: PetscCall(PetscMalloc1(bA->nr, &L));
902: /* Create the left vectors */
903: for (i = 0; i < bA->nr; i++) {
904: for (j = 0; j < bA->nc; j++) {
905: if (bA->m[i][j]) {
906: PetscCall(MatCreateVecs(bA->m[i][j], NULL, &L[i]));
907: break;
908: }
909: }
910: PetscCheck(j != bA->nc, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null row.");
911: }
913: PetscCall(VecCreateNest(comm, bA->nr, bA->isglobal.row, L, left));
914: for (i = 0; i < bA->nr; i++) PetscCall(VecDestroy(&L[i]));
916: PetscCall(PetscFree(L));
917: }
918: PetscFunctionReturn(PETSC_SUCCESS);
919: }
921: static PetscErrorCode MatView_Nest(Mat A, PetscViewer viewer)
922: {
923: Mat_Nest *bA = (Mat_Nest *)A->data;
924: PetscBool isascii, viewSub = PETSC_FALSE;
925: PetscInt i, j;
927: PetscFunctionBegin;
928: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
929: if (isascii) {
930: PetscCall(PetscOptionsGetBool(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_view_nest_sub", &viewSub, NULL));
931: PetscCall(PetscViewerASCIIPrintf(viewer, "Matrix object:\n"));
932: PetscCall(PetscViewerASCIIPushTab(viewer));
933: PetscCall(PetscViewerASCIIPrintf(viewer, "type=nest, rows=%" PetscInt_FMT ", cols=%" PetscInt_FMT "\n", bA->nr, bA->nc));
935: PetscCall(PetscViewerASCIIPrintf(viewer, "MatNest structure:\n"));
936: for (i = 0; i < bA->nr; i++) {
937: for (j = 0; j < bA->nc; j++) {
938: MatType type;
939: char name[256] = "", prefix[256] = "";
940: PetscInt NR, NC;
941: PetscBool isNest = PETSC_FALSE;
943: if (!bA->m[i][j]) {
944: PetscCall(PetscViewerASCIIPrintf(viewer, "(%" PetscInt_FMT ",%" PetscInt_FMT ") : NULL\n", i, j));
945: continue;
946: }
947: PetscCall(MatGetSize(bA->m[i][j], &NR, &NC));
948: PetscCall(MatGetType(bA->m[i][j], &type));
949: if (((PetscObject)bA->m[i][j])->name) PetscCall(PetscSNPrintf(name, sizeof(name), "name=\"%s\", ", ((PetscObject)bA->m[i][j])->name));
950: if (((PetscObject)bA->m[i][j])->prefix) PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "prefix=\"%s\", ", ((PetscObject)bA->m[i][j])->prefix));
951: PetscCall(PetscObjectTypeCompare((PetscObject)bA->m[i][j], MATNEST, &isNest));
953: PetscCall(PetscViewerASCIIPrintf(viewer, "(%" PetscInt_FMT ",%" PetscInt_FMT ") : %s%stype=%s, rows=%" PetscInt_FMT ", cols=%" PetscInt_FMT "\n", i, j, name, prefix, type, NR, NC));
955: if (isNest || viewSub) {
956: PetscCall(PetscViewerASCIIPushTab(viewer)); /* push1 */
957: PetscCall(MatView(bA->m[i][j], viewer));
958: PetscCall(PetscViewerASCIIPopTab(viewer)); /* pop1 */
959: }
960: }
961: }
962: PetscCall(PetscViewerASCIIPopTab(viewer)); /* pop0 */
963: }
964: PetscFunctionReturn(PETSC_SUCCESS);
965: }
967: static PetscErrorCode MatZeroEntries_Nest(Mat A)
968: {
969: Mat_Nest *bA = (Mat_Nest *)A->data;
970: PetscInt i, j;
972: PetscFunctionBegin;
973: for (i = 0; i < bA->nr; i++) {
974: for (j = 0; j < bA->nc; j++) {
975: if (!bA->m[i][j]) continue;
976: PetscCall(MatZeroEntries(bA->m[i][j]));
977: }
978: }
979: PetscFunctionReturn(PETSC_SUCCESS);
980: }
982: static PetscErrorCode MatCopy_Nest(Mat A, Mat B, MatStructure str)
983: {
984: Mat_Nest *bA = (Mat_Nest *)A->data, *bB = (Mat_Nest *)B->data;
985: PetscInt i, j, nr = bA->nr, nc = bA->nc;
986: PetscBool nnzstate = PETSC_FALSE;
988: PetscFunctionBegin;
989: PetscCheck(nr == bB->nr && nc == bB->nc, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Cannot copy a Mat_Nest of block size (%" PetscInt_FMT ",%" PetscInt_FMT ") to a Mat_Nest of block size (%" PetscInt_FMT ",%" PetscInt_FMT ")", bB->nr, bB->nc, nr, nc);
990: for (i = 0; i < nr; i++) {
991: for (j = 0; j < nc; j++) {
992: PetscObjectState subnnzstate = 0;
993: if (bA->m[i][j] && bB->m[i][j]) {
994: PetscCall(MatCopy(bA->m[i][j], bB->m[i][j], str));
995: } else PetscCheck(!bA->m[i][j] && !bB->m[i][j], PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Matrix block does not exist at %" PetscInt_FMT ",%" PetscInt_FMT, i, j);
996: PetscCall(MatGetNonzeroState(bB->m[i][j], &subnnzstate));
997: nnzstate = (PetscBool)(nnzstate || bB->nnzstate[i * nc + j] != subnnzstate);
998: bB->nnzstate[i * nc + j] = subnnzstate;
999: }
1000: }
1001: if (nnzstate) B->nonzerostate++;
1002: PetscFunctionReturn(PETSC_SUCCESS);
1003: }
1005: static PetscErrorCode MatAXPY_Nest(Mat Y, PetscScalar a, Mat X, MatStructure str)
1006: {
1007: Mat_Nest *bY = (Mat_Nest *)Y->data, *bX = (Mat_Nest *)X->data;
1008: PetscInt i, j, nr = bY->nr, nc = bY->nc;
1009: PetscBool nnzstate = PETSC_FALSE;
1011: PetscFunctionBegin;
1012: PetscCheck(nr == bX->nr && nc == bX->nc, PetscObjectComm((PetscObject)Y), PETSC_ERR_ARG_INCOMP, "Cannot AXPY a MatNest of block size (%" PetscInt_FMT ",%" PetscInt_FMT ") with a MatNest of block size (%" PetscInt_FMT ",%" PetscInt_FMT ")", bX->nr, bX->nc, nr, nc);
1013: for (i = 0; i < nr; i++) {
1014: for (j = 0; j < nc; j++) {
1015: PetscObjectState subnnzstate = 0;
1016: if (bY->m[i][j] && bX->m[i][j]) {
1017: PetscCall(MatAXPY(bY->m[i][j], a, bX->m[i][j], str));
1018: } else if (bX->m[i][j]) {
1019: Mat M;
1021: PetscCheck(str == DIFFERENT_NONZERO_PATTERN || str == UNKNOWN_NONZERO_PATTERN, PetscObjectComm((PetscObject)Y), PETSC_ERR_ARG_INCOMP, "Matrix block does not exist at %" PetscInt_FMT ",%" PetscInt_FMT ". Use DIFFERENT_NONZERO_PATTERN or UNKNOWN_NONZERO_PATTERN", i, j);
1022: PetscCall(MatDuplicate(bX->m[i][j], MAT_COPY_VALUES, &M));
1023: PetscCall(MatNestSetSubMat(Y, i, j, M));
1024: PetscCall(MatDestroy(&M));
1025: }
1026: if (bY->m[i][j]) PetscCall(MatGetNonzeroState(bY->m[i][j], &subnnzstate));
1027: nnzstate = (PetscBool)(nnzstate || bY->nnzstate[i * nc + j] != subnnzstate);
1028: bY->nnzstate[i * nc + j] = subnnzstate;
1029: }
1030: }
1031: if (nnzstate) Y->nonzerostate++;
1032: PetscFunctionReturn(PETSC_SUCCESS);
1033: }
1035: static PetscErrorCode MatDuplicate_Nest(Mat A, MatDuplicateOption op, Mat *B)
1036: {
1037: Mat_Nest *bA = (Mat_Nest *)A->data;
1038: Mat *b;
1039: PetscInt i, j, nr = bA->nr, nc = bA->nc;
1041: PetscFunctionBegin;
1042: PetscCall(PetscMalloc1(nr * nc, &b));
1043: for (i = 0; i < nr; i++) {
1044: for (j = 0; j < nc; j++) {
1045: if (bA->m[i][j]) {
1046: PetscCall(MatDuplicate(bA->m[i][j], op, &b[i * nc + j]));
1047: } else {
1048: b[i * nc + j] = NULL;
1049: }
1050: }
1051: }
1052: PetscCall(MatCreateNest(PetscObjectComm((PetscObject)A), nr, bA->isglobal.row, nc, bA->isglobal.col, b, B));
1053: /* Give the new MatNest exclusive ownership */
1054: for (i = 0; i < nr * nc; i++) PetscCall(MatDestroy(&b[i]));
1055: PetscCall(PetscFree(b));
1057: PetscCall(MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY));
1058: PetscCall(MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY));
1059: PetscFunctionReturn(PETSC_SUCCESS);
1060: }
1062: /* nest api */
1063: static PetscErrorCode MatNestGetSubMat_Nest(Mat A, PetscInt idxm, PetscInt jdxm, Mat *mat)
1064: {
1065: Mat_Nest *bA = (Mat_Nest *)A->data;
1067: PetscFunctionBegin;
1068: PetscCheck(idxm < bA->nr, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, idxm, bA->nr - 1);
1069: PetscCheck(jdxm < bA->nc, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Col too large: row %" PetscInt_FMT " max %" PetscInt_FMT, jdxm, bA->nc - 1);
1070: *mat = bA->m[idxm][jdxm];
1071: PetscFunctionReturn(PETSC_SUCCESS);
1072: }
1074: /*@
1075: MatNestGetSubMat - Returns a single, sub-matrix from a `MATNEST`
1077: Not Collective
1079: Input Parameters:
1080: + A - `MATNEST` matrix
1081: . idxm - index of the matrix within the nest matrix
1082: - jdxm - index of the matrix within the nest matrix
1084: Output Parameter:
1085: . sub - matrix at index `idxm`, `jdxm` within the nest matrix
1087: Level: developer
1089: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSize()`, `MatNestGetSubMats()`, `MatCreateNest()`, `MatNestSetSubMat()`,
1090: `MatNestGetLocalISs()`, `MatNestGetISs()`
1091: @*/
1092: PetscErrorCode MatNestGetSubMat(Mat A, PetscInt idxm, PetscInt jdxm, Mat *sub)
1093: {
1094: PetscFunctionBegin;
1098: PetscAssertPointer(sub, 4);
1099: PetscUseMethod(A, "MatNestGetSubMat_C", (Mat, PetscInt, PetscInt, Mat *), (A, idxm, jdxm, sub));
1100: PetscFunctionReturn(PETSC_SUCCESS);
1101: }
1103: static PetscErrorCode MatNestSetSubMat_Nest(Mat A, PetscInt idxm, PetscInt jdxm, Mat mat)
1104: {
1105: Mat_Nest *bA = (Mat_Nest *)A->data;
1106: PetscInt m, n, M, N, mi, ni, Mi, Ni;
1108: PetscFunctionBegin;
1109: PetscCheck(idxm < bA->nr, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, idxm, bA->nr - 1);
1110: PetscCheck(jdxm < bA->nc, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Col too large: row %" PetscInt_FMT " max %" PetscInt_FMT, jdxm, bA->nc - 1);
1111: if (mat) {
1112: PetscCall(MatGetLocalSize(mat, &m, &n));
1113: PetscCall(MatGetSize(mat, &M, &N));
1114: PetscCall(ISGetLocalSize(bA->isglobal.row[idxm], &mi));
1115: PetscCall(ISGetSize(bA->isglobal.row[idxm], &Mi));
1116: PetscCall(ISGetLocalSize(bA->isglobal.col[jdxm], &ni));
1117: PetscCall(ISGetSize(bA->isglobal.col[jdxm], &Ni));
1118: PetscCheck(M == Mi && N == Ni, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_INCOMP, "Submatrix dimension (%" PetscInt_FMT ",%" PetscInt_FMT ") incompatible with nest block (%" PetscInt_FMT ",%" PetscInt_FMT ")", M, N, Mi, Ni);
1119: PetscCheck(m == mi && n == ni, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_INCOMP, "Submatrix local dimension (%" PetscInt_FMT ",%" PetscInt_FMT ") incompatible with nest block (%" PetscInt_FMT ",%" PetscInt_FMT ")", m, n, mi, ni);
1120: }
1122: /* do not increase object state */
1123: if (mat == bA->m[idxm][jdxm]) PetscFunctionReturn(PETSC_SUCCESS);
1125: PetscCall(PetscObjectReference((PetscObject)mat));
1126: PetscCall(MatDestroy(&bA->m[idxm][jdxm]));
1127: bA->m[idxm][jdxm] = mat;
1128: PetscCall(PetscObjectStateIncrease((PetscObject)A));
1129: if (mat) PetscCall(MatGetNonzeroState(mat, &bA->nnzstate[idxm * bA->nc + jdxm]));
1130: else bA->nnzstate[idxm * bA->nc + jdxm] = 0;
1131: A->nonzerostate++;
1132: PetscFunctionReturn(PETSC_SUCCESS);
1133: }
1135: /*@
1136: MatNestSetSubMat - Set a single submatrix in the `MATNEST`
1138: Logically Collective
1140: Input Parameters:
1141: + A - `MATNEST` matrix
1142: . idxm - index of the matrix within the nest matrix
1143: . jdxm - index of the matrix within the nest matrix
1144: - sub - matrix at index `idxm`, `jdxm` within the nest matrix
1146: Level: developer
1148: Notes:
1149: The new submatrix must have the same size and communicator as that block of the nest.
1151: This increments the reference count of the submatrix.
1153: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestSetSubMats()`, `MatNestGetSubMats()`, `MatNestGetLocalISs()`, `MatCreateNest()`,
1154: `MatNestGetSubMat()`, `MatNestGetISs()`, `MatNestGetSize()`
1155: @*/
1156: PetscErrorCode MatNestSetSubMat(Mat A, PetscInt idxm, PetscInt jdxm, Mat sub)
1157: {
1158: PetscFunctionBegin;
1163: PetscTryMethod(A, "MatNestSetSubMat_C", (Mat, PetscInt, PetscInt, Mat), (A, idxm, jdxm, sub));
1164: PetscFunctionReturn(PETSC_SUCCESS);
1165: }
1167: static PetscErrorCode MatNestGetSubMats_Nest(Mat A, PetscInt *M, PetscInt *N, Mat ***mat)
1168: {
1169: Mat_Nest *bA = (Mat_Nest *)A->data;
1171: PetscFunctionBegin;
1172: if (M) *M = bA->nr;
1173: if (N) *N = bA->nc;
1174: if (mat) *mat = bA->m;
1175: PetscFunctionReturn(PETSC_SUCCESS);
1176: }
1178: /*@C
1179: MatNestGetSubMats - Returns the entire two dimensional array of matrices defining a `MATNEST` matrix.
1181: Not Collective
1183: Input Parameter:
1184: . A - nest matrix
1186: Output Parameters:
1187: + M - number of rows in the nest matrix
1188: . N - number of cols in the nest matrix
1189: - mat - array of matrices
1191: Level: developer
1193: Note:
1194: The user should not free the array `mat`.
1196: Fortran Notes:
1197: This routine has a calling sequence
1198: $ call MatNestGetSubMats(A, M, N, mat, ierr)
1199: where the space allocated for the optional argument `mat` is assumed large enough (if provided).
1200: Matrices in `mat` are returned in row-major order, see `MatCreateNest()` for an example.
1202: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSize()`, `MatNestGetSubMat()`, `MatNestGetLocalISs()`, `MatCreateNest()`,
1203: `MatNestSetSubMats()`, `MatNestGetISs()`, `MatNestSetSubMat()`
1204: @*/
1205: PetscErrorCode MatNestGetSubMats(Mat A, PetscInt *M, PetscInt *N, Mat ***mat)
1206: {
1207: PetscFunctionBegin;
1209: PetscUseMethod(A, "MatNestGetSubMats_C", (Mat, PetscInt *, PetscInt *, Mat ***), (A, M, N, mat));
1210: PetscFunctionReturn(PETSC_SUCCESS);
1211: }
1213: static PetscErrorCode MatNestGetSize_Nest(Mat A, PetscInt *M, PetscInt *N)
1214: {
1215: Mat_Nest *bA = (Mat_Nest *)A->data;
1217: PetscFunctionBegin;
1218: if (M) *M = bA->nr;
1219: if (N) *N = bA->nc;
1220: PetscFunctionReturn(PETSC_SUCCESS);
1221: }
1223: /*@
1224: MatNestGetSize - Returns the size of the `MATNEST` matrix.
1226: Not Collective
1228: Input Parameter:
1229: . A - `MATNEST` matrix
1231: Output Parameters:
1232: + M - number of rows in the nested mat
1233: - N - number of cols in the nested mat
1235: Level: developer
1237: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSubMat()`, `MatNestGetSubMats()`, `MatCreateNest()`, `MatNestGetLocalISs()`,
1238: `MatNestGetISs()`
1239: @*/
1240: PetscErrorCode MatNestGetSize(Mat A, PetscInt *M, PetscInt *N)
1241: {
1242: PetscFunctionBegin;
1244: PetscUseMethod(A, "MatNestGetSize_C", (Mat, PetscInt *, PetscInt *), (A, M, N));
1245: PetscFunctionReturn(PETSC_SUCCESS);
1246: }
1248: static PetscErrorCode MatNestGetISs_Nest(Mat A, IS rows[], IS cols[])
1249: {
1250: Mat_Nest *vs = (Mat_Nest *)A->data;
1251: PetscInt i;
1253: PetscFunctionBegin;
1254: if (rows)
1255: for (i = 0; i < vs->nr; i++) rows[i] = vs->isglobal.row[i];
1256: if (cols)
1257: for (i = 0; i < vs->nc; i++) cols[i] = vs->isglobal.col[i];
1258: PetscFunctionReturn(PETSC_SUCCESS);
1259: }
1261: /*@C
1262: MatNestGetISs - Returns the index sets partitioning the row and column spaces of a `MATNEST`
1264: Not Collective
1266: Input Parameter:
1267: . A - `MATNEST` matrix
1269: Output Parameters:
1270: + rows - array of row index sets
1271: - cols - array of column index sets
1273: Level: advanced
1275: Note:
1276: The user must have allocated arrays of the correct size. The reference count is not increased on the returned `IS`s.
1278: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSubMat()`, `MatNestGetSubMats()`, `MatNestGetSize()`, `MatNestGetLocalISs()`,
1279: `MatCreateNest()`, `MatNestSetSubMats()`
1280: @*/
1281: PetscErrorCode MatNestGetISs(Mat A, IS rows[], IS cols[])
1282: {
1283: PetscFunctionBegin;
1285: PetscUseMethod(A, "MatNestGetISs_C", (Mat, IS[], IS[]), (A, rows, cols));
1286: PetscFunctionReturn(PETSC_SUCCESS);
1287: }
1289: static PetscErrorCode MatNestGetLocalISs_Nest(Mat A, IS rows[], IS cols[])
1290: {
1291: Mat_Nest *vs = (Mat_Nest *)A->data;
1292: PetscInt i;
1294: PetscFunctionBegin;
1295: if (rows)
1296: for (i = 0; i < vs->nr; i++) rows[i] = vs->islocal.row[i];
1297: if (cols)
1298: for (i = 0; i < vs->nc; i++) cols[i] = vs->islocal.col[i];
1299: PetscFunctionReturn(PETSC_SUCCESS);
1300: }
1302: /*@C
1303: MatNestGetLocalISs - Returns the index sets partitioning the row and column spaces of a `MATNEST`
1305: Not Collective
1307: Input Parameter:
1308: . A - `MATNEST` matrix
1310: Output Parameters:
1311: + rows - array of row index sets (or `NULL` to ignore)
1312: - cols - array of column index sets (or `NULL` to ignore)
1314: Level: advanced
1316: Note:
1317: The user must have allocated arrays of the correct size. The reference count is not increased on the returned `IS`s.
1319: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSubMat()`, `MatNestGetSubMats()`, `MatNestGetSize()`, `MatNestGetISs()`, `MatCreateNest()`,
1320: `MatNestSetSubMats()`, `MatNestSetSubMat()`
1321: @*/
1322: PetscErrorCode MatNestGetLocalISs(Mat A, IS rows[], IS cols[])
1323: {
1324: PetscFunctionBegin;
1326: PetscUseMethod(A, "MatNestGetLocalISs_C", (Mat, IS[], IS[]), (A, rows, cols));
1327: PetscFunctionReturn(PETSC_SUCCESS);
1328: }
1330: static PetscErrorCode MatNestSetVecType_Nest(Mat A, VecType vtype)
1331: {
1332: PetscBool flg;
1334: PetscFunctionBegin;
1335: PetscCall(PetscStrcmp(vtype, VECNEST, &flg));
1336: /* In reality, this only distinguishes VECNEST and "other" */
1337: if (flg) A->ops->getvecs = MatCreateVecs_Nest;
1338: else A->ops->getvecs = (PetscErrorCode(*)(Mat, Vec *, Vec *))0;
1339: PetscFunctionReturn(PETSC_SUCCESS);
1340: }
1342: /*@C
1343: MatNestSetVecType - Sets the type of `Vec` returned by `MatCreateVecs()`
1345: Not Collective
1347: Input Parameters:
1348: + A - `MATNEST` matrix
1349: - vtype - `VecType` to use for creating vectors
1351: Level: developer
1353: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatCreateVecs()`, `MatCreateNest()`, `VecType`
1354: @*/
1355: PetscErrorCode MatNestSetVecType(Mat A, VecType vtype)
1356: {
1357: PetscFunctionBegin;
1359: PetscTryMethod(A, "MatNestSetVecType_C", (Mat, VecType), (A, vtype));
1360: PetscFunctionReturn(PETSC_SUCCESS);
1361: }
1363: static PetscErrorCode MatNestSetSubMats_Nest(Mat A, PetscInt nr, const IS is_row[], PetscInt nc, const IS is_col[], const Mat a[])
1364: {
1365: Mat_Nest *s = (Mat_Nest *)A->data;
1366: PetscInt i, j, m, n, M, N;
1367: PetscBool cong, isstd, sametype = PETSC_FALSE;
1368: VecType vtype, type;
1370: PetscFunctionBegin;
1371: PetscCall(MatReset_Nest(A));
1373: s->nr = nr;
1374: s->nc = nc;
1376: /* Create space for submatrices */
1377: PetscCall(PetscMalloc1(nr, &s->m));
1378: PetscCall(PetscMalloc1(nr * nc, &s->m[0]));
1379: for (i = 0; i < nr; i++) {
1380: s->m[i] = s->m[0] + i * nc;
1381: for (j = 0; j < nc; j++) {
1382: s->m[i][j] = a ? a[i * nc + j] : NULL;
1383: PetscCall(PetscObjectReference((PetscObject)s->m[i][j]));
1384: }
1385: }
1386: PetscCall(MatGetVecType(A, &vtype));
1387: PetscCall(PetscStrcmp(vtype, VECSTANDARD, &isstd));
1388: if (isstd) {
1389: /* check if all blocks have the same vectype */
1390: vtype = NULL;
1391: for (i = 0; i < nr; i++) {
1392: for (j = 0; j < nc; j++) {
1393: if (s->m[i][j]) {
1394: if (!vtype) { /* first visited block */
1395: PetscCall(MatGetVecType(s->m[i][j], &vtype));
1396: sametype = PETSC_TRUE;
1397: } else if (sametype) {
1398: PetscCall(MatGetVecType(s->m[i][j], &type));
1399: PetscCall(PetscStrcmp(vtype, type, &sametype));
1400: }
1401: }
1402: }
1403: }
1404: if (sametype) { /* propagate vectype */
1405: PetscCall(MatSetVecType(A, vtype));
1406: }
1407: }
1409: PetscCall(MatSetUp_NestIS_Private(A, nr, is_row, nc, is_col));
1411: PetscCall(PetscMalloc1(nr, &s->row_len));
1412: PetscCall(PetscMalloc1(nc, &s->col_len));
1413: for (i = 0; i < nr; i++) s->row_len[i] = -1;
1414: for (j = 0; j < nc; j++) s->col_len[j] = -1;
1416: PetscCall(PetscCalloc1(nr * nc, &s->nnzstate));
1417: for (i = 0; i < nr; i++) {
1418: for (j = 0; j < nc; j++) {
1419: if (s->m[i][j]) PetscCall(MatGetNonzeroState(s->m[i][j], &s->nnzstate[i * nc + j]));
1420: }
1421: }
1423: PetscCall(MatNestGetSizes_Private(A, &m, &n, &M, &N));
1425: PetscCall(PetscLayoutSetSize(A->rmap, M));
1426: PetscCall(PetscLayoutSetLocalSize(A->rmap, m));
1427: PetscCall(PetscLayoutSetSize(A->cmap, N));
1428: PetscCall(PetscLayoutSetLocalSize(A->cmap, n));
1430: PetscCall(PetscLayoutSetUp(A->rmap));
1431: PetscCall(PetscLayoutSetUp(A->cmap));
1433: /* disable operations that are not supported for non-square matrices,
1434: or matrices for which is_row != is_col */
1435: PetscCall(MatHasCongruentLayouts(A, &cong));
1436: if (cong && nr != nc) cong = PETSC_FALSE;
1437: if (cong) {
1438: for (i = 0; cong && i < nr; i++) PetscCall(ISEqualUnsorted(s->isglobal.row[i], s->isglobal.col[i], &cong));
1439: }
1440: if (!cong) {
1441: A->ops->missingdiagonal = NULL;
1442: A->ops->getdiagonal = NULL;
1443: A->ops->shift = NULL;
1444: A->ops->diagonalset = NULL;
1445: }
1447: PetscCall(PetscCalloc2(nr, &s->left, nc, &s->right));
1448: PetscCall(PetscObjectStateIncrease((PetscObject)A));
1449: A->nonzerostate++;
1450: PetscFunctionReturn(PETSC_SUCCESS);
1451: }
1453: /*@
1454: MatNestSetSubMats - Sets the nested submatrices in a `MATNEST`
1456: Collective
1458: Input Parameters:
1459: + A - `MATNEST` matrix
1460: . nr - number of nested row blocks
1461: . is_row - index sets for each nested row block, or `NULL` to make contiguous
1462: . nc - number of nested column blocks
1463: . is_col - index sets for each nested column block, or `NULL` to make contiguous
1464: - a - array of nr*nc submatrices, or `NULL`
1466: Level: advanced
1468: Notes:
1469: This always resets any block matrix information previously set.
1470: Pass `NULL` in the corresponding entry of `a` for an empty block.
1472: In both C and Fortran, `a` must be a row-major order array containing the matrices. See
1473: `MatCreateNest()` for an example.
1475: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatCreateNest()`, `MatNestSetSubMat()`, `MatNestGetSubMat()`, `MatNestGetSubMats()`
1476: @*/
1477: PetscErrorCode MatNestSetSubMats(Mat A, PetscInt nr, const IS is_row[], PetscInt nc, const IS is_col[], const Mat a[])
1478: {
1479: PetscFunctionBegin;
1482: PetscCheck(nr >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Number of rows cannot be negative");
1483: if (nr && is_row) {
1484: PetscAssertPointer(is_row, 3);
1486: }
1488: PetscCheck(nc >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Number of columns cannot be negative");
1489: if (nc && is_col) {
1490: PetscAssertPointer(is_col, 5);
1492: }
1493: PetscTryMethod(A, "MatNestSetSubMats_C", (Mat, PetscInt, const IS[], PetscInt, const IS[], const Mat[]), (A, nr, is_row, nc, is_col, a));
1494: PetscFunctionReturn(PETSC_SUCCESS);
1495: }
1497: static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A, PetscInt n, const IS islocal[], const IS isglobal[], PetscBool colflg, ISLocalToGlobalMapping *ltog)
1498: {
1499: PetscBool flg;
1500: PetscInt i, j, m, mi, *ix;
1502: PetscFunctionBegin;
1503: *ltog = NULL;
1504: for (i = 0, m = 0, flg = PETSC_FALSE; i < n; i++) {
1505: if (islocal[i]) {
1506: PetscCall(ISGetLocalSize(islocal[i], &mi));
1507: flg = PETSC_TRUE; /* We found a non-trivial entry */
1508: } else {
1509: PetscCall(ISGetLocalSize(isglobal[i], &mi));
1510: }
1511: m += mi;
1512: }
1513: if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
1515: PetscCall(PetscMalloc1(m, &ix));
1516: for (i = 0, m = 0; i < n; i++) {
1517: ISLocalToGlobalMapping smap = NULL;
1518: Mat sub = NULL;
1519: PetscSF sf;
1520: PetscLayout map;
1521: const PetscInt *ix2;
1523: if (!colflg) {
1524: PetscCall(MatNestFindNonzeroSubMatRow(A, i, &sub));
1525: } else {
1526: PetscCall(MatNestFindNonzeroSubMatCol(A, i, &sub));
1527: }
1528: if (sub) {
1529: if (!colflg) {
1530: PetscCall(MatGetLocalToGlobalMapping(sub, &smap, NULL));
1531: } else {
1532: PetscCall(MatGetLocalToGlobalMapping(sub, NULL, &smap));
1533: }
1534: }
1535: /*
1536: Now we need to extract the monolithic global indices that correspond to the given split global indices.
1537: In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces.
1538: */
1539: PetscCall(ISGetIndices(isglobal[i], &ix2));
1540: if (islocal[i]) {
1541: PetscInt *ilocal, *iremote;
1542: PetscInt mil, nleaves;
1544: PetscCall(ISGetLocalSize(islocal[i], &mi));
1545: PetscCheck(smap, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing local to global map");
1546: for (j = 0; j < mi; j++) ix[m + j] = j;
1547: PetscCall(ISLocalToGlobalMappingApply(smap, mi, ix + m, ix + m));
1549: /* PetscSFSetGraphLayout does not like negative indices */
1550: PetscCall(PetscMalloc2(mi, &ilocal, mi, &iremote));
1551: for (j = 0, nleaves = 0; j < mi; j++) {
1552: if (ix[m + j] < 0) continue;
1553: ilocal[nleaves] = j;
1554: iremote[nleaves] = ix[m + j];
1555: nleaves++;
1556: }
1557: PetscCall(ISGetLocalSize(isglobal[i], &mil));
1558: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)A), &sf));
1559: PetscCall(PetscLayoutCreate(PetscObjectComm((PetscObject)A), &map));
1560: PetscCall(PetscLayoutSetLocalSize(map, mil));
1561: PetscCall(PetscLayoutSetUp(map));
1562: PetscCall(PetscSFSetGraphLayout(sf, map, nleaves, ilocal, PETSC_USE_POINTER, iremote));
1563: PetscCall(PetscLayoutDestroy(&map));
1564: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, ix2, ix + m, MPI_REPLACE));
1565: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, ix2, ix + m, MPI_REPLACE));
1566: PetscCall(PetscSFDestroy(&sf));
1567: PetscCall(PetscFree2(ilocal, iremote));
1568: } else {
1569: PetscCall(ISGetLocalSize(isglobal[i], &mi));
1570: for (j = 0; j < mi; j++) ix[m + j] = ix2[i];
1571: }
1572: PetscCall(ISRestoreIndices(isglobal[i], &ix2));
1573: m += mi;
1574: }
1575: PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A), 1, m, ix, PETSC_OWN_POINTER, ltog));
1576: PetscFunctionReturn(PETSC_SUCCESS);
1577: }
1579: /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */
1580: /*
1581: nprocessors = NP
1582: Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1))
1583: proc 0: => (g_0,h_0,)
1584: proc 1: => (g_1,h_1,)
1585: ...
1586: proc nprocs-1: => (g_NP-1,h_NP-1,)
1588: proc 0: proc 1: proc nprocs-1:
1589: is[0] = (0,1,2,...,nlocal(g_0)-1) (0,1,...,nlocal(g_1)-1) (0,1,...,nlocal(g_NP-1))
1591: proc 0:
1592: is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1)
1593: proc 1:
1594: is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1)
1596: proc NP-1:
1597: is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1)
1598: */
1599: static PetscErrorCode MatSetUp_NestIS_Private(Mat A, PetscInt nr, const IS is_row[], PetscInt nc, const IS is_col[])
1600: {
1601: Mat_Nest *vs = (Mat_Nest *)A->data;
1602: PetscInt i, j, offset, n, nsum, bs;
1603: Mat sub = NULL;
1605: PetscFunctionBegin;
1606: PetscCall(PetscMalloc1(nr, &vs->isglobal.row));
1607: PetscCall(PetscMalloc1(nc, &vs->isglobal.col));
1608: if (is_row) { /* valid IS is passed in */
1609: /* refs on is[] are incremented */
1610: for (i = 0; i < vs->nr; i++) {
1611: PetscCall(PetscObjectReference((PetscObject)is_row[i]));
1613: vs->isglobal.row[i] = is_row[i];
1614: }
1615: } else { /* Create the ISs by inspecting sizes of a submatrix in each row */
1616: nsum = 0;
1617: for (i = 0; i < vs->nr; i++) { /* Add up the local sizes to compute the aggregate offset */
1618: PetscCall(MatNestFindNonzeroSubMatRow(A, i, &sub));
1619: PetscCheck(sub, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "No nonzero submatrix in row %" PetscInt_FMT, i);
1620: PetscCall(MatGetLocalSize(sub, &n, NULL));
1621: PetscCheck(n >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Sizes have not yet been set for submatrix");
1622: nsum += n;
1623: }
1624: PetscCallMPI(MPI_Scan(&nsum, &offset, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)A)));
1625: offset -= nsum;
1626: for (i = 0; i < vs->nr; i++) {
1627: PetscCall(MatNestFindNonzeroSubMatRow(A, i, &sub));
1628: PetscCall(MatGetLocalSize(sub, &n, NULL));
1629: PetscCall(MatGetBlockSizes(sub, &bs, NULL));
1630: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)sub), n, offset, 1, &vs->isglobal.row[i]));
1631: PetscCall(ISSetBlockSize(vs->isglobal.row[i], bs));
1632: offset += n;
1633: }
1634: }
1636: if (is_col) { /* valid IS is passed in */
1637: /* refs on is[] are incremented */
1638: for (j = 0; j < vs->nc; j++) {
1639: PetscCall(PetscObjectReference((PetscObject)is_col[j]));
1641: vs->isglobal.col[j] = is_col[j];
1642: }
1643: } else { /* Create the ISs by inspecting sizes of a submatrix in each column */
1644: offset = A->cmap->rstart;
1645: nsum = 0;
1646: for (j = 0; j < vs->nc; j++) {
1647: PetscCall(MatNestFindNonzeroSubMatCol(A, j, &sub));
1648: PetscCheck(sub, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "No nonzero submatrix in column %" PetscInt_FMT, i);
1649: PetscCall(MatGetLocalSize(sub, NULL, &n));
1650: PetscCheck(n >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Sizes have not yet been set for submatrix");
1651: nsum += n;
1652: }
1653: PetscCallMPI(MPI_Scan(&nsum, &offset, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)A)));
1654: offset -= nsum;
1655: for (j = 0; j < vs->nc; j++) {
1656: PetscCall(MatNestFindNonzeroSubMatCol(A, j, &sub));
1657: PetscCall(MatGetLocalSize(sub, NULL, &n));
1658: PetscCall(MatGetBlockSizes(sub, NULL, &bs));
1659: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)sub), n, offset, 1, &vs->isglobal.col[j]));
1660: PetscCall(ISSetBlockSize(vs->isglobal.col[j], bs));
1661: offset += n;
1662: }
1663: }
1665: /* Set up the local ISs */
1666: PetscCall(PetscMalloc1(vs->nr, &vs->islocal.row));
1667: PetscCall(PetscMalloc1(vs->nc, &vs->islocal.col));
1668: for (i = 0, offset = 0; i < vs->nr; i++) {
1669: IS isloc;
1670: ISLocalToGlobalMapping rmap = NULL;
1671: PetscInt nlocal, bs;
1672: PetscCall(MatNestFindNonzeroSubMatRow(A, i, &sub));
1673: if (sub) PetscCall(MatGetLocalToGlobalMapping(sub, &rmap, NULL));
1674: if (rmap) {
1675: PetscCall(MatGetBlockSizes(sub, &bs, NULL));
1676: PetscCall(ISLocalToGlobalMappingGetSize(rmap, &nlocal));
1677: PetscCall(ISCreateStride(PETSC_COMM_SELF, nlocal, offset, 1, &isloc));
1678: PetscCall(ISSetBlockSize(isloc, bs));
1679: } else {
1680: nlocal = 0;
1681: isloc = NULL;
1682: }
1683: vs->islocal.row[i] = isloc;
1684: offset += nlocal;
1685: }
1686: for (i = 0, offset = 0; i < vs->nc; i++) {
1687: IS isloc;
1688: ISLocalToGlobalMapping cmap = NULL;
1689: PetscInt nlocal, bs;
1690: PetscCall(MatNestFindNonzeroSubMatCol(A, i, &sub));
1691: if (sub) PetscCall(MatGetLocalToGlobalMapping(sub, NULL, &cmap));
1692: if (cmap) {
1693: PetscCall(MatGetBlockSizes(sub, NULL, &bs));
1694: PetscCall(ISLocalToGlobalMappingGetSize(cmap, &nlocal));
1695: PetscCall(ISCreateStride(PETSC_COMM_SELF, nlocal, offset, 1, &isloc));
1696: PetscCall(ISSetBlockSize(isloc, bs));
1697: } else {
1698: nlocal = 0;
1699: isloc = NULL;
1700: }
1701: vs->islocal.col[i] = isloc;
1702: offset += nlocal;
1703: }
1705: /* Set up the aggregate ISLocalToGlobalMapping */
1706: {
1707: ISLocalToGlobalMapping rmap, cmap;
1708: PetscCall(MatNestCreateAggregateL2G_Private(A, vs->nr, vs->islocal.row, vs->isglobal.row, PETSC_FALSE, &rmap));
1709: PetscCall(MatNestCreateAggregateL2G_Private(A, vs->nc, vs->islocal.col, vs->isglobal.col, PETSC_TRUE, &cmap));
1710: if (rmap && cmap) PetscCall(MatSetLocalToGlobalMapping(A, rmap, cmap));
1711: PetscCall(ISLocalToGlobalMappingDestroy(&rmap));
1712: PetscCall(ISLocalToGlobalMappingDestroy(&cmap));
1713: }
1715: if (PetscDefined(USE_DEBUG)) {
1716: for (i = 0; i < vs->nr; i++) {
1717: for (j = 0; j < vs->nc; j++) {
1718: PetscInt m, n, M, N, mi, ni, Mi, Ni;
1719: Mat B = vs->m[i][j];
1720: if (!B) continue;
1721: PetscCall(MatGetSize(B, &M, &N));
1722: PetscCall(MatGetLocalSize(B, &m, &n));
1723: PetscCall(ISGetSize(vs->isglobal.row[i], &Mi));
1724: PetscCall(ISGetSize(vs->isglobal.col[j], &Ni));
1725: PetscCall(ISGetLocalSize(vs->isglobal.row[i], &mi));
1726: PetscCall(ISGetLocalSize(vs->isglobal.col[j], &ni));
1727: PetscCheck(M == Mi && N == Ni, PetscObjectComm((PetscObject)sub), PETSC_ERR_ARG_INCOMP, "Global sizes (%" PetscInt_FMT ",%" PetscInt_FMT ") of nested submatrix (%" PetscInt_FMT ",%" PetscInt_FMT ") do not agree with space defined by index sets (%" PetscInt_FMT ",%" PetscInt_FMT ")", M, N, i, j, Mi, Ni);
1728: PetscCheck(m == mi && n == ni, PetscObjectComm((PetscObject)sub), PETSC_ERR_ARG_INCOMP, "Local sizes (%" PetscInt_FMT ",%" PetscInt_FMT ") of nested submatrix (%" PetscInt_FMT ",%" PetscInt_FMT ") do not agree with space defined by index sets (%" PetscInt_FMT ",%" PetscInt_FMT ")", m, n, i, j, mi, ni);
1729: }
1730: }
1731: }
1733: /* Set A->assembled if all non-null blocks are currently assembled */
1734: for (i = 0; i < vs->nr; i++) {
1735: for (j = 0; j < vs->nc; j++) {
1736: if (vs->m[i][j] && !vs->m[i][j]->assembled) PetscFunctionReturn(PETSC_SUCCESS);
1737: }
1738: }
1739: A->assembled = PETSC_TRUE;
1740: PetscFunctionReturn(PETSC_SUCCESS);
1741: }
1743: /*@C
1744: MatCreateNest - Creates a new `MATNEST` matrix containing several nested submatrices, each stored separately
1746: Collective
1748: Input Parameters:
1749: + comm - Communicator for the new `MATNEST`
1750: . nr - number of nested row blocks
1751: . is_row - index sets for each nested row block, or `NULL` to make contiguous
1752: . nc - number of nested column blocks
1753: . is_col - index sets for each nested column block, or `NULL` to make contiguous
1754: - a - array of nr*nc submatrices, empty submatrices can be passed using `NULL`
1756: Output Parameter:
1757: . B - new matrix
1759: Note:
1760: In both C and Fortran, `a` must be a row-major order array holding references to the matrices.
1761: For instance, to represent the matrix
1762: $\begin{bmatrix} A_{11} & A_{12} \\ A_{21} & A_{22}\end{bmatrix}$
1763: one should use `Mat a[4]={A11,A12,A21,A22}`.
1765: Level: advanced
1767: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatCreate()`, `VecCreateNest()`, `DMCreateMatrix()`, `MatNestSetSubMat()`,
1768: `MatNestGetSubMat()`, `MatNestGetLocalISs()`, `MatNestGetSize()`,
1769: `MatNestGetISs()`, `MatNestSetSubMats()`, `MatNestGetSubMats()`
1770: @*/
1771: PetscErrorCode MatCreateNest(MPI_Comm comm, PetscInt nr, const IS is_row[], PetscInt nc, const IS is_col[], const Mat a[], Mat *B)
1772: {
1773: PetscFunctionBegin;
1774: PetscCall(MatCreate(comm, B));
1775: PetscCall(MatSetType(*B, MATNEST));
1776: (*B)->preallocated = PETSC_TRUE;
1777: PetscCall(MatNestSetSubMats(*B, nr, is_row, nc, is_col, a));
1778: PetscFunctionReturn(PETSC_SUCCESS);
1779: }
1781: static PetscErrorCode MatConvert_Nest_SeqAIJ_fast(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
1782: {
1783: Mat_Nest *nest = (Mat_Nest *)A->data;
1784: Mat *trans;
1785: PetscScalar **avv;
1786: PetscScalar *vv;
1787: PetscInt **aii, **ajj;
1788: PetscInt *ii, *jj, *ci;
1789: PetscInt nr, nc, nnz, i, j;
1790: PetscBool done;
1792: PetscFunctionBegin;
1793: PetscCall(MatGetSize(A, &nr, &nc));
1794: if (reuse == MAT_REUSE_MATRIX) {
1795: PetscInt rnr;
1797: PetscCall(MatGetRowIJ(*newmat, 0, PETSC_FALSE, PETSC_FALSE, &rnr, (const PetscInt **)&ii, (const PetscInt **)&jj, &done));
1798: PetscCheck(done, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "MatGetRowIJ");
1799: PetscCheck(rnr == nr, PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "Cannot reuse matrix, wrong number of rows");
1800: PetscCall(MatSeqAIJGetArray(*newmat, &vv));
1801: }
1802: /* extract CSR for nested SeqAIJ matrices */
1803: nnz = 0;
1804: PetscCall(PetscCalloc4(nest->nr * nest->nc, &aii, nest->nr * nest->nc, &ajj, nest->nr * nest->nc, &avv, nest->nr * nest->nc, &trans));
1805: for (i = 0; i < nest->nr; ++i) {
1806: for (j = 0; j < nest->nc; ++j) {
1807: Mat B = nest->m[i][j];
1808: if (B) {
1809: PetscScalar *naa;
1810: PetscInt *nii, *njj, nnr;
1811: PetscBool istrans;
1813: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &istrans));
1814: if (istrans) {
1815: Mat Bt;
1817: PetscCall(MatTransposeGetMat(B, &Bt));
1818: PetscCall(MatTranspose(Bt, MAT_INITIAL_MATRIX, &trans[i * nest->nc + j]));
1819: B = trans[i * nest->nc + j];
1820: } else {
1821: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATHERMITIANTRANSPOSEVIRTUAL, &istrans));
1822: if (istrans) {
1823: Mat Bt;
1825: PetscCall(MatHermitianTransposeGetMat(B, &Bt));
1826: PetscCall(MatHermitianTranspose(Bt, MAT_INITIAL_MATRIX, &trans[i * nest->nc + j]));
1827: B = trans[i * nest->nc + j];
1828: }
1829: }
1830: PetscCall(MatGetRowIJ(B, 0, PETSC_FALSE, PETSC_FALSE, &nnr, (const PetscInt **)&nii, (const PetscInt **)&njj, &done));
1831: PetscCheck(done, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "MatGetRowIJ");
1832: PetscCall(MatSeqAIJGetArray(B, &naa));
1833: nnz += nii[nnr];
1835: aii[i * nest->nc + j] = nii;
1836: ajj[i * nest->nc + j] = njj;
1837: avv[i * nest->nc + j] = naa;
1838: }
1839: }
1840: }
1841: if (reuse != MAT_REUSE_MATRIX) {
1842: PetscCall(PetscMalloc1(nr + 1, &ii));
1843: PetscCall(PetscMalloc1(nnz, &jj));
1844: PetscCall(PetscMalloc1(nnz, &vv));
1845: } else {
1846: PetscCheck(nnz == ii[nr], PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "Cannot reuse matrix, wrong number of nonzeros");
1847: }
1849: /* new row pointer */
1850: PetscCall(PetscArrayzero(ii, nr + 1));
1851: for (i = 0; i < nest->nr; ++i) {
1852: PetscInt ncr, rst;
1854: PetscCall(ISStrideGetInfo(nest->isglobal.row[i], &rst, NULL));
1855: PetscCall(ISGetLocalSize(nest->isglobal.row[i], &ncr));
1856: for (j = 0; j < nest->nc; ++j) {
1857: if (aii[i * nest->nc + j]) {
1858: PetscInt *nii = aii[i * nest->nc + j];
1859: PetscInt ir;
1861: for (ir = rst; ir < ncr + rst; ++ir) {
1862: ii[ir + 1] += nii[1] - nii[0];
1863: nii++;
1864: }
1865: }
1866: }
1867: }
1868: for (i = 0; i < nr; i++) ii[i + 1] += ii[i];
1870: /* construct CSR for the new matrix */
1871: PetscCall(PetscCalloc1(nr, &ci));
1872: for (i = 0; i < nest->nr; ++i) {
1873: PetscInt ncr, rst;
1875: PetscCall(ISStrideGetInfo(nest->isglobal.row[i], &rst, NULL));
1876: PetscCall(ISGetLocalSize(nest->isglobal.row[i], &ncr));
1877: for (j = 0; j < nest->nc; ++j) {
1878: if (aii[i * nest->nc + j]) {
1879: PetscScalar *nvv = avv[i * nest->nc + j];
1880: PetscInt *nii = aii[i * nest->nc + j];
1881: PetscInt *njj = ajj[i * nest->nc + j];
1882: PetscInt ir, cst;
1884: PetscCall(ISStrideGetInfo(nest->isglobal.col[j], &cst, NULL));
1885: for (ir = rst; ir < ncr + rst; ++ir) {
1886: PetscInt ij, rsize = nii[1] - nii[0], ist = ii[ir] + ci[ir];
1888: for (ij = 0; ij < rsize; ij++) {
1889: jj[ist + ij] = *njj + cst;
1890: vv[ist + ij] = *nvv;
1891: njj++;
1892: nvv++;
1893: }
1894: ci[ir] += rsize;
1895: nii++;
1896: }
1897: }
1898: }
1899: }
1900: PetscCall(PetscFree(ci));
1902: /* restore info */
1903: for (i = 0; i < nest->nr; ++i) {
1904: for (j = 0; j < nest->nc; ++j) {
1905: Mat B = nest->m[i][j];
1906: if (B) {
1907: PetscInt nnr = 0, k = i * nest->nc + j;
1909: B = (trans[k] ? trans[k] : B);
1910: PetscCall(MatRestoreRowIJ(B, 0, PETSC_FALSE, PETSC_FALSE, &nnr, (const PetscInt **)&aii[k], (const PetscInt **)&ajj[k], &done));
1911: PetscCheck(done, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "MatRestoreRowIJ");
1912: PetscCall(MatSeqAIJRestoreArray(B, &avv[k]));
1913: PetscCall(MatDestroy(&trans[k]));
1914: }
1915: }
1916: }
1917: PetscCall(PetscFree4(aii, ajj, avv, trans));
1919: /* finalize newmat */
1920: if (reuse == MAT_INITIAL_MATRIX) {
1921: PetscCall(MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A), nr, nc, ii, jj, vv, newmat));
1922: } else if (reuse == MAT_INPLACE_MATRIX) {
1923: Mat B;
1925: PetscCall(MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A), nr, nc, ii, jj, vv, &B));
1926: PetscCall(MatHeaderReplace(A, &B));
1927: }
1928: PetscCall(MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY));
1929: PetscCall(MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY));
1930: {
1931: Mat_SeqAIJ *a = (Mat_SeqAIJ *)((*newmat)->data);
1932: a->free_a = PETSC_TRUE;
1933: a->free_ij = PETSC_TRUE;
1934: }
1935: PetscFunctionReturn(PETSC_SUCCESS);
1936: }
1938: PETSC_INTERN PetscErrorCode MatAXPY_Dense_Nest(Mat Y, PetscScalar a, Mat X)
1939: {
1940: Mat_Nest *nest = (Mat_Nest *)X->data;
1941: PetscInt i, j, k, rstart;
1942: PetscBool flg;
1944: PetscFunctionBegin;
1945: /* Fill by row */
1946: for (j = 0; j < nest->nc; ++j) {
1947: /* Using global column indices and ISAllGather() is not scalable. */
1948: IS bNis;
1949: PetscInt bN;
1950: const PetscInt *bNindices;
1951: PetscCall(ISAllGather(nest->isglobal.col[j], &bNis));
1952: PetscCall(ISGetSize(bNis, &bN));
1953: PetscCall(ISGetIndices(bNis, &bNindices));
1954: for (i = 0; i < nest->nr; ++i) {
1955: Mat B = nest->m[i][j], D = NULL;
1956: PetscInt bm, br;
1957: const PetscInt *bmindices;
1958: if (!B) continue;
1959: PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &flg, MATTRANSPOSEVIRTUAL, MATHERMITIANTRANSPOSEVIRTUAL, ""));
1960: if (flg) {
1961: PetscTryMethod(B, "MatTransposeGetMat_C", (Mat, Mat *), (B, &D));
1962: PetscTryMethod(B, "MatHermitianTransposeGetMat_C", (Mat, Mat *), (B, &D));
1963: PetscCall(MatConvert(B, ((PetscObject)D)->type_name, MAT_INITIAL_MATRIX, &D));
1964: B = D;
1965: }
1966: PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &flg, MATSEQSBAIJ, MATMPISBAIJ, ""));
1967: if (flg) {
1968: if (D) PetscCall(MatConvert(D, MATBAIJ, MAT_INPLACE_MATRIX, &D));
1969: else PetscCall(MatConvert(B, MATBAIJ, MAT_INITIAL_MATRIX, &D));
1970: B = D;
1971: }
1972: PetscCall(ISGetLocalSize(nest->isglobal.row[i], &bm));
1973: PetscCall(ISGetIndices(nest->isglobal.row[i], &bmindices));
1974: PetscCall(MatGetOwnershipRange(B, &rstart, NULL));
1975: for (br = 0; br < bm; ++br) {
1976: PetscInt row = bmindices[br], brncols, *cols;
1977: const PetscInt *brcols;
1978: const PetscScalar *brcoldata;
1979: PetscScalar *vals = NULL;
1980: PetscCall(MatGetRow(B, br + rstart, &brncols, &brcols, &brcoldata));
1981: PetscCall(PetscMalloc1(brncols, &cols));
1982: for (k = 0; k < brncols; k++) cols[k] = bNindices[brcols[k]];
1983: /*
1984: Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match.
1985: Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES.
1986: */
1987: if (a != 1.0) {
1988: PetscCall(PetscMalloc1(brncols, &vals));
1989: for (k = 0; k < brncols; k++) vals[k] = a * brcoldata[k];
1990: PetscCall(MatSetValues(Y, 1, &row, brncols, cols, vals, ADD_VALUES));
1991: PetscCall(PetscFree(vals));
1992: } else {
1993: PetscCall(MatSetValues(Y, 1, &row, brncols, cols, brcoldata, ADD_VALUES));
1994: }
1995: PetscCall(MatRestoreRow(B, br + rstart, &brncols, &brcols, &brcoldata));
1996: PetscCall(PetscFree(cols));
1997: }
1998: if (D) PetscCall(MatDestroy(&D));
1999: PetscCall(ISRestoreIndices(nest->isglobal.row[i], &bmindices));
2000: }
2001: PetscCall(ISRestoreIndices(bNis, &bNindices));
2002: PetscCall(ISDestroy(&bNis));
2003: }
2004: PetscCall(MatAssemblyBegin(Y, MAT_FINAL_ASSEMBLY));
2005: PetscCall(MatAssemblyEnd(Y, MAT_FINAL_ASSEMBLY));
2006: PetscFunctionReturn(PETSC_SUCCESS);
2007: }
2009: static PetscErrorCode MatConvert_Nest_AIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
2010: {
2011: Mat_Nest *nest = (Mat_Nest *)A->data;
2012: PetscInt m, n, M, N, i, j, k, *dnnz, *onnz = NULL, rstart, cstart, cend;
2013: PetscMPIInt size;
2014: Mat C;
2016: PetscFunctionBegin;
2017: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
2018: if (size == 1) { /* look for a special case with SeqAIJ matrices and strided-1, contiguous, blocks */
2019: PetscInt nf;
2020: PetscBool fast;
2022: PetscCall(PetscStrcmp(newtype, MATAIJ, &fast));
2023: if (!fast) PetscCall(PetscStrcmp(newtype, MATSEQAIJ, &fast));
2024: for (i = 0; i < nest->nr && fast; ++i) {
2025: for (j = 0; j < nest->nc && fast; ++j) {
2026: Mat B = nest->m[i][j];
2027: if (B) {
2028: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATSEQAIJ, &fast));
2029: if (!fast) {
2030: PetscBool istrans;
2032: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &istrans));
2033: if (istrans) {
2034: Mat Bt;
2036: PetscCall(MatTransposeGetMat(B, &Bt));
2037: PetscCall(PetscObjectTypeCompare((PetscObject)Bt, MATSEQAIJ, &fast));
2038: } else {
2039: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATHERMITIANTRANSPOSEVIRTUAL, &istrans));
2040: if (istrans) {
2041: Mat Bt;
2043: PetscCall(MatHermitianTransposeGetMat(B, &Bt));
2044: PetscCall(PetscObjectTypeCompare((PetscObject)Bt, MATSEQAIJ, &fast));
2045: }
2046: }
2047: }
2048: }
2049: }
2050: }
2051: for (i = 0, nf = 0; i < nest->nr && fast; ++i) {
2052: PetscCall(PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i], ISSTRIDE, &fast));
2053: if (fast) {
2054: PetscInt f, s;
2056: PetscCall(ISStrideGetInfo(nest->isglobal.row[i], &f, &s));
2057: if (f != nf || s != 1) {
2058: fast = PETSC_FALSE;
2059: } else {
2060: PetscCall(ISGetSize(nest->isglobal.row[i], &f));
2061: nf += f;
2062: }
2063: }
2064: }
2065: for (i = 0, nf = 0; i < nest->nc && fast; ++i) {
2066: PetscCall(PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i], ISSTRIDE, &fast));
2067: if (fast) {
2068: PetscInt f, s;
2070: PetscCall(ISStrideGetInfo(nest->isglobal.col[i], &f, &s));
2071: if (f != nf || s != 1) {
2072: fast = PETSC_FALSE;
2073: } else {
2074: PetscCall(ISGetSize(nest->isglobal.col[i], &f));
2075: nf += f;
2076: }
2077: }
2078: }
2079: if (fast) {
2080: PetscCall(MatConvert_Nest_SeqAIJ_fast(A, newtype, reuse, newmat));
2081: PetscFunctionReturn(PETSC_SUCCESS);
2082: }
2083: }
2084: PetscCall(MatGetSize(A, &M, &N));
2085: PetscCall(MatGetLocalSize(A, &m, &n));
2086: PetscCall(MatGetOwnershipRangeColumn(A, &cstart, &cend));
2087: if (reuse == MAT_REUSE_MATRIX) C = *newmat;
2088: else {
2089: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
2090: PetscCall(MatSetType(C, newtype));
2091: PetscCall(MatSetSizes(C, m, n, M, N));
2092: }
2093: PetscCall(PetscMalloc1(2 * m, &dnnz));
2094: if (m) {
2095: onnz = dnnz + m;
2096: for (k = 0; k < m; k++) {
2097: dnnz[k] = 0;
2098: onnz[k] = 0;
2099: }
2100: }
2101: for (j = 0; j < nest->nc; ++j) {
2102: IS bNis;
2103: PetscInt bN;
2104: const PetscInt *bNindices;
2105: PetscBool flg;
2106: /* Using global column indices and ISAllGather() is not scalable. */
2107: PetscCall(ISAllGather(nest->isglobal.col[j], &bNis));
2108: PetscCall(ISGetSize(bNis, &bN));
2109: PetscCall(ISGetIndices(bNis, &bNindices));
2110: for (i = 0; i < nest->nr; ++i) {
2111: PetscSF bmsf;
2112: PetscSFNode *iremote;
2113: Mat B = nest->m[i][j], D = NULL;
2114: PetscInt bm, *sub_dnnz, *sub_onnz, br;
2115: const PetscInt *bmindices;
2116: if (!B) continue;
2117: PetscCall(ISGetLocalSize(nest->isglobal.row[i], &bm));
2118: PetscCall(ISGetIndices(nest->isglobal.row[i], &bmindices));
2119: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf));
2120: PetscCall(PetscMalloc1(bm, &iremote));
2121: PetscCall(PetscMalloc1(bm, &sub_dnnz));
2122: PetscCall(PetscMalloc1(bm, &sub_onnz));
2123: for (k = 0; k < bm; ++k) {
2124: sub_dnnz[k] = 0;
2125: sub_onnz[k] = 0;
2126: }
2127: PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &flg, MATTRANSPOSEVIRTUAL, MATHERMITIANTRANSPOSEVIRTUAL, ""));
2128: if (flg) {
2129: PetscTryMethod(B, "MatTransposeGetMat_C", (Mat, Mat *), (B, &D));
2130: PetscTryMethod(B, "MatHermitianTransposeGetMat_C", (Mat, Mat *), (B, &D));
2131: PetscCall(MatConvert(B, ((PetscObject)D)->type_name, MAT_INITIAL_MATRIX, &D));
2132: B = D;
2133: }
2134: PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &flg, MATSEQSBAIJ, MATMPISBAIJ, ""));
2135: if (flg) {
2136: if (D) PetscCall(MatConvert(D, MATBAIJ, MAT_INPLACE_MATRIX, &D));
2137: else PetscCall(MatConvert(B, MATBAIJ, MAT_INITIAL_MATRIX, &D));
2138: B = D;
2139: }
2140: /*
2141: Locate the owners for all of the locally-owned global row indices for this row block.
2142: These determine the roots of PetscSF used to communicate preallocation data to row owners.
2143: The roots correspond to the dnnz and onnz entries; thus, there are two roots per row.
2144: */
2145: PetscCall(MatGetOwnershipRange(B, &rstart, NULL));
2146: for (br = 0; br < bm; ++br) {
2147: PetscInt row = bmindices[br], brncols, col;
2148: const PetscInt *brcols;
2149: PetscInt rowrel = 0; /* row's relative index on its owner rank */
2150: PetscMPIInt rowowner = 0;
2151: PetscCall(PetscLayoutFindOwnerIndex(A->rmap, row, &rowowner, &rowrel));
2152: /* how many roots */
2153: iremote[br].rank = rowowner;
2154: iremote[br].index = rowrel; /* edge from bmdnnz to dnnz */
2155: /* get nonzero pattern */
2156: PetscCall(MatGetRow(B, br + rstart, &brncols, &brcols, NULL));
2157: for (k = 0; k < brncols; k++) {
2158: col = bNindices[brcols[k]];
2159: if (col >= A->cmap->range[rowowner] && col < A->cmap->range[rowowner + 1]) {
2160: sub_dnnz[br]++;
2161: } else {
2162: sub_onnz[br]++;
2163: }
2164: }
2165: PetscCall(MatRestoreRow(B, br + rstart, &brncols, &brcols, NULL));
2166: }
2167: if (D) PetscCall(MatDestroy(&D));
2168: PetscCall(ISRestoreIndices(nest->isglobal.row[i], &bmindices));
2169: /* bsf will have to take care of disposing of bedges. */
2170: PetscCall(PetscSFSetGraph(bmsf, m, bm, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
2171: PetscCall(PetscSFReduceBegin(bmsf, MPIU_INT, sub_dnnz, dnnz, MPI_SUM));
2172: PetscCall(PetscSFReduceEnd(bmsf, MPIU_INT, sub_dnnz, dnnz, MPI_SUM));
2173: PetscCall(PetscSFReduceBegin(bmsf, MPIU_INT, sub_onnz, onnz, MPI_SUM));
2174: PetscCall(PetscSFReduceEnd(bmsf, MPIU_INT, sub_onnz, onnz, MPI_SUM));
2175: PetscCall(PetscFree(sub_dnnz));
2176: PetscCall(PetscFree(sub_onnz));
2177: PetscCall(PetscSFDestroy(&bmsf));
2178: }
2179: PetscCall(ISRestoreIndices(bNis, &bNindices));
2180: PetscCall(ISDestroy(&bNis));
2181: }
2182: /* Resize preallocation if overestimated */
2183: for (i = 0; i < m; i++) {
2184: dnnz[i] = PetscMin(dnnz[i], A->cmap->n);
2185: onnz[i] = PetscMin(onnz[i], A->cmap->N - A->cmap->n);
2186: }
2187: PetscCall(MatSeqAIJSetPreallocation(C, 0, dnnz));
2188: PetscCall(MatMPIAIJSetPreallocation(C, 0, dnnz, 0, onnz));
2189: PetscCall(PetscFree(dnnz));
2190: PetscCall(MatAXPY_Dense_Nest(C, 1.0, A));
2191: if (reuse == MAT_INPLACE_MATRIX) {
2192: PetscCall(MatHeaderReplace(A, &C));
2193: } else *newmat = C;
2194: PetscFunctionReturn(PETSC_SUCCESS);
2195: }
2197: static PetscErrorCode MatConvert_Nest_Dense(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
2198: {
2199: Mat B;
2200: PetscInt m, n, M, N;
2202: PetscFunctionBegin;
2203: PetscCall(MatGetSize(A, &M, &N));
2204: PetscCall(MatGetLocalSize(A, &m, &n));
2205: if (reuse == MAT_REUSE_MATRIX) {
2206: B = *newmat;
2207: PetscCall(MatZeroEntries(B));
2208: } else {
2209: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A), m, PETSC_DECIDE, M, N, NULL, &B));
2210: }
2211: PetscCall(MatAXPY_Dense_Nest(B, 1.0, A));
2212: if (reuse == MAT_INPLACE_MATRIX) {
2213: PetscCall(MatHeaderReplace(A, &B));
2214: } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
2215: PetscFunctionReturn(PETSC_SUCCESS);
2216: }
2218: static PetscErrorCode MatHasOperation_Nest(Mat mat, MatOperation op, PetscBool *has)
2219: {
2220: Mat_Nest *bA = (Mat_Nest *)mat->data;
2221: MatOperation opAdd;
2222: PetscInt i, j, nr = bA->nr, nc = bA->nc;
2223: PetscBool flg;
2225: PetscFunctionBegin;
2226: *has = PETSC_FALSE;
2227: if (op == MATOP_MULT || op == MATOP_MULT_ADD || op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) {
2228: opAdd = (op == MATOP_MULT || op == MATOP_MULT_ADD ? MATOP_MULT_ADD : MATOP_MULT_TRANSPOSE_ADD);
2229: for (j = 0; j < nc; j++) {
2230: for (i = 0; i < nr; i++) {
2231: if (!bA->m[i][j]) continue;
2232: PetscCall(MatHasOperation(bA->m[i][j], opAdd, &flg));
2233: if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
2234: }
2235: }
2236: }
2237: if (((void **)mat->ops)[op]) *has = PETSC_TRUE;
2238: PetscFunctionReturn(PETSC_SUCCESS);
2239: }
2241: /*MC
2242: MATNEST - "nest" - Matrix type consisting of nested submatrices, each stored separately.
2244: Level: intermediate
2246: Notes:
2247: This matrix type permits scalable use of `PCFIELDSPLIT` and avoids the large memory costs of extracting submatrices.
2248: It allows the use of symmetric and block formats for parts of multi-physics simulations.
2249: It is usually used with `DMCOMPOSITE` and `DMCreateMatrix()`
2251: Each of the submatrices lives on the same MPI communicator as the original nest matrix (though they can have zero
2252: rows/columns on some processes.) Thus this is not meant for cases where the submatrices live on far fewer processes
2253: than the nest matrix.
2255: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatCreate()`, `MatType`, `MatCreateNest()`, `MatNestSetSubMat()`, `MatNestGetSubMat()`,
2256: `VecCreateNest()`, `DMCreateMatrix()`, `DMCOMPOSITE`, `MatNestSetVecType()`, `MatNestGetLocalISs()`,
2257: `MatNestGetISs()`, `MatNestSetSubMats()`, `MatNestGetSubMats()`
2258: M*/
2259: PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A)
2260: {
2261: Mat_Nest *s;
2263: PetscFunctionBegin;
2264: PetscCall(PetscNew(&s));
2265: A->data = (void *)s;
2267: s->nr = -1;
2268: s->nc = -1;
2269: s->m = NULL;
2270: s->splitassembly = PETSC_FALSE;
2272: PetscCall(PetscMemzero(A->ops, sizeof(*A->ops)));
2274: A->ops->mult = MatMult_Nest;
2275: A->ops->multadd = MatMultAdd_Nest;
2276: A->ops->multtranspose = MatMultTranspose_Nest;
2277: A->ops->multtransposeadd = MatMultTransposeAdd_Nest;
2278: A->ops->transpose = MatTranspose_Nest;
2279: A->ops->multhermitiantranspose = MatMultHermitianTranspose_Nest;
2280: A->ops->multhermitiantransposeadd = MatMultHermitianTransposeAdd_Nest;
2281: A->ops->assemblybegin = MatAssemblyBegin_Nest;
2282: A->ops->assemblyend = MatAssemblyEnd_Nest;
2283: A->ops->zeroentries = MatZeroEntries_Nest;
2284: A->ops->copy = MatCopy_Nest;
2285: A->ops->axpy = MatAXPY_Nest;
2286: A->ops->duplicate = MatDuplicate_Nest;
2287: A->ops->createsubmatrix = MatCreateSubMatrix_Nest;
2288: A->ops->destroy = MatDestroy_Nest;
2289: A->ops->view = MatView_Nest;
2290: A->ops->getvecs = NULL; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */
2291: A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_Nest;
2292: A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest;
2293: A->ops->getdiagonal = MatGetDiagonal_Nest;
2294: A->ops->diagonalscale = MatDiagonalScale_Nest;
2295: A->ops->scale = MatScale_Nest;
2296: A->ops->shift = MatShift_Nest;
2297: A->ops->diagonalset = MatDiagonalSet_Nest;
2298: A->ops->setrandom = MatSetRandom_Nest;
2299: A->ops->hasoperation = MatHasOperation_Nest;
2300: A->ops->missingdiagonal = MatMissingDiagonal_Nest;
2302: A->spptr = NULL;
2303: A->assembled = PETSC_FALSE;
2305: /* expose Nest api's */
2306: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSubMat_C", MatNestGetSubMat_Nest));
2307: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetSubMat_C", MatNestSetSubMat_Nest));
2308: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSubMats_C", MatNestGetSubMats_Nest));
2309: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSize_C", MatNestGetSize_Nest));
2310: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetISs_C", MatNestGetISs_Nest));
2311: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetLocalISs_C", MatNestGetLocalISs_Nest));
2312: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetVecType_C", MatNestSetVecType_Nest));
2313: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetSubMats_C", MatNestSetSubMats_Nest));
2314: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_mpiaij_C", MatConvert_Nest_AIJ));
2315: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_seqaij_C", MatConvert_Nest_AIJ));
2316: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_aij_C", MatConvert_Nest_AIJ));
2317: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_is_C", MatConvert_Nest_IS));
2318: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_mpidense_C", MatConvert_Nest_Dense));
2319: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_seqdense_C", MatConvert_Nest_Dense));
2320: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_seqdense_C", MatProductSetFromOptions_Nest_Dense));
2321: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_mpidense_C", MatProductSetFromOptions_Nest_Dense));
2323: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATNEST));
2324: PetscFunctionReturn(PETSC_SUCCESS);
2325: }