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