Actual source code: matis.c

  1: /*
  2:     Creates a matrix class for using the Neumann-Neumann type preconditioners.
  3:     This stores the matrices in globally unassembled form. Each processor
  4:     assembles only its local Neumann problem and the parallel matrix vector
  5:     product is handled "implicitly".

  7:     Currently this allows for only one subdomain per processor.
  8: */

 10: #include <petsc/private/matisimpl.h>
 11: #include <../src/mat/impls/aij/mpi/mpiaij.h>
 12: #include <petsc/private/sfimpl.h>
 13: #include <petsc/private/vecimpl.h>
 14: #include <petsc/private/hashseti.h>

 16: #define MATIS_MAX_ENTRIES_INSERTION 2048
 17: static PetscErrorCode MatSetValuesLocal_IS(Mat, PetscInt, const PetscInt *, PetscInt, const PetscInt *, const PetscScalar *, InsertMode);
 18: static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat, PetscInt, const PetscInt *, PetscInt, const PetscInt *, const PetscScalar *, InsertMode);
 19: static PetscErrorCode MatISSetUpScatters_Private(Mat);

 21: static PetscErrorCode MatISContainerDestroyPtAP_Private(void *ptr)
 22: {
 23:   MatISPtAP ptap = (MatISPtAP)ptr;

 25:   PetscFunctionBegin;
 26:   PetscCall(MatDestroySubMatrices(ptap->ris1 ? 2 : 1, &ptap->lP));
 27:   PetscCall(ISDestroy(&ptap->cis0));
 28:   PetscCall(ISDestroy(&ptap->cis1));
 29:   PetscCall(ISDestroy(&ptap->ris0));
 30:   PetscCall(ISDestroy(&ptap->ris1));
 31:   PetscCall(PetscFree(ptap));
 32:   PetscFunctionReturn(PETSC_SUCCESS);
 33: }

 35: static PetscErrorCode MatPtAPNumeric_IS_XAIJ(Mat A, Mat P, Mat C)
 36: {
 37:   MatISPtAP      ptap;
 38:   Mat_IS        *matis = (Mat_IS *)A->data;
 39:   Mat            lA, lC;
 40:   MatReuse       reuse;
 41:   IS             ris[2], cis[2];
 42:   PetscContainer c;
 43:   PetscInt       n;

 45:   PetscFunctionBegin;
 46:   PetscCall(PetscObjectQuery((PetscObject)C, "_MatIS_PtAP", (PetscObject *)&c));
 47:   PetscCheck(c, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Missing PtAP information");
 48:   PetscCall(PetscContainerGetPointer(c, (void **)&ptap));
 49:   ris[0] = ptap->ris0;
 50:   ris[1] = ptap->ris1;
 51:   cis[0] = ptap->cis0;
 52:   cis[1] = ptap->cis1;
 53:   n      = ptap->ris1 ? 2 : 1;
 54:   reuse  = ptap->lP ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX;
 55:   PetscCall(MatCreateSubMatrices(P, n, ris, cis, reuse, &ptap->lP));

 57:   PetscCall(MatISGetLocalMat(A, &lA));
 58:   PetscCall(MatISGetLocalMat(C, &lC));
 59:   if (ptap->ris1) { /* unsymmetric A mapping */
 60:     Mat lPt;

 62:     PetscCall(MatTranspose(ptap->lP[1], MAT_INITIAL_MATRIX, &lPt));
 63:     PetscCall(MatMatMatMult(lPt, lA, ptap->lP[0], reuse, ptap->fill, &lC));
 64:     if (matis->storel2l) PetscCall(PetscObjectCompose((PetscObject)(A), "_MatIS_PtAP_l2l", (PetscObject)lPt));
 65:     PetscCall(MatDestroy(&lPt));
 66:   } else {
 67:     PetscCall(MatPtAP(lA, ptap->lP[0], reuse, ptap->fill, &lC));
 68:     if (matis->storel2l) PetscCall(PetscObjectCompose((PetscObject)C, "_MatIS_PtAP_l2l", (PetscObject)ptap->lP[0]));
 69:   }
 70:   if (reuse == MAT_INITIAL_MATRIX) {
 71:     PetscCall(MatISSetLocalMat(C, lC));
 72:     PetscCall(MatDestroy(&lC));
 73:   }
 74:   PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
 75:   PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
 76:   PetscFunctionReturn(PETSC_SUCCESS);
 77: }

 79: static PetscErrorCode MatGetNonzeroColumnsLocal_Private(Mat PT, IS *cis)
 80: {
 81:   Mat             Po, Pd;
 82:   IS              zd, zo;
 83:   const PetscInt *garray;
 84:   PetscInt       *aux, i, bs;
 85:   PetscInt        dc, stc, oc, ctd, cto;
 86:   PetscBool       ismpiaij, ismpibaij, isseqaij, isseqbaij;
 87:   MPI_Comm        comm;

 89:   PetscFunctionBegin;
 91:   PetscAssertPointer(cis, 2);
 92:   PetscCall(PetscObjectGetComm((PetscObject)PT, &comm));
 93:   bs = 1;
 94:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)PT, MATMPIAIJ, &ismpiaij));
 95:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)PT, MATMPIBAIJ, &ismpibaij));
 96:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)PT, MATSEQAIJ, &isseqaij));
 97:   PetscCall(PetscObjectTypeCompare((PetscObject)PT, MATSEQBAIJ, &isseqbaij));
 98:   if (isseqaij || isseqbaij) {
 99:     Pd     = PT;
100:     Po     = NULL;
101:     garray = NULL;
102:   } else if (ismpiaij) {
103:     PetscCall(MatMPIAIJGetSeqAIJ(PT, &Pd, &Po, &garray));
104:   } else if (ismpibaij) {
105:     PetscCall(MatMPIBAIJGetSeqBAIJ(PT, &Pd, &Po, &garray));
106:     PetscCall(MatGetBlockSize(PT, &bs));
107:   } else SETERRQ(comm, PETSC_ERR_SUP, "Not for matrix type %s", ((PetscObject)(PT))->type_name);

109:   /* identify any null columns in Pd or Po */
110:   /* We use a tolerance comparison since it may happen that, with geometric multigrid,
111:      some of the columns are not really zero, but very close to */
112:   zo = zd = NULL;
113:   if (Po) PetscCall(MatFindNonzeroRowsOrCols_Basic(Po, PETSC_TRUE, PETSC_SMALL, &zo));
114:   PetscCall(MatFindNonzeroRowsOrCols_Basic(Pd, PETSC_TRUE, PETSC_SMALL, &zd));

116:   PetscCall(MatGetLocalSize(PT, NULL, &dc));
117:   PetscCall(MatGetOwnershipRangeColumn(PT, &stc, NULL));
118:   if (Po) PetscCall(MatGetLocalSize(Po, NULL, &oc));
119:   else oc = 0;
120:   PetscCall(PetscMalloc1((dc + oc) / bs, &aux));
121:   if (zd) {
122:     const PetscInt *idxs;
123:     PetscInt        nz;

125:     /* this will throw an error if bs is not valid */
126:     PetscCall(ISSetBlockSize(zd, bs));
127:     PetscCall(ISGetLocalSize(zd, &nz));
128:     PetscCall(ISGetIndices(zd, &idxs));
129:     ctd = nz / bs;
130:     for (i = 0; i < ctd; i++) aux[i] = (idxs[bs * i] + stc) / bs;
131:     PetscCall(ISRestoreIndices(zd, &idxs));
132:   } else {
133:     ctd = dc / bs;
134:     for (i = 0; i < ctd; i++) aux[i] = i + stc / bs;
135:   }
136:   if (zo) {
137:     const PetscInt *idxs;
138:     PetscInt        nz;

140:     /* this will throw an error if bs is not valid */
141:     PetscCall(ISSetBlockSize(zo, bs));
142:     PetscCall(ISGetLocalSize(zo, &nz));
143:     PetscCall(ISGetIndices(zo, &idxs));
144:     cto = nz / bs;
145:     for (i = 0; i < cto; i++) aux[i + ctd] = garray[idxs[bs * i] / bs];
146:     PetscCall(ISRestoreIndices(zo, &idxs));
147:   } else {
148:     cto = oc / bs;
149:     for (i = 0; i < cto; i++) aux[i + ctd] = garray[i];
150:   }
151:   PetscCall(ISCreateBlock(comm, bs, ctd + cto, aux, PETSC_OWN_POINTER, cis));
152:   PetscCall(ISDestroy(&zd));
153:   PetscCall(ISDestroy(&zo));
154:   PetscFunctionReturn(PETSC_SUCCESS);
155: }

157: static PetscErrorCode MatPtAPSymbolic_IS_XAIJ(Mat A, Mat P, PetscReal fill, Mat C)
158: {
159:   Mat                    PT, lA;
160:   MatISPtAP              ptap;
161:   ISLocalToGlobalMapping Crl2g, Ccl2g, rl2g, cl2g;
162:   PetscContainer         c;
163:   MatType                lmtype;
164:   const PetscInt        *garray;
165:   PetscInt               ibs, N, dc;
166:   MPI_Comm               comm;

168:   PetscFunctionBegin;
169:   PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
170:   PetscCall(MatSetType(C, MATIS));
171:   PetscCall(MatISGetLocalMat(A, &lA));
172:   PetscCall(MatGetType(lA, &lmtype));
173:   PetscCall(MatISSetLocalMatType(C, lmtype));
174:   PetscCall(MatGetSize(P, NULL, &N));
175:   PetscCall(MatGetLocalSize(P, NULL, &dc));
176:   PetscCall(MatSetSizes(C, dc, dc, N, N));
177:   /* Not sure about this
178:   PetscCall(MatGetBlockSizes(P,NULL,&ibs));
179:   PetscCall(MatSetBlockSize(*C,ibs));
180: */

182:   PetscCall(PetscNew(&ptap));
183:   PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &c));
184:   PetscCall(PetscContainerSetPointer(c, ptap));
185:   PetscCall(PetscContainerSetUserDestroy(c, MatISContainerDestroyPtAP_Private));
186:   PetscCall(PetscObjectCompose((PetscObject)C, "_MatIS_PtAP", (PetscObject)c));
187:   PetscCall(PetscContainerDestroy(&c));
188:   ptap->fill = fill;

190:   PetscCall(MatISGetLocalToGlobalMapping(A, &rl2g, &cl2g));

192:   PetscCall(ISLocalToGlobalMappingGetBlockSize(cl2g, &ibs));
193:   PetscCall(ISLocalToGlobalMappingGetSize(cl2g, &N));
194:   PetscCall(ISLocalToGlobalMappingGetBlockIndices(cl2g, &garray));
195:   PetscCall(ISCreateBlock(comm, ibs, N / ibs, garray, PETSC_COPY_VALUES, &ptap->ris0));
196:   PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(cl2g, &garray));

198:   PetscCall(MatCreateSubMatrix(P, ptap->ris0, NULL, MAT_INITIAL_MATRIX, &PT));
199:   PetscCall(MatGetNonzeroColumnsLocal_Private(PT, &ptap->cis0));
200:   PetscCall(ISLocalToGlobalMappingCreateIS(ptap->cis0, &Ccl2g));
201:   PetscCall(MatDestroy(&PT));

203:   Crl2g = NULL;
204:   if (rl2g != cl2g) { /* unsymmetric A mapping */
205:     PetscBool same, lsame = PETSC_FALSE;
206:     PetscInt  N1, ibs1;

208:     PetscCall(ISLocalToGlobalMappingGetSize(rl2g, &N1));
209:     PetscCall(ISLocalToGlobalMappingGetBlockSize(rl2g, &ibs1));
210:     PetscCall(ISLocalToGlobalMappingGetBlockIndices(rl2g, &garray));
211:     PetscCall(ISCreateBlock(comm, ibs, N1 / ibs, garray, PETSC_COPY_VALUES, &ptap->ris1));
212:     PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(rl2g, &garray));
213:     if (ibs1 == ibs && N1 == N) { /* check if the l2gmaps are the same */
214:       const PetscInt *i1, *i2;

216:       PetscCall(ISBlockGetIndices(ptap->ris0, &i1));
217:       PetscCall(ISBlockGetIndices(ptap->ris1, &i2));
218:       PetscCall(PetscArraycmp(i1, i2, N, &lsame));
219:     }
220:     PetscCall(MPIU_Allreduce(&lsame, &same, 1, MPIU_BOOL, MPI_LAND, comm));
221:     if (same) {
222:       PetscCall(ISDestroy(&ptap->ris1));
223:     } else {
224:       PetscCall(MatCreateSubMatrix(P, ptap->ris1, NULL, MAT_INITIAL_MATRIX, &PT));
225:       PetscCall(MatGetNonzeroColumnsLocal_Private(PT, &ptap->cis1));
226:       PetscCall(ISLocalToGlobalMappingCreateIS(ptap->cis1, &Crl2g));
227:       PetscCall(MatDestroy(&PT));
228:     }
229:   }
230:   /* Not sure about this
231:   if (!Crl2g) {
232:     PetscCall(MatGetBlockSize(C,&ibs));
233:     PetscCall(ISLocalToGlobalMappingSetBlockSize(Ccl2g,ibs));
234:   }
235: */
236:   PetscCall(MatSetLocalToGlobalMapping(C, Crl2g ? Crl2g : Ccl2g, Ccl2g));
237:   PetscCall(ISLocalToGlobalMappingDestroy(&Crl2g));
238:   PetscCall(ISLocalToGlobalMappingDestroy(&Ccl2g));

240:   C->ops->ptapnumeric = MatPtAPNumeric_IS_XAIJ;
241:   PetscFunctionReturn(PETSC_SUCCESS);
242: }

244: static PetscErrorCode MatProductSymbolic_PtAP_IS_XAIJ(Mat C)
245: {
246:   Mat_Product *product = C->product;
247:   Mat          A = product->A, P = product->B;
248:   PetscReal    fill = product->fill;

250:   PetscFunctionBegin;
251:   PetscCall(MatPtAPSymbolic_IS_XAIJ(A, P, fill, C));
252:   C->ops->productnumeric = MatProductNumeric_PtAP;
253:   PetscFunctionReturn(PETSC_SUCCESS);
254: }

256: static PetscErrorCode MatProductSetFromOptions_IS_XAIJ_PtAP(Mat C)
257: {
258:   PetscFunctionBegin;
259:   C->ops->productsymbolic = MatProductSymbolic_PtAP_IS_XAIJ;
260:   PetscFunctionReturn(PETSC_SUCCESS);
261: }

263: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_IS_XAIJ(Mat C)
264: {
265:   Mat_Product *product = C->product;

267:   PetscFunctionBegin;
268:   if (product->type == MATPRODUCT_PtAP) PetscCall(MatProductSetFromOptions_IS_XAIJ_PtAP(C));
269:   PetscFunctionReturn(PETSC_SUCCESS);
270: }

272: static PetscErrorCode MatISContainerDestroyFields_Private(void *ptr)
273: {
274:   MatISLocalFields lf = (MatISLocalFields)ptr;
275:   PetscInt         i;

277:   PetscFunctionBegin;
278:   for (i = 0; i < lf->nr; i++) PetscCall(ISDestroy(&lf->rf[i]));
279:   for (i = 0; i < lf->nc; i++) PetscCall(ISDestroy(&lf->cf[i]));
280:   PetscCall(PetscFree2(lf->rf, lf->cf));
281:   PetscCall(PetscFree(lf));
282:   PetscFunctionReturn(PETSC_SUCCESS);
283: }

285: static PetscErrorCode MatConvert_SeqXAIJ_IS(Mat A, MatType type, MatReuse reuse, Mat *newmat)
286: {
287:   Mat B, lB;

289:   PetscFunctionBegin;
290:   if (reuse != MAT_REUSE_MATRIX) {
291:     ISLocalToGlobalMapping rl2g, cl2g;
292:     PetscInt               bs;
293:     IS                     is;

295:     PetscCall(MatGetBlockSize(A, &bs));
296:     PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), A->rmap->n / bs, 0, 1, &is));
297:     if (bs > 1) {
298:       IS       is2;
299:       PetscInt i, *aux;

301:       PetscCall(ISGetLocalSize(is, &i));
302:       PetscCall(ISGetIndices(is, (const PetscInt **)&aux));
303:       PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), bs, i, aux, PETSC_COPY_VALUES, &is2));
304:       PetscCall(ISRestoreIndices(is, (const PetscInt **)&aux));
305:       PetscCall(ISDestroy(&is));
306:       is = is2;
307:     }
308:     PetscCall(ISSetIdentity(is));
309:     PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g));
310:     PetscCall(ISDestroy(&is));
311:     PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), A->cmap->n / bs, 0, 1, &is));
312:     if (bs > 1) {
313:       IS       is2;
314:       PetscInt i, *aux;

316:       PetscCall(ISGetLocalSize(is, &i));
317:       PetscCall(ISGetIndices(is, (const PetscInt **)&aux));
318:       PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), bs, i, aux, PETSC_COPY_VALUES, &is2));
319:       PetscCall(ISRestoreIndices(is, (const PetscInt **)&aux));
320:       PetscCall(ISDestroy(&is));
321:       is = is2;
322:     }
323:     PetscCall(ISSetIdentity(is));
324:     PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g));
325:     PetscCall(ISDestroy(&is));
326:     PetscCall(MatCreateIS(PetscObjectComm((PetscObject)A), bs, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N, rl2g, cl2g, &B));
327:     PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
328:     PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
329:     PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &lB));
330:     if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
331:   } else {
332:     B = *newmat;
333:     PetscCall(PetscObjectReference((PetscObject)A));
334:     lB = A;
335:   }
336:   PetscCall(MatISSetLocalMat(B, lB));
337:   PetscCall(MatDestroy(&lB));
338:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
339:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
340:   if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &B));
341:   PetscFunctionReturn(PETSC_SUCCESS);
342: }

344: static PetscErrorCode MatISScaleDisassembling_Private(Mat A)
345: {
346:   Mat_IS         *matis = (Mat_IS *)A->data;
347:   PetscScalar    *aa;
348:   const PetscInt *ii, *jj;
349:   PetscInt        i, n, m;
350:   PetscInt       *ecount, **eneighs;
351:   PetscBool       flg;

353:   PetscFunctionBegin;
354:   PetscCall(MatGetRowIJ(matis->A, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg));
355:   PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure");
356:   PetscCall(ISLocalToGlobalMappingGetNodeInfo(matis->rmapping, &n, &ecount, &eneighs));
357:   PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected %" PetscInt_FMT " != %" PetscInt_FMT, m, n);
358:   PetscCall(MatSeqAIJGetArray(matis->A, &aa));
359:   for (i = 0; i < n; i++) {
360:     if (ecount[i] > 1) {
361:       PetscInt j;

363:       for (j = ii[i]; j < ii[i + 1]; j++) {
364:         PetscInt  i2   = jj[j], p, p2;
365:         PetscReal scal = 0.0;

367:         for (p = 0; p < ecount[i]; p++) {
368:           for (p2 = 0; p2 < ecount[i2]; p2++) {
369:             if (eneighs[i][p] == eneighs[i2][p2]) {
370:               scal += 1.0;
371:               break;
372:             }
373:           }
374:         }
375:         if (scal) aa[j] /= scal;
376:       }
377:     }
378:   }
379:   PetscCall(ISLocalToGlobalMappingRestoreNodeInfo(matis->rmapping, &n, &ecount, &eneighs));
380:   PetscCall(MatSeqAIJRestoreArray(matis->A, &aa));
381:   PetscCall(MatRestoreRowIJ(matis->A, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg));
382:   PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot restore IJ structure");
383:   PetscFunctionReturn(PETSC_SUCCESS);
384: }

386: typedef enum {
387:   MAT_IS_DISASSEMBLE_L2G_NATURAL,
388:   MAT_IS_DISASSEMBLE_L2G_MAT,
389:   MAT_IS_DISASSEMBLE_L2G_ND
390: } MatISDisassemblel2gType;

392: static PetscErrorCode MatMPIXAIJComputeLocalToGlobalMapping_Private(Mat A, ISLocalToGlobalMapping *l2g)
393: {
394:   Mat                     Ad, Ao;
395:   IS                      is, ndmap, ndsub;
396:   MPI_Comm                comm;
397:   const PetscInt         *garray, *ndmapi;
398:   PetscInt                bs, i, cnt, nl, *ncount, *ndmapc;
399:   PetscBool               ismpiaij, ismpibaij;
400:   const char *const       MatISDisassemblel2gTypes[] = {"NATURAL", "MAT", "ND", "MatISDisassemblel2gType", "MAT_IS_DISASSEMBLE_L2G_", NULL};
401:   MatISDisassemblel2gType mode                       = MAT_IS_DISASSEMBLE_L2G_NATURAL;
402:   MatPartitioning         part;
403:   PetscSF                 sf;
404:   PetscObject             dm;

406:   PetscFunctionBegin;
407:   PetscOptionsBegin(PetscObjectComm((PetscObject)A), ((PetscObject)A)->prefix, "MatIS l2g disassembling options", "Mat");
408:   PetscCall(PetscOptionsEnum("-mat_is_disassemble_l2g_type", "Type of local-to-global mapping to be used for disassembling", "MatISDisassemblel2gType", MatISDisassemblel2gTypes, (PetscEnum)mode, (PetscEnum *)&mode, NULL));
409:   PetscOptionsEnd();
410:   if (mode == MAT_IS_DISASSEMBLE_L2G_MAT) {
411:     PetscCall(MatGetLocalToGlobalMapping(A, l2g, NULL));
412:     PetscFunctionReturn(PETSC_SUCCESS);
413:   }
414:   PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
415:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij));
416:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIBAIJ, &ismpibaij));
417:   PetscCall(MatGetBlockSize(A, &bs));
418:   switch (mode) {
419:   case MAT_IS_DISASSEMBLE_L2G_ND:
420:     PetscCall(MatPartitioningCreate(comm, &part));
421:     PetscCall(MatPartitioningSetAdjacency(part, A));
422:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)part, ((PetscObject)A)->prefix));
423:     PetscCall(MatPartitioningSetFromOptions(part));
424:     PetscCall(MatPartitioningApplyND(part, &ndmap));
425:     PetscCall(MatPartitioningDestroy(&part));
426:     PetscCall(ISBuildTwoSided(ndmap, NULL, &ndsub));
427:     PetscCall(MatMPIAIJSetUseScalableIncreaseOverlap(A, PETSC_TRUE));
428:     PetscCall(MatIncreaseOverlap(A, 1, &ndsub, 1));
429:     PetscCall(ISLocalToGlobalMappingCreateIS(ndsub, l2g));

431:     /* it may happen that a separator node is not properly shared */
432:     PetscCall(ISLocalToGlobalMappingGetNodeInfo(*l2g, &nl, &ncount, NULL));
433:     PetscCall(PetscSFCreate(comm, &sf));
434:     PetscCall(ISLocalToGlobalMappingGetIndices(*l2g, &garray));
435:     PetscCall(PetscSFSetGraphLayout(sf, A->rmap, nl, NULL, PETSC_OWN_POINTER, garray));
436:     PetscCall(ISLocalToGlobalMappingRestoreIndices(*l2g, &garray));
437:     PetscCall(PetscCalloc1(A->rmap->n, &ndmapc));
438:     PetscCall(PetscSFReduceBegin(sf, MPIU_INT, ncount, ndmapc, MPI_REPLACE));
439:     PetscCall(PetscSFReduceEnd(sf, MPIU_INT, ncount, ndmapc, MPI_REPLACE));
440:     PetscCall(ISLocalToGlobalMappingRestoreNodeInfo(*l2g, NULL, &ncount, NULL));
441:     PetscCall(ISGetIndices(ndmap, &ndmapi));
442:     for (i = 0, cnt = 0; i < A->rmap->n; i++)
443:       if (ndmapi[i] < 0 && ndmapc[i] < 2) cnt++;

