Actual source code: matis.c

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

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

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

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

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

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

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

 43:   PetscObjectQuery((PetscObject)C, "_MatIS_PtAP", (PetscObject *)&c);
 45:   PetscContainerGetPointer(c, (void **)&ptap);
 46:   ris[0] = ptap->ris0;
 47:   ris[1] = ptap->ris1;
 48:   cis[0] = ptap->cis0;
 49:   cis[1] = ptap->cis1;
 50:   n      = ptap->ris1 ? 2 : 1;
 51:   reuse  = ptap->lP ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX;
 52:   MatCreateSubMatrices(P, n, ris, cis, reuse, &ptap->lP);

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

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

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

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

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

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

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

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

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

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

177:   PetscNew(&ptap);
178:   PetscContainerCreate(PETSC_COMM_SELF, &c);
179:   PetscContainerSetPointer(c, ptap);
180:   PetscContainerSetUserDestroy(c, MatISContainerDestroyPtAP_Private);
181:   PetscObjectCompose((PetscObject)C, "_MatIS_PtAP", (PetscObject)c);
182:   PetscContainerDestroy(&c);
183:   ptap->fill = fill;

185:   MatISGetLocalToGlobalMapping(A, &rl2g, &cl2g);

187:   ISLocalToGlobalMappingGetBlockSize(cl2g, &ibs);
188:   ISLocalToGlobalMappingGetSize(cl2g, &N);
189:   ISLocalToGlobalMappingGetBlockIndices(cl2g, &garray);
190:   ISCreateBlock(comm, ibs, N / ibs, garray, PETSC_COPY_VALUES, &ptap->ris0);
191:   ISLocalToGlobalMappingRestoreBlockIndices(cl2g, &garray);

193:   MatCreateSubMatrix(P, ptap->ris0, NULL, MAT_INITIAL_MATRIX, &PT);
194:   MatGetNonzeroColumnsLocal_Private(PT, &ptap->cis0);
195:   ISLocalToGlobalMappingCreateIS(ptap->cis0, &Ccl2g);
196:   MatDestroy(&PT);

198:   Crl2g = NULL;
199:   if (rl2g != cl2g) { /* unsymmetric A mapping */
200:     PetscBool same, lsame = PETSC_FALSE;
201:     PetscInt  N1, ibs1;

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

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

235:   C->ops->ptapnumeric = MatPtAPNumeric_IS_XAIJ;
236:   return 0;
237: }

239: /* ----------------------------------------- */
240: static PetscErrorCode MatProductSymbolic_PtAP_IS_XAIJ(Mat C)
241: {
242:   Mat_Product *product = C->product;
243:   Mat          A = product->A, P = product->B;
244:   PetscReal    fill = product->fill;

246:   MatPtAPSymbolic_IS_XAIJ(A, P, fill, C);
247:   C->ops->productnumeric = MatProductNumeric_PtAP;
248:   return 0;
249: }

251: static PetscErrorCode MatProductSetFromOptions_IS_XAIJ_PtAP(Mat C)
252: {
253:   C->ops->productsymbolic = MatProductSymbolic_PtAP_IS_XAIJ;
254:   return 0;
255: }

257: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_IS_XAIJ(Mat C)
258: {
259:   Mat_Product *product = C->product;

261:   if (product->type == MATPRODUCT_PtAP) MatProductSetFromOptions_IS_XAIJ_PtAP(C);
262:   return 0;
263: }

265: /* ----------------------------------------- */
266: static PetscErrorCode MatISContainerDestroyFields_Private(void *ptr)
267: {
268:   MatISLocalFields lf = (MatISLocalFields)ptr;
269:   PetscInt         i;

271:   for (i = 0; i < lf->nr; i++) ISDestroy(&lf->rf[i]);
272:   for (i = 0; i < lf->nc; i++) ISDestroy(&lf->cf[i]);
273:   PetscFree2(lf->rf, lf->cf);
274:   PetscFree(lf);
275:   return 0;
276: }

278: static PetscErrorCode MatConvert_SeqXAIJ_IS(Mat A, MatType type, MatReuse reuse, Mat *newmat)
279: {
280:   Mat B, lB;

282:   if (reuse != MAT_REUSE_MATRIX) {
283:     ISLocalToGlobalMapping rl2g, cl2g;
284:     PetscInt               bs;
285:     IS                     is;

287:     MatGetBlockSize(A, &bs);
288:     ISCreateStride(PetscObjectComm((PetscObject)A), A->rmap->n / bs, 0, 1, &is);
289:     if (bs > 1) {
290:       IS       is2;
291:       PetscInt i, *aux;

293:       ISGetLocalSize(is, &i);
294:       ISGetIndices(is, (const PetscInt **)&aux);
295:       ISCreateBlock(PetscObjectComm((PetscObject)A), bs, i, aux, PETSC_COPY_VALUES, &is2);
296:       ISRestoreIndices(is, (const PetscInt **)&aux);
297:       ISDestroy(&is);
298:       is = is2;
299:     }
300:     ISSetIdentity(is);
301:     ISLocalToGlobalMappingCreateIS(is, &rl2g);
302:     ISDestroy(&is);
303:     ISCreateStride(PetscObjectComm((PetscObject)A), A->cmap->n / bs, 0, 1, &is);
304:     if (bs > 1) {
305:       IS       is2;
306:       PetscInt i, *aux;

308:       ISGetLocalSize(is, &i);
309:       ISGetIndices(is, (const PetscInt **)&aux);
310:       ISCreateBlock(PetscObjectComm((PetscObject)A), bs, i, aux, PETSC_COPY_VALUES, &is2);
311:       ISRestoreIndices(is, (const PetscInt **)&aux);
312:       ISDestroy(&is);
313:       is = is2;
314:     }
315:     ISSetIdentity(is);
316:     ISLocalToGlobalMappingCreateIS(is, &cl2g);
317:     ISDestroy(&is);
318:     MatCreateIS(PetscObjectComm((PetscObject)A), bs, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N, rl2g, cl2g, &B);
319:     ISLocalToGlobalMappingDestroy(&rl2g);
320:     ISLocalToGlobalMappingDestroy(&cl2g);
321:     MatDuplicate(A, MAT_COPY_VALUES, &lB);
322:     if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
323:   } else {
324:     B = *newmat;
325:     PetscObjectReference((PetscObject)A);
326:     lB = A;
327:   }
328:   MatISSetLocalMat(B, lB);
329:   MatDestroy(&lB);
330:   MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
331:   MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
332:   if (reuse == MAT_INPLACE_MATRIX) MatHeaderReplace(A, &B);
333:   return 0;
334: }

336: static PetscErrorCode MatISScaleDisassembling_Private(Mat A)
337: {
338:   Mat_IS         *matis = (Mat_IS *)(A->data);
339:   PetscScalar    *aa;
340:   const PetscInt *ii, *jj;
341:   PetscInt        i, n, m;
342:   PetscInt       *ecount, **eneighs;
343:   PetscBool       flg;

345:   MatGetRowIJ(matis->A, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg);
347:   ISLocalToGlobalMappingGetNodeInfo(matis->rmapping, &n, &ecount, &eneighs);
349:   MatSeqAIJGetArray(matis->A, &aa);
350:   for (i = 0; i < n; i++) {
351:     if (ecount[i] > 1) {
352:       PetscInt j;

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

358:         for (p = 0; p < ecount[i]; p++) {
359:           for (p2 = 0; p2 < ecount[i2]; p2++) {
360:             if (eneighs[i][p] == eneighs[i2][p2]) {
361:               scal += 1.0;
362:               break;
363:             }
364:           }
365:         }
366:         if (scal) aa[j] /= scal;
367:       }
368:     }
369:   }
370:   ISLocalToGlobalMappingRestoreNodeInfo(matis->rmapping, &n, &ecount, &eneighs);
371:   MatSeqAIJRestoreArray(matis->A, &aa);
372:   MatRestoreRowIJ(matis->A, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg);
374:   return 0;
375: }

377: typedef enum {
378:   MAT_IS_DISASSEMBLE_L2G_NATURAL,
379:   MAT_IS_DISASSEMBLE_L2G_MAT,
380:   MAT_IS_DISASSEMBLE_L2G_ND
381: } MatISDisassemblel2gType;

383: static PetscErrorCode MatMPIXAIJComputeLocalToGlobalMapping_Private(Mat A, ISLocalToGlobalMapping *l2g)
384: {
385:   Mat                     Ad, Ao;
386:   IS                      is, ndmap, ndsub;
387:   MPI_Comm                comm;
388:   const PetscInt         *garray, *ndmapi;
389:   PetscInt                bs, i, cnt, nl, *ncount, *ndmapc;
390:   PetscBool               ismpiaij, ismpibaij;
391:   const char *const       MatISDisassemblel2gTypes[] = {"NATURAL", "MAT", "ND", "MatISDisassemblel2gType", "MAT_IS_DISASSEMBLE_L2G_", NULL};
392:   MatISDisassemblel2gType mode                       = MAT_IS_DISASSEMBLE_L2G_NATURAL;
393:   MatPartitioning         part;
394:   PetscSF                 sf;
395:   PetscObject             dm;

397:   PetscOptionsBegin(PetscObjectComm((PetscObject)A), ((PetscObject)A)->prefix, "MatIS l2g disassembling options", "Mat");
398:   PetscOptionsEnum("-mat_is_disassemble_l2g_type", "Type of local-to-global mapping to be used for disassembling", "MatISDisassemblel2gType", MatISDisassemblel2gTypes, (PetscEnum)mode, (PetscEnum *)&mode, NULL);
399:   PetscOptionsEnd();
400:   if (mode == MAT_IS_DISASSEMBLE_L2G_MAT) {
401:     MatGetLocalToGlobalMapping(A, l2g, NULL);
402:     return 0;
403:   }
404:   PetscObjectGetComm((PetscObject)A, &comm);
405:   PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij);
406:   PetscObjectBaseTypeCompare((PetscObject)A, MATMPIBAIJ, &ismpibaij);
407:   MatGetBlockSize(A, &bs);
408:   switch (mode) {
409:   case MAT_IS_DISASSEMBLE_L2G_ND:
410:     MatPartitioningCreate(comm, &part);
411:     MatPartitioningSetAdjacency(part, A);
412:     PetscObjectSetOptionsPrefix((PetscObject)part, ((PetscObject)A)->prefix);
413:     MatPartitioningSetFromOptions(part);
414:     MatPartitioningApplyND(part, &ndmap);
415:     MatPartitioningDestroy(&part);
416:     ISBuildTwoSided(ndmap, NULL, &ndsub);
417:     MatMPIAIJSetUseScalableIncreaseOverlap(A, PETSC_TRUE);
418:     MatIncreaseOverlap(A, 1, &ndsub, 1);
419:     ISLocalToGlobalMappingCreateIS(ndsub, l2g);

421:     /* it may happen that a separator node is not properly shared */
422:     ISLocalToGlobalMappingGetNodeInfo(*l2g, &nl, &ncount, NULL);
423:     PetscSFCreate(comm, &sf);
424:     ISLocalToGlobalMappingGetIndices(*l2g, &garray);
425:     PetscSFSetGraphLayout(sf, A->rmap, nl, NULL, PETSC_OWN_POINTER, garray);
426:     ISLocalToGlobalMappingRestoreIndices(*l2g, &garray);
427:     PetscCalloc1(A->rmap->n, &ndmapc);
428:     PetscSFReduceBegin(sf, MPIU_INT, ncount, ndmapc, MPI_REPLACE);
429:     PetscSFReduceEnd(sf, MPIU_INT, ncount, ndmapc, MPI_REPLACE);
430:     ISLocalToGlobalMappingRestoreNodeInfo(*l2g, NULL, &ncount, NULL);
431:     ISGetIndices(ndmap, &ndmapi);
432:     for (i = 0, cnt = 0; i < A->rmap->n; i++)
433:       if (ndmapi[i] < 0 && ndmapc[i] < 2) cnt++;

435:     MPIU_Allreduce(&cnt, &i, 1, MPIU_INT, MPI_MAX, comm);
436:     if (i) { /* we detected isolated separator nodes */
437:       Mat                    A2, A3;
438:       IS                    *workis, is2;
439:       PetscScalar           *vals;
440:       PetscInt               gcnt = i, *dnz, *onz, j, *lndmapi;
441:       ISLocalToGlobalMapping ll2g;
442:       PetscBool              flg;
443:       const PetscInt        *ii, *jj;

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

449:       PetscMalloc1(nl, &lndmapi);
450:       PetscSFBcastBegin(sf, MPIU_INT, dnz, lndmapi, MPI_REPLACE);

452:       /* compute adjacency of isolated separators node */
453:       PetscMalloc1(gcnt, &workis);
454:       for (i = 0, cnt = 0; i < A->rmap->n; i++) {
455:         if (ndmapi[i] < 0 && ndmapc[i] < 2) ISCreateStride(comm, 1, i + A->rmap->rstart, 1, &workis[cnt++]);
456:       }
457:       for (i = cnt; i < gcnt; i++) ISCreateStride(comm, 0, 0, 1, &workis[i]);
458:       for (i = 0; i < gcnt; i++) {
459:         PetscObjectSetName((PetscObject)workis[i], "ISOLATED");
460:         ISViewFromOptions(workis[i], NULL, "-view_isolated_separators");
461:       }

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

466:       /* end communicate global id of separators */
467:       PetscSFBcastEnd(sf, MPIU_INT, dnz, lndmapi, MPI_REPLACE);

469:       /* communicate new layers : create a matrix and transpose it */
470:       PetscArrayzero(dnz, A->rmap->n);
471:       PetscArrayzero(onz, A->rmap->n);
472:       for (i = 0, j = 0; i < A->rmap->n; i++) {
473:         if (ndmapi[i] < 0 && ndmapc[i] < 2) {
474:           const PetscInt *idxs;
475:           PetscInt        s;

477:           ISGetLocalSize(workis[j], &s);
478:           ISGetIndices(workis[j], &idxs);
479:           MatPreallocateSet(i + A->rmap->rstart, s, idxs, dnz, onz);
480:           j++;
481:         }
482:       }

485:       for (i = 0; i < gcnt; i++) {
486:         PetscObjectSetName((PetscObject)workis[i], "EXTENDED");
487:         ISViewFromOptions(workis[i], NULL, "-view_isolated_separators");
488:       }

490:       for (i = 0, j = 0; i < A->rmap->n; i++) j = PetscMax(j, dnz[i] + onz[i]);
491:       PetscMalloc1(j, &vals);
492:       for (i = 0; i < j; i++) vals[i] = 1.0;

494:       MatCreate(comm, &A2);
495:       MatSetType(A2, MATMPIAIJ);
496:       MatSetSizes(A2, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N);
497:       MatMPIAIJSetPreallocation(A2, 0, dnz, 0, onz);
498:       MatSetOption(A2, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE);
499:       for (i = 0, j = 0; i < A2->rmap->n; i++) {
500:         PetscInt        row = i + A2->rmap->rstart, s = dnz[i] + onz[i];
501:         const PetscInt *idxs;

503:         if (s) {
504:           ISGetIndices(workis[j], &idxs);
505:           MatSetValues(A2, 1, &row, s, idxs, vals, INSERT_VALUES);
506:           ISRestoreIndices(workis[j], &idxs);
507:           j++;
508:         }
509:       }
511:       PetscFree(vals);
512:       MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY);
513:       MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY);
514:       MatTranspose(A2, MAT_INPLACE_MATRIX, &A2);

516:       /* extract submatrix corresponding to the coupling "owned separators" x "isolated separators" */
517:       for (i = 0, j = 0; i < nl; i++)
518:         if (lndmapi[i] >= 0) lndmapi[j++] = lndmapi[i];
519:       ISCreateGeneral(comm, j, lndmapi, PETSC_USE_POINTER, &is);
520:       MatMPIAIJGetLocalMatCondensed(A2, MAT_INITIAL_MATRIX, &is, NULL, &A3);
521:       ISDestroy(&is);
522:       MatDestroy(&A2);

524:       /* extend local to global map to include connected isolated separators */
525:       PetscObjectQuery((PetscObject)A3, "_petsc_GetLocalMatCondensed_iscol", (PetscObject *)&is);
527:       ISLocalToGlobalMappingCreateIS(is, &ll2g);
528:       MatGetRowIJ(A3, 0, PETSC_FALSE, PETSC_FALSE, &i, &ii, &jj, &flg);
530:       ISCreateGeneral(PETSC_COMM_SELF, ii[i], jj, PETSC_COPY_VALUES, &is);
531:       MatRestoreRowIJ(A3, 0, PETSC_FALSE, PETSC_FALSE, &i, &ii, &jj, &flg);
533:       ISLocalToGlobalMappingApplyIS(ll2g, is, &is2);
534:       ISDestroy(&is);
535:       ISLocalToGlobalMappingDestroy(&ll2g);

537:       /* add new nodes to the local-to-global map */
538:       ISLocalToGlobalMappingDestroy(l2g);
539:       ISExpand(ndsub, is2, &is);
540:       ISDestroy(&is2);
541:       ISLocalToGlobalMappingCreateIS(is, l2g);
542:       ISDestroy(&is);

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

573:       MatGetLocalSize(Ad, NULL, &dc);
574:       MatGetLocalSize(Ao, NULL, &oc);
576:       MatGetOwnershipRangeColumn(A, &stc, NULL);
577:       PetscMalloc1((dc + oc) / bs, &aux);
578:       for (i = 0; i < dc / bs; i++) aux[i] = i + stc / bs;
579:       for (i = 0; i < oc / bs; i++) aux[i + dc / bs] = (ismpiaij ? garray[i * bs] / bs : garray[i]);
580:       ISCreateBlock(comm, bs, (dc + oc) / bs, aux, PETSC_OWN_POINTER, &is);
581:     } else {
582:       ISCreateBlock(comm, 1, 0, NULL, PETSC_OWN_POINTER, &is);
583:     }
584:     ISLocalToGlobalMappingCreateIS(is, l2g);
585:     ISDestroy(&is);
586:     break;
587:   default:
588:     SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unsupported l2g disassembling type %d", mode);
589:   }
590:   return 0;
591: }

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

