Actual source code: matis.c

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

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

 10: #include <petsc/private/matisimpl.h>
 11: #include <petsc/private/sfimpl.h>
 12: #include <petsc/private/vecimpl.h>
 13: #include <petsc/private/hashseti.h>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

696:   /* access relevant information from MPIAIJ */
697:   PetscCall(MatGetOwnershipRange(A, &str, NULL));
698:   PetscCall(MatGetOwnershipRangeColumn(A, &stc, NULL));
699:   PetscCall(MatGetLocalSize(A, &dr, &dc));
700:   PetscCall(MatGetLocalSize(Ao, NULL, &oc));
701:   PetscCall(MatGetRowIJ(Ad, 0, PETSC_FALSE, PETSC_FALSE, &i, &di, &dj, &flg));
702:   PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure");
703:   PetscCall(MatGetRowIJ(Ao, 0, PETSC_FALSE, PETSC_FALSE, &i, &oi, &oj, &flg));
704:   PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure");
705:   nnz = di[dr] + oi[dr];
706:   /* store original pointers to be restored later */
707:   odi = di;
708:   odj = dj;
709:   ooi = oi;
710:   ooj = oj;

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

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

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

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

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

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

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

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

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

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

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

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

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

847:       /* Nested matrices should be of type MATIS */
848:       PetscCall(PetscObjectTypeCompare((PetscObject)nest[i][j], MATTRANSPOSEVIRTUAL, &istrans[ij]));
849:       if (istrans[ij]) {
850:         Mat T, lT;
851:         PetscCall(MatTransposeGetMat(nest[i][j], &T));
852:         PetscCall(PetscObjectTypeCompare((PetscObject)T, MATIS, &ismatis));
853:         PetscCheck(ismatis, comm, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") (transposed) is not of type MATIS", i, j);
854:         PetscCall(MatISGetLocalMat(T, &lT));
855:         PetscCall(MatCreateTranspose(lT, &snest[ij]));
856:       } else {
857:         PetscCall(PetscObjectTypeCompare((PetscObject)nest[i][j], MATIS, &ismatis));
858:         PetscCheck(ismatis, comm, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") is not of type MATIS", i, j);
859:         PetscCall(MatISGetLocalMat(nest[i][j], &snest[ij]));
860:       }

862:       /* Check compatibility of local sizes */
863:       PetscCall(MatGetSize(snest[ij], &l1, &l2));
864:       PetscCall(MatGetBlockSizes(snest[ij], &lb1, &lb2));
865:       if (!l1 || !l2) continue;
866:       PetscCheck(!lr[i] || l1 == lr[i], PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid local size %" PetscInt_FMT " != %" PetscInt_FMT, i, j, lr[i], l1);
867:       PetscCheck(!lc[j] || l2 == lc[j], PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid local size %" PetscInt_FMT " != %" PetscInt_FMT, i, j, lc[j], l2);
868:       lr[i] = l1;
869:       lc[j] = l2;

871:       /* check compatibility for local matrix reusage */
872:       if (rnest && !rnest[i][j] != !snest[ij]) lreuse = PETSC_FALSE;
873:     }
874:   }

876:   if (PetscDefined(USE_DEBUG)) {
877:     /* Check compatibility of l2g maps for rows */
878:     for (i = 0; i < nr; i++) {
879:       rl2g = NULL;
880:       for (j = 0; j < nc; j++) {
881:         PetscInt n1, n2;

883:         if (!nest[i][j]) continue;
884:         if (istrans[i * nc + j]) {
885:           Mat T;

887:           PetscCall(MatTransposeGetMat(nest[i][j], &T));
888:           PetscCall(MatISGetLocalToGlobalMapping(T, NULL, &cl2g));
889:         } else {
890:           PetscCall(MatISGetLocalToGlobalMapping(nest[i][j], &cl2g, NULL));
891:         }
892:         PetscCall(ISLocalToGlobalMappingGetSize(cl2g, &n1));
893:         if (!n1) continue;
894:         if (!rl2g) {
895:           rl2g = cl2g;
896:         } else {
897:           const PetscInt *idxs1, *idxs2;
898:           PetscBool       same;

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

917:         if (!nest[j][i]) continue;
918:         if (istrans[j * nc + i]) {
919:           Mat T;

921:           PetscCall(MatTransposeGetMat(nest[j][i], &T));
922:           PetscCall(MatISGetLocalToGlobalMapping(T, &cl2g, NULL));
923:         } else {
924:           PetscCall(MatISGetLocalToGlobalMapping(nest[j][i], NULL, &cl2g));
925:         }
926:         PetscCall(ISLocalToGlobalMappingGetSize(cl2g, &n1));
927:         if (!n1) continue;
928:         if (!rl2g) {
929:           rl2g = cl2g;
930:         } else {
931:           const PetscInt *idxs1, *idxs2;
932:           PetscBool       same;

934:           PetscCall(ISLocalToGlobalMappingGetSize(rl2g, &n2));
935:           PetscCheck(n1 == n2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid column l2gmap size %" PetscInt_FMT " != %" PetscInt_FMT, j, i, n1, n2);
936:           PetscCall(ISLocalToGlobalMappingGetIndices(cl2g, &idxs1));
937:           PetscCall(ISLocalToGlobalMappingGetIndices(rl2g, &idxs2));
938:           PetscCall(PetscArraycmp(idxs1, idxs2, n1, &same));
939:           PetscCall(ISLocalToGlobalMappingRestoreIndices(cl2g, &idxs1));
940:           PetscCall(ISLocalToGlobalMappingRestoreIndices(rl2g, &idxs2));
941:           PetscCheck(same, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid column l2gmap", j, i);
942:         }
943:       }
944:     }
945:   }

947:   B = NULL;
948:   if (reuse != MAT_REUSE_MATRIX) {
949:     PetscInt stl;

951:     /* Create l2g map for the rows of the new matrix and index sets for the local MATNEST */
952:     for (i = 0, stl = 0; i < nr; i++) stl += lr[i];
953:     PetscCall(PetscMalloc1(stl, &l2gidxs));
954:     for (i = 0, stl = 0; i < nr; i++) {
955:       Mat             usedmat;
956:       Mat_IS         *matis;
957:       const PetscInt *idxs;

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

962:       /* l2gmap */
963:       j       = 0;
964:       usedmat = nest[i][j];
965:       while (!usedmat && j < nc - 1) usedmat = nest[i][++j];
966:       PetscCheck(usedmat, comm, PETSC_ERR_SUP, "Cannot find valid row mat");

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

987:     /* Create l2g map for columns of the new matrix and index sets for the local MATNEST */
988:     for (i = 0, stl = 0; i < nc; i++) stl += lc[i];
989:     PetscCall(PetscMalloc1(stl, &l2gidxs));
990:     for (i = 0, stl = 0; i < nc; i++) {
991:       Mat             usedmat;
992:       Mat_IS         *matis;
993:       const PetscInt *idxs;

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

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

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

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

1088:   /* Create local matrix in MATNEST format */
1089:   convert = PETSC_FALSE;
1090:   PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)A)->prefix, "-matis_convert_local_nest", &convert, NULL));
1091:   if (convert) {
1092:     Mat              M;
1093:     MatISLocalFields lf;
1094:     PetscContainer   c;

1096:     PetscCall(MatISGetLocalMat(*newmat, &lA));
1097:     PetscCall(MatConvert(lA, MATAIJ, MAT_INITIAL_MATRIX, &M));
1098:     PetscCall(MatISSetLocalMat(*newmat, M));
1099:     PetscCall(MatDestroy(&M));

1101:     /* attach local fields to the matrix */
1102:     PetscCall(PetscNew(&lf));
1103:     PetscCall(PetscMalloc2(nr, &lf->rf, nc, &lf->cf));
1104:     for (i = 0; i < nr; i++) {
1105:       PetscInt n, st;

1107:       PetscCall(ISGetLocalSize(islrow[i], &n));
1108:       PetscCall(ISStrideGetInfo(islrow[i], &st, NULL));
1109:       PetscCall(ISCreateStride(comm, n, st, 1, &lf->rf[i]));
1110:     }
1111:     for (i = 0; i < nc; i++) {
1112:       PetscInt n, st;

1114:       PetscCall(ISGetLocalSize(islcol[i], &n));
1115:       PetscCall(ISStrideGetInfo(islcol[i], &st, NULL));
1116:       PetscCall(ISCreateStride(comm, n, st, 1, &lf->cf[i]));
1117:     }
1118:     lf->nr = nr;
1119:     lf->nc = nc;
1120:     PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)(*newmat)), &c));
1121:     PetscCall(PetscContainerSetPointer(c, lf));
1122:     PetscCall(PetscContainerSetUserDestroy(c, MatISContainerDestroyFields_Private));
1123:     PetscCall(PetscObjectCompose((PetscObject)(*newmat), "_convert_nest_lfields", (PetscObject)c));
1124:     PetscCall(PetscContainerDestroy(&c));
1125:   }

1127:   /* Free workspace */
1128:   for (i = 0; i < nr; i++) PetscCall(ISDestroy(&islrow[i]));
1129:   for (i = 0; i < nc; i++) PetscCall(ISDestroy(&islcol[i]));
1130:   PetscCall(PetscFree6(isrow, iscol, islrow, islcol, snest, istrans));
1131:   PetscCall(PetscFree2(lr, lc));
1132:   PetscFunctionReturn(PETSC_SUCCESS);
1133: }

1135: static PetscErrorCode MatDiagonalScale_IS(Mat A, Vec l, Vec r)
1136: {
1137:   Mat_IS            *matis = (Mat_IS *)A->data;
1138:   Vec                ll, rr;
1139:   const PetscScalar *Y, *X;
1140:   PetscScalar       *x, *y;

1142:   PetscFunctionBegin;
1143:   if (l) {
1144:     ll = matis->y;
1145:     PetscCall(VecGetArrayRead(l, &Y));
1146:     PetscCall(VecGetArray(ll, &y));
1147:     PetscCall(PetscSFBcastBegin(matis->sf, MPIU_SCALAR, Y, y, MPI_REPLACE));
1148:   } else {
1149:     ll = NULL;
1150:   }
1151:   if (r) {
1152:     rr = matis->x;
1153:     PetscCall(VecGetArrayRead(r, &X));
1154:     PetscCall(VecGetArray(rr, &x));
1155:     PetscCall(PetscSFBcastBegin(matis->csf, MPIU_SCALAR, X, x, MPI_REPLACE));
1156:   } else {
1157:     rr = NULL;
1158:   }
1159:   if (ll) {
1160:     PetscCall(PetscSFBcastEnd(matis->sf, MPIU_SCALAR, Y, y, MPI_REPLACE));
1161:     PetscCall(VecRestoreArrayRead(l, &Y));
1162:     PetscCall(VecRestoreArray(ll, &y));
1163:   }
1164:   if (rr) {
1165:     PetscCall(PetscSFBcastEnd(matis->csf, MPIU_SCALAR, X, x, MPI_REPLACE));
1166:     PetscCall(VecRestoreArrayRead(r, &X));
1167:     PetscCall(VecRestoreArray(rr, &x));
1168:   }
1169:   PetscCall(MatDiagonalScale(matis->A, ll, rr));
1170:   PetscFunctionReturn(PETSC_SUCCESS);
1171: }

1173: static PetscErrorCode MatGetInfo_IS(Mat A, MatInfoType flag, MatInfo *ginfo)
1174: {
1175:   Mat_IS        *matis = (Mat_IS *)A->data;
1176:   MatInfo        info;
1177:   PetscLogDouble isend[6], irecv[6];
1178:   PetscInt       bs;

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

1207:     ginfo->nz_used      = irecv[0];
1208:     ginfo->nz_allocated = irecv[1];
1209:     ginfo->nz_unneeded  = irecv[2];
1210:     ginfo->memory       = irecv[3];
1211:     ginfo->mallocs      = irecv[4];
1212:     ginfo->assemblies   = irecv[5];
1213:   } else if (flag == MAT_GLOBAL_SUM) {
1214:     PetscCall(MPIU_Allreduce(isend, irecv, 5, MPIU_PETSCLOGDOUBLE, MPI_SUM, PetscObjectComm((PetscObject)A)));

1216:     ginfo->nz_used      = irecv[0];
1217:     ginfo->nz_allocated = irecv[1];
1218:     ginfo->nz_unneeded  = irecv[2];
1219:     ginfo->memory       = irecv[3];
1220:     ginfo->mallocs      = irecv[4];
1221:     ginfo->assemblies   = A->num_ass;
1222:   }
1223:   ginfo->block_size        = bs;
1224:   ginfo->fill_ratio_given  = 0;
1225:   ginfo->fill_ratio_needed = 0;
1226:   ginfo->factor_mallocs    = 0;
1227:   PetscFunctionReturn(PETSC_SUCCESS);
1228: }

1230: static PetscErrorCode MatTranspose_IS(Mat A, MatReuse reuse, Mat *B)
1231: {
1232:   Mat C, lC, lA;

1234:   PetscFunctionBegin;
1235:   if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *B));
1236:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
1237:     ISLocalToGlobalMapping rl2g, cl2g;
1238:     PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
1239:     PetscCall(MatSetSizes(C, A->cmap->n, A->rmap->n, A->cmap->N, A->rmap->N));
1240:     PetscCall(MatSetBlockSizes(C, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs)));
1241:     PetscCall(MatSetType(C, MATIS));
1242:     PetscCall(MatGetLocalToGlobalMapping(A, &rl2g, &cl2g));
1243:     PetscCall(MatSetLocalToGlobalMapping(C, cl2g, rl2g));
1244:   } else C = *B;

1246:   /* perform local transposition */
1247:   PetscCall(MatISGetLocalMat(A, &lA));
1248:   PetscCall(MatTranspose(lA, MAT_INITIAL_MATRIX, &lC));
1249:   PetscCall(MatSetLocalToGlobalMapping(lC, lA->cmap->mapping, lA->rmap->mapping));
1250:   PetscCall(MatISSetLocalMat(C, lC));
1251:   PetscCall(MatDestroy(&lC));

1253:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
1254:     *B = C;
1255:   } else {
1256:     PetscCall(MatHeaderMerge(A, &C));
1257:   }
1258:   PetscCall(MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY));
1259:   PetscCall(MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY));
1260:   PetscFunctionReturn(PETSC_SUCCESS);
1261: }

1263: static PetscErrorCode MatDiagonalSet_IS(Mat A, Vec D, InsertMode insmode)
1264: {
1265:   Mat_IS *is = (Mat_IS *)A->data;

1267:   PetscFunctionBegin;
1268:   if (D) { /* MatShift_IS pass D = NULL */
1269:     PetscCall(VecScatterBegin(is->rctx, D, is->y, INSERT_VALUES, SCATTER_FORWARD));
1270:     PetscCall(VecScatterEnd(is->rctx, D, is->y, INSERT_VALUES, SCATTER_FORWARD));
1271:   }
1272:   PetscCall(VecPointwiseDivide(is->y, is->y, is->counter));
1273:   PetscCall(MatDiagonalSet(is->A, is->y, insmode));
1274:   PetscFunctionReturn(PETSC_SUCCESS);
1275: }

1277: static PetscErrorCode MatShift_IS(Mat A, PetscScalar a)
1278: {
1279:   Mat_IS *is = (Mat_IS *)A->data;

1281:   PetscFunctionBegin;
1282:   PetscCall(VecSet(is->y, a));
1283:   PetscCall(MatDiagonalSet_IS(A, NULL, ADD_VALUES));
1284:   PetscFunctionReturn(PETSC_SUCCESS);
1285: }

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