445:     PetscCall(MPIU_Allreduce(&cnt, &i, 1, MPIU_INT, MPI_MAX, comm));
446:     if (i) { /* we detected isolated separator nodes */
447:       Mat                    A2, A3;
448:       IS                    *workis, is2;
449:       PetscScalar           *vals;
450:       PetscInt               gcnt = i, *dnz, *onz, j, *lndmapi;
451:       ISLocalToGlobalMapping ll2g;
452:       PetscBool              flg;
453:       const PetscInt        *ii, *jj;

455:       /* communicate global id of separators */
456:       MatPreallocateBegin(comm, A->rmap->n, A->cmap->n, dnz, onz);
457:       for (i = 0, cnt = 0; i < A->rmap->n; i++) dnz[i] = ndmapi[i] < 0 ? i + A->rmap->rstart : -1;

459:       PetscCall(PetscMalloc1(nl, &lndmapi));
460:       PetscCall(PetscSFBcastBegin(sf, MPIU_INT, dnz, lndmapi, MPI_REPLACE));

462:       /* compute adjacency of isolated separators node */
463:       PetscCall(PetscMalloc1(gcnt, &workis));
464:       for (i = 0, cnt = 0; i < A->rmap->n; i++) {
465:         if (ndmapi[i] < 0 && ndmapc[i] < 2) PetscCall(ISCreateStride(comm, 1, i + A->rmap->rstart, 1, &workis[cnt++]));
466:       }
467:       for (i = cnt; i < gcnt; i++) PetscCall(ISCreateStride(comm, 0, 0, 1, &workis[i]));
468:       for (i = 0; i < gcnt; i++) {
469:         PetscCall(PetscObjectSetName((PetscObject)workis[i], "ISOLATED"));
470:         PetscCall(ISViewFromOptions(workis[i], NULL, "-view_isolated_separators"));
471:       }

473:       /* no communications since all the ISes correspond to locally owned rows */
474:       PetscCall(MatIncreaseOverlap(A, gcnt, workis, 1));

476:       /* end communicate global id of separators */
477:       PetscCall(PetscSFBcastEnd(sf, MPIU_INT, dnz, lndmapi, MPI_REPLACE));

479:       /* communicate new layers : create a matrix and transpose it */
480:       PetscCall(PetscArrayzero(dnz, A->rmap->n));
481:       PetscCall(PetscArrayzero(onz, A->rmap->n));
482:       for (i = 0, j = 0; i < A->rmap->n; i++) {
483:         if (ndmapi[i] < 0 && ndmapc[i] < 2) {
484:           const PetscInt *idxs;
485:           PetscInt        s;

487:           PetscCall(ISGetLocalSize(workis[j], &s));
488:           PetscCall(ISGetIndices(workis[j], &idxs));
489:           PetscCall(MatPreallocateSet(i + A->rmap->rstart, s, idxs, dnz, onz));
490:           j++;
491:         }
492:       }
493:       PetscCheck(j == cnt, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected local count %" PetscInt_FMT " != %" PetscInt_FMT, j, cnt);

495:       for (i = 0; i < gcnt; i++) {
496:         PetscCall(PetscObjectSetName((PetscObject)workis[i], "EXTENDED"));
497:         PetscCall(ISViewFromOptions(workis[i], NULL, "-view_isolated_separators"));
498:       }

500:       for (i = 0, j = 0; i < A->rmap->n; i++) j = PetscMax(j, dnz[i] + onz[i]);
501:       PetscCall(PetscMalloc1(j, &vals));
502:       for (i = 0; i < j; i++) vals[i] = 1.0;

504:       PetscCall(MatCreate(comm, &A2));
505:       PetscCall(MatSetType(A2, MATMPIAIJ));
506:       PetscCall(MatSetSizes(A2, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
507:       PetscCall(MatMPIAIJSetPreallocation(A2, 0, dnz, 0, onz));
508:       PetscCall(MatSetOption(A2, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
509:       for (i = 0, j = 0; i < A2->rmap->n; i++) {
510:         PetscInt        row = i + A2->rmap->rstart, s = dnz[i] + onz[i];
511:         const PetscInt *idxs;

513:         if (s) {
514:           PetscCall(ISGetIndices(workis[j], &idxs));
515:           PetscCall(MatSetValues(A2, 1, &row, s, idxs, vals, INSERT_VALUES));
516:           PetscCall(ISRestoreIndices(workis[j], &idxs));
517:           j++;
518:         }
519:       }
520:       PetscCheck(j == cnt, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected local count %" PetscInt_FMT " != %" PetscInt_FMT, j, cnt);
521:       PetscCall(PetscFree(vals));
522:       PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY));
523:       PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY));
524:       PetscCall(MatTranspose(A2, MAT_INPLACE_MATRIX, &A2));

526:       /* extract submatrix corresponding to the coupling "owned separators" x "isolated separators" */
527:       for (i = 0, j = 0; i < nl; i++)
528:         if (lndmapi[i] >= 0) lndmapi[j++] = lndmapi[i];
529:       PetscCall(ISCreateGeneral(comm, j, lndmapi, PETSC_USE_POINTER, &is));
530:       PetscCall(MatMPIAIJGetLocalMatCondensed(A2, MAT_INITIAL_MATRIX, &is, NULL, &A3));
531:       PetscCall(ISDestroy(&is));
532:       PetscCall(MatDestroy(&A2));

534:       /* extend local to global map to include connected isolated separators */
535:       PetscCall(PetscObjectQuery((PetscObject)A3, "_petsc_GetLocalMatCondensed_iscol", (PetscObject *)&is));
536:       PetscCheck(is, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing column map");
537:       PetscCall(ISLocalToGlobalMappingCreateIS(is, &ll2g));
538:       PetscCall(MatGetRowIJ(A3, 0, PETSC_FALSE, PETSC_FALSE, &i, &ii, &jj, &flg));
539:       PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure");
540:       PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ii[i], jj, PETSC_COPY_VALUES, &is));
541:       PetscCall(MatRestoreRowIJ(A3, 0, PETSC_FALSE, PETSC_FALSE, &i, &ii, &jj, &flg));
542:       PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure");
543:       PetscCall(ISLocalToGlobalMappingApplyIS(ll2g, is, &is2));
544:       PetscCall(ISDestroy(&is));
545:       PetscCall(ISLocalToGlobalMappingDestroy(&ll2g));

547:       /* add new nodes to the local-to-global map */
548:       PetscCall(ISLocalToGlobalMappingDestroy(l2g));
549:       PetscCall(ISExpand(ndsub, is2, &is));
550:       PetscCall(ISDestroy(&is2));
551:       PetscCall(ISLocalToGlobalMappingCreateIS(is, l2g));
552:       PetscCall(ISDestroy(&is));

554:       PetscCall(MatDestroy(&A3));
555:       PetscCall(PetscFree(lndmapi));
556:       MatPreallocateEnd(dnz, onz);
557:       for (i = 0; i < gcnt; i++) PetscCall(ISDestroy(&workis[i]));
558:       PetscCall(PetscFree(workis));
559:     }
560:     PetscCall(ISRestoreIndices(ndmap, &ndmapi));
561:     PetscCall(PetscSFDestroy(&sf));
562:     PetscCall(PetscFree(ndmapc));
563:     PetscCall(ISDestroy(&ndmap));
564:     PetscCall(ISDestroy(&ndsub));
565:     PetscCall(ISLocalToGlobalMappingSetBlockSize(*l2g, bs));
566:     PetscCall(ISLocalToGlobalMappingViewFromOptions(*l2g, NULL, "-mat_is_nd_l2g_view"));
567:     break;
568:   case MAT_IS_DISASSEMBLE_L2G_NATURAL:
569:     PetscCall(PetscObjectQuery((PetscObject)A, "__PETSc_dm", (PetscObject *)&dm));
570:     if (dm) { /* if a matrix comes from a DM, most likely we can use the l2gmap if any */
571:       PetscCall(MatGetLocalToGlobalMapping(A, l2g, NULL));
572:       PetscCall(PetscObjectReference((PetscObject)*l2g));
573:       if (*l2g) PetscFunctionReturn(PETSC_SUCCESS);
574:     }
575:     if (ismpiaij) {
576:       PetscCall(MatMPIAIJGetSeqAIJ(A, &Ad, &Ao, &garray));
577:     } else if (ismpibaij) {
578:       PetscCall(MatMPIBAIJGetSeqBAIJ(A, &Ad, &Ao, &garray));
579:     } else SETERRQ(comm, PETSC_ERR_SUP, "Type %s", ((PetscObject)A)->type_name);
580:     if (A->rmap->n) {
581:       PetscInt dc, oc, stc, *aux;

583:       PetscCall(MatGetLocalSize(Ad, NULL, &dc));
584:       PetscCall(MatGetLocalSize(Ao, NULL, &oc));
585:       PetscCheck(!oc || garray, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "garray not present");
586:       PetscCall(MatGetOwnershipRangeColumn(A, &stc, NULL));
587:       PetscCall(PetscMalloc1((dc + oc) / bs, &aux));
588:       for (i = 0; i < dc / bs; i++) aux[i] = i + stc / bs;
589:       for (i = 0; i < oc / bs; i++) aux[i + dc / bs] = (ismpiaij ? garray[i * bs] / bs : garray[i]);
590:       PetscCall(ISCreateBlock(comm, bs, (dc + oc) / bs, aux, PETSC_OWN_POINTER, &is));
591:     } else {
592:       PetscCall(ISCreateBlock(comm, 1, 0, NULL, PETSC_OWN_POINTER, &is));
593:     }
594:     PetscCall(ISLocalToGlobalMappingCreateIS(is, l2g));
595:     PetscCall(ISDestroy(&is));
596:     break;
597:   default:
598:     SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unsupported l2g disassembling type %d", mode);
599:   }
600:   PetscFunctionReturn(PETSC_SUCCESS);
601: }

603: PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat A, MatType type, MatReuse reuse, Mat *newmat)
604: {
605:   Mat                    lA, Ad, Ao, B = NULL;
606:   ISLocalToGlobalMapping rl2g, cl2g;
607:   IS                     is;
608:   MPI_Comm               comm;
609:   void                  *ptrs[2];
610:   const char            *names[2] = {"_convert_csr_aux", "_convert_csr_data"};
611:   const PetscInt        *garray;
612:   PetscScalar           *dd, *od, *aa, *data;
613:   const PetscInt        *di, *dj, *oi, *oj;
614:   const PetscInt        *odi, *odj, *ooi, *ooj;
615:   PetscInt              *aux, *ii, *jj;
616:   PetscInt               bs, lc, dr, dc, oc, str, stc, nnz, i, jd, jo, cum;
617:   PetscBool              flg, ismpiaij, ismpibaij, was_inplace = PETSC_FALSE;
618:   PetscMPIInt            size;

620:   PetscFunctionBegin;
621:   PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
622:   PetscCallMPI(MPI_Comm_size(comm, &size));
623:   if (size == 1) {
624:     PetscCall(MatConvert_SeqXAIJ_IS(A, type, reuse, newmat));
625:     PetscFunctionReturn(PETSC_SUCCESS);
626:   }
627:   if (reuse != MAT_REUSE_MATRIX && A->cmap->N == A->rmap->N) {
628:     PetscCall(MatMPIXAIJComputeLocalToGlobalMapping_Private(A, &rl2g));
629:     PetscCall(MatCreate(comm, &B));
630:     PetscCall(MatSetType(B, MATIS));
631:     PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
632:     PetscCall(MatSetLocalToGlobalMapping(B, rl2g, rl2g));
633:     PetscCall(MatGetBlockSize(A, &bs));
634:     PetscCall(MatSetBlockSize(B, bs));
635:     PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
636:     if (reuse == MAT_INPLACE_MATRIX) was_inplace = PETSC_TRUE;
637:     reuse = MAT_REUSE_MATRIX;
638:   }
639:   if (reuse == MAT_REUSE_MATRIX) {
640:     Mat            *newlA, lA;
641:     IS              rows, cols;
642:     const PetscInt *ridx, *cidx;
643:     PetscInt        rbs, cbs, nr, nc;

645:     if (!B) B = *newmat;
646:     PetscCall(MatISGetLocalToGlobalMapping(B, &rl2g, &cl2g));
647:     PetscCall(ISLocalToGlobalMappingGetBlockIndices(rl2g, &ridx));
648:     PetscCall(ISLocalToGlobalMappingGetBlockIndices(cl2g, &cidx));
649:     PetscCall(ISLocalToGlobalMappingGetSize(rl2g, &nr));
650:     PetscCall(ISLocalToGlobalMappingGetSize(cl2g, &nc));
651:     PetscCall(ISLocalToGlobalMappingGetBlockSize(rl2g, &rbs));
652:     PetscCall(ISLocalToGlobalMappingGetBlockSize(cl2g, &cbs));
653:     PetscCall(ISCreateBlock(comm, rbs, nr / rbs, ridx, PETSC_USE_POINTER, &rows));
654:     if (rl2g != cl2g) {
655:       PetscCall(ISCreateBlock(comm, cbs, nc / cbs, cidx, PETSC_USE_POINTER, &cols));
656:     } else {
657:       PetscCall(PetscObjectReference((PetscObject)rows));
658:       cols = rows;
659:     }
660:     PetscCall(MatISGetLocalMat(B, &lA));
661:     PetscCall(MatCreateSubMatrices(A, 1, &rows, &cols, MAT_INITIAL_MATRIX, &newlA));
662:     PetscCall(MatConvert(newlA[0], MATSEQAIJ, MAT_INPLACE_MATRIX, &newlA[0]));
663:     PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(rl2g, &ridx));
664:     PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(cl2g, &cidx));
665:     PetscCall(ISDestroy(&rows));
666:     PetscCall(ISDestroy(&cols));
667:     if (!lA->preallocated) { /* first time */
668:       PetscCall(MatDuplicate(newlA[0], MAT_COPY_VALUES, &lA));
669:       PetscCall(MatISSetLocalMat(B, lA));
670:       PetscCall(PetscObjectDereference((PetscObject)lA));
671:     }
672:     PetscCall(MatCopy(newlA[0], lA, SAME_NONZERO_PATTERN));
673:     PetscCall(MatDestroySubMatrices(1, &newlA));
674:     PetscCall(MatISScaleDisassembling_Private(B));
675:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
676:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
677:     if (was_inplace) PetscCall(MatHeaderReplace(A, &B));
678:     else *newmat = B;
679:     PetscFunctionReturn(PETSC_SUCCESS);
680:   }
681:   /* rectangular case, just compress out the column space */
682:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij));
683:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIBAIJ, &ismpibaij));
684:   if (ismpiaij) {
685:     bs = 1;
686:     PetscCall(MatMPIAIJGetSeqAIJ(A, &Ad, &Ao, &garray));
687:   } else if (ismpibaij) {
688:     PetscCall(MatGetBlockSize(A, &bs));
689:     PetscCall(MatMPIBAIJGetSeqBAIJ(A, &Ad, &Ao, &garray));
690:     PetscCall(MatConvert(Ad, MATSEQAIJ, MAT_INITIAL_MATRIX, &Ad));
691:     PetscCall(MatConvert(Ao, MATSEQAIJ, MAT_INITIAL_MATRIX, &Ao));
692:   } else SETERRQ(comm, PETSC_ERR_SUP, "Type %s", ((PetscObject)A)->type_name);
693:   PetscCall(MatSeqAIJGetArray(Ad, &dd));
694:   PetscCall(MatSeqAIJGetArray(Ao, &od));

696:   /* access relevant information from MPIAIJ */
697:   PetscCall(MatGetOwnershipRange(A, &str, NULL));
698:   PetscCall(MatGetOwnershipRangeColumn(A, &stc, NULL));
699:   PetscCall(MatGetLocalSize(A, &dr, &dc));
700:   PetscCall(MatGetLocalSize(Ao, NULL, &oc));
701:   PetscCheck(!oc || garray, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "garray not present");

703:   PetscCall(MatGetRowIJ(Ad, 0, PETSC_FALSE, PETSC_FALSE, &i, &di, &dj, &flg));
704:   PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure");
705:   PetscCall(MatGetRowIJ(Ao, 0, PETSC_FALSE, PETSC_FALSE, &i, &oi, &oj, &flg));
706:   PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure");
707:   nnz = di[dr] + oi[dr];
708:   /* store original pointers to be restored later */
709:   odi = di;
710:   odj = dj;
711:   ooi = oi;
712:   ooj = oj;

714:   /* generate l2g maps for rows and cols */
715:   PetscCall(ISCreateStride(comm, dr / bs, str / bs, 1, &is));
716:   if (bs > 1) {
717:     IS is2;

719:     PetscCall(ISGetLocalSize(is, &i));
720:     PetscCall(ISGetIndices(is, (const PetscInt **)&aux));
721:     PetscCall(ISCreateBlock(comm, bs, i, aux, PETSC_COPY_VALUES, &is2));
722:     PetscCall(ISRestoreIndices(is, (const PetscInt **)&aux));
723:     PetscCall(ISDestroy(&is));
724:     is = is2;
725:   }
726:   PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g));
727:   PetscCall(ISDestroy(&is));
728:   if (dr) {
729:     PetscCall(PetscMalloc1((dc + oc) / bs, &aux));
730:     for (i = 0; i < dc / bs; i++) aux[i] = i + stc / bs;
731:     for (i = 0; i < oc / bs; i++) aux[i + dc / bs] = garray[i];
732:     PetscCall(ISCreateBlock(comm, bs, (dc + oc) / bs, aux, PETSC_OWN_POINTER, &is));
733:     lc = dc + oc;
734:   } else {
735:     PetscCall(ISCreateBlock(comm, bs, 0, NULL, PETSC_OWN_POINTER, &is));
736:     lc = 0;
737:   }
738:   PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g));
739:   PetscCall(ISDestroy(&is));

741:   /* create MATIS object */
742:   PetscCall(MatCreate(comm, &B));
743:   PetscCall(MatSetSizes(B, dr, dc, PETSC_DECIDE, PETSC_DECIDE));
744:   PetscCall(MatSetType(B, MATIS));
745:   PetscCall(MatSetBlockSize(B, bs));
746:   PetscCall(MatSetLocalToGlobalMapping(B, rl2g, cl2g));
747:   PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
748:   PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));

750:   /* merge local matrices */
751:   PetscCall(PetscMalloc1(nnz + dr + 1, &aux));
752:   PetscCall(PetscMalloc1(nnz, &data));
753:   ii  = aux;
754:   jj  = aux + dr + 1;
755:   aa  = data;
756:   *ii = *(di++) + *(oi++);
757:   for (jd = 0, jo = 0, cum = 0; *ii < nnz; cum++) {
758:     for (; jd < *di; jd++) {
759:       *jj++ = *dj++;
760:       *aa++ = *dd++;
761:     }
762:     for (; jo < *oi; jo++) {
763:       *jj++ = *oj++ + dc;
764:       *aa++ = *od++;
765:     }
766:     *(++ii) = *(di++) + *(oi++);
767:   }
768:   for (; cum < dr; cum++) *(++ii) = nnz;

770:   PetscCall(MatRestoreRowIJ(Ad, 0, PETSC_FALSE, PETSC_FALSE, &i, &odi, &odj, &flg));
771:   PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot restore IJ structure");
772:   PetscCall(MatRestoreRowIJ(Ao, 0, PETSC_FALSE, PETSC_FALSE, &i, &ooi, &ooj, &flg));
773:   PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot restore IJ structure");
774:   PetscCall(MatSeqAIJRestoreArray(Ad, &dd));
775:   PetscCall(MatSeqAIJRestoreArray(Ao, &od));

777:   ii = aux;
778:   jj = aux + dr + 1;
779:   aa = data;
780:   PetscCall(MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, dr, lc, ii, jj, aa, &lA));

782:   /* create containers to destroy the data */
783:   ptrs[0] = aux;
784:   ptrs[1] = data;
785:   for (i = 0; i < 2; i++) {
786:     PetscContainer c;

788:     PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &c));
789:     PetscCall(PetscContainerSetPointer(c, ptrs[i]));
790:     PetscCall(PetscContainerSetUserDestroy(c, PetscContainerUserDestroyDefault));
791:     PetscCall(PetscObjectCompose((PetscObject)lA, names[i], (PetscObject)c));
792:     PetscCall(PetscContainerDestroy(&c));
793:   }
794:   if (ismpibaij) { /* destroy converted local matrices */
795:     PetscCall(MatDestroy(&Ad));
796:     PetscCall(MatDestroy(&Ao));
797:   }

799:   /* finalize matrix */
800:   PetscCall(MatISSetLocalMat(B, lA));
801:   PetscCall(MatDestroy(&lA));
802:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
803:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
804:   if (reuse == MAT_INPLACE_MATRIX) {
805:     PetscCall(MatHeaderReplace(A, &B));
806:   } else *newmat = B;
807:   PetscFunctionReturn(PETSC_SUCCESS);
808: }