610:   PetscObjectGetComm((PetscObject)A, &comm);
611:   MPI_Comm_size(comm, &size);
612:   if (size == 1) {
613:     MatConvert_SeqXAIJ_IS(A, type, reuse, newmat);
614:     return 0;
615:   }
616:   if (reuse != MAT_REUSE_MATRIX && A->cmap->N == A->rmap->N) {
617:     MatMPIXAIJComputeLocalToGlobalMapping_Private(A, &rl2g);
618:     MatCreate(comm, &B);
619:     MatSetType(B, MATIS);
620:     MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N);
621:     MatSetLocalToGlobalMapping(B, rl2g, rl2g);
622:     MatGetBlockSize(A, &bs);
623:     MatSetBlockSize(B, bs);
624:     ISLocalToGlobalMappingDestroy(&rl2g);
625:     if (reuse == MAT_INPLACE_MATRIX) was_inplace = PETSC_TRUE;
626:     reuse = MAT_REUSE_MATRIX;
627:   }
628:   if (reuse == MAT_REUSE_MATRIX) {
629:     Mat            *newlA, lA;
630:     IS              rows, cols;
631:     const PetscInt *ridx, *cidx;
632:     PetscInt        rbs, cbs, nr, nc;

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

686:   /* access relevant information from MPIAIJ */
687:   MatGetOwnershipRange(A, &str, NULL);
688:   MatGetOwnershipRangeColumn(A, &stc, NULL);
689:   MatGetLocalSize(A, &dr, &dc);
690:   MatGetLocalSize(Ao, NULL, &oc);
691:   MatGetRowIJ(Ad, 0, PETSC_FALSE, PETSC_FALSE, &i, &di, &dj, &flg);
693:   MatGetRowIJ(Ao, 0, PETSC_FALSE, PETSC_FALSE, &i, &oi, &oj, &flg);
695:   nnz = di[dr] + oi[dr];
696:   /* store original pointers to be restored later */
697:   odi = di;
698:   odj = dj;
699:   ooi = oi;
700:   ooj = oj;

702:   /* generate l2g maps for rows and cols */
703:   ISCreateStride(comm, dr / bs, str / bs, 1, &is);
704:   if (bs > 1) {
705:     IS is2;

707:     ISGetLocalSize(is, &i);
708:     ISGetIndices(is, (const PetscInt **)&aux);
709:     ISCreateBlock(comm, bs, i, aux, PETSC_COPY_VALUES, &is2);
710:     ISRestoreIndices(is, (const PetscInt **)&aux);
711:     ISDestroy(&is);
712:     is = is2;
713:   }
714:   ISLocalToGlobalMappingCreateIS(is, &rl2g);
715:   ISDestroy(&is);
716:   if (dr) {
717:     PetscMalloc1((dc + oc) / bs, &aux);
718:     for (i = 0; i < dc / bs; i++) aux[i] = i + stc / bs;
719:     for (i = 0; i < oc / bs; i++) aux[i + dc / bs] = garray[i];
720:     ISCreateBlock(comm, bs, (dc + oc) / bs, aux, PETSC_OWN_POINTER, &is);
721:     lc = dc + oc;
722:   } else {
723:     ISCreateBlock(comm, bs, 0, NULL, PETSC_OWN_POINTER, &is);
724:     lc = 0;
725:   }
726:   ISLocalToGlobalMappingCreateIS(is, &cl2g);
727:   ISDestroy(&is);

729:   /* create MATIS object */
730:   MatCreate(comm, &B);
731:   MatSetSizes(B, dr, dc, PETSC_DECIDE, PETSC_DECIDE);
732:   MatSetType(B, MATIS);
733:   MatSetBlockSize(B, bs);
734:   MatSetLocalToGlobalMapping(B, rl2g, cl2g);
735:   ISLocalToGlobalMappingDestroy(&rl2g);
736:   ISLocalToGlobalMappingDestroy(&cl2g);

738:   /* merge local matrices */
739:   PetscMalloc1(nnz + dr + 1, &aux);
740:   PetscMalloc1(nnz, &data);
741:   ii  = aux;
742:   jj  = aux + dr + 1;
743:   aa  = data;
744:   *ii = *(di++) + *(oi++);
745:   for (jd = 0, jo = 0, cum = 0; *ii < nnz; cum++) {
746:     for (; jd < *di; jd++) {
747:       *jj++ = *dj++;
748:       *aa++ = *dd++;
749:     }
750:     for (; jo < *oi; jo++) {
751:       *jj++ = *oj++ + dc;
752:       *aa++ = *od++;
753:     }
754:     *(++ii) = *(di++) + *(oi++);
755:   }
756:   for (; cum < dr; cum++) *(++ii) = nnz;

758:   MatRestoreRowIJ(Ad, 0, PETSC_FALSE, PETSC_FALSE, &i, &odi, &odj, &flg);
760:   MatRestoreRowIJ(Ao, 0, PETSC_FALSE, PETSC_FALSE, &i, &ooi, &ooj, &flg);
762:   MatSeqAIJRestoreArray(Ad, &dd);
763:   MatSeqAIJRestoreArray(Ao, &od);

765:   ii = aux;
766:   jj = aux + dr + 1;
767:   aa = data;
768:   MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, dr, lc, ii, jj, aa, &lA);

770:   /* create containers to destroy the data */
771:   ptrs[0] = aux;
772:   ptrs[1] = data;
773:   for (i = 0; i < 2; i++) {
774:     PetscContainer c;

776:     PetscContainerCreate(PETSC_COMM_SELF, &c);
777:     PetscContainerSetPointer(c, ptrs[i]);
778:     PetscContainerSetUserDestroy(c, PetscContainerUserDestroyDefault);
779:     PetscObjectCompose((PetscObject)lA, names[i], (PetscObject)c);
780:     PetscContainerDestroy(&c);
781:   }
782:   if (ismpibaij) { /* destroy converted local matrices */
783:     MatDestroy(&Ad);
784:     MatDestroy(&Ao);
785:   }

787:   /* finalize matrix */
788:   MatISSetLocalMat(B, lA);
789:   MatDestroy(&lA);
790:   MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
791:   MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
792:   if (reuse == MAT_INPLACE_MATRIX) {
793:     MatHeaderReplace(A, &B);
794:   } else *newmat = B;
795:   return 0;
796: }

798: PETSC_INTERN PetscErrorCode MatConvert_Nest_IS(Mat A, MatType type, MatReuse reuse, Mat *newmat)
799: {
800:   Mat                  **nest, *snest, **rnest, lA, B;
801:   IS                    *iscol, *isrow, *islrow, *islcol;
802:   ISLocalToGlobalMapping rl2g, cl2g;
803:   MPI_Comm               comm;
804:   PetscInt              *lr, *lc, *l2gidxs;
805:   PetscInt               i, j, nr, nc, rbs, cbs;
806:   PetscBool              convert, lreuse, *istrans;

808:   MatNestGetSubMats(A, &nr, &nc, &nest);
809:   lreuse = PETSC_FALSE;
810:   rnest  = NULL;
811:   if (reuse == MAT_REUSE_MATRIX) {
812:     PetscBool ismatis, isnest;

814:     PetscObjectTypeCompare((PetscObject)*newmat, MATIS, &ismatis);
816:     MatISGetLocalMat(*newmat, &lA);
817:     PetscObjectTypeCompare((PetscObject)lA, MATNEST, &isnest);
818:     if (isnest) {
819:       MatNestGetSubMats(lA, &i, &j, &rnest);
820:       lreuse = (PetscBool)(i == nr && j == nc);
821:       if (!lreuse) rnest = NULL;
822:     }
823:   }
824:   PetscObjectGetComm((PetscObject)A, &comm);
825:   PetscCalloc2(nr, &lr, nc, &lc);
826:   PetscCalloc6(nr, &isrow, nc, &iscol, nr, &islrow, nc, &islcol, nr * nc, &snest, nr * nc, &istrans);
827:   MatNestGetISs(A, isrow, iscol);
828:   for (i = 0; i < nr; i++) {
829:     for (j = 0; j < nc; j++) {
830:       PetscBool ismatis;
831:       PetscInt  l1, l2, lb1, lb2, ij = i * nc + j;

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

836:       /* Nested matrices should be of type MATIS */
837:       PetscObjectTypeCompare((PetscObject)nest[i][j], MATTRANSPOSEVIRTUAL, &istrans[ij]);
838:       if (istrans[ij]) {
839:         Mat T, lT;
840:         MatTransposeGetMat(nest[i][j], &T);
841:         PetscObjectTypeCompare((PetscObject)T, MATIS, &ismatis);
843:         MatISGetLocalMat(T, &lT);
844:         MatCreateTranspose(lT, &snest[ij]);
845:       } else {
846:         PetscObjectTypeCompare((PetscObject)nest[i][j], MATIS, &ismatis);
848:         MatISGetLocalMat(nest[i][j], &snest[ij]);
849:       }

851:       /* Check compatibility of local sizes */
852:       MatGetSize(snest[ij], &l1, &l2);
853:       MatGetBlockSizes(snest[ij], &lb1, &lb2);
854:       if (!l1 || !l2) continue;
857:       lr[i] = l1;
858:       lc[j] = l2;

860:       /* check compatibility for local matrix reusage */
861:       if (rnest && !rnest[i][j] != !snest[ij]) lreuse = PETSC_FALSE;
862:     }
863:   }

865:   if (PetscDefined(USE_DEBUG)) {
866:     /* Check compatibility of l2g maps for rows */
867:     for (i = 0; i < nr; i++) {
868:       rl2g = NULL;
869:       for (j = 0; j < nc; j++) {
870:         PetscInt n1, n2;

872:         if (!nest[i][j]) continue;
873:         if (istrans[i * nc + j]) {
874:           Mat T;

876:           MatTransposeGetMat(nest[i][j], &T);
877:           MatISGetLocalToGlobalMapping(T, NULL, &cl2g);
878:         } else {
879:           MatISGetLocalToGlobalMapping(nest[i][j], &cl2g, NULL);
880:         }
881:         ISLocalToGlobalMappingGetSize(cl2g, &n1);
882:         if (!n1) continue;
883:         if (!rl2g) {
884:           rl2g = cl2g;
885:         } else {
886:           const PetscInt *idxs1, *idxs2;
887:           PetscBool       same;

889:           ISLocalToGlobalMappingGetSize(rl2g, &n2);
891:           ISLocalToGlobalMappingGetIndices(cl2g, &idxs1);
892:           ISLocalToGlobalMappingGetIndices(rl2g, &idxs2);
893:           PetscArraycmp(idxs1, idxs2, n1, &same);
894:           ISLocalToGlobalMappingRestoreIndices(cl2g, &idxs1);
895:           ISLocalToGlobalMappingRestoreIndices(rl2g, &idxs2);
897:         }
898:       }
899:     }
900:     /* Check compatibility of l2g maps for columns */
901:     for (i = 0; i < nc; i++) {
902:       rl2g = NULL;
903:       for (j = 0; j < nr; j++) {
904:         PetscInt n1, n2;

906:         if (!nest[j][i]) continue;
907:         if (istrans[j * nc + i]) {
908:           Mat T;

910:           MatTransposeGetMat(nest[j][i], &T);
911:           MatISGetLocalToGlobalMapping(T, &cl2g, NULL);
912:         } else {
913:           MatISGetLocalToGlobalMapping(nest[j][i], NULL, &cl2g);
914:         }
915:         ISLocalToGlobalMappingGetSize(cl2g, &n1);
916:         if (!n1) continue;
917:         if (!rl2g) {
918:           rl2g = cl2g;
919:         } else {
920:           const PetscInt *idxs1, *idxs2;
921:           PetscBool       same;

923:           ISLocalToGlobalMappingGetSize(rl2g, &n2);
925:           ISLocalToGlobalMappingGetIndices(cl2g, &idxs1);
926:           ISLocalToGlobalMappingGetIndices(rl2g, &idxs2);
927:           PetscArraycmp(idxs1, idxs2, n1, &same);
928:           ISLocalToGlobalMappingRestoreIndices(cl2g, &idxs1);
929:           ISLocalToGlobalMappingRestoreIndices(rl2g, &idxs2);
931:         }
932:       }
933:     }
934:   }

936:   B = NULL;
937:   if (reuse != MAT_REUSE_MATRIX) {
938:     PetscInt stl;

940:     /* Create l2g map for the rows of the new matrix and index sets for the local MATNEST */
941:     for (i = 0, stl = 0; i < nr; i++) stl += lr[i];
942:     PetscMalloc1(stl, &l2gidxs);
943:     for (i = 0, stl = 0; i < nr; i++) {
944:       Mat             usedmat;
945:       Mat_IS         *matis;
946:       const PetscInt *idxs;

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

951:       /* l2gmap */
952:       j       = 0;
953:       usedmat = nest[i][j];
954:       while (!usedmat && j < nc - 1) usedmat = nest[i][++j];

957:       if (istrans[i * nc + j]) {
958:         Mat T;
959:         MatTransposeGetMat(usedmat, &T);
960:         usedmat = T;
961:       }
962:       matis = (Mat_IS *)(usedmat->data);
963:       ISGetIndices(isrow[i], &idxs);
964:       if (istrans[i * nc + j]) {
965:         PetscSFBcastBegin(matis->csf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE);
966:         PetscSFBcastEnd(matis->csf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE);
967:       } else {
968:         PetscSFBcastBegin(matis->sf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE);
969:         PetscSFBcastEnd(matis->sf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE);
970:       }
971:       ISRestoreIndices(isrow[i], &idxs);
972:       stl += lr[i];
973:     }
974:     ISLocalToGlobalMappingCreate(comm, 1, stl, l2gidxs, PETSC_OWN_POINTER, &rl2g);

976:     /* Create l2g map for columns of the new matrix and index sets for the local MATNEST */
977:     for (i = 0, stl = 0; i < nc; i++) stl += lc[i];
978:     PetscMalloc1(stl, &l2gidxs);
979:     for (i = 0, stl = 0; i < nc; i++) {
980:       Mat             usedmat;
981:       Mat_IS         *matis;
982:       const PetscInt *idxs;

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

987:       /* l2gmap */
988:       j       = 0;
989:       usedmat = nest[j][i];
990:       while (!usedmat && j < nr - 1) usedmat = nest[++j][i];
992:       if (istrans[j * nc + i]) {
993:         Mat T;
994:         MatTransposeGetMat(usedmat, &T);
995:         usedmat = T;
996:       }
997:       matis = (Mat_IS *)(usedmat->data);
998:       ISGetIndices(iscol[i], &idxs);
999:       if (istrans[j * nc + i]) {
1000:         PetscSFBcastBegin(matis->sf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE);
1001:         PetscSFBcastEnd(matis->sf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE);
1002:       } else {
1003:         PetscSFBcastBegin(matis->csf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE);
1004:         PetscSFBcastEnd(matis->csf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE);
1005:       }
1006:       ISRestoreIndices(iscol[i], &idxs);
1007:       stl += lc[i];
1008:     }
1009:     ISLocalToGlobalMappingCreate(comm, 1, stl, l2gidxs, PETSC_OWN_POINTER, &cl2g);

1011:     /* Create MATIS */
1012:     MatCreate(comm, &B);
1013:     MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N);
1014:     MatGetBlockSizes(A, &rbs, &cbs);
1015:     MatSetBlockSizes(B, rbs, cbs);
1016:     MatSetType(B, MATIS);
1017:     MatISSetLocalMatType(B, MATNEST);
1018:     { /* hack : avoid setup of scatters */
1019:       Mat_IS *matis     = (Mat_IS *)(B->data);
1020:       matis->islocalref = PETSC_TRUE;
1021:     }
1022:     MatSetLocalToGlobalMapping(B, rl2g, cl2g);
1023:     ISLocalToGlobalMappingDestroy(&rl2g);
1024:     ISLocalToGlobalMappingDestroy(&cl2g);
1025:     MatCreateNest(PETSC_COMM_SELF, nr, islrow, nc, islcol, snest, &lA);
1026:     MatNestSetVecType(lA, VECNEST);
1027:     for (i = 0; i < nr * nc; i++) {
1028:       if (istrans[i]) MatDestroy(&snest[i]);
1029:     }
1030:     MatISSetLocalMat(B, lA);
1031:     MatDestroy(&lA);
1032:     { /* hack : setup of scatters done here */
1033:       Mat_IS *matis = (Mat_IS *)(B->data);

1035:       matis->islocalref = PETSC_FALSE;
1036:       MatISSetUpScatters_Private(B);
1037:     }
1038:     MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
1039:     MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
1040:     if (reuse == MAT_INPLACE_MATRIX) {
1041:       MatHeaderReplace(A, &B);
1042:     } else {
1043:       *newmat = B;
1044:     }
1045:   } else {
1046:     if (lreuse) {
1047:       MatISGetLocalMat(*newmat, &lA);
1048:       for (i = 0; i < nr; i++) {
1049:         for (j = 0; j < nc; j++) {
1050:           if (snest[i * nc + j]) {
1051:             MatNestSetSubMat(lA, i, j, snest[i * nc + j]);
1052:             if (istrans[i * nc + j]) MatDestroy(&snest[i * nc + j]);
1053:           }
1054:         }
1055:       }
1056:     } else {
1057:       PetscInt stl;
1058:       for (i = 0, stl = 0; i < nr; i++) {
1059:         ISCreateStride(PETSC_COMM_SELF, lr[i], stl, 1, &islrow[i]);
1060:         stl += lr[i];
1061:       }
1062:       for (i = 0, stl = 0; i < nc; i++) {
1063:         ISCreateStride(PETSC_COMM_SELF, lc[i], stl, 1, &islcol[i]);
1064:         stl += lc[i];
1065:       }
1066:       MatCreateNest(PETSC_COMM_SELF, nr, islrow, nc, islcol, snest, &lA);
1067:       for (i = 0; i < nr * nc; i++) {
1068:         if (istrans[i]) MatDestroy(&snest[i]);
1069:       }
1070:       MatISSetLocalMat(*newmat, lA);
1071:       MatDestroy(&lA);
1072:     }
1073:     MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY);
1074:     MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY);
1075:   }

