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_AB(Mat C)
251: {
252: PetscFunctionBegin;
253: C->ops->productsymbolic = MatProductSymbolic_Nest_Dense;
254: PetscFunctionReturn(PETSC_SUCCESS);
255: }
257: static PetscErrorCode MatProductSetFromOptions_Nest_Dense(Mat C)
258: {
259: Mat_Product *product = C->product;
261: PetscFunctionBegin;
262: if (product->type == MATPRODUCT_AB) PetscCall(MatProductSetFromOptions_Nest_Dense_AB(C));
263: PetscFunctionReturn(PETSC_SUCCESS);
264: }
266: static PetscErrorCode MatMultTranspose_Nest(Mat A, Vec x, Vec y)
267: {
268: Mat_Nest *bA = (Mat_Nest *)A->data;
269: Vec *bx = bA->left, *by = bA->right;
270: PetscInt i, j, nr = bA->nr, nc = bA->nc;
272: PetscFunctionBegin;
273: for (i = 0; i < nr; i++) PetscCall(VecGetSubVector(x, bA->isglobal.row[i], &bx[i]));
274: for (i = 0; i < nc; i++) PetscCall(VecGetSubVector(y, bA->isglobal.col[i], &by[i]));
275: for (j = 0; j < nc; j++) {
276: PetscCall(VecZeroEntries(by[j]));
277: for (i = 0; i < nr; i++) {
278: if (!bA->m[i][j]) continue;
279: /* y[j] <- y[j] + (A[i][j])^T * x[i] */
280: PetscCall(MatMultTransposeAdd(bA->m[i][j], bx[i], by[j], by[j]));
281: }
282: }
283: for (i = 0; i < nr; i++) PetscCall(VecRestoreSubVector(x, bA->isglobal.row[i], &bx[i]));
284: for (i = 0; i < nc; i++) PetscCall(VecRestoreSubVector(y, bA->isglobal.col[i], &by[i]));
285: PetscFunctionReturn(PETSC_SUCCESS);
286: }
288: static PetscErrorCode MatMultTransposeAdd_Nest(Mat A, Vec x, Vec y, Vec z)
289: {
290: Mat_Nest *bA = (Mat_Nest *)A->data;
291: Vec *bx = bA->left, *bz = bA->right;
292: PetscInt i, j, nr = bA->nr, nc = bA->nc;
294: PetscFunctionBegin;
295: for (i = 0; i < nr; i++) PetscCall(VecGetSubVector(x, bA->isglobal.row[i], &bx[i]));
296: for (i = 0; i < nc; i++) PetscCall(VecGetSubVector(z, bA->isglobal.col[i], &bz[i]));
297: for (j = 0; j < nc; j++) {
298: if (y != z) {
299: Vec by;
300: PetscCall(VecGetSubVector(y, bA->isglobal.col[j], &by));
301: PetscCall(VecCopy(by, bz[j]));
302: PetscCall(VecRestoreSubVector(y, bA->isglobal.col[j], &by));
303: }
304: for (i = 0; i < nr; i++) {
305: if (!bA->m[i][j]) continue;
306: /* z[j] <- y[j] + (A[i][j])^T * x[i] */
307: PetscCall(MatMultTransposeAdd(bA->m[i][j], bx[i], bz[j], bz[j]));
308: }
309: }
310: for (i = 0; i < nr; i++) PetscCall(VecRestoreSubVector(x, bA->isglobal.row[i], &bx[i]));
311: for (i = 0; i < nc; i++) PetscCall(VecRestoreSubVector(z, bA->isglobal.col[i], &bz[i]));
312: PetscFunctionReturn(PETSC_SUCCESS);
313: }
315: static PetscErrorCode MatTranspose_Nest(Mat A, MatReuse reuse, Mat *B)
316: {
317: Mat_Nest *bA = (Mat_Nest *)A->data, *bC;
318: Mat C;
319: PetscInt i, j, nr = bA->nr, nc = bA->nc;
321: PetscFunctionBegin;
322: if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *B));
323: PetscCheck(reuse != MAT_INPLACE_MATRIX || nr == nc, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_SIZ, "Square nested matrix only for in-place");
325: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
326: Mat *subs;
327: IS *is_row, *is_col;
329: PetscCall(PetscCalloc1(nr * nc, &subs));
330: PetscCall(PetscMalloc2(nr, &is_row, nc, &is_col));
331: PetscCall(MatNestGetISs(A, is_row, is_col));
332: if (reuse == MAT_INPLACE_MATRIX) {
333: for (i = 0; i < nr; i++) {
334: for (j = 0; j < nc; j++) subs[i + nr * j] = bA->m[i][j];
335: }
336: }
338: PetscCall(MatCreateNest(PetscObjectComm((PetscObject)A), nc, is_col, nr, is_row, subs, &C));
339: PetscCall(PetscFree(subs));
340: PetscCall(PetscFree2(is_row, is_col));
341: } else {
342: C = *B;
343: }
345: bC = (Mat_Nest *)C->data;
346: for (i = 0; i < nr; i++) {
347: for (j = 0; j < nc; j++) {
348: if (bA->m[i][j]) {
349: PetscCall(MatTranspose(bA->m[i][j], reuse, &(bC->m[j][i])));
350: } else {
351: bC->m[j][i] = NULL;
352: }
353: }
354: }
356: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
357: *B = C;
358: } else {
359: PetscCall(MatHeaderMerge(A, &C));
360: }
361: PetscFunctionReturn(PETSC_SUCCESS);
362: }
364: static PetscErrorCode MatNestDestroyISList(PetscInt n, IS **list)
365: {
366: IS *lst = *list;
367: PetscInt i;
369: PetscFunctionBegin;
370: if (!lst) PetscFunctionReturn(PETSC_SUCCESS);
371: for (i = 0; i < n; i++)
372: if (lst[i]) PetscCall(ISDestroy(&lst[i]));
373: PetscCall(PetscFree(lst));
374: *list = NULL;
375: PetscFunctionReturn(PETSC_SUCCESS);
376: }
378: static PetscErrorCode MatReset_Nest(Mat A)
379: {
380: Mat_Nest *vs = (Mat_Nest *)A->data;
381: PetscInt i, j;
383: PetscFunctionBegin;
384: /* release the matrices and the place holders */
385: PetscCall(MatNestDestroyISList(vs->nr, &vs->isglobal.row));
386: PetscCall(MatNestDestroyISList(vs->nc, &vs->isglobal.col));
387: PetscCall(MatNestDestroyISList(vs->nr, &vs->islocal.row));
388: PetscCall(MatNestDestroyISList(vs->nc, &vs->islocal.col));
390: PetscCall(PetscFree(vs->row_len));
391: PetscCall(PetscFree(vs->col_len));
392: PetscCall(PetscFree(vs->nnzstate));
394: PetscCall(PetscFree2(vs->left, vs->right));
396: /* release the matrices and the place holders */
397: if (vs->m) {
398: for (i = 0; i < vs->nr; i++) {
399: for (j = 0; j < vs->nc; j++) PetscCall(MatDestroy(&vs->m[i][j]));
400: PetscCall(PetscFree(vs->m[i]));
401: }
402: PetscCall(PetscFree(vs->m));
403: }
405: /* restore defaults */
406: vs->nr = 0;
407: vs->nc = 0;
408: vs->splitassembly = PETSC_FALSE;
409: PetscFunctionReturn(PETSC_SUCCESS);
410: }
412: static PetscErrorCode MatDestroy_Nest(Mat A)
413: {
414: PetscFunctionBegin;
415: PetscCall(MatReset_Nest(A));
416: PetscCall(PetscFree(A->data));
417: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSubMat_C", NULL));
418: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetSubMat_C", NULL));
419: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSubMats_C", NULL));
420: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSize_C", NULL));
421: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetISs_C", NULL));
422: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetLocalISs_C", NULL));
423: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetVecType_C", NULL));
424: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetSubMats_C", NULL));
425: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_mpiaij_C", NULL));
426: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_seqaij_C", NULL));
427: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_aij_C", NULL));
428: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_is_C", NULL));
429: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_mpidense_C", NULL));
430: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_seqdense_C", NULL));
431: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_seqdense_C", NULL));
432: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_mpidense_C", NULL));
433: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_dense_C", NULL));
434: PetscFunctionReturn(PETSC_SUCCESS);
435: }
437: static PetscErrorCode MatMissingDiagonal_Nest(Mat mat, PetscBool *missing, PetscInt *dd)
438: {
439: Mat_Nest *vs = (Mat_Nest *)mat->data;
440: PetscInt i;
442: PetscFunctionBegin;
443: if (dd) *dd = 0;
444: if (!vs->nr) {
445: *missing = PETSC_TRUE;
446: PetscFunctionReturn(PETSC_SUCCESS);
447: }
448: *missing = PETSC_FALSE;
449: for (i = 0; i < vs->nr && !(*missing); i++) {
450: *missing = PETSC_TRUE;
451: if (vs->m[i][i]) {
452: PetscCall(MatMissingDiagonal(vs->m[i][i], missing, NULL));
453: PetscCheck(!*missing || !dd, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "First missing entry not yet implemented");
454: }
455: }
456: PetscFunctionReturn(PETSC_SUCCESS);
457: }
459: static PetscErrorCode MatAssemblyBegin_Nest(Mat A, MatAssemblyType type)
460: {
461: Mat_Nest *vs = (Mat_Nest *)A->data;
462: PetscInt i, j;
463: PetscBool nnzstate = PETSC_FALSE;
465: PetscFunctionBegin;
466: for (i = 0; i < vs->nr; i++) {
467: for (j = 0; j < vs->nc; j++) {
468: PetscObjectState subnnzstate = 0;
469: if (vs->m[i][j]) {
470: PetscCall(MatAssemblyBegin(vs->m[i][j], type));
471: if (!vs->splitassembly) {
472: /* Note: split assembly will fail if the same block appears more than once (even indirectly through a nested
473: * sub-block). This could be fixed by adding a flag to Mat so that there was a way to check if a Mat was
474: * already performing an assembly, but the result would by more complicated and appears to offer less
475: * potential for diagnostics and correctness checking. Split assembly should be fixed once there is an
476: * interface for libraries to make asynchronous progress in "user-defined non-blocking collectives".
477: */
478: PetscCall(MatAssemblyEnd(vs->m[i][j], type));
479: PetscCall(MatGetNonzeroState(vs->m[i][j], &subnnzstate));
480: }
481: }
482: nnzstate = (PetscBool)(nnzstate || vs->nnzstate[i * vs->nc + j] != subnnzstate);
483: vs->nnzstate[i * vs->nc + j] = subnnzstate;
484: }
485: }
486: if (nnzstate) A->nonzerostate++;
487: PetscFunctionReturn(PETSC_SUCCESS);
488: }
490: static PetscErrorCode MatAssemblyEnd_Nest(Mat A, MatAssemblyType type)
491: {
492: Mat_Nest *vs = (Mat_Nest *)A->data;
493: PetscInt i, j;
495: PetscFunctionBegin;
496: for (i = 0; i < vs->nr; i++) {
497: for (j = 0; j < vs->nc; j++) {
498: if (vs->m[i][j]) {
499: if (vs->splitassembly) PetscCall(MatAssemblyEnd(vs->m[i][j], type));
500: }
501: }
502: }
503: PetscFunctionReturn(PETSC_SUCCESS);
504: }
506: static PetscErrorCode MatNestFindNonzeroSubMatRow(Mat A, PetscInt row, Mat *B)
507: {
508: Mat_Nest *vs = (Mat_Nest *)A->data;
509: PetscInt j;
510: Mat sub;
512: PetscFunctionBegin;
513: sub = (row < vs->nc) ? vs->m[row][row] : (Mat)NULL; /* Prefer to find on the diagonal */
514: for (j = 0; !sub && j < vs->nc; j++) sub = vs->m[row][j];
515: if (sub) PetscCall(MatSetUp(sub)); /* Ensure that the sizes are available */
516: *B = sub;
517: PetscFunctionReturn(PETSC_SUCCESS);
518: }
520: static PetscErrorCode MatNestFindNonzeroSubMatCol(Mat A, PetscInt col, Mat *B)
521: {
522: Mat_Nest *vs = (Mat_Nest *)A->data;
523: PetscInt i;
524: Mat sub;
526: PetscFunctionBegin;
527: sub = (col < vs->nr) ? vs->m[col][col] : (Mat)NULL; /* Prefer to find on the diagonal */
528: for (i = 0; !sub && i < vs->nr; i++) sub = vs->m[i][col];
529: if (sub) PetscCall(MatSetUp(sub)); /* Ensure that the sizes are available */
530: *B = sub;
531: PetscFunctionReturn(PETSC_SUCCESS);
532: }
534: static PetscErrorCode MatNestFindISRange(Mat A, PetscInt n, const IS list[], IS is, PetscInt *begin, PetscInt *end)
535: {
536: PetscInt i, j, size, m;
537: PetscBool flg;
538: IS out, concatenate[2];
540: PetscFunctionBegin;
541: PetscAssertPointer(list, 3);
543: if (begin) {
544: PetscAssertPointer(begin, 5);
545: *begin = -1;
546: }
547: if (end) {
548: PetscAssertPointer(end, 6);
549: *end = -1;
550: }
551: for (i = 0; i < n; i++) {
552: if (!list[i]) continue;
553: PetscCall(ISEqualUnsorted(list[i], is, &flg));
554: if (flg) {
555: if (begin) *begin = i;
556: if (end) *end = i + 1;
557: PetscFunctionReturn(PETSC_SUCCESS);
558: }
559: }
560: PetscCall(ISGetSize(is, &size));
561: for (i = 0; i < n - 1; i++) {
562: if (!list[i]) continue;
563: m = 0;
564: PetscCall(ISConcatenate(PetscObjectComm((PetscObject)A), 2, list + i, &out));
565: PetscCall(ISGetSize(out, &m));
566: for (j = i + 2; j < n && m < size; j++) {
567: if (list[j]) {
568: concatenate[0] = out;
569: concatenate[1] = list[j];
570: PetscCall(ISConcatenate(PetscObjectComm((PetscObject)A), 2, concatenate, &out));
571: PetscCall(ISDestroy(concatenate));
572: PetscCall(ISGetSize(out, &m));
573: }
574: }
575: if (m == size) {
576: PetscCall(ISEqualUnsorted(out, is, &flg));
577: if (flg) {
578: if (begin) *begin = i;
579: if (end) *end = j;
580: PetscCall(ISDestroy(&out));
581: PetscFunctionReturn(PETSC_SUCCESS);
582: }
583: }
584: PetscCall(ISDestroy(&out));
585: }
586: PetscFunctionReturn(PETSC_SUCCESS);
587: }
589: static PetscErrorCode MatNestFillEmptyMat_Private(Mat A, PetscInt i, PetscInt j, Mat *B)
590: {
591: Mat_Nest *vs = (Mat_Nest *)A->data;
592: PetscInt lr, lc;
594: PetscFunctionBegin;
595: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B));
596: PetscCall(ISGetLocalSize(vs->isglobal.row[i], &lr));
597: PetscCall(ISGetLocalSize(vs->isglobal.col[j], &lc));
598: PetscCall(MatSetSizes(*B, lr, lc, PETSC_DECIDE, PETSC_DECIDE));
599: PetscCall(MatSetType(*B, MATAIJ));
600: PetscCall(MatSeqAIJSetPreallocation(*B, 0, NULL));
601: PetscCall(MatMPIAIJSetPreallocation(*B, 0, NULL, 0, NULL));
602: PetscCall(MatSetUp(*B));
603: PetscCall(MatSetOption(*B, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
604: PetscCall(MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY));
605: PetscCall(MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY));
606: PetscFunctionReturn(PETSC_SUCCESS);
607: }
609: static PetscErrorCode MatNestGetBlock_Private(Mat A, PetscInt rbegin, PetscInt rend, PetscInt cbegin, PetscInt cend, Mat *B)
610: {
611: Mat_Nest *vs = (Mat_Nest *)A->data;
612: Mat *a;
613: PetscInt i, j, k, l, nr = rend - rbegin, nc = cend - cbegin;
614: char keyname[256];
615: PetscBool *b;
616: PetscBool flg;
618: PetscFunctionBegin;
619: *B = NULL;
620: PetscCall(PetscSNPrintf(keyname, sizeof(keyname), "NestBlock_%" PetscInt_FMT "-%" PetscInt_FMT "x%" PetscInt_FMT "-%" PetscInt_FMT, rbegin, rend, cbegin, cend));
621: PetscCall(PetscObjectQuery((PetscObject)A, keyname, (PetscObject *)B));
622: if (*B) PetscFunctionReturn(PETSC_SUCCESS);
624: PetscCall(PetscMalloc2(nr * nc, &a, nr * nc, &b));
625: for (i = 0; i < nr; i++) {
626: for (j = 0; j < nc; j++) {
627: a[i * nc + j] = vs->m[rbegin + i][cbegin + j];
628: b[i * nc + j] = PETSC_FALSE;
629: }
630: }
631: if (nc != vs->nc && nr != vs->nr) {
632: for (i = 0; i < nr; i++) {
633: for (j = 0; j < nc; j++) {
634: flg = PETSC_FALSE;
635: for (k = 0; (k < nr && !flg); k++) {
636: if (a[j + k * nc]) flg = PETSC_TRUE;
637: }
638: if (flg) {
639: flg = PETSC_FALSE;
640: for (l = 0; (l < nc && !flg); l++) {
641: if (a[i * nc + l]) flg = PETSC_TRUE;
642: }
643: }
644: if (!flg) {
645: b[i * nc + j] = PETSC_TRUE;
646: PetscCall(MatNestFillEmptyMat_Private(A, rbegin + i, cbegin + j, a + i * nc + j));
647: }
648: }
649: }
650: }
651: PetscCall(MatCreateNest(PetscObjectComm((PetscObject)A), nr, nr != vs->nr ? NULL : vs->isglobal.row, nc, nc != vs->nc ? NULL : vs->isglobal.col, a, B));
652: for (i = 0; i < nr; i++) {
653: for (j = 0; j < nc; j++) {
654: if (b[i * nc + j]) PetscCall(MatDestroy(a + i * nc + j));
655: }
656: }
657: PetscCall(PetscFree2(a, b));
658: (*B)->assembled = A->assembled;
659: PetscCall(PetscObjectCompose((PetscObject)A, keyname, (PetscObject)*B));
660: PetscCall(PetscObjectDereference((PetscObject)*B)); /* Leave the only remaining reference in the composition */
661: PetscFunctionReturn(PETSC_SUCCESS);
662: }
664: static PetscErrorCode MatNestFindSubMat(Mat A, struct MatNestISPair *is, IS isrow, IS iscol, Mat *B)
665: {
666: Mat_Nest *vs = (Mat_Nest *)A->data;
667: PetscInt rbegin, rend, cbegin, cend;
669: PetscFunctionBegin;
670: PetscCall(MatNestFindISRange(A, vs->nr, is->row, isrow, &rbegin, &rend));
671: PetscCall(MatNestFindISRange(A, vs->nc, is->col, iscol, &cbegin, &cend));
672: if (rend == rbegin + 1 && cend == cbegin + 1) {
673: if (!vs->m[rbegin][cbegin]) PetscCall(MatNestFillEmptyMat_Private(A, rbegin, cbegin, vs->m[rbegin] + cbegin));
674: *B = vs->m[rbegin][cbegin];
675: } else if (rbegin != -1 && cbegin != -1) {
676: PetscCall(MatNestGetBlock_Private(A, rbegin, rend, cbegin, cend, B));
677: } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Could not find index set");
678: PetscFunctionReturn(PETSC_SUCCESS);
679: }
681: /*
682: TODO: This does not actually returns a submatrix we can modify
683: */
684: static PetscErrorCode MatCreateSubMatrix_Nest(Mat A, IS isrow, IS iscol, MatReuse reuse, Mat *B)
685: {
686: Mat_Nest *vs = (Mat_Nest *)A->data;
687: Mat sub;
689: PetscFunctionBegin;
690: PetscCall(MatNestFindSubMat(A, &vs->isglobal, isrow, iscol, &sub));
691: switch (reuse) {
692: case MAT_INITIAL_MATRIX:
693: if (sub) PetscCall(PetscObjectReference((PetscObject)sub));
694: *B = sub;
695: break;
696: case MAT_REUSE_MATRIX:
697: PetscCheck(sub == *B, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Submatrix was not used before in this call");
698: break;
699: case MAT_IGNORE_MATRIX: /* Nothing to do */
700: break;
701: case MAT_INPLACE_MATRIX: /* Nothing to do */
702: SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MAT_INPLACE_MATRIX is not supported yet");
703: }
704: PetscFunctionReturn(PETSC_SUCCESS);
705: }
707: static PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A, IS isrow, IS iscol, Mat *B)
708: {
709: Mat_Nest *vs = (Mat_Nest *)A->data;
710: Mat sub;
712: PetscFunctionBegin;
713: PetscCall(MatNestFindSubMat(A, &vs->islocal, isrow, iscol, &sub));
714: /* We allow the submatrix to be NULL, perhaps it would be better for the user to return an empty matrix instead */
715: if (sub) PetscCall(PetscObjectReference((PetscObject)sub));
716: *B = sub;
717: PetscFunctionReturn(PETSC_SUCCESS);
718: }
720: static PetscErrorCode MatRestoreLocalSubMatrix_Nest(Mat A, IS isrow, IS iscol, Mat *B)
721: {
722: Mat_Nest *vs = (Mat_Nest *)A->data;
723: Mat sub;
725: PetscFunctionBegin;
726: PetscCall(MatNestFindSubMat(A, &vs->islocal, isrow, iscol, &sub));
727: PetscCheck(*B == sub, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Local submatrix has not been gotten");
728: if (sub) {
729: PetscCheck(((PetscObject)sub)->refct > 1, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Local submatrix has had reference count decremented too many times");
730: PetscCall(MatDestroy(B));
731: }
732: PetscFunctionReturn(PETSC_SUCCESS);
733: }
735: static PetscErrorCode MatGetDiagonal_Nest(Mat A, Vec v)
736: {
737: Mat_Nest *bA = (Mat_Nest *)A->data;
738: PetscInt i;
740: PetscFunctionBegin;
741: for (i = 0; i < bA->nr; i++) {
742: Vec bv;
743: PetscCall(VecGetSubVector(v, bA->isglobal.row[i], &bv));
744: if (bA->m[i][i]) {
745: PetscCall(MatGetDiagonal(bA->m[i][i], bv));
746: } else {
747: PetscCall(VecSet(bv, 0.0));
748: }
749: PetscCall(VecRestoreSubVector(v, bA->isglobal.row[i], &bv));
750: }
751: PetscFunctionReturn(PETSC_SUCCESS);
752: }
754: static PetscErrorCode MatDiagonalScale_Nest(Mat A, Vec l, Vec r)
755: {
756: Mat_Nest *bA = (Mat_Nest *)A->data;
757: Vec bl, *br;
758: PetscInt i, j;
760: PetscFunctionBegin;
761: PetscCall(PetscCalloc1(bA->nc, &br));
762: if (r) {
763: for (j = 0; j < bA->nc; j++) PetscCall(VecGetSubVector(r, bA->isglobal.col[j], &br[j]));
764: }
765: bl = NULL;
766: for (i = 0; i < bA->nr; i++) {
767: if (l) PetscCall(VecGetSubVector(l, bA->isglobal.row[i], &bl));
768: for (j = 0; j < bA->nc; j++) {
769: if (bA->m[i][j]) PetscCall(MatDiagonalScale(bA->m[i][j], bl, br[j]));
770: }
771: if (l) PetscCall(VecRestoreSubVector(l, bA->isglobal.row[i], &bl));
772: }
773: if (r) {
774: for (j = 0; j < bA->nc; j++) PetscCall(VecRestoreSubVector(r, bA->isglobal.col[j], &br[j]));
775: }
776: PetscCall(PetscFree(br));
777: PetscFunctionReturn(PETSC_SUCCESS);
778: }
780: static PetscErrorCode MatScale_Nest(Mat A, PetscScalar a)
781: {
782: Mat_Nest *bA = (Mat_Nest *)A->data;
783: PetscInt i, j;
785: PetscFunctionBegin;
786: for (i = 0; i < bA->nr; i++) {
787: for (j = 0; j < bA->nc; j++) {
788: if (bA->m[i][j]) PetscCall(MatScale(bA->m[i][j], a));
789: }
790: }
791: PetscFunctionReturn(PETSC_SUCCESS);
792: }
794: static PetscErrorCode MatShift_Nest(Mat A, PetscScalar a)
795: {
796: Mat_Nest *bA = (Mat_Nest *)A->data;
797: PetscInt i;
798: PetscBool nnzstate = PETSC_FALSE;
800: PetscFunctionBegin;
801: for (i = 0; i < bA->nr; i++) {
802: PetscObjectState subnnzstate = 0;
803: 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);
804: PetscCall(MatShift(bA->m[i][i], a));
805: PetscCall(MatGetNonzeroState(bA->m[i][i], &subnnzstate));
806: nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i * bA->nc + i] != subnnzstate);
807: bA->nnzstate[i * bA->nc + i] = subnnzstate;
808: }
809: if (nnzstate) A->nonzerostate++;
810: PetscFunctionReturn(PETSC_SUCCESS);
811: }
813: static PetscErrorCode MatDiagonalSet_Nest(Mat A, Vec D, InsertMode is)
814: {
815: Mat_Nest *bA = (Mat_Nest *)A->data;
816: PetscInt i;
817: PetscBool nnzstate = PETSC_FALSE;
819: PetscFunctionBegin;
820: for (i = 0; i < bA->nr; i++) {
821: PetscObjectState subnnzstate = 0;
822: Vec bv;
823: PetscCall(VecGetSubVector(D, bA->isglobal.row[i], &bv));
824: if (bA->m[i][i]) {
825: PetscCall(MatDiagonalSet(bA->m[i][i], bv, is));
826: PetscCall(MatGetNonzeroState(bA->m[i][i], &subnnzstate));
827: }
828: PetscCall(VecRestoreSubVector(D, bA->isglobal.row[i], &bv));
829: nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i * bA->nc + i] != subnnzstate);
830: bA->nnzstate[i * bA->nc + i] = subnnzstate;
831: }
832: if (nnzstate) A->nonzerostate++;
833: PetscFunctionReturn(PETSC_SUCCESS);
834: }
836: static PetscErrorCode MatSetRandom_Nest(Mat A, PetscRandom rctx)
837: {
838: Mat_Nest *bA = (Mat_Nest *)A->data;
839: PetscInt i, j;
841: PetscFunctionBegin;
842: for (i = 0; i < bA->nr; i++) {
843: for (j = 0; j < bA->nc; j++) {
844: if (bA->m[i][j]) PetscCall(MatSetRandom(bA->m[i][j], rctx));
845: }
846: }
847: PetscFunctionReturn(PETSC_SUCCESS);
848: }
850: static PetscErrorCode MatCreateVecs_Nest(Mat A, Vec *right, Vec *left)
851: {
852: Mat_Nest *bA = (Mat_Nest *)A->data;
853: Vec *L, *R;
854: MPI_Comm comm;
855: PetscInt i, j;
857: PetscFunctionBegin;
858: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
859: if (right) {
860: /* allocate R */
861: PetscCall(PetscMalloc1(bA->nc, &R));
862: /* Create the right vectors */
863: for (j = 0; j < bA->nc; j++) {
864: for (i = 0; i < bA->nr; i++) {
865: if (bA->m[i][j]) {
866: PetscCall(MatCreateVecs(bA->m[i][j], &R[j], NULL));
867: break;
868: }
869: }
870: PetscCheck(i != bA->nr, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null column.");
871: }
872: PetscCall(VecCreateNest(comm, bA->nc, bA->isglobal.col, R, right));
873: /* hand back control to the nest vector */
874: for (j = 0; j < bA->nc; j++) PetscCall(VecDestroy(&R[j]));
875: PetscCall(PetscFree(R));
876: }
878: if (left) {
879: /* allocate L */
880: PetscCall(PetscMalloc1(bA->nr, &L));
881: /* Create the left vectors */
882: for (i = 0; i < bA->nr; i++) {
883: for (j = 0; j < bA->nc; j++) {
884: if (bA->m[i][j]) {
885: PetscCall(MatCreateVecs(bA->m[i][j], NULL, &L[i]));
886: break;
887: }
888: }
889: PetscCheck(j != bA->nc, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null row.");
890: }
892: PetscCall(VecCreateNest(comm, bA->nr, bA->isglobal.row, L, left));
893: for (i = 0; i < bA->nr; i++) PetscCall(VecDestroy(&L[i]));
895: PetscCall(PetscFree(L));
896: }
897: PetscFunctionReturn(PETSC_SUCCESS);
898: }
900: static PetscErrorCode MatView_Nest(Mat A, PetscViewer viewer)
901: {
902: Mat_Nest *bA = (Mat_Nest *)A->data;
903: PetscBool isascii, viewSub = PETSC_FALSE;
904: PetscInt i, j;
906: PetscFunctionBegin;
907: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
908: if (isascii) {
909: PetscCall(PetscOptionsGetBool(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_view_nest_sub", &viewSub, NULL));
910: PetscCall(PetscViewerASCIIPrintf(viewer, "Matrix object:\n"));
911: PetscCall(PetscViewerASCIIPushTab(viewer));
912: PetscCall(PetscViewerASCIIPrintf(viewer, "type=nest, rows=%" PetscInt_FMT ", cols=%" PetscInt_FMT "\n", bA->nr, bA->nc));
914: PetscCall(PetscViewerASCIIPrintf(viewer, "MatNest structure:\n"));
915: for (i = 0; i < bA->nr; i++) {
916: for (j = 0; j < bA->nc; j++) {
917: MatType type;
918: char name[256] = "", prefix[256] = "";
919: PetscInt NR, NC;
920: PetscBool isNest = PETSC_FALSE;
922: if (!bA->m[i][j]) {
923: PetscCall(PetscViewerASCIIPrintf(viewer, "(%" PetscInt_FMT ",%" PetscInt_FMT ") : NULL\n", i, j));
924: continue;
925: }
926: PetscCall(MatGetSize(bA->m[i][j], &NR, &NC));
927: PetscCall(MatGetType(bA->m[i][j], &type));
928: if (((PetscObject)bA->m[i][j])->name) PetscCall(PetscSNPrintf(name, sizeof(name), "name=\"%s\", ", ((PetscObject)bA->m[i][j])->name));
929: if (((PetscObject)bA->m[i][j])->prefix) PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "prefix=\"%s\", ", ((PetscObject)bA->m[i][j])->prefix));
930: PetscCall(PetscObjectTypeCompare((PetscObject)bA->m[i][j], MATNEST, &isNest));
932: PetscCall(PetscViewerASCIIPrintf(viewer, "(%" PetscInt_FMT ",%" PetscInt_FMT ") : %s%stype=%s, rows=%" PetscInt_FMT ", cols=%" PetscInt_FMT "\n", i, j, name, prefix, type, NR, NC));
934: if (isNest || viewSub) {
935: PetscCall(PetscViewerASCIIPushTab(viewer)); /* push1 */
936: PetscCall(MatView(bA->m[i][j], viewer));
937: PetscCall(PetscViewerASCIIPopTab(viewer)); /* pop1 */
938: }
939: }
940: }
941: PetscCall(PetscViewerASCIIPopTab(viewer)); /* pop0 */
942: }
943: PetscFunctionReturn(PETSC_SUCCESS);
944: }
946: static PetscErrorCode MatZeroEntries_Nest(Mat A)
947: {
948: Mat_Nest *bA = (Mat_Nest *)A->data;
949: PetscInt i, j;
951: PetscFunctionBegin;
952: for (i = 0; i < bA->nr; i++) {
953: for (j = 0; j < bA->nc; j++) {
954: if (!bA->m[i][j]) continue;
955: PetscCall(MatZeroEntries(bA->m[i][j]));
956: }
957: }
958: PetscFunctionReturn(PETSC_SUCCESS);
959: }
961: static PetscErrorCode MatCopy_Nest(Mat A, Mat B, MatStructure str)
962: {
963: Mat_Nest *bA = (Mat_Nest *)A->data, *bB = (Mat_Nest *)B->data;
964: PetscInt i, j, nr = bA->nr, nc = bA->nc;
965: PetscBool nnzstate = PETSC_FALSE;
967: PetscFunctionBegin;
968: 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);
969: for (i = 0; i < nr; i++) {
970: for (j = 0; j < nc; j++) {
971: PetscObjectState subnnzstate = 0;
972: if (bA->m[i][j] && bB->m[i][j]) {
973: PetscCall(MatCopy(bA->m[i][j], bB->m[i][j], str));
974: } 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);
975: PetscCall(MatGetNonzeroState(bB->m[i][j], &subnnzstate));
976: nnzstate = (PetscBool)(nnzstate || bB->nnzstate[i * nc + j] != subnnzstate);
977: bB->nnzstate[i * nc + j] = subnnzstate;
978: }
979: }
980: if (nnzstate) B->nonzerostate++;
981: PetscFunctionReturn(PETSC_SUCCESS);
982: }
984: static PetscErrorCode MatAXPY_Nest(Mat Y, PetscScalar a, Mat X, MatStructure str)
985: {
986: Mat_Nest *bY = (Mat_Nest *)Y->data, *bX = (Mat_Nest *)X->data;
987: PetscInt i, j, nr = bY->nr, nc = bY->nc;
988: PetscBool nnzstate = PETSC_FALSE;
990: PetscFunctionBegin;
991: 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);
992: for (i = 0; i < nr; i++) {
993: for (j = 0; j < nc; j++) {
994: PetscObjectState subnnzstate = 0;
995: if (bY->m[i][j] && bX->m[i][j]) {
996: PetscCall(MatAXPY(bY->m[i][j], a, bX->m[i][j], str));
997: } else if (bX->m[i][j]) {
998: Mat M;
1000: 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);
1001: PetscCall(MatDuplicate(bX->m[i][j], MAT_COPY_VALUES, &M));
1002: PetscCall(MatNestSetSubMat(Y, i, j, M));
1003: PetscCall(MatDestroy(&M));
1004: }
1005: if (bY->m[i][j]) PetscCall(MatGetNonzeroState(bY->m[i][j], &subnnzstate));
1006: nnzstate = (PetscBool)(nnzstate || bY->nnzstate[i * nc + j] != subnnzstate);
1007: bY->nnzstate[i * nc + j] = subnnzstate;
1008: }
1009: }
1010: if (nnzstate) Y->nonzerostate++;
1011: PetscFunctionReturn(PETSC_SUCCESS);
1012: }
1014: static PetscErrorCode MatDuplicate_Nest(Mat A, MatDuplicateOption op, Mat *B)
1015: {
1016: Mat_Nest *bA = (Mat_Nest *)A->data;
1017: Mat *b;
1018: PetscInt i, j, nr = bA->nr, nc = bA->nc;
1020: PetscFunctionBegin;
1021: PetscCall(PetscMalloc1(nr * nc, &b));
1022: for (i = 0; i < nr; i++) {
1023: for (j = 0; j < nc; j++) {
1024: if (bA->m[i][j]) {
1025: PetscCall(MatDuplicate(bA->m[i][j], op, &b[i * nc + j]));
1026: } else {
1027: b[i * nc + j] = NULL;
1028: }
1029: }
1030: }
1031: PetscCall(MatCreateNest(PetscObjectComm((PetscObject)A), nr, bA->isglobal.row, nc, bA->isglobal.col, b, B));
1032: /* Give the new MatNest exclusive ownership */
1033: for (i = 0; i < nr * nc; i++) PetscCall(MatDestroy(&b[i]));
1034: PetscCall(PetscFree(b));
1036: PetscCall(MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY));
1037: PetscCall(MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY));
1038: PetscFunctionReturn(PETSC_SUCCESS);
1039: }
1041: /* nest api */
1042: static PetscErrorCode MatNestGetSubMat_Nest(Mat A, PetscInt idxm, PetscInt jdxm, Mat *mat)
1043: {
1044: Mat_Nest *bA = (Mat_Nest *)A->data;
1046: PetscFunctionBegin;
1047: PetscCheck(idxm < bA->nr, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, idxm, bA->nr - 1);
1048: PetscCheck(jdxm < bA->nc, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Col too large: row %" PetscInt_FMT " max %" PetscInt_FMT, jdxm, bA->nc - 1);
1049: *mat = bA->m[idxm][jdxm];
1050: PetscFunctionReturn(PETSC_SUCCESS);
1051: }
1053: /*@
1054: MatNestGetSubMat - Returns a single, sub-matrix from a `MATNEST`
1056: Not Collective
1058: Input Parameters:
1059: + A - `MATNEST` matrix
1060: . idxm - index of the matrix within the nest matrix
1061: - jdxm - index of the matrix within the nest matrix
1063: Output Parameter:
1064: . sub - matrix at index `idxm`, `jdxm` within the nest matrix
1066: Level: developer
1068: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSize()`, `MatNestGetSubMats()`, `MatCreateNest()`, `MatNestSetSubMat()`,
1069: `MatNestGetLocalISs()`, `MatNestGetISs()`
1070: @*/
1071: PetscErrorCode MatNestGetSubMat(Mat A, PetscInt idxm, PetscInt jdxm, Mat *sub)
1072: {
1073: PetscFunctionBegin;
1074: PetscUseMethod(A, "MatNestGetSubMat_C", (Mat, PetscInt, PetscInt, Mat *), (A, idxm, jdxm, sub));
1075: PetscFunctionReturn(PETSC_SUCCESS);
1076: }
1078: static PetscErrorCode MatNestSetSubMat_Nest(Mat A, PetscInt idxm, PetscInt jdxm, Mat mat)
1079: {
1080: Mat_Nest *bA = (Mat_Nest *)A->data;
1081: PetscInt m, n, M, N, mi, ni, Mi, Ni;
1083: PetscFunctionBegin;
1084: PetscCheck(idxm < bA->nr, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, idxm, bA->nr - 1);
1085: PetscCheck(jdxm < bA->nc, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Col too large: row %" PetscInt_FMT " max %" PetscInt_FMT, jdxm, bA->nc - 1);
1086: PetscCall(MatGetLocalSize(mat, &m, &n));
1087: PetscCall(MatGetSize(mat, &M, &N));
1088: PetscCall(ISGetLocalSize(bA->isglobal.row[idxm], &mi));
1089: PetscCall(ISGetSize(bA->isglobal.row[idxm], &Mi));
1090: PetscCall(ISGetLocalSize(bA->isglobal.col[jdxm], &ni));
1091: PetscCall(ISGetSize(bA->isglobal.col[jdxm], &Ni));
1092: 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);
1093: 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);
1095: /* do not increase object state */
1096: if (mat == bA->m[idxm][jdxm]) PetscFunctionReturn(PETSC_SUCCESS);
1098: PetscCall(PetscObjectReference((PetscObject)mat));
1099: PetscCall(MatDestroy(&bA->m[idxm][jdxm]));
1100: bA->m[idxm][jdxm] = mat;
1101: PetscCall(PetscObjectStateIncrease((PetscObject)A));
1102: PetscCall(MatGetNonzeroState(mat, &bA->nnzstate[idxm * bA->nc + jdxm]));
1103: A->nonzerostate++;
1104: PetscFunctionReturn(PETSC_SUCCESS);
1105: }
1107: /*@
1108: MatNestSetSubMat - Set a single submatrix in the `MATNEST`
1110: Logically Collective
1112: Input Parameters:
1113: + A - `MATNEST` matrix
1114: . idxm - index of the matrix within the nest matrix
1115: . jdxm - index of the matrix within the nest matrix
1116: - sub - matrix at index `idxm`, `jdxm` within the nest matrix
1118: Level: developer
1120: Notes:
1121: The new submatrix must have the same size and communicator as that block of the nest.
1123: This increments the reference count of the submatrix.
1125: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestSetSubMats()`, `MatNestGetSubMats()`, `MatNestGetLocalISs()`, `MatCreateNest()`,
1126: `MatNestGetSubMat()`, `MatNestGetISs()`, `MatNestGetSize()`
1127: @*/
1128: PetscErrorCode MatNestSetSubMat(Mat A, PetscInt idxm, PetscInt jdxm, Mat sub)
1129: {
1130: PetscFunctionBegin;
1131: PetscUseMethod(A, "MatNestSetSubMat_C", (Mat, PetscInt, PetscInt, Mat), (A, idxm, jdxm, sub));
1132: PetscFunctionReturn(PETSC_SUCCESS);
1133: }
1135: static PetscErrorCode MatNestGetSubMats_Nest(Mat A, PetscInt *M, PetscInt *N, Mat ***mat)
1136: {
1137: Mat_Nest *bA = (Mat_Nest *)A->data;
1139: PetscFunctionBegin;
1140: if (M) *M = bA->nr;
1141: if (N) *N = bA->nc;
1142: if (mat) *mat = bA->m;
1143: PetscFunctionReturn(PETSC_SUCCESS);
1144: }
1146: /*@C
1147: MatNestGetSubMats - Returns the entire two dimensional array of matrices defining a `MATNEST` matrix.
1149: Not Collective
1151: Input Parameter:
1152: . A - nest matrix
1154: Output Parameters:
1155: + M - number of rows in the nest matrix
1156: . N - number of cols in the nest matrix
1157: - mat - array of matrices
1159: Level: developer
1161: Note:
1162: The user should not free the array `mat`.
1164: Fortran Notes:
1165: This routine has a calling sequence
1166: $ call MatNestGetSubMats(A, M, N, mat, ierr)
1167: where the space allocated for the optional argument `mat` is assumed large enough (if provided).
1168: Matrices in `mat` are returned in row-major order, see `MatCreateNest()` for an example.
1170: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSize()`, `MatNestGetSubMat()`, `MatNestGetLocalISs()`, `MatCreateNest()`,
1171: `MatNestSetSubMats()`, `MatNestGetISs()`, `MatNestSetSubMat()`
1172: @*/
1173: PetscErrorCode MatNestGetSubMats(Mat A, PetscInt *M, PetscInt *N, Mat ***mat)
1174: {
1175: PetscFunctionBegin;
1176: PetscUseMethod(A, "MatNestGetSubMats_C", (Mat, PetscInt *, PetscInt *, Mat ***), (A, M, N, mat));
1177: PetscFunctionReturn(PETSC_SUCCESS);
1178: }
1180: static PetscErrorCode MatNestGetSize_Nest(Mat A, PetscInt *M, PetscInt *N)
1181: {
1182: Mat_Nest *bA = (Mat_Nest *)A->data;
1184: PetscFunctionBegin;
1185: if (M) *M = bA->nr;
1186: if (N) *N = bA->nc;
1187: PetscFunctionReturn(PETSC_SUCCESS);
1188: }
1190: /*@
1191: MatNestGetSize - Returns the size of the `MATNEST` matrix.
1193: Not Collective
1195: Input Parameter:
1196: . A - `MATNEST` matrix
1198: Output Parameters:
1199: + M - number of rows in the nested mat
1200: - N - number of cols in the nested mat
1202: Level: developer
1204: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSubMat()`, `MatNestGetSubMats()`, `MatCreateNest()`, `MatNestGetLocalISs()`,
1205: `MatNestGetISs()`
1206: @*/
1207: PetscErrorCode MatNestGetSize(Mat A, PetscInt *M, PetscInt *N)
1208: {
1209: PetscFunctionBegin;
1210: PetscUseMethod(A, "MatNestGetSize_C", (Mat, PetscInt *, PetscInt *), (A, M, N));
1211: PetscFunctionReturn(PETSC_SUCCESS);
1212: }
1214: static PetscErrorCode MatNestGetISs_Nest(Mat A, IS rows[], IS cols[])
1215: {
1216: Mat_Nest *vs = (Mat_Nest *)A->data;
1217: PetscInt i;
1219: PetscFunctionBegin;
1220: if (rows)
1221: for (i = 0; i < vs->nr; i++) rows[i] = vs->isglobal.row[i];
1222: if (cols)
1223: for (i = 0; i < vs->nc; i++) cols[i] = vs->isglobal.col[i];
1224: PetscFunctionReturn(PETSC_SUCCESS);
1225: }
1227: /*@C
1228: MatNestGetISs - Returns the index sets partitioning the row and column spaces of a `MATNEST`
1230: Not Collective
1232: Input Parameter:
1233: . A - `MATNEST` matrix
1235: Output Parameters:
1236: + rows - array of row index sets
1237: - cols - array of column index sets
1239: Level: advanced
1241: Note:
1242: The user must have allocated arrays of the correct size. The reference count is not increased on the returned `IS`s.
1244: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSubMat()`, `MatNestGetSubMats()`, `MatNestGetSize()`, `MatNestGetLocalISs()`,
1245: `MatCreateNest()`, `MatNestSetSubMats()`
1246: @*/
1247: PetscErrorCode MatNestGetISs(Mat A, IS rows[], IS cols[])
1248: {
1249: PetscFunctionBegin;
1251: PetscUseMethod(A, "MatNestGetISs_C", (Mat, IS[], IS[]), (A, rows, cols));
1252: PetscFunctionReturn(PETSC_SUCCESS);
1253: }
1255: static PetscErrorCode MatNestGetLocalISs_Nest(Mat A, IS rows[], IS cols[])
1256: {
1257: Mat_Nest *vs = (Mat_Nest *)A->data;
1258: PetscInt i;
1260: PetscFunctionBegin;
1261: if (rows)
1262: for (i = 0; i < vs->nr; i++) rows[i] = vs->islocal.row[i];
1263: if (cols)
1264: for (i = 0; i < vs->nc; i++) cols[i] = vs->islocal.col[i];
1265: PetscFunctionReturn(PETSC_SUCCESS);
1266: }
1268: /*@C
1269: MatNestGetLocalISs - Returns the index sets partitioning the row and column spaces of a `MATNEST`
1271: Not Collective
1273: Input Parameter:
1274: . A - `MATNEST` matrix
1276: Output Parameters:
1277: + rows - array of row index sets (or `NULL` to ignore)
1278: - cols - array of column index sets (or `NULL` to ignore)
1280: Level: advanced
1282: Note:
1283: The user must have allocated arrays of the correct size. The reference count is not increased on the returned `IS`s.
1285: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSubMat()`, `MatNestGetSubMats()`, `MatNestGetSize()`, `MatNestGetISs()`, `MatCreateNest()`,
1286: `MatNestSetSubMats()`, `MatNestSetSubMat()`
1287: @*/
1288: PetscErrorCode MatNestGetLocalISs(Mat A, IS rows[], IS cols[])
1289: {
1290: PetscFunctionBegin;
1292: PetscUseMethod(A, "MatNestGetLocalISs_C", (Mat, IS[], IS[]), (A, rows, cols));
1293: PetscFunctionReturn(PETSC_SUCCESS);
1294: }
1296: static PetscErrorCode MatNestSetVecType_Nest(Mat A, VecType vtype)
1297: {
1298: PetscBool flg;
1300: PetscFunctionBegin;
1301: PetscCall(PetscStrcmp(vtype, VECNEST, &flg));
1302: /* In reality, this only distinguishes VECNEST and "other" */
1303: if (flg) A->ops->getvecs = MatCreateVecs_Nest;
1304: else A->ops->getvecs = (PetscErrorCode(*)(Mat, Vec *, Vec *))0;
1305: PetscFunctionReturn(PETSC_SUCCESS);
1306: }
1308: /*@C
1309: MatNestSetVecType - Sets the type of `Vec` returned by `MatCreateVecs()`
1311: Not Collective
1313: Input Parameters:
1314: + A - `MATNEST` matrix
1315: - vtype - `VecType` to use for creating vectors
1317: Level: developer
1319: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatCreateVecs()`, `MatCreateNest()`, `VecType`
1320: @*/
1321: PetscErrorCode MatNestSetVecType(Mat A, VecType vtype)
1322: {
1323: PetscFunctionBegin;
1324: PetscTryMethod(A, "MatNestSetVecType_C", (Mat, VecType), (A, vtype));
1325: PetscFunctionReturn(PETSC_SUCCESS);
1326: }
1328: static PetscErrorCode MatNestSetSubMats_Nest(Mat A, PetscInt nr, const IS is_row[], PetscInt nc, const IS is_col[], const Mat a[])
1329: {
1330: Mat_Nest *s = (Mat_Nest *)A->data;
1331: PetscInt i, j, m, n, M, N;
1332: PetscBool cong, isstd, sametype = PETSC_FALSE;
1333: VecType vtype, type;
1335: PetscFunctionBegin;
1336: PetscCall(MatReset_Nest(A));
1338: s->nr = nr;
1339: s->nc = nc;
1341: /* Create space for submatrices */
1342: PetscCall(PetscMalloc1(nr, &s->m));
1343: for (i = 0; i < nr; i++) PetscCall(PetscMalloc1(nc, &s->m[i]));
1344: for (i = 0; i < nr; i++) {
1345: for (j = 0; j < nc; j++) {
1346: s->m[i][j] = a[i * nc + j];
1347: if (a[i * nc + j]) PetscCall(PetscObjectReference((PetscObject)a[i * nc + j]));
1348: }
1349: }
1350: PetscCall(MatGetVecType(A, &vtype));
1351: PetscCall(PetscStrcmp(vtype, VECSTANDARD, &isstd));
1352: if (isstd) {
1353: /* check if all blocks have the same vectype */
1354: vtype = NULL;
1355: for (i = 0; i < nr; i++) {
1356: for (j = 0; j < nc; j++) {
1357: if (a[i * nc + j]) {
1358: if (!vtype) { /* first visited block */
1359: PetscCall(MatGetVecType(a[i * nc + j], &vtype));
1360: sametype = PETSC_TRUE;
1361: } else if (sametype) {
1362: PetscCall(MatGetVecType(a[i * nc + j], &type));
1363: PetscCall(PetscStrcmp(vtype, type, &sametype));
1364: }
1365: }
1366: }
1367: }
1368: if (sametype) { /* propagate vectype */
1369: PetscCall(MatSetVecType(A, vtype));
1370: }
1371: }
1373: PetscCall(MatSetUp_NestIS_Private(A, nr, is_row, nc, is_col));
1375: PetscCall(PetscMalloc1(nr, &s->row_len));
1376: PetscCall(PetscMalloc1(nc, &s->col_len));
1377: for (i = 0; i < nr; i++) s->row_len[i] = -1;
1378: for (j = 0; j < nc; j++) s->col_len[j] = -1;
1380: PetscCall(PetscCalloc1(nr * nc, &s->nnzstate));
1381: for (i = 0; i < nr; i++) {
1382: for (j = 0; j < nc; j++) {
1383: if (s->m[i][j]) PetscCall(MatGetNonzeroState(s->m[i][j], &s->nnzstate[i * nc + j]));
1384: }
1385: }
1387: PetscCall(MatNestGetSizes_Private(A, &m, &n, &M, &N));
1389: PetscCall(PetscLayoutSetSize(A->rmap, M));
1390: PetscCall(PetscLayoutSetLocalSize(A->rmap, m));
1391: PetscCall(PetscLayoutSetSize(A->cmap, N));
1392: PetscCall(PetscLayoutSetLocalSize(A->cmap, n));
1394: PetscCall(PetscLayoutSetUp(A->rmap));
1395: PetscCall(PetscLayoutSetUp(A->cmap));
1397: /* disable operations that are not supported for non-square matrices,
1398: or matrices for which is_row != is_col */
1399: PetscCall(MatHasCongruentLayouts(A, &cong));
1400: if (cong && nr != nc) cong = PETSC_FALSE;
1401: if (cong) {
1402: for (i = 0; cong && i < nr; i++) PetscCall(ISEqualUnsorted(s->isglobal.row[i], s->isglobal.col[i], &cong));
1403: }
1404: if (!cong) {
1405: A->ops->missingdiagonal = NULL;
1406: A->ops->getdiagonal = NULL;
1407: A->ops->shift = NULL;
1408: A->ops->diagonalset = NULL;
1409: }
1411: PetscCall(PetscCalloc2(nr, &s->left, nc, &s->right));
1412: PetscCall(PetscObjectStateIncrease((PetscObject)A));
1413: A->nonzerostate++;
1414: PetscFunctionReturn(PETSC_SUCCESS);
1415: }
1417: /*@
1418: MatNestSetSubMats - Sets the nested submatrices in a `MATNEST`
1420: Collective
1422: Input Parameters:
1423: + A - `MATNEST` matrix
1424: . nr - number of nested row blocks
1425: . is_row - index sets for each nested row block, or `NULL` to make contiguous
1426: . nc - number of nested column blocks
1427: . is_col - index sets for each nested column block, or `NULL` to make contiguous
1428: - a - array of nr*nc submatrices, empty submatrices can be passed using `NULL`
1430: Level: advanced
1432: Notes:
1433: This always resets any submatrix information previously set
1435: In both C and Fortran, `a` must be a row-major order array containing the matrices. See
1436: `MatCreateNest()` for an example.
1438: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatCreateNest()`, `MatNestSetSubMat()`, `MatNestGetSubMat()`, `MatNestGetSubMats()`
1439: @*/
1440: PetscErrorCode MatNestSetSubMats(Mat A, PetscInt nr, const IS is_row[], PetscInt nc, const IS is_col[], const Mat a[])
1441: {
1442: PetscInt i;
1444: PetscFunctionBegin;
1446: PetscCheck(nr >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Number of rows cannot be negative");
1447: if (nr && is_row) {
1448: PetscAssertPointer(is_row, 3);
1450: }
1451: PetscCheck(nc >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Number of columns cannot be negative");
1452: if (nc && is_col) {
1453: PetscAssertPointer(is_col, 5);
1455: }
1456: if (nr * nc > 0) PetscAssertPointer(a, 6);
1457: PetscUseMethod(A, "MatNestSetSubMats_C", (Mat, PetscInt, const IS[], PetscInt, const IS[], const Mat[]), (A, nr, is_row, nc, is_col, a));
1458: PetscFunctionReturn(PETSC_SUCCESS);
1459: }
1461: static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A, PetscInt n, const IS islocal[], const IS isglobal[], PetscBool colflg, ISLocalToGlobalMapping *ltog)
1462: {
1463: PetscBool flg;
1464: PetscInt i, j, m, mi, *ix;
1466: PetscFunctionBegin;
1467: *ltog = NULL;
1468: for (i = 0, m = 0, flg = PETSC_FALSE; i < n; i++) {
1469: if (islocal[i]) {
1470: PetscCall(ISGetLocalSize(islocal[i], &mi));
1471: flg = PETSC_TRUE; /* We found a non-trivial entry */
1472: } else {
1473: PetscCall(ISGetLocalSize(isglobal[i], &mi));
1474: }
1475: m += mi;
1476: }
1477: if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
1479: PetscCall(PetscMalloc1(m, &ix));
1480: for (i = 0, m = 0; i < n; i++) {
1481: ISLocalToGlobalMapping smap = NULL;
1482: Mat sub = NULL;
1483: PetscSF sf;
1484: PetscLayout map;
1485: const PetscInt *ix2;
1487: if (!colflg) {
1488: PetscCall(MatNestFindNonzeroSubMatRow(A, i, &sub));
1489: } else {
1490: PetscCall(MatNestFindNonzeroSubMatCol(A, i, &sub));
1491: }
1492: if (sub) {
1493: if (!colflg) {
1494: PetscCall(MatGetLocalToGlobalMapping(sub, &smap, NULL));
1495: } else {
1496: PetscCall(MatGetLocalToGlobalMapping(sub, NULL, &smap));
1497: }
1498: }
1499: /*
1500: Now we need to extract the monolithic global indices that correspond to the given split global indices.
1501: In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces.
1502: */
1503: PetscCall(ISGetIndices(isglobal[i], &ix2));
1504: if (islocal[i]) {
1505: PetscInt *ilocal, *iremote;
1506: PetscInt mil, nleaves;
1508: PetscCall(ISGetLocalSize(islocal[i], &mi));
1509: PetscCheck(smap, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing local to global map");
1510: for (j = 0; j < mi; j++) ix[m + j] = j;
1511: PetscCall(ISLocalToGlobalMappingApply(smap, mi, ix + m, ix + m));
1513: /* PetscSFSetGraphLayout does not like negative indices */
1514: PetscCall(PetscMalloc2(mi, &ilocal, mi, &iremote));
1515: for (j = 0, nleaves = 0; j < mi; j++) {
1516: if (ix[m + j] < 0) continue;
1517: ilocal[nleaves] = j;
1518: iremote[nleaves] = ix[m + j];
1519: nleaves++;
1520: }
1521: PetscCall(ISGetLocalSize(isglobal[i], &mil));
1522: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)A), &sf));
1523: PetscCall(PetscLayoutCreate(PetscObjectComm((PetscObject)A), &map));
1524: PetscCall(PetscLayoutSetLocalSize(map, mil));
1525: PetscCall(PetscLayoutSetUp(map));
1526: PetscCall(PetscSFSetGraphLayout(sf, map, nleaves, ilocal, PETSC_USE_POINTER, iremote));
1527: PetscCall(PetscLayoutDestroy(&map));
1528: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, ix2, ix + m, MPI_REPLACE));
1529: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, ix2, ix + m, MPI_REPLACE));
1530: PetscCall(PetscSFDestroy(&sf));
1531: PetscCall(PetscFree2(ilocal, iremote));
1532: } else {
1533: PetscCall(ISGetLocalSize(isglobal[i], &mi));
1534: for (j = 0; j < mi; j++) ix[m + j] = ix2[i];
1535: }
1536: PetscCall(ISRestoreIndices(isglobal[i], &ix2));
1537: m += mi;
1538: }
1539: PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A), 1, m, ix, PETSC_OWN_POINTER, ltog));
1540: PetscFunctionReturn(PETSC_SUCCESS);
1541: }
1543: /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */
1544: /*
1545: nprocessors = NP
1546: Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1))
1547: proc 0: => (g_0,h_0,)
1548: proc 1: => (g_1,h_1,)
1549: ...
1550: proc nprocs-1: => (g_NP-1,h_NP-1,)
1552: proc 0: proc 1: proc nprocs-1:
1553: is[0] = (0,1,2,...,nlocal(g_0)-1) (0,1,...,nlocal(g_1)-1) (0,1,...,nlocal(g_NP-1))
1555: proc 0:
1556: is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1)
1557: proc 1:
1558: is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1)
1560: proc NP-1:
1561: is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1)
1562: */
1563: static PetscErrorCode MatSetUp_NestIS_Private(Mat A, PetscInt nr, const IS is_row[], PetscInt nc, const IS is_col[])
1564: {
1565: Mat_Nest *vs = (Mat_Nest *)A->data;
1566: PetscInt i, j, offset, n, nsum, bs;
1567: Mat sub = NULL;
1569: PetscFunctionBegin;
1570: PetscCall(PetscMalloc1(nr, &vs->isglobal.row));
1571: PetscCall(PetscMalloc1(nc, &vs->isglobal.col));
1572: if (is_row) { /* valid IS is passed in */
1573: /* refs on is[] are incremented */
1574: for (i = 0; i < vs->nr; i++) {
1575: PetscCall(PetscObjectReference((PetscObject)is_row[i]));
1577: vs->isglobal.row[i] = is_row[i];
1578: }
1579: } else { /* Create the ISs by inspecting sizes of a submatrix in each row */
1580: nsum = 0;
1581: for (i = 0; i < vs->nr; i++) { /* Add up the local sizes to compute the aggregate offset */
1582: PetscCall(MatNestFindNonzeroSubMatRow(A, i, &sub));
1583: PetscCheck(sub, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "No nonzero submatrix in row %" PetscInt_FMT, i);
1584: PetscCall(MatGetLocalSize(sub, &n, NULL));
1585: PetscCheck(n >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Sizes have not yet been set for submatrix");
1586: nsum += n;
1587: }
1588: PetscCallMPI(MPI_Scan(&nsum, &offset, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)A)));
1589: offset -= nsum;
1590: for (i = 0; i < vs->nr; i++) {
1591: PetscCall(MatNestFindNonzeroSubMatRow(A, i, &sub));
1592: PetscCall(MatGetLocalSize(sub, &n, NULL));
1593: PetscCall(MatGetBlockSizes(sub, &bs, NULL));
1594: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)sub), n, offset, 1, &vs->isglobal.row[i]));
1595: PetscCall(ISSetBlockSize(vs->isglobal.row[i], bs));
1596: offset += n;
1597: }
1598: }
1600: if (is_col) { /* valid IS is passed in */
1601: /* refs on is[] are incremented */
1602: for (j = 0; j < vs->nc; j++) {
1603: PetscCall(PetscObjectReference((PetscObject)is_col[j]));
1605: vs->isglobal.col[j] = is_col[j];
1606: }
1607: } else { /* Create the ISs by inspecting sizes of a submatrix in each column */
1608: offset = A->cmap->rstart;
1609: nsum = 0;
1610: for (j = 0; j < vs->nc; j++) {
1611: PetscCall(MatNestFindNonzeroSubMatCol(A, j, &sub));
1612: PetscCheck(sub, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "No nonzero submatrix in column %" PetscInt_FMT, i);
1613: PetscCall(MatGetLocalSize(sub, NULL, &n));
1614: PetscCheck(n >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Sizes have not yet been set for submatrix");
1615: nsum += n;
1616: }
1617: PetscCallMPI(MPI_Scan(&nsum, &offset, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)A)));
1618: offset -= nsum;
1619: for (j = 0; j < vs->nc; j++) {
1620: PetscCall(MatNestFindNonzeroSubMatCol(A, j, &sub));
1621: PetscCall(MatGetLocalSize(sub, NULL, &n));
1622: PetscCall(MatGetBlockSizes(sub, NULL, &bs));
1623: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)sub), n, offset, 1, &vs->isglobal.col[j]));
1624: PetscCall(ISSetBlockSize(vs->isglobal.col[j], bs));
1625: offset += n;
1626: }
1627: }
1629: /* Set up the local ISs */
1630: PetscCall(PetscMalloc1(vs->nr, &vs->islocal.row));
1631: PetscCall(PetscMalloc1(vs->nc, &vs->islocal.col));
1632: for (i = 0, offset = 0; i < vs->nr; i++) {
1633: IS isloc;
1634: ISLocalToGlobalMapping rmap = NULL;
1635: PetscInt nlocal, bs;
1636: PetscCall(MatNestFindNonzeroSubMatRow(A, i, &sub));
1637: if (sub) PetscCall(MatGetLocalToGlobalMapping(sub, &rmap, NULL));
1638: if (rmap) {
1639: PetscCall(MatGetBlockSizes(sub, &bs, NULL));
1640: PetscCall(ISLocalToGlobalMappingGetSize(rmap, &nlocal));
1641: PetscCall(ISCreateStride(PETSC_COMM_SELF, nlocal, offset, 1, &isloc));
1642: PetscCall(ISSetBlockSize(isloc, bs));
1643: } else {
1644: nlocal = 0;
1645: isloc = NULL;
1646: }
1647: vs->islocal.row[i] = isloc;
1648: offset += nlocal;
1649: }
1650: for (i = 0, offset = 0; i < vs->nc; i++) {
1651: IS isloc;
1652: ISLocalToGlobalMapping cmap = NULL;
1653: PetscInt nlocal, bs;
1654: PetscCall(MatNestFindNonzeroSubMatCol(A, i, &sub));
1655: if (sub) PetscCall(MatGetLocalToGlobalMapping(sub, NULL, &cmap));
1656: if (cmap) {
1657: PetscCall(MatGetBlockSizes(sub, NULL, &bs));
1658: PetscCall(ISLocalToGlobalMappingGetSize(cmap, &nlocal));
1659: PetscCall(ISCreateStride(PETSC_COMM_SELF, nlocal, offset, 1, &isloc));
1660: PetscCall(ISSetBlockSize(isloc, bs));
1661: } else {
1662: nlocal = 0;
1663: isloc = NULL;
1664: }
1665: vs->islocal.col[i] = isloc;
1666: offset += nlocal;
1667: }
1669: /* Set up the aggregate ISLocalToGlobalMapping */
1670: {
1671: ISLocalToGlobalMapping rmap, cmap;
1672: PetscCall(MatNestCreateAggregateL2G_Private(A, vs->nr, vs->islocal.row, vs->isglobal.row, PETSC_FALSE, &rmap));
1673: PetscCall(MatNestCreateAggregateL2G_Private(A, vs->nc, vs->islocal.col, vs->isglobal.col, PETSC_TRUE, &cmap));
1674: if (rmap && cmap) PetscCall(MatSetLocalToGlobalMapping(A, rmap, cmap));
1675: PetscCall(ISLocalToGlobalMappingDestroy(&rmap));
1676: PetscCall(ISLocalToGlobalMappingDestroy(&cmap));
1677: }
1679: if (PetscDefined(USE_DEBUG)) {
1680: for (i = 0; i < vs->nr; i++) {
1681: for (j = 0; j < vs->nc; j++) {
1682: PetscInt m, n, M, N, mi, ni, Mi, Ni;
1683: Mat B = vs->m[i][j];
1684: if (!B) continue;
1685: PetscCall(MatGetSize(B, &M, &N));
1686: PetscCall(MatGetLocalSize(B, &m, &n));
1687: PetscCall(ISGetSize(vs->isglobal.row[i], &Mi));
1688: PetscCall(ISGetSize(vs->isglobal.col[j], &Ni));
1689: PetscCall(ISGetLocalSize(vs->isglobal.row[i], &mi));
1690: PetscCall(ISGetLocalSize(vs->isglobal.col[j], &ni));
1691: 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);
1692: 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);
1693: }
1694: }
1695: }
1697: /* Set A->assembled if all non-null blocks are currently assembled */
1698: for (i = 0; i < vs->nr; i++) {
1699: for (j = 0; j < vs->nc; j++) {
1700: if (vs->m[i][j] && !vs->m[i][j]->assembled) PetscFunctionReturn(PETSC_SUCCESS);
1701: }
1702: }
1703: A->assembled = PETSC_TRUE;
1704: PetscFunctionReturn(PETSC_SUCCESS);
1705: }
1707: /*@C
1708: MatCreateNest - Creates a new `MATNEST` matrix containing several nested submatrices, each stored separately
1710: Collective
1712: Input Parameters:
1713: + comm - Communicator for the new `MATNEST`
1714: . nr - number of nested row blocks
1715: . is_row - index sets for each nested row block, or `NULL` to make contiguous
1716: . nc - number of nested column blocks
1717: . is_col - index sets for each nested column block, or `NULL` to make contiguous
1718: - a - array of nr*nc submatrices, empty submatrices can be passed using `NULL`
1720: Output Parameter:
1721: . B - new matrix
1723: Note:
1724: In both C and Fortran, `a` must be a row-major order array holding references to the matrices.
1725: For instance, to represent the matrix
1726: $\begin{bmatrix} A_{11} & A_{12} \\ A_{21} & A_{22}\end{bmatrix}$
1727: one should use `Mat a[4]={A11,A12,A21,A22}`.
1729: Level: advanced
1731: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatCreate()`, `VecCreateNest()`, `DMCreateMatrix()`, `MatNestSetSubMat()`,
1732: `MatNestGetSubMat()`, `MatNestGetLocalISs()`, `MatNestGetSize()`,
1733: `MatNestGetISs()`, `MatNestSetSubMats()`, `MatNestGetSubMats()`
1734: @*/
1735: PetscErrorCode MatCreateNest(MPI_Comm comm, PetscInt nr, const IS is_row[], PetscInt nc, const IS is_col[], const Mat a[], Mat *B)
1736: {
1737: Mat A;
1739: PetscFunctionBegin;
1740: *B = NULL;
1741: PetscCall(MatCreate(comm, &A));
1742: PetscCall(MatSetType(A, MATNEST));
1743: A->preallocated = PETSC_TRUE;
1744: PetscCall(MatNestSetSubMats(A, nr, is_row, nc, is_col, a));
1745: *B = A;
1746: PetscFunctionReturn(PETSC_SUCCESS);
1747: }
1749: static PetscErrorCode MatConvert_Nest_SeqAIJ_fast(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
1750: {
1751: Mat_Nest *nest = (Mat_Nest *)A->data;
1752: Mat *trans;
1753: PetscScalar **avv;
1754: PetscScalar *vv;
1755: PetscInt **aii, **ajj;
1756: PetscInt *ii, *jj, *ci;
1757: PetscInt nr, nc, nnz, i, j;
1758: PetscBool done;
1760: PetscFunctionBegin;
1761: PetscCall(MatGetSize(A, &nr, &nc));
1762: if (reuse == MAT_REUSE_MATRIX) {
1763: PetscInt rnr;
1765: PetscCall(MatGetRowIJ(*newmat, 0, PETSC_FALSE, PETSC_FALSE, &rnr, (const PetscInt **)&ii, (const PetscInt **)&jj, &done));
1766: PetscCheck(done, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "MatGetRowIJ");
1767: PetscCheck(rnr == nr, PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "Cannot reuse matrix, wrong number of rows");
1768: PetscCall(MatSeqAIJGetArray(*newmat, &vv));
1769: }
1770: /* extract CSR for nested SeqAIJ matrices */
1771: nnz = 0;
1772: PetscCall(PetscCalloc4(nest->nr * nest->nc, &aii, nest->nr * nest->nc, &ajj, nest->nr * nest->nc, &avv, nest->nr * nest->nc, &trans));
1773: for (i = 0; i < nest->nr; ++i) {
1774: for (j = 0; j < nest->nc; ++j) {
1775: Mat B = nest->m[i][j];
1776: if (B) {
1777: PetscScalar *naa;
1778: PetscInt *nii, *njj, nnr;
1779: PetscBool istrans;
1781: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &istrans));
1782: if (istrans) {
1783: Mat Bt;
1785: PetscCall(MatTransposeGetMat(B, &Bt));
1786: PetscCall(MatTranspose(Bt, MAT_INITIAL_MATRIX, &trans[i * nest->nc + j]));
1787: B = trans[i * nest->nc + j];
1788: } else {
1789: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATHERMITIANTRANSPOSEVIRTUAL, &istrans));
1790: if (istrans) {
1791: Mat Bt;
1793: PetscCall(MatHermitianTransposeGetMat(B, &Bt));
1794: PetscCall(MatHermitianTranspose(Bt, MAT_INITIAL_MATRIX, &trans[i * nest->nc + j]));
1795: B = trans[i * nest->nc + j];
1796: }
1797: }
1798: PetscCall(MatGetRowIJ(B, 0, PETSC_FALSE, PETSC_FALSE, &nnr, (const PetscInt **)&nii, (const PetscInt **)&njj, &done));
1799: PetscCheck(done, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "MatGetRowIJ");
1800: PetscCall(MatSeqAIJGetArray(B, &naa));
1801: nnz += nii[nnr];
1803: aii[i * nest->nc + j] = nii;
1804: ajj[i * nest->nc + j] = njj;
1805: avv[i * nest->nc + j] = naa;
1806: }
1807: }
1808: }
1809: if (reuse != MAT_REUSE_MATRIX) {
1810: PetscCall(PetscMalloc1(nr + 1, &ii));
1811: PetscCall(PetscMalloc1(nnz, &jj));
1812: PetscCall(PetscMalloc1(nnz, &vv));
1813: } else {
1814: PetscCheck(nnz == ii[nr], PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "Cannot reuse matrix, wrong number of nonzeros");
1815: }
1817: /* new row pointer */
1818: PetscCall(PetscArrayzero(ii, nr + 1));
1819: for (i = 0; i < nest->nr; ++i) {
1820: PetscInt ncr, rst;
1822: PetscCall(ISStrideGetInfo(nest->isglobal.row[i], &rst, NULL));
1823: PetscCall(ISGetLocalSize(nest->isglobal.row[i], &ncr));
1824: for (j = 0; j < nest->nc; ++j) {
1825: if (aii[i * nest->nc + j]) {
1826: PetscInt *nii = aii[i * nest->nc + j];
1827: PetscInt ir;
1829: for (ir = rst; ir < ncr + rst; ++ir) {
1830: ii[ir + 1] += nii[1] - nii[0];
1831: nii++;
1832: }
1833: }
1834: }
1835: }
1836: for (i = 0; i < nr; i++) ii[i + 1] += ii[i];
1838: /* construct CSR for the new matrix */
1839: PetscCall(PetscCalloc1(nr, &ci));
1840: for (i = 0; i < nest->nr; ++i) {
1841: PetscInt ncr, rst;
1843: PetscCall(ISStrideGetInfo(nest->isglobal.row[i], &rst, NULL));
1844: PetscCall(ISGetLocalSize(nest->isglobal.row[i], &ncr));
1845: for (j = 0; j < nest->nc; ++j) {
1846: if (aii[i * nest->nc + j]) {
1847: PetscScalar *nvv = avv[i * nest->nc + j];
1848: PetscInt *nii = aii[i * nest->nc + j];
1849: PetscInt *njj = ajj[i * nest->nc + j];
1850: PetscInt ir, cst;
1852: PetscCall(ISStrideGetInfo(nest->isglobal.col[j], &cst, NULL));
1853: for (ir = rst; ir < ncr + rst; ++ir) {
1854: PetscInt ij, rsize = nii[1] - nii[0], ist = ii[ir] + ci[ir];
1856: for (ij = 0; ij < rsize; ij++) {
1857: jj[ist + ij] = *njj + cst;
1858: vv[ist + ij] = *nvv;
1859: njj++;
1860: nvv++;
1861: }
1862: ci[ir] += rsize;
1863: nii++;
1864: }
1865: }
1866: }
1867: }
1868: PetscCall(PetscFree(ci));
1870: /* restore info */
1871: for (i = 0; i < nest->nr; ++i) {
1872: for (j = 0; j < nest->nc; ++j) {
1873: Mat B = nest->m[i][j];
1874: if (B) {
1875: PetscInt nnr = 0, k = i * nest->nc + j;
1877: B = (trans[k] ? trans[k] : B);
1878: PetscCall(MatRestoreRowIJ(B, 0, PETSC_FALSE, PETSC_FALSE, &nnr, (const PetscInt **)&aii[k], (const PetscInt **)&ajj[k], &done));
1879: PetscCheck(done, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "MatRestoreRowIJ");
1880: PetscCall(MatSeqAIJRestoreArray(B, &avv[k]));
1881: PetscCall(MatDestroy(&trans[k]));
1882: }
1883: }
1884: }
1885: PetscCall(PetscFree4(aii, ajj, avv, trans));
1887: /* finalize newmat */
1888: if (reuse == MAT_INITIAL_MATRIX) {
1889: PetscCall(MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A), nr, nc, ii, jj, vv, newmat));
1890: } else if (reuse == MAT_INPLACE_MATRIX) {
1891: Mat B;
1893: PetscCall(MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A), nr, nc, ii, jj, vv, &B));
1894: PetscCall(MatHeaderReplace(A, &B));
1895: }
1896: PetscCall(MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY));
1897: PetscCall(MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY));
1898: {
1899: Mat_SeqAIJ *a = (Mat_SeqAIJ *)((*newmat)->data);
1900: a->free_a = PETSC_TRUE;
1901: a->free_ij = PETSC_TRUE;
1902: }
1903: PetscFunctionReturn(PETSC_SUCCESS);
1904: }
1906: PETSC_INTERN PetscErrorCode MatAXPY_Dense_Nest(Mat Y, PetscScalar a, Mat X)
1907: {
1908: Mat_Nest *nest = (Mat_Nest *)X->data;
1909: PetscInt i, j, k, rstart;
1910: PetscBool flg;
1912: PetscFunctionBegin;
1913: /* Fill by row */
1914: for (j = 0; j < nest->nc; ++j) {
1915: /* Using global column indices and ISAllGather() is not scalable. */
1916: IS bNis;
1917: PetscInt bN;
1918: const PetscInt *bNindices;
1919: PetscCall(ISAllGather(nest->isglobal.col[j], &bNis));
1920: PetscCall(ISGetSize(bNis, &bN));
1921: PetscCall(ISGetIndices(bNis, &bNindices));
1922: for (i = 0; i < nest->nr; ++i) {
1923: Mat B = nest->m[i][j], D = NULL;
1924: PetscInt bm, br;
1925: const PetscInt *bmindices;
1926: if (!B) continue;
1927: PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &flg, MATTRANSPOSEVIRTUAL, MATHERMITIANTRANSPOSEVIRTUAL, ""));
1928: if (flg) {
1929: PetscTryMethod(B, "MatTransposeGetMat_C", (Mat, Mat *), (B, &D));
1930: PetscTryMethod(B, "MatHermitianTransposeGetMat_C", (Mat, Mat *), (B, &D));
1931: PetscCall(MatConvert(B, ((PetscObject)D)->type_name, MAT_INITIAL_MATRIX, &D));
1932: B = D;
1933: }
1934: PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &flg, MATSEQSBAIJ, MATMPISBAIJ, ""));
1935: if (flg) {
1936: if (D) PetscCall(MatConvert(D, MATBAIJ, MAT_INPLACE_MATRIX, &D));
1937: else PetscCall(MatConvert(B, MATBAIJ, MAT_INITIAL_MATRIX, &D));
1938: B = D;
1939: }
1940: PetscCall(ISGetLocalSize(nest->isglobal.row[i], &bm));
1941: PetscCall(ISGetIndices(nest->isglobal.row[i], &bmindices));
1942: PetscCall(MatGetOwnershipRange(B, &rstart, NULL));
1943: for (br = 0; br < bm; ++br) {
1944: PetscInt row = bmindices[br], brncols, *cols;
1945: const PetscInt *brcols;
1946: const PetscScalar *brcoldata;
1947: PetscScalar *vals = NULL;
1948: PetscCall(MatGetRow(B, br + rstart, &brncols, &brcols, &brcoldata));
1949: PetscCall(PetscMalloc1(brncols, &cols));
1950: for (k = 0; k < brncols; k++) cols[k] = bNindices[brcols[k]];
1951: /*
1952: Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match.
1953: Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES.
1954: */
1955: if (a != 1.0) {
1956: PetscCall(PetscMalloc1(brncols, &vals));
1957: for (k = 0; k < brncols; k++) vals[k] = a * brcoldata[k];
1958: PetscCall(MatSetValues(Y, 1, &row, brncols, cols, vals, ADD_VALUES));
1959: PetscCall(PetscFree(vals));
1960: } else {
1961: PetscCall(MatSetValues(Y, 1, &row, brncols, cols, brcoldata, ADD_VALUES));
1962: }
1963: PetscCall(MatRestoreRow(B, br + rstart, &brncols, &brcols, &brcoldata));
1964: PetscCall(PetscFree(cols));
1965: }
1966: if (D) PetscCall(MatDestroy(&D));
1967: PetscCall(ISRestoreIndices(nest->isglobal.row[i], &bmindices));
1968: }
1969: PetscCall(ISRestoreIndices(bNis, &bNindices));
1970: PetscCall(ISDestroy(&bNis));
1971: }
1972: PetscCall(MatAssemblyBegin(Y, MAT_FINAL_ASSEMBLY));
1973: PetscCall(MatAssemblyEnd(Y, MAT_FINAL_ASSEMBLY));
1974: PetscFunctionReturn(PETSC_SUCCESS);
1975: }
1977: static PetscErrorCode MatConvert_Nest_AIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
1978: {
1979: Mat_Nest *nest = (Mat_Nest *)A->data;
1980: PetscInt m, n, M, N, i, j, k, *dnnz, *onnz = NULL, rstart, cstart, cend;
1981: PetscMPIInt size;
1982: Mat C;
1984: PetscFunctionBegin;
1985: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1986: if (size == 1) { /* look for a special case with SeqAIJ matrices and strided-1, contiguous, blocks */
1987: PetscInt nf;
1988: PetscBool fast;
1990: PetscCall(PetscStrcmp(newtype, MATAIJ, &fast));
1991: if (!fast) PetscCall(PetscStrcmp(newtype, MATSEQAIJ, &fast));
1992: for (i = 0; i < nest->nr && fast; ++i) {
1993: for (j = 0; j < nest->nc && fast; ++j) {
1994: Mat B = nest->m[i][j];
1995: if (B) {
1996: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATSEQAIJ, &fast));
1997: if (!fast) {
1998: PetscBool istrans;
2000: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &istrans));
2001: if (istrans) {
2002: Mat Bt;
2004: PetscCall(MatTransposeGetMat(B, &Bt));
2005: PetscCall(PetscObjectTypeCompare((PetscObject)Bt, MATSEQAIJ, &fast));
2006: } else {
2007: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATHERMITIANTRANSPOSEVIRTUAL, &istrans));
2008: if (istrans) {
2009: Mat Bt;
2011: PetscCall(MatHermitianTransposeGetMat(B, &Bt));
2012: PetscCall(PetscObjectTypeCompare((PetscObject)Bt, MATSEQAIJ, &fast));
2013: }
2014: }
2015: }
2016: }
2017: }
2018: }
2019: for (i = 0, nf = 0; i < nest->nr && fast; ++i) {
2020: PetscCall(PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i], ISSTRIDE, &fast));
2021: if (fast) {
2022: PetscInt f, s;
2024: PetscCall(ISStrideGetInfo(nest->isglobal.row[i], &f, &s));
2025: if (f != nf || s != 1) {
2026: fast = PETSC_FALSE;
2027: } else {
2028: PetscCall(ISGetSize(nest->isglobal.row[i], &f));
2029: nf += f;
2030: }
2031: }
2032: }
2033: for (i = 0, nf = 0; i < nest->nc && fast; ++i) {
2034: PetscCall(PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i], ISSTRIDE, &fast));
2035: if (fast) {
2036: PetscInt f, s;
2038: PetscCall(ISStrideGetInfo(nest->isglobal.col[i], &f, &s));
2039: if (f != nf || s != 1) {
2040: fast = PETSC_FALSE;
2041: } else {
2042: PetscCall(ISGetSize(nest->isglobal.col[i], &f));
2043: nf += f;
2044: }
2045: }
2046: }
2047: if (fast) {
2048: PetscCall(MatConvert_Nest_SeqAIJ_fast(A, newtype, reuse, newmat));
2049: PetscFunctionReturn(PETSC_SUCCESS);
2050: }
2051: }
2052: PetscCall(MatGetSize(A, &M, &N));
2053: PetscCall(MatGetLocalSize(A, &m, &n));
2054: PetscCall(MatGetOwnershipRangeColumn(A, &cstart, &cend));
2055: if (reuse == MAT_REUSE_MATRIX) C = *newmat;
2056: else {
2057: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
2058: PetscCall(MatSetType(C, newtype));
2059: PetscCall(MatSetSizes(C, m, n, M, N));
2060: }
2061: PetscCall(PetscMalloc1(2 * m, &dnnz));
2062: if (m) {
2063: onnz = dnnz + m;
2064: for (k = 0; k < m; k++) {
2065: dnnz[k] = 0;
2066: onnz[k] = 0;
2067: }
2068: }
2069: for (j = 0; j < nest->nc; ++j) {
2070: IS bNis;
2071: PetscInt bN;
2072: const PetscInt *bNindices;
2073: PetscBool flg;
2074: /* Using global column indices and ISAllGather() is not scalable. */
2075: PetscCall(ISAllGather(nest->isglobal.col[j], &bNis));
2076: PetscCall(ISGetSize(bNis, &bN));
2077: PetscCall(ISGetIndices(bNis, &bNindices));
2078: for (i = 0; i < nest->nr; ++i) {
2079: PetscSF bmsf;
2080: PetscSFNode *iremote;
2081: Mat B = nest->m[i][j], D = NULL;
2082: PetscInt bm, *sub_dnnz, *sub_onnz, br;
2083: const PetscInt *bmindices;
2084: if (!B) continue;
2085: PetscCall(ISGetLocalSize(nest->isglobal.row[i], &bm));
2086: PetscCall(ISGetIndices(nest->isglobal.row[i], &bmindices));
2087: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf));
2088: PetscCall(PetscMalloc1(bm, &iremote));
2089: PetscCall(PetscMalloc1(bm, &sub_dnnz));
2090: PetscCall(PetscMalloc1(bm, &sub_onnz));
2091: for (k = 0; k < bm; ++k) {
2092: sub_dnnz[k] = 0;
2093: sub_onnz[k] = 0;
2094: }
2095: PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &flg, MATTRANSPOSEVIRTUAL, MATHERMITIANTRANSPOSEVIRTUAL, ""));
2096: if (flg) {
2097: PetscTryMethod(B, "MatTransposeGetMat_C", (Mat, Mat *), (B, &D));
2098: PetscTryMethod(B, "MatHermitianTransposeGetMat_C", (Mat, Mat *), (B, &D));
2099: PetscCall(MatConvert(B, ((PetscObject)D)->type_name, MAT_INITIAL_MATRIX, &D));
2100: B = D;
2101: }
2102: PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &flg, MATSEQSBAIJ, MATMPISBAIJ, ""));
2103: if (flg) {
2104: if (D) PetscCall(MatConvert(D, MATBAIJ, MAT_INPLACE_MATRIX, &D));
2105: else PetscCall(MatConvert(B, MATBAIJ, MAT_INITIAL_MATRIX, &D));
2106: B = D;
2107: }
2108: /*
2109: Locate the owners for all of the locally-owned global row indices for this row block.
2110: These determine the roots of PetscSF used to communicate preallocation data to row owners.
2111: The roots correspond to the dnnz and onnz entries; thus, there are two roots per row.
2112: */
2113: PetscCall(MatGetOwnershipRange(B, &rstart, NULL));
2114: for (br = 0; br < bm; ++br) {
2115: PetscInt row = bmindices[br], brncols, col;
2116: const PetscInt *brcols;
2117: PetscInt rowrel = 0; /* row's relative index on its owner rank */
2118: PetscMPIInt rowowner = 0;
2119: PetscCall(PetscLayoutFindOwnerIndex(A->rmap, row, &rowowner, &rowrel));
2120: /* how many roots */
2121: iremote[br].rank = rowowner;
2122: iremote[br].index = rowrel; /* edge from bmdnnz to dnnz */
2123: /* get nonzero pattern */
2124: PetscCall(MatGetRow(B, br + rstart, &brncols, &brcols, NULL));
2125: for (k = 0; k < brncols; k++) {
2126: col = bNindices[brcols[k]];
2127: if (col >= A->cmap->range[rowowner] && col < A->cmap->range[rowowner + 1]) {
2128: sub_dnnz[br]++;
2129: } else {
2130: sub_onnz[br]++;
2131: }
2132: }
2133: PetscCall(MatRestoreRow(B, br + rstart, &brncols, &brcols, NULL));
2134: }
2135: if (D) PetscCall(MatDestroy(&D));
2136: PetscCall(ISRestoreIndices(nest->isglobal.row[i], &bmindices));
2137: /* bsf will have to take care of disposing of bedges. */
2138: PetscCall(PetscSFSetGraph(bmsf, m, bm, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
2139: PetscCall(PetscSFReduceBegin(bmsf, MPIU_INT, sub_dnnz, dnnz, MPI_SUM));
2140: PetscCall(PetscSFReduceEnd(bmsf, MPIU_INT, sub_dnnz, dnnz, MPI_SUM));
2141: PetscCall(PetscSFReduceBegin(bmsf, MPIU_INT, sub_onnz, onnz, MPI_SUM));
2142: PetscCall(PetscSFReduceEnd(bmsf, MPIU_INT, sub_onnz, onnz, MPI_SUM));
2143: PetscCall(PetscFree(sub_dnnz));
2144: PetscCall(PetscFree(sub_onnz));
2145: PetscCall(PetscSFDestroy(&bmsf));
2146: }
2147: PetscCall(ISRestoreIndices(bNis, &bNindices));
2148: PetscCall(ISDestroy(&bNis));
2149: }
2150: /* Resize preallocation if overestimated */
2151: for (i = 0; i < m; i++) {
2152: dnnz[i] = PetscMin(dnnz[i], A->cmap->n);
2153: onnz[i] = PetscMin(onnz[i], A->cmap->N - A->cmap->n);
2154: }
2155: PetscCall(MatSeqAIJSetPreallocation(C, 0, dnnz));
2156: PetscCall(MatMPIAIJSetPreallocation(C, 0, dnnz, 0, onnz));
2157: PetscCall(PetscFree(dnnz));
2158: PetscCall(MatAXPY_Dense_Nest(C, 1.0, A));
2159: if (reuse == MAT_INPLACE_MATRIX) {
2160: PetscCall(MatHeaderReplace(A, &C));
2161: } else *newmat = C;
2162: PetscFunctionReturn(PETSC_SUCCESS);
2163: }
2165: static PetscErrorCode MatConvert_Nest_Dense(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
2166: {
2167: Mat B;
2168: PetscInt m, n, M, N;
2170: PetscFunctionBegin;
2171: PetscCall(MatGetSize(A, &M, &N));
2172: PetscCall(MatGetLocalSize(A, &m, &n));
2173: if (reuse == MAT_REUSE_MATRIX) {
2174: B = *newmat;
2175: PetscCall(MatZeroEntries(B));
2176: } else {
2177: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A), m, PETSC_DECIDE, M, N, NULL, &B));
2178: }
2179: PetscCall(MatAXPY_Dense_Nest(B, 1.0, A));
2180: if (reuse == MAT_INPLACE_MATRIX) {
2181: PetscCall(MatHeaderReplace(A, &B));
2182: } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
2183: PetscFunctionReturn(PETSC_SUCCESS);
2184: }
2186: static PetscErrorCode MatHasOperation_Nest(Mat mat, MatOperation op, PetscBool *has)
2187: {
2188: Mat_Nest *bA = (Mat_Nest *)mat->data;
2189: MatOperation opAdd;
2190: PetscInt i, j, nr = bA->nr, nc = bA->nc;
2191: PetscBool flg;
2192: PetscFunctionBegin;
2194: *has = PETSC_FALSE;
2195: if (op == MATOP_MULT || op == MATOP_MULT_ADD || op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) {
2196: opAdd = (op == MATOP_MULT || op == MATOP_MULT_ADD ? MATOP_MULT_ADD : MATOP_MULT_TRANSPOSE_ADD);
2197: for (j = 0; j < nc; j++) {
2198: for (i = 0; i < nr; i++) {
2199: if (!bA->m[i][j]) continue;
2200: PetscCall(MatHasOperation(bA->m[i][j], opAdd, &flg));
2201: if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
2202: }
2203: }
2204: }
2205: if (((void **)mat->ops)[op]) *has = PETSC_TRUE;
2206: PetscFunctionReturn(PETSC_SUCCESS);
2207: }
2209: /*MC
2210: MATNEST - "nest" - Matrix type consisting of nested submatrices, each stored separately.
2212: Level: intermediate
2214: Notes:
2215: This matrix type permits scalable use of `PCFIELDSPLIT` and avoids the large memory costs of extracting submatrices.
2216: It allows the use of symmetric and block formats for parts of multi-physics simulations.
2217: It is usually used with `DMCOMPOSITE` and `DMCreateMatrix()`
2219: Each of the submatrices lives on the same MPI communicator as the original nest matrix (though they can have zero
2220: rows/columns on some processes.) Thus this is not meant for cases where the submatrices live on far fewer processes
2221: than the nest matrix.
2223: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatCreate()`, `MatType`, `MatCreateNest()`, `MatNestSetSubMat()`, `MatNestGetSubMat()`,
2224: `VecCreateNest()`, `DMCreateMatrix()`, `DMCOMPOSITE`, `MatNestSetVecType()`, `MatNestGetLocalISs()`,
2225: `MatNestGetISs()`, `MatNestSetSubMats()`, `MatNestGetSubMats()`
2226: M*/
2227: PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A)
2228: {
2229: Mat_Nest *s;
2231: PetscFunctionBegin;
2232: PetscCall(PetscNew(&s));
2233: A->data = (void *)s;
2235: s->nr = -1;
2236: s->nc = -1;
2237: s->m = NULL;
2238: s->splitassembly = PETSC_FALSE;
2240: PetscCall(PetscMemzero(A->ops, sizeof(*A->ops)));
2242: A->ops->mult = MatMult_Nest;
2243: A->ops->multadd = MatMultAdd_Nest;
2244: A->ops->multtranspose = MatMultTranspose_Nest;
2245: A->ops->multtransposeadd = MatMultTransposeAdd_Nest;
2246: A->ops->transpose = MatTranspose_Nest;
2247: A->ops->assemblybegin = MatAssemblyBegin_Nest;
2248: A->ops->assemblyend = MatAssemblyEnd_Nest;
2249: A->ops->zeroentries = MatZeroEntries_Nest;
2250: A->ops->copy = MatCopy_Nest;
2251: A->ops->axpy = MatAXPY_Nest;
2252: A->ops->duplicate = MatDuplicate_Nest;
2253: A->ops->createsubmatrix = MatCreateSubMatrix_Nest;
2254: A->ops->destroy = MatDestroy_Nest;
2255: A->ops->view = MatView_Nest;
2256: A->ops->getvecs = NULL; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */
2257: A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_Nest;
2258: A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest;
2259: A->ops->getdiagonal = MatGetDiagonal_Nest;
2260: A->ops->diagonalscale = MatDiagonalScale_Nest;
2261: A->ops->scale = MatScale_Nest;
2262: A->ops->shift = MatShift_Nest;
2263: A->ops->diagonalset = MatDiagonalSet_Nest;
2264: A->ops->setrandom = MatSetRandom_Nest;
2265: A->ops->hasoperation = MatHasOperation_Nest;
2266: A->ops->missingdiagonal = MatMissingDiagonal_Nest;
2268: A->spptr = NULL;
2269: A->assembled = PETSC_FALSE;
2271: /* expose Nest api's */
2272: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSubMat_C", MatNestGetSubMat_Nest));
2273: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetSubMat_C", MatNestSetSubMat_Nest));
2274: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSubMats_C", MatNestGetSubMats_Nest));
2275: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSize_C", MatNestGetSize_Nest));
2276: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetISs_C", MatNestGetISs_Nest));
2277: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetLocalISs_C", MatNestGetLocalISs_Nest));
2278: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetVecType_C", MatNestSetVecType_Nest));
2279: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetSubMats_C", MatNestSetSubMats_Nest));
2280: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_mpiaij_C", MatConvert_Nest_AIJ));
2281: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_seqaij_C", MatConvert_Nest_AIJ));
2282: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_aij_C", MatConvert_Nest_AIJ));
2283: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_is_C", MatConvert_Nest_IS));
2284: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_mpidense_C", MatConvert_Nest_Dense));
2285: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_seqdense_C", MatConvert_Nest_Dense));
2286: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_seqdense_C", MatProductSetFromOptions_Nest_Dense));
2287: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_mpidense_C", MatProductSetFromOptions_Nest_Dense));
2288: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_dense_C", MatProductSetFromOptions_Nest_Dense));
2290: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATNEST));
2291: PetscFunctionReturn(PETSC_SUCCESS);
2292: }