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 <../src/mat/impls/is/matis.h>
     11: #include <petsc/private/sfimpl.h>
     12: #include <petsc/private/vecimpl.h>

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

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

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

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

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

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

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

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

     97:   PetscObjectGetComm((PetscObject)PT,&comm);
     98:   bs   = 1;
     99:   PetscObjectBaseTypeCompare((PetscObject)PT,MATMPIAIJ,&ismpiaij);
    100:   PetscObjectBaseTypeCompare((PetscObject)PT,MATMPIBAIJ,&ismpibaij);
    101:   PetscObjectBaseTypeCompare((PetscObject)PT,MATSEQAIJ,&isseqaij);
    102:   PetscObjectTypeCompare((PetscObject)PT,MATSEQBAIJ,&isseqbaij);
    103:   if (isseqaij || isseqbaij) {
    104:     Pd = PT;
    105:     Po = NULL;
    106:     garray = NULL;
    107:   } else if (ismpiaij) {
    108:     MatMPIAIJGetSeqAIJ(PT,&Pd,&Po,&garray);
    109:   } else if (ismpibaij) {
    110:     MatMPIBAIJGetSeqBAIJ(PT,&Pd,&Po,&garray);
    111:     MatGetBlockSize(PT,&bs);
    112:   } else SETERRQ1(comm,PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)(PT))->type_name);

    114:   /* identify any null columns in Pd or Po */
    115:   /* We use a tolerance comparison since it may happen that, with geometric multigrid,
    116:      some of the columns are not really zero, but very close to */
    117:   zo = zd = NULL;
    118:   if (Po) {
    119:     MatFindNonzeroRowsOrCols_Basic(Po,PETSC_TRUE,PETSC_SMALL,&zo);
    120:   }
    121:   MatFindNonzeroRowsOrCols_Basic(Pd,PETSC_TRUE,PETSC_SMALL,&zd);

    123:   MatGetLocalSize(PT,NULL,&dc);
    124:   MatGetOwnershipRangeColumn(PT,&stc,NULL);
    125:   if (Po) { MatGetLocalSize(Po,NULL,&oc); }
    126:   else oc = 0;
    127:   PetscMalloc1((dc+oc)/bs,&aux);
    128:   if (zd) {
    129:     const PetscInt *idxs;
    130:     PetscInt       nz;

    132:     /* this will throw an error if bs is not valid */
    133:     ISSetBlockSize(zd,bs);
    134:     ISGetLocalSize(zd,&nz);
    135:     ISGetIndices(zd,&idxs);
    136:     ctd  = nz/bs;
    137:     for (i=0; i<ctd; i++) aux[i] = (idxs[bs*i]+stc)/bs;
    138:     ISRestoreIndices(zd,&idxs);
    139:   } else {
    140:     ctd = dc/bs;
    141:     for (i=0; i<ctd; i++) aux[i] = i+stc/bs;
    142:   }
    143:   if (zo) {
    144:     const PetscInt *idxs;
    145:     PetscInt       nz;

    147:     /* this will throw an error if bs is not valid */
    148:     ISSetBlockSize(zo,bs);
    149:     ISGetLocalSize(zo,&nz);
    150:     ISGetIndices(zo,&idxs);
    151:     cto  = nz/bs;
    152:     for (i=0; i<cto; i++) aux[i+ctd] = garray[idxs[bs*i]/bs];
    153:     ISRestoreIndices(zo,&idxs);
    154:   } else {
    155:     cto = oc/bs;
    156:     for (i=0; i<cto; i++) aux[i+ctd] = garray[i];
    157:   }
    158:   ISCreateBlock(comm,bs,ctd+cto,aux,PETSC_OWN_POINTER,cis);
    159:   ISDestroy(&zd);
    160:   ISDestroy(&zo);
    161:   return(0);
    162: }

    164: static PetscErrorCode MatPtAPSymbolic_IS_XAIJ(Mat A,Mat P,PetscReal fill,Mat C)
    165: {
    166:   Mat                    PT,lA;
    167:   MatISPtAP              ptap;
    168:   ISLocalToGlobalMapping Crl2g,Ccl2g,rl2g,cl2g;
    169:   PetscContainer         c;
    170:   MatType                lmtype;
    171:   const PetscInt         *garray;
    172:   PetscInt               ibs,N,dc;
    173:   MPI_Comm               comm;
    174:   PetscErrorCode         ierr;

    177:   PetscObjectGetComm((PetscObject)A,&comm);
    178:   MatSetType(C,MATIS);
    179:   MatISGetLocalMat(A,&lA);
    180:   MatGetType(lA,&lmtype);
    181:   MatISSetLocalMatType(C,lmtype);
    182:   MatGetSize(P,NULL,&N);
    183:   MatGetLocalSize(P,NULL,&dc);
    184:   MatSetSizes(C,dc,dc,N,N);
    185: /* Not sure about this
    186:   MatGetBlockSizes(P,NULL,&ibs);
    187:   MatSetBlockSize(*C,ibs);
    188: */

    190:   PetscNew(&ptap);
    191:   PetscContainerCreate(PETSC_COMM_SELF,&c);
    192:   PetscContainerSetPointer(c,ptap);
    193:   PetscContainerSetUserDestroy(c,MatISContainerDestroyPtAP_Private);
    194:   PetscObjectCompose((PetscObject)C,"_MatIS_PtAP",(PetscObject)c);
    195:   PetscContainerDestroy(&c);
    196:   ptap->fill = fill;

    198:   MatGetLocalToGlobalMapping(A,&rl2g,&cl2g);

    200:   ISLocalToGlobalMappingGetBlockSize(cl2g,&ibs);
    201:   ISLocalToGlobalMappingGetSize(cl2g,&N);
    202:   ISLocalToGlobalMappingGetBlockIndices(cl2g,&garray);
    203:   ISCreateBlock(comm,ibs,N/ibs,garray,PETSC_COPY_VALUES,&ptap->ris0);
    204:   ISLocalToGlobalMappingRestoreBlockIndices(cl2g,&garray);

    206:   MatCreateSubMatrix(P,ptap->ris0,NULL,MAT_INITIAL_MATRIX,&PT);
    207:   MatGetNonzeroColumnsLocal_Private(PT,&ptap->cis0);
    208:   ISLocalToGlobalMappingCreateIS(ptap->cis0,&Ccl2g);
    209:   MatDestroy(&PT);

    211:   Crl2g = NULL;
    212:   if (rl2g != cl2g) { /* unsymmetric A mapping */
    213:     PetscBool same,lsame = PETSC_FALSE;
    214:     PetscInt  N1,ibs1;

    216:     ISLocalToGlobalMappingGetSize(rl2g,&N1);
    217:     ISLocalToGlobalMappingGetBlockSize(rl2g,&ibs1);
    218:     ISLocalToGlobalMappingGetBlockIndices(rl2g,&garray);
    219:     ISCreateBlock(comm,ibs,N/ibs,garray,PETSC_COPY_VALUES,&ptap->ris1);
    220:     ISLocalToGlobalMappingRestoreBlockIndices(rl2g,&garray);
    221:     if (ibs1 == ibs && N1 == N) { /* check if the l2gmaps are the same */
    222:       const PetscInt *i1,*i2;

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

    248:   C->ops->ptapnumeric = MatPtAPNumeric_IS_XAIJ;
    249:   return(0);
    250: }

    252: /* ----------------------------------------- */
    253: static PetscErrorCode MatProductSymbolic_PtAP_IS_XAIJ(Mat C)
    254: {
    256:   Mat_Product    *product = C->product;
    257:   Mat            A=product->A,P=product->B;
    258:   PetscReal      fill=product->fill;

    261:   MatPtAPSymbolic_IS_XAIJ(A,P,fill,C);
    262:   C->ops->productnumeric = MatProductNumeric_PtAP;
    263:   return(0);
    264: }

    266: static PetscErrorCode MatProductSetFromOptions_IS_XAIJ_PtAP(Mat C)
    267: {
    269:   C->ops->productsymbolic = MatProductSymbolic_PtAP_IS_XAIJ;
    270:   return(0);
    271: }

    273: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_IS_XAIJ(Mat C)
    274: {
    276:   Mat_Product    *product = C->product;

    279:   if (product->type == MATPRODUCT_PtAP) {
    280:     MatProductSetFromOptions_IS_XAIJ_PtAP(C);
    281:   }
    282:   return(0);
    283: }

    285: /* ----------------------------------------- */
    286: static PetscErrorCode MatISContainerDestroyFields_Private(void *ptr)
    287: {
    288:   MatISLocalFields lf = (MatISLocalFields)ptr;
    289:   PetscInt         i;
    290:   PetscErrorCode   ierr;

    293:   for (i=0;i<lf->nr;i++) {
    294:     ISDestroy(&lf->rf[i]);
    295:   }
    296:   for (i=0;i<lf->nc;i++) {
    297:     ISDestroy(&lf->cf[i]);
    298:   }
    299:   PetscFree2(lf->rf,lf->cf);
    300:   PetscFree(lf);
    301:   return(0);
    302: }

    304: static PetscErrorCode MatConvert_SeqXAIJ_IS(Mat A,MatType type,MatReuse reuse,Mat *newmat)
    305: {
    306:   Mat            B,lB;

    310:   if (reuse != MAT_REUSE_MATRIX) {
    311:     ISLocalToGlobalMapping rl2g,cl2g;
    312:     PetscInt               bs;
    313:     IS                     is;

    315:     MatGetBlockSize(A,&bs);
    316:     ISCreateStride(PetscObjectComm((PetscObject)A),A->rmap->n/bs,0,1,&is);
    317:     if (bs > 1) {
    318:       IS       is2;
    319:       PetscInt i,*aux;

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

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

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

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

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

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

    406: typedef enum {MAT_IS_DISASSEMBLE_L2G_NATURAL,MAT_IS_DISASSEMBLE_L2G_MAT, MAT_IS_DISASSEMBLE_L2G_ND} MatISDisassemblel2gType;

    408: static PetscErrorCode MatMPIXAIJComputeLocalToGlobalMapping_Private(Mat A, ISLocalToGlobalMapping *l2g)
    409: {
    410:   Mat                     Ad,Ao;
    411:   IS                      is,ndmap,ndsub;
    412:   MPI_Comm                comm;
    413:   const PetscInt          *garray,*ndmapi;
    414:   PetscInt                bs,i,cnt,nl,*ncount,*ndmapc;
    415:   PetscBool               ismpiaij,ismpibaij;
    416:   const char *const       MatISDisassemblel2gTypes[] = {"NATURAL","MAT","ND","MatISDisassemblel2gType","MAT_IS_DISASSEMBLE_L2G_",NULL};
    417:   MatISDisassemblel2gType mode = MAT_IS_DISASSEMBLE_L2G_NATURAL;
    418:   MatPartitioning         part;
    419:   PetscSF                 sf;
    420:   PetscErrorCode          ierr;

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

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

    462:     MPIU_Allreduce(&cnt,&i,1,MPIU_INT,MPI_MAX,comm);
    463:     if (i) { /* we detected isolated separator nodes */
    464:       Mat                    A2,A3;
    465:       IS                     *workis,is2;
    466:       PetscScalar            *vals;
    467:       PetscInt               gcnt = i,*dnz,*onz,j,*lndmapi;
    468:       ISLocalToGlobalMapping ll2g;
    469:       PetscBool              flg;
    470:       const PetscInt         *ii,*jj;

    472:       /* communicate global id of separators */
    473:       MatPreallocateInitialize(comm,A->rmap->n,A->cmap->n,dnz,onz);
    474:       for (i = 0, cnt = 0; i < A->rmap->n; i++)
    475:         dnz[i] = ndmapi[i] < 0 ? i + A->rmap->rstart : -1;

    477:       PetscMalloc1(nl,&lndmapi);
    478:       PetscSFBcastBegin(sf,MPIU_INT,dnz,lndmapi,MPI_REPLACE);

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

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

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

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

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

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

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

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

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

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

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

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

    576:       MatDestroy(&A3);
    577:       PetscFree(lndmapi);
    578:       MatPreallocateFinalize(dnz,onz);
    579:       for (i = 0; i < gcnt; i++) {
    580:         ISDestroy(&workis[i]);
    581:       }
    582:       PetscFree(workis);
    583:     }
    584:     ISRestoreIndices(ndmap,&ndmapi);
    585:     PetscSFDestroy(&sf);
    586:     PetscFree(ndmapc);
    587:     ISDestroy(&ndmap);
    588:     ISDestroy(&ndsub);
    589:     ISLocalToGlobalMappingSetBlockSize(*l2g,bs);
    590:     ISLocalToGlobalMappingViewFromOptions(*l2g,NULL,"-matis_nd_l2g_view");
    591:     break;
    592:   case MAT_IS_DISASSEMBLE_L2G_NATURAL:
    593:     if (ismpiaij) {
    594:       MatMPIAIJGetSeqAIJ(A,&Ad,&Ao,&garray);
    595:     } else if (ismpibaij) {
    596:       MatMPIBAIJGetSeqBAIJ(A,&Ad,&Ao,&garray);
    597:     } else SETERRQ1(comm,PETSC_ERR_SUP,"Type %s",((PetscObject)A)->type_name);
    598:     if (!garray) SETERRQ(comm,PETSC_ERR_ARG_WRONGSTATE,"garray not present");
    599:     if (A->rmap->n) {
    600:       PetscInt dc,oc,stc,*aux;

    602:       MatGetLocalSize(A,NULL,&dc);
    603:       MatGetLocalSize(Ao,NULL,&oc);
    604:       MatGetOwnershipRangeColumn(A,&stc,NULL);
    605:       PetscMalloc1((dc+oc)/bs,&aux);
    606:       for (i=0; i<dc/bs; i++) aux[i]       = i+stc/bs;
    607:       for (i=0; i<oc/bs; i++) aux[i+dc/bs] = garray[i];
    608:       ISCreateBlock(comm,bs,(dc+oc)/bs,aux,PETSC_OWN_POINTER,&is);
    609:     } else {
    610:       ISCreateBlock(comm,1,0,NULL,PETSC_OWN_POINTER,&is);
    611:     }
    612:     ISLocalToGlobalMappingCreateIS(is,l2g);
    613:     ISDestroy(&is);
    614:     break;
    615:   default:
    616:     SETERRQ1(comm,PETSC_ERR_ARG_WRONG,"Unsupported l2g disassembling type %D",mode);
    617:   }
    618:   return(0);
    619: }

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

    640:   PetscObjectGetComm((PetscObject)A,&comm);
    641:   MPI_Comm_size(comm,&size);
    642:   if (size == 1) {
    643:     MatConvert_SeqXAIJ_IS(A,type,reuse,newmat);
    644:     return(0);
    645:   }
    646:   if (reuse != MAT_REUSE_MATRIX && A->cmap->N == A->rmap->N) {
    647:     MatMPIXAIJComputeLocalToGlobalMapping_Private(A,&rl2g);
    648:     MatCreate(comm,&B);
    649:     MatSetType(B,MATIS);
    650:     MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
    651:     MatSetLocalToGlobalMapping(B,rl2g,rl2g);
    652:     MatGetBlockSize(A,&bs);
    653:     MatSetBlockSize(B,bs);
    654:     ISLocalToGlobalMappingDestroy(&rl2g);
    655:     if (reuse == MAT_INPLACE_MATRIX) was_inplace = PETSC_TRUE;
    656:     reuse = MAT_REUSE_MATRIX;
    657:   }
    658:   if (reuse == MAT_REUSE_MATRIX) {
    659:     Mat            *newlA, lA;
    660:     IS             rows, cols;
    661:     const PetscInt *ridx, *cidx;
    662:     PetscInt       rbs, cbs, nr, nc;

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

    716:   /* access relevant information from MPIAIJ */
    717:   MatGetOwnershipRange(A,&str,NULL);
    718:   MatGetOwnershipRangeColumn(A,&stc,NULL);
    719:   MatGetLocalSize(A,&dr,&dc);
    720:   MatGetLocalSize(Ao,NULL,&oc);
    721:   MatGetRowIJ(Ad,0,PETSC_FALSE,PETSC_FALSE,&i,&di,&dj,&flg);
    722:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get IJ structure");
    723:   MatGetRowIJ(Ao,0,PETSC_FALSE,PETSC_FALSE,&i,&oi,&oj,&flg);
    724:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get IJ structure");
    725:   nnz = di[dr] + oi[dr];
    726:   /* store original pointers to be restored later */
    727:   odi = di; odj = dj; ooi = oi; ooj = oj;

    729:   /* generate l2g maps for rows and cols */
    730:   ISCreateStride(comm,dr/bs,str/bs,1,&is);
    731:   if (bs > 1) {
    732:     IS is2;

    734:     ISGetLocalSize(is,&i);
    735:     ISGetIndices(is,(const PetscInt**)&aux);
    736:     ISCreateBlock(comm,bs,i,aux,PETSC_COPY_VALUES,&is2);
    737:     ISRestoreIndices(is,(const PetscInt**)&aux);
    738:     ISDestroy(&is);
    739:     is   = is2;
    740:   }
    741:   ISLocalToGlobalMappingCreateIS(is,&rl2g);
    742:   ISDestroy(&is);
    743:   if (dr) {
    744:     PetscMalloc1((dc+oc)/bs,&aux);
    745:     for (i=0; i<dc/bs; i++) aux[i]       = i+stc/bs;
    746:     for (i=0; i<oc/bs; i++) aux[i+dc/bs] = garray[i];
    747:     ISCreateBlock(comm,bs,(dc+oc)/bs,aux,PETSC_OWN_POINTER,&is);
    748:     lc   = dc+oc;
    749:   } else {
    750:     ISCreateBlock(comm,bs,0,NULL,PETSC_OWN_POINTER,&is);
    751:     lc   = 0;
    752:   }
    753:   ISLocalToGlobalMappingCreateIS(is,&cl2g);
    754:   ISDestroy(&is);

    756:   /* create MATIS object */
    757:   MatCreate(comm,&B);
    758:   MatSetSizes(B,dr,dc,PETSC_DECIDE,PETSC_DECIDE);
    759:   MatSetType(B,MATIS);
    760:   MatSetBlockSize(B,bs);
    761:   MatSetLocalToGlobalMapping(B,rl2g,cl2g);
    762:   ISLocalToGlobalMappingDestroy(&rl2g);
    763:   ISLocalToGlobalMappingDestroy(&cl2g);

    765:   /* merge local matrices */
    766:   PetscMalloc1(nnz+dr+1,&aux);
    767:   PetscMalloc1(nnz,&data);
    768:   ii   = aux;
    769:   jj   = aux+dr+1;
    770:   aa   = data;
    771:   *ii  = *(di++) + *(oi++);
    772:   for (jd=0,jo=0,cum=0;*ii<nnz;cum++)
    773:   {
    774:      for (;jd<*di;jd++) { *jj++ = *dj++;      *aa++ = *dd++; }
    775:      for (;jo<*oi;jo++) { *jj++ = *oj++ + dc; *aa++ = *od++; }
    776:      *(++ii) = *(di++) + *(oi++);
    777:   }
    778:   for (;cum<dr;cum++) *(++ii) = nnz;

    780:   MatRestoreRowIJ(Ad,0,PETSC_FALSE,PETSC_FALSE,&i,&odi,&odj,&flg);
    781:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot restore IJ structure");
    782:   MatRestoreRowIJ(Ao,0,PETSC_FALSE,PETSC_FALSE,&i,&ooi,&ooj,&flg);
    783:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot restore IJ structure");
    784:   MatSeqAIJRestoreArray(Ad,&dd);
    785:   MatSeqAIJRestoreArray(Ao,&od);

    787:   ii   = aux;
    788:   jj   = aux+dr+1;
    789:   aa   = data;
    790:   MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,dr,lc,ii,jj,aa,&lA);

    792:   /* create containers to destroy the data */
    793:   ptrs[0] = aux;
    794:   ptrs[1] = data;
    795:   for (i=0; i<2; i++) {
    796:     PetscContainer c;

    798:     PetscContainerCreate(PETSC_COMM_SELF,&c);
    799:     PetscContainerSetPointer(c,ptrs[i]);
    800:     PetscContainerSetUserDestroy(c,PetscContainerUserDestroyDefault);
    801:     PetscObjectCompose((PetscObject)lA,names[i],(PetscObject)c);
    802:     PetscContainerDestroy(&c);
    803:   }
    804:   if (ismpibaij) { /* destroy converted local matrices */
    805:     MatDestroy(&Ad);
    806:     MatDestroy(&Ao);
    807:   }

    809:   /* finalize matrix */
    810:   MatISSetLocalMat(B,lA);
    811:   MatDestroy(&lA);
    812:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    813:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
    814:   if (reuse == MAT_INPLACE_MATRIX) {
    815:     MatHeaderReplace(A,&B);
    816:   } else *newmat = B;
    817:   return(0);
    818: }

    820: PETSC_INTERN PetscErrorCode MatConvert_Nest_IS(Mat A,MatType type,MatReuse reuse,Mat *newmat)
    821: {
    822:   Mat                    **nest,*snest,**rnest,lA,B;
    823:   IS                     *iscol,*isrow,*islrow,*islcol;
    824:   ISLocalToGlobalMapping rl2g,cl2g;
    825:   MPI_Comm               comm;
    826:   PetscInt               *lr,*lc,*l2gidxs;
    827:   PetscInt               i,j,nr,nc,rbs,cbs;
    828:   PetscBool              convert,lreuse,*istrans;
    829:   PetscErrorCode         ierr;

    832:   MatNestGetSubMats(A,&nr,&nc,&nest);
    833:   lreuse = PETSC_FALSE;
    834:   rnest  = NULL;
    835:   if (reuse == MAT_REUSE_MATRIX) {
    836:     PetscBool ismatis,isnest;

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

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

    860:       /* Nested matrices should be of type MATIS */
    861:       PetscObjectTypeCompare((PetscObject)nest[i][j],MATTRANSPOSEMAT,&istrans[ij]);
    862:       if (istrans[ij]) {
    863:         Mat T,lT;
    864:         MatTransposeGetMat(nest[i][j],&T);
    865:         PetscObjectTypeCompare((PetscObject)T,MATIS,&ismatis);
    866:         if (!ismatis) SETERRQ2(comm,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) (transposed) is not of type MATIS",i,j);
    867:         MatISGetLocalMat(T,&lT);
    868:         MatCreateTranspose(lT,&snest[ij]);
    869:       } else {
    870:         PetscObjectTypeCompare((PetscObject)nest[i][j],MATIS,&ismatis);
    871:         if (!ismatis) SETERRQ2(comm,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) is not of type MATIS",i,j);
    872:         MatISGetLocalMat(nest[i][j],&snest[ij]);
    873:       }

    875:       /* Check compatibility of local sizes */
    876:       MatGetSize(snest[ij],&l1,&l2);
    877:       MatGetBlockSizes(snest[ij],&lb1,&lb2);
    878:       if (!l1 || !l2) continue;
    879:       if (lr[i] && l1 != lr[i]) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid local size %D != %D",i,j,lr[i],l1);
    880:       if (lc[j] && l2 != lc[j]) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid local size %D != %D",i,j,lc[j],l2);
    881:       lr[i] = l1;
    882:       lc[j] = l2;

    884:       /* check compatibilty for local matrix reusage */
    885:       if (rnest && !rnest[i][j] != !snest[ij]) lreuse = PETSC_FALSE;
    886:     }
    887:   }

    889:   if (PetscDefined (USE_DEBUG)) {
    890:     /* Check compatibility of l2g maps for rows */
    891:     for (i=0;i<nr;i++) {
    892:       rl2g = NULL;
    893:       for (j=0;j<nc;j++) {
    894:         PetscInt n1,n2;

    896:         if (!nest[i][j]) continue;
    897:         if (istrans[i*nc+j]) {
    898:           Mat T;

    900:           MatTransposeGetMat(nest[i][j],&T);
    901:           MatGetLocalToGlobalMapping(T,NULL,&cl2g);
    902:         } else {
    903:           MatGetLocalToGlobalMapping(nest[i][j],&cl2g,NULL);
    904:         }
    905:         ISLocalToGlobalMappingGetSize(cl2g,&n1);
    906:         if (!n1) continue;
    907:         if (!rl2g) {
    908:           rl2g = cl2g;
    909:         } else {
    910:           const PetscInt *idxs1,*idxs2;
    911:           PetscBool      same;

    913:           ISLocalToGlobalMappingGetSize(rl2g,&n2);
    914:           if (n1 != n2) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid row l2gmap size %D != %D",i,j,n1,n2);
    915:           ISLocalToGlobalMappingGetIndices(cl2g,&idxs1);
    916:           ISLocalToGlobalMappingGetIndices(rl2g,&idxs2);
    917:           PetscArraycmp(idxs1,idxs2,n1,&same);
    918:           ISLocalToGlobalMappingRestoreIndices(cl2g,&idxs1);
    919:           ISLocalToGlobalMappingRestoreIndices(rl2g,&idxs2);
    920:           if (!same) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid row l2gmap",i,j);
    921:         }
    922:       }
    923:     }
    924:     /* Check compatibility of l2g maps for columns */
    925:     for (i=0;i<nc;i++) {
    926:       rl2g = NULL;
    927:       for (j=0;j<nr;j++) {
    928:         PetscInt n1,n2;

    930:         if (!nest[j][i]) continue;
    931:         if (istrans[j*nc+i]) {
    932:           Mat T;

    934:           MatTransposeGetMat(nest[j][i],&T);
    935:           MatGetLocalToGlobalMapping(T,&cl2g,NULL);
    936:         } else {
    937:           MatGetLocalToGlobalMapping(nest[j][i],NULL,&cl2g);
    938:         }
    939:         ISLocalToGlobalMappingGetSize(cl2g,&n1);
    940:         if (!n1) continue;
    941:         if (!rl2g) {
    942:           rl2g = cl2g;
    943:         } else {
    944:           const PetscInt *idxs1,*idxs2;
    945:           PetscBool      same;

    947:           ISLocalToGlobalMappingGetSize(rl2g,&n2);
    948:           if (n1 != n2) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid column l2gmap size %D != %D",j,i,n1,n2);
    949:           ISLocalToGlobalMappingGetIndices(cl2g,&idxs1);
    950:           ISLocalToGlobalMappingGetIndices(rl2g,&idxs2);
    951:           PetscArraycmp(idxs1,idxs2,n1,&same);
    952:           ISLocalToGlobalMappingRestoreIndices(cl2g,&idxs1);
    953:           ISLocalToGlobalMappingRestoreIndices(rl2g,&idxs2);
    954:           if (!same) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid column l2gmap",j,i);
    955:         }
    956:       }
    957:     }
    958:   }

    960:   B = NULL;
    961:   if (reuse != MAT_REUSE_MATRIX) {
    962:     PetscInt stl;

    964:     /* Create l2g map for the rows of the new matrix and index sets for the local MATNEST */
    965:     for (i=0,stl=0;i<nr;i++) stl += lr[i];
    966:     PetscMalloc1(stl,&l2gidxs);
    967:     for (i=0,stl=0;i<nr;i++) {
    968:       Mat            usedmat;
    969:       Mat_IS         *matis;
    970:       const PetscInt *idxs;

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

    975:       /* l2gmap */
    976:       j = 0;
    977:       usedmat = nest[i][j];
    978:       while (!usedmat && j < nc-1) usedmat = nest[i][++j];
    979:       if (!usedmat) SETERRQ(comm,PETSC_ERR_SUP,"Cannot find valid row mat");

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

    1000:     /* Create l2g map for columns of the new matrix and index sets for the local MATNEST */
    1001:     for (i=0,stl=0;i<nc;i++) stl += lc[i];
    1002:     PetscMalloc1(stl,&l2gidxs);
    1003:     for (i=0,stl=0;i<nc;i++) {
    1004:       Mat            usedmat;
    1005:       Mat_IS         *matis;
    1006:       const PetscInt *idxs;

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

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

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

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

    1107:   /* Create local matrix in MATNEST format */
    1108:   convert = PETSC_FALSE;
    1109:   PetscOptionsGetBool(NULL,((PetscObject)A)->prefix,"-matis_convert_local_nest",&convert,NULL);
    1110:   if (convert) {
    1111:     Mat              M;
    1112:     MatISLocalFields lf;
    1113:     PetscContainer   c;

    1115:     MatISGetLocalMat(*newmat,&lA);
    1116:     MatConvert(lA,MATAIJ,MAT_INITIAL_MATRIX,&M);
    1117:     MatISSetLocalMat(*newmat,M);
    1118:     MatDestroy(&M);

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

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

    1133:       ISGetLocalSize(islcol[i],&n);
    1134:       ISStrideGetInfo(islcol[i],&st,NULL);
    1135:       ISCreateStride(comm,n,st,1,&lf->cf[i]);
    1136:     }
    1137:     lf->nr = nr;
    1138:     lf->nc = nc;
    1139:     PetscContainerCreate(PetscObjectComm((PetscObject)(*newmat)),&c);
    1140:     PetscContainerSetPointer(c,lf);
    1141:     PetscContainerSetUserDestroy(c,MatISContainerDestroyFields_Private);
    1142:     PetscObjectCompose((PetscObject)(*newmat),"_convert_nest_lfields",(PetscObject)c);
    1143:     PetscContainerDestroy(&c);
    1144:   }

    1146:   /* Free workspace */
    1147:   for (i=0;i<nr;i++) {
    1148:     ISDestroy(&islrow[i]);
    1149:   }
    1150:   for (i=0;i<nc;i++) {
    1151:     ISDestroy(&islcol[i]);
    1152:   }
    1153:   PetscFree6(isrow,iscol,islrow,islcol,snest,istrans);
    1154:   PetscFree2(lr,lc);
    1155:   return(0);
    1156: }

    1158: static PetscErrorCode MatDiagonalScale_IS(Mat A, Vec l, Vec r)
    1159: {
    1160:   Mat_IS            *matis = (Mat_IS*)A->data;
    1161:   Vec               ll,rr;
    1162:   const PetscScalar *Y,*X;
    1163:   PetscScalar       *x,*y;
    1164:   PetscErrorCode    ierr;

    1167:   if (l) {
    1168:     ll   = matis->y;
    1169:     VecGetArrayRead(l,&Y);
    1170:     VecGetArray(ll,&y);
    1171:     PetscSFBcastBegin(matis->sf,MPIU_SCALAR,Y,y,MPI_REPLACE);
    1172:   } else {
    1173:     ll = NULL;
    1174:   }
    1175:   if (r) {
    1176:     rr   = matis->x;
    1177:     VecGetArrayRead(r,&X);
    1178:     VecGetArray(rr,&x);
    1179:     PetscSFBcastBegin(matis->csf,MPIU_SCALAR,X,x,MPI_REPLACE);
    1180:   } else {
    1181:     rr = NULL;
    1182:   }
    1183:   if (ll) {
    1184:     PetscSFBcastEnd(matis->sf,MPIU_SCALAR,Y,y,MPI_REPLACE);
    1185:     VecRestoreArrayRead(l,&Y);
    1186:     VecRestoreArray(ll,&y);
    1187:   }
    1188:   if (rr) {
    1189:     PetscSFBcastEnd(matis->csf,MPIU_SCALAR,X,x,MPI_REPLACE);
    1190:     VecRestoreArrayRead(r,&X);
    1191:     VecRestoreArray(rr,&x);
    1192:   }
    1193:   MatDiagonalScale(matis->A,ll,rr);
    1194:   return(0);
    1195: }

    1197: static PetscErrorCode MatGetInfo_IS(Mat A,MatInfoType flag,MatInfo *ginfo)
    1198: {
    1199:   Mat_IS         *matis = (Mat_IS*)A->data;
    1200:   MatInfo        info;
    1201:   PetscLogDouble isend[6],irecv[6];
    1202:   PetscInt       bs;

    1206:   MatGetBlockSize(A,&bs);
    1207:   if (matis->A->ops->getinfo) {
    1208:     MatGetInfo(matis->A,MAT_LOCAL,&info);
    1209:     isend[0] = info.nz_used;
    1210:     isend[1] = info.nz_allocated;
    1211:     isend[2] = info.nz_unneeded;
    1212:     isend[3] = info.memory;
    1213:     isend[4] = info.mallocs;
    1214:   } else {
    1215:     isend[0] = 0.;
    1216:     isend[1] = 0.;
    1217:     isend[2] = 0.;
    1218:     isend[3] = 0.;
    1219:     isend[4] = 0.;
    1220:   }
    1221:   isend[5] = matis->A->num_ass;
    1222:   if (flag == MAT_LOCAL) {
    1223:     ginfo->nz_used      = isend[0];
    1224:     ginfo->nz_allocated = isend[1];
    1225:     ginfo->nz_unneeded  = isend[2];
    1226:     ginfo->memory       = isend[3];
    1227:     ginfo->mallocs      = isend[4];
    1228:     ginfo->assemblies   = isend[5];
    1229:   } else if (flag == MAT_GLOBAL_MAX) {
    1230:     MPIU_Allreduce(isend,irecv,6,MPIU_PETSCLOGDOUBLE,MPI_MAX,PetscObjectComm((PetscObject)A));

    1232:     ginfo->nz_used      = irecv[0];
    1233:     ginfo->nz_allocated = irecv[1];
    1234:     ginfo->nz_unneeded  = irecv[2];
    1235:     ginfo->memory       = irecv[3];
    1236:     ginfo->mallocs      = irecv[4];
    1237:     ginfo->assemblies   = irecv[5];
    1238:   } else if (flag == MAT_GLOBAL_SUM) {
    1239:     MPIU_Allreduce(isend,irecv,5,MPIU_PETSCLOGDOUBLE,MPI_SUM,PetscObjectComm((PetscObject)A));

    1241:     ginfo->nz_used      = irecv[0];
    1242:     ginfo->nz_allocated = irecv[1];
    1243:     ginfo->nz_unneeded  = irecv[2];
    1244:     ginfo->memory       = irecv[3];
    1245:     ginfo->mallocs      = irecv[4];
    1246:     ginfo->assemblies   = A->num_ass;
    1247:   }
    1248:   ginfo->block_size        = bs;
    1249:   ginfo->fill_ratio_given  = 0;
    1250:   ginfo->fill_ratio_needed = 0;
    1251:   ginfo->factor_mallocs    = 0;
    1252:   return(0);
    1253: }

    1255: PetscErrorCode MatTranspose_IS(Mat A,MatReuse reuse,Mat *B)
    1256: {
    1257:   Mat                    C,lC,lA;
    1258:   PetscErrorCode         ierr;

    1261:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
    1262:     ISLocalToGlobalMapping rl2g,cl2g;
    1263:     MatCreate(PetscObjectComm((PetscObject)A),&C);
    1264:     MatSetSizes(C,A->cmap->n,A->rmap->n,A->cmap->N,A->rmap->N);
    1265:     MatSetBlockSizes(C,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));
    1266:     MatSetType(C,MATIS);
    1267:     MatGetLocalToGlobalMapping(A,&rl2g,&cl2g);
    1268:     MatSetLocalToGlobalMapping(C,cl2g,rl2g);
    1269:   } else {
    1270:     C = *B;
    1271:   }

    1273:   /* perform local transposition */
    1274:   MatISGetLocalMat(A,&lA);
    1275:   MatTranspose(lA,MAT_INITIAL_MATRIX,&lC);
    1276:   MatISSetLocalMat(C,lC);
    1277:   MatDestroy(&lC);

    1279:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
    1280:     *B = C;
    1281:   } else {
    1282:     MatHeaderMerge(A,&C);
    1283:   }
    1284:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
    1285:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
    1286:   return(0);
    1287: }

    1289: PetscErrorCode  MatDiagonalSet_IS(Mat A,Vec D,InsertMode insmode)
    1290: {
    1291:   Mat_IS         *is = (Mat_IS*)A->data;

    1295:   if (D) { /* MatShift_IS pass D = NULL */
    1296:     VecScatterBegin(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD);
    1297:     VecScatterEnd(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD);
    1298:   }
    1299:   VecPointwiseDivide(is->y,is->y,is->counter);
    1300:   MatDiagonalSet(is->A,is->y,insmode);
    1301:   return(0);
    1302: }

    1304: PetscErrorCode  MatShift_IS(Mat A,PetscScalar a)
    1305: {
    1306:   Mat_IS         *is = (Mat_IS*)A->data;

    1310:   VecSet(is->y,a);
    1311:   MatDiagonalSet_IS(A,NULL,ADD_VALUES);
    1312:   return(0);
    1313: }

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

    1321:   if (PetscUnlikely(m > MATIS_MAX_ENTRIES_INSERTION || n > MATIS_MAX_ENTRIES_INSERTION)) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of row/column indices must be <= %D: they are %D %D",MATIS_MAX_ENTRIES_INSERTION,m,n);
    1322:   ISLocalToGlobalMappingApply(A->rmap->mapping,m,rows,rows_l);
    1323:   ISLocalToGlobalMappingApply(A->cmap->mapping,n,cols,cols_l);
    1324:   MatSetValuesLocal_IS(A,m,rows_l,n,cols_l,values,addv);
    1325:   return(0);
    1326: }

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

    1334:   if (PetscUnlikely(m > MATIS_MAX_ENTRIES_INSERTION || n > MATIS_MAX_ENTRIES_INSERTION)) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of row/column block indices must be <= %D: they are %D %D",MATIS_MAX_ENTRIES_INSERTION,m,n);
    1335:   ISLocalToGlobalMappingApplyBlock(A->rmap->mapping,m,rows,rows_l);
    1336:   ISLocalToGlobalMappingApplyBlock(A->cmap->mapping,n,cols,cols_l);
    1337:   MatSetValuesBlockedLocal_IS(A,m,rows_l,n,cols_l,values,addv);
    1338:   return(0);
    1339: }

    1341: static PetscErrorCode MatCreateSubMatrix_IS(Mat mat,IS irow,IS icol,MatReuse scall,Mat *newmat)
    1342: {
    1343:   Mat               locmat,newlocmat;
    1344:   Mat_IS            *newmatis;
    1345:   const PetscInt    *idxs;
    1346:   PetscInt          i,m,n;
    1347:   PetscErrorCode    ierr;

    1350:   if (scall == MAT_REUSE_MATRIX) {
    1351:     PetscBool ismatis;

    1353:     PetscObjectTypeCompare((PetscObject)*newmat,MATIS,&ismatis);
    1354:     if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Not of MATIS type");
    1355:     newmatis = (Mat_IS*)(*newmat)->data;
    1356:     if (!newmatis->getsub_ris) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Misses local row IS");
    1357:     if (!newmatis->getsub_cis) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Misses local col IS");
    1358:   }
    1359:   /* irow and icol may not have duplicate entries */
    1360:   if (PetscDefined(USE_DEBUG)) {
    1361:     Vec               rtest,ltest;
    1362:     const PetscScalar *array;

    1364:     MatCreateVecs(mat,&ltest,&rtest);
    1365:     ISGetLocalSize(irow,&n);
    1366:     ISGetIndices(irow,&idxs);
    1367:     for (i=0;i<n;i++) {
    1368:       VecSetValue(rtest,idxs[i],1.0,ADD_VALUES);
    1369:     }
    1370:     VecAssemblyBegin(rtest);
    1371:     VecAssemblyEnd(rtest);
    1372:     VecGetLocalSize(rtest,&n);
    1373:     VecGetOwnershipRange(rtest,&m,NULL);
    1374:     VecGetArrayRead(rtest,&array);
    1375:     for (i=0;i<n;i++) if (array[i] != 0. && array[i] != 1.) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Index %D counted %D times! Irow may not have duplicate entries",i+m,(PetscInt)PetscRealPart(array[i]));
    1376:     VecRestoreArrayRead(rtest,&array);
    1377:     ISRestoreIndices(irow,&idxs);
    1378:     ISGetLocalSize(icol,&n);
    1379:     ISGetIndices(icol,&idxs);
    1380:     for (i=0;i<n;i++) {
    1381:       VecSetValue(ltest,idxs[i],1.0,ADD_VALUES);
    1382:     }
    1383:     VecAssemblyBegin(ltest);
    1384:     VecAssemblyEnd(ltest);
    1385:     VecGetLocalSize(ltest,&n);
    1386:     VecGetOwnershipRange(ltest,&m,NULL);
    1387:     VecGetArrayRead(ltest,&array);
    1388:     for (i=0;i<n;i++) if (array[i] != 0. && array[i] != 1.) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Index %D counted %D times! Icol may not have duplicate entries",i+m,(PetscInt)PetscRealPart(array[i]));
    1389:     VecRestoreArrayRead(ltest,&array);
    1390:     ISRestoreIndices(icol,&idxs);
    1391:     VecDestroy(&rtest);
    1392:     VecDestroy(&ltest);
    1393:   }
    1394:   if (scall == MAT_INITIAL_MATRIX) {
    1395:     Mat_IS                 *matis = (Mat_IS*)mat->data;
    1396:     ISLocalToGlobalMapping rl2g;
    1397:     IS                     is;
    1398:     PetscInt               *lidxs,*lgidxs,*newgidxs;
    1399:     PetscInt               ll,newloc,irbs,icbs,arbs,acbs,rbs,cbs;
    1400:     PetscBool              cong;
    1401:     MPI_Comm               comm;

    1403:     PetscObjectGetComm((PetscObject)mat,&comm);
    1404:     MatGetBlockSizes(mat,&arbs,&acbs);
    1405:     ISGetBlockSize(irow,&irbs);
    1406:     ISGetBlockSize(icol,&icbs);
    1407:     rbs  = arbs == irbs ? irbs : 1;
    1408:     cbs  = acbs == icbs ? icbs : 1;
    1409:     ISGetLocalSize(irow,&m);
    1410:     ISGetLocalSize(icol,&n);
    1411:     MatCreate(comm,newmat);
    1412:     MatSetType(*newmat,MATIS);
    1413:     MatSetSizes(*newmat,m,n,PETSC_DECIDE,PETSC_DECIDE);
    1414:     MatSetBlockSizes(*newmat,rbs,cbs);
    1415:     /* communicate irow to their owners in the layout */
    1416:     ISGetIndices(irow,&idxs);
    1417:     PetscLayoutMapLocal(mat->rmap,m,idxs,&ll,&lidxs,&lgidxs);
    1418:     ISRestoreIndices(irow,&idxs);
    1419:     PetscArrayzero(matis->sf_rootdata,matis->sf->nroots);
    1420:     for (i=0;i<ll;i++) matis->sf_rootdata[lidxs[i]] = lgidxs[i]+1;
    1421:     PetscFree(lidxs);
    1422:     PetscFree(lgidxs);
    1423:     PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE);
    1424:     PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE);
    1425:     for (i=0,newloc=0;i<matis->sf->nleaves;i++) if (matis->sf_leafdata[i]) newloc++;
    1426:     PetscMalloc1(newloc,&newgidxs);
    1427:     PetscMalloc1(newloc,&lidxs);
    1428:     for (i=0,newloc=0;i<matis->sf->nleaves;i++)
    1429:       if (matis->sf_leafdata[i]) {
    1430:         lidxs[newloc] = i;
    1431:         newgidxs[newloc++] = matis->sf_leafdata[i]-1;
    1432:       }
    1433:     ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is);
    1434:     ISLocalToGlobalMappingCreateIS(is,&rl2g);
    1435:     ISLocalToGlobalMappingSetBlockSize(rl2g,rbs);
    1436:     ISDestroy(&is);
    1437:     /* local is to extract local submatrix */
    1438:     newmatis = (Mat_IS*)(*newmat)->data;
    1439:     ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_ris);
    1440:     MatHasCongruentLayouts(mat,&cong);
    1441:     if (cong && irow == icol && matis->csf == matis->sf) {
    1442:       MatSetLocalToGlobalMapping(*newmat,rl2g,rl2g);
    1443:       PetscObjectReference((PetscObject)newmatis->getsub_ris);
    1444:       newmatis->getsub_cis = newmatis->getsub_ris;
    1445:     } else {
    1446:       ISLocalToGlobalMapping cl2g;

    1448:       /* communicate icol to their owners in the layout */
    1449:       ISGetIndices(icol,&idxs);
    1450:       PetscLayoutMapLocal(mat->cmap,n,idxs,&ll,&lidxs,&lgidxs);
    1451:       ISRestoreIndices(icol,&idxs);
    1452:       PetscArrayzero(matis->csf_rootdata,matis->csf->nroots);
    1453:       for (i=0;i<ll;i++) matis->csf_rootdata[lidxs[i]] = lgidxs[i]+1;
    1454:       PetscFree(lidxs);
    1455:       PetscFree(lgidxs);
    1456:       PetscSFBcastBegin(matis->csf,MPIU_INT,matis->csf_rootdata,matis->csf_leafdata,MPI_REPLACE);
    1457:       PetscSFBcastEnd(matis->csf,MPIU_INT,matis->csf_rootdata,matis->csf_leafdata,MPI_REPLACE);
    1458:       for (i=0,newloc=0;i<matis->csf->nleaves;i++) if (matis->csf_leafdata[i]) newloc++;
    1459:       PetscMalloc1(newloc,&newgidxs);
    1460:       PetscMalloc1(newloc,&lidxs);
    1461:       for (i=0,newloc=0;i<matis->csf->nleaves;i++)
    1462:         if (matis->csf_leafdata[i]) {
    1463:           lidxs[newloc] = i;
    1464:           newgidxs[newloc++] = matis->csf_leafdata[i]-1;
    1465:         }
    1466:       ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is);
    1467:       ISLocalToGlobalMappingCreateIS(is,&cl2g);
    1468:       ISLocalToGlobalMappingSetBlockSize(cl2g,cbs);
    1469:       ISDestroy(&is);
    1470:       /* local is to extract local submatrix */
    1471:       ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_cis);
    1472:       MatSetLocalToGlobalMapping(*newmat,rl2g,cl2g);
    1473:       ISLocalToGlobalMappingDestroy(&cl2g);
    1474:     }
    1475:     ISLocalToGlobalMappingDestroy(&rl2g);
    1476:   } else {
    1477:     MatISGetLocalMat(*newmat,&newlocmat);
    1478:   }
    1479:   MatISGetLocalMat(mat,&locmat);
    1480:   newmatis = (Mat_IS*)(*newmat)->data;
    1481:   MatCreateSubMatrix(locmat,newmatis->getsub_ris,newmatis->getsub_cis,scall,&newlocmat);
    1482:   if (scall == MAT_INITIAL_MATRIX) {
    1483:     MatISSetLocalMat(*newmat,newlocmat);
    1484:     MatDestroy(&newlocmat);
    1485:   }
    1486:   MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
    1487:   MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
    1488:   return(0);
    1489: }

    1491: static PetscErrorCode MatCopy_IS(Mat A,Mat B,MatStructure str)
    1492: {
    1493:   Mat_IS         *a = (Mat_IS*)A->data,*b;
    1494:   PetscBool      ismatis;

    1498:   PetscObjectTypeCompare((PetscObject)B,MATIS,&ismatis);
    1499:   if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Need to be implemented");
    1500:   b = (Mat_IS*)B->data;
    1501:   MatCopy(a->A,b->A,str);
    1502:   PetscObjectStateIncrease((PetscObject)B);
    1503:   return(0);
    1504: }

    1506: static PetscErrorCode MatMissingDiagonal_IS(Mat A,PetscBool  *missing,PetscInt *d)
    1507: {
    1508:   Vec               v;
    1509:   const PetscScalar *array;
    1510:   PetscInt          i,n;
    1511:   PetscErrorCode    ierr;

    1514:   *missing = PETSC_FALSE;
    1515:   MatCreateVecs(A,NULL,&v);
    1516:   MatGetDiagonal(A,v);
    1517:   VecGetLocalSize(v,&n);
    1518:   VecGetArrayRead(v,&array);
    1519:   for (i=0;i<n;i++) if (array[i] == 0.) break;
    1520:   VecRestoreArrayRead(v,&array);
    1521:   VecDestroy(&v);
    1522:   if (i != n) *missing = PETSC_TRUE;
    1523:   if (d) {
    1524:     *d = -1;
    1525:     if (*missing) {
    1526:       PetscInt rstart;
    1527:       MatGetOwnershipRange(A,&rstart,NULL);
    1528:       *d = i+rstart;
    1529:     }
    1530:   }
    1531:   return(0);
    1532: }

    1534: static PetscErrorCode MatISSetUpSF_IS(Mat B)
    1535: {
    1536:   Mat_IS         *matis = (Mat_IS*)(B->data);
    1537:   const PetscInt *gidxs;
    1538:   PetscInt       nleaves;

    1542:   if (matis->sf) return(0);
    1543:   PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->sf);
    1544:   ISLocalToGlobalMappingGetIndices(B->rmap->mapping,&gidxs);
    1545:   ISLocalToGlobalMappingGetSize(B->rmap->mapping,&nleaves);
    1546:   PetscSFSetGraphLayout(matis->sf,B->rmap,nleaves,NULL,PETSC_OWN_POINTER,gidxs);
    1547:   ISLocalToGlobalMappingRestoreIndices(B->rmap->mapping,&gidxs);
    1548:   PetscMalloc2(matis->sf->nroots,&matis->sf_rootdata,matis->sf->nleaves,&matis->sf_leafdata);
    1549:   if (B->rmap->mapping != B->cmap->mapping) { /* setup SF for columns */
    1550:     ISLocalToGlobalMappingGetSize(B->cmap->mapping,&nleaves);
    1551:     PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->csf);
    1552:     ISLocalToGlobalMappingGetIndices(B->cmap->mapping,&gidxs);
    1553:     PetscSFSetGraphLayout(matis->csf,B->cmap,nleaves,NULL,PETSC_OWN_POINTER,gidxs);
    1554:     ISLocalToGlobalMappingRestoreIndices(B->cmap->mapping,&gidxs);
    1555:     PetscMalloc2(matis->csf->nroots,&matis->csf_rootdata,matis->csf->nleaves,&matis->csf_leafdata);
    1556:   } else {
    1557:     matis->csf = matis->sf;
    1558:     matis->csf_leafdata = matis->sf_leafdata;
    1559:     matis->csf_rootdata = matis->sf_rootdata;
    1560:   }
    1561:   return(0);
    1562: }

    1564: /*@
    1565:    MatISStoreL2L - Store local-to-local operators during the Galerkin process of MatPtAP.

    1567:    Collective

    1569:    Input Parameters:
    1570: +  A - the matrix
    1571: -  store - the boolean flag

    1573:    Level: advanced

    1575:    Notes:

    1577: .seealso: MatCreate(), MatCreateIS(), MatISSetPreallocation(), MatPtAP()
    1578: @*/
    1579: PetscErrorCode  MatISStoreL2L(Mat A, PetscBool store)
    1580: {

    1587:   PetscTryMethod(A,"MatISStoreL2L_C",(Mat,PetscBool),(A,store));
    1588:   return(0);
    1589: }

    1591: static PetscErrorCode MatISStoreL2L_IS(Mat A, PetscBool store)
    1592: {
    1593:   Mat_IS         *matis = (Mat_IS*)(A->data);

    1597:   matis->storel2l = store;
    1598:   if (!store) {
    1599:     PetscObjectCompose((PetscObject)(A),"_MatIS_PtAP_l2l",NULL);
    1600:   }
    1601:   return(0);
    1602: }

    1604: /*@
    1605:    MatISFixLocalEmpty - Compress out zero local rows from the local matrices

    1607:    Collective

    1609:    Input Parameters:
    1610: +  A - the matrix
    1611: -  fix - the boolean flag

    1613:    Level: advanced

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

    1617: .seealso: MatCreate(), MatCreateIS(), MatISSetPreallocation(), MatAssemblyEnd(), MAT_FINAL_ASSEMBLY
    1618: @*/
    1619: PetscErrorCode  MatISFixLocalEmpty(Mat A, PetscBool fix)
    1620: {

    1627:   PetscTryMethod(A,"MatISFixLocalEmpty_C",(Mat,PetscBool),(A,fix));
    1628:   return(0);
    1629: }

    1631: static PetscErrorCode MatISFixLocalEmpty_IS(Mat A, PetscBool fix)
    1632: {
    1633:   Mat_IS *matis = (Mat_IS*)(A->data);

    1636:   matis->locempty = fix;
    1637:   return(0);
    1638: }

    1640: /*@
    1641:    MatISSetPreallocation - Preallocates memory for a MATIS parallel matrix.

    1643:    Collective

    1645:    Input Parameters:
    1646: +  B - the matrix
    1647: .  d_nz  - number of nonzeros per row in DIAGONAL portion of local submatrix
    1648:            (same value is used for all local rows)
    1649: .  d_nnz - array containing the number of nonzeros in the various rows of the
    1650:            DIAGONAL portion of the local submatrix (possibly different for each row)
    1651:            or NULL, if d_nz is used to specify the nonzero structure.
    1652:            The size of this array is equal to the number of local rows, i.e 'm'.
    1653:            For matrices that will be factored, you must leave room for (and set)
    1654:            the diagonal entry even if it is zero.
    1655: .  o_nz  - number of nonzeros per row in the OFF-DIAGONAL portion of local
    1656:            submatrix (same value is used for all local rows).
    1657: -  o_nnz - array containing the number of nonzeros in the various rows of the
    1658:            OFF-DIAGONAL portion of the local submatrix (possibly different for
    1659:            each row) or NULL, if o_nz is used to specify the nonzero
    1660:            structure. The size of this array is equal to the number
    1661:            of local rows, i.e 'm'.

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

    1665:    Level: intermediate

    1667:    Notes:
    1668:     This function has the same interface as the MPIAIJ preallocation routine in order to simplify the transition
    1669:           from the asssembled format to the unassembled one. It overestimates the preallocation of MATIS local
    1670:           matrices; for exact preallocation, the user should set the preallocation directly on local matrix objects.

    1672: .seealso: MatCreate(), MatCreateIS(), MatMPIAIJSetPreallocation(), MatISGetLocalMat(), MATIS
    1673: @*/
    1674: PetscErrorCode  MatISSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
    1675: {

    1681:   PetscTryMethod(B,"MatISSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz));
    1682:   return(0);
    1683: }

    1685: /* this is used by DMDA */
    1686: PETSC_EXTERN PetscErrorCode  MatISSetPreallocation_IS(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
    1687: {
    1688:   Mat_IS         *matis = (Mat_IS*)(B->data);
    1689:   PetscInt       bs,i,nlocalcols;

    1693:   if (!matis->A) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"You should first call MatSetLocalToGlobalMapping");

    1695:   if (!d_nnz) for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] = d_nz;
    1696:   else for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] = d_nnz[i];

    1698:   if (!o_nnz) for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] += o_nz;
    1699:   else for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] += o_nnz[i];

    1701:   PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE);
    1702:   MatGetSize(matis->A,NULL,&nlocalcols);
    1703:   MatGetBlockSize(matis->A,&bs);
    1704:   PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE);

    1706:   for (i=0;i<matis->sf->nleaves;i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i],nlocalcols);
    1707:   MatSeqAIJSetPreallocation(matis->A,0,matis->sf_leafdata);
    1708: #if defined(PETSC_HAVE_HYPRE)
    1709:   MatHYPRESetPreallocation(matis->A,0,matis->sf_leafdata,0,NULL);
    1710: #endif

    1712:   for (i=0;i<matis->sf->nleaves/bs;i++) matis->sf_leafdata[i] = matis->sf_leafdata[i*bs]/bs;
    1713:   MatSeqBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata);

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

    1719:   /* for other matrix types */
    1720:   MatSetUp(matis->A);
    1721:   return(0);
    1722: }

    1724: PETSC_EXTERN PetscErrorCode MatISSetMPIXAIJPreallocation_Private(Mat A, Mat B, PetscBool maxreduce)
    1725: {
    1726:   Mat_IS          *matis = (Mat_IS*)(A->data);
    1727:   PetscInt        *my_dnz,*my_onz,*dnz,*onz,*mat_ranges,*row_ownership;
    1728:   const PetscInt  *global_indices_r,*global_indices_c;
    1729:   PetscInt        i,j,bs,rows,cols;
    1730:   PetscInt        lrows,lcols;
    1731:   PetscInt        local_rows,local_cols;
    1732:   PetscMPIInt     size;
    1733:   PetscBool       isdense,issbaij;
    1734:   PetscErrorCode  ierr;

    1737:   MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
    1738:   MatGetSize(A,&rows,&cols);
    1739:   MatGetBlockSize(A,&bs);
    1740:   MatGetSize(matis->A,&local_rows,&local_cols);
    1741:   PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isdense);
    1742:   PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);
    1743:   ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&global_indices_r);
    1744:   if (A->rmap->mapping != A->cmap->mapping) {
    1745:     ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&global_indices_c);
    1746:   } else {
    1747:     global_indices_c = global_indices_r;
    1748:   }

    1750:   if (issbaij) {
    1751:     MatGetRowUpperTriangular(matis->A);
    1752:   }
    1753:   /*
    1754:      An SF reduce is needed to sum up properly on shared rows.
    1755:      Note that generally preallocation is not exact, since it overestimates nonzeros
    1756:   */
    1757:   MatGetLocalSize(A,&lrows,&lcols);
    1758:   MatPreallocateInitialize(PetscObjectComm((PetscObject)A),lrows,lcols,dnz,onz);
    1759:   /* All processes need to compute entire row ownership */
    1760:   PetscMalloc1(rows,&row_ownership);
    1761:   MatGetOwnershipRanges(A,(const PetscInt**)&mat_ranges);
    1762:   for (i=0;i<size;i++) {
    1763:     for (j=mat_ranges[i];j<mat_ranges[i+1];j++) {
    1764:       row_ownership[j] = i;
    1765:     }
    1766:   }
    1767:   MatGetOwnershipRangesColumn(A,(const PetscInt**)&mat_ranges);

    1769:   /*
    1770:      my_dnz and my_onz contains exact contribution to preallocation from each local mat
    1771:      then, they will be summed up properly. This way, preallocation is always sufficient
    1772:   */
    1773:   PetscCalloc2(local_rows,&my_dnz,local_rows,&my_onz);
    1774:   /* preallocation as a MATAIJ */
    1775:   if (isdense) { /* special case for dense local matrices */
    1776:     for (i=0;i<local_rows;i++) {
    1777:       PetscInt owner = row_ownership[global_indices_r[i]];
    1778:       for (j=0;j<local_cols;j++) {
    1779:         PetscInt index_col = global_indices_c[j];
    1780:         if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1]) { /* diag block */
    1781:           my_dnz[i] += 1;
    1782:         } else { /* offdiag block */
    1783:           my_onz[i] += 1;
    1784:         }
    1785:       }
    1786:     }
    1787:   } else if (matis->A->ops->getrowij) {
    1788:     const PetscInt *ii,*jj,*jptr;
    1789:     PetscBool      done;
    1790:     MatGetRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done);
    1791:     if (!done) SETERRQ(PetscObjectComm((PetscObject)(matis->A)),PETSC_ERR_PLIB,"Error in MatGetRowIJ");
    1792:     jptr = jj;
    1793:     for (i=0;i<local_rows;i++) {
    1794:       PetscInt index_row = global_indices_r[i];
    1795:       for (j=0;j<ii[i+1]-ii[i];j++,jptr++) {
    1796:         PetscInt owner = row_ownership[index_row];
    1797:         PetscInt index_col = global_indices_c[*jptr];
    1798:         if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1]) { /* diag block */
    1799:           my_dnz[i] += 1;
    1800:         } else { /* offdiag block */
    1801:           my_onz[i] += 1;
    1802:         }
    1803:         /* same as before, interchanging rows and cols */
    1804:         if (issbaij && index_col != index_row) {
    1805:           owner = row_ownership[index_col];
    1806:           if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1]) {
    1807:             my_dnz[*jptr] += 1;
    1808:           } else {
    1809:             my_onz[*jptr] += 1;
    1810:           }
    1811:         }
    1812:       }
    1813:     }
    1814:     MatRestoreRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done);
    1815:     if (!done) SETERRQ(PetscObjectComm((PetscObject)(matis->A)),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ");
    1816:   } else { /* loop over rows and use MatGetRow */
    1817:     for (i=0;i<local_rows;i++) {
    1818:       const PetscInt *cols;
    1819:       PetscInt       ncols,index_row = global_indices_r[i];
    1820:       MatGetRow(matis->A,i,&ncols,&cols,NULL);
    1821:       for (j=0;j<ncols;j++) {
    1822:         PetscInt owner = row_ownership[index_row];
    1823:         PetscInt index_col = global_indices_c[cols[j]];
    1824:         if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1]) { /* diag block */
    1825:           my_dnz[i] += 1;
    1826:         } else { /* offdiag block */
    1827:           my_onz[i] += 1;
    1828:         }
    1829:         /* same as before, interchanging rows and cols */
    1830:         if (issbaij && index_col != index_row) {
    1831:           owner = row_ownership[index_col];
    1832:           if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1]) {
    1833:             my_dnz[cols[j]] += 1;
    1834:           } else {
    1835:             my_onz[cols[j]] += 1;
    1836:           }
    1837:         }
    1838:       }
    1839:       MatRestoreRow(matis->A,i,&ncols,&cols,NULL);
    1840:     }
    1841:   }
    1842:   if (global_indices_c != global_indices_r) {
    1843:     ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&global_indices_c);
    1844:   }
    1845:   ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&global_indices_r);
    1846:   PetscFree(row_ownership);

    1848:   /* Reduce my_dnz and my_onz */
    1849:   if (maxreduce) {
    1850:     PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);
    1851:     PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);
    1852:     PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);
    1853:     PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);
    1854:   } else {
    1855:     PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);
    1856:     PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);
    1857:     PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);
    1858:     PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);
    1859:   }
    1860:   PetscFree2(my_dnz,my_onz);

    1862:   /* Resize preallocation if overestimated */
    1863:   for (i=0;i<lrows;i++) {
    1864:     dnz[i] = PetscMin(dnz[i],lcols);
    1865:     onz[i] = PetscMin(onz[i],cols-lcols);
    1866:   }

    1868:   /* Set preallocation */
    1869:   MatSeqAIJSetPreallocation(B,0,dnz);
    1870:   MatMPIAIJSetPreallocation(B,0,dnz,0,onz);
    1871:   for (i=0;i<lrows;i+=bs) {
    1872:     PetscInt b, d = dnz[i],o = onz[i];

    1874:     for (b=1;b<bs;b++) {
    1875:       d = PetscMax(d,dnz[i+b]);
    1876:       o = PetscMax(o,onz[i+b]);
    1877:     }
    1878:     dnz[i/bs] = PetscMin(d/bs + d%bs,lcols/bs);
    1879:     onz[i/bs] = PetscMin(o/bs + o%bs,(cols-lcols)/bs);
    1880:   }
    1881:   MatSeqBAIJSetPreallocation(B,bs,0,dnz);
    1882:   MatMPIBAIJSetPreallocation(B,bs,0,dnz,0,onz);
    1883:   MatMPISBAIJSetPreallocation(B,bs,0,dnz,0,onz);
    1884:   MatPreallocateFinalize(dnz,onz);
    1885:   if (issbaij) {
    1886:     MatRestoreRowUpperTriangular(matis->A);
    1887:   }
    1888:   MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
    1889:   return(0);
    1890: }

    1892: PETSC_INTERN PetscErrorCode MatConvert_IS_XAIJ(Mat mat, MatType mtype, MatReuse reuse, Mat *M)
    1893: {
    1894:   Mat_IS            *matis = (Mat_IS*)(mat->data);
    1895:   Mat               local_mat,MT;
    1896:   PetscInt          rbs,cbs,rows,cols,lrows,lcols;
    1897:   PetscInt          local_rows,local_cols;
    1898:   PetscBool         isseqdense,isseqsbaij,isseqaij,isseqbaij;
    1899:   PetscMPIInt       size;
    1900:   const PetscScalar *array;
    1901:   PetscErrorCode    ierr;

    1904:   MPI_Comm_size(PetscObjectComm((PetscObject)mat),&size);
    1905:   if (size == 1 && mat->rmap->N == matis->A->rmap->N && mat->cmap->N == matis->A->cmap->N) {
    1906:     Mat      B;
    1907:     IS       irows = NULL,icols = NULL;
    1908:     PetscInt rbs,cbs;

    1910:     ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping,&rbs);
    1911:     ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping,&cbs);
    1912:     if (reuse != MAT_REUSE_MATRIX) { /* check if l2g maps are one-to-one */
    1913:       IS             rows,cols;
    1914:       const PetscInt *ridxs,*cidxs;
    1915:       PetscInt       i,nw,*work;

    1917:       ISLocalToGlobalMappingGetBlockIndices(mat->rmap->mapping,&ridxs);
    1918:       ISLocalToGlobalMappingGetSize(mat->rmap->mapping,&nw);
    1919:       nw   = nw/rbs;
    1920:       PetscCalloc1(nw,&work);
    1921:       for (i=0;i<nw;i++) work[ridxs[i]] += 1;
    1922:       for (i=0;i<nw;i++) if (!work[i] || work[i] > 1) break;
    1923:       if (i == nw) {
    1924:         ISCreateBlock(PETSC_COMM_SELF,rbs,nw,ridxs,PETSC_USE_POINTER,&rows);
    1925:         ISSetPermutation(rows);
    1926:         ISInvertPermutation(rows,PETSC_DECIDE,&irows);
    1927:         ISDestroy(&rows);
    1928:       }
    1929:       ISLocalToGlobalMappingRestoreBlockIndices(mat->rmap->mapping,&ridxs);
    1930:       PetscFree(work);
    1931:       if (irows && mat->rmap->mapping != mat->cmap->mapping) {
    1932:         ISLocalToGlobalMappingGetBlockIndices(mat->cmap->mapping,&cidxs);
    1933:         ISLocalToGlobalMappingGetSize(mat->cmap->mapping,&nw);
    1934:         nw   = nw/cbs;
    1935:         PetscCalloc1(nw,&work);
    1936:         for (i=0;i<nw;i++) work[cidxs[i]] += 1;
    1937:         for (i=0;i<nw;i++) if (!work[i] || work[i] > 1) break;
    1938:         if (i == nw) {
    1939:           ISCreateBlock(PETSC_COMM_SELF,cbs,nw,cidxs,PETSC_USE_POINTER,&cols);
    1940:           ISSetPermutation(cols);
    1941:           ISInvertPermutation(cols,PETSC_DECIDE,&icols);
    1942:           ISDestroy(&cols);
    1943:         }
    1944:         ISLocalToGlobalMappingRestoreBlockIndices(mat->cmap->mapping,&cidxs);
    1945:         PetscFree(work);
    1946:       } else if (irows) {
    1947:         PetscObjectReference((PetscObject)irows);
    1948:         icols = irows;
    1949:       }
    1950:     } else {
    1951:       PetscObjectQuery((PetscObject)(*M),"_MatIS_IS_XAIJ_irows",(PetscObject*)&irows);
    1952:       PetscObjectQuery((PetscObject)(*M),"_MatIS_IS_XAIJ_icols",(PetscObject*)&icols);
    1953:       if (irows) { PetscObjectReference((PetscObject)irows); }
    1954:       if (icols) { PetscObjectReference((PetscObject)icols); }
    1955:     }
    1956:     if (!irows || !icols) {
    1957:       ISDestroy(&icols);
    1958:       ISDestroy(&irows);
    1959:       goto general_assembly;
    1960:     }
    1961:     MatConvert(matis->A,mtype,MAT_INITIAL_MATRIX,&B);
    1962:     if (reuse != MAT_INPLACE_MATRIX) {
    1963:       MatCreateSubMatrix(B,irows,icols,reuse,M);
    1964:       PetscObjectCompose((PetscObject)(*M),"_MatIS_IS_XAIJ_irows",(PetscObject)irows);
    1965:       PetscObjectCompose((PetscObject)(*M),"_MatIS_IS_XAIJ_icols",(PetscObject)icols);
    1966:     } else {
    1967:       Mat C;

    1969:       MatCreateSubMatrix(B,irows,icols,MAT_INITIAL_MATRIX,&C);
    1970:       MatHeaderReplace(mat,&C);
    1971:     }
    1972:     MatDestroy(&B);
    1973:     ISDestroy(&icols);
    1974:     ISDestroy(&irows);
    1975:     return(0);
    1976:   }
    1977: general_assembly:
    1978:   MatGetSize(mat,&rows,&cols);
    1979:   ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping,&rbs);
    1980:   ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping,&cbs);
    1981:   MatGetLocalSize(mat,&lrows,&lcols);
    1982:   MatGetSize(matis->A,&local_rows,&local_cols);
    1983:   PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isseqdense);
    1984:   PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQAIJ,&isseqaij);
    1985:   PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQBAIJ,&isseqbaij);
    1986:   PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&isseqsbaij);
    1987:   if (!isseqdense && !isseqaij && !isseqbaij && !isseqsbaij) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)(matis->A))->type_name);
    1988:   if (PetscDefined (USE_DEBUG)) {
    1989:     PetscBool         lb[4],bb[4];

    1991:     lb[0] = isseqdense;
    1992:     lb[1] = isseqaij;
    1993:     lb[2] = isseqbaij;
    1994:     lb[3] = isseqsbaij;
    1995:     MPIU_Allreduce(lb,bb,4,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)mat));
    1996:     if (!bb[0] && !bb[1] && !bb[2] && !bb[3]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Local matrices must have the same type");
    1997:   }

    1999:   if (reuse != MAT_REUSE_MATRIX) {
    2000:     MatCreate(PetscObjectComm((PetscObject)mat),&MT);
    2001:     MatSetSizes(MT,lrows,lcols,rows,cols);
    2002:     MatSetType(MT,mtype);
    2003:     MatSetBlockSizes(MT,rbs,cbs);
    2004:     MatISSetMPIXAIJPreallocation_Private(mat,MT,PETSC_FALSE);
    2005:   } else {
    2006:     PetscInt mrbs,mcbs,mrows,mcols,mlrows,mlcols;

    2008:     /* some checks */
    2009:     MT   = *M;
    2010:     MatGetBlockSizes(MT,&mrbs,&mcbs);
    2011:     MatGetSize(MT,&mrows,&mcols);
    2012:     MatGetLocalSize(MT,&mlrows,&mlcols);
    2013:     if (mrows != rows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of rows (%d != %d)",rows,mrows);
    2014:     if (mcols != cols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of cols (%d != %d)",cols,mcols);
    2015:     if (mlrows != lrows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local rows (%d != %d)",lrows,mlrows);
    2016:     if (mlcols != lcols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local cols (%d != %d)",lcols,mlcols);
    2017:     if (mrbs != rbs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong row block size (%d != %d)",rbs,mrbs);
    2018:     if (mcbs != cbs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong col block size (%d != %d)",cbs,mcbs);
    2019:     MatZeroEntries(MT);
    2020:   }

    2022:   if (isseqsbaij || isseqbaij) {
    2023:     MatConvert(matis->A,MATSEQAIJ,MAT_INITIAL_MATRIX,&local_mat);
    2024:     isseqaij = PETSC_TRUE;
    2025:   } else {
    2026:     PetscObjectReference((PetscObject)matis->A);
    2027:     local_mat = matis->A;
    2028:   }

    2030:   /* Set values */
    2031:   MatSetLocalToGlobalMapping(MT,mat->rmap->mapping,mat->cmap->mapping);
    2032:   if (isseqdense) { /* special case for dense local matrices */
    2033:     PetscInt          i,*dummy;

    2035:     PetscMalloc1(PetscMax(local_rows,local_cols),&dummy);
    2036:     for (i=0;i<PetscMax(local_rows,local_cols);i++) dummy[i] = i;
    2037:     MatSetOption(MT,MAT_ROW_ORIENTED,PETSC_FALSE);
    2038:     MatDenseGetArrayRead(local_mat,&array);
    2039:     MatSetValuesLocal(MT,local_rows,dummy,local_cols,dummy,array,ADD_VALUES);
    2040:     MatDenseRestoreArrayRead(local_mat,&array);
    2041:     PetscFree(dummy);
    2042:   } else if (isseqaij) {
    2043:     const PetscInt *blocks;
    2044:     PetscInt       i,nvtxs,*xadj,*adjncy, nb;
    2045:     PetscBool      done;
    2046:     PetscScalar    *sarray;

    2048:     MatGetRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);
    2049:     if (!done) SETERRQ(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatGetRowIJ");
    2050:     MatSeqAIJGetArray(local_mat,&sarray);
    2051:     MatGetVariableBlockSizes(local_mat,&nb,&blocks);
    2052:     if (nb) { /* speed up assembly for special blocked matrices (used by BDDC) */
    2053:       PetscInt sum;

    2055:       for (i=0,sum=0;i<nb;i++) sum += blocks[i];
    2056:       if (sum == nvtxs) {
    2057:         PetscInt r;

    2059:         for (i=0,r=0;i<nb;i++) {
    2060:           if (PetscUnlikelyDebug(blocks[i] != xadj[r+1] - xadj[r])) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Invalid block sizes prescribed for block %D: expected %D, got %D",i,blocks[i],xadj[r+1] - xadj[r]);
    2061:           MatSetValuesLocal(MT,blocks[i],adjncy+xadj[r],blocks[i],adjncy+xadj[r],sarray+xadj[r],ADD_VALUES);
    2062:           r   += blocks[i];
    2063:         }
    2064:       } else {
    2065:         for (i=0;i<nvtxs;i++) {
    2066:           MatSetValuesLocal(MT,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],sarray+xadj[i],ADD_VALUES);
    2067:         }
    2068:       }
    2069:     } else {
    2070:       for (i=0;i<nvtxs;i++) {
    2071:         MatSetValuesLocal(MT,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],sarray+xadj[i],ADD_VALUES);
    2072:       }
    2073:     }
    2074:     MatRestoreRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);
    2075:     if (!done) SETERRQ(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ");
    2076:     MatSeqAIJRestoreArray(local_mat,&sarray);
    2077:   } else { /* very basic values insertion for all other matrix types */
    2078:     PetscInt i;

    2080:     for (i=0;i<local_rows;i++) {
    2081:       PetscInt       j;
    2082:       const PetscInt *local_indices_cols;

    2084:       MatGetRow(local_mat,i,&j,&local_indices_cols,&array);
    2085:       MatSetValuesLocal(MT,1,&i,j,local_indices_cols,array,ADD_VALUES);
    2086:       MatRestoreRow(local_mat,i,&j,&local_indices_cols,&array);
    2087:     }
    2088:   }
    2089:   MatAssemblyBegin(MT,MAT_FINAL_ASSEMBLY);
    2090:   MatDestroy(&local_mat);
    2091:   MatAssemblyEnd(MT,MAT_FINAL_ASSEMBLY);
    2092:   if (isseqdense) {
    2093:     MatSetOption(MT,MAT_ROW_ORIENTED,PETSC_TRUE);
    2094:   }
    2095:   if (reuse == MAT_INPLACE_MATRIX) {
    2096:     MatHeaderReplace(mat,&MT);
    2097:   } else if (reuse == MAT_INITIAL_MATRIX) {
    2098:     *M = MT;
    2099:   }
    2100:   return(0);
    2101: }

    2103: /*@
    2104:     MatISGetMPIXAIJ - Converts MATIS matrix into a parallel AIJ format

    2106:   Input Parameter:
    2107: +  mat - the matrix (should be of type MATIS)
    2108: -  reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX

    2110:   Output Parameter:
    2111: .  newmat - the matrix in AIJ format

    2113:   Level: developer

    2115:   Notes:
    2116:     This function has been deprecated and it will be removed in future releases. Update your code to use the MatConvert() interface.

    2118: .seealso: MATIS, MatConvert()
    2119: @*/
    2120: PetscErrorCode MatISGetMPIXAIJ(Mat mat, MatReuse reuse, Mat *newmat)
    2121: {

    2128:   if (reuse == MAT_REUSE_MATRIX) {
    2131:     if (mat == *newmat) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse the same matrix");
    2132:   }
    2133:   PetscUseMethod(mat,"MatISGetMPIXAIJ_C",(Mat,MatType,MatReuse,Mat*),(mat,MATAIJ,reuse,newmat));
    2134:   return(0);
    2135: }

    2137: PetscErrorCode MatDuplicate_IS(Mat mat,MatDuplicateOption op,Mat *newmat)
    2138: {
    2140:   Mat_IS         *matis = (Mat_IS*)(mat->data);
    2141:   PetscInt       rbs,cbs,m,n,M,N;
    2142:   Mat            B,localmat;

    2145:   ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping,&rbs);
    2146:   ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping,&cbs);
    2147:   MatGetSize(mat,&M,&N);
    2148:   MatGetLocalSize(mat,&m,&n);
    2149:   MatCreate(PetscObjectComm((PetscObject)mat),&B);
    2150:   MatSetSizes(B,m,n,M,N);
    2151:   MatSetBlockSize(B,rbs == cbs ? rbs : 1);
    2152:   MatSetType(B,MATIS);
    2153:   MatISSetLocalMatType(B,matis->lmattype);
    2154:   MatSetLocalToGlobalMapping(B,mat->rmap->mapping,mat->cmap->mapping);
    2155:   MatDuplicate(matis->A,op,&localmat);
    2156:   MatISSetLocalMat(B,localmat);
    2157:   MatDestroy(&localmat);
    2158:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    2159:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
    2160:   *newmat = B;
    2161:   return(0);
    2162: }

    2164: static PetscErrorCode MatIsHermitian_IS(Mat A,PetscReal tol,PetscBool  *flg)
    2165: {
    2167:   Mat_IS         *matis = (Mat_IS*)A->data;
    2168:   PetscBool      local_sym;

    2171:   MatIsHermitian(matis->A,tol,&local_sym);
    2172:   MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
    2173:   return(0);
    2174: }

    2176: static PetscErrorCode MatIsSymmetric_IS(Mat A,PetscReal tol,PetscBool  *flg)
    2177: {
    2179:   Mat_IS         *matis = (Mat_IS*)A->data;
    2180:   PetscBool      local_sym;

    2183:   MatIsSymmetric(matis->A,tol,&local_sym);
    2184:   MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
    2185:   return(0);
    2186: }

    2188: static PetscErrorCode MatIsStructurallySymmetric_IS(Mat A,PetscBool  *flg)
    2189: {
    2191:   Mat_IS         *matis = (Mat_IS*)A->data;
    2192:   PetscBool      local_sym;

    2195:   if (A->rmap->mapping != A->cmap->mapping) {
    2196:     *flg = PETSC_FALSE;
    2197:     return(0);
    2198:   }
    2199:   MatIsStructurallySymmetric(matis->A,&local_sym);
    2200:   MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
    2201:   return(0);
    2202: }

    2204: static PetscErrorCode MatDestroy_IS(Mat A)
    2205: {
    2207:   Mat_IS         *b = (Mat_IS*)A->data;

    2210:   PetscFree(b->bdiag);
    2211:   PetscFree(b->lmattype);
    2212:   MatDestroy(&b->A);
    2213:   VecScatterDestroy(&b->cctx);
    2214:   VecScatterDestroy(&b->rctx);
    2215:   VecDestroy(&b->x);
    2216:   VecDestroy(&b->y);
    2217:   VecDestroy(&b->counter);
    2218:   ISDestroy(&b->getsub_ris);
    2219:   ISDestroy(&b->getsub_cis);
    2220:   if (b->sf != b->csf) {
    2221:     PetscSFDestroy(&b->csf);
    2222:     PetscFree2(b->csf_rootdata,b->csf_leafdata);
    2223:   } else b->csf = NULL;
    2224:   PetscSFDestroy(&b->sf);
    2225:   PetscFree2(b->sf_rootdata,b->sf_leafdata);
    2226:   PetscFree(A->data);
    2227:   PetscObjectChangeTypeName((PetscObject)A,NULL);
    2228:   PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMatType_C",NULL);
    2229:   PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL);
    2230:   PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",NULL);
    2231:   PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",NULL);
    2232:   PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",NULL);
    2233:   PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",NULL);
    2234:   PetscObjectComposeFunction((PetscObject)A,"MatISFixLocalEmpty_C",NULL);
    2235:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpiaij_C",NULL);
    2236:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpibaij_C",NULL);
    2237:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpisbaij_C",NULL);
    2238:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqaij_C",NULL);
    2239:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqbaij_C",NULL);
    2240:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqsbaij_C",NULL);
    2241:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_aij_C",NULL);
    2242:   return(0);
    2243: }

    2245: static PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y)
    2246: {
    2248:   Mat_IS         *is  = (Mat_IS*)A->data;
    2249:   PetscScalar    zero = 0.0;

    2252:   /*  scatter the global vector x into the local work vector */
    2253:   VecScatterBegin(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
    2254:   VecScatterEnd(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);

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

    2259:   /* scatter product back into global memory */
    2260:   VecSet(y,zero);
    2261:   VecScatterBegin(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
    2262:   VecScatterEnd(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
    2263:   return(0);
    2264: }

    2266: static PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
    2267: {
    2268:   Vec            temp_vec;

    2272:   if (v3 != v2) {
    2273:     MatMult(A,v1,v3);
    2274:     VecAXPY(v3,1.0,v2);
    2275:   } else {
    2276:     VecDuplicate(v2,&temp_vec);
    2277:     MatMult(A,v1,temp_vec);
    2278:     VecAXPY(temp_vec,1.0,v2);
    2279:     VecCopy(temp_vec,v3);
    2280:     VecDestroy(&temp_vec);
    2281:   }
    2282:   return(0);
    2283: }

    2285: static PetscErrorCode MatMultTranspose_IS(Mat A,Vec y,Vec x)
    2286: {
    2287:   Mat_IS         *is = (Mat_IS*)A->data;

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

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

    2298:   /* scatter product back into global vector */
    2299:   VecSet(x,0);
    2300:   VecScatterBegin(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);
    2301:   VecScatterEnd(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);
    2302:   return(0);
    2303: }

    2305: static PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
    2306: {
    2307:   Vec            temp_vec;

    2311:   if (v3 != v2) {
    2312:     MatMultTranspose(A,v1,v3);
    2313:     VecAXPY(v3,1.0,v2);
    2314:   } else {
    2315:     VecDuplicate(v2,&temp_vec);
    2316:     MatMultTranspose(A,v1,temp_vec);
    2317:     VecAXPY(temp_vec,1.0,v2);
    2318:     VecCopy(temp_vec,v3);
    2319:     VecDestroy(&temp_vec);
    2320:   }
    2321:   return(0);
    2322: }

    2324: static PetscErrorCode MatView_IS(Mat A,PetscViewer viewer)
    2325: {
    2326:   Mat_IS         *a = (Mat_IS*)A->data;
    2328:   PetscViewer    sviewer;
    2329:   PetscBool      isascii,view = PETSC_TRUE;

    2332:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
    2333:   if (isascii)  {
    2334:     PetscViewerFormat format;

    2336:     PetscViewerGetFormat(viewer,&format);
    2337:     if (format == PETSC_VIEWER_ASCII_INFO) view = PETSC_FALSE;
    2338:   }
    2339:   if (!view) return(0);
    2340:   PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
    2341:   MatView(a->A,sviewer);
    2342:   PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
    2343:   PetscViewerFlush(viewer);
    2344:   return(0);
    2345: }

    2347: static PetscErrorCode MatInvertBlockDiagonal_IS(Mat mat,const PetscScalar **values)
    2348: {
    2349:   Mat_IS            *is = (Mat_IS*)mat->data;
    2350:   MPI_Datatype      nodeType;
    2351:   const PetscScalar *lv;
    2352:   PetscInt          bs;
    2353:   PetscErrorCode    ierr;

    2356:   MatGetBlockSize(mat,&bs);
    2357:   MatSetBlockSize(is->A,bs);
    2358:   MatInvertBlockDiagonal(is->A,&lv);
    2359:   if (!is->bdiag) {
    2360:     PetscMalloc1(bs*mat->rmap->n,&is->bdiag);
    2361:   }
    2362:   MPI_Type_contiguous(bs,MPIU_SCALAR,&nodeType);
    2363:   MPI_Type_commit(&nodeType);
    2364:   PetscSFReduceBegin(is->sf,nodeType,lv,is->bdiag,MPI_REPLACE);
    2365:   PetscSFReduceEnd(is->sf,nodeType,lv,is->bdiag,MPI_REPLACE);
    2366:   MPI_Type_free(&nodeType);
    2367:   if (values) *values = is->bdiag;
    2368:   return(0);
    2369: }

    2371: static PetscErrorCode MatISSetUpScatters_Private(Mat A)
    2372: {
    2373:   Vec            cglobal,rglobal;
    2374:   IS             from;
    2375:   Mat_IS         *is = (Mat_IS*)A->data;
    2376:   PetscScalar    sum;
    2377:   const PetscInt *garray;
    2378:   PetscInt       nr,rbs,nc,cbs;
    2379:   PetscBool      iscuda;

    2383:   ISLocalToGlobalMappingGetSize(A->rmap->mapping,&nr);
    2384:   ISLocalToGlobalMappingGetBlockSize(A->rmap->mapping,&rbs);
    2385:   ISLocalToGlobalMappingGetSize(A->cmap->mapping,&nc);
    2386:   ISLocalToGlobalMappingGetBlockSize(A->cmap->mapping,&cbs);
    2387:   VecDestroy(&is->x);
    2388:   VecDestroy(&is->y);
    2389:   VecDestroy(&is->counter);
    2390:   VecScatterDestroy(&is->rctx);
    2391:   VecScatterDestroy(&is->cctx);
    2392:   MatCreateVecs(is->A,&is->x,&is->y);
    2393:   VecBindToCPU(is->y,PETSC_TRUE);
    2394:   PetscObjectTypeCompare((PetscObject)is->y,VECSEQCUDA,&iscuda);
    2395:   if (iscuda) {
    2396:     PetscFree(A->defaultvectype);
    2397:     PetscStrallocpy(VECCUDA,&A->defaultvectype);
    2398:   }
    2399:   MatCreateVecs(A,&cglobal,&rglobal);
    2400:   VecBindToCPU(rglobal,PETSC_TRUE);
    2401:   ISLocalToGlobalMappingGetBlockIndices(A->rmap->mapping,&garray);
    2402:   ISCreateBlock(PetscObjectComm((PetscObject)A),rbs,nr/rbs,garray,PETSC_USE_POINTER,&from);
    2403:   VecScatterCreate(rglobal,from,is->y,NULL,&is->rctx);
    2404:   ISLocalToGlobalMappingRestoreBlockIndices(A->rmap->mapping,&garray);
    2405:   ISDestroy(&from);
    2406:   if (A->rmap->mapping != A->cmap->mapping) {
    2407:     ISLocalToGlobalMappingGetBlockIndices(A->cmap->mapping,&garray);
    2408:     ISCreateBlock(PetscObjectComm((PetscObject)A),cbs,nc/cbs,garray,PETSC_USE_POINTER,&from);
    2409:     VecScatterCreate(cglobal,from,is->x,NULL,&is->cctx);
    2410:     ISLocalToGlobalMappingRestoreBlockIndices(A->cmap->mapping,&garray);
    2411:     ISDestroy(&from);
    2412:   } else {
    2413:     PetscObjectReference((PetscObject)is->rctx);
    2414:     is->cctx = is->rctx;
    2415:   }
    2416:   VecDestroy(&cglobal);

    2418:   /* interface counter vector (local) */
    2419:   VecDuplicate(is->y,&is->counter);
    2420:   VecBindToCPU(is->counter,PETSC_TRUE);
    2421:   VecSet(is->y,1.);
    2422:   VecScatterBegin(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);
    2423:   VecScatterEnd(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);
    2424:   VecScatterBegin(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);
    2425:   VecScatterEnd(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);
    2426:   VecBindToCPU(is->y,PETSC_FALSE);
    2427:   VecBindToCPU(is->counter,PETSC_FALSE);

    2429:   /* special functions for block-diagonal matrices */
    2430:   VecSum(rglobal,&sum);
    2431:   if ((PetscInt)(PetscRealPart(sum)) == A->rmap->N && A->rmap->N == A->cmap->N && A->rmap->mapping == A->cmap->mapping) {
    2432:     A->ops->invertblockdiagonal = MatInvertBlockDiagonal_IS;
    2433:   } else {
    2434:     A->ops->invertblockdiagonal = NULL;
    2435:   }
    2436:   VecDestroy(&rglobal);

    2438:   /* setup SF for general purpose shared indices based communications */
    2439:   MatISSetUpSF_IS(A);
    2440:   return(0);
    2441: }

    2443: static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping)
    2444: {
    2446:   PetscInt       nr,rbs,nc,cbs;
    2447:   Mat_IS         *is = (Mat_IS*)A->data;

    2452:   MatDestroy(&is->A);
    2453:   if (is->csf != is->sf) {
    2454:     PetscSFDestroy(&is->csf);
    2455:     PetscFree2(is->csf_rootdata,is->csf_leafdata);
    2456:   } else is->csf = NULL;
    2457:   PetscSFDestroy(&is->sf);
    2458:   PetscFree2(is->sf_rootdata,is->sf_leafdata);
    2459:   PetscFree(is->bdiag);

    2461:   /* Setup Layout and set local to global maps */
    2462:   PetscLayoutSetUp(A->rmap);
    2463:   PetscLayoutSetUp(A->cmap);
    2464:   ISLocalToGlobalMappingGetSize(rmapping,&nr);
    2465:   ISLocalToGlobalMappingGetBlockSize(rmapping,&rbs);
    2466:   ISLocalToGlobalMappingGetSize(cmapping,&nc);
    2467:   ISLocalToGlobalMappingGetBlockSize(cmapping,&cbs);
    2468:   /* check if the two mappings are actually the same for square matrices (DOLFIN passes 2 different objects) */
    2469:   if (rmapping != cmapping && A->rmap->N == A->cmap->N) {
    2470:     PetscBool same,gsame;

    2472:     same = PETSC_FALSE;
    2473:     if (nr == nc && cbs == rbs) {
    2474:       const PetscInt *idxs1,*idxs2;

    2476:       ISLocalToGlobalMappingGetBlockIndices(rmapping,&idxs1);
    2477:       ISLocalToGlobalMappingGetBlockIndices(cmapping,&idxs2);
    2478:       PetscArraycmp(idxs1,idxs2,nr/rbs,&same);
    2479:       ISLocalToGlobalMappingRestoreBlockIndices(rmapping,&idxs1);
    2480:       ISLocalToGlobalMappingRestoreBlockIndices(cmapping,&idxs2);
    2481:     }
    2482:     MPIU_Allreduce(&same,&gsame,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
    2483:     if (gsame) cmapping = rmapping;
    2484:   }
    2485:   PetscLayoutSetBlockSize(A->rmap,rbs);
    2486:   PetscLayoutSetBlockSize(A->cmap,cbs);
    2487:   PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping);
    2488:   PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping);

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

    2500:   if (!is->islocalref) { /* setup scatters and local vectors for MatMult */
    2501:     MatISSetUpScatters_Private(A);
    2502:   }
    2503:   MatSetUp(A);
    2504:   return(0);
    2505: }

    2507: static PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv)
    2508: {
    2509:   Mat_IS         *is = (Mat_IS*)mat->data;
    2511:   PetscInt       i,zm,zn;
    2512:   PetscInt       rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION];

    2515:   if (PetscDefined(USE_DEBUG)) {
    2516:     if (m > MATIS_MAX_ENTRIES_INSERTION || n > MATIS_MAX_ENTRIES_INSERTION) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of row/column indices must be <= %D: they are %D %D",MATIS_MAX_ENTRIES_INSERTION,m,n);
    2517:     /* count negative indices */
    2518:     for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++;
    2519:     for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++;
    2520:   }
    2521:   ISGlobalToLocalMappingApply(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);
    2522:   ISGlobalToLocalMappingApply(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);
    2523:   if (PetscDefined(USE_DEBUG)) {
    2524:     /* count negative indices (should be the same as before) */
    2525:     for (i=0;i<m;i++) if (rows_l[i] < 0) zm--;
    2526:     for (i=0;i<n;i++) if (cols_l[i] < 0) zn--;
    2527:     if (!is->A->rmap->mapping && zm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the row indices can not be mapped! Maybe you should not use MATIS");
    2528:     if (!is->A->cmap->mapping && zn) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the column indices can not be mapped! Maybe you should not use MATIS");
    2529:   }
    2530:   MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);
    2531:   return(0);
    2532: }

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

    2542:   if (PetscDefined(USE_DEBUG)) {
    2543:     if (m > MATIS_MAX_ENTRIES_INSERTION || n > MATIS_MAX_ENTRIES_INSERTION) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of row/column block indices must be <= %D: they are %D %D",MATIS_MAX_ENTRIES_INSERTION,m,n);
    2544:     /* count negative indices */
    2545:     for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++;
    2546:     for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++;
    2547:   }
    2548:   ISGlobalToLocalMappingApplyBlock(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);
    2549:   ISGlobalToLocalMappingApplyBlock(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);
    2550:   if (PetscDefined(USE_DEBUG)) {
    2551:     /* count negative indices (should be the same as before) */
    2552:     for (i=0;i<m;i++) if (rows_l[i] < 0) zm--;
    2553:     for (i=0;i<n;i++) if (cols_l[i] < 0) zn--;
    2554:     if (!is->A->rmap->mapping && zm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the row indices can not be mapped! Maybe you should not use MATIS");
    2555:     if (!is->A->cmap->mapping && zn) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the column indices can not be mapped! Maybe you should not use MATIS");
    2556:   }
    2557:   MatSetValuesBlocked(is->A,m,rows_l,n,cols_l,values,addv);
    2558:   return(0);
    2559: }

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

    2567:   if (is->A->rmap->mapping) {
    2568:     MatSetValuesLocal(is->A,m,rows,n,cols,values,addv);
    2569:   } else {
    2570:     MatSetValues(is->A,m,rows,n,cols,values,addv);
    2571:   }
    2572:   return(0);
    2573: }

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

    2581:   if (is->A->rmap->mapping) {
    2582:     if (PetscDefined(USE_DEBUG)) {
    2583:       PetscInt ibs,bs;

    2585:       ISLocalToGlobalMappingGetBlockSize(is->A->rmap->mapping,&ibs);
    2586:       MatGetBlockSize(is->A,&bs);
    2587:       if (ibs != bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Different block sizes! local mat %D, local l2g map %D",bs,ibs);
    2588:     }
    2589:     MatSetValuesBlockedLocal(is->A,m,rows,n,cols,values,addv);
    2590:   } else {
    2591:     MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv);
    2592:   }
    2593:   return(0);
    2594: }

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

    2602:   if (!n) {
    2603:     is->pure_neumann = PETSC_TRUE;
    2604:   } else {
    2605:     PetscInt i;
    2606:     is->pure_neumann = PETSC_FALSE;

    2608:     if (columns) {
    2609:       MatZeroRowsColumns(is->A,n,rows,diag,NULL,NULL);
    2610:     } else {
    2611:       MatZeroRows(is->A,n,rows,diag,NULL,NULL);
    2612:     }
    2613:     if (diag != 0.) {
    2614:       const PetscScalar *array;
    2615:       VecGetArrayRead(is->counter,&array);
    2616:       for (i=0; i<n; i++) {
    2617:         MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);
    2618:       }
    2619:       VecRestoreArrayRead(is->counter,&array);
    2620:     }
    2621:     MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);
    2622:     MatAssemblyEnd(is->A,MAT_FINAL_ASSEMBLY);
    2623:   }
    2624:   return(0);
    2625: }

    2627: static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b,PetscBool columns)
    2628: {
    2629:   Mat_IS         *matis = (Mat_IS*)A->data;
    2630:   PetscInt       nr,nl,len,i;
    2631:   PetscInt       *lrows;

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

    2638:     PetscLayoutCompare(A->rmap,A->cmap,&cong);
    2639:     cong = (PetscBool)(cong && matis->sf == matis->csf);
    2640:     if (!cong && columns) SETERRQ(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");
    2641:     if (!cong && diag != 0.) SETERRQ(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");
    2642:     if (!cong && x && b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"A->rmap and A->cmap need to be congruent, and the l2g maps be the same");
    2643:   }
    2644:   /* get locally owned rows */
    2645:   PetscLayoutMapLocal(A->rmap,n,rows,&len,&lrows,NULL);
    2646:   /* fix right hand side if needed */
    2647:   if (x && b) {
    2648:     const PetscScalar *xx;
    2649:     PetscScalar       *bb;

    2651:     VecGetArrayRead(x, &xx);
    2652:     VecGetArray(b, &bb);
    2653:     for (i=0;i<len;++i) bb[lrows[i]] = diag*xx[lrows[i]];
    2654:     VecRestoreArrayRead(x, &xx);
    2655:     VecRestoreArray(b, &bb);
    2656:   }
    2657:   /* get rows associated to the local matrices */
    2658:   MatGetSize(matis->A,&nl,NULL);
    2659:   PetscArrayzero(matis->sf_leafdata,nl);
    2660:   PetscArrayzero(matis->sf_rootdata,A->rmap->n);
    2661:   for (i=0;i<len;i++) matis->sf_rootdata[lrows[i]] = 1;
    2662:   PetscFree(lrows);
    2663:   PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE);
    2664:   PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE);
    2665:   PetscMalloc1(nl,&lrows);
    2666:   for (i=0,nr=0;i<nl;i++) if (matis->sf_leafdata[i]) lrows[nr++] = i;
    2667:   MatISZeroRowsColumnsLocal_Private(A,nr,lrows,diag,columns);
    2668:   PetscFree(lrows);
    2669:   return(0);
    2670: }

    2672: static PetscErrorCode MatZeroRows_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
    2673: {

    2677:   MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_FALSE);
    2678:   return(0);
    2679: }

    2681: static PetscErrorCode MatZeroRowsColumns_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
    2682: {

    2686:   MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_TRUE);
    2687:   return(0);
    2688: }

    2690: static PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type)
    2691: {
    2692:   Mat_IS         *is = (Mat_IS*)A->data;

    2696:   MatAssemblyBegin(is->A,type);
    2697:   return(0);
    2698: }

    2700: static PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type)
    2701: {
    2702:   Mat_IS         *is = (Mat_IS*)A->data;

    2706:   MatAssemblyEnd(is->A,type);
    2707:   /* fix for local empty rows/cols */
    2708:   if (is->locempty && type == MAT_FINAL_ASSEMBLY) {
    2709:     Mat                    newlA;
    2710:     ISLocalToGlobalMapping rl2g,cl2g;
    2711:     IS                     nzr,nzc;
    2712:     PetscInt               nr,nc,nnzr,nnzc;
    2713:     PetscBool              lnewl2g,newl2g;

    2715:     MatGetSize(is->A,&nr,&nc);
    2716:     MatFindNonzeroRowsOrCols_Basic(is->A,PETSC_FALSE,PETSC_SMALL,&nzr);
    2717:     if (!nzr) {
    2718:       ISCreateStride(PetscObjectComm((PetscObject)is->A),nr,0,1,&nzr);
    2719:     }
    2720:     MatFindNonzeroRowsOrCols_Basic(is->A,PETSC_TRUE,PETSC_SMALL,&nzc);
    2721:     if (!nzc) {
    2722:       ISCreateStride(PetscObjectComm((PetscObject)is->A),nc,0,1,&nzc);
    2723:     }
    2724:     ISGetSize(nzr,&nnzr);
    2725:     ISGetSize(nzc,&nnzc);
    2726:     if (nnzr != nr || nnzc != nc) {
    2727:       ISLocalToGlobalMapping l2g;
    2728:       IS                     is1,is2;

    2730:       /* need new global l2g map */
    2731:       lnewl2g = PETSC_TRUE;
    2732:       MPI_Allreduce(&lnewl2g,&newl2g,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A));

    2734:       /* extract valid submatrix */
    2735:       MatCreateSubMatrix(is->A,nzr,nzc,MAT_INITIAL_MATRIX,&newlA);

    2737:       /* attach local l2g maps for successive calls of MatSetValues on the local matrix */
    2738:       ISLocalToGlobalMappingCreateIS(nzr,&l2g);
    2739:       ISCreateStride(PetscObjectComm((PetscObject)is->A),nr,0,1,&is1);
    2740:       ISGlobalToLocalMappingApplyIS(l2g,IS_GTOLM_MASK,is1,&is2);
    2741:       ISLocalToGlobalMappingDestroy(&l2g);
    2742:       if (is->A->rmap->mapping) { /* the matrix has a local-to-local map already attached (probably DMDA generated) */
    2743:         const PetscInt *idxs1,*idxs2;
    2744:         PetscInt       j,i,nl,*tidxs;

    2746:         ISLocalToGlobalMappingGetSize(is->A->rmap->mapping,&nl);
    2747:         ISLocalToGlobalMappingGetIndices(is->A->rmap->mapping,&idxs1);
    2748:         PetscMalloc1(nl,&tidxs);
    2749:         ISGetIndices(is2,&idxs2);
    2750:         for (i=0,j=0;i<nl;i++) tidxs[i] = idxs1[i] < 0 ? -1 : idxs2[j++];
    2751:         if (j != nr) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected count %D != %D",j,nr);
    2752:         ISLocalToGlobalMappingRestoreIndices(is->A->rmap->mapping,&idxs1);
    2753:         ISRestoreIndices(is2,&idxs2);
    2754:         ISDestroy(&is2);
    2755:         ISCreateGeneral(PetscObjectComm((PetscObject)is->A->rmap->mapping),nl,tidxs,PETSC_OWN_POINTER,&is2);
    2756:       }
    2757:       ISLocalToGlobalMappingCreateIS(is2,&rl2g);
    2758:       ISDestroy(&is1);
    2759:       ISDestroy(&is2);

    2761:       ISLocalToGlobalMappingCreateIS(nzc,&l2g);
    2762:       ISCreateStride(PetscObjectComm((PetscObject)is->A),nc,0,1,&is1);
    2763:       ISGlobalToLocalMappingApplyIS(l2g,IS_GTOLM_MASK,is1,&is2);
    2764:       ISLocalToGlobalMappingDestroy(&l2g);
    2765:       if (is->A->cmap->mapping) { /* the matrix has a local-to-local map already attached (probably DMDA generated) */
    2766:         const PetscInt *idxs1,*idxs2;
    2767:         PetscInt       j,i,nl,*tidxs;

    2769:         ISLocalToGlobalMappingGetSize(is->A->cmap->mapping,&nl);
    2770:         ISLocalToGlobalMappingGetIndices(is->A->cmap->mapping,&idxs1);
    2771:         PetscMalloc1(nl,&tidxs);
    2772:         ISGetIndices(is2,&idxs2);
    2773:         for (i=0,j=0;i<nl;i++) tidxs[i] = idxs1[i] < 0 ? -1 : idxs2[j++];
    2774:         if (j != nc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected count %D != %D",j,nc);
    2775:         ISLocalToGlobalMappingRestoreIndices(is->A->cmap->mapping,&idxs1);
    2776:         ISRestoreIndices(is2,&idxs2);
    2777:         ISDestroy(&is2);
    2778:         ISCreateGeneral(PetscObjectComm((PetscObject)is->A->cmap->mapping),nl,tidxs,PETSC_OWN_POINTER,&is2);
    2779:       }
    2780:       ISLocalToGlobalMappingCreateIS(is2,&cl2g);
    2781:       ISDestroy(&is1);
    2782:       ISDestroy(&is2);

    2784:       MatSetLocalToGlobalMapping(newlA,rl2g,cl2g);

    2786:       ISLocalToGlobalMappingDestroy(&rl2g);
    2787:       ISLocalToGlobalMappingDestroy(&cl2g);
    2788:     } else { /* local matrix fully populated */
    2789:       lnewl2g = PETSC_FALSE;
    2790:       MPI_Allreduce(&lnewl2g,&newl2g,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A));
    2791:       PetscObjectReference((PetscObject)is->A);
    2792:       newlA   = is->A;
    2793:     }
    2794:     /* attach new global l2g map if needed */
    2795:     if (newl2g) {
    2796:       IS             gnzr,gnzc;
    2797:       const PetscInt *grid,*gcid;

    2799:       ISLocalToGlobalMappingApplyIS(A->rmap->mapping,nzr,&gnzr);
    2800:       ISLocalToGlobalMappingApplyIS(A->cmap->mapping,nzc,&gnzc);
    2801:       ISGetIndices(gnzr,&grid);
    2802:       ISGetIndices(gnzc,&gcid);
    2803:       ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,nnzr,grid,PETSC_COPY_VALUES,&rl2g);
    2804:       ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,nnzc,gcid,PETSC_COPY_VALUES,&cl2g);
    2805:       ISRestoreIndices(gnzr,&grid);
    2806:       ISRestoreIndices(gnzc,&gcid);
    2807:       ISDestroy(&gnzr);
    2808:       ISDestroy(&gnzc);
    2809:       MatSetLocalToGlobalMapping(A,rl2g,cl2g);
    2810:       ISLocalToGlobalMappingDestroy(&rl2g);
    2811:       ISLocalToGlobalMappingDestroy(&cl2g);
    2812:     }
    2813:     MatISSetLocalMat(A,newlA);
    2814:     MatDestroy(&newlA);
    2815:     ISDestroy(&nzr);
    2816:     ISDestroy(&nzc);
    2817:     is->locempty = PETSC_FALSE;
    2818:   }
    2819:   return(0);
    2820: }

    2822: static PetscErrorCode MatISGetLocalMat_IS(Mat mat,Mat *local)
    2823: {
    2824:   Mat_IS *is = (Mat_IS*)mat->data;

    2827:   *local = is->A;
    2828:   return(0);
    2829: }

    2831: static PetscErrorCode MatISRestoreLocalMat_IS(Mat mat,Mat *local)
    2832: {
    2834:   *local = NULL;
    2835:   return(0);
    2836: }

    2838: /*@
    2839:     MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix.

    2841:   Input Parameter:
    2842: .  mat - the matrix

    2844:   Output Parameter:
    2845: .  local - the local matrix

    2847:   Level: advanced

    2849:   Notes:
    2850:     This can be called if you have precomputed the nonzero structure of the
    2851:   matrix and want to provide it to the inner matrix object to improve the performance
    2852:   of the MatSetValues() operation.

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

    2856: .seealso: MATIS
    2857: @*/
    2858: PetscErrorCode MatISGetLocalMat(Mat mat,Mat *local)
    2859: {

    2865:   PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat*),(mat,local));
    2866:   return(0);
    2867: }

    2869: /*@
    2870:     MatISRestoreLocalMat - Restores the local matrix obtained with MatISGetLocalMat()

    2872:   Input Parameter:
    2873: .  mat - the matrix

    2875:   Output Parameter:
    2876: .  local - the local matrix

    2878:   Level: advanced

    2880: .seealso: MATIS
    2881: @*/
    2882: PetscErrorCode MatISRestoreLocalMat(Mat mat,Mat *local)
    2883: {

    2889:   PetscUseMethod(mat,"MatISRestoreLocalMat_C",(Mat,Mat*),(mat,local));
    2890:   return(0);
    2891: }

    2893: static PetscErrorCode MatISSetLocalMatType_IS(Mat mat,MatType mtype)
    2894: {
    2895:   Mat_IS         *is = (Mat_IS*)mat->data;

    2899:   if (is->A) {
    2900:     MatSetType(is->A,mtype);
    2901:   }
    2902:   PetscFree(is->lmattype);
    2903:   PetscStrallocpy(mtype,&is->lmattype);
    2904:   return(0);
    2905: }

    2907: /*@
    2908:     MatISSetLocalMatType - Specifies the type of local matrix

    2910:   Input Parameter:
    2911: +  mat - the matrix
    2912: -  mtype - the local matrix type

    2914:   Output Parameter:

    2916:   Level: advanced

    2918: .seealso: MATIS, MatSetType(), MatType
    2919: @*/
    2920: PetscErrorCode MatISSetLocalMatType(Mat mat,MatType mtype)
    2921: {

    2926:   PetscUseMethod(mat,"MatISSetLocalMatType_C",(Mat,MatType),(mat,mtype));
    2927:   return(0);
    2928: }

    2930: static PetscErrorCode MatISSetLocalMat_IS(Mat mat,Mat local)
    2931: {
    2932:   Mat_IS         *is = (Mat_IS*)mat->data;
    2933:   PetscInt       nrows,ncols,orows,ocols;
    2935:   MatType        mtype,otype;
    2936:   PetscBool      sametype = PETSC_TRUE;

    2939:   if (is->A) {
    2940:     MatGetSize(is->A,&orows,&ocols);
    2941:     MatGetSize(local,&nrows,&ncols);
    2942:     if (orows != nrows || ocols != ncols) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local MATIS matrix should be of size %Dx%D (you passed a %Dx%D matrix)",orows,ocols,nrows,ncols);
    2943:     MatGetType(local,&mtype);
    2944:     MatGetType(is->A,&otype);
    2945:     PetscStrcmp(mtype,otype,&sametype);
    2946:   }
    2947:   PetscObjectReference((PetscObject)local);
    2948:   MatDestroy(&is->A);
    2949:   is->A = local;
    2950:   MatGetType(is->A,&mtype);
    2951:   MatISSetLocalMatType(mat,mtype);
    2952:   if (!sametype && !is->islocalref) {
    2953:     MatISSetUpScatters_Private(mat);
    2954:   }
    2955:   return(0);
    2956: }

    2958: /*@
    2959:     MatISSetLocalMat - Replace the local matrix stored inside a MATIS object.

    2961:   Collective on Mat

    2963:   Input Parameter:
    2964: +  mat - the matrix
    2965: -  local - the local matrix

    2967:   Output Parameter:

    2969:   Level: advanced

    2971:   Notes:
    2972:     This can be called if you have precomputed the local matrix and
    2973:   want to provide it to the matrix object MATIS.

    2975: .seealso: MATIS
    2976: @*/
    2977: PetscErrorCode MatISSetLocalMat(Mat mat,Mat local)
    2978: {

    2984:   PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local));
    2985:   return(0);
    2986: }

    2988: static PetscErrorCode MatZeroEntries_IS(Mat A)
    2989: {
    2990:   Mat_IS         *a = (Mat_IS*)A->data;

    2994:   MatZeroEntries(a->A);
    2995:   return(0);
    2996: }

    2998: static PetscErrorCode MatScale_IS(Mat A,PetscScalar a)
    2999: {
    3000:   Mat_IS         *is = (Mat_IS*)A->data;

    3004:   MatScale(is->A,a);
    3005:   return(0);
    3006: }

    3008: static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v)
    3009: {
    3010:   Mat_IS         *is = (Mat_IS*)A->data;

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

    3017:   /* scatter diagonal back into global vector */
    3018:   VecSet(v,0);
    3019:   VecScatterBegin(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);
    3020:   VecScatterEnd(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);
    3021:   return(0);
    3022: }

    3024: static PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg)
    3025: {
    3026:   Mat_IS         *a = (Mat_IS*)A->data;

    3030:   MatSetOption(a->A,op,flg);
    3031:   return(0);
    3032: }

    3034: static PetscErrorCode MatAXPY_IS(Mat Y,PetscScalar a,Mat X,MatStructure str)
    3035: {
    3036:   Mat_IS         *y = (Mat_IS*)Y->data;
    3037:   Mat_IS         *x;

    3041:   if (PetscDefined(USE_DEBUG)) {
    3042:     PetscBool      ismatis;
    3043:     PetscObjectTypeCompare((PetscObject)X,MATIS,&ismatis);
    3044:     if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot call MatAXPY(Y,a,X,str) with X not of type MATIS");
    3045:   }
    3046:   x = (Mat_IS*)X->data;
    3047:   MatAXPY(y->A,a,x->A,str);
    3048:   return(0);
    3049: }

    3051: static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A,IS row,IS col,Mat *submat)
    3052: {
    3053:   Mat                    lA;
    3054:   Mat_IS                 *matis;
    3055:   ISLocalToGlobalMapping rl2g,cl2g;
    3056:   IS                     is;
    3057:   const PetscInt         *rg,*rl;
    3058:   PetscInt               nrg;
    3059:   PetscInt               N,M,nrl,i,*idxs;
    3060:   PetscErrorCode         ierr;

    3063:   ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&rg);
    3064:   ISGetLocalSize(row,&nrl);
    3065:   ISGetIndices(row,&rl);
    3066:   ISLocalToGlobalMappingGetSize(A->rmap->mapping,&nrg);
    3067:   if (PetscDefined(USE_DEBUG)) {
    3068:     for (i=0; i<nrl; i++) if (rl[i]>=nrg) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row index %D -> %D greater then maximum possible %D",i,rl[i],nrg);
    3069:   }
    3070:   PetscMalloc1(nrg,&idxs);
    3071:   /* map from [0,nrl) to row */
    3072:   for (i=0;i<nrl;i++) idxs[i] = rl[i];
    3073:   for (i=nrl;i<nrg;i++) idxs[i] = -1;
    3074:   ISRestoreIndices(row,&rl);
    3075:   ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&rg);
    3076:   ISCreateGeneral(PetscObjectComm((PetscObject)A),nrg,idxs,PETSC_OWN_POINTER,&is);
    3077:   ISLocalToGlobalMappingCreateIS(is,&rl2g);
    3078:   ISDestroy(&is);
    3079:   /* compute new l2g map for columns */
    3080:   if (col != row || A->rmap->mapping != A->cmap->mapping) {
    3081:     const PetscInt *cg,*cl;
    3082:     PetscInt       ncg;
    3083:     PetscInt       ncl;

    3085:     ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&cg);
    3086:     ISGetLocalSize(col,&ncl);
    3087:     ISGetIndices(col,&cl);
    3088:     ISLocalToGlobalMappingGetSize(A->cmap->mapping,&ncg);
    3089:     if (PetscDefined(USE_DEBUG)) {
    3090:       for (i=0; i<ncl; i++) if (cl[i]>=ncg) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local column index %D -> %D greater then maximum possible %D",i,cl[i],ncg);
    3091:     }
    3092:     PetscMalloc1(ncg,&idxs);
    3093:     /* map from [0,ncl) to col */
    3094:     for (i=0;i<ncl;i++) idxs[i] = cl[i];
    3095:     for (i=ncl;i<ncg;i++) idxs[i] = -1;
    3096:     ISRestoreIndices(col,&cl);
    3097:     ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&cg);
    3098:     ISCreateGeneral(PetscObjectComm((PetscObject)A),ncg,idxs,PETSC_OWN_POINTER,&is);
    3099:     ISLocalToGlobalMappingCreateIS(is,&cl2g);
    3100:     ISDestroy(&is);
    3101:   } else {
    3102:     PetscObjectReference((PetscObject)rl2g);
    3103:     cl2g = rl2g;
    3104:   }
    3105:   /* create the MATIS submatrix */
    3106:   MatGetSize(A,&M,&N);
    3107:   MatCreate(PetscObjectComm((PetscObject)A),submat);
    3108:   MatSetSizes(*submat,PETSC_DECIDE,PETSC_DECIDE,M,N);
    3109:   MatSetType(*submat,MATIS);
    3110:   matis = (Mat_IS*)((*submat)->data);
    3111:   matis->islocalref = PETSC_TRUE;
    3112:   MatSetLocalToGlobalMapping(*submat,rl2g,cl2g);
    3113:   MatISGetLocalMat(A,&lA);
    3114:   MatISSetLocalMat(*submat,lA);
    3115:   MatSetUp(*submat);
    3116:   MatAssemblyBegin(*submat,MAT_FINAL_ASSEMBLY);
    3117:   MatAssemblyEnd(*submat,MAT_FINAL_ASSEMBLY);
    3118:   ISLocalToGlobalMappingDestroy(&rl2g);
    3119:   ISLocalToGlobalMappingDestroy(&cl2g);
    3120:   /* remove unsupported ops */
    3121:   PetscMemzero((*submat)->ops,sizeof(struct _MatOps));
    3122:   (*submat)->ops->destroy               = MatDestroy_IS;
    3123:   (*submat)->ops->setvalueslocal        = MatSetValuesLocal_SubMat_IS;
    3124:   (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS;
    3125:   (*submat)->ops->assemblybegin         = MatAssemblyBegin_IS;
    3126:   (*submat)->ops->assemblyend           = MatAssemblyEnd_IS;
    3127:   return(0);
    3128: }

    3130: static PetscErrorCode MatSetFromOptions_IS(PetscOptionItems *PetscOptionsObject, Mat A)
    3131: {
    3132:   Mat_IS         *a = (Mat_IS*)A->data;
    3133:   char           type[256];
    3134:   PetscBool      flg;

    3138:   PetscOptionsHead(PetscOptionsObject,"MATIS options");
    3139:   PetscOptionsBool("-matis_fixempty","Fix local matrices in case of empty local rows/columns","MatISFixLocalEmpty",a->locempty,&a->locempty,NULL);
    3140:   PetscOptionsBool("-matis_storel2l","Store local-to-local matrices generated from PtAP operations","MatISStoreL2L",a->storel2l,&a->storel2l,NULL);
    3141:   PetscOptionsFList("-matis_localmat_type","Matrix type","MatISSetLocalMatType",MatList,a->lmattype,type,256,&flg);
    3142:   if (flg) {
    3143:     MatISSetLocalMatType(A,type);
    3144:   }
    3145:   if (a->A) {
    3146:     MatSetFromOptions(a->A);
    3147:   }
    3148:   PetscOptionsTail();
    3149:   return(0);
    3150: }

    3152: /*@
    3153:     MatCreateIS - Creates a "process" unassembled matrix, assembled on each
    3154:        process but not across processes.

    3156:    Input Parameters:
    3157: +     comm    - MPI communicator that will share the matrix
    3158: .     bs      - block size of the matrix
    3159: .     m,n,M,N - local and/or global sizes of the left and right vector used in matrix vector products
    3160: .     rmap    - local to global map for rows
    3161: -     cmap    - local to global map for cols

    3163:    Output Parameter:
    3164: .    A - the resulting matrix

    3166:    Level: advanced

    3168:    Notes:
    3169:     See MATIS for more details.
    3170:           m and n are NOT related to the size of the map; they represent the size of the local parts of the vectors
    3171:           used in MatMult operations. The sizes of rmap and cmap define the size of the local matrices.
    3172:           If either rmap or cmap are NULL, then the matrix is assumed to be square.

    3174: .seealso: MATIS, MatSetLocalToGlobalMapping()
    3175: @*/
    3176: PetscErrorCode  MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping rmap,ISLocalToGlobalMapping cmap,Mat *A)
    3177: {

    3181:   if (!rmap && !cmap) SETERRQ(comm,PETSC_ERR_USER,"You need to provide at least one of the mappings");
    3182:   MatCreate(comm,A);
    3183:   MatSetSizes(*A,m,n,M,N);
    3184:   if (bs > 0) {
    3185:     MatSetBlockSize(*A,bs);
    3186:   }
    3187:   MatSetType(*A,MATIS);
    3188:   if (rmap && cmap) {
    3189:     MatSetLocalToGlobalMapping(*A,rmap,cmap);
    3190:   } else if (!rmap) {
    3191:     MatSetLocalToGlobalMapping(*A,cmap,cmap);
    3192:   } else {
    3193:     MatSetLocalToGlobalMapping(*A,rmap,rmap);
    3194:   }
    3195:   return(0);
    3196: }

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

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

    3208:    Notes:
    3209:     Options prefix for the inner matrix are given by -is_mat_xxx

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

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

    3216:   Level: advanced

    3218: .seealso: Mat, MatISGetLocalMat(), MatSetLocalToGlobalMapping(), MatISSetPreallocation(), MatCreateIS(), PCBDDC, KSPFETIDP

    3220: M*/

    3222: PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A)
    3223: {
    3225:   Mat_IS         *b;

    3228:   PetscNewLog(A,&b);
    3229:   PetscStrallocpy(MATAIJ,&b->lmattype);
    3230:   A->data = (void*)b;

    3232:   /* matrix ops */
    3233:   PetscMemzero(A->ops,sizeof(struct _MatOps));
    3234:   A->ops->mult                    = MatMult_IS;
    3235:   A->ops->multadd                 = MatMultAdd_IS;
    3236:   A->ops->multtranspose           = MatMultTranspose_IS;
    3237:   A->ops->multtransposeadd        = MatMultTransposeAdd_IS;
    3238:   A->ops->destroy                 = MatDestroy_IS;
    3239:   A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS;
    3240:   A->ops->setvalues               = MatSetValues_IS;
    3241:   A->ops->setvaluesblocked        = MatSetValuesBlocked_IS;
    3242:   A->ops->setvalueslocal          = MatSetValuesLocal_IS;
    3243:   A->ops->setvaluesblockedlocal   = MatSetValuesBlockedLocal_IS;
    3244:   A->ops->zerorows                = MatZeroRows_IS;
    3245:   A->ops->zerorowscolumns         = MatZeroRowsColumns_IS;
    3246:   A->ops->assemblybegin           = MatAssemblyBegin_IS;
    3247:   A->ops->assemblyend             = MatAssemblyEnd_IS;
    3248:   A->ops->view                    = MatView_IS;
    3249:   A->ops->zeroentries             = MatZeroEntries_IS;
    3250:   A->ops->scale                   = MatScale_IS;
    3251:   A->ops->getdiagonal             = MatGetDiagonal_IS;
    3252:   A->ops->setoption               = MatSetOption_IS;
    3253:   A->ops->ishermitian             = MatIsHermitian_IS;
    3254:   A->ops->issymmetric             = MatIsSymmetric_IS;
    3255:   A->ops->isstructurallysymmetric = MatIsStructurallySymmetric_IS;
    3256:   A->ops->duplicate               = MatDuplicate_IS;
    3257:   A->ops->missingdiagonal         = MatMissingDiagonal_IS;
    3258:   A->ops->copy                    = MatCopy_IS;
    3259:   A->ops->getlocalsubmatrix       = MatGetLocalSubMatrix_IS;
    3260:   A->ops->createsubmatrix         = MatCreateSubMatrix_IS;
    3261:   A->ops->axpy                    = MatAXPY_IS;
    3262:   A->ops->diagonalset             = MatDiagonalSet_IS;
    3263:   A->ops->shift                   = MatShift_IS;
    3264:   A->ops->transpose               = MatTranspose_IS;
    3265:   A->ops->getinfo                 = MatGetInfo_IS;
    3266:   A->ops->diagonalscale           = MatDiagonalScale_IS;
    3267:   A->ops->setfromoptions          = MatSetFromOptions_IS;

    3269:   /* special MATIS functions */
    3270:   PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMatType_C",MatISSetLocalMatType_IS);
    3271:   PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS);
    3272:   PetscObjectComposeFunction((PetscObject)A,"MatISRestoreLocalMat_C",MatISRestoreLocalMat_IS);
    3273:   PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS);
    3274:   PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",MatConvert_IS_XAIJ);
    3275:   PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",MatISSetPreallocation_IS);
    3276:   PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",MatISStoreL2L_IS);
    3277:   PetscObjectComposeFunction((PetscObject)A,"MatISFixLocalEmpty_C",MatISFixLocalEmpty_IS);
    3278:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpiaij_C",MatConvert_IS_XAIJ);
    3279:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpibaij_C",MatConvert_IS_XAIJ);
    3280:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpisbaij_C",MatConvert_IS_XAIJ);
    3281:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqaij_C",MatConvert_IS_XAIJ);
    3282:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqbaij_C",MatConvert_IS_XAIJ);
    3283:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqsbaij_C",MatConvert_IS_XAIJ);
    3284:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_aij_C",MatConvert_IS_XAIJ);
    3285:   PetscObjectChangeTypeName((PetscObject)A,MATIS);
    3286:   return(0);
    3287: }