1291:   PetscFunctionBegin;
1292:   PetscCheck(m <= MATIS_MAX_ENTRIES_INSERTION && n <= MATIS_MAX_ENTRIES_INSERTION, PETSC_COMM_SELF, PETSC_ERR_SUP, "Number of row/column indices must be <= %d: they are %" PetscInt_FMT " %" PetscInt_FMT, MATIS_MAX_ENTRIES_INSERTION, m, n);
1293:   PetscCall(ISLocalToGlobalMappingApply(A->rmap->mapping, m, rows, rows_l));
1294:   PetscCall(ISLocalToGlobalMappingApply(A->cmap->mapping, n, cols, cols_l));
1295:   PetscCall(MatSetValuesLocal_IS(A, m, rows_l, n, cols_l, values, addv));
1296:   PetscFunctionReturn(PETSC_SUCCESS);
1297: }

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

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

1311: static PetscErrorCode MatCreateSubMatrix_IS(Mat mat, IS irow, IS icol, MatReuse scall, Mat *newmat)
1312: {
1313:   Mat             locmat, newlocmat;
1314:   Mat_IS         *newmatis;
1315:   const PetscInt *idxs;
1316:   PetscInt        i, m, n;

1318:   PetscFunctionBegin;
1319:   if (scall == MAT_REUSE_MATRIX) {
1320:     PetscBool ismatis;

1322:     PetscCall(PetscObjectTypeCompare((PetscObject)*newmat, MATIS, &ismatis));
1323:     PetscCheck(ismatis, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_ARG_WRONG, "Cannot reuse matrix! Not of MATIS type");
1324:     newmatis = (Mat_IS *)(*newmat)->data;
1325:     PetscCheck(newmatis->getsub_ris, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_ARG_WRONG, "Cannot reuse matrix! Misses local row IS");
1326:     PetscCheck(newmatis->getsub_cis, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_ARG_WRONG, "Cannot reuse matrix! Misses local col IS");
1327:   }
1328:   /* irow and icol may not have duplicate entries */
1329:   if (PetscDefined(USE_DEBUG)) {
1330:     Vec                rtest, ltest;
1331:     const PetscScalar *array;

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

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

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

1458: static PetscErrorCode MatCopy_IS(Mat A, Mat B, MatStructure str)
1459: {
1460:   Mat_IS   *a = (Mat_IS *)A->data, *b;
1461:   PetscBool ismatis;

1463:   PetscFunctionBegin;
1464:   PetscCall(PetscObjectTypeCompare((PetscObject)B, MATIS, &ismatis));
1465:   PetscCheck(ismatis, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Need to be implemented");
1466:   b = (Mat_IS *)B->data;
1467:   PetscCall(MatCopy(a->A, b->A, str));
1468:   PetscCall(PetscObjectStateIncrease((PetscObject)B));
1469:   PetscFunctionReturn(PETSC_SUCCESS);
1470: }

1472: static PetscErrorCode MatMissingDiagonal_IS(Mat A, PetscBool *missing, PetscInt *d)
1473: {
1474:   Vec                v;
1475:   const PetscScalar *array;
1476:   PetscInt           i, n;

1478:   PetscFunctionBegin;
1479:   *missing = PETSC_FALSE;
1480:   PetscCall(MatCreateVecs(A, NULL, &v));
1481:   PetscCall(MatGetDiagonal(A, v));
1482:   PetscCall(VecGetLocalSize(v, &n));
1483:   PetscCall(VecGetArrayRead(v, &array));
1484:   for (i = 0; i < n; i++)
1485:     if (array[i] == 0.) break;
1486:   PetscCall(VecRestoreArrayRead(v, &array));
1487:   PetscCall(VecDestroy(&v));
1488:   if (i != n) *missing = PETSC_TRUE;
1489:   if (d) {
1490:     *d = -1;
1491:     if (*missing) {
1492:       PetscInt rstart;
1493:       PetscCall(MatGetOwnershipRange(A, &rstart, NULL));
1494:       *d = i + rstart;
1495:     }
1496:   }
1497:   PetscFunctionReturn(PETSC_SUCCESS);
1498: }

1500: static PetscErrorCode MatISSetUpSF_IS(Mat B)
1501: {
1502:   Mat_IS         *matis = (Mat_IS *)(B->data);
1503:   const PetscInt *gidxs;
1504:   PetscInt        nleaves;

1506:   PetscFunctionBegin;
1507:   if (matis->sf) PetscFunctionReturn(PETSC_SUCCESS);
1508:   PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)B), &matis->sf));
1509:   PetscCall(ISLocalToGlobalMappingGetIndices(matis->rmapping, &gidxs));
1510:   PetscCall(ISLocalToGlobalMappingGetSize(matis->rmapping, &nleaves));
1511:   PetscCall(PetscSFSetGraphLayout(matis->sf, B->rmap, nleaves, NULL, PETSC_OWN_POINTER, gidxs));
1512:   PetscCall(ISLocalToGlobalMappingRestoreIndices(matis->rmapping, &gidxs));
1513:   PetscCall(PetscMalloc2(matis->sf->nroots, &matis->sf_rootdata, matis->sf->nleaves, &matis->sf_leafdata));
1514:   if (matis->rmapping != matis->cmapping) { /* setup SF for columns */
1515:     PetscCall(ISLocalToGlobalMappingGetSize(matis->cmapping, &nleaves));
1516:     PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)B), &matis->csf));
1517:     PetscCall(ISLocalToGlobalMappingGetIndices(matis->cmapping, &gidxs));
1518:     PetscCall(PetscSFSetGraphLayout(matis->csf, B->cmap, nleaves, NULL, PETSC_OWN_POINTER, gidxs));
1519:     PetscCall(ISLocalToGlobalMappingRestoreIndices(matis->cmapping, &gidxs));
1520:     PetscCall(PetscMalloc2(matis->csf->nroots, &matis->csf_rootdata, matis->csf->nleaves, &matis->csf_leafdata));
1521:   } else {
1522:     matis->csf          = matis->sf;
1523:     matis->csf_leafdata = matis->sf_leafdata;
1524:     matis->csf_rootdata = matis->sf_rootdata;
1525:   }
1526:   PetscFunctionReturn(PETSC_SUCCESS);
1527: }

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

1532:   Collective

1534:   Input Parameters:
1535: + A     - the matrix
1536: - store - the boolean flag

1538:   Level: advanced

1540: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateIS()`, `MatISSetPreallocation()`, `MatPtAP()`
1541: @*/
1542: PetscErrorCode MatISStoreL2L(Mat A, PetscBool store)
1543: {
1544:   PetscFunctionBegin;
1548:   PetscTryMethod(A, "MatISStoreL2L_C", (Mat, PetscBool), (A, store));
1549:   PetscFunctionReturn(PETSC_SUCCESS);
1550: }

1552: static PetscErrorCode MatISStoreL2L_IS(Mat A, PetscBool store)
1553: {
1554:   Mat_IS *matis = (Mat_IS *)(A->data);

1556:   PetscFunctionBegin;
1557:   matis->storel2l = store;
1558:   if (!store) PetscCall(PetscObjectCompose((PetscObject)(A), "_MatIS_PtAP_l2l", NULL));
1559:   PetscFunctionReturn(PETSC_SUCCESS);
1560: }

1562: /*@
1563:   MatISFixLocalEmpty - Compress out zero local rows from the local matrices

1565:   Collective

1567:   Input Parameters:
1568: + A   - the matrix
1569: - fix - the boolean flag

1571:   Level: advanced

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

1576: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatCreate()`, `MatCreateIS()`, `MatISSetPreallocation()`, `MatAssemblyEnd()`, `MAT_FINAL_ASSEMBLY`
1577: @*/
1578: PetscErrorCode MatISFixLocalEmpty(Mat A, PetscBool fix)
1579: {
1580:   PetscFunctionBegin;
1584:   PetscTryMethod(A, "MatISFixLocalEmpty_C", (Mat, PetscBool), (A, fix));
1585:   PetscFunctionReturn(PETSC_SUCCESS);
1586: }

1588: static PetscErrorCode MatISFixLocalEmpty_IS(Mat A, PetscBool fix)
1589: {
1590:   Mat_IS *matis = (Mat_IS *)(A->data);

1592:   PetscFunctionBegin;
1593:   matis->locempty = fix;
1594:   PetscFunctionReturn(PETSC_SUCCESS);
1595: }

1597: /*@
1598:   MatISSetPreallocation - Preallocates memory for a `MATIS` parallel matrix.

1600:   Collective

1602:   Input Parameters:
1603: + B     - the matrix
1604: . d_nz  - number of nonzeros per row in DIAGONAL portion of local submatrix
1605:            (same value is used for all local rows)
1606: . d_nnz - array containing the number of nonzeros in the various rows of the
1607:            DIAGONAL portion of the local submatrix (possibly different for each row)
1608:            or `NULL`, if `d_nz` is used to specify the nonzero structure.
1609:            The size of this array is equal to the number of local rows, i.e `m`.
1610:            For matrices that will be factored, you must leave room for (and set)
1611:            the diagonal entry even if it is zero.
1612: . o_nz  - number of nonzeros per row in the OFF-DIAGONAL portion of local
1613:            submatrix (same value is used for all local rows).
1614: - o_nnz - array containing the number of nonzeros in the various rows of the
1615:            OFF-DIAGONAL portion of the local submatrix (possibly different for
1616:            each row) or `NULL`, if `o_nz` is used to specify the nonzero
1617:            structure. The size of this array is equal to the number
1618:            of local rows, i.e `m`.

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

1622:   Level: intermediate

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

1629: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateIS()`, `MatMPIAIJSetPreallocation()`, `MatISGetLocalMat()`, `MATIS`
1630: @*/
1631: PetscErrorCode MatISSetPreallocation(Mat B, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[])
1632: {
1633:   PetscFunctionBegin;
1636:   PetscTryMethod(B, "MatISSetPreallocation_C", (Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[]), (B, d_nz, d_nnz, o_nz, o_nnz));
1637:   PetscFunctionReturn(PETSC_SUCCESS);
1638: }

1640: /* this is used by DMDA */
1641: PETSC_EXTERN PetscErrorCode MatISSetPreallocation_IS(Mat B, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[])
1642: {
1643:   Mat_IS  *matis = (Mat_IS *)(B->data);
1644:   PetscInt bs, i, nlocalcols;

1646:   PetscFunctionBegin;
1647:   PetscCall(MatSetUp(B));
1648:   if (!d_nnz)
1649:     for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] = d_nz;
1650:   else
1651:     for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] = d_nnz[i];

1653:   if (!o_nnz)
1654:     for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] += o_nz;
1655:   else
1656:     for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] += o_nnz[i];

1658:   PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
1659:   PetscCall(MatGetSize(matis->A, NULL, &nlocalcols));
1660:   PetscCall(MatGetBlockSize(matis->A, &bs));
1661:   PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));

1663:   for (i = 0; i < matis->sf->nleaves; i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i], nlocalcols);
1664:   PetscCall(MatSeqAIJSetPreallocation(matis->A, 0, matis->sf_leafdata));
1665: #if defined(PETSC_HAVE_HYPRE)
1666:   PetscCall(MatHYPRESetPreallocation(matis->A, 0, matis->sf_leafdata, 0, NULL));
1667: #endif

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

1672:     matis->sf_leafdata[i] = matis->sf_leafdata[i * bs] / bs;
1673:     for (b = 1; b < bs; b++) matis->sf_leafdata[i] = PetscMax(matis->sf_leafdata[i], matis->sf_leafdata[i * bs + b] / bs);
1674:   }
1675:   PetscCall(MatSeqBAIJSetPreallocation(matis->A, bs, 0, matis->sf_leafdata));

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

1681:   /* for other matrix types */
1682:   PetscCall(MatSetUp(matis->A));
1683:   PetscFunctionReturn(PETSC_SUCCESS);
1684: }

1686: PETSC_EXTERN PetscErrorCode MatISSetMPIXAIJPreallocation_Private(Mat A, Mat B, PetscBool maxreduce)
1687: {
1688:   Mat_IS         *matis = (Mat_IS *)(A->data);
1689:   PetscInt       *my_dnz, *my_onz, *dnz, *onz, *mat_ranges, *row_ownership;
1690:   const PetscInt *global_indices_r, *global_indices_c;
1691:   PetscInt        i, j, bs, rows, cols;
1692:   PetscInt        lrows, lcols;
1693:   PetscInt        local_rows, local_cols;
1694:   PetscMPIInt     size;
1695:   PetscBool       isdense, issbaij;

1697:   PetscFunctionBegin;
1698:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1699:   PetscCall(MatGetSize(A, &rows, &cols));
1700:   PetscCall(MatGetBlockSize(A, &bs));
1701:   PetscCall(MatGetSize(matis->A, &local_rows, &local_cols));
1702:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQDENSE, &isdense));
1703:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQSBAIJ, &issbaij));
1704:   PetscCall(ISLocalToGlobalMappingGetIndices(matis->rmapping, &global_indices_r));
1705:   if (matis->rmapping != matis->cmapping) {
1706:     PetscCall(ISLocalToGlobalMappingGetIndices(matis->cmapping, &global_indices_c));
1707:   } else global_indices_c = global_indices_r;

1709:   if (issbaij) PetscCall(MatGetRowUpperTriangular(matis->A));
1710:   /*
1711:      An SF reduce is needed to sum up properly on shared rows.
1712:      Note that generally preallocation is not exact, since it overestimates nonzeros
1713:   */
1714:   PetscCall(MatGetLocalSize(A, &lrows, &lcols));
1715:   MatPreallocateBegin(PetscObjectComm((PetscObject)A), lrows, lcols, dnz, onz);
1716:   /* All processes need to compute entire row ownership */
1717:   PetscCall(PetscMalloc1(rows, &row_ownership));
1718:   PetscCall(MatGetOwnershipRanges(A, (const PetscInt **)&mat_ranges));
1719:   for (i = 0; i < size; i++) {
1720:     for (j = mat_ranges[i]; j < mat_ranges[i + 1]; j++) row_ownership[j] = i;
1721:   }
1722:   PetscCall(MatGetOwnershipRangesColumn(A, (const PetscInt **)&mat_ranges));

