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>
 13: #include <petsc/private/hashseti.h>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

241:   C->ops->ptapnumeric = MatPtAPNumeric_IS_XAIJ;
242:   return 0;
243: }

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

252:   MatPtAPSymbolic_IS_XAIJ(A,P,fill,C);
253:   C->ops->productnumeric = MatProductNumeric_PtAP;
254:   return 0;
255: }

257: static PetscErrorCode MatProductSetFromOptions_IS_XAIJ_PtAP(Mat C)
258: {
259:   C->ops->productsymbolic = MatProductSymbolic_PtAP_IS_XAIJ;
260:   return 0;
261: }

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

267:   if (product->type == MATPRODUCT_PtAP) {
268:     MatProductSetFromOptions_IS_XAIJ_PtAP(C);
269:   }
270:   return 0;
271: }

273: /* ----------------------------------------- */
274: static PetscErrorCode MatISContainerDestroyFields_Private(void *ptr)
275: {
276:   MatISLocalFields lf = (MatISLocalFields)ptr;
277:   PetscInt         i;

279:   for (i=0;i<lf->nr;i++) {
280:     ISDestroy(&lf->rf[i]);
281:   }
282:   for (i=0;i<lf->nc;i++) {
283:     ISDestroy(&lf->cf[i]);
284:   }
285:   PetscFree2(lf->rf,lf->cf);
286:   PetscFree(lf);
287:   return 0;
288: }

290: static PetscErrorCode MatConvert_SeqXAIJ_IS(Mat A,MatType type,MatReuse reuse,Mat *newmat)
291: {
292:   Mat            B,lB;

294:   if (reuse != MAT_REUSE_MATRIX) {
295:     ISLocalToGlobalMapping rl2g,cl2g;
296:     PetscInt               bs;
297:     IS                     is;

299:     MatGetBlockSize(A,&bs);
300:     ISCreateStride(PetscObjectComm((PetscObject)A),A->rmap->n/bs,0,1,&is);
301:     if (bs > 1) {
302:       IS       is2;
303:       PetscInt i,*aux;

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

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

350: static PetscErrorCode MatISScaleDisassembling_Private(Mat A)
351: {
352:   Mat_IS         *matis = (Mat_IS*)(A->data);
353:   PetscScalar    *aa;
354:   const PetscInt *ii,*jj;
355:   PetscInt       i,n,m;
356:   PetscInt       *ecount,**eneighs;
357:   PetscBool      flg;

359:   MatGetRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&m,&ii,&jj,&flg);
361:   ISLocalToGlobalMappingGetNodeInfo(matis->rmapping,&n,&ecount,&eneighs);
363:   MatSeqAIJGetArray(matis->A,&aa);
364:   for (i=0;i<n;i++) {
365:     if (ecount[i] > 1) {
366:       PetscInt j;

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

372:         for (p=0;p<ecount[i];p++) {
373:           for (p2=0;p2<ecount[i2];p2++) {
374:             if (eneighs[i][p] == eneighs[i2][p2]) { scal += 1.0; break; }
375:           }
376:         }
377:         if (scal) aa[j] /= scal;
378:       }
379:     }
380:   }
381:   ISLocalToGlobalMappingRestoreNodeInfo(matis->rmapping,&n,&ecount,&eneighs);
382:   MatSeqAIJRestoreArray(matis->A,&aa);
383:   MatRestoreRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&m,&ii,&jj,&flg);
385:   return 0;
386: }

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

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

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

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

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

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

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

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

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

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

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

490:           ISGetLocalSize(workis[j],&s);
491:           ISGetIndices(workis[j],&idxs);
492:           MatPreallocateSet(i+A->rmap->rstart,s,idxs,dnz,onz);
493:           j++;
494:         }
495:       }

498:       for (i = 0; i < gcnt; i++) {
499:         PetscObjectSetName((PetscObject)workis[i],"EXTENDED");
500:         ISViewFromOptions(workis[i],NULL,"-view_isolated_separators");
501:       }

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

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

516:         if (s) {
517:           ISGetIndices(workis[j],&idxs);
518:           MatSetValues(A2,1,&row,s,idxs,vals,INSERT_VALUES);
519:           ISRestoreIndices(workis[j],&idxs);
520:           j++;
521:         }
522:       }
524:       PetscFree(vals);
525:       MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY);
526:       MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY);
527:       MatTranspose(A2,MAT_INPLACE_MATRIX,&A2);

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

537:       /* extend local to global map to include connected isolated separators */
538:       PetscObjectQuery((PetscObject)A3,"_petsc_GetLocalMatCondensed_iscol",(PetscObject*)&is);
540:       ISLocalToGlobalMappingCreateIS(is,&ll2g);
541:       MatGetRowIJ(A3,0,PETSC_FALSE,PETSC_FALSE,&i,&ii,&jj,&flg);
543:       ISCreateGeneral(PETSC_COMM_SELF,ii[i],jj,PETSC_COPY_VALUES,&is);
544:       MatRestoreRowIJ(A3,0,PETSC_FALSE,PETSC_FALSE,&i,&ii,&jj,&flg);
546:       ISLocalToGlobalMappingApplyIS(ll2g,is,&is2);
547:       ISDestroy(&is);
548:       ISLocalToGlobalMappingDestroy(&ll2g);

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

557:       MatDestroy(&A3);
558:       PetscFree(lndmapi);
559:       MatPreallocateFinalize(dnz,onz);
560:       for (i = 0; i < gcnt; i++) {
561:         ISDestroy(&workis[i]);
562:       }
563:       PetscFree(workis);
564:     }
565:     ISRestoreIndices(ndmap,&ndmapi);
566:     PetscSFDestroy(&sf);
567:     PetscFree(ndmapc);
568:     ISDestroy(&ndmap);
569:     ISDestroy(&ndsub);
570:     ISLocalToGlobalMappingSetBlockSize(*l2g,bs);
571:     ISLocalToGlobalMappingViewFromOptions(*l2g,NULL,"-matis_nd_l2g_view");
572:     break;
573:   case MAT_IS_DISASSEMBLE_L2G_NATURAL:
574:     if (ismpiaij) {
575:       MatMPIAIJGetSeqAIJ(A,&Ad,&Ao,&garray);
576:     } else if (ismpibaij) {
577:       MatMPIBAIJGetSeqBAIJ(A,&Ad,&Ao,&garray);
578:     } else SETERRQ(comm,PETSC_ERR_SUP,"Type %s",((PetscObject)A)->type_name);
580:     if (A->rmap->n) {
581:       PetscInt dc,oc,stc,*aux;

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

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

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

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

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

708:   /* generate l2g maps for rows and cols */
709:   ISCreateStride(comm,dr/bs,str/bs,1,&is);
710:   if (bs > 1) {
711:     IS is2;

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

735:   /* create MATIS object */
736:   MatCreate(comm,&B);
737:   MatSetSizes(B,dr,dc,PETSC_DECIDE,PETSC_DECIDE);
738:   MatSetType(B,MATIS);
739:   MatSetBlockSize(B,bs);
740:   MatSetLocalToGlobalMapping(B,rl2g,cl2g);
741:   ISLocalToGlobalMappingDestroy(&rl2g);
742:   ISLocalToGlobalMappingDestroy(&cl2g);

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

759:   MatRestoreRowIJ(Ad,0,PETSC_FALSE,PETSC_FALSE,&i,&odi,&odj,&flg);
761:   MatRestoreRowIJ(Ao,0,PETSC_FALSE,PETSC_FALSE,&i,&ooi,&ooj,&flg);
763:   MatSeqAIJRestoreArray(Ad,&dd);
764:   MatSeqAIJRestoreArray(Ao,&od);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1092:     MatISGetLocalMat(*newmat,&lA);
1093:     MatConvert(lA,MATAIJ,MAT_INITIAL_MATRIX,&M);
1094:     MatISSetLocalMat(*newmat,M);
1095:     MatDestroy(&M);

1097:     /* attach local fields to the matrix */
1098:     PetscNew(&lf);
1099:     PetscMalloc2(nr,&lf->rf,nc,&lf->cf);
1100:     for (i=0;i<nr;i++) {
1101:       PetscInt n,st;

1103:       ISGetLocalSize(islrow[i],&n);
1104:       ISStrideGetInfo(islrow[i],&st,NULL);
1105:       ISCreateStride(comm,n,st,1,&lf->rf[i]);
1106:     }
1107:     for (i=0;i<nc;i++) {
1108:       PetscInt n,st;

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

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

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

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

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

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

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

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

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

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

1242:   /* perform local transposition */
1243:   MatISGetLocalMat(A,&lA);
1244:   MatTranspose(lA,MAT_INITIAL_MATRIX,&lC);
1245:   MatSetLocalToGlobalMapping(lC,lA->cmap->mapping,lA->rmap->mapping);
1246:   MatISSetLocalMat(C,lC);
1247:   MatDestroy(&lC);

1249:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
1250:     *B = C;
1251:   } else {
1252:     MatHeaderMerge(A,&C);
1253:   }
1254:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
1255:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
1256:   return 0;
1257: }

1259: static PetscErrorCode MatDiagonalSet_IS(Mat A,Vec D,InsertMode insmode)
1260: {
1261:   Mat_IS         *is = (Mat_IS*)A->data;

1263:   if (D) { /* MatShift_IS pass D = NULL */
1264:     VecScatterBegin(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD);
1265:     VecScatterEnd(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD);
1266:   }
1267:   VecPointwiseDivide(is->y,is->y,is->counter);
1268:   MatDiagonalSet(is->A,is->y,insmode);
1269:   return 0;
1270: }

1272: static PetscErrorCode MatShift_IS(Mat A,PetscScalar a)
1273: {
1274:   Mat_IS         *is = (Mat_IS*)A->data;

1276:   VecSet(is->y,a);
1277:   MatDiagonalSet_IS(A,NULL,ADD_VALUES);
1278:   return 0;
1279: }

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

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

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

1297:   ISLocalToGlobalMappingApplyBlock(A->rmap->mapping,m,rows,rows_l);
1298:   ISLocalToGlobalMappingApplyBlock(A->cmap->mapping,n,cols,cols_l);
1299:   MatSetValuesBlockedLocal_IS(A,m,rows_l,n,cols_l,values,addv);
1300:   return 0;
1301: }

