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