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