810: PETSC_INTERN PetscErrorCode MatConvert_Nest_IS(Mat A, MatType type, MatReuse reuse, Mat *newmat)
811: {
812:   Mat                  **nest, *snest, **rnest, lA, B;
813:   IS                    *iscol, *isrow, *islrow, *islcol;
814:   ISLocalToGlobalMapping rl2g, cl2g;
815:   MPI_Comm               comm;
816:   PetscInt              *lr, *lc, *l2gidxs;
817:   PetscInt               i, j, nr, nc, rbs, cbs;
818:   PetscBool              convert, lreuse, *istrans;
819:   PetscBool3             allow_repeated = PETSC_BOOL3_UNKNOWN;

821:   PetscFunctionBegin;
822:   PetscCall(MatNestGetSubMats(A, &nr, &nc, &nest));
823:   lreuse = PETSC_FALSE;
824:   rnest  = NULL;
825:   if (reuse == MAT_REUSE_MATRIX) {
826:     PetscBool ismatis, isnest;

828:     PetscCall(PetscObjectTypeCompare((PetscObject)*newmat, MATIS, &ismatis));
829:     PetscCheck(ismatis, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_USER, "Cannot reuse matrix of type %s", ((PetscObject)*newmat)->type_name);
830:     PetscCall(MatISGetLocalMat(*newmat, &lA));
831:     PetscCall(PetscObjectTypeCompare((PetscObject)lA, MATNEST, &isnest));
832:     if (isnest) {
833:       PetscCall(MatNestGetSubMats(lA, &i, &j, &rnest));
834:       lreuse = (PetscBool)(i == nr && j == nc);
835:       if (!lreuse) rnest = NULL;
836:     }
837:   }
838:   PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
839:   PetscCall(PetscCalloc2(nr, &lr, nc, &lc));
840:   PetscCall(PetscCalloc6(nr, &isrow, nc, &iscol, nr, &islrow, nc, &islcol, nr * nc, &snest, nr * nc, &istrans));
841:   PetscCall(MatNestGetISs(A, isrow, iscol));
842:   for (i = 0; i < nr; i++) {
843:     for (j = 0; j < nc; j++) {
844:       PetscBool ismatis, sallow;
845:       PetscInt  l1, l2, lb1, lb2, ij = i * nc + j;

847:       /* Null matrix pointers are allowed in MATNEST */
848:       if (!nest[i][j]) continue;

850:       /* Nested matrices should be of type MATIS */
851:       PetscCall(PetscObjectTypeCompare((PetscObject)nest[i][j], MATTRANSPOSEVIRTUAL, &istrans[ij]));
852:       if (istrans[ij]) {
853:         Mat T, lT;
854:         PetscCall(MatTransposeGetMat(nest[i][j], &T));
855:         PetscCall(PetscObjectTypeCompare((PetscObject)T, MATIS, &ismatis));
856:         PetscCheck(ismatis, comm, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") (transposed) is not of type MATIS", i, j);
857:         PetscCall(MatISGetAllowRepeated(T, &sallow));
858:         PetscCall(MatISGetLocalMat(T, &lT));
859:         PetscCall(MatCreateTranspose(lT, &snest[ij]));
860:       } else {
861:         PetscCall(PetscObjectTypeCompare((PetscObject)nest[i][j], MATIS, &ismatis));
862:         PetscCheck(ismatis, comm, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") is not of type MATIS", i, j);
863:         PetscCall(MatISGetAllowRepeated(nest[i][j], &sallow));
864:         PetscCall(MatISGetLocalMat(nest[i][j], &snest[ij]));
865:       }
866:       if (allow_repeated == PETSC_BOOL3_UNKNOWN) allow_repeated = PetscBoolToBool3(sallow);
867:       PetscCheck(sallow == PetscBool3ToBool(allow_repeated), comm, PETSC_ERR_SUP, "Cannot mix repeated and non repeated maps");

869:       /* Check compatibility of local sizes */
870:       PetscCall(MatGetSize(snest[ij], &l1, &l2));
871:       PetscCall(MatGetBlockSizes(snest[ij], &lb1, &lb2));
872:       if (!l1 || !l2) continue;
873:       PetscCheck(!lr[i] || l1 == lr[i], PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid local size %" PetscInt_FMT " != %" PetscInt_FMT, i, j, lr[i], l1);
874:       PetscCheck(!lc[j] || l2 == lc[j], PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid local size %" PetscInt_FMT " != %" PetscInt_FMT, i, j, lc[j], l2);
875:       lr[i] = l1;
876:       lc[j] = l2;

878:       /* check compatibility for local matrix reusage */
879:       if (rnest && !rnest[i][j] != !snest[ij]) lreuse = PETSC_FALSE;
880:     }
881:   }

883:   if (PetscDefined(USE_DEBUG)) {
884:     /* Check compatibility of l2g maps for rows */
885:     for (i = 0; i < nr; i++) {
886:       rl2g = NULL;
887:       for (j = 0; j < nc; j++) {
888:         PetscInt n1, n2;

890:         if (!nest[i][j]) continue;
891:         if (istrans[i * nc + j]) {
892:           Mat T;

894:           PetscCall(MatTransposeGetMat(nest[i][j], &T));
895:           PetscCall(MatISGetLocalToGlobalMapping(T, NULL, &cl2g));
896:         } else {
897:           PetscCall(MatISGetLocalToGlobalMapping(nest[i][j], &cl2g, NULL));
898:         }
899:         PetscCall(ISLocalToGlobalMappingGetSize(cl2g, &n1));
900:         if (!n1) continue;
901:         if (!rl2g) {
902:           rl2g = cl2g;
903:         } else {
904:           const PetscInt *idxs1, *idxs2;
905:           PetscBool       same;

907:           PetscCall(ISLocalToGlobalMappingGetSize(rl2g, &n2));
908:           PetscCheck(n1 == n2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid row l2gmap size %" PetscInt_FMT " != %" PetscInt_FMT, i, j, n1, n2);
909:           PetscCall(ISLocalToGlobalMappingGetIndices(cl2g, &idxs1));
910:           PetscCall(ISLocalToGlobalMappingGetIndices(rl2g, &idxs2));
911:           PetscCall(PetscArraycmp(idxs1, idxs2, n1, &same));
912:           PetscCall(ISLocalToGlobalMappingRestoreIndices(cl2g, &idxs1));
913:           PetscCall(ISLocalToGlobalMappingRestoreIndices(rl2g, &idxs2));
914:           PetscCheck(same, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid row l2gmap", i, j);
915:         }
916:       }
917:     }
918:     /* Check compatibility of l2g maps for columns */
919:     for (i = 0; i < nc; i++) {
920:       rl2g = NULL;
921:       for (j = 0; j < nr; j++) {
922:         PetscInt n1, n2;

924:         if (!nest[j][i]) continue;
925:         if (istrans[j * nc + i]) {
926:           Mat T;

928:           PetscCall(MatTransposeGetMat(nest[j][i], &T));
929:           PetscCall(MatISGetLocalToGlobalMapping(T, &cl2g, NULL));
930:         } else {
931:           PetscCall(MatISGetLocalToGlobalMapping(nest[j][i], NULL, &cl2g));
932:         }
933:         PetscCall(ISLocalToGlobalMappingGetSize(cl2g, &n1));
934:         if (!n1) continue;
935:         if (!rl2g) {
936:           rl2g = cl2g;
937:         } else {
938:           const PetscInt *idxs1, *idxs2;
939:           PetscBool       same;

941:           PetscCall(ISLocalToGlobalMappingGetSize(rl2g, &n2));
942:           PetscCheck(n1 == n2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid column l2gmap size %" PetscInt_FMT " != %" PetscInt_FMT, j, i, n1, n2);
943:           PetscCall(ISLocalToGlobalMappingGetIndices(cl2g, &idxs1));
944:           PetscCall(ISLocalToGlobalMappingGetIndices(rl2g, &idxs2));
945:           PetscCall(PetscArraycmp(idxs1, idxs2, n1, &same));
946:           PetscCall(ISLocalToGlobalMappingRestoreIndices(cl2g, &idxs1));
947:           PetscCall(ISLocalToGlobalMappingRestoreIndices(rl2g, &idxs2));
948:           PetscCheck(same, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid column l2gmap", j, i);
949:         }
950:       }
951:     }
952:   }

954:   B = NULL;
955:   if (reuse != MAT_REUSE_MATRIX) {
956:     PetscInt stl;

958:     /* Create l2g map for the rows of the new matrix and index sets for the local MATNEST */
959:     for (i = 0, stl = 0; i < nr; i++) stl += lr[i];
960:     PetscCall(PetscMalloc1(stl, &l2gidxs));
961:     for (i = 0, stl = 0; i < nr; i++) {
962:       Mat             usedmat;
963:       Mat_IS         *matis;
964:       const PetscInt *idxs;

966:       /* local IS for local NEST */
967:       PetscCall(ISCreateStride(PETSC_COMM_SELF, lr[i], stl, 1, &islrow[i]));

969:       /* l2gmap */
970:       j       = 0;
971:       usedmat = nest[i][j];
972:       while (!usedmat && j < nc - 1) usedmat = nest[i][++j];
973:       PetscCheck(usedmat, comm, PETSC_ERR_SUP, "Cannot find valid row mat");

975:       if (istrans[i * nc + j]) {
976:         Mat T;
977:         PetscCall(MatTransposeGetMat(usedmat, &T));
978:         usedmat = T;
979:       }
980:       matis = (Mat_IS *)usedmat->data;
981:       PetscCall(ISGetIndices(isrow[i], &idxs));
982:       if (istrans[i * nc + j]) {
983:         PetscCall(PetscSFBcastBegin(matis->csf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
984:         PetscCall(PetscSFBcastEnd(matis->csf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
985:       } else {
986:         PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
987:         PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
988:       }
989:       PetscCall(ISRestoreIndices(isrow[i], &idxs));
990:       stl += lr[i];
991:     }
992:     PetscCall(ISLocalToGlobalMappingCreate(comm, 1, stl, l2gidxs, PETSC_OWN_POINTER, &rl2g));

994:     /* Create l2g map for columns of the new matrix and index sets for the local MATNEST */
995:     for (i = 0, stl = 0; i < nc; i++) stl += lc[i];
996:     PetscCall(PetscMalloc1(stl, &l2gidxs));
997:     for (i = 0, stl = 0; i < nc; i++) {
998:       Mat             usedmat;
999:       Mat_IS         *matis;
1000:       const PetscInt *idxs;

1002:       /* local IS for local NEST */
1003:       PetscCall(ISCreateStride(PETSC_COMM_SELF, lc[i], stl, 1, &islcol[i]));

1005:       /* l2gmap */
1006:       j       = 0;
1007:       usedmat = nest[j][i];
1008:       while (!usedmat && j < nr - 1) usedmat = nest[++j][i];
1009:       PetscCheck(usedmat, comm, PETSC_ERR_SUP, "Cannot find valid column mat");
1010:       if (istrans[j * nc + i]) {
1011:         Mat T;
1012:         PetscCall(MatTransposeGetMat(usedmat, &T));
1013:         usedmat = T;
1014:       }
1015:       matis = (Mat_IS *)usedmat->data;
1016:       PetscCall(ISGetIndices(iscol[i], &idxs));
1017:       if (istrans[j * nc + i]) {
1018:         PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
1019:         PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
1020:       } else {
1021:         PetscCall(PetscSFBcastBegin(matis->csf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
1022:         PetscCall(PetscSFBcastEnd(matis->csf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
1023:       }
1024:       PetscCall(ISRestoreIndices(iscol[i], &idxs));
1025:       stl += lc[i];
1026:     }
1027:     PetscCall(ISLocalToGlobalMappingCreate(comm, 1, stl, l2gidxs, PETSC_OWN_POINTER, &cl2g));

1029:     /* Create MATIS */
1030:     PetscCall(MatCreate(comm, &B));
1031:     PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
1032:     PetscCall(MatGetBlockSizes(A, &rbs, &cbs));
1033:     PetscCall(MatSetBlockSizes(B, rbs, cbs));
1034:     PetscCall(MatSetType(B, MATIS));
1035:     PetscCall(MatISSetLocalMatType(B, MATNEST));
1036:     PetscCall(MatISSetAllowRepeated(B, PetscBool3ToBool(allow_repeated)));
1037:     { /* hack : avoid setup of scatters */
1038:       Mat_IS *matis     = (Mat_IS *)B->data;
1039:       matis->islocalref = PETSC_TRUE;
1040:     }
1041:     PetscCall(MatSetLocalToGlobalMapping(B, rl2g, cl2g));
1042:     PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
1043:     PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
1044:     PetscCall(MatCreateNest(PETSC_COMM_SELF, nr, islrow, nc, islcol, snest, &lA));
1045:     PetscCall(MatNestSetVecType(lA, VECNEST));
1046:     for (i = 0; i < nr * nc; i++) {
1047:       if (istrans[i]) PetscCall(MatDestroy(&snest[i]));
1048:     }
1049:     PetscCall(MatISSetLocalMat(B, lA));
1050:     PetscCall(MatDestroy(&lA));
1051:     { /* hack : setup of scatters done here */
1052:       Mat_IS *matis = (Mat_IS *)B->data;

1054:       matis->islocalref = PETSC_FALSE;
1055:       PetscCall(MatISSetUpScatters_Private(B));
1056:     }
1057:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
1058:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
1059:     if (reuse == MAT_INPLACE_MATRIX) {
1060:       PetscCall(MatHeaderReplace(A, &B));
1061:     } else {
1062:       *newmat = B;
1063:     }
1064:   } else {
1065:     if (lreuse) {
1066:       PetscCall(MatISGetLocalMat(*newmat, &lA));
1067:       for (i = 0; i < nr; i++) {
1068:         for (j = 0; j < nc; j++) {
1069:           if (snest[i * nc + j]) {
1070:             PetscCall(MatNestSetSubMat(lA, i, j, snest[i * nc + j]));
1071:             if (istrans[i * nc + j]) PetscCall(MatDestroy(&snest[i * nc + j]));
1072:           }
1073:         }
1074:       }
1075:     } else {
1076:       PetscInt stl;
1077:       for (i = 0, stl = 0; i < nr; i++) {
1078:         PetscCall(ISCreateStride(PETSC_COMM_SELF, lr[i], stl, 1, &islrow[i]));
1079:         stl += lr[i];
1080:       }
1081:       for (i = 0, stl = 0; i < nc; i++) {
1082:         PetscCall(ISCreateStride(PETSC_COMM_SELF, lc[i], stl, 1, &islcol[i]));
1083:         stl += lc[i];
1084:       }
1085:       PetscCall(MatCreateNest(PETSC_COMM_SELF, nr, islrow, nc, islcol, snest, &lA));
1086:       for (i = 0; i < nr * nc; i++) {
1087:         if (istrans[i]) PetscCall(MatDestroy(&snest[i]));
1088:       }
1089:       PetscCall(MatISSetLocalMat(*newmat, lA));
1090:       PetscCall(MatDestroy(&lA));
1091:     }
1092:     PetscCall(MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY));
1093:     PetscCall(MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY));
1094:   }

1096:   /* Create local matrix in MATNEST format */
1097:   convert = PETSC_FALSE;
1098:   PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)A)->prefix, "-mat_is_convert_local_nest", &convert, NULL));
1099:   if (convert) {
1100:     Mat              M;
1101:     MatISLocalFields lf;
1102:     PetscContainer   c;

1104:     PetscCall(MatISGetLocalMat(*newmat, &lA));
1105:     PetscCall(MatConvert(lA, MATAIJ, MAT_INITIAL_MATRIX, &M));
1106:     PetscCall(MatISSetLocalMat(*newmat, M));
1107:     PetscCall(MatDestroy(&M));

1109:     /* attach local fields to the matrix */
1110:     PetscCall(PetscNew(&lf));
1111:     PetscCall(PetscMalloc2(nr, &lf->rf, nc, &lf->cf));
1112:     for (i = 0; i < nr; i++) {
1113:       PetscInt n, st;

1115:       PetscCall(ISGetLocalSize(islrow[i], &n));
1116:       PetscCall(ISStrideGetInfo(islrow[i], &st, NULL));
1117:       PetscCall(ISCreateStride(comm, n, st, 1, &lf->rf[i]));
1118:     }
1119:     for (i = 0; i < nc; i++) {
1120:       PetscInt n, st;

1122:       PetscCall(ISGetLocalSize(islcol[i], &n));
1123:       PetscCall(ISStrideGetInfo(islcol[i], &st, NULL));
1124:       PetscCall(ISCreateStride(comm, n, st, 1, &lf->cf[i]));
1125:     }
1126:     lf->nr = nr;
1127:     lf->nc = nc;
1128:     PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)*newmat), &c));
1129:     PetscCall(PetscContainerSetPointer(c, lf));
1130:     PetscCall(PetscContainerSetUserDestroy(c, MatISContainerDestroyFields_Private));
1131:     PetscCall(PetscObjectCompose((PetscObject)*newmat, "_convert_nest_lfields", (PetscObject)c));
1132:     PetscCall(PetscContainerDestroy(&c));
1133:   }

1135:   /* Free workspace */
1136:   for (i = 0; i < nr; i++) PetscCall(ISDestroy(&islrow[i]));
1137:   for (i = 0; i < nc; i++) PetscCall(ISDestroy(&islcol[i]));
1138:   PetscCall(PetscFree6(isrow, iscol, islrow, islcol, snest, istrans));
1139:   PetscCall(PetscFree2(lr, lc));
1140:   PetscFunctionReturn(PETSC_SUCCESS);
1141: }

1143: static PetscErrorCode MatDiagonalScale_IS(Mat A, Vec l, Vec r)
1144: {
1145:   Mat_IS            *matis = (Mat_IS *)A->data;
1146:   Vec                ll, rr;
1147:   const PetscScalar *Y, *X;
1148:   PetscScalar       *x, *y;

1150:   PetscFunctionBegin;
1151:   if (l) {
1152:     ll = matis->y;
1153:     PetscCall(VecGetArrayRead(l, &Y));
1154:     PetscCall(VecGetArray(ll, &y));
1155:     PetscCall(PetscSFBcastBegin(matis->sf, MPIU_SCALAR, Y, y, MPI_REPLACE));
1156:   } else {
1157:     ll = NULL;
1158:   }
1159:   if (r) {
1160:     rr = matis->x;
1161:     PetscCall(VecGetArrayRead(r, &X));
1162:     PetscCall(VecGetArray(rr, &x));
1163:     PetscCall(PetscSFBcastBegin(matis->csf, MPIU_SCALAR, X, x, MPI_REPLACE));
1164:   } else {
1165:     rr = NULL;
1166:   }
1167:   if (ll) {
1168:     PetscCall(PetscSFBcastEnd(matis->sf, MPIU_SCALAR, Y, y, MPI_REPLACE));
1169:     PetscCall(VecRestoreArrayRead(l, &Y));
1170:     PetscCall(VecRestoreArray(ll, &y));
1171:   }
1172:   if (rr) {
1173:     PetscCall(PetscSFBcastEnd(matis->csf, MPIU_SCALAR, X, x, MPI_REPLACE));
1174:     PetscCall(VecRestoreArrayRead(r, &X));
1175:     PetscCall(VecRestoreArray(rr, &x));
1176:   }
1177:   PetscCall(MatDiagonalScale(matis->A, ll, rr));
1178:   PetscFunctionReturn(PETSC_SUCCESS);
1179: }

1181: static PetscErrorCode MatGetInfo_IS(Mat A, MatInfoType flag, MatInfo *ginfo)
1182: {
1183:   Mat_IS        *matis = (Mat_IS *)A->data;
1184:   MatInfo        info;
1185:   PetscLogDouble isend[6], irecv[6];
1186:   PetscInt       bs;

1188:   PetscFunctionBegin;
1189:   PetscCall(MatGetBlockSize(A, &bs));
1190:   if (matis->A->ops->getinfo) {
1191:     PetscCall(MatGetInfo(matis->A, MAT_LOCAL, &info));
1192:     isend[0] = info.nz_used;
1193:     isend[1] = info.nz_allocated;
1194:     isend[2] = info.nz_unneeded;
1195:     isend[3] = info.memory;
1196:     isend[4] = info.mallocs;
1197:   } else {
1198:     isend[0] = 0.;
1199:     isend[1] = 0.;
1200:     isend[2] = 0.;
1201:     isend[3] = 0.;
1202:     isend[4] = 0.;
1203:   }
1204:   isend[5] = matis->A->num_ass;
1205:   if (flag == MAT_LOCAL) {
1206:     ginfo->nz_used      = isend[0];
1207:     ginfo->nz_allocated = isend[1];
1208:     ginfo->nz_unneeded  = isend[2];
1209:     ginfo->memory       = isend[3];
1210:     ginfo->mallocs      = isend[4];
1211:     ginfo->assemblies   = isend[5];
1212:   } else if (flag == MAT_GLOBAL_MAX) {
1213:     PetscCall(MPIU_Allreduce(isend, irecv, 6, MPIU_PETSCLOGDOUBLE, MPI_MAX, PetscObjectComm((PetscObject)A)));

1215:     ginfo->nz_used      = irecv[0];
1216:     ginfo->nz_allocated = irecv[1];
1217:     ginfo->nz_unneeded  = irecv[2];
1218:     ginfo->memory       = irecv[3];
1219:     ginfo->mallocs      = irecv[4];
1220:     ginfo->assemblies   = irecv[5];
1221:   } else if (flag == MAT_GLOBAL_SUM) {
1222:     PetscCall(MPIU_Allreduce(isend, irecv, 5, MPIU_PETSCLOGDOUBLE, MPI_SUM, PetscObjectComm((PetscObject)A)));

1224:     ginfo->nz_used      = irecv[0];
1225:     ginfo->nz_allocated = irecv[1];
1226:     ginfo->nz_unneeded  = irecv[2];
1227:     ginfo->memory       = irecv[3];
1228:     ginfo->mallocs      = irecv[4];
1229:     ginfo->assemblies   = A->num_ass;
1230:   }
1231:   ginfo->block_size        = bs;
1232:   ginfo->fill_ratio_given  = 0;
1233:   ginfo->fill_ratio_needed = 0;
1234:   ginfo->factor_mallocs    = 0;
1235:   PetscFunctionReturn(PETSC_SUCCESS);
1236: }

1238: static PetscErrorCode MatTranspose_IS(Mat A, MatReuse reuse, Mat *B)
1239: {
1240:   Mat C, lC, lA;

1242:   PetscFunctionBegin;
1243:   if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *B));
1244:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
1245:     ISLocalToGlobalMapping rl2g, cl2g;
1246:     PetscBool              allow_repeated;

1248:     PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
1249:     PetscCall(MatSetSizes(C, A->cmap->n, A->rmap->n, A->cmap->N, A->rmap->N));
1250:     PetscCall(MatSetBlockSizes(C, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs)));
1251:     PetscCall(MatSetType(C, MATIS));
1252:     PetscCall(MatISGetAllowRepeated(A, &allow_repeated));
1253:     PetscCall(MatISSetAllowRepeated(C, allow_repeated));
1254:     PetscCall(MatGetLocalToGlobalMapping(A, &rl2g, &cl2g));
1255:     PetscCall(MatSetLocalToGlobalMapping(C, cl2g, rl2g));
1256:   } else C = *B;

1258:   /* perform local transposition */
1259:   PetscCall(MatISGetLocalMat(A, &lA));
1260:   PetscCall(MatTranspose(lA, MAT_INITIAL_MATRIX, &lC));
1261:   PetscCall(MatSetLocalToGlobalMapping(lC, lA->cmap->mapping, lA->rmap->mapping));
1262:   PetscCall(MatISSetLocalMat(C, lC));
1263:   PetscCall(MatDestroy(&lC));

1265:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
1266:     *B = C;
1267:   } else {
1268:     PetscCall(MatHeaderMerge(A, &C));
1269:   }
1270:   PetscCall(MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY));
1271:   PetscCall(MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY));
1272:   PetscFunctionReturn(PETSC_SUCCESS);
1273: }

1275: static PetscErrorCode MatDiagonalSet_IS(Mat A, Vec D, InsertMode insmode)
1276: {
1277:   Mat_IS *is = (Mat_IS *)A->data;

1279:   PetscFunctionBegin;
1280:   PetscCheck(!is->allow_repeated || insmode == ADD_VALUES, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "INSERT_VALUES with repeated entries not supported");
1281:   if (D) { /* MatShift_IS pass D = NULL */
1282:     PetscCall(VecScatterBegin(is->rctx, D, is->y, INSERT_VALUES, SCATTER_FORWARD));
1283:     PetscCall(VecScatterEnd(is->rctx, D, is->y, INSERT_VALUES, SCATTER_FORWARD));
1284:   }
1285:   PetscCall(VecPointwiseDivide(is->y, is->y, is->counter));
1286:   PetscCall(MatDiagonalSet(is->A, is->y, insmode));
1287:   PetscFunctionReturn(PETSC_SUCCESS);
1288: }

1290: static PetscErrorCode MatShift_IS(Mat A, PetscScalar a)
1291: {
1292:   Mat_IS *is = (Mat_IS *)A->data;

1294:   PetscFunctionBegin;
1295:   PetscCall(VecSet(is->y, a));
1296:   PetscCall(MatDiagonalSet_IS(A, NULL, ADD_VALUES));
1297:   PetscFunctionReturn(PETSC_SUCCESS);
1298: }

1300: static PetscErrorCode MatSetValuesLocal_SubMat_IS(Mat A, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
1301: {
1302:   PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION], cols_l[MATIS_MAX_ENTRIES_INSERTION];

1304:   PetscFunctionBegin;
1305:   PetscCheck(m <= MATIS_MAX_ENTRIES_INSERTION && n <= MATIS_MAX_ENTRIES_INSERTION, PETSC_COMM_SELF, PETSC_ERR_SUP, "Number of row/column indices must be <= %d: they are %" PetscInt_FMT " %" PetscInt_FMT, MATIS_MAX_ENTRIES_INSERTION, m, n);
1306:   PetscCall(ISLocalToGlobalMappingApply(A->rmap->mapping, m, rows, rows_l));
1307:   PetscCall(ISLocalToGlobalMappingApply(A->cmap->mapping, n, cols, cols_l));
1308:   PetscCall(MatSetValuesLocal_IS(A, m, rows_l, n, cols_l, values, addv));
1309:   PetscFunctionReturn(PETSC_SUCCESS);
1310: }

1312: static PetscErrorCode MatSetValuesBlockedLocal_SubMat_IS(Mat A, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
1313: {
1314:   PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION], cols_l[MATIS_MAX_ENTRIES_INSERTION];

1316:   PetscFunctionBegin;
1317:   PetscCheck(m <= MATIS_MAX_ENTRIES_INSERTION && n <= MATIS_MAX_ENTRIES_INSERTION, PETSC_COMM_SELF, PETSC_ERR_SUP, "Number of row/column block indices must be <= %d: they are %" PetscInt_FMT " %" PetscInt_FMT, MATIS_MAX_ENTRIES_INSERTION, m, n);
1318:   PetscCall(ISLocalToGlobalMappingApplyBlock(A->rmap->mapping, m, rows, rows_l));
1319:   PetscCall(ISLocalToGlobalMappingApplyBlock(A->cmap->mapping, n, cols, cols_l));
1320:   PetscCall(MatSetValuesBlockedLocal_IS(A, m, rows_l, n, cols_l, values, addv));
1321:   PetscFunctionReturn(PETSC_SUCCESS);
1322: }

1324: static PetscErrorCode MatCreateSubMatrix_IS(Mat mat, IS irow, IS icol, MatReuse scall, Mat *newmat)
1325: {
1326:   Mat             locmat, newlocmat;
1327:   Mat_IS         *newmatis;
1328:   const PetscInt *idxs;
1329:   PetscInt        i, m, n;

1331:   PetscFunctionBegin;
1332:   if (scall == MAT_REUSE_MATRIX) {
1333:     PetscBool ismatis;

1335:     PetscCall(PetscObjectTypeCompare((PetscObject)*newmat, MATIS, &ismatis));
1336:     PetscCheck(ismatis, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_ARG_WRONG, "Cannot reuse matrix! Not of MATIS type");
1337:     newmatis = (Mat_IS *)(*newmat)->data;
1338:     PetscCheck(newmatis->getsub_ris, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_ARG_WRONG, "Cannot reuse matrix! Misses local row IS");
1339:     PetscCheck(newmatis->getsub_cis, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_ARG_WRONG, "Cannot reuse matrix! Misses local col IS");
1340:   }
1341:   /* irow and icol may not have duplicate entries */
1342:   if (PetscDefined(USE_DEBUG)) {
1343:     Vec                rtest, ltest;
1344:     const PetscScalar *array;

1346:     PetscCall(MatCreateVecs(mat, &ltest, &rtest));
1347:     PetscCall(ISGetLocalSize(irow, &n));
1348:     PetscCall(ISGetIndices(irow, &idxs));
1349:     for (i = 0; i < n; i++) PetscCall(VecSetValue(rtest, idxs[i], 1.0, ADD_VALUES));
1350:     PetscCall(VecAssemblyBegin(rtest));
1351:     PetscCall(VecAssemblyEnd(rtest));
1352:     PetscCall(VecGetLocalSize(rtest, &n));
1353:     PetscCall(VecGetOwnershipRange(rtest, &m, NULL));
1354:     PetscCall(VecGetArrayRead(rtest, &array));
1355:     for (i = 0; i < n; i++) PetscCheck(array[i] == 0. || array[i] == 1., PETSC_COMM_SELF, PETSC_ERR_SUP, "Index %" PetscInt_FMT " counted %" PetscInt_FMT " times! Irow may not have duplicate entries", i + m, (PetscInt)PetscRealPart(array[i]));
1356:     PetscCall(VecRestoreArrayRead(rtest, &array));
1357:     PetscCall(ISRestoreIndices(irow, &idxs));
1358:     PetscCall(ISGetLocalSize(icol, &n));
1359:     PetscCall(ISGetIndices(icol, &idxs));
1360:     for (i = 0; i < n; i++) PetscCall(VecSetValue(ltest, idxs[i], 1.0, ADD_VALUES));
1361:     PetscCall(VecAssemblyBegin(ltest));
1362:     PetscCall(VecAssemblyEnd(ltest));
1363:     PetscCall(VecGetLocalSize(ltest, &n));
1364:     PetscCall(VecGetOwnershipRange(ltest, &m, NULL));
1365:     PetscCall(VecGetArrayRead(ltest, &array));
1366:     for (i = 0; i < n; i++) PetscCheck(array[i] == 0. || array[i] == 1., PETSC_COMM_SELF, PETSC_ERR_SUP, "Index %" PetscInt_FMT " counted %" PetscInt_FMT " times! Icol may not have duplicate entries", i + m, (PetscInt)PetscRealPart(array[i]));
1367:     PetscCall(VecRestoreArrayRead(ltest, &array));
1368:     PetscCall(ISRestoreIndices(icol, &idxs));
1369:     PetscCall(VecDestroy(&rtest));
1370:     PetscCall(VecDestroy(&ltest));
1371:   }
1372:   if (scall == MAT_INITIAL_MATRIX) {
1373:     Mat_IS                *matis = (Mat_IS *)mat->data;
1374:     ISLocalToGlobalMapping rl2g;
1375:     IS                     is;
1376:     PetscInt              *lidxs, *lgidxs, *newgidxs;
1377:     PetscInt               ll, newloc, irbs, icbs, arbs, acbs, rbs, cbs;
1378:     PetscBool              cong;
1379:     MPI_Comm               comm;

1381:     PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
1382:     PetscCall(MatGetBlockSizes(mat, &arbs, &acbs));
1383:     PetscCall(ISGetBlockSize(irow, &irbs));
1384:     PetscCall(ISGetBlockSize(icol, &icbs));
1385:     rbs = arbs == irbs ? irbs : 1;
1386:     cbs = acbs == icbs ? icbs : 1;
1387:     PetscCall(ISGetLocalSize(irow, &m));
1388:     PetscCall(ISGetLocalSize(icol, &n));
1389:     PetscCall(MatCreate(comm, newmat));
1390:     PetscCall(MatSetType(*newmat, MATIS));
1391:     PetscCall(MatISSetAllowRepeated(*newmat, matis->allow_repeated));
1392:     PetscCall(MatSetSizes(*newmat, m, n, PETSC_DECIDE, PETSC_DECIDE));
1393:     PetscCall(MatSetBlockSizes(*newmat, rbs, cbs));
1394:     /* communicate irow to their owners in the layout */
1395:     PetscCall(ISGetIndices(irow, &idxs));
1396:     PetscCall(PetscLayoutMapLocal(mat->rmap, m, idxs, &ll, &lidxs, &lgidxs));
1397:     PetscCall(ISRestoreIndices(irow, &idxs));
1398:     PetscCall(PetscArrayzero(matis->sf_rootdata, matis->sf->nroots));
1399:     for (i = 0; i < ll; i++) matis->sf_rootdata[lidxs[i]] = lgidxs[i] + 1;
1400:     PetscCall(PetscFree(lidxs));
1401:     PetscCall(PetscFree(lgidxs));
1402:     PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
1403:     PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
1404:     for (i = 0, newloc = 0; i < matis->sf->nleaves; i++)
1405:       if (matis->sf_leafdata[i]) newloc++;
1406:     PetscCall(PetscMalloc1(newloc, &newgidxs));
1407:     PetscCall(PetscMalloc1(newloc, &lidxs));
1408:     for (i = 0, newloc = 0; i < matis->sf->nleaves; i++)
1409:       if (matis->sf_leafdata[i]) {
1410:         lidxs[newloc]      = i;
1411:         newgidxs[newloc++] = matis->sf_leafdata[i] - 1;
1412:       }
1413:     PetscCall(ISCreateGeneral(comm, newloc, newgidxs, PETSC_OWN_POINTER, &is));
1414:     PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g));
1415:     PetscCall(ISLocalToGlobalMappingSetBlockSize(rl2g, rbs));
1416:     PetscCall(ISDestroy(&is));
1417:     /* local is to extract local submatrix */
1418:     newmatis = (Mat_IS *)(*newmat)->data;
1419:     PetscCall(ISCreateGeneral(comm, newloc, lidxs, PETSC_OWN_POINTER, &newmatis->getsub_ris));
1420:     PetscCall(MatHasCongruentLayouts(mat, &cong));
1421:     if (cong && irow == icol && matis->csf == matis->sf) {
1422:       PetscCall(MatSetLocalToGlobalMapping(*newmat, rl2g, rl2g));
1423:       PetscCall(PetscObjectReference((PetscObject)newmatis->getsub_ris));
1424:       newmatis->getsub_cis = newmatis->getsub_ris;
1425:     } else {
1426:       ISLocalToGlobalMapping cl2g;

1428:       /* communicate icol to their owners in the layout */
1429:       PetscCall(ISGetIndices(icol, &idxs));
1430:       PetscCall(PetscLayoutMapLocal(mat->cmap, n, idxs, &ll, &lidxs, &lgidxs));
1431:       PetscCall(ISRestoreIndices(icol, &idxs));
1432:       PetscCall(PetscArrayzero(matis->csf_rootdata, matis->csf->nroots));
1433:       for (i = 0; i < ll; i++) matis->csf_rootdata[lidxs[i]] = lgidxs[i] + 1;
1434:       PetscCall(PetscFree(lidxs));
1435:       PetscCall(PetscFree(lgidxs));
1436:       PetscCall(PetscSFBcastBegin(matis->csf, MPIU_INT, matis->csf_rootdata, matis->csf_leafdata, MPI_REPLACE));
1437:       PetscCall(PetscSFBcastEnd(matis->csf, MPIU_INT, matis->csf_rootdata, matis->csf_leafdata, MPI_REPLACE));
1438:       for (i = 0, newloc = 0; i < matis->csf->nleaves; i++)
1439:         if (matis->csf_leafdata[i]) newloc++;
1440:       PetscCall(PetscMalloc1(newloc, &newgidxs));
1441:       PetscCall(PetscMalloc1(newloc, &lidxs));
1442:       for (i = 0, newloc = 0; i < matis->csf->nleaves; i++)
1443:         if (matis->csf_leafdata[i]) {
1444:           lidxs[newloc]      = i;
1445:           newgidxs[newloc++] = matis->csf_leafdata[i] - 1;
1446:         }
1447:       PetscCall(ISCreateGeneral(comm, newloc, newgidxs, PETSC_OWN_POINTER, &is));
1448:       PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g));
1449:       PetscCall(ISLocalToGlobalMappingSetBlockSize(cl2g, cbs));
1450:       PetscCall(ISDestroy(&is));
1451:       /* local is to extract local submatrix */
1452:       PetscCall(ISCreateGeneral(comm, newloc, lidxs, PETSC_OWN_POINTER, &newmatis->getsub_cis));
1453:       PetscCall(MatSetLocalToGlobalMapping(*newmat, rl2g, cl2g));
1454:       PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
1455:     }
1456:     PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
1457:   } else {
1458:     PetscCall(MatISGetLocalMat(*newmat, &newlocmat));
1459:   }
1460:   PetscCall(MatISGetLocalMat(mat, &locmat));
1461:   newmatis = (Mat_IS *)(*newmat)->data;
1462:   PetscCall(MatCreateSubMatrix(locmat, newmatis->getsub_ris, newmatis->getsub_cis, scall, &newlocmat));
1463:   if (scall == MAT_INITIAL_MATRIX) {
1464:     PetscCall(MatISSetLocalMat(*newmat, newlocmat));
1465:     PetscCall(MatDestroy(&newlocmat));
1466:   }
1467:   PetscCall(MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY));
1468:   PetscCall(MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY));
1469:   PetscFunctionReturn(PETSC_SUCCESS);
1470: }