1077:   /* Create local matrix in MATNEST format */
1078:   convert = PETSC_FALSE;
1079:   PetscOptionsGetBool(NULL, ((PetscObject)A)->prefix, "-matis_convert_local_nest", &convert, NULL);
1080:   if (convert) {
1081:     Mat              M;
1082:     MatISLocalFields lf;
1083:     PetscContainer   c;

1085:     MatISGetLocalMat(*newmat, &lA);
1086:     MatConvert(lA, MATAIJ, MAT_INITIAL_MATRIX, &M);
1087:     MatISSetLocalMat(*newmat, M);
1088:     MatDestroy(&M);

1090:     /* attach local fields to the matrix */
1091:     PetscNew(&lf);
1092:     PetscMalloc2(nr, &lf->rf, nc, &lf->cf);
1093:     for (i = 0; i < nr; i++) {
1094:       PetscInt n, st;

1096:       ISGetLocalSize(islrow[i], &n);
1097:       ISStrideGetInfo(islrow[i], &st, NULL);
1098:       ISCreateStride(comm, n, st, 1, &lf->rf[i]);
1099:     }
1100:     for (i = 0; i < nc; i++) {
1101:       PetscInt n, st;

1103:       ISGetLocalSize(islcol[i], &n);
1104:       ISStrideGetInfo(islcol[i], &st, NULL);
1105:       ISCreateStride(comm, n, st, 1, &lf->cf[i]);
1106:     }
1107:     lf->nr = nr;
1108:     lf->nc = nc;
1109:     PetscContainerCreate(PetscObjectComm((PetscObject)(*newmat)), &c);
1110:     PetscContainerSetPointer(c, lf);
1111:     PetscContainerSetUserDestroy(c, MatISContainerDestroyFields_Private);
1112:     PetscObjectCompose((PetscObject)(*newmat), "_convert_nest_lfields", (PetscObject)c);
1113:     PetscContainerDestroy(&c);
1114:   }

1116:   /* Free workspace */
1117:   for (i = 0; i < nr; i++) ISDestroy(&islrow[i]);
1118:   for (i = 0; i < nc; i++) ISDestroy(&islcol[i]);
1119:   PetscFree6(isrow, iscol, islrow, islcol, snest, istrans);
1120:   PetscFree2(lr, lc);
1121:   return 0;
1122: }

1124: static PetscErrorCode MatDiagonalScale_IS(Mat A, Vec l, Vec r)
1125: {
1126:   Mat_IS            *matis = (Mat_IS *)A->data;
1127:   Vec                ll, rr;
1128:   const PetscScalar *Y, *X;
1129:   PetscScalar       *x, *y;

1131:   if (l) {
1132:     ll = matis->y;
1133:     VecGetArrayRead(l, &Y);
1134:     VecGetArray(ll, &y);
1135:     PetscSFBcastBegin(matis->sf, MPIU_SCALAR, Y, y, MPI_REPLACE);
1136:   } else {
1137:     ll = NULL;
1138:   }
1139:   if (r) {
1140:     rr = matis->x;
1141:     VecGetArrayRead(r, &X);
1142:     VecGetArray(rr, &x);
1143:     PetscSFBcastBegin(matis->csf, MPIU_SCALAR, X, x, MPI_REPLACE);
1144:   } else {
1145:     rr = NULL;
1146:   }
1147:   if (ll) {
1148:     PetscSFBcastEnd(matis->sf, MPIU_SCALAR, Y, y, MPI_REPLACE);
1149:     VecRestoreArrayRead(l, &Y);
1150:     VecRestoreArray(ll, &y);
1151:   }
1152:   if (rr) {
1153:     PetscSFBcastEnd(matis->csf, MPIU_SCALAR, X, x, MPI_REPLACE);
1154:     VecRestoreArrayRead(r, &X);
1155:     VecRestoreArray(rr, &x);
1156:   }
1157:   MatDiagonalScale(matis->A, ll, rr);
1158:   return 0;
1159: }

1161: static PetscErrorCode MatGetInfo_IS(Mat A, MatInfoType flag, MatInfo *ginfo)
1162: {
1163:   Mat_IS        *matis = (Mat_IS *)A->data;
1164:   MatInfo        info;
1165:   PetscLogDouble isend[6], irecv[6];
1166:   PetscInt       bs;

1168:   MatGetBlockSize(A, &bs);
1169:   if (matis->A->ops->getinfo) {
1170:     MatGetInfo(matis->A, MAT_LOCAL, &info);
1171:     isend[0] = info.nz_used;
1172:     isend[1] = info.nz_allocated;
1173:     isend[2] = info.nz_unneeded;
1174:     isend[3] = info.memory;
1175:     isend[4] = info.mallocs;
1176:   } else {
1177:     isend[0] = 0.;
1178:     isend[1] = 0.;
1179:     isend[2] = 0.;
1180:     isend[3] = 0.;
1181:     isend[4] = 0.;
1182:   }
1183:   isend[5] = matis->A->num_ass;
1184:   if (flag == MAT_LOCAL) {
1185:     ginfo->nz_used      = isend[0];
1186:     ginfo->nz_allocated = isend[1];
1187:     ginfo->nz_unneeded  = isend[2];
1188:     ginfo->memory       = isend[3];
1189:     ginfo->mallocs      = isend[4];
1190:     ginfo->assemblies   = isend[5];
1191:   } else if (flag == MAT_GLOBAL_MAX) {
1192:     MPIU_Allreduce(isend, irecv, 6, MPIU_PETSCLOGDOUBLE, MPI_MAX, PetscObjectComm((PetscObject)A));

1194:     ginfo->nz_used      = irecv[0];
1195:     ginfo->nz_allocated = irecv[1];
1196:     ginfo->nz_unneeded  = irecv[2];
1197:     ginfo->memory       = irecv[3];
1198:     ginfo->mallocs      = irecv[4];
1199:     ginfo->assemblies   = irecv[5];
1200:   } else if (flag == MAT_GLOBAL_SUM) {
1201:     MPIU_Allreduce(isend, irecv, 5, MPIU_PETSCLOGDOUBLE, MPI_SUM, PetscObjectComm((PetscObject)A));

1203:     ginfo->nz_used      = irecv[0];
1204:     ginfo->nz_allocated = irecv[1];
1205:     ginfo->nz_unneeded  = irecv[2];
1206:     ginfo->memory       = irecv[3];
1207:     ginfo->mallocs      = irecv[4];
1208:     ginfo->assemblies   = A->num_ass;
1209:   }
1210:   ginfo->block_size        = bs;
1211:   ginfo->fill_ratio_given  = 0;
1212:   ginfo->fill_ratio_needed = 0;
1213:   ginfo->factor_mallocs    = 0;
1214:   return 0;
1215: }

1217: static PetscErrorCode MatTranspose_IS(Mat A, MatReuse reuse, Mat *B)
1218: {
1219:   Mat C, lC, lA;

1221:   if (reuse == MAT_REUSE_MATRIX) MatTransposeCheckNonzeroState_Private(A, *B);
1222:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
1223:     ISLocalToGlobalMapping rl2g, cl2g;
1224:     MatCreate(PetscObjectComm((PetscObject)A), &C);
1225:     MatSetSizes(C, A->cmap->n, A->rmap->n, A->cmap->N, A->rmap->N);
1226:     MatSetBlockSizes(C, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs));
1227:     MatSetType(C, MATIS);
1228:     MatGetLocalToGlobalMapping(A, &rl2g, &cl2g);
1229:     MatSetLocalToGlobalMapping(C, cl2g, rl2g);
1230:   } else C = *B;

1232:   /* perform local transposition */
1233:   MatISGetLocalMat(A, &lA);
1234:   MatTranspose(lA, MAT_INITIAL_MATRIX, &lC);
1235:   MatSetLocalToGlobalMapping(lC, lA->cmap->mapping, lA->rmap->mapping);
1236:   MatISSetLocalMat(C, lC);
1237:   MatDestroy(&lC);

1239:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
1240:     *B = C;
1241:   } else {
1242:     MatHeaderMerge(A, &C);
1243:   }
1244:   MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY);
1245:   MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY);
1246:   return 0;
1247: }

1249: static PetscErrorCode MatDiagonalSet_IS(Mat A, Vec D, InsertMode insmode)
1250: {
1251:   Mat_IS *is = (Mat_IS *)A->data;

1253:   if (D) { /* MatShift_IS pass D = NULL */
1254:     VecScatterBegin(is->rctx, D, is->y, INSERT_VALUES, SCATTER_FORWARD);
1255:     VecScatterEnd(is->rctx, D, is->y, INSERT_VALUES, SCATTER_FORWARD);
1256:   }
1257:   VecPointwiseDivide(is->y, is->y, is->counter);
1258:   MatDiagonalSet(is->A, is->y, insmode);
1259:   return 0;
1260: }

1262: static PetscErrorCode MatShift_IS(Mat A, PetscScalar a)
1263: {
1264:   Mat_IS *is = (Mat_IS *)A->data;

1266:   VecSet(is->y, a);
1267:   MatDiagonalSet_IS(A, NULL, ADD_VALUES);
1268:   return 0;
1269: }

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

1276:   ISLocalToGlobalMappingApply(A->rmap->mapping, m, rows, rows_l);
1277:   ISLocalToGlobalMappingApply(A->cmap->mapping, n, cols, cols_l);
1278:   MatSetValuesLocal_IS(A, m, rows_l, n, cols_l, values, addv);
1279:   return 0;
1280: }

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

1287:   ISLocalToGlobalMappingApplyBlock(A->rmap->mapping, m, rows, rows_l);
1288:   ISLocalToGlobalMappingApplyBlock(A->cmap->mapping, n, cols, cols_l);
1289:   MatSetValuesBlockedLocal_IS(A, m, rows_l, n, cols_l, values, addv);
1290:   return 0;
1291: }

1293: static PetscErrorCode MatCreateSubMatrix_IS(Mat mat, IS irow, IS icol, MatReuse scall, Mat *newmat)
1294: {
1295:   Mat             locmat, newlocmat;
1296:   Mat_IS         *newmatis;
1297:   const PetscInt *idxs;
1298:   PetscInt        i, m, n;

1300:   if (scall == MAT_REUSE_MATRIX) {
1301:     PetscBool ismatis;

1303:     PetscObjectTypeCompare((PetscObject)*newmat, MATIS, &ismatis);
1305:     newmatis = (Mat_IS *)(*newmat)->data;
1308:   }
1309:   /* irow and icol may not have duplicate entries */
1310:   if (PetscDefined(USE_DEBUG)) {
1311:     Vec                rtest, ltest;
1312:     const PetscScalar *array;

1314:     MatCreateVecs(mat, &ltest, &rtest);
1315:     ISGetLocalSize(irow, &n);
1316:     ISGetIndices(irow, &idxs);
1317:     for (i = 0; i < n; i++) VecSetValue(rtest, idxs[i], 1.0, ADD_VALUES);
1318:     VecAssemblyBegin(rtest);
1319:     VecAssemblyEnd(rtest);
1320:     VecGetLocalSize(rtest, &n);
1321:     VecGetOwnershipRange(rtest, &m, NULL);
1322:     VecGetArrayRead(rtest, &array);
1324:     VecRestoreArrayRead(rtest, &array);
1325:     ISRestoreIndices(irow, &idxs);
1326:     ISGetLocalSize(icol, &n);
1327:     ISGetIndices(icol, &idxs);
1328:     for (i = 0; i < n; i++) VecSetValue(ltest, idxs[i], 1.0, ADD_VALUES);
1329:     VecAssemblyBegin(ltest);
1330:     VecAssemblyEnd(ltest);
1331:     VecGetLocalSize(ltest, &n);
1332:     VecGetOwnershipRange(ltest, &m, NULL);
1333:     VecGetArrayRead(ltest, &array);
1335:     VecRestoreArrayRead(ltest, &array);
1336:     ISRestoreIndices(icol, &idxs);
1337:     VecDestroy(&rtest);
1338:     VecDestroy(&ltest);
1339:   }
1340:   if (scall == MAT_INITIAL_MATRIX) {
1341:     Mat_IS                *matis = (Mat_IS *)mat->data;
1342:     ISLocalToGlobalMapping rl2g;
1343:     IS                     is;
1344:     PetscInt              *lidxs, *lgidxs, *newgidxs;
1345:     PetscInt               ll, newloc, irbs, icbs, arbs, acbs, rbs, cbs;
1346:     PetscBool              cong;
1347:     MPI_Comm               comm;

1349:     PetscObjectGetComm((PetscObject)mat, &comm);
1350:     MatGetBlockSizes(mat, &arbs, &acbs);
1351:     ISGetBlockSize(irow, &irbs);
1352:     ISGetBlockSize(icol, &icbs);
1353:     rbs = arbs == irbs ? irbs : 1;
1354:     cbs = acbs == icbs ? icbs : 1;
1355:     ISGetLocalSize(irow, &m);
1356:     ISGetLocalSize(icol, &n);
1357:     MatCreate(comm, newmat);
1358:     MatSetType(*newmat, MATIS);
1359:     MatSetSizes(*newmat, m, n, PETSC_DECIDE, PETSC_DECIDE);
1360:     MatSetBlockSizes(*newmat, rbs, cbs);
1361:     /* communicate irow to their owners in the layout */
1362:     ISGetIndices(irow, &idxs);
1363:     PetscLayoutMapLocal(mat->rmap, m, idxs, &ll, &lidxs, &lgidxs);
1364:     ISRestoreIndices(irow, &idxs);
1365:     PetscArrayzero(matis->sf_rootdata, matis->sf->nroots);
1366:     for (i = 0; i < ll; i++) matis->sf_rootdata[lidxs[i]] = lgidxs[i] + 1;
1367:     PetscFree(lidxs);
1368:     PetscFree(lgidxs);
1369:     PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE);
1370:     PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE);
1371:     for (i = 0, newloc = 0; i < matis->sf->nleaves; i++)
1372:       if (matis->sf_leafdata[i]) newloc++;
1373:     PetscMalloc1(newloc, &newgidxs);
1374:     PetscMalloc1(newloc, &lidxs);
1375:     for (i = 0, newloc = 0; i < matis->sf->nleaves; i++)
1376:       if (matis->sf_leafdata[i]) {
1377:         lidxs[newloc]      = i;
1378:         newgidxs[newloc++] = matis->sf_leafdata[i] - 1;
1379:       }
1380:     ISCreateGeneral(comm, newloc, newgidxs, PETSC_OWN_POINTER, &is);
1381:     ISLocalToGlobalMappingCreateIS(is, &rl2g);
1382:     ISLocalToGlobalMappingSetBlockSize(rl2g, rbs);
1383:     ISDestroy(&is);
1384:     /* local is to extract local submatrix */
1385:     newmatis = (Mat_IS *)(*newmat)->data;
1386:     ISCreateGeneral(comm, newloc, lidxs, PETSC_OWN_POINTER, &newmatis->getsub_ris);
1387:     MatHasCongruentLayouts(mat, &cong);
1388:     if (cong && irow == icol && matis->csf == matis->sf) {
1389:       MatSetLocalToGlobalMapping(*newmat, rl2g, rl2g);
1390:       PetscObjectReference((PetscObject)newmatis->getsub_ris);
1391:       newmatis->getsub_cis = newmatis->getsub_ris;
1392:     } else {
1393:       ISLocalToGlobalMapping cl2g;

1395:       /* communicate icol to their owners in the layout */
1396:       ISGetIndices(icol, &idxs);
1397:       PetscLayoutMapLocal(mat->cmap, n, idxs, &ll, &lidxs, &lgidxs);
1398:       ISRestoreIndices(icol, &idxs);
1399:       PetscArrayzero(matis->csf_rootdata, matis->csf->nroots);
1400:       for (i = 0; i < ll; i++) matis->csf_rootdata[lidxs[i]] = lgidxs[i] + 1;
1401:       PetscFree(lidxs);
1402:       PetscFree(lgidxs);
1403:       PetscSFBcastBegin(matis->csf, MPIU_INT, matis->csf_rootdata, matis->csf_leafdata, MPI_REPLACE);
1404:       PetscSFBcastEnd(matis->csf, MPIU_INT, matis->csf_rootdata, matis->csf_leafdata, MPI_REPLACE);
1405:       for (i = 0, newloc = 0; i < matis->csf->nleaves; i++)
1406:         if (matis->csf_leafdata[i]) newloc++;
1407:       PetscMalloc1(newloc, &newgidxs);
1408:       PetscMalloc1(newloc, &lidxs);
1409:       for (i = 0, newloc = 0; i < matis->csf->nleaves; i++)
1410:         if (matis->csf_leafdata[i]) {
1411:           lidxs[newloc]      = i;
1412:           newgidxs[newloc++] = matis->csf_leafdata[i] - 1;
1413:         }
1414:       ISCreateGeneral(comm, newloc, newgidxs, PETSC_OWN_POINTER, &is);
1415:       ISLocalToGlobalMappingCreateIS(is, &cl2g);
1416:       ISLocalToGlobalMappingSetBlockSize(cl2g, cbs);
1417:       ISDestroy(&is);
1418:       /* local is to extract local submatrix */
1419:       ISCreateGeneral(comm, newloc, lidxs, PETSC_OWN_POINTER, &newmatis->getsub_cis);
1420:       MatSetLocalToGlobalMapping(*newmat, rl2g, cl2g);
1421:       ISLocalToGlobalMappingDestroy(&cl2g);
1422:     }
1423:     ISLocalToGlobalMappingDestroy(&rl2g);
1424:   } else {
1425:     MatISGetLocalMat(*newmat, &newlocmat);
1426:   }
1427:   MatISGetLocalMat(mat, &locmat);
1428:   newmatis = (Mat_IS *)(*newmat)->data;
1429:   MatCreateSubMatrix(locmat, newmatis->getsub_ris, newmatis->getsub_cis, scall, &newlocmat);
1430:   if (scall == MAT_INITIAL_MATRIX) {
1431:     MatISSetLocalMat(*newmat, newlocmat);
1432:     MatDestroy(&newlocmat);
1433:   }
1434:   MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY);
1435:   MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY);
1436:   return 0;
1437: }

