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