1472: static PetscErrorCode MatCopy_IS(Mat A, Mat B, MatStructure str)
1473: {
1474:   Mat_IS   *a = (Mat_IS *)A->data, *b;
1475:   PetscBool ismatis;

1477:   PetscFunctionBegin;
1478:   PetscCall(PetscObjectTypeCompare((PetscObject)B, MATIS, &ismatis));
1479:   PetscCheck(ismatis, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Need to be implemented");
1480:   b = (Mat_IS *)B->data;
1481:   PetscCall(MatCopy(a->A, b->A, str));
1482:   PetscCall(PetscObjectStateIncrease((PetscObject)B));
1483:   PetscFunctionReturn(PETSC_SUCCESS);
1484: }

1486: static PetscErrorCode MatMissingDiagonal_IS(Mat A, PetscBool *missing, PetscInt *d)
1487: {
1488:   Vec                v;
1489:   const PetscScalar *array;
1490:   PetscInt           i, n;

1492:   PetscFunctionBegin;
1493:   *missing = PETSC_FALSE;
1494:   PetscCall(MatCreateVecs(A, NULL, &v));
1495:   PetscCall(MatGetDiagonal(A, v));
1496:   PetscCall(VecGetLocalSize(v, &n));
1497:   PetscCall(VecGetArrayRead(v, &array));
1498:   for (i = 0; i < n; i++)
1499:     if (array[i] == 0.) break;
1500:   PetscCall(VecRestoreArrayRead(v, &array));
1501:   PetscCall(VecDestroy(&v));
1502:   if (i != n) *missing = PETSC_TRUE;
1503:   if (d) {
1504:     *d = -1;
1505:     if (*missing) {
1506:       PetscInt rstart;
1507:       PetscCall(MatGetOwnershipRange(A, &rstart, NULL));
1508:       *d = i + rstart;
1509:     }
1510:   }
1511:   PetscFunctionReturn(PETSC_SUCCESS);
1512: }

1514: static PetscErrorCode MatISSetUpSF_IS(Mat B)
1515: {
1516:   Mat_IS         *matis = (Mat_IS *)B->data;
1517:   const PetscInt *gidxs;
1518:   PetscInt        nleaves;

1520:   PetscFunctionBegin;
1521:   if (matis->sf) PetscFunctionReturn(PETSC_SUCCESS);
1522:   PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)B), &matis->sf));
1523:   PetscCall(ISLocalToGlobalMappingGetIndices(matis->rmapping, &gidxs));
1524:   PetscCall(ISLocalToGlobalMappingGetSize(matis->rmapping, &nleaves));
1525:   PetscCall(PetscSFSetGraphLayout(matis->sf, B->rmap, nleaves, NULL, PETSC_OWN_POINTER, gidxs));
1526:   PetscCall(ISLocalToGlobalMappingRestoreIndices(matis->rmapping, &gidxs));
1527:   PetscCall(PetscMalloc2(matis->sf->nroots, &matis->sf_rootdata, matis->sf->nleaves, &matis->sf_leafdata));
1528:   if (matis->rmapping != matis->cmapping) { /* setup SF for columns */
1529:     PetscCall(ISLocalToGlobalMappingGetSize(matis->cmapping, &nleaves));
1530:     PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)B), &matis->csf));
1531:     PetscCall(ISLocalToGlobalMappingGetIndices(matis->cmapping, &gidxs));
1532:     PetscCall(PetscSFSetGraphLayout(matis->csf, B->cmap, nleaves, NULL, PETSC_OWN_POINTER, gidxs));
1533:     PetscCall(ISLocalToGlobalMappingRestoreIndices(matis->cmapping, &gidxs));
1534:     PetscCall(PetscMalloc2(matis->csf->nroots, &matis->csf_rootdata, matis->csf->nleaves, &matis->csf_leafdata));
1535:   } else {
1536:     matis->csf          = matis->sf;
1537:     matis->csf_leafdata = matis->sf_leafdata;
1538:     matis->csf_rootdata = matis->sf_rootdata;
1539:   }
1540:   PetscFunctionReturn(PETSC_SUCCESS);
1541: }

1543: /*@
1544:   MatISGetAllowRepeated - Get the flag to allow repeated entries in the local to global map

1546:   Not Collective

1548:   Input Parameter:
1549: . A - the matrix

1551:   Output Parameter:
1552: . flg - the boolean flag

1554:   Level: intermediate

1556: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateIS()`, `MatSetLocalToGlobalMapping()`, `MatISSetAllowRepeated()`
1557: @*/
1558: PetscErrorCode MatISGetAllowRepeated(Mat A, PetscBool *flg)
1559: {
1560:   PetscBool ismatis;

1562:   PetscFunctionBegin;
1564:   PetscAssertPointer(flg, 2);
1565:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATIS, &ismatis));
1566:   PetscCheck(ismatis, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not for matrix type %s", ((PetscObject)A)->type_name);
1567:   *flg = ((Mat_IS *)A->data)->allow_repeated;
1568:   PetscFunctionReturn(PETSC_SUCCESS);
1569: }

1571: /*@
1572:   MatISSetAllowRepeated - Set the flag to allow repeated entries in the local to global map

1574:   Logically Collective

1576:   Input Parameters:
1577: + A   - the matrix
1578: - flg - the boolean flag

1580:   Level: intermediate

1582:   Notes:
1583:   The default value is `PETSC_FALSE`.
1584:   When called AFTER calling `MatSetLocalToGlobalMapping()` it will recreate the local matrices
1585:   if `flg` is different from the previously set value.
1586:   Specifically, when `flg` is true it will just recreate the local matrices, while if
1587:   `flg` is false will assemble the local matrices summing up repeated entries.

1589: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateIS()`, `MatSetLocalToGlobalMapping()`, `MatISGetAllowRepeated()`
1590: @*/
1591: PetscErrorCode MatISSetAllowRepeated(Mat A, PetscBool flg)
1592: {
1593:   PetscFunctionBegin;
1597:   PetscTryMethod(A, "MatISSetAllowRepeated_C", (Mat, PetscBool), (A, flg));
1598:   PetscFunctionReturn(PETSC_SUCCESS);
1599: }

1601: static PetscErrorCode MatISSetAllowRepeated_IS(Mat A, PetscBool flg)
1602: {
1603:   Mat_IS                *matis = (Mat_IS *)A->data;
1604:   Mat                    lA    = NULL;
1605:   ISLocalToGlobalMapping lrmap, lcmap;

1607:   PetscFunctionBegin;
1608:   if (flg == matis->allow_repeated) PetscFunctionReturn(PETSC_SUCCESS);
1609:   if (!matis->A) { /* matrix has not been preallocated yet */
1610:     matis->allow_repeated = flg;
1611:     PetscFunctionReturn(PETSC_SUCCESS);
1612:   }
1613:   PetscCheck(!matis->islocalref, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not implemented for local references");
1614:   if (matis->allow_repeated) { /* we will assemble the old local matrix if needed */
1615:     lA = matis->A;
1616:     PetscCall(PetscObjectReference((PetscObject)lA));
1617:   }
1618:   /* In case flg is True, we only recreate the local matrix */
1619:   matis->allow_repeated = flg;
1620:   PetscCall(MatSetLocalToGlobalMapping(A, A->rmap->mapping, A->cmap->mapping));
1621:   if (lA) { /* assemble previous local matrix if needed */
1622:     Mat nA = matis->A;

1624:     PetscCall(MatGetLocalToGlobalMapping(nA, &lrmap, &lcmap));
1625:     if (!lrmap && !lcmap) {
1626:       PetscCall(MatISSetLocalMat(A, lA));
1627:     } else {
1628:       Mat            P = NULL, R = NULL;
1629:       MatProductType ptype;

1631:       if (lrmap == lcmap) {
1632:         ptype = MATPRODUCT_PtAP;
1633:         PetscCall(MatCreateFromISLocalToGlobalMapping(lcmap, nA, PETSC_TRUE, PETSC_FALSE, NULL, &P));
1634:         PetscCall(MatProductCreate(lA, P, NULL, &nA));
1635:       } else {
1636:         if (lcmap) PetscCall(MatCreateFromISLocalToGlobalMapping(lcmap, nA, PETSC_TRUE, PETSC_FALSE, NULL, &P));
1637:         if (lrmap) PetscCall(MatCreateFromISLocalToGlobalMapping(lrmap, nA, PETSC_FALSE, PETSC_TRUE, NULL, &R));
1638:         if (R && P) {
1639:           ptype = MATPRODUCT_ABC;
1640:           PetscCall(MatProductCreate(R, lA, P, &nA));
1641:         } else if (R) {
1642:           ptype = MATPRODUCT_AB;
1643:           PetscCall(MatProductCreate(R, lA, NULL, &nA));
1644:         } else {
1645:           ptype = MATPRODUCT_AB;
1646:           PetscCall(MatProductCreate(lA, P, NULL, &nA));
1647:         }
1648:       }
1649:       PetscCall(MatProductSetType(nA, ptype));
1650:       PetscCall(MatProductSetFromOptions(nA));
1651:       PetscCall(MatProductSymbolic(nA));
1652:       PetscCall(MatProductNumeric(nA));
1653:       PetscCall(MatProductClear(nA));
1654:       PetscCall(MatConvert(nA, matis->lmattype, MAT_INPLACE_MATRIX, &nA));
1655:       PetscCall(MatISSetLocalMat(A, nA));
1656:       PetscCall(MatDestroy(&nA));
1657:       PetscCall(MatDestroy(&P));
1658:       PetscCall(MatDestroy(&R));
1659:     }
1660:   }
1661:   PetscCall(MatDestroy(&lA));
1662:   PetscFunctionReturn(PETSC_SUCCESS);
1663: }

1665: /*@
1666:   MatISStoreL2L - Store local-to-local operators during the Galerkin process of computing `MatPtAP()`

1668:   Logically Collective

1670:   Input Parameters:
1671: + A     - the matrix
1672: - store - the boolean flag

1674:   Level: advanced

1676: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateIS()`, `MatISSetPreallocation()`, `MatPtAP()`
1677: @*/
1678: PetscErrorCode MatISStoreL2L(Mat A, PetscBool store)
1679: {
1680:   PetscFunctionBegin;
1684:   PetscTryMethod(A, "MatISStoreL2L_C", (Mat, PetscBool), (A, store));
1685:   PetscFunctionReturn(PETSC_SUCCESS);
1686: }

1688: static PetscErrorCode MatISStoreL2L_IS(Mat A, PetscBool store)
1689: {
1690:   Mat_IS *matis = (Mat_IS *)A->data;

1692:   PetscFunctionBegin;
1693:   matis->storel2l = store;
1694:   if (!store) PetscCall(PetscObjectCompose((PetscObject)(A), "_MatIS_PtAP_l2l", NULL));
1695:   PetscFunctionReturn(PETSC_SUCCESS);
1696: }

1698: /*@
1699:   MatISFixLocalEmpty - Compress out zero local rows from the local matrices

1701:   Logically Collective

1703:   Input Parameters:
1704: + A   - the matrix
1705: - fix - the boolean flag

1707:   Level: advanced

1709:   Note:
1710:   When `fix` is `PETSC_TRUE`, new local matrices and l2g maps are generated during the final assembly process.

1712: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatCreate()`, `MatCreateIS()`, `MatISSetPreallocation()`, `MatAssemblyEnd()`, `MAT_FINAL_ASSEMBLY`
1713: @*/
1714: PetscErrorCode MatISFixLocalEmpty(Mat A, PetscBool fix)
1715: {
1716:   PetscFunctionBegin;
1720:   PetscTryMethod(A, "MatISFixLocalEmpty_C", (Mat, PetscBool), (A, fix));
1721:   PetscFunctionReturn(PETSC_SUCCESS);
1722: }

1724: static PetscErrorCode MatISFixLocalEmpty_IS(Mat A, PetscBool fix)
1725: {
1726:   Mat_IS *matis = (Mat_IS *)A->data;

1728:   PetscFunctionBegin;
1729:   matis->locempty = fix;
1730:   PetscFunctionReturn(PETSC_SUCCESS);
1731: }

1733: /*@
1734:   MatISSetPreallocation - Preallocates memory for a `MATIS` parallel matrix.

1736:   Collective

1738:   Input Parameters:
1739: + B     - the matrix
1740: . d_nz  - number of nonzeros per row in DIAGONAL portion of local submatrix
1741:            (same value is used for all local rows)
1742: . d_nnz - array containing the number of nonzeros in the various rows of the
1743:            DIAGONAL portion of the local submatrix (possibly different for each row)
1744:            or `NULL`, if `d_nz` is used to specify the nonzero structure.
1745:            The size of this array is equal to the number of local rows, i.e `m`.
1746:            For matrices that will be factored, you must leave room for (and set)
1747:            the diagonal entry even if it is zero.
1748: . o_nz  - number of nonzeros per row in the OFF-DIAGONAL portion of local
1749:            submatrix (same value is used for all local rows).
1750: - o_nnz - array containing the number of nonzeros in the various rows of the
1751:            OFF-DIAGONAL portion of the local submatrix (possibly different for
1752:            each row) or `NULL`, if `o_nz` is used to specify the nonzero
1753:            structure. The size of this array is equal to the number
1754:            of local rows, i.e `m`.

1756:    If the *_nnz parameter is given then the *_nz parameter is ignored

1758:   Level: intermediate

1760:   Note:
1761:   This function has the same interface as the `MATMPIAIJ` preallocation routine in order to simplify the transition
1762:   from the asssembled format to the unassembled one. It overestimates the preallocation of `MATIS` local
1763:   matrices; for exact preallocation, the user should set the preallocation directly on local matrix objects.

1765: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateIS()`, `MatMPIAIJSetPreallocation()`, `MatISGetLocalMat()`, `MATIS`
1766: @*/
1767: PetscErrorCode MatISSetPreallocation(Mat B, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[])
1768: {
1769:   PetscFunctionBegin;
1772:   PetscTryMethod(B, "MatISSetPreallocation_C", (Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[]), (B, d_nz, d_nnz, o_nz, o_nnz));
1773:   PetscFunctionReturn(PETSC_SUCCESS);
1774: }

1776: /* this is used by DMDA */
1777: PETSC_EXTERN PetscErrorCode MatISSetPreallocation_IS(Mat B, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[])
1778: {
1779:   Mat_IS  *matis = (Mat_IS *)B->data;
1780:   PetscInt bs, i, nlocalcols;

1782:   PetscFunctionBegin;
1783:   PetscCall(MatSetUp(B));
1784:   if (!d_nnz)
1785:     for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] = d_nz;
1786:   else
1787:     for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] = d_nnz[i];

1789:   if (!o_nnz)
1790:     for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] += o_nz;
1791:   else
1792:     for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] += o_nnz[i];

1794:   PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
1795:   PetscCall(MatGetSize(matis->A, NULL, &nlocalcols));
1796:   PetscCall(MatGetBlockSize(matis->A, &bs));
1797:   PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));

1799:   for (i = 0; i < matis->sf->nleaves; i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i], nlocalcols);
1800:   PetscCall(MatSeqAIJSetPreallocation(matis->A, 0, matis->sf_leafdata));
1801: #if defined(PETSC_HAVE_HYPRE)
1802:   PetscCall(MatHYPRESetPreallocation(matis->A, 0, matis->sf_leafdata, 0, NULL));
1803: #endif

1805:   for (i = 0; i < matis->sf->nleaves / bs; i++) {
1806:     PetscInt b;

1808:     matis->sf_leafdata[i] = matis->sf_leafdata[i * bs] / bs;
1809:     for (b = 1; b < bs; b++) matis->sf_leafdata[i] = PetscMax(matis->sf_leafdata[i], matis->sf_leafdata[i * bs + b] / bs);
1810:   }
1811:   PetscCall(MatSeqBAIJSetPreallocation(matis->A, bs, 0, matis->sf_leafdata));

1813:   nlocalcols /= bs;
1814:   for (i = 0; i < matis->sf->nleaves / bs; i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i], nlocalcols - i);
1815:   PetscCall(MatSeqSBAIJSetPreallocation(matis->A, bs, 0, matis->sf_leafdata));

1817:   /* for other matrix types */
1818:   PetscCall(MatSetUp(matis->A));
1819:   PetscFunctionReturn(PETSC_SUCCESS);
1820: }

1822: PETSC_EXTERN PetscErrorCode MatISSetMPIXAIJPreallocation_Private(Mat A, Mat B, PetscBool maxreduce)
1823: {
1824:   Mat_IS         *matis = (Mat_IS *)A->data;
1825:   PetscInt       *my_dnz, *my_onz, *dnz, *onz, *mat_ranges, *row_ownership;
1826:   const PetscInt *global_indices_r, *global_indices_c;
1827:   PetscInt        i, j, bs, rows, cols;
1828:   PetscInt        lrows, lcols;
1829:   PetscInt        local_rows, local_cols;
1830:   PetscMPIInt     size;
1831:   PetscBool       isdense, issbaij;

1833:   PetscFunctionBegin;
1834:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1835:   PetscCall(MatGetSize(A, &rows, &cols));
1836:   PetscCall(MatGetBlockSize(A, &bs));
1837:   PetscCall(MatGetSize(matis->A, &local_rows, &local_cols));
1838:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQDENSE, &isdense));
1839:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQSBAIJ, &issbaij));
1840:   PetscCall(ISLocalToGlobalMappingGetIndices(matis->rmapping, &global_indices_r));
1841:   if (matis->rmapping != matis->cmapping) {
1842:     PetscCall(ISLocalToGlobalMappingGetIndices(matis->cmapping, &global_indices_c));
1843:   } else global_indices_c = global_indices_r;

1845:   if (issbaij) PetscCall(MatGetRowUpperTriangular(matis->A));
1846:   /*
1847:      An SF reduce is needed to sum up properly on shared rows.
1848:      Note that generally preallocation is not exact, since it overestimates nonzeros
1849:   */
1850:   PetscCall(MatGetLocalSize(A, &lrows, &lcols));
1851:   MatPreallocateBegin(PetscObjectComm((PetscObject)A), lrows, lcols, dnz, onz);
1852:   /* All processes need to compute entire row ownership */
1853:   PetscCall(PetscMalloc1(rows, &row_ownership));
1854:   PetscCall(MatGetOwnershipRanges(A, (const PetscInt **)&mat_ranges));
1855:   for (i = 0; i < size; i++) {
1856:     for (j = mat_ranges[i]; j < mat_ranges[i + 1]; j++) row_ownership[j] = i;
1857:   }
1858:   PetscCall(MatGetOwnershipRangesColumn(A, (const PetscInt **)&mat_ranges));

1860:   /*
1861:      my_dnz and my_onz contains exact contribution to preallocation from each local mat
1862:      then, they will be summed up properly. This way, preallocation is always sufficient
1863:   */
1864:   PetscCall(PetscCalloc2(local_rows, &my_dnz, local_rows, &my_onz));
1865:   /* preallocation as a MATAIJ */
1866:   if (isdense) { /* special case for dense local matrices */
1867:     for (i = 0; i < local_rows; i++) {
1868:       PetscInt owner = row_ownership[global_indices_r[i]];
1869:       for (j = 0; j < local_cols; j++) {
1870:         PetscInt index_col = global_indices_c[j];
1871:         if (index_col > mat_ranges[owner] - 1 && index_col < mat_ranges[owner + 1]) { /* diag block */
1872:           my_dnz[i] += 1;
1873:         } else { /* offdiag block */
1874:           my_onz[i] += 1;
1875:         }
1876:       }
1877:     }
1878:   } else if (matis->A->ops->getrowij) {
1879:     const PetscInt *ii, *jj, *jptr;
1880:     PetscBool       done;
1881:     PetscCall(MatGetRowIJ(matis->A, 0, PETSC_FALSE, PETSC_FALSE, &local_rows, &ii, &jj, &done));
1882:     PetscCheck(done, PetscObjectComm((PetscObject)matis->A), PETSC_ERR_PLIB, "Error in MatGetRowIJ");
1883:     jptr = jj;
1884:     for (i = 0; i < local_rows; i++) {
1885:       PetscInt index_row = global_indices_r[i];
1886:       for (j = 0; j < ii[i + 1] - ii[i]; j++, jptr++) {
1887:         PetscInt owner     = row_ownership[index_row];
1888:         PetscInt index_col = global_indices_c[*jptr];
1889:         if (index_col > mat_ranges[owner] - 1 && index_col < mat_ranges[owner + 1]) { /* diag block */
1890:           my_dnz[i] += 1;
1891:         } else { /* offdiag block */
1892:           my_onz[i] += 1;
1893:         }
1894:         /* same as before, interchanging rows and cols */
1895:         if (issbaij && index_col != index_row) {
1896:           owner = row_ownership[index_col];
1897:           if (index_row > mat_ranges[owner] - 1 && index_row < mat_ranges[owner + 1]) {
1898:             my_dnz[*jptr] += 1;
1899:           } else {
1900:             my_onz[*jptr] += 1;
1901:           }
1902:         }
1903:       }
1904:     }
1905:     PetscCall(MatRestoreRowIJ(matis->A, 0, PETSC_FALSE, PETSC_FALSE, &local_rows, &ii, &jj, &done));
1906:     PetscCheck(done, PetscObjectComm((PetscObject)matis->A), PETSC_ERR_PLIB, "Error in MatRestoreRowIJ");
1907:   } else { /* loop over rows and use MatGetRow */
1908:     for (i = 0; i < local_rows; i++) {
1909:       const PetscInt *cols;
1910:       PetscInt        ncols, index_row = global_indices_r[i];
1911:       PetscCall(MatGetRow(matis->A, i, &ncols, &cols, NULL));
1912:       for (j = 0; j < ncols; j++) {
1913:         PetscInt owner     = row_ownership[index_row];
1914:         PetscInt index_col = global_indices_c[cols[j]];
1915:         if (index_col > mat_ranges[owner] - 1 && index_col < mat_ranges[owner + 1]) { /* diag block */
1916:           my_dnz[i] += 1;
1917:         } else { /* offdiag block */
1918:           my_onz[i] += 1;
1919:         }
1920:         /* same as before, interchanging rows and cols */
1921:         if (issbaij && index_col != index_row) {
1922:           owner = row_ownership[index_col];
1923:           if (index_row > mat_ranges[owner] - 1 && index_row < mat_ranges[owner + 1]) {
1924:             my_dnz[cols[j]] += 1;
1925:           } else {
1926:             my_onz[cols[j]] += 1;
1927:           }
1928:         }
1929:       }
1930:       PetscCall(MatRestoreRow(matis->A, i, &ncols, &cols, NULL));
1931:     }
1932:   }
1933:   if (global_indices_c != global_indices_r) PetscCall(ISLocalToGlobalMappingRestoreIndices(matis->cmapping, &global_indices_c));
1934:   PetscCall(ISLocalToGlobalMappingRestoreIndices(matis->rmapping, &global_indices_r));
1935:   PetscCall(PetscFree(row_ownership));