1439: static PetscErrorCode MatCopy_IS(Mat A, Mat B, MatStructure str)
1440: {
1441:   Mat_IS   *a = (Mat_IS *)A->data, *b;
1442:   PetscBool ismatis;

1444:   PetscObjectTypeCompare((PetscObject)B, MATIS, &ismatis);
1446:   b = (Mat_IS *)B->data;
1447:   MatCopy(a->A, b->A, str);
1448:   PetscObjectStateIncrease((PetscObject)B);
1449:   return 0;
1450: }

1452: static PetscErrorCode MatMissingDiagonal_IS(Mat A, PetscBool *missing, PetscInt *d)
1453: {
1454:   Vec                v;
1455:   const PetscScalar *array;
1456:   PetscInt           i, n;

1458:   *missing = PETSC_FALSE;
1459:   MatCreateVecs(A, NULL, &v);
1460:   MatGetDiagonal(A, v);
1461:   VecGetLocalSize(v, &n);
1462:   VecGetArrayRead(v, &array);
1463:   for (i = 0; i < n; i++)
1464:     if (array[i] == 0.) break;
1465:   VecRestoreArrayRead(v, &array);
1466:   VecDestroy(&v);
1467:   if (i != n) *missing = PETSC_TRUE;
1468:   if (d) {
1469:     *d = -1;
1470:     if (*missing) {
1471:       PetscInt rstart;
1472:       MatGetOwnershipRange(A, &rstart, NULL);
1473:       *d = i + rstart;
1474:     }
1475:   }
1476:   return 0;
1477: }

1479: static PetscErrorCode MatISSetUpSF_IS(Mat B)
1480: {
1481:   Mat_IS         *matis = (Mat_IS *)(B->data);
1482:   const PetscInt *gidxs;
1483:   PetscInt        nleaves;

1485:   if (matis->sf) return 0;
1486:   PetscSFCreate(PetscObjectComm((PetscObject)B), &matis->sf);
1487:   ISLocalToGlobalMappingGetIndices(matis->rmapping, &gidxs);
1488:   ISLocalToGlobalMappingGetSize(matis->rmapping, &nleaves);
1489:   PetscSFSetGraphLayout(matis->sf, B->rmap, nleaves, NULL, PETSC_OWN_POINTER, gidxs);
1490:   ISLocalToGlobalMappingRestoreIndices(matis->rmapping, &gidxs);
1491:   PetscMalloc2(matis->sf->nroots, &matis->sf_rootdata, matis->sf->nleaves, &matis->sf_leafdata);
1492:   if (matis->rmapping != matis->cmapping) { /* setup SF for columns */
1493:     ISLocalToGlobalMappingGetSize(matis->cmapping, &nleaves);
1494:     PetscSFCreate(PetscObjectComm((PetscObject)B), &matis->csf);
1495:     ISLocalToGlobalMappingGetIndices(matis->cmapping, &gidxs);
1496:     PetscSFSetGraphLayout(matis->csf, B->cmap, nleaves, NULL, PETSC_OWN_POINTER, gidxs);
1497:     ISLocalToGlobalMappingRestoreIndices(matis->cmapping, &gidxs);
1498:     PetscMalloc2(matis->csf->nroots, &matis->csf_rootdata, matis->csf->nleaves, &matis->csf_leafdata);
1499:   } else {
1500:     matis->csf          = matis->sf;
1501:     matis->csf_leafdata = matis->sf_leafdata;
1502:     matis->csf_rootdata = matis->sf_rootdata;
1503:   }
1504:   return 0;
1505: }

1507: /*@
1508:    MatISStoreL2L - Store local-to-local operators during the Galerkin process of computing MatPtAP.

1510:    Collective

1512:    Input Parameters:
1513: +  A - the matrix
1514: -  store - the boolean flag

1516:    Level: advanced

1518: .seealso: `MatCreate()`, `MatCreateIS()`, `MatISSetPreallocation()`, `MatPtAP()`
1519: @*/
1520: PetscErrorCode MatISStoreL2L(Mat A, PetscBool store)
1521: {
1525:   PetscTryMethod(A, "MatISStoreL2L_C", (Mat, PetscBool), (A, store));
1526:   return 0;
1527: }

1529: static PetscErrorCode MatISStoreL2L_IS(Mat A, PetscBool store)
1530: {
1531:   Mat_IS *matis = (Mat_IS *)(A->data);

1533:   matis->storel2l = store;
1534:   if (!store) PetscObjectCompose((PetscObject)(A), "_MatIS_PtAP_l2l", NULL);
1535:   return 0;
1536: }

1538: /*@
1539:    MatISFixLocalEmpty - Compress out zero local rows from the local matrices

1541:    Collective

1543:    Input Parameters:
1544: +  A - the matrix
1545: -  fix - the boolean flag

1547:    Level: advanced

1549:    Note:
1550:    When fix is true, new local matrices and l2g maps are generated during the final assembly process.

1552: .seealso: `MATIS`, `MatCreate()`, `MatCreateIS()`, `MatISSetPreallocation()`, `MatAssemblyEnd()`, `MAT_FINAL_ASSEMBLY`
1553: @*/
1554: PetscErrorCode MatISFixLocalEmpty(Mat A, PetscBool fix)
1555: {
1559:   PetscTryMethod(A, "MatISFixLocalEmpty_C", (Mat, PetscBool), (A, fix));
1560:   return 0;
1561: }

1563: static PetscErrorCode MatISFixLocalEmpty_IS(Mat A, PetscBool fix)
1564: {
1565:   Mat_IS *matis = (Mat_IS *)(A->data);

1567:   matis->locempty = fix;
1568:   return 0;
1569: }

1571: /*@
1572:    MatISSetPreallocation - Preallocates memory for a `MATIS` parallel matrix.

1574:    Collective

1576:    Input Parameters:
1577: +  B - the matrix
1578: .  d_nz  - number of nonzeros per row in DIAGONAL portion of local submatrix
1579:            (same value is used for all local rows)
1580: .  d_nnz - array containing the number of nonzeros in the various rows of the
1581:            DIAGONAL portion of the local submatrix (possibly different for each row)
1582:            or NULL, if d_nz is used to specify the nonzero structure.
1583:            The size of this array is equal to the number of local rows, i.e 'm'.
1584:            For matrices that will be factored, you must leave room for (and set)
1585:            the diagonal entry even if it is zero.
1586: .  o_nz  - number of nonzeros per row in the OFF-DIAGONAL portion of local
1587:            submatrix (same value is used for all local rows).
1588: -  o_nnz - array containing the number of nonzeros in the various rows of the
1589:            OFF-DIAGONAL portion of the local submatrix (possibly different for
1590:            each row) or NULL, if o_nz is used to specify the nonzero
1591:            structure. The size of this array is equal to the number
1592:            of local rows, i.e 'm'.

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

1596:    Level: intermediate

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

1603: .seealso: `MatCreate()`, `MatCreateIS()`, `MatMPIAIJSetPreallocation()`, `MatISGetLocalMat()`, `MATIS`
1604: @*/
1605: PetscErrorCode MatISSetPreallocation(Mat B, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[])
1606: {
1609:   PetscTryMethod(B, "MatISSetPreallocation_C", (Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[]), (B, d_nz, d_nnz, o_nz, o_nnz));
1610:   return 0;
1611: }

1613: /* this is used by DMDA */
1614: PETSC_EXTERN PetscErrorCode MatISSetPreallocation_IS(Mat B, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[])
1615: {
1616:   Mat_IS  *matis = (Mat_IS *)(B->data);
1617:   PetscInt bs, i, nlocalcols;

1619:   MatSetUp(B);
1620:   if (!d_nnz)
1621:     for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] = d_nz;
1622:   else
1623:     for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] = d_nnz[i];

1625:   if (!o_nnz)
1626:     for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] += o_nz;
1627:   else
1628:     for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] += o_nnz[i];

1630:   PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE);
1631:   MatGetSize(matis->A, NULL, &nlocalcols);
1632:   MatGetBlockSize(matis->A, &bs);
1633:   PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE);

1635:   for (i = 0; i < matis->sf->nleaves; i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i], nlocalcols);
1636:   MatSeqAIJSetPreallocation(matis->A, 0, matis->sf_leafdata);
1637: #if defined(PETSC_HAVE_HYPRE)
1638:   MatHYPRESetPreallocation(matis->A, 0, matis->sf_leafdata, 0, NULL);
1639: #endif

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

1644:     matis->sf_leafdata[i] = matis->sf_leafdata[i * bs] / bs;
1645:     for (b = 1; b < bs; b++) matis->sf_leafdata[i] = PetscMax(matis->sf_leafdata[i], matis->sf_leafdata[i * bs + b] / bs);
1646:   }
1647:   MatSeqBAIJSetPreallocation(matis->A, bs, 0, matis->sf_leafdata);

1649:   nlocalcols /= bs;
1650:   for (i = 0; i < matis->sf->nleaves / bs; i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i], nlocalcols - i);
1651:   MatSeqSBAIJSetPreallocation(matis->A, bs, 0, matis->sf_leafdata);

1653:   /* for other matrix types */
1654:   MatSetUp(matis->A);
1655:   return 0;
1656: }

1658: PETSC_EXTERN PetscErrorCode MatISSetMPIXAIJPreallocation_Private(Mat A, Mat B, PetscBool maxreduce)
1659: {
1660:   Mat_IS         *matis = (Mat_IS *)(A->data);
1661:   PetscInt       *my_dnz, *my_onz, *dnz, *onz, *mat_ranges, *row_ownership;
1662:   const PetscInt *global_indices_r, *global_indices_c;
1663:   PetscInt        i, j, bs, rows, cols;
1664:   PetscInt        lrows, lcols;
1665:   PetscInt        local_rows, local_cols;
1666:   PetscMPIInt     size;
1667:   PetscBool       isdense, issbaij;

1669:   MPI_Comm_size(PetscObjectComm((PetscObject)A), &size);
1670:   MatGetSize(A, &rows, &cols);
1671:   MatGetBlockSize(A, &bs);
1672:   MatGetSize(matis->A, &local_rows, &local_cols);
1673:   PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQDENSE, &isdense);
1674:   PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQSBAIJ, &issbaij);
1675:   ISLocalToGlobalMappingGetIndices(matis->rmapping, &global_indices_r);
1676:   if (matis->rmapping != matis->cmapping) {
1677:     ISLocalToGlobalMappingGetIndices(matis->cmapping, &global_indices_c);
1678:   } else global_indices_c = global_indices_r;

1680:   if (issbaij) MatGetRowUpperTriangular(matis->A);
1681:   /*
1682:      An SF reduce is needed to sum up properly on shared rows.
1683:      Note that generally preallocation is not exact, since it overestimates nonzeros
1684:   */
1685:   MatGetLocalSize(A, &lrows, &lcols);
1686:   MatPreallocateBegin(PetscObjectComm((PetscObject)A), lrows, lcols, dnz, onz);
1687:   /* All processes need to compute entire row ownership */
1688:   PetscMalloc1(rows, &row_ownership);
1689:   MatGetOwnershipRanges(A, (const PetscInt **)&mat_ranges);
1690:   for (i = 0; i < size; i++) {
1691:     for (j = mat_ranges[i]; j < mat_ranges[i + 1]; j++) row_ownership[j] = i;
1692:   }
1693:   MatGetOwnershipRangesColumn(A, (const PetscInt **)&mat_ranges);

1695:   /*
1696:      my_dnz and my_onz contains exact contribution to preallocation from each local mat
1697:      then, they will be summed up properly. This way, preallocation is always sufficient
1698:   */
1699:   PetscCalloc2(local_rows, &my_dnz, local_rows, &my_onz);
1700:   /* preallocation as a MATAIJ */
1701:   if (isdense) { /* special case for dense local matrices */
1702:     for (i = 0; i < local_rows; i++) {
1703:       PetscInt owner = row_ownership[global_indices_r[i]];
1704:       for (j = 0; j < local_cols; j++) {
1705:         PetscInt index_col = global_indices_c[j];
1706:         if (index_col > mat_ranges[owner] - 1 && index_col < mat_ranges[owner + 1]) { /* diag block */
1707:           my_dnz[i] += 1;
1708:         } else { /* offdiag block */
1709:           my_onz[i] += 1;
1710:         }
1711:       }
1712:     }
1713:   } else if (matis->A->ops->getrowij) {
1714:     const PetscInt *ii, *jj, *jptr;
1715:     PetscBool       done;
1716:     MatGetRowIJ(matis->A, 0, PETSC_FALSE, PETSC_FALSE, &local_rows, &ii, &jj, &done);
1718:     jptr = jj;
1719:     for (i = 0; i < local_rows; i++) {
1720:       PetscInt index_row = global_indices_r[i];
1721:       for (j = 0; j < ii[i + 1] - ii[i]; j++, jptr++) {
1722:         PetscInt owner     = row_ownership[index_row];
1723:         PetscInt index_col = global_indices_c[*jptr];
1724:         if (index_col > mat_ranges[owner] - 1 && index_col < mat_ranges[owner + 1]) { /* diag block */
1725:           my_dnz[i] += 1;
1726:         } else { /* offdiag block */
1727:           my_onz[i] += 1;
1728:         }
1729:         /* same as before, interchanging rows and cols */
1730:         if (issbaij && index_col != index_row) {
1731:           owner = row_ownership[index_col];
1732:           if (index_row > mat_ranges[owner] - 1 && index_row < mat_ranges[owner + 1]) {
1733:             my_dnz[*jptr] += 1;
1734:           } else {
1735:             my_onz[*jptr] += 1;
1736:           }
1737:         }
1738:       }
1739:     }
1740:     MatRestoreRowIJ(matis->A, 0, PETSC_FALSE, PETSC_FALSE, &local_rows, &ii, &jj, &done);
1742:   } else { /* loop over rows and use MatGetRow */
1743:     for (i = 0; i < local_rows; i++) {
1744:       const PetscInt *cols;
1745:       PetscInt        ncols, index_row = global_indices_r[i];
1746:       MatGetRow(matis->A, i, &ncols, &cols, NULL);
1747:       for (j = 0; j < ncols; j++) {
1748:         PetscInt owner     = row_ownership[index_row];
1749:         PetscInt index_col = global_indices_c[cols[j]];
1750:         if (index_col > mat_ranges[owner] - 1 && index_col < mat_ranges[owner + 1]) { /* diag block */
1751:           my_dnz[i] += 1;
1752:         } else { /* offdiag block */
1753:           my_onz[i] += 1;
1754:         }
1755:         /* same as before, interchanging rows and cols */
1756:         if (issbaij && index_col != index_row) {
1757:           owner = row_ownership[index_col];
1758:           if (index_row > mat_ranges[owner] - 1 && index_row < mat_ranges[owner + 1]) {
1759:             my_dnz[cols[j]] += 1;
1760:           } else {
1761:             my_onz[cols[j]] += 1;
1762:           }
1763:         }
1764:       }
1765:       MatRestoreRow(matis->A, i, &ncols, &cols, NULL);
1766:     }
1767:   }
1768:   if (global_indices_c != global_indices_r) ISLocalToGlobalMappingRestoreIndices(matis->cmapping, &global_indices_c);
1769:   ISLocalToGlobalMappingRestoreIndices(matis->rmapping, &global_indices_r);
1770:   PetscFree(row_ownership);

1772:   /* Reduce my_dnz and my_onz */
1773:   if (maxreduce) {
1774:     PetscSFReduceBegin(matis->sf, MPIU_INT, my_dnz, dnz, MPI_MAX);
1775:     PetscSFReduceBegin(matis->sf, MPIU_INT, my_onz, onz, MPI_MAX);
1776:     PetscSFReduceEnd(matis->sf, MPIU_INT, my_dnz, dnz, MPI_MAX);
1777:     PetscSFReduceEnd(matis->sf, MPIU_INT, my_onz, onz, MPI_MAX);
1778:   } else {
1779:     PetscSFReduceBegin(matis->sf, MPIU_INT, my_dnz, dnz, MPI_SUM);
1780:     PetscSFReduceBegin(matis->sf, MPIU_INT, my_onz, onz, MPI_SUM);
1781:     PetscSFReduceEnd(matis->sf, MPIU_INT, my_dnz, dnz, MPI_SUM);
1782:     PetscSFReduceEnd(matis->sf, MPIU_INT, my_onz, onz, MPI_SUM);
1783:   }
1784:   PetscFree2(my_dnz, my_onz);

1786:   /* Resize preallocation if overestimated */
1787:   for (i = 0; i < lrows; i++) {
1788:     dnz[i] = PetscMin(dnz[i], lcols);
1789:     onz[i] = PetscMin(onz[i], cols - lcols);
1790:   }

1792:   /* Set preallocation */
1793:   MatSetBlockSizesFromMats(B, A, A);
1794:   MatSeqAIJSetPreallocation(B, 0, dnz);
1795:   MatMPIAIJSetPreallocation(B, 0, dnz, 0, onz);
1796:   for (i = 0; i < lrows; i += bs) {
1797:     PetscInt b, d = dnz[i], o = onz[i];

1799:     for (b = 1; b < bs; b++) {
1800:       d = PetscMax(d, dnz[i + b]);
1801:       o = PetscMax(o, onz[i + b]);
1802:     }
1803:     dnz[i / bs] = PetscMin(d / bs + d % bs, lcols / bs);
1804:     onz[i / bs] = PetscMin(o / bs + o % bs, (cols - lcols) / bs);
1805:   }
1806:   MatSeqBAIJSetPreallocation(B, bs, 0, dnz);
1807:   MatMPIBAIJSetPreallocation(B, bs, 0, dnz, 0, onz);
1808:   MatMPISBAIJSetPreallocation(B, bs, 0, dnz, 0, onz);
1809:   MatPreallocateEnd(dnz, onz);
1810:   if (issbaij) MatRestoreRowUpperTriangular(matis->A);
1811:   MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
1812:   return 0;
1813: }