1724:   /*
1725:      my_dnz and my_onz contains exact contribution to preallocation from each local mat
1726:      then, they will be summed up properly. This way, preallocation is always sufficient
1727:   */
1728:   PetscCall(PetscCalloc2(local_rows, &my_dnz, local_rows, &my_onz));
1729:   /* preallocation as a MATAIJ */
1730:   if (isdense) { /* special case for dense local matrices */
1731:     for (i = 0; i < local_rows; i++) {
1732:       PetscInt owner = row_ownership[global_indices_r[i]];
1733:       for (j = 0; j < local_cols; j++) {
1734:         PetscInt index_col = global_indices_c[j];
1735:         if (index_col > mat_ranges[owner] - 1 && index_col < mat_ranges[owner + 1]) { /* diag block */
1736:           my_dnz[i] += 1;
1737:         } else { /* offdiag block */
1738:           my_onz[i] += 1;
1739:         }
1740:       }
1741:     }
1742:   } else if (matis->A->ops->getrowij) {
1743:     const PetscInt *ii, *jj, *jptr;
1744:     PetscBool       done;
1745:     PetscCall(MatGetRowIJ(matis->A, 0, PETSC_FALSE, PETSC_FALSE, &local_rows, &ii, &jj, &done));
1746:     PetscCheck(done, PetscObjectComm((PetscObject)(matis->A)), PETSC_ERR_PLIB, "Error in MatGetRowIJ");
1747:     jptr = jj;
1748:     for (i = 0; i < local_rows; i++) {
1749:       PetscInt index_row = global_indices_r[i];
1750:       for (j = 0; j < ii[i + 1] - ii[i]; j++, jptr++) {
1751:         PetscInt owner     = row_ownership[index_row];
1752:         PetscInt index_col = global_indices_c[*jptr];
1753:         if (index_col > mat_ranges[owner] - 1 && index_col < mat_ranges[owner + 1]) { /* diag block */
1754:           my_dnz[i] += 1;
1755:         } else { /* offdiag block */
1756:           my_onz[i] += 1;
1757:         }
1758:         /* same as before, interchanging rows and cols */
1759:         if (issbaij && index_col != index_row) {
1760:           owner = row_ownership[index_col];
1761:           if (index_row > mat_ranges[owner] - 1 && index_row < mat_ranges[owner + 1]) {
1762:             my_dnz[*jptr] += 1;
1763:           } else {
1764:             my_onz[*jptr] += 1;
1765:           }
1766:         }
1767:       }
1768:     }
1769:     PetscCall(MatRestoreRowIJ(matis->A, 0, PETSC_FALSE, PETSC_FALSE, &local_rows, &ii, &jj, &done));
1770:     PetscCheck(done, PetscObjectComm((PetscObject)(matis->A)), PETSC_ERR_PLIB, "Error in MatRestoreRowIJ");
1771:   } else { /* loop over rows and use MatGetRow */
1772:     for (i = 0; i < local_rows; i++) {
1773:       const PetscInt *cols;
1774:       PetscInt        ncols, index_row = global_indices_r[i];
1775:       PetscCall(MatGetRow(matis->A, i, &ncols, &cols, NULL));
1776:       for (j = 0; j < ncols; j++) {
1777:         PetscInt owner     = row_ownership[index_row];
1778:         PetscInt index_col = global_indices_c[cols[j]];
1779:         if (index_col > mat_ranges[owner] - 1 && index_col < mat_ranges[owner + 1]) { /* diag block */
1780:           my_dnz[i] += 1;
1781:         } else { /* offdiag block */
1782:           my_onz[i] += 1;
1783:         }
1784:         /* same as before, interchanging rows and cols */
1785:         if (issbaij && index_col != index_row) {
1786:           owner = row_ownership[index_col];
1787:           if (index_row > mat_ranges[owner] - 1 && index_row < mat_ranges[owner + 1]) {
1788:             my_dnz[cols[j]] += 1;
1789:           } else {
1790:             my_onz[cols[j]] += 1;
1791:           }
1792:         }
1793:       }
1794:       PetscCall(MatRestoreRow(matis->A, i, &ncols, &cols, NULL));
1795:     }
1796:   }
1797:   if (global_indices_c != global_indices_r) PetscCall(ISLocalToGlobalMappingRestoreIndices(matis->cmapping, &global_indices_c));
1798:   PetscCall(ISLocalToGlobalMappingRestoreIndices(matis->rmapping, &global_indices_r));
1799:   PetscCall(PetscFree(row_ownership));

1801:   /* Reduce my_dnz and my_onz */
1802:   if (maxreduce) {
1803:     PetscCall(PetscSFReduceBegin(matis->sf, MPIU_INT, my_dnz, dnz, MPI_MAX));
1804:     PetscCall(PetscSFReduceBegin(matis->sf, MPIU_INT, my_onz, onz, MPI_MAX));
1805:     PetscCall(PetscSFReduceEnd(matis->sf, MPIU_INT, my_dnz, dnz, MPI_MAX));
1806:     PetscCall(PetscSFReduceEnd(matis->sf, MPIU_INT, my_onz, onz, MPI_MAX));
1807:   } else {
1808:     PetscCall(PetscSFReduceBegin(matis->sf, MPIU_INT, my_dnz, dnz, MPI_SUM));
1809:     PetscCall(PetscSFReduceBegin(matis->sf, MPIU_INT, my_onz, onz, MPI_SUM));
1810:     PetscCall(PetscSFReduceEnd(matis->sf, MPIU_INT, my_dnz, dnz, MPI_SUM));
1811:     PetscCall(PetscSFReduceEnd(matis->sf, MPIU_INT, my_onz, onz, MPI_SUM));
1812:   }
1813:   PetscCall(PetscFree2(my_dnz, my_onz));

1815:   /* Resize preallocation if overestimated */
1816:   for (i = 0; i < lrows; i++) {
1817:     dnz[i] = PetscMin(dnz[i], lcols);
1818:     onz[i] = PetscMin(onz[i], cols - lcols);
1819:   }

1821:   /* Set preallocation */
1822:   PetscCall(MatSetBlockSizesFromMats(B, A, A));
1823:   PetscCall(MatSeqAIJSetPreallocation(B, 0, dnz));
1824:   PetscCall(MatMPIAIJSetPreallocation(B, 0, dnz, 0, onz));
1825:   for (i = 0; i < lrows; i += bs) {
1826:     PetscInt b, d = dnz[i], o = onz[i];

1828:     for (b = 1; b < bs; b++) {
1829:       d = PetscMax(d, dnz[i + b]);
1830:       o = PetscMax(o, onz[i + b]);
1831:     }
1832:     dnz[i / bs] = PetscMin(d / bs + d % bs, lcols / bs);
1833:     onz[i / bs] = PetscMin(o / bs + o % bs, (cols - lcols) / bs);
1834:   }
1835:   PetscCall(MatSeqBAIJSetPreallocation(B, bs, 0, dnz));
1836:   PetscCall(MatMPIBAIJSetPreallocation(B, bs, 0, dnz, 0, onz));
1837:   PetscCall(MatMPISBAIJSetPreallocation(B, bs, 0, dnz, 0, onz));
1838:   MatPreallocateEnd(dnz, onz);
1839:   if (issbaij) PetscCall(MatRestoreRowUpperTriangular(matis->A));
1840:   PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
1841:   PetscFunctionReturn(PETSC_SUCCESS);
1842: }

1844: PETSC_INTERN PetscErrorCode MatConvert_IS_XAIJ(Mat mat, MatType mtype, MatReuse reuse, Mat *M)
1845: {
1846:   Mat_IS            *matis = (Mat_IS *)(mat->data);
1847:   Mat                local_mat, MT;
1848:   PetscInt           rbs, cbs, rows, cols, lrows, lcols;
1849:   PetscInt           local_rows, local_cols;
1850:   PetscBool          isseqdense, isseqsbaij, isseqaij, isseqbaij;
1851:   PetscMPIInt        size;
1852:   const PetscScalar *array;

1854:   PetscFunctionBegin;
1855:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)mat), &size));
1856:   if (size == 1 && mat->rmap->N == matis->A->rmap->N && mat->cmap->N == matis->A->cmap->N) {
1857:     Mat      B;
1858:     IS       irows = NULL, icols = NULL;
1859:     PetscInt rbs, cbs;

1861:     PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->rmapping, &rbs));
1862:     PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->cmapping, &cbs));
1863:     if (reuse != MAT_REUSE_MATRIX) { /* check if l2g maps are one-to-one */
1864:       IS              rows, cols;
1865:       const PetscInt *ridxs, *cidxs;
1866:       PetscInt        i, nw, *work;

1868:       PetscCall(ISLocalToGlobalMappingGetBlockIndices(matis->rmapping, &ridxs));
1869:       PetscCall(ISLocalToGlobalMappingGetSize(matis->rmapping, &nw));
1870:       nw = nw / rbs;
1871:       PetscCall(PetscCalloc1(nw, &work));
1872:       for (i = 0; i < nw; i++) work[ridxs[i]] += 1;
1873:       for (i = 0; i < nw; i++)
1874:         if (!work[i] || work[i] > 1) break;
1875:       if (i == nw) {
1876:         PetscCall(ISCreateBlock(PETSC_COMM_SELF, rbs, nw, ridxs, PETSC_USE_POINTER, &rows));
1877:         PetscCall(ISSetPermutation(rows));
1878:         PetscCall(ISInvertPermutation(rows, PETSC_DECIDE, &irows));
1879:         PetscCall(ISDestroy(&rows));
1880:       }
1881:       PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(matis->rmapping, &ridxs));
1882:       PetscCall(PetscFree(work));
1883:       if (irows && matis->rmapping != matis->cmapping) {
1884:         PetscCall(ISLocalToGlobalMappingGetBlockIndices(matis->cmapping, &cidxs));
1885:         PetscCall(ISLocalToGlobalMappingGetSize(matis->cmapping, &nw));
1886:         nw = nw / cbs;
1887:         PetscCall(PetscCalloc1(nw, &work));
1888:         for (i = 0; i < nw; i++) work[cidxs[i]] += 1;
1889:         for (i = 0; i < nw; i++)
1890:           if (!work[i] || work[i] > 1) break;
1891:         if (i == nw) {
1892:           PetscCall(ISCreateBlock(PETSC_COMM_SELF, cbs, nw, cidxs, PETSC_USE_POINTER, &cols));
1893:           PetscCall(ISSetPermutation(cols));
1894:           PetscCall(ISInvertPermutation(cols, PETSC_DECIDE, &icols));
1895:           PetscCall(ISDestroy(&cols));
1896:         }
1897:         PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(matis->cmapping, &cidxs));
1898:         PetscCall(PetscFree(work));
1899:       } else if (irows) {
1900:         PetscCall(PetscObjectReference((PetscObject)irows));
1901:         icols = irows;
1902:       }
1903:     } else {
1904:       PetscCall(PetscObjectQuery((PetscObject)(*M), "_MatIS_IS_XAIJ_irows", (PetscObject *)&irows));
1905:       PetscCall(PetscObjectQuery((PetscObject)(*M), "_MatIS_IS_XAIJ_icols", (PetscObject *)&icols));
1906:       if (irows) PetscCall(PetscObjectReference((PetscObject)irows));
1907:       if (icols) PetscCall(PetscObjectReference((PetscObject)icols));
1908:     }
1909:     if (!irows || !icols) {
1910:       PetscCall(ISDestroy(&icols));
1911:       PetscCall(ISDestroy(&irows));
1912:       goto general_assembly;
1913:     }
1914:     PetscCall(MatConvert(matis->A, mtype, MAT_INITIAL_MATRIX, &B));
1915:     if (reuse != MAT_INPLACE_MATRIX) {
1916:       PetscCall(MatCreateSubMatrix(B, irows, icols, reuse, M));
1917:       PetscCall(PetscObjectCompose((PetscObject)(*M), "_MatIS_IS_XAIJ_irows", (PetscObject)irows));
1918:       PetscCall(PetscObjectCompose((PetscObject)(*M), "_MatIS_IS_XAIJ_icols", (PetscObject)icols));
1919:     } else {
1920:       Mat C;

1922:       PetscCall(MatCreateSubMatrix(B, irows, icols, MAT_INITIAL_MATRIX, &C));
1923:       PetscCall(MatHeaderReplace(mat, &C));
1924:     }
1925:     PetscCall(MatDestroy(&B));
1926:     PetscCall(ISDestroy(&icols));
1927:     PetscCall(ISDestroy(&irows));
1928:     PetscFunctionReturn(PETSC_SUCCESS);
1929:   }
1930: general_assembly:
1931:   PetscCall(MatGetSize(mat, &rows, &cols));
1932:   PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->rmapping, &rbs));
1933:   PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->cmapping, &cbs));
1934:   PetscCall(MatGetLocalSize(mat, &lrows, &lcols));
1935:   PetscCall(MatGetSize(matis->A, &local_rows, &local_cols));
1936:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQDENSE, &isseqdense));
1937:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQAIJ, &isseqaij));
1938:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQBAIJ, &isseqbaij));
1939:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQSBAIJ, &isseqsbaij));
1940:   PetscCheck(isseqdense || isseqaij || isseqbaij || isseqsbaij, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for matrix type %s", ((PetscObject)(matis->A))->type_name);
1941:   if (PetscDefined(USE_DEBUG)) {
1942:     PetscBool lb[4], bb[4];

1944:     lb[0] = isseqdense;
1945:     lb[1] = isseqaij;
1946:     lb[2] = isseqbaij;
1947:     lb[3] = isseqsbaij;
1948:     PetscCall(MPIU_Allreduce(lb, bb, 4, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)mat)));
1949:     PetscCheck(bb[0] || bb[1] || bb[2] || bb[3], PETSC_COMM_SELF, PETSC_ERR_SUP, "Local matrices must have the same type");
1950:   }