1937:   /* Reduce my_dnz and my_onz */
1938:   if (maxreduce) {
1939:     PetscCall(PetscSFReduceBegin(matis->sf, MPIU_INT, my_dnz, dnz, MPI_MAX));
1940:     PetscCall(PetscSFReduceBegin(matis->sf, MPIU_INT, my_onz, onz, MPI_MAX));
1941:     PetscCall(PetscSFReduceEnd(matis->sf, MPIU_INT, my_dnz, dnz, MPI_MAX));
1942:     PetscCall(PetscSFReduceEnd(matis->sf, MPIU_INT, my_onz, onz, MPI_MAX));
1943:   } else {
1944:     PetscCall(PetscSFReduceBegin(matis->sf, MPIU_INT, my_dnz, dnz, MPI_SUM));
1945:     PetscCall(PetscSFReduceBegin(matis->sf, MPIU_INT, my_onz, onz, MPI_SUM));
1946:     PetscCall(PetscSFReduceEnd(matis->sf, MPIU_INT, my_dnz, dnz, MPI_SUM));
1947:     PetscCall(PetscSFReduceEnd(matis->sf, MPIU_INT, my_onz, onz, MPI_SUM));
1948:   }
1949:   PetscCall(PetscFree2(my_dnz, my_onz));

1951:   /* Resize preallocation if overestimated */
1952:   for (i = 0; i < lrows; i++) {
1953:     dnz[i] = PetscMin(dnz[i], lcols);
1954:     onz[i] = PetscMin(onz[i], cols - lcols);
1955:   }

1957:   /* Set preallocation */
1958:   PetscCall(MatSetBlockSizesFromMats(B, A, A));
1959:   PetscCall(MatSeqAIJSetPreallocation(B, 0, dnz));
1960:   PetscCall(MatMPIAIJSetPreallocation(B, 0, dnz, 0, onz));
1961:   for (i = 0; i < lrows; i += bs) {
1962:     PetscInt b, d = dnz[i], o = onz[i];

1964:     for (b = 1; b < bs; b++) {
1965:       d = PetscMax(d, dnz[i + b]);
1966:       o = PetscMax(o, onz[i + b]);
1967:     }
1968:     dnz[i / bs] = PetscMin(d / bs + d % bs, lcols / bs);
1969:     onz[i / bs] = PetscMin(o / bs + o % bs, (cols - lcols) / bs);
1970:   }
1971:   PetscCall(MatSeqBAIJSetPreallocation(B, bs, 0, dnz));
1972:   PetscCall(MatMPIBAIJSetPreallocation(B, bs, 0, dnz, 0, onz));
1973:   PetscCall(MatMPISBAIJSetPreallocation(B, bs, 0, dnz, 0, onz));
1974:   MatPreallocateEnd(dnz, onz);
1975:   if (issbaij) PetscCall(MatRestoreRowUpperTriangular(matis->A));
1976:   PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
1977:   PetscFunctionReturn(PETSC_SUCCESS);
1978: }

1980: PETSC_INTERN PetscErrorCode MatConvert_IS_XAIJ(Mat mat, MatType mtype, MatReuse reuse, Mat *M)
1981: {
1982:   Mat_IS            *matis = (Mat_IS *)mat->data;
1983:   Mat                local_mat, MT;
1984:   PetscInt           rbs, cbs, rows, cols, lrows, lcols;
1985:   PetscInt           local_rows, local_cols;
1986:   PetscBool          isseqdense, isseqsbaij, isseqaij, isseqbaij;
1987:   PetscMPIInt        size;
1988:   const PetscScalar *array;

1990:   PetscFunctionBegin;
1991:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)mat), &size));
1992:   if (size == 1 && mat->rmap->N == matis->A->rmap->N && mat->cmap->N == matis->A->cmap->N && !matis->allow_repeated) {
1993:     Mat      B;
1994:     IS       irows = NULL, icols = NULL;
1995:     PetscInt rbs, cbs;

1997:     PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->rmapping, &rbs));
1998:     PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->cmapping, &cbs));
1999:     if (reuse != MAT_REUSE_MATRIX) { /* check if l2g maps are one-to-one */
2000:       IS              rows, cols;
2001:       const PetscInt *ridxs, *cidxs;
2002:       PetscInt        i, nw;
2003:       PetscBT         work;

2005:       PetscCall(ISLocalToGlobalMappingGetBlockIndices(matis->rmapping, &ridxs));
2006:       PetscCall(ISLocalToGlobalMappingGetSize(matis->rmapping, &nw));
2007:       nw = nw / rbs;
2008:       PetscCall(PetscBTCreate(nw, &work));
2009:       for (i = 0; i < nw; i++) PetscCall(PetscBTSet(work, ridxs[i]));
2010:       for (i = 0; i < nw; i++)
2011:         if (!PetscBTLookup(work, i)) break;
2012:       if (i == nw) {
2013:         PetscCall(ISCreateBlock(PETSC_COMM_SELF, rbs, nw, ridxs, PETSC_USE_POINTER, &rows));
2014:         PetscCall(ISSetPermutation(rows));
2015:         PetscCall(ISInvertPermutation(rows, PETSC_DECIDE, &irows));
2016:         PetscCall(ISDestroy(&rows));
2017:       }
2018:       PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(matis->rmapping, &ridxs));
2019:       PetscCall(PetscBTDestroy(&work));
2020:       if (irows && matis->rmapping != matis->cmapping) {
2021:         PetscCall(ISLocalToGlobalMappingGetBlockIndices(matis->cmapping, &cidxs));
2022:         PetscCall(ISLocalToGlobalMappingGetSize(matis->cmapping, &nw));
2023:         nw = nw / cbs;
2024:         PetscCall(PetscBTCreate(nw, &work));
2025:         for (i = 0; i < nw; i++) PetscCall(PetscBTSet(work, cidxs[i]));
2026:         for (i = 0; i < nw; i++)
2027:           if (!PetscBTLookup(work, i)) break;
2028:         if (i == nw) {
2029:           PetscCall(ISCreateBlock(PETSC_COMM_SELF, cbs, nw, cidxs, PETSC_USE_POINTER, &cols));
2030:           PetscCall(ISSetPermutation(cols));
2031:           PetscCall(ISInvertPermutation(cols, PETSC_DECIDE, &icols));
2032:           PetscCall(ISDestroy(&cols));
2033:         }
2034:         PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(matis->cmapping, &cidxs));
2035:         PetscCall(PetscBTDestroy(&work));
2036:       } else if (irows) {
2037:         PetscCall(PetscObjectReference((PetscObject)irows));
2038:         icols = irows;
2039:       }
2040:     } else {
2041:       PetscCall(PetscObjectQuery((PetscObject)*M, "_MatIS_IS_XAIJ_irows", (PetscObject *)&irows));
2042:       PetscCall(PetscObjectQuery((PetscObject)*M, "_MatIS_IS_XAIJ_icols", (PetscObject *)&icols));
2043:       if (irows) PetscCall(PetscObjectReference((PetscObject)irows));
2044:       if (icols) PetscCall(PetscObjectReference((PetscObject)icols));
2045:     }
2046:     if (!irows || !icols) {
2047:       PetscCall(ISDestroy(&icols));
2048:       PetscCall(ISDestroy(&irows));
2049:       goto general_assembly;
2050:     }
2051:     PetscCall(MatConvert(matis->A, mtype, MAT_INITIAL_MATRIX, &B));
2052:     if (reuse != MAT_INPLACE_MATRIX) {
2053:       PetscCall(MatCreateSubMatrix(B, irows, icols, reuse, M));
2054:       PetscCall(PetscObjectCompose((PetscObject)*M, "_MatIS_IS_XAIJ_irows", (PetscObject)irows));
2055:       PetscCall(PetscObjectCompose((PetscObject)*M, "_MatIS_IS_XAIJ_icols", (PetscObject)icols));
2056:     } else {
2057:       Mat C;

2059:       PetscCall(MatCreateSubMatrix(B, irows, icols, MAT_INITIAL_MATRIX, &C));
2060:       PetscCall(MatHeaderReplace(mat, &C));
2061:     }
2062:     PetscCall(MatDestroy(&B));
2063:     PetscCall(ISDestroy(&icols));
2064:     PetscCall(ISDestroy(&irows));
2065:     PetscFunctionReturn(PETSC_SUCCESS);
2066:   }
2067: general_assembly:
2068:   PetscCall(MatGetSize(mat, &rows, &cols));
2069:   PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->rmapping, &rbs));
2070:   PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->cmapping, &cbs));
2071:   PetscCall(MatGetLocalSize(mat, &lrows, &lcols));
2072:   PetscCall(MatGetSize(matis->A, &local_rows, &local_cols));
2073:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQDENSE, &isseqdense));
2074:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQAIJ, &isseqaij));
2075:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQBAIJ, &isseqbaij));
2076:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQSBAIJ, &isseqsbaij));
2077:   PetscCheck(isseqdense || isseqaij || isseqbaij || isseqsbaij, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for matrix type %s", ((PetscObject)matis->A)->type_name);
2078:   if (PetscDefined(USE_DEBUG)) {
2079:     PetscBool lb[4], bb[4];

2081:     lb[0] = isseqdense;
2082:     lb[1] = isseqaij;
2083:     lb[2] = isseqbaij;
2084:     lb[3] = isseqsbaij;
2085:     PetscCall(MPIU_Allreduce(lb, bb, 4, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)mat)));
2086:     PetscCheck(bb[0] || bb[1] || bb[2] || bb[3], PETSC_COMM_SELF, PETSC_ERR_SUP, "Local matrices must have the same type");
2087:   }

2089:   if (reuse != MAT_REUSE_MATRIX) {
2090:     PetscCall(MatCreate(PetscObjectComm((PetscObject)mat), &MT));
2091:     PetscCall(MatSetSizes(MT, lrows, lcols, rows, cols));
2092:     PetscCall(MatSetType(MT, mtype));
2093:     PetscCall(MatSetBlockSizes(MT, rbs, cbs));
2094:     PetscCall(MatISSetMPIXAIJPreallocation_Private(mat, MT, PETSC_FALSE));
2095:   } else {
2096:     PetscInt mrbs, mcbs, mrows, mcols, mlrows, mlcols;

2098:     /* some checks */
2099:     MT = *M;
2100:     PetscCall(MatGetBlockSizes(MT, &mrbs, &mcbs));
2101:     PetscCall(MatGetSize(MT, &mrows, &mcols));
2102:     PetscCall(MatGetLocalSize(MT, &mlrows, &mlcols));
2103:     PetscCheck(mrows == rows, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of rows (%" PetscInt_FMT " != %" PetscInt_FMT ")", rows, mrows);
2104:     PetscCheck(mcols == cols, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of cols (%" PetscInt_FMT " != %" PetscInt_FMT ")", cols, mcols);
2105:     PetscCheck(mlrows == lrows, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of local rows (%" PetscInt_FMT " != %" PetscInt_FMT ")", lrows, mlrows);
2106:     PetscCheck(mlcols == lcols, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of local cols (%" PetscInt_FMT " != %" PetscInt_FMT ")", lcols, mlcols);
2107:     PetscCheck(mrbs == rbs, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong row block size (%" PetscInt_FMT " != %" PetscInt_FMT ")", rbs, mrbs);
2108:     PetscCheck(mcbs == cbs, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong col block size (%" PetscInt_FMT " != %" PetscInt_FMT ")", cbs, mcbs);
2109:     PetscCall(MatZeroEntries(MT));
2110:   }

2112:   if (isseqsbaij || isseqbaij) {
2113:     PetscCall(MatConvert(matis->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &local_mat));
2114:     isseqaij = PETSC_TRUE;
2115:   } else {
2116:     PetscCall(PetscObjectReference((PetscObject)matis->A));
2117:     local_mat = matis->A;
2118:   }

2120:   /* Set values */
2121:   PetscCall(MatSetLocalToGlobalMapping(MT, matis->rmapping, matis->cmapping));
2122:   if (isseqdense) { /* special case for dense local matrices */
2123:     PetscInt i, *dummy;

2125:     PetscCall(PetscMalloc1(PetscMax(local_rows, local_cols), &dummy));
2126:     for (i = 0; i < PetscMax(local_rows, local_cols); i++) dummy[i] = i;
2127:     PetscCall(MatSetOption(MT, MAT_ROW_ORIENTED, PETSC_FALSE));
2128:     PetscCall(MatDenseGetArrayRead(local_mat, &array));
2129:     PetscCall(MatSetValuesLocal(MT, local_rows, dummy, local_cols, dummy, array, ADD_VALUES));
2130:     PetscCall(MatDenseRestoreArrayRead(local_mat, &array));
2131:     PetscCall(PetscFree(dummy));
2132:   } else if (isseqaij) {
2133:     const PetscInt *blocks;
2134:     PetscInt        i, nvtxs, *xadj, *adjncy, nb;
2135:     PetscBool       done;
2136:     PetscScalar    *sarray;

2138:     PetscCall(MatGetRowIJ(local_mat, 0, PETSC_FALSE, PETSC_FALSE, &nvtxs, (const PetscInt **)&xadj, (const PetscInt **)&adjncy, &done));
2139:     PetscCheck(done, PetscObjectComm((PetscObject)local_mat), PETSC_ERR_PLIB, "Error in MatGetRowIJ");
2140:     PetscCall(MatSeqAIJGetArray(local_mat, &sarray));
2141:     PetscCall(MatGetVariableBlockSizes(local_mat, &nb, &blocks));
2142:     if (nb) { /* speed up assembly for special blocked matrices (used by BDDC) */
2143:       PetscInt sum;

2145:       for (i = 0, sum = 0; i < nb; i++) sum += blocks[i];
2146:       if (sum == nvtxs) {
2147:         PetscInt r;

2149:         for (i = 0, r = 0; i < nb; i++) {
2150:           PetscAssert(blocks[i] == xadj[r + 1] - xadj[r], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid block sizes prescribed for block %" PetscInt_FMT ": expected %" PetscInt_FMT ", got %" PetscInt_FMT, i, blocks[i], xadj[r + 1] - xadj[r]);
2151:           PetscCall(MatSetValuesLocal(MT, blocks[i], adjncy + xadj[r], blocks[i], adjncy + xadj[r], sarray + xadj[r], ADD_VALUES));
2152:           r += blocks[i];
2153:         }
2154:       } else {
2155:         for (i = 0; i < nvtxs; i++) PetscCall(MatSetValuesLocal(MT, 1, &i, xadj[i + 1] - xadj[i], adjncy + xadj[i], sarray + xadj[i], ADD_VALUES));
2156:       }
2157:     } else {
2158:       for (i = 0; i < nvtxs; i++) PetscCall(MatSetValuesLocal(MT, 1, &i, xadj[i + 1] - xadj[i], PetscSafePointerPlusOffset(adjncy, xadj[i]), PetscSafePointerPlusOffset(sarray, xadj[i]), ADD_VALUES));
2159:     }
2160:     PetscCall(MatRestoreRowIJ(local_mat, 0, PETSC_FALSE, PETSC_FALSE, &nvtxs, (const PetscInt **)&xadj, (const PetscInt **)&adjncy, &done));
2161:     PetscCheck(done, PetscObjectComm((PetscObject)local_mat), PETSC_ERR_PLIB, "Error in MatRestoreRowIJ");
2162:     PetscCall(MatSeqAIJRestoreArray(local_mat, &sarray));
2163:   } else { /* very basic values insertion for all other matrix types */
2164:     for (PetscInt i = 0; i < local_rows; i++) {
2165:       PetscInt        j;
2166:       const PetscInt *local_indices_cols;

2168:       PetscCall(MatGetRow(local_mat, i, &j, &local_indices_cols, &array));
2169:       PetscCall(MatSetValuesLocal(MT, 1, &i, j, local_indices_cols, array, ADD_VALUES));
2170:       PetscCall(MatRestoreRow(local_mat, i, &j, &local_indices_cols, &array));
2171:     }
2172:   }
2173:   PetscCall(MatDestroy(&local_mat));
2174:   PetscCall(MatAssemblyBegin(MT, MAT_FINAL_ASSEMBLY));
2175:   PetscCall(MatAssemblyEnd(MT, MAT_FINAL_ASSEMBLY));
2176:   if (isseqdense) PetscCall(MatSetOption(MT, MAT_ROW_ORIENTED, PETSC_TRUE));
2177:   if (reuse == MAT_INPLACE_MATRIX) {
2178:     PetscCall(MatHeaderReplace(mat, &MT));
2179:   } else if (reuse == MAT_INITIAL_MATRIX) {
2180:     *M = MT;
2181:   }
2182:   PetscFunctionReturn(PETSC_SUCCESS);
2183: }

2185: static PetscErrorCode MatDuplicate_IS(Mat mat, MatDuplicateOption op, Mat *newmat)
2186: {
2187:   Mat_IS  *matis = (Mat_IS *)mat->data;
2188:   PetscInt rbs, cbs, m, n, M, N;
2189:   Mat      B, localmat;

2191:   PetscFunctionBegin;
2192:   PetscCall(ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping, &rbs));
2193:   PetscCall(ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping, &cbs));
2194:   PetscCall(MatGetSize(mat, &M, &N));
2195:   PetscCall(MatGetLocalSize(mat, &m, &n));
2196:   PetscCall(MatCreate(PetscObjectComm((PetscObject)mat), &B));
2197:   PetscCall(MatSetSizes(B, m, n, M, N));
2198:   PetscCall(MatSetBlockSize(B, rbs == cbs ? rbs : 1));
2199:   PetscCall(MatSetType(B, MATIS));
2200:   PetscCall(MatISSetLocalMatType(B, matis->lmattype));
2201:   PetscCall(MatISSetAllowRepeated(B, matis->allow_repeated));
2202:   PetscCall(MatSetLocalToGlobalMapping(B, mat->rmap->mapping, mat->cmap->mapping));
2203:   PetscCall(MatDuplicate(matis->A, op, &localmat));
2204:   PetscCall(MatSetLocalToGlobalMapping(localmat, matis->A->rmap->mapping, matis->A->cmap->mapping));
2205:   PetscCall(MatISSetLocalMat(B, localmat));
2206:   PetscCall(MatDestroy(&localmat));
2207:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
2208:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
2209:   *newmat = B;
2210:   PetscFunctionReturn(PETSC_SUCCESS);
2211: }

2213: static PetscErrorCode MatIsHermitian_IS(Mat A, PetscReal tol, PetscBool *flg)
2214: {
2215:   Mat_IS   *matis = (Mat_IS *)A->data;
2216:   PetscBool local_sym;

2218:   PetscFunctionBegin;
2219:   PetscCall(MatIsHermitian(matis->A, tol, &local_sym));
2220:   PetscCall(MPIU_Allreduce(&local_sym, flg, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
2221:   PetscFunctionReturn(PETSC_SUCCESS);
2222: }

2224: static PetscErrorCode MatIsSymmetric_IS(Mat A, PetscReal tol, PetscBool *flg)
2225: {
2226:   Mat_IS   *matis = (Mat_IS *)A->data;
2227:   PetscBool local_sym;

2229:   PetscFunctionBegin;
2230:   if (matis->rmapping != matis->cmapping) {
2231:     *flg = PETSC_FALSE;
2232:     PetscFunctionReturn(PETSC_SUCCESS);
2233:   }
2234:   PetscCall(MatIsSymmetric(matis->A, tol, &local_sym));
2235:   PetscCall(MPIU_Allreduce(&local_sym, flg, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
2236:   PetscFunctionReturn(PETSC_SUCCESS);
2237: }

2239: static PetscErrorCode MatIsStructurallySymmetric_IS(Mat A, PetscBool *flg)
2240: {
2241:   Mat_IS   *matis = (Mat_IS *)A->data;
2242:   PetscBool local_sym;

2244:   PetscFunctionBegin;
2245:   if (matis->rmapping != matis->cmapping) {
2246:     *flg = PETSC_FALSE;
2247:     PetscFunctionReturn(PETSC_SUCCESS);
2248:   }
2249:   PetscCall(MatIsStructurallySymmetric(matis->A, &local_sym));
2250:   PetscCall(MPIU_Allreduce(&local_sym, flg, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
2251:   PetscFunctionReturn(PETSC_SUCCESS);
2252: }

2254: static PetscErrorCode MatDestroy_IS(Mat A)
2255: {
2256:   Mat_IS *b = (Mat_IS *)A->data;

2258:   PetscFunctionBegin;
2259:   PetscCall(PetscFree(b->bdiag));
2260:   PetscCall(PetscFree(b->lmattype));
2261:   PetscCall(MatDestroy(&b->A));
2262:   PetscCall(VecScatterDestroy(&b->cctx));
2263:   PetscCall(VecScatterDestroy(&b->rctx));
2264:   PetscCall(VecDestroy(&b->x));
2265:   PetscCall(VecDestroy(&b->y));
2266:   PetscCall(VecDestroy(&b->counter));
2267:   PetscCall(ISDestroy(&b->getsub_ris));
2268:   PetscCall(ISDestroy(&b->getsub_cis));
2269:   if (b->sf != b->csf) {
2270:     PetscCall(PetscSFDestroy(&b->csf));
2271:     PetscCall(PetscFree2(b->csf_rootdata, b->csf_leafdata));
2272:   } else b->csf = NULL;
2273:   PetscCall(PetscSFDestroy(&b->sf));
2274:   PetscCall(PetscFree2(b->sf_rootdata, b->sf_leafdata));
2275:   PetscCall(ISLocalToGlobalMappingDestroy(&b->rmapping));
2276:   PetscCall(ISLocalToGlobalMappingDestroy(&b->cmapping));
2277:   PetscCall(MatDestroy(&b->dA));
2278:   PetscCall(MatDestroy(&b->assembledA));
2279:   PetscCall(PetscFree(A->data));
2280:   PetscCall(PetscObjectChangeTypeName((PetscObject)A, NULL));
2281:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMatType_C", NULL));
2282:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalMat_C", NULL));
2283:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMat_C", NULL));
2284:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISRestoreLocalMat_C", NULL));
2285:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetPreallocation_C", NULL));
2286:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISStoreL2L_C", NULL));
2287:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISFixLocalEmpty_C", NULL));
2288:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalToGlobalMapping_C", NULL));
2289:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpiaij_C", NULL));
2290:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpibaij_C", NULL));
2291:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpisbaij_C", NULL));
2292:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqaij_C", NULL));
2293:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqbaij_C", NULL));
2294:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqsbaij_C", NULL));
2295:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_aij_C", NULL));
2296:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOOLocal_C", NULL));
2297:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOO_C", NULL));
2298:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", NULL));
2299:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetAllowRepeated_C", NULL));
2300:   PetscFunctionReturn(PETSC_SUCCESS);
2301: }

2303: static PetscErrorCode MatMult_IS(Mat A, Vec x, Vec y)
2304: {
2305:   Mat_IS     *is   = (Mat_IS *)A->data;
2306:   PetscScalar zero = 0.0;

2308:   PetscFunctionBegin;
2309:   /*  scatter the global vector x into the local work vector */
2310:   PetscCall(VecScatterBegin(is->cctx, x, is->x, INSERT_VALUES, SCATTER_FORWARD));
2311:   PetscCall(VecScatterEnd(is->cctx, x, is->x, INSERT_VALUES, SCATTER_FORWARD));

2313:   /* multiply the local matrix */
2314:   PetscCall(MatMult(is->A, is->x, is->y));

2316:   /* scatter product back into global memory */
2317:   PetscCall(VecSet(y, zero));
2318:   PetscCall(VecScatterBegin(is->rctx, is->y, y, ADD_VALUES, SCATTER_REVERSE));
2319:   PetscCall(VecScatterEnd(is->rctx, is->y, y, ADD_VALUES, SCATTER_REVERSE));
2320:   PetscFunctionReturn(PETSC_SUCCESS);
2321: }

2323: static PetscErrorCode MatMultAdd_IS(Mat A, Vec v1, Vec v2, Vec v3)
2324: {
2325:   Vec temp_vec;

2327:   PetscFunctionBegin; /*  v3 = v2 + A * v1.*/
2328:   if (v3 != v2) {
2329:     PetscCall(MatMult(A, v1, v3));
2330:     PetscCall(VecAXPY(v3, 1.0, v2));
2331:   } else {
2332:     PetscCall(VecDuplicate(v2, &temp_vec));
2333:     PetscCall(MatMult(A, v1, temp_vec));
2334:     PetscCall(VecAXPY(temp_vec, 1.0, v2));
2335:     PetscCall(VecCopy(temp_vec, v3));
2336:     PetscCall(VecDestroy(&temp_vec));
2337:   }
2338:   PetscFunctionReturn(PETSC_SUCCESS);
2339: }

2341: static PetscErrorCode MatMultTranspose_IS(Mat A, Vec y, Vec x)
2342: {
2343:   Mat_IS *is = (Mat_IS *)A->data;

2345:   PetscFunctionBegin;
2346:   /*  scatter the global vector x into the local work vector */
2347:   PetscCall(VecScatterBegin(is->rctx, y, is->y, INSERT_VALUES, SCATTER_FORWARD));
2348:   PetscCall(VecScatterEnd(is->rctx, y, is->y, INSERT_VALUES, SCATTER_FORWARD));

2350:   /* multiply the local matrix */
2351:   PetscCall(MatMultTranspose(is->A, is->y, is->x));

2353:   /* scatter product back into global vector */
2354:   PetscCall(VecSet(x, 0));
2355:   PetscCall(VecScatterBegin(is->cctx, is->x, x, ADD_VALUES, SCATTER_REVERSE));
2356:   PetscCall(VecScatterEnd(is->cctx, is->x, x, ADD_VALUES, SCATTER_REVERSE));
2357:   PetscFunctionReturn(PETSC_SUCCESS);
2358: }

2360: static PetscErrorCode MatMultTransposeAdd_IS(Mat A, Vec v1, Vec v2, Vec v3)
2361: {
2362:   Vec temp_vec;

2364:   PetscFunctionBegin; /*  v3 = v2 + A' * v1.*/
2365:   if (v3 != v2) {
2366:     PetscCall(MatMultTranspose(A, v1, v3));
2367:     PetscCall(VecAXPY(v3, 1.0, v2));
2368:   } else {
2369:     PetscCall(VecDuplicate(v2, &temp_vec));
2370:     PetscCall(MatMultTranspose(A, v1, temp_vec));
2371:     PetscCall(VecAXPY(temp_vec, 1.0, v2));
2372:     PetscCall(VecCopy(temp_vec, v3));
2373:     PetscCall(VecDestroy(&temp_vec));
2374:   }
2375:   PetscFunctionReturn(PETSC_SUCCESS);
2376: }

2378: static PetscErrorCode MatView_IS(Mat A, PetscViewer viewer)
2379: {
2380:   Mat_IS                *a = (Mat_IS *)A->data;
2381:   PetscViewer            sviewer;
2382:   PetscBool              isascii, isbinary, viewl2g = PETSC_FALSE, native;
2383:   PetscViewerFormat      format;
2384:   ISLocalToGlobalMapping rmap, cmap;

2386:   PetscFunctionBegin;
2387:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2388:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
2389:   PetscCall(PetscViewerGetFormat(viewer, &format));
2390:   native = (PetscBool)(format == PETSC_VIEWER_NATIVE);
2391:   if (native) {
2392:     rmap = A->rmap->mapping;
2393:     cmap = A->cmap->mapping;
2394:   } else {
2395:     rmap = a->rmapping;
2396:     cmap = a->cmapping;
2397:   }
2398:   if (isascii) {
2399:     if (format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(PETSC_SUCCESS);
2400:     if (format == PETSC_VIEWER_ASCII_INFO_DETAIL || format == PETSC_VIEWER_ASCII_MATLAB) viewl2g = PETSC_TRUE;
2401:   } else if (isbinary) {
2402:     PetscInt    tr[6], nr, nc;
2403:     char        lmattype[64] = {'\0'};
2404:     PetscMPIInt size;
2405:     PetscBool   skipHeader;
2406:     IS          is;

2408:     PetscCall(PetscViewerSetUp(viewer));
2409:     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)viewer), &size));
2410:     tr[0] = MAT_FILE_CLASSID;
2411:     tr[1] = A->rmap->N;
2412:     tr[2] = A->cmap->N;
2413:     tr[3] = -size; /* AIJ stores nnz here */
2414:     tr[4] = (PetscInt)(rmap == cmap);
2415:     tr[5] = a->allow_repeated;
2416:     PetscCall(PetscSNPrintf(lmattype, sizeof(lmattype), "%s", a->lmattype));

2418:     PetscCall(PetscViewerBinaryWrite(viewer, tr, PETSC_STATIC_ARRAY_LENGTH(tr), PETSC_INT));
2419:     PetscCall(PetscViewerBinaryWrite(viewer, lmattype, sizeof(lmattype), PETSC_CHAR));

2421:     /* first dump l2g info (we need the header for proper loading on different number of processes) */
2422:     PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));
2423:     PetscCall(PetscViewerBinarySetSkipHeader(viewer, PETSC_FALSE));
2424:     PetscCall(ISLocalToGlobalMappingView(rmap, viewer));
2425:     if (cmap != rmap) PetscCall(ISLocalToGlobalMappingView(cmap, viewer));

2427:     /* then the sizes of the local matrices */
2428:     PetscCall(MatGetSize(a->A, &nr, &nc));
2429:     PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)viewer), 1, &nr, PETSC_USE_POINTER, &is));
2430:     PetscCall(ISView(is, viewer));
2431:     PetscCall(ISDestroy(&is));
2432:     PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)viewer), 1, &nc, PETSC_USE_POINTER, &is));
2433:     PetscCall(ISView(is, viewer));
2434:     PetscCall(ISDestroy(&is));
2435:     PetscCall(PetscViewerBinarySetSkipHeader(viewer, skipHeader));
2436:   }
2437:   if (format == PETSC_VIEWER_ASCII_MATLAB) {
2438:     char        name[64];
2439:     PetscMPIInt size, rank;

2441:     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)viewer), &size));
2442:     PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
2443:     if (size > 1) PetscCall(PetscSNPrintf(name, sizeof(name), "lmat_%d", rank));
2444:     else PetscCall(PetscSNPrintf(name, sizeof(name), "lmat"));
2445:     PetscCall(PetscObjectSetName((PetscObject)a->A, name));
2446:   }