1815: PETSC_INTERN PetscErrorCode MatConvert_IS_XAIJ(Mat mat, MatType mtype, MatReuse reuse, Mat *M)
1816: {
1817:   Mat_IS            *matis = (Mat_IS *)(mat->data);
1818:   Mat                local_mat, MT;
1819:   PetscInt           rbs, cbs, rows, cols, lrows, lcols;
1820:   PetscInt           local_rows, local_cols;
1821:   PetscBool          isseqdense, isseqsbaij, isseqaij, isseqbaij;
1822:   PetscMPIInt        size;
1823:   const PetscScalar *array;

1825:   MPI_Comm_size(PetscObjectComm((PetscObject)mat), &size);
1826:   if (size == 1 && mat->rmap->N == matis->A->rmap->N && mat->cmap->N == matis->A->cmap->N) {
1827:     Mat      B;
1828:     IS       irows = NULL, icols = NULL;
1829:     PetscInt rbs, cbs;

1831:     ISLocalToGlobalMappingGetBlockSize(matis->rmapping, &rbs);
1832:     ISLocalToGlobalMappingGetBlockSize(matis->cmapping, &cbs);
1833:     if (reuse != MAT_REUSE_MATRIX) { /* check if l2g maps are one-to-one */
1834:       IS              rows, cols;
1835:       const PetscInt *ridxs, *cidxs;
1836:       PetscInt        i, nw, *work;

1838:       ISLocalToGlobalMappingGetBlockIndices(matis->rmapping, &ridxs);
1839:       ISLocalToGlobalMappingGetSize(matis->rmapping, &nw);
1840:       nw = nw / rbs;
1841:       PetscCalloc1(nw, &work);
1842:       for (i = 0; i < nw; i++) work[ridxs[i]] += 1;
1843:       for (i = 0; i < nw; i++)
1844:         if (!work[i] || work[i] > 1) break;
1845:       if (i == nw) {
1846:         ISCreateBlock(PETSC_COMM_SELF, rbs, nw, ridxs, PETSC_USE_POINTER, &rows);
1847:         ISSetPermutation(rows);
1848:         ISInvertPermutation(rows, PETSC_DECIDE, &irows);
1849:         ISDestroy(&rows);
1850:       }
1851:       ISLocalToGlobalMappingRestoreBlockIndices(matis->rmapping, &ridxs);
1852:       PetscFree(work);
1853:       if (irows && matis->rmapping != matis->cmapping) {
1854:         ISLocalToGlobalMappingGetBlockIndices(matis->cmapping, &cidxs);
1855:         ISLocalToGlobalMappingGetSize(matis->cmapping, &nw);
1856:         nw = nw / cbs;
1857:         PetscCalloc1(nw, &work);
1858:         for (i = 0; i < nw; i++) work[cidxs[i]] += 1;
1859:         for (i = 0; i < nw; i++)
1860:           if (!work[i] || work[i] > 1) break;
1861:         if (i == nw) {
1862:           ISCreateBlock(PETSC_COMM_SELF, cbs, nw, cidxs, PETSC_USE_POINTER, &cols);
1863:           ISSetPermutation(cols);
1864:           ISInvertPermutation(cols, PETSC_DECIDE, &icols);
1865:           ISDestroy(&cols);
1866:         }
1867:         ISLocalToGlobalMappingRestoreBlockIndices(matis->cmapping, &cidxs);
1868:         PetscFree(work);
1869:       } else if (irows) {
1870:         PetscObjectReference((PetscObject)irows);
1871:         icols = irows;
1872:       }
1873:     } else {
1874:       PetscObjectQuery((PetscObject)(*M), "_MatIS_IS_XAIJ_irows", (PetscObject *)&irows);
1875:       PetscObjectQuery((PetscObject)(*M), "_MatIS_IS_XAIJ_icols", (PetscObject *)&icols);
1876:       if (irows) PetscObjectReference((PetscObject)irows);
1877:       if (icols) PetscObjectReference((PetscObject)icols);
1878:     }
1879:     if (!irows || !icols) {
1880:       ISDestroy(&icols);
1881:       ISDestroy(&irows);
1882:       goto general_assembly;
1883:     }
1884:     MatConvert(matis->A, mtype, MAT_INITIAL_MATRIX, &B);
1885:     if (reuse != MAT_INPLACE_MATRIX) {
1886:       MatCreateSubMatrix(B, irows, icols, reuse, M);
1887:       PetscObjectCompose((PetscObject)(*M), "_MatIS_IS_XAIJ_irows", (PetscObject)irows);
1888:       PetscObjectCompose((PetscObject)(*M), "_MatIS_IS_XAIJ_icols", (PetscObject)icols);
1889:     } else {
1890:       Mat C;

1892:       MatCreateSubMatrix(B, irows, icols, MAT_INITIAL_MATRIX, &C);
1893:       MatHeaderReplace(mat, &C);
1894:     }
1895:     MatDestroy(&B);
1896:     ISDestroy(&icols);
1897:     ISDestroy(&irows);
1898:     return 0;
1899:   }
1900: general_assembly:
1901:   MatGetSize(mat, &rows, &cols);
1902:   ISLocalToGlobalMappingGetBlockSize(matis->rmapping, &rbs);
1903:   ISLocalToGlobalMappingGetBlockSize(matis->cmapping, &cbs);
1904:   MatGetLocalSize(mat, &lrows, &lcols);
1905:   MatGetSize(matis->A, &local_rows, &local_cols);
1906:   PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQDENSE, &isseqdense);
1907:   PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQAIJ, &isseqaij);
1908:   PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQBAIJ, &isseqbaij);
1909:   PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQSBAIJ, &isseqsbaij);
1911:   if (PetscDefined(USE_DEBUG)) {
1912:     PetscBool lb[4], bb[4];

1914:     lb[0] = isseqdense;
1915:     lb[1] = isseqaij;
1916:     lb[2] = isseqbaij;
1917:     lb[3] = isseqsbaij;
1918:     MPIU_Allreduce(lb, bb, 4, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)mat));
1920:   }

1922:   if (reuse != MAT_REUSE_MATRIX) {
1923:     MatCreate(PetscObjectComm((PetscObject)mat), &MT);
1924:     MatSetSizes(MT, lrows, lcols, rows, cols);
1925:     MatSetType(MT, mtype);
1926:     MatSetBlockSizes(MT, rbs, cbs);
1927:     MatISSetMPIXAIJPreallocation_Private(mat, MT, PETSC_FALSE);
1928:   } else {
1929:     PetscInt mrbs, mcbs, mrows, mcols, mlrows, mlcols;

1931:     /* some checks */
1932:     MT = *M;
1933:     MatGetBlockSizes(MT, &mrbs, &mcbs);
1934:     MatGetSize(MT, &mrows, &mcols);
1935:     MatGetLocalSize(MT, &mlrows, &mlcols);
1942:     MatZeroEntries(MT);
1943:   }

1945:   if (isseqsbaij || isseqbaij) {
1946:     MatConvert(matis->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &local_mat);
1947:     isseqaij = PETSC_TRUE;
1948:   } else {
1949:     PetscObjectReference((PetscObject)matis->A);
1950:     local_mat = matis->A;
1951:   }

1953:   /* Set values */
1954:   MatSetLocalToGlobalMapping(MT, matis->rmapping, matis->cmapping);
1955:   if (isseqdense) { /* special case for dense local matrices */
1956:     PetscInt i, *dummy;

1958:     PetscMalloc1(PetscMax(local_rows, local_cols), &dummy);
1959:     for (i = 0; i < PetscMax(local_rows, local_cols); i++) dummy[i] = i;
1960:     MatSetOption(MT, MAT_ROW_ORIENTED, PETSC_FALSE);
1961:     MatDenseGetArrayRead(local_mat, &array);
1962:     MatSetValuesLocal(MT, local_rows, dummy, local_cols, dummy, array, ADD_VALUES);
1963:     MatDenseRestoreArrayRead(local_mat, &array);
1964:     PetscFree(dummy);
1965:   } else if (isseqaij) {
1966:     const PetscInt *blocks;
1967:     PetscInt        i, nvtxs, *xadj, *adjncy, nb;
1968:     PetscBool       done;
1969:     PetscScalar    *sarray;

1971:     MatGetRowIJ(local_mat, 0, PETSC_FALSE, PETSC_FALSE, &nvtxs, (const PetscInt **)&xadj, (const PetscInt **)&adjncy, &done);
1973:     MatSeqAIJGetArray(local_mat, &sarray);
1974:     MatGetVariableBlockSizes(local_mat, &nb, &blocks);
1975:     if (nb) { /* speed up assembly for special blocked matrices (used by BDDC) */
1976:       PetscInt sum;

1978:       for (i = 0, sum = 0; i < nb; i++) sum += blocks[i];
1979:       if (sum == nvtxs) {
1980:         PetscInt r;

1982:         for (i = 0, r = 0; i < nb; i++) {
1983:           PetscAssert(blocks[i] == xadj[r + 1] - xadj[r], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid block sizes prescribed for block %" PetscInt_FMT ": expected %" PetscInt_FMT ", got %" PetscInt_FMT, i, blocks[i], xadj[r + 1] - xadj[r]);
1984:           MatSetValuesLocal(MT, blocks[i], adjncy + xadj[r], blocks[i], adjncy + xadj[r], sarray + xadj[r], ADD_VALUES);
1985:           r += blocks[i];
1986:         }
1987:       } else {
1988:         for (i = 0; i < nvtxs; i++) MatSetValuesLocal(MT, 1, &i, xadj[i + 1] - xadj[i], adjncy + xadj[i], sarray + xadj[i], ADD_VALUES);
1989:       }
1990:     } else {
1991:       for (i = 0; i < nvtxs; i++) MatSetValuesLocal(MT, 1, &i, xadj[i + 1] - xadj[i], adjncy + xadj[i], sarray + xadj[i], ADD_VALUES);
1992:     }
1993:     MatRestoreRowIJ(local_mat, 0, PETSC_FALSE, PETSC_FALSE, &nvtxs, (const PetscInt **)&xadj, (const PetscInt **)&adjncy, &done);
1995:     MatSeqAIJRestoreArray(local_mat, &sarray);
1996:   } else { /* very basic values insertion for all other matrix types */
1997:     PetscInt i;

1999:     for (i = 0; i < local_rows; i++) {
2000:       PetscInt        j;
2001:       const PetscInt *local_indices_cols;

2003:       MatGetRow(local_mat, i, &j, &local_indices_cols, &array);
2004:       MatSetValuesLocal(MT, 1, &i, j, local_indices_cols, array, ADD_VALUES);
2005:       MatRestoreRow(local_mat, i, &j, &local_indices_cols, &array);
2006:     }
2007:   }
2008:   MatAssemblyBegin(MT, MAT_FINAL_ASSEMBLY);
2009:   MatDestroy(&local_mat);
2010:   MatAssemblyEnd(MT, MAT_FINAL_ASSEMBLY);
2011:   if (isseqdense) MatSetOption(MT, MAT_ROW_ORIENTED, PETSC_TRUE);
2012:   if (reuse == MAT_INPLACE_MATRIX) {
2013:     MatHeaderReplace(mat, &MT);
2014:   } else if (reuse == MAT_INITIAL_MATRIX) {
2015:     *M = MT;
2016:   }
2017:   return 0;
2018: }

2020: /*@
2021:     MatISGetMPIXAIJ - Converts `MATIS` matrix into a parallel `MATAIJ` format

2023:   Input Parameters:
2024: +  mat - the matrix (should be of type `MATIS`)
2025: -  reuse - either `MAT_INITIAL_MATRIX` or `MAT_REUSE_MATRIX`

2027:   Output Parameter:
2028: .  newmat - the matrix in `MATAIJ` format

2030:   Level: developer

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

2035: .seealso: `MATIS`, `MatConvert()`
2036: @*/
2037: PetscErrorCode MatISGetMPIXAIJ(Mat mat, MatReuse reuse, Mat *newmat)
2038: {
2042:   if (reuse == MAT_REUSE_MATRIX) {
2046:   }
2047:   PetscUseMethod(mat, "MatISGetMPIXAIJ_C", (Mat, MatType, MatReuse, Mat *), (mat, MATAIJ, reuse, newmat));
2048:   return 0;
2049: }

2051: static PetscErrorCode MatDuplicate_IS(Mat mat, MatDuplicateOption op, Mat *newmat)
2052: {
2053:   Mat_IS  *matis = (Mat_IS *)(mat->data);
2054:   PetscInt rbs, cbs, m, n, M, N;
2055:   Mat      B, localmat;

2057:   ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping, &rbs);
2058:   ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping, &cbs);
2059:   MatGetSize(mat, &M, &N);
2060:   MatGetLocalSize(mat, &m, &n);
2061:   MatCreate(PetscObjectComm((PetscObject)mat), &B);
2062:   MatSetSizes(B, m, n, M, N);
2063:   MatSetBlockSize(B, rbs == cbs ? rbs : 1);
2064:   MatSetType(B, MATIS);
2065:   MatISSetLocalMatType(B, matis->lmattype);
2066:   MatSetLocalToGlobalMapping(B, mat->rmap->mapping, mat->cmap->mapping);
2067:   MatDuplicate(matis->A, op, &localmat);
2068:   MatSetLocalToGlobalMapping(localmat, matis->A->rmap->mapping, matis->A->cmap->mapping);
2069:   MatISSetLocalMat(B, localmat);
2070:   MatDestroy(&localmat);
2071:   MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
2072:   MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
2073:   *newmat = B;
2074:   return 0;
2075: }

2077: static PetscErrorCode MatIsHermitian_IS(Mat A, PetscReal tol, PetscBool *flg)
2078: {
2079:   Mat_IS   *matis = (Mat_IS *)A->data;
2080:   PetscBool local_sym;

2082:   MatIsHermitian(matis->A, tol, &local_sym);
2083:   MPIU_Allreduce(&local_sym, flg, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A));
2084:   return 0;
2085: }

2087: static PetscErrorCode MatIsSymmetric_IS(Mat A, PetscReal tol, PetscBool *flg)
2088: {
2089:   Mat_IS   *matis = (Mat_IS *)A->data;
2090:   PetscBool local_sym;

2092:   if (matis->rmapping != matis->cmapping) {
2093:     *flg = PETSC_FALSE;
2094:     return 0;
2095:   }
2096:   MatIsSymmetric(matis->A, tol, &local_sym);
2097:   MPIU_Allreduce(&local_sym, flg, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A));
2098:   return 0;
2099: }

2101: static PetscErrorCode MatIsStructurallySymmetric_IS(Mat A, PetscBool *flg)
2102: {
2103:   Mat_IS   *matis = (Mat_IS *)A->data;
2104:   PetscBool local_sym;

2106:   if (matis->rmapping != matis->cmapping) {
2107:     *flg = PETSC_FALSE;
2108:     return 0;
2109:   }
2110:   MatIsStructurallySymmetric(matis->A, &local_sym);
2111:   MPIU_Allreduce(&local_sym, flg, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A));
2112:   return 0;
2113: }

2115: static PetscErrorCode MatDestroy_IS(Mat A)
2116: {
2117:   Mat_IS *b = (Mat_IS *)A->data;

2119:   PetscFree(b->bdiag);
2120:   PetscFree(b->lmattype);
2121:   MatDestroy(&b->A);
2122:   VecScatterDestroy(&b->cctx);
2123:   VecScatterDestroy(&b->rctx);
2124:   VecDestroy(&b->x);
2125:   VecDestroy(&b->y);
2126:   VecDestroy(&b->counter);
2127:   ISDestroy(&b->getsub_ris);
2128:   ISDestroy(&b->getsub_cis);
2129:   if (b->sf != b->csf) {
2130:     PetscSFDestroy(&b->csf);
2131:     PetscFree2(b->csf_rootdata, b->csf_leafdata);
2132:   } else b->csf = NULL;
2133:   PetscSFDestroy(&b->sf);
2134:   PetscFree2(b->sf_rootdata, b->sf_leafdata);
2135:   ISLocalToGlobalMappingDestroy(&b->rmapping);
2136:   ISLocalToGlobalMappingDestroy(&b->cmapping);
2137:   MatDestroy(&b->dA);
2138:   MatDestroy(&b->assembledA);
2139:   PetscFree(A->data);
2140:   PetscObjectChangeTypeName((PetscObject)A, NULL);
2141:   PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMatType_C", NULL);
2142:   PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalMat_C", NULL);
2143:   PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMat_C", NULL);
2144:   PetscObjectComposeFunction((PetscObject)A, "MatISRestoreLocalMat_C", NULL);
2145:   PetscObjectComposeFunction((PetscObject)A, "MatISGetMPIXAIJ_C", NULL);
2146:   PetscObjectComposeFunction((PetscObject)A, "MatISSetPreallocation_C", NULL);
2147:   PetscObjectComposeFunction((PetscObject)A, "MatISStoreL2L_C", NULL);
2148:   PetscObjectComposeFunction((PetscObject)A, "MatISFixLocalEmpty_C", NULL);
2149:   PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalToGlobalMapping_C", NULL);
2150:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpiaij_C", NULL);
2151:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpibaij_C", NULL);
2152:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpisbaij_C", NULL);
2153:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqaij_C", NULL);
2154:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqbaij_C", NULL);
2155:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqsbaij_C", NULL);
2156:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_aij_C", NULL);
2157:   PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOOLocal_C", NULL);
2158:   PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOO_C", NULL);
2159:   PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", NULL);
2160:   return 0;
2161: }

2163: static PetscErrorCode MatMult_IS(Mat A, Vec x, Vec y)
2164: {
2165:   Mat_IS     *is   = (Mat_IS *)A->data;
2166:   PetscScalar zero = 0.0;

2168:   /*  scatter the global vector x into the local work vector */
2169:   VecScatterBegin(is->cctx, x, is->x, INSERT_VALUES, SCATTER_FORWARD);
2170:   VecScatterEnd(is->cctx, x, is->x, INSERT_VALUES, SCATTER_FORWARD);

2172:   /* multiply the local matrix */
2173:   MatMult(is->A, is->x, is->y);

2175:   /* scatter product back into global memory */
2176:   VecSet(y, zero);
2177:   VecScatterBegin(is->rctx, is->y, y, ADD_VALUES, SCATTER_REVERSE);
2178:   VecScatterEnd(is->rctx, is->y, y, ADD_VALUES, SCATTER_REVERSE);
2179:   return 0;
2180: }