1303: static PetscErrorCode MatCreateSubMatrix_IS(Mat mat,IS irow,IS icol,MatReuse scall,Mat *newmat)
1304: {
1305:   Mat               locmat,newlocmat;
1306:   Mat_IS            *newmatis;
1307:   const PetscInt    *idxs;
1308:   PetscInt          i,m,n;

1310:   if (scall == MAT_REUSE_MATRIX) {
1311:     PetscBool ismatis;

1313:     PetscObjectTypeCompare((PetscObject)*newmat,MATIS,&ismatis);
1315:     newmatis = (Mat_IS*)(*newmat)->data;
1318:   }
1319:   /* irow and icol may not have duplicate entries */
1320:   if (PetscDefined(USE_DEBUG)) {
1321:     Vec               rtest,ltest;
1322:     const PetscScalar *array;

1324:     MatCreateVecs(mat,&ltest,&rtest);
1325:     ISGetLocalSize(irow,&n);
1326:     ISGetIndices(irow,&idxs);
1327:     for (i=0;i<n;i++) {
1328:       VecSetValue(rtest,idxs[i],1.0,ADD_VALUES);
1329:     }
1330:     VecAssemblyBegin(rtest);
1331:     VecAssemblyEnd(rtest);
1332:     VecGetLocalSize(rtest,&n);
1333:     VecGetOwnershipRange(rtest,&m,NULL);
1334:     VecGetArrayRead(rtest,&array);
1336:     VecRestoreArrayRead(rtest,&array);
1337:     ISRestoreIndices(irow,&idxs);
1338:     ISGetLocalSize(icol,&n);
1339:     ISGetIndices(icol,&idxs);
1340:     for (i=0;i<n;i++) {
1341:       VecSetValue(ltest,idxs[i],1.0,ADD_VALUES);
1342:     }
1343:     VecAssemblyBegin(ltest);
1344:     VecAssemblyEnd(ltest);
1345:     VecGetLocalSize(ltest,&n);
1346:     VecGetOwnershipRange(ltest,&m,NULL);
1347:     VecGetArrayRead(ltest,&array);
1349:     VecRestoreArrayRead(ltest,&array);
1350:     ISRestoreIndices(icol,&idxs);
1351:     VecDestroy(&rtest);
1352:     VecDestroy(&ltest);
1353:   }
1354:   if (scall == MAT_INITIAL_MATRIX) {
1355:     Mat_IS                 *matis = (Mat_IS*)mat->data;
1356:     ISLocalToGlobalMapping rl2g;
1357:     IS                     is;
1358:     PetscInt               *lidxs,*lgidxs,*newgidxs;
1359:     PetscInt               ll,newloc,irbs,icbs,arbs,acbs,rbs,cbs;
1360:     PetscBool              cong;
1361:     MPI_Comm               comm;

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

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

1451: static PetscErrorCode MatCopy_IS(Mat A,Mat B,MatStructure str)
1452: {
1453:   Mat_IS         *a = (Mat_IS*)A->data,*b;
1454:   PetscBool      ismatis;

1456:   PetscObjectTypeCompare((PetscObject)B,MATIS,&ismatis);
1458:   b = (Mat_IS*)B->data;
1459:   MatCopy(a->A,b->A,str);
1460:   PetscObjectStateIncrease((PetscObject)B);
1461:   return 0;
1462: }

1464: static PetscErrorCode MatMissingDiagonal_IS(Mat A,PetscBool  *missing,PetscInt *d)
1465: {
1466:   Vec               v;
1467:   const PetscScalar *array;
1468:   PetscInt          i,n;

1470:   *missing = PETSC_FALSE;
1471:   MatCreateVecs(A,NULL,&v);
1472:   MatGetDiagonal(A,v);
1473:   VecGetLocalSize(v,&n);
1474:   VecGetArrayRead(v,&array);
1475:   for (i=0;i<n;i++) if (array[i] == 0.) break;
1476:   VecRestoreArrayRead(v,&array);
1477:   VecDestroy(&v);
1478:   if (i != n) *missing = PETSC_TRUE;
1479:   if (d) {
1480:     *d = -1;
1481:     if (*missing) {
1482:       PetscInt rstart;
1483:       MatGetOwnershipRange(A,&rstart,NULL);
1484:       *d = i+rstart;
1485:     }
1486:   }
1487:   return 0;
1488: }

1490: static PetscErrorCode MatISSetUpSF_IS(Mat B)
1491: {
1492:   Mat_IS         *matis = (Mat_IS*)(B->data);
1493:   const PetscInt *gidxs;
1494:   PetscInt       nleaves;

1496:   if (matis->sf) return 0;
1497:   PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->sf);
1498:   ISLocalToGlobalMappingGetIndices(matis->rmapping,&gidxs);
1499:   ISLocalToGlobalMappingGetSize(matis->rmapping,&nleaves);
1500:   PetscSFSetGraphLayout(matis->sf,B->rmap,nleaves,NULL,PETSC_OWN_POINTER,gidxs);
1501:   ISLocalToGlobalMappingRestoreIndices(matis->rmapping,&gidxs);
1502:   PetscMalloc2(matis->sf->nroots,&matis->sf_rootdata,matis->sf->nleaves,&matis->sf_leafdata);
1503:   if (matis->rmapping != matis->cmapping) { /* setup SF for columns */
1504:     ISLocalToGlobalMappingGetSize(matis->cmapping,&nleaves);
1505:     PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->csf);
1506:     ISLocalToGlobalMappingGetIndices(matis->cmapping,&gidxs);
1507:     PetscSFSetGraphLayout(matis->csf,B->cmap,nleaves,NULL,PETSC_OWN_POINTER,gidxs);
1508:     ISLocalToGlobalMappingRestoreIndices(matis->cmapping,&gidxs);
1509:     PetscMalloc2(matis->csf->nroots,&matis->csf_rootdata,matis->csf->nleaves,&matis->csf_leafdata);
1510:   } else {
1511:     matis->csf = matis->sf;
1512:     matis->csf_leafdata = matis->sf_leafdata;
1513:     matis->csf_rootdata = matis->sf_rootdata;
1514:   }
1515:   return 0;
1516: }

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

1521:    Collective

1523:    Input Parameters:
1524: +  A - the matrix
1525: -  store - the boolean flag

1527:    Level: advanced

1529:    Notes:

1531: .seealso: MatCreate(), MatCreateIS(), MatISSetPreallocation(), MatPtAP()
1532: @*/
1533: PetscErrorCode MatISStoreL2L(Mat A, PetscBool store)
1534: {
1538:   PetscTryMethod(A,"MatISStoreL2L_C",(Mat,PetscBool),(A,store));
1539:   return 0;
1540: }

1542: static PetscErrorCode MatISStoreL2L_IS(Mat A, PetscBool store)
1543: {
1544:   Mat_IS         *matis = (Mat_IS*)(A->data);

1546:   matis->storel2l = store;
1547:   if (!store) {
1548:     PetscObjectCompose((PetscObject)(A),"_MatIS_PtAP_l2l",NULL);
1549:   }
1550:   return 0;
1551: }

1553: /*@
1554:    MatISFixLocalEmpty - Compress out zero local rows from the local matrices

1556:    Collective

1558:    Input Parameters:
1559: +  A - the matrix
1560: -  fix - the boolean flag

1562:    Level: advanced

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

1566: .seealso: MatCreate(), MatCreateIS(), MatISSetPreallocation(), MatAssemblyEnd(), MAT_FINAL_ASSEMBLY
1567: @*/
1568: PetscErrorCode MatISFixLocalEmpty(Mat A, PetscBool fix)
1569: {
1573:   PetscTryMethod(A,"MatISFixLocalEmpty_C",(Mat,PetscBool),(A,fix));
1574:   return 0;
1575: }

1577: static PetscErrorCode MatISFixLocalEmpty_IS(Mat A, PetscBool fix)
1578: {
1579:   Mat_IS *matis = (Mat_IS*)(A->data);

1581:   matis->locempty = fix;
1582:   return 0;
1583: }

1585: /*@
1586:    MatISSetPreallocation - Preallocates memory for a MATIS parallel matrix.

1588:    Collective

1590:    Input Parameters:
1591: +  B - the matrix
1592: .  d_nz  - number of nonzeros per row in DIAGONAL portion of local submatrix
1593:            (same value is used for all local rows)
1594: .  d_nnz - array containing the number of nonzeros in the various rows of the
1595:            DIAGONAL portion of the local submatrix (possibly different for each row)
1596:            or NULL, if d_nz is used to specify the nonzero structure.
1597:            The size of this array is equal to the number of local rows, i.e 'm'.
1598:            For matrices that will be factored, you must leave room for (and set)
1599:            the diagonal entry even if it is zero.
1600: .  o_nz  - number of nonzeros per row in the OFF-DIAGONAL portion of local
1601:            submatrix (same value is used for all local rows).
1602: -  o_nnz - array containing the number of nonzeros in the various rows of the
1603:            OFF-DIAGONAL portion of the local submatrix (possibly different for
1604:            each row) or NULL, if o_nz is used to specify the nonzero
1605:            structure. The size of this array is equal to the number
1606:            of local rows, i.e 'm'.

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

1610:    Level: intermediate

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

1617: .seealso: MatCreate(), MatCreateIS(), MatMPIAIJSetPreallocation(), MatISGetLocalMat(), MATIS
1618: @*/
1619: PetscErrorCode MatISSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
1620: {
1623:   PetscTryMethod(B,"MatISSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz));
1624:   return 0;
1625: }

1627: /* this is used by DMDA */
1628: PETSC_EXTERN PetscErrorCode MatISSetPreallocation_IS(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
1629: {
1630:   Mat_IS         *matis = (Mat_IS*)(B->data);
1631:   PetscInt       bs,i,nlocalcols;

1633:   MatSetUp(B);
1634:   if (!d_nnz) for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] = d_nz;
1635:   else for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] = d_nnz[i];

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

1640:   PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE);
1641:   MatGetSize(matis->A,NULL,&nlocalcols);
1642:   MatGetBlockSize(matis->A,&bs);
1643:   PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE);

1645:   for (i=0;i<matis->sf->nleaves;i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i],nlocalcols);
1646:   MatSeqAIJSetPreallocation(matis->A,0,matis->sf_leafdata);
1647: #if defined(PETSC_HAVE_HYPRE)
1648:   MatHYPRESetPreallocation(matis->A,0,matis->sf_leafdata,0,NULL);
1649: #endif

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