2448:   /* Dump the local matrices */
2449:   if (isbinary) { /* ViewerGetSubViewer does not work in parallel */
2450:     PetscBool   isaij;
2451:     PetscInt    nr, nc;
2452:     Mat         lA, B;
2453:     Mat_MPIAIJ *b;

2455:     /* We create a temporary MPIAIJ matrix that stores the unassembled operator */
2456:     PetscCall(PetscObjectBaseTypeCompare((PetscObject)a->A, MATAIJ, &isaij));
2457:     if (!isaij) PetscCall(MatConvert(a->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &lA));
2458:     else {
2459:       PetscCall(PetscObjectReference((PetscObject)a->A));
2460:       lA = a->A;
2461:     }
2462:     PetscCall(MatCreate(PetscObjectComm((PetscObject)viewer), &B));
2463:     PetscCall(MatSetType(B, MATMPIAIJ));
2464:     PetscCall(MatGetSize(lA, &nr, &nc));
2465:     PetscCall(MatSetSizes(B, nr, nc, PETSC_DECIDE, PETSC_DECIDE));
2466:     PetscCall(MatMPIAIJSetPreallocation(B, 0, NULL, 0, NULL));

2468:     b = (Mat_MPIAIJ *)B->data;
2469:     PetscCall(MatDestroy(&b->A));
2470:     b->A = lA;

2472:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
2473:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
2474:     PetscCall(MatView(B, viewer));
2475:     PetscCall(MatDestroy(&B));
2476:   } else {
2477:     PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
2478:     PetscCall(MatView(a->A, sviewer));
2479:     PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
2480:   }

2482:   /* with ASCII, we dump the l2gmaps at the end */
2483:   if (viewl2g) {
2484:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
2485:       PetscCall(PetscObjectSetName((PetscObject)rmap, "row"));
2486:       PetscCall(ISLocalToGlobalMappingView(rmap, viewer));
2487:       PetscCall(PetscObjectSetName((PetscObject)cmap, "col"));
2488:       PetscCall(ISLocalToGlobalMappingView(cmap, viewer));
2489:     } else {
2490:       PetscCall(ISLocalToGlobalMappingView(rmap, viewer));
2491:       PetscCall(ISLocalToGlobalMappingView(cmap, viewer));
2492:     }
2493:   }
2494:   PetscFunctionReturn(PETSC_SUCCESS);
2495: }

2497: static PetscErrorCode MatLoad_IS(Mat A, PetscViewer viewer)
2498: {
2499:   ISLocalToGlobalMapping rmap, cmap;
2500:   MPI_Comm               comm = PetscObjectComm((PetscObject)A);
2501:   PetscBool              isbinary, samel, allow, isbaij;
2502:   PetscInt               tr[6], M, N, nr, nc, Asize, isn;
2503:   const PetscInt        *idx;
2504:   PetscMPIInt            size;
2505:   char                   lmattype[64];
2506:   Mat                    dA, lA;
2507:   IS                     is;

2509:   PetscFunctionBegin;
2510:   PetscCheckSameComm(A, 1, viewer, 2);
2511:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
2512:   PetscCheck(isbinary, PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Invalid viewer of type %s", ((PetscObject)viewer)->type_name);

2514:   PetscCall(PetscViewerBinaryRead(viewer, tr, PETSC_STATIC_ARRAY_LENGTH(tr), NULL, PETSC_INT));
2515:   PetscCheck(tr[0] == MAT_FILE_CLASSID, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a matrix next in file");
2516:   PetscCheck(tr[1] >= 0, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a IS matrix next in file");
2517:   PetscCheck(tr[2] >= 0, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a IS matrix next in file");
2518:   PetscCheck(tr[3] < 0, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a IS matrix next in file");
2519:   PetscCheck(tr[4] == 0 || tr[4] == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a IS matrix next in file");
2520:   PetscCheck(tr[5] == 0 || tr[5] == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a IS matrix next in file");
2521:   M     = tr[1];
2522:   N     = tr[2];
2523:   Asize = -tr[3];
2524:   samel = (PetscBool)tr[4];
2525:   allow = (PetscBool)tr[5];
2526:   PetscCall(PetscViewerBinaryRead(viewer, lmattype, sizeof(lmattype), NULL, PETSC_CHAR));

2528:   /* if we are loading from a larger set of processes, allow repeated entries */
2529:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)viewer), &size));
2530:   if (Asize > size) allow = PETSC_TRUE;

2532:   /* set global sizes if not set already */
2533:   if (A->rmap->N < 0) A->rmap->N = M;
2534:   if (A->cmap->N < 0) A->cmap->N = N;
2535:   PetscCall(PetscLayoutSetUp(A->rmap));
2536:   PetscCall(PetscLayoutSetUp(A->cmap));
2537:   PetscCheck(M == A->rmap->N, comm, PETSC_ERR_ARG_SIZ, "Matrix rows should be %" PetscInt_FMT ", found %" PetscInt_FMT, M, A->rmap->N);
2538:   PetscCheck(N == A->cmap->N, comm, PETSC_ERR_ARG_SIZ, "Matrix columns should be %" PetscInt_FMT ", found %" PetscInt_FMT, N, A->cmap->N);

2540:   /* load l2g maps */
2541:   PetscCall(ISLocalToGlobalMappingCreate(comm, 0, 0, NULL, PETSC_USE_POINTER, &rmap));
2542:   PetscCall(ISLocalToGlobalMappingLoad(rmap, viewer));
2543:   if (!samel) {
2544:     PetscCall(ISLocalToGlobalMappingCreate(comm, 0, 0, NULL, PETSC_USE_POINTER, &cmap));
2545:     PetscCall(ISLocalToGlobalMappingLoad(cmap, viewer));
2546:   } else {
2547:     PetscCall(PetscObjectReference((PetscObject)rmap));
2548:     cmap = rmap;
2549:   }

2551:   /* load sizes of local matrices */
2552:   PetscCall(ISCreate(comm, &is));
2553:   PetscCall(ISSetType(is, ISGENERAL));
2554:   PetscCall(ISLoad(is, viewer));
2555:   PetscCall(ISGetLocalSize(is, &isn));
2556:   PetscCall(ISGetIndices(is, &idx));
2557:   nr = 0;
2558:   for (PetscInt i = 0; i < isn; i++) nr += idx[i];
2559:   PetscCall(ISRestoreIndices(is, &idx));
2560:   PetscCall(ISDestroy(&is));
2561:   PetscCall(ISCreate(comm, &is));
2562:   PetscCall(ISSetType(is, ISGENERAL));
2563:   PetscCall(ISLoad(is, viewer));
2564:   PetscCall(ISGetLocalSize(is, &isn));
2565:   PetscCall(ISGetIndices(is, &idx));
2566:   nc = 0;
2567:   for (PetscInt i = 0; i < isn; i++) nc += idx[i];
2568:   PetscCall(ISRestoreIndices(is, &idx));
2569:   PetscCall(ISDestroy(&is));

2571:   /* now load the unassembled operator */
2572:   PetscCall(MatCreate(comm, &dA));
2573:   PetscCall(MatSetType(dA, MATMPIAIJ));
2574:   PetscCall(MatSetSizes(dA, nr, nc, PETSC_DECIDE, PETSC_DECIDE));
2575:   PetscCall(MatLoad(dA, viewer));
2576:   PetscCall(MatMPIAIJGetSeqAIJ(dA, &lA, NULL, NULL));
2577:   PetscCall(PetscObjectReference((PetscObject)lA));
2578:   PetscCall(MatDestroy(&dA));

2580:   /* and convert to the desired format */
2581:   PetscCall(PetscStrcmpAny(lmattype, &isbaij, MATSBAIJ, MATSEQSBAIJ, ""));
2582:   if (isbaij) PetscCall(MatSetOption(lA, MAT_SYMMETRIC, PETSC_TRUE));
2583:   PetscCall(MatConvert(lA, lmattype, MAT_INPLACE_MATRIX, &lA));

2585:   /* assemble the MATIS object */
2586:   PetscCall(MatISSetAllowRepeated(A, allow));
2587:   PetscCall(MatSetLocalToGlobalMapping(A, rmap, cmap));
2588:   PetscCall(MatISSetLocalMat(A, lA));
2589:   PetscCall(MatDestroy(&lA));
2590:   PetscCall(ISLocalToGlobalMappingDestroy(&rmap));
2591:   PetscCall(ISLocalToGlobalMappingDestroy(&cmap));
2592:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
2593:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
2594:   PetscFunctionReturn(PETSC_SUCCESS);
2595: }

2597: static PetscErrorCode MatInvertBlockDiagonal_IS(Mat mat, const PetscScalar **values)
2598: {
2599:   Mat_IS            *is = (Mat_IS *)mat->data;
2600:   MPI_Datatype       nodeType;
2601:   const PetscScalar *lv;
2602:   PetscInt           bs;

2604:   PetscFunctionBegin;
2605:   PetscCall(MatGetBlockSize(mat, &bs));
2606:   PetscCall(MatSetBlockSize(is->A, bs));
2607:   PetscCall(MatInvertBlockDiagonal(is->A, &lv));
2608:   if (!is->bdiag) PetscCall(PetscMalloc1(bs * mat->rmap->n, &is->bdiag));
2609:   PetscCallMPI(MPI_Type_contiguous(bs, MPIU_SCALAR, &nodeType));
2610:   PetscCallMPI(MPI_Type_commit(&nodeType));
2611:   PetscCall(PetscSFReduceBegin(is->sf, nodeType, lv, is->bdiag, MPI_REPLACE));
2612:   PetscCall(PetscSFReduceEnd(is->sf, nodeType, lv, is->bdiag, MPI_REPLACE));
2613:   PetscCallMPI(MPI_Type_free(&nodeType));
2614:   if (values) *values = is->bdiag;
2615:   PetscFunctionReturn(PETSC_SUCCESS);
2616: }

2618: static PetscErrorCode MatISSetUpScatters_Private(Mat A)
2619: {
2620:   Vec             cglobal, rglobal;
2621:   IS              from;
2622:   Mat_IS         *is = (Mat_IS *)A->data;
2623:   PetscScalar     sum;
2624:   const PetscInt *garray;
2625:   PetscInt        nr, rbs, nc, cbs;
2626:   VecType         rtype;

2628:   PetscFunctionBegin;
2629:   PetscCall(ISLocalToGlobalMappingGetSize(is->rmapping, &nr));
2630:   PetscCall(ISLocalToGlobalMappingGetBlockSize(is->rmapping, &rbs));
2631:   PetscCall(ISLocalToGlobalMappingGetSize(is->cmapping, &nc));
2632:   PetscCall(ISLocalToGlobalMappingGetBlockSize(is->cmapping, &cbs));
2633:   PetscCall(VecDestroy(&is->x));
2634:   PetscCall(VecDestroy(&is->y));
2635:   PetscCall(VecDestroy(&is->counter));
2636:   PetscCall(VecScatterDestroy(&is->rctx));
2637:   PetscCall(VecScatterDestroy(&is->cctx));
2638:   PetscCall(MatCreateVecs(is->A, &is->x, &is->y));
2639:   PetscCall(VecBindToCPU(is->y, PETSC_TRUE));
2640:   PetscCall(VecGetRootType_Private(is->y, &rtype));
2641:   PetscCall(PetscFree(A->defaultvectype));
2642:   PetscCall(PetscStrallocpy(rtype, &A->defaultvectype));
2643:   PetscCall(MatCreateVecs(A, &cglobal, &rglobal));
2644:   PetscCall(VecBindToCPU(rglobal, PETSC_TRUE));
2645:   PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->rmapping, &garray));
2646:   PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), rbs, nr / rbs, garray, PETSC_USE_POINTER, &from));
2647:   PetscCall(VecScatterCreate(rglobal, from, is->y, NULL, &is->rctx));
2648:   PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->rmapping, &garray));
2649:   PetscCall(ISDestroy(&from));
2650:   if (is->rmapping != is->cmapping) {
2651:     PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->cmapping, &garray));
2652:     PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), cbs, nc / cbs, garray, PETSC_USE_POINTER, &from));
2653:     PetscCall(VecScatterCreate(cglobal, from, is->x, NULL, &is->cctx));
2654:     PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->cmapping, &garray));
2655:     PetscCall(ISDestroy(&from));
2656:   } else {
2657:     PetscCall(PetscObjectReference((PetscObject)is->rctx));
2658:     is->cctx = is->rctx;
2659:   }
2660:   PetscCall(VecDestroy(&cglobal));

2662:   /* interface counter vector (local) */
2663:   PetscCall(VecDuplicate(is->y, &is->counter));
2664:   PetscCall(VecBindToCPU(is->counter, PETSC_TRUE));
2665:   PetscCall(VecSet(is->y, 1.));
2666:   PetscCall(VecScatterBegin(is->rctx, is->y, rglobal, ADD_VALUES, SCATTER_REVERSE));
2667:   PetscCall(VecScatterEnd(is->rctx, is->y, rglobal, ADD_VALUES, SCATTER_REVERSE));
2668:   PetscCall(VecScatterBegin(is->rctx, rglobal, is->counter, INSERT_VALUES, SCATTER_FORWARD));
2669:   PetscCall(VecScatterEnd(is->rctx, rglobal, is->counter, INSERT_VALUES, SCATTER_FORWARD));
2670:   PetscCall(VecBindToCPU(is->y, PETSC_FALSE));
2671:   PetscCall(VecBindToCPU(is->counter, PETSC_FALSE));

2673:   /* special functions for block-diagonal matrices */
2674:   PetscCall(VecSum(rglobal, &sum));
2675:   A->ops->invertblockdiagonal = NULL;
2676:   if ((PetscInt)(PetscRealPart(sum)) == A->rmap->N && A->rmap->N == A->cmap->N && is->rmapping == is->cmapping) A->ops->invertblockdiagonal = MatInvertBlockDiagonal_IS;
2677:   PetscCall(VecDestroy(&rglobal));

2679:   /* setup SF for general purpose shared indices based communications */
2680:   PetscCall(MatISSetUpSF_IS(A));
2681:   PetscFunctionReturn(PETSC_SUCCESS);
2682: }

2684: static PetscErrorCode MatISFilterL2GMap(Mat A, ISLocalToGlobalMapping map, ISLocalToGlobalMapping *nmap, ISLocalToGlobalMapping *lmap)
2685: {
2686:   Mat_IS                    *matis = (Mat_IS *)A->data;
2687:   IS                         is;
2688:   ISLocalToGlobalMappingType l2gtype;
2689:   const PetscInt            *idxs;
2690:   PetscHSetI                 ht;
2691:   PetscInt                  *nidxs;
2692:   PetscInt                   i, n, bs, c;
2693:   PetscBool                  flg[] = {PETSC_FALSE, PETSC_FALSE};

2695:   PetscFunctionBegin;
2696:   PetscCall(ISLocalToGlobalMappingGetSize(map, &n));
2697:   PetscCall(ISLocalToGlobalMappingGetBlockSize(map, &bs));
2698:   PetscCall(ISLocalToGlobalMappingGetBlockIndices(map, &idxs));
2699:   PetscCall(PetscHSetICreate(&ht));
2700:   PetscCall(PetscMalloc1(n / bs, &nidxs));
2701:   for (i = 0, c = 0; i < n / bs; i++) {
2702:     PetscBool missing = PETSC_TRUE;
2703:     if (idxs[i] < 0) {
2704:       flg[0] = PETSC_TRUE;
2705:       continue;
2706:     }
2707:     if (!matis->allow_repeated) PetscCall(PetscHSetIQueryAdd(ht, idxs[i], &missing));
2708:     if (!missing) flg[1] = PETSC_TRUE;
2709:     else nidxs[c++] = idxs[i];
2710:   }
2711:   PetscCall(PetscHSetIDestroy(&ht));
2712:   PetscCall(MPIU_Allreduce(MPI_IN_PLACE, flg, 2, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A)));
2713:   if (!flg[0] && !flg[1]) { /* Entries are all non negative and unique */
2714:     *nmap = NULL;
2715:     *lmap = NULL;
2716:     PetscCall(PetscFree(nidxs));
2717:     PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(map, &idxs));
2718:     PetscFunctionReturn(PETSC_SUCCESS);
2719:   }

2721:   /* New l2g map without negative indices (and repeated indices if not allowed) */
2722:   PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), bs, c, nidxs, PETSC_USE_POINTER, &is));
2723:   PetscCall(ISLocalToGlobalMappingCreateIS(is, nmap));
2724:   PetscCall(ISDestroy(&is));
2725:   PetscCall(ISLocalToGlobalMappingGetType(map, &l2gtype));
2726:   PetscCall(ISLocalToGlobalMappingSetType(*nmap, l2gtype));

2728:   /* New local l2g map for repeated indices if not allowed */
2729:   PetscCall(ISGlobalToLocalMappingApplyBlock(*nmap, IS_GTOLM_MASK, n / bs, idxs, NULL, nidxs));
2730:   PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, n / bs, nidxs, PETSC_USE_POINTER, &is));
2731:   PetscCall(ISLocalToGlobalMappingCreateIS(is, lmap));
2732:   PetscCall(ISDestroy(&is));
2733:   PetscCall(PetscFree(nidxs));
2734:   PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(map, &idxs));
2735:   PetscFunctionReturn(PETSC_SUCCESS);
2736: }

2738: static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A, ISLocalToGlobalMapping rmapping, ISLocalToGlobalMapping cmapping)
2739: {
2740:   Mat_IS                *is            = (Mat_IS *)A->data;
2741:   ISLocalToGlobalMapping localrmapping = NULL, localcmapping = NULL;
2742:   PetscInt               nr, rbs, nc, cbs;
2743:   PetscBool              cong, freem[] = {PETSC_FALSE, PETSC_FALSE};

2745:   PetscFunctionBegin;
2746:   if (rmapping) PetscCheckSameComm(A, 1, rmapping, 2);
2747:   if (cmapping) PetscCheckSameComm(A, 1, cmapping, 3);

2749:   PetscCall(ISLocalToGlobalMappingDestroy(&is->rmapping));
2750:   PetscCall(ISLocalToGlobalMappingDestroy(&is->cmapping));
2751:   PetscCall(PetscLayoutSetUp(A->rmap));
2752:   PetscCall(PetscLayoutSetUp(A->cmap));
2753:   PetscCall(MatHasCongruentLayouts(A, &cong));

2755:   /* If NULL, local space matches global space */
2756:   if (!rmapping) {
2757:     IS is;

2759:     PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), A->rmap->N, 0, 1, &is));
2760:     PetscCall(ISLocalToGlobalMappingCreateIS(is, &rmapping));
2761:     if (A->rmap->bs > 0) PetscCall(ISLocalToGlobalMappingSetBlockSize(rmapping, A->rmap->bs));
2762:     PetscCall(ISDestroy(&is));
2763:     freem[0] = PETSC_TRUE;
2764:     if (!cmapping && cong && A->rmap->bs == A->cmap->bs) cmapping = rmapping;
2765:   } else if (!is->islocalref) { /* check if the l2g map has negative or repeated entries */
2766:     PetscCall(MatISFilterL2GMap(A, rmapping, &is->rmapping, &localrmapping));
2767:     if (rmapping == cmapping) {
2768:       PetscCall(PetscObjectReference((PetscObject)is->rmapping));
2769:       is->cmapping = is->rmapping;
2770:       PetscCall(PetscObjectReference((PetscObject)localrmapping));
2771:       localcmapping = localrmapping;
2772:     }
2773:   }
2774:   if (!cmapping) {
2775:     IS is;

2777:     PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), A->cmap->N, 0, 1, &is));
2778:     PetscCall(ISLocalToGlobalMappingCreateIS(is, &cmapping));
2779:     if (A->cmap->bs > 0) PetscCall(ISLocalToGlobalMappingSetBlockSize(cmapping, A->cmap->bs));
2780:     PetscCall(ISDestroy(&is));
2781:     freem[1] = PETSC_TRUE;
2782:   } else if (cmapping != rmapping && !is->islocalref) { /* check if the l2g map has negative or repeated entries */
2783:     PetscCall(MatISFilterL2GMap(A, cmapping, &is->cmapping, &localcmapping));
2784:   }
2785:   if (!is->rmapping) {
2786:     PetscCall(PetscObjectReference((PetscObject)rmapping));
2787:     is->rmapping = rmapping;
2788:   }
2789:   if (!is->cmapping) {
2790:     PetscCall(PetscObjectReference((PetscObject)cmapping));
2791:     is->cmapping = cmapping;
2792:   }