2182: static PetscErrorCode MatMultAdd_IS(Mat A, Vec v1, Vec v2, Vec v3)
2183: {
2184:   Vec temp_vec;

2186:   if (v3 != v2) {
2187:     MatMult(A, v1, v3);
2188:     VecAXPY(v3, 1.0, v2);
2189:   } else {
2190:     VecDuplicate(v2, &temp_vec);
2191:     MatMult(A, v1, temp_vec);
2192:     VecAXPY(temp_vec, 1.0, v2);
2193:     VecCopy(temp_vec, v3);
2194:     VecDestroy(&temp_vec);
2195:   }
2196:   return 0;
2197: }

2199: static PetscErrorCode MatMultTranspose_IS(Mat A, Vec y, Vec x)
2200: {
2201:   Mat_IS *is = (Mat_IS *)A->data;

2203:   /*  scatter the global vector x into the local work vector */
2204:   VecScatterBegin(is->rctx, y, is->y, INSERT_VALUES, SCATTER_FORWARD);
2205:   VecScatterEnd(is->rctx, y, is->y, INSERT_VALUES, SCATTER_FORWARD);

2207:   /* multiply the local matrix */
2208:   MatMultTranspose(is->A, is->y, is->x);

2210:   /* scatter product back into global vector */
2211:   VecSet(x, 0);
2212:   VecScatterBegin(is->cctx, is->x, x, ADD_VALUES, SCATTER_REVERSE);
2213:   VecScatterEnd(is->cctx, is->x, x, ADD_VALUES, SCATTER_REVERSE);
2214:   return 0;
2215: }

2217: static PetscErrorCode MatMultTransposeAdd_IS(Mat A, Vec v1, Vec v2, Vec v3)
2218: {
2219:   Vec temp_vec;

2221:   if (v3 != v2) {
2222:     MatMultTranspose(A, v1, v3);
2223:     VecAXPY(v3, 1.0, v2);
2224:   } else {
2225:     VecDuplicate(v2, &temp_vec);
2226:     MatMultTranspose(A, v1, temp_vec);
2227:     VecAXPY(temp_vec, 1.0, v2);
2228:     VecCopy(temp_vec, v3);
2229:     VecDestroy(&temp_vec);
2230:   }
2231:   return 0;
2232: }

2234: static PetscErrorCode MatView_IS(Mat A, PetscViewer viewer)
2235: {
2236:   Mat_IS     *a = (Mat_IS *)A->data;
2237:   PetscViewer sviewer;
2238:   PetscBool   isascii, view = PETSC_TRUE;

2240:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
2241:   if (isascii) {
2242:     PetscViewerFormat format;

2244:     PetscViewerGetFormat(viewer, &format);
2245:     if (format == PETSC_VIEWER_ASCII_INFO) view = PETSC_FALSE;
2246:   }
2247:   if (!view) return 0;
2248:   PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
2249:   MatView(a->A, sviewer);
2250:   PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
2251:   PetscViewerFlush(viewer);
2252:   return 0;
2253: }

2255: static PetscErrorCode MatInvertBlockDiagonal_IS(Mat mat, const PetscScalar **values)
2256: {
2257:   Mat_IS            *is = (Mat_IS *)mat->data;
2258:   MPI_Datatype       nodeType;
2259:   const PetscScalar *lv;
2260:   PetscInt           bs;

2262:   MatGetBlockSize(mat, &bs);
2263:   MatSetBlockSize(is->A, bs);
2264:   MatInvertBlockDiagonal(is->A, &lv);
2265:   if (!is->bdiag) PetscMalloc1(bs * mat->rmap->n, &is->bdiag);
2266:   MPI_Type_contiguous(bs, MPIU_SCALAR, &nodeType);
2267:   MPI_Type_commit(&nodeType);
2268:   PetscSFReduceBegin(is->sf, nodeType, lv, is->bdiag, MPI_REPLACE);
2269:   PetscSFReduceEnd(is->sf, nodeType, lv, is->bdiag, MPI_REPLACE);
2270:   MPI_Type_free(&nodeType);
2271:   if (values) *values = is->bdiag;
2272:   return 0;
2273: }

2275: static PetscErrorCode MatISSetUpScatters_Private(Mat A)
2276: {
2277:   Vec             cglobal, rglobal;
2278:   IS              from;
2279:   Mat_IS         *is = (Mat_IS *)A->data;
2280:   PetscScalar     sum;
2281:   const PetscInt *garray;
2282:   PetscInt        nr, rbs, nc, cbs;
2283:   VecType         rtype;

2285:   ISLocalToGlobalMappingGetSize(is->rmapping, &nr);
2286:   ISLocalToGlobalMappingGetBlockSize(is->rmapping, &rbs);
2287:   ISLocalToGlobalMappingGetSize(is->cmapping, &nc);
2288:   ISLocalToGlobalMappingGetBlockSize(is->cmapping, &cbs);
2289:   VecDestroy(&is->x);
2290:   VecDestroy(&is->y);
2291:   VecDestroy(&is->counter);
2292:   VecScatterDestroy(&is->rctx);
2293:   VecScatterDestroy(&is->cctx);
2294:   MatCreateVecs(is->A, &is->x, &is->y);
2295:   VecBindToCPU(is->y, PETSC_TRUE);
2296:   VecGetRootType_Private(is->y, &rtype);
2297:   PetscFree(A->defaultvectype);
2298:   PetscStrallocpy(rtype, &A->defaultvectype);
2299:   MatCreateVecs(A, &cglobal, &rglobal);
2300:   VecBindToCPU(rglobal, PETSC_TRUE);
2301:   ISLocalToGlobalMappingGetBlockIndices(is->rmapping, &garray);
2302:   ISCreateBlock(PetscObjectComm((PetscObject)A), rbs, nr / rbs, garray, PETSC_USE_POINTER, &from);
2303:   VecScatterCreate(rglobal, from, is->y, NULL, &is->rctx);
2304:   ISLocalToGlobalMappingRestoreBlockIndices(is->rmapping, &garray);
2305:   ISDestroy(&from);
2306:   if (is->rmapping != is->cmapping) {
2307:     ISLocalToGlobalMappingGetBlockIndices(is->cmapping, &garray);
2308:     ISCreateBlock(PetscObjectComm((PetscObject)A), cbs, nc / cbs, garray, PETSC_USE_POINTER, &from);
2309:     VecScatterCreate(cglobal, from, is->x, NULL, &is->cctx);
2310:     ISLocalToGlobalMappingRestoreBlockIndices(is->cmapping, &garray);
2311:     ISDestroy(&from);
2312:   } else {
2313:     PetscObjectReference((PetscObject)is->rctx);
2314:     is->cctx = is->rctx;
2315:   }
2316:   VecDestroy(&cglobal);

2318:   /* interface counter vector (local) */
2319:   VecDuplicate(is->y, &is->counter);
2320:   VecBindToCPU(is->counter, PETSC_TRUE);
2321:   VecSet(is->y, 1.);
2322:   VecScatterBegin(is->rctx, is->y, rglobal, ADD_VALUES, SCATTER_REVERSE);
2323:   VecScatterEnd(is->rctx, is->y, rglobal, ADD_VALUES, SCATTER_REVERSE);
2324:   VecScatterBegin(is->rctx, rglobal, is->counter, INSERT_VALUES, SCATTER_FORWARD);
2325:   VecScatterEnd(is->rctx, rglobal, is->counter, INSERT_VALUES, SCATTER_FORWARD);
2326:   VecBindToCPU(is->y, PETSC_FALSE);
2327:   VecBindToCPU(is->counter, PETSC_FALSE);

2329:   /* special functions for block-diagonal matrices */
2330:   VecSum(rglobal, &sum);
2331:   A->ops->invertblockdiagonal = NULL;
2332:   if ((PetscInt)(PetscRealPart(sum)) == A->rmap->N && A->rmap->N == A->cmap->N && is->rmapping == is->cmapping) A->ops->invertblockdiagonal = MatInvertBlockDiagonal_IS;
2333:   VecDestroy(&rglobal);

2335:   /* setup SF for general purpose shared indices based communications */
2336:   MatISSetUpSF_IS(A);
2337:   return 0;
2338: }

2340: static PetscErrorCode MatISFilterL2GMap(Mat A, ISLocalToGlobalMapping map, ISLocalToGlobalMapping *nmap, ISLocalToGlobalMapping *lmap)
2341: {
2342:   IS                         is;
2343:   ISLocalToGlobalMappingType l2gtype;
2344:   const PetscInt            *idxs;
2345:   PetscHSetI                 ht;
2346:   PetscInt                  *nidxs;
2347:   PetscInt                   i, n, bs, c;
2348:   PetscBool                  flg[] = {PETSC_FALSE, PETSC_FALSE};

2350:   ISLocalToGlobalMappingGetSize(map, &n);
2351:   ISLocalToGlobalMappingGetBlockSize(map, &bs);
2352:   ISLocalToGlobalMappingGetBlockIndices(map, &idxs);
2353:   PetscHSetICreate(&ht);
2354:   PetscMalloc1(n / bs, &nidxs);
2355:   for (i = 0, c = 0; i < n / bs; i++) {
2356:     PetscBool missing;
2357:     if (idxs[i] < 0) {
2358:       flg[0] = PETSC_TRUE;
2359:       continue;
2360:     }
2361:     PetscHSetIQueryAdd(ht, idxs[i], &missing);
2362:     if (!missing) flg[1] = PETSC_TRUE;
2363:     else nidxs[c++] = idxs[i];
2364:   }
2365:   PetscHSetIDestroy(&ht);
2366:   MPIU_Allreduce(MPI_IN_PLACE, flg, 2, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A));
2367:   if (!flg[0] && !flg[1]) { /* Entries are all non negative and unique */
2368:     *nmap = NULL;
2369:     *lmap = NULL;
2370:     PetscFree(nidxs);
2371:     ISLocalToGlobalMappingRestoreBlockIndices(map, &idxs);
2372:     return 0;
2373:   }

2375:   /* New l2g map without negative or repeated indices */
2376:   ISCreateBlock(PetscObjectComm((PetscObject)A), bs, c, nidxs, PETSC_USE_POINTER, &is);
2377:   ISLocalToGlobalMappingCreateIS(is, nmap);
2378:   ISDestroy(&is);
2379:   ISLocalToGlobalMappingGetType(map, &l2gtype);
2380:   ISLocalToGlobalMappingSetType(*nmap, l2gtype);

2382:   /* New local l2g map for repeated indices */
2383:   ISGlobalToLocalMappingApplyBlock(*nmap, IS_GTOLM_MASK, n / bs, idxs, NULL, nidxs);
2384:   ISCreateBlock(PETSC_COMM_SELF, bs, n / bs, nidxs, PETSC_USE_POINTER, &is);
2385:   ISLocalToGlobalMappingCreateIS(is, lmap);
2386:   ISDestroy(&is);

2388:   PetscFree(nidxs);
2389:   ISLocalToGlobalMappingRestoreBlockIndices(map, &idxs);
2390:   return 0;
2391: }

2393: static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A, ISLocalToGlobalMapping rmapping, ISLocalToGlobalMapping cmapping)
2394: {
2395:   Mat_IS                *is            = (Mat_IS *)A->data;
2396:   ISLocalToGlobalMapping localrmapping = NULL, localcmapping = NULL;
2397:   PetscBool              cong, freem[]                       = {PETSC_FALSE, PETSC_FALSE};
2398:   PetscInt               nr, rbs, nc, cbs;


2403:   ISLocalToGlobalMappingDestroy(&is->rmapping);
2404:   ISLocalToGlobalMappingDestroy(&is->cmapping);
2405:   PetscLayoutSetUp(A->rmap);
2406:   PetscLayoutSetUp(A->cmap);
2407:   MatHasCongruentLayouts(A, &cong);

2409:   /* If NULL, local space matches global space */
2410:   if (!rmapping) {
2411:     IS is;

2413:     ISCreateStride(PetscObjectComm((PetscObject)A), A->rmap->N, 0, 1, &is);
2414:     ISLocalToGlobalMappingCreateIS(is, &rmapping);
2415:     if (A->rmap->bs > 0) ISLocalToGlobalMappingSetBlockSize(rmapping, A->rmap->bs);
2416:     ISDestroy(&is);
2417:     freem[0] = PETSC_TRUE;
2418:     if (!cmapping && cong && A->rmap->bs == A->cmap->bs) cmapping = rmapping;
2419:   } else if (!is->islocalref) { /* check if the l2g map has negative or repeated entries */
2420:     MatISFilterL2GMap(A, rmapping, &is->rmapping, &localrmapping);
2421:     if (rmapping == cmapping) {
2422:       PetscObjectReference((PetscObject)is->rmapping);
2423:       is->cmapping = is->rmapping;
2424:       PetscObjectReference((PetscObject)localrmapping);
2425:       localcmapping = localrmapping;
2426:     }
2427:   }
2428:   if (!cmapping) {
2429:     IS is;

2431:     ISCreateStride(PetscObjectComm((PetscObject)A), A->cmap->N, 0, 1, &is);
2432:     ISLocalToGlobalMappingCreateIS(is, &cmapping);
2433:     if (A->cmap->bs > 0) ISLocalToGlobalMappingSetBlockSize(cmapping, A->cmap->bs);
2434:     ISDestroy(&is);
2435:     freem[1] = PETSC_TRUE;
2436:   } else if (cmapping != rmapping && !is->islocalref) { /* check if the l2g map has negative or repeated entries */
2437:     MatISFilterL2GMap(A, cmapping, &is->cmapping, &localcmapping);
2438:   }
2439:   if (!is->rmapping) {
2440:     PetscObjectReference((PetscObject)rmapping);
2441:     is->rmapping = rmapping;
2442:   }
2443:   if (!is->cmapping) {
2444:     PetscObjectReference((PetscObject)cmapping);
2445:     is->cmapping = cmapping;
2446:   }

2448:   /* Clean up */
2449:   MatDestroy(&is->A);
2450:   if (is->csf != is->sf) {
2451:     PetscSFDestroy(&is->csf);
2452:     PetscFree2(is->csf_rootdata, is->csf_leafdata);
2453:   } else is->csf = NULL;
2454:   PetscSFDestroy(&is->sf);
2455:   PetscFree2(is->sf_rootdata, is->sf_leafdata);
2456:   PetscFree(is->bdiag);

2458:   /* check if the two mappings are actually the same for square matrices since MATIS has some optimization for this case
2459:      (DOLFIN passes 2 different objects) */
2460:   ISLocalToGlobalMappingGetSize(is->rmapping, &nr);
2461:   ISLocalToGlobalMappingGetBlockSize(is->rmapping, &rbs);
2462:   ISLocalToGlobalMappingGetSize(is->cmapping, &nc);
2463:   ISLocalToGlobalMappingGetBlockSize(is->cmapping, &cbs);
2464:   if (is->rmapping != is->cmapping && cong) {
2465:     PetscBool same = PETSC_FALSE;
2466:     if (nr == nc && cbs == rbs) {
2467:       const PetscInt *idxs1, *idxs2;

2469:       ISLocalToGlobalMappingGetBlockIndices(is->rmapping, &idxs1);
2470:       ISLocalToGlobalMappingGetBlockIndices(is->cmapping, &idxs2);
2471:       PetscArraycmp(idxs1, idxs2, nr / rbs, &same);
2472:       ISLocalToGlobalMappingRestoreBlockIndices(is->rmapping, &idxs1);
2473:       ISLocalToGlobalMappingRestoreBlockIndices(is->cmapping, &idxs2);
2474:     }
2475:     MPIU_Allreduce(MPI_IN_PLACE, &same, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A));
2476:     if (same) {
2477:       ISLocalToGlobalMappingDestroy(&is->cmapping);
2478:       PetscObjectReference((PetscObject)is->rmapping);
2479:       is->cmapping = is->rmapping;
2480:     }
2481:   }
2482:   PetscLayoutSetBlockSize(A->rmap, rbs);
2483:   PetscLayoutSetBlockSize(A->cmap, cbs);
2484:   /* Pass the user defined maps to the layout */
2485:   PetscLayoutSetISLocalToGlobalMapping(A->rmap, rmapping);
2486:   PetscLayoutSetISLocalToGlobalMapping(A->cmap, cmapping);
2487:   if (freem[0]) ISLocalToGlobalMappingDestroy(&rmapping);
2488:   if (freem[1]) ISLocalToGlobalMappingDestroy(&cmapping);

2490:   /* Create the local matrix A */
2491:   MatCreate(PETSC_COMM_SELF, &is->A);
2492:   MatSetType(is->A, is->lmattype);
2493:   MatSetSizes(is->A, nr, nc, nr, nc);
2494:   MatSetBlockSizes(is->A, rbs, cbs);
2495:   MatSetOptionsPrefix(is->A, "is_");
2496:   MatAppendOptionsPrefix(is->A, ((PetscObject)A)->prefix);
2497:   PetscLayoutSetUp(is->A->rmap);
2498:   PetscLayoutSetUp(is->A->cmap);
2499:   MatSetLocalToGlobalMapping(is->A, localrmapping, localcmapping);
2500:   ISLocalToGlobalMappingDestroy(&localrmapping);
2501:   ISLocalToGlobalMappingDestroy(&localcmapping);

2503:   /* setup scatters and local vectors for MatMult */
2504:   if (!is->islocalref) MatISSetUpScatters_Private(A);
2505:   A->preallocated = PETSC_TRUE;
2506:   return 0;
2507: }

2509: static PetscErrorCode MatSetUp_IS(Mat A)
2510: {
2511:   ISLocalToGlobalMapping rmap, cmap;

2513:   MatGetLocalToGlobalMapping(A, &rmap, &cmap);
2514:   if (!rmap && !cmap) MatSetLocalToGlobalMapping(A, NULL, NULL);
2515:   return 0;
2516: }

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

2523:   ISGlobalToLocalMappingApply(is->rmapping, IS_GTOLM_MASK, m, rows, &m, rows_l);
2524:   if (m != n || rows != cols || is->cmapping != is->rmapping) {
2525:     ISGlobalToLocalMappingApply(is->cmapping, IS_GTOLM_MASK, n, cols, &n, cols_l);
2526:     MatSetValues(is->A, m, rows_l, n, cols_l, values, addv);
2527:   } else {
2528:     MatSetValues(is->A, m, rows_l, m, rows_l, values, addv);
2529:   }
2530:   return 0;
2531: }

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

