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