2794:   /* Clean up */
2795:   is->lnnzstate = 0;
2796:   PetscCall(MatDestroy(&is->dA));
2797:   PetscCall(MatDestroy(&is->assembledA));
2798:   PetscCall(MatDestroy(&is->A));
2799:   if (is->csf != is->sf) {
2800:     PetscCall(PetscSFDestroy(&is->csf));
2801:     PetscCall(PetscFree2(is->csf_rootdata, is->csf_leafdata));
2802:   } else is->csf = NULL;
2803:   PetscCall(PetscSFDestroy(&is->sf));
2804:   PetscCall(PetscFree2(is->sf_rootdata, is->sf_leafdata));
2805:   PetscCall(PetscFree(is->bdiag));

2807:   /* check if the two mappings are actually the same for square matrices since MATIS has some optimization for this case
2808:      (DOLFIN passes 2 different objects) */
2809:   PetscCall(ISLocalToGlobalMappingGetSize(is->rmapping, &nr));
2810:   PetscCall(ISLocalToGlobalMappingGetBlockSize(is->rmapping, &rbs));
2811:   PetscCall(ISLocalToGlobalMappingGetSize(is->cmapping, &nc));
2812:   PetscCall(ISLocalToGlobalMappingGetBlockSize(is->cmapping, &cbs));
2813:   if (is->rmapping != is->cmapping && cong) {
2814:     PetscBool same = PETSC_FALSE;
2815:     if (nr == nc && cbs == rbs) {
2816:       const PetscInt *idxs1, *idxs2;

2818:       PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->rmapping, &idxs1));
2819:       PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->cmapping, &idxs2));
2820:       PetscCall(PetscArraycmp(idxs1, idxs2, nr / rbs, &same));
2821:       PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->rmapping, &idxs1));
2822:       PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->cmapping, &idxs2));
2823:     }
2824:     PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &same, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
2825:     if (same) {
2826:       PetscCall(ISLocalToGlobalMappingDestroy(&is->cmapping));
2827:       PetscCall(PetscObjectReference((PetscObject)is->rmapping));
2828:       is->cmapping = is->rmapping;
2829:     }
2830:   }
2831:   PetscCall(PetscLayoutSetBlockSize(A->rmap, rbs));
2832:   PetscCall(PetscLayoutSetBlockSize(A->cmap, cbs));
2833:   /* Pass the user defined maps to the layout */
2834:   PetscCall(PetscLayoutSetISLocalToGlobalMapping(A->rmap, rmapping));
2835:   PetscCall(PetscLayoutSetISLocalToGlobalMapping(A->cmap, cmapping));
2836:   if (freem[0]) PetscCall(ISLocalToGlobalMappingDestroy(&rmapping));
2837:   if (freem[1]) PetscCall(ISLocalToGlobalMappingDestroy(&cmapping));

2839:   /* Create the local matrix A */
2840:   PetscCall(MatCreate(PETSC_COMM_SELF, &is->A));
2841:   PetscCall(MatSetType(is->A, is->lmattype));
2842:   PetscCall(MatSetSizes(is->A, nr, nc, nr, nc));
2843:   PetscCall(MatSetBlockSizes(is->A, rbs, cbs));
2844:   PetscCall(MatSetOptionsPrefix(is->A, "is_"));
2845:   PetscCall(MatAppendOptionsPrefix(is->A, ((PetscObject)A)->prefix));
2846:   PetscCall(PetscLayoutSetUp(is->A->rmap));
2847:   PetscCall(PetscLayoutSetUp(is->A->cmap));
2848:   PetscCall(MatSetLocalToGlobalMapping(is->A, localrmapping, localcmapping));
2849:   PetscCall(ISLocalToGlobalMappingDestroy(&localrmapping));
2850:   PetscCall(ISLocalToGlobalMappingDestroy(&localcmapping));

2852:   /* setup scatters and local vectors for MatMult */
2853:   if (!is->islocalref) PetscCall(MatISSetUpScatters_Private(A));
2854:   A->preallocated = PETSC_TRUE;
2855:   PetscFunctionReturn(PETSC_SUCCESS);
2856: }

2858: static PetscErrorCode MatSetUp_IS(Mat A)
2859: {
2860:   ISLocalToGlobalMapping rmap, cmap;

2862:   PetscFunctionBegin;
2863:   PetscCall(MatGetLocalToGlobalMapping(A, &rmap, &cmap));
2864:   if (!rmap && !cmap) PetscCall(MatSetLocalToGlobalMapping(A, NULL, NULL));
2865:   PetscFunctionReturn(PETSC_SUCCESS);
2866: }

2868: static PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
2869: {
2870:   Mat_IS  *is = (Mat_IS *)mat->data;
2871:   PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION], cols_l[MATIS_MAX_ENTRIES_INSERTION];

2873:   PetscFunctionBegin;
2874:   PetscCall(ISGlobalToLocalMappingApply(is->rmapping, IS_GTOLM_MASK, m, rows, &m, rows_l));
2875:   if (m != n || rows != cols || is->cmapping != is->rmapping) {
2876:     PetscCall(ISGlobalToLocalMappingApply(is->cmapping, IS_GTOLM_MASK, n, cols, &n, cols_l));
2877:     PetscCall(MatSetValues(is->A, m, rows_l, n, cols_l, values, addv));
2878:   } else {
2879:     PetscCall(MatSetValues(is->A, m, rows_l, m, rows_l, values, addv));
2880:   }
2881:   PetscFunctionReturn(PETSC_SUCCESS);
2882: }

2884: static PetscErrorCode MatSetValuesBlocked_IS(Mat mat, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
2885: {
2886:   Mat_IS  *is = (Mat_IS *)mat->data;
2887:   PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION], cols_l[MATIS_MAX_ENTRIES_INSERTION];

2889:   PetscFunctionBegin;
2890:   PetscCall(ISGlobalToLocalMappingApplyBlock(is->rmapping, IS_GTOLM_MASK, m, rows, &m, rows_l));
2891:   if (m != n || rows != cols || is->cmapping != is->rmapping) {
2892:     PetscCall(ISGlobalToLocalMappingApply(is->cmapping, IS_GTOLM_MASK, n, cols, &n, cols_l));
2893:     PetscCall(MatSetValuesBlocked(is->A, m, rows_l, n, cols_l, values, addv));
2894:   } else {
2895:     PetscCall(MatSetValuesBlocked(is->A, m, rows_l, n, rows_l, values, addv));
2896:   }
2897:   PetscFunctionReturn(PETSC_SUCCESS);
2898: }

2900: static PetscErrorCode MatSetValuesLocal_IS(Mat A, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
2901: {
2902:   Mat_IS *is = (Mat_IS *)A->data;

2904:   PetscFunctionBegin;
2905:   if (is->A->rmap->mapping || is->A->cmap->mapping) {
2906:     PetscCall(MatSetValuesLocal(is->A, m, rows, n, cols, values, addv));
2907:   } else {
2908:     PetscCall(MatSetValues(is->A, m, rows, n, cols, values, addv));
2909:   }
2910:   PetscFunctionReturn(PETSC_SUCCESS);
2911: }

2913: static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat A, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
2914: {
2915:   Mat_IS *is = (Mat_IS *)A->data;

2917:   PetscFunctionBegin;
2918:   if (is->A->rmap->mapping || is->A->cmap->mapping) {
2919:     PetscCall(MatSetValuesBlockedLocal(is->A, m, rows, n, cols, values, addv));
2920:   } else {
2921:     PetscCall(MatSetValuesBlocked(is->A, m, rows, n, cols, values, addv));
2922:   }
2923:   PetscFunctionReturn(PETSC_SUCCESS);
2924: }

2926: static PetscErrorCode MatISZeroRowsColumnsLocal_Private(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, PetscBool columns)
2927: {
2928:   Mat_IS *is = (Mat_IS *)A->data;

2930:   PetscFunctionBegin;
2931:   if (!n) {
2932:     is->pure_neumann = PETSC_TRUE;
2933:   } else {
2934:     PetscInt i;
2935:     is->pure_neumann = PETSC_FALSE;

2937:     if (columns) {
2938:       PetscCall(MatZeroRowsColumns(is->A, n, rows, diag, NULL, NULL));
2939:     } else {
2940:       PetscCall(MatZeroRows(is->A, n, rows, diag, NULL, NULL));
2941:     }
2942:     if (diag != 0.) {
2943:       const PetscScalar *array;
2944:       PetscCall(VecGetArrayRead(is->counter, &array));
2945:       for (i = 0; i < n; i++) PetscCall(MatSetValue(is->A, rows[i], rows[i], diag / (array[rows[i]]), INSERT_VALUES));
2946:       PetscCall(VecRestoreArrayRead(is->counter, &array));
2947:     }
2948:     PetscCall(MatAssemblyBegin(is->A, MAT_FINAL_ASSEMBLY));
2949:     PetscCall(MatAssemblyEnd(is->A, MAT_FINAL_ASSEMBLY));
2950:   }
2951:   PetscFunctionReturn(PETSC_SUCCESS);
2952: }

2954: static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b, PetscBool columns)
2955: {
2956:   Mat_IS   *matis = (Mat_IS *)A->data;
2957:   PetscInt  nr, nl, len, i;
2958:   PetscInt *lrows;

2960:   PetscFunctionBegin;
2961:   if (PetscUnlikelyDebug(columns || diag != 0. || (x && b))) {
2962:     PetscBool cong;

2964:     PetscCall(PetscLayoutCompare(A->rmap, A->cmap, &cong));
2965:     cong = (PetscBool)(cong && matis->sf == matis->csf);
2966:     PetscCheck(cong || !columns, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Columns can be zeroed if and only if A->rmap and A->cmap are congruent and the l2g maps are the same for MATIS");
2967:     PetscCheck(cong || diag == 0., PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Nonzero diagonal value supported if and only if A->rmap and A->cmap are congruent and the l2g maps are the same for MATIS");
2968:     PetscCheck(cong || !x || !b, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "A->rmap and A->cmap need to be congruent, and the l2g maps be the same");
2969:   }
2970:   /* get locally owned rows */
2971:   PetscCall(PetscLayoutMapLocal(A->rmap, n, rows, &len, &lrows, NULL));
2972:   /* fix right hand side if needed */
2973:   if (x && b) {
2974:     const PetscScalar *xx;
2975:     PetscScalar       *bb;

2977:     PetscCall(VecGetArrayRead(x, &xx));
2978:     PetscCall(VecGetArray(b, &bb));
2979:     for (i = 0; i < len; ++i) bb[lrows[i]] = diag * xx[lrows[i]];
2980:     PetscCall(VecRestoreArrayRead(x, &xx));
2981:     PetscCall(VecRestoreArray(b, &bb));
2982:   }
2983:   /* get rows associated to the local matrices */
2984:   PetscCall(MatGetSize(matis->A, &nl, NULL));
2985:   PetscCall(PetscArrayzero(matis->sf_leafdata, nl));
2986:   PetscCall(PetscArrayzero(matis->sf_rootdata, A->rmap->n));
2987:   for (i = 0; i < len; i++) matis->sf_rootdata[lrows[i]] = 1;
2988:   PetscCall(PetscFree(lrows));
2989:   PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
2990:   PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
2991:   PetscCall(PetscMalloc1(nl, &lrows));
2992:   for (i = 0, nr = 0; i < nl; i++)
2993:     if (matis->sf_leafdata[i]) lrows[nr++] = i;
2994:   PetscCall(MatISZeroRowsColumnsLocal_Private(A, nr, lrows, diag, columns));
2995:   PetscCall(PetscFree(lrows));
2996:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
2997:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
2998:   PetscFunctionReturn(PETSC_SUCCESS);
2999: }

3001: static PetscErrorCode MatZeroRows_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
3002: {
3003:   PetscFunctionBegin;
3004:   PetscCall(MatZeroRowsColumns_Private_IS(A, n, rows, diag, x, b, PETSC_FALSE));
3005:   PetscFunctionReturn(PETSC_SUCCESS);
3006: }

3008: static PetscErrorCode MatZeroRowsColumns_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
3009: {
3010:   PetscFunctionBegin;
3011:   PetscCall(MatZeroRowsColumns_Private_IS(A, n, rows, diag, x, b, PETSC_TRUE));
3012:   PetscFunctionReturn(PETSC_SUCCESS);
3013: }

3015: static PetscErrorCode MatAssemblyBegin_IS(Mat A, MatAssemblyType type)
3016: {
3017:   Mat_IS *is = (Mat_IS *)A->data;

3019:   PetscFunctionBegin;
3020:   PetscCall(MatAssemblyBegin(is->A, type));
3021:   PetscFunctionReturn(PETSC_SUCCESS);
3022: }

3024: static PetscErrorCode MatAssemblyEnd_IS(Mat A, MatAssemblyType type)
3025: {
3026:   Mat_IS   *is = (Mat_IS *)A->data;
3027:   PetscBool lnnz;

3029:   PetscFunctionBegin;
3030:   PetscCall(MatAssemblyEnd(is->A, type));
3031:   /* fix for local empty rows/cols */
3032:   if (is->locempty && type == MAT_FINAL_ASSEMBLY) {
3033:     Mat                    newlA;
3034:     ISLocalToGlobalMapping rl2g, cl2g;
3035:     IS                     nzr, nzc;
3036:     PetscInt               nr, nc, nnzr, nnzc;
3037:     PetscBool              lnewl2g, newl2g;

3039:     PetscCall(MatGetSize(is->A, &nr, &nc));
3040:     PetscCall(MatFindNonzeroRowsOrCols_Basic(is->A, PETSC_FALSE, PETSC_SMALL, &nzr));
3041:     if (!nzr) PetscCall(ISCreateStride(PetscObjectComm((PetscObject)is->A), nr, 0, 1, &nzr));
3042:     PetscCall(MatFindNonzeroRowsOrCols_Basic(is->A, PETSC_TRUE, PETSC_SMALL, &nzc));
3043:     if (!nzc) PetscCall(ISCreateStride(PetscObjectComm((PetscObject)is->A), nc, 0, 1, &nzc));
3044:     PetscCall(ISGetSize(nzr, &nnzr));
3045:     PetscCall(ISGetSize(nzc, &nnzc));
3046:     if (nnzr != nr || nnzc != nc) { /* need new global l2g map */
3047:       lnewl2g = PETSC_TRUE;
3048:       PetscCall(MPIU_Allreduce(&lnewl2g, &newl2g, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A)));

3050:       /* extract valid submatrix */
3051:       PetscCall(MatCreateSubMatrix(is->A, nzr, nzc, MAT_INITIAL_MATRIX, &newlA));
3052:     } else { /* local matrix fully populated */
3053:       lnewl2g = PETSC_FALSE;
3054:       PetscCall(MPIU_Allreduce(&lnewl2g, &newl2g, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A)));
3055:       PetscCall(PetscObjectReference((PetscObject)is->A));
3056:       newlA = is->A;
3057:     }

3059:     /* attach new global l2g map if needed */
3060:     if (newl2g) {
3061:       IS              zr, zc;
3062:       const PetscInt *ridxs, *cidxs, *zridxs, *zcidxs;
3063:       PetscInt       *nidxs, i;

3065:       PetscCall(ISComplement(nzr, 0, nr, &zr));
3066:       PetscCall(ISComplement(nzc, 0, nc, &zc));
3067:       PetscCall(PetscMalloc1(PetscMax(nr, nc), &nidxs));
3068:       PetscCall(ISLocalToGlobalMappingGetIndices(is->rmapping, &ridxs));
3069:       PetscCall(ISLocalToGlobalMappingGetIndices(is->cmapping, &cidxs));
3070:       PetscCall(ISGetIndices(zr, &zridxs));
3071:       PetscCall(ISGetIndices(zc, &zcidxs));
3072:       PetscCall(ISGetLocalSize(zr, &nnzr));
3073:       PetscCall(ISGetLocalSize(zc, &nnzc));

3075:       PetscCall(PetscArraycpy(nidxs, ridxs, nr));
3076:       for (i = 0; i < nnzr; i++) nidxs[zridxs[i]] = -1;
3077:       PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A), 1, nr, nidxs, PETSC_COPY_VALUES, &rl2g));
3078:       PetscCall(PetscArraycpy(nidxs, cidxs, nc));
3079:       for (i = 0; i < nnzc; i++) nidxs[zcidxs[i]] = -1;
3080:       PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A), 1, nc, nidxs, PETSC_COPY_VALUES, &cl2g));

3082:       PetscCall(ISRestoreIndices(zr, &zridxs));
3083:       PetscCall(ISRestoreIndices(zc, &zcidxs));
3084:       PetscCall(ISLocalToGlobalMappingRestoreIndices(is->rmapping, &ridxs));
3085:       PetscCall(ISLocalToGlobalMappingRestoreIndices(is->cmapping, &cidxs));
3086:       PetscCall(ISDestroy(&nzr));
3087:       PetscCall(ISDestroy(&nzc));
3088:       PetscCall(ISDestroy(&zr));
3089:       PetscCall(ISDestroy(&zc));
3090:       PetscCall(PetscFree(nidxs));
3091:       PetscCall(MatSetLocalToGlobalMapping(A, rl2g, cl2g));
3092:       PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
3093:       PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
3094:     }
3095:     PetscCall(MatISSetLocalMat(A, newlA));
3096:     PetscCall(MatDestroy(&newlA));
3097:     PetscCall(ISDestroy(&nzr));
3098:     PetscCall(ISDestroy(&nzc));
3099:     is->locempty = PETSC_FALSE;
3100:   }
3101:   lnnz          = (PetscBool)(is->A->nonzerostate == is->lnnzstate);
3102:   is->lnnzstate = is->A->nonzerostate;
3103:   PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &lnnz, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
3104:   if (!lnnz) A->nonzerostate++;
3105:   PetscFunctionReturn(PETSC_SUCCESS);
3106: }

3108: static PetscErrorCode MatISGetLocalMat_IS(Mat mat, Mat *local)
3109: {
3110:   Mat_IS *is = (Mat_IS *)mat->data;

3112:   PetscFunctionBegin;
3113:   *local = is->A;
3114:   PetscFunctionReturn(PETSC_SUCCESS);
3115: }

3117: static PetscErrorCode MatISRestoreLocalMat_IS(Mat mat, Mat *local)
3118: {
3119:   PetscFunctionBegin;
3120:   *local = NULL;
3121:   PetscFunctionReturn(PETSC_SUCCESS);
3122: }

3124: /*@
3125:   MatISGetLocalMat - Gets the local matrix stored inside a `MATIS` matrix.

3127:   Not Collective.

3129:   Input Parameter:
3130: . mat - the matrix

3132:   Output Parameter:
3133: . local - the local matrix

3135:   Level: intermediate

3137:   Notes:
3138:   This can be called if you have precomputed the nonzero structure of the
3139:   matrix and want to provide it to the inner matrix object to improve the performance
3140:   of the `MatSetValues()` operation.

3142:   Call `MatISRestoreLocalMat()` when finished with the local matrix.

3144: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatISRestoreLocalMat()`
3145: @*/
3146: PetscErrorCode MatISGetLocalMat(Mat mat, Mat *local)
3147: {
3148:   PetscFunctionBegin;
3150:   PetscAssertPointer(local, 2);
3151:   PetscUseMethod(mat, "MatISGetLocalMat_C", (Mat, Mat *), (mat, local));
3152:   PetscFunctionReturn(PETSC_SUCCESS);
3153: }

3155: /*@
3156:   MatISRestoreLocalMat - Restores the local matrix obtained with `MatISGetLocalMat()`

3158:   Not Collective.

3160:   Input Parameters:
3161: + mat   - the matrix
3162: - local - the local matrix

3164:   Level: intermediate

3166: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatISGetLocalMat()`
3167: @*/
3168: PetscErrorCode MatISRestoreLocalMat(Mat mat, Mat *local)
3169: {
3170:   PetscFunctionBegin;
3172:   PetscAssertPointer(local, 2);
3173:   PetscUseMethod(mat, "MatISRestoreLocalMat_C", (Mat, Mat *), (mat, local));
3174:   PetscFunctionReturn(PETSC_SUCCESS);
3175: }

3177: static PetscErrorCode MatISSetLocalMatType_IS(Mat mat, MatType mtype)
3178: {
3179:   Mat_IS *is = (Mat_IS *)mat->data;

3181:   PetscFunctionBegin;
3182:   if (is->A) PetscCall(MatSetType(is->A, mtype));
3183:   PetscCall(PetscFree(is->lmattype));
3184:   PetscCall(PetscStrallocpy(mtype, &is->lmattype));
3185:   PetscFunctionReturn(PETSC_SUCCESS);
3186: }

3188: /*@C
3189:   MatISSetLocalMatType - Specifies the type of local matrix inside the `MATIS`

3191:   Logically Collective.

3193:   Input Parameters:
3194: + mat   - the matrix
3195: - mtype - the local matrix type

3197:   Level: intermediate

3199: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatSetType()`, `MatType`
3200: @*/
3201: PetscErrorCode MatISSetLocalMatType(Mat mat, MatType mtype)
3202: {
3203:   PetscFunctionBegin;
3205:   PetscUseMethod(mat, "MatISSetLocalMatType_C", (Mat, MatType), (mat, mtype));
3206:   PetscFunctionReturn(PETSC_SUCCESS);
3207: }

3209: static PetscErrorCode MatISSetLocalMat_IS(Mat mat, Mat local)
3210: {
3211:   Mat_IS   *is = (Mat_IS *)mat->data;
3212:   PetscInt  nrows, ncols, orows, ocols;
3213:   MatType   mtype, otype;
3214:   PetscBool sametype = PETSC_TRUE;

3216:   PetscFunctionBegin;
3217:   if (is->A && !is->islocalref) {
3218:     PetscCall(MatGetSize(is->A, &orows, &ocols));
3219:     PetscCall(MatGetSize(local, &nrows, &ncols));
3220:     PetscCheck(orows == nrows && ocols == ncols, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local MATIS matrix should be of size %" PetscInt_FMT "x%" PetscInt_FMT " (passed a %" PetscInt_FMT "x%" PetscInt_FMT " matrix)", orows, ocols, nrows, ncols);
3221:     PetscCall(MatGetType(local, &mtype));
3222:     PetscCall(MatGetType(is->A, &otype));
3223:     PetscCall(PetscStrcmp(mtype, otype, &sametype));
3224:   }
3225:   PetscCall(PetscObjectReference((PetscObject)local));
3226:   PetscCall(MatDestroy(&is->A));
3227:   is->A = local;
3228:   PetscCall(MatGetType(is->A, &mtype));
3229:   PetscCall(MatISSetLocalMatType(mat, mtype));
3230:   if (!sametype && !is->islocalref) PetscCall(MatISSetUpScatters_Private(mat));
3231:   is->lnnzstate = 0;
3232:   PetscFunctionReturn(PETSC_SUCCESS);
3233: }

3235: /*@
3236:   MatISSetLocalMat - Replace the local matrix stored inside a `MATIS` object.

3238:   Not Collective

3240:   Input Parameters:
3241: + mat   - the matrix
3242: - local - the local matrix

3244:   Level: intermediate

3246: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatISSetLocalMatType`, `MatISGetLocalMat()`
3247: @*/
3248: PetscErrorCode MatISSetLocalMat(Mat mat, Mat local)
3249: {
3250:   PetscFunctionBegin;
3253:   PetscUseMethod(mat, "MatISSetLocalMat_C", (Mat, Mat), (mat, local));
3254:   PetscFunctionReturn(PETSC_SUCCESS);
3255: }

3257: static PetscErrorCode MatZeroEntries_IS(Mat A)
3258: {
3259:   Mat_IS *a = (Mat_IS *)A->data;

3261:   PetscFunctionBegin;
3262:   PetscCall(MatZeroEntries(a->A));
3263:   PetscFunctionReturn(PETSC_SUCCESS);
3264: }

3266: static PetscErrorCode MatScale_IS(Mat A, PetscScalar a)
3267: {
3268:   Mat_IS *is = (Mat_IS *)A->data;

3270:   PetscFunctionBegin;
3271:   PetscCall(MatScale(is->A, a));
3272:   PetscFunctionReturn(PETSC_SUCCESS);
3273: }

3275: static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v)
3276: {
3277:   Mat_IS *is = (Mat_IS *)A->data;

3279:   PetscFunctionBegin;
3280:   /* get diagonal of the local matrix */
3281:   PetscCall(MatGetDiagonal(is->A, is->y));

3283:   /* scatter diagonal back into global vector */
3284:   PetscCall(VecSet(v, 0));
3285:   PetscCall(VecScatterBegin(is->rctx, is->y, v, ADD_VALUES, SCATTER_REVERSE));
3286:   PetscCall(VecScatterEnd(is->rctx, is->y, v, ADD_VALUES, SCATTER_REVERSE));
3287:   PetscFunctionReturn(PETSC_SUCCESS);
3288: }

3290: static PetscErrorCode MatSetOption_IS(Mat A, MatOption op, PetscBool flg)
3291: {
3292:   Mat_IS *a = (Mat_IS *)A->data;

3294:   PetscFunctionBegin;
3295:   PetscCall(MatSetOption(a->A, op, flg));
3296:   PetscFunctionReturn(PETSC_SUCCESS);
3297: }

3299: static PetscErrorCode MatAXPY_IS(Mat Y, PetscScalar a, Mat X, MatStructure str)
3300: {
3301:   Mat_IS *y = (Mat_IS *)Y->data;
3302:   Mat_IS *x;

3304:   PetscFunctionBegin;
3305:   if (PetscDefined(USE_DEBUG)) {
3306:     PetscBool ismatis;
3307:     PetscCall(PetscObjectTypeCompare((PetscObject)X, MATIS, &ismatis));
3308:     PetscCheck(ismatis, PetscObjectComm((PetscObject)Y), PETSC_ERR_SUP, "Cannot call MatAXPY(Y,a,X,str) with X not of type MATIS");
3309:   }
3310:   x = (Mat_IS *)X->data;
3311:   PetscCall(MatAXPY(y->A, a, x->A, str));
3312:   PetscFunctionReturn(PETSC_SUCCESS);
3313: }