1952:   if (reuse != MAT_REUSE_MATRIX) {
1953:     PetscCall(MatCreate(PetscObjectComm((PetscObject)mat), &MT));
1954:     PetscCall(MatSetSizes(MT, lrows, lcols, rows, cols));
1955:     PetscCall(MatSetType(MT, mtype));
1956:     PetscCall(MatSetBlockSizes(MT, rbs, cbs));
1957:     PetscCall(MatISSetMPIXAIJPreallocation_Private(mat, MT, PETSC_FALSE));
1958:   } else {
1959:     PetscInt mrbs, mcbs, mrows, mcols, mlrows, mlcols;

1961:     /* some checks */
1962:     MT = *M;
1963:     PetscCall(MatGetBlockSizes(MT, &mrbs, &mcbs));
1964:     PetscCall(MatGetSize(MT, &mrows, &mcols));
1965:     PetscCall(MatGetLocalSize(MT, &mlrows, &mlcols));
1966:     PetscCheck(mrows == rows, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of rows (%" PetscInt_FMT " != %" PetscInt_FMT ")", rows, mrows);
1967:     PetscCheck(mcols == cols, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of cols (%" PetscInt_FMT " != %" PetscInt_FMT ")", cols, mcols);
1968:     PetscCheck(mlrows == lrows, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of local rows (%" PetscInt_FMT " != %" PetscInt_FMT ")", lrows, mlrows);
1969:     PetscCheck(mlcols == lcols, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of local cols (%" PetscInt_FMT " != %" PetscInt_FMT ")", lcols, mlcols);
1970:     PetscCheck(mrbs == rbs, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong row block size (%" PetscInt_FMT " != %" PetscInt_FMT ")", rbs, mrbs);
1971:     PetscCheck(mcbs == cbs, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong col block size (%" PetscInt_FMT " != %" PetscInt_FMT ")", cbs, mcbs);
1972:     PetscCall(MatZeroEntries(MT));
1973:   }

1975:   if (isseqsbaij || isseqbaij) {
1976:     PetscCall(MatConvert(matis->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &local_mat));
1977:     isseqaij = PETSC_TRUE;
1978:   } else {
1979:     PetscCall(PetscObjectReference((PetscObject)matis->A));
1980:     local_mat = matis->A;
1981:   }

1983:   /* Set values */
1984:   PetscCall(MatSetLocalToGlobalMapping(MT, matis->rmapping, matis->cmapping));
1985:   if (isseqdense) { /* special case for dense local matrices */
1986:     PetscInt i, *dummy;

1988:     PetscCall(PetscMalloc1(PetscMax(local_rows, local_cols), &dummy));
1989:     for (i = 0; i < PetscMax(local_rows, local_cols); i++) dummy[i] = i;
1990:     PetscCall(MatSetOption(MT, MAT_ROW_ORIENTED, PETSC_FALSE));
1991:     PetscCall(MatDenseGetArrayRead(local_mat, &array));
1992:     PetscCall(MatSetValuesLocal(MT, local_rows, dummy, local_cols, dummy, array, ADD_VALUES));
1993:     PetscCall(MatDenseRestoreArrayRead(local_mat, &array));
1994:     PetscCall(PetscFree(dummy));
1995:   } else if (isseqaij) {
1996:     const PetscInt *blocks;
1997:     PetscInt        i, nvtxs, *xadj, *adjncy, nb;
1998:     PetscBool       done;
1999:     PetscScalar    *sarray;

2001:     PetscCall(MatGetRowIJ(local_mat, 0, PETSC_FALSE, PETSC_FALSE, &nvtxs, (const PetscInt **)&xadj, (const PetscInt **)&adjncy, &done));
2002:     PetscCheck(done, PetscObjectComm((PetscObject)local_mat), PETSC_ERR_PLIB, "Error in MatGetRowIJ");
2003:     PetscCall(MatSeqAIJGetArray(local_mat, &sarray));
2004:     PetscCall(MatGetVariableBlockSizes(local_mat, &nb, &blocks));
2005:     if (nb) { /* speed up assembly for special blocked matrices (used by BDDC) */
2006:       PetscInt sum;

2008:       for (i = 0, sum = 0; i < nb; i++) sum += blocks[i];
2009:       if (sum == nvtxs) {
2010:         PetscInt r;

2012:         for (i = 0, r = 0; i < nb; i++) {
2013:           PetscAssert(blocks[i] == xadj[r + 1] - xadj[r], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid block sizes prescribed for block %" PetscInt_FMT ": expected %" PetscInt_FMT ", got %" PetscInt_FMT, i, blocks[i], xadj[r + 1] - xadj[r]);
2014:           PetscCall(MatSetValuesLocal(MT, blocks[i], adjncy + xadj[r], blocks[i], adjncy + xadj[r], sarray + xadj[r], ADD_VALUES));
2015:           r += blocks[i];
2016:         }
2017:       } else {
2018:         for (i = 0; i < nvtxs; i++) PetscCall(MatSetValuesLocal(MT, 1, &i, xadj[i + 1] - xadj[i], adjncy + xadj[i], sarray + xadj[i], ADD_VALUES));
2019:       }
2020:     } else {
2021:       for (i = 0; i < nvtxs; i++) PetscCall(MatSetValuesLocal(MT, 1, &i, xadj[i + 1] - xadj[i], adjncy + xadj[i], sarray + xadj[i], ADD_VALUES));
2022:     }
2023:     PetscCall(MatRestoreRowIJ(local_mat, 0, PETSC_FALSE, PETSC_FALSE, &nvtxs, (const PetscInt **)&xadj, (const PetscInt **)&adjncy, &done));
2024:     PetscCheck(done, PetscObjectComm((PetscObject)local_mat), PETSC_ERR_PLIB, "Error in MatRestoreRowIJ");
2025:     PetscCall(MatSeqAIJRestoreArray(local_mat, &sarray));
2026:   } else { /* very basic values insertion for all other matrix types */
2027:     PetscInt i;

2029:     for (i = 0; i < local_rows; i++) {
2030:       PetscInt        j;
2031:       const PetscInt *local_indices_cols;

2033:       PetscCall(MatGetRow(local_mat, i, &j, &local_indices_cols, &array));
2034:       PetscCall(MatSetValuesLocal(MT, 1, &i, j, local_indices_cols, array, ADD_VALUES));
2035:       PetscCall(MatRestoreRow(local_mat, i, &j, &local_indices_cols, &array));
2036:     }
2037:   }
2038:   PetscCall(MatAssemblyBegin(MT, MAT_FINAL_ASSEMBLY));
2039:   PetscCall(MatDestroy(&local_mat));
2040:   PetscCall(MatAssemblyEnd(MT, MAT_FINAL_ASSEMBLY));
2041:   if (isseqdense) PetscCall(MatSetOption(MT, MAT_ROW_ORIENTED, PETSC_TRUE));
2042:   if (reuse == MAT_INPLACE_MATRIX) {
2043:     PetscCall(MatHeaderReplace(mat, &MT));
2044:   } else if (reuse == MAT_INITIAL_MATRIX) {
2045:     *M = MT;
2046:   }
2047:   PetscFunctionReturn(PETSC_SUCCESS);
2048: }

2050: /*@
2051:   MatISGetMPIXAIJ - Converts `MATIS` matrix into a parallel `MATAIJ` format

2053:   Input Parameters:
2054: + mat   - the matrix (should be of type `MATIS`)
2055: - reuse - either `MAT_INITIAL_MATRIX` or `MAT_REUSE_MATRIX`

2057:   Output Parameter:
2058: . newmat - the matrix in `MATAIJ` format

2060:   Level: deprecated

2062:   Note:
2063:   This function has been deprecated and it will be removed in future releases. Update your code to use the `MatConvert()` interface.

2065: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatConvert()`
2066: @*/
2067: PetscErrorCode MatISGetMPIXAIJ(Mat mat, MatReuse reuse, Mat *newmat)
2068: {
2069:   PetscFunctionBegin;
2072:   PetscAssertPointer(newmat, 3);
2073:   if (reuse == MAT_REUSE_MATRIX) {
2075:     PetscCheckSameComm(mat, 1, *newmat, 3);
2076:     PetscCheck(mat != *newmat, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse the same matrix");
2077:   }
2078:   PetscUseMethod(mat, "MatISGetMPIXAIJ_C", (Mat, MatType, MatReuse, Mat *), (mat, MATAIJ, reuse, newmat));
2079:   PetscFunctionReturn(PETSC_SUCCESS);
2080: }

2082: static PetscErrorCode MatDuplicate_IS(Mat mat, MatDuplicateOption op, Mat *newmat)
2083: {
2084:   Mat_IS  *matis = (Mat_IS *)(mat->data);
2085:   PetscInt rbs, cbs, m, n, M, N;
2086:   Mat      B, localmat;

2088:   PetscFunctionBegin;
2089:   PetscCall(ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping, &rbs));
2090:   PetscCall(ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping, &cbs));
2091:   PetscCall(MatGetSize(mat, &M, &N));
2092:   PetscCall(MatGetLocalSize(mat, &m, &n));
2093:   PetscCall(MatCreate(PetscObjectComm((PetscObject)mat), &B));
2094:   PetscCall(MatSetSizes(B, m, n, M, N));
2095:   PetscCall(MatSetBlockSize(B, rbs == cbs ? rbs : 1));
2096:   PetscCall(MatSetType(B, MATIS));
2097:   PetscCall(MatISSetLocalMatType(B, matis->lmattype));
2098:   PetscCall(MatSetLocalToGlobalMapping(B, mat->rmap->mapping, mat->cmap->mapping));
2099:   PetscCall(MatDuplicate(matis->A, op, &localmat));
2100:   PetscCall(MatSetLocalToGlobalMapping(localmat, matis->A->rmap->mapping, matis->A->cmap->mapping));
2101:   PetscCall(MatISSetLocalMat(B, localmat));
2102:   PetscCall(MatDestroy(&localmat));
2103:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
2104:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
2105:   *newmat = B;
2106:   PetscFunctionReturn(PETSC_SUCCESS);
2107: }

2109: static PetscErrorCode MatIsHermitian_IS(Mat A, PetscReal tol, PetscBool *flg)
2110: {
2111:   Mat_IS   *matis = (Mat_IS *)A->data;
2112:   PetscBool local_sym;

2114:   PetscFunctionBegin;
2115:   PetscCall(MatIsHermitian(matis->A, tol, &local_sym));
2116:   PetscCall(MPIU_Allreduce(&local_sym, flg, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
2117:   PetscFunctionReturn(PETSC_SUCCESS);
2118: }

2120: static PetscErrorCode MatIsSymmetric_IS(Mat A, PetscReal tol, PetscBool *flg)
2121: {
2122:   Mat_IS   *matis = (Mat_IS *)A->data;
2123:   PetscBool local_sym;

2125:   PetscFunctionBegin;
2126:   if (matis->rmapping != matis->cmapping) {
2127:     *flg = PETSC_FALSE;
2128:     PetscFunctionReturn(PETSC_SUCCESS);
2129:   }
2130:   PetscCall(MatIsSymmetric(matis->A, tol, &local_sym));
2131:   PetscCall(MPIU_Allreduce(&local_sym, flg, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
2132:   PetscFunctionReturn(PETSC_SUCCESS);
2133: }

2135: static PetscErrorCode MatIsStructurallySymmetric_IS(Mat A, PetscBool *flg)
2136: {
2137:   Mat_IS   *matis = (Mat_IS *)A->data;
2138:   PetscBool local_sym;

2140:   PetscFunctionBegin;
2141:   if (matis->rmapping != matis->cmapping) {
2142:     *flg = PETSC_FALSE;
2143:     PetscFunctionReturn(PETSC_SUCCESS);
2144:   }
2145:   PetscCall(MatIsStructurallySymmetric(matis->A, &local_sym));
2146:   PetscCall(MPIU_Allreduce(&local_sym, flg, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
2147:   PetscFunctionReturn(PETSC_SUCCESS);
2148: }

2150: static PetscErrorCode MatDestroy_IS(Mat A)
2151: {
2152:   Mat_IS *b = (Mat_IS *)A->data;

2154:   PetscFunctionBegin;
2155:   PetscCall(PetscFree(b->bdiag));
2156:   PetscCall(PetscFree(b->lmattype));
2157:   PetscCall(MatDestroy(&b->A));
2158:   PetscCall(VecScatterDestroy(&b->cctx));
2159:   PetscCall(VecScatterDestroy(&b->rctx));
2160:   PetscCall(VecDestroy(&b->x));
2161:   PetscCall(VecDestroy(&b->y));
2162:   PetscCall(VecDestroy(&b->counter));
2163:   PetscCall(ISDestroy(&b->getsub_ris));
2164:   PetscCall(ISDestroy(&b->getsub_cis));
2165:   if (b->sf != b->csf) {
2166:     PetscCall(PetscSFDestroy(&b->csf));
2167:     PetscCall(PetscFree2(b->csf_rootdata, b->csf_leafdata));
2168:   } else b->csf = NULL;
2169:   PetscCall(PetscSFDestroy(&b->sf));
2170:   PetscCall(PetscFree2(b->sf_rootdata, b->sf_leafdata));
2171:   PetscCall(ISLocalToGlobalMappingDestroy(&b->rmapping));
2172:   PetscCall(ISLocalToGlobalMappingDestroy(&b->cmapping));
2173:   PetscCall(MatDestroy(&b->dA));
2174:   PetscCall(MatDestroy(&b->assembledA));
2175:   PetscCall(PetscFree(A->data));
2176:   PetscCall(PetscObjectChangeTypeName((PetscObject)A, NULL));
2177:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMatType_C", NULL));
2178:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalMat_C", NULL));
2179:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMat_C", NULL));
2180:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISRestoreLocalMat_C", NULL));
2181:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetMPIXAIJ_C", NULL));
2182:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetPreallocation_C", NULL));
2183:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISStoreL2L_C", NULL));
2184:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISFixLocalEmpty_C", NULL));
2185:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalToGlobalMapping_C", NULL));
2186:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpiaij_C", NULL));
2187:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpibaij_C", NULL));
2188:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpisbaij_C", NULL));
2189:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqaij_C", NULL));
2190:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqbaij_C", NULL));
2191:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqsbaij_C", NULL));
2192:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_aij_C", NULL));
2193:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOOLocal_C", NULL));
2194:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOO_C", NULL));
2195:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", NULL));
2196:   PetscFunctionReturn(PETSC_SUCCESS);
2197: }

2199: static PetscErrorCode MatMult_IS(Mat A, Vec x, Vec y)
2200: {
2201:   Mat_IS     *is   = (Mat_IS *)A->data;
2202:   PetscScalar zero = 0.0;

2204:   PetscFunctionBegin;
2205:   /*  scatter the global vector x into the local work vector */
2206:   PetscCall(VecScatterBegin(is->cctx, x, is->x, INSERT_VALUES, SCATTER_FORWARD));
2207:   PetscCall(VecScatterEnd(is->cctx, x, is->x, INSERT_VALUES, SCATTER_FORWARD));

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

2212:   /* scatter product back into global memory */
2213:   PetscCall(VecSet(y, zero));
2214:   PetscCall(VecScatterBegin(is->rctx, is->y, y, ADD_VALUES, SCATTER_REVERSE));
2215:   PetscCall(VecScatterEnd(is->rctx, is->y, y, ADD_VALUES, SCATTER_REVERSE));
2216:   PetscFunctionReturn(PETSC_SUCCESS);
2217: }

2219: static PetscErrorCode MatMultAdd_IS(Mat A, Vec v1, Vec v2, Vec v3)
2220: {
2221:   Vec temp_vec;

2223:   PetscFunctionBegin; /*  v3 = v2 + A * v1.*/
2224:   if (v3 != v2) {
2225:     PetscCall(MatMult(A, v1, v3));
2226:     PetscCall(VecAXPY(v3, 1.0, v2));
2227:   } else {
2228:     PetscCall(VecDuplicate(v2, &temp_vec));
2229:     PetscCall(MatMult(A, v1, temp_vec));
2230:     PetscCall(VecAXPY(temp_vec, 1.0, v2));
2231:     PetscCall(VecCopy(temp_vec, v3));
2232:     PetscCall(VecDestroy(&temp_vec));
2233:   }
2234:   PetscFunctionReturn(PETSC_SUCCESS);
2235: }

2237: static PetscErrorCode MatMultTranspose_IS(Mat A, Vec y, Vec x)
2238: {
2239:   Mat_IS *is = (Mat_IS *)A->data;

2241:   PetscFunctionBegin;
2242:   /*  scatter the global vector x into the local work vector */
2243:   PetscCall(VecScatterBegin(is->rctx, y, is->y, INSERT_VALUES, SCATTER_FORWARD));
2244:   PetscCall(VecScatterEnd(is->rctx, y, is->y, INSERT_VALUES, SCATTER_FORWARD));

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

2249:   /* scatter product back into global vector */
2250:   PetscCall(VecSet(x, 0));
2251:   PetscCall(VecScatterBegin(is->cctx, is->x, x, ADD_VALUES, SCATTER_REVERSE));
2252:   PetscCall(VecScatterEnd(is->cctx, is->x, x, ADD_VALUES, SCATTER_REVERSE));
2253:   PetscFunctionReturn(PETSC_SUCCESS);
2254: }

2256: static PetscErrorCode MatMultTransposeAdd_IS(Mat A, Vec v1, Vec v2, Vec v3)
2257: {
2258:   Vec temp_vec;

2260:   PetscFunctionBegin; /*  v3 = v2 + A' * v1.*/
2261:   if (v3 != v2) {
2262:     PetscCall(MatMultTranspose(A, v1, v3));
2263:     PetscCall(VecAXPY(v3, 1.0, v2));
2264:   } else {
2265:     PetscCall(VecDuplicate(v2, &temp_vec));
2266:     PetscCall(MatMultTranspose(A, v1, temp_vec));
2267:     PetscCall(VecAXPY(temp_vec, 1.0, v2));
2268:     PetscCall(VecCopy(temp_vec, v3));
2269:     PetscCall(VecDestroy(&temp_vec));
2270:   }
2271:   PetscFunctionReturn(PETSC_SUCCESS);
2272: }

2274: static PetscErrorCode MatView_IS(Mat A, PetscViewer viewer)
2275: {
2276:   Mat_IS     *a = (Mat_IS *)A->data;
2277:   PetscViewer sviewer;
2278:   PetscBool   isascii, view = PETSC_TRUE;

2280:   PetscFunctionBegin;
2281:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2282:   if (isascii) {
2283:     PetscViewerFormat format;

2285:     PetscCall(PetscViewerGetFormat(viewer, &format));
2286:     if (format == PETSC_VIEWER_ASCII_INFO) view = PETSC_FALSE;
2287:   }
2288:   if (!view) PetscFunctionReturn(PETSC_SUCCESS);
2289:   PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
2290:   PetscCall(MatView(a->A, sviewer));
2291:   PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
2292:   PetscCall(PetscViewerFlush(viewer));
2293:   PetscFunctionReturn(PETSC_SUCCESS);
2294: }

2296: static PetscErrorCode MatInvertBlockDiagonal_IS(Mat mat, const PetscScalar **values)
2297: {
2298:   Mat_IS            *is = (Mat_IS *)mat->data;
2299:   MPI_Datatype       nodeType;
2300:   const PetscScalar *lv;
2301:   PetscInt           bs;

2303:   PetscFunctionBegin;
2304:   PetscCall(MatGetBlockSize(mat, &bs));
2305:   PetscCall(MatSetBlockSize(is->A, bs));
2306:   PetscCall(MatInvertBlockDiagonal(is->A, &lv));
2307:   if (!is->bdiag) PetscCall(PetscMalloc1(bs * mat->rmap->n, &is->bdiag));
2308:   PetscCallMPI(MPI_Type_contiguous(bs, MPIU_SCALAR, &nodeType));
2309:   PetscCallMPI(MPI_Type_commit(&nodeType));
2310:   PetscCall(PetscSFReduceBegin(is->sf, nodeType, lv, is->bdiag, MPI_REPLACE));
2311:   PetscCall(PetscSFReduceEnd(is->sf, nodeType, lv, is->bdiag, MPI_REPLACE));
2312:   PetscCallMPI(MPI_Type_free(&nodeType));
2313:   if (values) *values = is->bdiag;
2314:   PetscFunctionReturn(PETSC_SUCCESS);
2315: }

2317: static PetscErrorCode MatISSetUpScatters_Private(Mat A)
2318: {
2319:   Vec             cglobal, rglobal;
2320:   IS              from;
2321:   Mat_IS         *is = (Mat_IS *)A->data;
2322:   PetscScalar     sum;
2323:   const PetscInt *garray;
2324:   PetscInt        nr, rbs, nc, cbs;
2325:   VecType         rtype;

2327:   PetscFunctionBegin;
2328:   PetscCall(ISLocalToGlobalMappingGetSize(is->rmapping, &nr));
2329:   PetscCall(ISLocalToGlobalMappingGetBlockSize(is->rmapping, &rbs));
2330:   PetscCall(ISLocalToGlobalMappingGetSize(is->cmapping, &nc));
2331:   PetscCall(ISLocalToGlobalMappingGetBlockSize(is->cmapping, &cbs));
2332:   PetscCall(VecDestroy(&is->x));
2333:   PetscCall(VecDestroy(&is->y));
2334:   PetscCall(VecDestroy(&is->counter));
2335:   PetscCall(VecScatterDestroy(&is->rctx));
2336:   PetscCall(VecScatterDestroy(&is->cctx));
2337:   PetscCall(MatCreateVecs(is->A, &is->x, &is->y));
2338:   PetscCall(VecBindToCPU(is->y, PETSC_TRUE));
2339:   PetscCall(VecGetRootType_Private(is->y, &rtype));
2340:   PetscCall(PetscFree(A->defaultvectype));
2341:   PetscCall(PetscStrallocpy(rtype, &A->defaultvectype));
2342:   PetscCall(MatCreateVecs(A, &cglobal, &rglobal));
2343:   PetscCall(VecBindToCPU(rglobal, PETSC_TRUE));
2344:   PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->rmapping, &garray));
2345:   PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), rbs, nr / rbs, garray, PETSC_USE_POINTER, &from));
2346:   PetscCall(VecScatterCreate(rglobal, from, is->y, NULL, &is->rctx));
2347:   PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->rmapping, &garray));
2348:   PetscCall(ISDestroy(&from));
2349:   if (is->rmapping != is->cmapping) {
2350:     PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->cmapping, &garray));
2351:     PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), cbs, nc / cbs, garray, PETSC_USE_POINTER, &from));
2352:     PetscCall(VecScatterCreate(cglobal, from, is->x, NULL, &is->cctx));
2353:     PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->cmapping, &garray));
2354:     PetscCall(ISDestroy(&from));
2355:   } else {
2356:     PetscCall(PetscObjectReference((PetscObject)is->rctx));
2357:     is->cctx = is->rctx;
2358:   }
2359:   PetscCall(VecDestroy(&cglobal));