1654:     matis->sf_leafdata[i] = matis->sf_leafdata[i*bs]/bs;
1655:     for (b=1;b<bs;b++) {
1656:       matis->sf_leafdata[i] = PetscMax(matis->sf_leafdata[i],matis->sf_leafdata[i*bs+b]/bs);
1657:     }
1658:   }
1659:   MatSeqBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata);

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

1665:   /* for other matrix types */
1666:   MatSetUp(matis->A);
1667:   return 0;
1668: }

1670: PETSC_EXTERN PetscErrorCode MatISSetMPIXAIJPreallocation_Private(Mat A, Mat B, PetscBool maxreduce)
1671: {
1672:   Mat_IS          *matis = (Mat_IS*)(A->data);
1673:   PetscInt        *my_dnz,*my_onz,*dnz,*onz,*mat_ranges,*row_ownership;
1674:   const PetscInt  *global_indices_r,*global_indices_c;
1675:   PetscInt        i,j,bs,rows,cols;
1676:   PetscInt        lrows,lcols;
1677:   PetscInt        local_rows,local_cols;
1678:   PetscMPIInt     size;
1679:   PetscBool       isdense,issbaij;
1680:   PetscErrorCode  ierr;

1682:   MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
1683:   MatGetSize(A,&rows,&cols);
1684:   MatGetBlockSize(A,&bs);
1685:   MatGetSize(matis->A,&local_rows,&local_cols);
1686:   PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isdense);
1687:   PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);
1688:   ISLocalToGlobalMappingGetIndices(matis->rmapping,&global_indices_r);
1689:   if (matis->rmapping != matis->cmapping) {
1690:     ISLocalToGlobalMappingGetIndices(matis->cmapping,&global_indices_c);
1691:   } else global_indices_c = global_indices_r;

1693:   if (issbaij) MatGetRowUpperTriangular(matis->A);
1694:   /*
1695:      An SF reduce is needed to sum up properly on shared rows.
1696:      Note that generally preallocation is not exact, since it overestimates nonzeros
1697:   */
1698:   MatGetLocalSize(A,&lrows,&lcols);
1699:   MatPreallocateInitialize(PetscObjectComm((PetscObject)A),lrows,lcols,dnz,onz);
1700:   /* All processes need to compute entire row ownership */
1701:   PetscMalloc1(rows,&row_ownership);
1702:   MatGetOwnershipRanges(A,(const PetscInt**)&mat_ranges);
1703:   for (i=0;i<size;i++) {
1704:     for (j=mat_ranges[i];j<mat_ranges[i+1];j++) row_ownership[j] = i;
1705:   }
1706:   MatGetOwnershipRangesColumn(A,(const PetscInt**)&mat_ranges);

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

1787:   /* Reduce my_dnz and my_onz */
1788:   if (maxreduce) {
1789:     PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);
1790:     PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);
1791:     PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);
1792:     PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);
1793:   } else {
1794:     PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);
1795:     PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);
1796:     PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);
1797:     PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);
1798:   }
1799:   PetscFree2(my_dnz,my_onz);

1801:   /* Resize preallocation if overestimated */
1802:   for (i=0;i<lrows;i++) {
1803:     dnz[i] = PetscMin(dnz[i],lcols);
1804:     onz[i] = PetscMin(onz[i],cols-lcols);
1805:   }

1807:   /* Set preallocation */
1808:   MatSetBlockSizesFromMats(B,A,A);
1809:   MatSeqAIJSetPreallocation(B,0,dnz);
1810:   MatMPIAIJSetPreallocation(B,0,dnz,0,onz);
1811:   for (i=0;i<lrows;i+=bs) {
1812:     PetscInt b, d = dnz[i],o = onz[i];

1814:     for (b=1;b<bs;b++) {
1815:       d = PetscMax(d,dnz[i+b]);
1816:       o = PetscMax(o,onz[i+b]);
1817:     }
1818:     dnz[i/bs] = PetscMin(d/bs + d%bs,lcols/bs);
1819:     onz[i/bs] = PetscMin(o/bs + o%bs,(cols-lcols)/bs);
1820:   }
1821:   MatSeqBAIJSetPreallocation(B,bs,0,dnz);
1822:   MatMPIBAIJSetPreallocation(B,bs,0,dnz,0,onz);
1823:   MatMPISBAIJSetPreallocation(B,bs,0,dnz,0,onz);
1824:   MatPreallocateFinalize(dnz,onz);
1825:   if (issbaij) MatRestoreRowUpperTriangular(matis->A);
1826:   MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
1827:   return 0;
1828: }

1830: PETSC_INTERN PetscErrorCode MatConvert_IS_XAIJ(Mat mat, MatType mtype, MatReuse reuse, Mat *M)
1831: {
1832:   Mat_IS            *matis = (Mat_IS*)(mat->data);
1833:   Mat               local_mat,MT;
1834:   PetscInt          rbs,cbs,rows,cols,lrows,lcols;
1835:   PetscInt          local_rows,local_cols;
1836:   PetscBool         isseqdense,isseqsbaij,isseqaij,isseqbaij;
1837:   PetscMPIInt       size;
1838:   const PetscScalar *array;

1840:   MPI_Comm_size(PetscObjectComm((PetscObject)mat),&size);
1841:   if (size == 1 && mat->rmap->N == matis->A->rmap->N && mat->cmap->N == matis->A->cmap->N) {
1842:     Mat      B;
1843:     IS       irows = NULL,icols = NULL;
1844:     PetscInt rbs,cbs;

1846:     ISLocalToGlobalMappingGetBlockSize(matis->rmapping,&rbs);
1847:     ISLocalToGlobalMappingGetBlockSize(matis->cmapping,&cbs);
1848:     if (reuse != MAT_REUSE_MATRIX) { /* check if l2g maps are one-to-one */
1849:       IS             rows,cols;
1850:       const PetscInt *ridxs,*cidxs;
1851:       PetscInt       i,nw,*work;

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

1905:       MatCreateSubMatrix(B,irows,icols,MAT_INITIAL_MATRIX,&C);
1906:       MatHeaderReplace(mat,&C);
1907:     }
1908:     MatDestroy(&B);
1909:     ISDestroy(&icols);
1910:     ISDestroy(&irows);
1911:     return 0;
1912:   }
1913: general_assembly:
1914:   MatGetSize(mat,&rows,&cols);
1915:   ISLocalToGlobalMappingGetBlockSize(matis->rmapping,&rbs);
1916:   ISLocalToGlobalMappingGetBlockSize(matis->cmapping,&cbs);
1917:   MatGetLocalSize(mat,&lrows,&lcols);
1918:   MatGetSize(matis->A,&local_rows,&local_cols);
1919:   PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isseqdense);
1920:   PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQAIJ,&isseqaij);
1921:   PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQBAIJ,&isseqbaij);
1922:   PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&isseqsbaij);
1924:   if (PetscDefined (USE_DEBUG)) {
1925:     PetscBool         lb[4],bb[4];

1927:     lb[0] = isseqdense;
1928:     lb[1] = isseqaij;
1929:     lb[2] = isseqbaij;
1930:     lb[3] = isseqsbaij;
1931:     MPIU_Allreduce(lb,bb,4,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)mat));
1933:   }

1935:   if (reuse != MAT_REUSE_MATRIX) {
1936:     MatCreate(PetscObjectComm((PetscObject)mat),&MT);
1937:     MatSetSizes(MT,lrows,lcols,rows,cols);
1938:     MatSetType(MT,mtype);
1939:     MatSetBlockSizes(MT,rbs,cbs);
1940:     MatISSetMPIXAIJPreallocation_Private(mat,MT,PETSC_FALSE);
1941:   } else {
1942:     PetscInt mrbs,mcbs,mrows,mcols,mlrows,mlcols;

1944:     /* some checks */
1945:     MT   = *M;
1946:     MatGetBlockSizes(MT,&mrbs,&mcbs);
1947:     MatGetSize(MT,&mrows,&mcols);
1948:     MatGetLocalSize(MT,&mlrows,&mlcols);
1955:     MatZeroEntries(MT);
1956:   }

1958:   if (isseqsbaij || isseqbaij) {
1959:     MatConvert(matis->A,MATSEQAIJ,MAT_INITIAL_MATRIX,&local_mat);
1960:     isseqaij = PETSC_TRUE;
1961:   } else {
1962:     PetscObjectReference((PetscObject)matis->A);
1963:     local_mat = matis->A;
1964:   }

1966:   /* Set values */
1967:   MatSetLocalToGlobalMapping(MT,matis->rmapping,matis->cmapping);
1968:   if (isseqdense) { /* special case for dense local matrices */
1969:     PetscInt          i,*dummy;

1971:     PetscMalloc1(PetscMax(local_rows,local_cols),&dummy);
1972:     for (i=0;i<PetscMax(local_rows,local_cols);i++) dummy[i] = i;
1973:     MatSetOption(MT,MAT_ROW_ORIENTED,PETSC_FALSE);
1974:     MatDenseGetArrayRead(local_mat,&array);
1975:     MatSetValuesLocal(MT,local_rows,dummy,local_cols,dummy,array,ADD_VALUES);
1976:     MatDenseRestoreArrayRead(local_mat,&array);
1977:     PetscFree(dummy);
1978:   } else if (isseqaij) {
1979:     const PetscInt *blocks;
1980:     PetscInt       i,nvtxs,*xadj,*adjncy, nb;
1981:     PetscBool      done;
1982:     PetscScalar    *sarray;

1984:     MatGetRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);
1986:     MatSeqAIJGetArray(local_mat,&sarray);
1987:     MatGetVariableBlockSizes(local_mat,&nb,&blocks);
1988:     if (nb) { /* speed up assembly for special blocked matrices (used by BDDC) */
1989:       PetscInt sum;

1991:       for (i=0,sum=0;i<nb;i++) sum += blocks[i];
1992:       if (sum == nvtxs) {
1993:         PetscInt r;

1995:         for (i=0,r=0;i<nb;i++) {
1996:           PetscAssert(blocks[i] == xadj[r+1] - xadj[r],PETSC_COMM_SELF,PETSC_ERR_PLIB,"Invalid block sizes prescribed for block %" PetscInt_FMT ": expected %" PetscInt_FMT ", got %" PetscInt_FMT,i,blocks[i],xadj[r+1] - xadj[r]);
1997:           MatSetValuesLocal(MT,blocks[i],adjncy+xadj[r],blocks[i],adjncy+xadj[r],sarray+xadj[r],ADD_VALUES);
1998:           r   += blocks[i];
1999:         }
2000:       } else {
2001:         for (i=0;i<nvtxs;i++) {
2002:           MatSetValuesLocal(MT,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],sarray+xadj[i],ADD_VALUES);
2003:         }
2004:       }
2005:     } else {
2006:       for (i=0;i<nvtxs;i++) {
2007:         MatSetValuesLocal(MT,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],sarray+xadj[i],ADD_VALUES);
2008:       }
2009:     }
2010:     MatRestoreRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);
2012:     MatSeqAIJRestoreArray(local_mat,&sarray);
2013:   } else { /* very basic values insertion for all other matrix types */
2014:     PetscInt i;

2016:     for (i=0;i<local_rows;i++) {
2017:       PetscInt       j;
2018:       const PetscInt *local_indices_cols;

2020:       MatGetRow(local_mat,i,&j,&local_indices_cols,&array);
2021:       MatSetValuesLocal(MT,1,&i,j,local_indices_cols,array,ADD_VALUES);
2022:       MatRestoreRow(local_mat,i,&j,&local_indices_cols,&array);
2023:     }
2024:   }
2025:   MatAssemblyBegin(MT,MAT_FINAL_ASSEMBLY);
2026:   MatDestroy(&local_mat);
2027:   MatAssemblyEnd(MT,MAT_FINAL_ASSEMBLY);
2028:   if (isseqdense) {
2029:     MatSetOption(MT,MAT_ROW_ORIENTED,PETSC_TRUE);
2030:   }
2031:   if (reuse == MAT_INPLACE_MATRIX) {
2032:     MatHeaderReplace(mat,&MT);
2033:   } else if (reuse == MAT_INITIAL_MATRIX) {
2034:     *M = MT;
2035:   }
2036:   return 0;
2037: }