3315: static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A, IS row, IS col, Mat *submat)
3316: {
3317:   Mat                    lA;
3318:   Mat_IS                *matis = (Mat_IS *)A->data;
3319:   ISLocalToGlobalMapping rl2g, cl2g;
3320:   IS                     is;
3321:   const PetscInt        *rg, *rl;
3322:   PetscInt               nrg;
3323:   PetscInt               N, M, nrl, i, *idxs;

3325:   PetscFunctionBegin;
3326:   PetscCall(ISLocalToGlobalMappingGetIndices(A->rmap->mapping, &rg));
3327:   PetscCall(ISGetLocalSize(row, &nrl));
3328:   PetscCall(ISGetIndices(row, &rl));
3329:   PetscCall(ISLocalToGlobalMappingGetSize(A->rmap->mapping, &nrg));
3330:   if (PetscDefined(USE_DEBUG)) {
3331:     for (i = 0; i < nrl; i++) PetscCheck(rl[i] < nrg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Local row index %" PetscInt_FMT " -> %" PetscInt_FMT " greater then maximum possible %" PetscInt_FMT, i, rl[i], nrg);
3332:   }
3333:   PetscCall(PetscMalloc1(nrg, &idxs));
3334:   /* map from [0,nrl) to row */
3335:   for (i = 0; i < nrl; i++) idxs[i] = rl[i];
3336:   for (i = nrl; i < nrg; i++) idxs[i] = -1;
3337:   PetscCall(ISRestoreIndices(row, &rl));
3338:   PetscCall(ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping, &rg));
3339:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A), nrg, idxs, PETSC_OWN_POINTER, &is));
3340:   PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g));
3341:   PetscCall(ISDestroy(&is));
3342:   /* compute new l2g map for columns */
3343:   if (col != row || matis->rmapping != matis->cmapping || matis->A->rmap->mapping != matis->A->cmap->mapping) {
3344:     const PetscInt *cg, *cl;
3345:     PetscInt        ncg;
3346:     PetscInt        ncl;

3348:     PetscCall(ISLocalToGlobalMappingGetIndices(A->cmap->mapping, &cg));
3349:     PetscCall(ISGetLocalSize(col, &ncl));
3350:     PetscCall(ISGetIndices(col, &cl));
3351:     PetscCall(ISLocalToGlobalMappingGetSize(A->cmap->mapping, &ncg));
3352:     if (PetscDefined(USE_DEBUG)) {
3353:       for (i = 0; i < ncl; i++) PetscCheck(cl[i] < ncg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Local column index %" PetscInt_FMT " -> %" PetscInt_FMT " greater then maximum possible %" PetscInt_FMT, i, cl[i], ncg);
3354:     }
3355:     PetscCall(PetscMalloc1(ncg, &idxs));
3356:     /* map from [0,ncl) to col */
3357:     for (i = 0; i < ncl; i++) idxs[i] = cl[i];
3358:     for (i = ncl; i < ncg; i++) idxs[i] = -1;
3359:     PetscCall(ISRestoreIndices(col, &cl));
3360:     PetscCall(ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping, &cg));
3361:     PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A), ncg, idxs, PETSC_OWN_POINTER, &is));
3362:     PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g));
3363:     PetscCall(ISDestroy(&is));
3364:   } else {
3365:     PetscCall(PetscObjectReference((PetscObject)rl2g));
3366:     cl2g = rl2g;
3367:   }
3368:   /* create the MATIS submatrix */
3369:   PetscCall(MatGetSize(A, &M, &N));
3370:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), submat));
3371:   PetscCall(MatSetSizes(*submat, PETSC_DECIDE, PETSC_DECIDE, M, N));
3372:   PetscCall(MatSetType(*submat, MATIS));
3373:   matis             = (Mat_IS *)((*submat)->data);
3374:   matis->islocalref = PETSC_TRUE;
3375:   PetscCall(MatSetLocalToGlobalMapping(*submat, rl2g, cl2g));
3376:   PetscCall(MatISGetLocalMat(A, &lA));
3377:   PetscCall(MatISSetLocalMat(*submat, lA));
3378:   PetscCall(MatSetUp(*submat));
3379:   PetscCall(MatAssemblyBegin(*submat, MAT_FINAL_ASSEMBLY));
3380:   PetscCall(MatAssemblyEnd(*submat, MAT_FINAL_ASSEMBLY));
3381:   PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
3382:   PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));

3384:   /* remove unsupported ops */
3385:   PetscCall(PetscMemzero((*submat)->ops, sizeof(struct _MatOps)));
3386:   (*submat)->ops->destroy               = MatDestroy_IS;
3387:   (*submat)->ops->setvalueslocal        = MatSetValuesLocal_SubMat_IS;
3388:   (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS;
3389:   (*submat)->ops->assemblybegin         = MatAssemblyBegin_IS;
3390:   (*submat)->ops->assemblyend           = MatAssemblyEnd_IS;
3391:   PetscFunctionReturn(PETSC_SUCCESS);
3392: }

3394: static PetscErrorCode MatSetFromOptions_IS(Mat A, PetscOptionItems *PetscOptionsObject)
3395: {
3396:   Mat_IS   *a = (Mat_IS *)A->data;
3397:   char      type[256];
3398:   PetscBool flg;

3400:   PetscFunctionBegin;
3401:   PetscOptionsHeadBegin(PetscOptionsObject, "MATIS options");
3402:   PetscCall(PetscOptionsDeprecated("-matis_keepassembled", "-mat_is_keepassembled", "3.21", NULL));
3403:   PetscCall(PetscOptionsDeprecated("-matis_fixempty", "-mat_is_fixempty", "3.21", NULL));
3404:   PetscCall(PetscOptionsDeprecated("-matis_storel2l", "-mat_is_storel2l", "3.21", NULL));
3405:   PetscCall(PetscOptionsDeprecated("-matis_localmat_type", "-mat_is_localmat_type", "3.21", NULL));
3406:   PetscCall(PetscOptionsBool("-mat_is_keepassembled", "Store an assembled version if needed", NULL, a->keepassembled, &a->keepassembled, NULL));
3407:   PetscCall(PetscOptionsBool("-mat_is_fixempty", "Fix local matrices in case of empty local rows/columns", "MatISFixLocalEmpty", a->locempty, &a->locempty, NULL));
3408:   PetscCall(PetscOptionsBool("-mat_is_storel2l", "Store local-to-local matrices generated from PtAP operations", "MatISStoreL2L", a->storel2l, &a->storel2l, NULL));
3409:   PetscCall(PetscOptionsBool("-mat_is_allow_repeated", "Allow local repeated entries", "MatISSetAllowRepeated", a->allow_repeated, &a->allow_repeated, NULL));
3410:   PetscCall(PetscOptionsFList("-mat_is_localmat_type", "Matrix type", "MatISSetLocalMatType", MatList, a->lmattype, type, 256, &flg));
3411:   if (flg) PetscCall(MatISSetLocalMatType(A, type));
3412:   if (a->A) PetscCall(MatSetFromOptions(a->A));
3413:   PetscOptionsHeadEnd();
3414:   PetscFunctionReturn(PETSC_SUCCESS);
3415: }

3417: /*@
3418:   MatCreateIS - Creates a "process" unassembled matrix.

3420:   Collective.

3422:   Input Parameters:
3423: + comm - MPI communicator that will share the matrix
3424: . bs   - block size of the matrix
3425: . m    - local size of left vector used in matrix vector products
3426: . n    - local size of right vector used in matrix vector products
3427: . M    - global size of left vector used in matrix vector products
3428: . N    - global size of right vector used in matrix vector products
3429: . rmap - local to global map for rows
3430: - cmap - local to global map for cols

3432:   Output Parameter:
3433: . A - the resulting matrix

3435:   Level: intermediate

3437:   Notes:
3438:   `m` and `n` are NOT related to the size of the map; they represent the size of the local parts of the distributed vectors
3439:   used in `MatMult()` operations. The local sizes of `rmap` and `cmap` define the size of the local matrices.

3441:   If `rmap` (`cmap`) is `NULL`, then the local row (column) spaces matches the global space.

3443: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatSetLocalToGlobalMapping()`
3444: @*/
3445: PetscErrorCode MatCreateIS(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt M, PetscInt N, ISLocalToGlobalMapping rmap, ISLocalToGlobalMapping cmap, Mat *A)
3446: {
3447:   PetscFunctionBegin;
3448:   PetscCall(MatCreate(comm, A));
3449:   PetscCall(MatSetSizes(*A, m, n, M, N));
3450:   if (bs > 0) PetscCall(MatSetBlockSize(*A, bs));
3451:   PetscCall(MatSetType(*A, MATIS));
3452:   PetscCall(MatSetLocalToGlobalMapping(*A, rmap, cmap));
3453:   PetscFunctionReturn(PETSC_SUCCESS);
3454: }

3456: static PetscErrorCode MatHasOperation_IS(Mat A, MatOperation op, PetscBool *has)
3457: {
3458:   Mat_IS      *a              = (Mat_IS *)A->data;
3459:   MatOperation tobefiltered[] = {MATOP_MULT_ADD, MATOP_MULT_TRANSPOSE_ADD, MATOP_GET_DIAGONAL_BLOCK, MATOP_INCREASE_OVERLAP};

3461:   PetscFunctionBegin;
3462:   *has = PETSC_FALSE;
3463:   if (!((void **)A->ops)[op] || !a->A) PetscFunctionReturn(PETSC_SUCCESS);
3464:   *has = PETSC_TRUE;
3465:   for (PetscInt i = 0; i < (PetscInt)PETSC_STATIC_ARRAY_LENGTH(tobefiltered); i++)
3466:     if (op == tobefiltered[i]) PetscFunctionReturn(PETSC_SUCCESS);
3467:   PetscCall(MatHasOperation(a->A, op, has));
3468:   PetscFunctionReturn(PETSC_SUCCESS);
3469: }

3471: static PetscErrorCode MatSetValuesCOO_IS(Mat A, const PetscScalar v[], InsertMode imode)
3472: {
3473:   Mat_IS *a = (Mat_IS *)A->data;

3475:   PetscFunctionBegin;
3476:   PetscCall(MatSetValuesCOO(a->A, v, imode));
3477:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
3478:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
3479:   PetscFunctionReturn(PETSC_SUCCESS);
3480: }

3482: static PetscErrorCode MatSetPreallocationCOOLocal_IS(Mat A, PetscCount ncoo, PetscInt coo_i[], PetscInt coo_j[])
3483: {
3484:   Mat_IS *a = (Mat_IS *)A->data;

3486:   PetscFunctionBegin;
3487:   PetscCheck(a->A, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to provide l2g map first via MatSetLocalToGlobalMapping");
3488:   if (a->A->rmap->mapping || a->A->cmap->mapping) {
3489:     PetscCall(MatSetPreallocationCOOLocal(a->A, ncoo, coo_i, coo_j));
3490:   } else {
3491:     PetscCall(MatSetPreallocationCOO(a->A, ncoo, coo_i, coo_j));
3492:   }
3493:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", MatSetValuesCOO_IS));
3494:   A->preallocated = PETSC_TRUE;
3495:   PetscFunctionReturn(PETSC_SUCCESS);
3496: }

3498: static PetscErrorCode MatSetPreallocationCOO_IS(Mat A, PetscCount ncoo, PetscInt coo_i[], PetscInt coo_j[])
3499: {
3500:   Mat_IS *a = (Mat_IS *)A->data;

3502:   PetscFunctionBegin;
3503:   PetscCheck(a->A, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to provide l2g map first via MatSetLocalToGlobalMapping");
3504:   PetscCheck(ncoo <= PETSC_MAX_INT, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "ncoo %" PetscCount_FMT " overflowed PetscInt; configure --with-64-bit-indices or request support", ncoo);
3505:   PetscCall(ISGlobalToLocalMappingApply(a->rmapping, IS_GTOLM_MASK, ncoo, coo_i, NULL, coo_i));
3506:   PetscCall(ISGlobalToLocalMappingApply(a->cmapping, IS_GTOLM_MASK, ncoo, coo_j, NULL, coo_j));
3507:   PetscCall(MatSetPreallocationCOO(a->A, ncoo, coo_i, coo_j));
3508:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", MatSetValuesCOO_IS));
3509:   A->preallocated = PETSC_TRUE;
3510:   PetscFunctionReturn(PETSC_SUCCESS);
3511: }

3513: static PetscErrorCode MatISGetAssembled_Private(Mat A, Mat *tA)
3514: {
3515:   Mat_IS          *a = (Mat_IS *)A->data;
3516:   PetscObjectState Astate, aAstate       = PETSC_MIN_INT;
3517:   PetscObjectState Annzstate, aAnnzstate = PETSC_MIN_INT;

3519:   PetscFunctionBegin;
3520:   PetscCall(PetscObjectStateGet((PetscObject)A, &Astate));
3521:   Annzstate = A->nonzerostate;
3522:   if (a->assembledA) {
3523:     PetscCall(PetscObjectStateGet((PetscObject)a->assembledA, &aAstate));
3524:     aAnnzstate = a->assembledA->nonzerostate;
3525:   }
3526:   if (aAnnzstate != Annzstate) PetscCall(MatDestroy(&a->assembledA));
3527:   if (Astate != aAstate || !a->assembledA) {
3528:     MatType     aAtype;
3529:     PetscMPIInt size;
3530:     PetscInt    rbs, cbs, bs;

3532:     /* the assembled form is used as temporary storage for parallel operations
3533:        like createsubmatrices and the like, do not waste device memory */
3534:     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
3535:     PetscCall(ISLocalToGlobalMappingGetBlockSize(a->cmapping, &cbs));
3536:     PetscCall(ISLocalToGlobalMappingGetBlockSize(a->rmapping, &rbs));
3537:     bs = rbs == cbs ? rbs : 1;
3538:     if (a->assembledA) PetscCall(MatGetType(a->assembledA, &aAtype));
3539:     else if (size > 1) aAtype = bs > 1 ? MATMPIBAIJ : MATMPIAIJ;
3540:     else aAtype = bs > 1 ? MATSEQBAIJ : MATSEQAIJ;

3542:     PetscCall(MatConvert(A, aAtype, a->assembledA ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX, &a->assembledA));
3543:     PetscCall(PetscObjectStateSet((PetscObject)a->assembledA, Astate));
3544:     a->assembledA->nonzerostate = Annzstate;
3545:   }
3546:   PetscCall(PetscObjectReference((PetscObject)a->assembledA));
3547:   *tA = a->assembledA;
3548:   if (!a->keepassembled) PetscCall(MatDestroy(&a->assembledA));
3549:   PetscFunctionReturn(PETSC_SUCCESS);
3550: }

3552: static PetscErrorCode MatISRestoreAssembled_Private(Mat A, Mat *tA)
3553: {
3554:   PetscFunctionBegin;
3555:   PetscCall(MatDestroy(tA));
3556:   *tA = NULL;
3557:   PetscFunctionReturn(PETSC_SUCCESS);
3558: }

3560: static PetscErrorCode MatGetDiagonalBlock_IS(Mat A, Mat *dA)
3561: {
3562:   Mat_IS          *a = (Mat_IS *)A->data;
3563:   PetscObjectState Astate, dAstate = PETSC_MIN_INT;

3565:   PetscFunctionBegin;
3566:   PetscCall(PetscObjectStateGet((PetscObject)A, &Astate));
3567:   if (a->dA) PetscCall(PetscObjectStateGet((PetscObject)a->dA, &dAstate));
3568:   if (Astate != dAstate) {
3569:     Mat     tA;
3570:     MatType ltype;

3572:     PetscCall(MatDestroy(&a->dA));
3573:     PetscCall(MatISGetAssembled_Private(A, &tA));
3574:     PetscCall(MatGetDiagonalBlock(tA, &a->dA));
3575:     PetscCall(MatPropagateSymmetryOptions(tA, a->dA));
3576:     PetscCall(MatGetType(a->A, &ltype));
3577:     PetscCall(MatConvert(a->dA, ltype, MAT_INPLACE_MATRIX, &a->dA));
3578:     PetscCall(PetscObjectReference((PetscObject)a->dA));
3579:     PetscCall(MatISRestoreAssembled_Private(A, &tA));
3580:     PetscCall(PetscObjectStateSet((PetscObject)a->dA, Astate));
3581:   }
3582:   *dA = a->dA;
3583:   PetscFunctionReturn(PETSC_SUCCESS);
3584: }

3586: static PetscErrorCode MatCreateSubMatrices_IS(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse reuse, Mat *submat[])
3587: {
3588:   Mat tA;

3590:   PetscFunctionBegin;
3591:   PetscCall(MatISGetAssembled_Private(A, &tA));
3592:   PetscCall(MatCreateSubMatrices(tA, n, irow, icol, reuse, submat));
3593:   /* MatCreateSubMatrices_MPIAIJ is a mess at the moment */
3594: #if 0
3595:   {
3596:     Mat_IS    *a = (Mat_IS*)A->data;
3597:     MatType   ltype;
3598:     VecType   vtype;
3599:     char      *flg;

3601:     PetscCall(MatGetType(a->A,&ltype));
3602:     PetscCall(MatGetVecType(a->A,&vtype));
3603:     PetscCall(PetscStrstr(vtype,"cuda",&flg));
3604:     if (!flg) PetscCall(PetscStrstr(vtype,"hip",&flg));
3605:     if (!flg) PetscCall(PetscStrstr(vtype,"kokkos",&flg));
3606:     if (flg) {
3607:       for (PetscInt i = 0; i < n; i++) {
3608:         Mat sA = (*submat)[i];

3610:         PetscCall(MatConvert(sA,ltype,MAT_INPLACE_MATRIX,&sA));
3611:         (*submat)[i] = sA;
3612:       }
3613:     }
3614:   }
3615: #endif
3616:   PetscCall(MatISRestoreAssembled_Private(A, &tA));
3617:   PetscFunctionReturn(PETSC_SUCCESS);
3618: }

3620: static PetscErrorCode MatIncreaseOverlap_IS(Mat A, PetscInt n, IS is[], PetscInt ov)
3621: {
3622:   Mat tA;

3624:   PetscFunctionBegin;
3625:   PetscCall(MatISGetAssembled_Private(A, &tA));
3626:   PetscCall(MatIncreaseOverlap(tA, n, is, ov));
3627:   PetscCall(MatISRestoreAssembled_Private(A, &tA));
3628:   PetscFunctionReturn(PETSC_SUCCESS);
3629: }

3631: /*@
3632:   MatISGetLocalToGlobalMapping - Gets the local-to-global numbering of the `MATIS` object

3634:   Not Collective

3636:   Input Parameter:
3637: . A - the matrix

3639:   Output Parameters:
3640: + rmapping - row mapping
3641: - cmapping - column mapping

3643:   Level: advanced

3645:   Note:
3646:   The returned map can be different from the one used to construct the `MATIS` object, since it will not contain negative or repeated indices.

3648: .seealso: [](ch_matrices), `Mat`, `MatIS`, `MatSetLocalToGlobalMapping()`
3649: @*/
3650: PetscErrorCode MatISGetLocalToGlobalMapping(Mat A, ISLocalToGlobalMapping *rmapping, ISLocalToGlobalMapping *cmapping)
3651: {
3652:   PetscFunctionBegin;
3655:   if (rmapping) PetscAssertPointer(rmapping, 2);
3656:   if (cmapping) PetscAssertPointer(cmapping, 3);
3657:   PetscUseMethod(A, "MatISGetLocalToGlobalMapping_C", (Mat, ISLocalToGlobalMapping *, ISLocalToGlobalMapping *), (A, rmapping, cmapping));
3658:   PetscFunctionReturn(PETSC_SUCCESS);
3659: }

3661: static PetscErrorCode MatISGetLocalToGlobalMapping_IS(Mat A, ISLocalToGlobalMapping *r, ISLocalToGlobalMapping *c)
3662: {
3663:   Mat_IS *a = (Mat_IS *)A->data;

3665:   PetscFunctionBegin;
3666:   if (r) *r = a->rmapping;
3667:   if (c) *c = a->cmapping;
3668:   PetscFunctionReturn(PETSC_SUCCESS);
3669: }

3671: /*MC
3672:   MATIS - MATIS = "is" - A matrix type to be used for non-overlapping domain decomposition methods (e.g. `PCBDDC` or `KSPFETIDP`).
3673:   This stores the matrices in globally unassembled form and the parallel matrix vector product is handled "implicitly".

3675:   Options Database Keys:
3676: + -mat_type is           - Set the matrix type to `MATIS`.
3677: . -mat_is_allow_repeated - Allow repeated entries in the local part of the local to global maps.
3678: . -mat_is_fixempty       - Fix local matrices in case of empty local rows/columns.
3679: - -mat_is_storel2l       - Store the local-to-local operators generated by the Galerkin process of `MatPtAP()`.

3681:   Level: intermediate

3683:   Notes:
3684:   Options prefix for the inner matrix are given by `-is_mat_xxx`

3686:   You must call `MatSetLocalToGlobalMapping()` before using this matrix type.

3688:   You can do matrix preallocation on the local matrix after you obtain it with
3689:   `MatISGetLocalMat()`; otherwise, you could use `MatISSetPreallocation()` or `MatXAIJSetPreallocation()`

3691: .seealso: [](ch_matrices), `Mat`, `MATIS`, `Mat`, `MatISGetLocalMat()`, `MatSetLocalToGlobalMapping()`, `MatISSetPreallocation()`, `MatCreateIS()`, `PCBDDC`, `KSPFETIDP`
3692: M*/
3693: PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A)
3694: {
3695:   Mat_IS *a;

3697:   PetscFunctionBegin;
3698:   PetscCall(PetscNew(&a));
3699:   PetscCall(PetscStrallocpy(MATAIJ, &a->lmattype));
3700:   A->data = (void *)a;

3702:   /* matrix ops */
3703:   PetscCall(PetscMemzero(A->ops, sizeof(struct _MatOps)));
3704:   A->ops->mult                    = MatMult_IS;
3705:   A->ops->multadd                 = MatMultAdd_IS;
3706:   A->ops->multtranspose           = MatMultTranspose_IS;
3707:   A->ops->multtransposeadd        = MatMultTransposeAdd_IS;
3708:   A->ops->destroy                 = MatDestroy_IS;
3709:   A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS;
3710:   A->ops->setvalues               = MatSetValues_IS;
3711:   A->ops->setvaluesblocked        = MatSetValuesBlocked_IS;
3712:   A->ops->setvalueslocal          = MatSetValuesLocal_IS;
3713:   A->ops->setvaluesblockedlocal   = MatSetValuesBlockedLocal_IS;
3714:   A->ops->zerorows                = MatZeroRows_IS;
3715:   A->ops->zerorowscolumns         = MatZeroRowsColumns_IS;
3716:   A->ops->assemblybegin           = MatAssemblyBegin_IS;
3717:   A->ops->assemblyend             = MatAssemblyEnd_IS;
3718:   A->ops->view                    = MatView_IS;
3719:   A->ops->load                    = MatLoad_IS;
3720:   A->ops->zeroentries             = MatZeroEntries_IS;
3721:   A->ops->scale                   = MatScale_IS;
3722:   A->ops->getdiagonal             = MatGetDiagonal_IS;
3723:   A->ops->setoption               = MatSetOption_IS;
3724:   A->ops->ishermitian             = MatIsHermitian_IS;
3725:   A->ops->issymmetric             = MatIsSymmetric_IS;
3726:   A->ops->isstructurallysymmetric = MatIsStructurallySymmetric_IS;
3727:   A->ops->duplicate               = MatDuplicate_IS;
3728:   A->ops->missingdiagonal         = MatMissingDiagonal_IS;
3729:   A->ops->copy                    = MatCopy_IS;
3730:   A->ops->getlocalsubmatrix       = MatGetLocalSubMatrix_IS;
3731:   A->ops->createsubmatrix         = MatCreateSubMatrix_IS;
3732:   A->ops->axpy                    = MatAXPY_IS;
3733:   A->ops->diagonalset             = MatDiagonalSet_IS;
3734:   A->ops->shift                   = MatShift_IS;
3735:   A->ops->transpose               = MatTranspose_IS;
3736:   A->ops->getinfo                 = MatGetInfo_IS;
3737:   A->ops->diagonalscale           = MatDiagonalScale_IS;
3738:   A->ops->setfromoptions          = MatSetFromOptions_IS;
3739:   A->ops->setup                   = MatSetUp_IS;
3740:   A->ops->hasoperation            = MatHasOperation_IS;
3741:   A->ops->getdiagonalblock        = MatGetDiagonalBlock_IS;
3742:   A->ops->createsubmatrices       = MatCreateSubMatrices_IS;
3743:   A->ops->increaseoverlap         = MatIncreaseOverlap_IS;

3745:   /* special MATIS functions */
3746:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMatType_C", MatISSetLocalMatType_IS));
3747:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalMat_C", MatISGetLocalMat_IS));
3748:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISRestoreLocalMat_C", MatISRestoreLocalMat_IS));
3749:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMat_C", MatISSetLocalMat_IS));
3750:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetPreallocation_C", MatISSetPreallocation_IS));
3751:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetAllowRepeated_C", MatISSetAllowRepeated_IS));
3752:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISStoreL2L_C", MatISStoreL2L_IS));
3753:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISFixLocalEmpty_C", MatISFixLocalEmpty_IS));
3754:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalToGlobalMapping_C", MatISGetLocalToGlobalMapping_IS));
3755:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpiaij_C", MatConvert_IS_XAIJ));
3756:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpibaij_C", MatConvert_IS_XAIJ));
3757:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpisbaij_C", MatConvert_IS_XAIJ));
3758:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqaij_C", MatConvert_IS_XAIJ));
3759:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqbaij_C", MatConvert_IS_XAIJ));
3760:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqsbaij_C", MatConvert_IS_XAIJ));
3761:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_aij_C", MatConvert_IS_XAIJ));
3762:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOOLocal_C", MatSetPreallocationCOOLocal_IS));
3763:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOO_C", MatSetPreallocationCOO_IS));
3764:   PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATIS));
3765:   PetscFunctionReturn(PETSC_SUCCESS);
3766: }