2361:   /* interface counter vector (local) */
2362:   PetscCall(VecDuplicate(is->y, &is->counter));
2363:   PetscCall(VecBindToCPU(is->counter, PETSC_TRUE));
2364:   PetscCall(VecSet(is->y, 1.));
2365:   PetscCall(VecScatterBegin(is->rctx, is->y, rglobal, ADD_VALUES, SCATTER_REVERSE));
2366:   PetscCall(VecScatterEnd(is->rctx, is->y, rglobal, ADD_VALUES, SCATTER_REVERSE));
2367:   PetscCall(VecScatterBegin(is->rctx, rglobal, is->counter, INSERT_VALUES, SCATTER_FORWARD));
2368:   PetscCall(VecScatterEnd(is->rctx, rglobal, is->counter, INSERT_VALUES, SCATTER_FORWARD));
2369:   PetscCall(VecBindToCPU(is->y, PETSC_FALSE));
2370:   PetscCall(VecBindToCPU(is->counter, PETSC_FALSE));

2372:   /* special functions for block-diagonal matrices */
2373:   PetscCall(VecSum(rglobal, &sum));
2374:   A->ops->invertblockdiagonal = NULL;
2375:   if ((PetscInt)(PetscRealPart(sum)) == A->rmap->N && A->rmap->N == A->cmap->N && is->rmapping == is->cmapping) A->ops->invertblockdiagonal = MatInvertBlockDiagonal_IS;
2376:   PetscCall(VecDestroy(&rglobal));

2378:   /* setup SF for general purpose shared indices based communications */
2379:   PetscCall(MatISSetUpSF_IS(A));
2380:   PetscFunctionReturn(PETSC_SUCCESS);
2381: }

2383: static PetscErrorCode MatISFilterL2GMap(Mat A, ISLocalToGlobalMapping map, ISLocalToGlobalMapping *nmap, ISLocalToGlobalMapping *lmap)
2384: {
2385:   IS                         is;
2386:   ISLocalToGlobalMappingType l2gtype;
2387:   const PetscInt            *idxs;
2388:   PetscHSetI                 ht;
2389:   PetscInt                  *nidxs;
2390:   PetscInt                   i, n, bs, c;
2391:   PetscBool                  flg[] = {PETSC_FALSE, PETSC_FALSE};

2393:   PetscFunctionBegin;
2394:   PetscCall(ISLocalToGlobalMappingGetSize(map, &n));
2395:   PetscCall(ISLocalToGlobalMappingGetBlockSize(map, &bs));
2396:   PetscCall(ISLocalToGlobalMappingGetBlockIndices(map, &idxs));
2397:   PetscCall(PetscHSetICreate(&ht));
2398:   PetscCall(PetscMalloc1(n / bs, &nidxs));
2399:   for (i = 0, c = 0; i < n / bs; i++) {
2400:     PetscBool missing;
2401:     if (idxs[i] < 0) {
2402:       flg[0] = PETSC_TRUE;
2403:       continue;
2404:     }
2405:     PetscCall(PetscHSetIQueryAdd(ht, idxs[i], &missing));
2406:     if (!missing) flg[1] = PETSC_TRUE;
2407:     else nidxs[c++] = idxs[i];
2408:   }
2409:   PetscCall(PetscHSetIDestroy(&ht));
2410:   PetscCall(MPIU_Allreduce(MPI_IN_PLACE, flg, 2, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A)));
2411:   if (!flg[0] && !flg[1]) { /* Entries are all non negative and unique */
2412:     *nmap = NULL;
2413:     *lmap = NULL;
2414:     PetscCall(PetscFree(nidxs));
2415:     PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(map, &idxs));
2416:     PetscFunctionReturn(PETSC_SUCCESS);
2417:   }

2419:   /* New l2g map without negative or repeated indices */
2420:   PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), bs, c, nidxs, PETSC_USE_POINTER, &is));
2421:   PetscCall(ISLocalToGlobalMappingCreateIS(is, nmap));
2422:   PetscCall(ISDestroy(&is));
2423:   PetscCall(ISLocalToGlobalMappingGetType(map, &l2gtype));
2424:   PetscCall(ISLocalToGlobalMappingSetType(*nmap, l2gtype));

2426:   /* New local l2g map for repeated indices */
2427:   PetscCall(ISGlobalToLocalMappingApplyBlock(*nmap, IS_GTOLM_MASK, n / bs, idxs, NULL, nidxs));
2428:   PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, n / bs, nidxs, PETSC_USE_POINTER, &is));
2429:   PetscCall(ISLocalToGlobalMappingCreateIS(is, lmap));
2430:   PetscCall(ISDestroy(&is));

2432:   PetscCall(PetscFree(nidxs));
2433:   PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(map, &idxs));
2434:   PetscFunctionReturn(PETSC_SUCCESS);
2435: }

2437: static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A, ISLocalToGlobalMapping rmapping, ISLocalToGlobalMapping cmapping)
2438: {
2439:   Mat_IS                *is            = (Mat_IS *)A->data;
2440:   ISLocalToGlobalMapping localrmapping = NULL, localcmapping = NULL;
2441:   PetscBool              cong, freem[]                       = {PETSC_FALSE, PETSC_FALSE};
2442:   PetscInt               nr, rbs, nc, cbs;

2444:   PetscFunctionBegin;
2445:   if (rmapping) PetscCheckSameComm(A, 1, rmapping, 2);
2446:   if (cmapping) PetscCheckSameComm(A, 1, cmapping, 3);

2448:   PetscCall(ISLocalToGlobalMappingDestroy(&is->rmapping));
2449:   PetscCall(ISLocalToGlobalMappingDestroy(&is->cmapping));
2450:   PetscCall(PetscLayoutSetUp(A->rmap));
2451:   PetscCall(PetscLayoutSetUp(A->cmap));
2452:   PetscCall(MatHasCongruentLayouts(A, &cong));

2454:   /* If NULL, local space matches global space */
2455:   if (!rmapping) {
2456:     IS is;

2458:     PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), A->rmap->N, 0, 1, &is));
2459:     PetscCall(ISLocalToGlobalMappingCreateIS(is, &rmapping));
2460:     if (A->rmap->bs > 0) PetscCall(ISLocalToGlobalMappingSetBlockSize(rmapping, A->rmap->bs));
2461:     PetscCall(ISDestroy(&is));
2462:     freem[0] = PETSC_TRUE;
2463:     if (!cmapping && cong && A->rmap->bs == A->cmap->bs) cmapping = rmapping;
2464:   } else if (!is->islocalref) { /* check if the l2g map has negative or repeated entries */
2465:     PetscCall(MatISFilterL2GMap(A, rmapping, &is->rmapping, &localrmapping));
2466:     if (rmapping == cmapping) {
2467:       PetscCall(PetscObjectReference((PetscObject)is->rmapping));
2468:       is->cmapping = is->rmapping;
2469:       PetscCall(PetscObjectReference((PetscObject)localrmapping));
2470:       localcmapping = localrmapping;
2471:     }
2472:   }
2473:   if (!cmapping) {
2474:     IS is;

2476:     PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), A->cmap->N, 0, 1, &is));
2477:     PetscCall(ISLocalToGlobalMappingCreateIS(is, &cmapping));
2478:     if (A->cmap->bs > 0) PetscCall(ISLocalToGlobalMappingSetBlockSize(cmapping, A->cmap->bs));
2479:     PetscCall(ISDestroy(&is));
2480:     freem[1] = PETSC_TRUE;
2481:   } else if (cmapping != rmapping && !is->islocalref) { /* check if the l2g map has negative or repeated entries */
2482:     PetscCall(MatISFilterL2GMap(A, cmapping, &is->cmapping, &localcmapping));
2483:   }
2484:   if (!is->rmapping) {
2485:     PetscCall(PetscObjectReference((PetscObject)rmapping));
2486:     is->rmapping = rmapping;
2487:   }
2488:   if (!is->cmapping) {
2489:     PetscCall(PetscObjectReference((PetscObject)cmapping));
2490:     is->cmapping = cmapping;
2491:   }

2493:   /* Clean up */
2494:   PetscCall(MatDestroy(&is->A));
2495:   if (is->csf != is->sf) {
2496:     PetscCall(PetscSFDestroy(&is->csf));
2497:     PetscCall(PetscFree2(is->csf_rootdata, is->csf_leafdata));
2498:   } else is->csf = NULL;
2499:   PetscCall(PetscSFDestroy(&is->sf));
2500:   PetscCall(PetscFree2(is->sf_rootdata, is->sf_leafdata));
2501:   PetscCall(PetscFree(is->bdiag));

2503:   /* check if the two mappings are actually the same for square matrices since MATIS has some optimization for this case
2504:      (DOLFIN passes 2 different objects) */
2505:   PetscCall(ISLocalToGlobalMappingGetSize(is->rmapping, &nr));
2506:   PetscCall(ISLocalToGlobalMappingGetBlockSize(is->rmapping, &rbs));
2507:   PetscCall(ISLocalToGlobalMappingGetSize(is->cmapping, &nc));
2508:   PetscCall(ISLocalToGlobalMappingGetBlockSize(is->cmapping, &cbs));
2509:   if (is->rmapping != is->cmapping && cong) {
2510:     PetscBool same = PETSC_FALSE;
2511:     if (nr == nc && cbs == rbs) {
2512:       const PetscInt *idxs1, *idxs2;

2514:       PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->rmapping, &idxs1));
2515:       PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->cmapping, &idxs2));
2516:       PetscCall(PetscArraycmp(idxs1, idxs2, nr / rbs, &same));
2517:       PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->rmapping, &idxs1));
2518:       PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->cmapping, &idxs2));
2519:     }
2520:     PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &same, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
2521:     if (same) {
2522:       PetscCall(ISLocalToGlobalMappingDestroy(&is->cmapping));
2523:       PetscCall(PetscObjectReference((PetscObject)is->rmapping));
2524:       is->cmapping = is->rmapping;
2525:     }
2526:   }
2527:   PetscCall(PetscLayoutSetBlockSize(A->rmap, rbs));
2528:   PetscCall(PetscLayoutSetBlockSize(A->cmap, cbs));
2529:   /* Pass the user defined maps to the layout */
2530:   PetscCall(PetscLayoutSetISLocalToGlobalMapping(A->rmap, rmapping));
2531:   PetscCall(PetscLayoutSetISLocalToGlobalMapping(A->cmap, cmapping));
2532:   if (freem[0]) PetscCall(ISLocalToGlobalMappingDestroy(&rmapping));
2533:   if (freem[1]) PetscCall(ISLocalToGlobalMappingDestroy(&cmapping));

