Actual source code: matis.c

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

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

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

     16: #define MATIS_MAX_ENTRIES_INSERTION 2048

     18: /* copied from src/mat/impls/localref/mlocalref.c */
     19: #define IndexSpaceGet(buf, nrow, ncol, irowm, icolm) \
     20:   do { \
     21:     if (nrow + ncol > (PetscInt)PETSC_STATIC_ARRAY_LENGTH(buf)) { \
     22:       PetscCall(PetscMalloc2(nrow, &irowm, ncol, &icolm)); \
     23:     } else { \
     24:       irowm = &buf[0]; \
     25:       icolm = &buf[nrow]; \
     26:     } \
     27:   } while (0)

     29: #define IndexSpaceRestore(buf, nrow, ncol, irowm, icolm) \
     30:   do { \
     31:     if (nrow + ncol > (PetscInt)PETSC_STATIC_ARRAY_LENGTH(buf)) PetscCall(PetscFree2(irowm, icolm)); \
     32:   } while (0)

     34: static void BlockIndicesExpand(PetscInt n, const PetscInt idx[], PetscInt bs, PetscInt idxm[])
     35: {
     36:   for (PetscInt i = 0; i < n; i++) {
     37:     for (PetscInt j = 0; j < bs; j++) idxm[i * bs + j] = idx[i] * bs + j;
     38:   }
     39: }

     41: static PetscErrorCode MatSetValuesLocal_IS(Mat, PetscInt, const PetscInt *, PetscInt, const PetscInt *, const PetscScalar *, InsertMode);
     42: static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat, PetscInt, const PetscInt *, PetscInt, const PetscInt *, const PetscScalar *, InsertMode);
     43: static PetscErrorCode MatISSetUpScatters_Private(Mat);

     45: static PetscErrorCode MatISContainerDestroyPtAP_Private(PetscCtxRt ptr)
     46: {
     47:   MatISPtAP ptap = *(MatISPtAP *)ptr;

     49:   PetscFunctionBegin;
     50:   PetscCall(MatDestroySubMatrices(ptap->ris1 ? 2 : 1, &ptap->lP));
     51:   PetscCall(ISDestroy(&ptap->cis0));
     52:   PetscCall(ISDestroy(&ptap->cis1));
     53:   PetscCall(ISDestroy(&ptap->ris0));
     54:   PetscCall(ISDestroy(&ptap->ris1));
     55:   PetscCall(PetscFree(ptap));
     56:   PetscFunctionReturn(PETSC_SUCCESS);
     57: }

     59: static PetscErrorCode MatPtAPNumeric_IS_XAIJ(Mat A, Mat P, Mat C)
     60: {
     61:   MatISPtAP      ptap;
     62:   Mat_IS        *matis = (Mat_IS *)A->data;
     63:   Mat            lA, lC;
     64:   MatReuse       reuse;
     65:   IS             ris[2], cis[2];
     66:   PetscContainer c;
     67:   PetscInt       n;

     69:   PetscFunctionBegin;
     70:   PetscCall(PetscObjectQuery((PetscObject)C, "_MatIS_PtAP", (PetscObject *)&c));
     71:   PetscCheck(c, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Missing PtAP information");
     72:   PetscCall(PetscContainerGetPointer(c, &ptap));
     73:   ris[0] = ptap->ris0;
     74:   ris[1] = ptap->ris1;
     75:   cis[0] = ptap->cis0;
     76:   cis[1] = ptap->cis1;
     77:   n      = ptap->ris1 ? 2 : 1;
     78:   reuse  = ptap->lP ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX;
     79:   PetscCall(MatCreateSubMatrices(P, n, ris, cis, reuse, &ptap->lP));

     81:   PetscCall(MatISGetLocalMat(A, &lA));
     82:   PetscCall(MatISGetLocalMat(C, &lC));
     83:   if (ptap->ris1) { /* unsymmetric A mapping */
     84:     Mat lPt;

     86:     PetscCall(MatTranspose(ptap->lP[1], MAT_INITIAL_MATRIX, &lPt));
     87:     PetscCall(MatMatMatMult(lPt, lA, ptap->lP[0], reuse, ptap->fill, &lC));
     88:     if (matis->storel2l) PetscCall(PetscObjectCompose((PetscObject)A, "_MatIS_PtAP_l2l", (PetscObject)lPt));
     89:     PetscCall(MatDestroy(&lPt));
     90:   } else {
     91:     PetscCall(MatPtAP(lA, ptap->lP[0], reuse, ptap->fill, &lC));
     92:     if (matis->storel2l) PetscCall(PetscObjectCompose((PetscObject)C, "_MatIS_PtAP_l2l", (PetscObject)ptap->lP[0]));
     93:   }
     94:   if (reuse == MAT_INITIAL_MATRIX) {
     95:     PetscCall(MatISSetLocalMat(C, lC));
     96:     PetscCall(MatDestroy(&lC));
     97:   }
     98:   PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
     99:   PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
    100:   PetscFunctionReturn(PETSC_SUCCESS);
    101: }

    103: static PetscErrorCode MatGetNonzeroColumnsLocal_Private(Mat PT, IS *cis)
    104: {
    105:   Mat             Po, Pd;
    106:   IS              zd, zo;
    107:   const PetscInt *garray;
    108:   PetscInt       *aux, i, bs;
    109:   PetscInt        dc, stc, oc, ctd, cto;
    110:   PetscBool       ismpiaij, ismpibaij, isseqaij, isseqbaij;
    111:   MPI_Comm        comm;

    113:   PetscFunctionBegin;
    115:   PetscAssertPointer(cis, 2);
    116:   PetscCall(PetscObjectGetComm((PetscObject)PT, &comm));
    117:   bs = 1;
    118:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)PT, MATMPIAIJ, &ismpiaij));
    119:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)PT, MATMPIBAIJ, &ismpibaij));
    120:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)PT, MATSEQAIJ, &isseqaij));
    121:   PetscCall(PetscObjectTypeCompare((PetscObject)PT, MATSEQBAIJ, &isseqbaij));
    122:   if (isseqaij || isseqbaij) {
    123:     Pd     = PT;
    124:     Po     = NULL;
    125:     garray = NULL;
    126:   } else if (ismpiaij) {
    127:     PetscCall(MatMPIAIJGetSeqAIJ(PT, &Pd, &Po, &garray));
    128:   } else if (ismpibaij) {
    129:     PetscCall(MatMPIBAIJGetSeqBAIJ(PT, &Pd, &Po, &garray));
    130:     PetscCall(MatGetBlockSize(PT, &bs));
    131:   } else SETERRQ(comm, PETSC_ERR_SUP, "Not for matrix type %s", ((PetscObject)PT)->type_name);

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

    140:   PetscCall(MatGetLocalSize(PT, NULL, &dc));
    141:   PetscCall(MatGetOwnershipRangeColumn(PT, &stc, NULL));
    142:   if (Po) PetscCall(MatGetLocalSize(Po, NULL, &oc));
    143:   else oc = 0;
    144:   PetscCall(PetscMalloc1((dc + oc) / bs, &aux));
    145:   if (zd) {
    146:     const PetscInt *idxs;
    147:     PetscInt        nz;

    149:     /* this will throw an error if bs is not valid */
    150:     PetscCall(ISSetBlockSize(zd, bs));
    151:     PetscCall(ISGetLocalSize(zd, &nz));
    152:     PetscCall(ISGetIndices(zd, &idxs));
    153:     ctd = nz / bs;
    154:     for (i = 0; i < ctd; i++) aux[i] = (idxs[bs * i] + stc) / bs;
    155:     PetscCall(ISRestoreIndices(zd, &idxs));
    156:   } else {
    157:     ctd = dc / bs;
    158:     for (i = 0; i < ctd; i++) aux[i] = i + stc / bs;
    159:   }
    160:   if (zo) {
    161:     const PetscInt *idxs;
    162:     PetscInt        nz;

    164:     /* this will throw an error if bs is not valid */
    165:     PetscCall(ISSetBlockSize(zo, bs));
    166:     PetscCall(ISGetLocalSize(zo, &nz));
    167:     PetscCall(ISGetIndices(zo, &idxs));
    168:     cto = nz / bs;
    169:     for (i = 0; i < cto; i++) aux[i + ctd] = garray[idxs[bs * i] / bs];
    170:     PetscCall(ISRestoreIndices(zo, &idxs));
    171:   } else {
    172:     cto = oc / bs;
    173:     for (i = 0; i < cto; i++) aux[i + ctd] = garray[i];
    174:   }
    175:   PetscCall(ISCreateBlock(comm, bs, ctd + cto, aux, PETSC_OWN_POINTER, cis));
    176:   PetscCall(ISDestroy(&zd));
    177:   PetscCall(ISDestroy(&zo));
    178:   PetscFunctionReturn(PETSC_SUCCESS);
    179: }

    181: static PetscErrorCode MatPtAPSymbolic_IS_XAIJ(Mat A, Mat P, PetscReal fill, Mat C)
    182: {
    183:   Mat                    PT, lA;
    184:   MatISPtAP              ptap;
    185:   ISLocalToGlobalMapping Crl2g, Ccl2g, rl2g, cl2g;
    186:   PetscContainer         c;
    187:   MatType                lmtype;
    188:   const PetscInt        *garray;
    189:   PetscInt               ibs, N, dc;
    190:   MPI_Comm               comm;

    192:   PetscFunctionBegin;
    193:   PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
    194:   PetscCall(MatSetType(C, MATIS));
    195:   PetscCall(MatISGetLocalMat(A, &lA));
    196:   PetscCall(MatGetType(lA, &lmtype));
    197:   PetscCall(MatISSetLocalMatType(C, lmtype));
    198:   PetscCall(MatGetSize(P, NULL, &N));
    199:   PetscCall(MatGetLocalSize(P, NULL, &dc));
    200:   PetscCall(MatSetSizes(C, dc, dc, N, N));
    201:   /* Not sure about this
    202:   PetscCall(MatGetBlockSizes(P,NULL,&ibs));
    203:   PetscCall(MatSetBlockSize(*C,ibs));
    204: */

    206:   PetscCall(PetscNew(&ptap));
    207:   PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &c));
    208:   PetscCall(PetscContainerSetPointer(c, ptap));
    209:   PetscCall(PetscContainerSetCtxDestroy(c, MatISContainerDestroyPtAP_Private));
    210:   PetscCall(PetscObjectCompose((PetscObject)C, "_MatIS_PtAP", (PetscObject)c));
    211:   PetscCall(PetscContainerDestroy(&c));
    212:   ptap->fill = fill;

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

    216:   PetscCall(ISLocalToGlobalMappingGetBlockSize(cl2g, &ibs));
    217:   PetscCall(ISLocalToGlobalMappingGetSize(cl2g, &N));
    218:   PetscCall(ISLocalToGlobalMappingGetBlockIndices(cl2g, &garray));
    219:   PetscCall(ISCreateBlock(comm, ibs, N / ibs, garray, PETSC_COPY_VALUES, &ptap->ris0));
    220:   PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(cl2g, &garray));

    222:   PetscCall(MatCreateSubMatrix(P, ptap->ris0, NULL, MAT_INITIAL_MATRIX, &PT));
    223:   PetscCall(MatGetNonzeroColumnsLocal_Private(PT, &ptap->cis0));
    224:   PetscCall(ISLocalToGlobalMappingCreateIS(ptap->cis0, &Ccl2g));
    225:   PetscCall(MatDestroy(&PT));

    227:   Crl2g = NULL;
    228:   if (rl2g != cl2g) { /* unsymmetric A mapping */
    229:     PetscBool same, lsame = PETSC_FALSE;
    230:     PetscInt  N1, ibs1;

    232:     PetscCall(ISLocalToGlobalMappingGetSize(rl2g, &N1));
    233:     PetscCall(ISLocalToGlobalMappingGetBlockSize(rl2g, &ibs1));
    234:     PetscCall(ISLocalToGlobalMappingGetBlockIndices(rl2g, &garray));
    235:     PetscCall(ISCreateBlock(comm, ibs, N1 / ibs, garray, PETSC_COPY_VALUES, &ptap->ris1));
    236:     PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(rl2g, &garray));
    237:     if (ibs1 == ibs && N1 == N) { /* check if the l2gmaps are the same */
    238:       const PetscInt *i1, *i2;

    240:       PetscCall(ISBlockGetIndices(ptap->ris0, &i1));
    241:       PetscCall(ISBlockGetIndices(ptap->ris1, &i2));
    242:       PetscCall(PetscArraycmp(i1, i2, N, &lsame));
    243:     }
    244:     PetscCallMPI(MPIU_Allreduce(&lsame, &same, 1, MPI_C_BOOL, MPI_LAND, comm));
    245:     if (same) {
    246:       PetscCall(ISDestroy(&ptap->ris1));
    247:     } else {
    248:       PetscCall(MatCreateSubMatrix(P, ptap->ris1, NULL, MAT_INITIAL_MATRIX, &PT));
    249:       PetscCall(MatGetNonzeroColumnsLocal_Private(PT, &ptap->cis1));
    250:       PetscCall(ISLocalToGlobalMappingCreateIS(ptap->cis1, &Crl2g));
    251:       PetscCall(MatDestroy(&PT));
    252:     }
    253:   }
    254:   /* Not sure about this
    255:   if (!Crl2g) {
    256:     PetscCall(MatGetBlockSize(C,&ibs));
    257:     PetscCall(ISLocalToGlobalMappingSetBlockSize(Ccl2g,ibs));
    258:   }
    259: */
    260:   PetscCall(MatSetLocalToGlobalMapping(C, Crl2g ? Crl2g : Ccl2g, Ccl2g));
    261:   PetscCall(ISLocalToGlobalMappingDestroy(&Crl2g));
    262:   PetscCall(ISLocalToGlobalMappingDestroy(&Ccl2g));

    264:   C->ops->ptapnumeric = MatPtAPNumeric_IS_XAIJ;
    265:   PetscFunctionReturn(PETSC_SUCCESS);
    266: }

    268: static PetscErrorCode MatProductSymbolic_PtAP_IS_XAIJ(Mat C)
    269: {
    270:   Mat_Product *product = C->product;
    271:   Mat          A = product->A, P = product->B;
    272:   PetscReal    fill = product->fill;

    274:   PetscFunctionBegin;
    275:   PetscCall(MatPtAPSymbolic_IS_XAIJ(A, P, fill, C));
    276:   C->ops->productnumeric = MatProductNumeric_PtAP;
    277:   PetscFunctionReturn(PETSC_SUCCESS);
    278: }

    280: static PetscErrorCode MatProductSetFromOptions_IS_XAIJ_PtAP(Mat C)
    281: {
    282:   PetscFunctionBegin;
    283:   C->ops->productsymbolic = MatProductSymbolic_PtAP_IS_XAIJ;
    284:   PetscFunctionReturn(PETSC_SUCCESS);
    285: }

    287: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_IS_XAIJ(Mat C)
    288: {
    289:   Mat_Product *product = C->product;

    291:   PetscFunctionBegin;
    292:   if (product->type == MATPRODUCT_PtAP) PetscCall(MatProductSetFromOptions_IS_XAIJ_PtAP(C));
    293:   PetscFunctionReturn(PETSC_SUCCESS);
    294: }

    296: static PetscErrorCode MatISContainerDestroyFields_Private(PetscCtxRt ptr)
    297: {
    298:   MatISLocalFields lf = *(MatISLocalFields *)ptr;
    299:   PetscInt         i;

    301:   PetscFunctionBegin;
    302:   for (i = 0; i < lf->nr; i++) PetscCall(ISDestroy(&lf->rf[i]));
    303:   for (i = 0; i < lf->nc; i++) PetscCall(ISDestroy(&lf->cf[i]));
    304:   PetscCall(PetscFree2(lf->rf, lf->cf));
    305:   PetscCall(PetscFree(lf));
    306:   PetscFunctionReturn(PETSC_SUCCESS);
    307: }

    309: static PetscErrorCode MatConvert_SeqXAIJ_IS(Mat A, MatType type, MatReuse reuse, Mat *newmat)
    310: {
    311:   Mat B, lB;

    313:   PetscFunctionBegin;
    314:   if (reuse != MAT_REUSE_MATRIX) {
    315:     ISLocalToGlobalMapping rl2g, cl2g;
    316:     PetscInt               bs;
    317:     IS                     is;

    319:     PetscCall(MatGetBlockSize(A, &bs));
    320:     PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), A->rmap->n / bs, 0, 1, &is));
    321:     if (bs > 1) {
    322:       IS       is2;
    323:       PetscInt i, *aux;

    325:       PetscCall(ISGetLocalSize(is, &i));
    326:       PetscCall(ISGetIndices(is, (const PetscInt **)&aux));
    327:       PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), bs, i, aux, PETSC_COPY_VALUES, &is2));
    328:       PetscCall(ISRestoreIndices(is, (const PetscInt **)&aux));
    329:       PetscCall(ISDestroy(&is));
    330:       is = is2;
    331:     }
    332:     PetscCall(ISSetIdentity(is));
    333:     PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g));
    334:     PetscCall(ISDestroy(&is));
    335:     PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), A->cmap->n / bs, 0, 1, &is));
    336:     if (bs > 1) {
    337:       IS       is2;
    338:       PetscInt i, *aux;

    340:       PetscCall(ISGetLocalSize(is, &i));
    341:       PetscCall(ISGetIndices(is, (const PetscInt **)&aux));
    342:       PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), bs, i, aux, PETSC_COPY_VALUES, &is2));
    343:       PetscCall(ISRestoreIndices(is, (const PetscInt **)&aux));
    344:       PetscCall(ISDestroy(&is));
    345:       is = is2;
    346:     }
    347:     PetscCall(ISSetIdentity(is));
    348:     PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g));
    349:     PetscCall(ISDestroy(&is));
    350:     PetscCall(MatCreateIS(PetscObjectComm((PetscObject)A), bs, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N, rl2g, cl2g, &B));
    351:     PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
    352:     PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
    353:     PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &lB));
    354:     if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
    355:   } else {
    356:     B = *newmat;
    357:     PetscCall(PetscObjectReference((PetscObject)A));
    358:     lB = A;
    359:   }
    360:   PetscCall(MatISSetLocalMat(B, lB));
    361:   PetscCall(MatDestroy(&lB));
    362:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
    363:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
    364:   if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &B));
    365:   PetscFunctionReturn(PETSC_SUCCESS);
    366: }

    368: static PetscErrorCode MatISScaleDisassembling_Private(Mat A)
    369: {
    370:   Mat_IS         *matis = (Mat_IS *)A->data;
    371:   PetscScalar    *aa;
    372:   const PetscInt *ii, *jj;
    373:   PetscInt        i, n, m;
    374:   PetscInt       *ecount, **eneighs;
    375:   PetscBool       flg;

    377:   PetscFunctionBegin;
    378:   PetscCall(MatGetRowIJ(matis->A, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg));
    379:   PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure");
    380:   PetscCall(ISLocalToGlobalMappingGetNodeInfo(matis->rmapping, &n, &ecount, &eneighs));
    381:   PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected %" PetscInt_FMT " != %" PetscInt_FMT, m, n);
    382:   PetscCall(MatSeqAIJGetArray(matis->A, &aa));
    383:   for (i = 0; i < n; i++) {
    384:     if (ecount[i] > 1) {
    385:       PetscInt j;

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

    391:         for (p = 0; p < ecount[i]; p++) {
    392:           for (p2 = 0; p2 < ecount[i2]; p2++) {
    393:             if (eneighs[i][p] == eneighs[i2][p2]) {
    394:               scal += 1.0;
    395:               break;
    396:             }
    397:           }
    398:         }
    399:         if (scal) aa[j] /= scal;
    400:       }
    401:     }
    402:   }
    403:   PetscCall(ISLocalToGlobalMappingRestoreNodeInfo(matis->rmapping, &n, &ecount, &eneighs));
    404:   PetscCall(MatSeqAIJRestoreArray(matis->A, &aa));
    405:   PetscCall(MatRestoreRowIJ(matis->A, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg));
    406:   PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot restore IJ structure");
    407:   PetscFunctionReturn(PETSC_SUCCESS);
    408: }

    410: typedef enum {
    411:   MAT_IS_DISASSEMBLE_L2G_NATURAL,
    412:   MAT_IS_DISASSEMBLE_L2G_MAT,
    413:   MAT_IS_DISASSEMBLE_L2G_ND
    414: } MatISDisassemblel2gType;

    416: static PetscErrorCode MatMPIXAIJComputeLocalToGlobalMapping_Private(Mat A, ISLocalToGlobalMapping *l2g)
    417: {
    418:   Mat                     Ad, Ao;
    419:   IS                      is, ndmap, ndsub;
    420:   MPI_Comm                comm;
    421:   const PetscInt         *garray, *ndmapi;
    422:   PetscInt                bs, i, cnt, nl, *ncount, *ndmapc;
    423:   PetscBool               ismpiaij, ismpibaij;
    424:   const char *const       MatISDisassemblel2gTypes[] = {"NATURAL", "MAT", "ND", "MatISDisassemblel2gType", "MAT_IS_DISASSEMBLE_L2G_", NULL};
    425:   MatISDisassemblel2gType mode                       = MAT_IS_DISASSEMBLE_L2G_NATURAL;
    426:   MatPartitioning         part;
    427:   PetscSF                 sf;
    428:   PetscObject             dm;

    430:   PetscFunctionBegin;
    431:   PetscOptionsBegin(PetscObjectComm((PetscObject)A), ((PetscObject)A)->prefix, "MatIS l2g disassembling options", "Mat");
    432:   PetscCall(PetscOptionsEnum("-mat_is_disassemble_l2g_type", "Type of local-to-global mapping to be used for disassembling", "MatISDisassemblel2gType", MatISDisassemblel2gTypes, (PetscEnum)mode, (PetscEnum *)&mode, NULL));
    433:   PetscOptionsEnd();
    434:   if (mode == MAT_IS_DISASSEMBLE_L2G_MAT) {
    435:     PetscCall(MatGetLocalToGlobalMapping(A, l2g, NULL));
    436:     PetscFunctionReturn(PETSC_SUCCESS);
    437:   }
    438:   PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
    439:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij));
    440:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIBAIJ, &ismpibaij));
    441:   PetscCall(MatGetBlockSize(A, &bs));
    442:   switch (mode) {
    443:   case MAT_IS_DISASSEMBLE_L2G_ND:
    444:     PetscCall(MatPartitioningCreate(comm, &part));
    445:     PetscCall(MatPartitioningSetAdjacency(part, A));
    446:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)part, ((PetscObject)A)->prefix));
    447:     PetscCall(MatPartitioningSetFromOptions(part));
    448:     PetscCall(MatPartitioningApplyND(part, &ndmap));
    449:     PetscCall(MatPartitioningDestroy(&part));
    450:     PetscCall(ISBuildTwoSided(ndmap, NULL, &ndsub));
    451:     PetscCall(MatMPIAIJSetUseScalableIncreaseOverlap(A, PETSC_TRUE));
    452:     PetscCall(MatIncreaseOverlap(A, 1, &ndsub, 1));
    453:     PetscCall(ISLocalToGlobalMappingCreateIS(ndsub, l2g));

    455:     /* it may happen that a separator node is not properly shared */
    456:     PetscCall(ISLocalToGlobalMappingGetNodeInfo(*l2g, &nl, &ncount, NULL));
    457:     PetscCall(PetscSFCreate(comm, &sf));
    458:     PetscCall(ISLocalToGlobalMappingGetIndices(*l2g, &garray));
    459:     PetscCall(PetscSFSetGraphLayout(sf, A->rmap, nl, NULL, PETSC_OWN_POINTER, garray));
    460:     PetscCall(ISLocalToGlobalMappingRestoreIndices(*l2g, &garray));
    461:     PetscCall(PetscCalloc1(A->rmap->n, &ndmapc));
    462:     PetscCall(PetscSFReduceBegin(sf, MPIU_INT, ncount, ndmapc, MPI_REPLACE));
    463:     PetscCall(PetscSFReduceEnd(sf, MPIU_INT, ncount, ndmapc, MPI_REPLACE));
    464:     PetscCall(ISLocalToGlobalMappingRestoreNodeInfo(*l2g, NULL, &ncount, NULL));
    465:     PetscCall(ISGetIndices(ndmap, &ndmapi));
    466:     for (i = 0, cnt = 0; i < A->rmap->n; i++)
    467:       if (ndmapi[i] < 0 && ndmapc[i] < 2) cnt++;

    469:     PetscCallMPI(MPIU_Allreduce(&cnt, &i, 1, MPIU_INT, MPI_MAX, comm));
    470:     if (i) { /* we detected isolated separator nodes */
    471:       Mat                    A2, A3;
    472:       IS                    *workis, is2;
    473:       PetscScalar           *vals;
    474:       PetscInt               gcnt = i, *dnz, *onz, j, *lndmapi;
    475:       ISLocalToGlobalMapping ll2g;
    476:       PetscBool              flg;
    477:       const PetscInt        *ii, *jj;

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

    483:       PetscCall(PetscMalloc1(nl, &lndmapi));
    484:       PetscCall(PetscSFBcastBegin(sf, MPIU_INT, dnz, lndmapi, MPI_REPLACE));

    486:       /* compute adjacency of isolated separators node */
    487:       PetscCall(PetscMalloc1(gcnt, &workis));
    488:       for (i = 0, cnt = 0; i < A->rmap->n; i++) {
    489:         if (ndmapi[i] < 0 && ndmapc[i] < 2) PetscCall(ISCreateStride(comm, 1, i + A->rmap->rstart, 1, &workis[cnt++]));
    490:       }
    491:       for (i = cnt; i < gcnt; i++) PetscCall(ISCreateStride(comm, 0, 0, 1, &workis[i]));
    492:       for (i = 0; i < gcnt; i++) {
    493:         PetscCall(PetscObjectSetName((PetscObject)workis[i], "ISOLATED"));
    494:         PetscCall(ISViewFromOptions(workis[i], NULL, "-view_isolated_separators"));
    495:       }

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

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

    503:       /* communicate new layers : create a matrix and transpose it */
    504:       PetscCall(PetscArrayzero(dnz, A->rmap->n));
    505:       PetscCall(PetscArrayzero(onz, A->rmap->n));
    506:       for (i = 0, j = 0; i < A->rmap->n; i++) {
    507:         if (ndmapi[i] < 0 && ndmapc[i] < 2) {
    508:           const PetscInt *idxs;
    509:           PetscInt        s;

    511:           PetscCall(ISGetLocalSize(workis[j], &s));
    512:           PetscCall(ISGetIndices(workis[j], &idxs));
    513:           PetscCall(MatPreallocateSet(i + A->rmap->rstart, s, idxs, dnz, onz));
    514:           j++;
    515:         }
    516:       }
    517:       PetscCheck(j == cnt, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected local count %" PetscInt_FMT " != %" PetscInt_FMT, j, cnt);

    519:       for (i = 0; i < gcnt; i++) {
    520:         PetscCall(PetscObjectSetName((PetscObject)workis[i], "EXTENDED"));
    521:         PetscCall(ISViewFromOptions(workis[i], NULL, "-view_isolated_separators"));
    522:       }

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

    528:       PetscCall(MatCreate(comm, &A2));
    529:       PetscCall(MatSetType(A2, MATMPIAIJ));
    530:       PetscCall(MatSetSizes(A2, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
    531:       PetscCall(MatMPIAIJSetPreallocation(A2, 0, dnz, 0, onz));
    532:       PetscCall(MatSetOption(A2, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
    533:       for (i = 0, j = 0; i < A2->rmap->n; i++) {
    534:         PetscInt        row = i + A2->rmap->rstart, s = dnz[i] + onz[i];
    535:         const PetscInt *idxs;

    537:         if (s) {
    538:           PetscCall(ISGetIndices(workis[j], &idxs));
    539:           PetscCall(MatSetValues(A2, 1, &row, s, idxs, vals, INSERT_VALUES));
    540:           PetscCall(ISRestoreIndices(workis[j], &idxs));
    541:           j++;
    542:         }
    543:       }
    544:       PetscCheck(j == cnt, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected local count %" PetscInt_FMT " != %" PetscInt_FMT, j, cnt);
    545:       PetscCall(PetscFree(vals));
    546:       PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY));
    547:       PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY));
    548:       PetscCall(MatTranspose(A2, MAT_INPLACE_MATRIX, &A2));

    550:       /* extract submatrix corresponding to the coupling "owned separators" x "isolated separators" */
    551:       for (i = 0, j = 0; i < nl; i++)
    552:         if (lndmapi[i] >= 0) lndmapi[j++] = lndmapi[i];
    553:       PetscCall(ISCreateGeneral(comm, j, lndmapi, PETSC_USE_POINTER, &is));
    554:       PetscCall(MatMPIAIJGetLocalMatCondensed(A2, MAT_INITIAL_MATRIX, &is, NULL, &A3));
    555:       PetscCall(ISDestroy(&is));
    556:       PetscCall(MatDestroy(&A2));

    558:       /* extend local to global map to include connected isolated separators */
    559:       PetscCall(PetscObjectQuery((PetscObject)A3, "_petsc_GetLocalMatCondensed_iscol", (PetscObject *)&is));
    560:       PetscCheck(is, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing column map");
    561:       PetscCall(ISLocalToGlobalMappingCreateIS(is, &ll2g));
    562:       PetscCall(MatGetRowIJ(A3, 0, PETSC_FALSE, PETSC_FALSE, &i, &ii, &jj, &flg));
    563:       PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure");
    564:       PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ii[i], jj, PETSC_COPY_VALUES, &is));
    565:       PetscCall(MatRestoreRowIJ(A3, 0, PETSC_FALSE, PETSC_FALSE, &i, &ii, &jj, &flg));
    566:       PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure");
    567:       PetscCall(ISLocalToGlobalMappingApplyIS(ll2g, is, &is2));
    568:       PetscCall(ISDestroy(&is));
    569:       PetscCall(ISLocalToGlobalMappingDestroy(&ll2g));

    571:       /* add new nodes to the local-to-global map */
    572:       PetscCall(ISLocalToGlobalMappingDestroy(l2g));
    573:       PetscCall(ISExpand(ndsub, is2, &is));
    574:       PetscCall(ISDestroy(&is2));
    575:       PetscCall(ISLocalToGlobalMappingCreateIS(is, l2g));
    576:       PetscCall(ISDestroy(&is));

    578:       PetscCall(MatDestroy(&A3));
    579:       PetscCall(PetscFree(lndmapi));
    580:       MatPreallocateEnd(dnz, onz);
    581:       for (i = 0; i < gcnt; i++) PetscCall(ISDestroy(&workis[i]));
    582:       PetscCall(PetscFree(workis));
    583:     }
    584:     PetscCall(ISRestoreIndices(ndmap, &ndmapi));
    585:     PetscCall(PetscSFDestroy(&sf));
    586:     PetscCall(PetscFree(ndmapc));
    587:     PetscCall(ISDestroy(&ndmap));
    588:     PetscCall(ISDestroy(&ndsub));
    589:     PetscCall(ISLocalToGlobalMappingSetBlockSize(*l2g, bs));
    590:     PetscCall(ISLocalToGlobalMappingViewFromOptions(*l2g, NULL, "-mat_is_nd_l2g_view"));
    591:     break;
    592:   case MAT_IS_DISASSEMBLE_L2G_NATURAL:
    593:     PetscCall(PetscObjectQuery((PetscObject)A, "__PETSc_dm", &dm));
    594:     if (dm) { /* if a matrix comes from a DM, most likely we can use the l2gmap if any */
    595:       PetscCall(MatGetLocalToGlobalMapping(A, l2g, NULL));
    596:       PetscCall(PetscObjectReference((PetscObject)*l2g));
    597:       if (*l2g) PetscFunctionReturn(PETSC_SUCCESS);
    598:     }
    599:     if (ismpiaij) {
    600:       PetscCall(MatMPIAIJGetSeqAIJ(A, &Ad, &Ao, &garray));
    601:     } else if (ismpibaij) {
    602:       PetscCall(MatMPIBAIJGetSeqBAIJ(A, &Ad, &Ao, &garray));
    603:     } else SETERRQ(comm, PETSC_ERR_SUP, "Type %s", ((PetscObject)A)->type_name);
    604:     if (A->rmap->n) {
    605:       PetscInt dc, oc, stc, *aux;

    607:       PetscCall(MatGetLocalSize(Ad, NULL, &dc));
    608:       PetscCall(MatGetLocalSize(Ao, NULL, &oc));
    609:       PetscCheck(!oc || garray, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "garray not present");
    610:       PetscCall(MatGetOwnershipRangeColumn(A, &stc, NULL));
    611:       PetscCall(PetscMalloc1((dc + oc) / bs, &aux));
    612:       for (i = 0; i < dc / bs; i++) aux[i] = i + stc / bs;
    613:       for (i = 0; i < oc / bs; i++) aux[i + dc / bs] = (ismpiaij ? garray[i * bs] / bs : garray[i]);
    614:       PetscCall(ISCreateBlock(comm, bs, (dc + oc) / bs, aux, PETSC_OWN_POINTER, &is));
    615:     } else {
    616:       PetscCall(ISCreateBlock(comm, 1, 0, NULL, PETSC_OWN_POINTER, &is));
    617:     }
    618:     PetscCall(ISLocalToGlobalMappingCreateIS(is, l2g));
    619:     PetscCall(ISDestroy(&is));
    620:     break;
    621:   default:
    622:     SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unsupported l2g disassembling type %d", mode);
    623:   }
    624:   PetscFunctionReturn(PETSC_SUCCESS);
    625: }

    627: PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat A, MatType type, MatReuse reuse, Mat *newmat)
    628: {
    629:   Mat                    lA, Ad, Ao, B = NULL;
    630:   ISLocalToGlobalMapping rl2g, cl2g;
    631:   IS                     is;
    632:   MPI_Comm               comm;
    633:   void                  *ptrs[2];
    634:   const char            *names[2] = {"_convert_csr_aux", "_convert_csr_data"};
    635:   const PetscInt        *garray;
    636:   PetscScalar           *dd, *od, *aa, *data;
    637:   const PetscInt        *di, *dj, *oi, *oj;
    638:   const PetscInt        *odi, *odj, *ooi, *ooj;
    639:   PetscInt              *aux, *ii, *jj;
    640:   PetscInt               rbs, cbs, lc, dr, dc, oc, str, stc, nnz, i, jd, jo, cum;
    641:   PetscBool              flg, ismpiaij, ismpibaij, was_inplace = PETSC_FALSE, cong;
    642:   PetscMPIInt            size;

    644:   PetscFunctionBegin;
    645:   PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
    646:   PetscCallMPI(MPI_Comm_size(comm, &size));
    647:   if (size == 1) {
    648:     PetscCall(MatConvert_SeqXAIJ_IS(A, type, reuse, newmat));
    649:     PetscFunctionReturn(PETSC_SUCCESS);
    650:   }
    651:   PetscCall(MatGetBlockSizes(A, &rbs, &cbs));
    652:   PetscCall(MatHasCongruentLayouts(A, &cong));
    653:   if (reuse != MAT_REUSE_MATRIX && cong && rbs == cbs) {
    654:     PetscCall(MatMPIXAIJComputeLocalToGlobalMapping_Private(A, &rl2g));
    655:     PetscCall(MatCreate(comm, &B));
    656:     PetscCall(MatSetType(B, MATIS));
    657:     PetscCall(MatSetSizes(B, A->rmap->n, A->rmap->n, A->rmap->N, A->rmap->N));
    658:     PetscCall(MatSetLocalToGlobalMapping(B, rl2g, rl2g));
    659:     PetscCall(MatSetBlockSizes(B, rbs, rbs));
    660:     PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
    661:     if (reuse == MAT_INPLACE_MATRIX) was_inplace = PETSC_TRUE;
    662:     reuse = MAT_REUSE_MATRIX;
    663:   }
    664:   if (reuse == MAT_REUSE_MATRIX) {
    665:     Mat            *newlA, lA;
    666:     IS              rows, cols;
    667:     const PetscInt *ridx, *cidx;
    668:     PetscInt        nr, nc;

    670:     if (!B) B = *newmat;
    671:     PetscCall(MatISGetLocalToGlobalMapping(B, &rl2g, &cl2g));
    672:     PetscCall(ISLocalToGlobalMappingGetBlockIndices(rl2g, &ridx));
    673:     PetscCall(ISLocalToGlobalMappingGetBlockIndices(cl2g, &cidx));
    674:     PetscCall(ISLocalToGlobalMappingGetSize(rl2g, &nr));
    675:     PetscCall(ISLocalToGlobalMappingGetSize(cl2g, &nc));
    676:     PetscCall(ISLocalToGlobalMappingGetBlockSize(rl2g, &rbs));
    677:     PetscCall(ISLocalToGlobalMappingGetBlockSize(cl2g, &cbs));
    678:     PetscCall(ISCreateBlock(comm, rbs, nr / rbs, ridx, PETSC_USE_POINTER, &rows));
    679:     if (rl2g != cl2g) {
    680:       PetscCall(ISCreateBlock(comm, cbs, nc / cbs, cidx, PETSC_USE_POINTER, &cols));
    681:     } else {
    682:       PetscCall(PetscObjectReference((PetscObject)rows));
    683:       cols = rows;
    684:     }
    685:     PetscCall(MatISGetLocalMat(B, &lA));
    686:     PetscCall(MatCreateSubMatrices(A, 1, &rows, &cols, MAT_INITIAL_MATRIX, &newlA));
    687:     PetscCall(MatConvert(newlA[0], MATSEQAIJ, MAT_INPLACE_MATRIX, &newlA[0]));
    688:     PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(rl2g, &ridx));
    689:     PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(cl2g, &cidx));
    690:     PetscCall(ISDestroy(&rows));
    691:     PetscCall(ISDestroy(&cols));
    692:     if (!lA->preallocated) { /* first time */
    693:       PetscCall(MatDuplicate(newlA[0], MAT_COPY_VALUES, &lA));
    694:       PetscCall(MatISSetLocalMat(B, lA));
    695:       PetscCall(PetscObjectDereference((PetscObject)lA));
    696:     }
    697:     PetscCall(MatCopy(newlA[0], lA, SAME_NONZERO_PATTERN));
    698:     PetscCall(MatDestroySubMatrices(1, &newlA));
    699:     PetscCall(MatISScaleDisassembling_Private(B));
    700:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
    701:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
    702:     if (was_inplace) PetscCall(MatHeaderReplace(A, &B));
    703:     else *newmat = B;
    704:     PetscFunctionReturn(PETSC_SUCCESS);
    705:   }
    706:   /* general case, just compress out the column space */
    707:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij));
    708:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIBAIJ, &ismpibaij));
    709:   if (ismpiaij) {
    710:     cbs = 1; /* We cannot guarantee the off-process matrix will respect the column block size */
    711:     PetscCall(MatMPIAIJGetSeqAIJ(A, &Ad, &Ao, &garray));
    712:   } else if (ismpibaij) {
    713:     PetscCall(MatMPIBAIJGetSeqBAIJ(A, &Ad, &Ao, &garray));
    714:     PetscCall(MatConvert(Ad, MATSEQAIJ, MAT_INITIAL_MATRIX, &Ad));
    715:     PetscCall(MatConvert(Ao, MATSEQAIJ, MAT_INITIAL_MATRIX, &Ao));
    716:   } else SETERRQ(comm, PETSC_ERR_SUP, "Type %s", ((PetscObject)A)->type_name);
    717:   PetscCall(MatSeqAIJGetArray(Ad, &dd));
    718:   PetscCall(MatSeqAIJGetArray(Ao, &od));

    720:   /* access relevant information from MPIAIJ */
    721:   PetscCall(MatGetOwnershipRange(A, &str, NULL));
    722:   PetscCall(MatGetOwnershipRangeColumn(A, &stc, NULL));
    723:   PetscCall(MatGetLocalSize(Ad, &dr, &dc));
    724:   PetscCall(MatGetLocalSize(Ao, NULL, &oc));
    725:   PetscCheck(!oc || garray, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "garray not present");

    727:   PetscCall(MatGetRowIJ(Ad, 0, PETSC_FALSE, PETSC_FALSE, &i, &di, &dj, &flg));
    728:   PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure");
    729:   PetscCall(MatGetRowIJ(Ao, 0, PETSC_FALSE, PETSC_FALSE, &i, &oi, &oj, &flg));
    730:   PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure");
    731:   nnz = di[dr] + oi[dr];
    732:   /* store original pointers to be restored later */
    733:   odi = di;
    734:   odj = dj;
    735:   ooi = oi;
    736:   ooj = oj;

    738:   /* generate l2g maps for rows and cols */
    739:   PetscCall(ISCreateStride(comm, dr / rbs, str / rbs, 1, &is));
    740:   if (rbs > 1) {
    741:     IS is2;

    743:     PetscCall(ISGetLocalSize(is, &i));
    744:     PetscCall(ISGetIndices(is, (const PetscInt **)&aux));
    745:     PetscCall(ISCreateBlock(comm, rbs, i, aux, PETSC_COPY_VALUES, &is2));
    746:     PetscCall(ISRestoreIndices(is, (const PetscInt **)&aux));
    747:     PetscCall(ISDestroy(&is));
    748:     is = is2;
    749:   }
    750:   PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g));
    751:   PetscCall(ISDestroy(&is));
    752:   if (dr) {
    753:     PetscCall(PetscMalloc1((dc + oc) / cbs, &aux));
    754:     for (i = 0; i < dc / cbs; i++) aux[i] = i + stc / cbs;
    755:     for (i = 0; i < oc / cbs; i++) aux[i + dc / cbs] = garray[i];
    756:     PetscCall(ISCreateBlock(comm, cbs, (dc + oc) / cbs, aux, PETSC_OWN_POINTER, &is));
    757:     lc = dc + oc;
    758:   } else {
    759:     PetscCall(ISCreateBlock(comm, cbs, 0, NULL, PETSC_OWN_POINTER, &is));
    760:     lc = 0;
    761:   }
    762:   PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g));
    763:   PetscCall(ISDestroy(&is));

    765:   /* create MATIS object */
    766:   PetscCall(MatCreate(comm, &B));
    767:   PetscCall(MatSetSizes(B, dr, dc, PETSC_DECIDE, PETSC_DECIDE));
    768:   PetscCall(MatSetType(B, MATIS));
    769:   PetscCall(MatSetBlockSizes(B, rbs, cbs));
    770:   PetscCall(MatSetLocalToGlobalMapping(B, rl2g, cl2g));
    771:   PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
    772:   PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));

    774:   /* merge local matrices */
    775:   PetscCall(PetscMalloc1(nnz + dr + 1, &aux));
    776:   PetscCall(PetscMalloc1(nnz, &data));
    777:   ii  = aux;
    778:   jj  = aux + dr + 1;
    779:   aa  = data;
    780:   *ii = *(di++) + *(oi++);
    781:   for (jd = 0, jo = 0, cum = 0; *ii < nnz; cum++) {
    782:     for (; jd < *di; jd++) {
    783:       *jj++ = *dj++;
    784:       *aa++ = *dd++;
    785:     }
    786:     for (; jo < *oi; jo++) {
    787:       *jj++ = *oj++ + dc;
    788:       *aa++ = *od++;
    789:     }
    790:     *(++ii) = *(di++) + *(oi++);
    791:   }
    792:   for (; cum < dr; cum++) *(++ii) = nnz;

    794:   PetscCall(MatRestoreRowIJ(Ad, 0, PETSC_FALSE, PETSC_FALSE, &i, &odi, &odj, &flg));
    795:   PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot restore IJ structure");
    796:   PetscCall(MatRestoreRowIJ(Ao, 0, PETSC_FALSE, PETSC_FALSE, &i, &ooi, &ooj, &flg));
    797:   PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot restore IJ structure");
    798:   PetscCall(MatSeqAIJRestoreArray(Ad, &dd));
    799:   PetscCall(MatSeqAIJRestoreArray(Ao, &od));

    801:   ii = aux;
    802:   jj = aux + dr + 1;
    803:   aa = data;
    804:   PetscCall(MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, dr, lc, ii, jj, aa, &lA));

    806:   /* create containers to destroy the data */
    807:   ptrs[0] = aux;
    808:   ptrs[1] = data;
    809:   for (i = 0; i < 2; i++) PetscCall(PetscObjectContainerCompose((PetscObject)lA, names[i], ptrs[i], PetscCtxDestroyDefault));
    810:   if (ismpibaij) { /* destroy converted local matrices */
    811:     PetscCall(MatDestroy(&Ad));
    812:     PetscCall(MatDestroy(&Ao));
    813:   }

    815:   /* finalize matrix */
    816:   PetscCall(MatISSetLocalMat(B, lA));
    817:   PetscCall(MatDestroy(&lA));
    818:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
    819:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
    820:   if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &B));
    821:   else *newmat = B;
    822:   PetscFunctionReturn(PETSC_SUCCESS);
    823: }

    825: PETSC_INTERN PetscErrorCode MatConvert_Nest_IS(Mat A, MatType type, MatReuse reuse, Mat *newmat)
    826: {
    827:   Mat                  **nest, *snest, **rnest, lA, B;
    828:   IS                    *iscol, *isrow, *islrow, *islcol;
    829:   ISLocalToGlobalMapping rl2g, cl2g;
    830:   MPI_Comm               comm;
    831:   PetscInt              *lr, *lc, *l2gidxs;
    832:   PetscInt               i, j, nr, nc, rbs, cbs;
    833:   PetscBool              convert, lreuse, *istrans;
    834:   PetscBool3             allow_repeated = PETSC_BOOL3_UNKNOWN;

    836:   PetscFunctionBegin;
    837:   PetscCall(MatNestGetSubMats(A, &nr, &nc, &nest));
    838:   lreuse = PETSC_FALSE;
    839:   rnest  = NULL;
    840:   if (reuse == MAT_REUSE_MATRIX) {
    841:     PetscBool ismatis, isnest;

    843:     PetscCall(PetscObjectTypeCompare((PetscObject)*newmat, MATIS, &ismatis));
    844:     PetscCheck(ismatis, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_USER, "Cannot reuse matrix of type %s", ((PetscObject)*newmat)->type_name);
    845:     PetscCall(MatISGetLocalMat(*newmat, &lA));
    846:     PetscCall(PetscObjectTypeCompare((PetscObject)lA, MATNEST, &isnest));
    847:     if (isnest) {
    848:       PetscCall(MatNestGetSubMats(lA, &i, &j, &rnest));
    849:       lreuse = (PetscBool)(i == nr && j == nc);
    850:       if (!lreuse) rnest = NULL;
    851:     }
    852:   }
    853:   PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
    854:   PetscCall(PetscCalloc2(nr, &lr, nc, &lc));
    855:   PetscCall(PetscCalloc6(nr, &isrow, nc, &iscol, nr, &islrow, nc, &islcol, nr * nc, &snest, nr * nc, &istrans));
    856:   PetscCall(MatNestGetISs(A, isrow, iscol));
    857:   for (i = 0; i < nr; i++) {
    858:     for (j = 0; j < nc; j++) {
    859:       PetscBool ismatis, sallow;
    860:       PetscInt  l1, l2, lb1, lb2, ij = i * nc + j;

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

    865:       /* Nested matrices should be of type MATIS */
    866:       PetscCall(PetscObjectTypeCompare((PetscObject)nest[i][j], MATTRANSPOSEVIRTUAL, &istrans[ij]));
    867:       if (istrans[ij]) {
    868:         Mat T, lT;
    869:         PetscCall(MatTransposeGetMat(nest[i][j], &T));
    870:         PetscCall(PetscObjectTypeCompare((PetscObject)T, MATIS, &ismatis));
    871:         PetscCheck(ismatis, comm, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") (transposed) is not of type MATIS", i, j);
    872:         PetscCall(MatISGetAllowRepeated(T, &sallow));
    873:         PetscCall(MatISGetLocalMat(T, &lT));
    874:         PetscCall(MatCreateTranspose(lT, &snest[ij]));
    875:       } else {
    876:         PetscCall(PetscObjectTypeCompare((PetscObject)nest[i][j], MATIS, &ismatis));
    877:         PetscCheck(ismatis, comm, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") is not of type MATIS", i, j);
    878:         PetscCall(MatISGetAllowRepeated(nest[i][j], &sallow));
    879:         PetscCall(MatISGetLocalMat(nest[i][j], &snest[ij]));
    880:       }
    881:       if (allow_repeated == PETSC_BOOL3_UNKNOWN) allow_repeated = PetscBoolToBool3(sallow);
    882:       PetscCheck(sallow == PetscBool3ToBool(allow_repeated), comm, PETSC_ERR_SUP, "Cannot mix repeated and non repeated maps");

    884:       /* Check compatibility of local sizes */
    885:       PetscCall(MatGetSize(snest[ij], &l1, &l2));
    886:       PetscCall(MatGetBlockSizes(snest[ij], &lb1, &lb2));
    887:       if (!l1 || !l2) continue;
    888:       PetscCheck(!lr[i] || l1 == lr[i], PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid local size %" PetscInt_FMT " != %" PetscInt_FMT, i, j, lr[i], l1);
    889:       PetscCheck(!lc[j] || l2 == lc[j], PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid local size %" PetscInt_FMT " != %" PetscInt_FMT, i, j, lc[j], l2);
    890:       lr[i] = l1;
    891:       lc[j] = l2;

    893:       /* check compatibility for local matrix reusage */
    894:       if (rnest && !rnest[i][j] != !snest[ij]) lreuse = PETSC_FALSE;
    895:     }
    896:   }

    898:   if (PetscDefined(USE_DEBUG)) {
    899:     /* Check compatibility of l2g maps for rows */
    900:     for (i = 0; i < nr; i++) {
    901:       rl2g = NULL;
    902:       for (j = 0; j < nc; j++) {
    903:         PetscInt n1, n2;

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

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

    922:           PetscCall(ISLocalToGlobalMappingGetSize(rl2g, &n2));
    923:           PetscCheck(n1 == n2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid row l2gmap size %" PetscInt_FMT " != %" PetscInt_FMT, i, j, n1, n2);
    924:           PetscCall(ISLocalToGlobalMappingGetIndices(cl2g, &idxs1));
    925:           PetscCall(ISLocalToGlobalMappingGetIndices(rl2g, &idxs2));
    926:           PetscCall(PetscArraycmp(idxs1, idxs2, n1, &same));
    927:           PetscCall(ISLocalToGlobalMappingRestoreIndices(cl2g, &idxs1));
    928:           PetscCall(ISLocalToGlobalMappingRestoreIndices(rl2g, &idxs2));
    929:           PetscCheck(same, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid row l2gmap", i, j);
    930:         }
    931:       }
    932:     }
    933:     /* Check compatibility of l2g maps for columns */
    934:     for (i = 0; i < nc; i++) {
    935:       rl2g = NULL;
    936:       for (j = 0; j < nr; j++) {
    937:         PetscInt n1, n2;

    939:         if (!nest[j][i]) continue;
    940:         if (istrans[j * nc + i]) {
    941:           Mat T;

    943:           PetscCall(MatTransposeGetMat(nest[j][i], &T));
    944:           PetscCall(MatISGetLocalToGlobalMapping(T, &cl2g, NULL));
    945:         } else {
    946:           PetscCall(MatISGetLocalToGlobalMapping(nest[j][i], NULL, &cl2g));
    947:         }
    948:         PetscCall(ISLocalToGlobalMappingGetSize(cl2g, &n1));
    949:         if (!n1) continue;
    950:         if (!rl2g) {
    951:           rl2g = cl2g;
    952:         } else {
    953:           const PetscInt *idxs1, *idxs2;
    954:           PetscBool       same;

    956:           PetscCall(ISLocalToGlobalMappingGetSize(rl2g, &n2));
    957:           PetscCheck(n1 == n2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid column l2gmap size %" PetscInt_FMT " != %" PetscInt_FMT, j, i, n1, n2);
    958:           PetscCall(ISLocalToGlobalMappingGetIndices(cl2g, &idxs1));
    959:           PetscCall(ISLocalToGlobalMappingGetIndices(rl2g, &idxs2));
    960:           PetscCall(PetscArraycmp(idxs1, idxs2, n1, &same));
    961:           PetscCall(ISLocalToGlobalMappingRestoreIndices(cl2g, &idxs1));
    962:           PetscCall(ISLocalToGlobalMappingRestoreIndices(rl2g, &idxs2));
    963:           PetscCheck(same, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid column l2gmap", j, i);
    964:         }
    965:       }
    966:     }
    967:   }

    969:   B = NULL;
    970:   if (reuse != MAT_REUSE_MATRIX) {
    971:     PetscInt stl;

    973:     /* Create l2g map for the rows of the new matrix and index sets for the local MATNEST */
    974:     for (i = 0, stl = 0; i < nr; i++) stl += lr[i];
    975:     PetscCall(PetscMalloc1(stl, &l2gidxs));
    976:     for (i = 0, stl = 0; i < nr; i++) {
    977:       Mat             usedmat;
    978:       Mat_IS         *matis;
    979:       const PetscInt *idxs;

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

    984:       /* l2gmap */
    985:       j       = 0;
    986:       usedmat = nest[i][j];
    987:       while (!usedmat && j < nc - 1) usedmat = nest[i][++j];
    988:       PetscCheck(usedmat, comm, PETSC_ERR_SUP, "Cannot find valid row mat");

    990:       if (istrans[i * nc + j]) {
    991:         Mat T;
    992:         PetscCall(MatTransposeGetMat(usedmat, &T));
    993:         usedmat = T;
    994:       }
    995:       matis = (Mat_IS *)usedmat->data;
    996:       PetscCall(ISGetIndices(isrow[i], &idxs));
    997:       if (istrans[i * nc + j]) {
    998:         PetscCall(PetscSFBcastBegin(matis->csf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
    999:         PetscCall(PetscSFBcastEnd(matis->csf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
    1000:       } else {
    1001:         PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
    1002:         PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
    1003:       }
    1004:       PetscCall(ISRestoreIndices(isrow[i], &idxs));
    1005:       stl += lr[i];
    1006:     }
    1007:     PetscCall(ISLocalToGlobalMappingCreate(comm, 1, stl, l2gidxs, PETSC_OWN_POINTER, &rl2g));

    1009:     /* Create l2g map for columns of the new matrix and index sets for the local MATNEST */
    1010:     for (i = 0, stl = 0; i < nc; i++) stl += lc[i];
    1011:     PetscCall(PetscMalloc1(stl, &l2gidxs));
    1012:     for (i = 0, stl = 0; i < nc; i++) {
    1013:       Mat             usedmat;
    1014:       Mat_IS         *matis;
    1015:       const PetscInt *idxs;

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

    1020:       /* l2gmap */
    1021:       j       = 0;
    1022:       usedmat = nest[j][i];
    1023:       while (!usedmat && j < nr - 1) usedmat = nest[++j][i];
    1024:       PetscCheck(usedmat, comm, PETSC_ERR_SUP, "Cannot find valid column mat");
    1025:       if (istrans[j * nc + i]) {
    1026:         Mat T;
    1027:         PetscCall(MatTransposeGetMat(usedmat, &T));
    1028:         usedmat = T;
    1029:       }
    1030:       matis = (Mat_IS *)usedmat->data;
    1031:       PetscCall(ISGetIndices(iscol[i], &idxs));
    1032:       if (istrans[j * nc + i]) {
    1033:         PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
    1034:         PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
    1035:       } else {
    1036:         PetscCall(PetscSFBcastBegin(matis->csf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
    1037:         PetscCall(PetscSFBcastEnd(matis->csf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE));
    1038:       }
    1039:       PetscCall(ISRestoreIndices(iscol[i], &idxs));
    1040:       stl += lc[i];
    1041:     }
    1042:     PetscCall(ISLocalToGlobalMappingCreate(comm, 1, stl, l2gidxs, PETSC_OWN_POINTER, &cl2g));

    1044:     /* Create MATIS */
    1045:     PetscCall(MatCreate(comm, &B));
    1046:     PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
    1047:     PetscCall(MatGetBlockSizes(A, &rbs, &cbs));
    1048:     PetscCall(MatSetBlockSizes(B, rbs, cbs));
    1049:     PetscCall(MatSetType(B, MATIS));
    1050:     PetscCall(MatISSetLocalMatType(B, MATNEST));
    1051:     PetscCall(MatISSetAllowRepeated(B, PetscBool3ToBool(allow_repeated)));
    1052:     { /* hack : avoid setup of scatters */
    1053:       Mat_IS *matis     = (Mat_IS *)B->data;
    1054:       matis->islocalref = B;
    1055:     }
    1056:     PetscCall(MatSetLocalToGlobalMapping(B, rl2g, cl2g));
    1057:     PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
    1058:     PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
    1059:     PetscCall(MatCreateNest(PETSC_COMM_SELF, nr, islrow, nc, islcol, snest, &lA));
    1060:     PetscCall(MatNestSetVecType(lA, VECNEST));
    1061:     for (i = 0; i < nr * nc; i++) {
    1062:       if (istrans[i]) PetscCall(MatDestroy(&snest[i]));
    1063:     }
    1064:     PetscCall(MatISSetLocalMat(B, lA));
    1065:     PetscCall(MatDestroy(&lA));
    1066:     { /* hack : setup of scatters done here */
    1067:       Mat_IS *matis = (Mat_IS *)B->data;

    1069:       matis->islocalref = NULL;
    1070:       PetscCall(MatISSetUpScatters_Private(B));
    1071:     }
    1072:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
    1073:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
    1074:     if (reuse == MAT_INPLACE_MATRIX) {
    1075:       PetscCall(MatHeaderReplace(A, &B));
    1076:     } else {
    1077:       *newmat = B;
    1078:     }
    1079:   } else {
    1080:     if (lreuse) {
    1081:       PetscCall(MatISGetLocalMat(*newmat, &lA));
    1082:       for (i = 0; i < nr; i++) {
    1083:         for (j = 0; j < nc; j++) {
    1084:           if (snest[i * nc + j]) {
    1085:             PetscCall(MatNestSetSubMat(lA, i, j, snest[i * nc + j]));
    1086:             if (istrans[i * nc + j]) PetscCall(MatDestroy(&snest[i * nc + j]));
    1087:           }
    1088:         }
    1089:       }
    1090:     } else {
    1091:       PetscInt stl;
    1092:       for (i = 0, stl = 0; i < nr; i++) {
    1093:         PetscCall(ISCreateStride(PETSC_COMM_SELF, lr[i], stl, 1, &islrow[i]));
    1094:         stl += lr[i];
    1095:       }
    1096:       for (i = 0, stl = 0; i < nc; i++) {
    1097:         PetscCall(ISCreateStride(PETSC_COMM_SELF, lc[i], stl, 1, &islcol[i]));
    1098:         stl += lc[i];
    1099:       }
    1100:       PetscCall(MatCreateNest(PETSC_COMM_SELF, nr, islrow, nc, islcol, snest, &lA));
    1101:       for (i = 0; i < nr * nc; i++) {
    1102:         if (istrans[i]) PetscCall(MatDestroy(&snest[i]));
    1103:       }
    1104:       PetscCall(MatISSetLocalMat(*newmat, lA));
    1105:       PetscCall(MatDestroy(&lA));
    1106:     }
    1107:     PetscCall(MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY));
    1108:     PetscCall(MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY));
    1109:   }

    1111:   /* Create local matrix in MATNEST format */
    1112:   convert = PETSC_FALSE;
    1113:   PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)A)->prefix, "-mat_is_convert_local_nest", &convert, NULL));
    1114:   if (convert) {
    1115:     Mat              M;
    1116:     MatISLocalFields lf;

    1118:     PetscCall(MatISGetLocalMat(*newmat, &lA));
    1119:     PetscCall(MatConvert(lA, MATAIJ, MAT_INITIAL_MATRIX, &M));
    1120:     PetscCall(MatISSetLocalMat(*newmat, M));
    1121:     PetscCall(MatDestroy(&M));

    1123:     /* attach local fields to the matrix */
    1124:     PetscCall(PetscNew(&lf));
    1125:     PetscCall(PetscMalloc2(nr, &lf->rf, nc, &lf->cf));
    1126:     for (i = 0; i < nr; i++) {
    1127:       PetscInt n, st;

    1129:       PetscCall(ISGetLocalSize(islrow[i], &n));
    1130:       PetscCall(ISStrideGetInfo(islrow[i], &st, NULL));
    1131:       PetscCall(ISCreateStride(comm, n, st, 1, &lf->rf[i]));
    1132:     }
    1133:     for (i = 0; i < nc; i++) {
    1134:       PetscInt n, st;

    1136:       PetscCall(ISGetLocalSize(islcol[i], &n));
    1137:       PetscCall(ISStrideGetInfo(islcol[i], &st, NULL));
    1138:       PetscCall(ISCreateStride(comm, n, st, 1, &lf->cf[i]));
    1139:     }
    1140:     lf->nr = nr;
    1141:     lf->nc = nc;
    1142:     PetscCall(PetscObjectContainerCompose((PetscObject)*newmat, "_convert_nest_lfields", lf, MatISContainerDestroyFields_Private));
    1143:   }

    1145:   /* Free workspace */
    1146:   for (i = 0; i < nr; i++) PetscCall(ISDestroy(&islrow[i]));
    1147:   for (i = 0; i < nc; i++) PetscCall(ISDestroy(&islcol[i]));
    1148:   PetscCall(PetscFree6(isrow, iscol, islrow, islcol, snest, istrans));
    1149:   PetscCall(PetscFree2(lr, lc));
    1150:   PetscFunctionReturn(PETSC_SUCCESS);
    1151: }

    1153: static PetscErrorCode MatDiagonalScale_IS(Mat A, Vec l, Vec r)
    1154: {
    1155:   Mat_IS            *matis = (Mat_IS *)A->data;
    1156:   Vec                ll, rr;
    1157:   const PetscScalar *Y, *X;
    1158:   PetscScalar       *x, *y;

    1160:   PetscFunctionBegin;
    1161:   if (l) {
    1162:     ll = matis->y;
    1163:     PetscCall(VecGetArrayRead(l, &Y));
    1164:     PetscCall(VecGetArray(ll, &y));
    1165:     PetscCall(PetscSFBcastBegin(matis->sf, MPIU_SCALAR, Y, y, MPI_REPLACE));
    1166:   } else {
    1167:     ll = NULL;
    1168:   }
    1169:   if (r) {
    1170:     rr = matis->x;
    1171:     PetscCall(VecGetArrayRead(r, &X));
    1172:     PetscCall(VecGetArray(rr, &x));
    1173:     PetscCall(PetscSFBcastBegin(matis->csf, MPIU_SCALAR, X, x, MPI_REPLACE));
    1174:   } else {
    1175:     rr = NULL;
    1176:   }
    1177:   if (ll) {
    1178:     PetscCall(PetscSFBcastEnd(matis->sf, MPIU_SCALAR, Y, y, MPI_REPLACE));
    1179:     PetscCall(VecRestoreArrayRead(l, &Y));
    1180:     PetscCall(VecRestoreArray(ll, &y));
    1181:   }
    1182:   if (rr) {
    1183:     PetscCall(PetscSFBcastEnd(matis->csf, MPIU_SCALAR, X, x, MPI_REPLACE));
    1184:     PetscCall(VecRestoreArrayRead(r, &X));
    1185:     PetscCall(VecRestoreArray(rr, &x));
    1186:   }
    1187:   PetscCall(MatDiagonalScale(matis->A, ll, rr));
    1188:   PetscFunctionReturn(PETSC_SUCCESS);
    1189: }

    1191: static PetscErrorCode MatGetInfo_IS(Mat A, MatInfoType flag, MatInfo *ginfo)
    1192: {
    1193:   Mat_IS        *matis = (Mat_IS *)A->data;
    1194:   MatInfo        info;
    1195:   PetscLogDouble isend[6], irecv[6];
    1196:   PetscInt       bs;

    1198:   PetscFunctionBegin;
    1199:   PetscCall(MatGetBlockSize(A, &bs));
    1200:   if (matis->A->ops->getinfo) {
    1201:     PetscCall(MatGetInfo(matis->A, MAT_LOCAL, &info));
    1202:     isend[0] = info.nz_used;
    1203:     isend[1] = info.nz_allocated;
    1204:     isend[2] = info.nz_unneeded;
    1205:     isend[3] = info.memory;
    1206:     isend[4] = info.mallocs;
    1207:   } else {
    1208:     isend[0] = 0.;
    1209:     isend[1] = 0.;
    1210:     isend[2] = 0.;
    1211:     isend[3] = 0.;
    1212:     isend[4] = 0.;
    1213:   }
    1214:   isend[5] = matis->A->num_ass;
    1215:   if (flag == MAT_LOCAL) {
    1216:     ginfo->nz_used      = isend[0];
    1217:     ginfo->nz_allocated = isend[1];
    1218:     ginfo->nz_unneeded  = isend[2];
    1219:     ginfo->memory       = isend[3];
    1220:     ginfo->mallocs      = isend[4];
    1221:     ginfo->assemblies   = isend[5];
    1222:   } else if (flag == MAT_GLOBAL_MAX) {
    1223:     PetscCallMPI(MPIU_Allreduce(isend, irecv, 6, MPIU_PETSCLOGDOUBLE, MPI_MAX, PetscObjectComm((PetscObject)A)));

    1225:     ginfo->nz_used      = irecv[0];
    1226:     ginfo->nz_allocated = irecv[1];
    1227:     ginfo->nz_unneeded  = irecv[2];
    1228:     ginfo->memory       = irecv[3];
    1229:     ginfo->mallocs      = irecv[4];
    1230:     ginfo->assemblies   = irecv[5];
    1231:   } else if (flag == MAT_GLOBAL_SUM) {
    1232:     PetscCallMPI(MPIU_Allreduce(isend, irecv, 5, MPIU_PETSCLOGDOUBLE, MPI_SUM, PetscObjectComm((PetscObject)A)));

    1234:     ginfo->nz_used      = irecv[0];
    1235:     ginfo->nz_allocated = irecv[1];
    1236:     ginfo->nz_unneeded  = irecv[2];
    1237:     ginfo->memory       = irecv[3];
    1238:     ginfo->mallocs      = irecv[4];
    1239:     ginfo->assemblies   = A->num_ass;
    1240:   }
    1241:   ginfo->block_size        = bs;
    1242:   ginfo->fill_ratio_given  = 0;
    1243:   ginfo->fill_ratio_needed = 0;
    1244:   ginfo->factor_mallocs    = 0;
    1245:   PetscFunctionReturn(PETSC_SUCCESS);
    1246: }

    1248: static PetscErrorCode MatTranspose_IS(Mat A, MatReuse reuse, Mat *B)
    1249: {
    1250:   Mat C, lC, lA;

    1252:   PetscFunctionBegin;
    1253:   if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *B));
    1254:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
    1255:     ISLocalToGlobalMapping rl2g, cl2g;
    1256:     PetscBool              allow_repeated;

    1258:     PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
    1259:     PetscCall(MatSetSizes(C, A->cmap->n, A->rmap->n, A->cmap->N, A->rmap->N));
    1260:     PetscCall(MatSetBlockSizes(C, A->cmap->bs, A->rmap->bs));
    1261:     PetscCall(MatSetType(C, MATIS));
    1262:     PetscCall(MatISGetAllowRepeated(A, &allow_repeated));
    1263:     PetscCall(MatISSetAllowRepeated(C, allow_repeated));
    1264:     PetscCall(MatGetLocalToGlobalMapping(A, &rl2g, &cl2g));
    1265:     PetscCall(MatSetLocalToGlobalMapping(C, cl2g, rl2g));
    1266:   } else C = *B;

    1268:   /* perform local transposition */
    1269:   PetscCall(MatISGetLocalMat(A, &lA));
    1270:   PetscCall(MatTranspose(lA, MAT_INITIAL_MATRIX, &lC));
    1271:   PetscCall(MatSetLocalToGlobalMapping(lC, lA->cmap->mapping, lA->rmap->mapping));
    1272:   PetscCall(MatISSetLocalMat(C, lC));
    1273:   PetscCall(MatDestroy(&lC));

    1275:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
    1276:     *B = C;
    1277:   } else {
    1278:     PetscCall(MatHeaderMerge(A, &C));
    1279:   }
    1280:   PetscCall(MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY));
    1281:   PetscCall(MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY));
    1282:   PetscFunctionReturn(PETSC_SUCCESS);
    1283: }

    1285: static PetscErrorCode MatDiagonalSet_IS(Mat A, Vec D, InsertMode insmode)
    1286: {
    1287:   Mat_IS *is = (Mat_IS *)A->data;

    1289:   PetscFunctionBegin;
    1290:   PetscCheck(!is->allow_repeated || insmode == ADD_VALUES, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "INSERT_VALUES with repeated entries not supported");
    1291:   if (D) { /* MatShift_IS pass D = NULL */
    1292:     PetscCall(VecScatterBegin(is->rctx, D, is->y, INSERT_VALUES, SCATTER_FORWARD));
    1293:     PetscCall(VecScatterEnd(is->rctx, D, is->y, INSERT_VALUES, SCATTER_FORWARD));
    1294:   }
    1295:   PetscCall(VecPointwiseDivide(is->y, is->y, is->counter));
    1296:   PetscCall(MatDiagonalSet(is->A, is->y, insmode));
    1297:   PetscFunctionReturn(PETSC_SUCCESS);
    1298: }

    1300: static PetscErrorCode MatShift_IS(Mat A, PetscScalar a)
    1301: {
    1302:   Mat_IS *is = (Mat_IS *)A->data;

    1304:   PetscFunctionBegin;
    1305:   PetscCall(VecSet(is->y, a));
    1306:   PetscCall(MatDiagonalSet_IS(A, NULL, ADD_VALUES));
    1307:   PetscFunctionReturn(PETSC_SUCCESS);
    1308: }

    1310: static PetscErrorCode MatSetValuesLocal_SubMat_IS(Mat A, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
    1311: {
    1312:   PetscInt buf[2 * MATIS_MAX_ENTRIES_INSERTION], *rows_l = NULL, *cols_l = NULL;

    1314:   PetscFunctionBegin;
    1315:   IndexSpaceGet(buf, m, n, rows_l, cols_l);
    1316:   PetscCall(ISLocalToGlobalMappingApply(A->rmap->mapping, m, rows, rows_l));
    1317:   PetscCall(ISLocalToGlobalMappingApply(A->cmap->mapping, n, cols, cols_l));
    1318:   PetscCall(MatSetValuesLocal_IS(A, m, rows_l, n, cols_l, values, addv));
    1319:   IndexSpaceRestore(buf, m, n, rows_l, cols_l);
    1320:   PetscFunctionReturn(PETSC_SUCCESS);
    1321: }

    1323: static PetscErrorCode MatSetValuesBlockedLocal_SubMat_IS(Mat A, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
    1324: {
    1325:   PetscInt buf[2 * MATIS_MAX_ENTRIES_INSERTION], *rows_l = NULL, *cols_l = NULL, rbs, cbs;

    1327:   PetscFunctionBegin;
    1328:   /* We cannot guarantee the local matrix will have the same block size of the original matrix */
    1329:   PetscCall(ISLocalToGlobalMappingGetBlockSize(A->rmap->mapping, &rbs));
    1330:   PetscCall(ISLocalToGlobalMappingGetBlockSize(A->cmap->mapping, &cbs));
    1331:   IndexSpaceGet(buf, m * rbs, n * cbs, rows_l, cols_l);
    1332:   BlockIndicesExpand(m, rows, rbs, rows_l);
    1333:   BlockIndicesExpand(n, cols, cbs, cols_l);
    1334:   PetscCall(ISLocalToGlobalMappingApply(A->rmap->mapping, m * rbs, rows_l, rows_l));
    1335:   PetscCall(ISLocalToGlobalMappingApply(A->cmap->mapping, n * cbs, cols_l, cols_l));
    1336:   PetscCall(MatSetValuesLocal_IS(A, m * rbs, rows_l, n * cbs, cols_l, values, addv));
    1337:   IndexSpaceRestore(buf, m * rbs, n * cbs, rows_l, cols_l);
    1338:   PetscFunctionReturn(PETSC_SUCCESS);
    1339: }

    1341: static PetscErrorCode MatZeroRowsLocal_SubMat_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
    1342: {
    1343:   PetscInt *rows_l;
    1344:   Mat_IS   *is = (Mat_IS *)A->data;

    1346:   PetscFunctionBegin;
    1347:   PetscCall(PetscMalloc1(n, &rows_l));
    1348:   PetscCall(ISLocalToGlobalMappingApply(A->rmap->mapping, n, rows, rows_l));
    1349:   PetscCall(MatZeroRowsLocal(is->islocalref, n, rows_l, diag, x, b));
    1350:   PetscCall(PetscFree(rows_l));
    1351:   PetscFunctionReturn(PETSC_SUCCESS);
    1352: }

    1354: static PetscErrorCode MatZeroRowsColumnsLocal_SubMat_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
    1355: {
    1356:   PetscInt *rows_l;
    1357:   Mat_IS   *is = (Mat_IS *)A->data;

    1359:   PetscFunctionBegin;
    1360:   PetscCall(PetscMalloc1(n, &rows_l));
    1361:   PetscCall(ISLocalToGlobalMappingApply(A->rmap->mapping, n, rows, rows_l));
    1362:   PetscCall(MatZeroRowsColumnsLocal(is->islocalref, n, rows_l, diag, x, b));
    1363:   PetscCall(PetscFree(rows_l));
    1364:   PetscFunctionReturn(PETSC_SUCCESS);
    1365: }

    1367: static PetscErrorCode MatCreateSubMatrix_IS(Mat mat, IS irow, IS icol, MatReuse scall, Mat *newmat)
    1368: {
    1369:   Mat             locmat, newlocmat;
    1370:   Mat_IS         *newmatis;
    1371:   const PetscInt *idxs;
    1372:   PetscInt        i, m, n;

    1374:   PetscFunctionBegin;
    1375:   if (scall == MAT_REUSE_MATRIX) {
    1376:     PetscBool ismatis;

    1378:     PetscCall(PetscObjectTypeCompare((PetscObject)*newmat, MATIS, &ismatis));
    1379:     PetscCheck(ismatis, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_ARG_WRONG, "Cannot reuse matrix! Not of MATIS type");
    1380:     newmatis = (Mat_IS *)(*newmat)->data;
    1381:     PetscCheck(newmatis->getsub_ris, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_ARG_WRONG, "Cannot reuse matrix! Misses local row IS");
    1382:     PetscCheck(newmatis->getsub_cis, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_ARG_WRONG, "Cannot reuse matrix! Misses local col IS");
    1383:   }
    1384:   /* irow and icol may not have duplicate entries */
    1385:   if (PetscDefined(USE_DEBUG)) {
    1386:     Vec                rtest, ltest;
    1387:     const PetscScalar *array;

    1389:     PetscCall(MatCreateVecs(mat, &ltest, &rtest));
    1390:     PetscCall(ISGetLocalSize(irow, &n));
    1391:     PetscCall(ISGetIndices(irow, &idxs));
    1392:     for (i = 0; i < n; i++) PetscCall(VecSetValue(rtest, idxs[i], 1.0, ADD_VALUES));
    1393:     PetscCall(VecAssemblyBegin(rtest));
    1394:     PetscCall(VecAssemblyEnd(rtest));
    1395:     PetscCall(VecGetLocalSize(rtest, &n));
    1396:     PetscCall(VecGetOwnershipRange(rtest, &m, NULL));
    1397:     PetscCall(VecGetArrayRead(rtest, &array));
    1398:     for (i = 0; i < n; i++) PetscCheck(array[i] == 0. || array[i] == 1., PETSC_COMM_SELF, PETSC_ERR_SUP, "Index %" PetscInt_FMT " counted %" PetscInt_FMT " times! Irow may not have duplicate entries", i + m, (PetscInt)PetscRealPart(array[i]));
    1399:     PetscCall(VecRestoreArrayRead(rtest, &array));
    1400:     PetscCall(ISRestoreIndices(irow, &idxs));
    1401:     PetscCall(ISGetLocalSize(icol, &n));
    1402:     PetscCall(ISGetIndices(icol, &idxs));
    1403:     for (i = 0; i < n; i++) PetscCall(VecSetValue(ltest, idxs[i], 1.0, ADD_VALUES));
    1404:     PetscCall(VecAssemblyBegin(ltest));
    1405:     PetscCall(VecAssemblyEnd(ltest));
    1406:     PetscCall(VecGetLocalSize(ltest, &n));
    1407:     PetscCall(VecGetOwnershipRange(ltest, &m, NULL));
    1408:     PetscCall(VecGetArrayRead(ltest, &array));
    1409:     for (i = 0; i < n; i++) PetscCheck(array[i] == 0. || array[i] == 1., PETSC_COMM_SELF, PETSC_ERR_SUP, "Index %" PetscInt_FMT " counted %" PetscInt_FMT " times! Icol may not have duplicate entries", i + m, (PetscInt)PetscRealPart(array[i]));
    1410:     PetscCall(VecRestoreArrayRead(ltest, &array));
    1411:     PetscCall(ISRestoreIndices(icol, &idxs));
    1412:     PetscCall(VecDestroy(&rtest));
    1413:     PetscCall(VecDestroy(&ltest));
    1414:   }
    1415:   if (scall == MAT_INITIAL_MATRIX) {
    1416:     Mat_IS                *matis = (Mat_IS *)mat->data;
    1417:     ISLocalToGlobalMapping rl2g;
    1418:     IS                     is;
    1419:     PetscInt              *lidxs, *lgidxs, *newgidxs;
    1420:     PetscInt               ll, newloc, irbs, icbs, arbs, acbs, rbs, cbs;
    1421:     PetscBool              cong;
    1422:     MPI_Comm               comm;

    1424:     PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
    1425:     PetscCall(MatGetBlockSizes(mat, &arbs, &acbs));
    1426:     PetscCall(ISGetBlockSize(irow, &irbs));
    1427:     PetscCall(ISGetBlockSize(icol, &icbs));
    1428:     rbs = arbs == irbs ? irbs : 1;
    1429:     cbs = acbs == icbs ? icbs : 1;
    1430:     PetscCall(ISGetLocalSize(irow, &m));
    1431:     PetscCall(ISGetLocalSize(icol, &n));
    1432:     PetscCall(MatCreate(comm, newmat));
    1433:     PetscCall(MatSetType(*newmat, MATIS));
    1434:     PetscCall(MatISSetAllowRepeated(*newmat, matis->allow_repeated));
    1435:     PetscCall(MatSetSizes(*newmat, m, n, PETSC_DECIDE, PETSC_DECIDE));
    1436:     PetscCall(MatSetBlockSizes(*newmat, rbs, cbs));
    1437:     /* communicate irow to their owners in the layout */
    1438:     PetscCall(ISGetIndices(irow, &idxs));
    1439:     PetscCall(PetscLayoutMapLocal(mat->rmap, m, idxs, &ll, &lidxs, &lgidxs));
    1440:     PetscCall(ISRestoreIndices(irow, &idxs));
    1441:     PetscCall(PetscArrayzero(matis->sf_rootdata, matis->sf->nroots));
    1442:     for (i = 0; i < ll; i++) matis->sf_rootdata[lidxs[i]] = lgidxs[i] + 1;
    1443:     PetscCall(PetscFree(lidxs));
    1444:     PetscCall(PetscFree(lgidxs));
    1445:     PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
    1446:     PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
    1447:     for (i = 0, newloc = 0; i < matis->sf->nleaves; i++)
    1448:       if (matis->sf_leafdata[i]) newloc++;
    1449:     PetscCall(PetscMalloc1(newloc, &newgidxs));
    1450:     PetscCall(PetscMalloc1(newloc, &lidxs));
    1451:     for (i = 0, newloc = 0; i < matis->sf->nleaves; i++)
    1452:       if (matis->sf_leafdata[i]) {
    1453:         lidxs[newloc]      = i;
    1454:         newgidxs[newloc++] = matis->sf_leafdata[i] - 1;
    1455:       }
    1456:     PetscCall(ISCreateGeneral(comm, newloc, newgidxs, PETSC_OWN_POINTER, &is));
    1457:     PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g));
    1458:     PetscCall(ISLocalToGlobalMappingSetBlockSize(rl2g, rbs));
    1459:     PetscCall(ISDestroy(&is));
    1460:     /* local is to extract local submatrix */
    1461:     newmatis = (Mat_IS *)(*newmat)->data;
    1462:     PetscCall(ISCreateGeneral(comm, newloc, lidxs, PETSC_OWN_POINTER, &newmatis->getsub_ris));
    1463:     PetscCall(MatHasCongruentLayouts(mat, &cong));
    1464:     if (cong && irow == icol && matis->csf == matis->sf) {
    1465:       PetscCall(MatSetLocalToGlobalMapping(*newmat, rl2g, rl2g));
    1466:       PetscCall(PetscObjectReference((PetscObject)newmatis->getsub_ris));
    1467:       newmatis->getsub_cis = newmatis->getsub_ris;
    1468:     } else {
    1469:       ISLocalToGlobalMapping cl2g;

    1471:       /* communicate icol to their owners in the layout */
    1472:       PetscCall(ISGetIndices(icol, &idxs));
    1473:       PetscCall(PetscLayoutMapLocal(mat->cmap, n, idxs, &ll, &lidxs, &lgidxs));
    1474:       PetscCall(ISRestoreIndices(icol, &idxs));
    1475:       PetscCall(PetscArrayzero(matis->csf_rootdata, matis->csf->nroots));
    1476:       for (i = 0; i < ll; i++) matis->csf_rootdata[lidxs[i]] = lgidxs[i] + 1;
    1477:       PetscCall(PetscFree(lidxs));
    1478:       PetscCall(PetscFree(lgidxs));
    1479:       PetscCall(PetscSFBcastBegin(matis->csf, MPIU_INT, matis->csf_rootdata, matis->csf_leafdata, MPI_REPLACE));
    1480:       PetscCall(PetscSFBcastEnd(matis->csf, MPIU_INT, matis->csf_rootdata, matis->csf_leafdata, MPI_REPLACE));
    1481:       for (i = 0, newloc = 0; i < matis->csf->nleaves; i++)
    1482:         if (matis->csf_leafdata[i]) newloc++;
    1483:       PetscCall(PetscMalloc1(newloc, &newgidxs));
    1484:       PetscCall(PetscMalloc1(newloc, &lidxs));
    1485:       for (i = 0, newloc = 0; i < matis->csf->nleaves; i++)
    1486:         if (matis->csf_leafdata[i]) {
    1487:           lidxs[newloc]      = i;
    1488:           newgidxs[newloc++] = matis->csf_leafdata[i] - 1;
    1489:         }
    1490:       PetscCall(ISCreateGeneral(comm, newloc, newgidxs, PETSC_OWN_POINTER, &is));
    1491:       PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g));
    1492:       PetscCall(ISLocalToGlobalMappingSetBlockSize(cl2g, cbs));
    1493:       PetscCall(ISDestroy(&is));
    1494:       /* local is to extract local submatrix */
    1495:       PetscCall(ISCreateGeneral(comm, newloc, lidxs, PETSC_OWN_POINTER, &newmatis->getsub_cis));
    1496:       PetscCall(MatSetLocalToGlobalMapping(*newmat, rl2g, cl2g));
    1497:       PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
    1498:     }
    1499:     PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
    1500:   } else {
    1501:     PetscCall(MatISGetLocalMat(*newmat, &newlocmat));
    1502:   }
    1503:   PetscCall(MatISGetLocalMat(mat, &locmat));
    1504:   newmatis = (Mat_IS *)(*newmat)->data;
    1505:   PetscCall(MatCreateSubMatrix(locmat, newmatis->getsub_ris, newmatis->getsub_cis, scall, &newlocmat));
    1506:   if (scall == MAT_INITIAL_MATRIX) {
    1507:     PetscCall(MatISSetLocalMat(*newmat, newlocmat));
    1508:     PetscCall(MatDestroy(&newlocmat));
    1509:   }
    1510:   PetscCall(MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY));
    1511:   PetscCall(MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY));
    1512:   PetscFunctionReturn(PETSC_SUCCESS);
    1513: }

    1515: static PetscErrorCode MatCopy_IS(Mat A, Mat B, MatStructure str)
    1516: {
    1517:   Mat_IS   *a = (Mat_IS *)A->data, *b;
    1518:   PetscBool ismatis;

    1520:   PetscFunctionBegin;
    1521:   PetscCall(PetscObjectTypeCompare((PetscObject)B, MATIS, &ismatis));
    1522:   PetscCheck(ismatis, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Need to be implemented");
    1523:   b = (Mat_IS *)B->data;
    1524:   PetscCall(MatCopy(a->A, b->A, str));
    1525:   PetscCall(PetscObjectStateIncrease((PetscObject)B));
    1526:   PetscFunctionReturn(PETSC_SUCCESS);
    1527: }

    1529: static PetscErrorCode MatISSetUpSF_IS(Mat B)
    1530: {
    1531:   Mat_IS         *matis = (Mat_IS *)B->data;
    1532:   const PetscInt *gidxs;
    1533:   PetscInt        nleaves;

    1535:   PetscFunctionBegin;
    1536:   if (matis->sf) PetscFunctionReturn(PETSC_SUCCESS);
    1537:   PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)B), &matis->sf));
    1538:   PetscCall(ISLocalToGlobalMappingGetIndices(matis->rmapping, &gidxs));
    1539:   PetscCall(ISLocalToGlobalMappingGetSize(matis->rmapping, &nleaves));
    1540:   PetscCall(PetscSFSetGraphLayout(matis->sf, B->rmap, nleaves, NULL, PETSC_OWN_POINTER, gidxs));
    1541:   PetscCall(ISLocalToGlobalMappingRestoreIndices(matis->rmapping, &gidxs));
    1542:   PetscCall(PetscMalloc2(matis->sf->nroots, &matis->sf_rootdata, matis->sf->nleaves, &matis->sf_leafdata));
    1543:   if (matis->rmapping != matis->cmapping) { /* setup SF for columns */
    1544:     PetscCall(ISLocalToGlobalMappingGetSize(matis->cmapping, &nleaves));
    1545:     PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)B), &matis->csf));
    1546:     PetscCall(ISLocalToGlobalMappingGetIndices(matis->cmapping, &gidxs));
    1547:     PetscCall(PetscSFSetGraphLayout(matis->csf, B->cmap, nleaves, NULL, PETSC_OWN_POINTER, gidxs));
    1548:     PetscCall(ISLocalToGlobalMappingRestoreIndices(matis->cmapping, &gidxs));
    1549:     PetscCall(PetscMalloc2(matis->csf->nroots, &matis->csf_rootdata, matis->csf->nleaves, &matis->csf_leafdata));
    1550:   } else {
    1551:     matis->csf          = matis->sf;
    1552:     matis->csf_leafdata = matis->sf_leafdata;
    1553:     matis->csf_rootdata = matis->sf_rootdata;
    1554:   }
    1555:   PetscFunctionReturn(PETSC_SUCCESS);
    1556: }

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

    1561:   Not Collective

    1563:   Input Parameter:
    1564: . A - the matrix

    1566:   Output Parameter:
    1567: . flg - the boolean flag

    1569:   Level: intermediate

    1571: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateIS()`, `MatSetLocalToGlobalMapping()`, `MatISSetAllowRepeated()`
    1572: @*/
    1573: PetscErrorCode MatISGetAllowRepeated(Mat A, PetscBool *flg)
    1574: {
    1575:   PetscBool ismatis;

    1577:   PetscFunctionBegin;
    1579:   PetscAssertPointer(flg, 2);
    1580:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATIS, &ismatis));
    1581:   PetscCheck(ismatis, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not for matrix type %s", ((PetscObject)A)->type_name);
    1582:   *flg = ((Mat_IS *)A->data)->allow_repeated;
    1583:   PetscFunctionReturn(PETSC_SUCCESS);
    1584: }

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

    1589:   Logically Collective

    1591:   Input Parameters:
    1592: + A   - the matrix
    1593: - flg - the boolean flag

    1595:   Level: intermediate

    1597:   Notes:
    1598:   The default value is `PETSC_FALSE`.
    1599:   When called AFTER calling `MatSetLocalToGlobalMapping()` it will recreate the local matrices
    1600:   if `flg` is different from the previously set value.
    1601:   Specifically, when `flg` is true it will just recreate the local matrices, while if
    1602:   `flg` is false will assemble the local matrices summing up repeated entries.

    1604: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateIS()`, `MatSetLocalToGlobalMapping()`, `MatISGetAllowRepeated()`
    1605: @*/
    1606: PetscErrorCode MatISSetAllowRepeated(Mat A, PetscBool flg)
    1607: {
    1608:   PetscFunctionBegin;
    1612:   PetscTryMethod(A, "MatISSetAllowRepeated_C", (Mat, PetscBool), (A, flg));
    1613:   PetscFunctionReturn(PETSC_SUCCESS);
    1614: }

    1616: static PetscErrorCode MatISSetAllowRepeated_IS(Mat A, PetscBool flg)
    1617: {
    1618:   Mat_IS                *matis = (Mat_IS *)A->data;
    1619:   Mat                    lA    = NULL;
    1620:   ISLocalToGlobalMapping lrmap, lcmap;

    1622:   PetscFunctionBegin;
    1623:   if (flg == matis->allow_repeated) PetscFunctionReturn(PETSC_SUCCESS);
    1624:   if (!matis->A) { /* matrix has not been preallocated yet */
    1625:     matis->allow_repeated = flg;
    1626:     PetscFunctionReturn(PETSC_SUCCESS);
    1627:   }
    1628:   PetscCheck(!matis->islocalref, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not implemented for local references");
    1629:   if (matis->allow_repeated) { /* we will assemble the old local matrix if needed */
    1630:     lA = matis->A;
    1631:     PetscCall(PetscObjectReference((PetscObject)lA));
    1632:   }
    1633:   /* In case flg is True, we only recreate the local matrix */
    1634:   matis->allow_repeated = flg;
    1635:   PetscCall(MatSetLocalToGlobalMapping(A, A->rmap->mapping, A->cmap->mapping));
    1636:   if (lA) { /* assemble previous local matrix if needed */
    1637:     Mat nA = matis->A;

    1639:     PetscCall(MatGetLocalToGlobalMapping(nA, &lrmap, &lcmap));
    1640:     if (!lrmap && !lcmap) {
    1641:       PetscCall(MatISSetLocalMat(A, lA));
    1642:     } else {
    1643:       Mat            P = NULL, R = NULL;
    1644:       MatProductType ptype;

    1646:       if (lrmap == lcmap) {
    1647:         ptype = MATPRODUCT_PtAP;
    1648:         PetscCall(MatCreateFromISLocalToGlobalMapping(lcmap, nA, PETSC_TRUE, PETSC_FALSE, NULL, &P));
    1649:         PetscCall(MatProductCreate(lA, P, NULL, &nA));
    1650:       } else {
    1651:         if (lcmap) PetscCall(MatCreateFromISLocalToGlobalMapping(lcmap, nA, PETSC_TRUE, PETSC_FALSE, NULL, &P));
    1652:         if (lrmap) PetscCall(MatCreateFromISLocalToGlobalMapping(lrmap, nA, PETSC_FALSE, PETSC_TRUE, NULL, &R));
    1653:         if (R && P) {
    1654:           ptype = MATPRODUCT_ABC;
    1655:           PetscCall(MatProductCreate(R, lA, P, &nA));
    1656:         } else if (R) {
    1657:           ptype = MATPRODUCT_AB;
    1658:           PetscCall(MatProductCreate(R, lA, NULL, &nA));
    1659:         } else {
    1660:           ptype = MATPRODUCT_AB;
    1661:           PetscCall(MatProductCreate(lA, P, NULL, &nA));
    1662:         }
    1663:       }
    1664:       PetscCall(MatProductSetType(nA, ptype));
    1665:       PetscCall(MatProductSetFromOptions(nA));
    1666:       PetscCall(MatProductSymbolic(nA));
    1667:       PetscCall(MatProductNumeric(nA));
    1668:       PetscCall(MatProductClear(nA));
    1669:       PetscCall(MatConvert(nA, matis->lmattype, MAT_INPLACE_MATRIX, &nA));
    1670:       PetscCall(MatISSetLocalMat(A, nA));
    1671:       PetscCall(MatDestroy(&nA));
    1672:       PetscCall(MatDestroy(&P));
    1673:       PetscCall(MatDestroy(&R));
    1674:     }
    1675:   }
    1676:   PetscCall(MatDestroy(&lA));
    1677:   PetscFunctionReturn(PETSC_SUCCESS);
    1678: }

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

    1683:   Logically Collective

    1685:   Input Parameters:
    1686: + A     - the matrix
    1687: - store - the boolean flag

    1689:   Level: advanced

    1691: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateIS()`, `MatISSetPreallocation()`, `MatPtAP()`
    1692: @*/
    1693: PetscErrorCode MatISStoreL2L(Mat A, PetscBool store)
    1694: {
    1695:   PetscFunctionBegin;
    1699:   PetscTryMethod(A, "MatISStoreL2L_C", (Mat, PetscBool), (A, store));
    1700:   PetscFunctionReturn(PETSC_SUCCESS);
    1701: }

    1703: static PetscErrorCode MatISStoreL2L_IS(Mat A, PetscBool store)
    1704: {
    1705:   Mat_IS *matis = (Mat_IS *)A->data;

    1707:   PetscFunctionBegin;
    1708:   matis->storel2l = store;
    1709:   if (!store) PetscCall(PetscObjectCompose((PetscObject)A, "_MatIS_PtAP_l2l", NULL));
    1710:   PetscFunctionReturn(PETSC_SUCCESS);
    1711: }

    1713: /*@
    1714:   MatISFixLocalEmpty - Compress out zero local rows from the local matrices

    1716:   Logically Collective

    1718:   Input Parameters:
    1719: + A   - the matrix
    1720: - fix - the boolean flag

    1722:   Level: advanced

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

    1727: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatCreate()`, `MatCreateIS()`, `MatISSetPreallocation()`, `MatAssemblyEnd()`, `MAT_FINAL_ASSEMBLY`
    1728: @*/
    1729: PetscErrorCode MatISFixLocalEmpty(Mat A, PetscBool fix)
    1730: {
    1731:   PetscFunctionBegin;
    1735:   PetscTryMethod(A, "MatISFixLocalEmpty_C", (Mat, PetscBool), (A, fix));
    1736:   PetscFunctionReturn(PETSC_SUCCESS);
    1737: }

    1739: static PetscErrorCode MatISFixLocalEmpty_IS(Mat A, PetscBool fix)
    1740: {
    1741:   Mat_IS *matis = (Mat_IS *)A->data;

    1743:   PetscFunctionBegin;
    1744:   matis->locempty = fix;
    1745:   PetscFunctionReturn(PETSC_SUCCESS);
    1746: }

    1748: /*@
    1749:   MatISSetPreallocation - Preallocates memory for a `MATIS` parallel matrix.

    1751:   Collective

    1753:   Input Parameters:
    1754: + B     - the matrix
    1755: . d_nz  - number of nonzeros per row in DIAGONAL portion of local submatrix
    1756:            (same value is used for all local rows)
    1757: . d_nnz - array containing the number of nonzeros in the various rows of the
    1758:            DIAGONAL portion of the local submatrix (possibly different for each row)
    1759:            or `NULL`, if `d_nz` is used to specify the nonzero structure.
    1760:            The size of this array is equal to the number of local rows, i.e `m`.
    1761:            For matrices that will be factored, you must leave room for (and set)
    1762:            the diagonal entry even if it is zero.
    1763: . o_nz  - number of nonzeros per row in the OFF-DIAGONAL portion of local
    1764:            submatrix (same value is used for all local rows).
    1765: - o_nnz - array containing the number of nonzeros in the various rows of the
    1766:            OFF-DIAGONAL portion of the local submatrix (possibly different for
    1767:            each row) or `NULL`, if `o_nz` is used to specify the nonzero
    1768:            structure. The size of this array is equal to the number
    1769:            of local rows, i.e `m`.

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

    1773:   Level: intermediate

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

    1780: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateIS()`, `MatMPIAIJSetPreallocation()`, `MatISGetLocalMat()`, `MATIS`
    1781: @*/
    1782: PetscErrorCode MatISSetPreallocation(Mat B, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[])
    1783: {
    1784:   PetscFunctionBegin;
    1787:   PetscTryMethod(B, "MatISSetPreallocation_C", (Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[]), (B, d_nz, d_nnz, o_nz, o_nnz));
    1788:   PetscFunctionReturn(PETSC_SUCCESS);
    1789: }

    1791: static PetscErrorCode MatISSetPreallocation_IS(Mat B, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[])
    1792: {
    1793:   Mat_IS  *matis = (Mat_IS *)B->data;
    1794:   PetscInt bs, i, nlocalcols;

    1796:   PetscFunctionBegin;
    1797:   PetscCall(MatSetUp(B));
    1798:   if (!d_nnz)
    1799:     for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] = d_nz;
    1800:   else
    1801:     for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] = d_nnz[i];

    1803:   if (!o_nnz)
    1804:     for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] += o_nz;
    1805:   else
    1806:     for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] += o_nnz[i];

    1808:   PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
    1809:   PetscCall(MatGetSize(matis->A, NULL, &nlocalcols));
    1810:   PetscCall(MatGetBlockSize(matis->A, &bs));
    1811:   PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));

    1813:   for (i = 0; i < matis->sf->nleaves; i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i], nlocalcols);
    1814:   PetscCall(MatSeqAIJSetPreallocation(matis->A, 0, matis->sf_leafdata));
    1815: #if defined(PETSC_HAVE_HYPRE)
    1816:   PetscCall(MatHYPRESetPreallocation(matis->A, 0, matis->sf_leafdata, 0, NULL));
    1817: #endif

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

    1822:     matis->sf_leafdata[i] = matis->sf_leafdata[i * bs] / bs;
    1823:     for (b = 1; b < bs; b++) matis->sf_leafdata[i] = PetscMax(matis->sf_leafdata[i], matis->sf_leafdata[i * bs + b] / bs);
    1824:   }
    1825:   PetscCall(MatSeqBAIJSetPreallocation(matis->A, bs, 0, matis->sf_leafdata));

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

    1831:   /* for other matrix types */
    1832:   PetscCall(MatSetUp(matis->A));
    1833:   PetscFunctionReturn(PETSC_SUCCESS);
    1834: }

    1836: PETSC_INTERN PetscErrorCode MatConvert_IS_XAIJ(Mat mat, MatType mtype, MatReuse reuse, Mat *M)
    1837: {
    1838:   Mat_IS            *matis     = (Mat_IS *)mat->data;
    1839:   Mat                local_mat = NULL, MT;
    1840:   PetscInt           rbs, cbs, rows, cols, lrows, lcols;
    1841:   PetscInt           local_rows, local_cols;
    1842:   PetscBool          isseqdense, isseqsbaij, isseqaij, isseqbaij;
    1843:   PetscMPIInt        size;
    1844:   const PetscScalar *array;

    1846:   PetscFunctionBegin;
    1847:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)mat), &size));
    1848:   if (size == 1 && mat->rmap->N == matis->A->rmap->N && mat->cmap->N == matis->A->cmap->N && !matis->allow_repeated) {
    1849:     Mat      B;
    1850:     IS       irows = NULL, icols = NULL;
    1851:     PetscInt rbs, cbs;

    1853:     PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->rmapping, &rbs));
    1854:     PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->cmapping, &cbs));
    1855:     if (reuse != MAT_REUSE_MATRIX) { /* check if l2g maps are one-to-one */
    1856:       IS              rows, cols;
    1857:       const PetscInt *ridxs, *cidxs;
    1858:       PetscInt        i, nw;
    1859:       PetscBT         work;

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

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

    1937:     lb[0] = isseqdense;
    1938:     lb[1] = isseqaij;
    1939:     lb[2] = isseqbaij;
    1940:     lb[3] = isseqsbaij;
    1941:     PetscCallMPI(MPIU_Allreduce(lb, bb, 4, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)mat)));
    1942:     PetscCheck(bb[0] || bb[1] || bb[2] || bb[3], PETSC_COMM_SELF, PETSC_ERR_SUP, "Local matrices must have the same type");
    1943:   }

    1945:   if (reuse != MAT_REUSE_MATRIX) {
    1946:     PetscCount ncoo;
    1947:     PetscInt  *coo_i, *coo_j;

    1949:     PetscCall(MatCreate(PetscObjectComm((PetscObject)mat), &MT));
    1950:     PetscCall(MatSetSizes(MT, lrows, lcols, rows, cols));
    1951:     PetscCall(MatSetType(MT, mtype));
    1952:     PetscCall(MatSetBlockSizes(MT, rbs, cbs));
    1953:     if (!isseqaij && !isseqdense) {
    1954:       PetscCall(MatConvert(matis->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &local_mat));
    1955:     } else {
    1956:       PetscCall(PetscObjectReference((PetscObject)matis->A));
    1957:       local_mat = matis->A;
    1958:     }
    1959:     PetscCall(MatSetLocalToGlobalMapping(MT, matis->rmapping, matis->cmapping));
    1960:     if (isseqdense) {
    1961:       PetscInt nr, nc;

    1963:       PetscCall(MatGetSize(local_mat, &nr, &nc));
    1964:       ncoo = nr * nc;
    1965:       PetscCall(PetscMalloc2(ncoo, &coo_i, ncoo, &coo_j));
    1966:       for (PetscInt j = 0; j < nc; j++) {
    1967:         for (PetscInt i = 0; i < nr; i++) {
    1968:           coo_i[j * nr + i] = i;
    1969:           coo_j[j * nr + i] = j;
    1970:         }
    1971:       }
    1972:     } else {
    1973:       const PetscInt *ii, *jj;
    1974:       PetscInt        nr;
    1975:       PetscBool       done;

    1977:       PetscCall(MatGetRowIJ(local_mat, 0, PETSC_FALSE, PETSC_FALSE, &nr, &ii, &jj, &done));
    1978:       PetscCheck(done, PetscObjectComm((PetscObject)local_mat), PETSC_ERR_PLIB, "Error in MatGetRowIJ");
    1979:       ncoo = ii[nr];
    1980:       PetscCall(PetscMalloc2(ncoo, &coo_i, ncoo, &coo_j));
    1981:       PetscCall(PetscArraycpy(coo_j, jj, ncoo));
    1982:       for (PetscInt i = 0; i < nr; i++) {
    1983:         for (PetscInt j = ii[i]; j < ii[i + 1]; j++) coo_i[j] = i;
    1984:       }
    1985:       PetscCall(MatRestoreRowIJ(local_mat, 0, PETSC_FALSE, PETSC_FALSE, &nr, &ii, &jj, &done));
    1986:       PetscCheck(done, PetscObjectComm((PetscObject)local_mat), PETSC_ERR_PLIB, "Error in MatRestoreRowIJ");
    1987:     }
    1988:     PetscCall(MatSetPreallocationCOOLocal(MT, ncoo, coo_i, coo_j));
    1989:     PetscCall(PetscFree2(coo_i, coo_j));
    1990:   } else {
    1991:     PetscInt mrbs, mcbs, mrows, mcols, mlrows, mlcols;

    1993:     /* some checks */
    1994:     MT = *M;
    1995:     PetscCall(MatGetBlockSizes(MT, &mrbs, &mcbs));
    1996:     PetscCall(MatGetSize(MT, &mrows, &mcols));
    1997:     PetscCall(MatGetLocalSize(MT, &mlrows, &mlcols));
    1998:     PetscCheck(mrows == rows, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of rows (%" PetscInt_FMT " != %" PetscInt_FMT ")", rows, mrows);
    1999:     PetscCheck(mcols == cols, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of cols (%" PetscInt_FMT " != %" PetscInt_FMT ")", cols, mcols);
    2000:     PetscCheck(mlrows == lrows, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of local rows (%" PetscInt_FMT " != %" PetscInt_FMT ")", lrows, mlrows);
    2001:     PetscCheck(mlcols == lcols, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of local cols (%" PetscInt_FMT " != %" PetscInt_FMT ")", lcols, mlcols);
    2002:     PetscCheck(mrbs == rbs, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong row block size (%" PetscInt_FMT " != %" PetscInt_FMT ")", rbs, mrbs);
    2003:     PetscCheck(mcbs == cbs, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong col block size (%" PetscInt_FMT " != %" PetscInt_FMT ")", cbs, mcbs);
    2004:     PetscCall(MatZeroEntries(MT));
    2005:     if (!isseqaij && !isseqdense) {
    2006:       PetscCall(MatConvert(matis->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &local_mat));
    2007:     } else {
    2008:       PetscCall(PetscObjectReference((PetscObject)matis->A));
    2009:       local_mat = matis->A;
    2010:     }
    2011:   }

    2013:   /* Set values */
    2014:   if (isseqdense) {
    2015:     PetscCall(MatDenseGetArrayRead(local_mat, &array));
    2016:     PetscCall(MatSetValuesCOO(MT, array, INSERT_VALUES));
    2017:     PetscCall(MatDenseRestoreArrayRead(local_mat, &array));
    2018:   } else {
    2019:     PetscCall(MatSeqAIJGetArrayRead(local_mat, &array));
    2020:     PetscCall(MatSetValuesCOO(MT, array, INSERT_VALUES));
    2021:     PetscCall(MatSeqAIJRestoreArrayRead(local_mat, &array));
    2022:   }
    2023:   PetscCall(MatDestroy(&local_mat));
    2024:   PetscCall(MatAssemblyBegin(MT, MAT_FINAL_ASSEMBLY));
    2025:   PetscCall(MatAssemblyEnd(MT, MAT_FINAL_ASSEMBLY));
    2026:   if (reuse == MAT_INPLACE_MATRIX) {
    2027:     PetscCall(MatHeaderReplace(mat, &MT));
    2028:   } else if (reuse == MAT_INITIAL_MATRIX) {
    2029:     *M = MT;
    2030:   }
    2031:   PetscFunctionReturn(PETSC_SUCCESS);
    2032: }

    2034: static PetscErrorCode MatDuplicate_IS(Mat mat, MatDuplicateOption op, Mat *newmat)
    2035: {
    2036:   Mat_IS  *matis = (Mat_IS *)mat->data;
    2037:   PetscInt rbs, cbs, m, n, M, N;
    2038:   Mat      B, localmat;

    2040:   PetscFunctionBegin;
    2041:   PetscCall(ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping, &rbs));
    2042:   PetscCall(ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping, &cbs));
    2043:   PetscCall(MatGetSize(mat, &M, &N));
    2044:   PetscCall(MatGetLocalSize(mat, &m, &n));
    2045:   PetscCall(MatCreate(PetscObjectComm((PetscObject)mat), &B));
    2046:   PetscCall(MatSetSizes(B, m, n, M, N));
    2047:   PetscCall(MatSetBlockSize(B, rbs == cbs ? rbs : 1));
    2048:   PetscCall(MatSetType(B, MATIS));
    2049:   PetscCall(MatISSetLocalMatType(B, matis->lmattype));
    2050:   PetscCall(MatISSetAllowRepeated(B, matis->allow_repeated));
    2051:   PetscCall(MatSetLocalToGlobalMapping(B, mat->rmap->mapping, mat->cmap->mapping));
    2052:   PetscCall(MatDuplicate(matis->A, op, &localmat));
    2053:   PetscCall(MatSetLocalToGlobalMapping(localmat, matis->A->rmap->mapping, matis->A->cmap->mapping));
    2054:   PetscCall(MatISSetLocalMat(B, localmat));
    2055:   PetscCall(MatDestroy(&localmat));
    2056:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
    2057:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
    2058:   *newmat = B;
    2059:   PetscFunctionReturn(PETSC_SUCCESS);
    2060: }

    2062: static PetscErrorCode MatIsHermitian_IS(Mat A, PetscReal tol, PetscBool *flg)
    2063: {
    2064:   Mat_IS   *matis = (Mat_IS *)A->data;
    2065:   PetscBool local_sym;

    2067:   PetscFunctionBegin;
    2068:   PetscCall(MatIsHermitian(matis->A, tol, &local_sym));
    2069:   PetscCallMPI(MPIU_Allreduce(&local_sym, flg, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
    2070:   PetscFunctionReturn(PETSC_SUCCESS);
    2071: }

    2073: static PetscErrorCode MatIsSymmetric_IS(Mat A, PetscReal tol, PetscBool *flg)
    2074: {
    2075:   Mat_IS   *matis = (Mat_IS *)A->data;
    2076:   PetscBool local_sym;

    2078:   PetscFunctionBegin;
    2079:   if (matis->rmapping != matis->cmapping) {
    2080:     *flg = PETSC_FALSE;
    2081:     PetscFunctionReturn(PETSC_SUCCESS);
    2082:   }
    2083:   PetscCall(MatIsSymmetric(matis->A, tol, &local_sym));
    2084:   PetscCallMPI(MPIU_Allreduce(&local_sym, flg, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
    2085:   PetscFunctionReturn(PETSC_SUCCESS);
    2086: }

    2088: static PetscErrorCode MatIsStructurallySymmetric_IS(Mat A, PetscBool *flg)
    2089: {
    2090:   Mat_IS   *matis = (Mat_IS *)A->data;
    2091:   PetscBool local_sym;

    2093:   PetscFunctionBegin;
    2094:   if (matis->rmapping != matis->cmapping) {
    2095:     *flg = PETSC_FALSE;
    2096:     PetscFunctionReturn(PETSC_SUCCESS);
    2097:   }
    2098:   PetscCall(MatIsStructurallySymmetric(matis->A, &local_sym));
    2099:   PetscCallMPI(MPIU_Allreduce(&local_sym, flg, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
    2100:   PetscFunctionReturn(PETSC_SUCCESS);
    2101: }

    2103: static PetscErrorCode MatDestroy_IS(Mat A)
    2104: {
    2105:   Mat_IS *b = (Mat_IS *)A->data;

    2107:   PetscFunctionBegin;
    2108:   PetscCall(PetscFree(b->bdiag));
    2109:   PetscCall(PetscFree(b->lmattype));
    2110:   PetscCall(MatDestroy(&b->A));
    2111:   PetscCall(VecScatterDestroy(&b->cctx));
    2112:   PetscCall(VecScatterDestroy(&b->rctx));
    2113:   PetscCall(VecDestroy(&b->x));
    2114:   PetscCall(VecDestroy(&b->y));
    2115:   PetscCall(VecDestroy(&b->counter));
    2116:   PetscCall(ISDestroy(&b->getsub_ris));
    2117:   PetscCall(ISDestroy(&b->getsub_cis));
    2118:   if (b->sf != b->csf) {
    2119:     PetscCall(PetscSFDestroy(&b->csf));
    2120:     PetscCall(PetscFree2(b->csf_rootdata, b->csf_leafdata));
    2121:   } else b->csf = NULL;
    2122:   PetscCall(PetscSFDestroy(&b->sf));
    2123:   PetscCall(PetscFree2(b->sf_rootdata, b->sf_leafdata));
    2124:   PetscCall(ISLocalToGlobalMappingDestroy(&b->rmapping));
    2125:   PetscCall(ISLocalToGlobalMappingDestroy(&b->cmapping));
    2126:   PetscCall(MatDestroy(&b->dA));
    2127:   PetscCall(MatDestroy(&b->assembledA));
    2128:   PetscCall(PetscFree(A->data));
    2129:   PetscCall(PetscObjectChangeTypeName((PetscObject)A, NULL));
    2130:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMatType_C", NULL));
    2131:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalMat_C", NULL));
    2132:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMat_C", NULL));
    2133:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISRestoreLocalMat_C", NULL));
    2134:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetPreallocation_C", NULL));
    2135:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISStoreL2L_C", NULL));
    2136:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISFixLocalEmpty_C", NULL));
    2137:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalToGlobalMapping_C", NULL));
    2138:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpiaij_C", NULL));
    2139:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpibaij_C", NULL));
    2140:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpisbaij_C", NULL));
    2141:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqaij_C", NULL));
    2142:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqbaij_C", NULL));
    2143:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqsbaij_C", NULL));
    2144:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_aij_C", NULL));
    2145:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOOLocal_C", NULL));
    2146:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOO_C", NULL));
    2147:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", NULL));
    2148:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetAllowRepeated_C", NULL));
    2149:   PetscFunctionReturn(PETSC_SUCCESS);
    2150: }

    2152: static PetscErrorCode MatMult_IS(Mat A, Vec x, Vec y)
    2153: {
    2154:   Mat_IS     *is   = (Mat_IS *)A->data;
    2155:   PetscScalar zero = 0.0;

    2157:   PetscFunctionBegin;
    2158:   /*  scatter the global vector x into the local work vector */
    2159:   PetscCall(VecScatterBegin(is->cctx, x, is->x, INSERT_VALUES, SCATTER_FORWARD));
    2160:   PetscCall(VecScatterEnd(is->cctx, x, is->x, INSERT_VALUES, SCATTER_FORWARD));

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

    2165:   /* scatter product back into global memory */
    2166:   PetscCall(VecSet(y, zero));
    2167:   PetscCall(VecScatterBegin(is->rctx, is->y, y, ADD_VALUES, SCATTER_REVERSE));
    2168:   PetscCall(VecScatterEnd(is->rctx, is->y, y, ADD_VALUES, SCATTER_REVERSE));
    2169:   PetscFunctionReturn(PETSC_SUCCESS);
    2170: }

    2172: static PetscErrorCode MatMultAdd_IS(Mat A, Vec v1, Vec v2, Vec v3)
    2173: {
    2174:   Vec temp_vec;

    2176:   PetscFunctionBegin; /*  v3 = v2 + A * v1.*/
    2177:   if (v3 != v2) {
    2178:     PetscCall(MatMult(A, v1, v3));
    2179:     PetscCall(VecAXPY(v3, 1.0, v2));
    2180:   } else {
    2181:     PetscCall(VecDuplicate(v2, &temp_vec));
    2182:     PetscCall(MatMult(A, v1, temp_vec));
    2183:     PetscCall(VecAXPY(temp_vec, 1.0, v2));
    2184:     PetscCall(VecCopy(temp_vec, v3));
    2185:     PetscCall(VecDestroy(&temp_vec));
    2186:   }
    2187:   PetscFunctionReturn(PETSC_SUCCESS);
    2188: }

    2190: static PetscErrorCode MatMultTranspose_IS(Mat A, Vec y, Vec x)
    2191: {
    2192:   Mat_IS *is = (Mat_IS *)A->data;

    2194:   PetscFunctionBegin;
    2195:   /*  scatter the global vector x into the local work vector */
    2196:   PetscCall(VecScatterBegin(is->rctx, y, is->y, INSERT_VALUES, SCATTER_FORWARD));
    2197:   PetscCall(VecScatterEnd(is->rctx, y, is->y, INSERT_VALUES, SCATTER_FORWARD));

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

    2202:   /* scatter product back into global vector */
    2203:   PetscCall(VecSet(x, 0));
    2204:   PetscCall(VecScatterBegin(is->cctx, is->x, x, ADD_VALUES, SCATTER_REVERSE));
    2205:   PetscCall(VecScatterEnd(is->cctx, is->x, x, ADD_VALUES, SCATTER_REVERSE));
    2206:   PetscFunctionReturn(PETSC_SUCCESS);
    2207: }

    2209: static PetscErrorCode MatMultTransposeAdd_IS(Mat A, Vec v1, Vec v2, Vec v3)
    2210: {
    2211:   Vec temp_vec;

    2213:   PetscFunctionBegin; /*  v3 = v2 + A' * v1.*/
    2214:   if (v3 != v2) {
    2215:     PetscCall(MatMultTranspose(A, v1, v3));
    2216:     PetscCall(VecAXPY(v3, 1.0, v2));
    2217:   } else {
    2218:     PetscCall(VecDuplicate(v2, &temp_vec));
    2219:     PetscCall(MatMultTranspose(A, v1, temp_vec));
    2220:     PetscCall(VecAXPY(temp_vec, 1.0, v2));
    2221:     PetscCall(VecCopy(temp_vec, v3));
    2222:     PetscCall(VecDestroy(&temp_vec));
    2223:   }
    2224:   PetscFunctionReturn(PETSC_SUCCESS);
    2225: }

    2227: static PetscErrorCode ISLocalToGlobalMappingView_Multi(ISLocalToGlobalMapping mapping, PetscInt lsize, PetscInt gsize, const PetscInt vblocks[], PetscViewer viewer)
    2228: {
    2229:   PetscInt        tr[3], n;
    2230:   const PetscInt *indices;

    2232:   PetscFunctionBegin;
    2233:   tr[0] = IS_LTOGM_FILE_CLASSID;
    2234:   tr[1] = 1;
    2235:   tr[2] = gsize;
    2236:   PetscCall(PetscViewerBinaryWrite(viewer, tr, 3, PETSC_INT));
    2237:   PetscCall(PetscViewerBinaryWriteAll(viewer, vblocks, lsize, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
    2238:   PetscCall(ISLocalToGlobalMappingGetSize(mapping, &n));
    2239:   PetscCall(ISLocalToGlobalMappingGetIndices(mapping, &indices));
    2240:   PetscCall(PetscViewerBinaryWriteAll(viewer, indices, n, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
    2241:   PetscCall(ISLocalToGlobalMappingRestoreIndices(mapping, &indices));
    2242:   PetscFunctionReturn(PETSC_SUCCESS);
    2243: }

    2245: static PetscErrorCode MatView_IS(Mat A, PetscViewer viewer)
    2246: {
    2247:   Mat_IS                *a = (Mat_IS *)A->data;
    2248:   PetscViewer            sviewer;
    2249:   PetscBool              isascii, isbinary, viewl2g = PETSC_FALSE, native;
    2250:   PetscViewerFormat      format;
    2251:   ISLocalToGlobalMapping rmap, cmap;

    2253:   PetscFunctionBegin;
    2254:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    2255:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
    2256:   PetscCall(PetscViewerGetFormat(viewer, &format));
    2257:   native = (PetscBool)(format == PETSC_VIEWER_NATIVE);
    2258:   if (native) {
    2259:     rmap = A->rmap->mapping;
    2260:     cmap = A->cmap->mapping;
    2261:   } else {
    2262:     rmap = a->rmapping;
    2263:     cmap = a->cmapping;
    2264:   }
    2265:   if (isascii) {
    2266:     if (format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(PETSC_SUCCESS);
    2267:     if (format == PETSC_VIEWER_ASCII_INFO_DETAIL || format == PETSC_VIEWER_ASCII_MATLAB) viewl2g = PETSC_TRUE;
    2268:   } else if (isbinary) {
    2269:     PetscInt        tr[6], nr, nc, lsize = 0;
    2270:     char            lmattype[64] = {'\0'};
    2271:     PetscMPIInt     size;
    2272:     PetscBool       skipHeader, vbs = PETSC_FALSE;
    2273:     IS              is;
    2274:     const PetscInt *vblocks = NULL;

    2276:     PetscCall(PetscViewerSetUp(viewer));
    2277:     PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)A)->prefix, "-mat_is_view_variableblocksizes", &vbs, NULL));
    2278:     if (vbs) {
    2279:       PetscCall(MatGetVariableBlockSizes(a->A, &lsize, &vblocks));
    2280:       PetscCall(PetscMPIIntCast(lsize, &size));
    2281:       PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &size, 1, MPI_INT, MPI_SUM, PetscObjectComm((PetscObject)viewer)));
    2282:     } else {
    2283:       PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)viewer), &size));
    2284:     }
    2285:     tr[0] = MAT_FILE_CLASSID;
    2286:     tr[1] = A->rmap->N;
    2287:     tr[2] = A->cmap->N;
    2288:     tr[3] = -size; /* AIJ stores nnz here */
    2289:     tr[4] = (PetscInt)(rmap == cmap);
    2290:     tr[5] = a->allow_repeated;
    2291:     PetscCall(PetscSNPrintf(lmattype, sizeof(lmattype), "%s", a->lmattype));

    2293:     PetscCall(PetscViewerBinaryWrite(viewer, tr, PETSC_STATIC_ARRAY_LENGTH(tr), PETSC_INT));
    2294:     PetscCall(PetscViewerBinaryWrite(viewer, lmattype, sizeof(lmattype), PETSC_CHAR));

    2296:     /* first dump l2g info (we need the header for proper loading on different number of processes) */
    2297:     PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));
    2298:     PetscCall(PetscViewerBinarySetSkipHeader(viewer, PETSC_FALSE));
    2299:     if (vbs) {
    2300:       PetscCall(ISLocalToGlobalMappingView_Multi(rmap, lsize, size, vblocks, viewer));
    2301:       if (cmap != rmap) PetscCall(ISLocalToGlobalMappingView_Multi(cmap, lsize, size, vblocks, viewer));
    2302:       PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)viewer), lsize, vblocks, PETSC_USE_POINTER, &is));
    2303:       PetscCall(ISView(is, viewer));
    2304:       PetscCall(ISView(is, viewer));
    2305:       PetscCall(ISDestroy(&is));
    2306:     } else {
    2307:       PetscCall(ISLocalToGlobalMappingView(rmap, viewer));
    2308:       if (cmap != rmap) PetscCall(ISLocalToGlobalMappingView(cmap, viewer));

    2310:       /* then the sizes of the local matrices */
    2311:       PetscCall(MatGetSize(a->A, &nr, &nc));
    2312:       PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)viewer), 1, &nr, PETSC_USE_POINTER, &is));
    2313:       PetscCall(ISView(is, viewer));
    2314:       PetscCall(ISDestroy(&is));
    2315:       PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)viewer), 1, &nc, PETSC_USE_POINTER, &is));
    2316:       PetscCall(ISView(is, viewer));
    2317:       PetscCall(ISDestroy(&is));
    2318:     }
    2319:     PetscCall(PetscViewerBinarySetSkipHeader(viewer, skipHeader));
    2320:   }
    2321:   if (format == PETSC_VIEWER_ASCII_MATLAB) {
    2322:     char        name[64];
    2323:     PetscMPIInt size, rank;

    2325:     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)viewer), &size));
    2326:     PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
    2327:     if (size > 1) PetscCall(PetscSNPrintf(name, sizeof(name), "lmat_%d", rank));
    2328:     else PetscCall(PetscSNPrintf(name, sizeof(name), "lmat"));
    2329:     PetscCall(PetscObjectSetName((PetscObject)a->A, name));
    2330:   }

    2332:   /* Dump the local matrices */
    2333:   if (isbinary) { /* ViewerGetSubViewer does not work in parallel */
    2334:     PetscBool   isaij;
    2335:     PetscInt    nr, nc;
    2336:     Mat         lA, B;
    2337:     Mat_MPIAIJ *b;

    2339:     /* We create a temporary MPIAIJ matrix that stores the unassembled operator */
    2340:     PetscCall(PetscObjectBaseTypeCompare((PetscObject)a->A, MATAIJ, &isaij));
    2341:     if (!isaij) PetscCall(MatConvert(a->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &lA));
    2342:     else {
    2343:       PetscCall(PetscObjectReference((PetscObject)a->A));
    2344:       lA = a->A;
    2345:     }
    2346:     PetscCall(MatCreate(PetscObjectComm((PetscObject)viewer), &B));
    2347:     PetscCall(MatSetType(B, MATMPIAIJ));
    2348:     PetscCall(MatGetSize(lA, &nr, &nc));
    2349:     PetscCall(MatSetSizes(B, nr, nc, PETSC_DECIDE, PETSC_DECIDE));
    2350:     PetscCall(MatMPIAIJSetPreallocation(B, 0, NULL, 0, NULL));

    2352:     b = (Mat_MPIAIJ *)B->data;
    2353:     PetscCall(MatDestroy(&b->A));
    2354:     b->A = lA;

    2356:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
    2357:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
    2358:     PetscCall(MatView(B, viewer));
    2359:     PetscCall(MatDestroy(&B));
    2360:   } else {
    2361:     PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
    2362:     PetscCall(MatView(a->A, sviewer));
    2363:     PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
    2364:   }

    2366:   /* with ASCII, we dump the l2gmaps at the end */
    2367:   if (viewl2g) {
    2368:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
    2369:       PetscCall(PetscObjectSetName((PetscObject)rmap, "row"));
    2370:       PetscCall(ISLocalToGlobalMappingView(rmap, viewer));
    2371:       PetscCall(PetscObjectSetName((PetscObject)cmap, "col"));
    2372:       PetscCall(ISLocalToGlobalMappingView(cmap, viewer));
    2373:     } else {
    2374:       PetscCall(ISLocalToGlobalMappingView(rmap, viewer));
    2375:       PetscCall(ISLocalToGlobalMappingView(cmap, viewer));
    2376:     }
    2377:   }
    2378:   PetscFunctionReturn(PETSC_SUCCESS);
    2379: }

    2381: static PetscErrorCode ISLocalToGlobalMappingHasRepeatedLocal_Private(ISLocalToGlobalMapping map, PetscBool *has)
    2382: {
    2383:   const PetscInt *idxs;
    2384:   PetscHSetI      ht;
    2385:   PetscInt        n, bs;

    2387:   PetscFunctionBegin;
    2388:   PetscCall(ISLocalToGlobalMappingGetSize(map, &n));
    2389:   PetscCall(ISLocalToGlobalMappingGetBlockSize(map, &bs));
    2390:   PetscCall(ISLocalToGlobalMappingGetBlockIndices(map, &idxs));
    2391:   PetscCall(PetscHSetICreate(&ht));
    2392:   *has = PETSC_FALSE;
    2393:   for (PetscInt i = 0; i < n / bs; i++) {
    2394:     PetscBool missing = PETSC_TRUE;
    2395:     if (idxs[i] < 0) continue;
    2396:     PetscCall(PetscHSetIQueryAdd(ht, idxs[i], &missing));
    2397:     if (!missing) {
    2398:       *has = PETSC_TRUE;
    2399:       break;
    2400:     }
    2401:   }
    2402:   PetscCall(PetscHSetIDestroy(&ht));
    2403:   PetscFunctionReturn(PETSC_SUCCESS);
    2404: }

    2406: static PetscErrorCode MatLoad_IS(Mat A, PetscViewer viewer)
    2407: {
    2408:   ISLocalToGlobalMapping rmap, cmap;
    2409:   MPI_Comm               comm = PetscObjectComm((PetscObject)A);
    2410:   PetscBool              isbinary, samel, allow, isbaij;
    2411:   PetscInt               tr[6], M, N, nr, nc, Asize, isn;
    2412:   const PetscInt        *idx;
    2413:   PetscMPIInt            size;
    2414:   char                   lmattype[64];
    2415:   Mat                    dA, lA;
    2416:   IS                     is;

    2418:   PetscFunctionBegin;
    2419:   PetscCheckSameComm(A, 1, viewer, 2);
    2420:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
    2421:   PetscCheck(isbinary, PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Invalid viewer of type %s", ((PetscObject)viewer)->type_name);

    2423:   PetscCall(PetscViewerBinaryRead(viewer, tr, PETSC_STATIC_ARRAY_LENGTH(tr), NULL, PETSC_INT));
    2424:   PetscCheck(tr[0] == MAT_FILE_CLASSID, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a matrix next in file");
    2425:   PetscCheck(tr[1] >= 0, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a IS matrix next in file");
    2426:   PetscCheck(tr[2] >= 0, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a IS matrix next in file");
    2427:   PetscCheck(tr[3] < 0, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a IS matrix next in file");
    2428:   PetscCheck(tr[4] == 0 || tr[4] == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a IS matrix next in file");
    2429:   PetscCheck(tr[5] == 0 || tr[5] == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a IS matrix next in file");
    2430:   M     = tr[1];
    2431:   N     = tr[2];
    2432:   Asize = -tr[3];
    2433:   samel = (PetscBool)tr[4];
    2434:   allow = (PetscBool)tr[5];
    2435:   PetscCall(PetscViewerBinaryRead(viewer, lmattype, sizeof(lmattype), NULL, PETSC_CHAR));

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

    2441:   /* set global sizes if not set already */
    2442:   if (A->rmap->N < 0) A->rmap->N = M;
    2443:   if (A->cmap->N < 0) A->cmap->N = N;
    2444:   PetscCall(PetscLayoutSetUp(A->rmap));
    2445:   PetscCall(PetscLayoutSetUp(A->cmap));
    2446:   PetscCheck(M == A->rmap->N, comm, PETSC_ERR_ARG_SIZ, "Matrix rows should be %" PetscInt_FMT ", found %" PetscInt_FMT, M, A->rmap->N);
    2447:   PetscCheck(N == A->cmap->N, comm, PETSC_ERR_ARG_SIZ, "Matrix columns should be %" PetscInt_FMT ", found %" PetscInt_FMT, N, A->cmap->N);

    2449:   /* load l2g maps */
    2450:   PetscCall(ISLocalToGlobalMappingCreate(comm, 0, 0, NULL, PETSC_USE_POINTER, &rmap));
    2451:   PetscCall(ISLocalToGlobalMappingLoad(rmap, viewer));
    2452:   if (!samel) {
    2453:     PetscCall(ISLocalToGlobalMappingCreate(comm, 0, 0, NULL, PETSC_USE_POINTER, &cmap));
    2454:     PetscCall(ISLocalToGlobalMappingLoad(cmap, viewer));
    2455:   } else {
    2456:     PetscCall(PetscObjectReference((PetscObject)rmap));
    2457:     cmap = rmap;
    2458:   }

    2460:   /* load sizes of local matrices */
    2461:   PetscCall(ISCreate(comm, &is));
    2462:   PetscCall(ISSetType(is, ISGENERAL));
    2463:   PetscCall(ISLoad(is, viewer));
    2464:   PetscCall(ISGetLocalSize(is, &isn));
    2465:   PetscCall(ISGetIndices(is, &idx));
    2466:   nr = 0;
    2467:   for (PetscInt i = 0; i < isn; i++) nr += idx[i];
    2468:   PetscCall(ISRestoreIndices(is, &idx));
    2469:   PetscCall(ISDestroy(&is));
    2470:   PetscCall(ISCreate(comm, &is));
    2471:   PetscCall(ISSetType(is, ISGENERAL));
    2472:   PetscCall(ISLoad(is, viewer));
    2473:   PetscCall(ISGetLocalSize(is, &isn));
    2474:   PetscCall(ISGetIndices(is, &idx));
    2475:   nc = 0;
    2476:   for (PetscInt i = 0; i < isn; i++) nc += idx[i];
    2477:   PetscCall(ISRestoreIndices(is, &idx));
    2478:   PetscCall(ISDestroy(&is));

    2480:   /* now load the unassembled operator */
    2481:   PetscCall(MatCreate(comm, &dA));
    2482:   PetscCall(MatSetType(dA, MATMPIAIJ));
    2483:   PetscCall(MatSetSizes(dA, nr, nc, PETSC_DECIDE, PETSC_DECIDE));
    2484:   PetscCall(MatLoad(dA, viewer));
    2485:   PetscCall(MatMPIAIJGetSeqAIJ(dA, &lA, NULL, NULL));
    2486:   PetscCall(PetscObjectReference((PetscObject)lA));
    2487:   PetscCall(MatDestroy(&dA));

    2489:   /* and convert to the desired format */
    2490:   PetscCall(PetscStrcmpAny(lmattype, &isbaij, MATSBAIJ, MATSEQSBAIJ, ""));
    2491:   if (isbaij) PetscCall(MatSetOption(lA, MAT_SYMMETRIC, PETSC_TRUE));
    2492:   PetscCall(MatConvert(lA, lmattype, MAT_INPLACE_MATRIX, &lA));

    2494:   /* check if we actually have repeated entries */
    2495:   if (allow) {
    2496:     PetscBool rhas, chas, hasrepeated;

    2498:     PetscCall(ISLocalToGlobalMappingHasRepeatedLocal_Private(rmap, &rhas));
    2499:     if (rmap != cmap) PetscCall(ISLocalToGlobalMappingHasRepeatedLocal_Private(cmap, &chas));
    2500:     else chas = rhas;
    2501:     hasrepeated = (PetscBool)(rhas || chas);
    2502:     PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &hasrepeated, 1, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A)));
    2503:     if (!hasrepeated) allow = PETSC_FALSE;
    2504:   }

    2506:   /* assemble the MATIS object */
    2507:   PetscCall(MatISSetAllowRepeated(A, allow));
    2508:   PetscCall(MatSetLocalToGlobalMapping(A, rmap, cmap));
    2509:   PetscCall(MatISSetLocalMat(A, lA));
    2510:   PetscCall(MatDestroy(&lA));
    2511:   PetscCall(ISLocalToGlobalMappingDestroy(&rmap));
    2512:   PetscCall(ISLocalToGlobalMappingDestroy(&cmap));
    2513:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    2514:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    2515:   PetscFunctionReturn(PETSC_SUCCESS);
    2516: }

    2518: static PetscErrorCode MatInvertBlockDiagonal_IS(Mat mat, const PetscScalar **values)
    2519: {
    2520:   Mat_IS            *is = (Mat_IS *)mat->data;
    2521:   MPI_Datatype       nodeType;
    2522:   const PetscScalar *lv;
    2523:   PetscInt           bs;
    2524:   PetscMPIInt        mbs;

    2526:   PetscFunctionBegin;
    2527:   PetscCall(MatGetBlockSize(mat, &bs));
    2528:   PetscCall(MatSetBlockSize(is->A, bs));
    2529:   PetscCall(MatInvertBlockDiagonal(is->A, &lv));
    2530:   if (!is->bdiag) PetscCall(PetscMalloc1(bs * mat->rmap->n, &is->bdiag));
    2531:   PetscCall(PetscMPIIntCast(bs, &mbs));
    2532:   PetscCallMPI(MPI_Type_contiguous(mbs, MPIU_SCALAR, &nodeType));
    2533:   PetscCallMPI(MPI_Type_commit(&nodeType));
    2534:   PetscCall(PetscSFReduceBegin(is->sf, nodeType, lv, is->bdiag, MPI_REPLACE));
    2535:   PetscCall(PetscSFReduceEnd(is->sf, nodeType, lv, is->bdiag, MPI_REPLACE));
    2536:   PetscCallMPI(MPI_Type_free(&nodeType));
    2537:   if (values) *values = is->bdiag;
    2538:   PetscFunctionReturn(PETSC_SUCCESS);
    2539: }

    2541: static PetscErrorCode MatISSetUpScatters_Private(Mat A)
    2542: {
    2543:   Vec             cglobal, rglobal;
    2544:   IS              from;
    2545:   Mat_IS         *is = (Mat_IS *)A->data;
    2546:   PetscScalar     sum;
    2547:   const PetscInt *garray;
    2548:   PetscInt        nr, rbs, nc, cbs;
    2549:   VecType         rtype;

    2551:   PetscFunctionBegin;
    2552:   PetscCall(ISLocalToGlobalMappingGetSize(is->rmapping, &nr));
    2553:   PetscCall(ISLocalToGlobalMappingGetBlockSize(is->rmapping, &rbs));
    2554:   PetscCall(ISLocalToGlobalMappingGetSize(is->cmapping, &nc));
    2555:   PetscCall(ISLocalToGlobalMappingGetBlockSize(is->cmapping, &cbs));
    2556:   PetscCall(VecDestroy(&is->x));
    2557:   PetscCall(VecDestroy(&is->y));
    2558:   PetscCall(VecDestroy(&is->counter));
    2559:   PetscCall(VecScatterDestroy(&is->rctx));
    2560:   PetscCall(VecScatterDestroy(&is->cctx));
    2561:   PetscCall(MatCreateVecs(is->A, &is->x, &is->y));
    2562:   PetscCall(VecBindToCPU(is->y, PETSC_TRUE));
    2563:   PetscCall(VecGetRootType_Private(is->y, &rtype));
    2564:   PetscCall(PetscFree(A->defaultvectype));
    2565:   PetscCall(PetscStrallocpy(rtype, &A->defaultvectype));
    2566:   PetscCall(MatCreateVecs(A, &cglobal, &rglobal));
    2567:   PetscCall(VecBindToCPU(rglobal, PETSC_TRUE));
    2568:   PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->rmapping, &garray));
    2569:   PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), rbs, nr / rbs, garray, PETSC_USE_POINTER, &from));
    2570:   PetscCall(VecScatterCreate(rglobal, from, is->y, NULL, &is->rctx));
    2571:   PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->rmapping, &garray));
    2572:   PetscCall(ISDestroy(&from));
    2573:   if (is->rmapping != is->cmapping) {
    2574:     PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->cmapping, &garray));
    2575:     PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), cbs, nc / cbs, garray, PETSC_USE_POINTER, &from));
    2576:     PetscCall(VecScatterCreate(cglobal, from, is->x, NULL, &is->cctx));
    2577:     PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->cmapping, &garray));
    2578:     PetscCall(ISDestroy(&from));
    2579:   } else {
    2580:     PetscCall(PetscObjectReference((PetscObject)is->rctx));
    2581:     is->cctx = is->rctx;
    2582:   }
    2583:   PetscCall(VecDestroy(&cglobal));

    2585:   /* interface counter vector (local) */
    2586:   PetscCall(VecDuplicate(is->y, &is->counter));
    2587:   PetscCall(VecBindToCPU(is->counter, PETSC_TRUE));
    2588:   PetscCall(VecSet(is->y, 1.));
    2589:   PetscCall(VecScatterBegin(is->rctx, is->y, rglobal, ADD_VALUES, SCATTER_REVERSE));
    2590:   PetscCall(VecScatterEnd(is->rctx, is->y, rglobal, ADD_VALUES, SCATTER_REVERSE));
    2591:   PetscCall(VecScatterBegin(is->rctx, rglobal, is->counter, INSERT_VALUES, SCATTER_FORWARD));
    2592:   PetscCall(VecScatterEnd(is->rctx, rglobal, is->counter, INSERT_VALUES, SCATTER_FORWARD));
    2593:   PetscCall(VecBindToCPU(is->y, PETSC_FALSE));
    2594:   PetscCall(VecBindToCPU(is->counter, PETSC_FALSE));

    2596:   /* special functions for block-diagonal matrices */
    2597:   PetscCall(VecSum(rglobal, &sum));
    2598:   A->ops->invertblockdiagonal = NULL;
    2599:   if ((PetscInt)(PetscRealPart(sum)) == A->rmap->N && A->rmap->N == A->cmap->N && is->rmapping == is->cmapping) A->ops->invertblockdiagonal = MatInvertBlockDiagonal_IS;
    2600:   PetscCall(VecDestroy(&rglobal));

    2602:   /* setup SF for general purpose shared indices based communications */
    2603:   PetscCall(MatISSetUpSF_IS(A));
    2604:   PetscFunctionReturn(PETSC_SUCCESS);
    2605: }

    2607: static PetscErrorCode MatISFilterL2GMap(Mat A, ISLocalToGlobalMapping map, ISLocalToGlobalMapping *nmap, ISLocalToGlobalMapping *lmap)
    2608: {
    2609:   Mat_IS                    *matis = (Mat_IS *)A->data;
    2610:   IS                         is;
    2611:   ISLocalToGlobalMappingType l2gtype;
    2612:   const PetscInt            *idxs;
    2613:   PetscHSetI                 ht;
    2614:   PetscInt                  *nidxs;
    2615:   PetscInt                   i, n, bs, c;
    2616:   PetscBool                  flg[] = {PETSC_FALSE, PETSC_FALSE};

    2618:   PetscFunctionBegin;
    2619:   PetscCall(ISLocalToGlobalMappingGetSize(map, &n));
    2620:   PetscCall(ISLocalToGlobalMappingGetBlockSize(map, &bs));
    2621:   PetscCall(ISLocalToGlobalMappingGetBlockIndices(map, &idxs));
    2622:   PetscCall(PetscHSetICreate(&ht));
    2623:   PetscCall(PetscMalloc1(n / bs, &nidxs));
    2624:   for (i = 0, c = 0; i < n / bs; i++) {
    2625:     PetscBool missing = PETSC_TRUE;
    2626:     if (idxs[i] < 0) {
    2627:       flg[0] = PETSC_TRUE;
    2628:       continue;
    2629:     }
    2630:     if (!matis->allow_repeated) PetscCall(PetscHSetIQueryAdd(ht, idxs[i], &missing));
    2631:     if (!missing) flg[1] = PETSC_TRUE;
    2632:     else nidxs[c++] = idxs[i];
    2633:   }
    2634:   PetscCall(PetscHSetIDestroy(&ht));
    2635:   PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, flg, 2, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A)));
    2636:   if (!flg[0] && !flg[1]) { /* Entries are all non negative and unique */
    2637:     *nmap = NULL;
    2638:     *lmap = NULL;
    2639:     PetscCall(PetscFree(nidxs));
    2640:     PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(map, &idxs));
    2641:     PetscFunctionReturn(PETSC_SUCCESS);
    2642:   }

    2644:   /* New l2g map without negative indices (and repeated indices if not allowed) */
    2645:   PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), bs, c, nidxs, PETSC_USE_POINTER, &is));
    2646:   PetscCall(ISLocalToGlobalMappingCreateIS(is, nmap));
    2647:   PetscCall(ISDestroy(&is));
    2648:   PetscCall(ISLocalToGlobalMappingGetType(map, &l2gtype));
    2649:   PetscCall(ISLocalToGlobalMappingSetType(*nmap, l2gtype));

    2651:   /* New local l2g map for repeated indices if not allowed */
    2652:   PetscCall(ISGlobalToLocalMappingApplyBlock(*nmap, IS_GTOLM_MASK, n / bs, idxs, NULL, nidxs));
    2653:   PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, n / bs, nidxs, PETSC_USE_POINTER, &is));
    2654:   PetscCall(ISLocalToGlobalMappingCreateIS(is, lmap));
    2655:   PetscCall(ISDestroy(&is));
    2656:   PetscCall(PetscFree(nidxs));
    2657:   PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(map, &idxs));
    2658:   PetscFunctionReturn(PETSC_SUCCESS);
    2659: }

    2661: static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A, ISLocalToGlobalMapping rmapping, ISLocalToGlobalMapping cmapping)
    2662: {
    2663:   Mat_IS                *is            = (Mat_IS *)A->data;
    2664:   ISLocalToGlobalMapping localrmapping = NULL, localcmapping = NULL;
    2665:   PetscInt               nr, rbs, nc, cbs;
    2666:   PetscBool              cong, freem[] = {PETSC_FALSE, PETSC_FALSE};

    2668:   PetscFunctionBegin;
    2669:   if (rmapping) PetscCheckSameComm(A, 1, rmapping, 2);
    2670:   if (cmapping) PetscCheckSameComm(A, 1, cmapping, 3);

    2672:   PetscCall(ISLocalToGlobalMappingDestroy(&is->rmapping));
    2673:   PetscCall(ISLocalToGlobalMappingDestroy(&is->cmapping));
    2674:   PetscCall(PetscLayoutSetUp(A->rmap));
    2675:   PetscCall(PetscLayoutSetUp(A->cmap));
    2676:   PetscCall(MatHasCongruentLayouts(A, &cong));

    2678:   /* If NULL, local space matches global space */
    2679:   if (!rmapping) {
    2680:     IS is;

    2682:     PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), A->rmap->N, 0, 1, &is));
    2683:     PetscCall(ISLocalToGlobalMappingCreateIS(is, &rmapping));
    2684:     PetscCall(ISLocalToGlobalMappingSetBlockSize(rmapping, A->rmap->bs));
    2685:     PetscCall(ISDestroy(&is));
    2686:     freem[0] = PETSC_TRUE;
    2687:     if (!cmapping && cong && A->rmap->bs == A->cmap->bs) cmapping = rmapping;
    2688:   } else if (!is->islocalref) { /* check if the l2g map has negative or repeated entries */
    2689:     PetscCall(MatISFilterL2GMap(A, rmapping, &is->rmapping, &localrmapping));
    2690:     if (rmapping == cmapping) {
    2691:       PetscCall(PetscObjectReference((PetscObject)is->rmapping));
    2692:       is->cmapping = is->rmapping;
    2693:       PetscCall(PetscObjectReference((PetscObject)localrmapping));
    2694:       localcmapping = localrmapping;
    2695:     }
    2696:   }
    2697:   if (!cmapping) {
    2698:     IS is;

    2700:     PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), A->cmap->N, 0, 1, &is));
    2701:     PetscCall(ISLocalToGlobalMappingCreateIS(is, &cmapping));
    2702:     PetscCall(ISLocalToGlobalMappingSetBlockSize(cmapping, A->cmap->bs));
    2703:     PetscCall(ISDestroy(&is));
    2704:     freem[1] = PETSC_TRUE;
    2705:   } else if (cmapping != rmapping && !is->islocalref) { /* check if the l2g map has negative or repeated entries */
    2706:     PetscCall(MatISFilterL2GMap(A, cmapping, &is->cmapping, &localcmapping));
    2707:   }
    2708:   if (!is->rmapping) {
    2709:     PetscCall(PetscObjectReference((PetscObject)rmapping));
    2710:     is->rmapping = rmapping;
    2711:   }
    2712:   if (!is->cmapping) {
    2713:     PetscCall(PetscObjectReference((PetscObject)cmapping));
    2714:     is->cmapping = cmapping;
    2715:   }

    2717:   /* Clean up */
    2718:   is->lnnzstate = 0;
    2719:   PetscCall(MatDestroy(&is->dA));
    2720:   PetscCall(MatDestroy(&is->assembledA));
    2721:   PetscCall(MatDestroy(&is->A));
    2722:   if (is->csf != is->sf) {
    2723:     PetscCall(PetscSFDestroy(&is->csf));
    2724:     PetscCall(PetscFree2(is->csf_rootdata, is->csf_leafdata));
    2725:   } else is->csf = NULL;
    2726:   PetscCall(PetscSFDestroy(&is->sf));
    2727:   PetscCall(PetscFree2(is->sf_rootdata, is->sf_leafdata));
    2728:   PetscCall(PetscFree(is->bdiag));

    2730:   /* check if the two mappings are actually the same for square matrices since MATIS has some optimization for this case
    2731:      (DOLFIN passes 2 different objects) */
    2732:   PetscCall(ISLocalToGlobalMappingGetSize(is->rmapping, &nr));
    2733:   PetscCall(ISLocalToGlobalMappingGetBlockSize(is->rmapping, &rbs));
    2734:   PetscCall(ISLocalToGlobalMappingGetSize(is->cmapping, &nc));
    2735:   PetscCall(ISLocalToGlobalMappingGetBlockSize(is->cmapping, &cbs));
    2736:   if (is->rmapping != is->cmapping && cong) {
    2737:     PetscBool same = PETSC_FALSE;
    2738:     if (nr == nc && cbs == rbs) {
    2739:       const PetscInt *idxs1, *idxs2;

    2741:       PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->rmapping, &idxs1));
    2742:       PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->cmapping, &idxs2));
    2743:       PetscCall(PetscArraycmp(idxs1, idxs2, nr / rbs, &same));
    2744:       PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->rmapping, &idxs1));
    2745:       PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->cmapping, &idxs2));
    2746:     }
    2747:     PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &same, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
    2748:     if (same) {
    2749:       PetscCall(ISLocalToGlobalMappingDestroy(&is->cmapping));
    2750:       PetscCall(PetscObjectReference((PetscObject)is->rmapping));
    2751:       is->cmapping = is->rmapping;
    2752:     }
    2753:   }
    2754:   PetscCall(PetscLayoutSetBlockSize(A->rmap, rbs));
    2755:   PetscCall(PetscLayoutSetBlockSize(A->cmap, cbs));
    2756:   /* Pass the user defined maps to the layout */
    2757:   PetscCall(PetscLayoutSetISLocalToGlobalMapping(A->rmap, rmapping));
    2758:   PetscCall(PetscLayoutSetISLocalToGlobalMapping(A->cmap, cmapping));
    2759:   if (freem[0]) PetscCall(ISLocalToGlobalMappingDestroy(&rmapping));
    2760:   if (freem[1]) PetscCall(ISLocalToGlobalMappingDestroy(&cmapping));

    2762:   if (!is->islocalref) {
    2763:     /* Create the local matrix A */
    2764:     PetscCall(MatCreate(PETSC_COMM_SELF, &is->A));
    2765:     PetscCall(MatSetType(is->A, is->lmattype));
    2766:     PetscCall(MatSetSizes(is->A, nr, nc, nr, nc));
    2767:     PetscCall(MatSetBlockSizes(is->A, rbs, cbs));
    2768:     PetscCall(MatSetOptionsPrefix(is->A, "is_"));
    2769:     PetscCall(MatAppendOptionsPrefix(is->A, ((PetscObject)A)->prefix));
    2770:     PetscCall(PetscLayoutSetUp(is->A->rmap));
    2771:     PetscCall(PetscLayoutSetUp(is->A->cmap));
    2772:     PetscCall(MatSetLocalToGlobalMapping(is->A, localrmapping, localcmapping));
    2773:     PetscCall(ISLocalToGlobalMappingDestroy(&localrmapping));
    2774:     PetscCall(ISLocalToGlobalMappingDestroy(&localcmapping));
    2775:     /* setup scatters and local vectors for MatMult */
    2776:     PetscCall(MatISSetUpScatters_Private(A));
    2777:   }
    2778:   A->preallocated = PETSC_TRUE;
    2779:   PetscFunctionReturn(PETSC_SUCCESS);
    2780: }

    2782: static PetscErrorCode MatSetUp_IS(Mat A)
    2783: {
    2784:   Mat_IS                *is = (Mat_IS *)A->data;
    2785:   ISLocalToGlobalMapping rmap, cmap;

    2787:   PetscFunctionBegin;
    2788:   if (!is->sf) {
    2789:     PetscCall(MatGetLocalToGlobalMapping(A, &rmap, &cmap));
    2790:     PetscCall(MatSetLocalToGlobalMapping(A, rmap, cmap));
    2791:   }
    2792:   PetscFunctionReturn(PETSC_SUCCESS);
    2793: }

    2795: static PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
    2796: {
    2797:   Mat_IS  *is = (Mat_IS *)mat->data;
    2798:   PetscInt buf[2 * MATIS_MAX_ENTRIES_INSERTION], *rows_l = NULL, *cols_l = NULL;

    2800:   PetscFunctionBegin;
    2801:   IndexSpaceGet(buf, m, n, rows_l, cols_l);
    2802:   PetscCall(ISGlobalToLocalMappingApply(is->rmapping, IS_GTOLM_MASK, m, rows, &m, rows_l));
    2803:   if (m != n || rows != cols || is->cmapping != is->rmapping) {
    2804:     PetscCall(ISGlobalToLocalMappingApply(is->cmapping, IS_GTOLM_MASK, n, cols, &n, cols_l));
    2805:     PetscCall(MatSetValues(is->A, m, rows_l, n, cols_l, values, addv));
    2806:   } else {
    2807:     PetscCall(MatSetValues(is->A, m, rows_l, m, rows_l, values, addv));
    2808:   }
    2809:   IndexSpaceRestore(buf, m, n, rows_l, cols_l);
    2810:   PetscFunctionReturn(PETSC_SUCCESS);
    2811: }

    2813: static PetscErrorCode MatSetValuesBlocked_IS(Mat mat, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
    2814: {
    2815:   Mat_IS  *is = (Mat_IS *)mat->data;
    2816:   PetscInt buf[2 * MATIS_MAX_ENTRIES_INSERTION], *rows_l = NULL, *cols_l = NULL;

    2818:   PetscFunctionBegin;
    2819:   IndexSpaceGet(buf, m, n, rows_l, cols_l);
    2820:   PetscCall(ISGlobalToLocalMappingApplyBlock(is->rmapping, IS_GTOLM_MASK, m, rows, &m, rows_l));
    2821:   if (m != n || rows != cols || is->cmapping != is->rmapping) {
    2822:     PetscCall(ISGlobalToLocalMappingApplyBlock(is->cmapping, IS_GTOLM_MASK, n, cols, &n, cols_l));
    2823:     PetscCall(MatSetValuesBlocked(is->A, m, rows_l, n, cols_l, values, addv));
    2824:   } else {
    2825:     PetscCall(MatSetValuesBlocked(is->A, m, rows_l, m, rows_l, values, addv));
    2826:   }
    2827:   IndexSpaceRestore(buf, m, n, rows_l, cols_l);
    2828:   PetscFunctionReturn(PETSC_SUCCESS);
    2829: }

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

    2835:   PetscFunctionBegin;
    2836:   if (is->A->rmap->mapping || is->A->cmap->mapping) {
    2837:     PetscCall(MatSetValuesLocal(is->A, m, rows, n, cols, values, addv));
    2838:   } else {
    2839:     PetscCall(MatSetValues(is->A, m, rows, n, cols, values, addv));
    2840:   }
    2841:   PetscFunctionReturn(PETSC_SUCCESS);
    2842: }

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

    2848:   PetscFunctionBegin;
    2849:   if (is->A->rmap->mapping || is->A->cmap->mapping) {
    2850:     PetscCall(MatSetValuesBlockedLocal(is->A, m, rows, n, cols, values, addv));
    2851:   } else {
    2852:     PetscCall(MatSetValuesBlocked(is->A, m, rows, n, cols, values, addv));
    2853:   }
    2854:   PetscFunctionReturn(PETSC_SUCCESS);
    2855: }

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

    2861:   PetscFunctionBegin;
    2862:   if (!n) PetscFunctionReturn(PETSC_SUCCESS);
    2863:   is->pure_neumann = PETSC_FALSE;

    2865:   if (columns) {
    2866:     PetscCall(MatZeroRowsColumns(is->A, n, rows, diag, NULL, NULL));
    2867:   } else {
    2868:     PetscCall(MatZeroRows(is->A, n, rows, diag, NULL, NULL));
    2869:   }
    2870:   if (diag != 0.) {
    2871:     const PetscScalar *array;

    2873:     PetscCall(VecGetArrayRead(is->counter, &array));
    2874:     for (PetscInt i = 0; i < n; i++) PetscCall(MatSetValue(is->A, rows[i], rows[i], diag / (array[rows[i]]), INSERT_VALUES));
    2875:     PetscCall(VecRestoreArrayRead(is->counter, &array));
    2876:     PetscCall(MatAssemblyBegin(is->A, MAT_FINAL_ASSEMBLY));
    2877:     PetscCall(MatAssemblyEnd(is->A, MAT_FINAL_ASSEMBLY));
    2878:   }
    2879:   PetscFunctionReturn(PETSC_SUCCESS);
    2880: }

    2882: static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b, PetscBool columns)
    2883: {
    2884:   Mat_IS   *matis = (Mat_IS *)A->data;
    2885:   PetscInt  nr, nl, len;
    2886:   PetscInt *lrows = NULL;

    2888:   PetscFunctionBegin;
    2889:   if (PetscUnlikelyDebug(columns || diag != 0. || (x && b))) {
    2890:     PetscBool cong;

    2892:     PetscCall(PetscLayoutCompare(A->rmap, A->cmap, &cong));
    2893:     cong = (PetscBool)(cong && matis->sf == matis->csf);
    2894:     PetscCheck(cong || !columns, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Columns can be zeroed if and only if A->rmap and A->cmap are congruent and the l2g maps are the same for MATIS");
    2895:     PetscCheck(cong || diag == 0., PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Nonzero diagonal value supported if and only if A->rmap and A->cmap are congruent and the l2g maps are the same for MATIS");
    2896:     PetscCheck(cong || !x || !b, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "A->rmap and A->cmap need to be congruent, and the l2g maps be the same");
    2897:   }
    2898:   PetscCall(MatGetSize(matis->A, &nl, NULL));
    2899:   /* get locally owned rows */
    2900:   PetscCall(PetscLayoutMapLocal(A->rmap, n, rows, &len, &lrows, NULL));
    2901:   /* fix right-hand side if needed */
    2902:   if (x && b) {
    2903:     const PetscScalar *xx;
    2904:     PetscScalar       *bb;

    2906:     if (columns) {
    2907:       /* Subtract the column contributions: b[i] -= A[i,r] * x[r] for non-zeroed rows i.
    2908:          Build x_zeroed = x restricted to the zeroed (locally owned) rows, 0 elsewhere,
    2909:          then compute b -= A * x_zeroed using the original (unmodified) local matrices.
    2910:          The zeroed rows of b are overwritten below with diag * x[r], so no special
    2911:          treatment is needed for them in the MatMult() output. */
    2912:       Vec          x_zeroed, temp;
    2913:       PetscScalar *xz, *saved;

    2915:       PetscCall(PetscMalloc1(len, &saved));
    2916:       PetscCall(VecDuplicate(x, &x_zeroed));
    2917:       PetscCall(VecGetArrayRead(x, &xx));
    2918:       PetscCall(VecGetArray(x_zeroed, &xz));
    2919:       for (PetscInt i = 0; i < len; i++) {
    2920:         xz[lrows[i]] = xx[lrows[i]];
    2921:         saved[i]     = diag * xx[lrows[i]];
    2922:       }
    2923:       PetscCall(VecRestoreArray(x_zeroed, &xz));
    2924:       PetscCall(VecRestoreArrayRead(x, &xx));
    2925:       PetscCall(VecDuplicate(b, &temp));
    2926:       PetscCall(MatMult(A, x_zeroed, temp));
    2927:       PetscCall(VecAXPY(b, -1.0, temp));
    2928:       PetscCall(VecDestroy(&temp));
    2929:       PetscCall(VecDestroy(&x_zeroed));
    2930:       /* Overwrite zeroed rows: b[r] = diag * x[r] (after the MatMult() so it is not clobbered) */
    2931:       PetscCall(VecGetArray(b, &bb));
    2932:       for (PetscInt i = 0; i < len; i++) bb[lrows[i]] = saved[i];
    2933:       PetscCall(VecRestoreArray(b, &bb));
    2934:       PetscCall(PetscFree(saved));
    2935:     } else {
    2936:       /* MatZeroRows(): only set b[r] = diag * x[r] for the zeroed rows */
    2937:       PetscCall(VecGetArrayRead(x, &xx));
    2938:       PetscCall(VecGetArray(b, &bb));
    2939:       for (PetscInt i = 0; i < len; ++i) bb[lrows[i]] = diag * xx[lrows[i]];
    2940:       PetscCall(VecRestoreArrayRead(x, &xx));
    2941:       PetscCall(VecRestoreArray(b, &bb));
    2942:     }
    2943:   }
    2944:   /* get rows associated to the local matrices */
    2945:   PetscCall(PetscArrayzero(matis->sf_leafdata, nl));
    2946:   PetscCall(PetscArrayzero(matis->sf_rootdata, A->rmap->n));
    2947:   for (PetscInt i = 0; i < len; i++) matis->sf_rootdata[lrows[i]] = 1;
    2948:   PetscCall(PetscFree(lrows));
    2949:   PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
    2950:   PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
    2951:   PetscCall(PetscMalloc1(nl, &lrows));
    2952:   nr = 0;
    2953:   for (PetscInt i = 0; i < nl; i++)
    2954:     if (matis->sf_leafdata[i]) lrows[nr++] = i;
    2955:   PetscCall(MatISZeroRowsColumnsLocal_Private(A, nr, lrows, diag, columns));
    2956:   PetscCall(PetscFree(lrows));
    2957:   PetscFunctionReturn(PETSC_SUCCESS);
    2958: }

    2960: static PetscErrorCode MatZeroRows_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
    2961: {
    2962:   PetscFunctionBegin;
    2963:   PetscCall(MatZeroRowsColumns_Private_IS(A, n, rows, diag, x, b, PETSC_FALSE));
    2964:   PetscFunctionReturn(PETSC_SUCCESS);
    2965: }

    2967: static PetscErrorCode MatZeroRowsColumns_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
    2968: {
    2969:   PetscFunctionBegin;
    2970:   PetscCall(MatZeroRowsColumns_Private_IS(A, n, rows, diag, x, b, PETSC_TRUE));
    2971:   PetscFunctionReturn(PETSC_SUCCESS);
    2972: }

    2974: static PetscErrorCode MatAssemblyBegin_IS(Mat A, MatAssemblyType type)
    2975: {
    2976:   Mat_IS *is = (Mat_IS *)A->data;

    2978:   PetscFunctionBegin;
    2979:   PetscCall(MatAssemblyBegin(is->A, type));
    2980:   PetscFunctionReturn(PETSC_SUCCESS);
    2981: }

    2983: static PetscErrorCode MatAssemblyEnd_IS(Mat A, MatAssemblyType type)
    2984: {
    2985:   Mat_IS   *is = (Mat_IS *)A->data;
    2986:   PetscBool lnnz;

    2988:   PetscFunctionBegin;
    2989:   PetscCall(MatAssemblyEnd(is->A, type));
    2990:   /* fix for local empty rows/cols */
    2991:   if (is->locempty && type == MAT_FINAL_ASSEMBLY) {
    2992:     Mat                    newlA;
    2993:     ISLocalToGlobalMapping rl2g, cl2g;
    2994:     IS                     nzr, nzc;
    2995:     PetscInt               nr, nc, nnzr, nnzc;
    2996:     PetscBool              lnewl2g, newl2g;

    2998:     PetscCall(MatGetSize(is->A, &nr, &nc));
    2999:     PetscCall(MatFindNonzeroRowsOrCols_Basic(is->A, PETSC_FALSE, PETSC_SMALL, &nzr));
    3000:     if (!nzr) PetscCall(ISCreateStride(PetscObjectComm((PetscObject)is->A), nr, 0, 1, &nzr));
    3001:     PetscCall(MatFindNonzeroRowsOrCols_Basic(is->A, PETSC_TRUE, PETSC_SMALL, &nzc));
    3002:     if (!nzc) PetscCall(ISCreateStride(PetscObjectComm((PetscObject)is->A), nc, 0, 1, &nzc));
    3003:     PetscCall(ISGetSize(nzr, &nnzr));
    3004:     PetscCall(ISGetSize(nzc, &nnzc));
    3005:     if (nnzr != nr || nnzc != nc) { /* need new global l2g map */
    3006:       lnewl2g = PETSC_TRUE;
    3007:       PetscCallMPI(MPIU_Allreduce(&lnewl2g, &newl2g, 1, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A)));

    3009:       /* extract valid submatrix */
    3010:       PetscCall(MatCreateSubMatrix(is->A, nzr, nzc, MAT_INITIAL_MATRIX, &newlA));
    3011:     } else { /* local matrix fully populated */
    3012:       lnewl2g = PETSC_FALSE;
    3013:       PetscCallMPI(MPIU_Allreduce(&lnewl2g, &newl2g, 1, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A)));
    3014:       PetscCall(PetscObjectReference((PetscObject)is->A));
    3015:       newlA = is->A;
    3016:     }

    3018:     /* attach new global l2g map if needed */
    3019:     if (newl2g) {
    3020:       IS              zr, zc;
    3021:       const PetscInt *ridxs, *cidxs, *zridxs, *zcidxs;
    3022:       PetscInt       *nidxs, i;

    3024:       PetscCall(ISComplement(nzr, 0, nr, &zr));
    3025:       PetscCall(ISComplement(nzc, 0, nc, &zc));
    3026:       PetscCall(PetscMalloc1(PetscMax(nr, nc), &nidxs));
    3027:       PetscCall(ISLocalToGlobalMappingGetIndices(is->rmapping, &ridxs));
    3028:       PetscCall(ISLocalToGlobalMappingGetIndices(is->cmapping, &cidxs));
    3029:       PetscCall(ISGetIndices(zr, &zridxs));
    3030:       PetscCall(ISGetIndices(zc, &zcidxs));
    3031:       PetscCall(ISGetLocalSize(zr, &nnzr));
    3032:       PetscCall(ISGetLocalSize(zc, &nnzc));

    3034:       PetscCall(PetscArraycpy(nidxs, ridxs, nr));
    3035:       for (i = 0; i < nnzr; i++) nidxs[zridxs[i]] = -1;
    3036:       PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A), 1, nr, nidxs, PETSC_COPY_VALUES, &rl2g));
    3037:       PetscCall(PetscArraycpy(nidxs, cidxs, nc));
    3038:       for (i = 0; i < nnzc; i++) nidxs[zcidxs[i]] = -1;
    3039:       PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A), 1, nc, nidxs, PETSC_COPY_VALUES, &cl2g));

    3041:       PetscCall(ISRestoreIndices(zr, &zridxs));
    3042:       PetscCall(ISRestoreIndices(zc, &zcidxs));
    3043:       PetscCall(ISLocalToGlobalMappingRestoreIndices(is->rmapping, &ridxs));
    3044:       PetscCall(ISLocalToGlobalMappingRestoreIndices(is->cmapping, &cidxs));
    3045:       PetscCall(ISDestroy(&nzr));
    3046:       PetscCall(ISDestroy(&nzc));
    3047:       PetscCall(ISDestroy(&zr));
    3048:       PetscCall(ISDestroy(&zc));
    3049:       PetscCall(PetscFree(nidxs));
    3050:       PetscCall(MatSetLocalToGlobalMapping(A, rl2g, cl2g));
    3051:       PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
    3052:       PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
    3053:     }
    3054:     PetscCall(MatISSetLocalMat(A, newlA));
    3055:     PetscCall(MatDestroy(&newlA));
    3056:     PetscCall(ISDestroy(&nzr));
    3057:     PetscCall(ISDestroy(&nzc));
    3058:     is->locempty = PETSC_FALSE;
    3059:   }
    3060:   lnnz          = (PetscBool)(is->A->nonzerostate == is->lnnzstate);
    3061:   is->lnnzstate = is->A->nonzerostate;
    3062:   PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &lnnz, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
    3063:   if (!lnnz) A->nonzerostate++;
    3064:   PetscFunctionReturn(PETSC_SUCCESS);
    3065: }

    3067: static PetscErrorCode MatISGetLocalMat_IS(Mat mat, Mat *local)
    3068: {
    3069:   Mat_IS *is = (Mat_IS *)mat->data;

    3071:   PetscFunctionBegin;
    3072:   *local = is->A;
    3073:   PetscFunctionReturn(PETSC_SUCCESS);
    3074: }

    3076: static PetscErrorCode MatISRestoreLocalMat_IS(Mat mat, Mat *local)
    3077: {
    3078:   PetscFunctionBegin;
    3079:   *local = NULL;
    3080:   PetscFunctionReturn(PETSC_SUCCESS);
    3081: }

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

    3086:   Not Collective.

    3088:   Input Parameter:
    3089: . mat - the matrix

    3091:   Output Parameter:
    3092: . local - the local matrix

    3094:   Level: intermediate

    3096:   Notes:
    3097:   This can be called if you have precomputed the nonzero structure of the
    3098:   matrix and want to provide it to the inner matrix object to improve the performance
    3099:   of the `MatSetValues()` operation.

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

    3103: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatISRestoreLocalMat()`
    3104: @*/
    3105: PetscErrorCode MatISGetLocalMat(Mat mat, Mat *local)
    3106: {
    3107:   PetscFunctionBegin;
    3109:   PetscAssertPointer(local, 2);
    3110:   PetscUseMethod(mat, "MatISGetLocalMat_C", (Mat, Mat *), (mat, local));
    3111:   PetscFunctionReturn(PETSC_SUCCESS);
    3112: }

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

    3117:   Not Collective.

    3119:   Input Parameters:
    3120: + mat   - the matrix
    3121: - local - the local matrix

    3123:   Level: intermediate

    3125: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatISGetLocalMat()`
    3126: @*/
    3127: PetscErrorCode MatISRestoreLocalMat(Mat mat, Mat *local)
    3128: {
    3129:   PetscFunctionBegin;
    3131:   PetscAssertPointer(local, 2);
    3132:   PetscUseMethod(mat, "MatISRestoreLocalMat_C", (Mat, Mat *), (mat, local));
    3133:   PetscFunctionReturn(PETSC_SUCCESS);
    3134: }

    3136: static PetscErrorCode MatISSetLocalMatType_IS(Mat mat, MatType mtype)
    3137: {
    3138:   Mat_IS *is = (Mat_IS *)mat->data;

    3140:   PetscFunctionBegin;
    3141:   if (is->A) PetscCall(MatSetType(is->A, mtype));
    3142:   PetscCall(PetscFree(is->lmattype));
    3143:   PetscCall(PetscStrallocpy(mtype, &is->lmattype));
    3144:   PetscFunctionReturn(PETSC_SUCCESS);
    3145: }

    3147: /*@
    3148:   MatISSetLocalMatType - Specifies the type of local matrix inside the `MATIS`

    3150:   Logically Collective.

    3152:   Input Parameters:
    3153: + mat   - the matrix
    3154: - mtype - the local matrix type

    3156:   Level: intermediate

    3158: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatSetType()`, `MatType`
    3159: @*/
    3160: PetscErrorCode MatISSetLocalMatType(Mat mat, MatType mtype)
    3161: {
    3162:   PetscFunctionBegin;
    3164:   PetscUseMethod(mat, "MatISSetLocalMatType_C", (Mat, MatType), (mat, mtype));
    3165:   PetscFunctionReturn(PETSC_SUCCESS);
    3166: }

    3168: static PetscErrorCode MatISSetLocalMat_IS(Mat mat, Mat local)
    3169: {
    3170:   Mat_IS   *is = (Mat_IS *)mat->data;
    3171:   PetscInt  nrows, ncols, orows, ocols;
    3172:   MatType   mtype, otype;
    3173:   PetscBool sametype = PETSC_TRUE;

    3175:   PetscFunctionBegin;
    3176:   if (is->A && !is->islocalref) {
    3177:     PetscCall(MatGetSize(is->A, &orows, &ocols));
    3178:     PetscCall(MatGetSize(local, &nrows, &ncols));
    3179:     PetscCheck(orows == nrows && ocols == ncols, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local MATIS matrix should be of size %" PetscInt_FMT "x%" PetscInt_FMT " (passed a %" PetscInt_FMT "x%" PetscInt_FMT " matrix)", orows, ocols, nrows, ncols);
    3180:     PetscCall(MatGetType(local, &mtype));
    3181:     PetscCall(MatGetType(is->A, &otype));
    3182:     PetscCall(PetscStrcmp(mtype, otype, &sametype));
    3183:   }
    3184:   PetscCall(PetscObjectReference((PetscObject)local));
    3185:   PetscCall(MatDestroy(&is->A));
    3186:   is->A = local;
    3187:   PetscCall(MatGetType(is->A, &mtype));
    3188:   PetscCall(MatISSetLocalMatType(mat, mtype));
    3189:   if (!sametype && !is->islocalref) PetscCall(MatISSetUpScatters_Private(mat));
    3190:   is->lnnzstate = 0;
    3191:   PetscFunctionReturn(PETSC_SUCCESS);
    3192: }

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

    3197:   Not Collective

    3199:   Input Parameters:
    3200: + mat   - the matrix
    3201: - local - the local matrix

    3203:   Level: intermediate

    3205: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatISSetLocalMatType`, `MatISGetLocalMat()`
    3206: @*/
    3207: PetscErrorCode MatISSetLocalMat(Mat mat, Mat local)
    3208: {
    3209:   PetscFunctionBegin;
    3212:   PetscUseMethod(mat, "MatISSetLocalMat_C", (Mat, Mat), (mat, local));
    3213:   PetscFunctionReturn(PETSC_SUCCESS);
    3214: }

    3216: static PetscErrorCode MatZeroEntries_IS(Mat A)
    3217: {
    3218:   Mat_IS *a = (Mat_IS *)A->data;

    3220:   PetscFunctionBegin;
    3221:   PetscCall(MatZeroEntries(a->A));
    3222:   PetscFunctionReturn(PETSC_SUCCESS);
    3223: }

    3225: static PetscErrorCode MatScale_IS(Mat A, PetscScalar a)
    3226: {
    3227:   Mat_IS *is = (Mat_IS *)A->data;

    3229:   PetscFunctionBegin;
    3230:   PetscCall(MatScale(is->A, a));
    3231:   PetscFunctionReturn(PETSC_SUCCESS);
    3232: }

    3234: static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v)
    3235: {
    3236:   Mat_IS *is = (Mat_IS *)A->data;

    3238:   PetscFunctionBegin;
    3239:   /* get diagonal of the local matrix */
    3240:   PetscCall(MatGetDiagonal(is->A, is->y));

    3242:   /* scatter diagonal back into global vector */
    3243:   PetscCall(VecSet(v, 0));
    3244:   PetscCall(VecScatterBegin(is->rctx, is->y, v, ADD_VALUES, SCATTER_REVERSE));
    3245:   PetscCall(VecScatterEnd(is->rctx, is->y, v, ADD_VALUES, SCATTER_REVERSE));
    3246:   PetscFunctionReturn(PETSC_SUCCESS);
    3247: }

    3249: static PetscErrorCode MatSetOption_IS(Mat A, MatOption op, PetscBool flg)
    3250: {
    3251:   Mat_IS *a = (Mat_IS *)A->data;

    3253:   PetscFunctionBegin;
    3254:   PetscCall(MatSetOption(a->A, op, flg));
    3255:   PetscFunctionReturn(PETSC_SUCCESS);
    3256: }

    3258: static PetscErrorCode MatAXPY_IS(Mat Y, PetscScalar a, Mat X, MatStructure str)
    3259: {
    3260:   Mat_IS *y = (Mat_IS *)Y->data;
    3261:   Mat_IS *x;

    3263:   PetscFunctionBegin;
    3264:   if (PetscDefined(USE_DEBUG)) {
    3265:     PetscBool ismatis;
    3266:     PetscCall(PetscObjectTypeCompare((PetscObject)X, MATIS, &ismatis));
    3267:     PetscCheck(ismatis, PetscObjectComm((PetscObject)Y), PETSC_ERR_SUP, "Cannot call MatAXPY(Y,a,X,str) with X not of type MATIS");
    3268:   }
    3269:   x = (Mat_IS *)X->data;
    3270:   PetscCall(MatAXPY(y->A, a, x->A, str));
    3271:   PetscFunctionReturn(PETSC_SUCCESS);
    3272: }

    3274: static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A, IS row, IS col, Mat *submat)
    3275: {
    3276:   Mat                    lA;
    3277:   Mat_IS                *matis = (Mat_IS *)A->data;
    3278:   ISLocalToGlobalMapping rl2g, cl2g;
    3279:   IS                     is;
    3280:   const PetscInt        *rg, *rl;
    3281:   PetscInt               nrg, rbs, cbs;
    3282:   PetscInt               N, M, nrl, i, *idxs;

    3284:   PetscFunctionBegin;
    3285:   PetscCall(ISGetBlockSize(row, &rbs));
    3286:   PetscCall(ISGetBlockSize(col, &cbs));
    3287:   PetscCall(ISLocalToGlobalMappingGetIndices(A->rmap->mapping, &rg));
    3288:   PetscCall(ISGetLocalSize(row, &nrl));
    3289:   PetscCall(ISGetIndices(row, &rl));
    3290:   PetscCall(ISLocalToGlobalMappingGetSize(A->rmap->mapping, &nrg));
    3291:   if (PetscDefined(USE_DEBUG)) {
    3292:     for (i = 0; i < nrl; i++) PetscCheck(rl[i] < nrg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Local row index %" PetscInt_FMT " -> %" PetscInt_FMT " greater than maximum possible %" PetscInt_FMT, i, rl[i], nrg);
    3293:   }
    3294:   if (nrg % rbs) nrg = rbs * (nrg / rbs + 1);
    3295:   PetscCall(PetscMalloc1(nrg, &idxs));
    3296:   /* map from [0,nrl) to row */
    3297:   for (i = 0; i < nrl; i++) idxs[i] = rl[i];
    3298:   for (i = nrl; i < nrg; i++) idxs[i] = -1;
    3299:   PetscCall(ISRestoreIndices(row, &rl));
    3300:   PetscCall(ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping, &rg));
    3301:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A), nrg, idxs, PETSC_OWN_POINTER, &is));
    3302:   PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g));
    3303:   PetscCall(ISLocalToGlobalMappingSetBlockSize(rl2g, rbs));
    3304:   PetscCall(ISDestroy(&is));
    3305:   /* compute new l2g map for columns */
    3306:   if (col != row || matis->rmapping != matis->cmapping || matis->A->rmap->mapping != matis->A->cmap->mapping) {
    3307:     const PetscInt *cg, *cl;
    3308:     PetscInt        ncg;
    3309:     PetscInt        ncl;

    3311:     PetscCall(ISLocalToGlobalMappingGetIndices(A->cmap->mapping, &cg));
    3312:     PetscCall(ISGetLocalSize(col, &ncl));
    3313:     PetscCall(ISGetIndices(col, &cl));
    3314:     PetscCall(ISLocalToGlobalMappingGetSize(A->cmap->mapping, &ncg));
    3315:     if (PetscDefined(USE_DEBUG)) {
    3316:       for (i = 0; i < ncl; i++) PetscCheck(cl[i] < ncg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Local column index %" PetscInt_FMT " -> %" PetscInt_FMT " greater than maximum possible %" PetscInt_FMT, i, cl[i], ncg);
    3317:     }
    3318:     if (ncg % cbs) ncg = cbs * (ncg / cbs + 1);
    3319:     PetscCall(PetscMalloc1(ncg, &idxs));
    3320:     /* map from [0,ncl) to col */
    3321:     for (i = 0; i < ncl; i++) idxs[i] = cl[i];
    3322:     for (i = ncl; i < ncg; i++) idxs[i] = -1;
    3323:     PetscCall(ISRestoreIndices(col, &cl));
    3324:     PetscCall(ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping, &cg));
    3325:     PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A), ncg, idxs, PETSC_OWN_POINTER, &is));
    3326:     PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g));
    3327:     PetscCall(ISLocalToGlobalMappingSetBlockSize(cl2g, cbs));
    3328:     PetscCall(ISDestroy(&is));
    3329:   } else {
    3330:     PetscCall(PetscObjectReference((PetscObject)rl2g));
    3331:     cl2g = rl2g;
    3332:   }

    3334:   /* create the MATIS submatrix */
    3335:   PetscCall(MatGetSize(A, &M, &N));
    3336:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), submat));
    3337:   PetscCall(MatSetSizes(*submat, PETSC_DECIDE, PETSC_DECIDE, M, N));
    3338:   PetscCall(MatSetType(*submat, MATIS));
    3339:   matis             = (Mat_IS *)((*submat)->data);
    3340:   matis->islocalref = A;
    3341:   PetscCall(MatSetLocalToGlobalMapping(*submat, rl2g, cl2g));
    3342:   PetscCall(MatISGetLocalMat(A, &lA));
    3343:   PetscCall(MatISSetLocalMat(*submat, lA));
    3344:   PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
    3345:   PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));

    3347:   /* remove unsupported ops */
    3348:   PetscCall(PetscMemzero((*submat)->ops, sizeof(struct _MatOps)));
    3349:   (*submat)->ops->destroy               = MatDestroy_IS;
    3350:   (*submat)->ops->setvalueslocal        = MatSetValuesLocal_SubMat_IS;
    3351:   (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS;
    3352:   (*submat)->ops->zerorowslocal         = MatZeroRowsLocal_SubMat_IS;
    3353:   (*submat)->ops->zerorowscolumnslocal  = MatZeroRowsColumnsLocal_SubMat_IS;
    3354:   (*submat)->ops->getlocalsubmatrix     = MatGetLocalSubMatrix_IS;
    3355:   PetscFunctionReturn(PETSC_SUCCESS);
    3356: }

    3358: static PetscErrorCode MatSetFromOptions_IS(Mat A, PetscOptionItems PetscOptionsObject)
    3359: {
    3360:   Mat_IS   *a = (Mat_IS *)A->data;
    3361:   char      type[256];
    3362:   PetscBool flg;

    3364:   PetscFunctionBegin;
    3365:   PetscOptionsHeadBegin(PetscOptionsObject, "MATIS options");
    3366:   PetscCall(PetscOptionsDeprecated("-matis_keepassembled", "-mat_is_keepassembled", "3.21", NULL));
    3367:   PetscCall(PetscOptionsDeprecated("-matis_fixempty", "-mat_is_fixempty", "3.21", NULL));
    3368:   PetscCall(PetscOptionsDeprecated("-matis_storel2l", "-mat_is_storel2l", "3.21", NULL));
    3369:   PetscCall(PetscOptionsDeprecated("-matis_localmat_type", "-mat_is_localmat_type", "3.21", NULL));
    3370:   PetscCall(PetscOptionsBool("-mat_is_keepassembled", "Store an assembled version if needed", NULL, a->keepassembled, &a->keepassembled, NULL));
    3371:   PetscCall(PetscOptionsBool("-mat_is_fixempty", "Fix local matrices in case of empty local rows/columns", "MatISFixLocalEmpty", a->locempty, &a->locempty, NULL));
    3372:   PetscCall(PetscOptionsBool("-mat_is_storel2l", "Store local-to-local matrices generated from PtAP operations", "MatISStoreL2L", a->storel2l, &a->storel2l, NULL));
    3373:   PetscCall(PetscOptionsBool("-mat_is_allow_repeated", "Allow local repeated entries", "MatISSetAllowRepeated", a->allow_repeated, &a->allow_repeated, NULL));
    3374:   PetscCall(PetscOptionsFList("-mat_is_localmat_type", "Matrix type", "MatISSetLocalMatType", MatList, a->lmattype, type, 256, &flg));
    3375:   if (flg) PetscCall(MatISSetLocalMatType(A, type));
    3376:   if (a->A) PetscCall(MatSetFromOptions(a->A));
    3377:   PetscOptionsHeadEnd();
    3378:   PetscFunctionReturn(PETSC_SUCCESS);
    3379: }

    3381: /*@
    3382:   MatCreateIS - Creates a "process" unassembled matrix.

    3384:   Collective.

    3386:   Input Parameters:
    3387: + comm - MPI communicator that will share the matrix
    3388: . bs   - block size of the matrix
    3389: . m    - local size of left vector used in matrix vector products
    3390: . n    - local size of right vector used in matrix vector products
    3391: . M    - global size of left vector used in matrix vector products
    3392: . N    - global size of right vector used in matrix vector products
    3393: . rmap - local to global map for rows
    3394: - cmap - local to global map for cols

    3396:   Output Parameter:
    3397: . A - the resulting matrix

    3399:   Level: intermediate

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

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

    3407: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatSetLocalToGlobalMapping()`
    3408: @*/
    3409: PetscErrorCode MatCreateIS(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt M, PetscInt N, ISLocalToGlobalMapping rmap, ISLocalToGlobalMapping cmap, Mat *A)
    3410: {
    3411:   PetscFunctionBegin;
    3412:   PetscCall(MatCreate(comm, A));
    3413:   PetscCall(MatSetSizes(*A, m, n, M, N));
    3414:   if (bs > 0) PetscCall(MatSetBlockSize(*A, bs));
    3415:   PetscCall(MatSetType(*A, MATIS));
    3416:   PetscCall(MatSetLocalToGlobalMapping(*A, rmap, cmap));
    3417:   PetscFunctionReturn(PETSC_SUCCESS);
    3418: }

    3420: static PetscErrorCode MatHasOperation_IS(Mat A, MatOperation op, PetscBool *has)
    3421: {
    3422:   Mat_IS      *a              = (Mat_IS *)A->data;
    3423:   MatOperation tobefiltered[] = {MATOP_MULT_ADD, MATOP_MULT_TRANSPOSE_ADD, MATOP_GET_DIAGONAL_BLOCK, MATOP_INCREASE_OVERLAP};

    3425:   PetscFunctionBegin;
    3426:   *has = PETSC_FALSE;
    3427:   if (!((void **)A->ops)[op] || !a->A) PetscFunctionReturn(PETSC_SUCCESS);
    3428:   *has = PETSC_TRUE;
    3429:   for (PetscInt i = 0; i < (PetscInt)PETSC_STATIC_ARRAY_LENGTH(tobefiltered); i++)
    3430:     if (op == tobefiltered[i]) PetscFunctionReturn(PETSC_SUCCESS);
    3431:   PetscCall(MatHasOperation(a->A, op, has));
    3432:   PetscFunctionReturn(PETSC_SUCCESS);
    3433: }

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

    3439:   PetscFunctionBegin;
    3440:   PetscCall(MatSetValuesCOO(a->A, v, imode));
    3441:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    3442:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    3443:   PetscFunctionReturn(PETSC_SUCCESS);
    3444: }

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

    3450:   PetscFunctionBegin;
    3451:   PetscCheck(a->A, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to provide l2g map first via MatSetLocalToGlobalMapping");
    3452:   if (a->A->rmap->mapping || a->A->cmap->mapping) {
    3453:     PetscCall(MatSetPreallocationCOOLocal(a->A, ncoo, coo_i, coo_j));
    3454:   } else {
    3455:     PetscCall(MatSetPreallocationCOO(a->A, ncoo, coo_i, coo_j));
    3456:   }
    3457:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", MatSetValuesCOO_IS));
    3458:   A->preallocated = PETSC_TRUE;
    3459:   PetscFunctionReturn(PETSC_SUCCESS);
    3460: }

    3462: static PetscErrorCode MatSetPreallocationCOO_IS(Mat A, PetscCount ncoo, PetscInt coo_i[], PetscInt coo_j[])
    3463: {
    3464:   Mat_IS  *a = (Mat_IS *)A->data;
    3465:   PetscInt ncoo_i;

    3467:   PetscFunctionBegin;
    3468:   PetscCheck(a->A, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to provide l2g map first via MatSetLocalToGlobalMapping");
    3469:   PetscCall(PetscIntCast(ncoo, &ncoo_i));
    3470:   PetscCall(ISGlobalToLocalMappingApply(a->rmapping, IS_GTOLM_MASK, ncoo_i, coo_i, NULL, coo_i));
    3471:   PetscCall(ISGlobalToLocalMappingApply(a->cmapping, IS_GTOLM_MASK, ncoo_i, coo_j, NULL, coo_j));
    3472:   PetscCall(MatSetPreallocationCOO(a->A, ncoo, coo_i, coo_j));
    3473:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", MatSetValuesCOO_IS));
    3474:   A->preallocated = PETSC_TRUE;
    3475:   PetscFunctionReturn(PETSC_SUCCESS);
    3476: }

    3478: static PetscErrorCode MatISGetAssembled_Private(Mat A, Mat *tA)
    3479: {
    3480:   Mat_IS          *a = (Mat_IS *)A->data;
    3481:   PetscObjectState Astate, aAstate       = PETSC_INT_MIN;
    3482:   PetscObjectState Annzstate, aAnnzstate = PETSC_INT_MIN;

    3484:   PetscFunctionBegin;
    3485:   PetscCall(PetscObjectStateGet((PetscObject)A, &Astate));
    3486:   Annzstate = A->nonzerostate;
    3487:   if (a->assembledA) {
    3488:     PetscCall(PetscObjectStateGet((PetscObject)a->assembledA, &aAstate));
    3489:     aAnnzstate = a->assembledA->nonzerostate;
    3490:   }
    3491:   if (aAnnzstate != Annzstate) PetscCall(MatDestroy(&a->assembledA));
    3492:   if (Astate != aAstate || !a->assembledA) {
    3493:     MatType     aAtype;
    3494:     PetscMPIInt size;
    3495:     PetscInt    rbs, cbs, bs;

    3497:     /* the assembled form is used as temporary storage for parallel operations
    3498:        like createsubmatrices and the like, do not waste device memory */
    3499:     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
    3500:     PetscCall(ISLocalToGlobalMappingGetBlockSize(a->cmapping, &cbs));
    3501:     PetscCall(ISLocalToGlobalMappingGetBlockSize(a->rmapping, &rbs));
    3502:     bs = rbs == cbs ? rbs : 1;
    3503:     if (a->assembledA) PetscCall(MatGetType(a->assembledA, &aAtype));
    3504:     else if (size > 1) aAtype = bs > 1 ? MATMPIBAIJ : MATMPIAIJ;
    3505:     else aAtype = bs > 1 ? MATSEQBAIJ : MATSEQAIJ;

    3507:     PetscCall(MatConvert(A, aAtype, a->assembledA ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX, &a->assembledA));
    3508:     PetscCall(PetscObjectStateSet((PetscObject)a->assembledA, Astate));
    3509:     a->assembledA->nonzerostate = Annzstate;
    3510:   }
    3511:   PetscCall(PetscObjectReference((PetscObject)a->assembledA));
    3512:   *tA = a->assembledA;
    3513:   if (!a->keepassembled) PetscCall(MatDestroy(&a->assembledA));
    3514:   PetscFunctionReturn(PETSC_SUCCESS);
    3515: }

    3517: static PetscErrorCode MatISRestoreAssembled_Private(Mat A, Mat *tA)
    3518: {
    3519:   PetscFunctionBegin;
    3520:   PetscCall(MatDestroy(tA));
    3521:   *tA = NULL;
    3522:   PetscFunctionReturn(PETSC_SUCCESS);
    3523: }

    3525: static PetscErrorCode MatGetDiagonalBlock_IS(Mat A, Mat *dA)
    3526: {
    3527:   Mat_IS          *a = (Mat_IS *)A->data;
    3528:   PetscObjectState Astate, dAstate = PETSC_INT_MIN;

    3530:   PetscFunctionBegin;
    3531:   PetscCall(PetscObjectStateGet((PetscObject)A, &Astate));
    3532:   if (a->dA) PetscCall(PetscObjectStateGet((PetscObject)a->dA, &dAstate));
    3533:   if (Astate != dAstate) {
    3534:     Mat     tA;
    3535:     MatType ltype;

    3537:     PetscCall(MatDestroy(&a->dA));
    3538:     PetscCall(MatISGetAssembled_Private(A, &tA));
    3539:     PetscCall(MatGetDiagonalBlock(tA, &a->dA));
    3540:     PetscCall(MatPropagateSymmetryOptions(tA, a->dA));
    3541:     PetscCall(MatGetType(a->A, &ltype));
    3542:     PetscCall(MatConvert(a->dA, ltype, MAT_INPLACE_MATRIX, &a->dA));
    3543:     PetscCall(PetscObjectReference((PetscObject)a->dA));
    3544:     PetscCall(MatISRestoreAssembled_Private(A, &tA));
    3545:     PetscCall(PetscObjectStateSet((PetscObject)a->dA, Astate));
    3546:   }
    3547:   *dA = a->dA;
    3548:   PetscFunctionReturn(PETSC_SUCCESS);
    3549: }

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

    3555:   PetscFunctionBegin;
    3556:   PetscCall(MatISGetAssembled_Private(A, &tA));
    3557:   PetscCall(MatCreateSubMatrices(tA, n, irow, icol, reuse, submat));
    3558:   /* MatCreateSubMatrices_MPIAIJ is a mess at the moment */
    3559: #if 0
    3560:   {
    3561:     Mat_IS    *a = (Mat_IS*)A->data;
    3562:     MatType   ltype;
    3563:     VecType   vtype;
    3564:     char      *flg;

    3566:     PetscCall(MatGetType(a->A,&ltype));
    3567:     PetscCall(MatGetVecType(a->A,&vtype));
    3568:     PetscCall(PetscStrstr(vtype,"cuda",&flg));
    3569:     if (!flg) PetscCall(PetscStrstr(vtype,"hip",&flg));
    3570:     if (!flg) PetscCall(PetscStrstr(vtype,"kokkos",&flg));
    3571:     if (flg) {
    3572:       for (PetscInt i = 0; i < n; i++) {
    3573:         Mat sA = (*submat)[i];

    3575:         PetscCall(MatConvert(sA,ltype,MAT_INPLACE_MATRIX,&sA));
    3576:         (*submat)[i] = sA;
    3577:       }
    3578:     }
    3579:   }
    3580: #endif
    3581:   PetscCall(MatISRestoreAssembled_Private(A, &tA));
    3582:   PetscFunctionReturn(PETSC_SUCCESS);
    3583: }

    3585: static PetscErrorCode MatIncreaseOverlap_IS(Mat A, PetscInt n, IS is[], PetscInt ov)
    3586: {
    3587:   Mat tA;

    3589:   PetscFunctionBegin;
    3590:   PetscCall(MatISGetAssembled_Private(A, &tA));
    3591:   PetscCall(MatIncreaseOverlap(tA, n, is, ov));
    3592:   PetscCall(MatISRestoreAssembled_Private(A, &tA));
    3593:   PetscFunctionReturn(PETSC_SUCCESS);
    3594: }

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

    3599:   Not Collective

    3601:   Input Parameter:
    3602: . A - the matrix

    3604:   Output Parameters:
    3605: + rmapping - row mapping
    3606: - cmapping - column mapping

    3608:   Level: advanced

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

    3613: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatSetLocalToGlobalMapping()`
    3614: @*/
    3615: PetscErrorCode MatISGetLocalToGlobalMapping(Mat A, ISLocalToGlobalMapping *rmapping, ISLocalToGlobalMapping *cmapping)
    3616: {
    3617:   PetscFunctionBegin;
    3620:   if (rmapping) PetscAssertPointer(rmapping, 2);
    3621:   if (cmapping) PetscAssertPointer(cmapping, 3);
    3622:   PetscUseMethod(A, "MatISGetLocalToGlobalMapping_C", (Mat, ISLocalToGlobalMapping *, ISLocalToGlobalMapping *), (A, rmapping, cmapping));
    3623:   PetscFunctionReturn(PETSC_SUCCESS);
    3624: }

    3626: static PetscErrorCode MatISGetLocalToGlobalMapping_IS(Mat A, ISLocalToGlobalMapping *r, ISLocalToGlobalMapping *c)
    3627: {
    3628:   Mat_IS *a = (Mat_IS *)A->data;

    3630:   PetscFunctionBegin;
    3631:   if (r) *r = a->rmapping;
    3632:   if (c) *c = a->cmapping;
    3633:   PetscFunctionReturn(PETSC_SUCCESS);
    3634: }

    3636: static PetscErrorCode MatSetBlockSizes_IS(Mat A, PetscInt rbs, PetscInt cbs)
    3637: {
    3638:   Mat_IS *a = (Mat_IS *)A->data;

    3640:   PetscFunctionBegin;
    3641:   if (a->A) PetscCall(MatSetBlockSizes(a->A, rbs, cbs));
    3642:   PetscFunctionReturn(PETSC_SUCCESS);
    3643: }

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

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

    3655:   Level: intermediate

    3657:   Notes:
    3658:   Options prefix for the inner matrix are given by `-is_mat_xxx`

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

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

    3665: .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatISGetLocalMat()`, `MatSetLocalToGlobalMapping()`, `MatISSetPreallocation()`, `MatCreateIS()`, `PCBDDC`, `KSPFETIDP`
    3666: M*/
    3667: PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A)
    3668: {
    3669:   Mat_IS *a;

    3671:   PetscFunctionBegin;
    3672:   PetscCall(PetscNew(&a));
    3673:   PetscCall(PetscStrallocpy(MATAIJ, &a->lmattype));
    3674:   A->data = (void *)a;

    3676:   /* matrix ops */
    3677:   PetscCall(PetscMemzero(A->ops, sizeof(struct _MatOps)));
    3678:   A->ops->mult                    = MatMult_IS;
    3679:   A->ops->multadd                 = MatMultAdd_IS;
    3680:   A->ops->multtranspose           = MatMultTranspose_IS;
    3681:   A->ops->multtransposeadd        = MatMultTransposeAdd_IS;
    3682:   A->ops->destroy                 = MatDestroy_IS;
    3683:   A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS;
    3684:   A->ops->setvalues               = MatSetValues_IS;
    3685:   A->ops->setvaluesblocked        = MatSetValuesBlocked_IS;
    3686:   A->ops->setvalueslocal          = MatSetValuesLocal_IS;
    3687:   A->ops->setvaluesblockedlocal   = MatSetValuesBlockedLocal_IS;
    3688:   A->ops->zerorows                = MatZeroRows_IS;
    3689:   A->ops->zerorowscolumns         = MatZeroRowsColumns_IS;
    3690:   A->ops->assemblybegin           = MatAssemblyBegin_IS;
    3691:   A->ops->assemblyend             = MatAssemblyEnd_IS;
    3692:   A->ops->view                    = MatView_IS;
    3693:   A->ops->load                    = MatLoad_IS;
    3694:   A->ops->zeroentries             = MatZeroEntries_IS;
    3695:   A->ops->scale                   = MatScale_IS;
    3696:   A->ops->getdiagonal             = MatGetDiagonal_IS;
    3697:   A->ops->setoption               = MatSetOption_IS;
    3698:   A->ops->ishermitian             = MatIsHermitian_IS;
    3699:   A->ops->issymmetric             = MatIsSymmetric_IS;
    3700:   A->ops->isstructurallysymmetric = MatIsStructurallySymmetric_IS;
    3701:   A->ops->duplicate               = MatDuplicate_IS;
    3702:   A->ops->copy                    = MatCopy_IS;
    3703:   A->ops->getlocalsubmatrix       = MatGetLocalSubMatrix_IS;
    3704:   A->ops->createsubmatrix         = MatCreateSubMatrix_IS;
    3705:   A->ops->axpy                    = MatAXPY_IS;
    3706:   A->ops->diagonalset             = MatDiagonalSet_IS;
    3707:   A->ops->shift                   = MatShift_IS;
    3708:   A->ops->transpose               = MatTranspose_IS;
    3709:   A->ops->getinfo                 = MatGetInfo_IS;
    3710:   A->ops->diagonalscale           = MatDiagonalScale_IS;
    3711:   A->ops->setfromoptions          = MatSetFromOptions_IS;
    3712:   A->ops->setup                   = MatSetUp_IS;
    3713:   A->ops->hasoperation            = MatHasOperation_IS;
    3714:   A->ops->getdiagonalblock        = MatGetDiagonalBlock_IS;
    3715:   A->ops->createsubmatrices       = MatCreateSubMatrices_IS;
    3716:   A->ops->increaseoverlap         = MatIncreaseOverlap_IS;
    3717:   A->ops->setblocksizes           = MatSetBlockSizes_IS;

    3719:   /* special MATIS functions */
    3720:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMatType_C", MatISSetLocalMatType_IS));
    3721:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalMat_C", MatISGetLocalMat_IS));
    3722:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISRestoreLocalMat_C", MatISRestoreLocalMat_IS));
    3723:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMat_C", MatISSetLocalMat_IS));
    3724:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetPreallocation_C", MatISSetPreallocation_IS));
    3725:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetAllowRepeated_C", MatISSetAllowRepeated_IS));
    3726:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISStoreL2L_C", MatISStoreL2L_IS));
    3727:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISFixLocalEmpty_C", MatISFixLocalEmpty_IS));
    3728:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalToGlobalMapping_C", MatISGetLocalToGlobalMapping_IS));
    3729:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpiaij_C", MatConvert_IS_XAIJ));
    3730:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpibaij_C", MatConvert_IS_XAIJ));
    3731:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpisbaij_C", MatConvert_IS_XAIJ));
    3732:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqaij_C", MatConvert_IS_XAIJ));
    3733:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqbaij_C", MatConvert_IS_XAIJ));
    3734:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqsbaij_C", MatConvert_IS_XAIJ));
    3735:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_aij_C", MatConvert_IS_XAIJ));
    3736:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOOLocal_C", MatSetPreallocationCOOLocal_IS));
    3737:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOO_C", MatSetPreallocationCOO_IS));
    3738:   PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATIS));
    3739:   PetscFunctionReturn(PETSC_SUCCESS);
    3740: }