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