2535:   /* Create the local matrix A */
2536:   PetscCall(MatCreate(PETSC_COMM_SELF, &is->A));
2537:   PetscCall(MatSetType(is->A, is->lmattype));
2538:   PetscCall(MatSetSizes(is->A, nr, nc, nr, nc));
2539:   PetscCall(MatSetBlockSizes(is->A, rbs, cbs));
2540:   PetscCall(MatSetOptionsPrefix(is->A, "is_"));
2541:   PetscCall(MatAppendOptionsPrefix(is->A, ((PetscObject)A)->prefix));
2542:   PetscCall(PetscLayoutSetUp(is->A->rmap));
2543:   PetscCall(PetscLayoutSetUp(is->A->cmap));
2544:   PetscCall(MatSetLocalToGlobalMapping(is->A, localrmapping, localcmapping));
2545:   PetscCall(ISLocalToGlobalMappingDestroy(&localrmapping));
2546:   PetscCall(ISLocalToGlobalMappingDestroy(&localcmapping));

2548:   /* setup scatters and local vectors for MatMult */
2549:   if (!is->islocalref) PetscCall(MatISSetUpScatters_Private(A));
2550:   A->preallocated = PETSC_TRUE;
2551:   PetscFunctionReturn(PETSC_SUCCESS);
2552: }

2554: static PetscErrorCode MatSetUp_IS(Mat A)
2555: {
2556:   ISLocalToGlobalMapping rmap, cmap;

2558:   PetscFunctionBegin;
2559:   PetscCall(MatGetLocalToGlobalMapping(A, &rmap, &cmap));
2560:   if (!rmap && !cmap) PetscCall(MatSetLocalToGlobalMapping(A, NULL, NULL));
2561:   PetscFunctionReturn(PETSC_SUCCESS);
2562: }

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

2569:   PetscFunctionBegin;
2570:   PetscCall(ISGlobalToLocalMappingApply(is->rmapping, IS_GTOLM_MASK, m, rows, &m, rows_l));
2571:   if (m != n || rows != cols || is->cmapping != is->rmapping) {
2572:     PetscCall(ISGlobalToLocalMappingApply(is->cmapping, IS_GTOLM_MASK, n, cols, &n, cols_l));
2573:     PetscCall(MatSetValues(is->A, m, rows_l, n, cols_l, values, addv));
2574:   } else {
2575:     PetscCall(MatSetValues(is->A, m, rows_l, m, rows_l, values, addv));
2576:   }
2577:   PetscFunctionReturn(PETSC_SUCCESS);
2578: }

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

2585:   PetscFunctionBegin;
2586:   PetscCall(ISGlobalToLocalMappingApplyBlock(is->rmapping, IS_GTOLM_MASK, m, rows, &m, rows_l));
2587:   if (m != n || rows != cols || is->cmapping != is->rmapping) {
2588:     PetscCall(ISGlobalToLocalMappingApply(is->cmapping, IS_GTOLM_MASK, n, cols, &n, cols_l));
2589:     PetscCall(MatSetValuesBlocked(is->A, m, rows_l, n, cols_l, values, addv));
2590:   } else {
2591:     PetscCall(MatSetValuesBlocked(is->A, m, rows_l, n, rows_l, values, addv));
2592:   }
2593:   PetscFunctionReturn(PETSC_SUCCESS);
2594: }

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

2600:   PetscFunctionBegin;
2601:   if (is->A->rmap->mapping || is->A->cmap->mapping) {
2602:     PetscCall(MatSetValuesLocal(is->A, m, rows, n, cols, values, addv));
2603:   } else {
2604:     PetscCall(MatSetValues(is->A, m, rows, n, cols, values, addv));
2605:   }
2606:   PetscFunctionReturn(PETSC_SUCCESS);
2607: }

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

2613:   PetscFunctionBegin;
2614:   if (is->A->rmap->mapping || is->A->cmap->mapping) {
2615:     PetscCall(MatSetValuesBlockedLocal(is->A, m, rows, n, cols, values, addv));
2616:   } else {
2617:     PetscCall(MatSetValuesBlocked(is->A, m, rows, n, cols, values, addv));
2618:   }
2619:   PetscFunctionReturn(PETSC_SUCCESS);
2620: }

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

2626:   PetscFunctionBegin;
2627:   if (!n) {
2628:     is->pure_neumann = PETSC_TRUE;
2629:   } else {
2630:     PetscInt i;
2631:     is->pure_neumann = PETSC_FALSE;

2633:     if (columns) {
2634:       PetscCall(MatZeroRowsColumns(is->A, n, rows, diag, NULL, NULL));
2635:     } else {
2636:       PetscCall(MatZeroRows(is->A, n, rows, diag, NULL, NULL));
2637:     }
2638:     if (diag != 0.) {
2639:       const PetscScalar *array;
2640:       PetscCall(VecGetArrayRead(is->counter, &array));
2641:       for (i = 0; i < n; i++) PetscCall(MatSetValue(is->A, rows[i], rows[i], diag / (array[rows[i]]), INSERT_VALUES));
2642:       PetscCall(VecRestoreArrayRead(is->counter, &array));
2643:     }
2644:     PetscCall(MatAssemblyBegin(is->A, MAT_FINAL_ASSEMBLY));
2645:     PetscCall(MatAssemblyEnd(is->A, MAT_FINAL_ASSEMBLY));
2646:   }
2647:   PetscFunctionReturn(PETSC_SUCCESS);
2648: }

2650: static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b, PetscBool columns)
2651: {
2652:   Mat_IS   *matis = (Mat_IS *)A->data;
2653:   PetscInt  nr, nl, len, i;
2654:   PetscInt *lrows;

2656:   PetscFunctionBegin;
2657:   if (PetscUnlikelyDebug(columns || diag != 0. || (x && b))) {
2658:     PetscBool cong;

2660:     PetscCall(PetscLayoutCompare(A->rmap, A->cmap, &cong));
2661:     cong = (PetscBool)(cong && matis->sf == matis->csf);
2662:     PetscCheck(cong || !columns, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Columns can be zeroed if and only if A->rmap and A->cmap are congruent and the l2g maps are the same for MATIS");
2663:     PetscCheck(cong || diag == 0., PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Nonzero diagonal value supported if and only if A->rmap and A->cmap are congruent and the l2g maps are the same for MATIS");
2664:     PetscCheck(cong || !x || !b, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "A->rmap and A->cmap need to be congruent, and the l2g maps be the same");
2665:   }
2666:   /* get locally owned rows */
2667:   PetscCall(PetscLayoutMapLocal(A->rmap, n, rows, &len, &lrows, NULL));
2668:   /* fix right hand side if needed */
2669:   if (x && b) {
2670:     const PetscScalar *xx;
2671:     PetscScalar       *bb;

2673:     PetscCall(VecGetArrayRead(x, &xx));
2674:     PetscCall(VecGetArray(b, &bb));
2675:     for (i = 0; i < len; ++i) bb[lrows[i]] = diag * xx[lrows[i]];
2676:     PetscCall(VecRestoreArrayRead(x, &xx));
2677:     PetscCall(VecRestoreArray(b, &bb));
2678:   }
2679:   /* get rows associated to the local matrices */
2680:   PetscCall(MatGetSize(matis->A, &nl, NULL));
2681:   PetscCall(PetscArrayzero(matis->sf_leafdata, nl));
2682:   PetscCall(PetscArrayzero(matis->sf_rootdata, A->rmap->n));
2683:   for (i = 0; i < len; i++) matis->sf_rootdata[lrows[i]] = 1;
2684:   PetscCall(PetscFree(lrows));
2685:   PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
2686:   PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
2687:   PetscCall(PetscMalloc1(nl, &lrows));
2688:   for (i = 0, nr = 0; i < nl; i++)
2689:     if (matis->sf_leafdata[i]) lrows[nr++] = i;
2690:   PetscCall(MatISZeroRowsColumnsLocal_Private(A, nr, lrows, diag, columns));
2691:   PetscCall(PetscFree(lrows));
2692:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
2693:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
2694:   PetscFunctionReturn(PETSC_SUCCESS);
2695: }

2697: static PetscErrorCode MatZeroRows_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
2698: {
2699:   PetscFunctionBegin;
2700:   PetscCall(MatZeroRowsColumns_Private_IS(A, n, rows, diag, x, b, PETSC_FALSE));
2701:   PetscFunctionReturn(PETSC_SUCCESS);
2702: }

2704: static PetscErrorCode MatZeroRowsColumns_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
2705: {
2706:   PetscFunctionBegin;
2707:   PetscCall(MatZeroRowsColumns_Private_IS(A, n, rows, diag, x, b, PETSC_TRUE));
2708:   PetscFunctionReturn(PETSC_SUCCESS);
2709: }

2711: static PetscErrorCode MatAssemblyBegin_IS(Mat A, MatAssemblyType type)
2712: {
2713:   Mat_IS *is = (Mat_IS *)A->data;

2715:   PetscFunctionBegin;
2716:   PetscCall(MatAssemblyBegin(is->A, type));
2717:   PetscFunctionReturn(PETSC_SUCCESS);
2718: }

2720: static PetscErrorCode MatAssemblyEnd_IS(Mat A, MatAssemblyType type)
2721: {
2722:   Mat_IS   *is = (Mat_IS *)A->data;
2723:   PetscBool lnnz;

2725:   PetscFunctionBegin;
2726:   PetscCall(MatAssemblyEnd(is->A, type));
2727:   /* fix for local empty rows/cols */
2728:   if (is->locempty && type == MAT_FINAL_ASSEMBLY) {
2729:     Mat                    newlA;
2730:     ISLocalToGlobalMapping rl2g, cl2g;
2731:     IS                     nzr, nzc;
2732:     PetscInt               nr, nc, nnzr, nnzc;
2733:     PetscBool              lnewl2g, newl2g;

2735:     PetscCall(MatGetSize(is->A, &nr, &nc));
2736:     PetscCall(MatFindNonzeroRowsOrCols_Basic(is->A, PETSC_FALSE, PETSC_SMALL, &nzr));
2737:     if (!nzr) PetscCall(ISCreateStride(PetscObjectComm((PetscObject)is->A), nr, 0, 1, &nzr));
2738:     PetscCall(MatFindNonzeroRowsOrCols_Basic(is->A, PETSC_TRUE, PETSC_SMALL, &nzc));
2739:     if (!nzc) PetscCall(ISCreateStride(PetscObjectComm((PetscObject)is->A), nc, 0, 1, &nzc));
2740:     PetscCall(ISGetSize(nzr, &nnzr));
2741:     PetscCall(ISGetSize(nzc, &nnzc));
2742:     if (nnzr != nr || nnzc != nc) { /* need new global l2g map */
2743:       lnewl2g = PETSC_TRUE;
2744:       PetscCall(MPIU_Allreduce(&lnewl2g, &newl2g, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A)));

2746:       /* extract valid submatrix */
2747:       PetscCall(MatCreateSubMatrix(is->A, nzr, nzc, MAT_INITIAL_MATRIX, &newlA));
2748:     } else { /* local matrix fully populated */
2749:       lnewl2g = PETSC_FALSE;
2750:       PetscCall(MPIU_Allreduce(&lnewl2g, &newl2g, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A)));
2751:       PetscCall(PetscObjectReference((PetscObject)is->A));
2752:       newlA = is->A;
2753:     }

2755:     /* attach new global l2g map if needed */
2756:     if (newl2g) {
2757:       IS              zr, zc;
2758:       const PetscInt *ridxs, *cidxs, *zridxs, *zcidxs;
2759:       PetscInt       *nidxs, i;

2761:       PetscCall(ISComplement(nzr, 0, nr, &zr));
2762:       PetscCall(ISComplement(nzc, 0, nc, &zc));
2763:       PetscCall(PetscMalloc1(PetscMax(nr, nc), &nidxs));
2764:       PetscCall(ISLocalToGlobalMappingGetIndices(is->rmapping, &ridxs));
2765:       PetscCall(ISLocalToGlobalMappingGetIndices(is->cmapping, &cidxs));
2766:       PetscCall(ISGetIndices(zr, &zridxs));
2767:       PetscCall(ISGetIndices(zc, &zcidxs));
2768:       PetscCall(ISGetLocalSize(zr, &nnzr));
2769:       PetscCall(ISGetLocalSize(zc, &nnzc));

2771:       PetscCall(PetscArraycpy(nidxs, ridxs, nr));
2772:       for (i = 0; i < nnzr; i++) nidxs[zridxs[i]] = -1;
2773:       PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A), 1, nr, nidxs, PETSC_COPY_VALUES, &rl2g));
2774:       PetscCall(PetscArraycpy(nidxs, cidxs, nc));
2775:       for (i = 0; i < nnzc; i++) nidxs[zcidxs[i]] = -1;
2776:       PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A), 1, nc, nidxs, PETSC_COPY_VALUES, &cl2g));

2778:       PetscCall(ISRestoreIndices(zr, &zridxs));
2779:       PetscCall(ISRestoreIndices(zc, &zcidxs));
2780:       PetscCall(ISLocalToGlobalMappingRestoreIndices(is->rmapping, &ridxs));
2781:       PetscCall(ISLocalToGlobalMappingRestoreIndices(is->cmapping, &cidxs));
2782:       PetscCall(ISDestroy(&nzr));
2783:       PetscCall(ISDestroy(&nzc));
2784:       PetscCall(ISDestroy(&zr));
2785:       PetscCall(ISDestroy(&zc));
2786:       PetscCall(PetscFree(nidxs));
2787:       PetscCall(MatSetLocalToGlobalMapping(A, rl2g, cl2g));
2788:       PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
2789:       PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
2790:     }
2791:     PetscCall(MatISSetLocalMat(A, newlA));
2792:     PetscCall(MatDestroy(&newlA));
2793:     PetscCall(ISDestroy(&nzr));
2794:     PetscCall(ISDestroy(&nzc));
2795:     is->locempty = PETSC_FALSE;
2796:   }
2797:   lnnz          = (PetscBool)(is->A->nonzerostate == is->lnnzstate);
2798:   is->lnnzstate = is->A->nonzerostate;
2799:   PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &lnnz, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A)));
2800:   if (lnnz) A->nonzerostate++;
2801:   PetscFunctionReturn(PETSC_SUCCESS);
2802: }

2804: static PetscErrorCode MatISGetLocalMat_IS(Mat mat, Mat *local)
2805: {
2806:   Mat_IS *is = (Mat_IS *)mat->data;

2808:   PetscFunctionBegin;
2809:   *local = is->A;
2810:   PetscFunctionReturn(PETSC_SUCCESS);
2811: }

2813: static PetscErrorCode MatISRestoreLocalMat_IS(Mat mat, Mat *local)
2814: {
2815:   PetscFunctionBegin;
2816:   *local = NULL;
2817:   PetscFunctionReturn(PETSC_SUCCESS);
2818: }

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

2823:   Input Parameter:
2824: . mat - the matrix

2826:   Output Parameter:
2827: . local - the local matrix

2829:   Level: advanced

2831:   Notes:
2832:   This can be called if you have precomputed the nonzero structure of the
2833:   matrix and want to provide it to the inner matrix object to improve the performance
2834:   of the `MatSetValues()` operation.

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

