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 <petsc/private/sfimpl.h>
12: #include <petsc/private/vecimpl.h>
13: #include <petsc/private/hashseti.h>
15: #define MATIS_MAX_ENTRIES_INSERTION 2048
16: static PetscErrorCode MatSetValuesLocal_IS(Mat, PetscInt, const PetscInt *, PetscInt, const PetscInt *, const PetscScalar *, InsertMode);
17: static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat, PetscInt, const PetscInt *, PetscInt, const PetscInt *, const PetscScalar *, InsertMode);
18: static PetscErrorCode MatISSetUpScatters_Private(Mat);
20: static PetscErrorCode MatISContainerDestroyPtAP_Private(void *ptr)
21: {
22: MatISPtAP ptap = (MatISPtAP)ptr;
24: PetscFunctionBegin;
25: PetscCall(MatDestroySubMatrices(ptap->ris1 ? 2 : 1, &ptap->lP));
26: PetscCall(ISDestroy(&ptap->cis0));
27: PetscCall(ISDestroy(&ptap->cis1));
28: PetscCall(ISDestroy(&ptap->ris0));
29: PetscCall(ISDestroy(&ptap->ris1));
30: PetscCall(PetscFree(ptap));
31: PetscFunctionReturn(PETSC_SUCCESS);
32: }
34: static PetscErrorCode MatPtAPNumeric_IS_XAIJ(Mat A, Mat P, Mat C)
35: {
36: MatISPtAP ptap;
37: Mat_IS *matis = (Mat_IS *)A->data;
38: Mat lA, lC;
39: MatReuse reuse;
40: IS ris[2], cis[2];
41: PetscContainer c;
42: PetscInt n;
44: PetscFunctionBegin;
45: PetscCall(PetscObjectQuery((PetscObject)C, "_MatIS_PtAP", (PetscObject *)&c));
46: PetscCheck(c, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Missing PtAP information");
47: PetscCall(PetscContainerGetPointer(c, (void **)&ptap));
48: ris[0] = ptap->ris0;
49: ris[1] = ptap->ris1;
50: cis[0] = ptap->cis0;
51: cis[1] = ptap->cis1;
52: n = ptap->ris1 ? 2 : 1;
53: reuse = ptap->lP ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX;
54: PetscCall(MatCreateSubMatrices(P, n, ris, cis, reuse, &ptap->lP));
56: PetscCall(MatISGetLocalMat(A, &lA));
57: PetscCall(MatISGetLocalMat(C, &lC));
58: if (ptap->ris1) { /* unsymmetric A mapping */
59: Mat lPt;
61: PetscCall(MatTranspose(ptap->lP[1], MAT_INITIAL_MATRIX, &lPt));
62: PetscCall(MatMatMatMult(lPt, lA, ptap->lP[0], reuse, ptap->fill, &lC));
63: if (matis->storel2l) PetscCall(PetscObjectCompose((PetscObject)(A), "_MatIS_PtAP_l2l", (PetscObject)lPt));
64: PetscCall(MatDestroy(&lPt));
65: } else {
66: PetscCall(MatPtAP(lA, ptap->lP[0], reuse, ptap->fill, &lC));
67: if (matis->storel2l) PetscCall(PetscObjectCompose((PetscObject)C, "_MatIS_PtAP_l2l", (PetscObject)ptap->lP[0]));
68: }
69: if (reuse == MAT_INITIAL_MATRIX) {
70: PetscCall(MatISSetLocalMat(C, lC));
71: PetscCall(MatDestroy(&lC));
72: }
73: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
74: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
75: PetscFunctionReturn(PETSC_SUCCESS);
76: }
78: static PetscErrorCode MatGetNonzeroColumnsLocal_Private(Mat PT, IS *cis)
79: {
80: Mat Po, Pd;
81: IS zd, zo;
82: const PetscInt *garray;
83: PetscInt *aux, i, bs;
84: PetscInt dc, stc, oc, ctd, cto;
85: PetscBool ismpiaij, ismpibaij, isseqaij, isseqbaij;
86: MPI_Comm comm;
88: PetscFunctionBegin;
90: PetscAssertPointer(cis, 2);
91: PetscCall(PetscObjectGetComm((PetscObject)PT, &comm));
92: bs = 1;
93: PetscCall(PetscObjectBaseTypeCompare((PetscObject)PT, MATMPIAIJ, &ismpiaij));
94: PetscCall(PetscObjectBaseTypeCompare((PetscObject)PT, MATMPIBAIJ, &ismpibaij));
95: PetscCall(PetscObjectBaseTypeCompare((PetscObject)PT, MATSEQAIJ, &isseqaij));
96: PetscCall(PetscObjectTypeCompare((PetscObject)PT, MATSEQBAIJ, &isseqbaij));
97: if (isseqaij || isseqbaij) {
98: Pd = PT;
99: Po = NULL;
100: garray = NULL;
101: } else if (ismpiaij) {
102: PetscCall(MatMPIAIJGetSeqAIJ(PT, &Pd, &Po, &garray));
103: } else if (ismpibaij) {
104: PetscCall(MatMPIBAIJGetSeqBAIJ(PT, &Pd, &Po, &garray));
105: PetscCall(MatGetBlockSize(PT, &bs));
106: } else SETERRQ(comm, PETSC_ERR_SUP, "Not for matrix type %s", ((PetscObject)(PT))->type_name);
108: /* identify any null columns in Pd or Po */
109: /* We use a tolerance comparison since it may happen that, with geometric multigrid,
110: some of the columns are not really zero, but very close to */
111: zo = zd = NULL;
112: if (Po) PetscCall(MatFindNonzeroRowsOrCols_Basic(Po, PETSC_TRUE, PETSC_SMALL, &zo));
113: PetscCall(MatFindNonzeroRowsOrCols_Basic(Pd, PETSC_TRUE, PETSC_SMALL, &zd));
115: PetscCall(MatGetLocalSize(PT, NULL, &dc));
116: PetscCall(MatGetOwnershipRangeColumn(PT, &stc, NULL));
117: if (Po) PetscCall(MatGetLocalSize(Po, NULL, &oc));
118: else oc = 0;
119: PetscCall(PetscMalloc1((dc + oc) / bs, &aux));
120: if (zd) {
121: const PetscInt *idxs;
122: PetscInt nz;
124: /* this will throw an error if bs is not valid */
125: PetscCall(ISSetBlockSize(zd, bs));
126: PetscCall(ISGetLocalSize(zd, &nz));
127: PetscCall(ISGetIndices(zd, &idxs));
128: ctd = nz / bs;
129: for (i = 0; i < ctd; i++) aux[i] = (idxs[bs * i] + stc) / bs;
130: PetscCall(ISRestoreIndices(zd, &idxs));
131: } else {
132: ctd = dc / bs;
133: for (i = 0; i < ctd; i++) aux[i] = i + stc / bs;
134: }
135: if (zo) {
136: const PetscInt *idxs;
137: PetscInt nz;
139: /* this will throw an error if bs is not valid */
140: PetscCall(ISSetBlockSize(zo, bs));
141: PetscCall(ISGetLocalSize(zo, &nz));
142: PetscCall(ISGetIndices(zo, &idxs));
143: cto = nz / bs;
144: for (i = 0; i < cto; i++) aux[i + ctd] = garray[idxs[bs * i] / bs];
145: PetscCall(ISRestoreIndices(zo, &idxs));
146: } else {
147: cto = oc / bs;
148: for (i = 0; i < cto; i++) aux[i + ctd] = garray[i];
149: }
150: PetscCall(ISCreateBlock(comm, bs, ctd + cto, aux, PETSC_OWN_POINTER, cis));
151: PetscCall(ISDestroy(&zd));
152: PetscCall(ISDestroy(&zo));
153: PetscFunctionReturn(PETSC_SUCCESS);
154: }
156: static PetscErrorCode MatPtAPSymbolic_IS_XAIJ(Mat A, Mat P, PetscReal fill, Mat C)
157: {
158: Mat PT, lA;
159: MatISPtAP ptap;
160: ISLocalToGlobalMapping Crl2g, Ccl2g, rl2g, cl2g;
161: PetscContainer c;
162: MatType lmtype;
163: const PetscInt *garray;
164: PetscInt ibs, N, dc;
165: MPI_Comm comm;
167: PetscFunctionBegin;
168: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
169: PetscCall(MatSetType(C, MATIS));
170: PetscCall(MatISGetLocalMat(A, &lA));
171: PetscCall(MatGetType(lA, &lmtype));
172: PetscCall(MatISSetLocalMatType(C, lmtype));
173: PetscCall(MatGetSize(P, NULL, &N));
174: PetscCall(MatGetLocalSize(P, NULL, &dc));
175: PetscCall(MatSetSizes(C, dc, dc, N, N));
176: /* Not sure about this
177: PetscCall(MatGetBlockSizes(P,NULL,&ibs));
178: PetscCall(MatSetBlockSize(*C,ibs));
179: */
181: PetscCall(PetscNew(&ptap));
182: PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &c));
183: PetscCall(PetscContainerSetPointer(c, ptap));
184: PetscCall(PetscContainerSetUserDestroy(c, MatISContainerDestroyPtAP_Private));
185: PetscCall(PetscObjectCompose((PetscObject)C, "_MatIS_PtAP", (PetscObject)c));
186: PetscCall(PetscContainerDestroy(&c));
187: ptap->fill = fill;
189: PetscCall(MatISGetLocalToGlobalMapping(A, &rl2g, &cl2g));
191: PetscCall(ISLocalToGlobalMappingGetBlockSize(cl2g, &ibs));
192: PetscCall(ISLocalToGlobalMappingGetSize(cl2g, &N));
193: PetscCall(ISLocalToGlobalMappingGetBlockIndices(cl2g, &garray));
194: PetscCall(ISCreateBlock(comm, ibs, N / ibs, garray, PETSC_COPY_VALUES, &ptap->ris0));
195: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(cl2g, &garray));
197: PetscCall(MatCreateSubMatrix(P, ptap->ris0, NULL, MAT_INITIAL_MATRIX, &PT));
198: PetscCall(MatGetNonzeroColumnsLocal_Private(PT, &ptap->cis0));
199: PetscCall(ISLocalToGlobalMappingCreateIS(ptap->cis0, &Ccl2g));
200: PetscCall(MatDestroy(&PT));
202: Crl2g = NULL;
203: if (rl2g != cl2g) { /* unsymmetric A mapping */
204: PetscBool same, lsame = PETSC_FALSE;
205: PetscInt N1, ibs1;
207: PetscCall(ISLocalToGlobalMappingGetSize(rl2g, &N1));
208: PetscCall(ISLocalToGlobalMappingGetBlockSize(rl2g, &ibs1));
209: PetscCall(ISLocalToGlobalMappingGetBlockIndices(rl2g, &garray));
210: PetscCall(ISCreateBlock(comm, ibs, N / ibs, garray, PETSC_COPY_VALUES, &ptap->ris1));
211: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(rl2g, &garray));
212: if (ibs1 == ibs && N1 == N) { /* check if the l2gmaps are the same */
213: const PetscInt *i1, *i2;
215: PetscCall(ISBlockGetIndices(ptap->ris0, &i1));
216: PetscCall(ISBlockGetIndices(ptap->ris1, &i2));
217: PetscCall(PetscArraycmp(i1, i2, N, &lsame));
218: }
219: PetscCall(MPIU_Allreduce(&lsame, &same, 1, MPIU_BOOL, MPI_LAND, comm));
220: if (same) {
221: PetscCall(ISDestroy(&ptap->ris1));
222: } else {
223: PetscCall(MatCreateSubMatrix(P, ptap->ris1, NULL, MAT_INITIAL_MATRIX, &PT));
224: PetscCall(MatGetNonzeroColumnsLocal_Private(PT, &ptap->cis1));
225: PetscCall(ISLocalToGlobalMappingCreateIS(ptap->cis1, &Crl2g));
226: PetscCall(MatDestroy(&PT));
227: }
228: }
229: /* Not sure about this
230: if (!Crl2g) {
231: PetscCall(MatGetBlockSize(C,&ibs));
232: PetscCall(ISLocalToGlobalMappingSetBlockSize(Ccl2g,ibs));
233: }
234: */
235: PetscCall(MatSetLocalToGlobalMapping(C, Crl2g ? Crl2g : Ccl2g, Ccl2g));
236: PetscCall(ISLocalToGlobalMappingDestroy(&Crl2g));
237: PetscCall(ISLocalToGlobalMappingDestroy(&Ccl2g));
239: C->ops->ptapnumeric = MatPtAPNumeric_IS_XAIJ;
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: static PetscErrorCode MatProductSymbolic_PtAP_IS_XAIJ(Mat C)
244: {
245: Mat_Product *product = C->product;
246: Mat A = product->A, P = product->B;
247: PetscReal fill = product->fill;
249: PetscFunctionBegin;
250: PetscCall(MatPtAPSymbolic_IS_XAIJ(A, P, fill, C));
251: C->ops->productnumeric = MatProductNumeric_PtAP;
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }
255: static PetscErrorCode MatProductSetFromOptions_IS_XAIJ_PtAP(Mat C)
256: {
257: PetscFunctionBegin;
258: C->ops->productsymbolic = MatProductSymbolic_PtAP_IS_XAIJ;
259: PetscFunctionReturn(PETSC_SUCCESS);
260: }
262: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_IS_XAIJ(Mat C)
263: {
264: Mat_Product *product = C->product;
266: PetscFunctionBegin;
267: if (product->type == MATPRODUCT_PtAP) PetscCall(MatProductSetFromOptions_IS_XAIJ_PtAP(C));
268: PetscFunctionReturn(PETSC_SUCCESS);
269: }
271: static PetscErrorCode MatISContainerDestroyFields_Private(void *ptr)
272: {
273: MatISLocalFields lf = (MatISLocalFields)ptr;
274: PetscInt i;
276: PetscFunctionBegin;
277: for (i = 0; i < lf->nr; i++) PetscCall(ISDestroy(&lf->rf[i]));
278: for (i = 0; i < lf->nc; i++) PetscCall(ISDestroy(&lf->cf[i]));
279: PetscCall(PetscFree2(lf->rf, lf->cf));
280: PetscCall(PetscFree(lf));
281: PetscFunctionReturn(PETSC_SUCCESS);
282: }
284: static PetscErrorCode MatConvert_SeqXAIJ_IS(Mat A, MatType type, MatReuse reuse, Mat *newmat)
285: {
286: Mat B, lB;
288: PetscFunctionBegin;
289: if (reuse != MAT_REUSE_MATRIX) {
290: ISLocalToGlobalMapping rl2g, cl2g;
291: PetscInt bs;
292: IS is;
294: PetscCall(MatGetBlockSize(A, &bs));
295: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), A->rmap->n / bs, 0, 1, &is));
296: if (bs > 1) {
297: IS is2;
298: PetscInt i, *aux;
300: PetscCall(ISGetLocalSize(is, &i));
301: PetscCall(ISGetIndices(is, (const PetscInt **)&aux));
302: PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), bs, i, aux, PETSC_COPY_VALUES, &is2));
303: PetscCall(ISRestoreIndices(is, (const PetscInt **)&aux));
304: PetscCall(ISDestroy(&is));
305: is = is2;
306: }
307: PetscCall(ISSetIdentity(is));
308: PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g));
309: PetscCall(ISDestroy(&is));
310: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), A->cmap->n / bs, 0, 1, &is));
311: if (bs > 1) {
312: IS is2;
313: PetscInt i, *aux;
315: PetscCall(ISGetLocalSize(is, &i));
316: PetscCall(ISGetIndices(is, (const PetscInt **)&aux));
317: PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), bs, i, aux, PETSC_COPY_VALUES, &is2));
318: PetscCall(ISRestoreIndices(is, (const PetscInt **)&aux));
319: PetscCall(ISDestroy(&is));
320: is = is2;
321: }
322: PetscCall(ISSetIdentity(is));
323: PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g));
324: PetscCall(ISDestroy(&is));
325: PetscCall(MatCreateIS(PetscObjectComm((PetscObject)A), bs, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N, rl2g, cl2g, &B));
326: PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
327: PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
328: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &lB));
329: if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
330: } else {
331: B = *newmat;
332: PetscCall(PetscObjectReference((PetscObject)A));
333: lB = A;
334: }
335: PetscCall(MatISSetLocalMat(B, lB));
336: PetscCall(MatDestroy(&lB));
337: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
338: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
339: if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &B));
340: PetscFunctionReturn(PETSC_SUCCESS);
341: }
343: static PetscErrorCode MatISScaleDisassembling_Private(Mat A)
344: {
345: Mat_IS *matis = (Mat_IS *)(A->data);
346: PetscScalar *aa;
347: const PetscInt *ii, *jj;
348: PetscInt i, n, m;
349: PetscInt *ecount, **eneighs;
350: PetscBool flg;
352: PetscFunctionBegin;
353: PetscCall(MatGetRowIJ(matis->A, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg));
354: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure");
355: PetscCall(ISLocalToGlobalMappingGetNodeInfo(matis->rmapping, &n, &ecount, &eneighs));
356: PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected %" PetscInt_FMT " != %" PetscInt_FMT, m, n);
357: PetscCall(MatSeqAIJGetArray(matis->A, &aa));
358: for (i = 0; i < n; i++) {
359: if (ecount[i] > 1) {
360: PetscInt j;
362: for (j = ii[i]; j < ii[i + 1]; j++) {
363: PetscInt i2 = jj[j], p, p2;
364: PetscReal scal = 0.0;
366: for (p = 0; p < ecount[i]; p++) {
367: for (p2 = 0; p2 < ecount[i2]; p2++) {
368: if (eneighs[i][p] == eneighs[i2][p2]) {
369: scal += 1.0;
370: break;
371: }
372: }
373: }
374: if (scal) aa[j] /= scal;
375: }
376: }
377: }
378: PetscCall(ISLocalToGlobalMappingRestoreNodeInfo(matis->rmapping, &n, &ecount, &eneighs));
379: PetscCall(MatSeqAIJRestoreArray(matis->A, &aa));
380: PetscCall(MatRestoreRowIJ(matis->A, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg));
381: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot restore IJ structure");
382: PetscFunctionReturn(PETSC_SUCCESS);
383: }
385: typedef enum {
386: MAT_IS_DISASSEMBLE_L2G_NATURAL,
387: MAT_IS_DISASSEMBLE_L2G_MAT,
388: MAT_IS_DISASSEMBLE_L2G_ND
389: } MatISDisassemblel2gType;
391: static PetscErrorCode MatMPIXAIJComputeLocalToGlobalMapping_Private(Mat A, ISLocalToGlobalMapping *l2g)
392: {
393: Mat Ad, Ao;
394: IS is, ndmap, ndsub;
395: MPI_Comm comm;
396: const PetscInt *garray, *ndmapi;
397: PetscInt bs, i, cnt, nl, *ncount, *ndmapc;
398: PetscBool ismpiaij, ismpibaij;
399: const char *const MatISDisassemblel2gTypes[] = {"NATURAL", "MAT", "ND", "MatISDisassemblel2gType", "MAT_IS_DISASSEMBLE_L2G_", NULL};
400: MatISDisassemblel2gType mode = MAT_IS_DISASSEMBLE_L2G_NATURAL;
401: MatPartitioning part;
402: PetscSF sf;
403: PetscObject dm;
405: PetscFunctionBegin;
406: PetscOptionsBegin(PetscObjectComm((PetscObject)A), ((PetscObject)A)->prefix, "MatIS l2g disassembling options", "Mat");
407: 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));
408: PetscOptionsEnd();
409: if (mode == MAT_IS_DISASSEMBLE_L2G_MAT) {
410: PetscCall(MatGetLocalToGlobalMapping(A, l2g, NULL));
411: PetscFunctionReturn(PETSC_SUCCESS);
412: }
413: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
414: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij));
415: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIBAIJ, &ismpibaij));
416: PetscCall(MatGetBlockSize(A, &bs));
417: switch (mode) {
418: case MAT_IS_DISASSEMBLE_L2G_ND:
419: PetscCall(MatPartitioningCreate(comm, &part));
420: PetscCall(MatPartitioningSetAdjacency(part, A));
421: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)part, ((PetscObject)A)->prefix));
422: PetscCall(MatPartitioningSetFromOptions(part));
423: PetscCall(MatPartitioningApplyND(part, &ndmap));
424: PetscCall(MatPartitioningDestroy(&part));
425: PetscCall(ISBuildTwoSided(ndmap, NULL, &ndsub));
426: PetscCall(MatMPIAIJSetUseScalableIncreaseOverlap(A, PETSC_TRUE));
427: PetscCall(MatIncreaseOverlap(A, 1, &ndsub, 1));
428: PetscCall(ISLocalToGlobalMappingCreateIS(ndsub, l2g));
430: /* it may happen that a separator node is not properly shared */
431: PetscCall(ISLocalToGlobalMappingGetNodeInfo(*l2g, &nl, &ncount, NULL));
432: PetscCall(PetscSFCreate(comm, &sf));
433: PetscCall(ISLocalToGlobalMappingGetIndices(*l2g, &garray));
434: PetscCall(PetscSFSetGraphLayout(sf, A->rmap, nl, NULL, PETSC_OWN_POINTER, garray));
435: PetscCall(ISLocalToGlobalMappingRestoreIndices(*l2g, &garray));
436: PetscCall(PetscCalloc1(A->rmap->n, &ndmapc));
437: PetscCall(PetscSFReduceBegin(sf, MPIU_INT, ncount, ndmapc, MPI_REPLACE));
438: PetscCall(PetscSFReduceEnd(sf, MPIU_INT, ncount, ndmapc, MPI_REPLACE));
439: PetscCall(ISLocalToGlobalMappingRestoreNodeInfo(*l2g, NULL, &ncount, NULL));
440: PetscCall(ISGetIndices(ndmap, &ndmapi));
441: for (i = 0, cnt = 0; i < A->rmap->n; i++)
442: if (ndmapi[i] < 0 && ndmapc[i] < 2) cnt++;
444: PetscCall(MPIU_Allreduce(&cnt, &i, 1, MPIU_INT, MPI_MAX, comm));
445: if (i) { /* we detected isolated separator nodes */
446: Mat A2, A3;
447: IS *workis, is2;
448: PetscScalar *vals;
449: PetscInt gcnt = i, *dnz, *onz, j, *lndmapi;
450: ISLocalToGlobalMapping ll2g;
451: PetscBool flg;
452: const PetscInt *ii, *jj;
454: /* communicate global id of separators */
455: MatPreallocateBegin(comm, A->rmap->n, A->cmap->n, dnz, onz);
456: for (i = 0, cnt = 0; i < A->rmap->n; i++) dnz[i] = ndmapi[i] < 0 ? i + A->rmap->rstart : -1;
458: PetscCall(PetscMalloc1(nl, &lndmapi));
459: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, dnz, lndmapi, MPI_REPLACE));
461: /* compute adjacency of isolated separators node */
462: PetscCall(PetscMalloc1(gcnt, &workis));
463: for (i = 0, cnt = 0; i < A->rmap->n; i++) {
464: if (ndmapi[i] < 0 && ndmapc[i] < 2) PetscCall(ISCreateStride(comm, 1, i + A->rmap->rstart, 1, &workis[cnt++]));
465: }
466: for (i = cnt; i < gcnt; i++) PetscCall(ISCreateStride(comm, 0, 0, 1, &workis[i]));
467: for (i = 0; i < gcnt; i++) {
468: PetscCall(PetscObjectSetName((PetscObject)workis[i], "ISOLATED"));
469: PetscCall(ISViewFromOptions(workis[i], NULL, "-view_isolated_separators"));
470: }
472: /* no communications since all the ISes correspond to locally owned rows */
473: PetscCall(MatIncreaseOverlap(A, gcnt, workis, 1));
475: /* end communicate global id of separators */
476: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, dnz, lndmapi, MPI_REPLACE));
478: /* communicate new layers : create a matrix and transpose it */
479: PetscCall(PetscArrayzero(dnz, A->rmap->n));
480: PetscCall(PetscArrayzero(onz, A->rmap->n));
481: for (i = 0, j = 0; i < A->rmap->n; i++) {
482: if (ndmapi[i] < 0 && ndmapc[i] < 2) {
483: const PetscInt *idxs;
484: PetscInt s;
486: PetscCall(ISGetLocalSize(workis[j], &s));
487: PetscCall(ISGetIndices(workis[j], &idxs));
488: PetscCall(MatPreallocateSet(i + A->rmap->rstart, s, idxs, dnz, onz));
489: j++;
490: }
491: }
492: PetscCheck(j == cnt, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected local count %" PetscInt_FMT " != %" PetscInt_FMT, j, cnt);
494: for (i = 0; i < gcnt; i++) {
495: PetscCall(PetscObjectSetName((PetscObject)workis[i], "EXTENDED"));
496: PetscCall(ISViewFromOptions(workis[i], NULL, "-view_isolated_separators"));
497: }
499: for (i = 0, j = 0; i < A->rmap->n; i++) j = PetscMax(j, dnz[i] + onz[i]);
500: PetscCall(PetscMalloc1(j, &vals));
501: for (i = 0; i < j; i++) vals[i] = 1.0;
503: PetscCall(MatCreate(comm, &A2));
504: PetscCall(MatSetType(A2, MATMPIAIJ));
505: PetscCall(MatSetSizes(A2, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
506: PetscCall(MatMPIAIJSetPreallocation(A2, 0, dnz, 0, onz));
507: PetscCall(MatSetOption(A2, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
508: for (i = 0, j = 0; i < A2->rmap->n; i++) {
509: PetscInt row = i + A2->rmap->rstart, s = dnz[i] + onz[i];
510: const PetscInt *idxs;
512: if (s) {
513: PetscCall(ISGetIndices(workis[j], &idxs));
514: PetscCall(MatSetValues(A2, 1, &row, s, idxs, vals, INSERT_VALUES));
515: PetscCall(ISRestoreIndices(workis[j], &idxs));
516: j++;
517: }
518: }
519: PetscCheck(j == cnt, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected local count %" PetscInt_FMT " != %" PetscInt_FMT, j, cnt);
520: PetscCall(PetscFree(vals));
521: PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY));
522: PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY));
523: PetscCall(MatTranspose(A2, MAT_INPLACE_MATRIX, &A2));
525: /* extract submatrix corresponding to the coupling "owned separators" x "isolated separators" */
526: for (i = 0, j = 0; i < nl; i++)
527: if (lndmapi[i] >= 0) lndmapi[j++] = lndmapi[i];
528: PetscCall(ISCreateGeneral(comm, j, lndmapi, PETSC_USE_POINTER, &is));
529: PetscCall(MatMPIAIJGetLocalMatCondensed(A2, MAT_INITIAL_MATRIX, &is, NULL, &A3));
530: PetscCall(ISDestroy(&is));
531: PetscCall(MatDestroy(&A2));
533: /* extend local to global map to include connected isolated separators */
534: PetscCall(PetscObjectQuery((PetscObject)A3, "_petsc_GetLocalMatCondensed_iscol", (PetscObject *)&is));
535: PetscCheck(is, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing column map");
536: PetscCall(ISLocalToGlobalMappingCreateIS(is, &ll2g));
537: PetscCall(MatGetRowIJ(A3, 0, PETSC_FALSE, PETSC_FALSE, &i, &ii, &jj, &flg));
538: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure");
539: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ii[i], jj, PETSC_COPY_VALUES, &is));
540: PetscCall(MatRestoreRowIJ(A3, 0, PETSC_FALSE, PETSC_FALSE, &i, &ii, &jj, &flg));
541: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure");
542: PetscCall(ISLocalToGlobalMappingApplyIS(ll2g, is, &is2));
543: PetscCall(ISDestroy(&is));
544: PetscCall(ISLocalToGlobalMappingDestroy(&ll2g));
546: /* add new nodes to the local-to-global map */
547: PetscCall(ISLocalToGlobalMappingDestroy(l2g));
548: PetscCall(ISExpand(ndsub, is2, &is));
549: PetscCall(ISDestroy(&is2));
550: PetscCall(ISLocalToGlobalMappingCreateIS(is, l2g));
551: PetscCall(ISDestroy(&is));
553: PetscCall(MatDestroy(&A3));
554: PetscCall(PetscFree(lndmapi));
555: MatPreallocateEnd(dnz, onz);
556: for (i = 0; i < gcnt; i++) PetscCall(ISDestroy(&workis[i]));
557: PetscCall(PetscFree(workis));
558: }
559: PetscCall(ISRestoreIndices(ndmap, &ndmapi));
560: PetscCall(PetscSFDestroy(&sf));
561: PetscCall(PetscFree(ndmapc));
562: PetscCall(ISDestroy(&ndmap));
563: PetscCall(ISDestroy(&ndsub));
564: PetscCall(ISLocalToGlobalMappingSetBlockSize(*l2g, bs));
565: PetscCall(ISLocalToGlobalMappingViewFromOptions(*l2g, NULL, "-matis_nd_l2g_view"));
566: break;
567: case MAT_IS_DISASSEMBLE_L2G_NATURAL:
568: PetscCall(PetscObjectQuery((PetscObject)A, "__PETSc_dm", (PetscObject *)&dm));
569: if (dm) { /* if a matrix comes from a DM, most likely we can use the l2gmap if any */
570: PetscCall(MatGetLocalToGlobalMapping(A, l2g, NULL));
571: PetscCall(PetscObjectReference((PetscObject)*l2g));
572: if (*l2g) PetscFunctionReturn(PETSC_SUCCESS);
573: }
574: if (ismpiaij) {
575: PetscCall(MatMPIAIJGetSeqAIJ(A, &Ad, &Ao, &garray));
576: } else if (ismpibaij) {
577: PetscCall(MatMPIBAIJGetSeqBAIJ(A, &Ad, &Ao, &garray));
578: } else SETERRQ(comm, PETSC_ERR_SUP, "Type %s", ((PetscObject)A)->type_name);
579: if (A->rmap->n) {
580: PetscInt dc, oc, stc, *aux;
582: PetscCall(MatGetLocalSize(Ad, NULL, &dc));
583: PetscCall(MatGetLocalSize(Ao, NULL, &oc));
584: PetscCheck(!oc || garray, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "garray not present");
585: PetscCall(MatGetOwnershipRangeColumn(A, &stc, NULL));
586: PetscCall(PetscMalloc1((dc + oc) / bs, &aux));
587: for (i = 0; i < dc / bs; i++) aux[i] = i + stc / bs;
588: for (i = 0; i < oc / bs; i++) aux[i + dc / bs] = (ismpiaij ? garray[i * bs] / bs : garray[i]);
589: PetscCall(ISCreateBlock(comm, bs, (dc + oc) / bs, aux, PETSC_OWN_POINTER, &is));
590: } else {
591: PetscCall(ISCreateBlock(comm, 1, 0, NULL, PETSC_OWN_POINTER, &is));
592: }
593: PetscCall(ISLocalToGlobalMappingCreateIS(is, l2g));
594: PetscCall(ISDestroy(&is));
595: break;
596: default:
597: SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unsupported l2g disassembling type %d", mode);
598: }
599: PetscFunctionReturn(PETSC_SUCCESS);
600: }
602: PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat A, MatType type, MatReuse reuse, Mat *newmat)
603: {
604: Mat lA, Ad, Ao, B = NULL;
605: ISLocalToGlobalMapping rl2g, cl2g;
606: IS is;
607: MPI_Comm comm;
608: void *ptrs[2];
609: const char *names[2] = {"_convert_csr_aux", "_convert_csr_data"};
610: const PetscInt *garray;
611: PetscScalar *dd, *od, *aa, *data;
612: const PetscInt *di, *dj, *oi, *oj;
613: const PetscInt *odi, *odj, *ooi, *ooj;
614: PetscInt *aux, *ii, *jj;
615: PetscInt bs, lc, dr, dc, oc, str, stc, nnz, i, jd, jo, cum;
616: PetscBool flg, ismpiaij, ismpibaij, was_inplace = PETSC_FALSE;
617: PetscMPIInt size;
619: PetscFunctionBegin;
620: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
621: PetscCallMPI(MPI_Comm_size(comm, &size));
622: if (size == 1) {
623: PetscCall(MatConvert_SeqXAIJ_IS(A, type, reuse, newmat));
624: PetscFunctionReturn(PETSC_SUCCESS);
625: }
626: if (reuse != MAT_REUSE_MATRIX && A->cmap->N == A->rmap->N) {
627: PetscCall(MatMPIXAIJComputeLocalToGlobalMapping_Private(A, &rl2g));
628: PetscCall(MatCreate(comm, &B));
629: PetscCall(MatSetType(B, MATIS));
630: PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
631: PetscCall(MatSetLocalToGlobalMapping(B, rl2g, rl2g));
632: PetscCall(MatGetBlockSize(A, &bs));
633: PetscCall(MatSetBlockSize(B, bs));
634: PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
635: if (reuse == MAT_INPLACE_MATRIX) was_inplace = PETSC_TRUE;
636: reuse = MAT_REUSE_MATRIX;
637: }
638: if (reuse == MAT_REUSE_MATRIX) {
639: Mat *newlA, lA;
640: IS rows, cols;
641: const PetscInt *ridx, *cidx;
642: PetscInt rbs, cbs, nr, nc;
644: if (!B) B = *newmat;
645: PetscCall(MatISGetLocalToGlobalMapping(B, &rl2g, &cl2g));
646: PetscCall(ISLocalToGlobalMappingGetBlockIndices(rl2g, &ridx));
647: PetscCall(ISLocalToGlobalMappingGetBlockIndices(cl2g, &cidx));
648: PetscCall(ISLocalToGlobalMappingGetSize(rl2g, &nr));
649: PetscCall(ISLocalToGlobalMappingGetSize(cl2g, &nc));
650: PetscCall(ISLocalToGlobalMappingGetBlockSize(rl2g, &rbs));
651: PetscCall(ISLocalToGlobalMappingGetBlockSize(cl2g, &cbs));
652: PetscCall(ISCreateBlock(comm, rbs, nr / rbs, ridx, PETSC_USE_POINTER, &rows));
653: if (rl2g != cl2g) {
654: PetscCall(ISCreateBlock(comm, cbs, nc / cbs, cidx, PETSC_USE_POINTER, &cols));
655: } else {
656: PetscCall(PetscObjectReference((PetscObject)rows));
657: cols = rows;
658: }
659: PetscCall(MatISGetLocalMat(B, &lA));
660: PetscCall(MatCreateSubMatrices(A, 1, &rows, &cols, MAT_INITIAL_MATRIX, &newlA));
661: PetscCall(MatConvert(newlA[0], MATSEQAIJ, MAT_INPLACE_MATRIX, &newlA[0]));
662: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(rl2g, &ridx));
663: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(cl2g, &cidx));
664: PetscCall(ISDestroy(&rows));
665: PetscCall(ISDestroy(&cols));
666: if (!lA->preallocated) { /* first time */
667: PetscCall(MatDuplicate(newlA[0], MAT_COPY_VALUES, &lA));
668: PetscCall(MatISSetLocalMat(B, lA));
669: PetscCall(PetscObjectDereference((PetscObject)lA));
670: }
671: PetscCall(MatCopy(newlA[0], lA, SAME_NONZERO_PATTERN));
672: PetscCall(MatDestroySubMatrices(1, &newlA));
673: PetscCall(MatISScaleDisassembling_Private(B));
674: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
675: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
676: if (was_inplace) PetscCall(MatHeaderReplace(A, &B));
677: else *newmat = B;
678: PetscFunctionReturn(PETSC_SUCCESS);
679: }
680: /* rectangular case, just compress out the column space */
681: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij));
682: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIBAIJ, &ismpibaij));
683: if (ismpiaij) {
684: bs = 1;
685: PetscCall(MatMPIAIJGetSeqAIJ(A, &Ad, &Ao, &garray));
686: } else if (ismpibaij) {
687: PetscCall(MatGetBlockSize(A, &bs));
688: PetscCall(MatMPIBAIJGetSeqBAIJ(A, &Ad, &Ao, &garray));
689: PetscCall(MatConvert(Ad, MATSEQAIJ, MAT_INITIAL_MATRIX, &Ad));
690: PetscCall(MatConvert(Ao, MATSEQAIJ, MAT_INITIAL_MATRIX, &Ao));
691: } else SETERRQ(comm, PETSC_ERR_SUP, "Type %s", ((PetscObject)A)->type_name);
692: PetscCall(MatSeqAIJGetArray(Ad, &dd));
693: PetscCall(MatSeqAIJGetArray(Ao, &od));
694: PetscCheck(garray, comm, PETSC_ERR_ARG_WRONGSTATE, "garray not present");
696: /* access relevant information from MPIAIJ */
697: PetscCall(MatGetOwnershipRange(A, &str, NULL));
698: PetscCall(MatGetOwnershipRangeColumn(A, &stc, NULL));
699: PetscCall(MatGetLocalSize(A, &dr, &dc));
700: PetscCall(MatGetLocalSize(Ao, NULL, &oc));
701: PetscCall(MatGetRowIJ(Ad, 0, PETSC_FALSE, PETSC_FALSE, &i, &di, &dj, &flg));
702: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure");
703: PetscCall(MatGetRowIJ(Ao, 0, PETSC_FALSE, PETSC_FALSE, &i, &oi, &oj, &flg));
704: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure");
705: nnz = di[dr] + oi[dr];
706: /* store original pointers to be restored later */
707: odi = di;
708: odj = dj;
709: ooi = oi;
710: ooj = oj;
712: /* generate l2g maps for rows and cols */
713: PetscCall(ISCreateStride(comm, dr / bs, str / bs, 1, &is));
714: if (bs > 1) {
715: IS is2;
717: PetscCall(ISGetLocalSize(is, &i));
718: PetscCall(ISGetIndices(is, (const PetscInt **)&aux));
719: PetscCall(ISCreateBlock(comm, bs, i, aux, PETSC_COPY_VALUES, &is2));
720: PetscCall(ISRestoreIndices(is, (const PetscInt **)&aux));
721: PetscCall(ISDestroy(&is));
722: is = is2;
723: }
724: PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g));
725: PetscCall(ISDestroy(&is));
726: if (dr) {
727: PetscCall(PetscMalloc1((dc + oc) / bs, &aux));
728: for (i = 0; i < dc / bs; i++) aux[i] = i + stc / bs;
729: for (i = 0; i < oc / bs; i++) aux[i + dc / bs] = garray[i];
730: PetscCall(ISCreateBlock(comm, bs, (dc + oc) / bs, aux, PETSC_OWN_POINTER, &is));
731: lc = dc + oc;
732: } else {
733: PetscCall(ISCreateBlock(comm, bs, 0, NULL, PETSC_OWN_POINTER, &is));
734: lc = 0;
735: }
736: PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g));
737: PetscCall(ISDestroy(&is));
739: /* create MATIS object */
740: PetscCall(MatCreate(comm, &B));
741: PetscCall(MatSetSizes(B, dr, dc, PETSC_DECIDE, PETSC_DECIDE));
742: PetscCall(MatSetType(B, MATIS));
743: PetscCall(MatSetBlockSize(B, bs));
744: PetscCall(MatSetLocalToGlobalMapping(B, rl2g, cl2g));
745: PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
746: PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
748: /* merge local matrices */
749: PetscCall(PetscMalloc1(nnz + dr + 1, &aux));
750: PetscCall(PetscMalloc1(nnz, &data));
751: ii = aux;
752: jj = aux + dr + 1;
753: aa = data;
754: *ii = *(di++) + *(oi++);
755: for (jd = 0, jo = 0, cum = 0; *ii < nnz; cum++) {
756: for (; jd < *di; jd++) {
757: *jj++ = *dj++;
758: *aa++ = *dd++;
759: }
760: for (; jo < *oi; jo++) {
761: *jj++ = *oj++ + dc;
762: *aa++ = *od++;
763: }
764: *(++ii) = *(di++) + *(oi++);
765: }
766: for (; cum < dr; cum++) *(++ii) = nnz;
768: PetscCall(MatRestoreRowIJ(Ad, 0, PETSC_FALSE, PETSC_FALSE, &i, &odi, &odj, &flg));
769: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot restore IJ structure");
770: PetscCall(MatRestoreRowIJ(Ao, 0, PETSC_FALSE, PETSC_FALSE, &i, &ooi, &ooj, &flg));
771: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot restore IJ structure");
772: PetscCall(MatSeqAIJRestoreArray(Ad, &dd));
773: PetscCall(MatSeqAIJRestoreArray(Ao, &od));
775: ii = aux;
776: jj = aux + dr + 1;
777: aa = data;
778: PetscCall(MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, dr, lc, ii, jj, aa, &lA));
780: /* create containers to destroy the data */
781: ptrs[0] = aux;
782: ptrs[1] = data;
783: for (i = 0; i < 2; i++) {
784: PetscContainer c;
786: PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &c));
787: PetscCall(PetscContainerSetPointer(c, ptrs[i]));
788: PetscCall(PetscContainerSetUserDestroy(c, PetscContainerUserDestroyDefault));
789: PetscCall(PetscObjectCompose((PetscObject)lA, names[i], (PetscObject)c));
790: PetscCall(PetscContainerDestroy(&c));
791: }
792: if (ismpibaij) { /* destroy converted local matrices */
793: PetscCall(MatDestroy(&Ad));
794: PetscCall(MatDestroy(&Ao));
795: }
797: /* finalize matrix */
798: PetscCall(MatISSetLocalMat(B, lA));
799: PetscCall(MatDestroy(&lA));
800: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
801: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
802: if (reuse == MAT_INPLACE_MATRIX) {
803: PetscCall(MatHeaderReplace(A, &B));
804: } else *newmat = B;
805: PetscFunctionReturn(PETSC_SUCCESS);
806: }
808: PETSC_INTERN PetscErrorCode MatConvert_Nest_IS(Mat A, MatType type, MatReuse reuse, Mat *newmat)
809: {
810: Mat **nest, *snest, **rnest, lA, B;
811: IS *iscol, *isrow, *islrow, *islcol;
812: ISLocalToGlobalMapping rl2g, cl2g;
813: MPI_Comm comm;
814: PetscInt *lr, *lc, *l2gidxs;
815: PetscInt i, j, nr, nc, rbs, cbs;
816: PetscBool convert, lreuse, *istrans;
818: PetscFunctionBegin;
819: PetscCall(MatNestGetSubMats(A, &nr, &nc, &nest));
820: lreuse = PETSC_FALSE;
821: rnest = NULL;
822: if (reuse == MAT_REUSE_MATRIX) {
823: PetscBool ismatis, isnest;
825: PetscCall(PetscObjectTypeCompare((PetscObject)*newmat, MATIS, &ismatis));
826: PetscCheck(ismatis, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_USER, "Cannot reuse matrix of type %s", ((PetscObject)(*newmat))->type_name);
827: PetscCall(MatISGetLocalMat(*newmat, &lA));
828: PetscCall(PetscObjectTypeCompare((PetscObject)lA, MATNEST, &isnest));
829: if (isnest) {
830: PetscCall(MatNestGetSubMats(lA, &i, &j, &rnest));
831: lreuse = (PetscBool)(i == nr && j == nc);
832: if (!lreuse) rnest = NULL;
833: }
834: }
835: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
836: PetscCall(PetscCalloc2(nr, &lr, nc, &lc));
837: PetscCall(PetscCalloc6(nr, &isrow, nc, &iscol, nr, &islrow, nc, &islcol, nr * nc, &snest, nr * nc, &istrans));
838: PetscCall(MatNestGetISs(A, isrow, iscol));
839: for (i = 0; i < nr; i++) {
840: for (j = 0; j < nc; j++) {
841: PetscBool ismatis;
842: PetscInt l1, l2, lb1, lb2, ij = i * nc + j;
844: /* Null matrix pointers are allowed in MATNEST */
845: if (!nest[i][j]) continue;
847: /* Nested matrices should be of type MATIS */
848: PetscCall(PetscObjectTypeCompare((PetscObject)nest[i][j], MATTRANSPOSEVIRTUAL, &istrans[ij]));
849: if (istrans[ij]) {
850: Mat T, lT;
851: PetscCall(MatTransposeGetMat(nest[i][j], &T));
852: PetscCall(PetscObjectTypeCompare((PetscObject)T, MATIS, &ismatis));
853: 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);
854: PetscCall(MatISGetLocalMat(T, &lT));
855: PetscCall(MatCreateTranspose(lT, &snest[ij]));
856: } else {
857: PetscCall(PetscObjectTypeCompare((PetscObject)nest[i][j], MATIS, &ismatis));
858: 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);
859: PetscCall(MatISGetLocalMat(nest[i][j], &snest[ij]));
860: }
862: /* Check compatibility of local sizes */
863: PetscCall(MatGetSize(snest[ij], &l1, &l2));
864: PetscCall(MatGetBlockSizes(snest[ij], &lb1, &lb2));
865: if (!l1 || !l2) continue;
866: 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);
867: 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);
868: lr[i] = l1;
869: lc[j] = l2;
871: /* check compatibility for local matrix reusage */
872: if (rnest && !rnest[i][j] != !snest[ij]) lreuse = PETSC_FALSE;
873: }
874: }
876: if (PetscDefined(USE_DEBUG)) {
877: /* Check compatibility of l2g maps for rows */
878: for (i = 0; i < nr; i++) {
879: rl2g = NULL;
880: for (j = 0; j < nc; j++) {
881: PetscInt n1, n2;
883: if (!nest[i][j]) continue;
884: if (istrans[i * nc + j]) {
885: Mat T;
887: PetscCall(MatTransposeGetMat(nest[i][j], &T));
888: PetscCall(MatISGetLocalToGlobalMapping(T, NULL, &cl2g));
889: } else {
890: PetscCall(MatISGetLocalToGlobalMapping(nest[i][j], &cl2g, NULL));
891: }
892: PetscCall(ISLocalToGlobalMappingGetSize(cl2g, &n1));
893: if (!n1) continue;
894: if (!rl2g) {
895: rl2g = cl2g;
896: } else {
897: const PetscInt *idxs1, *idxs2;
898: PetscBool same;
900: PetscCall(ISLocalToGlobalMappingGetSize(rl2g, &n2));
901: 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);
902: PetscCall(ISLocalToGlobalMappingGetIndices(cl2g, &idxs1));
903: PetscCall(ISLocalToGlobalMappingGetIndices(rl2g, &idxs2));
904: PetscCall(PetscArraycmp(idxs1, idxs2, n1, &same));
905: PetscCall(ISLocalToGlobalMappingRestoreIndices(cl2g, &idxs1));
906: PetscCall(ISLocalToGlobalMappingRestoreIndices(rl2g, &idxs2));
907: 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);
908: }
909: }
910: }
911: /* Check compatibility of l2g maps for columns */
912: for (i = 0; i < nc; i++) {
913: rl2g = NULL;
914: for (j = 0; j < nr; j++) {
915: PetscInt n1, n2;
917: if (!nest[j][i]) continue;
918: if (istrans[j * nc + i]) {
919: Mat T;
921: PetscCall(MatTransposeGetMat(nest[j][i], &T));
922: PetscCall(MatISGetLocalToGlobalMapping(T, &cl2g, NULL));
923: } else {
924: PetscCall(MatISGetLocalToGlobalMapping(nest[j][i], NULL, &cl2g));
925: }
926: PetscCall(ISLocalToGlobalMappingGetSize(cl2g, &n1));
927: if (!n1) continue;
928: if (!rl2g) {
929: rl2g = cl2g;
930: } else {
931: const PetscInt *idxs1, *idxs2;
932: PetscBool same;
934: PetscCall(ISLocalToGlobalMappingGetSize(rl2g, &n2));
935: 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);
936: PetscCall(ISLocalToGlobalMappingGetIndices(cl2g, &idxs1));
937: PetscCall(ISLocalToGlobalMappingGetIndices(rl2g, &idxs2));
938: PetscCall(PetscArraycmp(idxs1, idxs2, n1, &same));
939: PetscCall(ISLocalToGlobalMappingRestoreIndices(cl2g, &idxs1));
940: PetscCall(ISLocalToGlobalMappingRestoreIndices(rl2g, &idxs2));
941: 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);
942: }
943: }
944: }
945: }
947: B = NULL;
948: if (reuse != MAT_REUSE_MATRIX) {
949: PetscInt stl;
951: /* Create l2g map for the rows of the new matrix and index sets for the local MATNEST */
952: for (i = 0, stl = 0; i < nr; i++) stl += lr[i];
953: PetscCall(PetscMalloc1(stl, &l2gidxs));
954: for (i = 0, stl = 0; i < nr; i++) {
955: Mat usedmat;
956: Mat_IS *matis;
957: const PetscInt *idxs;
959: /* local IS for local NEST */
960: PetscCall(ISCreateStride(PETSC_COMM_SELF, lr[i], stl, 1, &islrow[i]));
962: /* l2gmap */
963: j = 0;
964: usedmat = nest[i][j];
965: while (!usedmat && j < nc - 1) usedmat = nest[i][++j];
966: PetscCheck(usedmat, comm, PETSC_ERR_SUP, "Cannot find valid row mat");
968: if (istrans[i * nc + j]) {
969: Mat T;
970: PetscCall(MatTransposeGetMat(usedmat, &T));
971: usedmat = T;
972: }
973: matis = (Mat_IS *)(usedmat->data);
974: PetscCall(ISGetIndices(isrow[i], &idxs));
975: if (istrans[i * nc + j]) {
976: PetscCall(PetscSFBcastBegin(matis->csf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
977: PetscCall(PetscSFBcastEnd(matis->csf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
978: } else {
979: PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
980: PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
981: }
982: PetscCall(ISRestoreIndices(isrow[i], &idxs));
983: stl += lr[i];
984: }
985: PetscCall(ISLocalToGlobalMappingCreate(comm, 1, stl, l2gidxs, PETSC_OWN_POINTER, &rl2g));
987: /* Create l2g map for columns of the new matrix and index sets for the local MATNEST */
988: for (i = 0, stl = 0; i < nc; i++) stl += lc[i];
989: PetscCall(PetscMalloc1(stl, &l2gidxs));
990: for (i = 0, stl = 0; i < nc; i++) {
991: Mat usedmat;
992: Mat_IS *matis;
993: const PetscInt *idxs;
995: /* local IS for local NEST */
996: PetscCall(ISCreateStride(PETSC_COMM_SELF, lc[i], stl, 1, &islcol[i]));
998: /* l2gmap */
999: j = 0;
1000: usedmat = nest[j][i];
1001: while (!usedmat && j < nr - 1) usedmat = nest[++j][i];
1002: PetscCheck(usedmat, comm, PETSC_ERR_SUP, "Cannot find valid column mat");
1003: if (istrans[j * nc + i]) {
1004: Mat T;
1005: PetscCall(MatTransposeGetMat(usedmat, &T));
1006: usedmat = T;
1007: }
1008: matis = (Mat_IS *)(usedmat->data);
1009: PetscCall(ISGetIndices(iscol[i], &idxs));
1010: if (istrans[j * nc + i]) {
1011: PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
1012: PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
1013: } else {
1014: PetscCall(PetscSFBcastBegin(matis->csf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
1015: PetscCall(PetscSFBcastEnd(matis->csf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
1016: }
1017: PetscCall(ISRestoreIndices(iscol[i], &idxs));
1018: stl += lc[i];
1019: }
1020: PetscCall(ISLocalToGlobalMappingCreate(comm, 1, stl, l2gidxs, PETSC_OWN_POINTER, &cl2g));
1022: /* Create MATIS */
1023: PetscCall(MatCreate(comm, &B));
1024: PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
1025: PetscCall(MatGetBlockSizes(A, &rbs, &cbs));
1026: PetscCall(MatSetBlockSizes(B, rbs, cbs));
1027: PetscCall(MatSetType(B, MATIS));
1028: PetscCall(MatISSetLocalMatType(B, MATNEST));
1029: { /* hack : avoid setup of scatters */
1030: Mat_IS *matis = (Mat_IS *)(B->data);
1031: matis->islocalref = PETSC_TRUE;
1032: }
1033: PetscCall(MatSetLocalToGlobalMapping(B, rl2g, cl2g));
1034: PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
1035: PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
1036: PetscCall(MatCreateNest(PETSC_COMM_SELF, nr, islrow, nc, islcol, snest, &lA));
1037: PetscCall(MatNestSetVecType(lA, VECNEST));
1038: for (i = 0; i < nr * nc; i++) {
1039: if (istrans[i]) PetscCall(MatDestroy(&snest[i]));
1040: }
1041: PetscCall(MatISSetLocalMat(B, lA));
1042: PetscCall(MatDestroy(&lA));
1043: { /* hack : setup of scatters done here */
1044: Mat_IS *matis = (Mat_IS *)(B->data);
1046: matis->islocalref = PETSC_FALSE;
1047: PetscCall(MatISSetUpScatters_Private(B));
1048: }
1049: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
1050: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
1051: if (reuse == MAT_INPLACE_MATRIX) {
1052: PetscCall(MatHeaderReplace(A, &B));
1053: } else {
1054: *newmat = B;
1055: }
1056: } else {
1057: if (lreuse) {
1058: PetscCall(MatISGetLocalMat(*newmat, &lA));
1059: for (i = 0; i < nr; i++) {
1060: for (j = 0; j < nc; j++) {
1061: if (snest[i * nc + j]) {
1062: PetscCall(MatNestSetSubMat(lA, i, j, snest[i * nc + j]));
1063: if (istrans[i * nc + j]) PetscCall(MatDestroy(&snest[i * nc + j]));
1064: }
1065: }
1066: }
1067: } else {
1068: PetscInt stl;
1069: for (i = 0, stl = 0; i < nr; i++) {
1070: PetscCall(ISCreateStride(PETSC_COMM_SELF, lr[i], stl, 1, &islrow[i]));
1071: stl += lr[i];
1072: }
1073: for (i = 0, stl = 0; i < nc; i++) {
1074: PetscCall(ISCreateStride(PETSC_COMM_SELF, lc[i], stl, 1, &islcol[i]));
1075: stl += lc[i];
1076: }
1077: PetscCall(MatCreateNest(PETSC_COMM_SELF, nr, islrow, nc, islcol, snest, &lA));
1078: for (i = 0; i < nr * nc; i++) {
1079: if (istrans[i]) PetscCall(MatDestroy(&snest[i]));
1080: }
1081: PetscCall(MatISSetLocalMat(*newmat, lA));
1082: PetscCall(MatDestroy(&lA));
1083: }
1084: PetscCall(MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY));
1085: PetscCall(MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY));
1086: }
1088: /* Create local matrix in MATNEST format */
1089: convert = PETSC_FALSE;
1090: PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)A)->prefix, "-matis_convert_local_nest", &convert, NULL));
1091: if (convert) {
1092: Mat M;
1093: MatISLocalFields lf;
1094: PetscContainer c;
1096: PetscCall(MatISGetLocalMat(*newmat, &lA));
1097: PetscCall(MatConvert(lA, MATAIJ, MAT_INITIAL_MATRIX, &M));
1098: PetscCall(MatISSetLocalMat(*newmat, M));
1099: PetscCall(MatDestroy(&M));
1101: /* attach local fields to the matrix */
1102: PetscCall(PetscNew(&lf));
1103: PetscCall(PetscMalloc2(nr, &lf->rf, nc, &lf->cf));
1104: for (i = 0; i < nr; i++) {
1105: PetscInt n, st;
1107: PetscCall(ISGetLocalSize(islrow[i], &n));
1108: PetscCall(ISStrideGetInfo(islrow[i], &st, NULL));
1109: PetscCall(ISCreateStride(comm, n, st, 1, &lf->rf[i]));
1110: }
1111: for (i = 0; i < nc; i++) {
1112: PetscInt n, st;
1114: PetscCall(ISGetLocalSize(islcol[i], &n));
1115: PetscCall(ISStrideGetInfo(islcol[i], &st, NULL));
1116: PetscCall(ISCreateStride(comm, n, st, 1, &lf->cf[i]));
1117: }
1118: lf->nr = nr;
1119: lf->nc = nc;
1120: PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)(*newmat)), &c));
1121: PetscCall(PetscContainerSetPointer(c, lf));
1122: PetscCall(PetscContainerSetUserDestroy(c, MatISContainerDestroyFields_Private));
1123: PetscCall(PetscObjectCompose((PetscObject)(*newmat), "_convert_nest_lfields", (PetscObject)c));
1124: PetscCall(PetscContainerDestroy(&c));
1125: }
1127: /* Free workspace */
1128: for (i = 0; i < nr; i++) PetscCall(ISDestroy(&islrow[i]));
1129: for (i = 0; i < nc; i++) PetscCall(ISDestroy(&islcol[i]));
1130: PetscCall(PetscFree6(isrow, iscol, islrow, islcol, snest, istrans));
1131: PetscCall(PetscFree2(lr, lc));
1132: PetscFunctionReturn(PETSC_SUCCESS);
1133: }
1135: static PetscErrorCode MatDiagonalScale_IS(Mat A, Vec l, Vec r)
1136: {
1137: Mat_IS *matis = (Mat_IS *)A->data;
1138: Vec ll, rr;
1139: const PetscScalar *Y, *X;
1140: PetscScalar *x, *y;
1142: PetscFunctionBegin;
1143: if (l) {
1144: ll = matis->y;
1145: PetscCall(VecGetArrayRead(l, &Y));
1146: PetscCall(VecGetArray(ll, &y));
1147: PetscCall(PetscSFBcastBegin(matis->sf, MPIU_SCALAR, Y, y, MPI_REPLACE));
1148: } else {
1149: ll = NULL;
1150: }
1151: if (r) {
1152: rr = matis->x;
1153: PetscCall(VecGetArrayRead(r, &X));
1154: PetscCall(VecGetArray(rr, &x));
1155: PetscCall(PetscSFBcastBegin(matis->csf, MPIU_SCALAR, X, x, MPI_REPLACE));
1156: } else {
1157: rr = NULL;
1158: }
1159: if (ll) {
1160: PetscCall(PetscSFBcastEnd(matis->sf, MPIU_SCALAR, Y, y, MPI_REPLACE));
1161: PetscCall(VecRestoreArrayRead(l, &Y));
1162: PetscCall(VecRestoreArray(ll, &y));
1163: }
1164: if (rr) {
1165: PetscCall(PetscSFBcastEnd(matis->csf, MPIU_SCALAR, X, x, MPI_REPLACE));
1166: PetscCall(VecRestoreArrayRead(r, &X));
1167: PetscCall(VecRestoreArray(rr, &x));
1168: }
1169: PetscCall(MatDiagonalScale(matis->A, ll, rr));
1170: PetscFunctionReturn(PETSC_SUCCESS);
1171: }
1173: static PetscErrorCode MatGetInfo_IS(Mat A, MatInfoType flag, MatInfo *ginfo)
1174: {
1175: Mat_IS *matis = (Mat_IS *)A->data;
1176: MatInfo info;
1177: PetscLogDouble isend[6], irecv[6];
1178: PetscInt bs;
1180: PetscFunctionBegin;
1181: PetscCall(MatGetBlockSize(A, &bs));
1182: if (matis->A->ops->getinfo) {
1183: PetscCall(MatGetInfo(matis->A, MAT_LOCAL, &info));
1184: isend[0] = info.nz_used;
1185: isend[1] = info.nz_allocated;
1186: isend[2] = info.nz_unneeded;
1187: isend[3] = info.memory;
1188: isend[4] = info.mallocs;
1189: } else {
1190: isend[0] = 0.;
1191: isend[1] = 0.;
1192: isend[2] = 0.;
1193: isend[3] = 0.;
1194: isend[4] = 0.;
1195: }
1196: isend[5] = matis->A->num_ass;
1197: if (flag == MAT_LOCAL) {
1198: ginfo->nz_used = isend[0];
1199: ginfo->nz_allocated = isend[1];
1200: ginfo->nz_unneeded = isend[2];
1201: ginfo->memory = isend[3];
1202: ginfo->mallocs = isend[4];
1203: ginfo->assemblies = isend[5];
1204: } else if (flag == MAT_GLOBAL_MAX) {
1205: PetscCall(MPIU_Allreduce(isend, irecv, 6, MPIU_PETSCLOGDOUBLE, MPI_MAX, PetscObjectComm((PetscObject)A)));
1207: ginfo->nz_used = irecv[0];
1208: ginfo->nz_allocated = irecv[1];
1209: ginfo->nz_unneeded = irecv[2];
1210: ginfo->memory = irecv[3];
1211: ginfo->mallocs = irecv[4];
1212: ginfo->assemblies = irecv[5];
1213: } else if (flag == MAT_GLOBAL_SUM) {
1214: PetscCall(MPIU_Allreduce(isend, irecv, 5, MPIU_PETSCLOGDOUBLE, MPI_SUM, PetscObjectComm((PetscObject)A)));
1216: ginfo->nz_used = irecv[0];
1217: ginfo->nz_allocated = irecv[1];
1218: ginfo->nz_unneeded = irecv[2];
1219: ginfo->memory = irecv[3];
1220: ginfo->mallocs = irecv[4];
1221: ginfo->assemblies = A->num_ass;
1222: }
1223: ginfo->block_size = bs;
1224: ginfo->fill_ratio_given = 0;
1225: ginfo->fill_ratio_needed = 0;
1226: ginfo->factor_mallocs = 0;
1227: PetscFunctionReturn(PETSC_SUCCESS);
1228: }
1230: static PetscErrorCode MatTranspose_IS(Mat A, MatReuse reuse, Mat *B)
1231: {
1232: Mat C, lC, lA;
1234: PetscFunctionBegin;
1235: if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *B));
1236: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
1237: ISLocalToGlobalMapping rl2g, cl2g;
1238: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
1239: PetscCall(MatSetSizes(C, A->cmap->n, A->rmap->n, A->cmap->N, A->rmap->N));
1240: PetscCall(MatSetBlockSizes(C, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs)));
1241: PetscCall(MatSetType(C, MATIS));
1242: PetscCall(MatGetLocalToGlobalMapping(A, &rl2g, &cl2g));
1243: PetscCall(MatSetLocalToGlobalMapping(C, cl2g, rl2g));
1244: } else C = *B;
1246: /* perform local transposition */
1247: PetscCall(MatISGetLocalMat(A, &lA));
1248: PetscCall(MatTranspose(lA, MAT_INITIAL_MATRIX, &lC));
1249: PetscCall(MatSetLocalToGlobalMapping(lC, lA->cmap->mapping, lA->rmap->mapping));
1250: PetscCall(MatISSetLocalMat(C, lC));
1251: PetscCall(MatDestroy(&lC));
1253: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
1254: *B = C;
1255: } else {
1256: PetscCall(MatHeaderMerge(A, &C));
1257: }
1258: PetscCall(MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY));
1259: PetscCall(MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY));
1260: PetscFunctionReturn(PETSC_SUCCESS);
1261: }
1263: static PetscErrorCode MatDiagonalSet_IS(Mat A, Vec D, InsertMode insmode)
1264: {
1265: Mat_IS *is = (Mat_IS *)A->data;
1267: PetscFunctionBegin;
1268: if (D) { /* MatShift_IS pass D = NULL */
1269: PetscCall(VecScatterBegin(is->rctx, D, is->y, INSERT_VALUES, SCATTER_FORWARD));
1270: PetscCall(VecScatterEnd(is->rctx, D, is->y, INSERT_VALUES, SCATTER_FORWARD));
1271: }
1272: PetscCall(VecPointwiseDivide(is->y, is->y, is->counter));
1273: PetscCall(MatDiagonalSet(is->A, is->y, insmode));
1274: PetscFunctionReturn(PETSC_SUCCESS);
1275: }
1277: static PetscErrorCode MatShift_IS(Mat A, PetscScalar a)
1278: {
1279: Mat_IS *is = (Mat_IS *)A->data;
1281: PetscFunctionBegin;
1282: PetscCall(VecSet(is->y, a));
1283: PetscCall(MatDiagonalSet_IS(A, NULL, ADD_VALUES));
1284: PetscFunctionReturn(PETSC_SUCCESS);
1285: }
1287: static PetscErrorCode MatSetValuesLocal_SubMat_IS(Mat A, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
1288: {
1289: PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION], cols_l[MATIS_MAX_ENTRIES_INSERTION];
1291: PetscFunctionBegin;
1292: PetscCheck(m <= MATIS_MAX_ENTRIES_INSERTION && n <= MATIS_MAX_ENTRIES_INSERTION, PETSC_COMM_SELF, PETSC_ERR_SUP, "Number of row/column indices must be <= %d: they are %" PetscInt_FMT " %" PetscInt_FMT, MATIS_MAX_ENTRIES_INSERTION, m, n);
1293: PetscCall(ISLocalToGlobalMappingApply(A->rmap->mapping, m, rows, rows_l));
1294: PetscCall(ISLocalToGlobalMappingApply(A->cmap->mapping, n, cols, cols_l));
1295: PetscCall(MatSetValuesLocal_IS(A, m, rows_l, n, cols_l, values, addv));
1296: PetscFunctionReturn(PETSC_SUCCESS);
1297: }
1299: static PetscErrorCode MatSetValuesBlockedLocal_SubMat_IS(Mat A, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
1300: {
1301: PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION], cols_l[MATIS_MAX_ENTRIES_INSERTION];
1303: PetscFunctionBegin;
1304: PetscCheck(m <= MATIS_MAX_ENTRIES_INSERTION && n <= MATIS_MAX_ENTRIES_INSERTION, PETSC_COMM_SELF, PETSC_ERR_SUP, "Number of row/column block indices must be <= %d: they are %" PetscInt_FMT " %" PetscInt_FMT, MATIS_MAX_ENTRIES_INSERTION, m, n);
1305: PetscCall(ISLocalToGlobalMappingApplyBlock(A->rmap->mapping, m, rows, rows_l));
1306: PetscCall(ISLocalToGlobalMappingApplyBlock(A->cmap->mapping, n, cols, cols_l));
1307: PetscCall(MatSetValuesBlockedLocal_IS(A, m, rows_l, n, cols_l, values, addv));
1308: PetscFunctionReturn(PETSC_SUCCESS);
1309: }
1311: static PetscErrorCode MatCreateSubMatrix_IS(Mat mat, IS irow, IS icol, MatReuse scall, Mat *newmat)
1312: {
1313: Mat locmat, newlocmat;
1314: Mat_IS *newmatis;
1315: const PetscInt *idxs;
1316: PetscInt i, m, n;
1318: PetscFunctionBegin;
1319: if (scall == MAT_REUSE_MATRIX) {
1320: PetscBool ismatis;
1322: PetscCall(PetscObjectTypeCompare((PetscObject)*newmat, MATIS, &ismatis));
1323: PetscCheck(ismatis, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_ARG_WRONG, "Cannot reuse matrix! Not of MATIS type");
1324: newmatis = (Mat_IS *)(*newmat)->data;
1325: PetscCheck(newmatis->getsub_ris, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_ARG_WRONG, "Cannot reuse matrix! Misses local row IS");
1326: PetscCheck(newmatis->getsub_cis, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_ARG_WRONG, "Cannot reuse matrix! Misses local col IS");
1327: }
1328: /* irow and icol may not have duplicate entries */
1329: if (PetscDefined(USE_DEBUG)) {
1330: Vec rtest, ltest;
1331: const PetscScalar *array;
1333: PetscCall(MatCreateVecs(mat, <est, &rtest));
1334: PetscCall(ISGetLocalSize(irow, &n));
1335: PetscCall(ISGetIndices(irow, &idxs));
1336: for (i = 0; i < n; i++) PetscCall(VecSetValue(rtest, idxs[i], 1.0, ADD_VALUES));
1337: PetscCall(VecAssemblyBegin(rtest));
1338: PetscCall(VecAssemblyEnd(rtest));
1339: PetscCall(VecGetLocalSize(rtest, &n));
1340: PetscCall(VecGetOwnershipRange(rtest, &m, NULL));
1341: PetscCall(VecGetArrayRead(rtest, &array));
1342: 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]));
1343: PetscCall(VecRestoreArrayRead(rtest, &array));
1344: PetscCall(ISRestoreIndices(irow, &idxs));
1345: PetscCall(ISGetLocalSize(icol, &n));
1346: PetscCall(ISGetIndices(icol, &idxs));
1347: for (i = 0; i < n; i++) PetscCall(VecSetValue(ltest, idxs[i], 1.0, ADD_VALUES));
1348: PetscCall(VecAssemblyBegin(ltest));
1349: PetscCall(VecAssemblyEnd(ltest));
1350: PetscCall(VecGetLocalSize(ltest, &n));
1351: PetscCall(VecGetOwnershipRange(ltest, &m, NULL));
1352: PetscCall(VecGetArrayRead(ltest, &array));
1353: 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]));
1354: PetscCall(VecRestoreArrayRead(ltest, &array));
1355: PetscCall(ISRestoreIndices(icol, &idxs));
1356: PetscCall(VecDestroy(&rtest));
1357: PetscCall(VecDestroy(<est));
1358: }
1359: if (scall == MAT_INITIAL_MATRIX) {
1360: Mat_IS *matis = (Mat_IS *)mat->data;
1361: ISLocalToGlobalMapping rl2g;
1362: IS is;
1363: PetscInt *lidxs, *lgidxs, *newgidxs;
1364: PetscInt ll, newloc, irbs, icbs, arbs, acbs, rbs, cbs;
1365: PetscBool cong;
1366: MPI_Comm comm;
1368: PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
1369: PetscCall(MatGetBlockSizes(mat, &arbs, &acbs));
1370: PetscCall(ISGetBlockSize(irow, &irbs));
1371: PetscCall(ISGetBlockSize(icol, &icbs));
1372: rbs = arbs == irbs ? irbs : 1;
1373: cbs = acbs == icbs ? icbs : 1;
1374: PetscCall(ISGetLocalSize(irow, &m));
1375: PetscCall(ISGetLocalSize(icol, &n));
1376: PetscCall(MatCreate(comm, newmat));
1377: PetscCall(MatSetType(*newmat, MATIS));
1378: PetscCall(MatSetSizes(*newmat, m, n, PETSC_DECIDE, PETSC_DECIDE));
1379: PetscCall(MatSetBlockSizes(*newmat, rbs, cbs));
1380: /* communicate irow to their owners in the layout */
1381: PetscCall(ISGetIndices(irow, &idxs));
1382: PetscCall(PetscLayoutMapLocal(mat->rmap, m, idxs, &ll, &lidxs, &lgidxs));
1383: PetscCall(ISRestoreIndices(irow, &idxs));
1384: PetscCall(PetscArrayzero(matis->sf_rootdata, matis->sf->nroots));
1385: for (i = 0; i < ll; i++) matis->sf_rootdata[lidxs[i]] = lgidxs[i] + 1;
1386: PetscCall(PetscFree(lidxs));
1387: PetscCall(PetscFree(lgidxs));
1388: PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
1389: PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
1390: for (i = 0, newloc = 0; i < matis->sf->nleaves; i++)
1391: if (matis->sf_leafdata[i]) newloc++;
1392: PetscCall(PetscMalloc1(newloc, &newgidxs));
1393: PetscCall(PetscMalloc1(newloc, &lidxs));
1394: for (i = 0, newloc = 0; i < matis->sf->nleaves; i++)
1395: if (matis->sf_leafdata[i]) {
1396: lidxs[newloc] = i;
1397: newgidxs[newloc++] = matis->sf_leafdata[i] - 1;
1398: }
1399: PetscCall(ISCreateGeneral(comm, newloc, newgidxs, PETSC_OWN_POINTER, &is));
1400: PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g));
1401: PetscCall(ISLocalToGlobalMappingSetBlockSize(rl2g, rbs));
1402: PetscCall(ISDestroy(&is));
1403: /* local is to extract local submatrix */
1404: newmatis = (Mat_IS *)(*newmat)->data;
1405: PetscCall(ISCreateGeneral(comm, newloc, lidxs, PETSC_OWN_POINTER, &newmatis->getsub_ris));
1406: PetscCall(MatHasCongruentLayouts(mat, &cong));
1407: if (cong && irow == icol && matis->csf == matis->sf) {
1408: PetscCall(MatSetLocalToGlobalMapping(*newmat, rl2g, rl2g));
1409: PetscCall(PetscObjectReference((PetscObject)newmatis->getsub_ris));
1410: newmatis->getsub_cis = newmatis->getsub_ris;
1411: } else {
1412: ISLocalToGlobalMapping cl2g;
1414: /* communicate icol to their owners in the layout */
1415: PetscCall(ISGetIndices(icol, &idxs));
1416: PetscCall(PetscLayoutMapLocal(mat->cmap, n, idxs, &ll, &lidxs, &lgidxs));
1417: PetscCall(ISRestoreIndices(icol, &idxs));
1418: PetscCall(PetscArrayzero(matis->csf_rootdata, matis->csf->nroots));
1419: for (i = 0; i < ll; i++) matis->csf_rootdata[lidxs[i]] = lgidxs[i] + 1;
1420: PetscCall(PetscFree(lidxs));
1421: PetscCall(PetscFree(lgidxs));
1422: PetscCall(PetscSFBcastBegin(matis->csf, MPIU_INT, matis->csf_rootdata, matis->csf_leafdata, MPI_REPLACE));
1423: PetscCall(PetscSFBcastEnd(matis->csf, MPIU_INT, matis->csf_rootdata, matis->csf_leafdata, MPI_REPLACE));
1424: for (i = 0, newloc = 0; i < matis->csf->nleaves; i++)
1425: if (matis->csf_leafdata[i]) newloc++;
1426: PetscCall(PetscMalloc1(newloc, &newgidxs));
1427: PetscCall(PetscMalloc1(newloc, &lidxs));
1428: for (i = 0, newloc = 0; i < matis->csf->nleaves; i++)
1429: if (matis->csf_leafdata[i]) {
1430: lidxs[newloc] = i;
1431: newgidxs[newloc++] = matis->csf_leafdata[i] - 1;
1432: }
1433: PetscCall(ISCreateGeneral(comm, newloc, newgidxs, PETSC_OWN_POINTER, &is));
1434: PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g));
1435: PetscCall(ISLocalToGlobalMappingSetBlockSize(cl2g, cbs));
1436: PetscCall(ISDestroy(&is));
1437: /* local is to extract local submatrix */
1438: PetscCall(ISCreateGeneral(comm, newloc, lidxs, PETSC_OWN_POINTER, &newmatis->getsub_cis));
1439: PetscCall(MatSetLocalToGlobalMapping(*newmat, rl2g, cl2g));
1440: PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
1441: }
1442: PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
1443: } else {
1444: PetscCall(MatISGetLocalMat(*newmat, &newlocmat));
1445: }
1446: PetscCall(MatISGetLocalMat(mat, &locmat));
1447: newmatis = (Mat_IS *)(*newmat)->data;
1448: PetscCall(MatCreateSubMatrix(locmat, newmatis->getsub_ris, newmatis->getsub_cis, scall, &newlocmat));
1449: if (scall == MAT_INITIAL_MATRIX) {
1450: PetscCall(MatISSetLocalMat(*newmat, newlocmat));
1451: PetscCall(MatDestroy(&newlocmat));
1452: }
1453: PetscCall(MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY));
1454: PetscCall(MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY));
1455: PetscFunctionReturn(PETSC_SUCCESS);
1456: }
1458: static PetscErrorCode MatCopy_IS(Mat A, Mat B, MatStructure str)
1459: {
1460: Mat_IS *a = (Mat_IS *)A->data, *b;
1461: PetscBool ismatis;
1463: PetscFunctionBegin;
1464: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATIS, &ismatis));
1465: PetscCheck(ismatis, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Need to be implemented");
1466: b = (Mat_IS *)B->data;
1467: PetscCall(MatCopy(a->A, b->A, str));
1468: PetscCall(PetscObjectStateIncrease((PetscObject)B));
1469: PetscFunctionReturn(PETSC_SUCCESS);
1470: }
1472: static PetscErrorCode MatMissingDiagonal_IS(Mat A, PetscBool *missing, PetscInt *d)
1473: {
1474: Vec v;
1475: const PetscScalar *array;
1476: PetscInt i, n;
1478: PetscFunctionBegin;
1479: *missing = PETSC_FALSE;
1480: PetscCall(MatCreateVecs(A, NULL, &v));
1481: PetscCall(MatGetDiagonal(A, v));
1482: PetscCall(VecGetLocalSize(v, &n));
1483: PetscCall(VecGetArrayRead(v, &array));
1484: for (i = 0; i < n; i++)
1485: if (array[i] == 0.) break;
1486: PetscCall(VecRestoreArrayRead(v, &array));
1487: PetscCall(VecDestroy(&v));
1488: if (i != n) *missing = PETSC_TRUE;
1489: if (d) {
1490: *d = -1;
1491: if (*missing) {
1492: PetscInt rstart;
1493: PetscCall(MatGetOwnershipRange(A, &rstart, NULL));
1494: *d = i + rstart;
1495: }
1496: }
1497: PetscFunctionReturn(PETSC_SUCCESS);
1498: }
1500: static PetscErrorCode MatISSetUpSF_IS(Mat B)
1501: {
1502: Mat_IS *matis = (Mat_IS *)(B->data);
1503: const PetscInt *gidxs;
1504: PetscInt nleaves;
1506: PetscFunctionBegin;
1507: if (matis->sf) PetscFunctionReturn(PETSC_SUCCESS);
1508: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)B), &matis->sf));
1509: PetscCall(ISLocalToGlobalMappingGetIndices(matis->rmapping, &gidxs));
1510: PetscCall(ISLocalToGlobalMappingGetSize(matis->rmapping, &nleaves));
1511: PetscCall(PetscSFSetGraphLayout(matis->sf, B->rmap, nleaves, NULL, PETSC_OWN_POINTER, gidxs));
1512: PetscCall(ISLocalToGlobalMappingRestoreIndices(matis->rmapping, &gidxs));
1513: PetscCall(PetscMalloc2(matis->sf->nroots, &matis->sf_rootdata, matis->sf->nleaves, &matis->sf_leafdata));
1514: if (matis->rmapping != matis->cmapping) { /* setup SF for columns */
1515: PetscCall(ISLocalToGlobalMappingGetSize(matis->cmapping, &nleaves));
1516: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)B), &matis->csf));
1517: PetscCall(ISLocalToGlobalMappingGetIndices(matis->cmapping, &gidxs));
1518: PetscCall(PetscSFSetGraphLayout(matis->csf, B->cmap, nleaves, NULL, PETSC_OWN_POINTER, gidxs));
1519: PetscCall(ISLocalToGlobalMappingRestoreIndices(matis->cmapping, &gidxs));
1520: PetscCall(PetscMalloc2(matis->csf->nroots, &matis->csf_rootdata, matis->csf->nleaves, &matis->csf_leafdata));
1521: } else {
1522: matis->csf = matis->sf;
1523: matis->csf_leafdata = matis->sf_leafdata;
1524: matis->csf_rootdata = matis->sf_rootdata;
1525: }
1526: PetscFunctionReturn(PETSC_SUCCESS);
1527: }
1529: /*@
1530: MatISStoreL2L - Store local-to-local operators during the Galerkin process of computing `MatPtAP()`
1532: Collective
1534: Input Parameters:
1535: + A - the matrix
1536: - store - the boolean flag
1538: Level: advanced
1540: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateIS()`, `MatISSetPreallocation()`, `MatPtAP()`
1541: @*/
1542: PetscErrorCode MatISStoreL2L(Mat A, PetscBool store)
1543: {
1544: PetscFunctionBegin;
1548: PetscTryMethod(A, "MatISStoreL2L_C", (Mat, PetscBool), (A, store));
1549: PetscFunctionReturn(PETSC_SUCCESS);
1550: }
1552: static PetscErrorCode MatISStoreL2L_IS(Mat A, PetscBool store)
1553: {
1554: Mat_IS *matis = (Mat_IS *)(A->data);
1556: PetscFunctionBegin;
1557: matis->storel2l = store;
1558: if (!store) PetscCall(PetscObjectCompose((PetscObject)(A), "_MatIS_PtAP_l2l", NULL));
1559: PetscFunctionReturn(PETSC_SUCCESS);
1560: }
1562: /*@
1563: MatISFixLocalEmpty - Compress out zero local rows from the local matrices
1565: Collective
1567: Input Parameters:
1568: + A - the matrix
1569: - fix - the boolean flag
1571: Level: advanced
1573: Note:
1574: When `fix` is `PETSC_TRUE`, new local matrices and l2g maps are generated during the final assembly process.
1576: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatCreate()`, `MatCreateIS()`, `MatISSetPreallocation()`, `MatAssemblyEnd()`, `MAT_FINAL_ASSEMBLY`
1577: @*/
1578: PetscErrorCode MatISFixLocalEmpty(Mat A, PetscBool fix)
1579: {
1580: PetscFunctionBegin;
1584: PetscTryMethod(A, "MatISFixLocalEmpty_C", (Mat, PetscBool), (A, fix));
1585: PetscFunctionReturn(PETSC_SUCCESS);
1586: }
1588: static PetscErrorCode MatISFixLocalEmpty_IS(Mat A, PetscBool fix)
1589: {
1590: Mat_IS *matis = (Mat_IS *)(A->data);
1592: PetscFunctionBegin;
1593: matis->locempty = fix;
1594: PetscFunctionReturn(PETSC_SUCCESS);
1595: }
1597: /*@
1598: MatISSetPreallocation - Preallocates memory for a `MATIS` parallel matrix.
1600: Collective
1602: Input Parameters:
1603: + B - the matrix
1604: . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix
1605: (same value is used for all local rows)
1606: . d_nnz - array containing the number of nonzeros in the various rows of the
1607: DIAGONAL portion of the local submatrix (possibly different for each row)
1608: or `NULL`, if `d_nz` is used to specify the nonzero structure.
1609: The size of this array is equal to the number of local rows, i.e `m`.
1610: For matrices that will be factored, you must leave room for (and set)
1611: the diagonal entry even if it is zero.
1612: . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local
1613: submatrix (same value is used for all local rows).
1614: - o_nnz - array containing the number of nonzeros in the various rows of the
1615: OFF-DIAGONAL portion of the local submatrix (possibly different for
1616: each row) or `NULL`, if `o_nz` is used to specify the nonzero
1617: structure. The size of this array is equal to the number
1618: of local rows, i.e `m`.
1620: If the *_nnz parameter is given then the *_nz parameter is ignored
1622: Level: intermediate
1624: Note:
1625: This function has the same interface as the `MATMPIAIJ` preallocation routine in order to simplify the transition
1626: from the asssembled format to the unassembled one. It overestimates the preallocation of `MATIS` local
1627: matrices; for exact preallocation, the user should set the preallocation directly on local matrix objects.
1629: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateIS()`, `MatMPIAIJSetPreallocation()`, `MatISGetLocalMat()`, `MATIS`
1630: @*/
1631: PetscErrorCode MatISSetPreallocation(Mat B, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[])
1632: {
1633: PetscFunctionBegin;
1636: PetscTryMethod(B, "MatISSetPreallocation_C", (Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[]), (B, d_nz, d_nnz, o_nz, o_nnz));
1637: PetscFunctionReturn(PETSC_SUCCESS);
1638: }
1640: /* this is used by DMDA */
1641: PETSC_EXTERN PetscErrorCode MatISSetPreallocation_IS(Mat B, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[])
1642: {
1643: Mat_IS *matis = (Mat_IS *)(B->data);
1644: PetscInt bs, i, nlocalcols;
1646: PetscFunctionBegin;
1647: PetscCall(MatSetUp(B));
1648: if (!d_nnz)
1649: for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] = d_nz;
1650: else
1651: for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] = d_nnz[i];
1653: if (!o_nnz)
1654: for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] += o_nz;
1655: else
1656: for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] += o_nnz[i];
1658: PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
1659: PetscCall(MatGetSize(matis->A, NULL, &nlocalcols));
1660: PetscCall(MatGetBlockSize(matis->A, &bs));
1661: PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
1663: for (i = 0; i < matis->sf->nleaves; i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i], nlocalcols);
1664: PetscCall(MatSeqAIJSetPreallocation(matis->A, 0, matis->sf_leafdata));
1665: #if defined(PETSC_HAVE_HYPRE)
1666: PetscCall(MatHYPRESetPreallocation(matis->A, 0, matis->sf_leafdata, 0, NULL));
1667: #endif
1669: for (i = 0; i < matis->sf->nleaves / bs; i++) {
1670: PetscInt b;
1672: matis->sf_leafdata[i] = matis->sf_leafdata[i * bs] / bs;
1673: for (b = 1; b < bs; b++) matis->sf_leafdata[i] = PetscMax(matis->sf_leafdata[i], matis->sf_leafdata[i * bs + b] / bs);
1674: }
1675: PetscCall(MatSeqBAIJSetPreallocation(matis->A, bs, 0, matis->sf_leafdata));
1677: nlocalcols /= bs;
1678: for (i = 0; i < matis->sf->nleaves / bs; i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i], nlocalcols - i);
1679: PetscCall(MatSeqSBAIJSetPreallocation(matis->A, bs, 0, matis->sf_leafdata));
1681: /* for other matrix types */
1682: PetscCall(MatSetUp(matis->A));
1683: PetscFunctionReturn(PETSC_SUCCESS);
1684: }
1686: PETSC_EXTERN PetscErrorCode MatISSetMPIXAIJPreallocation_Private(Mat A, Mat B, PetscBool maxreduce)
1687: {
1688: Mat_IS *matis = (Mat_IS *)(A->data);
1689: PetscInt *my_dnz, *my_onz, *dnz, *onz, *mat_ranges, *row_ownership;
1690: const PetscInt *global_indices_r, *global_indices_c;
1691: PetscInt i, j, bs, rows, cols;
1692: PetscInt lrows, lcols;
1693: PetscInt local_rows, local_cols;
1694: PetscMPIInt size;
1695: PetscBool isdense, issbaij;
1697: PetscFunctionBegin;
1698: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1699: PetscCall(MatGetSize(A, &rows, &cols));
1700: PetscCall(MatGetBlockSize(A, &bs));
1701: PetscCall(MatGetSize(matis->A, &local_rows, &local_cols));
1702: PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQDENSE, &isdense));
1703: PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQSBAIJ, &issbaij));
1704: PetscCall(ISLocalToGlobalMappingGetIndices(matis->rmapping, &global_indices_r));
1705: if (matis->rmapping != matis->cmapping) {
1706: PetscCall(ISLocalToGlobalMappingGetIndices(matis->cmapping, &global_indices_c));
1707: } else global_indices_c = global_indices_r;
1709: if (issbaij) PetscCall(MatGetRowUpperTriangular(matis->A));
1710: /*
1711: An SF reduce is needed to sum up properly on shared rows.
1712: Note that generally preallocation is not exact, since it overestimates nonzeros
1713: */
1714: PetscCall(MatGetLocalSize(A, &lrows, &lcols));
1715: MatPreallocateBegin(PetscObjectComm((PetscObject)A), lrows, lcols, dnz, onz);
1716: /* All processes need to compute entire row ownership */
1717: PetscCall(PetscMalloc1(rows, &row_ownership));
1718: PetscCall(MatGetOwnershipRanges(A, (const PetscInt **)&mat_ranges));
1719: for (i = 0; i < size; i++) {
1720: for (j = mat_ranges[i]; j < mat_ranges[i + 1]; j++) row_ownership[j] = i;
1721: }
1722: PetscCall(MatGetOwnershipRangesColumn(A, (const PetscInt **)&mat_ranges));
1724: /*
1725: my_dnz and my_onz contains exact contribution to preallocation from each local mat
1726: then, they will be summed up properly. This way, preallocation is always sufficient
1727: */
1728: PetscCall(PetscCalloc2(local_rows, &my_dnz, local_rows, &my_onz));
1729: /* preallocation as a MATAIJ */
1730: if (isdense) { /* special case for dense local matrices */
1731: for (i = 0; i < local_rows; i++) {
1732: PetscInt owner = row_ownership[global_indices_r[i]];
1733: for (j = 0; j < local_cols; j++) {
1734: PetscInt index_col = global_indices_c[j];
1735: if (index_col > mat_ranges[owner] - 1 && index_col < mat_ranges[owner + 1]) { /* diag block */
1736: my_dnz[i] += 1;
1737: } else { /* offdiag block */
1738: my_onz[i] += 1;
1739: }
1740: }
1741: }
1742: } else if (matis->A->ops->getrowij) {
1743: const PetscInt *ii, *jj, *jptr;
1744: PetscBool done;
1745: PetscCall(MatGetRowIJ(matis->A, 0, PETSC_FALSE, PETSC_FALSE, &local_rows, &ii, &jj, &done));
1746: PetscCheck(done, PetscObjectComm((PetscObject)(matis->A)), PETSC_ERR_PLIB, "Error in MatGetRowIJ");
1747: jptr = jj;
1748: for (i = 0; i < local_rows; i++) {
1749: PetscInt index_row = global_indices_r[i];
1750: for (j = 0; j < ii[i + 1] - ii[i]; j++, jptr++) {
1751: PetscInt owner = row_ownership[index_row];
1752: PetscInt index_col = global_indices_c[*jptr];
1753: if (index_col > mat_ranges[owner] - 1 && index_col < mat_ranges[owner + 1]) { /* diag block */
1754: my_dnz[i] += 1;
1755: } else { /* offdiag block */
1756: my_onz[i] += 1;
1757: }
1758: /* same as before, interchanging rows and cols */
1759: if (issbaij && index_col != index_row) {
1760: owner = row_ownership[index_col];
1761: if (index_row > mat_ranges[owner] - 1 && index_row < mat_ranges[owner + 1]) {
1762: my_dnz[*jptr] += 1;
1763: } else {
1764: my_onz[*jptr] += 1;
1765: }
1766: }
1767: }
1768: }
1769: PetscCall(MatRestoreRowIJ(matis->A, 0, PETSC_FALSE, PETSC_FALSE, &local_rows, &ii, &jj, &done));
1770: PetscCheck(done, PetscObjectComm((PetscObject)(matis->A)), PETSC_ERR_PLIB, "Error in MatRestoreRowIJ");
1771: } else { /* loop over rows and use MatGetRow */
1772: for (i = 0; i < local_rows; i++) {
1773: const PetscInt *cols;
1774: PetscInt ncols, index_row = global_indices_r[i];
1775: PetscCall(MatGetRow(matis->A, i, &ncols, &cols, NULL));
1776: for (j = 0; j < ncols; j++) {
1777: PetscInt owner = row_ownership[index_row];
1778: PetscInt index_col = global_indices_c[cols[j]];
1779: if (index_col > mat_ranges[owner] - 1 && index_col < mat_ranges[owner + 1]) { /* diag block */
1780: my_dnz[i] += 1;
1781: } else { /* offdiag block */
1782: my_onz[i] += 1;
1783: }
1784: /* same as before, interchanging rows and cols */
1785: if (issbaij && index_col != index_row) {
1786: owner = row_ownership[index_col];
1787: if (index_row > mat_ranges[owner] - 1 && index_row < mat_ranges[owner + 1]) {
1788: my_dnz[cols[j]] += 1;
1789: } else {
1790: my_onz[cols[j]] += 1;
1791: }
1792: }
1793: }
1794: PetscCall(MatRestoreRow(matis->A, i, &ncols, &cols, NULL));
1795: }
1796: }
1797: if (global_indices_c != global_indices_r) PetscCall(ISLocalToGlobalMappingRestoreIndices(matis->cmapping, &global_indices_c));
1798: PetscCall(ISLocalToGlobalMappingRestoreIndices(matis->rmapping, &global_indices_r));
1799: PetscCall(PetscFree(row_ownership));
1801: /* Reduce my_dnz and my_onz */
1802: if (maxreduce) {
1803: PetscCall(PetscSFReduceBegin(matis->sf, MPIU_INT, my_dnz, dnz, MPI_MAX));
1804: PetscCall(PetscSFReduceBegin(matis->sf, MPIU_INT, my_onz, onz, MPI_MAX));
1805: PetscCall(PetscSFReduceEnd(matis->sf, MPIU_INT, my_dnz, dnz, MPI_MAX));
1806: PetscCall(PetscSFReduceEnd(matis->sf, MPIU_INT, my_onz, onz, MPI_MAX));
1807: } else {
1808: PetscCall(PetscSFReduceBegin(matis->sf, MPIU_INT, my_dnz, dnz, MPI_SUM));
1809: PetscCall(PetscSFReduceBegin(matis->sf, MPIU_INT, my_onz, onz, MPI_SUM));
1810: PetscCall(PetscSFReduceEnd(matis->sf, MPIU_INT, my_dnz, dnz, MPI_SUM));
1811: PetscCall(PetscSFReduceEnd(matis->sf, MPIU_INT, my_onz, onz, MPI_SUM));
1812: }
1813: PetscCall(PetscFree2(my_dnz, my_onz));
1815: /* Resize preallocation if overestimated */
1816: for (i = 0; i < lrows; i++) {
1817: dnz[i] = PetscMin(dnz[i], lcols);
1818: onz[i] = PetscMin(onz[i], cols - lcols);
1819: }
1821: /* Set preallocation */
1822: PetscCall(MatSetBlockSizesFromMats(B, A, A));
1823: PetscCall(MatSeqAIJSetPreallocation(B, 0, dnz));
1824: PetscCall(MatMPIAIJSetPreallocation(B, 0, dnz, 0, onz));
1825: for (i = 0; i < lrows; i += bs) {
1826: PetscInt b, d = dnz[i], o = onz[i];
1828: for (b = 1; b < bs; b++) {
1829: d = PetscMax(d, dnz[i + b]);
1830: o = PetscMax(o, onz[i + b]);
1831: }
1832: dnz[i / bs] = PetscMin(d / bs + d % bs, lcols / bs);
1833: onz[i / bs] = PetscMin(o / bs + o % bs, (cols - lcols) / bs);
1834: }
1835: PetscCall(MatSeqBAIJSetPreallocation(B, bs, 0, dnz));
1836: PetscCall(MatMPIBAIJSetPreallocation(B, bs, 0, dnz, 0, onz));
1837: PetscCall(MatMPISBAIJSetPreallocation(B, bs, 0, dnz, 0, onz));
1838: MatPreallocateEnd(dnz, onz);
1839: if (issbaij) PetscCall(MatRestoreRowUpperTriangular(matis->A));
1840: PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
1841: PetscFunctionReturn(PETSC_SUCCESS);
1842: }
1844: PETSC_INTERN PetscErrorCode MatConvert_IS_XAIJ(Mat mat, MatType mtype, MatReuse reuse, Mat *M)
1845: {
1846: Mat_IS *matis = (Mat_IS *)(mat->data);
1847: Mat local_mat, MT;
1848: PetscInt rbs, cbs, rows, cols, lrows, lcols;
1849: PetscInt local_rows, local_cols;
1850: PetscBool isseqdense, isseqsbaij, isseqaij, isseqbaij;
1851: PetscMPIInt size;
1852: const PetscScalar *array;
1854: PetscFunctionBegin;
1855: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)mat), &size));
1856: if (size == 1 && mat->rmap->N == matis->A->rmap->N && mat->cmap->N == matis->A->cmap->N) {
1857: Mat B;
1858: IS irows = NULL, icols = NULL;
1859: PetscInt rbs, cbs;
1861: PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->rmapping, &rbs));
1862: PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->cmapping, &cbs));
1863: if (reuse != MAT_REUSE_MATRIX) { /* check if l2g maps are one-to-one */
1864: IS rows, cols;
1865: const PetscInt *ridxs, *cidxs;
1866: PetscInt i, nw, *work;
1868: PetscCall(ISLocalToGlobalMappingGetBlockIndices(matis->rmapping, &ridxs));
1869: PetscCall(ISLocalToGlobalMappingGetSize(matis->rmapping, &nw));
1870: nw = nw / rbs;
1871: PetscCall(PetscCalloc1(nw, &work));
1872: for (i = 0; i < nw; i++) work[ridxs[i]] += 1;
1873: for (i = 0; i < nw; i++)
1874: if (!work[i] || work[i] > 1) break;
1875: if (i == nw) {
1876: PetscCall(ISCreateBlock(PETSC_COMM_SELF, rbs, nw, ridxs, PETSC_USE_POINTER, &rows));
1877: PetscCall(ISSetPermutation(rows));
1878: PetscCall(ISInvertPermutation(rows, PETSC_DECIDE, &irows));
1879: PetscCall(ISDestroy(&rows));
1880: }
1881: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(matis->rmapping, &ridxs));
1882: PetscCall(PetscFree(work));
1883: if (irows && matis->rmapping != matis->cmapping) {
1884: PetscCall(ISLocalToGlobalMappingGetBlockIndices(matis->cmapping, &cidxs));
1885: PetscCall(ISLocalToGlobalMappingGetSize(matis->cmapping, &nw));
1886: nw = nw / cbs;
1887: PetscCall(PetscCalloc1(nw, &work));
1888: for (i = 0; i < nw; i++) work[cidxs[i]] += 1;
1889: for (i = 0; i < nw; i++)
1890: if (!work[i] || work[i] > 1) break;
1891: if (i == nw) {
1892: PetscCall(ISCreateBlock(PETSC_COMM_SELF, cbs, nw, cidxs, PETSC_USE_POINTER, &cols));
1893: PetscCall(ISSetPermutation(cols));
1894: PetscCall(ISInvertPermutation(cols, PETSC_DECIDE, &icols));
1895: PetscCall(ISDestroy(&cols));
1896: }
1897: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(matis->cmapping, &cidxs));
1898: PetscCall(PetscFree(work));
1899: } else if (irows) {
1900: PetscCall(PetscObjectReference((PetscObject)irows));
1901: icols = irows;
1902: }
1903: } else {
1904: PetscCall(PetscObjectQuery((PetscObject)(*M), "_MatIS_IS_XAIJ_irows", (PetscObject *)&irows));
1905: PetscCall(PetscObjectQuery((PetscObject)(*M), "_MatIS_IS_XAIJ_icols", (PetscObject *)&icols));
1906: if (irows) PetscCall(PetscObjectReference((PetscObject)irows));
1907: if (icols) PetscCall(PetscObjectReference((PetscObject)icols));
1908: }
1909: if (!irows || !icols) {
1910: PetscCall(ISDestroy(&icols));
1911: PetscCall(ISDestroy(&irows));
1912: goto general_assembly;
1913: }
1914: PetscCall(MatConvert(matis->A, mtype, MAT_INITIAL_MATRIX, &B));
1915: if (reuse != MAT_INPLACE_MATRIX) {
1916: PetscCall(MatCreateSubMatrix(B, irows, icols, reuse, M));
1917: PetscCall(PetscObjectCompose((PetscObject)(*M), "_MatIS_IS_XAIJ_irows", (PetscObject)irows));
1918: PetscCall(PetscObjectCompose((PetscObject)(*M), "_MatIS_IS_XAIJ_icols", (PetscObject)icols));
1919: } else {
1920: Mat C;
1922: PetscCall(MatCreateSubMatrix(B, irows, icols, MAT_INITIAL_MATRIX, &C));
1923: PetscCall(MatHeaderReplace(mat, &C));
1924: }
1925: PetscCall(MatDestroy(&B));
1926: PetscCall(ISDestroy(&icols));
1927: PetscCall(ISDestroy(&irows));
1928: PetscFunctionReturn(PETSC_SUCCESS);
1929: }
1930: general_assembly:
1931: PetscCall(MatGetSize(mat, &rows, &cols));
1932: PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->rmapping, &rbs));
1933: PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->cmapping, &cbs));
1934: PetscCall(MatGetLocalSize(mat, &lrows, &lcols));
1935: PetscCall(MatGetSize(matis->A, &local_rows, &local_cols));
1936: PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQDENSE, &isseqdense));
1937: PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQAIJ, &isseqaij));
1938: PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQBAIJ, &isseqbaij));
1939: PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQSBAIJ, &isseqsbaij));
1940: PetscCheck(isseqdense || isseqaij || isseqbaij || isseqsbaij, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for matrix type %s", ((PetscObject)(matis->A))->type_name);
1941: if (PetscDefined(USE_DEBUG)) {
1942: PetscBool lb[4], bb[4];
1944: lb[0] = isseqdense;
1945: lb[1] = isseqaij;
1946: lb[2] = isseqbaij;
1947: lb[3] = isseqsbaij;
1948: PetscCall(MPIU_Allreduce(lb, bb, 4, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)mat)));
1949: PetscCheck(bb[0] || bb[1] || bb[2] || bb[3], PETSC_COMM_SELF, PETSC_ERR_SUP, "Local matrices must have the same type");
1950: }
1952: if (reuse != MAT_REUSE_MATRIX) {
1953: PetscCall(MatCreate(PetscObjectComm((PetscObject)mat), &MT));
1954: PetscCall(MatSetSizes(MT, lrows, lcols, rows, cols));
1955: PetscCall(MatSetType(MT, mtype));
1956: PetscCall(MatSetBlockSizes(MT, rbs, cbs));
1957: PetscCall(MatISSetMPIXAIJPreallocation_Private(mat, MT, PETSC_FALSE));
1958: } else {
1959: PetscInt mrbs, mcbs, mrows, mcols, mlrows, mlcols;
1961: /* some checks */
1962: MT = *M;
1963: PetscCall(MatGetBlockSizes(MT, &mrbs, &mcbs));
1964: PetscCall(MatGetSize(MT, &mrows, &mcols));
1965: PetscCall(MatGetLocalSize(MT, &mlrows, &mlcols));
1966: PetscCheck(mrows == rows, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of rows (%" PetscInt_FMT " != %" PetscInt_FMT ")", rows, mrows);
1967: PetscCheck(mcols == cols, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of cols (%" PetscInt_FMT " != %" PetscInt_FMT ")", cols, mcols);
1968: PetscCheck(mlrows == lrows, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of local rows (%" PetscInt_FMT " != %" PetscInt_FMT ")", lrows, mlrows);
1969: PetscCheck(mlcols == lcols, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of local cols (%" PetscInt_FMT " != %" PetscInt_FMT ")", lcols, mlcols);
1970: PetscCheck(mrbs == rbs, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong row block size (%" PetscInt_FMT " != %" PetscInt_FMT ")", rbs, mrbs);
1971: PetscCheck(mcbs == cbs, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong col block size (%" PetscInt_FMT " != %" PetscInt_FMT ")", cbs, mcbs);
1972: PetscCall(MatZeroEntries(MT));
1973: }
1975: if (isseqsbaij || isseqbaij) {
1976: PetscCall(MatConvert(matis->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &local_mat));
1977: isseqaij = PETSC_TRUE;
1978: } else {
1979: PetscCall(PetscObjectReference((PetscObject)matis->A));
1980: local_mat = matis->A;
1981: }
1983: /* Set values */
1984: PetscCall(MatSetLocalToGlobalMapping(MT, matis->rmapping, matis->cmapping));
1985: if (isseqdense) { /* special case for dense local matrices */
1986: PetscInt i, *dummy;
1988: PetscCall(PetscMalloc1(PetscMax(local_rows, local_cols), &dummy));
1989: for (i = 0; i < PetscMax(local_rows, local_cols); i++) dummy[i] = i;
1990: PetscCall(MatSetOption(MT, MAT_ROW_ORIENTED, PETSC_FALSE));
1991: PetscCall(MatDenseGetArrayRead(local_mat, &array));
1992: PetscCall(MatSetValuesLocal(MT, local_rows, dummy, local_cols, dummy, array, ADD_VALUES));
1993: PetscCall(MatDenseRestoreArrayRead(local_mat, &array));
1994: PetscCall(PetscFree(dummy));
1995: } else if (isseqaij) {
1996: const PetscInt *blocks;
1997: PetscInt i, nvtxs, *xadj, *adjncy, nb;
1998: PetscBool done;
1999: PetscScalar *sarray;
2001: PetscCall(MatGetRowIJ(local_mat, 0, PETSC_FALSE, PETSC_FALSE, &nvtxs, (const PetscInt **)&xadj, (const PetscInt **)&adjncy, &done));
2002: PetscCheck(done, PetscObjectComm((PetscObject)local_mat), PETSC_ERR_PLIB, "Error in MatGetRowIJ");
2003: PetscCall(MatSeqAIJGetArray(local_mat, &sarray));
2004: PetscCall(MatGetVariableBlockSizes(local_mat, &nb, &blocks));
2005: if (nb) { /* speed up assembly for special blocked matrices (used by BDDC) */
2006: PetscInt sum;
2008: for (i = 0, sum = 0; i < nb; i++) sum += blocks[i];
2009: if (sum == nvtxs) {
2010: PetscInt r;
2012: for (i = 0, r = 0; i < nb; i++) {
2013: PetscAssert(blocks[i] == xadj[r + 1] - xadj[r], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid block sizes prescribed for block %" PetscInt_FMT ": expected %" PetscInt_FMT ", got %" PetscInt_FMT, i, blocks[i], xadj[r + 1] - xadj[r]);
2014: PetscCall(MatSetValuesLocal(MT, blocks[i], adjncy + xadj[r], blocks[i], adjncy + xadj[r], sarray + xadj[r], ADD_VALUES));
2015: r += blocks[i];
2016: }
2017: } else {
2018: for (i = 0; i < nvtxs; i++) PetscCall(MatSetValuesLocal(MT, 1, &i, xadj[i + 1] - xadj[i], adjncy + xadj[i], sarray + xadj[i], ADD_VALUES));
2019: }
2020: } else {
2021: for (i = 0; i < nvtxs; i++) PetscCall(MatSetValuesLocal(MT, 1, &i, xadj[i + 1] - xadj[i], adjncy + xadj[i], sarray + xadj[i], ADD_VALUES));
2022: }
2023: PetscCall(MatRestoreRowIJ(local_mat, 0, PETSC_FALSE, PETSC_FALSE, &nvtxs, (const PetscInt **)&xadj, (const PetscInt **)&adjncy, &done));
2024: PetscCheck(done, PetscObjectComm((PetscObject)local_mat), PETSC_ERR_PLIB, "Error in MatRestoreRowIJ");
2025: PetscCall(MatSeqAIJRestoreArray(local_mat, &sarray));
2026: } else { /* very basic values insertion for all other matrix types */
2027: PetscInt i;
2029: for (i = 0; i < local_rows; i++) {
2030: PetscInt j;
2031: const PetscInt *local_indices_cols;
2033: PetscCall(MatGetRow(local_mat, i, &j, &local_indices_cols, &array));
2034: PetscCall(MatSetValuesLocal(MT, 1, &i, j, local_indices_cols, array, ADD_VALUES));
2035: PetscCall(MatRestoreRow(local_mat, i, &j, &local_indices_cols, &array));
2036: }
2037: }
2038: PetscCall(MatAssemblyBegin(MT, MAT_FINAL_ASSEMBLY));
2039: PetscCall(MatDestroy(&local_mat));
2040: PetscCall(MatAssemblyEnd(MT, MAT_FINAL_ASSEMBLY));
2041: if (isseqdense) PetscCall(MatSetOption(MT, MAT_ROW_ORIENTED, PETSC_TRUE));
2042: if (reuse == MAT_INPLACE_MATRIX) {
2043: PetscCall(MatHeaderReplace(mat, &MT));
2044: } else if (reuse == MAT_INITIAL_MATRIX) {
2045: *M = MT;
2046: }
2047: PetscFunctionReturn(PETSC_SUCCESS);
2048: }
2050: /*@
2051: MatISGetMPIXAIJ - Converts `MATIS` matrix into a parallel `MATAIJ` format
2053: Input Parameters:
2054: + mat - the matrix (should be of type `MATIS`)
2055: - reuse - either `MAT_INITIAL_MATRIX` or `MAT_REUSE_MATRIX`
2057: Output Parameter:
2058: . newmat - the matrix in `MATAIJ` format
2060: Level: deprecated
2062: Note:
2063: This function has been deprecated and it will be removed in future releases. Update your code to use the `MatConvert()` interface.
2065: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatConvert()`
2066: @*/
2067: PetscErrorCode MatISGetMPIXAIJ(Mat mat, MatReuse reuse, Mat *newmat)
2068: {
2069: PetscFunctionBegin;
2072: PetscAssertPointer(newmat, 3);
2073: if (reuse == MAT_REUSE_MATRIX) {
2075: PetscCheckSameComm(mat, 1, *newmat, 3);
2076: PetscCheck(mat != *newmat, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse the same matrix");
2077: }
2078: PetscUseMethod(mat, "MatISGetMPIXAIJ_C", (Mat, MatType, MatReuse, Mat *), (mat, MATAIJ, reuse, newmat));
2079: PetscFunctionReturn(PETSC_SUCCESS);
2080: }
2082: static PetscErrorCode MatDuplicate_IS(Mat mat, MatDuplicateOption op, Mat *newmat)
2083: {
2084: Mat_IS *matis = (Mat_IS *)(mat->data);
2085: PetscInt rbs, cbs, m, n, M, N;
2086: Mat B, localmat;
2088: PetscFunctionBegin;
2089: PetscCall(ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping, &rbs));
2090: PetscCall(ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping, &cbs));
2091: PetscCall(MatGetSize(mat, &M, &N));
2092: PetscCall(MatGetLocalSize(mat, &m, &n));
2093: PetscCall(MatCreate(PetscObjectComm((PetscObject)mat), &B));
2094: PetscCall(MatSetSizes(B, m, n, M, N));
2095: PetscCall(MatSetBlockSize(B, rbs == cbs ? rbs : 1));
2096: PetscCall(MatSetType(B, MATIS));
2097: PetscCall(MatISSetLocalMatType(B, matis->lmattype));
2098: PetscCall(MatSetLocalToGlobalMapping(B, mat->rmap->mapping, mat->cmap->mapping));
2099: PetscCall(MatDuplicate(matis->A, op, &localmat));
2100: PetscCall(MatSetLocalToGlobalMapping(localmat, matis->A->rmap->mapping, matis->A->cmap->mapping));
2101: PetscCall(MatISSetLocalMat(B, localmat));
2102: PetscCall(MatDestroy(&localmat));
2103: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
2104: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
2105: *newmat = B;
2106: PetscFunctionReturn(PETSC_SUCCESS);
2107: }
2109: static PetscErrorCode MatIsHermitian_IS(Mat A, PetscReal tol, PetscBool *flg)
2110: {
2111: Mat_IS *matis = (Mat_IS *)A->data;
2112: PetscBool local_sym;
2114: PetscFunctionBegin;
2115: PetscCall(MatIsHermitian(matis->A, tol, &local_sym));
2116: PetscCall(MPIU_Allreduce(&local_sym, flg, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
2117: PetscFunctionReturn(PETSC_SUCCESS);
2118: }
2120: static PetscErrorCode MatIsSymmetric_IS(Mat A, PetscReal tol, PetscBool *flg)
2121: {
2122: Mat_IS *matis = (Mat_IS *)A->data;
2123: PetscBool local_sym;
2125: PetscFunctionBegin;
2126: if (matis->rmapping != matis->cmapping) {
2127: *flg = PETSC_FALSE;
2128: PetscFunctionReturn(PETSC_SUCCESS);
2129: }
2130: PetscCall(MatIsSymmetric(matis->A, tol, &local_sym));
2131: PetscCall(MPIU_Allreduce(&local_sym, flg, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
2132: PetscFunctionReturn(PETSC_SUCCESS);
2133: }
2135: static PetscErrorCode MatIsStructurallySymmetric_IS(Mat A, PetscBool *flg)
2136: {
2137: Mat_IS *matis = (Mat_IS *)A->data;
2138: PetscBool local_sym;
2140: PetscFunctionBegin;
2141: if (matis->rmapping != matis->cmapping) {
2142: *flg = PETSC_FALSE;
2143: PetscFunctionReturn(PETSC_SUCCESS);
2144: }
2145: PetscCall(MatIsStructurallySymmetric(matis->A, &local_sym));
2146: PetscCall(MPIU_Allreduce(&local_sym, flg, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
2147: PetscFunctionReturn(PETSC_SUCCESS);
2148: }
2150: static PetscErrorCode MatDestroy_IS(Mat A)
2151: {
2152: Mat_IS *b = (Mat_IS *)A->data;
2154: PetscFunctionBegin;
2155: PetscCall(PetscFree(b->bdiag));
2156: PetscCall(PetscFree(b->lmattype));
2157: PetscCall(MatDestroy(&b->A));
2158: PetscCall(VecScatterDestroy(&b->cctx));
2159: PetscCall(VecScatterDestroy(&b->rctx));
2160: PetscCall(VecDestroy(&b->x));
2161: PetscCall(VecDestroy(&b->y));
2162: PetscCall(VecDestroy(&b->counter));
2163: PetscCall(ISDestroy(&b->getsub_ris));
2164: PetscCall(ISDestroy(&b->getsub_cis));
2165: if (b->sf != b->csf) {
2166: PetscCall(PetscSFDestroy(&b->csf));
2167: PetscCall(PetscFree2(b->csf_rootdata, b->csf_leafdata));
2168: } else b->csf = NULL;
2169: PetscCall(PetscSFDestroy(&b->sf));
2170: PetscCall(PetscFree2(b->sf_rootdata, b->sf_leafdata));
2171: PetscCall(ISLocalToGlobalMappingDestroy(&b->rmapping));
2172: PetscCall(ISLocalToGlobalMappingDestroy(&b->cmapping));
2173: PetscCall(MatDestroy(&b->dA));
2174: PetscCall(MatDestroy(&b->assembledA));
2175: PetscCall(PetscFree(A->data));
2176: PetscCall(PetscObjectChangeTypeName((PetscObject)A, NULL));
2177: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMatType_C", NULL));
2178: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalMat_C", NULL));
2179: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMat_C", NULL));
2180: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISRestoreLocalMat_C", NULL));
2181: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetMPIXAIJ_C", NULL));
2182: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetPreallocation_C", NULL));
2183: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISStoreL2L_C", NULL));
2184: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISFixLocalEmpty_C", NULL));
2185: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalToGlobalMapping_C", NULL));
2186: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpiaij_C", NULL));
2187: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpibaij_C", NULL));
2188: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpisbaij_C", NULL));
2189: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqaij_C", NULL));
2190: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqbaij_C", NULL));
2191: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqsbaij_C", NULL));
2192: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_aij_C", NULL));
2193: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOOLocal_C", NULL));
2194: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOO_C", NULL));
2195: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", NULL));
2196: PetscFunctionReturn(PETSC_SUCCESS);
2197: }
2199: static PetscErrorCode MatMult_IS(Mat A, Vec x, Vec y)
2200: {
2201: Mat_IS *is = (Mat_IS *)A->data;
2202: PetscScalar zero = 0.0;
2204: PetscFunctionBegin;
2205: /* scatter the global vector x into the local work vector */
2206: PetscCall(VecScatterBegin(is->cctx, x, is->x, INSERT_VALUES, SCATTER_FORWARD));
2207: PetscCall(VecScatterEnd(is->cctx, x, is->x, INSERT_VALUES, SCATTER_FORWARD));
2209: /* multiply the local matrix */
2210: PetscCall(MatMult(is->A, is->x, is->y));
2212: /* scatter product back into global memory */
2213: PetscCall(VecSet(y, zero));
2214: PetscCall(VecScatterBegin(is->rctx, is->y, y, ADD_VALUES, SCATTER_REVERSE));
2215: PetscCall(VecScatterEnd(is->rctx, is->y, y, ADD_VALUES, SCATTER_REVERSE));
2216: PetscFunctionReturn(PETSC_SUCCESS);
2217: }
2219: static PetscErrorCode MatMultAdd_IS(Mat A, Vec v1, Vec v2, Vec v3)
2220: {
2221: Vec temp_vec;
2223: PetscFunctionBegin; /* v3 = v2 + A * v1.*/
2224: if (v3 != v2) {
2225: PetscCall(MatMult(A, v1, v3));
2226: PetscCall(VecAXPY(v3, 1.0, v2));
2227: } else {
2228: PetscCall(VecDuplicate(v2, &temp_vec));
2229: PetscCall(MatMult(A, v1, temp_vec));
2230: PetscCall(VecAXPY(temp_vec, 1.0, v2));
2231: PetscCall(VecCopy(temp_vec, v3));
2232: PetscCall(VecDestroy(&temp_vec));
2233: }
2234: PetscFunctionReturn(PETSC_SUCCESS);
2235: }
2237: static PetscErrorCode MatMultTranspose_IS(Mat A, Vec y, Vec x)
2238: {
2239: Mat_IS *is = (Mat_IS *)A->data;
2241: PetscFunctionBegin;
2242: /* scatter the global vector x into the local work vector */
2243: PetscCall(VecScatterBegin(is->rctx, y, is->y, INSERT_VALUES, SCATTER_FORWARD));
2244: PetscCall(VecScatterEnd(is->rctx, y, is->y, INSERT_VALUES, SCATTER_FORWARD));
2246: /* multiply the local matrix */
2247: PetscCall(MatMultTranspose(is->A, is->y, is->x));
2249: /* scatter product back into global vector */
2250: PetscCall(VecSet(x, 0));
2251: PetscCall(VecScatterBegin(is->cctx, is->x, x, ADD_VALUES, SCATTER_REVERSE));
2252: PetscCall(VecScatterEnd(is->cctx, is->x, x, ADD_VALUES, SCATTER_REVERSE));
2253: PetscFunctionReturn(PETSC_SUCCESS);
2254: }
2256: static PetscErrorCode MatMultTransposeAdd_IS(Mat A, Vec v1, Vec v2, Vec v3)
2257: {
2258: Vec temp_vec;
2260: PetscFunctionBegin; /* v3 = v2 + A' * v1.*/
2261: if (v3 != v2) {
2262: PetscCall(MatMultTranspose(A, v1, v3));
2263: PetscCall(VecAXPY(v3, 1.0, v2));
2264: } else {
2265: PetscCall(VecDuplicate(v2, &temp_vec));
2266: PetscCall(MatMultTranspose(A, v1, temp_vec));
2267: PetscCall(VecAXPY(temp_vec, 1.0, v2));
2268: PetscCall(VecCopy(temp_vec, v3));
2269: PetscCall(VecDestroy(&temp_vec));
2270: }
2271: PetscFunctionReturn(PETSC_SUCCESS);
2272: }
2274: static PetscErrorCode MatView_IS(Mat A, PetscViewer viewer)
2275: {
2276: Mat_IS *a = (Mat_IS *)A->data;
2277: PetscViewer sviewer;
2278: PetscBool isascii, view = PETSC_TRUE;
2280: PetscFunctionBegin;
2281: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2282: if (isascii) {
2283: PetscViewerFormat format;
2285: PetscCall(PetscViewerGetFormat(viewer, &format));
2286: if (format == PETSC_VIEWER_ASCII_INFO) view = PETSC_FALSE;
2287: }
2288: if (!view) PetscFunctionReturn(PETSC_SUCCESS);
2289: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
2290: PetscCall(MatView(a->A, sviewer));
2291: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
2292: PetscCall(PetscViewerFlush(viewer));
2293: PetscFunctionReturn(PETSC_SUCCESS);
2294: }
2296: static PetscErrorCode MatInvertBlockDiagonal_IS(Mat mat, const PetscScalar **values)
2297: {
2298: Mat_IS *is = (Mat_IS *)mat->data;
2299: MPI_Datatype nodeType;
2300: const PetscScalar *lv;
2301: PetscInt bs;
2303: PetscFunctionBegin;
2304: PetscCall(MatGetBlockSize(mat, &bs));
2305: PetscCall(MatSetBlockSize(is->A, bs));
2306: PetscCall(MatInvertBlockDiagonal(is->A, &lv));
2307: if (!is->bdiag) PetscCall(PetscMalloc1(bs * mat->rmap->n, &is->bdiag));
2308: PetscCallMPI(MPI_Type_contiguous(bs, MPIU_SCALAR, &nodeType));
2309: PetscCallMPI(MPI_Type_commit(&nodeType));
2310: PetscCall(PetscSFReduceBegin(is->sf, nodeType, lv, is->bdiag, MPI_REPLACE));
2311: PetscCall(PetscSFReduceEnd(is->sf, nodeType, lv, is->bdiag, MPI_REPLACE));
2312: PetscCallMPI(MPI_Type_free(&nodeType));
2313: if (values) *values = is->bdiag;
2314: PetscFunctionReturn(PETSC_SUCCESS);
2315: }
2317: static PetscErrorCode MatISSetUpScatters_Private(Mat A)
2318: {
2319: Vec cglobal, rglobal;
2320: IS from;
2321: Mat_IS *is = (Mat_IS *)A->data;
2322: PetscScalar sum;
2323: const PetscInt *garray;
2324: PetscInt nr, rbs, nc, cbs;
2325: VecType rtype;
2327: PetscFunctionBegin;
2328: PetscCall(ISLocalToGlobalMappingGetSize(is->rmapping, &nr));
2329: PetscCall(ISLocalToGlobalMappingGetBlockSize(is->rmapping, &rbs));
2330: PetscCall(ISLocalToGlobalMappingGetSize(is->cmapping, &nc));
2331: PetscCall(ISLocalToGlobalMappingGetBlockSize(is->cmapping, &cbs));
2332: PetscCall(VecDestroy(&is->x));
2333: PetscCall(VecDestroy(&is->y));
2334: PetscCall(VecDestroy(&is->counter));
2335: PetscCall(VecScatterDestroy(&is->rctx));
2336: PetscCall(VecScatterDestroy(&is->cctx));
2337: PetscCall(MatCreateVecs(is->A, &is->x, &is->y));
2338: PetscCall(VecBindToCPU(is->y, PETSC_TRUE));
2339: PetscCall(VecGetRootType_Private(is->y, &rtype));
2340: PetscCall(PetscFree(A->defaultvectype));
2341: PetscCall(PetscStrallocpy(rtype, &A->defaultvectype));
2342: PetscCall(MatCreateVecs(A, &cglobal, &rglobal));
2343: PetscCall(VecBindToCPU(rglobal, PETSC_TRUE));
2344: PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->rmapping, &garray));
2345: PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), rbs, nr / rbs, garray, PETSC_USE_POINTER, &from));
2346: PetscCall(VecScatterCreate(rglobal, from, is->y, NULL, &is->rctx));
2347: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->rmapping, &garray));
2348: PetscCall(ISDestroy(&from));
2349: if (is->rmapping != is->cmapping) {
2350: PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->cmapping, &garray));
2351: PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), cbs, nc / cbs, garray, PETSC_USE_POINTER, &from));
2352: PetscCall(VecScatterCreate(cglobal, from, is->x, NULL, &is->cctx));
2353: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->cmapping, &garray));
2354: PetscCall(ISDestroy(&from));
2355: } else {
2356: PetscCall(PetscObjectReference((PetscObject)is->rctx));
2357: is->cctx = is->rctx;
2358: }
2359: PetscCall(VecDestroy(&cglobal));
2361: /* interface counter vector (local) */
2362: PetscCall(VecDuplicate(is->y, &is->counter));
2363: PetscCall(VecBindToCPU(is->counter, PETSC_TRUE));
2364: PetscCall(VecSet(is->y, 1.));
2365: PetscCall(VecScatterBegin(is->rctx, is->y, rglobal, ADD_VALUES, SCATTER_REVERSE));
2366: PetscCall(VecScatterEnd(is->rctx, is->y, rglobal, ADD_VALUES, SCATTER_REVERSE));
2367: PetscCall(VecScatterBegin(is->rctx, rglobal, is->counter, INSERT_VALUES, SCATTER_FORWARD));
2368: PetscCall(VecScatterEnd(is->rctx, rglobal, is->counter, INSERT_VALUES, SCATTER_FORWARD));
2369: PetscCall(VecBindToCPU(is->y, PETSC_FALSE));
2370: PetscCall(VecBindToCPU(is->counter, PETSC_FALSE));
2372: /* special functions for block-diagonal matrices */
2373: PetscCall(VecSum(rglobal, &sum));
2374: A->ops->invertblockdiagonal = NULL;
2375: if ((PetscInt)(PetscRealPart(sum)) == A->rmap->N && A->rmap->N == A->cmap->N && is->rmapping == is->cmapping) A->ops->invertblockdiagonal = MatInvertBlockDiagonal_IS;
2376: PetscCall(VecDestroy(&rglobal));
2378: /* setup SF for general purpose shared indices based communications */
2379: PetscCall(MatISSetUpSF_IS(A));
2380: PetscFunctionReturn(PETSC_SUCCESS);
2381: }
2383: static PetscErrorCode MatISFilterL2GMap(Mat A, ISLocalToGlobalMapping map, ISLocalToGlobalMapping *nmap, ISLocalToGlobalMapping *lmap)
2384: {
2385: IS is;
2386: ISLocalToGlobalMappingType l2gtype;
2387: const PetscInt *idxs;
2388: PetscHSetI ht;
2389: PetscInt *nidxs;
2390: PetscInt i, n, bs, c;
2391: PetscBool flg[] = {PETSC_FALSE, PETSC_FALSE};
2393: PetscFunctionBegin;
2394: PetscCall(ISLocalToGlobalMappingGetSize(map, &n));
2395: PetscCall(ISLocalToGlobalMappingGetBlockSize(map, &bs));
2396: PetscCall(ISLocalToGlobalMappingGetBlockIndices(map, &idxs));
2397: PetscCall(PetscHSetICreate(&ht));
2398: PetscCall(PetscMalloc1(n / bs, &nidxs));
2399: for (i = 0, c = 0; i < n / bs; i++) {
2400: PetscBool missing;
2401: if (idxs[i] < 0) {
2402: flg[0] = PETSC_TRUE;
2403: continue;
2404: }
2405: PetscCall(PetscHSetIQueryAdd(ht, idxs[i], &missing));
2406: if (!missing) flg[1] = PETSC_TRUE;
2407: else nidxs[c++] = idxs[i];
2408: }
2409: PetscCall(PetscHSetIDestroy(&ht));
2410: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, flg, 2, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A)));
2411: if (!flg[0] && !flg[1]) { /* Entries are all non negative and unique */
2412: *nmap = NULL;
2413: *lmap = NULL;
2414: PetscCall(PetscFree(nidxs));
2415: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(map, &idxs));
2416: PetscFunctionReturn(PETSC_SUCCESS);
2417: }
2419: /* New l2g map without negative or repeated indices */
2420: PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), bs, c, nidxs, PETSC_USE_POINTER, &is));
2421: PetscCall(ISLocalToGlobalMappingCreateIS(is, nmap));
2422: PetscCall(ISDestroy(&is));
2423: PetscCall(ISLocalToGlobalMappingGetType(map, &l2gtype));
2424: PetscCall(ISLocalToGlobalMappingSetType(*nmap, l2gtype));
2426: /* New local l2g map for repeated indices */
2427: PetscCall(ISGlobalToLocalMappingApplyBlock(*nmap, IS_GTOLM_MASK, n / bs, idxs, NULL, nidxs));
2428: PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, n / bs, nidxs, PETSC_USE_POINTER, &is));
2429: PetscCall(ISLocalToGlobalMappingCreateIS(is, lmap));
2430: PetscCall(ISDestroy(&is));
2432: PetscCall(PetscFree(nidxs));
2433: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(map, &idxs));
2434: PetscFunctionReturn(PETSC_SUCCESS);
2435: }
2437: static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A, ISLocalToGlobalMapping rmapping, ISLocalToGlobalMapping cmapping)
2438: {
2439: Mat_IS *is = (Mat_IS *)A->data;
2440: ISLocalToGlobalMapping localrmapping = NULL, localcmapping = NULL;
2441: PetscBool cong, freem[] = {PETSC_FALSE, PETSC_FALSE};
2442: PetscInt nr, rbs, nc, cbs;
2444: PetscFunctionBegin;
2445: if (rmapping) PetscCheckSameComm(A, 1, rmapping, 2);
2446: if (cmapping) PetscCheckSameComm(A, 1, cmapping, 3);
2448: PetscCall(ISLocalToGlobalMappingDestroy(&is->rmapping));
2449: PetscCall(ISLocalToGlobalMappingDestroy(&is->cmapping));
2450: PetscCall(PetscLayoutSetUp(A->rmap));
2451: PetscCall(PetscLayoutSetUp(A->cmap));
2452: PetscCall(MatHasCongruentLayouts(A, &cong));
2454: /* If NULL, local space matches global space */
2455: if (!rmapping) {
2456: IS is;
2458: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), A->rmap->N, 0, 1, &is));
2459: PetscCall(ISLocalToGlobalMappingCreateIS(is, &rmapping));
2460: if (A->rmap->bs > 0) PetscCall(ISLocalToGlobalMappingSetBlockSize(rmapping, A->rmap->bs));
2461: PetscCall(ISDestroy(&is));
2462: freem[0] = PETSC_TRUE;
2463: if (!cmapping && cong && A->rmap->bs == A->cmap->bs) cmapping = rmapping;
2464: } else if (!is->islocalref) { /* check if the l2g map has negative or repeated entries */
2465: PetscCall(MatISFilterL2GMap(A, rmapping, &is->rmapping, &localrmapping));
2466: if (rmapping == cmapping) {
2467: PetscCall(PetscObjectReference((PetscObject)is->rmapping));
2468: is->cmapping = is->rmapping;
2469: PetscCall(PetscObjectReference((PetscObject)localrmapping));
2470: localcmapping = localrmapping;
2471: }
2472: }
2473: if (!cmapping) {
2474: IS is;
2476: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), A->cmap->N, 0, 1, &is));
2477: PetscCall(ISLocalToGlobalMappingCreateIS(is, &cmapping));
2478: if (A->cmap->bs > 0) PetscCall(ISLocalToGlobalMappingSetBlockSize(cmapping, A->cmap->bs));
2479: PetscCall(ISDestroy(&is));
2480: freem[1] = PETSC_TRUE;
2481: } else if (cmapping != rmapping && !is->islocalref) { /* check if the l2g map has negative or repeated entries */
2482: PetscCall(MatISFilterL2GMap(A, cmapping, &is->cmapping, &localcmapping));
2483: }
2484: if (!is->rmapping) {
2485: PetscCall(PetscObjectReference((PetscObject)rmapping));
2486: is->rmapping = rmapping;
2487: }
2488: if (!is->cmapping) {
2489: PetscCall(PetscObjectReference((PetscObject)cmapping));
2490: is->cmapping = cmapping;
2491: }
2493: /* Clean up */
2494: PetscCall(MatDestroy(&is->A));
2495: if (is->csf != is->sf) {
2496: PetscCall(PetscSFDestroy(&is->csf));
2497: PetscCall(PetscFree2(is->csf_rootdata, is->csf_leafdata));
2498: } else is->csf = NULL;
2499: PetscCall(PetscSFDestroy(&is->sf));
2500: PetscCall(PetscFree2(is->sf_rootdata, is->sf_leafdata));
2501: PetscCall(PetscFree(is->bdiag));
2503: /* check if the two mappings are actually the same for square matrices since MATIS has some optimization for this case
2504: (DOLFIN passes 2 different objects) */
2505: PetscCall(ISLocalToGlobalMappingGetSize(is->rmapping, &nr));
2506: PetscCall(ISLocalToGlobalMappingGetBlockSize(is->rmapping, &rbs));
2507: PetscCall(ISLocalToGlobalMappingGetSize(is->cmapping, &nc));
2508: PetscCall(ISLocalToGlobalMappingGetBlockSize(is->cmapping, &cbs));
2509: if (is->rmapping != is->cmapping && cong) {
2510: PetscBool same = PETSC_FALSE;
2511: if (nr == nc && cbs == rbs) {
2512: const PetscInt *idxs1, *idxs2;
2514: PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->rmapping, &idxs1));
2515: PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->cmapping, &idxs2));
2516: PetscCall(PetscArraycmp(idxs1, idxs2, nr / rbs, &same));
2517: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->rmapping, &idxs1));
2518: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->cmapping, &idxs2));
2519: }
2520: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &same, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
2521: if (same) {
2522: PetscCall(ISLocalToGlobalMappingDestroy(&is->cmapping));
2523: PetscCall(PetscObjectReference((PetscObject)is->rmapping));
2524: is->cmapping = is->rmapping;
2525: }
2526: }
2527: PetscCall(PetscLayoutSetBlockSize(A->rmap, rbs));
2528: PetscCall(PetscLayoutSetBlockSize(A->cmap, cbs));
2529: /* Pass the user defined maps to the layout */
2530: PetscCall(PetscLayoutSetISLocalToGlobalMapping(A->rmap, rmapping));
2531: PetscCall(PetscLayoutSetISLocalToGlobalMapping(A->cmap, cmapping));
2532: if (freem[0]) PetscCall(ISLocalToGlobalMappingDestroy(&rmapping));
2533: if (freem[1]) PetscCall(ISLocalToGlobalMappingDestroy(&cmapping));
2535: /* Create the local matrix A */
2536: PetscCall(MatCreate(PETSC_COMM_SELF, &is->A));
2537: PetscCall(MatSetType(is->A, is->lmattype));
2538: PetscCall(MatSetSizes(is->A, nr, nc, nr, nc));
2539: PetscCall(MatSetBlockSizes(is->A, rbs, cbs));
2540: PetscCall(MatSetOptionsPrefix(is->A, "is_"));
2541: PetscCall(MatAppendOptionsPrefix(is->A, ((PetscObject)A)->prefix));
2542: PetscCall(PetscLayoutSetUp(is->A->rmap));
2543: PetscCall(PetscLayoutSetUp(is->A->cmap));
2544: PetscCall(MatSetLocalToGlobalMapping(is->A, localrmapping, localcmapping));
2545: PetscCall(ISLocalToGlobalMappingDestroy(&localrmapping));
2546: PetscCall(ISLocalToGlobalMappingDestroy(&localcmapping));
2548: /* setup scatters and local vectors for MatMult */
2549: if (!is->islocalref) PetscCall(MatISSetUpScatters_Private(A));
2550: A->preallocated = PETSC_TRUE;
2551: PetscFunctionReturn(PETSC_SUCCESS);
2552: }
2554: static PetscErrorCode MatSetUp_IS(Mat A)
2555: {
2556: ISLocalToGlobalMapping rmap, cmap;
2558: PetscFunctionBegin;
2559: PetscCall(MatGetLocalToGlobalMapping(A, &rmap, &cmap));
2560: if (!rmap && !cmap) PetscCall(MatSetLocalToGlobalMapping(A, NULL, NULL));
2561: PetscFunctionReturn(PETSC_SUCCESS);
2562: }
2564: static PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
2565: {
2566: Mat_IS *is = (Mat_IS *)mat->data;
2567: PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION], cols_l[MATIS_MAX_ENTRIES_INSERTION];
2569: PetscFunctionBegin;
2570: PetscCall(ISGlobalToLocalMappingApply(is->rmapping, IS_GTOLM_MASK, m, rows, &m, rows_l));
2571: if (m != n || rows != cols || is->cmapping != is->rmapping) {
2572: PetscCall(ISGlobalToLocalMappingApply(is->cmapping, IS_GTOLM_MASK, n, cols, &n, cols_l));
2573: PetscCall(MatSetValues(is->A, m, rows_l, n, cols_l, values, addv));
2574: } else {
2575: PetscCall(MatSetValues(is->A, m, rows_l, m, rows_l, values, addv));
2576: }
2577: PetscFunctionReturn(PETSC_SUCCESS);
2578: }
2580: static PetscErrorCode MatSetValuesBlocked_IS(Mat mat, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
2581: {
2582: Mat_IS *is = (Mat_IS *)mat->data;
2583: PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION], cols_l[MATIS_MAX_ENTRIES_INSERTION];
2585: PetscFunctionBegin;
2586: PetscCall(ISGlobalToLocalMappingApplyBlock(is->rmapping, IS_GTOLM_MASK, m, rows, &m, rows_l));
2587: if (m != n || rows != cols || is->cmapping != is->rmapping) {
2588: PetscCall(ISGlobalToLocalMappingApply(is->cmapping, IS_GTOLM_MASK, n, cols, &n, cols_l));
2589: PetscCall(MatSetValuesBlocked(is->A, m, rows_l, n, cols_l, values, addv));
2590: } else {
2591: PetscCall(MatSetValuesBlocked(is->A, m, rows_l, n, rows_l, values, addv));
2592: }
2593: PetscFunctionReturn(PETSC_SUCCESS);
2594: }
2596: static PetscErrorCode MatSetValuesLocal_IS(Mat A, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
2597: {
2598: Mat_IS *is = (Mat_IS *)A->data;
2600: PetscFunctionBegin;
2601: if (is->A->rmap->mapping || is->A->cmap->mapping) {
2602: PetscCall(MatSetValuesLocal(is->A, m, rows, n, cols, values, addv));
2603: } else {
2604: PetscCall(MatSetValues(is->A, m, rows, n, cols, values, addv));
2605: }
2606: PetscFunctionReturn(PETSC_SUCCESS);
2607: }
2609: static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat A, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
2610: {
2611: Mat_IS *is = (Mat_IS *)A->data;
2613: PetscFunctionBegin;
2614: if (is->A->rmap->mapping || is->A->cmap->mapping) {
2615: PetscCall(MatSetValuesBlockedLocal(is->A, m, rows, n, cols, values, addv));
2616: } else {
2617: PetscCall(MatSetValuesBlocked(is->A, m, rows, n, cols, values, addv));
2618: }
2619: PetscFunctionReturn(PETSC_SUCCESS);
2620: }
2622: static PetscErrorCode MatISZeroRowsColumnsLocal_Private(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, PetscBool columns)
2623: {
2624: Mat_IS *is = (Mat_IS *)A->data;
2626: PetscFunctionBegin;
2627: if (!n) {
2628: is->pure_neumann = PETSC_TRUE;
2629: } else {
2630: PetscInt i;
2631: is->pure_neumann = PETSC_FALSE;
2633: if (columns) {
2634: PetscCall(MatZeroRowsColumns(is->A, n, rows, diag, NULL, NULL));
2635: } else {
2636: PetscCall(MatZeroRows(is->A, n, rows, diag, NULL, NULL));
2637: }
2638: if (diag != 0.) {
2639: const PetscScalar *array;
2640: PetscCall(VecGetArrayRead(is->counter, &array));
2641: for (i = 0; i < n; i++) PetscCall(MatSetValue(is->A, rows[i], rows[i], diag / (array[rows[i]]), INSERT_VALUES));
2642: PetscCall(VecRestoreArrayRead(is->counter, &array));
2643: }
2644: PetscCall(MatAssemblyBegin(is->A, MAT_FINAL_ASSEMBLY));
2645: PetscCall(MatAssemblyEnd(is->A, MAT_FINAL_ASSEMBLY));
2646: }
2647: PetscFunctionReturn(PETSC_SUCCESS);
2648: }
2650: static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b, PetscBool columns)
2651: {
2652: Mat_IS *matis = (Mat_IS *)A->data;
2653: PetscInt nr, nl, len, i;
2654: PetscInt *lrows;
2656: PetscFunctionBegin;
2657: if (PetscUnlikelyDebug(columns || diag != 0. || (x && b))) {
2658: PetscBool cong;
2660: PetscCall(PetscLayoutCompare(A->rmap, A->cmap, &cong));
2661: cong = (PetscBool)(cong && matis->sf == matis->csf);
2662: 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");
2663: 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");
2664: 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");
2665: }
2666: /* get locally owned rows */
2667: PetscCall(PetscLayoutMapLocal(A->rmap, n, rows, &len, &lrows, NULL));
2668: /* fix right hand side if needed */
2669: if (x && b) {
2670: const PetscScalar *xx;
2671: PetscScalar *bb;
2673: PetscCall(VecGetArrayRead(x, &xx));
2674: PetscCall(VecGetArray(b, &bb));
2675: for (i = 0; i < len; ++i) bb[lrows[i]] = diag * xx[lrows[i]];
2676: PetscCall(VecRestoreArrayRead(x, &xx));
2677: PetscCall(VecRestoreArray(b, &bb));
2678: }
2679: /* get rows associated to the local matrices */
2680: PetscCall(MatGetSize(matis->A, &nl, NULL));
2681: PetscCall(PetscArrayzero(matis->sf_leafdata, nl));
2682: PetscCall(PetscArrayzero(matis->sf_rootdata, A->rmap->n));
2683: for (i = 0; i < len; i++) matis->sf_rootdata[lrows[i]] = 1;
2684: PetscCall(PetscFree(lrows));
2685: PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
2686: PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
2687: PetscCall(PetscMalloc1(nl, &lrows));
2688: for (i = 0, nr = 0; i < nl; i++)
2689: if (matis->sf_leafdata[i]) lrows[nr++] = i;
2690: PetscCall(MatISZeroRowsColumnsLocal_Private(A, nr, lrows, diag, columns));
2691: PetscCall(PetscFree(lrows));
2692: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
2693: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
2694: PetscFunctionReturn(PETSC_SUCCESS);
2695: }
2697: static PetscErrorCode MatZeroRows_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
2698: {
2699: PetscFunctionBegin;
2700: PetscCall(MatZeroRowsColumns_Private_IS(A, n, rows, diag, x, b, PETSC_FALSE));
2701: PetscFunctionReturn(PETSC_SUCCESS);
2702: }
2704: static PetscErrorCode MatZeroRowsColumns_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
2705: {
2706: PetscFunctionBegin;
2707: PetscCall(MatZeroRowsColumns_Private_IS(A, n, rows, diag, x, b, PETSC_TRUE));
2708: PetscFunctionReturn(PETSC_SUCCESS);
2709: }
2711: static PetscErrorCode MatAssemblyBegin_IS(Mat A, MatAssemblyType type)
2712: {
2713: Mat_IS *is = (Mat_IS *)A->data;
2715: PetscFunctionBegin;
2716: PetscCall(MatAssemblyBegin(is->A, type));
2717: PetscFunctionReturn(PETSC_SUCCESS);
2718: }
2720: static PetscErrorCode MatAssemblyEnd_IS(Mat A, MatAssemblyType type)
2721: {
2722: Mat_IS *is = (Mat_IS *)A->data;
2723: PetscBool lnnz;
2725: PetscFunctionBegin;
2726: PetscCall(MatAssemblyEnd(is->A, type));
2727: /* fix for local empty rows/cols */
2728: if (is->locempty && type == MAT_FINAL_ASSEMBLY) {
2729: Mat newlA;
2730: ISLocalToGlobalMapping rl2g, cl2g;
2731: IS nzr, nzc;
2732: PetscInt nr, nc, nnzr, nnzc;
2733: PetscBool lnewl2g, newl2g;
2735: PetscCall(MatGetSize(is->A, &nr, &nc));
2736: PetscCall(MatFindNonzeroRowsOrCols_Basic(is->A, PETSC_FALSE, PETSC_SMALL, &nzr));
2737: if (!nzr) PetscCall(ISCreateStride(PetscObjectComm((PetscObject)is->A), nr, 0, 1, &nzr));
2738: PetscCall(MatFindNonzeroRowsOrCols_Basic(is->A, PETSC_TRUE, PETSC_SMALL, &nzc));
2739: if (!nzc) PetscCall(ISCreateStride(PetscObjectComm((PetscObject)is->A), nc, 0, 1, &nzc));
2740: PetscCall(ISGetSize(nzr, &nnzr));
2741: PetscCall(ISGetSize(nzc, &nnzc));
2742: if (nnzr != nr || nnzc != nc) { /* need new global l2g map */
2743: lnewl2g = PETSC_TRUE;
2744: PetscCall(MPIU_Allreduce(&lnewl2g, &newl2g, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A)));
2746: /* extract valid submatrix */
2747: PetscCall(MatCreateSubMatrix(is->A, nzr, nzc, MAT_INITIAL_MATRIX, &newlA));
2748: } else { /* local matrix fully populated */
2749: lnewl2g = PETSC_FALSE;
2750: PetscCall(MPIU_Allreduce(&lnewl2g, &newl2g, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A)));
2751: PetscCall(PetscObjectReference((PetscObject)is->A));
2752: newlA = is->A;
2753: }
2755: /* attach new global l2g map if needed */
2756: if (newl2g) {
2757: IS zr, zc;
2758: const PetscInt *ridxs, *cidxs, *zridxs, *zcidxs;
2759: PetscInt *nidxs, i;
2761: PetscCall(ISComplement(nzr, 0, nr, &zr));
2762: PetscCall(ISComplement(nzc, 0, nc, &zc));
2763: PetscCall(PetscMalloc1(PetscMax(nr, nc), &nidxs));
2764: PetscCall(ISLocalToGlobalMappingGetIndices(is->rmapping, &ridxs));
2765: PetscCall(ISLocalToGlobalMappingGetIndices(is->cmapping, &cidxs));
2766: PetscCall(ISGetIndices(zr, &zridxs));
2767: PetscCall(ISGetIndices(zc, &zcidxs));
2768: PetscCall(ISGetLocalSize(zr, &nnzr));
2769: PetscCall(ISGetLocalSize(zc, &nnzc));
2771: PetscCall(PetscArraycpy(nidxs, ridxs, nr));
2772: for (i = 0; i < nnzr; i++) nidxs[zridxs[i]] = -1;
2773: PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A), 1, nr, nidxs, PETSC_COPY_VALUES, &rl2g));
2774: PetscCall(PetscArraycpy(nidxs, cidxs, nc));
2775: for (i = 0; i < nnzc; i++) nidxs[zcidxs[i]] = -1;
2776: PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A), 1, nc, nidxs, PETSC_COPY_VALUES, &cl2g));
2778: PetscCall(ISRestoreIndices(zr, &zridxs));
2779: PetscCall(ISRestoreIndices(zc, &zcidxs));
2780: PetscCall(ISLocalToGlobalMappingRestoreIndices(is->rmapping, &ridxs));
2781: PetscCall(ISLocalToGlobalMappingRestoreIndices(is->cmapping, &cidxs));
2782: PetscCall(ISDestroy(&nzr));
2783: PetscCall(ISDestroy(&nzc));
2784: PetscCall(ISDestroy(&zr));
2785: PetscCall(ISDestroy(&zc));
2786: PetscCall(PetscFree(nidxs));
2787: PetscCall(MatSetLocalToGlobalMapping(A, rl2g, cl2g));
2788: PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
2789: PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
2790: }
2791: PetscCall(MatISSetLocalMat(A, newlA));
2792: PetscCall(MatDestroy(&newlA));
2793: PetscCall(ISDestroy(&nzr));
2794: PetscCall(ISDestroy(&nzc));
2795: is->locempty = PETSC_FALSE;
2796: }
2797: lnnz = (PetscBool)(is->A->nonzerostate == is->lnnzstate);
2798: is->lnnzstate = is->A->nonzerostate;
2799: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &lnnz, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A)));
2800: if (lnnz) A->nonzerostate++;
2801: PetscFunctionReturn(PETSC_SUCCESS);
2802: }
2804: static PetscErrorCode MatISGetLocalMat_IS(Mat mat, Mat *local)
2805: {
2806: Mat_IS *is = (Mat_IS *)mat->data;
2808: PetscFunctionBegin;
2809: *local = is->A;
2810: PetscFunctionReturn(PETSC_SUCCESS);
2811: }
2813: static PetscErrorCode MatISRestoreLocalMat_IS(Mat mat, Mat *local)
2814: {
2815: PetscFunctionBegin;
2816: *local = NULL;
2817: PetscFunctionReturn(PETSC_SUCCESS);
2818: }
2820: /*@
2821: MatISGetLocalMat - Gets the local matrix stored inside a `MATIS` matrix.
2823: Input Parameter:
2824: . mat - the matrix
2826: Output Parameter:
2827: . local - the local matrix
2829: Level: advanced
2831: Notes:
2832: This can be called if you have precomputed the nonzero structure of the
2833: matrix and want to provide it to the inner matrix object to improve the performance
2834: of the `MatSetValues()` operation.
2836: Call `MatISRestoreLocalMat()` when finished with the local matrix.
2838: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatISRestoreLocalMat()`
2839: @*/
2840: PetscErrorCode MatISGetLocalMat(Mat mat, Mat *local)
2841: {
2842: PetscFunctionBegin;
2844: PetscAssertPointer(local, 2);
2845: PetscUseMethod(mat, "MatISGetLocalMat_C", (Mat, Mat *), (mat, local));
2846: PetscFunctionReturn(PETSC_SUCCESS);
2847: }
2849: /*@
2850: MatISRestoreLocalMat - Restores the local matrix obtained with `MatISGetLocalMat()`
2852: Input Parameters:
2853: + mat - the matrix
2854: - local - the local matrix
2856: Level: advanced
2858: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatISGetLocalMat()`
2859: @*/
2860: PetscErrorCode MatISRestoreLocalMat(Mat mat, Mat *local)
2861: {
2862: PetscFunctionBegin;
2864: PetscAssertPointer(local, 2);
2865: PetscUseMethod(mat, "MatISRestoreLocalMat_C", (Mat, Mat *), (mat, local));
2866: PetscFunctionReturn(PETSC_SUCCESS);
2867: }
2869: static PetscErrorCode MatISSetLocalMatType_IS(Mat mat, MatType mtype)
2870: {
2871: Mat_IS *is = (Mat_IS *)mat->data;
2873: PetscFunctionBegin;
2874: if (is->A) PetscCall(MatSetType(is->A, mtype));
2875: PetscCall(PetscFree(is->lmattype));
2876: PetscCall(PetscStrallocpy(mtype, &is->lmattype));
2877: PetscFunctionReturn(PETSC_SUCCESS);
2878: }
2880: /*@C
2881: MatISSetLocalMatType - Specifies the type of local matrix inside the `MATIS`
2883: Input Parameters:
2884: + mat - the matrix
2885: - mtype - the local matrix type
2887: Level: advanced
2889: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatSetType()`, `MatType`
2890: @*/
2891: PetscErrorCode MatISSetLocalMatType(Mat mat, MatType mtype)
2892: {
2893: PetscFunctionBegin;
2895: PetscUseMethod(mat, "MatISSetLocalMatType_C", (Mat, MatType), (mat, mtype));
2896: PetscFunctionReturn(PETSC_SUCCESS);
2897: }
2899: static PetscErrorCode MatISSetLocalMat_IS(Mat mat, Mat local)
2900: {
2901: Mat_IS *is = (Mat_IS *)mat->data;
2902: PetscInt nrows, ncols, orows, ocols;
2903: MatType mtype, otype;
2904: PetscBool sametype = PETSC_TRUE;
2906: PetscFunctionBegin;
2907: if (is->A && !is->islocalref) {
2908: PetscCall(MatGetSize(is->A, &orows, &ocols));
2909: PetscCall(MatGetSize(local, &nrows, &ncols));
2910: PetscCheck(orows == nrows && ocols == ncols, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local MATIS matrix should be of size %" PetscInt_FMT "x%" PetscInt_FMT " (you passed a %" PetscInt_FMT "x%" PetscInt_FMT " matrix)", orows, ocols, nrows, ncols);
2911: PetscCall(MatGetType(local, &mtype));
2912: PetscCall(MatGetType(is->A, &otype));
2913: PetscCall(PetscStrcmp(mtype, otype, &sametype));
2914: }
2915: PetscCall(PetscObjectReference((PetscObject)local));
2916: PetscCall(MatDestroy(&is->A));
2917: is->A = local;
2918: PetscCall(MatGetType(is->A, &mtype));
2919: PetscCall(MatISSetLocalMatType(mat, mtype));
2920: if (!sametype && !is->islocalref) PetscCall(MatISSetUpScatters_Private(mat));
2921: PetscFunctionReturn(PETSC_SUCCESS);
2922: }
2924: /*@
2925: MatISSetLocalMat - Replace the local matrix stored inside a `MATIS` object.
2927: Collective
2929: Input Parameters:
2930: + mat - the matrix
2931: - local - the local matrix
2933: Level: advanced
2935: Notes:
2936: Any previous matrix within the `MATIS` has its reference count decreased by one.
2938: This can be called if you have precomputed the local matrix and
2939: want to provide it to the matrix object `MATIS`.
2941: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatISSetLocalMatType`, `MatISGetLocalMat()`
2942: @*/
2943: PetscErrorCode MatISSetLocalMat(Mat mat, Mat local)
2944: {
2945: PetscFunctionBegin;
2948: PetscUseMethod(mat, "MatISSetLocalMat_C", (Mat, Mat), (mat, local));
2949: PetscFunctionReturn(PETSC_SUCCESS);
2950: }
2952: static PetscErrorCode MatZeroEntries_IS(Mat A)
2953: {
2954: Mat_IS *a = (Mat_IS *)A->data;
2956: PetscFunctionBegin;
2957: PetscCall(MatZeroEntries(a->A));
2958: PetscFunctionReturn(PETSC_SUCCESS);
2959: }
2961: static PetscErrorCode MatScale_IS(Mat A, PetscScalar a)
2962: {
2963: Mat_IS *is = (Mat_IS *)A->data;
2965: PetscFunctionBegin;
2966: PetscCall(MatScale(is->A, a));
2967: PetscFunctionReturn(PETSC_SUCCESS);
2968: }
2970: static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v)
2971: {
2972: Mat_IS *is = (Mat_IS *)A->data;
2974: PetscFunctionBegin;
2975: /* get diagonal of the local matrix */
2976: PetscCall(MatGetDiagonal(is->A, is->y));
2978: /* scatter diagonal back into global vector */
2979: PetscCall(VecSet(v, 0));
2980: PetscCall(VecScatterBegin(is->rctx, is->y, v, ADD_VALUES, SCATTER_REVERSE));
2981: PetscCall(VecScatterEnd(is->rctx, is->y, v, ADD_VALUES, SCATTER_REVERSE));
2982: PetscFunctionReturn(PETSC_SUCCESS);
2983: }
2985: static PetscErrorCode MatSetOption_IS(Mat A, MatOption op, PetscBool flg)
2986: {
2987: Mat_IS *a = (Mat_IS *)A->data;
2989: PetscFunctionBegin;
2990: PetscCall(MatSetOption(a->A, op, flg));
2991: PetscFunctionReturn(PETSC_SUCCESS);
2992: }
2994: static PetscErrorCode MatAXPY_IS(Mat Y, PetscScalar a, Mat X, MatStructure str)
2995: {
2996: Mat_IS *y = (Mat_IS *)Y->data;
2997: Mat_IS *x;
2999: PetscFunctionBegin;
3000: if (PetscDefined(USE_DEBUG)) {
3001: PetscBool ismatis;
3002: PetscCall(PetscObjectTypeCompare((PetscObject)X, MATIS, &ismatis));
3003: PetscCheck(ismatis, PetscObjectComm((PetscObject)Y), PETSC_ERR_SUP, "Cannot call MatAXPY(Y,a,X,str) with X not of type MATIS");
3004: }
3005: x = (Mat_IS *)X->data;
3006: PetscCall(MatAXPY(y->A, a, x->A, str));
3007: PetscFunctionReturn(PETSC_SUCCESS);
3008: }
3010: static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A, IS row, IS col, Mat *submat)
3011: {
3012: Mat lA;
3013: Mat_IS *matis = (Mat_IS *)(A->data);
3014: ISLocalToGlobalMapping rl2g, cl2g;
3015: IS is;
3016: const PetscInt *rg, *rl;
3017: PetscInt nrg;
3018: PetscInt N, M, nrl, i, *idxs;
3020: PetscFunctionBegin;
3021: PetscCall(ISLocalToGlobalMappingGetIndices(A->rmap->mapping, &rg));
3022: PetscCall(ISGetLocalSize(row, &nrl));
3023: PetscCall(ISGetIndices(row, &rl));
3024: PetscCall(ISLocalToGlobalMappingGetSize(A->rmap->mapping, &nrg));
3025: if (PetscDefined(USE_DEBUG)) {
3026: 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 then maximum possible %" PetscInt_FMT, i, rl[i], nrg);
3027: }
3028: PetscCall(PetscMalloc1(nrg, &idxs));
3029: /* map from [0,nrl) to row */
3030: for (i = 0; i < nrl; i++) idxs[i] = rl[i];
3031: for (i = nrl; i < nrg; i++) idxs[i] = -1;
3032: PetscCall(ISRestoreIndices(row, &rl));
3033: PetscCall(ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping, &rg));
3034: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A), nrg, idxs, PETSC_OWN_POINTER, &is));
3035: PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g));
3036: PetscCall(ISDestroy(&is));
3037: /* compute new l2g map for columns */
3038: if (col != row || matis->rmapping != matis->cmapping || matis->A->rmap->mapping != matis->A->cmap->mapping) {
3039: const PetscInt *cg, *cl;
3040: PetscInt ncg;
3041: PetscInt ncl;
3043: PetscCall(ISLocalToGlobalMappingGetIndices(A->cmap->mapping, &cg));
3044: PetscCall(ISGetLocalSize(col, &ncl));
3045: PetscCall(ISGetIndices(col, &cl));
3046: PetscCall(ISLocalToGlobalMappingGetSize(A->cmap->mapping, &ncg));
3047: if (PetscDefined(USE_DEBUG)) {
3048: 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 then maximum possible %" PetscInt_FMT, i, cl[i], ncg);
3049: }
3050: PetscCall(PetscMalloc1(ncg, &idxs));
3051: /* map from [0,ncl) to col */
3052: for (i = 0; i < ncl; i++) idxs[i] = cl[i];
3053: for (i = ncl; i < ncg; i++) idxs[i] = -1;
3054: PetscCall(ISRestoreIndices(col, &cl));
3055: PetscCall(ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping, &cg));
3056: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A), ncg, idxs, PETSC_OWN_POINTER, &is));
3057: PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g));
3058: PetscCall(ISDestroy(&is));
3059: } else {
3060: PetscCall(PetscObjectReference((PetscObject)rl2g));
3061: cl2g = rl2g;
3062: }
3063: /* create the MATIS submatrix */
3064: PetscCall(MatGetSize(A, &M, &N));
3065: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), submat));
3066: PetscCall(MatSetSizes(*submat, PETSC_DECIDE, PETSC_DECIDE, M, N));
3067: PetscCall(MatSetType(*submat, MATIS));
3068: matis = (Mat_IS *)((*submat)->data);
3069: matis->islocalref = PETSC_TRUE;
3070: PetscCall(MatSetLocalToGlobalMapping(*submat, rl2g, cl2g));
3071: PetscCall(MatISGetLocalMat(A, &lA));
3072: PetscCall(MatISSetLocalMat(*submat, lA));
3073: PetscCall(MatSetUp(*submat));
3074: PetscCall(MatAssemblyBegin(*submat, MAT_FINAL_ASSEMBLY));
3075: PetscCall(MatAssemblyEnd(*submat, MAT_FINAL_ASSEMBLY));
3076: PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
3077: PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
3079: /* remove unsupported ops */
3080: PetscCall(PetscMemzero((*submat)->ops, sizeof(struct _MatOps)));
3081: (*submat)->ops->destroy = MatDestroy_IS;
3082: (*submat)->ops->setvalueslocal = MatSetValuesLocal_SubMat_IS;
3083: (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS;
3084: (*submat)->ops->assemblybegin = MatAssemblyBegin_IS;
3085: (*submat)->ops->assemblyend = MatAssemblyEnd_IS;
3086: PetscFunctionReturn(PETSC_SUCCESS);
3087: }
3089: static PetscErrorCode MatSetFromOptions_IS(Mat A, PetscOptionItems *PetscOptionsObject)
3090: {
3091: Mat_IS *a = (Mat_IS *)A->data;
3092: char type[256];
3093: PetscBool flg;
3095: PetscFunctionBegin;
3096: PetscOptionsHeadBegin(PetscOptionsObject, "MATIS options");
3097: PetscCall(PetscOptionsBool("-matis_keepassembled", "Store an assembled version if needed", "MatISKeepAssembled", a->keepassembled, &a->keepassembled, NULL));
3098: PetscCall(PetscOptionsBool("-matis_fixempty", "Fix local matrices in case of empty local rows/columns", "MatISFixLocalEmpty", a->locempty, &a->locempty, NULL));
3099: PetscCall(PetscOptionsBool("-matis_storel2l", "Store local-to-local matrices generated from PtAP operations", "MatISStoreL2L", a->storel2l, &a->storel2l, NULL));
3100: PetscCall(PetscOptionsFList("-matis_localmat_type", "Matrix type", "MatISSetLocalMatType", MatList, a->lmattype, type, 256, &flg));
3101: if (flg) PetscCall(MatISSetLocalMatType(A, type));
3102: if (a->A) PetscCall(MatSetFromOptions(a->A));
3103: PetscOptionsHeadEnd();
3104: PetscFunctionReturn(PETSC_SUCCESS);
3105: }
3107: /*@
3108: MatCreateIS - Creates a "process" unassembled matrix, `MATIS`, assembled on each
3109: process but not across processes.
3111: Input Parameters:
3112: + comm - MPI communicator that will share the matrix
3113: . bs - block size of the matrix
3114: . m - local size of left vector used in matrix vector products
3115: . n - local size of right vector used in matrix vector products
3116: . M - global size of left vector used in matrix vector products
3117: . N - global size of right vector used in matrix vector products
3118: . rmap - local to global map for rows
3119: - cmap - local to global map for cols
3121: Output Parameter:
3122: . A - the resulting matrix
3124: Level: advanced
3126: Notes:
3127: `m` and `n` are NOT related to the size of the map; they represent the size of the local parts of the distributed vectors
3128: used in `MatMult()` operations. The sizes of rmap and cmap define the size of the local matrices.
3130: If `rmap` (`cmap`) is `NULL`, then the local row (column) spaces matches the global space.
3132: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatSetLocalToGlobalMapping()`
3133: @*/
3134: PetscErrorCode MatCreateIS(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt M, PetscInt N, ISLocalToGlobalMapping rmap, ISLocalToGlobalMapping cmap, Mat *A)
3135: {
3136: PetscFunctionBegin;
3137: PetscCall(MatCreate(comm, A));
3138: PetscCall(MatSetSizes(*A, m, n, M, N));
3139: if (bs > 0) PetscCall(MatSetBlockSize(*A, bs));
3140: PetscCall(MatSetType(*A, MATIS));
3141: PetscCall(MatSetLocalToGlobalMapping(*A, rmap, cmap));
3142: PetscFunctionReturn(PETSC_SUCCESS);
3143: }
3145: static PetscErrorCode MatHasOperation_IS(Mat A, MatOperation op, PetscBool *has)
3146: {
3147: Mat_IS *a = (Mat_IS *)A->data;
3148: MatOperation tobefiltered[] = {MATOP_MULT_ADD, MATOP_MULT_TRANSPOSE_ADD, MATOP_GET_DIAGONAL_BLOCK, MATOP_INCREASE_OVERLAP};
3150: PetscFunctionBegin;
3151: *has = PETSC_FALSE;
3152: if (!((void **)A->ops)[op] || !a->A) PetscFunctionReturn(PETSC_SUCCESS);
3153: *has = PETSC_TRUE;
3154: for (PetscInt i = 0; i < (PetscInt)PETSC_STATIC_ARRAY_LENGTH(tobefiltered); i++)
3155: if (op == tobefiltered[i]) PetscFunctionReturn(PETSC_SUCCESS);
3156: PetscCall(MatHasOperation(a->A, op, has));
3157: PetscFunctionReturn(PETSC_SUCCESS);
3158: }
3160: static PetscErrorCode MatSetValuesCOO_IS(Mat A, const PetscScalar v[], InsertMode imode)
3161: {
3162: Mat_IS *a = (Mat_IS *)A->data;
3164: PetscFunctionBegin;
3165: PetscCall(MatSetValuesCOO(a->A, v, imode));
3166: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
3167: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
3168: PetscFunctionReturn(PETSC_SUCCESS);
3169: }
3171: static PetscErrorCode MatSetPreallocationCOOLocal_IS(Mat A, PetscCount ncoo, PetscInt coo_i[], PetscInt coo_j[])
3172: {
3173: Mat_IS *a = (Mat_IS *)A->data;
3175: PetscFunctionBegin;
3176: PetscCheck(a->A, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to provide l2g map first via MatSetLocalToGlobalMapping");
3177: if (a->A->rmap->mapping || a->A->cmap->mapping) {
3178: PetscCall(MatSetPreallocationCOOLocal(a->A, ncoo, coo_i, coo_j));
3179: } else {
3180: PetscCall(MatSetPreallocationCOO(a->A, ncoo, coo_i, coo_j));
3181: }
3182: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", MatSetValuesCOO_IS));
3183: A->preallocated = PETSC_TRUE;
3184: PetscFunctionReturn(PETSC_SUCCESS);
3185: }
3187: static PetscErrorCode MatSetPreallocationCOO_IS(Mat A, PetscCount ncoo, PetscInt coo_i[], PetscInt coo_j[])
3188: {
3189: Mat_IS *a = (Mat_IS *)A->data;
3191: PetscFunctionBegin;
3192: PetscCheck(a->A, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to provide l2g map first via MatSetLocalToGlobalMapping");
3193: PetscCheck(ncoo <= PETSC_MAX_INT, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "ncoo %" PetscCount_FMT " overflowed PetscInt; configure --with-64-bit-indices or request support", ncoo);
3194: PetscCall(ISGlobalToLocalMappingApply(a->rmapping, IS_GTOLM_MASK, ncoo, coo_i, NULL, coo_i));
3195: PetscCall(ISGlobalToLocalMappingApply(a->cmapping, IS_GTOLM_MASK, ncoo, coo_j, NULL, coo_j));
3196: PetscCall(MatSetPreallocationCOO(a->A, ncoo, coo_i, coo_j));
3197: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", MatSetValuesCOO_IS));
3198: A->preallocated = PETSC_TRUE;
3199: PetscFunctionReturn(PETSC_SUCCESS);
3200: }
3202: static PetscErrorCode MatISGetAssembled_Private(Mat A, Mat *tA)
3203: {
3204: Mat_IS *a = (Mat_IS *)A->data;
3205: PetscObjectState Astate, aAstate = PETSC_MIN_INT;
3206: PetscObjectState Annzstate, aAnnzstate = PETSC_MIN_INT;
3208: PetscFunctionBegin;
3209: PetscCall(PetscObjectStateGet((PetscObject)A, &Astate));
3210: Annzstate = A->nonzerostate;
3211: if (a->assembledA) {
3212: PetscCall(PetscObjectStateGet((PetscObject)a->assembledA, &aAstate));
3213: aAnnzstate = a->assembledA->nonzerostate;
3214: }
3215: if (aAnnzstate != Annzstate) PetscCall(MatDestroy(&a->assembledA));
3216: if (Astate != aAstate || !a->assembledA) {
3217: MatType aAtype;
3218: PetscMPIInt size;
3219: PetscInt rbs, cbs, bs;
3221: /* the assembled form is used as temporary storage for parallel operations
3222: like createsubmatrices and the like, do not waste device memory */
3223: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
3224: PetscCall(ISLocalToGlobalMappingGetBlockSize(a->cmapping, &cbs));
3225: PetscCall(ISLocalToGlobalMappingGetBlockSize(a->rmapping, &rbs));
3226: bs = rbs == cbs ? rbs : 1;
3227: if (a->assembledA) PetscCall(MatGetType(a->assembledA, &aAtype));
3228: else if (size > 1) aAtype = bs > 1 ? MATMPIBAIJ : MATMPIAIJ;
3229: else aAtype = bs > 1 ? MATSEQBAIJ : MATSEQAIJ;
3231: PetscCall(MatConvert(A, aAtype, a->assembledA ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX, &a->assembledA));
3232: PetscCall(PetscObjectStateSet((PetscObject)a->assembledA, Astate));
3233: a->assembledA->nonzerostate = Annzstate;
3234: }
3235: PetscCall(PetscObjectReference((PetscObject)a->assembledA));
3236: *tA = a->assembledA;
3237: if (!a->keepassembled) PetscCall(MatDestroy(&a->assembledA));
3238: PetscFunctionReturn(PETSC_SUCCESS);
3239: }
3241: static PetscErrorCode MatISRestoreAssembled_Private(Mat A, Mat *tA)
3242: {
3243: PetscFunctionBegin;
3244: PetscCall(MatDestroy(tA));
3245: *tA = NULL;
3246: PetscFunctionReturn(PETSC_SUCCESS);
3247: }
3249: static PetscErrorCode MatGetDiagonalBlock_IS(Mat A, Mat *dA)
3250: {
3251: Mat_IS *a = (Mat_IS *)A->data;
3252: PetscObjectState Astate, dAstate = PETSC_MIN_INT;
3254: PetscFunctionBegin;
3255: PetscCall(PetscObjectStateGet((PetscObject)A, &Astate));
3256: if (a->dA) PetscCall(PetscObjectStateGet((PetscObject)a->dA, &dAstate));
3257: if (Astate != dAstate) {
3258: Mat tA;
3259: MatType ltype;
3261: PetscCall(MatDestroy(&a->dA));
3262: PetscCall(MatISGetAssembled_Private(A, &tA));
3263: PetscCall(MatGetDiagonalBlock(tA, &a->dA));
3264: PetscCall(MatPropagateSymmetryOptions(tA, a->dA));
3265: PetscCall(MatGetType(a->A, <ype));
3266: PetscCall(MatConvert(a->dA, ltype, MAT_INPLACE_MATRIX, &a->dA));
3267: PetscCall(PetscObjectReference((PetscObject)a->dA));
3268: PetscCall(MatISRestoreAssembled_Private(A, &tA));
3269: PetscCall(PetscObjectStateSet((PetscObject)a->dA, Astate));
3270: }
3271: *dA = a->dA;
3272: PetscFunctionReturn(PETSC_SUCCESS);
3273: }
3275: static PetscErrorCode MatCreateSubMatrices_IS(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse reuse, Mat *submat[])
3276: {
3277: Mat tA;
3279: PetscFunctionBegin;
3280: PetscCall(MatISGetAssembled_Private(A, &tA));
3281: PetscCall(MatCreateSubMatrices(tA, n, irow, icol, reuse, submat));
3282: /* MatCreateSubMatrices_MPIAIJ is a mess at the moment */
3283: #if 0
3284: {
3285: Mat_IS *a = (Mat_IS*)A->data;
3286: MatType ltype;
3287: VecType vtype;
3288: char *flg;
3290: PetscCall(MatGetType(a->A,<ype));
3291: PetscCall(MatGetVecType(a->A,&vtype));
3292: PetscCall(PetscStrstr(vtype,"cuda",&flg));
3293: if (!flg) PetscCall(PetscStrstr(vtype,"hip",&flg));
3294: if (!flg) PetscCall(PetscStrstr(vtype,"kokkos",&flg));
3295: if (flg) {
3296: for (PetscInt i = 0; i < n; i++) {
3297: Mat sA = (*submat)[i];
3299: PetscCall(MatConvert(sA,ltype,MAT_INPLACE_MATRIX,&sA));
3300: (*submat)[i] = sA;
3301: }
3302: }
3303: }
3304: #endif
3305: PetscCall(MatISRestoreAssembled_Private(A, &tA));
3306: PetscFunctionReturn(PETSC_SUCCESS);
3307: }
3309: static PetscErrorCode MatIncreaseOverlap_IS(Mat A, PetscInt n, IS is[], PetscInt ov)
3310: {
3311: Mat tA;
3313: PetscFunctionBegin;
3314: PetscCall(MatISGetAssembled_Private(A, &tA));
3315: PetscCall(MatIncreaseOverlap(tA, n, is, ov));
3316: PetscCall(MatISRestoreAssembled_Private(A, &tA));
3317: PetscFunctionReturn(PETSC_SUCCESS);
3318: }
3320: /*@
3321: MatISGetLocalToGlobalMapping - Gets the local-to-global numbering of the `MATIS` object
3323: Not Collective
3325: Input Parameter:
3326: . A - the matrix
3328: Output Parameters:
3329: + rmapping - row mapping
3330: - cmapping - column mapping
3332: Level: advanced
3334: Note:
3335: The returned map can be different from the one used to construct the `MATIS` object, since it will not contain negative or repeated indices.
3337: .seealso: [](ch_matrices), `Mat`, `MatIS`, `MatSetLocalToGlobalMapping()`
3338: @*/
3339: PetscErrorCode MatISGetLocalToGlobalMapping(Mat A, ISLocalToGlobalMapping *rmapping, ISLocalToGlobalMapping *cmapping)
3340: {
3341: PetscFunctionBegin;
3344: if (rmapping) PetscAssertPointer(rmapping, 2);
3345: if (cmapping) PetscAssertPointer(cmapping, 3);
3346: PetscUseMethod(A, "MatISGetLocalToGlobalMapping_C", (Mat, ISLocalToGlobalMapping *, ISLocalToGlobalMapping *), (A, rmapping, cmapping));
3347: PetscFunctionReturn(PETSC_SUCCESS);
3348: }
3350: static PetscErrorCode MatISGetLocalToGlobalMapping_IS(Mat A, ISLocalToGlobalMapping *r, ISLocalToGlobalMapping *c)
3351: {
3352: Mat_IS *a = (Mat_IS *)A->data;
3354: PetscFunctionBegin;
3355: if (r) *r = a->rmapping;
3356: if (c) *c = a->cmapping;
3357: PetscFunctionReturn(PETSC_SUCCESS);
3358: }
3360: /*MC
3361: MATIS - MATIS = "is" - A matrix type to be used for using the non-overlapping domain decomposition methods (e.g. `PCBDDC` or `KSPFETIDP`).
3362: This stores the matrices in globally unassembled form. Each processor assembles only its local Neumann problem and the parallel matrix vector
3363: product is handled "implicitly".
3365: Options Database Keys:
3366: + -mat_type is - sets the matrix type to `MATIS`.
3367: . -matis_fixempty - Fixes local matrices in case of empty local rows/columns.
3368: - -matis_storel2l - stores the local-to-local operators generated by the Galerkin process of `MatPtAP()`.
3370: Level: advanced
3372: Notes:
3373: Options prefix for the inner matrix are given by `-is_mat_xxx`
3375: You must call `MatSetLocalToGlobalMapping()` before using this matrix type.
3377: You can do matrix preallocation on the local matrix after you obtain it with
3378: `MatISGetLocalMat()`; otherwise, you could use `MatISSetPreallocation()`
3380: .seealso: [](ch_matrices), `Mat`, `MATIS`, `Mat`, `MatISGetLocalMat()`, `MatSetLocalToGlobalMapping()`, `MatISSetPreallocation()`, `MatCreateIS()`, `PCBDDC`, `KSPFETIDP`
3381: M*/
3382: PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A)
3383: {
3384: Mat_IS *a;
3386: PetscFunctionBegin;
3387: PetscCall(PetscNew(&a));
3388: PetscCall(PetscStrallocpy(MATAIJ, &a->lmattype));
3389: A->data = (void *)a;
3391: /* matrix ops */
3392: PetscCall(PetscMemzero(A->ops, sizeof(struct _MatOps)));
3393: A->ops->mult = MatMult_IS;
3394: A->ops->multadd = MatMultAdd_IS;
3395: A->ops->multtranspose = MatMultTranspose_IS;
3396: A->ops->multtransposeadd = MatMultTransposeAdd_IS;
3397: A->ops->destroy = MatDestroy_IS;
3398: A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS;
3399: A->ops->setvalues = MatSetValues_IS;
3400: A->ops->setvaluesblocked = MatSetValuesBlocked_IS;
3401: A->ops->setvalueslocal = MatSetValuesLocal_IS;
3402: A->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_IS;
3403: A->ops->zerorows = MatZeroRows_IS;
3404: A->ops->zerorowscolumns = MatZeroRowsColumns_IS;
3405: A->ops->assemblybegin = MatAssemblyBegin_IS;
3406: A->ops->assemblyend = MatAssemblyEnd_IS;
3407: A->ops->view = MatView_IS;
3408: A->ops->zeroentries = MatZeroEntries_IS;
3409: A->ops->scale = MatScale_IS;
3410: A->ops->getdiagonal = MatGetDiagonal_IS;
3411: A->ops->setoption = MatSetOption_IS;
3412: A->ops->ishermitian = MatIsHermitian_IS;
3413: A->ops->issymmetric = MatIsSymmetric_IS;
3414: A->ops->isstructurallysymmetric = MatIsStructurallySymmetric_IS;
3415: A->ops->duplicate = MatDuplicate_IS;
3416: A->ops->missingdiagonal = MatMissingDiagonal_IS;
3417: A->ops->copy = MatCopy_IS;
3418: A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_IS;
3419: A->ops->createsubmatrix = MatCreateSubMatrix_IS;
3420: A->ops->axpy = MatAXPY_IS;
3421: A->ops->diagonalset = MatDiagonalSet_IS;
3422: A->ops->shift = MatShift_IS;
3423: A->ops->transpose = MatTranspose_IS;
3424: A->ops->getinfo = MatGetInfo_IS;
3425: A->ops->diagonalscale = MatDiagonalScale_IS;
3426: A->ops->setfromoptions = MatSetFromOptions_IS;
3427: A->ops->setup = MatSetUp_IS;
3428: A->ops->hasoperation = MatHasOperation_IS;
3429: A->ops->getdiagonalblock = MatGetDiagonalBlock_IS;
3430: A->ops->createsubmatrices = MatCreateSubMatrices_IS;
3431: A->ops->increaseoverlap = MatIncreaseOverlap_IS;
3433: /* special MATIS functions */
3434: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMatType_C", MatISSetLocalMatType_IS));
3435: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalMat_C", MatISGetLocalMat_IS));
3436: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISRestoreLocalMat_C", MatISRestoreLocalMat_IS));
3437: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMat_C", MatISSetLocalMat_IS));
3438: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetMPIXAIJ_C", MatConvert_IS_XAIJ));
3439: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetPreallocation_C", MatISSetPreallocation_IS));
3440: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISStoreL2L_C", MatISStoreL2L_IS));
3441: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISFixLocalEmpty_C", MatISFixLocalEmpty_IS));
3442: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalToGlobalMapping_C", MatISGetLocalToGlobalMapping_IS));
3443: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpiaij_C", MatConvert_IS_XAIJ));
3444: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpibaij_C", MatConvert_IS_XAIJ));
3445: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpisbaij_C", MatConvert_IS_XAIJ));
3446: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqaij_C", MatConvert_IS_XAIJ));
3447: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqbaij_C", MatConvert_IS_XAIJ));
3448: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqsbaij_C", MatConvert_IS_XAIJ));
3449: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_aij_C", MatConvert_IS_XAIJ));
3450: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOOLocal_C", MatSetPreallocationCOOLocal_IS));
3451: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOO_C", MatSetPreallocationCOO_IS));
3452: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATIS));
3453: PetscFunctionReturn(PETSC_SUCCESS);
3454: }