2538:   ISGlobalToLocalMappingApplyBlock(is->rmapping, IS_GTOLM_MASK, m, rows, &m, rows_l);
2539:   if (m != n || rows != cols || is->cmapping != is->rmapping) {
2540:     ISGlobalToLocalMappingApply(is->cmapping, IS_GTOLM_MASK, n, cols, &n, cols_l);
2541:     MatSetValuesBlocked(is->A, m, rows_l, n, cols_l, values, addv);
2542:   } else {
2543:     MatSetValuesBlocked(is->A, m, rows_l, n, rows_l, values, addv);
2544:   }
2545:   return 0;
2546: }

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

2552:   if (is->A->rmap->mapping || is->A->cmap->mapping) {
2553:     MatSetValuesLocal(is->A, m, rows, n, cols, values, addv);
2554:   } else {
2555:     MatSetValues(is->A, m, rows, n, cols, values, addv);
2556:   }
2557:   return 0;
2558: }

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

2564:   if (is->A->rmap->mapping || is->A->cmap->mapping) {
2565:     MatSetValuesBlockedLocal(is->A, m, rows, n, cols, values, addv);
2566:   } else {
2567:     MatSetValuesBlocked(is->A, m, rows, n, cols, values, addv);
2568:   }
2569:   return 0;
2570: }

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

2576:   if (!n) {
2577:     is->pure_neumann = PETSC_TRUE;
2578:   } else {
2579:     PetscInt i;
2580:     is->pure_neumann = PETSC_FALSE;

2582:     if (columns) {
2583:       MatZeroRowsColumns(is->A, n, rows, diag, NULL, NULL);
2584:     } else {
2585:       MatZeroRows(is->A, n, rows, diag, NULL, NULL);
2586:     }
2587:     if (diag != 0.) {
2588:       const PetscScalar *array;
2589:       VecGetArrayRead(is->counter, &array);
2590:       for (i = 0; i < n; i++) MatSetValue(is->A, rows[i], rows[i], diag / (array[rows[i]]), INSERT_VALUES);
2591:       VecRestoreArrayRead(is->counter, &array);
2592:     }
2593:     MatAssemblyBegin(is->A, MAT_FINAL_ASSEMBLY);
2594:     MatAssemblyEnd(is->A, MAT_FINAL_ASSEMBLY);
2595:   }
2596:   return 0;
2597: }

2599: static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b, PetscBool columns)
2600: {
2601:   Mat_IS   *matis = (Mat_IS *)A->data;
2602:   PetscInt  nr, nl, len, i;
2603:   PetscInt *lrows;

2605:   if (PetscUnlikelyDebug(columns || diag != 0. || (x && b))) {
2606:     PetscBool cong;

2608:     PetscLayoutCompare(A->rmap, A->cmap, &cong);
2609:     cong = (PetscBool)(cong && matis->sf == matis->csf);
2613:   }
2614:   /* get locally owned rows */
2615:   PetscLayoutMapLocal(A->rmap, n, rows, &len, &lrows, NULL);
2616:   /* fix right hand side if needed */
2617:   if (x && b) {
2618:     const PetscScalar *xx;
2619:     PetscScalar       *bb;

2621:     VecGetArrayRead(x, &xx);
2622:     VecGetArray(b, &bb);
2623:     for (i = 0; i < len; ++i) bb[lrows[i]] = diag * xx[lrows[i]];
2624:     VecRestoreArrayRead(x, &xx);
2625:     VecRestoreArray(b, &bb);
2626:   }
2627:   /* get rows associated to the local matrices */
2628:   MatGetSize(matis->A, &nl, NULL);
2629:   PetscArrayzero(matis->sf_leafdata, nl);
2630:   PetscArrayzero(matis->sf_rootdata, A->rmap->n);
2631:   for (i = 0; i < len; i++) matis->sf_rootdata[lrows[i]] = 1;
2632:   PetscFree(lrows);
2633:   PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE);
2634:   PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE);
2635:   PetscMalloc1(nl, &lrows);
2636:   for (i = 0, nr = 0; i < nl; i++)
2637:     if (matis->sf_leafdata[i]) lrows[nr++] = i;
2638:   MatISZeroRowsColumnsLocal_Private(A, nr, lrows, diag, columns);
2639:   PetscFree(lrows);
2640:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
2641:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
2642:   return 0;
2643: }

2645: static PetscErrorCode MatZeroRows_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
2646: {
2647:   MatZeroRowsColumns_Private_IS(A, n, rows, diag, x, b, PETSC_FALSE);
2648:   return 0;
2649: }

2651: static PetscErrorCode MatZeroRowsColumns_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
2652: {
2653:   MatZeroRowsColumns_Private_IS(A, n, rows, diag, x, b, PETSC_TRUE);
2654:   return 0;
2655: }

2657: static PetscErrorCode MatAssemblyBegin_IS(Mat A, MatAssemblyType type)
2658: {
2659:   Mat_IS *is = (Mat_IS *)A->data;

2661:   MatAssemblyBegin(is->A, type);
2662:   return 0;
2663: }

2665: static PetscErrorCode MatAssemblyEnd_IS(Mat A, MatAssemblyType type)
2666: {
2667:   Mat_IS   *is = (Mat_IS *)A->data;
2668:   PetscBool lnnz;

2670:   MatAssemblyEnd(is->A, type);
2671:   /* fix for local empty rows/cols */
2672:   if (is->locempty && type == MAT_FINAL_ASSEMBLY) {
2673:     Mat                    newlA;
2674:     ISLocalToGlobalMapping rl2g, cl2g;
2675:     IS                     nzr, nzc;
2676:     PetscInt               nr, nc, nnzr, nnzc;
2677:     PetscBool              lnewl2g, newl2g;

2679:     MatGetSize(is->A, &nr, &nc);
2680:     MatFindNonzeroRowsOrCols_Basic(is->A, PETSC_FALSE, PETSC_SMALL, &nzr);
2681:     if (!nzr) ISCreateStride(PetscObjectComm((PetscObject)is->A), nr, 0, 1, &nzr);
2682:     MatFindNonzeroRowsOrCols_Basic(is->A, PETSC_TRUE, PETSC_SMALL, &nzc);
2683:     if (!nzc) ISCreateStride(PetscObjectComm((PetscObject)is->A), nc, 0, 1, &nzc);
2684:     ISGetSize(nzr, &nnzr);
2685:     ISGetSize(nzc, &nnzc);
2686:     if (nnzr != nr || nnzc != nc) { /* need new global l2g map */
2687:       lnewl2g = PETSC_TRUE;
2688:       MPI_Allreduce(&lnewl2g, &newl2g, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A));

2690:       /* extract valid submatrix */
2691:       MatCreateSubMatrix(is->A, nzr, nzc, MAT_INITIAL_MATRIX, &newlA);
2692:     } else { /* local matrix fully populated */
2693:       lnewl2g = PETSC_FALSE;
2694:       MPI_Allreduce(&lnewl2g, &newl2g, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A));
2695:       PetscObjectReference((PetscObject)is->A);
2696:       newlA = is->A;
2697:     }

2699:     /* attach new global l2g map if needed */
2700:     if (newl2g) {
2701:       IS              zr, zc;
2702:       const PetscInt *ridxs, *cidxs, *zridxs, *zcidxs;
2703:       PetscInt       *nidxs, i;

2705:       ISComplement(nzr, 0, nr, &zr);
2706:       ISComplement(nzc, 0, nc, &zc);
2707:       PetscMalloc1(PetscMax(nr, nc), &nidxs);
2708:       ISLocalToGlobalMappingGetIndices(is->rmapping, &ridxs);
2709:       ISLocalToGlobalMappingGetIndices(is->cmapping, &cidxs);
2710:       ISGetIndices(zr, &zridxs);
2711:       ISGetIndices(zc, &zcidxs);
2712:       ISGetLocalSize(zr, &nnzr);
2713:       ISGetLocalSize(zc, &nnzc);

2715:       PetscArraycpy(nidxs, ridxs, nr);
2716:       for (i = 0; i < nnzr; i++) nidxs[zridxs[i]] = -1;
2717:       ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A), 1, nr, nidxs, PETSC_COPY_VALUES, &rl2g);
2718:       PetscArraycpy(nidxs, cidxs, nc);
2719:       for (i = 0; i < nnzc; i++) nidxs[zcidxs[i]] = -1;
2720:       ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A), 1, nc, nidxs, PETSC_COPY_VALUES, &cl2g);

2722:       ISRestoreIndices(zr, &zridxs);
2723:       ISRestoreIndices(zc, &zcidxs);
2724:       ISLocalToGlobalMappingRestoreIndices(is->rmapping, &ridxs);
2725:       ISLocalToGlobalMappingRestoreIndices(is->cmapping, &cidxs);
2726:       ISDestroy(&nzr);
2727:       ISDestroy(&nzc);
2728:       ISDestroy(&zr);
2729:       ISDestroy(&zc);
2730:       PetscFree(nidxs);
2731:       MatSetLocalToGlobalMapping(A, rl2g, cl2g);
2732:       ISLocalToGlobalMappingDestroy(&rl2g);
2733:       ISLocalToGlobalMappingDestroy(&cl2g);
2734:     }
2735:     MatISSetLocalMat(A, newlA);
2736:     MatDestroy(&newlA);
2737:     ISDestroy(&nzr);
2738:     ISDestroy(&nzc);
2739:     is->locempty = PETSC_FALSE;
2740:   }
2741:   lnnz          = (PetscBool)(is->A->nonzerostate == is->lnnzstate);
2742:   is->lnnzstate = is->A->nonzerostate;
2743:   MPIU_Allreduce(MPI_IN_PLACE, &lnnz, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A));
2744:   if (lnnz) A->nonzerostate++;
2745:   return 0;
2746: }

2748: static PetscErrorCode MatISGetLocalMat_IS(Mat mat, Mat *local)
2749: {
2750:   Mat_IS *is = (Mat_IS *)mat->data;

2752:   *local = is->A;
2753:   return 0;
2754: }

2756: static PetscErrorCode MatISRestoreLocalMat_IS(Mat mat, Mat *local)
2757: {
2758:   *local = NULL;
2759:   return 0;
2760: }

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

2765:   Input Parameter:
2766: .  mat - the matrix

2768:   Output Parameter:
2769: .  local - the local matrix

2771:   Level: advanced

2773:   Notes:
2774:     This can be called if you have precomputed the nonzero structure of the
2775:   matrix and want to provide it to the inner matrix object to improve the performance
2776:   of the `MatSetValues()` operation.

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

2780: .seealso: `MATIS`, `MatISRestoreLocalMat()`
2781: @*/
2782: PetscErrorCode MatISGetLocalMat(Mat mat, Mat *local)
2783: {
2786:   PetscUseMethod(mat, "MatISGetLocalMat_C", (Mat, Mat *), (mat, local));
2787:   return 0;
2788: }

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

2793:   Input Parameter:
2794: .  mat - the matrix

2796:   Output Parameter:
2797: .  local - the local matrix

2799:   Level: advanced

2801: .seealso: `MATIS`, `MatISGetLocalMat()`
2802: @*/
2803: PetscErrorCode MatISRestoreLocalMat(Mat mat, Mat *local)
2804: {
2807:   PetscUseMethod(mat, "MatISRestoreLocalMat_C", (Mat, Mat *), (mat, local));
2808:   return 0;
2809: }

2811: static PetscErrorCode MatISSetLocalMatType_IS(Mat mat, MatType mtype)
2812: {
2813:   Mat_IS *is = (Mat_IS *)mat->data;

2815:   if (is->A) MatSetType(is->A, mtype);
2816:   PetscFree(is->lmattype);
2817:   PetscStrallocpy(mtype, &is->lmattype);
2818:   return 0;
2819: }

2821: /*@
2822:     MatISSetLocalMatType - Specifies the type of local matrix inside the `MATIS`

2824:   Input Parameters:
2825: +  mat - the matrix
2826: -  mtype - the local matrix type

2828:   Level: advanced

2830: .seealso: `MATIS`, `MatSetType()`, `MatType`
2831: @*/
2832: PetscErrorCode MatISSetLocalMatType(Mat mat, MatType mtype)
2833: {
2835:   PetscUseMethod(mat, "MatISSetLocalMatType_C", (Mat, MatType), (mat, mtype));
2836:   return 0;
2837: }

2839: static PetscErrorCode MatISSetLocalMat_IS(Mat mat, Mat local)
2840: {
2841:   Mat_IS   *is = (Mat_IS *)mat->data;
2842:   PetscInt  nrows, ncols, orows, ocols;
2843:   MatType   mtype, otype;
2844:   PetscBool sametype = PETSC_TRUE;

2846:   if (is->A && !is->islocalref) {
2847:     MatGetSize(is->A, &orows, &ocols);
2848:     MatGetSize(local, &nrows, &ncols);
2850:     MatGetType(local, &mtype);
2851:     MatGetType(is->A, &otype);
2852:     PetscStrcmp(mtype, otype, &sametype);
2853:   }
2854:   PetscObjectReference((PetscObject)local);
2855:   MatDestroy(&is->A);
2856:   is->A = local;
2857:   MatGetType(is->A, &mtype);
2858:   MatISSetLocalMatType(mat, mtype);
2859:   if (!sametype && !is->islocalref) MatISSetUpScatters_Private(mat);
2860:   return 0;
2861: }

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

2866:   Collective

2868:   Input Parameters:
2869: +  mat - the matrix
2870: -  local - the local matrix

2872:   Level: advanced

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

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

2880: .seealso: `MATIS`, `MatISSetLocalMatType`, `MatISGetLocalMat()`
2881: @*/
2882: PetscErrorCode MatISSetLocalMat(Mat mat, Mat local)
2883: {
2886:   PetscUseMethod(mat, "MatISSetLocalMat_C", (Mat, Mat), (mat, local));
2887:   return 0;
2888: }

2890: static PetscErrorCode MatZeroEntries_IS(Mat A)
2891: {
2892:   Mat_IS *a = (Mat_IS *)A->data;

2894:   MatZeroEntries(a->A);
2895:   return 0;
2896: }

2898: static PetscErrorCode MatScale_IS(Mat A, PetscScalar a)
2899: {
2900:   Mat_IS *is = (Mat_IS *)A->data;

2902:   MatScale(is->A, a);
2903:   return 0;
2904: }

2906: static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v)
2907: {
2908:   Mat_IS *is = (Mat_IS *)A->data;

2910:   /* get diagonal of the local matrix */
2911:   MatGetDiagonal(is->A, is->y);

2913:   /* scatter diagonal back into global vector */
2914:   VecSet(v, 0);
2915:   VecScatterBegin(is->rctx, is->y, v, ADD_VALUES, SCATTER_REVERSE);
2916:   VecScatterEnd(is->rctx, is->y, v, ADD_VALUES, SCATTER_REVERSE);
2917:   return 0;
2918: }

2920: static PetscErrorCode MatSetOption_IS(Mat A, MatOption op, PetscBool flg)
2921: {
2922:   Mat_IS *a = (Mat_IS *)A->data;

2924:   MatSetOption(a->A, op, flg);
2925:   return 0;
2926: }

2928: static PetscErrorCode MatAXPY_IS(Mat Y, PetscScalar a, Mat X, MatStructure str)
2929: {
2930:   Mat_IS *y = (Mat_IS *)Y->data;
2931:   Mat_IS *x;

2933:   if (PetscDefined(USE_DEBUG)) {
2934:     PetscBool ismatis;
2935:     PetscObjectTypeCompare((PetscObject)X, MATIS, &ismatis);
2937:   }
2938:   x = (Mat_IS *)X->data;
2939:   MatAXPY(y->A, a, x->A, str);
2940:   return 0;
2941: }

2943: static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A, IS row, IS col, Mat *submat)
2944: {
2945:   Mat                    lA;
2946:   Mat_IS                *matis = (Mat_IS *)(A->data);
2947:   ISLocalToGlobalMapping rl2g, cl2g;
2948:   IS                     is;
2949:   const PetscInt        *rg, *rl;
2950:   PetscInt               nrg;
2951:   PetscInt               N, M, nrl, i, *idxs;

2953:   ISLocalToGlobalMappingGetIndices(A->rmap->mapping, &rg);
2954:   ISGetLocalSize(row, &nrl);
2955:   ISGetIndices(row, &rl);
2956:   ISLocalToGlobalMappingGetSize(A->rmap->mapping, &nrg);
2957:   if (PetscDefined(USE_DEBUG)) {
2959:   }
2960:   PetscMalloc1(nrg, &idxs);
2961:   /* map from [0,nrl) to row */
2962:   for (i = 0; i < nrl; i++) idxs[i] = rl[i];
2963:   for (i = nrl; i < nrg; i++) idxs[i] = -1;
2964:   ISRestoreIndices(row, &rl);
2965:   ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping, &rg);
2966:   ISCreateGeneral(PetscObjectComm((PetscObject)A), nrg, idxs, PETSC_OWN_POINTER, &is);
2967:   ISLocalToGlobalMappingCreateIS(is, &rl2g);
2968:   ISDestroy(&is);
2969:   /* compute new l2g map for columns */
2970:   if (col != row || matis->rmapping != matis->cmapping || matis->A->rmap->mapping != matis->A->cmap->mapping) {
2971:     const PetscInt *cg, *cl;
2972:     PetscInt        ncg;
2973:     PetscInt        ncl;

2975:     ISLocalToGlobalMappingGetIndices(A->cmap->mapping, &cg);
2976:     ISGetLocalSize(col, &ncl);
2977:     ISGetIndices(col, &cl);
2978:     ISLocalToGlobalMappingGetSize(A->cmap->mapping, &ncg);
2979:     if (PetscDefined(USE_DEBUG)) {
2981:     }
2982:     PetscMalloc1(ncg, &idxs);
2983:     /* map from [0,ncl) to col */
2984:     for (i = 0; i < ncl; i++) idxs[i] = cl[i];
2985:     for (i = ncl; i < ncg; i++) idxs[i] = -1;
2986:     ISRestoreIndices(col, &cl);
2987:     ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping, &cg);
2988:     ISCreateGeneral(PetscObjectComm((PetscObject)A), ncg, idxs, PETSC_OWN_POINTER, &is);
2989:     ISLocalToGlobalMappingCreateIS(is, &cl2g);
2990:     ISDestroy(&is);
2991:   } else {
2992:     PetscObjectReference((PetscObject)rl2g);
2993:     cl2g = rl2g;
2994:   }
2995:   /* create the MATIS submatrix */
2996:   MatGetSize(A, &M, &N);
2997:   MatCreate(PetscObjectComm((PetscObject)A), submat);
2998:   MatSetSizes(*submat, PETSC_DECIDE, PETSC_DECIDE, M, N);
2999:   MatSetType(*submat, MATIS);
3000:   matis             = (Mat_IS *)((*submat)->data);
3001:   matis->islocalref = PETSC_TRUE;
3002:   MatSetLocalToGlobalMapping(*submat, rl2g, cl2g);
3003:   MatISGetLocalMat(A, &lA);
3004:   MatISSetLocalMat(*submat, lA);
3005:   MatSetUp(*submat);
3006:   MatAssemblyBegin(*submat, MAT_FINAL_ASSEMBLY);
3007:   MatAssemblyEnd(*submat, MAT_FINAL_ASSEMBLY);
3008:   ISLocalToGlobalMappingDestroy(&rl2g);
3009:   ISLocalToGlobalMappingDestroy(&cl2g);

