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: *m = *n = *M = *N = 0;
18: for (i = 0; i < bA->nr; i++) { /* rows */
19: PetscInt sm, sM;
20: ISGetLocalSize(bA->isglobal.row[i], &sm);
21: ISGetSize(bA->isglobal.row[i], &sM);
22: *m += sm;
23: *M += sM;
24: }
25: for (j = 0; j < bA->nc; j++) { /* cols */
26: PetscInt sn, sN;
27: ISGetLocalSize(bA->isglobal.col[j], &sn);
28: ISGetSize(bA->isglobal.col[j], &sN);
29: *n += sn;
30: *N += sN;
31: }
32: return 0;
33: }
35: /* operations */
36: static PetscErrorCode MatMult_Nest(Mat A, Vec x, Vec y)
37: {
38: Mat_Nest *bA = (Mat_Nest *)A->data;
39: Vec *bx = bA->right, *by = bA->left;
40: PetscInt i, j, nr = bA->nr, nc = bA->nc;
42: for (i = 0; i < nr; i++) VecGetSubVector(y, bA->isglobal.row[i], &by[i]);
43: for (i = 0; i < nc; i++) VecGetSubVector(x, bA->isglobal.col[i], &bx[i]);
44: for (i = 0; i < nr; i++) {
45: VecZeroEntries(by[i]);
46: for (j = 0; j < nc; j++) {
47: if (!bA->m[i][j]) continue;
48: /* y[i] <- y[i] + A[i][j] * x[j] */
49: MatMultAdd(bA->m[i][j], bx[j], by[i], by[i]);
50: }
51: }
52: for (i = 0; i < nr; i++) VecRestoreSubVector(y, bA->isglobal.row[i], &by[i]);
53: for (i = 0; i < nc; i++) VecRestoreSubVector(x, bA->isglobal.col[i], &bx[i]);
54: return 0;
55: }
57: static PetscErrorCode MatMultAdd_Nest(Mat A, Vec x, Vec y, Vec z)
58: {
59: Mat_Nest *bA = (Mat_Nest *)A->data;
60: Vec *bx = bA->right, *bz = bA->left;
61: PetscInt i, j, nr = bA->nr, nc = bA->nc;
63: for (i = 0; i < nr; i++) VecGetSubVector(z, bA->isglobal.row[i], &bz[i]);
64: for (i = 0; i < nc; i++) VecGetSubVector(x, bA->isglobal.col[i], &bx[i]);
65: for (i = 0; i < nr; i++) {
66: if (y != z) {
67: Vec by;
68: VecGetSubVector(y, bA->isglobal.row[i], &by);
69: VecCopy(by, bz[i]);
70: VecRestoreSubVector(y, bA->isglobal.row[i], &by);
71: }
72: for (j = 0; j < nc; j++) {
73: if (!bA->m[i][j]) continue;
74: /* y[i] <- y[i] + A[i][j] * x[j] */
75: MatMultAdd(bA->m[i][j], bx[j], bz[i], bz[i]);
76: }
77: }
78: for (i = 0; i < nr; i++) VecRestoreSubVector(z, bA->isglobal.row[i], &bz[i]);
79: for (i = 0; i < nc; i++) VecRestoreSubVector(x, bA->isglobal.col[i], &bx[i]);
80: return 0;
81: }
83: typedef struct {
84: Mat *workC; /* array of Mat with specific containers depending on the underlying MatMatMult implementation */
85: PetscScalar *tarray; /* buffer for storing all temporary products A[i][j] B[j] */
86: PetscInt *dm, *dn, k; /* displacements and number of submatrices */
87: } Nest_Dense;
89: PETSC_INTERN PetscErrorCode MatProductNumeric_Nest_Dense(Mat C)
90: {
91: Mat_Nest *bA;
92: Nest_Dense *contents;
93: Mat viewB, viewC, productB, workC;
94: const PetscScalar *barray;
95: PetscScalar *carray;
96: PetscInt i, j, M, N, nr, nc, ldb, ldc;
97: Mat A, B;
99: MatCheckProduct(C, 3);
100: A = C->product->A;
101: B = C->product->B;
102: MatGetSize(B, NULL, &N);
103: if (!N) {
104: MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
105: MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
106: return 0;
107: }
108: contents = (Nest_Dense *)C->product->data;
110: bA = (Mat_Nest *)A->data;
111: nr = bA->nr;
112: nc = bA->nc;
113: MatDenseGetLDA(B, &ldb);
114: MatDenseGetLDA(C, &ldc);
115: MatZeroEntries(C);
116: MatDenseGetArrayRead(B, &barray);
117: MatDenseGetArray(C, &carray);
118: for (i = 0; i < nr; i++) {
119: ISGetSize(bA->isglobal.row[i], &M);
120: MatCreateDense(PetscObjectComm((PetscObject)A), contents->dm[i + 1] - contents->dm[i], PETSC_DECIDE, M, N, carray + contents->dm[i], &viewC);
121: MatDenseSetLDA(viewC, ldc);
122: for (j = 0; j < nc; j++) {
123: if (!bA->m[i][j]) continue;
124: ISGetSize(bA->isglobal.col[j], &M);
125: MatCreateDense(PetscObjectComm((PetscObject)A), contents->dn[j + 1] - contents->dn[j], PETSC_DECIDE, M, N, (PetscScalar *)(barray + contents->dn[j]), &viewB);
126: MatDenseSetLDA(viewB, ldb);
128: /* MatMatMultNumeric(bA->m[i][j],viewB,contents->workC[i*nc + j]); */
129: workC = contents->workC[i * nc + j];
130: productB = workC->product->B;
131: workC->product->B = viewB; /* use newly created dense matrix viewB */
132: MatProductNumeric(workC);
133: MatDestroy(&viewB);
134: workC->product->B = productB; /* resume original B */
136: /* C[i] <- workC + C[i] */
137: MatAXPY(viewC, 1.0, contents->workC[i * nc + j], SAME_NONZERO_PATTERN);
138: }
139: MatDestroy(&viewC);
140: }
141: MatDenseRestoreArray(C, &carray);
142: MatDenseRestoreArrayRead(B, &barray);
144: MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
145: MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
146: return 0;
147: }
149: PetscErrorCode MatNest_DenseDestroy(void *ctx)
150: {
151: Nest_Dense *contents = (Nest_Dense *)ctx;
152: PetscInt i;
154: PetscFree(contents->tarray);
155: for (i = 0; i < contents->k; i++) MatDestroy(contents->workC + i);
156: PetscFree3(contents->dm, contents->dn, contents->workC);
157: PetscFree(contents);
158: return 0;
159: }
161: PETSC_INTERN PetscErrorCode MatProductSymbolic_Nest_Dense(Mat C)
162: {
163: Mat_Nest *bA;
164: Mat viewB, workC;
165: const PetscScalar *barray;
166: PetscInt i, j, M, N, m, n, nr, nc, maxm = 0, ldb;
167: Nest_Dense *contents = NULL;
168: PetscBool cisdense;
169: Mat A, B;
170: PetscReal fill;
172: MatCheckProduct(C, 4);
174: A = C->product->A;
175: B = C->product->B;
176: fill = C->product->fill;
177: bA = (Mat_Nest *)A->data;
178: nr = bA->nr;
179: nc = bA->nc;
180: MatGetLocalSize(C, &m, &n);
181: MatGetSize(C, &M, &N);
182: if (m == PETSC_DECIDE || n == PETSC_DECIDE || M == PETSC_DECIDE || N == PETSC_DECIDE) {
183: MatGetLocalSize(B, NULL, &n);
184: MatGetSize(B, NULL, &N);
185: MatGetLocalSize(A, &m, NULL);
186: MatGetSize(A, &M, NULL);
187: MatSetSizes(C, m, n, M, N);
188: }
189: PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATMPIDENSE, MATSEQDENSECUDA, MATMPIDENSECUDA, "");
190: if (!cisdense) MatSetType(C, ((PetscObject)B)->type_name);
191: MatSetUp(C);
192: if (!N) {
193: C->ops->productnumeric = MatProductNumeric_Nest_Dense;
194: return 0;
195: }
197: PetscNew(&contents);
198: C->product->data = contents;
199: C->product->destroy = MatNest_DenseDestroy;
200: PetscCalloc3(nr + 1, &contents->dm, nc + 1, &contents->dn, nr * nc, &contents->workC);
201: contents->k = nr * nc;
202: for (i = 0; i < nr; i++) {
203: ISGetLocalSize(bA->isglobal.row[i], contents->dm + i + 1);
204: maxm = PetscMax(maxm, contents->dm[i + 1]);
205: contents->dm[i + 1] += contents->dm[i];
206: }
207: for (i = 0; i < nc; i++) {
208: ISGetLocalSize(bA->isglobal.col[i], contents->dn + i + 1);
209: contents->dn[i + 1] += contents->dn[i];
210: }
211: PetscMalloc1(maxm * N, &contents->tarray);
212: MatDenseGetLDA(B, &ldb);
213: MatGetSize(B, NULL, &N);
214: MatDenseGetArrayRead(B, &barray);
215: /* loops are permuted compared to MatMatMultNumeric so that viewB is created only once per column of A */
216: for (j = 0; j < nc; j++) {
217: ISGetSize(bA->isglobal.col[j], &M);
218: MatCreateDense(PetscObjectComm((PetscObject)A), contents->dn[j + 1] - contents->dn[j], PETSC_DECIDE, M, N, (PetscScalar *)(barray + contents->dn[j]), &viewB);
219: MatDenseSetLDA(viewB, ldb);
220: for (i = 0; i < nr; i++) {
221: if (!bA->m[i][j]) continue;
222: /* MatMatMultSymbolic may attach a specific container (depending on MatType of bA->m[i][j]) to workC[i][j] */
224: MatProductCreate(bA->m[i][j], viewB, NULL, &contents->workC[i * nc + j]);
225: workC = contents->workC[i * nc + j];
226: MatProductSetType(workC, MATPRODUCT_AB);
227: MatProductSetAlgorithm(workC, "default");
228: MatProductSetFill(workC, fill);
229: MatProductSetFromOptions(workC);
230: MatProductSymbolic(workC);
232: /* since tarray will be shared by all Mat */
233: MatSeqDenseSetPreallocation(workC, contents->tarray);
234: MatMPIDenseSetPreallocation(workC, contents->tarray);
235: }
236: MatDestroy(&viewB);
237: }
238: MatDenseRestoreArrayRead(B, &barray);
240: C->ops->productnumeric = MatProductNumeric_Nest_Dense;
241: return 0;
242: }
244: /* --------------------------------------------------------- */
245: static PetscErrorCode MatProductSetFromOptions_Nest_Dense_AB(Mat C)
246: {
247: C->ops->productsymbolic = MatProductSymbolic_Nest_Dense;
248: return 0;
249: }
251: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Nest_Dense(Mat C)
252: {
253: Mat_Product *product = C->product;
255: if (product->type == MATPRODUCT_AB) MatProductSetFromOptions_Nest_Dense_AB(C);
256: return 0;
257: }
258: /* --------------------------------------------------------- */
260: static PetscErrorCode MatMultTranspose_Nest(Mat A, Vec x, Vec y)
261: {
262: Mat_Nest *bA = (Mat_Nest *)A->data;
263: Vec *bx = bA->left, *by = bA->right;
264: PetscInt i, j, nr = bA->nr, nc = bA->nc;
266: for (i = 0; i < nr; i++) VecGetSubVector(x, bA->isglobal.row[i], &bx[i]);
267: for (i = 0; i < nc; i++) VecGetSubVector(y, bA->isglobal.col[i], &by[i]);
268: for (j = 0; j < nc; j++) {
269: 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: MatMultTransposeAdd(bA->m[i][j], bx[i], by[j], by[j]);
274: }
275: }
276: for (i = 0; i < nr; i++) VecRestoreSubVector(x, bA->isglobal.row[i], &bx[i]);
277: for (i = 0; i < nc; i++) VecRestoreSubVector(y, bA->isglobal.col[i], &by[i]);
278: return 0;
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: for (i = 0; i < nr; i++) VecGetSubVector(x, bA->isglobal.row[i], &bx[i]);
288: for (i = 0; i < nc; i++) VecGetSubVector(z, bA->isglobal.col[i], &bz[i]);
289: for (j = 0; j < nc; j++) {
290: if (y != z) {
291: Vec by;
292: VecGetSubVector(y, bA->isglobal.col[j], &by);
293: VecCopy(by, bz[j]);
294: VecRestoreSubVector(y, bA->isglobal.col[j], &by);
295: }
296: for (i = 0; i < nr; i++) {
297: if (!bA->m[i][j]) continue;
298: /* z[j] <- y[j] + (A[i][j])^T * x[i] */
299: MatMultTransposeAdd(bA->m[i][j], bx[i], bz[j], bz[j]);
300: }
301: }
302: for (i = 0; i < nr; i++) VecRestoreSubVector(x, bA->isglobal.row[i], &bx[i]);
303: for (i = 0; i < nc; i++) VecRestoreSubVector(z, bA->isglobal.col[i], &bz[i]);
304: return 0;
305: }
307: static PetscErrorCode MatTranspose_Nest(Mat A, MatReuse reuse, Mat *B)
308: {
309: Mat_Nest *bA = (Mat_Nest *)A->data, *bC;
310: Mat C;
311: PetscInt i, j, nr = bA->nr, nc = bA->nc;
313: if (reuse == MAT_REUSE_MATRIX) MatTransposeCheckNonzeroState_Private(A, *B);
316: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
317: Mat *subs;
318: IS *is_row, *is_col;
320: PetscCalloc1(nr * nc, &subs);
321: PetscMalloc2(nr, &is_row, nc, &is_col);
322: MatNestGetISs(A, is_row, is_col);
323: if (reuse == MAT_INPLACE_MATRIX) {
324: for (i = 0; i < nr; i++) {
325: for (j = 0; j < nc; j++) subs[i + nr * j] = bA->m[i][j];
326: }
327: }
329: MatCreateNest(PetscObjectComm((PetscObject)A), nc, is_col, nr, is_row, subs, &C);
330: PetscFree(subs);
331: PetscFree2(is_row, is_col);
332: } else {
333: C = *B;
334: }
336: bC = (Mat_Nest *)C->data;
337: for (i = 0; i < nr; i++) {
338: for (j = 0; j < nc; j++) {
339: if (bA->m[i][j]) {
340: MatTranspose(bA->m[i][j], reuse, &(bC->m[j][i]));
341: } else {
342: bC->m[j][i] = NULL;
343: }
344: }
345: }
347: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
348: *B = C;
349: } else {
350: MatHeaderMerge(A, &C);
351: }
352: return 0;
353: }
355: static PetscErrorCode MatNestDestroyISList(PetscInt n, IS **list)
356: {
357: IS *lst = *list;
358: PetscInt i;
360: if (!lst) return 0;
361: for (i = 0; i < n; i++)
362: if (lst[i]) ISDestroy(&lst[i]);
363: PetscFree(lst);
364: *list = NULL;
365: return 0;
366: }
368: static PetscErrorCode MatReset_Nest(Mat A)
369: {
370: Mat_Nest *vs = (Mat_Nest *)A->data;
371: PetscInt i, j;
373: /* release the matrices and the place holders */
374: MatNestDestroyISList(vs->nr, &vs->isglobal.row);
375: MatNestDestroyISList(vs->nc, &vs->isglobal.col);
376: MatNestDestroyISList(vs->nr, &vs->islocal.row);
377: MatNestDestroyISList(vs->nc, &vs->islocal.col);
379: PetscFree(vs->row_len);
380: PetscFree(vs->col_len);
381: PetscFree(vs->nnzstate);
383: PetscFree2(vs->left, vs->right);
385: /* release the matrices and the place holders */
386: if (vs->m) {
387: for (i = 0; i < vs->nr; i++) {
388: for (j = 0; j < vs->nc; j++) MatDestroy(&vs->m[i][j]);
389: PetscFree(vs->m[i]);
390: }
391: PetscFree(vs->m);
392: }
394: /* restore defaults */
395: vs->nr = 0;
396: vs->nc = 0;
397: vs->splitassembly = PETSC_FALSE;
398: return 0;
399: }
401: static PetscErrorCode MatDestroy_Nest(Mat A)
402: {
403: MatReset_Nest(A);
404: PetscFree(A->data);
405: PetscObjectComposeFunction((PetscObject)A, "MatNestGetSubMat_C", NULL);
406: PetscObjectComposeFunction((PetscObject)A, "MatNestSetSubMat_C", NULL);
407: PetscObjectComposeFunction((PetscObject)A, "MatNestGetSubMats_C", NULL);
408: PetscObjectComposeFunction((PetscObject)A, "MatNestGetSize_C", NULL);
409: PetscObjectComposeFunction((PetscObject)A, "MatNestGetISs_C", NULL);
410: PetscObjectComposeFunction((PetscObject)A, "MatNestGetLocalISs_C", NULL);
411: PetscObjectComposeFunction((PetscObject)A, "MatNestSetVecType_C", NULL);
412: PetscObjectComposeFunction((PetscObject)A, "MatNestSetSubMats_C", NULL);
413: PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_mpiaij_C", NULL);
414: PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_seqaij_C", NULL);
415: PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_aij_C", NULL);
416: PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_is_C", NULL);
417: PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_mpidense_C", NULL);
418: PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_seqdense_C", NULL);
419: PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_seqdense_C", NULL);
420: PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_mpidense_C", NULL);
421: PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_dense_C", NULL);
422: return 0;
423: }
425: static PetscErrorCode MatMissingDiagonal_Nest(Mat mat, PetscBool *missing, PetscInt *dd)
426: {
427: Mat_Nest *vs = (Mat_Nest *)mat->data;
428: PetscInt i;
430: if (dd) *dd = 0;
431: if (!vs->nr) {
432: *missing = PETSC_TRUE;
433: return 0;
434: }
435: *missing = PETSC_FALSE;
436: for (i = 0; i < vs->nr && !(*missing); i++) {
437: *missing = PETSC_TRUE;
438: if (vs->m[i][i]) {
439: MatMissingDiagonal(vs->m[i][i], missing, NULL);
441: }
442: }
443: return 0;
444: }
446: static PetscErrorCode MatAssemblyBegin_Nest(Mat A, MatAssemblyType type)
447: {
448: Mat_Nest *vs = (Mat_Nest *)A->data;
449: PetscInt i, j;
450: PetscBool nnzstate = PETSC_FALSE;
452: for (i = 0; i < vs->nr; i++) {
453: for (j = 0; j < vs->nc; j++) {
454: PetscObjectState subnnzstate = 0;
455: if (vs->m[i][j]) {
456: MatAssemblyBegin(vs->m[i][j], type);
457: if (!vs->splitassembly) {
458: /* Note: split assembly will fail if the same block appears more than once (even indirectly through a nested
459: * sub-block). This could be fixed by adding a flag to Mat so that there was a way to check if a Mat was
460: * already performing an assembly, but the result would by more complicated and appears to offer less
461: * potential for diagnostics and correctness checking. Split assembly should be fixed once there is an
462: * interface for libraries to make asynchronous progress in "user-defined non-blocking collectives".
463: */
464: MatAssemblyEnd(vs->m[i][j], type);
465: MatGetNonzeroState(vs->m[i][j], &subnnzstate);
466: }
467: }
468: nnzstate = (PetscBool)(nnzstate || vs->nnzstate[i * vs->nc + j] != subnnzstate);
469: vs->nnzstate[i * vs->nc + j] = subnnzstate;
470: }
471: }
472: if (nnzstate) A->nonzerostate++;
473: return 0;
474: }
476: static PetscErrorCode MatAssemblyEnd_Nest(Mat A, MatAssemblyType type)
477: {
478: Mat_Nest *vs = (Mat_Nest *)A->data;
479: PetscInt i, j;
481: for (i = 0; i < vs->nr; i++) {
482: for (j = 0; j < vs->nc; j++) {
483: if (vs->m[i][j]) {
484: if (vs->splitassembly) MatAssemblyEnd(vs->m[i][j], type);
485: }
486: }
487: }
488: return 0;
489: }
491: static PetscErrorCode MatNestFindNonzeroSubMatRow(Mat A, PetscInt row, Mat *B)
492: {
493: Mat_Nest *vs = (Mat_Nest *)A->data;
494: PetscInt j;
495: Mat sub;
497: sub = (row < vs->nc) ? vs->m[row][row] : (Mat)NULL; /* Prefer to find on the diagonal */
498: for (j = 0; !sub && j < vs->nc; j++) sub = vs->m[row][j];
499: if (sub) MatSetUp(sub); /* Ensure that the sizes are available */
500: *B = sub;
501: return 0;
502: }
504: static PetscErrorCode MatNestFindNonzeroSubMatCol(Mat A, PetscInt col, Mat *B)
505: {
506: Mat_Nest *vs = (Mat_Nest *)A->data;
507: PetscInt i;
508: Mat sub;
510: sub = (col < vs->nr) ? vs->m[col][col] : (Mat)NULL; /* Prefer to find on the diagonal */
511: for (i = 0; !sub && i < vs->nr; i++) sub = vs->m[i][col];
512: if (sub) MatSetUp(sub); /* Ensure that the sizes are available */
513: *B = sub;
514: return 0;
515: }
517: static PetscErrorCode MatNestFindISRange(Mat A, PetscInt n, const IS list[], IS is, PetscInt *begin, PetscInt *end)
518: {
519: PetscInt i, j, size, m;
520: PetscBool flg;
521: IS out, concatenate[2];
525: if (begin) {
527: *begin = -1;
528: }
529: if (end) {
531: *end = -1;
532: }
533: for (i = 0; i < n; i++) {
534: if (!list[i]) continue;
535: ISEqualUnsorted(list[i], is, &flg);
536: if (flg) {
537: if (begin) *begin = i;
538: if (end) *end = i + 1;
539: return 0;
540: }
541: }
542: ISGetSize(is, &size);
543: for (i = 0; i < n - 1; i++) {
544: if (!list[i]) continue;
545: m = 0;
546: ISConcatenate(PetscObjectComm((PetscObject)A), 2, list + i, &out);
547: ISGetSize(out, &m);
548: for (j = i + 2; j < n && m < size; j++) {
549: if (list[j]) {
550: concatenate[0] = out;
551: concatenate[1] = list[j];
552: ISConcatenate(PetscObjectComm((PetscObject)A), 2, concatenate, &out);
553: ISDestroy(concatenate);
554: ISGetSize(out, &m);
555: }
556: }
557: if (m == size) {
558: ISEqualUnsorted(out, is, &flg);
559: if (flg) {
560: if (begin) *begin = i;
561: if (end) *end = j;
562: ISDestroy(&out);
563: return 0;
564: }
565: }
566: ISDestroy(&out);
567: }
568: return 0;
569: }
571: static PetscErrorCode MatNestFillEmptyMat_Private(Mat A, PetscInt i, PetscInt j, Mat *B)
572: {
573: Mat_Nest *vs = (Mat_Nest *)A->data;
574: PetscInt lr, lc;
576: MatCreate(PetscObjectComm((PetscObject)A), B);
577: ISGetLocalSize(vs->isglobal.row[i], &lr);
578: ISGetLocalSize(vs->isglobal.col[j], &lc);
579: MatSetSizes(*B, lr, lc, PETSC_DECIDE, PETSC_DECIDE);
580: MatSetType(*B, MATAIJ);
581: MatSeqAIJSetPreallocation(*B, 0, NULL);
582: MatMPIAIJSetPreallocation(*B, 0, NULL, 0, NULL);
583: MatSetUp(*B);
584: MatSetOption(*B, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE);
585: MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY);
586: MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY);
587: return 0;
588: }
590: static PetscErrorCode MatNestGetBlock_Private(Mat A, PetscInt rbegin, PetscInt rend, PetscInt cbegin, PetscInt cend, Mat *B)
591: {
592: Mat_Nest *vs = (Mat_Nest *)A->data;
593: Mat *a;
594: PetscInt i, j, k, l, nr = rend - rbegin, nc = cend - cbegin;
595: char keyname[256];
596: PetscBool *b;
597: PetscBool flg;
599: *B = NULL;
600: PetscSNPrintf(keyname, sizeof(keyname), "NestBlock_%" PetscInt_FMT "-%" PetscInt_FMT "x%" PetscInt_FMT "-%" PetscInt_FMT, rbegin, rend, cbegin, cend);
601: PetscObjectQuery((PetscObject)A, keyname, (PetscObject *)B);
602: if (*B) return 0;
604: PetscMalloc2(nr * nc, &a, nr * nc, &b);
605: for (i = 0; i < nr; i++) {
606: for (j = 0; j < nc; j++) {
607: a[i * nc + j] = vs->m[rbegin + i][cbegin + j];
608: b[i * nc + j] = PETSC_FALSE;
609: }
610: }
611: if (nc != vs->nc && nr != vs->nr) {
612: for (i = 0; i < nr; i++) {
613: for (j = 0; j < nc; j++) {
614: flg = PETSC_FALSE;
615: for (k = 0; (k < nr && !flg); k++) {
616: if (a[j + k * nc]) flg = PETSC_TRUE;
617: }
618: if (flg) {
619: flg = PETSC_FALSE;
620: for (l = 0; (l < nc && !flg); l++) {
621: if (a[i * nc + l]) flg = PETSC_TRUE;
622: }
623: }
624: if (!flg) {
625: b[i * nc + j] = PETSC_TRUE;
626: MatNestFillEmptyMat_Private(A, rbegin + i, cbegin + j, a + i * nc + j);
627: }
628: }
629: }
630: }
631: MatCreateNest(PetscObjectComm((PetscObject)A), nr, nr != vs->nr ? NULL : vs->isglobal.row, nc, nc != vs->nc ? NULL : vs->isglobal.col, a, B);
632: for (i = 0; i < nr; i++) {
633: for (j = 0; j < nc; j++) {
634: if (b[i * nc + j]) MatDestroy(a + i * nc + j);
635: }
636: }
637: PetscFree2(a, b);
638: (*B)->assembled = A->assembled;
639: PetscObjectCompose((PetscObject)A, keyname, (PetscObject)*B);
640: PetscObjectDereference((PetscObject)*B); /* Leave the only remaining reference in the composition */
641: return 0;
642: }
644: static PetscErrorCode MatNestFindSubMat(Mat A, struct MatNestISPair *is, IS isrow, IS iscol, Mat *B)
645: {
646: Mat_Nest *vs = (Mat_Nest *)A->data;
647: PetscInt rbegin, rend, cbegin, cend;
649: MatNestFindISRange(A, vs->nr, is->row, isrow, &rbegin, &rend);
650: MatNestFindISRange(A, vs->nc, is->col, iscol, &cbegin, &cend);
651: if (rend == rbegin + 1 && cend == cbegin + 1) {
652: if (!vs->m[rbegin][cbegin]) MatNestFillEmptyMat_Private(A, rbegin, cbegin, vs->m[rbegin] + cbegin);
653: *B = vs->m[rbegin][cbegin];
654: } else if (rbegin != -1 && cbegin != -1) {
655: MatNestGetBlock_Private(A, rbegin, rend, cbegin, cend, B);
656: } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Could not find index set");
657: return 0;
658: }
660: /*
661: TODO: This does not actually returns a submatrix we can modify
662: */
663: static PetscErrorCode MatCreateSubMatrix_Nest(Mat A, IS isrow, IS iscol, MatReuse reuse, Mat *B)
664: {
665: Mat_Nest *vs = (Mat_Nest *)A->data;
666: Mat sub;
668: MatNestFindSubMat(A, &vs->isglobal, isrow, iscol, &sub);
669: switch (reuse) {
670: case MAT_INITIAL_MATRIX:
671: if (sub) PetscObjectReference((PetscObject)sub);
672: *B = sub;
673: break;
674: case MAT_REUSE_MATRIX:
676: break;
677: case MAT_IGNORE_MATRIX: /* Nothing to do */
678: break;
679: case MAT_INPLACE_MATRIX: /* Nothing to do */
680: SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MAT_INPLACE_MATRIX is not supported yet");
681: }
682: return 0;
683: }
685: PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A, IS isrow, IS iscol, Mat *B)
686: {
687: Mat_Nest *vs = (Mat_Nest *)A->data;
688: Mat sub;
690: MatNestFindSubMat(A, &vs->islocal, isrow, iscol, &sub);
691: /* We allow the submatrix to be NULL, perhaps it would be better for the user to return an empty matrix instead */
692: if (sub) PetscObjectReference((PetscObject)sub);
693: *B = sub;
694: return 0;
695: }
697: static PetscErrorCode MatRestoreLocalSubMatrix_Nest(Mat A, IS isrow, IS iscol, Mat *B)
698: {
699: Mat_Nest *vs = (Mat_Nest *)A->data;
700: Mat sub;
702: MatNestFindSubMat(A, &vs->islocal, isrow, iscol, &sub);
704: if (sub) {
706: MatDestroy(B);
707: }
708: return 0;
709: }
711: static PetscErrorCode MatGetDiagonal_Nest(Mat A, Vec v)
712: {
713: Mat_Nest *bA = (Mat_Nest *)A->data;
714: PetscInt i;
716: for (i = 0; i < bA->nr; i++) {
717: Vec bv;
718: VecGetSubVector(v, bA->isglobal.row[i], &bv);
719: if (bA->m[i][i]) {
720: MatGetDiagonal(bA->m[i][i], bv);
721: } else {
722: VecSet(bv, 0.0);
723: }
724: VecRestoreSubVector(v, bA->isglobal.row[i], &bv);
725: }
726: return 0;
727: }
729: static PetscErrorCode MatDiagonalScale_Nest(Mat A, Vec l, Vec r)
730: {
731: Mat_Nest *bA = (Mat_Nest *)A->data;
732: Vec bl, *br;
733: PetscInt i, j;
735: PetscCalloc1(bA->nc, &br);
736: if (r) {
737: for (j = 0; j < bA->nc; j++) VecGetSubVector(r, bA->isglobal.col[j], &br[j]);
738: }
739: bl = NULL;
740: for (i = 0; i < bA->nr; i++) {
741: if (l) VecGetSubVector(l, bA->isglobal.row[i], &bl);
742: for (j = 0; j < bA->nc; j++) {
743: if (bA->m[i][j]) MatDiagonalScale(bA->m[i][j], bl, br[j]);
744: }
745: if (l) VecRestoreSubVector(l, bA->isglobal.row[i], &bl);
746: }
747: if (r) {
748: for (j = 0; j < bA->nc; j++) VecRestoreSubVector(r, bA->isglobal.col[j], &br[j]);
749: }
750: PetscFree(br);
751: return 0;
752: }
754: static PetscErrorCode MatScale_Nest(Mat A, PetscScalar a)
755: {
756: Mat_Nest *bA = (Mat_Nest *)A->data;
757: PetscInt i, j;
759: for (i = 0; i < bA->nr; i++) {
760: for (j = 0; j < bA->nc; j++) {
761: if (bA->m[i][j]) MatScale(bA->m[i][j], a);
762: }
763: }
764: return 0;
765: }
767: static PetscErrorCode MatShift_Nest(Mat A, PetscScalar a)
768: {
769: Mat_Nest *bA = (Mat_Nest *)A->data;
770: PetscInt i;
771: PetscBool nnzstate = PETSC_FALSE;
773: for (i = 0; i < bA->nr; i++) {
774: PetscObjectState subnnzstate = 0;
776: MatShift(bA->m[i][i], a);
777: MatGetNonzeroState(bA->m[i][i], &subnnzstate);
778: nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i * bA->nc + i] != subnnzstate);
779: bA->nnzstate[i * bA->nc + i] = subnnzstate;
780: }
781: if (nnzstate) A->nonzerostate++;
782: return 0;
783: }
785: static PetscErrorCode MatDiagonalSet_Nest(Mat A, Vec D, InsertMode is)
786: {
787: Mat_Nest *bA = (Mat_Nest *)A->data;
788: PetscInt i;
789: PetscBool nnzstate = PETSC_FALSE;
791: for (i = 0; i < bA->nr; i++) {
792: PetscObjectState subnnzstate = 0;
793: Vec bv;
794: VecGetSubVector(D, bA->isglobal.row[i], &bv);
795: if (bA->m[i][i]) {
796: MatDiagonalSet(bA->m[i][i], bv, is);
797: MatGetNonzeroState(bA->m[i][i], &subnnzstate);
798: }
799: VecRestoreSubVector(D, bA->isglobal.row[i], &bv);
800: nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i * bA->nc + i] != subnnzstate);
801: bA->nnzstate[i * bA->nc + i] = subnnzstate;
802: }
803: if (nnzstate) A->nonzerostate++;
804: return 0;
805: }
807: static PetscErrorCode MatSetRandom_Nest(Mat A, PetscRandom rctx)
808: {
809: Mat_Nest *bA = (Mat_Nest *)A->data;
810: PetscInt i, j;
812: for (i = 0; i < bA->nr; i++) {
813: for (j = 0; j < bA->nc; j++) {
814: if (bA->m[i][j]) MatSetRandom(bA->m[i][j], rctx);
815: }
816: }
817: return 0;
818: }
820: static PetscErrorCode MatCreateVecs_Nest(Mat A, Vec *right, Vec *left)
821: {
822: Mat_Nest *bA = (Mat_Nest *)A->data;
823: Vec *L, *R;
824: MPI_Comm comm;
825: PetscInt i, j;
827: PetscObjectGetComm((PetscObject)A, &comm);
828: if (right) {
829: /* allocate R */
830: PetscMalloc1(bA->nc, &R);
831: /* Create the right vectors */
832: for (j = 0; j < bA->nc; j++) {
833: for (i = 0; i < bA->nr; i++) {
834: if (bA->m[i][j]) {
835: MatCreateVecs(bA->m[i][j], &R[j], NULL);
836: break;
837: }
838: }
840: }
841: VecCreateNest(comm, bA->nc, bA->isglobal.col, R, right);
842: /* hand back control to the nest vector */
843: for (j = 0; j < bA->nc; j++) VecDestroy(&R[j]);
844: PetscFree(R);
845: }
847: if (left) {
848: /* allocate L */
849: PetscMalloc1(bA->nr, &L);
850: /* Create the left vectors */
851: for (i = 0; i < bA->nr; i++) {
852: for (j = 0; j < bA->nc; j++) {
853: if (bA->m[i][j]) {
854: MatCreateVecs(bA->m[i][j], NULL, &L[i]);
855: break;
856: }
857: }
859: }
861: VecCreateNest(comm, bA->nr, bA->isglobal.row, L, left);
862: for (i = 0; i < bA->nr; i++) VecDestroy(&L[i]);
864: PetscFree(L);
865: }
866: return 0;
867: }
869: static PetscErrorCode MatView_Nest(Mat A, PetscViewer viewer)
870: {
871: Mat_Nest *bA = (Mat_Nest *)A->data;
872: PetscBool isascii, viewSub = PETSC_FALSE;
873: PetscInt i, j;
875: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
876: if (isascii) {
877: PetscOptionsGetBool(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_view_nest_sub", &viewSub, NULL);
878: PetscViewerASCIIPrintf(viewer, "Matrix object: \n");
879: PetscViewerASCIIPushTab(viewer);
880: PetscViewerASCIIPrintf(viewer, "type=nest, rows=%" PetscInt_FMT ", cols=%" PetscInt_FMT " \n", bA->nr, bA->nc);
882: PetscViewerASCIIPrintf(viewer, "MatNest structure: \n");
883: for (i = 0; i < bA->nr; i++) {
884: for (j = 0; j < bA->nc; j++) {
885: MatType type;
886: char name[256] = "", prefix[256] = "";
887: PetscInt NR, NC;
888: PetscBool isNest = PETSC_FALSE;
890: if (!bA->m[i][j]) {
891: PetscViewerASCIIPrintf(viewer, "(%" PetscInt_FMT ",%" PetscInt_FMT ") : NULL \n", i, j);
892: continue;
893: }
894: MatGetSize(bA->m[i][j], &NR, &NC);
895: MatGetType(bA->m[i][j], &type);
896: if (((PetscObject)bA->m[i][j])->name) PetscSNPrintf(name, sizeof(name), "name=\"%s\", ", ((PetscObject)bA->m[i][j])->name);
897: if (((PetscObject)bA->m[i][j])->prefix) PetscSNPrintf(prefix, sizeof(prefix), "prefix=\"%s\", ", ((PetscObject)bA->m[i][j])->prefix);
898: PetscObjectTypeCompare((PetscObject)bA->m[i][j], MATNEST, &isNest);
900: PetscViewerASCIIPrintf(viewer, "(%" PetscInt_FMT ",%" PetscInt_FMT ") : %s%stype=%s, rows=%" PetscInt_FMT ", cols=%" PetscInt_FMT " \n", i, j, name, prefix, type, NR, NC);
902: if (isNest || viewSub) {
903: PetscViewerASCIIPushTab(viewer); /* push1 */
904: MatView(bA->m[i][j], viewer);
905: PetscViewerASCIIPopTab(viewer); /* pop1 */
906: }
907: }
908: }
909: PetscViewerASCIIPopTab(viewer); /* pop0 */
910: }
911: return 0;
912: }
914: static PetscErrorCode MatZeroEntries_Nest(Mat A)
915: {
916: Mat_Nest *bA = (Mat_Nest *)A->data;
917: PetscInt i, j;
919: for (i = 0; i < bA->nr; i++) {
920: for (j = 0; j < bA->nc; j++) {
921: if (!bA->m[i][j]) continue;
922: MatZeroEntries(bA->m[i][j]);
923: }
924: }
925: return 0;
926: }
928: static PetscErrorCode MatCopy_Nest(Mat A, Mat B, MatStructure str)
929: {
930: Mat_Nest *bA = (Mat_Nest *)A->data, *bB = (Mat_Nest *)B->data;
931: PetscInt i, j, nr = bA->nr, nc = bA->nc;
932: PetscBool nnzstate = PETSC_FALSE;
935: for (i = 0; i < nr; i++) {
936: for (j = 0; j < nc; j++) {
937: PetscObjectState subnnzstate = 0;
938: if (bA->m[i][j] && bB->m[i][j]) {
939: MatCopy(bA->m[i][j], bB->m[i][j], str);
941: MatGetNonzeroState(bB->m[i][j], &subnnzstate);
942: nnzstate = (PetscBool)(nnzstate || bB->nnzstate[i * nc + j] != subnnzstate);
943: bB->nnzstate[i * nc + j] = subnnzstate;
944: }
945: }
946: if (nnzstate) B->nonzerostate++;
947: return 0;
948: }
950: static PetscErrorCode MatAXPY_Nest(Mat Y, PetscScalar a, Mat X, MatStructure str)
951: {
952: Mat_Nest *bY = (Mat_Nest *)Y->data, *bX = (Mat_Nest *)X->data;
953: PetscInt i, j, nr = bY->nr, nc = bY->nc;
954: PetscBool nnzstate = PETSC_FALSE;
957: for (i = 0; i < nr; i++) {
958: for (j = 0; j < nc; j++) {
959: PetscObjectState subnnzstate = 0;
960: if (bY->m[i][j] && bX->m[i][j]) {
961: MatAXPY(bY->m[i][j], a, bX->m[i][j], str);
962: } else if (bX->m[i][j]) {
963: Mat M;
966: MatDuplicate(bX->m[i][j], MAT_COPY_VALUES, &M);
967: MatNestSetSubMat(Y, i, j, M);
968: MatDestroy(&M);
969: }
970: if (bY->m[i][j]) MatGetNonzeroState(bY->m[i][j], &subnnzstate);
971: nnzstate = (PetscBool)(nnzstate || bY->nnzstate[i * nc + j] != subnnzstate);
972: bY->nnzstate[i * nc + j] = subnnzstate;
973: }
974: }
975: if (nnzstate) Y->nonzerostate++;
976: return 0;
977: }
979: static PetscErrorCode MatDuplicate_Nest(Mat A, MatDuplicateOption op, Mat *B)
980: {
981: Mat_Nest *bA = (Mat_Nest *)A->data;
982: Mat *b;
983: PetscInt i, j, nr = bA->nr, nc = bA->nc;
985: PetscMalloc1(nr * nc, &b);
986: for (i = 0; i < nr; i++) {
987: for (j = 0; j < nc; j++) {
988: if (bA->m[i][j]) {
989: MatDuplicate(bA->m[i][j], op, &b[i * nc + j]);
990: } else {
991: b[i * nc + j] = NULL;
992: }
993: }
994: }
995: MatCreateNest(PetscObjectComm((PetscObject)A), nr, bA->isglobal.row, nc, bA->isglobal.col, b, B);
996: /* Give the new MatNest exclusive ownership */
997: for (i = 0; i < nr * nc; i++) MatDestroy(&b[i]);
998: PetscFree(b);
1000: MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY);
1001: MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY);
1002: return 0;
1003: }
1005: /* nest api */
1006: PetscErrorCode MatNestGetSubMat_Nest(Mat A, PetscInt idxm, PetscInt jdxm, Mat *mat)
1007: {
1008: Mat_Nest *bA = (Mat_Nest *)A->data;
1012: *mat = bA->m[idxm][jdxm];
1013: return 0;
1014: }
1016: /*@
1017: MatNestGetSubMat - Returns a single, sub-matrix from a `MATNEST`
1019: Not collective
1021: Input Parameters:
1022: + A - `MATNEST` matrix
1023: . idxm - index of the matrix within the nest matrix
1024: - jdxm - index of the matrix within the nest matrix
1026: Output Parameter:
1027: . sub - matrix at index idxm,jdxm within the nest matrix
1029: Level: developer
1031: .seealso: `MATNEST`, `MatNestGetSize()`, `MatNestGetSubMats()`, `MatCreateNest()`, `MATNEST`, `MatNestSetSubMat()`,
1032: `MatNestGetLocalISs()`, `MatNestGetISs()`
1033: @*/
1034: PetscErrorCode MatNestGetSubMat(Mat A, PetscInt idxm, PetscInt jdxm, Mat *sub)
1035: {
1036: PetscUseMethod(A, "MatNestGetSubMat_C", (Mat, PetscInt, PetscInt, Mat *), (A, idxm, jdxm, sub));
1037: return 0;
1038: }
1040: PetscErrorCode MatNestSetSubMat_Nest(Mat A, PetscInt idxm, PetscInt jdxm, Mat mat)
1041: {
1042: Mat_Nest *bA = (Mat_Nest *)A->data;
1043: PetscInt m, n, M, N, mi, ni, Mi, Ni;
1047: MatGetLocalSize(mat, &m, &n);
1048: MatGetSize(mat, &M, &N);
1049: ISGetLocalSize(bA->isglobal.row[idxm], &mi);
1050: ISGetSize(bA->isglobal.row[idxm], &Mi);
1051: ISGetLocalSize(bA->isglobal.col[jdxm], &ni);
1052: ISGetSize(bA->isglobal.col[jdxm], &Ni);
1056: /* do not increase object state */
1057: if (mat == bA->m[idxm][jdxm]) return 0;
1059: PetscObjectReference((PetscObject)mat);
1060: MatDestroy(&bA->m[idxm][jdxm]);
1061: bA->m[idxm][jdxm] = mat;
1062: PetscObjectStateIncrease((PetscObject)A);
1063: MatGetNonzeroState(mat, &bA->nnzstate[idxm * bA->nc + jdxm]);
1064: A->nonzerostate++;
1065: return 0;
1066: }
1068: /*@
1069: MatNestSetSubMat - Set a single submatrix in the `MATNEST`
1071: Logically collective on the submatrix communicator
1073: Input Parameters:
1074: + A - `MATNEST` matrix
1075: . idxm - index of the matrix within the nest matrix
1076: . jdxm - index of the matrix within the nest matrix
1077: - sub - matrix at index idxm,jdxm within the nest matrix
1079: Notes:
1080: The new submatrix must have the same size and communicator as that block of the nest.
1082: This increments the reference count of the submatrix.
1084: Level: developer
1086: .seealso: `MATNEST`, `MatNestSetSubMats()`, `MatNestGetSubMats()`, `MatNestGetLocalISs()`, `MATNEST`, `MatCreateNest()`,
1087: `MatNestGetSubMat()`, `MatNestGetISs()`, `MatNestGetSize()`
1088: @*/
1089: PetscErrorCode MatNestSetSubMat(Mat A, PetscInt idxm, PetscInt jdxm, Mat sub)
1090: {
1091: PetscUseMethod(A, "MatNestSetSubMat_C", (Mat, PetscInt, PetscInt, Mat), (A, idxm, jdxm, sub));
1092: return 0;
1093: }
1095: PetscErrorCode MatNestGetSubMats_Nest(Mat A, PetscInt *M, PetscInt *N, Mat ***mat)
1096: {
1097: Mat_Nest *bA = (Mat_Nest *)A->data;
1099: if (M) *M = bA->nr;
1100: if (N) *N = bA->nc;
1101: if (mat) *mat = bA->m;
1102: return 0;
1103: }
1105: /*@C
1106: MatNestGetSubMats - Returns the entire two dimensional array of matrices defining a `MATNEST` matrix.
1108: Not collective
1110: Input Parameter:
1111: . A - nest matrix
1113: Output Parameters:
1114: + M - number of rows in the nest matrix
1115: . N - number of cols in the nest matrix
1116: - mat - 2d array of matrices
1118: Note:
1119: The user should not free the array mat.
1121: Fortran Note:
1122: This routine has a calling sequence
1123: $ call MatNestGetSubMats(A, M, N, mat, ierr)
1124: where the space allocated for the optional argument mat is assumed large enough (if provided).
1126: Level: developer
1128: .seealso: `MATNEST`, `MatNestGetSize()`, `MatNestGetSubMat()`, `MatNestGetLocalISs()`, `MATNEST`, `MatCreateNest()`,
1129: `MatNestSetSubMats()`, `MatNestGetISs()`, `MatNestSetSubMat()`
1130: @*/
1131: PetscErrorCode MatNestGetSubMats(Mat A, PetscInt *M, PetscInt *N, Mat ***mat)
1132: {
1133: PetscUseMethod(A, "MatNestGetSubMats_C", (Mat, PetscInt *, PetscInt *, Mat ***), (A, M, N, mat));
1134: return 0;
1135: }
1137: PetscErrorCode MatNestGetSize_Nest(Mat A, PetscInt *M, PetscInt *N)
1138: {
1139: Mat_Nest *bA = (Mat_Nest *)A->data;
1141: if (M) *M = bA->nr;
1142: if (N) *N = bA->nc;
1143: return 0;
1144: }
1146: /*@
1147: MatNestGetSize - Returns the size of the `MATNEST` matrix.
1149: Not collective
1151: Input Parameter:
1152: . A - `MATNEST` matrix
1154: Output Parameters:
1155: + M - number of rows in the nested mat
1156: - N - number of cols in the nested mat
1158: Level: developer
1160: .seealso: `MATNEST`, `MatNestGetSubMat()`, `MatNestGetSubMats()`, `MATNEST`, `MatCreateNest()`, `MatNestGetLocalISs()`,
1161: `MatNestGetISs()`
1162: @*/
1163: PetscErrorCode MatNestGetSize(Mat A, PetscInt *M, PetscInt *N)
1164: {
1165: PetscUseMethod(A, "MatNestGetSize_C", (Mat, PetscInt *, PetscInt *), (A, M, N));
1166: return 0;
1167: }
1169: static PetscErrorCode MatNestGetISs_Nest(Mat A, IS rows[], IS cols[])
1170: {
1171: Mat_Nest *vs = (Mat_Nest *)A->data;
1172: PetscInt i;
1174: if (rows)
1175: for (i = 0; i < vs->nr; i++) rows[i] = vs->isglobal.row[i];
1176: if (cols)
1177: for (i = 0; i < vs->nc; i++) cols[i] = vs->isglobal.col[i];
1178: return 0;
1179: }
1181: /*@C
1182: MatNestGetISs - Returns the index sets partitioning the row and column spaces of a `MATNEST`
1184: Not collective
1186: Input Parameter:
1187: . A - `MATNEST` matrix
1189: Output Parameters:
1190: + rows - array of row index sets
1191: - cols - array of column index sets
1193: Level: advanced
1195: Note:
1196: The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs.
1198: .seealso: `MATNEST`, `MatNestGetSubMat()`, `MatNestGetSubMats()`, `MatNestGetSize()`, `MatNestGetLocalISs()`, `MATNEST`,
1199: `MatCreateNest()`, `MatNestGetSubMats()`, `MatNestSetSubMats()`
1200: @*/
1201: PetscErrorCode MatNestGetISs(Mat A, IS rows[], IS cols[])
1202: {
1204: PetscUseMethod(A, "MatNestGetISs_C", (Mat, IS[], IS[]), (A, rows, cols));
1205: return 0;
1206: }
1208: static PetscErrorCode MatNestGetLocalISs_Nest(Mat A, IS rows[], IS cols[])
1209: {
1210: Mat_Nest *vs = (Mat_Nest *)A->data;
1211: PetscInt i;
1213: if (rows)
1214: for (i = 0; i < vs->nr; i++) rows[i] = vs->islocal.row[i];
1215: if (cols)
1216: for (i = 0; i < vs->nc; i++) cols[i] = vs->islocal.col[i];
1217: return 0;
1218: }
1220: /*@C
1221: MatNestGetLocalISs - Returns the index sets partitioning the row and column spaces of a `MATNEST`
1223: Not collective
1225: Input Parameter:
1226: . A - `MATNEST` matrix
1228: Output Parameters:
1229: + rows - array of row index sets (or NULL to ignore)
1230: - cols - array of column index sets (or NULL to ignore)
1232: Level: advanced
1234: Note:
1235: The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs.
1237: .seealso: `MATNEST`, `MatNestGetSubMat()`, `MatNestGetSubMats()`, `MatNestGetSize()`, `MatNestGetISs()`, `MatCreateNest()`,
1238: `MATNEST`, `MatNestSetSubMats()`, `MatNestSetSubMat()`
1239: @*/
1240: PetscErrorCode MatNestGetLocalISs(Mat A, IS rows[], IS cols[])
1241: {
1243: PetscUseMethod(A, "MatNestGetLocalISs_C", (Mat, IS[], IS[]), (A, rows, cols));
1244: return 0;
1245: }
1247: PetscErrorCode MatNestSetVecType_Nest(Mat A, VecType vtype)
1248: {
1249: PetscBool flg;
1251: PetscStrcmp(vtype, VECNEST, &flg);
1252: /* In reality, this only distinguishes VECNEST and "other" */
1253: if (flg) A->ops->getvecs = MatCreateVecs_Nest;
1254: else A->ops->getvecs = (PetscErrorCode(*)(Mat, Vec *, Vec *))0;
1255: return 0;
1256: }
1258: /*@C
1259: MatNestSetVecType - Sets the type of `Vec` returned by `MatCreateVecs()`
1261: Not collective
1263: Input Parameters:
1264: + A - `MATNEST` matrix
1265: - vtype - `VecType` to use for creating vectors
1267: Level: developer
1269: .seealso: `MATNEST`, `MatCreateVecs()`, `MATNEST`, `MatCreateNest()`, `VecType`
1270: @*/
1271: PetscErrorCode MatNestSetVecType(Mat A, VecType vtype)
1272: {
1273: PetscTryMethod(A, "MatNestSetVecType_C", (Mat, VecType), (A, vtype));
1274: return 0;
1275: }
1277: PetscErrorCode MatNestSetSubMats_Nest(Mat A, PetscInt nr, const IS is_row[], PetscInt nc, const IS is_col[], const Mat a[])
1278: {
1279: Mat_Nest *s = (Mat_Nest *)A->data;
1280: PetscInt i, j, m, n, M, N;
1281: PetscBool cong, isstd, sametype = PETSC_FALSE;
1282: VecType vtype, type;
1284: MatReset_Nest(A);
1286: s->nr = nr;
1287: s->nc = nc;
1289: /* Create space for submatrices */
1290: PetscMalloc1(nr, &s->m);
1291: for (i = 0; i < nr; i++) PetscMalloc1(nc, &s->m[i]);
1292: for (i = 0; i < nr; i++) {
1293: for (j = 0; j < nc; j++) {
1294: s->m[i][j] = a[i * nc + j];
1295: if (a[i * nc + j]) PetscObjectReference((PetscObject)a[i * nc + j]);
1296: }
1297: }
1298: MatGetVecType(A, &vtype);
1299: PetscStrcmp(vtype, VECSTANDARD, &isstd);
1300: if (isstd) {
1301: /* check if all blocks have the same vectype */
1302: vtype = NULL;
1303: for (i = 0; i < nr; i++) {
1304: for (j = 0; j < nc; j++) {
1305: if (a[i * nc + j]) {
1306: if (!vtype) { /* first visited block */
1307: MatGetVecType(a[i * nc + j], &vtype);
1308: sametype = PETSC_TRUE;
1309: } else if (sametype) {
1310: MatGetVecType(a[i * nc + j], &type);
1311: PetscStrcmp(vtype, type, &sametype);
1312: }
1313: }
1314: }
1315: }
1316: if (sametype) { /* propagate vectype */
1317: MatSetVecType(A, vtype);
1318: }
1319: }
1321: MatSetUp_NestIS_Private(A, nr, is_row, nc, is_col);
1323: PetscMalloc1(nr, &s->row_len);
1324: PetscMalloc1(nc, &s->col_len);
1325: for (i = 0; i < nr; i++) s->row_len[i] = -1;
1326: for (j = 0; j < nc; j++) s->col_len[j] = -1;
1328: PetscCalloc1(nr * nc, &s->nnzstate);
1329: for (i = 0; i < nr; i++) {
1330: for (j = 0; j < nc; j++) {
1331: if (s->m[i][j]) MatGetNonzeroState(s->m[i][j], &s->nnzstate[i * nc + j]);
1332: }
1333: }
1335: MatNestGetSizes_Private(A, &m, &n, &M, &N);
1337: PetscLayoutSetSize(A->rmap, M);
1338: PetscLayoutSetLocalSize(A->rmap, m);
1339: PetscLayoutSetSize(A->cmap, N);
1340: PetscLayoutSetLocalSize(A->cmap, n);
1342: PetscLayoutSetUp(A->rmap);
1343: PetscLayoutSetUp(A->cmap);
1345: /* disable operations that are not supported for non-square matrices,
1346: or matrices for which is_row != is_col */
1347: MatHasCongruentLayouts(A, &cong);
1348: if (cong && nr != nc) cong = PETSC_FALSE;
1349: if (cong) {
1350: for (i = 0; cong && i < nr; i++) ISEqualUnsorted(s->isglobal.row[i], s->isglobal.col[i], &cong);
1351: }
1352: if (!cong) {
1353: A->ops->missingdiagonal = NULL;
1354: A->ops->getdiagonal = NULL;
1355: A->ops->shift = NULL;
1356: A->ops->diagonalset = NULL;
1357: }
1359: PetscCalloc2(nr, &s->left, nc, &s->right);
1360: PetscObjectStateIncrease((PetscObject)A);
1361: A->nonzerostate++;
1362: return 0;
1363: }
1365: /*@
1366: MatNestSetSubMats - Sets the nested submatrices in a `MATNEST`
1368: Collective
1370: Input Parameters:
1371: + A - `MATNEST` matrix
1372: . nr - number of nested row blocks
1373: . is_row - index sets for each nested row block, or NULL to make contiguous
1374: . nc - number of nested column blocks
1375: . is_col - index sets for each nested column block, or NULL to make contiguous
1376: - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL
1378: Note:
1379: This always resets any submatrix information previously set
1381: Level: advanced
1383: .seealso: `MATNEST`, `MatCreateNest()`, `MATNEST`, `MatNestSetSubMat()`, `MatNestGetSubMat()`, `MatNestGetSubMats()`
1384: @*/
1385: PetscErrorCode MatNestSetSubMats(Mat A, PetscInt nr, const IS is_row[], PetscInt nc, const IS is_col[], const Mat a[])
1386: {
1387: PetscInt i;
1391: if (nr && is_row) {
1394: }
1396: if (nc && is_col) {
1399: }
1401: PetscUseMethod(A, "MatNestSetSubMats_C", (Mat, PetscInt, const IS[], PetscInt, const IS[], const Mat[]), (A, nr, is_row, nc, is_col, a));
1402: return 0;
1403: }
1405: static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A, PetscInt n, const IS islocal[], const IS isglobal[], PetscBool colflg, ISLocalToGlobalMapping *ltog)
1406: {
1407: PetscBool flg;
1408: PetscInt i, j, m, mi, *ix;
1410: *ltog = NULL;
1411: for (i = 0, m = 0, flg = PETSC_FALSE; i < n; i++) {
1412: if (islocal[i]) {
1413: ISGetLocalSize(islocal[i], &mi);
1414: flg = PETSC_TRUE; /* We found a non-trivial entry */
1415: } else {
1416: ISGetLocalSize(isglobal[i], &mi);
1417: }
1418: m += mi;
1419: }
1420: if (!flg) return 0;
1422: PetscMalloc1(m, &ix);
1423: for (i = 0, m = 0; i < n; i++) {
1424: ISLocalToGlobalMapping smap = NULL;
1425: Mat sub = NULL;
1426: PetscSF sf;
1427: PetscLayout map;
1428: const PetscInt *ix2;
1430: if (!colflg) {
1431: MatNestFindNonzeroSubMatRow(A, i, &sub);
1432: } else {
1433: MatNestFindNonzeroSubMatCol(A, i, &sub);
1434: }
1435: if (sub) {
1436: if (!colflg) {
1437: MatGetLocalToGlobalMapping(sub, &smap, NULL);
1438: } else {
1439: MatGetLocalToGlobalMapping(sub, NULL, &smap);
1440: }
1441: }
1442: /*
1443: Now we need to extract the monolithic global indices that correspond to the given split global indices.
1444: In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces.
1445: */
1446: ISGetIndices(isglobal[i], &ix2);
1447: if (islocal[i]) {
1448: PetscInt *ilocal, *iremote;
1449: PetscInt mil, nleaves;
1451: ISGetLocalSize(islocal[i], &mi);
1453: for (j = 0; j < mi; j++) ix[m + j] = j;
1454: ISLocalToGlobalMappingApply(smap, mi, ix + m, ix + m);
1456: /* PetscSFSetGraphLayout does not like negative indices */
1457: PetscMalloc2(mi, &ilocal, mi, &iremote);
1458: for (j = 0, nleaves = 0; j < mi; j++) {
1459: if (ix[m + j] < 0) continue;
1460: ilocal[nleaves] = j;
1461: iremote[nleaves] = ix[m + j];
1462: nleaves++;
1463: }
1464: ISGetLocalSize(isglobal[i], &mil);
1465: PetscSFCreate(PetscObjectComm((PetscObject)A), &sf);
1466: PetscLayoutCreate(PetscObjectComm((PetscObject)A), &map);
1467: PetscLayoutSetLocalSize(map, mil);
1468: PetscLayoutSetUp(map);
1469: PetscSFSetGraphLayout(sf, map, nleaves, ilocal, PETSC_USE_POINTER, iremote);
1470: PetscLayoutDestroy(&map);
1471: PetscSFBcastBegin(sf, MPIU_INT, ix2, ix + m, MPI_REPLACE);
1472: PetscSFBcastEnd(sf, MPIU_INT, ix2, ix + m, MPI_REPLACE);
1473: PetscSFDestroy(&sf);
1474: PetscFree2(ilocal, iremote);
1475: } else {
1476: ISGetLocalSize(isglobal[i], &mi);
1477: for (j = 0; j < mi; j++) ix[m + j] = ix2[i];
1478: }
1479: ISRestoreIndices(isglobal[i], &ix2);
1480: m += mi;
1481: }
1482: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A), 1, m, ix, PETSC_OWN_POINTER, ltog);
1483: return 0;
1484: }
1486: /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */
1487: /*
1488: nprocessors = NP
1489: Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1))
1490: proc 0: => (g_0,h_0,)
1491: proc 1: => (g_1,h_1,)
1492: ...
1493: proc nprocs-1: => (g_NP-1,h_NP-1,)
1495: proc 0: proc 1: proc nprocs-1:
1496: is[0] = (0,1,2,...,nlocal(g_0)-1) (0,1,...,nlocal(g_1)-1) (0,1,...,nlocal(g_NP-1))
1498: proc 0:
1499: is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1)
1500: proc 1:
1501: is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1)
1503: proc NP-1:
1504: is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1)
1505: */
1506: static PetscErrorCode MatSetUp_NestIS_Private(Mat A, PetscInt nr, const IS is_row[], PetscInt nc, const IS is_col[])
1507: {
1508: Mat_Nest *vs = (Mat_Nest *)A->data;
1509: PetscInt i, j, offset, n, nsum, bs;
1510: Mat sub = NULL;
1512: PetscMalloc1(nr, &vs->isglobal.row);
1513: PetscMalloc1(nc, &vs->isglobal.col);
1514: if (is_row) { /* valid IS is passed in */
1515: /* refs on is[] are incremented */
1516: for (i = 0; i < vs->nr; i++) {
1517: PetscObjectReference((PetscObject)is_row[i]);
1519: vs->isglobal.row[i] = is_row[i];
1520: }
1521: } else { /* Create the ISs by inspecting sizes of a submatrix in each row */
1522: nsum = 0;
1523: for (i = 0; i < vs->nr; i++) { /* Add up the local sizes to compute the aggregate offset */
1524: MatNestFindNonzeroSubMatRow(A, i, &sub);
1526: MatGetLocalSize(sub, &n, NULL);
1528: nsum += n;
1529: }
1530: MPI_Scan(&nsum, &offset, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)A));
1531: offset -= nsum;
1532: for (i = 0; i < vs->nr; i++) {
1533: MatNestFindNonzeroSubMatRow(A, i, &sub);
1534: MatGetLocalSize(sub, &n, NULL);
1535: MatGetBlockSizes(sub, &bs, NULL);
1536: ISCreateStride(PetscObjectComm((PetscObject)sub), n, offset, 1, &vs->isglobal.row[i]);
1537: ISSetBlockSize(vs->isglobal.row[i], bs);
1538: offset += n;
1539: }
1540: }
1542: if (is_col) { /* valid IS is passed in */
1543: /* refs on is[] are incremented */
1544: for (j = 0; j < vs->nc; j++) {
1545: PetscObjectReference((PetscObject)is_col[j]);
1547: vs->isglobal.col[j] = is_col[j];
1548: }
1549: } else { /* Create the ISs by inspecting sizes of a submatrix in each column */
1550: offset = A->cmap->rstart;
1551: nsum = 0;
1552: for (j = 0; j < vs->nc; j++) {
1553: MatNestFindNonzeroSubMatCol(A, j, &sub);
1555: MatGetLocalSize(sub, NULL, &n);
1557: nsum += n;
1558: }
1559: MPI_Scan(&nsum, &offset, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)A));
1560: offset -= nsum;
1561: for (j = 0; j < vs->nc; j++) {
1562: MatNestFindNonzeroSubMatCol(A, j, &sub);
1563: MatGetLocalSize(sub, NULL, &n);
1564: MatGetBlockSizes(sub, NULL, &bs);
1565: ISCreateStride(PetscObjectComm((PetscObject)sub), n, offset, 1, &vs->isglobal.col[j]);
1566: ISSetBlockSize(vs->isglobal.col[j], bs);
1567: offset += n;
1568: }
1569: }
1571: /* Set up the local ISs */
1572: PetscMalloc1(vs->nr, &vs->islocal.row);
1573: PetscMalloc1(vs->nc, &vs->islocal.col);
1574: for (i = 0, offset = 0; i < vs->nr; i++) {
1575: IS isloc;
1576: ISLocalToGlobalMapping rmap = NULL;
1577: PetscInt nlocal, bs;
1578: MatNestFindNonzeroSubMatRow(A, i, &sub);
1579: if (sub) MatGetLocalToGlobalMapping(sub, &rmap, NULL);
1580: if (rmap) {
1581: MatGetBlockSizes(sub, &bs, NULL);
1582: ISLocalToGlobalMappingGetSize(rmap, &nlocal);
1583: ISCreateStride(PETSC_COMM_SELF, nlocal, offset, 1, &isloc);
1584: ISSetBlockSize(isloc, bs);
1585: } else {
1586: nlocal = 0;
1587: isloc = NULL;
1588: }
1589: vs->islocal.row[i] = isloc;
1590: offset += nlocal;
1591: }
1592: for (i = 0, offset = 0; i < vs->nc; i++) {
1593: IS isloc;
1594: ISLocalToGlobalMapping cmap = NULL;
1595: PetscInt nlocal, bs;
1596: MatNestFindNonzeroSubMatCol(A, i, &sub);
1597: if (sub) MatGetLocalToGlobalMapping(sub, NULL, &cmap);
1598: if (cmap) {
1599: MatGetBlockSizes(sub, NULL, &bs);
1600: ISLocalToGlobalMappingGetSize(cmap, &nlocal);
1601: ISCreateStride(PETSC_COMM_SELF, nlocal, offset, 1, &isloc);
1602: ISSetBlockSize(isloc, bs);
1603: } else {
1604: nlocal = 0;
1605: isloc = NULL;
1606: }
1607: vs->islocal.col[i] = isloc;
1608: offset += nlocal;
1609: }
1611: /* Set up the aggregate ISLocalToGlobalMapping */
1612: {
1613: ISLocalToGlobalMapping rmap, cmap;
1614: MatNestCreateAggregateL2G_Private(A, vs->nr, vs->islocal.row, vs->isglobal.row, PETSC_FALSE, &rmap);
1615: MatNestCreateAggregateL2G_Private(A, vs->nc, vs->islocal.col, vs->isglobal.col, PETSC_TRUE, &cmap);
1616: if (rmap && cmap) MatSetLocalToGlobalMapping(A, rmap, cmap);
1617: ISLocalToGlobalMappingDestroy(&rmap);
1618: ISLocalToGlobalMappingDestroy(&cmap);
1619: }
1621: if (PetscDefined(USE_DEBUG)) {
1622: for (i = 0; i < vs->nr; i++) {
1623: for (j = 0; j < vs->nc; j++) {
1624: PetscInt m, n, M, N, mi, ni, Mi, Ni;
1625: Mat B = vs->m[i][j];
1626: if (!B) continue;
1627: MatGetSize(B, &M, &N);
1628: MatGetLocalSize(B, &m, &n);
1629: ISGetSize(vs->isglobal.row[i], &Mi);
1630: ISGetSize(vs->isglobal.col[j], &Ni);
1631: ISGetLocalSize(vs->isglobal.row[i], &mi);
1632: ISGetLocalSize(vs->isglobal.col[j], &ni);
1635: }
1636: }
1637: }
1639: /* Set A->assembled if all non-null blocks are currently assembled */
1640: for (i = 0; i < vs->nr; i++) {
1641: for (j = 0; j < vs->nc; j++) {
1642: if (vs->m[i][j] && !vs->m[i][j]->assembled) return 0;
1643: }
1644: }
1645: A->assembled = PETSC_TRUE;
1646: return 0;
1647: }
1649: /*@C
1650: MatCreateNest - Creates a new `MATNEST` matrix containing several nested submatrices, each stored separately
1652: Collective
1654: Input Parameters:
1655: + comm - Communicator for the new `MATNEST`
1656: . nr - number of nested row blocks
1657: . is_row - index sets for each nested row block, or NULL to make contiguous
1658: . nc - number of nested column blocks
1659: . is_col - index sets for each nested column block, or NULL to make contiguous
1660: - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL
1662: Output Parameter:
1663: . B - new matrix
1665: Level: advanced
1667: .seealso: `MATNEST`, `MatCreate()`, `VecCreateNest()`, `DMCreateMatrix()`, `MATNEST`, `MatNestSetSubMat()`,
1668: `MatNestGetSubMat()`, `MatNestGetLocalISs()`, `MatNestGetSize()`,
1669: `MatNestGetISs()`, `MatNestSetSubMats()`, `MatNestGetSubMats()`
1670: @*/
1671: PetscErrorCode MatCreateNest(MPI_Comm comm, PetscInt nr, const IS is_row[], PetscInt nc, const IS is_col[], const Mat a[], Mat *B)
1672: {
1673: Mat A;
1675: *B = NULL;
1676: MatCreate(comm, &A);
1677: MatSetType(A, MATNEST);
1678: A->preallocated = PETSC_TRUE;
1679: MatNestSetSubMats(A, nr, is_row, nc, is_col, a);
1680: *B = A;
1681: return 0;
1682: }
1684: PetscErrorCode MatConvert_Nest_SeqAIJ_fast(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
1685: {
1686: Mat_Nest *nest = (Mat_Nest *)A->data;
1687: Mat *trans;
1688: PetscScalar **avv;
1689: PetscScalar *vv;
1690: PetscInt **aii, **ajj;
1691: PetscInt *ii, *jj, *ci;
1692: PetscInt nr, nc, nnz, i, j;
1693: PetscBool done;
1695: MatGetSize(A, &nr, &nc);
1696: if (reuse == MAT_REUSE_MATRIX) {
1697: PetscInt rnr;
1699: MatGetRowIJ(*newmat, 0, PETSC_FALSE, PETSC_FALSE, &rnr, (const PetscInt **)&ii, (const PetscInt **)&jj, &done);
1702: MatSeqAIJGetArray(*newmat, &vv);
1703: }
1704: /* extract CSR for nested SeqAIJ matrices */
1705: nnz = 0;
1706: PetscCalloc4(nest->nr * nest->nc, &aii, nest->nr * nest->nc, &ajj, nest->nr * nest->nc, &avv, nest->nr * nest->nc, &trans);
1707: for (i = 0; i < nest->nr; ++i) {
1708: for (j = 0; j < nest->nc; ++j) {
1709: Mat B = nest->m[i][j];
1710: if (B) {
1711: PetscScalar *naa;
1712: PetscInt *nii, *njj, nnr;
1713: PetscBool istrans;
1715: PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &istrans);
1716: if (istrans) {
1717: Mat Bt;
1719: MatTransposeGetMat(B, &Bt);
1720: MatTranspose(Bt, MAT_INITIAL_MATRIX, &trans[i * nest->nc + j]);
1721: B = trans[i * nest->nc + j];
1722: } else {
1723: PetscObjectTypeCompare((PetscObject)B, MATHERMITIANTRANSPOSEVIRTUAL, &istrans);
1724: if (istrans) {
1725: Mat Bt;
1727: MatHermitianTransposeGetMat(B, &Bt);
1728: MatHermitianTranspose(Bt, MAT_INITIAL_MATRIX, &trans[i * nest->nc + j]);
1729: B = trans[i * nest->nc + j];
1730: }
1731: }
1732: MatGetRowIJ(B, 0, PETSC_FALSE, PETSC_FALSE, &nnr, (const PetscInt **)&nii, (const PetscInt **)&njj, &done);
1734: MatSeqAIJGetArray(B, &naa);
1735: nnz += nii[nnr];
1737: aii[i * nest->nc + j] = nii;
1738: ajj[i * nest->nc + j] = njj;
1739: avv[i * nest->nc + j] = naa;
1740: }
1741: }
1742: }
1743: if (reuse != MAT_REUSE_MATRIX) {
1744: PetscMalloc1(nr + 1, &ii);
1745: PetscMalloc1(nnz, &jj);
1746: PetscMalloc1(nnz, &vv);
1747: } else {
1749: }
1751: /* new row pointer */
1752: PetscArrayzero(ii, nr + 1);
1753: for (i = 0; i < nest->nr; ++i) {
1754: PetscInt ncr, rst;
1756: ISStrideGetInfo(nest->isglobal.row[i], &rst, NULL);
1757: ISGetLocalSize(nest->isglobal.row[i], &ncr);
1758: for (j = 0; j < nest->nc; ++j) {
1759: if (aii[i * nest->nc + j]) {
1760: PetscInt *nii = aii[i * nest->nc + j];
1761: PetscInt ir;
1763: for (ir = rst; ir < ncr + rst; ++ir) {
1764: ii[ir + 1] += nii[1] - nii[0];
1765: nii++;
1766: }
1767: }
1768: }
1769: }
1770: for (i = 0; i < nr; i++) ii[i + 1] += ii[i];
1772: /* construct CSR for the new matrix */
1773: PetscCalloc1(nr, &ci);
1774: for (i = 0; i < nest->nr; ++i) {
1775: PetscInt ncr, rst;
1777: ISStrideGetInfo(nest->isglobal.row[i], &rst, NULL);
1778: ISGetLocalSize(nest->isglobal.row[i], &ncr);
1779: for (j = 0; j < nest->nc; ++j) {
1780: if (aii[i * nest->nc + j]) {
1781: PetscScalar *nvv = avv[i * nest->nc + j];
1782: PetscInt *nii = aii[i * nest->nc + j];
1783: PetscInt *njj = ajj[i * nest->nc + j];
1784: PetscInt ir, cst;
1786: ISStrideGetInfo(nest->isglobal.col[j], &cst, NULL);
1787: for (ir = rst; ir < ncr + rst; ++ir) {
1788: PetscInt ij, rsize = nii[1] - nii[0], ist = ii[ir] + ci[ir];
1790: for (ij = 0; ij < rsize; ij++) {
1791: jj[ist + ij] = *njj + cst;
1792: vv[ist + ij] = *nvv;
1793: njj++;
1794: nvv++;
1795: }
1796: ci[ir] += rsize;
1797: nii++;
1798: }
1799: }
1800: }
1801: }
1802: PetscFree(ci);
1804: /* restore info */
1805: for (i = 0; i < nest->nr; ++i) {
1806: for (j = 0; j < nest->nc; ++j) {
1807: Mat B = nest->m[i][j];
1808: if (B) {
1809: PetscInt nnr = 0, k = i * nest->nc + j;
1811: B = (trans[k] ? trans[k] : B);
1812: MatRestoreRowIJ(B, 0, PETSC_FALSE, PETSC_FALSE, &nnr, (const PetscInt **)&aii[k], (const PetscInt **)&ajj[k], &done);
1814: MatSeqAIJRestoreArray(B, &avv[k]);
1815: MatDestroy(&trans[k]);
1816: }
1817: }
1818: }
1819: PetscFree4(aii, ajj, avv, trans);
1821: /* finalize newmat */
1822: if (reuse == MAT_INITIAL_MATRIX) {
1823: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A), nr, nc, ii, jj, vv, newmat);
1824: } else if (reuse == MAT_INPLACE_MATRIX) {
1825: Mat B;
1827: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A), nr, nc, ii, jj, vv, &B);
1828: MatHeaderReplace(A, &B);
1829: }
1830: MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY);
1831: MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY);
1832: {
1833: Mat_SeqAIJ *a = (Mat_SeqAIJ *)((*newmat)->data);
1834: a->free_a = PETSC_TRUE;
1835: a->free_ij = PETSC_TRUE;
1836: }
1837: return 0;
1838: }
1840: PETSC_INTERN PetscErrorCode MatAXPY_Dense_Nest(Mat Y, PetscScalar a, Mat X)
1841: {
1842: Mat_Nest *nest = (Mat_Nest *)X->data;
1843: PetscInt i, j, k, rstart;
1844: PetscBool flg;
1846: /* Fill by row */
1847: for (j = 0; j < nest->nc; ++j) {
1848: /* Using global column indices and ISAllGather() is not scalable. */
1849: IS bNis;
1850: PetscInt bN;
1851: const PetscInt *bNindices;
1852: ISAllGather(nest->isglobal.col[j], &bNis);
1853: ISGetSize(bNis, &bN);
1854: ISGetIndices(bNis, &bNindices);
1855: for (i = 0; i < nest->nr; ++i) {
1856: Mat B = nest->m[i][j], D = NULL;
1857: PetscInt bm, br;
1858: const PetscInt *bmindices;
1859: if (!B) continue;
1860: PetscObjectTypeCompareAny((PetscObject)B, &flg, MATTRANSPOSEVIRTUAL, MATHERMITIANTRANSPOSEVIRTUAL, "");
1861: if (flg) {
1862: PetscTryMethod(B, "MatTransposeGetMat_C", (Mat, Mat *), (B, &D));
1863: PetscTryMethod(B, "MatHermitianTransposeGetMat_C", (Mat, Mat *), (B, &D));
1864: MatConvert(B, ((PetscObject)D)->type_name, MAT_INITIAL_MATRIX, &D);
1865: B = D;
1866: }
1867: PetscObjectTypeCompareAny((PetscObject)B, &flg, MATSEQSBAIJ, MATMPISBAIJ, "");
1868: if (flg) {
1869: if (D) MatConvert(D, MATBAIJ, MAT_INPLACE_MATRIX, &D);
1870: else MatConvert(B, MATBAIJ, MAT_INITIAL_MATRIX, &D);
1871: B = D;
1872: }
1873: ISGetLocalSize(nest->isglobal.row[i], &bm);
1874: ISGetIndices(nest->isglobal.row[i], &bmindices);
1875: MatGetOwnershipRange(B, &rstart, NULL);
1876: for (br = 0; br < bm; ++br) {
1877: PetscInt row = bmindices[br], brncols, *cols;
1878: const PetscInt *brcols;
1879: const PetscScalar *brcoldata;
1880: PetscScalar *vals = NULL;
1881: MatGetRow(B, br + rstart, &brncols, &brcols, &brcoldata);
1882: PetscMalloc1(brncols, &cols);
1883: for (k = 0; k < brncols; k++) cols[k] = bNindices[brcols[k]];
1884: /*
1885: Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match.
1886: Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES.
1887: */
1888: if (a != 1.0) {
1889: PetscMalloc1(brncols, &vals);
1890: for (k = 0; k < brncols; k++) vals[k] = a * brcoldata[k];
1891: MatSetValues(Y, 1, &row, brncols, cols, vals, ADD_VALUES);
1892: PetscFree(vals);
1893: } else {
1894: MatSetValues(Y, 1, &row, brncols, cols, brcoldata, ADD_VALUES);
1895: }
1896: MatRestoreRow(B, br + rstart, &brncols, &brcols, &brcoldata);
1897: PetscFree(cols);
1898: }
1899: if (D) MatDestroy(&D);
1900: ISRestoreIndices(nest->isglobal.row[i], &bmindices);
1901: }
1902: ISRestoreIndices(bNis, &bNindices);
1903: ISDestroy(&bNis);
1904: }
1905: MatAssemblyBegin(Y, MAT_FINAL_ASSEMBLY);
1906: MatAssemblyEnd(Y, MAT_FINAL_ASSEMBLY);
1907: return 0;
1908: }
1910: PetscErrorCode MatConvert_Nest_AIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
1911: {
1912: Mat_Nest *nest = (Mat_Nest *)A->data;
1913: PetscInt m, n, M, N, i, j, k, *dnnz, *onnz, rstart, cstart, cend;
1914: PetscMPIInt size;
1915: Mat C;
1917: MPI_Comm_size(PetscObjectComm((PetscObject)A), &size);
1918: if (size == 1) { /* look for a special case with SeqAIJ matrices and strided-1, contiguous, blocks */
1919: PetscInt nf;
1920: PetscBool fast;
1922: PetscStrcmp(newtype, MATAIJ, &fast);
1923: if (!fast) PetscStrcmp(newtype, MATSEQAIJ, &fast);
1924: for (i = 0; i < nest->nr && fast; ++i) {
1925: for (j = 0; j < nest->nc && fast; ++j) {
1926: Mat B = nest->m[i][j];
1927: if (B) {
1928: PetscObjectTypeCompare((PetscObject)B, MATSEQAIJ, &fast);
1929: if (!fast) {
1930: PetscBool istrans;
1932: PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &istrans);
1933: if (istrans) {
1934: Mat Bt;
1936: MatTransposeGetMat(B, &Bt);
1937: PetscObjectTypeCompare((PetscObject)Bt, MATSEQAIJ, &fast);
1938: } else {
1939: PetscObjectTypeCompare((PetscObject)B, MATHERMITIANTRANSPOSEVIRTUAL, &istrans);
1940: if (istrans) {
1941: Mat Bt;
1943: MatHermitianTransposeGetMat(B, &Bt);
1944: PetscObjectTypeCompare((PetscObject)Bt, MATSEQAIJ, &fast);
1945: }
1946: }
1947: }
1948: }
1949: }
1950: }
1951: for (i = 0, nf = 0; i < nest->nr && fast; ++i) {
1952: PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i], ISSTRIDE, &fast);
1953: if (fast) {
1954: PetscInt f, s;
1956: ISStrideGetInfo(nest->isglobal.row[i], &f, &s);
1957: if (f != nf || s != 1) {
1958: fast = PETSC_FALSE;
1959: } else {
1960: ISGetSize(nest->isglobal.row[i], &f);
1961: nf += f;
1962: }
1963: }
1964: }
1965: for (i = 0, nf = 0; i < nest->nc && fast; ++i) {
1966: PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i], ISSTRIDE, &fast);
1967: if (fast) {
1968: PetscInt f, s;
1970: ISStrideGetInfo(nest->isglobal.col[i], &f, &s);
1971: if (f != nf || s != 1) {
1972: fast = PETSC_FALSE;
1973: } else {
1974: ISGetSize(nest->isglobal.col[i], &f);
1975: nf += f;
1976: }
1977: }
1978: }
1979: if (fast) {
1980: MatConvert_Nest_SeqAIJ_fast(A, newtype, reuse, newmat);
1981: return 0;
1982: }
1983: }
1984: MatGetSize(A, &M, &N);
1985: MatGetLocalSize(A, &m, &n);
1986: MatGetOwnershipRangeColumn(A, &cstart, &cend);
1987: if (reuse == MAT_REUSE_MATRIX) C = *newmat;
1988: else {
1989: MatCreate(PetscObjectComm((PetscObject)A), &C);
1990: MatSetType(C, newtype);
1991: MatSetSizes(C, m, n, M, N);
1992: }
1993: PetscMalloc1(2 * m, &dnnz);
1994: onnz = dnnz + m;
1995: for (k = 0; k < m; k++) {
1996: dnnz[k] = 0;
1997: onnz[k] = 0;
1998: }
1999: for (j = 0; j < nest->nc; ++j) {
2000: IS bNis;
2001: PetscInt bN;
2002: const PetscInt *bNindices;
2003: PetscBool flg;
2004: /* Using global column indices and ISAllGather() is not scalable. */
2005: ISAllGather(nest->isglobal.col[j], &bNis);
2006: ISGetSize(bNis, &bN);
2007: ISGetIndices(bNis, &bNindices);
2008: for (i = 0; i < nest->nr; ++i) {
2009: PetscSF bmsf;
2010: PetscSFNode *iremote;
2011: Mat B = nest->m[i][j], D = NULL;
2012: PetscInt bm, *sub_dnnz, *sub_onnz, br;
2013: const PetscInt *bmindices;
2014: if (!B) continue;
2015: ISGetLocalSize(nest->isglobal.row[i], &bm);
2016: ISGetIndices(nest->isglobal.row[i], &bmindices);
2017: PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf);
2018: PetscMalloc1(bm, &iremote);
2019: PetscMalloc1(bm, &sub_dnnz);
2020: PetscMalloc1(bm, &sub_onnz);
2021: for (k = 0; k < bm; ++k) {
2022: sub_dnnz[k] = 0;
2023: sub_onnz[k] = 0;
2024: }
2025: PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &flg);
2026: if (flg) {
2027: PetscTryMethod(B, "MatTransposeGetMat_C", (Mat, Mat *), (B, &D));
2028: PetscTryMethod(B, "MatHermitianTransposeGetMat_C", (Mat, Mat *), (B, &D));
2029: MatConvert(B, ((PetscObject)D)->type_name, MAT_INITIAL_MATRIX, &D);
2030: B = D;
2031: }
2032: PetscObjectTypeCompareAny((PetscObject)B, &flg, MATSEQSBAIJ, MATMPISBAIJ, "");
2033: if (flg) {
2034: if (D) MatConvert(D, MATBAIJ, MAT_INPLACE_MATRIX, &D);
2035: else MatConvert(B, MATBAIJ, MAT_INITIAL_MATRIX, &D);
2036: B = D;
2037: }
2038: /*
2039: Locate the owners for all of the locally-owned global row indices for this row block.
2040: These determine the roots of PetscSF used to communicate preallocation data to row owners.
2041: The roots correspond to the dnnz and onnz entries; thus, there are two roots per row.
2042: */
2043: MatGetOwnershipRange(B, &rstart, NULL);
2044: for (br = 0; br < bm; ++br) {
2045: PetscInt row = bmindices[br], brncols, col;
2046: const PetscInt *brcols;
2047: PetscInt rowrel = 0; /* row's relative index on its owner rank */
2048: PetscMPIInt rowowner = 0;
2049: PetscLayoutFindOwnerIndex(A->rmap, row, &rowowner, &rowrel);
2050: /* how many roots */
2051: iremote[br].rank = rowowner;
2052: iremote[br].index = rowrel; /* edge from bmdnnz to dnnz */
2053: /* get nonzero pattern */
2054: MatGetRow(B, br + rstart, &brncols, &brcols, NULL);
2055: for (k = 0; k < brncols; k++) {
2056: col = bNindices[brcols[k]];
2057: if (col >= A->cmap->range[rowowner] && col < A->cmap->range[rowowner + 1]) {
2058: sub_dnnz[br]++;
2059: } else {
2060: sub_onnz[br]++;
2061: }
2062: }
2063: MatRestoreRow(B, br + rstart, &brncols, &brcols, NULL);
2064: }
2065: if (D) MatDestroy(&D);
2066: ISRestoreIndices(nest->isglobal.row[i], &bmindices);
2067: /* bsf will have to take care of disposing of bedges. */
2068: PetscSFSetGraph(bmsf, m, bm, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);
2069: PetscSFReduceBegin(bmsf, MPIU_INT, sub_dnnz, dnnz, MPI_SUM);
2070: PetscSFReduceEnd(bmsf, MPIU_INT, sub_dnnz, dnnz, MPI_SUM);
2071: PetscSFReduceBegin(bmsf, MPIU_INT, sub_onnz, onnz, MPI_SUM);
2072: PetscSFReduceEnd(bmsf, MPIU_INT, sub_onnz, onnz, MPI_SUM);
2073: PetscFree(sub_dnnz);
2074: PetscFree(sub_onnz);
2075: PetscSFDestroy(&bmsf);
2076: }
2077: ISRestoreIndices(bNis, &bNindices);
2078: ISDestroy(&bNis);
2079: }
2080: /* Resize preallocation if overestimated */
2081: for (i = 0; i < m; i++) {
2082: dnnz[i] = PetscMin(dnnz[i], A->cmap->n);
2083: onnz[i] = PetscMin(onnz[i], A->cmap->N - A->cmap->n);
2084: }
2085: MatSeqAIJSetPreallocation(C, 0, dnnz);
2086: MatMPIAIJSetPreallocation(C, 0, dnnz, 0, onnz);
2087: PetscFree(dnnz);
2088: MatAXPY_Dense_Nest(C, 1.0, A);
2089: if (reuse == MAT_INPLACE_MATRIX) {
2090: MatHeaderReplace(A, &C);
2091: } else *newmat = C;
2092: return 0;
2093: }
2095: PetscErrorCode MatConvert_Nest_Dense(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
2096: {
2097: Mat B;
2098: PetscInt m, n, M, N;
2100: MatGetSize(A, &M, &N);
2101: MatGetLocalSize(A, &m, &n);
2102: if (reuse == MAT_REUSE_MATRIX) {
2103: B = *newmat;
2104: MatZeroEntries(B);
2105: } else {
2106: MatCreateDense(PetscObjectComm((PetscObject)A), m, PETSC_DECIDE, M, N, NULL, &B);
2107: }
2108: MatAXPY_Dense_Nest(B, 1.0, A);
2109: if (reuse == MAT_INPLACE_MATRIX) {
2110: MatHeaderReplace(A, &B);
2111: } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
2112: return 0;
2113: }
2115: PetscErrorCode MatHasOperation_Nest(Mat mat, MatOperation op, PetscBool *has)
2116: {
2117: Mat_Nest *bA = (Mat_Nest *)mat->data;
2118: MatOperation opAdd;
2119: PetscInt i, j, nr = bA->nr, nc = bA->nc;
2120: PetscBool flg;
2122: *has = PETSC_FALSE;
2123: if (op == MATOP_MULT || op == MATOP_MULT_ADD || op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) {
2124: opAdd = (op == MATOP_MULT || op == MATOP_MULT_ADD ? MATOP_MULT_ADD : MATOP_MULT_TRANSPOSE_ADD);
2125: for (j = 0; j < nc; j++) {
2126: for (i = 0; i < nr; i++) {
2127: if (!bA->m[i][j]) continue;
2128: MatHasOperation(bA->m[i][j], opAdd, &flg);
2129: if (!flg) return 0;
2130: }
2131: }
2132: }
2133: if (((void **)mat->ops)[op]) *has = PETSC_TRUE;
2134: return 0;
2135: }
2137: /*MC
2138: MATNEST - MATNEST = "nest" - Matrix type consisting of nested submatrices, each stored separately.
2140: Level: intermediate
2142: Notes:
2143: This matrix type permits scalable use of `PCFIELDSPLIT` and avoids the large memory costs of extracting submatrices.
2144: It allows the use of symmetric and block formats for parts of multi-physics simulations.
2145: It is usually used with `DMCOMPOSITE` and `DMCreateMatrix()`
2147: Each of the submatrices lives on the same MPI communicator as the original nest matrix (though they can have zero
2148: rows/columns on some processes.) Thus this is not meant for cases where the submatrices live on far fewer processes
2149: than the nest matrix.
2151: .seealso: `MATNEST`, `MatCreate()`, `MatType`, `MatCreateNest()`, `MatNestSetSubMat()`, `MatNestGetSubMat()`,
2152: `VecCreateNest()`, `DMCreateMatrix()`, `DMCOMPOSITE`, `MatNestSetVecType()`, `MatNestGetLocalISs()`,
2153: `MatNestGetISs()`, `MatNestSetSubMats()`, `MatNestGetSubMats()`
2154: M*/
2155: PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A)
2156: {
2157: Mat_Nest *s;
2159: PetscNew(&s);
2160: A->data = (void *)s;
2162: s->nr = -1;
2163: s->nc = -1;
2164: s->m = NULL;
2165: s->splitassembly = PETSC_FALSE;
2167: PetscMemzero(A->ops, sizeof(*A->ops));
2169: A->ops->mult = MatMult_Nest;
2170: A->ops->multadd = MatMultAdd_Nest;
2171: A->ops->multtranspose = MatMultTranspose_Nest;
2172: A->ops->multtransposeadd = MatMultTransposeAdd_Nest;
2173: A->ops->transpose = MatTranspose_Nest;
2174: A->ops->assemblybegin = MatAssemblyBegin_Nest;
2175: A->ops->assemblyend = MatAssemblyEnd_Nest;
2176: A->ops->zeroentries = MatZeroEntries_Nest;
2177: A->ops->copy = MatCopy_Nest;
2178: A->ops->axpy = MatAXPY_Nest;
2179: A->ops->duplicate = MatDuplicate_Nest;
2180: A->ops->createsubmatrix = MatCreateSubMatrix_Nest;
2181: A->ops->destroy = MatDestroy_Nest;
2182: A->ops->view = MatView_Nest;
2183: A->ops->getvecs = NULL; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */
2184: A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_Nest;
2185: A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest;
2186: A->ops->getdiagonal = MatGetDiagonal_Nest;
2187: A->ops->diagonalscale = MatDiagonalScale_Nest;
2188: A->ops->scale = MatScale_Nest;
2189: A->ops->shift = MatShift_Nest;
2190: A->ops->diagonalset = MatDiagonalSet_Nest;
2191: A->ops->setrandom = MatSetRandom_Nest;
2192: A->ops->hasoperation = MatHasOperation_Nest;
2193: A->ops->missingdiagonal = MatMissingDiagonal_Nest;
2195: A->spptr = NULL;
2196: A->assembled = PETSC_FALSE;
2198: /* expose Nest api's */
2199: PetscObjectComposeFunction((PetscObject)A, "MatNestGetSubMat_C", MatNestGetSubMat_Nest);
2200: PetscObjectComposeFunction((PetscObject)A, "MatNestSetSubMat_C", MatNestSetSubMat_Nest);
2201: PetscObjectComposeFunction((PetscObject)A, "MatNestGetSubMats_C", MatNestGetSubMats_Nest);
2202: PetscObjectComposeFunction((PetscObject)A, "MatNestGetSize_C", MatNestGetSize_Nest);
2203: PetscObjectComposeFunction((PetscObject)A, "MatNestGetISs_C", MatNestGetISs_Nest);
2204: PetscObjectComposeFunction((PetscObject)A, "MatNestGetLocalISs_C", MatNestGetLocalISs_Nest);
2205: PetscObjectComposeFunction((PetscObject)A, "MatNestSetVecType_C", MatNestSetVecType_Nest);
2206: PetscObjectComposeFunction((PetscObject)A, "MatNestSetSubMats_C", MatNestSetSubMats_Nest);
2207: PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_mpiaij_C", MatConvert_Nest_AIJ);
2208: PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_seqaij_C", MatConvert_Nest_AIJ);
2209: PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_aij_C", MatConvert_Nest_AIJ);
2210: PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_is_C", MatConvert_Nest_IS);
2211: PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_mpidense_C", MatConvert_Nest_Dense);
2212: PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_seqdense_C", MatConvert_Nest_Dense);
2213: PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_seqdense_C", MatProductSetFromOptions_Nest_Dense);
2214: PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_mpidense_C", MatProductSetFromOptions_Nest_Dense);
2215: PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_dense_C", MatProductSetFromOptions_Nest_Dense);
2217: PetscObjectChangeTypeName((PetscObject)A, MATNEST);
2218: return 0;
2219: }