2838: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatISRestoreLocalMat()`
2839: @*/
2840: PetscErrorCode MatISGetLocalMat(Mat mat, Mat *local)
2841: {
2842:   PetscFunctionBegin;
2844:   PetscAssertPointer(local, 2);
2845:   PetscUseMethod(mat, "MatISGetLocalMat_C", (Mat, Mat *), (mat, local));
2846:   PetscFunctionReturn(PETSC_SUCCESS);
2847: }

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

2852:   Input Parameters:
2853: + mat   - the matrix
2854: - local - the local matrix

2856:   Level: advanced

2858: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatISGetLocalMat()`
2859: @*/
2860: PetscErrorCode MatISRestoreLocalMat(Mat mat, Mat *local)
2861: {
2862:   PetscFunctionBegin;
2864:   PetscAssertPointer(local, 2);
2865:   PetscUseMethod(mat, "MatISRestoreLocalMat_C", (Mat, Mat *), (mat, local));
2866:   PetscFunctionReturn(PETSC_SUCCESS);
2867: }

2869: static PetscErrorCode MatISSetLocalMatType_IS(Mat mat, MatType mtype)
2870: {
2871:   Mat_IS *is = (Mat_IS *)mat->data;

2873:   PetscFunctionBegin;
2874:   if (is->A) PetscCall(MatSetType(is->A, mtype));
2875:   PetscCall(PetscFree(is->lmattype));
2876:   PetscCall(PetscStrallocpy(mtype, &is->lmattype));
2877:   PetscFunctionReturn(PETSC_SUCCESS);
2878: }

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

2883:   Input Parameters:
2884: + mat   - the matrix
2885: - mtype - the local matrix type

2887:   Level: advanced

2889: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatSetType()`, `MatType`
2890: @*/
2891: PetscErrorCode MatISSetLocalMatType(Mat mat, MatType mtype)
2892: {
2893:   PetscFunctionBegin;
2895:   PetscUseMethod(mat, "MatISSetLocalMatType_C", (Mat, MatType), (mat, mtype));
2896:   PetscFunctionReturn(PETSC_SUCCESS);
2897: }

2899: static PetscErrorCode MatISSetLocalMat_IS(Mat mat, Mat local)
2900: {
2901:   Mat_IS   *is = (Mat_IS *)mat->data;
2902:   PetscInt  nrows, ncols, orows, ocols;
2903:   MatType   mtype, otype;
2904:   PetscBool sametype = PETSC_TRUE;

2906:   PetscFunctionBegin;
2907:   if (is->A && !is->islocalref) {
2908:     PetscCall(MatGetSize(is->A, &orows, &ocols));
2909:     PetscCall(MatGetSize(local, &nrows, &ncols));
2910:     PetscCheck(orows == nrows && ocols == ncols, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local MATIS matrix should be of size %" PetscInt_FMT "x%" PetscInt_FMT " (you passed a %" PetscInt_FMT "x%" PetscInt_FMT " matrix)", orows, ocols, nrows, ncols);
2911:     PetscCall(MatGetType(local, &mtype));
2912:     PetscCall(MatGetType(is->A, &otype));
2913:     PetscCall(PetscStrcmp(mtype, otype, &sametype));
2914:   }
2915:   PetscCall(PetscObjectReference((PetscObject)local));
2916:   PetscCall(MatDestroy(&is->A));
2917:   is->A = local;
2918:   PetscCall(MatGetType(is->A, &mtype));
2919:   PetscCall(MatISSetLocalMatType(mat, mtype));
2920:   if (!sametype && !is->islocalref) PetscCall(MatISSetUpScatters_Private(mat));
2921:   PetscFunctionReturn(PETSC_SUCCESS);
2922: }

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

2927:   Collective

2929:   Input Parameters:
2930: + mat   - the matrix
2931: - local - the local matrix

2933:   Level: advanced

2935:   Notes:
2936:   Any previous matrix within the `MATIS` has its reference count decreased by one.

2938:   This can be called if you have precomputed the local matrix and
2939:   want to provide it to the matrix object `MATIS`.

2941: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatISSetLocalMatType`, `MatISGetLocalMat()`
2942: @*/
2943: PetscErrorCode MatISSetLocalMat(Mat mat, Mat local)
2944: {
2945:   PetscFunctionBegin;
2948:   PetscUseMethod(mat, "MatISSetLocalMat_C", (Mat, Mat), (mat, local));
2949:   PetscFunctionReturn(PETSC_SUCCESS);
2950: }

2952: static PetscErrorCode MatZeroEntries_IS(Mat A)
2953: {
2954:   Mat_IS *a = (Mat_IS *)A->data;

2956:   PetscFunctionBegin;
2957:   PetscCall(MatZeroEntries(a->A));
2958:   PetscFunctionReturn(PETSC_SUCCESS);
2959: }

2961: static PetscErrorCode MatScale_IS(Mat A, PetscScalar a)
2962: {
2963:   Mat_IS *is = (Mat_IS *)A->data;

2965:   PetscFunctionBegin;
2966:   PetscCall(MatScale(is->A, a));
2967:   PetscFunctionReturn(PETSC_SUCCESS);
2968: }

2970: static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v)
2971: {
2972:   Mat_IS *is = (Mat_IS *)A->data;

2974:   PetscFunctionBegin;
2975:   /* get diagonal of the local matrix */
2976:   PetscCall(MatGetDiagonal(is->A, is->y));

2978:   /* scatter diagonal back into global vector */
2979:   PetscCall(VecSet(v, 0));
2980:   PetscCall(VecScatterBegin(is->rctx, is->y, v, ADD_VALUES, SCATTER_REVERSE));
2981:   PetscCall(VecScatterEnd(is->rctx, is->y, v, ADD_VALUES, SCATTER_REVERSE));
2982:   PetscFunctionReturn(PETSC_SUCCESS);
2983: }

2985: static PetscErrorCode MatSetOption_IS(Mat A, MatOption op, PetscBool flg)
2986: {
2987:   Mat_IS *a = (Mat_IS *)A->data;

2989:   PetscFunctionBegin;
2990:   PetscCall(MatSetOption(a->A, op, flg));
2991:   PetscFunctionReturn(PETSC_SUCCESS);
2992: }

2994: static PetscErrorCode MatAXPY_IS(Mat Y, PetscScalar a, Mat X, MatStructure str)
2995: {
2996:   Mat_IS *y = (Mat_IS *)Y->data;
2997:   Mat_IS *x;

2999:   PetscFunctionBegin;
3000:   if (PetscDefined(USE_DEBUG)) {
3001:     PetscBool ismatis;
3002:     PetscCall(PetscObjectTypeCompare((PetscObject)X, MATIS, &ismatis));
3003:     PetscCheck(ismatis, PetscObjectComm((PetscObject)Y), PETSC_ERR_SUP, "Cannot call MatAXPY(Y,a,X,str) with X not of type MATIS");
3004:   }
3005:   x = (Mat_IS *)X->data;
3006:   PetscCall(MatAXPY(y->A, a, x->A, str));
3007:   PetscFunctionReturn(PETSC_SUCCESS);
3008: }

3010: static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A, IS row, IS col, Mat *submat)
3011: {
3012:   Mat                    lA;
3013:   Mat_IS                *matis = (Mat_IS *)(A->data);
3014:   ISLocalToGlobalMapping rl2g, cl2g;
3015:   IS                     is;
3016:   const PetscInt        *rg, *rl;
3017:   PetscInt               nrg;
3018:   PetscInt               N, M, nrl, i, *idxs;

3020:   PetscFunctionBegin;
3021:   PetscCall(ISLocalToGlobalMappingGetIndices(A->rmap->mapping, &rg));
3022:   PetscCall(ISGetLocalSize(row, &nrl));
3023:   PetscCall(ISGetIndices(row, &rl));
3024:   PetscCall(ISLocalToGlobalMappingGetSize(A->rmap->mapping, &nrg));
3025:   if (PetscDefined(USE_DEBUG)) {
3026:     for (i = 0; i < nrl; i++) PetscCheck(rl[i] < nrg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Local row index %" PetscInt_FMT " -> %" PetscInt_FMT " greater then maximum possible %" PetscInt_FMT, i, rl[i], nrg);
3027:   }
3028:   PetscCall(PetscMalloc1(nrg, &idxs));
3029:   /* map from [0,nrl) to row */
3030:   for (i = 0; i < nrl; i++) idxs[i] = rl[i];
3031:   for (i = nrl; i < nrg; i++) idxs[i] = -1;
3032:   PetscCall(ISRestoreIndices(row, &rl));
3033:   PetscCall(ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping, &rg));
3034:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A), nrg, idxs, PETSC_OWN_POINTER, &is));
3035:   PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g));
3036:   PetscCall(ISDestroy(&is));
3037:   /* compute new l2g map for columns */
3038:   if (col != row || matis->rmapping != matis->cmapping || matis->A->rmap->mapping != matis->A->cmap->mapping) {
3039:     const PetscInt *cg, *cl;
3040:     PetscInt        ncg;
3041:     PetscInt        ncl;

3043:     PetscCall(ISLocalToGlobalMappingGetIndices(A->cmap->mapping, &cg));
3044:     PetscCall(ISGetLocalSize(col, &ncl));
3045:     PetscCall(ISGetIndices(col, &cl));
3046:     PetscCall(ISLocalToGlobalMappingGetSize(A->cmap->mapping, &ncg));
3047:     if (PetscDefined(USE_DEBUG)) {
3048:       for (i = 0; i < ncl; i++) PetscCheck(cl[i] < ncg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Local column index %" PetscInt_FMT " -> %" PetscInt_FMT " greater then maximum possible %" PetscInt_FMT, i, cl[i], ncg);
3049:     }
3050:     PetscCall(PetscMalloc1(ncg, &idxs));
3051:     /* map from [0,ncl) to col */
3052:     for (i = 0; i < ncl; i++) idxs[i] = cl[i];
3053:     for (i = ncl; i < ncg; i++) idxs[i] = -1;
3054:     PetscCall(ISRestoreIndices(col, &cl));
3055:     PetscCall(ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping, &cg));
3056:     PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A), ncg, idxs, PETSC_OWN_POINTER, &is));
3057:     PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g));
3058:     PetscCall(ISDestroy(&is));
3059:   } else {
3060:     PetscCall(PetscObjectReference((PetscObject)rl2g));
3061:     cl2g = rl2g;
3062:   }
3063:   /* create the MATIS submatrix */
3064:   PetscCall(MatGetSize(A, &M, &N));
3065:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), submat));
3066:   PetscCall(MatSetSizes(*submat, PETSC_DECIDE, PETSC_DECIDE, M, N));
3067:   PetscCall(MatSetType(*submat, MATIS));
3068:   matis             = (Mat_IS *)((*submat)->data);
3069:   matis->islocalref = PETSC_TRUE;
3070:   PetscCall(MatSetLocalToGlobalMapping(*submat, rl2g, cl2g));
3071:   PetscCall(MatISGetLocalMat(A, &lA));
3072:   PetscCall(MatISSetLocalMat(*submat, lA));
3073:   PetscCall(MatSetUp(*submat));
3074:   PetscCall(MatAssemblyBegin(*submat, MAT_FINAL_ASSEMBLY));
3075:   PetscCall(MatAssemblyEnd(*submat, MAT_FINAL_ASSEMBLY));
3076:   PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
3077:   PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));

3079:   /* remove unsupported ops */
3080:   PetscCall(PetscMemzero((*submat)->ops, sizeof(struct _MatOps)));
3081:   (*submat)->ops->destroy               = MatDestroy_IS;
3082:   (*submat)->ops->setvalueslocal        = MatSetValuesLocal_SubMat_IS;
3083:   (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS;
3084:   (*submat)->ops->assemblybegin         = MatAssemblyBegin_IS;
3085:   (*submat)->ops->assemblyend           = MatAssemblyEnd_IS;
3086:   PetscFunctionReturn(PETSC_SUCCESS);
3087: }

3089: static PetscErrorCode MatSetFromOptions_IS(Mat A, PetscOptionItems *PetscOptionsObject)
3090: {
3091:   Mat_IS   *a = (Mat_IS *)A->data;
3092:   char      type[256];
3093:   PetscBool flg;

3095:   PetscFunctionBegin;
3096:   PetscOptionsHeadBegin(PetscOptionsObject, "MATIS options");
3097:   PetscCall(PetscOptionsBool("-matis_keepassembled", "Store an assembled version if needed", "MatISKeepAssembled", a->keepassembled, &a->keepassembled, NULL));
3098:   PetscCall(PetscOptionsBool("-matis_fixempty", "Fix local matrices in case of empty local rows/columns", "MatISFixLocalEmpty", a->locempty, &a->locempty, NULL));
3099:   PetscCall(PetscOptionsBool("-matis_storel2l", "Store local-to-local matrices generated from PtAP operations", "MatISStoreL2L", a->storel2l, &a->storel2l, NULL));
3100:   PetscCall(PetscOptionsFList("-matis_localmat_type", "Matrix type", "MatISSetLocalMatType", MatList, a->lmattype, type, 256, &flg));
3101:   if (flg) PetscCall(MatISSetLocalMatType(A, type));
3102:   if (a->A) PetscCall(MatSetFromOptions(a->A));
3103:   PetscOptionsHeadEnd();
3104:   PetscFunctionReturn(PETSC_SUCCESS);
3105: }

3107: /*@
3108:   MatCreateIS - Creates a "process" unassembled matrix, `MATIS`, assembled on each
3109:   process but not across processes.

3111:   Input Parameters:
3112: + comm - MPI communicator that will share the matrix
3113: . bs   - block size of the matrix
3114: . m    - local size of left vector used in matrix vector products
3115: . n    - local size of right vector used in matrix vector products
3116: . M    - global size of left vector used in matrix vector products
3117: . N    - global size of right vector used in matrix vector products
3118: . rmap - local to global map for rows
3119: - cmap - local to global map for cols

3121:   Output Parameter:
3122: . A - the resulting matrix

3124:   Level: advanced

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

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

3132: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatSetLocalToGlobalMapping()`
3133: @*/
3134: PetscErrorCode MatCreateIS(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt M, PetscInt N, ISLocalToGlobalMapping rmap, ISLocalToGlobalMapping cmap, Mat *A)
3135: {
3136:   PetscFunctionBegin;
3137:   PetscCall(MatCreate(comm, A));
3138:   PetscCall(MatSetSizes(*A, m, n, M, N));
3139:   if (bs > 0) PetscCall(MatSetBlockSize(*A, bs));
3140:   PetscCall(MatSetType(*A, MATIS));
3141:   PetscCall(MatSetLocalToGlobalMapping(*A, rmap, cmap));
3142:   PetscFunctionReturn(PETSC_SUCCESS);
3143: }

3145: static PetscErrorCode MatHasOperation_IS(Mat A, MatOperation op, PetscBool *has)
3146: {
3147:   Mat_IS      *a              = (Mat_IS *)A->data;
3148:   MatOperation tobefiltered[] = {MATOP_MULT_ADD, MATOP_MULT_TRANSPOSE_ADD, MATOP_GET_DIAGONAL_BLOCK, MATOP_INCREASE_OVERLAP};

3150:   PetscFunctionBegin;
3151:   *has = PETSC_FALSE;
3152:   if (!((void **)A->ops)[op] || !a->A) PetscFunctionReturn(PETSC_SUCCESS);
3153:   *has = PETSC_TRUE;
3154:   for (PetscInt i = 0; i < (PetscInt)PETSC_STATIC_ARRAY_LENGTH(tobefiltered); i++)
3155:     if (op == tobefiltered[i]) PetscFunctionReturn(PETSC_SUCCESS);
3156:   PetscCall(MatHasOperation(a->A, op, has));
3157:   PetscFunctionReturn(PETSC_SUCCESS);
3158: }

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