2039: /*@
2040:     MatISGetMPIXAIJ - Converts MATIS matrix into a parallel AIJ format

2042:   Input Parameters:
2043: +  mat - the matrix (should be of type MATIS)
2044: -  reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX

2046:   Output Parameter:
2047: .  newmat - the matrix in AIJ format

2049:   Level: developer

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

2054: .seealso: MATIS, MatConvert()
2055: @*/
2056: PetscErrorCode MatISGetMPIXAIJ(Mat mat, MatReuse reuse, Mat *newmat)
2057: {
2061:   if (reuse == MAT_REUSE_MATRIX) {
2065:   }
2066:   PetscUseMethod(mat,"MatISGetMPIXAIJ_C",(Mat,MatType,MatReuse,Mat*),(mat,MATAIJ,reuse,newmat));
2067:   return 0;
2068: }

2070: static PetscErrorCode MatDuplicate_IS(Mat mat,MatDuplicateOption op,Mat *newmat)
2071: {
2072:   Mat_IS         *matis = (Mat_IS*)(mat->data);
2073:   PetscInt       rbs,cbs,m,n,M,N;
2074:   Mat            B,localmat;

2076:   ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping,&rbs);
2077:   ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping,&cbs);
2078:   MatGetSize(mat,&M,&N);
2079:   MatGetLocalSize(mat,&m,&n);
2080:   MatCreate(PetscObjectComm((PetscObject)mat),&B);
2081:   MatSetSizes(B,m,n,M,N);
2082:   MatSetBlockSize(B,rbs == cbs ? rbs : 1);
2083:   MatSetType(B,MATIS);
2084:   MatISSetLocalMatType(B,matis->lmattype);
2085:   MatSetLocalToGlobalMapping(B,mat->rmap->mapping,mat->cmap->mapping);
2086:   MatDuplicate(matis->A,op,&localmat);
2087:   MatSetLocalToGlobalMapping(localmat,matis->A->rmap->mapping,matis->A->cmap->mapping);
2088:   MatISSetLocalMat(B,localmat);
2089:   MatDestroy(&localmat);
2090:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
2091:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
2092:   *newmat = B;
2093:   return 0;
2094: }

2096: static PetscErrorCode MatIsHermitian_IS(Mat A,PetscReal tol,PetscBool  *flg)
2097: {
2098:   Mat_IS         *matis = (Mat_IS*)A->data;
2099:   PetscBool      local_sym;

2101:   MatIsHermitian(matis->A,tol,&local_sym);
2102:   MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
2103:   return 0;
2104: }

2106: static PetscErrorCode MatIsSymmetric_IS(Mat A,PetscReal tol,PetscBool *flg)
2107: {
2108:   Mat_IS         *matis = (Mat_IS*)A->data;
2109:   PetscBool      local_sym;

2111:   if (matis->rmapping != matis->cmapping) {
2112:     *flg = PETSC_FALSE;
2113:     return 0;
2114:   }
2115:   MatIsSymmetric(matis->A,tol,&local_sym);
2116:   MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
2117:   return 0;
2118: }

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

2125:   if (matis->rmapping != matis->cmapping) {
2126:     *flg = PETSC_FALSE;
2127:     return 0;
2128:   }
2129:   MatIsStructurallySymmetric(matis->A,&local_sym);
2130:   MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
2131:   return 0;
2132: }

2134: static PetscErrorCode MatDestroy_IS(Mat A)
2135: {
2136:   Mat_IS         *b = (Mat_IS*)A->data;

2138:   PetscFree(b->bdiag);
2139:   PetscFree(b->lmattype);
2140:   MatDestroy(&b->A);
2141:   VecScatterDestroy(&b->cctx);
2142:   VecScatterDestroy(&b->rctx);
2143:   VecDestroy(&b->x);
2144:   VecDestroy(&b->y);
2145:   VecDestroy(&b->counter);
2146:   ISDestroy(&b->getsub_ris);
2147:   ISDestroy(&b->getsub_cis);
2148:   if (b->sf != b->csf) {
2149:     PetscSFDestroy(&b->csf);
2150:     PetscFree2(b->csf_rootdata,b->csf_leafdata);
2151:   } else b->csf = NULL;
2152:   PetscSFDestroy(&b->sf);
2153:   PetscFree2(b->sf_rootdata,b->sf_leafdata);
2154:   ISLocalToGlobalMappingDestroy(&b->rmapping);
2155:   ISLocalToGlobalMappingDestroy(&b->cmapping);
2156:   PetscFree(A->data);
2157:   PetscObjectChangeTypeName((PetscObject)A,NULL);
2158:   PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMatType_C",NULL);
2159:   PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL);
2160:   PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",NULL);
2161:   PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",NULL);
2162:   PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",NULL);
2163:   PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",NULL);
2164:   PetscObjectComposeFunction((PetscObject)A,"MatISFixLocalEmpty_C",NULL);
2165:   PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalToGlobalMapping_C",NULL);
2166:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpiaij_C",NULL);
2167:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpibaij_C",NULL);
2168:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpisbaij_C",NULL);
2169:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqaij_C",NULL);
2170:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqbaij_C",NULL);
2171:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqsbaij_C",NULL);
2172:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_aij_C",NULL);
2173:   PetscObjectComposeFunction((PetscObject)A,"MatSetPreallocationCOOLocal_C",NULL);
2174:   PetscObjectComposeFunction((PetscObject)A,"MatSetPreallocationCOO_C",NULL);
2175:   PetscObjectComposeFunction((PetscObject)A,"MatSetValuesCOO_C",NULL);
2176:   return 0;
2177: }