3011:   /* remove unsupported ops */
3012:   PetscMemzero((*submat)->ops, sizeof(struct _MatOps));
3013:   (*submat)->ops->destroy               = MatDestroy_IS;
3014:   (*submat)->ops->setvalueslocal        = MatSetValuesLocal_SubMat_IS;
3015:   (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS;
3016:   (*submat)->ops->assemblybegin         = MatAssemblyBegin_IS;
3017:   (*submat)->ops->assemblyend           = MatAssemblyEnd_IS;
3018:   return 0;
3019: }

3021: static PetscErrorCode MatSetFromOptions_IS(Mat A, PetscOptionItems *PetscOptionsObject)
3022: {
3023:   Mat_IS   *a = (Mat_IS *)A->data;
3024:   char      type[256];
3025:   PetscBool flg;

3027:   PetscOptionsHeadBegin(PetscOptionsObject, "MATIS options");
3028:   PetscOptionsBool("-matis_keepassembled", "Store an assembled version if needed", "MatISKeepAssembled", a->keepassembled, &a->keepassembled, NULL);
3029:   PetscOptionsBool("-matis_fixempty", "Fix local matrices in case of empty local rows/columns", "MatISFixLocalEmpty", a->locempty, &a->locempty, NULL);
3030:   PetscOptionsBool("-matis_storel2l", "Store local-to-local matrices generated from PtAP operations", "MatISStoreL2L", a->storel2l, &a->storel2l, NULL);
3031:   PetscOptionsFList("-matis_localmat_type", "Matrix type", "MatISSetLocalMatType", MatList, a->lmattype, type, 256, &flg);
3032:   if (flg) MatISSetLocalMatType(A, type);
3033:   if (a->A) MatSetFromOptions(a->A);
3034:   PetscOptionsHeadEnd();
3035:   return 0;
3036: }

3038: /*@
3039:     MatCreateIS - Creates a "process" unassembled matrix, `MATIS`, assembled on each
3040:        process but not across processes.

3042:    Input Parameters:
3043: +     comm    - MPI communicator that will share the matrix
3044: .     bs      - block size of the matrix
3045: .     m,n,M,N - local and/or global sizes of the left and right vector used in matrix vector products
3046: .     rmap    - local to global map for rows
3047: -     cmap    - local to global map for cols

3049:    Output Parameter:
3050: .    A - the resulting matrix

3052:    Level: advanced

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

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

3060: .seealso: `MATIS`, `MatSetLocalToGlobalMapping()`
3061: @*/
3062: PetscErrorCode MatCreateIS(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt M, PetscInt N, ISLocalToGlobalMapping rmap, ISLocalToGlobalMapping cmap, Mat *A)
3063: {
3064:   MatCreate(comm, A);
3065:   MatSetSizes(*A, m, n, M, N);
3066:   if (bs > 0) MatSetBlockSize(*A, bs);
3067:   MatSetType(*A, MATIS);
3068:   MatSetLocalToGlobalMapping(*A, rmap, cmap);
3069:   return 0;
3070: }

3072: static PetscErrorCode MatHasOperation_IS(Mat A, MatOperation op, PetscBool *has)
3073: {
3074:   Mat_IS      *a              = (Mat_IS *)A->data;
3075:   MatOperation tobefiltered[] = {MATOP_MULT_ADD, MATOP_MULT_TRANSPOSE_ADD, MATOP_GET_DIAGONAL_BLOCK, MATOP_INCREASE_OVERLAP};

3077:   *has = PETSC_FALSE;
3078:   if (!((void **)A->ops)[op] || !a->A) return 0;
3079:   *has = PETSC_TRUE;
3080:   for (PetscInt i = 0; i < (PetscInt)PETSC_STATIC_ARRAY_LENGTH(tobefiltered); i++)
3081:     if (op == tobefiltered[i]) return 0;
3082:   MatHasOperation(a->A, op, has);
3083:   return 0;
3084: }

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

3090:   MatSetValuesCOO(a->A, v, imode);
3091:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
3092:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
3093:   return 0;
3094: }

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

3101:   if (a->A->rmap->mapping || a->A->cmap->mapping) {
3102:     MatSetPreallocationCOOLocal(a->A, ncoo, coo_i, coo_j);
3103:   } else {
3104:     MatSetPreallocationCOO(a->A, ncoo, coo_i, coo_j);
3105:   }
3106:   PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", MatSetValuesCOO_IS);
3107:   A->preallocated = PETSC_TRUE;
3108:   return 0;
3109: }

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

3117:   ISGlobalToLocalMappingApply(a->rmapping, IS_GTOLM_MASK, ncoo, coo_i, NULL, coo_i);
3118:   ISGlobalToLocalMappingApply(a->cmapping, IS_GTOLM_MASK, ncoo, coo_j, NULL, coo_j);
3119:   MatSetPreallocationCOO(a->A, ncoo, coo_i, coo_j);
3120:   PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", MatSetValuesCOO_IS);
3121:   A->preallocated = PETSC_TRUE;
3122:   return 0;
3123: }

3125: static PetscErrorCode MatISGetAssembled_Private(Mat A, Mat *tA)
3126: {
3127:   Mat_IS          *a = (Mat_IS *)A->data;
3128:   PetscObjectState Astate, aAstate       = PETSC_MIN_INT;
3129:   PetscObjectState Annzstate, aAnnzstate = PETSC_MIN_INT;

3131:   PetscObjectStateGet((PetscObject)A, &Astate);
3132:   Annzstate = A->nonzerostate;
3133:   if (a->assembledA) {
3134:     PetscObjectStateGet((PetscObject)a->assembledA, &aAstate);
3135:     aAnnzstate = a->assembledA->nonzerostate;
3136:   }
3137:   if (aAnnzstate != Annzstate) MatDestroy(&a->assembledA);
3138:   if (Astate != aAstate || !a->assembledA) {
3139:     MatType     aAtype;
3140:     PetscMPIInt size;
3141:     PetscInt    rbs, cbs, bs;

3143:     /* the assembled form is used as temporary storage for parallel operations
3144:        like createsubmatrices and the like, do not waste device memory */
3145:     MPI_Comm_size(PetscObjectComm((PetscObject)A), &size);
3146:     ISLocalToGlobalMappingGetBlockSize(a->cmapping, &cbs);
3147:     ISLocalToGlobalMappingGetBlockSize(a->rmapping, &rbs);
3148:     bs = rbs == cbs ? rbs : 1;
3149:     if (a->assembledA) MatGetType(a->assembledA, &aAtype);
3150:     else if (size > 1) aAtype = bs > 1 ? MATMPIBAIJ : MATMPIAIJ;
3151:     else aAtype = bs > 1 ? MATSEQBAIJ : MATSEQAIJ;

3153:     MatConvert(A, aAtype, a->assembledA ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX, &a->assembledA);
3154:     PetscObjectStateSet((PetscObject)a->assembledA, Astate);
3155:     a->assembledA->nonzerostate = Annzstate;
3156:   }
3157:   PetscObjectReference((PetscObject)a->assembledA);
3158:   *tA = a->assembledA;
3159:   if (!a->keepassembled) MatDestroy(&a->assembledA);
3160:   return 0;
3161: }

3163: static PetscErrorCode MatISRestoreAssembled_Private(Mat A, Mat *tA)
3164: {
3165:   MatDestroy(tA);
3166:   *tA = NULL;
3167:   return 0;
3168: }

3170: static PetscErrorCode MatGetDiagonalBlock_IS(Mat A, Mat *dA)
3171: {
3172:   Mat_IS          *a = (Mat_IS *)A->data;
3173:   PetscObjectState Astate, dAstate = PETSC_MIN_INT;

3175:   PetscObjectStateGet((PetscObject)A, &Astate);
3176:   if (a->dA) PetscObjectStateGet((PetscObject)a->dA, &dAstate);
3177:   if (Astate != dAstate) {
3178:     Mat     tA;
3179:     MatType ltype;

3181:     MatDestroy(&a->dA);
3182:     MatISGetAssembled_Private(A, &tA);
3183:     MatGetDiagonalBlock(tA, &a->dA);
3184:     MatPropagateSymmetryOptions(tA, a->dA);
3185:     MatGetType(a->A, &ltype);
3186:     MatConvert(a->dA, ltype, MAT_INPLACE_MATRIX, &a->dA);
3187:     PetscObjectReference((PetscObject)a->dA);
3188:     MatISRestoreAssembled_Private(A, &tA);
3189:     PetscObjectStateSet((PetscObject)a->dA, Astate);
3190:   }
3191:   *dA = a->dA;
3192:   return 0;
3193: }

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

3199:   MatISGetAssembled_Private(A, &tA);
3200:   MatCreateSubMatrices(tA, n, irow, icol, reuse, submat);
3201:   /* MatCreateSubMatrices_MPIAIJ is a mess at the moment */
3202: #if 0
3203:   {
3204:     Mat_IS    *a = (Mat_IS*)A->data;
3205:     MatType   ltype;
3206:     VecType   vtype;
3207:     char      *flg;

3209:     MatGetType(a->A,&ltype);
3210:     MatGetVecType(a->A,&vtype);
3211:     PetscStrstr(vtype,"cuda",&flg);
3212:     if (!flg) PetscStrstr(vtype,"hip",&flg);
3213:     if (!flg) PetscStrstr(vtype,"kokkos",&flg);
3214:     if (flg) {
3215:       for (PetscInt i = 0; i < n; i++) {
3216:         Mat sA = (*submat)[i];

3218:         MatConvert(sA,ltype,MAT_INPLACE_MATRIX,&sA);
3219:         (*submat)[i] = sA;
3220:       }
3221:     }
3222:   }
3223: #endif
3224:   MatISRestoreAssembled_Private(A, &tA);
3225:   return 0;
3226: }

3228: static PetscErrorCode MatIncreaseOverlap_IS(Mat A, PetscInt n, IS is[], PetscInt ov)
3229: {
3230:   Mat tA;

3232:   MatISGetAssembled_Private(A, &tA);
3233:   MatIncreaseOverlap(tA, n, is, ov);
3234:   MatISRestoreAssembled_Private(A, &tA);
3235:   return 0;
3236: }

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

3241:    Not Collective

3243:    Input Parameter:
3244: .  A - the matrix

3246:    Output Parameters:
3247: +  rmapping - row mapping
3248: -  cmapping - column mapping

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

3253:    Level: advanced

3255: .seealso: `MatIS`, `MatSetLocalToGlobalMapping()`
3256: @*/
3257: PetscErrorCode MatISGetLocalToGlobalMapping(Mat A, ISLocalToGlobalMapping *rmapping, ISLocalToGlobalMapping *cmapping)
3258: {
3263:   PetscUseMethod(A, "MatISGetLocalToGlobalMapping_C", (Mat, ISLocalToGlobalMapping *, ISLocalToGlobalMapping *), (A, rmapping, cmapping));
3264:   return 0;
3265: }

3267: static PetscErrorCode MatISGetLocalToGlobalMapping_IS(Mat A, ISLocalToGlobalMapping *r, ISLocalToGlobalMapping *c)
3268: {
3269:   Mat_IS *a = (Mat_IS *)A->data;

3271:   if (r) *r = a->rmapping;
3272:   if (c) *c = a->cmapping;
3273:   return 0;
3274: }

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

3281:    Options Database Keys:
3282: + -mat_type is - sets the matrix type to "is" during a call to MatSetFromOptions()
3283: . -matis_fixempty - Fixes local matrices in case of empty local rows/columns.
3284: - -matis_storel2l - stores the local-to-local operators generated by the Galerkin process of MatPtAP().

3286:    Notes:
3287:    Options prefix for the inner matrix are given by -is_mat_xxx

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

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

3294:   Level: advanced

3296: .seealso: `MATIS`, `Mat`, `MatISGetLocalMat()`, `MatSetLocalToGlobalMapping()`, `MatISSetPreallocation()`, `MatCreateIS()`, `PCBDDC`, `KSPFETIDP`
3297: M*/
3298: PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A)
3299: {
3300:   Mat_IS *a;

3302:   PetscNew(&a);
3303:   PetscStrallocpy(MATAIJ, &a->lmattype);
3304:   A->data = (void *)a;

3306:   /* matrix ops */
3307:   PetscMemzero(A->ops, sizeof(struct _MatOps));
3308:   A->ops->mult                    = MatMult_IS;
3309:   A->ops->multadd                 = MatMultAdd_IS;
3310:   A->ops->multtranspose           = MatMultTranspose_IS;
3311:   A->ops->multtransposeadd        = MatMultTransposeAdd_IS;
3312:   A->ops->destroy                 = MatDestroy_IS;
3313:   A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS;
3314:   A->ops->setvalues               = MatSetValues_IS;
3315:   A->ops->setvaluesblocked        = MatSetValuesBlocked_IS;
3316:   A->ops->setvalueslocal          = MatSetValuesLocal_IS;
3317:   A->ops->setvaluesblockedlocal   = MatSetValuesBlockedLocal_IS;
3318:   A->ops->zerorows                = MatZeroRows_IS;
3319:   A->ops->zerorowscolumns         = MatZeroRowsColumns_IS;
3320:   A->ops->assemblybegin           = MatAssemblyBegin_IS;
3321:   A->ops->assemblyend             = MatAssemblyEnd_IS;
3322:   A->ops->view                    = MatView_IS;
3323:   A->ops->zeroentries             = MatZeroEntries_IS;
3324:   A->ops->scale                   = MatScale_IS;
3325:   A->ops->getdiagonal             = MatGetDiagonal_IS;
3326:   A->ops->setoption               = MatSetOption_IS;
3327:   A->ops->ishermitian             = MatIsHermitian_IS;
3328:   A->ops->issymmetric             = MatIsSymmetric_IS;
3329:   A->ops->isstructurallysymmetric = MatIsStructurallySymmetric_IS;
3330:   A->ops->duplicate               = MatDuplicate_IS;
3331:   A->ops->missingdiagonal         = MatMissingDiagonal_IS;
3332:   A->ops->copy                    = MatCopy_IS;
3333:   A->ops->getlocalsubmatrix       = MatGetLocalSubMatrix_IS;
3334:   A->ops->createsubmatrix         = MatCreateSubMatrix_IS;
3335:   A->ops->axpy                    = MatAXPY_IS;
3336:   A->ops->diagonalset             = MatDiagonalSet_IS;
3337:   A->ops->shift                   = MatShift_IS;
3338:   A->ops->transpose               = MatTranspose_IS;
3339:   A->ops->getinfo                 = MatGetInfo_IS;
3340:   A->ops->diagonalscale           = MatDiagonalScale_IS;
3341:   A->ops->setfromoptions          = MatSetFromOptions_IS;
3342:   A->ops->setup                   = MatSetUp_IS;
3343:   A->ops->hasoperation            = MatHasOperation_IS;
3344:   A->ops->getdiagonalblock        = MatGetDiagonalBlock_IS;
3345:   A->ops->createsubmatrices       = MatCreateSubMatrices_IS;
3346:   A->ops->increaseoverlap         = MatIncreaseOverlap_IS;

3348:   /* special MATIS functions */
3349:   PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMatType_C", MatISSetLocalMatType_IS);
3350:   PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalMat_C", MatISGetLocalMat_IS);
3351:   PetscObjectComposeFunction((PetscObject)A, "MatISRestoreLocalMat_C", MatISRestoreLocalMat_IS);
3352:   PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMat_C", MatISSetLocalMat_IS);
3353:   PetscObjectComposeFunction((PetscObject)A, "MatISGetMPIXAIJ_C", MatConvert_IS_XAIJ);
3354:   PetscObjectComposeFunction((PetscObject)A, "MatISSetPreallocation_C", MatISSetPreallocation_IS);
3355:   PetscObjectComposeFunction((PetscObject)A, "MatISStoreL2L_C", MatISStoreL2L_IS);
3356:   PetscObjectComposeFunction((PetscObject)A, "MatISFixLocalEmpty_C", MatISFixLocalEmpty_IS);
3357:   PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalToGlobalMapping_C", MatISGetLocalToGlobalMapping_IS);
3358:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpiaij_C", MatConvert_IS_XAIJ);
3359:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpibaij_C", MatConvert_IS_XAIJ);
3360:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpisbaij_C", MatConvert_IS_XAIJ);
3361:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqaij_C", MatConvert_IS_XAIJ);
3362:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqbaij_C", MatConvert_IS_XAIJ);
3363:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqsbaij_C", MatConvert_IS_XAIJ);
3364:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_aij_C", MatConvert_IS_XAIJ);
3365:   PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOOLocal_C", MatSetPreallocationCOOLocal_IS);
3366:   PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOO_C", MatSetPreallocationCOO_IS);
3367:   PetscObjectChangeTypeName((PetscObject)A, MATIS);
3368:   return 0;
3369: }