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