Actual source code: matis.c
1: /*
2: Creates a matrix class for using the Neumann-Neumann type preconditioners.
3: This stores the matrices in globally unassembled form. Each processor
4: assembles only its local Neumann problem and the parallel matrix vector
5: product is handled "implicitly".
7: Currently this allows for only one subdomain per processor.
8: */
10: #include <petsc/private/matisimpl.h>
11: #include <../src/mat/impls/aij/mpi/mpiaij.h>
12: #include <petsc/private/sfimpl.h>
13: #include <petsc/private/vecimpl.h>
14: #include <petsc/private/hashseti.h>
16: #define MATIS_MAX_ENTRIES_INSERTION 2048
18: /* copied from src/mat/impls/localref/mlocalref.c */
19: #define IndexSpaceGet(buf, nrow, ncol, irowm, icolm) \
20: do { \
21: if (nrow + ncol > (PetscInt)PETSC_STATIC_ARRAY_LENGTH(buf)) { \
22: PetscCall(PetscMalloc2(nrow, &irowm, ncol, &icolm)); \
23: } else { \
24: irowm = &buf[0]; \
25: icolm = &buf[nrow]; \
26: } \
27: } while (0)
29: #define IndexSpaceRestore(buf, nrow, ncol, irowm, icolm) \
30: do { \
31: if (nrow + ncol > (PetscInt)PETSC_STATIC_ARRAY_LENGTH(buf)) PetscCall(PetscFree2(irowm, icolm)); \
32: } while (0)
34: static void BlockIndicesExpand(PetscInt n, const PetscInt idx[], PetscInt bs, PetscInt idxm[])
35: {
36: for (PetscInt i = 0; i < n; i++) {
37: for (PetscInt j = 0; j < bs; j++) idxm[i * bs + j] = idx[i] * bs + j;
38: }
39: }
41: static PetscErrorCode MatSetValuesLocal_IS(Mat, PetscInt, const PetscInt *, PetscInt, const PetscInt *, const PetscScalar *, InsertMode);
42: static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat, PetscInt, const PetscInt *, PetscInt, const PetscInt *, const PetscScalar *, InsertMode);
43: static PetscErrorCode MatISSetUpScatters_Private(Mat);
45: static PetscErrorCode MatISContainerDestroyPtAP_Private(PetscCtxRt ptr)
46: {
47: MatISPtAP ptap = *(MatISPtAP *)ptr;
49: PetscFunctionBegin;
50: PetscCall(MatDestroySubMatrices(ptap->ris1 ? 2 : 1, &ptap->lP));
51: PetscCall(ISDestroy(&ptap->cis0));
52: PetscCall(ISDestroy(&ptap->cis1));
53: PetscCall(ISDestroy(&ptap->ris0));
54: PetscCall(ISDestroy(&ptap->ris1));
55: PetscCall(PetscFree(ptap));
56: PetscFunctionReturn(PETSC_SUCCESS);
57: }
59: static PetscErrorCode MatPtAPNumeric_IS_XAIJ(Mat A, Mat P, Mat C)
60: {
61: MatISPtAP ptap;
62: Mat_IS *matis = (Mat_IS *)A->data;
63: Mat lA, lC;
64: MatReuse reuse;
65: IS ris[2], cis[2];
66: PetscContainer c;
67: PetscInt n;
69: PetscFunctionBegin;
70: PetscCall(PetscObjectQuery((PetscObject)C, "_MatIS_PtAP", (PetscObject *)&c));
71: PetscCheck(c, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Missing PtAP information");
72: PetscCall(PetscContainerGetPointer(c, &ptap));
73: ris[0] = ptap->ris0;
74: ris[1] = ptap->ris1;
75: cis[0] = ptap->cis0;
76: cis[1] = ptap->cis1;
77: n = ptap->ris1 ? 2 : 1;
78: reuse = ptap->lP ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX;
79: PetscCall(MatCreateSubMatrices(P, n, ris, cis, reuse, &ptap->lP));
81: PetscCall(MatISGetLocalMat(A, &lA));
82: PetscCall(MatISGetLocalMat(C, &lC));
83: if (ptap->ris1) { /* unsymmetric A mapping */
84: Mat lPt;
86: PetscCall(MatTranspose(ptap->lP[1], MAT_INITIAL_MATRIX, &lPt));
87: PetscCall(MatMatMatMult(lPt, lA, ptap->lP[0], reuse, ptap->fill, &lC));
88: if (matis->storel2l) PetscCall(PetscObjectCompose((PetscObject)A, "_MatIS_PtAP_l2l", (PetscObject)lPt));
89: PetscCall(MatDestroy(&lPt));
90: } else {
91: PetscCall(MatPtAP(lA, ptap->lP[0], reuse, ptap->fill, &lC));
92: if (matis->storel2l) PetscCall(PetscObjectCompose((PetscObject)C, "_MatIS_PtAP_l2l", (PetscObject)ptap->lP[0]));
93: }
94: if (reuse == MAT_INITIAL_MATRIX) {
95: PetscCall(MatISSetLocalMat(C, lC));
96: PetscCall(MatDestroy(&lC));
97: }
98: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
99: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
100: PetscFunctionReturn(PETSC_SUCCESS);
101: }
103: static PetscErrorCode MatGetNonzeroColumnsLocal_Private(Mat PT, IS *cis)
104: {
105: Mat Po, Pd;
106: IS zd, zo;
107: const PetscInt *garray;
108: PetscInt *aux, i, bs;
109: PetscInt dc, stc, oc, ctd, cto;
110: PetscBool ismpiaij, ismpibaij, isseqaij, isseqbaij;
111: MPI_Comm comm;
113: PetscFunctionBegin;
115: PetscAssertPointer(cis, 2);
116: PetscCall(PetscObjectGetComm((PetscObject)PT, &comm));
117: bs = 1;
118: PetscCall(PetscObjectBaseTypeCompare((PetscObject)PT, MATMPIAIJ, &ismpiaij));
119: PetscCall(PetscObjectBaseTypeCompare((PetscObject)PT, MATMPIBAIJ, &ismpibaij));
120: PetscCall(PetscObjectBaseTypeCompare((PetscObject)PT, MATSEQAIJ, &isseqaij));
121: PetscCall(PetscObjectTypeCompare((PetscObject)PT, MATSEQBAIJ, &isseqbaij));
122: if (isseqaij || isseqbaij) {
123: Pd = PT;
124: Po = NULL;
125: garray = NULL;
126: } else if (ismpiaij) {
127: PetscCall(MatMPIAIJGetSeqAIJ(PT, &Pd, &Po, &garray));
128: } else if (ismpibaij) {
129: PetscCall(MatMPIBAIJGetSeqBAIJ(PT, &Pd, &Po, &garray));
130: PetscCall(MatGetBlockSize(PT, &bs));
131: } else SETERRQ(comm, PETSC_ERR_SUP, "Not for matrix type %s", ((PetscObject)PT)->type_name);
133: /* identify any null columns in Pd or Po */
134: /* We use a tolerance comparison since it may happen that, with geometric multigrid,
135: some of the columns are not really zero, but very close to */
136: zo = zd = NULL;
137: if (Po) PetscCall(MatFindNonzeroRowsOrCols_Basic(Po, PETSC_TRUE, PETSC_SMALL, &zo));
138: PetscCall(MatFindNonzeroRowsOrCols_Basic(Pd, PETSC_TRUE, PETSC_SMALL, &zd));
140: PetscCall(MatGetLocalSize(PT, NULL, &dc));
141: PetscCall(MatGetOwnershipRangeColumn(PT, &stc, NULL));
142: if (Po) PetscCall(MatGetLocalSize(Po, NULL, &oc));
143: else oc = 0;
144: PetscCall(PetscMalloc1((dc + oc) / bs, &aux));
145: if (zd) {
146: const PetscInt *idxs;
147: PetscInt nz;
149: /* this will throw an error if bs is not valid */
150: PetscCall(ISSetBlockSize(zd, bs));
151: PetscCall(ISGetLocalSize(zd, &nz));
152: PetscCall(ISGetIndices(zd, &idxs));
153: ctd = nz / bs;
154: for (i = 0; i < ctd; i++) aux[i] = (idxs[bs * i] + stc) / bs;
155: PetscCall(ISRestoreIndices(zd, &idxs));
156: } else {
157: ctd = dc / bs;
158: for (i = 0; i < ctd; i++) aux[i] = i + stc / bs;
159: }
160: if (zo) {
161: const PetscInt *idxs;
162: PetscInt nz;
164: /* this will throw an error if bs is not valid */
165: PetscCall(ISSetBlockSize(zo, bs));
166: PetscCall(ISGetLocalSize(zo, &nz));
167: PetscCall(ISGetIndices(zo, &idxs));
168: cto = nz / bs;
169: for (i = 0; i < cto; i++) aux[i + ctd] = garray[idxs[bs * i] / bs];
170: PetscCall(ISRestoreIndices(zo, &idxs));
171: } else {
172: cto = oc / bs;
173: for (i = 0; i < cto; i++) aux[i + ctd] = garray[i];
174: }
175: PetscCall(ISCreateBlock(comm, bs, ctd + cto, aux, PETSC_OWN_POINTER, cis));
176: PetscCall(ISDestroy(&zd));
177: PetscCall(ISDestroy(&zo));
178: PetscFunctionReturn(PETSC_SUCCESS);
179: }
181: static PetscErrorCode MatPtAPSymbolic_IS_XAIJ(Mat A, Mat P, PetscReal fill, Mat C)
182: {
183: Mat PT, lA;
184: MatISPtAP ptap;
185: ISLocalToGlobalMapping Crl2g, Ccl2g, rl2g, cl2g;
186: PetscContainer c;
187: MatType lmtype;
188: const PetscInt *garray;
189: PetscInt ibs, N, dc;
190: MPI_Comm comm;
192: PetscFunctionBegin;
193: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
194: PetscCall(MatSetType(C, MATIS));
195: PetscCall(MatISGetLocalMat(A, &lA));
196: PetscCall(MatGetType(lA, &lmtype));
197: PetscCall(MatISSetLocalMatType(C, lmtype));
198: PetscCall(MatGetSize(P, NULL, &N));
199: PetscCall(MatGetLocalSize(P, NULL, &dc));
200: PetscCall(MatSetSizes(C, dc, dc, N, N));
201: /* Not sure about this
202: PetscCall(MatGetBlockSizes(P,NULL,&ibs));
203: PetscCall(MatSetBlockSize(*C,ibs));
204: */
206: PetscCall(PetscNew(&ptap));
207: PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &c));
208: PetscCall(PetscContainerSetPointer(c, ptap));
209: PetscCall(PetscContainerSetCtxDestroy(c, MatISContainerDestroyPtAP_Private));
210: PetscCall(PetscObjectCompose((PetscObject)C, "_MatIS_PtAP", (PetscObject)c));
211: PetscCall(PetscContainerDestroy(&c));
212: ptap->fill = fill;
214: PetscCall(MatISGetLocalToGlobalMapping(A, &rl2g, &cl2g));
216: PetscCall(ISLocalToGlobalMappingGetBlockSize(cl2g, &ibs));
217: PetscCall(ISLocalToGlobalMappingGetSize(cl2g, &N));
218: PetscCall(ISLocalToGlobalMappingGetBlockIndices(cl2g, &garray));
219: PetscCall(ISCreateBlock(comm, ibs, N / ibs, garray, PETSC_COPY_VALUES, &ptap->ris0));
220: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(cl2g, &garray));
222: PetscCall(MatCreateSubMatrix(P, ptap->ris0, NULL, MAT_INITIAL_MATRIX, &PT));
223: PetscCall(MatGetNonzeroColumnsLocal_Private(PT, &ptap->cis0));
224: PetscCall(ISLocalToGlobalMappingCreateIS(ptap->cis0, &Ccl2g));
225: PetscCall(MatDestroy(&PT));
227: Crl2g = NULL;
228: if (rl2g != cl2g) { /* unsymmetric A mapping */
229: PetscBool same, lsame = PETSC_FALSE;
230: PetscInt N1, ibs1;
232: PetscCall(ISLocalToGlobalMappingGetSize(rl2g, &N1));
233: PetscCall(ISLocalToGlobalMappingGetBlockSize(rl2g, &ibs1));
234: PetscCall(ISLocalToGlobalMappingGetBlockIndices(rl2g, &garray));
235: PetscCall(ISCreateBlock(comm, ibs, N1 / ibs, garray, PETSC_COPY_VALUES, &ptap->ris1));
236: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(rl2g, &garray));
237: if (ibs1 == ibs && N1 == N) { /* check if the l2gmaps are the same */
238: const PetscInt *i1, *i2;
240: PetscCall(ISBlockGetIndices(ptap->ris0, &i1));
241: PetscCall(ISBlockGetIndices(ptap->ris1, &i2));
242: PetscCall(PetscArraycmp(i1, i2, N, &lsame));
243: }
244: PetscCallMPI(MPIU_Allreduce(&lsame, &same, 1, MPI_C_BOOL, MPI_LAND, comm));
245: if (same) {
246: PetscCall(ISDestroy(&ptap->ris1));
247: } else {
248: PetscCall(MatCreateSubMatrix(P, ptap->ris1, NULL, MAT_INITIAL_MATRIX, &PT));
249: PetscCall(MatGetNonzeroColumnsLocal_Private(PT, &ptap->cis1));
250: PetscCall(ISLocalToGlobalMappingCreateIS(ptap->cis1, &Crl2g));
251: PetscCall(MatDestroy(&PT));
252: }
253: }
254: /* Not sure about this
255: if (!Crl2g) {
256: PetscCall(MatGetBlockSize(C,&ibs));
257: PetscCall(ISLocalToGlobalMappingSetBlockSize(Ccl2g,ibs));
258: }
259: */
260: PetscCall(MatSetLocalToGlobalMapping(C, Crl2g ? Crl2g : Ccl2g, Ccl2g));
261: PetscCall(ISLocalToGlobalMappingDestroy(&Crl2g));
262: PetscCall(ISLocalToGlobalMappingDestroy(&Ccl2g));
264: C->ops->ptapnumeric = MatPtAPNumeric_IS_XAIJ;
265: PetscFunctionReturn(PETSC_SUCCESS);
266: }
268: static PetscErrorCode MatProductSymbolic_PtAP_IS_XAIJ(Mat C)
269: {
270: Mat_Product *product = C->product;
271: Mat A = product->A, P = product->B;
272: PetscReal fill = product->fill;
274: PetscFunctionBegin;
275: PetscCall(MatPtAPSymbolic_IS_XAIJ(A, P, fill, C));
276: C->ops->productnumeric = MatProductNumeric_PtAP;
277: PetscFunctionReturn(PETSC_SUCCESS);
278: }
280: static PetscErrorCode MatProductSetFromOptions_IS_XAIJ_PtAP(Mat C)
281: {
282: PetscFunctionBegin;
283: C->ops->productsymbolic = MatProductSymbolic_PtAP_IS_XAIJ;
284: PetscFunctionReturn(PETSC_SUCCESS);
285: }
287: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_IS_XAIJ(Mat C)
288: {
289: Mat_Product *product = C->product;
291: PetscFunctionBegin;
292: if (product->type == MATPRODUCT_PtAP) PetscCall(MatProductSetFromOptions_IS_XAIJ_PtAP(C));
293: PetscFunctionReturn(PETSC_SUCCESS);
294: }
296: static PetscErrorCode MatISContainerDestroyFields_Private(PetscCtxRt ptr)
297: {
298: MatISLocalFields lf = *(MatISLocalFields *)ptr;
299: PetscInt i;
301: PetscFunctionBegin;
302: for (i = 0; i < lf->nr; i++) PetscCall(ISDestroy(&lf->rf[i]));
303: for (i = 0; i < lf->nc; i++) PetscCall(ISDestroy(&lf->cf[i]));
304: PetscCall(PetscFree2(lf->rf, lf->cf));
305: PetscCall(PetscFree(lf));
306: PetscFunctionReturn(PETSC_SUCCESS);
307: }
309: static PetscErrorCode MatConvert_SeqXAIJ_IS(Mat A, MatType type, MatReuse reuse, Mat *newmat)
310: {
311: Mat B, lB;
313: PetscFunctionBegin;
314: if (reuse != MAT_REUSE_MATRIX) {
315: ISLocalToGlobalMapping rl2g, cl2g;
316: PetscInt bs;
317: IS is;
319: PetscCall(MatGetBlockSize(A, &bs));
320: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), A->rmap->n / bs, 0, 1, &is));
321: if (bs > 1) {
322: IS is2;
323: PetscInt i, *aux;
325: PetscCall(ISGetLocalSize(is, &i));
326: PetscCall(ISGetIndices(is, (const PetscInt **)&aux));
327: PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), bs, i, aux, PETSC_COPY_VALUES, &is2));
328: PetscCall(ISRestoreIndices(is, (const PetscInt **)&aux));
329: PetscCall(ISDestroy(&is));
330: is = is2;
331: }
332: PetscCall(ISSetIdentity(is));
333: PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g));
334: PetscCall(ISDestroy(&is));
335: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), A->cmap->n / bs, 0, 1, &is));
336: if (bs > 1) {
337: IS is2;
338: PetscInt i, *aux;
340: PetscCall(ISGetLocalSize(is, &i));
341: PetscCall(ISGetIndices(is, (const PetscInt **)&aux));
342: PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), bs, i, aux, PETSC_COPY_VALUES, &is2));
343: PetscCall(ISRestoreIndices(is, (const PetscInt **)&aux));
344: PetscCall(ISDestroy(&is));
345: is = is2;
346: }
347: PetscCall(ISSetIdentity(is));
348: PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g));
349: PetscCall(ISDestroy(&is));
350: PetscCall(MatCreateIS(PetscObjectComm((PetscObject)A), bs, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N, rl2g, cl2g, &B));
351: PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
352: PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
353: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &lB));
354: if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
355: } else {
356: B = *newmat;
357: PetscCall(PetscObjectReference((PetscObject)A));
358: lB = A;
359: }
360: PetscCall(MatISSetLocalMat(B, lB));
361: PetscCall(MatDestroy(&lB));
362: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
363: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
364: if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &B));
365: PetscFunctionReturn(PETSC_SUCCESS);
366: }
368: static PetscErrorCode MatISScaleDisassembling_Private(Mat A)
369: {
370: Mat_IS *matis = (Mat_IS *)A->data;
371: PetscScalar *aa;
372: const PetscInt *ii, *jj;
373: PetscInt i, n, m;
374: PetscInt *ecount, **eneighs;
375: PetscBool flg;
377: PetscFunctionBegin;
378: PetscCall(MatGetRowIJ(matis->A, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg));
379: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure");
380: PetscCall(ISLocalToGlobalMappingGetNodeInfo(matis->rmapping, &n, &ecount, &eneighs));
381: PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected %" PetscInt_FMT " != %" PetscInt_FMT, m, n);
382: PetscCall(MatSeqAIJGetArray(matis->A, &aa));
383: for (i = 0; i < n; i++) {
384: if (ecount[i] > 1) {
385: for (PetscInt j = ii[i]; j < ii[i + 1]; j++) {
386: PetscInt i2 = jj[j], p, p2;
387: PetscReal scal = 0.0;
389: for (p = 0; p < ecount[i]; p++) {
390: for (p2 = 0; p2 < ecount[i2]; p2++) {
391: if (eneighs[i][p] == eneighs[i2][p2]) {
392: scal += 1.0;
393: break;
394: }
395: }
396: }
397: if (scal) aa[j] /= scal;
398: }
399: }
400: }
401: PetscCall(ISLocalToGlobalMappingRestoreNodeInfo(matis->rmapping, &n, &ecount, &eneighs));
402: PetscCall(MatSeqAIJRestoreArray(matis->A, &aa));
403: PetscCall(MatRestoreRowIJ(matis->A, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg));
404: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot restore IJ structure");
405: PetscFunctionReturn(PETSC_SUCCESS);
406: }
408: typedef enum {
409: MAT_IS_DISASSEMBLE_L2G_NATURAL,
410: MAT_IS_DISASSEMBLE_L2G_MAT,
411: MAT_IS_DISASSEMBLE_L2G_ND
412: } MatISDisassemblel2gType;
414: static PetscErrorCode MatMPIXAIJComputeLocalToGlobalMapping_Private(Mat A, ISLocalToGlobalMapping *l2g)
415: {
416: Mat Ad, Ao;
417: IS is, ndmap, ndsub;
418: MPI_Comm comm;
419: const PetscInt *garray, *ndmapi;
420: PetscInt bs, i, cnt, nl, *ncount, *ndmapc;
421: PetscBool ismpiaij, ismpibaij;
422: const char *const MatISDisassemblel2gTypes[] = {"NATURAL", "MAT", "ND", "MatISDisassemblel2gType", "MAT_IS_DISASSEMBLE_L2G_", NULL};
423: MatISDisassemblel2gType mode = MAT_IS_DISASSEMBLE_L2G_NATURAL;
424: MatPartitioning part;
425: PetscSF sf;
426: PetscObject dm;
428: PetscFunctionBegin;
429: PetscOptionsBegin(PetscObjectComm((PetscObject)A), ((PetscObject)A)->prefix, "MatIS l2g disassembling options", "Mat");
430: PetscCall(PetscOptionsEnum("-mat_is_disassemble_l2g_type", "Type of local-to-global mapping to be used for disassembling", "MatISDisassemblel2gType", MatISDisassemblel2gTypes, (PetscEnum)mode, (PetscEnum *)&mode, NULL));
431: PetscOptionsEnd();
432: if (mode == MAT_IS_DISASSEMBLE_L2G_MAT) {
433: PetscCall(MatGetLocalToGlobalMapping(A, l2g, NULL));
434: PetscFunctionReturn(PETSC_SUCCESS);
435: }
436: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
437: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij));
438: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIBAIJ, &ismpibaij));
439: PetscCall(MatGetBlockSize(A, &bs));
440: switch (mode) {
441: case MAT_IS_DISASSEMBLE_L2G_ND:
442: PetscCall(MatPartitioningCreate(comm, &part));
443: PetscCall(MatPartitioningSetAdjacency(part, A));
444: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)part, ((PetscObject)A)->prefix));
445: PetscCall(MatPartitioningSetFromOptions(part));
446: PetscCall(MatPartitioningApplyND(part, &ndmap));
447: PetscCall(MatPartitioningDestroy(&part));
448: PetscCall(ISBuildTwoSided(ndmap, NULL, &ndsub));
449: PetscCall(MatMPIAIJSetUseScalableIncreaseOverlap(A, PETSC_TRUE));
450: PetscCall(MatIncreaseOverlap(A, 1, &ndsub, 1));
451: PetscCall(ISLocalToGlobalMappingCreateIS(ndsub, l2g));
453: /* it may happen that a separator node is not properly shared */
454: PetscCall(ISLocalToGlobalMappingGetNodeInfo(*l2g, &nl, &ncount, NULL));
455: PetscCall(PetscSFCreate(comm, &sf));
456: PetscCall(ISLocalToGlobalMappingGetIndices(*l2g, &garray));
457: PetscCall(PetscSFSetGraphLayout(sf, A->rmap, nl, NULL, PETSC_OWN_POINTER, garray));
458: PetscCall(ISLocalToGlobalMappingRestoreIndices(*l2g, &garray));
459: PetscCall(PetscCalloc1(A->rmap->n, &ndmapc));
460: PetscCall(PetscSFReduceBegin(sf, MPIU_INT, ncount, ndmapc, MPI_REPLACE));
461: PetscCall(PetscSFReduceEnd(sf, MPIU_INT, ncount, ndmapc, MPI_REPLACE));
462: PetscCall(ISLocalToGlobalMappingRestoreNodeInfo(*l2g, NULL, &ncount, NULL));
463: PetscCall(ISGetIndices(ndmap, &ndmapi));
464: for (i = 0, cnt = 0; i < A->rmap->n; i++)
465: if (ndmapi[i] < 0 && ndmapc[i] < 2) cnt++;
467: PetscCallMPI(MPIU_Allreduce(&cnt, &i, 1, MPIU_INT, MPI_MAX, comm));
468: if (i) { /* we detected isolated separator nodes */
469: Mat A2, A3;
470: IS *workis, is2;
471: PetscScalar *vals;
472: PetscInt gcnt = i, *dnz, *onz, j, *lndmapi;
473: ISLocalToGlobalMapping ll2g;
474: PetscBool flg;
475: const PetscInt *ii, *jj;
477: /* communicate global id of separators */
478: MatPreallocateBegin(comm, A->rmap->n, A->cmap->n, dnz, onz);
479: for (i = 0, cnt = 0; i < A->rmap->n; i++) dnz[i] = ndmapi[i] < 0 ? i + A->rmap->rstart : -1;
481: PetscCall(PetscMalloc1(nl, &lndmapi));
482: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, dnz, lndmapi, MPI_REPLACE));
484: /* compute adjacency of isolated separators node */
485: PetscCall(PetscMalloc1(gcnt, &workis));
486: for (i = 0, cnt = 0; i < A->rmap->n; i++) {
487: if (ndmapi[i] < 0 && ndmapc[i] < 2) PetscCall(ISCreateStride(comm, 1, i + A->rmap->rstart, 1, &workis[cnt++]));
488: }
489: for (i = cnt; i < gcnt; i++) PetscCall(ISCreateStride(comm, 0, 0, 1, &workis[i]));
490: for (i = 0; i < gcnt; i++) {
491: PetscCall(PetscObjectSetName((PetscObject)workis[i], "ISOLATED"));
492: PetscCall(ISViewFromOptions(workis[i], NULL, "-view_isolated_separators"));
493: }
495: /* no communications since all the ISes correspond to locally owned rows */
496: PetscCall(MatIncreaseOverlap(A, gcnt, workis, 1));
498: /* end communicate global id of separators */
499: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, dnz, lndmapi, MPI_REPLACE));
501: /* communicate new layers : create a matrix and transpose it */
502: PetscCall(PetscArrayzero(dnz, A->rmap->n));
503: PetscCall(PetscArrayzero(onz, A->rmap->n));
504: for (i = 0, j = 0; i < A->rmap->n; i++) {
505: if (ndmapi[i] < 0 && ndmapc[i] < 2) {
506: const PetscInt *idxs;
507: PetscInt s;
509: PetscCall(ISGetLocalSize(workis[j], &s));
510: PetscCall(ISGetIndices(workis[j], &idxs));
511: PetscCall(MatPreallocateSet(i + A->rmap->rstart, s, idxs, dnz, onz));
512: j++;
513: }
514: }
515: PetscCheck(j == cnt, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected local count %" PetscInt_FMT " != %" PetscInt_FMT, j, cnt);
517: for (i = 0; i < gcnt; i++) {
518: PetscCall(PetscObjectSetName((PetscObject)workis[i], "EXTENDED"));
519: PetscCall(ISViewFromOptions(workis[i], NULL, "-view_isolated_separators"));
520: }
522: for (i = 0, j = 0; i < A->rmap->n; i++) j = PetscMax(j, dnz[i] + onz[i]);
523: PetscCall(PetscMalloc1(j, &vals));
524: for (i = 0; i < j; i++) vals[i] = 1.0;
526: PetscCall(MatCreate(comm, &A2));
527: PetscCall(MatSetType(A2, MATMPIAIJ));
528: PetscCall(MatSetSizes(A2, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
529: PetscCall(MatMPIAIJSetPreallocation(A2, 0, dnz, 0, onz));
530: PetscCall(MatSetOption(A2, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
531: for (i = 0, j = 0; i < A2->rmap->n; i++) {
532: PetscInt row = i + A2->rmap->rstart, s = dnz[i] + onz[i];
533: const PetscInt *idxs;
535: if (s) {
536: PetscCall(ISGetIndices(workis[j], &idxs));
537: PetscCall(MatSetValues(A2, 1, &row, s, idxs, vals, INSERT_VALUES));
538: PetscCall(ISRestoreIndices(workis[j], &idxs));
539: j++;
540: }
541: }
542: PetscCheck(j == cnt, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected local count %" PetscInt_FMT " != %" PetscInt_FMT, j, cnt);
543: PetscCall(PetscFree(vals));
544: PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY));
545: PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY));
546: PetscCall(MatTranspose(A2, MAT_INPLACE_MATRIX, &A2));
548: /* extract submatrix corresponding to the coupling "owned separators" x "isolated separators" */
549: for (i = 0, j = 0; i < nl; i++)
550: if (lndmapi[i] >= 0) lndmapi[j++] = lndmapi[i];
551: PetscCall(ISCreateGeneral(comm, j, lndmapi, PETSC_USE_POINTER, &is));
552: PetscCall(MatMPIAIJGetLocalMatCondensed(A2, MAT_INITIAL_MATRIX, &is, NULL, &A3));
553: PetscCall(ISDestroy(&is));
554: PetscCall(MatDestroy(&A2));
556: /* extend local to global map to include connected isolated separators */
557: PetscCall(PetscObjectQuery((PetscObject)A3, "_petsc_GetLocalMatCondensed_iscol", (PetscObject *)&is));
558: PetscCheck(is, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing column map");
559: PetscCall(ISLocalToGlobalMappingCreateIS(is, &ll2g));
560: PetscCall(MatGetRowIJ(A3, 0, PETSC_FALSE, PETSC_FALSE, &i, &ii, &jj, &flg));
561: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure");
562: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ii[i], jj, PETSC_COPY_VALUES, &is));
563: PetscCall(MatRestoreRowIJ(A3, 0, PETSC_FALSE, PETSC_FALSE, &i, &ii, &jj, &flg));
564: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure");
565: PetscCall(ISLocalToGlobalMappingApplyIS(ll2g, is, &is2));
566: PetscCall(ISDestroy(&is));
567: PetscCall(ISLocalToGlobalMappingDestroy(&ll2g));
569: /* add new nodes to the local-to-global map */
570: PetscCall(ISLocalToGlobalMappingDestroy(l2g));
571: PetscCall(ISExpand(ndsub, is2, &is));
572: PetscCall(ISDestroy(&is2));
573: PetscCall(ISLocalToGlobalMappingCreateIS(is, l2g));
574: PetscCall(ISDestroy(&is));
576: PetscCall(MatDestroy(&A3));
577: PetscCall(PetscFree(lndmapi));
578: MatPreallocateEnd(dnz, onz);
579: for (i = 0; i < gcnt; i++) PetscCall(ISDestroy(&workis[i]));
580: PetscCall(PetscFree(workis));
581: }
582: PetscCall(ISRestoreIndices(ndmap, &ndmapi));
583: PetscCall(PetscSFDestroy(&sf));
584: PetscCall(PetscFree(ndmapc));
585: PetscCall(ISDestroy(&ndmap));
586: PetscCall(ISDestroy(&ndsub));
587: PetscCall(ISLocalToGlobalMappingSetBlockSize(*l2g, bs));
588: PetscCall(ISLocalToGlobalMappingViewFromOptions(*l2g, NULL, "-mat_is_nd_l2g_view"));
589: break;
590: case MAT_IS_DISASSEMBLE_L2G_NATURAL:
591: PetscCall(PetscObjectQuery((PetscObject)A, "__PETSc_dm", &dm));
592: if (dm) { /* if a matrix comes from a DM, most likely we can use the l2gmap if any */
593: PetscCall(MatGetLocalToGlobalMapping(A, l2g, NULL));
594: PetscCall(PetscObjectReference((PetscObject)*l2g));
595: if (*l2g) PetscFunctionReturn(PETSC_SUCCESS);
596: }
597: if (ismpiaij) {
598: PetscCall(MatMPIAIJGetSeqAIJ(A, &Ad, &Ao, &garray));
599: } else if (ismpibaij) {
600: PetscCall(MatMPIBAIJGetSeqBAIJ(A, &Ad, &Ao, &garray));
601: } else SETERRQ(comm, PETSC_ERR_SUP, "Type %s", ((PetscObject)A)->type_name);
602: if (A->rmap->n) {
603: PetscInt dc, oc, stc, *aux;
605: PetscCall(MatGetLocalSize(Ad, NULL, &dc));
606: PetscCall(MatGetLocalSize(Ao, NULL, &oc));
607: PetscCheck(!oc || garray, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "garray not present");
608: PetscCall(MatGetOwnershipRangeColumn(A, &stc, NULL));
609: PetscCall(PetscMalloc1((dc + oc) / bs, &aux));
610: for (i = 0; i < dc / bs; i++) aux[i] = i + stc / bs;
611: for (i = 0; i < oc / bs; i++) aux[i + dc / bs] = (ismpiaij ? garray[i * bs] / bs : garray[i]);
612: PetscCall(ISCreateBlock(comm, bs, (dc + oc) / bs, aux, PETSC_OWN_POINTER, &is));
613: } else {
614: PetscCall(ISCreateBlock(comm, 1, 0, NULL, PETSC_OWN_POINTER, &is));
615: }
616: PetscCall(ISLocalToGlobalMappingCreateIS(is, l2g));
617: PetscCall(ISDestroy(&is));
618: break;
619: default:
620: SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unsupported l2g disassembling type %d", mode);
621: }
622: PetscFunctionReturn(PETSC_SUCCESS);
623: }
625: PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat A, MatType type, MatReuse reuse, Mat *newmat)
626: {
627: Mat lA, Ad, Ao, B = NULL;
628: ISLocalToGlobalMapping rl2g, cl2g;
629: IS is;
630: MPI_Comm comm;
631: void *ptrs[2];
632: const char *names[2] = {"_convert_csr_aux", "_convert_csr_data"};
633: const PetscInt *garray;
634: PetscScalar *dd, *od, *aa, *data;
635: const PetscInt *di, *dj, *oi, *oj;
636: const PetscInt *odi, *odj, *ooi, *ooj;
637: PetscInt *aux, *ii, *jj;
638: PetscInt rbs, cbs, lc, dr, dc, oc, str, stc, nnz, i, jd, jo, cum;
639: PetscBool flg, ismpiaij, ismpibaij, was_inplace = PETSC_FALSE, cong;
640: PetscMPIInt size;
642: PetscFunctionBegin;
643: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
644: PetscCallMPI(MPI_Comm_size(comm, &size));
645: if (size == 1) {
646: PetscCall(MatConvert_SeqXAIJ_IS(A, type, reuse, newmat));
647: PetscFunctionReturn(PETSC_SUCCESS);
648: }
649: PetscCall(MatGetBlockSizes(A, &rbs, &cbs));
650: PetscCall(MatHasCongruentLayouts(A, &cong));
651: if (reuse != MAT_REUSE_MATRIX && cong && rbs == cbs) {
652: PetscCall(MatMPIXAIJComputeLocalToGlobalMapping_Private(A, &rl2g));
653: PetscCall(MatCreate(comm, &B));
654: PetscCall(MatSetType(B, MATIS));
655: PetscCall(MatSetSizes(B, A->rmap->n, A->rmap->n, A->rmap->N, A->rmap->N));
656: PetscCall(MatSetLocalToGlobalMapping(B, rl2g, rl2g));
657: PetscCall(MatSetBlockSizes(B, rbs, rbs));
658: PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
659: if (reuse == MAT_INPLACE_MATRIX) was_inplace = PETSC_TRUE;
660: reuse = MAT_REUSE_MATRIX;
661: }
662: if (reuse == MAT_REUSE_MATRIX) {
663: Mat *newlA, lA;
664: IS rows, cols;
665: const PetscInt *ridx, *cidx;
666: PetscInt nr, nc;
668: if (!B) B = *newmat;
669: PetscCall(MatISGetLocalToGlobalMapping(B, &rl2g, &cl2g));
670: PetscCall(ISLocalToGlobalMappingGetBlockIndices(rl2g, &ridx));
671: PetscCall(ISLocalToGlobalMappingGetBlockIndices(cl2g, &cidx));
672: PetscCall(ISLocalToGlobalMappingGetSize(rl2g, &nr));
673: PetscCall(ISLocalToGlobalMappingGetSize(cl2g, &nc));
674: PetscCall(ISLocalToGlobalMappingGetBlockSize(rl2g, &rbs));
675: PetscCall(ISLocalToGlobalMappingGetBlockSize(cl2g, &cbs));
676: PetscCall(ISCreateBlock(comm, rbs, nr / rbs, ridx, PETSC_USE_POINTER, &rows));
677: if (rl2g != cl2g) {
678: PetscCall(ISCreateBlock(comm, cbs, nc / cbs, cidx, PETSC_USE_POINTER, &cols));
679: } else {
680: PetscCall(PetscObjectReference((PetscObject)rows));
681: cols = rows;
682: }
683: PetscCall(MatISGetLocalMat(B, &lA));
684: PetscCall(MatCreateSubMatrices(A, 1, &rows, &cols, MAT_INITIAL_MATRIX, &newlA));
685: PetscCall(MatConvert(newlA[0], MATSEQAIJ, MAT_INPLACE_MATRIX, &newlA[0]));
686: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(rl2g, &ridx));
687: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(cl2g, &cidx));
688: PetscCall(ISDestroy(&rows));
689: PetscCall(ISDestroy(&cols));
690: if (!lA->preallocated) { /* first time */
691: PetscCall(MatDuplicate(newlA[0], MAT_COPY_VALUES, &lA));
692: PetscCall(MatISSetLocalMat(B, lA));
693: PetscCall(PetscObjectDereference((PetscObject)lA));
694: }
695: PetscCall(MatCopy(newlA[0], lA, SAME_NONZERO_PATTERN));
696: PetscCall(MatDestroySubMatrices(1, &newlA));
697: PetscCall(MatISScaleDisassembling_Private(B));
698: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
699: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
700: if (was_inplace) PetscCall(MatHeaderReplace(A, &B));
701: else *newmat = B;
702: PetscFunctionReturn(PETSC_SUCCESS);
703: }
704: /* general case, just compress out the column space */
705: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij));
706: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIBAIJ, &ismpibaij));
707: if (ismpiaij) {
708: cbs = 1; /* We cannot guarantee the off-process matrix will respect the column block size */
709: PetscCall(MatMPIAIJGetSeqAIJ(A, &Ad, &Ao, &garray));
710: } else if (ismpibaij) {
711: PetscCall(MatMPIBAIJGetSeqBAIJ(A, &Ad, &Ao, &garray));
712: PetscCall(MatConvert(Ad, MATSEQAIJ, MAT_INITIAL_MATRIX, &Ad));
713: PetscCall(MatConvert(Ao, MATSEQAIJ, MAT_INITIAL_MATRIX, &Ao));
714: } else SETERRQ(comm, PETSC_ERR_SUP, "Type %s", ((PetscObject)A)->type_name);
715: PetscCall(MatSeqAIJGetArray(Ad, &dd));
716: PetscCall(MatSeqAIJGetArray(Ao, &od));
718: /* access relevant information from MPIAIJ */
719: PetscCall(MatGetOwnershipRange(A, &str, NULL));
720: PetscCall(MatGetOwnershipRangeColumn(A, &stc, NULL));
721: PetscCall(MatGetLocalSize(Ad, &dr, &dc));
722: PetscCall(MatGetLocalSize(Ao, NULL, &oc));
723: PetscCheck(!oc || garray, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "garray not present");
725: PetscCall(MatGetRowIJ(Ad, 0, PETSC_FALSE, PETSC_FALSE, &i, &di, &dj, &flg));
726: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure");
727: PetscCall(MatGetRowIJ(Ao, 0, PETSC_FALSE, PETSC_FALSE, &i, &oi, &oj, &flg));
728: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure");
729: nnz = di[dr] + oi[dr];
730: /* store original pointers to be restored later */
731: odi = di;
732: odj = dj;
733: ooi = oi;
734: ooj = oj;
736: /* generate l2g maps for rows and cols */
737: PetscCall(ISCreateStride(comm, dr / rbs, str / rbs, 1, &is));
738: if (rbs > 1) {
739: IS is2;
741: PetscCall(ISGetLocalSize(is, &i));
742: PetscCall(ISGetIndices(is, (const PetscInt **)&aux));
743: PetscCall(ISCreateBlock(comm, rbs, i, aux, PETSC_COPY_VALUES, &is2));
744: PetscCall(ISRestoreIndices(is, (const PetscInt **)&aux));
745: PetscCall(ISDestroy(&is));
746: is = is2;
747: }
748: PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g));
749: PetscCall(ISDestroy(&is));
750: if (dr) {
751: PetscCall(PetscMalloc1((dc + oc) / cbs, &aux));
752: for (i = 0; i < dc / cbs; i++) aux[i] = i + stc / cbs;
753: for (i = 0; i < oc / cbs; i++) aux[i + dc / cbs] = garray[i];
754: PetscCall(ISCreateBlock(comm, cbs, (dc + oc) / cbs, aux, PETSC_OWN_POINTER, &is));
755: lc = dc + oc;
756: } else {
757: PetscCall(ISCreateBlock(comm, cbs, 0, NULL, PETSC_OWN_POINTER, &is));
758: lc = 0;
759: }
760: PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g));
761: PetscCall(ISDestroy(&is));
763: /* create MATIS object */
764: PetscCall(MatCreate(comm, &B));
765: PetscCall(MatSetSizes(B, dr, dc, PETSC_DECIDE, PETSC_DECIDE));
766: PetscCall(MatSetType(B, MATIS));
767: PetscCall(MatSetBlockSizes(B, rbs, cbs));
768: PetscCall(MatSetLocalToGlobalMapping(B, rl2g, cl2g));
769: PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
770: PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
772: /* merge local matrices */
773: PetscCall(PetscMalloc1(nnz + dr + 1, &aux));
774: PetscCall(PetscMalloc1(nnz, &data));
775: ii = aux;
776: jj = aux + dr + 1;
777: aa = data;
778: *ii = *(di++) + *(oi++);
779: for (jd = 0, jo = 0, cum = 0; *ii < nnz; cum++) {
780: for (; jd < *di; jd++) {
781: *jj++ = *dj++;
782: *aa++ = *dd++;
783: }
784: for (; jo < *oi; jo++) {
785: *jj++ = *oj++ + dc;
786: *aa++ = *od++;
787: }
788: *(++ii) = *(di++) + *(oi++);
789: }
790: for (; cum < dr; cum++) *(++ii) = nnz;
792: PetscCall(MatRestoreRowIJ(Ad, 0, PETSC_FALSE, PETSC_FALSE, &i, &odi, &odj, &flg));
793: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot restore IJ structure");
794: PetscCall(MatRestoreRowIJ(Ao, 0, PETSC_FALSE, PETSC_FALSE, &i, &ooi, &ooj, &flg));
795: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot restore IJ structure");
796: PetscCall(MatSeqAIJRestoreArray(Ad, &dd));
797: PetscCall(MatSeqAIJRestoreArray(Ao, &od));
799: ii = aux;
800: jj = aux + dr + 1;
801: aa = data;
802: PetscCall(MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, dr, lc, ii, jj, aa, &lA));
804: /* create containers to destroy the data */
805: ptrs[0] = aux;
806: ptrs[1] = data;
807: for (i = 0; i < 2; i++) PetscCall(PetscObjectContainerCompose((PetscObject)lA, names[i], ptrs[i], PetscCtxDestroyDefault));
808: if (ismpibaij) { /* destroy converted local matrices */
809: PetscCall(MatDestroy(&Ad));
810: PetscCall(MatDestroy(&Ao));
811: }
813: /* finalize matrix */
814: PetscCall(MatISSetLocalMat(B, lA));
815: PetscCall(MatDestroy(&lA));
816: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
817: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
818: if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &B));
819: else *newmat = B;
820: PetscFunctionReturn(PETSC_SUCCESS);
821: }
823: PETSC_INTERN PetscErrorCode MatConvert_Nest_IS(Mat A, MatType type, MatReuse reuse, Mat *newmat)
824: {
825: Mat **nest, *snest, **rnest, lA, B;
826: IS *iscol, *isrow, *islrow, *islcol;
827: ISLocalToGlobalMapping rl2g, cl2g;
828: MPI_Comm comm;
829: PetscInt *lr, *lc, *l2gidxs;
830: PetscInt i, j, nr, nc, rbs, cbs;
831: PetscBool convert, lreuse, *istrans;
832: PetscBool3 allow_repeated = PETSC_BOOL3_UNKNOWN;
834: PetscFunctionBegin;
835: PetscCall(MatNestGetSubMats(A, &nr, &nc, &nest));
836: lreuse = PETSC_FALSE;
837: rnest = NULL;
838: if (reuse == MAT_REUSE_MATRIX) {
839: PetscBool ismatis, isnest;
841: PetscCall(PetscObjectTypeCompare((PetscObject)*newmat, MATIS, &ismatis));
842: PetscCheck(ismatis, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_USER, "Cannot reuse matrix of type %s", ((PetscObject)*newmat)->type_name);
843: PetscCall(MatISGetLocalMat(*newmat, &lA));
844: PetscCall(PetscObjectTypeCompare((PetscObject)lA, MATNEST, &isnest));
845: if (isnest) {
846: PetscCall(MatNestGetSubMats(lA, &i, &j, &rnest));
847: lreuse = (PetscBool)(i == nr && j == nc);
848: if (!lreuse) rnest = NULL;
849: }
850: }
851: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
852: PetscCall(PetscCalloc2(nr, &lr, nc, &lc));
853: PetscCall(PetscCalloc6(nr, &isrow, nc, &iscol, nr, &islrow, nc, &islcol, nr * nc, &snest, nr * nc, &istrans));
854: PetscCall(MatNestGetISs(A, isrow, iscol));
855: for (i = 0; i < nr; i++) {
856: for (j = 0; j < nc; j++) {
857: PetscBool ismatis, sallow;
858: PetscInt l1, l2, lb1, lb2, ij = i * nc + j;
860: /* Null matrix pointers are allowed in MATNEST */
861: if (!nest[i][j]) continue;
863: /* Nested matrices should be of type MATIS */
864: PetscCall(PetscObjectTypeCompare((PetscObject)nest[i][j], MATTRANSPOSEVIRTUAL, &istrans[ij]));
865: if (istrans[ij]) {
866: Mat T, lT;
867: PetscCall(MatTransposeGetMat(nest[i][j], &T));
868: PetscCall(PetscObjectTypeCompare((PetscObject)T, MATIS, &ismatis));
869: PetscCheck(ismatis, comm, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") (transposed) is not of type MATIS", i, j);
870: PetscCall(MatISGetAllowRepeated(T, &sallow));
871: PetscCall(MatISGetLocalMat(T, &lT));
872: PetscCall(MatCreateTranspose(lT, &snest[ij]));
873: } else {
874: PetscCall(PetscObjectTypeCompare((PetscObject)nest[i][j], MATIS, &ismatis));
875: PetscCheck(ismatis, comm, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") is not of type MATIS", i, j);
876: PetscCall(MatISGetAllowRepeated(nest[i][j], &sallow));
877: PetscCall(MatISGetLocalMat(nest[i][j], &snest[ij]));
878: }
879: if (allow_repeated == PETSC_BOOL3_UNKNOWN) allow_repeated = PetscBoolToBool3(sallow);
880: PetscCheck(sallow == PetscBool3ToBool(allow_repeated), comm, PETSC_ERR_SUP, "Cannot mix repeated and non repeated maps");
882: /* Check compatibility of local sizes */
883: PetscCall(MatGetSize(snest[ij], &l1, &l2));
884: PetscCall(MatGetBlockSizes(snest[ij], &lb1, &lb2));
885: if (!l1 || !l2) continue;
886: PetscCheck(!lr[i] || l1 == lr[i], PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid local size %" PetscInt_FMT " != %" PetscInt_FMT, i, j, lr[i], l1);
887: PetscCheck(!lc[j] || l2 == lc[j], PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid local size %" PetscInt_FMT " != %" PetscInt_FMT, i, j, lc[j], l2);
888: lr[i] = l1;
889: lc[j] = l2;
891: /* check compatibility for local matrix reusage */
892: if (rnest && !rnest[i][j] != !snest[ij]) lreuse = PETSC_FALSE;
893: }
894: }
896: if (PetscDefined(USE_DEBUG)) {
897: /* Check compatibility of l2g maps for rows */
898: for (i = 0; i < nr; i++) {
899: rl2g = NULL;
900: for (j = 0; j < nc; j++) {
901: PetscInt n1, n2;
903: if (!nest[i][j]) continue;
904: if (istrans[i * nc + j]) {
905: Mat T;
907: PetscCall(MatTransposeGetMat(nest[i][j], &T));
908: PetscCall(MatISGetLocalToGlobalMapping(T, NULL, &cl2g));
909: } else {
910: PetscCall(MatISGetLocalToGlobalMapping(nest[i][j], &cl2g, NULL));
911: }
912: PetscCall(ISLocalToGlobalMappingGetSize(cl2g, &n1));
913: if (!n1) continue;
914: if (!rl2g) {
915: rl2g = cl2g;
916: } else {
917: const PetscInt *idxs1, *idxs2;
918: PetscBool same;
920: PetscCall(ISLocalToGlobalMappingGetSize(rl2g, &n2));
921: PetscCheck(n1 == n2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid row l2gmap size %" PetscInt_FMT " != %" PetscInt_FMT, i, j, n1, n2);
922: PetscCall(ISLocalToGlobalMappingGetIndices(cl2g, &idxs1));
923: PetscCall(ISLocalToGlobalMappingGetIndices(rl2g, &idxs2));
924: PetscCall(PetscArraycmp(idxs1, idxs2, n1, &same));
925: PetscCall(ISLocalToGlobalMappingRestoreIndices(cl2g, &idxs1));
926: PetscCall(ISLocalToGlobalMappingRestoreIndices(rl2g, &idxs2));
927: PetscCheck(same, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid row l2gmap", i, j);
928: }
929: }
930: }
931: /* Check compatibility of l2g maps for columns */
932: for (i = 0; i < nc; i++) {
933: rl2g = NULL;
934: for (j = 0; j < nr; j++) {
935: PetscInt n1, n2;
937: if (!nest[j][i]) continue;
938: if (istrans[j * nc + i]) {
939: Mat T;
941: PetscCall(MatTransposeGetMat(nest[j][i], &T));
942: PetscCall(MatISGetLocalToGlobalMapping(T, &cl2g, NULL));
943: } else {
944: PetscCall(MatISGetLocalToGlobalMapping(nest[j][i], NULL, &cl2g));
945: }
946: PetscCall(ISLocalToGlobalMappingGetSize(cl2g, &n1));
947: if (!n1) continue;
948: if (!rl2g) {
949: rl2g = cl2g;
950: } else {
951: const PetscInt *idxs1, *idxs2;
952: PetscBool same;
954: PetscCall(ISLocalToGlobalMappingGetSize(rl2g, &n2));
955: PetscCheck(n1 == n2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid column l2gmap size %" PetscInt_FMT " != %" PetscInt_FMT, j, i, n1, n2);
956: PetscCall(ISLocalToGlobalMappingGetIndices(cl2g, &idxs1));
957: PetscCall(ISLocalToGlobalMappingGetIndices(rl2g, &idxs2));
958: PetscCall(PetscArraycmp(idxs1, idxs2, n1, &same));
959: PetscCall(ISLocalToGlobalMappingRestoreIndices(cl2g, &idxs1));
960: PetscCall(ISLocalToGlobalMappingRestoreIndices(rl2g, &idxs2));
961: PetscCheck(same, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid column l2gmap", j, i);
962: }
963: }
964: }
965: }
967: B = NULL;
968: if (reuse != MAT_REUSE_MATRIX) {
969: PetscInt stl;
971: /* Create l2g map for the rows of the new matrix and index sets for the local MATNEST */
972: for (i = 0, stl = 0; i < nr; i++) stl += lr[i];
973: PetscCall(PetscMalloc1(stl, &l2gidxs));
974: for (i = 0, stl = 0; i < nr; i++) {
975: Mat usedmat;
976: Mat_IS *matis;
977: const PetscInt *idxs;
979: /* local IS for local NEST */
980: PetscCall(ISCreateStride(PETSC_COMM_SELF, lr[i], stl, 1, &islrow[i]));
982: /* l2gmap */
983: j = 0;
984: usedmat = nest[i][j];
985: while (!usedmat && j < nc - 1) usedmat = nest[i][++j];
986: PetscCheck(usedmat, comm, PETSC_ERR_SUP, "Cannot find valid row mat");
988: if (istrans[i * nc + j]) {
989: Mat T;
990: PetscCall(MatTransposeGetMat(usedmat, &T));
991: usedmat = T;
992: }
993: matis = (Mat_IS *)usedmat->data;
994: PetscCall(ISGetIndices(isrow[i], &idxs));
995: if (istrans[i * nc + j]) {
996: PetscCall(PetscSFBcastBegin(matis->csf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
997: PetscCall(PetscSFBcastEnd(matis->csf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
998: } else {
999: PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
1000: PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
1001: }
1002: PetscCall(ISRestoreIndices(isrow[i], &idxs));
1003: stl += lr[i];
1004: }
1005: PetscCall(ISLocalToGlobalMappingCreate(comm, 1, stl, l2gidxs, PETSC_OWN_POINTER, &rl2g));
1007: /* Create l2g map for columns of the new matrix and index sets for the local MATNEST */
1008: for (i = 0, stl = 0; i < nc; i++) stl += lc[i];
1009: PetscCall(PetscMalloc1(stl, &l2gidxs));
1010: for (i = 0, stl = 0; i < nc; i++) {
1011: Mat usedmat;
1012: Mat_IS *matis;
1013: const PetscInt *idxs;
1015: /* local IS for local NEST */
1016: PetscCall(ISCreateStride(PETSC_COMM_SELF, lc[i], stl, 1, &islcol[i]));
1018: /* l2gmap */
1019: j = 0;
1020: usedmat = nest[j][i];
1021: while (!usedmat && j < nr - 1) usedmat = nest[++j][i];
1022: PetscCheck(usedmat, comm, PETSC_ERR_SUP, "Cannot find valid column mat");
1023: if (istrans[j * nc + i]) {
1024: Mat T;
1025: PetscCall(MatTransposeGetMat(usedmat, &T));
1026: usedmat = T;
1027: }
1028: matis = (Mat_IS *)usedmat->data;
1029: PetscCall(ISGetIndices(iscol[i], &idxs));
1030: if (istrans[j * nc + i]) {
1031: PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
1032: PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
1033: } else {
1034: PetscCall(PetscSFBcastBegin(matis->csf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
1035: PetscCall(PetscSFBcastEnd(matis->csf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
1036: }
1037: PetscCall(ISRestoreIndices(iscol[i], &idxs));
1038: stl += lc[i];
1039: }
1040: PetscCall(ISLocalToGlobalMappingCreate(comm, 1, stl, l2gidxs, PETSC_OWN_POINTER, &cl2g));
1042: /* Create MATIS */
1043: PetscCall(MatCreate(comm, &B));
1044: PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
1045: PetscCall(MatGetBlockSizes(A, &rbs, &cbs));
1046: PetscCall(MatSetBlockSizes(B, rbs, cbs));
1047: PetscCall(MatSetType(B, MATIS));
1048: PetscCall(MatISSetLocalMatType(B, MATNEST));
1049: PetscCall(MatISSetAllowRepeated(B, PetscBool3ToBool(allow_repeated)));
1050: { /* hack : avoid setup of scatters */
1051: Mat_IS *matis = (Mat_IS *)B->data;
1052: matis->islocalref = B;
1053: }
1054: PetscCall(MatSetLocalToGlobalMapping(B, rl2g, cl2g));
1055: PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
1056: PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
1057: PetscCall(MatCreateNest(PETSC_COMM_SELF, nr, islrow, nc, islcol, snest, &lA));
1058: PetscCall(MatNestSetVecType(lA, VECNEST));
1059: for (i = 0; i < nr * nc; i++) {
1060: if (istrans[i]) PetscCall(MatDestroy(&snest[i]));
1061: }
1062: PetscCall(MatISSetLocalMat(B, lA));
1063: PetscCall(MatDestroy(&lA));
1064: { /* hack : setup of scatters done here */
1065: Mat_IS *matis = (Mat_IS *)B->data;
1067: matis->islocalref = NULL;
1068: PetscCall(MatISSetUpScatters_Private(B));
1069: }
1070: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
1071: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
1072: if (reuse == MAT_INPLACE_MATRIX) {
1073: PetscCall(MatHeaderReplace(A, &B));
1074: } else {
1075: *newmat = B;
1076: }
1077: } else {
1078: if (lreuse) {
1079: PetscCall(MatISGetLocalMat(*newmat, &lA));
1080: for (i = 0; i < nr; i++) {
1081: for (j = 0; j < nc; j++) {
1082: if (snest[i * nc + j]) {
1083: PetscCall(MatNestSetSubMat(lA, i, j, snest[i * nc + j]));
1084: if (istrans[i * nc + j]) PetscCall(MatDestroy(&snest[i * nc + j]));
1085: }
1086: }
1087: }
1088: } else {
1089: PetscInt stl;
1090: for (i = 0, stl = 0; i < nr; i++) {
1091: PetscCall(ISCreateStride(PETSC_COMM_SELF, lr[i], stl, 1, &islrow[i]));
1092: stl += lr[i];
1093: }
1094: for (i = 0, stl = 0; i < nc; i++) {
1095: PetscCall(ISCreateStride(PETSC_COMM_SELF, lc[i], stl, 1, &islcol[i]));
1096: stl += lc[i];
1097: }
1098: PetscCall(MatCreateNest(PETSC_COMM_SELF, nr, islrow, nc, islcol, snest, &lA));
1099: for (i = 0; i < nr * nc; i++) {
1100: if (istrans[i]) PetscCall(MatDestroy(&snest[i]));
1101: }
1102: PetscCall(MatISSetLocalMat(*newmat, lA));
1103: PetscCall(MatDestroy(&lA));
1104: }
1105: PetscCall(MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY));
1106: PetscCall(MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY));
1107: }
1109: /* Create local matrix in MATNEST format */
1110: convert = PETSC_FALSE;
1111: PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)A)->prefix, "-mat_is_convert_local_nest", &convert, NULL));
1112: if (convert) {
1113: Mat M;
1114: MatISLocalFields lf;
1116: PetscCall(MatISGetLocalMat(*newmat, &lA));
1117: PetscCall(MatConvert(lA, MATAIJ, MAT_INITIAL_MATRIX, &M));
1118: PetscCall(MatISSetLocalMat(*newmat, M));
1119: PetscCall(MatDestroy(&M));
1121: /* attach local fields to the matrix */
1122: PetscCall(PetscNew(&lf));
1123: PetscCall(PetscMalloc2(nr, &lf->rf, nc, &lf->cf));
1124: for (i = 0; i < nr; i++) {
1125: PetscInt n, st;
1127: PetscCall(ISGetLocalSize(islrow[i], &n));
1128: PetscCall(ISStrideGetInfo(islrow[i], &st, NULL));
1129: PetscCall(ISCreateStride(comm, n, st, 1, &lf->rf[i]));
1130: }
1131: for (i = 0; i < nc; i++) {
1132: PetscInt n, st;
1134: PetscCall(ISGetLocalSize(islcol[i], &n));
1135: PetscCall(ISStrideGetInfo(islcol[i], &st, NULL));
1136: PetscCall(ISCreateStride(comm, n, st, 1, &lf->cf[i]));
1137: }
1138: lf->nr = nr;
1139: lf->nc = nc;
1140: PetscCall(PetscObjectContainerCompose((PetscObject)*newmat, "_convert_nest_lfields", lf, MatISContainerDestroyFields_Private));
1141: }
1143: /* Free workspace */
1144: for (i = 0; i < nr; i++) PetscCall(ISDestroy(&islrow[i]));
1145: for (i = 0; i < nc; i++) PetscCall(ISDestroy(&islcol[i]));
1146: PetscCall(PetscFree6(isrow, iscol, islrow, islcol, snest, istrans));
1147: PetscCall(PetscFree2(lr, lc));
1148: PetscFunctionReturn(PETSC_SUCCESS);
1149: }
1151: static PetscErrorCode MatDiagonalScale_IS(Mat A, Vec l, Vec r)
1152: {
1153: Mat_IS *matis = (Mat_IS *)A->data;
1154: Vec ll, rr;
1155: const PetscScalar *Y, *X;
1156: PetscScalar *x, *y;
1158: PetscFunctionBegin;
1159: if (l) {
1160: ll = matis->y;
1161: PetscCall(VecGetArrayRead(l, &Y));
1162: PetscCall(VecGetArray(ll, &y));
1163: PetscCall(PetscSFBcastBegin(matis->sf, MPIU_SCALAR, Y, y, MPI_REPLACE));
1164: } else {
1165: ll = NULL;
1166: }
1167: if (r) {
1168: rr = matis->x;
1169: PetscCall(VecGetArrayRead(r, &X));
1170: PetscCall(VecGetArray(rr, &x));
1171: PetscCall(PetscSFBcastBegin(matis->csf, MPIU_SCALAR, X, x, MPI_REPLACE));
1172: } else {
1173: rr = NULL;
1174: }
1175: if (ll) {
1176: PetscCall(PetscSFBcastEnd(matis->sf, MPIU_SCALAR, Y, y, MPI_REPLACE));
1177: PetscCall(VecRestoreArrayRead(l, &Y));
1178: PetscCall(VecRestoreArray(ll, &y));
1179: }
1180: if (rr) {
1181: PetscCall(PetscSFBcastEnd(matis->csf, MPIU_SCALAR, X, x, MPI_REPLACE));
1182: PetscCall(VecRestoreArrayRead(r, &X));
1183: PetscCall(VecRestoreArray(rr, &x));
1184: }
1185: PetscCall(MatDiagonalScale(matis->A, ll, rr));
1186: PetscFunctionReturn(PETSC_SUCCESS);
1187: }
1189: static PetscErrorCode MatGetInfo_IS(Mat A, MatInfoType flag, MatInfo *ginfo)
1190: {
1191: Mat_IS *matis = (Mat_IS *)A->data;
1192: MatInfo info;
1193: PetscLogDouble isend[6], irecv[6];
1194: PetscInt bs;
1196: PetscFunctionBegin;
1197: PetscCall(MatGetBlockSize(A, &bs));
1198: if (matis->A->ops->getinfo) {
1199: PetscCall(MatGetInfo(matis->A, MAT_LOCAL, &info));
1200: isend[0] = info.nz_used;
1201: isend[1] = info.nz_allocated;
1202: isend[2] = info.nz_unneeded;
1203: isend[3] = info.memory;
1204: isend[4] = info.mallocs;
1205: } else {
1206: isend[0] = 0.;
1207: isend[1] = 0.;
1208: isend[2] = 0.;
1209: isend[3] = 0.;
1210: isend[4] = 0.;
1211: }
1212: isend[5] = matis->A->num_ass;
1213: if (flag == MAT_LOCAL) {
1214: ginfo->nz_used = isend[0];
1215: ginfo->nz_allocated = isend[1];
1216: ginfo->nz_unneeded = isend[2];
1217: ginfo->memory = isend[3];
1218: ginfo->mallocs = isend[4];
1219: ginfo->assemblies = isend[5];
1220: } else if (flag == MAT_GLOBAL_MAX) {
1221: PetscCallMPI(MPIU_Allreduce(isend, irecv, 6, MPIU_PETSCLOGDOUBLE, MPI_MAX, PetscObjectComm((PetscObject)A)));
1223: ginfo->nz_used = irecv[0];
1224: ginfo->nz_allocated = irecv[1];
1225: ginfo->nz_unneeded = irecv[2];
1226: ginfo->memory = irecv[3];
1227: ginfo->mallocs = irecv[4];
1228: ginfo->assemblies = irecv[5];
1229: } else if (flag == MAT_GLOBAL_SUM) {
1230: PetscCallMPI(MPIU_Allreduce(isend, irecv, 5, MPIU_PETSCLOGDOUBLE, MPI_SUM, PetscObjectComm((PetscObject)A)));
1232: ginfo->nz_used = irecv[0];
1233: ginfo->nz_allocated = irecv[1];
1234: ginfo->nz_unneeded = irecv[2];
1235: ginfo->memory = irecv[3];
1236: ginfo->mallocs = irecv[4];
1237: ginfo->assemblies = A->num_ass;
1238: }
1239: ginfo->block_size = bs;
1240: ginfo->fill_ratio_given = 0;
1241: ginfo->fill_ratio_needed = 0;
1242: ginfo->factor_mallocs = 0;
1243: PetscFunctionReturn(PETSC_SUCCESS);
1244: }
1246: static PetscErrorCode MatTranspose_IS(Mat A, MatReuse reuse, Mat *B)
1247: {
1248: Mat C, lC, lA;
1250: PetscFunctionBegin;
1251: if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *B));
1252: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
1253: ISLocalToGlobalMapping rl2g, cl2g;
1254: PetscBool allow_repeated;
1256: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
1257: PetscCall(MatSetSizes(C, A->cmap->n, A->rmap->n, A->cmap->N, A->rmap->N));
1258: PetscCall(MatSetBlockSizes(C, A->cmap->bs, A->rmap->bs));
1259: PetscCall(MatSetType(C, MATIS));
1260: PetscCall(MatISGetAllowRepeated(A, &allow_repeated));
1261: PetscCall(MatISSetAllowRepeated(C, allow_repeated));
1262: PetscCall(MatGetLocalToGlobalMapping(A, &rl2g, &cl2g));
1263: PetscCall(MatSetLocalToGlobalMapping(C, cl2g, rl2g));
1264: } else C = *B;
1266: /* perform local transposition */
1267: PetscCall(MatISGetLocalMat(A, &lA));
1268: PetscCall(MatTranspose(lA, MAT_INITIAL_MATRIX, &lC));
1269: PetscCall(MatSetLocalToGlobalMapping(lC, lA->cmap->mapping, lA->rmap->mapping));
1270: PetscCall(MatISSetLocalMat(C, lC));
1271: PetscCall(MatDestroy(&lC));
1273: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
1274: *B = C;
1275: } else {
1276: PetscCall(MatHeaderMerge(A, &C));
1277: }
1278: PetscCall(MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY));
1279: PetscCall(MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY));
1280: PetscFunctionReturn(PETSC_SUCCESS);
1281: }
1283: static PetscErrorCode MatDiagonalSet_IS(Mat A, Vec D, InsertMode insmode)
1284: {
1285: Mat_IS *is = (Mat_IS *)A->data;
1287: PetscFunctionBegin;
1288: PetscCheck(!is->allow_repeated || insmode == ADD_VALUES, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "INSERT_VALUES with repeated entries not supported");
1289: if (D) { /* MatShift_IS pass D = NULL */
1290: PetscCall(VecScatterBegin(is->rctx, D, is->y, INSERT_VALUES, SCATTER_FORWARD));
1291: PetscCall(VecScatterEnd(is->rctx, D, is->y, INSERT_VALUES, SCATTER_FORWARD));
1292: }
1293: PetscCall(VecPointwiseDivide(is->y, is->y, is->counter));
1294: PetscCall(MatDiagonalSet(is->A, is->y, insmode));
1295: PetscFunctionReturn(PETSC_SUCCESS);
1296: }
1298: static PetscErrorCode MatShift_IS(Mat A, PetscScalar a)
1299: {
1300: Mat_IS *is = (Mat_IS *)A->data;
1302: PetscFunctionBegin;
1303: PetscCall(VecSet(is->y, a));
1304: PetscCall(MatDiagonalSet_IS(A, NULL, ADD_VALUES));
1305: PetscFunctionReturn(PETSC_SUCCESS);
1306: }
1308: static PetscErrorCode MatSetValuesLocal_SubMat_IS(Mat A, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
1309: {
1310: PetscInt buf[2 * MATIS_MAX_ENTRIES_INSERTION], *rows_l = NULL, *cols_l = NULL;
1312: PetscFunctionBegin;
1313: IndexSpaceGet(buf, m, n, rows_l, cols_l);
1314: PetscCall(ISLocalToGlobalMappingApply(A->rmap->mapping, m, rows, rows_l));
1315: PetscCall(ISLocalToGlobalMappingApply(A->cmap->mapping, n, cols, cols_l));
1316: PetscCall(MatSetValuesLocal_IS(A, m, rows_l, n, cols_l, values, addv));
1317: IndexSpaceRestore(buf, m, n, rows_l, cols_l);
1318: PetscFunctionReturn(PETSC_SUCCESS);
1319: }
1321: static PetscErrorCode MatSetValuesBlockedLocal_SubMat_IS(Mat A, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
1322: {
1323: PetscInt buf[2 * MATIS_MAX_ENTRIES_INSERTION], *rows_l = NULL, *cols_l = NULL, rbs, cbs;
1325: PetscFunctionBegin;
1326: /* We cannot guarantee the local matrix will have the same block size of the original matrix */
1327: PetscCall(ISLocalToGlobalMappingGetBlockSize(A->rmap->mapping, &rbs));
1328: PetscCall(ISLocalToGlobalMappingGetBlockSize(A->cmap->mapping, &cbs));
1329: IndexSpaceGet(buf, m * rbs, n * cbs, rows_l, cols_l);
1330: BlockIndicesExpand(m, rows, rbs, rows_l);
1331: BlockIndicesExpand(n, cols, cbs, cols_l);
1332: PetscCall(ISLocalToGlobalMappingApply(A->rmap->mapping, m * rbs, rows_l, rows_l));
1333: PetscCall(ISLocalToGlobalMappingApply(A->cmap->mapping, n * cbs, cols_l, cols_l));
1334: PetscCall(MatSetValuesLocal_IS(A, m * rbs, rows_l, n * cbs, cols_l, values, addv));
1335: IndexSpaceRestore(buf, m * rbs, n * cbs, rows_l, cols_l);
1336: PetscFunctionReturn(PETSC_SUCCESS);
1337: }
1339: static PetscErrorCode MatZeroRowsLocal_SubMat_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
1340: {
1341: PetscInt *rows_l;
1342: Mat_IS *is = (Mat_IS *)A->data;
1344: PetscFunctionBegin;
1345: PetscCall(PetscMalloc1(n, &rows_l));
1346: PetscCall(ISLocalToGlobalMappingApply(A->rmap->mapping, n, rows, rows_l));
1347: PetscCall(MatZeroRowsLocal(is->islocalref, n, rows_l, diag, x, b));
1348: PetscCall(PetscFree(rows_l));
1349: PetscFunctionReturn(PETSC_SUCCESS);
1350: }
1352: static PetscErrorCode MatZeroRowsColumnsLocal_SubMat_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
1353: {
1354: PetscInt *rows_l;
1355: Mat_IS *is = (Mat_IS *)A->data;
1357: PetscFunctionBegin;
1358: PetscCall(PetscMalloc1(n, &rows_l));
1359: PetscCall(ISLocalToGlobalMappingApply(A->rmap->mapping, n, rows, rows_l));
1360: PetscCall(MatZeroRowsColumnsLocal(is->islocalref, n, rows_l, diag, x, b));
1361: PetscCall(PetscFree(rows_l));
1362: PetscFunctionReturn(PETSC_SUCCESS);
1363: }
1365: static PetscErrorCode MatCreateSubMatrix_IS(Mat mat, IS irow, IS icol, MatReuse scall, Mat *newmat)
1366: {
1367: Mat locmat, newlocmat;
1368: Mat_IS *newmatis;
1369: const PetscInt *idxs;
1370: PetscInt i, m, n;
1372: PetscFunctionBegin;
1373: if (scall == MAT_REUSE_MATRIX) {
1374: PetscBool ismatis;
1376: PetscCall(PetscObjectTypeCompare((PetscObject)*newmat, MATIS, &ismatis));
1377: PetscCheck(ismatis, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_ARG_WRONG, "Cannot reuse matrix! Not of MATIS type");
1378: newmatis = (Mat_IS *)(*newmat)->data;
1379: PetscCheck(newmatis->getsub_ris, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_ARG_WRONG, "Cannot reuse matrix! Misses local row IS");
1380: PetscCheck(newmatis->getsub_cis, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_ARG_WRONG, "Cannot reuse matrix! Misses local col IS");
1381: }
1382: /* irow and icol may not have duplicate entries */
1383: if (PetscDefined(USE_DEBUG)) {
1384: Vec rtest, ltest;
1385: const PetscScalar *array;
1387: PetscCall(MatCreateVecs(mat, <est, &rtest));
1388: PetscCall(ISGetLocalSize(irow, &n));
1389: PetscCall(ISGetIndices(irow, &idxs));
1390: for (i = 0; i < n; i++) PetscCall(VecSetValue(rtest, idxs[i], 1.0, ADD_VALUES));
1391: PetscCall(VecAssemblyBegin(rtest));
1392: PetscCall(VecAssemblyEnd(rtest));
1393: PetscCall(VecGetLocalSize(rtest, &n));
1394: PetscCall(VecGetOwnershipRange(rtest, &m, NULL));
1395: PetscCall(VecGetArrayRead(rtest, &array));
1396: for (i = 0; i < n; i++) PetscCheck(array[i] == 0. || array[i] == 1., PETSC_COMM_SELF, PETSC_ERR_SUP, "Index %" PetscInt_FMT " counted %" PetscInt_FMT " times! Irow may not have duplicate entries", i + m, (PetscInt)PetscRealPart(array[i]));
1397: PetscCall(VecRestoreArrayRead(rtest, &array));
1398: PetscCall(ISRestoreIndices(irow, &idxs));
1399: PetscCall(ISGetLocalSize(icol, &n));
1400: PetscCall(ISGetIndices(icol, &idxs));
1401: for (i = 0; i < n; i++) PetscCall(VecSetValue(ltest, idxs[i], 1.0, ADD_VALUES));
1402: PetscCall(VecAssemblyBegin(ltest));
1403: PetscCall(VecAssemblyEnd(ltest));
1404: PetscCall(VecGetLocalSize(ltest, &n));
1405: PetscCall(VecGetOwnershipRange(ltest, &m, NULL));
1406: PetscCall(VecGetArrayRead(ltest, &array));
1407: for (i = 0; i < n; i++) PetscCheck(array[i] == 0. || array[i] == 1., PETSC_COMM_SELF, PETSC_ERR_SUP, "Index %" PetscInt_FMT " counted %" PetscInt_FMT " times! Icol may not have duplicate entries", i + m, (PetscInt)PetscRealPart(array[i]));
1408: PetscCall(VecRestoreArrayRead(ltest, &array));
1409: PetscCall(ISRestoreIndices(icol, &idxs));
1410: PetscCall(VecDestroy(&rtest));
1411: PetscCall(VecDestroy(<est));
1412: }
1413: if (scall == MAT_INITIAL_MATRIX) {
1414: Mat_IS *matis = (Mat_IS *)mat->data;
1415: ISLocalToGlobalMapping rl2g;
1416: IS is;
1417: PetscInt *lidxs, *lgidxs, *newgidxs;
1418: PetscInt ll, newloc, irbs, icbs, arbs, acbs, rbs, cbs;
1419: PetscBool cong;
1420: MPI_Comm comm;
1422: PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
1423: PetscCall(MatGetBlockSizes(mat, &arbs, &acbs));
1424: PetscCall(ISGetBlockSize(irow, &irbs));
1425: PetscCall(ISGetBlockSize(icol, &icbs));
1426: rbs = arbs == irbs ? irbs : 1;
1427: cbs = acbs == icbs ? icbs : 1;
1428: PetscCall(ISGetLocalSize(irow, &m));
1429: PetscCall(ISGetLocalSize(icol, &n));
1430: PetscCall(MatCreate(comm, newmat));
1431: PetscCall(MatSetType(*newmat, MATIS));
1432: PetscCall(MatISSetAllowRepeated(*newmat, matis->allow_repeated));
1433: PetscCall(MatSetSizes(*newmat, m, n, PETSC_DECIDE, PETSC_DECIDE));
1434: PetscCall(MatSetBlockSizes(*newmat, rbs, cbs));
1435: /* communicate irow to their owners in the layout */
1436: PetscCall(ISGetIndices(irow, &idxs));
1437: PetscCall(PetscLayoutMapLocal(mat->rmap, m, idxs, &ll, &lidxs, &lgidxs));
1438: PetscCall(ISRestoreIndices(irow, &idxs));
1439: PetscCall(PetscArrayzero(matis->sf_rootdata, matis->sf->nroots));
1440: for (i = 0; i < ll; i++) matis->sf_rootdata[lidxs[i]] = lgidxs[i] + 1;
1441: PetscCall(PetscFree(lidxs));
1442: PetscCall(PetscFree(lgidxs));
1443: PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
1444: PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
1445: for (i = 0, newloc = 0; i < matis->sf->nleaves; i++)
1446: if (matis->sf_leafdata[i]) newloc++;
1447: PetscCall(PetscMalloc1(newloc, &newgidxs));
1448: PetscCall(PetscMalloc1(newloc, &lidxs));
1449: for (i = 0, newloc = 0; i < matis->sf->nleaves; i++)
1450: if (matis->sf_leafdata[i]) {
1451: lidxs[newloc] = i;
1452: newgidxs[newloc++] = matis->sf_leafdata[i] - 1;
1453: }
1454: PetscCall(ISCreateGeneral(comm, newloc, newgidxs, PETSC_OWN_POINTER, &is));
1455: PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g));
1456: PetscCall(ISLocalToGlobalMappingSetBlockSize(rl2g, rbs));
1457: PetscCall(ISDestroy(&is));
1458: /* local is to extract local submatrix */
1459: newmatis = (Mat_IS *)(*newmat)->data;
1460: PetscCall(ISCreateGeneral(comm, newloc, lidxs, PETSC_OWN_POINTER, &newmatis->getsub_ris));
1461: PetscCall(MatHasCongruentLayouts(mat, &cong));
1462: if (cong && irow == icol && matis->csf == matis->sf) {
1463: PetscCall(MatSetLocalToGlobalMapping(*newmat, rl2g, rl2g));
1464: PetscCall(PetscObjectReference((PetscObject)newmatis->getsub_ris));
1465: newmatis->getsub_cis = newmatis->getsub_ris;
1466: } else {
1467: ISLocalToGlobalMapping cl2g;
1469: /* communicate icol to their owners in the layout */
1470: PetscCall(ISGetIndices(icol, &idxs));
1471: PetscCall(PetscLayoutMapLocal(mat->cmap, n, idxs, &ll, &lidxs, &lgidxs));
1472: PetscCall(ISRestoreIndices(icol, &idxs));
1473: PetscCall(PetscArrayzero(matis->csf_rootdata, matis->csf->nroots));
1474: for (i = 0; i < ll; i++) matis->csf_rootdata[lidxs[i]] = lgidxs[i] + 1;
1475: PetscCall(PetscFree(lidxs));
1476: PetscCall(PetscFree(lgidxs));
1477: PetscCall(PetscSFBcastBegin(matis->csf, MPIU_INT, matis->csf_rootdata, matis->csf_leafdata, MPI_REPLACE));
1478: PetscCall(PetscSFBcastEnd(matis->csf, MPIU_INT, matis->csf_rootdata, matis->csf_leafdata, MPI_REPLACE));
1479: for (i = 0, newloc = 0; i < matis->csf->nleaves; i++)
1480: if (matis->csf_leafdata[i]) newloc++;
1481: PetscCall(PetscMalloc1(newloc, &newgidxs));
1482: PetscCall(PetscMalloc1(newloc, &lidxs));
1483: for (i = 0, newloc = 0; i < matis->csf->nleaves; i++)
1484: if (matis->csf_leafdata[i]) {
1485: lidxs[newloc] = i;
1486: newgidxs[newloc++] = matis->csf_leafdata[i] - 1;
1487: }
1488: PetscCall(ISCreateGeneral(comm, newloc, newgidxs, PETSC_OWN_POINTER, &is));
1489: PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g));
1490: PetscCall(ISLocalToGlobalMappingSetBlockSize(cl2g, cbs));
1491: PetscCall(ISDestroy(&is));
1492: /* local is to extract local submatrix */
1493: PetscCall(ISCreateGeneral(comm, newloc, lidxs, PETSC_OWN_POINTER, &newmatis->getsub_cis));
1494: PetscCall(MatSetLocalToGlobalMapping(*newmat, rl2g, cl2g));
1495: PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
1496: }
1497: PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
1498: } else {
1499: PetscCall(MatISGetLocalMat(*newmat, &newlocmat));
1500: }
1501: PetscCall(MatISGetLocalMat(mat, &locmat));
1502: newmatis = (Mat_IS *)(*newmat)->data;
1503: PetscCall(MatCreateSubMatrix(locmat, newmatis->getsub_ris, newmatis->getsub_cis, scall, &newlocmat));
1504: if (scall == MAT_INITIAL_MATRIX) {
1505: PetscCall(MatISSetLocalMat(*newmat, newlocmat));
1506: PetscCall(MatDestroy(&newlocmat));
1507: }
1508: PetscCall(MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY));
1509: PetscCall(MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY));
1510: PetscFunctionReturn(PETSC_SUCCESS);
1511: }
1513: static PetscErrorCode MatCopy_IS(Mat A, Mat B, MatStructure str)
1514: {
1515: Mat_IS *a = (Mat_IS *)A->data, *b;
1516: PetscBool ismatis;
1518: PetscFunctionBegin;
1519: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATIS, &ismatis));
1520: PetscCheck(ismatis, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Need to be implemented");
1521: b = (Mat_IS *)B->data;
1522: PetscCall(MatCopy(a->A, b->A, str));
1523: PetscCall(PetscObjectStateIncrease((PetscObject)B));
1524: PetscFunctionReturn(PETSC_SUCCESS);
1525: }
1527: static PetscErrorCode MatISSetUpSF_IS(Mat B)
1528: {
1529: Mat_IS *matis = (Mat_IS *)B->data;
1530: const PetscInt *gidxs;
1531: PetscInt nleaves;
1533: PetscFunctionBegin;
1534: if (matis->sf) PetscFunctionReturn(PETSC_SUCCESS);
1535: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)B), &matis->sf));
1536: PetscCall(ISLocalToGlobalMappingGetIndices(matis->rmapping, &gidxs));
1537: PetscCall(ISLocalToGlobalMappingGetSize(matis->rmapping, &nleaves));
1538: PetscCall(PetscSFSetGraphLayout(matis->sf, B->rmap, nleaves, NULL, PETSC_OWN_POINTER, gidxs));
1539: PetscCall(ISLocalToGlobalMappingRestoreIndices(matis->rmapping, &gidxs));
1540: PetscCall(PetscMalloc2(matis->sf->nroots, &matis->sf_rootdata, matis->sf->nleaves, &matis->sf_leafdata));
1541: if (matis->rmapping != matis->cmapping) { /* setup SF for columns */
1542: PetscCall(ISLocalToGlobalMappingGetSize(matis->cmapping, &nleaves));
1543: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)B), &matis->csf));
1544: PetscCall(ISLocalToGlobalMappingGetIndices(matis->cmapping, &gidxs));
1545: PetscCall(PetscSFSetGraphLayout(matis->csf, B->cmap, nleaves, NULL, PETSC_OWN_POINTER, gidxs));
1546: PetscCall(ISLocalToGlobalMappingRestoreIndices(matis->cmapping, &gidxs));
1547: PetscCall(PetscMalloc2(matis->csf->nroots, &matis->csf_rootdata, matis->csf->nleaves, &matis->csf_leafdata));
1548: } else {
1549: matis->csf = matis->sf;
1550: matis->csf_leafdata = matis->sf_leafdata;
1551: matis->csf_rootdata = matis->sf_rootdata;
1552: }
1553: PetscFunctionReturn(PETSC_SUCCESS);
1554: }
1556: /*@
1557: MatISGetAllowRepeated - Get the flag to allow repeated entries in the local to global map
1559: Not Collective
1561: Input Parameter:
1562: . A - the matrix
1564: Output Parameter:
1565: . flg - the boolean flag
1567: Level: intermediate
1569: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateIS()`, `MatSetLocalToGlobalMapping()`, `MatISSetAllowRepeated()`
1570: @*/
1571: PetscErrorCode MatISGetAllowRepeated(Mat A, PetscBool *flg)
1572: {
1573: PetscBool ismatis;
1575: PetscFunctionBegin;
1577: PetscAssertPointer(flg, 2);
1578: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATIS, &ismatis));
1579: PetscCheck(ismatis, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not for matrix type %s", ((PetscObject)A)->type_name);
1580: *flg = ((Mat_IS *)A->data)->allow_repeated;
1581: PetscFunctionReturn(PETSC_SUCCESS);
1582: }
1584: /*@
1585: MatISSetAllowRepeated - Set the flag to allow repeated entries in the local to global map
1587: Logically Collective
1589: Input Parameters:
1590: + A - the matrix
1591: - flg - the boolean flag
1593: Level: intermediate
1595: Notes:
1596: The default value is `PETSC_FALSE`.
1597: When called AFTER calling `MatSetLocalToGlobalMapping()` it will recreate the local matrices
1598: if `flg` is different from the previously set value.
1599: Specifically, when `flg` is true it will just recreate the local matrices, while if
1600: `flg` is false will assemble the local matrices summing up repeated entries.
1602: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateIS()`, `MatSetLocalToGlobalMapping()`, `MatISGetAllowRepeated()`
1603: @*/
1604: PetscErrorCode MatISSetAllowRepeated(Mat A, PetscBool flg)
1605: {
1606: PetscFunctionBegin;
1610: PetscTryMethod(A, "MatISSetAllowRepeated_C", (Mat, PetscBool), (A, flg));
1611: PetscFunctionReturn(PETSC_SUCCESS);
1612: }
1614: static PetscErrorCode MatISSetAllowRepeated_IS(Mat A, PetscBool flg)
1615: {
1616: Mat_IS *matis = (Mat_IS *)A->data;
1617: Mat lA = NULL;
1618: ISLocalToGlobalMapping lrmap, lcmap;
1620: PetscFunctionBegin;
1621: if (flg == matis->allow_repeated) PetscFunctionReturn(PETSC_SUCCESS);
1622: if (!matis->A) { /* matrix has not been preallocated yet */
1623: matis->allow_repeated = flg;
1624: PetscFunctionReturn(PETSC_SUCCESS);
1625: }
1626: PetscCheck(!matis->islocalref, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not implemented for local references");
1627: if (matis->allow_repeated) { /* we will assemble the old local matrix if needed */
1628: lA = matis->A;
1629: PetscCall(PetscObjectReference((PetscObject)lA));
1630: }
1631: /* In case flg is True, we only recreate the local matrix */
1632: matis->allow_repeated = flg;
1633: PetscCall(MatSetLocalToGlobalMapping(A, A->rmap->mapping, A->cmap->mapping));
1634: if (lA) { /* assemble previous local matrix if needed */
1635: Mat nA = matis->A;
1637: PetscCall(MatGetLocalToGlobalMapping(nA, &lrmap, &lcmap));
1638: if (!lrmap && !lcmap) {
1639: PetscCall(MatISSetLocalMat(A, lA));
1640: } else {
1641: Mat P = NULL, R = NULL;
1642: MatProductType ptype;
1644: if (lrmap == lcmap) {
1645: ptype = MATPRODUCT_PtAP;
1646: PetscCall(MatCreateFromISLocalToGlobalMapping(lcmap, nA, PETSC_TRUE, PETSC_FALSE, NULL, &P));
1647: PetscCall(MatProductCreate(lA, P, NULL, &nA));
1648: } else {
1649: if (lcmap) PetscCall(MatCreateFromISLocalToGlobalMapping(lcmap, nA, PETSC_TRUE, PETSC_FALSE, NULL, &P));
1650: if (lrmap) PetscCall(MatCreateFromISLocalToGlobalMapping(lrmap, nA, PETSC_FALSE, PETSC_TRUE, NULL, &R));
1651: if (R && P) {
1652: ptype = MATPRODUCT_ABC;
1653: PetscCall(MatProductCreate(R, lA, P, &nA));
1654: } else if (R) {
1655: ptype = MATPRODUCT_AB;
1656: PetscCall(MatProductCreate(R, lA, NULL, &nA));
1657: } else {
1658: ptype = MATPRODUCT_AB;
1659: PetscCall(MatProductCreate(lA, P, NULL, &nA));
1660: }
1661: }
1662: PetscCall(MatProductSetType(nA, ptype));
1663: PetscCall(MatProductSetFromOptions(nA));
1664: PetscCall(MatProductSymbolic(nA));
1665: PetscCall(MatProductNumeric(nA));
1666: PetscCall(MatProductClear(nA));
1667: PetscCall(MatConvert(nA, matis->lmattype, MAT_INPLACE_MATRIX, &nA));
1668: PetscCall(MatISSetLocalMat(A, nA));
1669: PetscCall(MatDestroy(&nA));
1670: PetscCall(MatDestroy(&P));
1671: PetscCall(MatDestroy(&R));
1672: }
1673: }
1674: PetscCall(MatDestroy(&lA));
1675: PetscFunctionReturn(PETSC_SUCCESS);
1676: }
1678: /*@
1679: MatISStoreL2L - Store local-to-local operators during the Galerkin process of computing `MatPtAP()`
1681: Logically Collective
1683: Input Parameters:
1684: + A - the matrix
1685: - store - the boolean flag
1687: Level: advanced
1689: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateIS()`, `MatISSetPreallocation()`, `MatPtAP()`
1690: @*/
1691: PetscErrorCode MatISStoreL2L(Mat A, PetscBool store)
1692: {
1693: PetscFunctionBegin;
1697: PetscTryMethod(A, "MatISStoreL2L_C", (Mat, PetscBool), (A, store));
1698: PetscFunctionReturn(PETSC_SUCCESS);
1699: }
1701: static PetscErrorCode MatISStoreL2L_IS(Mat A, PetscBool store)
1702: {
1703: Mat_IS *matis = (Mat_IS *)A->data;
1705: PetscFunctionBegin;
1706: matis->storel2l = store;
1707: if (!store) PetscCall(PetscObjectCompose((PetscObject)A, "_MatIS_PtAP_l2l", NULL));
1708: PetscFunctionReturn(PETSC_SUCCESS);
1709: }
1711: /*@
1712: MatISFixLocalEmpty - Compress out zero local rows from the local matrices
1714: Logically Collective
1716: Input Parameters:
1717: + A - the matrix
1718: - fix - the boolean flag
1720: Level: advanced
1722: Note:
1723: When `fix` is `PETSC_TRUE`, new local matrices and l2g maps are generated during the final assembly process.
1725: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatCreate()`, `MatCreateIS()`, `MatISSetPreallocation()`, `MatAssemblyEnd()`, `MAT_FINAL_ASSEMBLY`
1726: @*/
1727: PetscErrorCode MatISFixLocalEmpty(Mat A, PetscBool fix)
1728: {
1729: PetscFunctionBegin;
1733: PetscTryMethod(A, "MatISFixLocalEmpty_C", (Mat, PetscBool), (A, fix));
1734: PetscFunctionReturn(PETSC_SUCCESS);
1735: }
1737: static PetscErrorCode MatISFixLocalEmpty_IS(Mat A, PetscBool fix)
1738: {
1739: Mat_IS *matis = (Mat_IS *)A->data;
1741: PetscFunctionBegin;
1742: matis->locempty = fix;
1743: PetscFunctionReturn(PETSC_SUCCESS);
1744: }
1746: /*@
1747: MatISSetPreallocation - Preallocates memory for a `MATIS` parallel matrix.
1749: Collective
1751: Input Parameters:
1752: + B - the matrix
1753: . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix
1754: (same value is used for all local rows)
1755: . d_nnz - array containing the number of nonzeros in the various rows of the
1756: DIAGONAL portion of the local submatrix (possibly different for each row)
1757: or `NULL`, if `d_nz` is used to specify the nonzero structure.
1758: The size of this array is equal to the number of local rows, i.e `m`.
1759: For matrices that will be factored, you must leave room for (and set)
1760: the diagonal entry even if it is zero.
1761: . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local
1762: submatrix (same value is used for all local rows).
1763: - o_nnz - array containing the number of nonzeros in the various rows of the
1764: OFF-DIAGONAL portion of the local submatrix (possibly different for
1765: each row) or `NULL`, if `o_nz` is used to specify the nonzero
1766: structure. The size of this array is equal to the number
1767: of local rows, i.e `m`.
1769: If the *_nnz parameter is given then the *_nz parameter is ignored
1771: Level: intermediate
1773: Note:
1774: This function has the same interface as the `MATMPIAIJ` preallocation routine in order to simplify the transition
1775: from the asssembled format to the unassembled one. It overestimates the preallocation of `MATIS` local
1776: matrices; for exact preallocation, the user should set the preallocation directly on local matrix objects.
1778: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateIS()`, `MatMPIAIJSetPreallocation()`, `MatISGetLocalMat()`, `MATIS`
1779: @*/
1780: PetscErrorCode MatISSetPreallocation(Mat B, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[])
1781: {
1782: PetscFunctionBegin;
1785: PetscTryMethod(B, "MatISSetPreallocation_C", (Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[]), (B, d_nz, d_nnz, o_nz, o_nnz));
1786: PetscFunctionReturn(PETSC_SUCCESS);
1787: }
1789: static PetscErrorCode MatISSetPreallocation_IS(Mat B, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[])
1790: {
1791: Mat_IS *matis = (Mat_IS *)B->data;
1792: PetscInt bs, i, nlocalcols;
1794: PetscFunctionBegin;
1795: PetscCall(MatSetUp(B));
1796: if (!d_nnz)
1797: for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] = d_nz;
1798: else
1799: for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] = d_nnz[i];
1801: if (!o_nnz)
1802: for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] += o_nz;
1803: else
1804: for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] += o_nnz[i];
1806: PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
1807: PetscCall(MatGetSize(matis->A, NULL, &nlocalcols));
1808: PetscCall(MatGetBlockSize(matis->A, &bs));
1809: PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
1811: for (i = 0; i < matis->sf->nleaves; i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i], nlocalcols);
1812: PetscCall(MatSeqAIJSetPreallocation(matis->A, 0, matis->sf_leafdata));
1813: #if defined(PETSC_HAVE_HYPRE)
1814: PetscCall(MatHYPRESetPreallocation(matis->A, 0, matis->sf_leafdata, 0, NULL));
1815: #endif
1817: for (i = 0; i < matis->sf->nleaves / bs; i++) {
1818: matis->sf_leafdata[i] = matis->sf_leafdata[i * bs] / bs;
1819: for (PetscInt b = 1; b < bs; b++) matis->sf_leafdata[i] = PetscMax(matis->sf_leafdata[i], matis->sf_leafdata[i * bs + b] / bs);
1820: }
1821: PetscCall(MatSeqBAIJSetPreallocation(matis->A, bs, 0, matis->sf_leafdata));
1823: nlocalcols /= bs;
1824: for (i = 0; i < matis->sf->nleaves / bs; i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i], nlocalcols - i);
1825: PetscCall(MatSeqSBAIJSetPreallocation(matis->A, bs, 0, matis->sf_leafdata));
1827: /* for other matrix types */
1828: PetscCall(MatSetUp(matis->A));
1829: PetscFunctionReturn(PETSC_SUCCESS);
1830: }
1832: PETSC_INTERN PetscErrorCode MatConvert_IS_XAIJ(Mat mat, MatType mtype, MatReuse reuse, Mat *M)
1833: {
1834: Mat_IS *matis = (Mat_IS *)mat->data;
1835: Mat local_mat = NULL, MT;
1836: PetscInt rbs, cbs, rows, cols, lrows, lcols;
1837: PetscInt local_rows, local_cols;
1838: PetscBool isseqdense, isseqsbaij, isseqaij, isseqbaij;
1839: PetscMPIInt size;
1840: const PetscScalar *array;
1842: PetscFunctionBegin;
1843: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)mat), &size));
1844: if (size == 1 && mat->rmap->N == matis->A->rmap->N && mat->cmap->N == matis->A->cmap->N && !matis->allow_repeated) {
1845: Mat B;
1846: IS irows = NULL, icols = NULL;
1847: PetscInt rbs, cbs;
1849: PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->rmapping, &rbs));
1850: PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->cmapping, &cbs));
1851: if (reuse != MAT_REUSE_MATRIX) { /* check if l2g maps are one-to-one */
1852: IS rows, cols;
1853: const PetscInt *ridxs, *cidxs;
1854: PetscInt i, nw;
1855: PetscBT work;
1857: PetscCall(ISLocalToGlobalMappingGetBlockIndices(matis->rmapping, &ridxs));
1858: PetscCall(ISLocalToGlobalMappingGetSize(matis->rmapping, &nw));
1859: nw = nw / rbs;
1860: PetscCall(PetscBTCreate(nw, &work));
1861: for (i = 0; i < nw; i++) PetscCall(PetscBTSet(work, ridxs[i]));
1862: for (i = 0; i < nw; i++)
1863: if (!PetscBTLookup(work, i)) break;
1864: if (i == nw) {
1865: PetscCall(ISCreateBlock(PETSC_COMM_SELF, rbs, nw, ridxs, PETSC_USE_POINTER, &rows));
1866: PetscCall(ISSetPermutation(rows));
1867: PetscCall(ISInvertPermutation(rows, PETSC_DECIDE, &irows));
1868: PetscCall(ISDestroy(&rows));
1869: }
1870: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(matis->rmapping, &ridxs));
1871: PetscCall(PetscBTDestroy(&work));
1872: if (irows && matis->rmapping != matis->cmapping) {
1873: PetscCall(ISLocalToGlobalMappingGetBlockIndices(matis->cmapping, &cidxs));
1874: PetscCall(ISLocalToGlobalMappingGetSize(matis->cmapping, &nw));
1875: nw = nw / cbs;
1876: PetscCall(PetscBTCreate(nw, &work));
1877: for (i = 0; i < nw; i++) PetscCall(PetscBTSet(work, cidxs[i]));
1878: for (i = 0; i < nw; i++)
1879: if (!PetscBTLookup(work, i)) break;
1880: if (i == nw) {
1881: PetscCall(ISCreateBlock(PETSC_COMM_SELF, cbs, nw, cidxs, PETSC_USE_POINTER, &cols));
1882: PetscCall(ISSetPermutation(cols));
1883: PetscCall(ISInvertPermutation(cols, PETSC_DECIDE, &icols));
1884: PetscCall(ISDestroy(&cols));
1885: }
1886: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(matis->cmapping, &cidxs));
1887: PetscCall(PetscBTDestroy(&work));
1888: } else if (irows) {
1889: PetscCall(PetscObjectReference((PetscObject)irows));
1890: icols = irows;
1891: }
1892: } else {
1893: PetscCall(PetscObjectQuery((PetscObject)*M, "_MatIS_IS_XAIJ_irows", (PetscObject *)&irows));
1894: PetscCall(PetscObjectQuery((PetscObject)*M, "_MatIS_IS_XAIJ_icols", (PetscObject *)&icols));
1895: PetscCall(PetscObjectReference((PetscObject)irows));
1896: PetscCall(PetscObjectReference((PetscObject)icols));
1897: }
1898: if (!irows || !icols) {
1899: PetscCall(ISDestroy(&icols));
1900: PetscCall(ISDestroy(&irows));
1901: goto general_assembly;
1902: }
1903: PetscCall(MatConvert(matis->A, mtype, MAT_INITIAL_MATRIX, &B));
1904: if (reuse != MAT_INPLACE_MATRIX) {
1905: PetscCall(MatCreateSubMatrix(B, irows, icols, reuse, M));
1906: PetscCall(PetscObjectCompose((PetscObject)*M, "_MatIS_IS_XAIJ_irows", (PetscObject)irows));
1907: PetscCall(PetscObjectCompose((PetscObject)*M, "_MatIS_IS_XAIJ_icols", (PetscObject)icols));
1908: } else {
1909: Mat C;
1911: PetscCall(MatCreateSubMatrix(B, irows, icols, MAT_INITIAL_MATRIX, &C));
1912: PetscCall(MatHeaderReplace(mat, &C));
1913: }
1914: PetscCall(MatDestroy(&B));
1915: PetscCall(ISDestroy(&icols));
1916: PetscCall(ISDestroy(&irows));
1917: PetscFunctionReturn(PETSC_SUCCESS);
1918: }
1919: general_assembly:
1920: PetscCall(MatGetSize(mat, &rows, &cols));
1921: PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->rmapping, &rbs));
1922: PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->cmapping, &cbs));
1923: PetscCall(MatGetLocalSize(mat, &lrows, &lcols));
1924: PetscCall(MatGetSize(matis->A, &local_rows, &local_cols));
1925: PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQDENSE, &isseqdense));
1926: PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQAIJ, &isseqaij));
1927: PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQBAIJ, &isseqbaij));
1928: PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQSBAIJ, &isseqsbaij));
1929: PetscCheck(isseqdense || isseqaij || isseqbaij || isseqsbaij, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for matrix type %s", ((PetscObject)matis->A)->type_name);
1930: if (PetscDefined(USE_DEBUG)) {
1931: PetscBool lb[4], bb[4];
1933: lb[0] = isseqdense;
1934: lb[1] = isseqaij;
1935: lb[2] = isseqbaij;
1936: lb[3] = isseqsbaij;
1937: PetscCallMPI(MPIU_Allreduce(lb, bb, 4, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)mat)));
1938: PetscCheck(bb[0] || bb[1] || bb[2] || bb[3], PETSC_COMM_SELF, PETSC_ERR_SUP, "Local matrices must have the same type");
1939: }
1941: if (reuse != MAT_REUSE_MATRIX) {
1942: PetscCount ncoo;
1943: PetscInt *coo_i, *coo_j;
1945: PetscCall(MatCreate(PetscObjectComm((PetscObject)mat), &MT));
1946: PetscCall(MatSetSizes(MT, lrows, lcols, rows, cols));
1947: PetscCall(MatSetType(MT, mtype));
1948: PetscCall(MatSetBlockSizes(MT, rbs, cbs));
1949: if (!isseqaij && !isseqdense) {
1950: PetscCall(MatConvert(matis->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &local_mat));
1951: } else {
1952: PetscCall(PetscObjectReference((PetscObject)matis->A));
1953: local_mat = matis->A;
1954: }
1955: PetscCall(MatSetLocalToGlobalMapping(MT, matis->rmapping, matis->cmapping));
1956: if (isseqdense) {
1957: PetscInt nr, nc;
1959: PetscCall(MatGetSize(local_mat, &nr, &nc));
1960: ncoo = nr * nc;
1961: PetscCall(PetscMalloc2(ncoo, &coo_i, ncoo, &coo_j));
1962: for (PetscInt j = 0; j < nc; j++) {
1963: for (PetscInt i = 0; i < nr; i++) {
1964: coo_i[j * nr + i] = i;
1965: coo_j[j * nr + i] = j;
1966: }
1967: }
1968: } else {
1969: const PetscInt *ii, *jj;
1970: PetscInt nr;
1971: PetscBool done;
1973: PetscCall(MatGetRowIJ(local_mat, 0, PETSC_FALSE, PETSC_FALSE, &nr, &ii, &jj, &done));
1974: PetscCheck(done, PetscObjectComm((PetscObject)local_mat), PETSC_ERR_PLIB, "Error in MatGetRowIJ");
1975: ncoo = ii[nr];
1976: PetscCall(PetscMalloc2(ncoo, &coo_i, ncoo, &coo_j));
1977: PetscCall(PetscArraycpy(coo_j, jj, ncoo));
1978: for (PetscInt i = 0; i < nr; i++) {
1979: for (PetscInt j = ii[i]; j < ii[i + 1]; j++) coo_i[j] = i;
1980: }
1981: PetscCall(MatRestoreRowIJ(local_mat, 0, PETSC_FALSE, PETSC_FALSE, &nr, &ii, &jj, &done));
1982: PetscCheck(done, PetscObjectComm((PetscObject)local_mat), PETSC_ERR_PLIB, "Error in MatRestoreRowIJ");
1983: }
1984: PetscCall(MatSetPreallocationCOOLocal(MT, ncoo, coo_i, coo_j));
1985: PetscCall(PetscFree2(coo_i, coo_j));
1986: } else {
1987: PetscInt mrbs, mcbs, mrows, mcols, mlrows, mlcols;
1989: /* some checks */
1990: MT = *M;
1991: PetscCall(MatGetBlockSizes(MT, &mrbs, &mcbs));
1992: PetscCall(MatGetSize(MT, &mrows, &mcols));
1993: PetscCall(MatGetLocalSize(MT, &mlrows, &mlcols));
1994: PetscCheck(mrows == rows, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of rows (%" PetscInt_FMT " != %" PetscInt_FMT ")", rows, mrows);
1995: PetscCheck(mcols == cols, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of cols (%" PetscInt_FMT " != %" PetscInt_FMT ")", cols, mcols);
1996: PetscCheck(mlrows == lrows, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of local rows (%" PetscInt_FMT " != %" PetscInt_FMT ")", lrows, mlrows);
1997: PetscCheck(mlcols == lcols, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of local cols (%" PetscInt_FMT " != %" PetscInt_FMT ")", lcols, mlcols);
1998: PetscCheck(mrbs == rbs, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong row block size (%" PetscInt_FMT " != %" PetscInt_FMT ")", rbs, mrbs);
1999: PetscCheck(mcbs == cbs, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong col block size (%" PetscInt_FMT " != %" PetscInt_FMT ")", cbs, mcbs);
2000: PetscCall(MatZeroEntries(MT));
2001: if (!isseqaij && !isseqdense) {
2002: PetscCall(MatConvert(matis->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &local_mat));
2003: } else {
2004: PetscCall(PetscObjectReference((PetscObject)matis->A));
2005: local_mat = matis->A;
2006: }
2007: }
2009: /* Set values */
2010: if (isseqdense) {
2011: PetscCall(MatDenseGetArrayRead(local_mat, &array));
2012: PetscCall(MatSetValuesCOO(MT, array, INSERT_VALUES));
2013: PetscCall(MatDenseRestoreArrayRead(local_mat, &array));
2014: } else {
2015: PetscCall(MatSeqAIJGetArrayRead(local_mat, &array));
2016: PetscCall(MatSetValuesCOO(MT, array, INSERT_VALUES));
2017: PetscCall(MatSeqAIJRestoreArrayRead(local_mat, &array));
2018: }
2019: PetscCall(MatDestroy(&local_mat));
2020: PetscCall(MatAssemblyBegin(MT, MAT_FINAL_ASSEMBLY));
2021: PetscCall(MatAssemblyEnd(MT, MAT_FINAL_ASSEMBLY));
2022: if (reuse == MAT_INPLACE_MATRIX) {
2023: PetscCall(MatHeaderReplace(mat, &MT));
2024: } else if (reuse == MAT_INITIAL_MATRIX) {
2025: *M = MT;
2026: }
2027: PetscFunctionReturn(PETSC_SUCCESS);
2028: }
2030: static PetscErrorCode MatDuplicate_IS(Mat mat, MatDuplicateOption op, Mat *newmat)
2031: {
2032: Mat_IS *matis = (Mat_IS *)mat->data;
2033: PetscInt rbs, cbs, m, n, M, N;
2034: Mat B, localmat;
2036: PetscFunctionBegin;
2037: PetscCall(ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping, &rbs));
2038: PetscCall(ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping, &cbs));
2039: PetscCall(MatGetSize(mat, &M, &N));
2040: PetscCall(MatGetLocalSize(mat, &m, &n));
2041: PetscCall(MatCreate(PetscObjectComm((PetscObject)mat), &B));
2042: PetscCall(MatSetSizes(B, m, n, M, N));
2043: PetscCall(MatSetBlockSize(B, rbs == cbs ? rbs : 1));
2044: PetscCall(MatSetType(B, MATIS));
2045: PetscCall(MatISSetLocalMatType(B, matis->lmattype));
2046: PetscCall(MatISSetAllowRepeated(B, matis->allow_repeated));
2047: PetscCall(MatSetLocalToGlobalMapping(B, mat->rmap->mapping, mat->cmap->mapping));
2048: PetscCall(MatDuplicate(matis->A, op, &localmat));
2049: PetscCall(MatSetLocalToGlobalMapping(localmat, matis->A->rmap->mapping, matis->A->cmap->mapping));
2050: PetscCall(MatISSetLocalMat(B, localmat));
2051: PetscCall(MatDestroy(&localmat));
2052: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
2053: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
2054: *newmat = B;
2055: PetscFunctionReturn(PETSC_SUCCESS);
2056: }
2058: static PetscErrorCode MatIsHermitian_IS(Mat A, PetscReal tol, PetscBool *flg)
2059: {
2060: Mat_IS *matis = (Mat_IS *)A->data;
2061: PetscBool local_sym;
2063: PetscFunctionBegin;
2064: PetscCall(MatIsHermitian(matis->A, tol, &local_sym));
2065: PetscCallMPI(MPIU_Allreduce(&local_sym, flg, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
2066: PetscFunctionReturn(PETSC_SUCCESS);
2067: }
2069: static PetscErrorCode MatIsSymmetric_IS(Mat A, PetscReal tol, PetscBool *flg)
2070: {
2071: Mat_IS *matis = (Mat_IS *)A->data;
2072: PetscBool local_sym;
2074: PetscFunctionBegin;
2075: if (matis->rmapping != matis->cmapping) {
2076: *flg = PETSC_FALSE;
2077: PetscFunctionReturn(PETSC_SUCCESS);
2078: }
2079: PetscCall(MatIsSymmetric(matis->A, tol, &local_sym));
2080: PetscCallMPI(MPIU_Allreduce(&local_sym, flg, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
2081: PetscFunctionReturn(PETSC_SUCCESS);
2082: }
2084: static PetscErrorCode MatIsStructurallySymmetric_IS(Mat A, PetscBool *flg)
2085: {
2086: Mat_IS *matis = (Mat_IS *)A->data;
2087: PetscBool local_sym;
2089: PetscFunctionBegin;
2090: if (matis->rmapping != matis->cmapping) {
2091: *flg = PETSC_FALSE;
2092: PetscFunctionReturn(PETSC_SUCCESS);
2093: }
2094: PetscCall(MatIsStructurallySymmetric(matis->A, &local_sym));
2095: PetscCallMPI(MPIU_Allreduce(&local_sym, flg, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
2096: PetscFunctionReturn(PETSC_SUCCESS);
2097: }
2099: static PetscErrorCode MatDestroy_IS(Mat A)
2100: {
2101: Mat_IS *b = (Mat_IS *)A->data;
2103: PetscFunctionBegin;
2104: PetscCall(PetscFree(b->bdiag));
2105: PetscCall(PetscFree(b->lmattype));
2106: PetscCall(MatDestroy(&b->A));
2107: PetscCall(VecScatterDestroy(&b->cctx));
2108: PetscCall(VecScatterDestroy(&b->rctx));
2109: PetscCall(VecDestroy(&b->x));
2110: PetscCall(VecDestroy(&b->y));
2111: PetscCall(VecDestroy(&b->counter));
2112: PetscCall(ISDestroy(&b->getsub_ris));
2113: PetscCall(ISDestroy(&b->getsub_cis));
2114: if (b->sf != b->csf) {
2115: PetscCall(PetscSFDestroy(&b->csf));
2116: PetscCall(PetscFree2(b->csf_rootdata, b->csf_leafdata));
2117: } else b->csf = NULL;
2118: PetscCall(PetscSFDestroy(&b->sf));
2119: PetscCall(PetscFree2(b->sf_rootdata, b->sf_leafdata));
2120: PetscCall(ISLocalToGlobalMappingDestroy(&b->rmapping));
2121: PetscCall(ISLocalToGlobalMappingDestroy(&b->cmapping));
2122: PetscCall(MatDestroy(&b->dA));
2123: PetscCall(MatDestroy(&b->assembledA));
2124: PetscCall(PetscFree(A->data));
2125: PetscCall(PetscObjectChangeTypeName((PetscObject)A, NULL));
2126: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMatType_C", NULL));
2127: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalMat_C", NULL));
2128: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMat_C", NULL));
2129: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISRestoreLocalMat_C", NULL));
2130: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetPreallocation_C", NULL));
2131: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISStoreL2L_C", NULL));
2132: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISFixLocalEmpty_C", NULL));
2133: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalToGlobalMapping_C", NULL));
2134: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpiaij_C", NULL));
2135: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpibaij_C", NULL));
2136: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpisbaij_C", NULL));
2137: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqaij_C", NULL));
2138: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqbaij_C", NULL));
2139: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqsbaij_C", NULL));
2140: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_aij_C", NULL));
2141: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOOLocal_C", NULL));
2142: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOO_C", NULL));
2143: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", NULL));
2144: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetAllowRepeated_C", NULL));
2145: PetscFunctionReturn(PETSC_SUCCESS);
2146: }
2148: static PetscErrorCode MatMult_IS(Mat A, Vec x, Vec y)
2149: {
2150: Mat_IS *is = (Mat_IS *)A->data;
2151: PetscScalar zero = 0.0;
2153: PetscFunctionBegin;
2154: /* scatter the global vector x into the local work vector */
2155: PetscCall(VecScatterBegin(is->cctx, x, is->x, INSERT_VALUES, SCATTER_FORWARD));
2156: PetscCall(VecScatterEnd(is->cctx, x, is->x, INSERT_VALUES, SCATTER_FORWARD));
2158: /* multiply the local matrix */
2159: PetscCall(MatMult(is->A, is->x, is->y));
2161: /* scatter product back into global memory */
2162: PetscCall(VecSet(y, zero));
2163: PetscCall(VecScatterBegin(is->rctx, is->y, y, ADD_VALUES, SCATTER_REVERSE));
2164: PetscCall(VecScatterEnd(is->rctx, is->y, y, ADD_VALUES, SCATTER_REVERSE));
2165: PetscFunctionReturn(PETSC_SUCCESS);
2166: }
2168: static PetscErrorCode MatMultAdd_IS(Mat A, Vec v1, Vec v2, Vec v3)
2169: {
2170: Vec temp_vec;
2172: PetscFunctionBegin; /* v3 = v2 + A * v1.*/
2173: if (v3 != v2) {
2174: PetscCall(MatMult(A, v1, v3));
2175: PetscCall(VecAXPY(v3, 1.0, v2));
2176: } else {
2177: PetscCall(VecDuplicate(v2, &temp_vec));
2178: PetscCall(MatMult(A, v1, temp_vec));
2179: PetscCall(VecAXPY(temp_vec, 1.0, v2));
2180: PetscCall(VecCopy(temp_vec, v3));
2181: PetscCall(VecDestroy(&temp_vec));
2182: }
2183: PetscFunctionReturn(PETSC_SUCCESS);
2184: }
2186: static PetscErrorCode MatMultTranspose_IS(Mat A, Vec y, Vec x)
2187: {
2188: Mat_IS *is = (Mat_IS *)A->data;
2190: PetscFunctionBegin;
2191: /* scatter the global vector x into the local work vector */
2192: PetscCall(VecScatterBegin(is->rctx, y, is->y, INSERT_VALUES, SCATTER_FORWARD));
2193: PetscCall(VecScatterEnd(is->rctx, y, is->y, INSERT_VALUES, SCATTER_FORWARD));
2195: /* multiply the local matrix */
2196: PetscCall(MatMultTranspose(is->A, is->y, is->x));
2198: /* scatter product back into global vector */
2199: PetscCall(VecSet(x, 0));
2200: PetscCall(VecScatterBegin(is->cctx, is->x, x, ADD_VALUES, SCATTER_REVERSE));
2201: PetscCall(VecScatterEnd(is->cctx, is->x, x, ADD_VALUES, SCATTER_REVERSE));
2202: PetscFunctionReturn(PETSC_SUCCESS);
2203: }
2205: static PetscErrorCode MatMultTransposeAdd_IS(Mat A, Vec v1, Vec v2, Vec v3)
2206: {
2207: Vec temp_vec;
2209: PetscFunctionBegin; /* v3 = v2 + A' * v1.*/
2210: if (v3 != v2) {
2211: PetscCall(MatMultTranspose(A, v1, v3));
2212: PetscCall(VecAXPY(v3, 1.0, v2));
2213: } else {
2214: PetscCall(VecDuplicate(v2, &temp_vec));
2215: PetscCall(MatMultTranspose(A, v1, temp_vec));
2216: PetscCall(VecAXPY(temp_vec, 1.0, v2));
2217: PetscCall(VecCopy(temp_vec, v3));
2218: PetscCall(VecDestroy(&temp_vec));
2219: }
2220: PetscFunctionReturn(PETSC_SUCCESS);
2221: }
2223: static PetscErrorCode ISLocalToGlobalMappingView_Multi(ISLocalToGlobalMapping mapping, PetscInt lsize, PetscInt gsize, const PetscInt vblocks[], PetscViewer viewer)
2224: {
2225: PetscInt tr[3], n;
2226: const PetscInt *indices;
2228: PetscFunctionBegin;
2229: tr[0] = IS_LTOGM_FILE_CLASSID;
2230: tr[1] = 1;
2231: tr[2] = gsize;
2232: PetscCall(PetscViewerBinaryWrite(viewer, tr, 3, PETSC_INT));
2233: PetscCall(PetscViewerBinaryWriteAll(viewer, vblocks, lsize, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
2234: PetscCall(ISLocalToGlobalMappingGetSize(mapping, &n));
2235: PetscCall(ISLocalToGlobalMappingGetIndices(mapping, &indices));
2236: PetscCall(PetscViewerBinaryWriteAll(viewer, indices, n, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
2237: PetscCall(ISLocalToGlobalMappingRestoreIndices(mapping, &indices));
2238: PetscFunctionReturn(PETSC_SUCCESS);
2239: }
2241: static PetscErrorCode MatView_IS(Mat A, PetscViewer viewer)
2242: {
2243: Mat_IS *a = (Mat_IS *)A->data;
2244: PetscViewer sviewer;
2245: PetscBool isascii, isbinary, viewl2g = PETSC_FALSE, native;
2246: PetscViewerFormat format;
2247: ISLocalToGlobalMapping rmap, cmap;
2249: PetscFunctionBegin;
2250: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2251: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
2252: PetscCall(PetscViewerGetFormat(viewer, &format));
2253: native = (PetscBool)(format == PETSC_VIEWER_NATIVE);
2254: if (native) {
2255: rmap = A->rmap->mapping;
2256: cmap = A->cmap->mapping;
2257: } else {
2258: rmap = a->rmapping;
2259: cmap = a->cmapping;
2260: }
2261: if (isascii) {
2262: if (format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(PETSC_SUCCESS);
2263: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL || format == PETSC_VIEWER_ASCII_MATLAB) viewl2g = PETSC_TRUE;
2264: } else if (isbinary) {
2265: PetscInt tr[6], nr, nc, lsize = 0;
2266: char lmattype[64] = {'\0'};
2267: PetscMPIInt size;
2268: PetscBool skipHeader, vbs = PETSC_FALSE;
2269: IS is;
2270: const PetscInt *vblocks = NULL;
2272: PetscCall(PetscViewerSetUp(viewer));
2273: PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)A)->prefix, "-mat_is_view_variableblocksizes", &vbs, NULL));
2274: if (vbs) {
2275: PetscCall(MatGetVariableBlockSizes(a->A, &lsize, &vblocks));
2276: PetscCall(PetscMPIIntCast(lsize, &size));
2277: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &size, 1, MPI_INT, MPI_SUM, PetscObjectComm((PetscObject)viewer)));
2278: } else {
2279: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)viewer), &size));
2280: }
2281: tr[0] = MAT_FILE_CLASSID;
2282: tr[1] = A->rmap->N;
2283: tr[2] = A->cmap->N;
2284: tr[3] = -size; /* AIJ stores nnz here */
2285: tr[4] = (PetscInt)(rmap == cmap);
2286: tr[5] = a->allow_repeated;
2287: PetscCall(PetscSNPrintf(lmattype, sizeof(lmattype), "%s", a->lmattype));
2289: PetscCall(PetscViewerBinaryWrite(viewer, tr, PETSC_STATIC_ARRAY_LENGTH(tr), PETSC_INT));
2290: PetscCall(PetscViewerBinaryWrite(viewer, lmattype, sizeof(lmattype), PETSC_CHAR));
2292: /* first dump l2g info (we need the header for proper loading on different number of processes) */
2293: PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));
2294: PetscCall(PetscViewerBinarySetSkipHeader(viewer, PETSC_FALSE));
2295: if (vbs) {
2296: PetscCall(ISLocalToGlobalMappingView_Multi(rmap, lsize, size, vblocks, viewer));
2297: if (cmap != rmap) PetscCall(ISLocalToGlobalMappingView_Multi(cmap, lsize, size, vblocks, viewer));
2298: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)viewer), lsize, vblocks, PETSC_USE_POINTER, &is));
2299: PetscCall(ISView(is, viewer));
2300: PetscCall(ISView(is, viewer));
2301: PetscCall(ISDestroy(&is));
2302: } else {
2303: PetscCall(ISLocalToGlobalMappingView(rmap, viewer));
2304: if (cmap != rmap) PetscCall(ISLocalToGlobalMappingView(cmap, viewer));
2306: /* then the sizes of the local matrices */
2307: PetscCall(MatGetSize(a->A, &nr, &nc));
2308: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)viewer), 1, &nr, PETSC_USE_POINTER, &is));
2309: PetscCall(ISView(is, viewer));
2310: PetscCall(ISDestroy(&is));
2311: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)viewer), 1, &nc, PETSC_USE_POINTER, &is));
2312: PetscCall(ISView(is, viewer));
2313: PetscCall(ISDestroy(&is));
2314: }
2315: PetscCall(PetscViewerBinarySetSkipHeader(viewer, skipHeader));
2316: }
2317: if (format == PETSC_VIEWER_ASCII_MATLAB) {
2318: char name[64];
2319: PetscMPIInt size, rank;
2321: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)viewer), &size));
2322: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
2323: if (size > 1) PetscCall(PetscSNPrintf(name, sizeof(name), "lmat_%d", rank));
2324: else PetscCall(PetscSNPrintf(name, sizeof(name), "lmat"));
2325: PetscCall(PetscObjectSetName((PetscObject)a->A, name));
2326: }
2328: /* Dump the local matrices */
2329: if (isbinary) { /* ViewerGetSubViewer does not work in parallel */
2330: PetscBool isaij;
2331: PetscInt nr, nc;
2332: Mat lA, B;
2333: Mat_MPIAIJ *b;
2335: /* We create a temporary MPIAIJ matrix that stores the unassembled operator */
2336: PetscCall(PetscObjectBaseTypeCompare((PetscObject)a->A, MATAIJ, &isaij));
2337: if (!isaij) PetscCall(MatConvert(a->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &lA));
2338: else {
2339: PetscCall(PetscObjectReference((PetscObject)a->A));
2340: lA = a->A;
2341: }
2342: PetscCall(MatCreate(PetscObjectComm((PetscObject)viewer), &B));
2343: PetscCall(MatSetType(B, MATMPIAIJ));
2344: PetscCall(MatGetSize(lA, &nr, &nc));
2345: PetscCall(MatSetSizes(B, nr, nc, PETSC_DECIDE, PETSC_DECIDE));
2346: PetscCall(MatMPIAIJSetPreallocation(B, 0, NULL, 0, NULL));
2348: b = (Mat_MPIAIJ *)B->data;
2349: PetscCall(MatDestroy(&b->A));
2350: b->A = lA;
2352: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
2353: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
2354: PetscCall(MatView(B, viewer));
2355: PetscCall(MatDestroy(&B));
2356: } else {
2357: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
2358: PetscCall(MatView(a->A, sviewer));
2359: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
2360: }
2362: /* with ASCII, we dump the l2gmaps at the end */
2363: if (viewl2g) {
2364: if (format == PETSC_VIEWER_ASCII_MATLAB) {
2365: PetscCall(PetscObjectSetName((PetscObject)rmap, "row"));
2366: PetscCall(ISLocalToGlobalMappingView(rmap, viewer));
2367: PetscCall(PetscObjectSetName((PetscObject)cmap, "col"));
2368: PetscCall(ISLocalToGlobalMappingView(cmap, viewer));
2369: } else {
2370: PetscCall(ISLocalToGlobalMappingView(rmap, viewer));
2371: PetscCall(ISLocalToGlobalMappingView(cmap, viewer));
2372: }
2373: }
2374: PetscFunctionReturn(PETSC_SUCCESS);
2375: }
2377: static PetscErrorCode ISLocalToGlobalMappingHasRepeatedLocal_Private(ISLocalToGlobalMapping map, PetscBool *has)
2378: {
2379: const PetscInt *idxs;
2380: PetscHSetI ht;
2381: PetscInt n, bs;
2383: PetscFunctionBegin;
2384: PetscCall(ISLocalToGlobalMappingGetSize(map, &n));
2385: PetscCall(ISLocalToGlobalMappingGetBlockSize(map, &bs));
2386: PetscCall(ISLocalToGlobalMappingGetBlockIndices(map, &idxs));
2387: PetscCall(PetscHSetICreate(&ht));
2388: *has = PETSC_FALSE;
2389: for (PetscInt i = 0; i < n / bs; i++) {
2390: PetscBool missing = PETSC_TRUE;
2391: if (idxs[i] < 0) continue;
2392: PetscCall(PetscHSetIQueryAdd(ht, idxs[i], &missing));
2393: if (!missing) {
2394: *has = PETSC_TRUE;
2395: break;
2396: }
2397: }
2398: PetscCall(PetscHSetIDestroy(&ht));
2399: PetscFunctionReturn(PETSC_SUCCESS);
2400: }
2402: static PetscErrorCode MatLoad_IS(Mat A, PetscViewer viewer)
2403: {
2404: ISLocalToGlobalMapping rmap, cmap;
2405: MPI_Comm comm = PetscObjectComm((PetscObject)A);
2406: PetscBool isbinary, samel, allow, isbaij;
2407: PetscInt tr[6], M, N, nr, nc, Asize, isn;
2408: const PetscInt *idx;
2409: PetscMPIInt size;
2410: char lmattype[64];
2411: Mat dA, lA;
2412: IS is;
2414: PetscFunctionBegin;
2415: PetscCheckSameComm(A, 1, viewer, 2);
2416: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
2417: PetscCheck(isbinary, PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Invalid viewer of type %s", ((PetscObject)viewer)->type_name);
2419: PetscCall(PetscViewerBinaryRead(viewer, tr, PETSC_STATIC_ARRAY_LENGTH(tr), NULL, PETSC_INT));
2420: PetscCheck(tr[0] == MAT_FILE_CLASSID, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a matrix next in file");
2421: PetscCheck(tr[1] >= 0, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a IS matrix next in file");
2422: PetscCheck(tr[2] >= 0, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a IS matrix next in file");
2423: PetscCheck(tr[3] < 0, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a IS matrix next in file");
2424: PetscCheck(tr[4] == 0 || tr[4] == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a IS matrix next in file");
2425: PetscCheck(tr[5] == 0 || tr[5] == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a IS matrix next in file");
2426: M = tr[1];
2427: N = tr[2];
2428: Asize = -tr[3];
2429: samel = (PetscBool)tr[4];
2430: allow = (PetscBool)tr[5];
2431: PetscCall(PetscViewerBinaryRead(viewer, lmattype, sizeof(lmattype), NULL, PETSC_CHAR));
2433: /* if we are loading from a larger set of processes, allow repeated entries */
2434: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)viewer), &size));
2435: if (Asize > size) allow = PETSC_TRUE;
2437: /* set global sizes if not set already */
2438: if (A->rmap->N < 0) A->rmap->N = M;
2439: if (A->cmap->N < 0) A->cmap->N = N;
2440: PetscCall(PetscLayoutSetUp(A->rmap));
2441: PetscCall(PetscLayoutSetUp(A->cmap));
2442: PetscCheck(M == A->rmap->N, comm, PETSC_ERR_ARG_SIZ, "Matrix rows should be %" PetscInt_FMT ", found %" PetscInt_FMT, M, A->rmap->N);
2443: PetscCheck(N == A->cmap->N, comm, PETSC_ERR_ARG_SIZ, "Matrix columns should be %" PetscInt_FMT ", found %" PetscInt_FMT, N, A->cmap->N);
2445: /* load l2g maps */
2446: PetscCall(ISLocalToGlobalMappingCreate(comm, 0, 0, NULL, PETSC_USE_POINTER, &rmap));
2447: PetscCall(ISLocalToGlobalMappingLoad(rmap, viewer));
2448: if (!samel) {
2449: PetscCall(ISLocalToGlobalMappingCreate(comm, 0, 0, NULL, PETSC_USE_POINTER, &cmap));
2450: PetscCall(ISLocalToGlobalMappingLoad(cmap, viewer));
2451: } else {
2452: PetscCall(PetscObjectReference((PetscObject)rmap));
2453: cmap = rmap;
2454: }
2456: /* load sizes of local matrices */
2457: PetscCall(ISCreate(comm, &is));
2458: PetscCall(ISSetType(is, ISGENERAL));
2459: PetscCall(ISLoad(is, viewer));
2460: PetscCall(ISGetLocalSize(is, &isn));
2461: PetscCall(ISGetIndices(is, &idx));
2462: nr = 0;
2463: for (PetscInt i = 0; i < isn; i++) nr += idx[i];
2464: PetscCall(ISRestoreIndices(is, &idx));
2465: PetscCall(ISDestroy(&is));
2466: PetscCall(ISCreate(comm, &is));
2467: PetscCall(ISSetType(is, ISGENERAL));
2468: PetscCall(ISLoad(is, viewer));
2469: PetscCall(ISGetLocalSize(is, &isn));
2470: PetscCall(ISGetIndices(is, &idx));
2471: nc = 0;
2472: for (PetscInt i = 0; i < isn; i++) nc += idx[i];
2473: PetscCall(ISRestoreIndices(is, &idx));
2474: PetscCall(ISDestroy(&is));
2476: /* now load the unassembled operator */
2477: PetscCall(MatCreate(comm, &dA));
2478: PetscCall(MatSetType(dA, MATMPIAIJ));
2479: PetscCall(MatSetSizes(dA, nr, nc, PETSC_DECIDE, PETSC_DECIDE));
2480: PetscCall(MatLoad(dA, viewer));
2481: PetscCall(MatMPIAIJGetSeqAIJ(dA, &lA, NULL, NULL));
2482: PetscCall(PetscObjectReference((PetscObject)lA));
2483: PetscCall(MatDestroy(&dA));
2485: /* and convert to the desired format */
2486: PetscCall(PetscStrcmpAny(lmattype, &isbaij, MATSBAIJ, MATSEQSBAIJ, ""));
2487: if (isbaij) PetscCall(MatSetOption(lA, MAT_SYMMETRIC, PETSC_TRUE));
2488: PetscCall(MatConvert(lA, lmattype, MAT_INPLACE_MATRIX, &lA));
2490: /* check if we actually have repeated entries */
2491: if (allow) {
2492: PetscBool rhas, chas, hasrepeated;
2494: PetscCall(ISLocalToGlobalMappingHasRepeatedLocal_Private(rmap, &rhas));
2495: if (rmap != cmap) PetscCall(ISLocalToGlobalMappingHasRepeatedLocal_Private(cmap, &chas));
2496: else chas = rhas;
2497: hasrepeated = (PetscBool)(rhas || chas);
2498: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &hasrepeated, 1, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A)));
2499: if (!hasrepeated) allow = PETSC_FALSE;
2500: }
2502: /* assemble the MATIS object */
2503: PetscCall(MatISSetAllowRepeated(A, allow));
2504: PetscCall(MatSetLocalToGlobalMapping(A, rmap, cmap));
2505: PetscCall(MatISSetLocalMat(A, lA));
2506: PetscCall(MatDestroy(&lA));
2507: PetscCall(ISLocalToGlobalMappingDestroy(&rmap));
2508: PetscCall(ISLocalToGlobalMappingDestroy(&cmap));
2509: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
2510: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
2511: PetscFunctionReturn(PETSC_SUCCESS);
2512: }
2514: static PetscErrorCode MatInvertBlockDiagonal_IS(Mat mat, const PetscScalar **values)
2515: {
2516: Mat_IS *is = (Mat_IS *)mat->data;
2517: MPI_Datatype nodeType;
2518: const PetscScalar *lv;
2519: PetscInt bs;
2520: PetscMPIInt mbs;
2522: PetscFunctionBegin;
2523: PetscCall(MatGetBlockSize(mat, &bs));
2524: PetscCall(MatSetBlockSize(is->A, bs));
2525: PetscCall(MatInvertBlockDiagonal(is->A, &lv));
2526: if (!is->bdiag) PetscCall(PetscMalloc1(bs * mat->rmap->n, &is->bdiag));
2527: PetscCall(PetscMPIIntCast(bs, &mbs));
2528: PetscCallMPI(MPI_Type_contiguous(mbs, MPIU_SCALAR, &nodeType));
2529: PetscCallMPI(MPI_Type_commit(&nodeType));
2530: PetscCall(PetscSFReduceBegin(is->sf, nodeType, lv, is->bdiag, MPI_REPLACE));
2531: PetscCall(PetscSFReduceEnd(is->sf, nodeType, lv, is->bdiag, MPI_REPLACE));
2532: PetscCallMPI(MPI_Type_free(&nodeType));
2533: if (values) *values = is->bdiag;
2534: PetscFunctionReturn(PETSC_SUCCESS);
2535: }
2537: static PetscErrorCode MatISSetUpScatters_Private(Mat A)
2538: {
2539: Vec cglobal, rglobal;
2540: IS from;
2541: Mat_IS *is = (Mat_IS *)A->data;
2542: PetscScalar sum;
2543: const PetscInt *garray;
2544: PetscInt nr, rbs, nc, cbs;
2545: VecType rtype;
2547: PetscFunctionBegin;
2548: PetscCall(ISLocalToGlobalMappingGetSize(is->rmapping, &nr));
2549: PetscCall(ISLocalToGlobalMappingGetBlockSize(is->rmapping, &rbs));
2550: PetscCall(ISLocalToGlobalMappingGetSize(is->cmapping, &nc));
2551: PetscCall(ISLocalToGlobalMappingGetBlockSize(is->cmapping, &cbs));
2552: PetscCall(VecDestroy(&is->x));
2553: PetscCall(VecDestroy(&is->y));
2554: PetscCall(VecDestroy(&is->counter));
2555: PetscCall(VecScatterDestroy(&is->rctx));
2556: PetscCall(VecScatterDestroy(&is->cctx));
2557: PetscCall(MatCreateVecs(is->A, &is->x, &is->y));
2558: PetscCall(VecBindToCPU(is->y, PETSC_TRUE));
2559: PetscCall(VecGetRootType_Private(is->y, &rtype));
2560: PetscCall(PetscFree(A->defaultvectype));
2561: PetscCall(PetscStrallocpy(rtype, &A->defaultvectype));
2562: PetscCall(MatCreateVecs(A, &cglobal, &rglobal));
2563: PetscCall(VecBindToCPU(rglobal, PETSC_TRUE));
2564: PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->rmapping, &garray));
2565: PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), rbs, nr / rbs, garray, PETSC_USE_POINTER, &from));
2566: PetscCall(VecScatterCreate(rglobal, from, is->y, NULL, &is->rctx));
2567: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->rmapping, &garray));
2568: PetscCall(ISDestroy(&from));
2569: if (is->rmapping != is->cmapping) {
2570: PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->cmapping, &garray));
2571: PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), cbs, nc / cbs, garray, PETSC_USE_POINTER, &from));
2572: PetscCall(VecScatterCreate(cglobal, from, is->x, NULL, &is->cctx));
2573: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->cmapping, &garray));
2574: PetscCall(ISDestroy(&from));
2575: } else {
2576: PetscCall(PetscObjectReference((PetscObject)is->rctx));
2577: is->cctx = is->rctx;
2578: }
2579: PetscCall(VecDestroy(&cglobal));
2581: /* interface counter vector (local) */
2582: PetscCall(VecDuplicate(is->y, &is->counter));
2583: PetscCall(VecBindToCPU(is->counter, PETSC_TRUE));
2584: PetscCall(VecSet(is->y, 1.));
2585: PetscCall(VecScatterBegin(is->rctx, is->y, rglobal, ADD_VALUES, SCATTER_REVERSE));
2586: PetscCall(VecScatterEnd(is->rctx, is->y, rglobal, ADD_VALUES, SCATTER_REVERSE));
2587: PetscCall(VecScatterBegin(is->rctx, rglobal, is->counter, INSERT_VALUES, SCATTER_FORWARD));
2588: PetscCall(VecScatterEnd(is->rctx, rglobal, is->counter, INSERT_VALUES, SCATTER_FORWARD));
2589: PetscCall(VecBindToCPU(is->y, PETSC_FALSE));
2590: PetscCall(VecBindToCPU(is->counter, PETSC_FALSE));
2592: /* special functions for block-diagonal matrices */
2593: PetscCall(VecSum(rglobal, &sum));
2594: A->ops->invertblockdiagonal = NULL;
2595: if ((PetscInt)(PetscRealPart(sum)) == A->rmap->N && A->rmap->N == A->cmap->N && is->rmapping == is->cmapping) A->ops->invertblockdiagonal = MatInvertBlockDiagonal_IS;
2596: PetscCall(VecDestroy(&rglobal));
2598: /* setup SF for general purpose shared indices based communications */
2599: PetscCall(MatISSetUpSF_IS(A));
2600: PetscFunctionReturn(PETSC_SUCCESS);
2601: }
2603: static PetscErrorCode MatISFilterL2GMap(Mat A, ISLocalToGlobalMapping map, ISLocalToGlobalMapping *nmap, ISLocalToGlobalMapping *lmap)
2604: {
2605: Mat_IS *matis = (Mat_IS *)A->data;
2606: IS is;
2607: ISLocalToGlobalMappingType l2gtype;
2608: const PetscInt *idxs;
2609: PetscHSetI ht;
2610: PetscInt *nidxs;
2611: PetscInt i, n, bs, c;
2612: PetscBool flg[] = {PETSC_FALSE, PETSC_FALSE};
2614: PetscFunctionBegin;
2615: PetscCall(ISLocalToGlobalMappingGetSize(map, &n));
2616: PetscCall(ISLocalToGlobalMappingGetBlockSize(map, &bs));
2617: PetscCall(ISLocalToGlobalMappingGetBlockIndices(map, &idxs));
2618: PetscCall(PetscHSetICreate(&ht));
2619: PetscCall(PetscMalloc1(n / bs, &nidxs));
2620: for (i = 0, c = 0; i < n / bs; i++) {
2621: PetscBool missing = PETSC_TRUE;
2622: if (idxs[i] < 0) {
2623: flg[0] = PETSC_TRUE;
2624: continue;
2625: }
2626: if (!matis->allow_repeated) PetscCall(PetscHSetIQueryAdd(ht, idxs[i], &missing));
2627: if (!missing) flg[1] = PETSC_TRUE;
2628: else nidxs[c++] = idxs[i];
2629: }
2630: PetscCall(PetscHSetIDestroy(&ht));
2631: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, flg, 2, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A)));
2632: if (!flg[0] && !flg[1]) { /* Entries are all non negative and unique */
2633: *nmap = NULL;
2634: *lmap = NULL;
2635: PetscCall(PetscFree(nidxs));
2636: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(map, &idxs));
2637: PetscFunctionReturn(PETSC_SUCCESS);
2638: }
2640: /* New l2g map without negative indices (and repeated indices if not allowed) */
2641: PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), bs, c, nidxs, PETSC_USE_POINTER, &is));
2642: PetscCall(ISLocalToGlobalMappingCreateIS(is, nmap));
2643: PetscCall(ISDestroy(&is));
2644: PetscCall(ISLocalToGlobalMappingGetType(map, &l2gtype));
2645: PetscCall(ISLocalToGlobalMappingSetType(*nmap, l2gtype));
2647: /* New local l2g map for repeated indices if not allowed */
2648: PetscCall(ISGlobalToLocalMappingApplyBlock(*nmap, IS_GTOLM_MASK, n / bs, idxs, NULL, nidxs));
2649: PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, n / bs, nidxs, PETSC_USE_POINTER, &is));
2650: PetscCall(ISLocalToGlobalMappingCreateIS(is, lmap));
2651: PetscCall(ISDestroy(&is));
2652: PetscCall(PetscFree(nidxs));
2653: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(map, &idxs));
2654: PetscFunctionReturn(PETSC_SUCCESS);
2655: }
2657: static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A, ISLocalToGlobalMapping rmapping, ISLocalToGlobalMapping cmapping)
2658: {
2659: Mat_IS *is = (Mat_IS *)A->data;
2660: ISLocalToGlobalMapping localrmapping = NULL, localcmapping = NULL;
2661: PetscInt nr, rbs, nc, cbs;
2662: PetscBool cong, freem[] = {PETSC_FALSE, PETSC_FALSE};
2664: PetscFunctionBegin;
2665: if (rmapping) PetscCheckSameComm(A, 1, rmapping, 2);
2666: if (cmapping) PetscCheckSameComm(A, 1, cmapping, 3);
2668: PetscCall(ISLocalToGlobalMappingDestroy(&is->rmapping));
2669: PetscCall(ISLocalToGlobalMappingDestroy(&is->cmapping));
2670: PetscCall(PetscLayoutSetUp(A->rmap));
2671: PetscCall(PetscLayoutSetUp(A->cmap));
2672: PetscCall(MatHasCongruentLayouts(A, &cong));
2674: /* If NULL, local space matches global space */
2675: if (!rmapping) {
2676: IS is;
2678: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), A->rmap->N, 0, 1, &is));
2679: PetscCall(ISLocalToGlobalMappingCreateIS(is, &rmapping));
2680: PetscCall(ISLocalToGlobalMappingSetBlockSize(rmapping, A->rmap->bs));
2681: PetscCall(ISDestroy(&is));
2682: freem[0] = PETSC_TRUE;
2683: if (!cmapping && cong && A->rmap->bs == A->cmap->bs) cmapping = rmapping;
2684: } else if (!is->islocalref) { /* check if the l2g map has negative or repeated entries */
2685: PetscCall(MatISFilterL2GMap(A, rmapping, &is->rmapping, &localrmapping));
2686: if (rmapping == cmapping) {
2687: PetscCall(PetscObjectReference((PetscObject)is->rmapping));
2688: is->cmapping = is->rmapping;
2689: PetscCall(PetscObjectReference((PetscObject)localrmapping));
2690: localcmapping = localrmapping;
2691: }
2692: }
2693: if (!cmapping) {
2694: IS is;
2696: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), A->cmap->N, 0, 1, &is));
2697: PetscCall(ISLocalToGlobalMappingCreateIS(is, &cmapping));
2698: PetscCall(ISLocalToGlobalMappingSetBlockSize(cmapping, A->cmap->bs));
2699: PetscCall(ISDestroy(&is));
2700: freem[1] = PETSC_TRUE;
2701: } else if (cmapping != rmapping && !is->islocalref) { /* check if the l2g map has negative or repeated entries */
2702: PetscCall(MatISFilterL2GMap(A, cmapping, &is->cmapping, &localcmapping));
2703: }
2704: if (!is->rmapping) {
2705: PetscCall(PetscObjectReference((PetscObject)rmapping));
2706: is->rmapping = rmapping;
2707: }
2708: if (!is->cmapping) {
2709: PetscCall(PetscObjectReference((PetscObject)cmapping));
2710: is->cmapping = cmapping;
2711: }
2713: /* Clean up */
2714: is->lnnzstate = 0;
2715: PetscCall(MatDestroy(&is->dA));
2716: PetscCall(MatDestroy(&is->assembledA));
2717: PetscCall(MatDestroy(&is->A));
2718: if (is->csf != is->sf) {
2719: PetscCall(PetscSFDestroy(&is->csf));
2720: PetscCall(PetscFree2(is->csf_rootdata, is->csf_leafdata));
2721: } else is->csf = NULL;
2722: PetscCall(PetscSFDestroy(&is->sf));
2723: PetscCall(PetscFree2(is->sf_rootdata, is->sf_leafdata));
2724: PetscCall(PetscFree(is->bdiag));
2726: /* check if the two mappings are actually the same for square matrices since MATIS has some optimization for this case
2727: (DOLFIN passes 2 different objects) */
2728: PetscCall(ISLocalToGlobalMappingGetSize(is->rmapping, &nr));
2729: PetscCall(ISLocalToGlobalMappingGetBlockSize(is->rmapping, &rbs));
2730: PetscCall(ISLocalToGlobalMappingGetSize(is->cmapping, &nc));
2731: PetscCall(ISLocalToGlobalMappingGetBlockSize(is->cmapping, &cbs));
2732: if (is->rmapping != is->cmapping && cong) {
2733: PetscBool same = PETSC_FALSE;
2734: if (nr == nc && cbs == rbs) {
2735: const PetscInt *idxs1, *idxs2;
2737: PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->rmapping, &idxs1));
2738: PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->cmapping, &idxs2));
2739: PetscCall(PetscArraycmp(idxs1, idxs2, nr / rbs, &same));
2740: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->rmapping, &idxs1));
2741: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->cmapping, &idxs2));
2742: }
2743: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &same, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
2744: if (same) {
2745: PetscCall(ISLocalToGlobalMappingDestroy(&is->cmapping));
2746: PetscCall(PetscObjectReference((PetscObject)is->rmapping));
2747: is->cmapping = is->rmapping;
2748: }
2749: }
2750: PetscCall(PetscLayoutSetBlockSize(A->rmap, rbs));
2751: PetscCall(PetscLayoutSetBlockSize(A->cmap, cbs));
2752: /* Pass the user defined maps to the layout */
2753: PetscCall(PetscLayoutSetISLocalToGlobalMapping(A->rmap, rmapping));
2754: PetscCall(PetscLayoutSetISLocalToGlobalMapping(A->cmap, cmapping));
2755: if (freem[0]) PetscCall(ISLocalToGlobalMappingDestroy(&rmapping));
2756: if (freem[1]) PetscCall(ISLocalToGlobalMappingDestroy(&cmapping));
2758: if (!is->islocalref) {
2759: /* Create the local matrix A */
2760: PetscCall(MatCreate(PETSC_COMM_SELF, &is->A));
2761: PetscCall(MatSetType(is->A, is->lmattype));
2762: PetscCall(MatSetSizes(is->A, nr, nc, nr, nc));
2763: PetscCall(MatSetBlockSizes(is->A, rbs, cbs));
2764: PetscCall(MatSetOptionsPrefix(is->A, "is_"));
2765: PetscCall(MatAppendOptionsPrefix(is->A, ((PetscObject)A)->prefix));
2766: PetscCall(PetscLayoutSetUp(is->A->rmap));
2767: PetscCall(PetscLayoutSetUp(is->A->cmap));
2768: PetscCall(MatSetLocalToGlobalMapping(is->A, localrmapping, localcmapping));
2769: PetscCall(ISLocalToGlobalMappingDestroy(&localrmapping));
2770: PetscCall(ISLocalToGlobalMappingDestroy(&localcmapping));
2771: /* setup scatters and local vectors for MatMult */
2772: PetscCall(MatISSetUpScatters_Private(A));
2773: }
2774: A->preallocated = PETSC_TRUE;
2775: PetscFunctionReturn(PETSC_SUCCESS);
2776: }
2778: static PetscErrorCode MatSetUp_IS(Mat A)
2779: {
2780: Mat_IS *is = (Mat_IS *)A->data;
2781: ISLocalToGlobalMapping rmap, cmap;
2783: PetscFunctionBegin;
2784: if (!is->sf) {
2785: PetscCall(MatGetLocalToGlobalMapping(A, &rmap, &cmap));
2786: PetscCall(MatSetLocalToGlobalMapping(A, rmap, cmap));
2787: }
2788: PetscFunctionReturn(PETSC_SUCCESS);
2789: }
2791: static PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
2792: {
2793: Mat_IS *is = (Mat_IS *)mat->data;
2794: PetscInt buf[2 * MATIS_MAX_ENTRIES_INSERTION], *rows_l = NULL, *cols_l = NULL;
2796: PetscFunctionBegin;
2797: IndexSpaceGet(buf, m, n, rows_l, cols_l);
2798: PetscCall(ISGlobalToLocalMappingApply(is->rmapping, IS_GTOLM_MASK, m, rows, &m, rows_l));
2799: if (m != n || rows != cols || is->cmapping != is->rmapping) {
2800: PetscCall(ISGlobalToLocalMappingApply(is->cmapping, IS_GTOLM_MASK, n, cols, &n, cols_l));
2801: PetscCall(MatSetValues(is->A, m, rows_l, n, cols_l, values, addv));
2802: } else {
2803: PetscCall(MatSetValues(is->A, m, rows_l, m, rows_l, values, addv));
2804: }
2805: IndexSpaceRestore(buf, m, n, rows_l, cols_l);
2806: PetscFunctionReturn(PETSC_SUCCESS);
2807: }
2809: static PetscErrorCode MatSetValuesBlocked_IS(Mat mat, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
2810: {
2811: Mat_IS *is = (Mat_IS *)mat->data;
2812: PetscInt buf[2 * MATIS_MAX_ENTRIES_INSERTION], *rows_l = NULL, *cols_l = NULL;
2814: PetscFunctionBegin;
2815: IndexSpaceGet(buf, m, n, rows_l, cols_l);
2816: PetscCall(ISGlobalToLocalMappingApplyBlock(is->rmapping, IS_GTOLM_MASK, m, rows, &m, rows_l));
2817: if (m != n || rows != cols || is->cmapping != is->rmapping) {
2818: PetscCall(ISGlobalToLocalMappingApplyBlock(is->cmapping, IS_GTOLM_MASK, n, cols, &n, cols_l));
2819: PetscCall(MatSetValuesBlocked(is->A, m, rows_l, n, cols_l, values, addv));
2820: } else {
2821: PetscCall(MatSetValuesBlocked(is->A, m, rows_l, m, rows_l, values, addv));
2822: }
2823: IndexSpaceRestore(buf, m, n, rows_l, cols_l);
2824: PetscFunctionReturn(PETSC_SUCCESS);
2825: }
2827: static PetscErrorCode MatSetValuesLocal_IS(Mat A, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
2828: {
2829: Mat_IS *is = (Mat_IS *)A->data;
2831: PetscFunctionBegin;
2832: if (is->A->rmap->mapping || is->A->cmap->mapping) {
2833: PetscCall(MatSetValuesLocal(is->A, m, rows, n, cols, values, addv));
2834: } else {
2835: PetscCall(MatSetValues(is->A, m, rows, n, cols, values, addv));
2836: }
2837: PetscFunctionReturn(PETSC_SUCCESS);
2838: }
2840: static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat A, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
2841: {
2842: Mat_IS *is = (Mat_IS *)A->data;
2844: PetscFunctionBegin;
2845: if (is->A->rmap->mapping || is->A->cmap->mapping) {
2846: PetscCall(MatSetValuesBlockedLocal(is->A, m, rows, n, cols, values, addv));
2847: } else {
2848: PetscCall(MatSetValuesBlocked(is->A, m, rows, n, cols, values, addv));
2849: }
2850: PetscFunctionReturn(PETSC_SUCCESS);
2851: }
2853: static PetscErrorCode MatISZeroRowsColumnsLocal_Private(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, PetscBool columns)
2854: {
2855: Mat_IS *is = (Mat_IS *)A->data;
2857: PetscFunctionBegin;
2858: if (!n) PetscFunctionReturn(PETSC_SUCCESS);
2859: is->pure_neumann = PETSC_FALSE;
2861: if (columns) {
2862: PetscCall(MatZeroRowsColumns(is->A, n, rows, diag, NULL, NULL));
2863: } else {
2864: PetscCall(MatZeroRows(is->A, n, rows, diag, NULL, NULL));
2865: }
2866: if (diag != 0.) {
2867: const PetscScalar *array;
2869: PetscCall(VecGetArrayRead(is->counter, &array));
2870: for (PetscInt i = 0; i < n; i++) PetscCall(MatSetValue(is->A, rows[i], rows[i], diag / (array[rows[i]]), INSERT_VALUES));
2871: PetscCall(VecRestoreArrayRead(is->counter, &array));
2872: PetscCall(MatAssemblyBegin(is->A, MAT_FINAL_ASSEMBLY));
2873: PetscCall(MatAssemblyEnd(is->A, MAT_FINAL_ASSEMBLY));
2874: }
2875: PetscFunctionReturn(PETSC_SUCCESS);
2876: }
2878: static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b, PetscBool columns)
2879: {
2880: Mat_IS *matis = (Mat_IS *)A->data;
2881: PetscInt nr, nl, len;
2882: PetscInt *lrows = NULL;
2884: PetscFunctionBegin;
2885: if (PetscUnlikelyDebug(columns || diag != 0. || (x && b))) {
2886: PetscBool cong;
2888: PetscCall(PetscLayoutCompare(A->rmap, A->cmap, &cong));
2889: cong = (PetscBool)(cong && matis->sf == matis->csf);
2890: PetscCheck(cong || !columns, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Columns can be zeroed if and only if A->rmap and A->cmap are congruent and the l2g maps are the same for MATIS");
2891: PetscCheck(cong || diag == 0., PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Nonzero diagonal value supported if and only if A->rmap and A->cmap are congruent and the l2g maps are the same for MATIS");
2892: PetscCheck(cong || !x || !b, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "A->rmap and A->cmap need to be congruent, and the l2g maps be the same");
2893: }
2894: PetscCall(MatGetSize(matis->A, &nl, NULL));
2895: /* get locally owned rows */
2896: PetscCall(PetscLayoutMapLocal(A->rmap, n, rows, &len, &lrows, NULL));
2897: /* fix right-hand side if needed */
2898: if (x && b) {
2899: const PetscScalar *xx;
2900: PetscScalar *bb;
2902: if (columns) {
2903: /* Subtract the column contributions: b[i] -= A[i,r] * x[r] for non-zeroed rows i.
2904: Build x_zeroed = x restricted to the zeroed (locally owned) rows, 0 elsewhere,
2905: then compute b -= A * x_zeroed using the original (unmodified) local matrices.
2906: The zeroed rows of b are overwritten below with diag * x[r], so no special
2907: treatment is needed for them in the MatMult() output. */
2908: Vec x_zeroed, temp;
2909: PetscScalar *xz, *saved;
2911: PetscCall(PetscMalloc1(len, &saved));
2912: PetscCall(VecDuplicate(x, &x_zeroed));
2913: PetscCall(VecGetArrayRead(x, &xx));
2914: PetscCall(VecGetArray(x_zeroed, &xz));
2915: for (PetscInt i = 0; i < len; i++) {
2916: xz[lrows[i]] = xx[lrows[i]];
2917: saved[i] = diag * xx[lrows[i]];
2918: }
2919: PetscCall(VecRestoreArray(x_zeroed, &xz));
2920: PetscCall(VecRestoreArrayRead(x, &xx));
2921: PetscCall(VecDuplicate(b, &temp));
2922: PetscCall(MatMult(A, x_zeroed, temp));
2923: PetscCall(VecAXPY(b, -1.0, temp));
2924: PetscCall(VecDestroy(&temp));
2925: PetscCall(VecDestroy(&x_zeroed));
2926: /* Overwrite zeroed rows: b[r] = diag * x[r] (after the MatMult() so it is not clobbered) */
2927: PetscCall(VecGetArray(b, &bb));
2928: for (PetscInt i = 0; i < len; i++) bb[lrows[i]] = saved[i];
2929: PetscCall(VecRestoreArray(b, &bb));
2930: PetscCall(PetscFree(saved));
2931: } else {
2932: /* MatZeroRows(): only set b[r] = diag * x[r] for the zeroed rows */
2933: PetscCall(VecGetArrayRead(x, &xx));
2934: PetscCall(VecGetArray(b, &bb));
2935: for (PetscInt i = 0; i < len; ++i) bb[lrows[i]] = diag * xx[lrows[i]];
2936: PetscCall(VecRestoreArrayRead(x, &xx));
2937: PetscCall(VecRestoreArray(b, &bb));
2938: }
2939: }
2940: /* get rows associated to the local matrices */
2941: PetscCall(PetscArrayzero(matis->sf_leafdata, nl));
2942: PetscCall(PetscArrayzero(matis->sf_rootdata, A->rmap->n));
2943: for (PetscInt i = 0; i < len; i++) matis->sf_rootdata[lrows[i]] = 1;
2944: PetscCall(PetscFree(lrows));
2945: PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
2946: PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
2947: PetscCall(PetscMalloc1(nl, &lrows));
2948: nr = 0;
2949: for (PetscInt i = 0; i < nl; i++)
2950: if (matis->sf_leafdata[i]) lrows[nr++] = i;
2951: PetscCall(MatISZeroRowsColumnsLocal_Private(A, nr, lrows, diag, columns));
2952: PetscCall(PetscFree(lrows));
2953: PetscFunctionReturn(PETSC_SUCCESS);
2954: }
2956: static PetscErrorCode MatZeroRows_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
2957: {
2958: PetscFunctionBegin;
2959: PetscCall(MatZeroRowsColumns_Private_IS(A, n, rows, diag, x, b, PETSC_FALSE));
2960: PetscFunctionReturn(PETSC_SUCCESS);
2961: }
2963: static PetscErrorCode MatZeroRowsColumns_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
2964: {
2965: PetscFunctionBegin;
2966: PetscCall(MatZeroRowsColumns_Private_IS(A, n, rows, diag, x, b, PETSC_TRUE));
2967: PetscFunctionReturn(PETSC_SUCCESS);
2968: }
2970: static PetscErrorCode MatAssemblyBegin_IS(Mat A, MatAssemblyType type)
2971: {
2972: Mat_IS *is = (Mat_IS *)A->data;
2974: PetscFunctionBegin;
2975: PetscCall(MatAssemblyBegin(is->A, type));
2976: PetscFunctionReturn(PETSC_SUCCESS);
2977: }
2979: static PetscErrorCode MatAssemblyEnd_IS(Mat A, MatAssemblyType type)
2980: {
2981: Mat_IS *is = (Mat_IS *)A->data;
2982: PetscBool lnnz;
2984: PetscFunctionBegin;
2985: PetscCall(MatAssemblyEnd(is->A, type));
2986: /* fix for local empty rows/cols */
2987: if (is->locempty && type == MAT_FINAL_ASSEMBLY) {
2988: Mat newlA;
2989: ISLocalToGlobalMapping rl2g, cl2g;
2990: IS nzr, nzc;
2991: PetscInt nr, nc, nnzr, nnzc;
2992: PetscBool lnewl2g, newl2g;
2994: PetscCall(MatGetSize(is->A, &nr, &nc));
2995: PetscCall(MatFindNonzeroRowsOrCols_Basic(is->A, PETSC_FALSE, PETSC_SMALL, &nzr));
2996: if (!nzr) PetscCall(ISCreateStride(PetscObjectComm((PetscObject)is->A), nr, 0, 1, &nzr));
2997: PetscCall(MatFindNonzeroRowsOrCols_Basic(is->A, PETSC_TRUE, PETSC_SMALL, &nzc));
2998: if (!nzc) PetscCall(ISCreateStride(PetscObjectComm((PetscObject)is->A), nc, 0, 1, &nzc));
2999: PetscCall(ISGetSize(nzr, &nnzr));
3000: PetscCall(ISGetSize(nzc, &nnzc));
3001: if (nnzr != nr || nnzc != nc) { /* need new global l2g map */
3002: lnewl2g = PETSC_TRUE;
3003: PetscCallMPI(MPIU_Allreduce(&lnewl2g, &newl2g, 1, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A)));
3005: /* extract valid submatrix */
3006: PetscCall(MatCreateSubMatrix(is->A, nzr, nzc, MAT_INITIAL_MATRIX, &newlA));
3007: } else { /* local matrix fully populated */
3008: lnewl2g = PETSC_FALSE;
3009: PetscCallMPI(MPIU_Allreduce(&lnewl2g, &newl2g, 1, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A)));
3010: PetscCall(PetscObjectReference((PetscObject)is->A));
3011: newlA = is->A;
3012: }
3014: /* attach new global l2g map if needed */
3015: if (newl2g) {
3016: IS zr, zc;
3017: const PetscInt *ridxs, *cidxs, *zridxs, *zcidxs;
3018: PetscInt *nidxs, i;
3020: PetscCall(ISComplement(nzr, 0, nr, &zr));
3021: PetscCall(ISComplement(nzc, 0, nc, &zc));
3022: PetscCall(PetscMalloc1(PetscMax(nr, nc), &nidxs));
3023: PetscCall(ISLocalToGlobalMappingGetIndices(is->rmapping, &ridxs));
3024: PetscCall(ISLocalToGlobalMappingGetIndices(is->cmapping, &cidxs));
3025: PetscCall(ISGetIndices(zr, &zridxs));
3026: PetscCall(ISGetIndices(zc, &zcidxs));
3027: PetscCall(ISGetLocalSize(zr, &nnzr));
3028: PetscCall(ISGetLocalSize(zc, &nnzc));
3030: PetscCall(PetscArraycpy(nidxs, ridxs, nr));
3031: for (i = 0; i < nnzr; i++) nidxs[zridxs[i]] = -1;
3032: PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A), 1, nr, nidxs, PETSC_COPY_VALUES, &rl2g));
3033: PetscCall(PetscArraycpy(nidxs, cidxs, nc));
3034: for (i = 0; i < nnzc; i++) nidxs[zcidxs[i]] = -1;
3035: PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A), 1, nc, nidxs, PETSC_COPY_VALUES, &cl2g));
3037: PetscCall(ISRestoreIndices(zr, &zridxs));
3038: PetscCall(ISRestoreIndices(zc, &zcidxs));
3039: PetscCall(ISLocalToGlobalMappingRestoreIndices(is->rmapping, &ridxs));
3040: PetscCall(ISLocalToGlobalMappingRestoreIndices(is->cmapping, &cidxs));
3041: PetscCall(ISDestroy(&nzr));
3042: PetscCall(ISDestroy(&nzc));
3043: PetscCall(ISDestroy(&zr));
3044: PetscCall(ISDestroy(&zc));
3045: PetscCall(PetscFree(nidxs));
3046: PetscCall(MatSetLocalToGlobalMapping(A, rl2g, cl2g));
3047: PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
3048: PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
3049: }
3050: PetscCall(MatISSetLocalMat(A, newlA));
3051: PetscCall(MatDestroy(&newlA));
3052: PetscCall(ISDestroy(&nzr));
3053: PetscCall(ISDestroy(&nzc));
3054: is->locempty = PETSC_FALSE;
3055: }
3056: lnnz = (PetscBool)(is->A->nonzerostate == is->lnnzstate);
3057: is->lnnzstate = is->A->nonzerostate;
3058: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &lnnz, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
3059: if (!lnnz) A->nonzerostate++;
3060: PetscFunctionReturn(PETSC_SUCCESS);
3061: }
3063: static PetscErrorCode MatISGetLocalMat_IS(Mat mat, Mat *local)
3064: {
3065: Mat_IS *is = (Mat_IS *)mat->data;
3067: PetscFunctionBegin;
3068: *local = is->A;
3069: PetscFunctionReturn(PETSC_SUCCESS);
3070: }
3072: static PetscErrorCode MatISRestoreLocalMat_IS(Mat mat, Mat *local)
3073: {
3074: PetscFunctionBegin;
3075: *local = NULL;
3076: PetscFunctionReturn(PETSC_SUCCESS);
3077: }
3079: /*@
3080: MatISGetLocalMat - Gets the local matrix stored inside a `MATIS` matrix.
3082: Not Collective.
3084: Input Parameter:
3085: . mat - the matrix
3087: Output Parameter:
3088: . local - the local matrix
3090: Level: intermediate
3092: Notes:
3093: This can be called if you have precomputed the nonzero structure of the
3094: matrix and want to provide it to the inner matrix object to improve the performance
3095: of the `MatSetValues()` operation.
3097: Call `MatISRestoreLocalMat()` when finished with the local matrix.
3099: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatISRestoreLocalMat()`
3100: @*/
3101: PetscErrorCode MatISGetLocalMat(Mat mat, Mat *local)
3102: {
3103: PetscFunctionBegin;
3105: PetscAssertPointer(local, 2);
3106: PetscUseMethod(mat, "MatISGetLocalMat_C", (Mat, Mat *), (mat, local));
3107: PetscFunctionReturn(PETSC_SUCCESS);
3108: }
3110: /*@
3111: MatISRestoreLocalMat - Restores the local matrix obtained with `MatISGetLocalMat()`
3113: Not Collective.
3115: Input Parameters:
3116: + mat - the matrix
3117: - local - the local matrix
3119: Level: intermediate
3121: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatISGetLocalMat()`
3122: @*/
3123: PetscErrorCode MatISRestoreLocalMat(Mat mat, Mat *local)
3124: {
3125: PetscFunctionBegin;
3127: PetscAssertPointer(local, 2);
3128: PetscUseMethod(mat, "MatISRestoreLocalMat_C", (Mat, Mat *), (mat, local));
3129: PetscFunctionReturn(PETSC_SUCCESS);
3130: }
3132: static PetscErrorCode MatISSetLocalMatType_IS(Mat mat, MatType mtype)
3133: {
3134: Mat_IS *is = (Mat_IS *)mat->data;
3136: PetscFunctionBegin;
3137: if (is->A) PetscCall(MatSetType(is->A, mtype));
3138: PetscCall(PetscFree(is->lmattype));
3139: PetscCall(PetscStrallocpy(mtype, &is->lmattype));
3140: PetscFunctionReturn(PETSC_SUCCESS);
3141: }
3143: /*@
3144: MatISSetLocalMatType - Specifies the type of local matrix inside the `MATIS`
3146: Logically Collective.
3148: Input Parameters:
3149: + mat - the matrix
3150: - mtype - the local matrix type
3152: Level: intermediate
3154: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatSetType()`, `MatType`
3155: @*/
3156: PetscErrorCode MatISSetLocalMatType(Mat mat, MatType mtype)
3157: {
3158: PetscFunctionBegin;
3160: PetscUseMethod(mat, "MatISSetLocalMatType_C", (Mat, MatType), (mat, mtype));
3161: PetscFunctionReturn(PETSC_SUCCESS);
3162: }
3164: static PetscErrorCode MatISSetLocalMat_IS(Mat mat, Mat local)
3165: {
3166: Mat_IS *is = (Mat_IS *)mat->data;
3167: PetscInt nrows, ncols, orows, ocols;
3168: MatType mtype, otype;
3169: PetscBool sametype = PETSC_TRUE;
3171: PetscFunctionBegin;
3172: if (is->A && !is->islocalref) {
3173: PetscCall(MatGetSize(is->A, &orows, &ocols));
3174: PetscCall(MatGetSize(local, &nrows, &ncols));
3175: PetscCheck(orows == nrows && ocols == ncols, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local MATIS matrix should be of size %" PetscInt_FMT "x%" PetscInt_FMT " (passed a %" PetscInt_FMT "x%" PetscInt_FMT " matrix)", orows, ocols, nrows, ncols);
3176: PetscCall(MatGetType(local, &mtype));
3177: PetscCall(MatGetType(is->A, &otype));
3178: PetscCall(PetscStrcmp(mtype, otype, &sametype));
3179: }
3180: PetscCall(PetscObjectReference((PetscObject)local));
3181: PetscCall(MatDestroy(&is->A));
3182: is->A = local;
3183: PetscCall(MatGetType(is->A, &mtype));
3184: PetscCall(MatISSetLocalMatType(mat, mtype));
3185: if (!sametype && !is->islocalref) PetscCall(MatISSetUpScatters_Private(mat));
3186: is->lnnzstate = 0;
3187: PetscFunctionReturn(PETSC_SUCCESS);
3188: }
3190: /*@
3191: MatISSetLocalMat - Replace the local matrix stored inside a `MATIS` object.
3193: Not Collective
3195: Input Parameters:
3196: + mat - the matrix
3197: - local - the local matrix
3199: Level: intermediate
3201: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatISSetLocalMatType`, `MatISGetLocalMat()`
3202: @*/
3203: PetscErrorCode MatISSetLocalMat(Mat mat, Mat local)
3204: {
3205: PetscFunctionBegin;
3208: PetscUseMethod(mat, "MatISSetLocalMat_C", (Mat, Mat), (mat, local));
3209: PetscFunctionReturn(PETSC_SUCCESS);
3210: }
3212: static PetscErrorCode MatZeroEntries_IS(Mat A)
3213: {
3214: Mat_IS *a = (Mat_IS *)A->data;
3216: PetscFunctionBegin;
3217: PetscCall(MatZeroEntries(a->A));
3218: PetscFunctionReturn(PETSC_SUCCESS);
3219: }
3221: static PetscErrorCode MatScale_IS(Mat A, PetscScalar a)
3222: {
3223: Mat_IS *is = (Mat_IS *)A->data;
3225: PetscFunctionBegin;
3226: PetscCall(MatScale(is->A, a));
3227: PetscFunctionReturn(PETSC_SUCCESS);
3228: }
3230: static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v)
3231: {
3232: Mat_IS *is = (Mat_IS *)A->data;
3234: PetscFunctionBegin;
3235: /* get diagonal of the local matrix */
3236: PetscCall(MatGetDiagonal(is->A, is->y));
3238: /* scatter diagonal back into global vector */
3239: PetscCall(VecSet(v, 0));
3240: PetscCall(VecScatterBegin(is->rctx, is->y, v, ADD_VALUES, SCATTER_REVERSE));
3241: PetscCall(VecScatterEnd(is->rctx, is->y, v, ADD_VALUES, SCATTER_REVERSE));
3242: PetscFunctionReturn(PETSC_SUCCESS);
3243: }
3245: static PetscErrorCode MatSetOption_IS(Mat A, MatOption op, PetscBool flg)
3246: {
3247: Mat_IS *a = (Mat_IS *)A->data;
3249: PetscFunctionBegin;
3250: PetscCall(MatSetOption(a->A, op, flg));
3251: PetscFunctionReturn(PETSC_SUCCESS);
3252: }
3254: static PetscErrorCode MatAXPY_IS(Mat Y, PetscScalar a, Mat X, MatStructure str)
3255: {
3256: Mat_IS *y = (Mat_IS *)Y->data;
3257: Mat_IS *x;
3259: PetscFunctionBegin;
3260: if (PetscDefined(USE_DEBUG)) {
3261: PetscBool ismatis;
3262: PetscCall(PetscObjectTypeCompare((PetscObject)X, MATIS, &ismatis));
3263: PetscCheck(ismatis, PetscObjectComm((PetscObject)Y), PETSC_ERR_SUP, "Cannot call MatAXPY(Y,a,X,str) with X not of type MATIS");
3264: }
3265: x = (Mat_IS *)X->data;
3266: PetscCall(MatAXPY(y->A, a, x->A, str));
3267: PetscFunctionReturn(PETSC_SUCCESS);
3268: }
3270: static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A, IS row, IS col, Mat *submat)
3271: {
3272: Mat lA;
3273: Mat_IS *matis = (Mat_IS *)A->data;
3274: ISLocalToGlobalMapping rl2g, cl2g;
3275: IS is;
3276: const PetscInt *rg, *rl;
3277: PetscInt nrg, rbs, cbs;
3278: PetscInt N, M, nrl, i, *idxs;
3280: PetscFunctionBegin;
3281: PetscCall(ISGetBlockSize(row, &rbs));
3282: PetscCall(ISGetBlockSize(col, &cbs));
3283: PetscCall(ISLocalToGlobalMappingGetIndices(A->rmap->mapping, &rg));
3284: PetscCall(ISGetLocalSize(row, &nrl));
3285: PetscCall(ISGetIndices(row, &rl));
3286: PetscCall(ISLocalToGlobalMappingGetSize(A->rmap->mapping, &nrg));
3287: if (PetscDefined(USE_DEBUG)) {
3288: for (i = 0; i < nrl; i++) PetscCheck(rl[i] < nrg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Local row index %" PetscInt_FMT " -> %" PetscInt_FMT " greater than maximum possible %" PetscInt_FMT, i, rl[i], nrg);
3289: }
3290: if (nrg % rbs) nrg = rbs * (nrg / rbs + 1);
3291: PetscCall(PetscMalloc1(nrg, &idxs));
3292: /* map from [0,nrl) to row */
3293: for (i = 0; i < nrl; i++) idxs[i] = rl[i];
3294: for (i = nrl; i < nrg; i++) idxs[i] = -1;
3295: PetscCall(ISRestoreIndices(row, &rl));
3296: PetscCall(ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping, &rg));
3297: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A), nrg, idxs, PETSC_OWN_POINTER, &is));
3298: PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g));
3299: PetscCall(ISLocalToGlobalMappingSetBlockSize(rl2g, rbs));
3300: PetscCall(ISDestroy(&is));
3301: /* compute new l2g map for columns */
3302: if (col != row || matis->rmapping != matis->cmapping || matis->A->rmap->mapping != matis->A->cmap->mapping) {
3303: const PetscInt *cg, *cl;
3304: PetscInt ncg;
3305: PetscInt ncl;
3307: PetscCall(ISLocalToGlobalMappingGetIndices(A->cmap->mapping, &cg));
3308: PetscCall(ISGetLocalSize(col, &ncl));
3309: PetscCall(ISGetIndices(col, &cl));
3310: PetscCall(ISLocalToGlobalMappingGetSize(A->cmap->mapping, &ncg));
3311: if (PetscDefined(USE_DEBUG)) {
3312: for (i = 0; i < ncl; i++) PetscCheck(cl[i] < ncg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Local column index %" PetscInt_FMT " -> %" PetscInt_FMT " greater than maximum possible %" PetscInt_FMT, i, cl[i], ncg);
3313: }
3314: if (ncg % cbs) ncg = cbs * (ncg / cbs + 1);
3315: PetscCall(PetscMalloc1(ncg, &idxs));
3316: /* map from [0,ncl) to col */
3317: for (i = 0; i < ncl; i++) idxs[i] = cl[i];
3318: for (i = ncl; i < ncg; i++) idxs[i] = -1;
3319: PetscCall(ISRestoreIndices(col, &cl));
3320: PetscCall(ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping, &cg));
3321: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A), ncg, idxs, PETSC_OWN_POINTER, &is));
3322: PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g));
3323: PetscCall(ISLocalToGlobalMappingSetBlockSize(cl2g, cbs));
3324: PetscCall(ISDestroy(&is));
3325: } else {
3326: PetscCall(PetscObjectReference((PetscObject)rl2g));
3327: cl2g = rl2g;
3328: }
3330: /* create the MATIS submatrix */
3331: PetscCall(MatGetSize(A, &M, &N));
3332: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), submat));
3333: PetscCall(MatSetSizes(*submat, PETSC_DECIDE, PETSC_DECIDE, M, N));
3334: PetscCall(MatSetType(*submat, MATIS));
3335: matis = (Mat_IS *)((*submat)->data);
3336: matis->islocalref = A;
3337: PetscCall(MatSetLocalToGlobalMapping(*submat, rl2g, cl2g));
3338: PetscCall(MatISGetLocalMat(A, &lA));
3339: PetscCall(MatISSetLocalMat(*submat, lA));
3340: PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
3341: PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
3343: /* remove unsupported ops */
3344: PetscCall(PetscMemzero((*submat)->ops, sizeof(struct _MatOps)));
3345: (*submat)->ops->destroy = MatDestroy_IS;
3346: (*submat)->ops->setvalueslocal = MatSetValuesLocal_SubMat_IS;
3347: (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS;
3348: (*submat)->ops->zerorowslocal = MatZeroRowsLocal_SubMat_IS;
3349: (*submat)->ops->zerorowscolumnslocal = MatZeroRowsColumnsLocal_SubMat_IS;
3350: (*submat)->ops->getlocalsubmatrix = MatGetLocalSubMatrix_IS;
3351: PetscFunctionReturn(PETSC_SUCCESS);
3352: }
3354: static PetscErrorCode MatSetFromOptions_IS(Mat A, PetscOptionItems PetscOptionsObject)
3355: {
3356: Mat_IS *a = (Mat_IS *)A->data;
3357: char type[256];
3358: PetscBool flg;
3360: PetscFunctionBegin;
3361: PetscOptionsHeadBegin(PetscOptionsObject, "MATIS options");
3362: PetscCall(PetscOptionsDeprecated("-matis_keepassembled", "-mat_is_keepassembled", "3.21", NULL));
3363: PetscCall(PetscOptionsDeprecated("-matis_fixempty", "-mat_is_fixempty", "3.21", NULL));
3364: PetscCall(PetscOptionsDeprecated("-matis_storel2l", "-mat_is_storel2l", "3.21", NULL));
3365: PetscCall(PetscOptionsDeprecated("-matis_localmat_type", "-mat_is_localmat_type", "3.21", NULL));
3366: PetscCall(PetscOptionsBool("-mat_is_keepassembled", "Store an assembled version if needed", NULL, a->keepassembled, &a->keepassembled, NULL));
3367: PetscCall(PetscOptionsBool("-mat_is_fixempty", "Fix local matrices in case of empty local rows/columns", "MatISFixLocalEmpty", a->locempty, &a->locempty, NULL));
3368: PetscCall(PetscOptionsBool("-mat_is_storel2l", "Store local-to-local matrices generated from PtAP operations", "MatISStoreL2L", a->storel2l, &a->storel2l, NULL));
3369: PetscCall(PetscOptionsBool("-mat_is_allow_repeated", "Allow local repeated entries", "MatISSetAllowRepeated", a->allow_repeated, &a->allow_repeated, NULL));
3370: PetscCall(PetscOptionsFList("-mat_is_localmat_type", "Matrix type", "MatISSetLocalMatType", MatList, a->lmattype, type, 256, &flg));
3371: if (flg) PetscCall(MatISSetLocalMatType(A, type));
3372: if (a->A) PetscCall(MatSetFromOptions(a->A));
3373: PetscOptionsHeadEnd();
3374: PetscFunctionReturn(PETSC_SUCCESS);
3375: }
3377: /*@
3378: MatCreateIS - Creates a "process" unassembled matrix.
3380: Collective.
3382: Input Parameters:
3383: + comm - MPI communicator that will share the matrix
3384: . bs - block size of the matrix
3385: . m - local size of left vector used in matrix vector products
3386: . n - local size of right vector used in matrix vector products
3387: . M - global size of left vector used in matrix vector products
3388: . N - global size of right vector used in matrix vector products
3389: . rmap - local to global map for rows
3390: - cmap - local to global map for cols
3392: Output Parameter:
3393: . A - the resulting matrix
3395: Level: intermediate
3397: Notes:
3398: `m` and `n` are NOT related to the size of the map; they represent the size of the local parts of the distributed vectors
3399: used in `MatMult()` operations. The local sizes of `rmap` and `cmap` define the size of the local matrices.
3401: If `rmap` (`cmap`) is `NULL`, then the local row (column) spaces matches the global space.
3403: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatSetLocalToGlobalMapping()`
3404: @*/
3405: PetscErrorCode MatCreateIS(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt M, PetscInt N, ISLocalToGlobalMapping rmap, ISLocalToGlobalMapping cmap, Mat *A)
3406: {
3407: PetscFunctionBegin;
3408: PetscCall(MatCreate(comm, A));
3409: PetscCall(MatSetSizes(*A, m, n, M, N));
3410: if (bs > 0) PetscCall(MatSetBlockSize(*A, bs));
3411: PetscCall(MatSetType(*A, MATIS));
3412: PetscCall(MatSetLocalToGlobalMapping(*A, rmap, cmap));
3413: PetscFunctionReturn(PETSC_SUCCESS);
3414: }
3416: static PetscErrorCode MatHasOperation_IS(Mat A, MatOperation op, PetscBool *has)
3417: {
3418: Mat_IS *a = (Mat_IS *)A->data;
3419: MatOperation tobefiltered[] = {MATOP_MULT_ADD, MATOP_MULT_TRANSPOSE_ADD, MATOP_GET_DIAGONAL_BLOCK, MATOP_INCREASE_OVERLAP};
3421: PetscFunctionBegin;
3422: *has = PETSC_FALSE;
3423: if (!((void **)A->ops)[op] || !a->A) PetscFunctionReturn(PETSC_SUCCESS);
3424: *has = PETSC_TRUE;
3425: for (PetscInt i = 0; i < (PetscInt)PETSC_STATIC_ARRAY_LENGTH(tobefiltered); i++)
3426: if (op == tobefiltered[i]) PetscFunctionReturn(PETSC_SUCCESS);
3427: PetscCall(MatHasOperation(a->A, op, has));
3428: PetscFunctionReturn(PETSC_SUCCESS);
3429: }
3431: static PetscErrorCode MatSetValuesCOO_IS(Mat A, const PetscScalar v[], InsertMode imode)
3432: {
3433: Mat_IS *a = (Mat_IS *)A->data;
3435: PetscFunctionBegin;
3436: PetscCall(MatSetValuesCOO(a->A, v, imode));
3437: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
3438: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
3439: PetscFunctionReturn(PETSC_SUCCESS);
3440: }
3442: static PetscErrorCode MatSetPreallocationCOOLocal_IS(Mat A, PetscCount ncoo, PetscInt coo_i[], PetscInt coo_j[])
3443: {
3444: Mat_IS *a = (Mat_IS *)A->data;
3446: PetscFunctionBegin;
3447: PetscCheck(a->A, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to provide l2g map first via MatSetLocalToGlobalMapping");
3448: if (a->A->rmap->mapping || a->A->cmap->mapping) {
3449: PetscCall(MatSetPreallocationCOOLocal(a->A, ncoo, coo_i, coo_j));
3450: } else {
3451: PetscCall(MatSetPreallocationCOO(a->A, ncoo, coo_i, coo_j));
3452: }
3453: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", MatSetValuesCOO_IS));
3454: A->preallocated = PETSC_TRUE;
3455: PetscFunctionReturn(PETSC_SUCCESS);
3456: }
3458: static PetscErrorCode MatSetPreallocationCOO_IS(Mat A, PetscCount ncoo, PetscInt coo_i[], PetscInt coo_j[])
3459: {
3460: Mat_IS *a = (Mat_IS *)A->data;
3461: PetscInt ncoo_i;
3463: PetscFunctionBegin;
3464: PetscCheck(a->A, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to provide l2g map first via MatSetLocalToGlobalMapping");
3465: PetscCall(PetscIntCast(ncoo, &ncoo_i));
3466: PetscCall(ISGlobalToLocalMappingApply(a->rmapping, IS_GTOLM_MASK, ncoo_i, coo_i, NULL, coo_i));
3467: PetscCall(ISGlobalToLocalMappingApply(a->cmapping, IS_GTOLM_MASK, ncoo_i, coo_j, NULL, coo_j));
3468: PetscCall(MatSetPreallocationCOO(a->A, ncoo, coo_i, coo_j));
3469: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", MatSetValuesCOO_IS));
3470: A->preallocated = PETSC_TRUE;
3471: PetscFunctionReturn(PETSC_SUCCESS);
3472: }
3474: static PetscErrorCode MatISGetAssembled_Private(Mat A, Mat *tA)
3475: {
3476: Mat_IS *a = (Mat_IS *)A->data;
3477: PetscObjectState Astate, aAstate = PETSC_INT_MIN;
3478: PetscObjectState Annzstate, aAnnzstate = PETSC_INT_MIN;
3480: PetscFunctionBegin;
3481: PetscCall(PetscObjectStateGet((PetscObject)A, &Astate));
3482: Annzstate = A->nonzerostate;
3483: if (a->assembledA) {
3484: PetscCall(PetscObjectStateGet((PetscObject)a->assembledA, &aAstate));
3485: aAnnzstate = a->assembledA->nonzerostate;
3486: }
3487: if (aAnnzstate != Annzstate) PetscCall(MatDestroy(&a->assembledA));
3488: if (Astate != aAstate || !a->assembledA) {
3489: MatType aAtype;
3490: PetscMPIInt size;
3491: PetscInt rbs, cbs, bs;
3493: /* the assembled form is used as temporary storage for parallel operations
3494: like createsubmatrices and the like, do not waste device memory */
3495: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
3496: PetscCall(ISLocalToGlobalMappingGetBlockSize(a->cmapping, &cbs));
3497: PetscCall(ISLocalToGlobalMappingGetBlockSize(a->rmapping, &rbs));
3498: bs = rbs == cbs ? rbs : 1;
3499: if (a->assembledA) PetscCall(MatGetType(a->assembledA, &aAtype));
3500: else if (size > 1) aAtype = bs > 1 ? MATMPIBAIJ : MATMPIAIJ;
3501: else aAtype = bs > 1 ? MATSEQBAIJ : MATSEQAIJ;
3503: PetscCall(MatConvert(A, aAtype, a->assembledA ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX, &a->assembledA));
3504: PetscCall(PetscObjectStateSet((PetscObject)a->assembledA, Astate));
3505: a->assembledA->nonzerostate = Annzstate;
3506: }
3507: PetscCall(PetscObjectReference((PetscObject)a->assembledA));
3508: *tA = a->assembledA;
3509: if (!a->keepassembled) PetscCall(MatDestroy(&a->assembledA));
3510: PetscFunctionReturn(PETSC_SUCCESS);
3511: }
3513: static PetscErrorCode MatISRestoreAssembled_Private(Mat A, Mat *tA)
3514: {
3515: PetscFunctionBegin;
3516: PetscCall(MatDestroy(tA));
3517: *tA = NULL;
3518: PetscFunctionReturn(PETSC_SUCCESS);
3519: }
3521: static PetscErrorCode MatGetDiagonalBlock_IS(Mat A, Mat *dA)
3522: {
3523: Mat_IS *a = (Mat_IS *)A->data;
3524: PetscObjectState Astate, dAstate = PETSC_INT_MIN;
3526: PetscFunctionBegin;
3527: PetscCall(PetscObjectStateGet((PetscObject)A, &Astate));
3528: if (a->dA) PetscCall(PetscObjectStateGet((PetscObject)a->dA, &dAstate));
3529: if (Astate != dAstate) {
3530: Mat tA;
3531: MatType ltype;
3533: PetscCall(MatDestroy(&a->dA));
3534: PetscCall(MatISGetAssembled_Private(A, &tA));
3535: PetscCall(MatGetDiagonalBlock(tA, &a->dA));
3536: PetscCall(MatPropagateSymmetryOptions(tA, a->dA));
3537: PetscCall(MatGetType(a->A, <ype));
3538: PetscCall(MatConvert(a->dA, ltype, MAT_INPLACE_MATRIX, &a->dA));
3539: PetscCall(PetscObjectReference((PetscObject)a->dA));
3540: PetscCall(MatISRestoreAssembled_Private(A, &tA));
3541: PetscCall(PetscObjectStateSet((PetscObject)a->dA, Astate));
3542: }
3543: *dA = a->dA;
3544: PetscFunctionReturn(PETSC_SUCCESS);
3545: }
3547: static PetscErrorCode MatCreateSubMatrices_IS(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse reuse, Mat *submat[])
3548: {
3549: Mat tA;
3551: PetscFunctionBegin;
3552: PetscCall(MatISGetAssembled_Private(A, &tA));
3553: PetscCall(MatCreateSubMatrices(tA, n, irow, icol, reuse, submat));
3554: /* MatCreateSubMatrices_MPIAIJ is a mess at the moment */
3555: #if 0
3556: {
3557: Mat_IS *a = (Mat_IS*)A->data;
3558: MatType ltype;
3559: VecType vtype;
3560: char *flg;
3562: PetscCall(MatGetType(a->A,<ype));
3563: PetscCall(MatGetVecType(a->A,&vtype));
3564: PetscCall(PetscStrstr(vtype,"cuda",&flg));
3565: if (!flg) PetscCall(PetscStrstr(vtype,"hip",&flg));
3566: if (!flg) PetscCall(PetscStrstr(vtype,"kokkos",&flg));
3567: if (flg) {
3568: for (PetscInt i = 0; i < n; i++) {
3569: Mat sA = (*submat)[i];
3571: PetscCall(MatConvert(sA,ltype,MAT_INPLACE_MATRIX,&sA));
3572: (*submat)[i] = sA;
3573: }
3574: }
3575: }
3576: #endif
3577: PetscCall(MatISRestoreAssembled_Private(A, &tA));
3578: PetscFunctionReturn(PETSC_SUCCESS);
3579: }
3581: static PetscErrorCode MatIncreaseOverlap_IS(Mat A, PetscInt n, IS is[], PetscInt ov)
3582: {
3583: Mat tA;
3585: PetscFunctionBegin;
3586: PetscCall(MatISGetAssembled_Private(A, &tA));
3587: PetscCall(MatIncreaseOverlap(tA, n, is, ov));
3588: PetscCall(MatISRestoreAssembled_Private(A, &tA));
3589: PetscFunctionReturn(PETSC_SUCCESS);
3590: }
3592: /*@
3593: MatISGetLocalToGlobalMapping - Gets the local-to-global numbering of the `MATIS` object
3595: Not Collective
3597: Input Parameter:
3598: . A - the matrix
3600: Output Parameters:
3601: + rmapping - row mapping
3602: - cmapping - column mapping
3604: Level: advanced
3606: Note:
3607: The returned map can be different from the one used to construct the `MATIS` object, since it will not contain negative or repeated indices.
3609: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatSetLocalToGlobalMapping()`
3610: @*/
3611: PetscErrorCode MatISGetLocalToGlobalMapping(Mat A, ISLocalToGlobalMapping *rmapping, ISLocalToGlobalMapping *cmapping)
3612: {
3613: PetscFunctionBegin;
3616: if (rmapping) PetscAssertPointer(rmapping, 2);
3617: if (cmapping) PetscAssertPointer(cmapping, 3);
3618: PetscUseMethod(A, "MatISGetLocalToGlobalMapping_C", (Mat, ISLocalToGlobalMapping *, ISLocalToGlobalMapping *), (A, rmapping, cmapping));
3619: PetscFunctionReturn(PETSC_SUCCESS);
3620: }
3622: static PetscErrorCode MatISGetLocalToGlobalMapping_IS(Mat A, ISLocalToGlobalMapping *r, ISLocalToGlobalMapping *c)
3623: {
3624: Mat_IS *a = (Mat_IS *)A->data;
3626: PetscFunctionBegin;
3627: if (r) *r = a->rmapping;
3628: if (c) *c = a->cmapping;
3629: PetscFunctionReturn(PETSC_SUCCESS);
3630: }
3632: static PetscErrorCode MatSetBlockSizes_IS(Mat A, PetscInt rbs, PetscInt cbs)
3633: {
3634: Mat_IS *a = (Mat_IS *)A->data;
3636: PetscFunctionBegin;
3637: if (a->A) PetscCall(MatSetBlockSizes(a->A, rbs, cbs));
3638: PetscFunctionReturn(PETSC_SUCCESS);
3639: }
3641: /*MC
3642: MATIS - MATIS = "is" - A matrix type to be used for non-overlapping domain decomposition methods (e.g. `PCBDDC` or `KSPFETIDP`).
3643: This stores the matrices in globally unassembled form and the parallel matrix vector product is handled "implicitly".
3645: Options Database Keys:
3646: + -mat_type is - Set the matrix type to `MATIS`.
3647: . -mat_is_allow_repeated - Allow repeated entries in the local part of the local to global maps.
3648: . -mat_is_fixempty - Fix local matrices in case of empty local rows/columns.
3649: - -mat_is_storel2l - Store the local-to-local operators generated by the Galerkin process of `MatPtAP()`.
3651: Level: intermediate
3653: Notes:
3654: Options prefix for the inner matrix are given by `-is_mat_xxx`
3656: You must call `MatSetLocalToGlobalMapping()` before using this matrix type.
3658: You can do matrix preallocation on the local matrix after you obtain it with
3659: `MatISGetLocalMat()`; otherwise, you could use `MatISSetPreallocation()` or `MatXAIJSetPreallocation()`
3661: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatISGetLocalMat()`, `MatSetLocalToGlobalMapping()`, `MatISSetPreallocation()`, `MatCreateIS()`, `PCBDDC`, `KSPFETIDP`
3662: M*/
3663: PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A)
3664: {
3665: Mat_IS *a;
3667: PetscFunctionBegin;
3668: PetscCall(PetscNew(&a));
3669: PetscCall(PetscStrallocpy(MATAIJ, &a->lmattype));
3670: A->data = (void *)a;
3672: /* matrix ops */
3673: PetscCall(PetscMemzero(A->ops, sizeof(struct _MatOps)));
3674: A->ops->mult = MatMult_IS;
3675: A->ops->multadd = MatMultAdd_IS;
3676: A->ops->multtranspose = MatMultTranspose_IS;
3677: A->ops->multtransposeadd = MatMultTransposeAdd_IS;
3678: A->ops->destroy = MatDestroy_IS;
3679: A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS;
3680: A->ops->setvalues = MatSetValues_IS;
3681: A->ops->setvaluesblocked = MatSetValuesBlocked_IS;
3682: A->ops->setvalueslocal = MatSetValuesLocal_IS;
3683: A->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_IS;
3684: A->ops->zerorows = MatZeroRows_IS;
3685: A->ops->zerorowscolumns = MatZeroRowsColumns_IS;
3686: A->ops->assemblybegin = MatAssemblyBegin_IS;
3687: A->ops->assemblyend = MatAssemblyEnd_IS;
3688: A->ops->view = MatView_IS;
3689: A->ops->load = MatLoad_IS;
3690: A->ops->zeroentries = MatZeroEntries_IS;
3691: A->ops->scale = MatScale_IS;
3692: A->ops->getdiagonal = MatGetDiagonal_IS;
3693: A->ops->setoption = MatSetOption_IS;
3694: A->ops->ishermitian = MatIsHermitian_IS;
3695: A->ops->issymmetric = MatIsSymmetric_IS;
3696: A->ops->isstructurallysymmetric = MatIsStructurallySymmetric_IS;
3697: A->ops->duplicate = MatDuplicate_IS;
3698: A->ops->copy = MatCopy_IS;
3699: A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_IS;
3700: A->ops->createsubmatrix = MatCreateSubMatrix_IS;
3701: A->ops->axpy = MatAXPY_IS;
3702: A->ops->diagonalset = MatDiagonalSet_IS;
3703: A->ops->shift = MatShift_IS;
3704: A->ops->transpose = MatTranspose_IS;
3705: A->ops->getinfo = MatGetInfo_IS;
3706: A->ops->diagonalscale = MatDiagonalScale_IS;
3707: A->ops->setfromoptions = MatSetFromOptions_IS;
3708: A->ops->setup = MatSetUp_IS;
3709: A->ops->hasoperation = MatHasOperation_IS;
3710: A->ops->getdiagonalblock = MatGetDiagonalBlock_IS;
3711: A->ops->createsubmatrices = MatCreateSubMatrices_IS;
3712: A->ops->increaseoverlap = MatIncreaseOverlap_IS;
3713: A->ops->setblocksizes = MatSetBlockSizes_IS;
3715: /* special MATIS functions */
3716: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMatType_C", MatISSetLocalMatType_IS));
3717: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalMat_C", MatISGetLocalMat_IS));
3718: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISRestoreLocalMat_C", MatISRestoreLocalMat_IS));
3719: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMat_C", MatISSetLocalMat_IS));
3720: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetPreallocation_C", MatISSetPreallocation_IS));
3721: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetAllowRepeated_C", MatISSetAllowRepeated_IS));
3722: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISStoreL2L_C", MatISStoreL2L_IS));
3723: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISFixLocalEmpty_C", MatISFixLocalEmpty_IS));
3724: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalToGlobalMapping_C", MatISGetLocalToGlobalMapping_IS));
3725: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpiaij_C", MatConvert_IS_XAIJ));
3726: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpibaij_C", MatConvert_IS_XAIJ));
3727: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpisbaij_C", MatConvert_IS_XAIJ));
3728: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqaij_C", MatConvert_IS_XAIJ));
3729: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqbaij_C", MatConvert_IS_XAIJ));
3730: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqsbaij_C", MatConvert_IS_XAIJ));
3731: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_aij_C", MatConvert_IS_XAIJ));
3732: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOOLocal_C", MatSetPreallocationCOOLocal_IS));
3733: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOO_C", MatSetPreallocationCOO_IS));
3734: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATIS));
3735: PetscFunctionReturn(PETSC_SUCCESS);
3736: }