2179: static PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y)
2180: {
2181:   Mat_IS         *is  = (Mat_IS*)A->data;
2182:   PetscScalar    zero = 0.0;

2184:   /*  scatter the global vector x into the local work vector */
2185:   VecScatterBegin(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
2186:   VecScatterEnd(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);

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

2191:   /* scatter product back into global memory */
2192:   VecSet(y,zero);
2193:   VecScatterBegin(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
2194:   VecScatterEnd(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
2195:   return 0;
2196: }

2198: static PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
2199: {
2200:   Vec            temp_vec;

2202:   if (v3 != v2) {
2203:     MatMult(A,v1,v3);
2204:     VecAXPY(v3,1.0,v2);
2205:   } else {
2206:     VecDuplicate(v2,&temp_vec);
2207:     MatMult(A,v1,temp_vec);
2208:     VecAXPY(temp_vec,1.0,v2);
2209:     VecCopy(temp_vec,v3);
2210:     VecDestroy(&temp_vec);
2211:   }
2212:   return 0;
2213: }

2215: static PetscErrorCode MatMultTranspose_IS(Mat A,Vec y,Vec x)
2216: {
2217:   Mat_IS         *is = (Mat_IS*)A->data;

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

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

2226:   /* scatter product back into global vector */
2227:   VecSet(x,0);
2228:   VecScatterBegin(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);
2229:   VecScatterEnd(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);
2230:   return 0;
2231: }

2233: static PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
2234: {
2235:   Vec            temp_vec;

2237:   if (v3 != v2) {
2238:     MatMultTranspose(A,v1,v3);
2239:     VecAXPY(v3,1.0,v2);
2240:   } else {
2241:     VecDuplicate(v2,&temp_vec);
2242:     MatMultTranspose(A,v1,temp_vec);
2243:     VecAXPY(temp_vec,1.0,v2);
2244:     VecCopy(temp_vec,v3);
2245:     VecDestroy(&temp_vec);
2246:   }
2247:   return 0;
2248: }

2250: static PetscErrorCode MatView_IS(Mat A,PetscViewer viewer)
2251: {
2252:   Mat_IS         *a = (Mat_IS*)A->data;
2253:   PetscViewer    sviewer;
2254:   PetscBool      isascii,view = PETSC_TRUE;

2256:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
2257:   if (isascii)  {
2258:     PetscViewerFormat format;

2260:     PetscViewerGetFormat(viewer,&format);
2261:     if (format == PETSC_VIEWER_ASCII_INFO) view = PETSC_FALSE;
2262:   }
2263:   if (!view) return 0;
2264:   PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
2265:   MatView(a->A,sviewer);
2266:   PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
2267:   PetscViewerFlush(viewer);
2268:   return 0;
2269: }

2271: static PetscErrorCode MatInvertBlockDiagonal_IS(Mat mat,const PetscScalar **values)
2272: {
2273:   Mat_IS            *is = (Mat_IS*)mat->data;
2274:   MPI_Datatype      nodeType;
2275:   const PetscScalar *lv;
2276:   PetscInt          bs;

2278:   MatGetBlockSize(mat,&bs);
2279:   MatSetBlockSize(is->A,bs);
2280:   MatInvertBlockDiagonal(is->A,&lv);
2281:   if (!is->bdiag) {
2282:     PetscMalloc1(bs*mat->rmap->n,&is->bdiag);
2283:   }
2284:   MPI_Type_contiguous(bs,MPIU_SCALAR,&nodeType);
2285:   MPI_Type_commit(&nodeType);
2286:   PetscSFReduceBegin(is->sf,nodeType,lv,is->bdiag,MPI_REPLACE);
2287:   PetscSFReduceEnd(is->sf,nodeType,lv,is->bdiag,MPI_REPLACE);
2288:   MPI_Type_free(&nodeType);
2289:   if (values) *values = is->bdiag;
2290:   return 0;
2291: }

2293: static PetscErrorCode MatISSetUpScatters_Private(Mat A)
2294: {
2295:   Vec            cglobal,rglobal;
2296:   IS             from;
2297:   Mat_IS         *is = (Mat_IS*)A->data;
2298:   PetscScalar    sum;
2299:   const PetscInt *garray;
2300:   PetscInt       nr,rbs,nc,cbs;
2301:   VecType        rtype;

2303:   ISLocalToGlobalMappingGetSize(is->rmapping,&nr);
2304:   ISLocalToGlobalMappingGetBlockSize(is->rmapping,&rbs);
2305:   ISLocalToGlobalMappingGetSize(is->cmapping,&nc);
2306:   ISLocalToGlobalMappingGetBlockSize(is->cmapping,&cbs);
2307:   VecDestroy(&is->x);
2308:   VecDestroy(&is->y);
2309:   VecDestroy(&is->counter);
2310:   VecScatterDestroy(&is->rctx);
2311:   VecScatterDestroy(&is->cctx);
2312:   MatCreateVecs(is->A,&is->x,&is->y);
2313:   VecBindToCPU(is->y,PETSC_TRUE);
2314:   VecGetRootType_Private(is->y,&rtype);
2315:   PetscFree(A->defaultvectype);
2316:   PetscStrallocpy(rtype,&A->defaultvectype);
2317:   MatCreateVecs(A,&cglobal,&rglobal);
2318:   VecBindToCPU(rglobal,PETSC_TRUE);
2319:   ISLocalToGlobalMappingGetBlockIndices(is->rmapping,&garray);
2320:   ISCreateBlock(PetscObjectComm((PetscObject)A),rbs,nr/rbs,garray,PETSC_USE_POINTER,&from);
2321:   VecScatterCreate(rglobal,from,is->y,NULL,&is->rctx);
2322:   ISLocalToGlobalMappingRestoreBlockIndices(is->rmapping,&garray);
2323:   ISDestroy(&from);
2324:   if (is->rmapping != is->cmapping) {
2325:     ISLocalToGlobalMappingGetBlockIndices(is->cmapping,&garray);
2326:     ISCreateBlock(PetscObjectComm((PetscObject)A),cbs,nc/cbs,garray,PETSC_USE_POINTER,&from);
2327:     VecScatterCreate(cglobal,from,is->x,NULL,&is->cctx);
2328:     ISLocalToGlobalMappingRestoreBlockIndices(is->cmapping,&garray);
2329:     ISDestroy(&from);
2330:   } else {
2331:     PetscObjectReference((PetscObject)is->rctx);
2332:     is->cctx = is->rctx;
2333:   }
2334:   VecDestroy(&cglobal);

2336:   /* interface counter vector (local) */
2337:   VecDuplicate(is->y,&is->counter);
2338:   VecBindToCPU(is->counter,PETSC_TRUE);
2339:   VecSet(is->y,1.);
2340:   VecScatterBegin(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);
2341:   VecScatterEnd(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);
2342:   VecScatterBegin(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);
2343:   VecScatterEnd(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);
2344:   VecBindToCPU(is->y,PETSC_FALSE);
2345:   VecBindToCPU(is->counter,PETSC_FALSE);

2347:   /* special functions for block-diagonal matrices */
2348:   VecSum(rglobal,&sum);
2349:   A->ops->invertblockdiagonal = NULL;
2350:   if ((PetscInt)(PetscRealPart(sum)) == A->rmap->N && A->rmap->N == A->cmap->N && is->rmapping == is->cmapping) A->ops->invertblockdiagonal = MatInvertBlockDiagonal_IS;
2351:   VecDestroy(&rglobal);

2353:   /* setup SF for general purpose shared indices based communications */
2354:   MatISSetUpSF_IS(A);
2355:   return 0;
2356: }

2358: static PetscErrorCode MatISFilterL2GMap(Mat A, ISLocalToGlobalMapping map, ISLocalToGlobalMapping *nmap, ISLocalToGlobalMapping *lmap)
2359: {
2360:   IS                         is;
2361:   ISLocalToGlobalMappingType l2gtype;
2362:   const PetscInt             *idxs;
2363:   PetscHSetI                 ht;
2364:   PetscInt                   *nidxs;
2365:   PetscInt                   i,n,bs,c;
2366:   PetscBool                  flg[] = {PETSC_FALSE,PETSC_FALSE};

2368:   ISLocalToGlobalMappingGetSize(map,&n);
2369:   ISLocalToGlobalMappingGetBlockSize(map,&bs);
2370:   ISLocalToGlobalMappingGetBlockIndices(map,&idxs);
2371:   PetscHSetICreate(&ht);
2372:   PetscMalloc1(n/bs,&nidxs);
2373:   for (i=0,c=0;i<n/bs;i++) {
2374:     PetscBool missing;
2375:     if (idxs[i] < 0) { flg[0] = PETSC_TRUE; continue; }
2376:     PetscHSetIQueryAdd(ht,idxs[i],&missing);
2377:     if (!missing) flg[1] = PETSC_TRUE;
2378:     else nidxs[c++] = idxs[i];
2379:   }
2380:   PetscHSetIDestroy(&ht);
2381:   MPIU_Allreduce(MPI_IN_PLACE,flg,2,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A));
2382:   if (!flg[0] && !flg[1]) { /* Entries are all non negative and unique */
2383:     *nmap = NULL;
2384:     *lmap = NULL;
2385:     PetscFree(nidxs);
2386:     ISLocalToGlobalMappingRestoreBlockIndices(map,&idxs);
2387:     return 0;
2388:   }

2390:   /* New l2g map without negative or repeated indices */
2391:   ISCreateBlock(PetscObjectComm((PetscObject)A),bs,c,nidxs,PETSC_USE_POINTER,&is);
2392:   ISLocalToGlobalMappingCreateIS(is,nmap);
2393:   ISDestroy(&is);
2394:   ISLocalToGlobalMappingGetType(map,&l2gtype);
2395:   ISLocalToGlobalMappingSetType(*nmap,l2gtype);

2397:   /* New local l2g map for repeated indices */
2398:   ISGlobalToLocalMappingApplyBlock(*nmap,IS_GTOLM_MASK,n/bs,idxs,NULL,nidxs);
2399:   ISCreateBlock(PETSC_COMM_SELF,bs,n/bs,nidxs,PETSC_USE_POINTER,&is);
2400:   ISLocalToGlobalMappingCreateIS(is,lmap);
2401:   ISDestroy(&is);

2403:   PetscFree(nidxs);
2404:   ISLocalToGlobalMappingRestoreBlockIndices(map,&idxs);
2405:   return 0;
2406: }

2408: static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping)
2409: {
2410:   Mat_IS                 *is = (Mat_IS*)A->data;
2411:   ISLocalToGlobalMapping localrmapping = NULL, localcmapping = NULL;
2412:   PetscBool              cong, freem[] = { PETSC_FALSE, PETSC_FALSE };
2413:   PetscInt               nr,rbs,nc,cbs;


2418:   ISLocalToGlobalMappingDestroy(&is->rmapping);
2419:   ISLocalToGlobalMappingDestroy(&is->cmapping);
2420:   PetscLayoutSetUp(A->rmap);
2421:   PetscLayoutSetUp(A->cmap);
2422:   MatHasCongruentLayouts(A,&cong);

2424:   /* If NULL, local space matches global space */
2425:   if (!rmapping) {
2426:     IS is;

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

2446:     ISCreateStride(PetscObjectComm((PetscObject)A),A->cmap->N,0,1,&is);
2447:     ISLocalToGlobalMappingCreateIS(is,&cmapping);
2448:     if (A->cmap->bs > 0) ISLocalToGlobalMappingSetBlockSize(cmapping,A->cmap->bs);
2449:     ISDestroy(&is);
2450:     freem[1] = PETSC_TRUE;
2451:   } else if (cmapping != rmapping && !is->islocalref) { /* check if the l2g map has negative or repeated entries */
2452:     MatISFilterL2GMap(A,cmapping,&is->cmapping,&localcmapping);
2453:   }
2454:   if (!is->rmapping) {
2455:     PetscObjectReference((PetscObject)rmapping);
2456:     is->rmapping = rmapping;
2457:   }
2458:   if (!is->cmapping) {
2459:     PetscObjectReference((PetscObject)cmapping);
2460:     is->cmapping = cmapping;
2461:   }

2463:   /* Clean up */
2464:   MatDestroy(&is->A);
2465:   if (is->csf != is->sf) {
2466:     PetscSFDestroy(&is->csf);
2467:     PetscFree2(is->csf_rootdata,is->csf_leafdata);
2468:   } else is->csf = NULL;
2469:   PetscSFDestroy(&is->sf);
2470:   PetscFree2(is->sf_rootdata,is->sf_leafdata);
2471:   PetscFree(is->bdiag);

2473:   /* check if the two mappings are actually the same for square matrices since MATIS has some optimization for this case
2474:      (DOLFIN passes 2 different objects) */
2475:   ISLocalToGlobalMappingGetSize(is->rmapping,&nr);
2476:   ISLocalToGlobalMappingGetBlockSize(is->rmapping,&rbs);
2477:   ISLocalToGlobalMappingGetSize(is->cmapping,&nc);
2478:   ISLocalToGlobalMappingGetBlockSize(is->cmapping,&cbs);
2479:   if (is->rmapping != is->cmapping && cong) {
2480:     PetscBool same = PETSC_FALSE;
2481:     if (nr == nc && cbs == rbs) {
2482:       const PetscInt *idxs1,*idxs2;

2484:       ISLocalToGlobalMappingGetBlockIndices(is->rmapping,&idxs1);
2485:       ISLocalToGlobalMappingGetBlockIndices(is->cmapping,&idxs2);
2486:       PetscArraycmp(idxs1,idxs2,nr/rbs,&same);
2487:       ISLocalToGlobalMappingRestoreBlockIndices(is->rmapping,&idxs1);
2488:       ISLocalToGlobalMappingRestoreBlockIndices(is->cmapping,&idxs2);
2489:     }
2490:     MPIU_Allreduce(MPI_IN_PLACE,&same,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
2491:     if (same) {
2492:       ISLocalToGlobalMappingDestroy(&is->cmapping);
2493:       PetscObjectReference((PetscObject)is->rmapping);
2494:       is->cmapping = is->rmapping;
2495:     }
2496:   }
2497:   PetscLayoutSetBlockSize(A->rmap,rbs);
2498:   PetscLayoutSetBlockSize(A->cmap,cbs);
2499:   /* Pass the user defined maps to the layout */
2500:   PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping);
2501:   PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping);
2502:   if (freem[0]) ISLocalToGlobalMappingDestroy(&rmapping);
2503:   if (freem[1]) ISLocalToGlobalMappingDestroy(&cmapping);

2505:   /* Create the local matrix A */
2506:   MatCreate(PETSC_COMM_SELF,&is->A);
2507:   MatSetType(is->A,is->lmattype);
2508:   MatSetSizes(is->A,nr,nc,nr,nc);
2509:   MatSetBlockSizes(is->A,rbs,cbs);
2510:   MatSetOptionsPrefix(is->A,"is_");
2511:   MatAppendOptionsPrefix(is->A,((PetscObject)A)->prefix);
2512:   PetscLayoutSetUp(is->A->rmap);
2513:   PetscLayoutSetUp(is->A->cmap);
2514:   MatSetLocalToGlobalMapping(is->A,localrmapping,localcmapping);
2515:   ISLocalToGlobalMappingDestroy(&localrmapping);
2516:   ISLocalToGlobalMappingDestroy(&localcmapping);

2518:   /* setup scatters and local vectors for MatMult */
2519:   if (!is->islocalref) MatISSetUpScatters_Private(A);
2520:   A->preallocated = PETSC_TRUE;
2521:   return 0;
2522: }