3164:   PetscFunctionBegin;
3165:   PetscCall(MatSetValuesCOO(a->A, v, imode));
3166:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
3167:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
3168:   PetscFunctionReturn(PETSC_SUCCESS);
3169: }

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

3175:   PetscFunctionBegin;
3176:   PetscCheck(a->A, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to provide l2g map first via MatSetLocalToGlobalMapping");
3177:   if (a->A->rmap->mapping || a->A->cmap->mapping) {
3178:     PetscCall(MatSetPreallocationCOOLocal(a->A, ncoo, coo_i, coo_j));
3179:   } else {
3180:     PetscCall(MatSetPreallocationCOO(a->A, ncoo, coo_i, coo_j));
3181:   }
3182:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", MatSetValuesCOO_IS));
3183:   A->preallocated = PETSC_TRUE;
3184:   PetscFunctionReturn(PETSC_SUCCESS);
3185: }

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

3191:   PetscFunctionBegin;
3192:   PetscCheck(a->A, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to provide l2g map first via MatSetLocalToGlobalMapping");
3193:   PetscCheck(ncoo <= PETSC_MAX_INT, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "ncoo %" PetscCount_FMT " overflowed PetscInt; configure --with-64-bit-indices or request support", ncoo);
3194:   PetscCall(ISGlobalToLocalMappingApply(a->rmapping, IS_GTOLM_MASK, ncoo, coo_i, NULL, coo_i));
3195:   PetscCall(ISGlobalToLocalMappingApply(a->cmapping, IS_GTOLM_MASK, ncoo, coo_j, NULL, coo_j));
3196:   PetscCall(MatSetPreallocationCOO(a->A, ncoo, coo_i, coo_j));
3197:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", MatSetValuesCOO_IS));
3198:   A->preallocated = PETSC_TRUE;
3199:   PetscFunctionReturn(PETSC_SUCCESS);
3200: }

3202: static PetscErrorCode MatISGetAssembled_Private(Mat A, Mat *tA)
3203: {
3204:   Mat_IS          *a = (Mat_IS *)A->data;
3205:   PetscObjectState Astate, aAstate       = PETSC_MIN_INT;
3206:   PetscObjectState Annzstate, aAnnzstate = PETSC_MIN_INT;

3208:   PetscFunctionBegin;
3209:   PetscCall(PetscObjectStateGet((PetscObject)A, &Astate));
3210:   Annzstate = A->nonzerostate;
3211:   if (a->assembledA) {
3212:     PetscCall(PetscObjectStateGet((PetscObject)a->assembledA, &aAstate));
3213:     aAnnzstate = a->assembledA->nonzerostate;
3214:   }
3215:   if (aAnnzstate != Annzstate) PetscCall(MatDestroy(&a->assembledA));
3216:   if (Astate != aAstate || !a->assembledA) {
3217:     MatType     aAtype;
3218:     PetscMPIInt size;
3219:     PetscInt    rbs, cbs, bs;

3221:     /* the assembled form is used as temporary storage for parallel operations
3222:        like createsubmatrices and the like, do not waste device memory */
3223:     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
3224:     PetscCall(ISLocalToGlobalMappingGetBlockSize(a->cmapping, &cbs));
3225:     PetscCall(ISLocalToGlobalMappingGetBlockSize(a->rmapping, &rbs));
3226:     bs = rbs == cbs ? rbs : 1;
3227:     if (a->assembledA) PetscCall(MatGetType(a->assembledA, &aAtype));
3228:     else if (size > 1) aAtype = bs > 1 ? MATMPIBAIJ : MATMPIAIJ;
3229:     else aAtype = bs > 1 ? MATSEQBAIJ : MATSEQAIJ;

3231:     PetscCall(MatConvert(A, aAtype, a->assembledA ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX, &a->assembledA));
3232:     PetscCall(PetscObjectStateSet((PetscObject)a->assembledA, Astate));
3233:     a->assembledA->nonzerostate = Annzstate;
3234:   }
3235:   PetscCall(PetscObjectReference((PetscObject)a->assembledA));
3236:   *tA = a->assembledA;
3237:   if (!a->keepassembled) PetscCall(MatDestroy(&a->assembledA));
3238:   PetscFunctionReturn(PETSC_SUCCESS);
3239: }

3241: static PetscErrorCode MatISRestoreAssembled_Private(Mat A, Mat *tA)
3242: {
3243:   PetscFunctionBegin;
3244:   PetscCall(MatDestroy(tA));
3245:   *tA = NULL;
3246:   PetscFunctionReturn(PETSC_SUCCESS);
3247: }

3249: static PetscErrorCode MatGetDiagonalBlock_IS(Mat A, Mat *dA)
3250: {
3251:   Mat_IS          *a = (Mat_IS *)A->data;
3252:   PetscObjectState Astate, dAstate = PETSC_MIN_INT;

3254:   PetscFunctionBegin;
3255:   PetscCall(PetscObjectStateGet((PetscObject)A, &Astate));
3256:   if (a->dA) PetscCall(PetscObjectStateGet((PetscObject)a->dA, &dAstate));
3257:   if (Astate != dAstate) {
3258:     Mat     tA;
3259:     MatType ltype;

3261:     PetscCall(MatDestroy(&a->dA));
3262:     PetscCall(MatISGetAssembled_Private(A, &tA));
3263:     PetscCall(MatGetDiagonalBlock(tA, &a->dA));
3264:     PetscCall(MatPropagateSymmetryOptions(tA, a->dA));
3265:     PetscCall(MatGetType(a->A, &ltype));
3266:     PetscCall(MatConvert(a->dA, ltype, MAT_INPLACE_MATRIX, &a->dA));
3267:     PetscCall(PetscObjectReference((PetscObject)a->dA));
3268:     PetscCall(MatISRestoreAssembled_Private(A, &tA));
3269:     PetscCall(PetscObjectStateSet((PetscObject)a->dA, Astate));
3270:   }
3271:   *dA = a->dA;
3272:   PetscFunctionReturn(PETSC_SUCCESS);
3273: }

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

3279:   PetscFunctionBegin;
3280:   PetscCall(MatISGetAssembled_Private(A, &tA));
3281:   PetscCall(MatCreateSubMatrices(tA, n, irow, icol, reuse, submat));
3282:   /* MatCreateSubMatrices_MPIAIJ is a mess at the moment */
3283: #if 0
3284:   {
3285:     Mat_IS    *a = (Mat_IS*)A->data;
3286:     MatType   ltype;
3287:     VecType   vtype;
3288:     char      *flg;

3290:     PetscCall(MatGetType(a->A,&ltype));
3291:     PetscCall(MatGetVecType(a->A,&vtype));
3292:     PetscCall(PetscStrstr(vtype,"cuda",&flg));
3293:     if (!flg) PetscCall(PetscStrstr(vtype,"hip",&flg));
3294:     if (!flg) PetscCall(PetscStrstr(vtype,"kokkos",&flg));
3295:     if (flg) {
3296:       for (PetscInt i = 0; i < n; i++) {
3297:         Mat sA = (*submat)[i];

3299:         PetscCall(MatConvert(sA,ltype,MAT_INPLACE_MATRIX,&sA));
3300:         (*submat)[i] = sA;
3301:       }
3302:     }
3303:   }
3304: #endif
3305:   PetscCall(MatISRestoreAssembled_Private(A, &tA));
3306:   PetscFunctionReturn(PETSC_SUCCESS);
3307: }

3309: static PetscErrorCode MatIncreaseOverlap_IS(Mat A, PetscInt n, IS is[], PetscInt ov)
3310: {
3311:   Mat tA;

3313:   PetscFunctionBegin;
3314:   PetscCall(MatISGetAssembled_Private(A, &tA));
3315:   PetscCall(MatIncreaseOverlap(tA, n, is, ov));
3316:   PetscCall(MatISRestoreAssembled_Private(A, &tA));
3317:   PetscFunctionReturn(PETSC_SUCCESS);
3318: }

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

3323:   Not Collective

3325:   Input Parameter:
3326: . A - the matrix

3328:   Output Parameters:
3329: + rmapping - row mapping
3330: - cmapping - column mapping

3332:   Level: advanced

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

3337: .seealso: [](ch_matrices), `Mat`, `MatIS`, `MatSetLocalToGlobalMapping()`
3338: @*/
3339: PetscErrorCode MatISGetLocalToGlobalMapping(Mat A, ISLocalToGlobalMapping *rmapping, ISLocalToGlobalMapping *cmapping)
3340: {
3341:   PetscFunctionBegin;
3344:   if (rmapping) PetscAssertPointer(rmapping, 2);
3345:   if (cmapping) PetscAssertPointer(cmapping, 3);
3346:   PetscUseMethod(A, "MatISGetLocalToGlobalMapping_C", (Mat, ISLocalToGlobalMapping *, ISLocalToGlobalMapping *), (A, rmapping, cmapping));
3347:   PetscFunctionReturn(PETSC_SUCCESS);
3348: }

3350: static PetscErrorCode MatISGetLocalToGlobalMapping_IS(Mat A, ISLocalToGlobalMapping *r, ISLocalToGlobalMapping *c)
3351: {
3352:   Mat_IS *a = (Mat_IS *)A->data;

3354:   PetscFunctionBegin;
3355:   if (r) *r = a->rmapping;
3356:   if (c) *c = a->cmapping;
3357:   PetscFunctionReturn(PETSC_SUCCESS);
3358: }

3360: /*MC
3361:    MATIS - MATIS = "is" - A matrix type to be used for using the non-overlapping domain decomposition methods (e.g. `PCBDDC` or `KSPFETIDP`).
3362:    This stores the matrices in globally unassembled form. Each processor assembles only its local Neumann problem and the parallel matrix vector
3363:    product is handled "implicitly".

3365:    Options Database Keys:
3366: + -mat_type is - sets the matrix type to `MATIS`.
3367: . -matis_fixempty - Fixes local matrices in case of empty local rows/columns.
3368: - -matis_storel2l - stores the local-to-local operators generated by the Galerkin process of `MatPtAP()`.

3370:   Level: advanced

3372:    Notes:
3373:    Options prefix for the inner matrix are given by `-is_mat_xxx`

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

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

3380: .seealso: [](ch_matrices), `Mat`, `MATIS`, `Mat`, `MatISGetLocalMat()`, `MatSetLocalToGlobalMapping()`, `MatISSetPreallocation()`, `MatCreateIS()`, `PCBDDC`, `KSPFETIDP`
3381: M*/
3382: PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A)
3383: {
3384:   Mat_IS *a;

3386:   PetscFunctionBegin;
3387:   PetscCall(PetscNew(&a));
3388:   PetscCall(PetscStrallocpy(MATAIJ, &a->lmattype));
3389:   A->data = (void *)a;

3391:   /* matrix ops */
3392:   PetscCall(PetscMemzero(A->ops, sizeof(struct _MatOps)));
3393:   A->ops->mult                    = MatMult_IS;
3394:   A->ops->multadd                 = MatMultAdd_IS;
3395:   A->ops->multtranspose           = MatMultTranspose_IS;
3396:   A->ops->multtransposeadd        = MatMultTransposeAdd_IS;
3397:   A->ops->destroy                 = MatDestroy_IS;
3398:   A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS;
3399:   A->ops->setvalues               = MatSetValues_IS;
3400:   A->ops->setvaluesblocked        = MatSetValuesBlocked_IS;
3401:   A->ops->setvalueslocal          = MatSetValuesLocal_IS;
3402:   A->ops->setvaluesblockedlocal   = MatSetValuesBlockedLocal_IS;
3403:   A->ops->zerorows                = MatZeroRows_IS;
3404:   A->ops->zerorowscolumns         = MatZeroRowsColumns_IS;
3405:   A->ops->assemblybegin           = MatAssemblyBegin_IS;
3406:   A->ops->assemblyend             = MatAssemblyEnd_IS;
3407:   A->ops->view                    = MatView_IS;
3408:   A->ops->zeroentries             = MatZeroEntries_IS;
3409:   A->ops->scale                   = MatScale_IS;
3410:   A->ops->getdiagonal             = MatGetDiagonal_IS;
3411:   A->ops->setoption               = MatSetOption_IS;
3412:   A->ops->ishermitian             = MatIsHermitian_IS;
3413:   A->ops->issymmetric             = MatIsSymmetric_IS;
3414:   A->ops->isstructurallysymmetric = MatIsStructurallySymmetric_IS;
3415:   A->ops->duplicate               = MatDuplicate_IS;
3416:   A->ops->missingdiagonal         = MatMissingDiagonal_IS;
3417:   A->ops->copy                    = MatCopy_IS;
3418:   A->ops->getlocalsubmatrix       = MatGetLocalSubMatrix_IS;
3419:   A->ops->createsubmatrix         = MatCreateSubMatrix_IS;
3420:   A->ops->axpy                    = MatAXPY_IS;
3421:   A->ops->diagonalset             = MatDiagonalSet_IS;
3422:   A->ops->shift                   = MatShift_IS;
3423:   A->ops->transpose               = MatTranspose_IS;
3424:   A->ops->getinfo                 = MatGetInfo_IS;
3425:   A->ops->diagonalscale           = MatDiagonalScale_IS;
3426:   A->ops->setfromoptions          = MatSetFromOptions_IS;
3427:   A->ops->setup                   = MatSetUp_IS;
3428:   A->ops->hasoperation            = MatHasOperation_IS;
3429:   A->ops->getdiagonalblock        = MatGetDiagonalBlock_IS;
3430:   A->ops->createsubmatrices       = MatCreateSubMatrices_IS;
3431:   A->ops->increaseoverlap         = MatIncreaseOverlap_IS;

3433:   /* special MATIS functions */
3434:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMatType_C", MatISSetLocalMatType_IS));
3435:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalMat_C", MatISGetLocalMat_IS));
3436:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISRestoreLocalMat_C", MatISRestoreLocalMat_IS));
3437:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMat_C", MatISSetLocalMat_IS));
3438:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetMPIXAIJ_C", MatConvert_IS_XAIJ));
3439:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetPreallocation_C", MatISSetPreallocation_IS));
3440:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISStoreL2L_C", MatISStoreL2L_IS));
3441:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISFixLocalEmpty_C", MatISFixLocalEmpty_IS));
3442:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalToGlobalMapping_C", MatISGetLocalToGlobalMapping_IS));
3443:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpiaij_C", MatConvert_IS_XAIJ));
3444:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpibaij_C", MatConvert_IS_XAIJ));
3445:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpisbaij_C", MatConvert_IS_XAIJ));
3446:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqaij_C", MatConvert_IS_XAIJ));
3447:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqbaij_C", MatConvert_IS_XAIJ));
3448:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqsbaij_C", MatConvert_IS_XAIJ));
3449:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_aij_C", MatConvert_IS_XAIJ));
3450:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOOLocal_C", MatSetPreallocationCOOLocal_IS));
3451:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOO_C", MatSetPreallocationCOO_IS));
3452:   PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATIS));
3453:   PetscFunctionReturn(PETSC_SUCCESS);
3454: }