2524: static PetscErrorCode MatSetUp_IS(Mat A)
2525: {
2526:   ISLocalToGlobalMapping rmap, cmap;

2528:   MatGetLocalToGlobalMapping(A,&rmap,&cmap);
2529:   if (!rmap && !cmap) {
2530:     MatSetLocalToGlobalMapping(A,NULL,NULL);
2531:   }
2532:   return 0;
2533: }

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

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

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

2555:   ISGlobalToLocalMappingApplyBlock(is->rmapping,IS_GTOLM_MASK,m,rows,&m,rows_l);
2556:   if (m != n || rows != cols || is->cmapping != is->rmapping) {
2557:     ISGlobalToLocalMappingApply(is->cmapping,IS_GTOLM_MASK,n,cols,&n,cols_l);
2558:     MatSetValuesBlocked(is->A,m,rows_l,n,cols_l,values,addv);
2559:   } else {
2560:     MatSetValuesBlocked(is->A,m,rows_l,n,rows_l,values,addv);
2561:   }
2562:   return 0;
2563: }

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

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

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

2581:   if (is->A->rmap->mapping || is->A->cmap->mapping) {
2582:     MatSetValuesBlockedLocal(is->A,m,rows,n,cols,values,addv);
2583:   } else {
2584:     MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv);
2585:   }
2586:   return 0;
2587: }

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

2593:   if (!n) {
2594:     is->pure_neumann = PETSC_TRUE;
2595:   } else {
2596:     PetscInt i;
2597:     is->pure_neumann = PETSC_FALSE;

2599:     if (columns) {
2600:       MatZeroRowsColumns(is->A,n,rows,diag,NULL,NULL);
2601:     } else {
2602:       MatZeroRows(is->A,n,rows,diag,NULL,NULL);
2603:     }
2604:     if (diag != 0.) {
2605:       const PetscScalar *array;
2606:       VecGetArrayRead(is->counter,&array);
2607:       for (i=0; i<n; i++) {
2608:         MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);
2609:       }
2610:       VecRestoreArrayRead(is->counter,&array);
2611:     }
2612:     MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);
2613:     MatAssemblyEnd(is->A,MAT_FINAL_ASSEMBLY);
2614:   }
2615:   return 0;
2616: }

2618: static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b,PetscBool columns)
2619: {
2620:   Mat_IS         *matis = (Mat_IS*)A->data;
2621:   PetscInt       nr,nl,len,i;
2622:   PetscInt       *lrows;

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

2627:     PetscLayoutCompare(A->rmap,A->cmap,&cong);
2628:     cong = (PetscBool)(cong && matis->sf == matis->csf);
2632:   }
2633:   /* get locally owned rows */
2634:   PetscLayoutMapLocal(A->rmap,n,rows,&len,&lrows,NULL);
2635:   /* fix right hand side if needed */
2636:   if (x && b) {
2637:     const PetscScalar *xx;
2638:     PetscScalar       *bb;

2640:     VecGetArrayRead(x, &xx);
2641:     VecGetArray(b, &bb);
2642:     for (i=0;i<len;++i) bb[lrows[i]] = diag*xx[lrows[i]];
2643:     VecRestoreArrayRead(x, &xx);
2644:     VecRestoreArray(b, &bb);
2645:   }
2646:   /* get rows associated to the local matrices */
2647:   MatGetSize(matis->A,&nl,NULL);
2648:   PetscArrayzero(matis->sf_leafdata,nl);
2649:   PetscArrayzero(matis->sf_rootdata,A->rmap->n);
2650:   for (i=0;i<len;i++) matis->sf_rootdata[lrows[i]] = 1;
2651:   PetscFree(lrows);
2652:   PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE);
2653:   PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE);
2654:   PetscMalloc1(nl,&lrows);
2655:   for (i=0,nr=0;i<nl;i++) if (matis->sf_leafdata[i]) lrows[nr++] = i;
2656:   MatISZeroRowsColumnsLocal_Private(A,nr,lrows,diag,columns);
2657:   PetscFree(lrows);
2658:   return 0;
2659: }

2661: static PetscErrorCode MatZeroRows_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
2662: {
2663:   MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_FALSE);
2664:   return 0;
2665: }

2667: static PetscErrorCode MatZeroRowsColumns_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
2668: {
2669:   MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_TRUE);
2670:   return 0;
2671: }

2673: static PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type)
2674: {
2675:   Mat_IS         *is = (Mat_IS*)A->data;

2677:   MatAssemblyBegin(is->A,type);
2678:   return 0;
2679: }

2681: static PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type)
2682: {
2683:   Mat_IS         *is = (Mat_IS*)A->data;

2685:   MatAssemblyEnd(is->A,type);
2686:   /* fix for local empty rows/cols */
2687:   if (is->locempty && type == MAT_FINAL_ASSEMBLY) {
2688:     Mat                    newlA;
2689:     ISLocalToGlobalMapping rl2g,cl2g;
2690:     IS                     nzr,nzc;
2691:     PetscInt               nr,nc,nnzr,nnzc;
2692:     PetscBool              lnewl2g,newl2g;

2694:     MatGetSize(is->A,&nr,&nc);
2695:     MatFindNonzeroRowsOrCols_Basic(is->A,PETSC_FALSE,PETSC_SMALL,&nzr);
2696:     if (!nzr) {
2697:       ISCreateStride(PetscObjectComm((PetscObject)is->A),nr,0,1,&nzr);
2698:     }
2699:     MatFindNonzeroRowsOrCols_Basic(is->A,PETSC_TRUE,PETSC_SMALL,&nzc);
2700:     if (!nzc) {
2701:       ISCreateStride(PetscObjectComm((PetscObject)is->A),nc,0,1,&nzc);
2702:     }
2703:     ISGetSize(nzr,&nnzr);
2704:     ISGetSize(nzc,&nnzc);
2705:     if (nnzr != nr || nnzc != nc) { /* need new global l2g map */
2706:       lnewl2g = PETSC_TRUE;
2707:       MPI_Allreduce(&lnewl2g,&newl2g,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A));

2709:       /* extract valid submatrix */
2710:       MatCreateSubMatrix(is->A,nzr,nzc,MAT_INITIAL_MATRIX,&newlA);
2711:     } else { /* local matrix fully populated */
2712:       lnewl2g = PETSC_FALSE;
2713:       MPI_Allreduce(&lnewl2g,&newl2g,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A));
2714:       PetscObjectReference((PetscObject)is->A);
2715:       newlA   = is->A;
2716:     }

2718:     /* attach new global l2g map if needed */
2719:     if (newl2g) {
2720:       IS              zr,zc;
2721:       const  PetscInt *ridxs,*cidxs,*zridxs,*zcidxs;
2722:       PetscInt        *nidxs,i;

2724:       ISComplement(nzr,0,nr,&zr);
2725:       ISComplement(nzc,0,nc,&zc);
2726:       PetscMalloc1(PetscMax(nr,nc),&nidxs);
2727:       ISLocalToGlobalMappingGetIndices(is->rmapping,&ridxs);
2728:       ISLocalToGlobalMappingGetIndices(is->cmapping,&cidxs);
2729:       ISGetIndices(zr,&zridxs);
2730:       ISGetIndices(zc,&zcidxs);
2731:       ISGetLocalSize(zr,&nnzr);
2732:       ISGetLocalSize(zc,&nnzc);

2734:       PetscArraycpy(nidxs,ridxs,nr);
2735:       for (i = 0; i < nnzr; i++) nidxs[zridxs[i]] = -1;
2736:       ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,nr,nidxs,PETSC_COPY_VALUES,&rl2g);
2737:       PetscArraycpy(nidxs,cidxs,nc);
2738:       for (i = 0; i < nnzc; i++) nidxs[zcidxs[i]] = -1;
2739:       ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,nc,nidxs,PETSC_COPY_VALUES,&cl2g);

2741:       ISRestoreIndices(zr,&zridxs);
2742:       ISRestoreIndices(zc,&zcidxs);
2743:       ISLocalToGlobalMappingRestoreIndices(is->rmapping,&ridxs);
2744:       ISLocalToGlobalMappingRestoreIndices(is->cmapping,&cidxs);
2745:       ISDestroy(&nzr);
2746:       ISDestroy(&nzc);
2747:       ISDestroy(&zr);
2748:       ISDestroy(&zc);
2749:       PetscFree(nidxs);
2750:       MatSetLocalToGlobalMapping(A,rl2g,cl2g);
2751:       ISLocalToGlobalMappingDestroy(&rl2g);
2752:       ISLocalToGlobalMappingDestroy(&cl2g);
2753:     }
2754:     MatISSetLocalMat(A,newlA);
2755:     MatDestroy(&newlA);
2756:     ISDestroy(&nzr);
2757:     ISDestroy(&nzc);
2758:     is->locempty = PETSC_FALSE;
2759:   }
2760:   return 0;
2761: }

2763: static PetscErrorCode MatISGetLocalMat_IS(Mat mat,Mat *local)
2764: {
2765:   Mat_IS *is = (Mat_IS*)mat->data;

2767:   *local = is->A;
2768:   return 0;
2769: }

2771: static PetscErrorCode MatISRestoreLocalMat_IS(Mat mat,Mat *local)
2772: {
2773:   *local = NULL;
2774:   return 0;
2775: }

2777: /*@
2778:     MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix.

2780:   Input Parameter:
2781: .  mat - the matrix

2783:   Output Parameter:
2784: .  local - the local matrix

2786:   Level: advanced

2788:   Notes:
2789:     This can be called if you have precomputed the nonzero structure of the
2790:   matrix and want to provide it to the inner matrix object to improve the performance
2791:   of the MatSetValues() operation.

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

2795: .seealso: MATIS
2796: @*/
2797: PetscErrorCode MatISGetLocalMat(Mat mat,Mat *local)
2798: {
2801:   PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat*),(mat,local));
2802:   return 0;
2803: }

2805: /*@
2806:     MatISRestoreLocalMat - Restores the local matrix obtained with MatISGetLocalMat()

2808:   Input Parameter:
2809: .  mat - the matrix

2811:   Output Parameter:
2812: .  local - the local matrix

2814:   Level: advanced

2816: .seealso: MATIS
2817: @*/
2818: PetscErrorCode MatISRestoreLocalMat(Mat mat,Mat *local)
2819: {
2822:   PetscUseMethod(mat,"MatISRestoreLocalMat_C",(Mat,Mat*),(mat,local));
2823:   return 0;
2824: }

2826: static PetscErrorCode MatISSetLocalMatType_IS(Mat mat,MatType mtype)
2827: {
2828:   Mat_IS         *is = (Mat_IS*)mat->data;

2830:   if (is->A) {
2831:     MatSetType(is->A,mtype);
2832:   }
2833:   PetscFree(is->lmattype);
2834:   PetscStrallocpy(mtype,&is->lmattype);
2835:   return 0;
2836: }

2838: /*@
2839:     MatISSetLocalMatType - Specifies the type of local matrix

2841:   Input Parameters:
2842: +  mat - the matrix
2843: -  mtype - the local matrix type

2845:   Output Parameter:

2847:   Level: advanced

2849: .seealso: MATIS, MatSetType(), MatType
2850: @*/
2851: PetscErrorCode MatISSetLocalMatType(Mat mat,MatType mtype)
2852: {
2854:   PetscUseMethod(mat,"MatISSetLocalMatType_C",(Mat,MatType),(mat,mtype));
2855:   return 0;
2856: }

2858: static PetscErrorCode MatISSetLocalMat_IS(Mat mat,Mat local)
2859: {
2860:   Mat_IS         *is = (Mat_IS*)mat->data;
2861:   PetscInt       nrows,ncols,orows,ocols;
2862:   MatType        mtype,otype;
2863:   PetscBool      sametype = PETSC_TRUE;

2865:   if (is->A && !is->islocalref) {
2866:     MatGetSize(is->A,&orows,&ocols);
2867:     MatGetSize(local,&nrows,&ncols);
2869:     MatGetType(local,&mtype);
2870:     MatGetType(is->A,&otype);
2871:     PetscStrcmp(mtype,otype,&sametype);
2872:   }
2873:   PetscObjectReference((PetscObject)local);
2874:   MatDestroy(&is->A);
2875:   is->A = local;
2876:   MatGetType(is->A,&mtype);
2877:   MatISSetLocalMatType(mat,mtype);
2878:   if (!sametype && !is->islocalref) {
2879:     MatISSetUpScatters_Private(mat);
2880:   }
2881:   return 0;
2882: }

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

2887:   Collective on Mat

2889:   Input Parameters:
2890: +  mat - the matrix
2891: -  local - the local matrix

2893:   Output Parameter:

2895:   Level: advanced

2897:   Notes:
2898:     This can be called if you have precomputed the local matrix and
2899:   want to provide it to the matrix object MATIS.

2901: .seealso: MATIS
2902: @*/
2903: PetscErrorCode MatISSetLocalMat(Mat mat,Mat local)
2904: {
2907:   PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local));
2908:   return 0;
2909: }

2911: static PetscErrorCode MatZeroEntries_IS(Mat A)
2912: {
2913:   Mat_IS         *a = (Mat_IS*)A->data;

2915:   MatZeroEntries(a->A);
2916:   return 0;
2917: }

2919: static PetscErrorCode MatScale_IS(Mat A,PetscScalar a)
2920: {
2921:   Mat_IS         *is = (Mat_IS*)A->data;

2923:   MatScale(is->A,a);
2924:   return 0;
2925: }

2927: static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v)
2928: {
2929:   Mat_IS         *is = (Mat_IS*)A->data;

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

2934:   /* scatter diagonal back into global vector */
2935:   VecSet(v,0);
2936:   VecScatterBegin(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);
2937:   VecScatterEnd(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);
2938:   return 0;
2939: }

2941: static PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg)
2942: {
2943:   Mat_IS         *a = (Mat_IS*)A->data;

2945:   MatSetOption(a->A,op,flg);
2946:   return 0;
2947: }

2949: static PetscErrorCode MatAXPY_IS(Mat Y,PetscScalar a,Mat X,MatStructure str)
2950: {
2951:   Mat_IS         *y = (Mat_IS*)Y->data;
2952:   Mat_IS         *x;

2954:   if (PetscDefined(USE_DEBUG)) {
2955:     PetscBool      ismatis;
2956:     PetscObjectTypeCompare((PetscObject)X,MATIS,&ismatis);
2958:   }
2959:   x = (Mat_IS*)X->data;
2960:   MatAXPY(y->A,a,x->A,str);
2961:   return 0;
2962: }

2964: static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A,IS row,IS col,Mat *submat)
2965: {
2966:   Mat                    lA;
2967:   Mat_IS                 *matis = (Mat_IS*)(A->data);
2968:   ISLocalToGlobalMapping rl2g,cl2g;
2969:   IS                     is;
2970:   const PetscInt         *rg,*rl;
2971:   PetscInt               nrg;
2972:   PetscInt               N,M,nrl,i,*idxs;

2974:   ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&rg);
2975:   ISGetLocalSize(row,&nrl);
2976:   ISGetIndices(row,&rl);
2977:   ISLocalToGlobalMappingGetSize(A->rmap->mapping,&nrg);
2978:   if (PetscDefined(USE_DEBUG)) {
2980:   }
2981:   PetscMalloc1(nrg,&idxs);
2982:   /* map from [0,nrl) to row */
2983:   for (i=0;i<nrl;i++) idxs[i] = rl[i];
2984:   for (i=nrl;i<nrg;i++) idxs[i] = -1;
2985:   ISRestoreIndices(row,&rl);
2986:   ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&rg);
2987:   ISCreateGeneral(PetscObjectComm((PetscObject)A),nrg,idxs,PETSC_OWN_POINTER,&is);
2988:   ISLocalToGlobalMappingCreateIS(is,&rl2g);
2989:   ISDestroy(&is);
2990:   /* compute new l2g map for columns */
2991:   if (col != row || matis->rmapping != matis->cmapping || matis->A->rmap->mapping != matis->A->cmap->mapping) {
2992:     const PetscInt *cg,*cl;
2993:     PetscInt       ncg;
2994:     PetscInt       ncl;

2996:     ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&cg);
2997:     ISGetLocalSize(col,&ncl);
2998:     ISGetIndices(col,&cl);
2999:     ISLocalToGlobalMappingGetSize(A->cmap->mapping,&ncg);
3000:     if (PetscDefined(USE_DEBUG)) {
3002:     }
3003:     PetscMalloc1(ncg,&idxs);
3004:     /* map from [0,ncl) to col */
3005:     for (i=0;i<ncl;i++) idxs[i] = cl[i];
3006:     for (i=ncl;i<ncg;i++) idxs[i] = -1;
3007:     ISRestoreIndices(col,&cl);
3008:     ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&cg);
3009:     ISCreateGeneral(PetscObjectComm((PetscObject)A),ncg,idxs,PETSC_OWN_POINTER,&is);
3010:     ISLocalToGlobalMappingCreateIS(is,&cl2g);
3011:     ISDestroy(&is);
3012:   } else {
3013:     PetscObjectReference((PetscObject)rl2g);
3014:     cl2g = rl2g;
3015:   }
3016:   /* create the MATIS submatrix */
3017:   MatGetSize(A,&M,&N);
3018:   MatCreate(PetscObjectComm((PetscObject)A),submat);
3019:   MatSetSizes(*submat,PETSC_DECIDE,PETSC_DECIDE,M,N);
3020:   MatSetType(*submat,MATIS);
3021:   matis = (Mat_IS*)((*submat)->data);
3022:   matis->islocalref = PETSC_TRUE;
3023:   MatSetLocalToGlobalMapping(*submat,rl2g,cl2g);
3024:   MatISGetLocalMat(A,&lA);
3025:   MatISSetLocalMat(*submat,lA);
3026:   MatSetUp(*submat);
3027:   MatAssemblyBegin(*submat,MAT_FINAL_ASSEMBLY);
3028:   MatAssemblyEnd(*submat,MAT_FINAL_ASSEMBLY);
3029:   ISLocalToGlobalMappingDestroy(&rl2g);
3030:   ISLocalToGlobalMappingDestroy(&cl2g);

3032:   /* remove unsupported ops */
3033:   PetscMemzero((*submat)->ops,sizeof(struct _MatOps));
3034:   (*submat)->ops->destroy               = MatDestroy_IS;
3035:   (*submat)->ops->setvalueslocal        = MatSetValuesLocal_SubMat_IS;
3036:   (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS;
3037:   (*submat)->ops->assemblybegin         = MatAssemblyBegin_IS;
3038:   (*submat)->ops->assemblyend           = MatAssemblyEnd_IS;
3039:   return 0;
3040: }

3042: static PetscErrorCode MatSetFromOptions_IS(PetscOptionItems *PetscOptionsObject, Mat A)
3043: {
3044:   Mat_IS         *a = (Mat_IS*)A->data;
3045:   char           type[256];
3046:   PetscBool      flg;

3048:   PetscOptionsHead(PetscOptionsObject,"MATIS options");
3049:   PetscOptionsBool("-matis_fixempty","Fix local matrices in case of empty local rows/columns","MatISFixLocalEmpty",a->locempty,&a->locempty,NULL);
3050:   PetscOptionsBool("-matis_storel2l","Store local-to-local matrices generated from PtAP operations","MatISStoreL2L",a->storel2l,&a->storel2l,NULL);
3051:   PetscOptionsFList("-matis_localmat_type","Matrix type","MatISSetLocalMatType",MatList,a->lmattype,type,256,&flg);
3052:   if (flg) {
3053:     MatISSetLocalMatType(A,type);
3054:   }
3055:   if (a->A) {
3056:     MatSetFromOptions(a->A);
3057:   }
3058:   PetscOptionsTail();
3059:   return 0;
3060: }

3062: /*@
3063:     MatCreateIS - Creates a "process" unassembled matrix, assembled on each
3064:        process but not across processes.

3066:    Input Parameters:
3067: +     comm    - MPI communicator that will share the matrix
3068: .     bs      - block size of the matrix
3069: .     m,n,M,N - local and/or global sizes of the left and right vector used in matrix vector products
3070: .     rmap    - local to global map for rows
3071: -     cmap    - local to global map for cols

3073:    Output Parameter:
3074: .    A - the resulting matrix

3076:    Level: advanced

3078:    Notes:
3079:     See MATIS for more details.
3080:     m and n are NOT related to the size of the map; they represent the size of the local parts of the distributed vectors
3081:     used in MatMult operations. The sizes of rmap and cmap define the size of the local matrices.
3082:     If rmap (cmap) is NULL, then the local row (column) spaces matches the global space.

3084: .seealso: MATIS, MatSetLocalToGlobalMapping()
3085: @*/
3086: PetscErrorCode MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping rmap,ISLocalToGlobalMapping cmap,Mat *A)
3087: {
3088:   MatCreate(comm,A);
3089:   MatSetSizes(*A,m,n,M,N);
3090:   if (bs > 0) {
3091:     MatSetBlockSize(*A,bs);
3092:   }
3093:   MatSetType(*A,MATIS);
3094:   MatSetLocalToGlobalMapping(*A,rmap,cmap);
3095:   return 0;
3096: }

3098: static PetscErrorCode MatHasOperation_IS(Mat A, MatOperation op, PetscBool *has)
3099: {
3100:   Mat_IS         *a = (Mat_IS*)A->data;

3102:   *has = PETSC_FALSE;
3103:   if (!((void**)A->ops)[op]) return 0;
3104:   MatHasOperation(a->A,op,has);
3105:   return 0;
3106: }

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

3112:   MatSetValuesCOO(a->A,v,imode);
3113:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
3114:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
3115:   return 0;
3116: }

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

3123:   if (a->A->rmap->mapping || a->A->cmap->mapping) {
3124:     MatSetPreallocationCOOLocal(a->A,ncoo,coo_i,coo_j);
3125:   } else {
3126:     MatSetPreallocationCOO(a->A,ncoo,coo_i,coo_j);
3127:   }
3128:   PetscObjectComposeFunction((PetscObject)A,"MatSetValuesCOO_C",MatSetValuesCOO_IS);
3129:   A->preallocated = PETSC_TRUE;
3130:   return 0;
3131: }

3133: static PetscErrorCode MatSetPreallocationCOO_IS(Mat A,PetscCount ncoo,const PetscInt coo_i[],const PetscInt coo_j[])
3134: {
3135:   Mat_IS         *a = (Mat_IS*)A->data;
3136:   PetscInt       *coo_il, *coo_jl, incoo;

3140:   PetscMalloc2(ncoo,&coo_il,ncoo,&coo_jl);
3141:   ISGlobalToLocalMappingApply(a->rmapping,IS_GTOLM_MASK,ncoo,coo_i,&incoo,coo_il);
3142:   ISGlobalToLocalMappingApply(a->cmapping,IS_GTOLM_MASK,ncoo,coo_j,&incoo,coo_jl);
3143:   MatSetPreallocationCOO(a->A,ncoo,coo_il,coo_jl);
3144:   PetscFree2(coo_il,coo_jl);
3145:   PetscObjectComposeFunction((PetscObject)A,"MatSetValuesCOO_C",MatSetValuesCOO_IS);
3146:   A->preallocated = PETSC_TRUE;
3147:   return 0;
3148: }

3150: /*@
3151:    MatISGetLocalToGlobalMapping - Gets the local-to-global numbering of the MATIS object

3153:    Not Collective

3155:    Input Parameter:
3156: .  A - the matrix

3158:    Output Parameters:
3159: +  rmapping - row mapping
3160: -  cmapping - column mapping

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

3164:    Level: advanced

3166: .seealso:  MatSetLocalToGlobalMapping()
3167: @*/
3168: PetscErrorCode MatISGetLocalToGlobalMapping(Mat A,ISLocalToGlobalMapping *rmapping,ISLocalToGlobalMapping *cmapping)
3169: {
3174:   PetscUseMethod(A,"MatISGetLocalToGlobalMapping_C",(Mat,ISLocalToGlobalMapping*,ISLocalToGlobalMapping*),(A,rmapping,cmapping));
3175:   return 0;
3176: }

3178: static PetscErrorCode MatISGetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping *r, ISLocalToGlobalMapping *c)
3179: {
3180:   Mat_IS *a = (Mat_IS*)A->data;

3182:   if (r) *r = a->rmapping;
3183:   if (c) *c = a->cmapping;
3184:   return 0;
3185: }

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

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

3197:    Notes:
3198:     Options prefix for the inner matrix are given by -is_mat_xxx

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

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

3205:   Level: advanced

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

3209: M*/
3210: PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A)
3211: {
3212:   Mat_IS         *a;

3214:   PetscNewLog(A,&a);
3215:   PetscStrallocpy(MATAIJ,&a->lmattype);
3216:   A->data = (void*)a;

3218:   /* matrix ops */
3219:   PetscMemzero(A->ops,sizeof(struct _MatOps));
3220:   A->ops->mult                    = MatMult_IS;
3221:   A->ops->multadd                 = MatMultAdd_IS;
3222:   A->ops->multtranspose           = MatMultTranspose_IS;
3223:   A->ops->multtransposeadd        = MatMultTransposeAdd_IS;
3224:   A->ops->destroy                 = MatDestroy_IS;
3225:   A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS;
3226:   A->ops->setvalues               = MatSetValues_IS;
3227:   A->ops->setvaluesblocked        = MatSetValuesBlocked_IS;
3228:   A->ops->setvalueslocal          = MatSetValuesLocal_IS;
3229:   A->ops->setvaluesblockedlocal   = MatSetValuesBlockedLocal_IS;
3230:   A->ops->zerorows                = MatZeroRows_IS;
3231:   A->ops->zerorowscolumns         = MatZeroRowsColumns_IS;
3232:   A->ops->assemblybegin           = MatAssemblyBegin_IS;
3233:   A->ops->assemblyend             = MatAssemblyEnd_IS;
3234:   A->ops->view                    = MatView_IS;
3235:   A->ops->zeroentries             = MatZeroEntries_IS;
3236:   A->ops->scale                   = MatScale_IS;
3237:   A->ops->getdiagonal             = MatGetDiagonal_IS;
3238:   A->ops->setoption               = MatSetOption_IS;
3239:   A->ops->ishermitian             = MatIsHermitian_IS;
3240:   A->ops->issymmetric             = MatIsSymmetric_IS;
3241:   A->ops->isstructurallysymmetric = MatIsStructurallySymmetric_IS;
3242:   A->ops->duplicate               = MatDuplicate_IS;
3243:   A->ops->missingdiagonal         = MatMissingDiagonal_IS;
3244:   A->ops->copy                    = MatCopy_IS;
3245:   A->ops->getlocalsubmatrix       = MatGetLocalSubMatrix_IS;
3246:   A->ops->createsubmatrix         = MatCreateSubMatrix_IS;
3247:   A->ops->axpy                    = MatAXPY_IS;
3248:   A->ops->diagonalset             = MatDiagonalSet_IS;
3249:   A->ops->shift                   = MatShift_IS;
3250:   A->ops->transpose               = MatTranspose_IS;
3251:   A->ops->getinfo                 = MatGetInfo_IS;
3252:   A->ops->diagonalscale           = MatDiagonalScale_IS;
3253:   A->ops->setfromoptions          = MatSetFromOptions_IS;
3254:   A->ops->setup                   = MatSetUp_IS;
3255:   A->ops->hasoperation            = MatHasOperation_IS;

3257:   /* special MATIS functions */
3258:   PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMatType_C",MatISSetLocalMatType_IS);
3259:   PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS);
3260:   PetscObjectComposeFunction((PetscObject)A,"MatISRestoreLocalMat_C",MatISRestoreLocalMat_IS);
3261:   PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS);
3262:   PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",MatConvert_IS_XAIJ);
3263:   PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",MatISSetPreallocation_IS);
3264:   PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",MatISStoreL2L_IS);
3265:   PetscObjectComposeFunction((PetscObject)A,"MatISFixLocalEmpty_C",MatISFixLocalEmpty_IS);
3266:   PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalToGlobalMapping_C",MatISGetLocalToGlobalMapping_IS);
3267:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpiaij_C",MatConvert_IS_XAIJ);
3268:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpibaij_C",MatConvert_IS_XAIJ);
3269:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpisbaij_C",MatConvert_IS_XAIJ);
3270:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqaij_C",MatConvert_IS_XAIJ);
3271:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqbaij_C",MatConvert_IS_XAIJ);
3272:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqsbaij_C",MatConvert_IS_XAIJ);
3273:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_aij_C",MatConvert_IS_XAIJ);
3274:   PetscObjectComposeFunction((PetscObject)A,"MatSetPreallocationCOOLocal_C",MatSetPreallocationCOOLocal_IS);
3275:   PetscObjectComposeFunction((PetscObject)A,"MatSetPreallocationCOO_C",MatSetPreallocationCOO_IS);
3276:   PetscObjectChangeTypeName((PetscObject)A,MATIS);
3277:   return 0;
3278: }