Actual source code: matis.c

petsc-3.10.5 2019-03-28
Report Typos and Errors
  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>

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

 18: static PetscErrorCode MatISContainerDestroyPtAP_Private(void *ptr)
 19: {
 20:   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;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

252: PETSC_INTERN PetscErrorCode MatPtAP_IS_XAIJ(Mat A,Mat P,MatReuse scall,PetscReal fill,Mat *C)
253: {

257:   if (scall == MAT_INITIAL_MATRIX) {
258:     PetscLogEventBegin(MAT_PtAPSymbolic,A,P,0,0);
259:     MatPtAPSymbolic_IS_XAIJ(A,P,fill,C);
260:     PetscLogEventEnd(MAT_PtAPSymbolic,A,P,0,0);
261:   }
262:   PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);
263:   ((*C)->ops->ptapnumeric)(A,P,*C);
264:   PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);
265:   return(0);
266: }

268: static PetscErrorCode MatISContainerDestroyFields_Private(void *ptr)
269: {
270:   MatISLocalFields lf = (MatISLocalFields)ptr;
271:   PetscInt         i;
272:   PetscErrorCode   ierr;

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

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

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

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

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

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

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

359:   MatGetRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&m,&ii,&jj,&flg);
360:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get IJ structure");
361:   ISLocalToGlobalMappingGetNodeInfo(A->rmap->mapping,&n,&ecount,&eneighs);
362:   if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected %D != %D",m,n);
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(A->rmap->mapping,&n,&ecount,&eneighs);
382:   MatSeqAIJRestoreArray(matis->A,&aa);
383:   MatRestoreRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&m,&ii,&jj,&flg);
384:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot restore IJ structure");
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_",0};
399:   MatISDisassemblel2gType mode = MAT_IS_DISASSEMBLE_L2G_NATURAL;
400:   MatPartitioning         part;
401:   PetscSF                 sf;
402:   PetscErrorCode          ierr;

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

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

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

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

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

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

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

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

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

491:           ISGetLocalSize(workis[j],&s);
492:           ISGetIndices(workis[j],&idxs);
493:           MatPreallocateSet(i+A->rmap->rstart,s,idxs,dnz,onz);
494:           j++;
495:         }
496:       }
497:       if (j != cnt) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected local count %D != %D",j,cnt);

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

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

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

517:         if (s) {
518:           ISGetIndices(workis[j],&idxs);
519:           MatSetValues(A2,1,&row,s,idxs,vals,INSERT_VALUES);
520:           ISRestoreIndices(workis[j],&idxs);
521:           j++;
522:         }
523:       }
524:       if (j != cnt) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected local count %D != %D",j,cnt);
525:       PetscFree(vals);
526:       MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY);
527:       MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY);
528:       MatTranspose(A2,MAT_INPLACE_MATRIX,&A2);

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

538:       /* extend local to global map to include connected isolated separators */
539:       PetscObjectQuery((PetscObject)A3,"_petsc_GetLocalMatCondensed_iscol",(PetscObject*)&is);
540:       if (!is) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing column map");
541:       ISLocalToGlobalMappingCreateIS(is,&ll2g);
542:       MatGetRowIJ(A3,0,PETSC_FALSE,PETSC_FALSE,&i,&ii,&jj,&flg);
543:       if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get IJ structure");
544:       ISCreateGeneral(PETSC_COMM_SELF,ii[i],jj,PETSC_COPY_VALUES,&is);
545:       MatRestoreRowIJ(A3,0,PETSC_FALSE,PETSC_FALSE,&i,&ii,&jj,&flg);
546:       if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get IJ structure");
547:       ISLocalToGlobalMappingApplyIS(ll2g,is,&is2);
548:       ISDestroy(&is);
549:       ISLocalToGlobalMappingDestroy(&ll2g);

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

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

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

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

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

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

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

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

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

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

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

762:   MatRestoreRowIJ(Ad,0,PETSC_FALSE,PETSC_FALSE,&i,&odi,&odj,&flg);
763:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot restore IJ structure");
764:   MatRestoreRowIJ(Ao,0,PETSC_FALSE,PETSC_FALSE,&i,&ooi,&ooj,&flg);
765:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot restore IJ structure");
766:   MatSeqAIJRestoreArray(Ad,&dd);
767:   MatSeqAIJRestoreArray(Ao,&od);

769:   ii   = aux;
770:   jj   = aux+dr+1;
771:   aa   = data;
772:   MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,dr,lc,ii,jj,aa,&lA);

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

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

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

802: /*@
803:    MatISSetUpSF - Setup star forest objects used by MatIS.

805:    Collective on Mat

807:    Input Parameters:
808: +  A - the matrix

810:    Level: advanced

812:    Notes:
813:     This function does not need to be called by the user.

815: .keywords: matrix

817: .seealso: MatCreate(), MatCreateIS(), MatISSetPreallocation(), MatISGetLocalMat()
818: @*/
819: PetscErrorCode  MatISSetUpSF(Mat A)
820: {

826:   PetscTryMethod(A,"MatISSetUpSF_C",(Mat),(A));
827:   return(0);
828: }

830: PETSC_INTERN PetscErrorCode MatConvert_Nest_IS(Mat A,MatType type,MatReuse reuse,Mat *newmat)
831: {
832:   Mat                    **nest,*snest,**rnest,lA,B;
833:   IS                     *iscol,*isrow,*islrow,*islcol;
834:   ISLocalToGlobalMapping rl2g,cl2g;
835:   MPI_Comm               comm;
836:   PetscInt               *lr,*lc,*l2gidxs;
837:   PetscInt               i,j,nr,nc,rbs,cbs;
838:   PetscBool              convert,lreuse,*istrans;
839:   PetscErrorCode         ierr;

842:   MatNestGetSubMats(A,&nr,&nc,&nest);
843:   lreuse = PETSC_FALSE;
844:   rnest  = NULL;
845:   if (reuse == MAT_REUSE_MATRIX) {
846:     PetscBool ismatis,isnest;

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

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

872:       /* Nested matrices should be of type MATIS */
873:       PetscObjectTypeCompare((PetscObject)nest[i][j],MATTRANSPOSEMAT,&istrans[ij]);
874:       if (istrans[ij]) {
875:         Mat T,lT;
876:         MatTransposeGetMat(nest[i][j],&T);
877:         PetscObjectTypeCompare((PetscObject)T,MATIS,&ismatis);
878:         if (!ismatis) SETERRQ2(comm,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) (transposed) is not of type MATIS",i,j);
879:         MatISGetLocalMat(T,&lT);
880:         MatCreateTranspose(lT,&snest[ij]);
881:       } else {
882:         PetscObjectTypeCompare((PetscObject)nest[i][j],MATIS,&ismatis);
883:         if (!ismatis) SETERRQ2(comm,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) is not of type MATIS",i,j);
884:         MatISGetLocalMat(nest[i][j],&snest[ij]);
885:       }

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

896:       /* check compatibilty for local matrix reusage */
897:       if (rnest && !rnest[i][j] != !snest[ij]) lreuse = PETSC_FALSE;
898:     }
899:   }

901: #if defined (PETSC_USE_DEBUG)
902:   /* Check compatibility of l2g maps for rows */
903:   for (i=0;i<nr;i++) {
904:     rl2g = NULL;
905:     for (j=0;j<nc;j++) {
906:       PetscInt n1,n2;

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

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

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

942:       if (!nest[j][i]) continue;
943:       if (istrans[j*nc+i]) {
944:         Mat T;

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

959:         ISLocalToGlobalMappingGetSize(rl2g,&n2);
960:         if (n1 != n2) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid column l2gmap size %D != %D",j,i,n1,n2);
961:         ISLocalToGlobalMappingGetIndices(cl2g,&idxs1);
962:         ISLocalToGlobalMappingGetIndices(rl2g,&idxs2);
963:         PetscMemcmp(idxs1,idxs2,n1*sizeof(PetscInt),&same);
964:         ISLocalToGlobalMappingRestoreIndices(cl2g,&idxs1);
965:         ISLocalToGlobalMappingRestoreIndices(rl2g,&idxs2);
966:         if (!same) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid column l2gmap",j,i);
967:       }
968:     }
969:   }
970: #endif

972:   B = NULL;
973:   if (reuse != MAT_REUSE_MATRIX) {
974:     PetscInt stl;

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

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

987:       /* l2gmap */
988:       j = 0;
989:       usedmat = nest[i][j];
990:       while (!usedmat && j < nc-1) usedmat = nest[i][++j];
991:       if (!usedmat) SETERRQ(comm,PETSC_ERR_SUP,"Cannot find valid row mat");

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

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

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

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

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

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

1121:   /* Create local matrix in MATNEST format */
1122:   convert = PETSC_FALSE;
1123:   PetscOptionsGetBool(NULL,((PetscObject)A)->prefix,"-matis_convert_local_nest",&convert,NULL);
1124:   if (convert) {
1125:     Mat              M;
1126:     MatISLocalFields lf;
1127:     PetscContainer   c;

1129:     MatISGetLocalMat(*newmat,&lA);
1130:     MatConvert(lA,MATAIJ,MAT_INITIAL_MATRIX,&M);
1131:     MatISSetLocalMat(*newmat,M);
1132:     MatDestroy(&M);

1134:     /* attach local fields to the matrix */
1135:     PetscNew(&lf);
1136:     PetscCalloc2(nr,&lf->rf,nc,&lf->cf);
1137:     for (i=0;i<nr;i++) {
1138:       PetscInt n,st;

1140:       ISGetLocalSize(islrow[i],&n);
1141:       ISStrideGetInfo(islrow[i],&st,NULL);
1142:       ISCreateStride(comm,n,st,1,&lf->rf[i]);
1143:     }
1144:     for (i=0;i<nc;i++) {
1145:       PetscInt n,st;

1147:       ISGetLocalSize(islcol[i],&n);
1148:       ISStrideGetInfo(islcol[i],&st,NULL);
1149:       ISCreateStride(comm,n,st,1,&lf->cf[i]);
1150:     }
1151:     lf->nr = nr;
1152:     lf->nc = nc;
1153:     PetscContainerCreate(PetscObjectComm((PetscObject)(*newmat)),&c);
1154:     PetscContainerSetPointer(c,lf);
1155:     PetscContainerSetUserDestroy(c,MatISContainerDestroyFields_Private);
1156:     PetscObjectCompose((PetscObject)(*newmat),"_convert_nest_lfields",(PetscObject)c);
1157:     PetscContainerDestroy(&c);
1158:   }

1160:   /* Free workspace */
1161:   for (i=0;i<nr;i++) {
1162:     ISDestroy(&islrow[i]);
1163:   }
1164:   for (i=0;i<nc;i++) {
1165:     ISDestroy(&islcol[i]);
1166:   }
1167:   PetscFree6(isrow,iscol,islrow,islcol,snest,istrans);
1168:   PetscFree2(lr,lc);
1169:   return(0);
1170: }

1172: static PetscErrorCode MatDiagonalScale_IS(Mat A, Vec l, Vec r)
1173: {
1174:   Mat_IS            *matis = (Mat_IS*)A->data;
1175:   Vec               ll,rr;
1176:   const PetscScalar *Y,*X;
1177:   PetscScalar       *x,*y;
1178:   PetscErrorCode    ierr;

1181:   MatISSetUpSF(A);
1182:   if (l) {
1183:     ll   = matis->y;
1184:     VecGetArrayRead(l,&Y);
1185:     VecGetArray(ll,&y);
1186:     PetscSFBcastBegin(matis->sf,MPIU_SCALAR,Y,y);
1187:   } else {
1188:     ll = NULL;
1189:   }
1190:   if (r) {
1191:     rr   = matis->x;
1192:     VecGetArrayRead(r,&X);
1193:     VecGetArray(rr,&x);
1194:     PetscSFBcastBegin(matis->csf,MPIU_SCALAR,X,x);
1195:   } else {
1196:     rr = NULL;
1197:   }
1198:   if (ll) {
1199:     PetscSFBcastEnd(matis->sf,MPIU_SCALAR,Y,y);
1200:     VecRestoreArrayRead(l,&Y);
1201:     VecRestoreArray(ll,&y);
1202:   }
1203:   if (rr) {
1204:     PetscSFBcastEnd(matis->csf,MPIU_SCALAR,X,x);
1205:     VecRestoreArrayRead(r,&X);
1206:     VecRestoreArray(rr,&x);
1207:   }
1208:   MatDiagonalScale(matis->A,ll,rr);
1209:   return(0);
1210: }

1212: static PetscErrorCode MatGetInfo_IS(Mat A,MatInfoType flag,MatInfo *ginfo)
1213: {
1214:   Mat_IS         *matis = (Mat_IS*)A->data;
1215:   MatInfo        info;
1216:   PetscReal      isend[6],irecv[6];
1217:   PetscInt       bs;

1221:   MatGetBlockSize(A,&bs);
1222:   if (matis->A->ops->getinfo) {
1223:     MatGetInfo(matis->A,MAT_LOCAL,&info);
1224:     isend[0] = info.nz_used;
1225:     isend[1] = info.nz_allocated;
1226:     isend[2] = info.nz_unneeded;
1227:     isend[3] = info.memory;
1228:     isend[4] = info.mallocs;
1229:   } else {
1230:     isend[0] = 0.;
1231:     isend[1] = 0.;
1232:     isend[2] = 0.;
1233:     isend[3] = 0.;
1234:     isend[4] = 0.;
1235:   }
1236:   isend[5] = matis->A->num_ass;
1237:   if (flag == MAT_LOCAL) {
1238:     ginfo->nz_used      = isend[0];
1239:     ginfo->nz_allocated = isend[1];
1240:     ginfo->nz_unneeded  = isend[2];
1241:     ginfo->memory       = isend[3];
1242:     ginfo->mallocs      = isend[4];
1243:     ginfo->assemblies   = isend[5];
1244:   } else if (flag == MAT_GLOBAL_MAX) {
1245:     MPIU_Allreduce(isend,irecv,6,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)A));

1247:     ginfo->nz_used      = irecv[0];
1248:     ginfo->nz_allocated = irecv[1];
1249:     ginfo->nz_unneeded  = irecv[2];
1250:     ginfo->memory       = irecv[3];
1251:     ginfo->mallocs      = irecv[4];
1252:     ginfo->assemblies   = irecv[5];
1253:   } else if (flag == MAT_GLOBAL_SUM) {
1254:     MPIU_Allreduce(isend,irecv,5,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)A));

1256:     ginfo->nz_used      = irecv[0];
1257:     ginfo->nz_allocated = irecv[1];
1258:     ginfo->nz_unneeded  = irecv[2];
1259:     ginfo->memory       = irecv[3];
1260:     ginfo->mallocs      = irecv[4];
1261:     ginfo->assemblies   = A->num_ass;
1262:   }
1263:   ginfo->block_size        = bs;
1264:   ginfo->fill_ratio_given  = 0;
1265:   ginfo->fill_ratio_needed = 0;
1266:   ginfo->factor_mallocs    = 0;
1267:   return(0);
1268: }

1270: PetscErrorCode MatTranspose_IS(Mat A,MatReuse reuse,Mat *B)
1271: {
1272:   Mat                    C,lC,lA;
1273:   PetscErrorCode         ierr;

1276:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
1277:     ISLocalToGlobalMapping rl2g,cl2g;
1278:     MatCreate(PetscObjectComm((PetscObject)A),&C);
1279:     MatSetSizes(C,A->cmap->n,A->rmap->n,A->cmap->N,A->rmap->N);
1280:     MatSetBlockSizes(C,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));
1281:     MatSetType(C,MATIS);
1282:     MatGetLocalToGlobalMapping(A,&rl2g,&cl2g);
1283:     MatSetLocalToGlobalMapping(C,cl2g,rl2g);
1284:   } else {
1285:     C = *B;
1286:   }

1288:   /* perform local transposition */
1289:   MatISGetLocalMat(A,&lA);
1290:   MatTranspose(lA,MAT_INITIAL_MATRIX,&lC);
1291:   MatISSetLocalMat(C,lC);
1292:   MatDestroy(&lC);

1294:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
1295:     *B = C;
1296:   } else {
1297:     MatHeaderMerge(A,&C);
1298:   }
1299:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
1300:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
1301:   return(0);
1302: }

1304: PetscErrorCode  MatDiagonalSet_IS(Mat A,Vec D,InsertMode insmode)
1305: {
1306:   Mat_IS         *is = (Mat_IS*)A->data;

1310:   if (D) { /* MatShift_IS pass D = NULL */
1311:     VecScatterBegin(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD);
1312:     VecScatterEnd(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD);
1313:   }
1314:   VecPointwiseDivide(is->y,is->y,is->counter);
1315:   MatDiagonalSet(is->A,is->y,insmode);
1316:   return(0);
1317: }

1319: PetscErrorCode  MatShift_IS(Mat A,PetscScalar a)
1320: {
1321:   Mat_IS         *is = (Mat_IS*)A->data;

1325:   VecSet(is->y,a);
1326:   MatDiagonalSet_IS(A,NULL,ADD_VALUES);
1327:   return(0);
1328: }

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

1336: #if defined(PETSC_USE_DEBUG)
1337:   if (m > MATIS_MAX_ENTRIES_INSERTION || n > MATIS_MAX_ENTRIES_INSERTION) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of row/column indices must be <= %D: they are %D %D",MATIS_MAX_ENTRIES_INSERTION,m,n);
1338: #endif
1339:   ISLocalToGlobalMappingApply(A->rmap->mapping,m,rows,rows_l);
1340:   ISLocalToGlobalMappingApply(A->cmap->mapping,n,cols,cols_l);
1341:   MatSetValuesLocal_IS(A,m,rows_l,n,cols_l,values,addv);
1342:   return(0);
1343: }

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

1351: #if defined(PETSC_USE_DEBUG)
1352:   if (m > MATIS_MAX_ENTRIES_INSERTION || n > MATIS_MAX_ENTRIES_INSERTION) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of row/column block indices must be <= %D: they are %D %D",MATIS_MAX_ENTRIES_INSERTION,m,n);
1353: #endif
1354:   ISLocalToGlobalMappingApplyBlock(A->rmap->mapping,m,rows,rows_l);
1355:   ISLocalToGlobalMappingApplyBlock(A->cmap->mapping,n,cols,cols_l);
1356:   MatSetValuesBlockedLocal_IS(A,m,rows_l,n,cols_l,values,addv);
1357:   return(0);
1358: }

1360: static PetscErrorCode PetscLayoutMapLocal_Private(PetscLayout map,PetscInt N,const PetscInt idxs[], PetscInt *on,PetscInt **oidxs,PetscInt **ogidxs)
1361: {
1362:   PetscInt      *owners = map->range;
1363:   PetscInt       n      = map->n;
1364:   PetscSF        sf;
1365:   PetscInt      *lidxs,*work = NULL;
1366:   PetscSFNode   *ridxs;
1367:   PetscMPIInt    rank;
1368:   PetscInt       r, p = 0, len = 0;

1372:   if (on) *on = 0;              /* squelch -Wmaybe-uninitialized */
1373:   /* Create SF where leaves are input idxs and roots are owned idxs (code adapted from MatZeroRowsMapLocal_Private) */
1374:   MPI_Comm_rank(map->comm,&rank);
1375:   PetscMalloc1(n,&lidxs);
1376:   for (r = 0; r < n; ++r) lidxs[r] = -1;
1377:   PetscMalloc1(N,&ridxs);
1378:   for (r = 0; r < N; ++r) {
1379:     const PetscInt idx = idxs[r];
1380:     if (idx < 0 || map->N <= idx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Index %D out of range [0,%D)",idx,map->N);
1381:     if (idx < owners[p] || owners[p+1] <= idx) { /* short-circuit the search if the last p owns this idx too */
1382:       PetscLayoutFindOwner(map,idx,&p);
1383:     }
1384:     ridxs[r].rank = p;
1385:     ridxs[r].index = idxs[r] - owners[p];
1386:   }
1387:   PetscSFCreate(map->comm,&sf);
1388:   PetscSFSetGraph(sf,n,N,NULL,PETSC_OWN_POINTER,ridxs,PETSC_OWN_POINTER);
1389:   PetscSFReduceBegin(sf,MPIU_INT,(PetscInt*)idxs,lidxs,MPI_LOR);
1390:   PetscSFReduceEnd(sf,MPIU_INT,(PetscInt*)idxs,lidxs,MPI_LOR);
1391:   if (ogidxs) { /* communicate global idxs */
1392:     PetscInt cum = 0,start,*work2;

1394:     PetscMalloc1(n,&work);
1395:     PetscCalloc1(N,&work2);
1396:     for (r = 0; r < N; ++r) if (idxs[r] >=0) cum++;
1397:     MPI_Scan(&cum,&start,1,MPIU_INT,MPI_SUM,map->comm);
1398:     start -= cum;
1399:     cum = 0;
1400:     for (r = 0; r < N; ++r) if (idxs[r] >=0) work2[r] = start+cum++;
1401:     PetscSFReduceBegin(sf,MPIU_INT,work2,work,MPIU_REPLACE);
1402:     PetscSFReduceEnd(sf,MPIU_INT,work2,work,MPIU_REPLACE);
1403:     PetscFree(work2);
1404:   }
1405:   PetscSFDestroy(&sf);
1406:   /* Compress and put in indices */
1407:   for (r = 0; r < n; ++r)
1408:     if (lidxs[r] >= 0) {
1409:       if (work) work[len] = work[r];
1410:       lidxs[len++] = r;
1411:     }
1412:   if (on) *on = len;
1413:   if (oidxs) *oidxs = lidxs;
1414:   if (ogidxs) *ogidxs = work;
1415:   return(0);
1416: }

1418: static PetscErrorCode MatCreateSubMatrix_IS(Mat mat,IS irow,IS icol,MatReuse scall,Mat *newmat)
1419: {
1420:   Mat               locmat,newlocmat;
1421:   Mat_IS            *newmatis;
1422: #if defined(PETSC_USE_DEBUG)
1423:   Vec               rtest,ltest;
1424:   const PetscScalar *array;
1425: #endif
1426:   const PetscInt    *idxs;
1427:   PetscInt          i,m,n;
1428:   PetscErrorCode    ierr;

1431:   if (scall == MAT_REUSE_MATRIX) {
1432:     PetscBool ismatis;

1434:     PetscObjectTypeCompare((PetscObject)*newmat,MATIS,&ismatis);
1435:     if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Not of MATIS type");
1436:     newmatis = (Mat_IS*)(*newmat)->data;
1437:     if (!newmatis->getsub_ris) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Misses local row IS");
1438:     if (!newmatis->getsub_cis) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Misses local col IS");
1439:   }
1440:   /* irow and icol may not have duplicate entries */
1441: #if defined(PETSC_USE_DEBUG)
1442:   MatCreateVecs(mat,&ltest,&rtest);
1443:   ISGetLocalSize(irow,&n);
1444:   ISGetIndices(irow,&idxs);
1445:   for (i=0;i<n;i++) {
1446:     VecSetValue(rtest,idxs[i],1.0,ADD_VALUES);
1447:   }
1448:   VecAssemblyBegin(rtest);
1449:   VecAssemblyEnd(rtest);
1450:   VecGetLocalSize(rtest,&n);
1451:   VecGetOwnershipRange(rtest,&m,NULL);
1452:   VecGetArrayRead(rtest,&array);
1453:   for (i=0;i<n;i++) if (array[i] != 0. && array[i] != 1.) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Index %D counted %D times! Irow may not have duplicate entries",i+m,(PetscInt)PetscRealPart(array[i]));
1454:   VecRestoreArrayRead(rtest,&array);
1455:   ISRestoreIndices(irow,&idxs);
1456:   ISGetLocalSize(icol,&n);
1457:   ISGetIndices(icol,&idxs);
1458:   for (i=0;i<n;i++) {
1459:     VecSetValue(ltest,idxs[i],1.0,ADD_VALUES);
1460:   }
1461:   VecAssemblyBegin(ltest);
1462:   VecAssemblyEnd(ltest);
1463:   VecGetLocalSize(ltest,&n);
1464:   VecGetOwnershipRange(ltest,&m,NULL);
1465:   VecGetArrayRead(ltest,&array);
1466:   for (i=0;i<n;i++) if (array[i] != 0. && array[i] != 1.) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Index %D counted %D times! Icol may not have duplicate entries",i+m,(PetscInt)PetscRealPart(array[i]));
1467:   VecRestoreArrayRead(ltest,&array);
1468:   ISRestoreIndices(icol,&idxs);
1469:   VecDestroy(&rtest);
1470:   VecDestroy(&ltest);
1471: #endif
1472:   if (scall == MAT_INITIAL_MATRIX) {
1473:     Mat_IS                 *matis = (Mat_IS*)mat->data;
1474:     ISLocalToGlobalMapping rl2g;
1475:     IS                     is;
1476:     PetscInt               *lidxs,*lgidxs,*newgidxs;
1477:     PetscInt               ll,newloc,irbs,icbs,arbs,acbs,rbs,cbs;
1478:     PetscBool              cong;
1479:     MPI_Comm               comm;

1481:     PetscObjectGetComm((PetscObject)mat,&comm);
1482:     MatGetBlockSizes(mat,&arbs,&acbs);
1483:     ISGetBlockSize(irow,&irbs);
1484:     ISGetBlockSize(icol,&icbs);
1485:     rbs  = arbs == irbs ? irbs : 1;
1486:     cbs  = acbs == icbs ? icbs : 1;
1487:     ISGetLocalSize(irow,&m);
1488:     ISGetLocalSize(icol,&n);
1489:     MatCreate(comm,newmat);
1490:     MatSetType(*newmat,MATIS);
1491:     MatSetSizes(*newmat,m,n,PETSC_DECIDE,PETSC_DECIDE);
1492:     MatSetBlockSizes(*newmat,rbs,cbs);
1493:     /* communicate irow to their owners in the layout */
1494:     ISGetIndices(irow,&idxs);
1495:     PetscLayoutMapLocal_Private(mat->rmap,m,idxs,&ll,&lidxs,&lgidxs);
1496:     ISRestoreIndices(irow,&idxs);
1497:     MatISSetUpSF(mat);
1498:     PetscMemzero(matis->sf_rootdata,matis->sf->nroots*sizeof(PetscInt));
1499:     for (i=0;i<ll;i++) matis->sf_rootdata[lidxs[i]] = lgidxs[i]+1;
1500:     PetscFree(lidxs);
1501:     PetscFree(lgidxs);
1502:     PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
1503:     PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
1504:     for (i=0,newloc=0;i<matis->sf->nleaves;i++) if (matis->sf_leafdata[i]) newloc++;
1505:     PetscMalloc1(newloc,&newgidxs);
1506:     PetscMalloc1(newloc,&lidxs);
1507:     for (i=0,newloc=0;i<matis->sf->nleaves;i++)
1508:       if (matis->sf_leafdata[i]) {
1509:         lidxs[newloc] = i;
1510:         newgidxs[newloc++] = matis->sf_leafdata[i]-1;
1511:       }
1512:     ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is);
1513:     ISLocalToGlobalMappingCreateIS(is,&rl2g);
1514:     ISLocalToGlobalMappingSetBlockSize(rl2g,rbs);
1515:     ISDestroy(&is);
1516:     /* local is to extract local submatrix */
1517:     newmatis = (Mat_IS*)(*newmat)->data;
1518:     ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_ris);
1519:     MatHasCongruentLayouts(mat,&cong);
1520:     if (cong && irow == icol && matis->csf == matis->sf) {
1521:       MatSetLocalToGlobalMapping(*newmat,rl2g,rl2g);
1522:       PetscObjectReference((PetscObject)newmatis->getsub_ris);
1523:       newmatis->getsub_cis = newmatis->getsub_ris;
1524:     } else {
1525:       ISLocalToGlobalMapping cl2g;

1527:       /* communicate icol to their owners in the layout */
1528:       ISGetIndices(icol,&idxs);
1529:       PetscLayoutMapLocal_Private(mat->cmap,n,idxs,&ll,&lidxs,&lgidxs);
1530:       ISRestoreIndices(icol,&idxs);
1531:       PetscMemzero(matis->csf_rootdata,matis->csf->nroots*sizeof(PetscInt));
1532:       for (i=0;i<ll;i++) matis->csf_rootdata[lidxs[i]] = lgidxs[i]+1;
1533:       PetscFree(lidxs);
1534:       PetscFree(lgidxs);
1535:       PetscSFBcastBegin(matis->csf,MPIU_INT,matis->csf_rootdata,matis->csf_leafdata);
1536:       PetscSFBcastEnd(matis->csf,MPIU_INT,matis->csf_rootdata,matis->csf_leafdata);
1537:       for (i=0,newloc=0;i<matis->csf->nleaves;i++) if (matis->csf_leafdata[i]) newloc++;
1538:       PetscMalloc1(newloc,&newgidxs);
1539:       PetscMalloc1(newloc,&lidxs);
1540:       for (i=0,newloc=0;i<matis->csf->nleaves;i++)
1541:         if (matis->csf_leafdata[i]) {
1542:           lidxs[newloc] = i;
1543:           newgidxs[newloc++] = matis->csf_leafdata[i]-1;
1544:         }
1545:       ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is);
1546:       ISLocalToGlobalMappingCreateIS(is,&cl2g);
1547:       ISLocalToGlobalMappingSetBlockSize(cl2g,cbs);
1548:       ISDestroy(&is);
1549:       /* local is to extract local submatrix */
1550:       ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_cis);
1551:       MatSetLocalToGlobalMapping(*newmat,rl2g,cl2g);
1552:       ISLocalToGlobalMappingDestroy(&cl2g);
1553:     }
1554:     ISLocalToGlobalMappingDestroy(&rl2g);
1555:   } else {
1556:     MatISGetLocalMat(*newmat,&newlocmat);
1557:   }
1558:   MatISGetLocalMat(mat,&locmat);
1559:   newmatis = (Mat_IS*)(*newmat)->data;
1560:   MatCreateSubMatrix(locmat,newmatis->getsub_ris,newmatis->getsub_cis,scall,&newlocmat);
1561:   if (scall == MAT_INITIAL_MATRIX) {
1562:     MatISSetLocalMat(*newmat,newlocmat);
1563:     MatDestroy(&newlocmat);
1564:   }
1565:   MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
1566:   MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
1567:   return(0);
1568: }

1570: static PetscErrorCode MatCopy_IS(Mat A,Mat B,MatStructure str)
1571: {
1572:   Mat_IS         *a = (Mat_IS*)A->data,*b;
1573:   PetscBool      ismatis;

1577:   PetscObjectTypeCompare((PetscObject)B,MATIS,&ismatis);
1578:   if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Need to be implemented");
1579:   b = (Mat_IS*)B->data;
1580:   MatCopy(a->A,b->A,str);
1581:   PetscObjectStateIncrease((PetscObject)B);
1582:   return(0);
1583: }

1585: static PetscErrorCode MatMissingDiagonal_IS(Mat A,PetscBool  *missing,PetscInt *d)
1586: {
1587:   Vec               v;
1588:   const PetscScalar *array;
1589:   PetscInt          i,n;
1590:   PetscErrorCode    ierr;

1593:   *missing = PETSC_FALSE;
1594:   MatCreateVecs(A,NULL,&v);
1595:   MatGetDiagonal(A,v);
1596:   VecGetLocalSize(v,&n);
1597:   VecGetArrayRead(v,&array);
1598:   for (i=0;i<n;i++) if (array[i] == 0.) break;
1599:   VecRestoreArrayRead(v,&array);
1600:   VecDestroy(&v);
1601:   if (i != n) *missing = PETSC_TRUE;
1602:   if (d) {
1603:     *d = -1;
1604:     if (*missing) {
1605:       PetscInt rstart;
1606:       MatGetOwnershipRange(A,&rstart,NULL);
1607:       *d = i+rstart;
1608:     }
1609:   }
1610:   return(0);
1611: }

1613: static PetscErrorCode MatISSetUpSF_IS(Mat B)
1614: {
1615:   Mat_IS         *matis = (Mat_IS*)(B->data);
1616:   const PetscInt *gidxs;
1617:   PetscInt       nleaves;

1621:   if (matis->sf) return(0);
1622:   PetscLayoutSetUp(B->rmap);
1623:   PetscLayoutSetUp(B->cmap);
1624:   PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->sf);
1625:   ISLocalToGlobalMappingGetIndices(B->rmap->mapping,&gidxs);
1626:   ISLocalToGlobalMappingGetSize(B->rmap->mapping,&nleaves);
1627:   PetscSFSetGraphLayout(matis->sf,B->rmap,nleaves,NULL,PETSC_OWN_POINTER,gidxs);
1628:   ISLocalToGlobalMappingRestoreIndices(B->rmap->mapping,&gidxs);
1629:   PetscMalloc2(matis->sf->nroots,&matis->sf_rootdata,matis->sf->nleaves,&matis->sf_leafdata);
1630:   if (B->rmap->mapping != B->cmap->mapping) { /* setup SF for columns */
1631:     ISLocalToGlobalMappingGetSize(B->cmap->mapping,&nleaves);
1632:     PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->csf);
1633:     ISLocalToGlobalMappingGetIndices(B->cmap->mapping,&gidxs);
1634:     PetscSFSetGraphLayout(matis->csf,B->cmap,nleaves,NULL,PETSC_OWN_POINTER,gidxs);
1635:     ISLocalToGlobalMappingRestoreIndices(B->cmap->mapping,&gidxs);
1636:     PetscMalloc2(matis->csf->nroots,&matis->csf_rootdata,matis->csf->nleaves,&matis->csf_leafdata);
1637:   } else {
1638:     matis->csf = matis->sf;
1639:     matis->csf_leafdata = matis->sf_leafdata;
1640:     matis->csf_rootdata = matis->sf_rootdata;
1641:   }
1642:   return(0);
1643: }

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

1648:    Collective on MPI_Comm

1650:    Input Parameters:
1651: +  A - the matrix
1652: -  store - the boolean flag

1654:    Level: advanced

1656:    Notes:

1658: .keywords: matrix

1660: .seealso: MatCreate(), MatCreateIS(), MatISSetPreallocation(), MatPtAP()
1661: @*/
1662: PetscErrorCode  MatISStoreL2L(Mat A, PetscBool store)
1663: {

1670:   PetscTryMethod(A,"MatISStoreL2L_C",(Mat,PetscBool),(A,store));
1671:   return(0);
1672: }

1674: static PetscErrorCode MatISStoreL2L_IS(Mat A, PetscBool store)
1675: {
1676:   Mat_IS         *matis = (Mat_IS*)(A->data);

1680:   matis->storel2l = store;
1681:   if (!store) {
1682:     PetscObjectCompose((PetscObject)(A),"_MatIS_PtAP_l2l",NULL);
1683:   }
1684:   return(0);
1685: }

1687: /*@
1688:    MatISFixLocalEmpty - Compress out zero local rows from the local matrices

1690:    Collective on MPI_Comm

1692:    Input Parameters:
1693: +  A - the matrix
1694: -  fix - the boolean flag

1696:    Level: advanced

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

1700: .keywords: matrix

1702: .seealso: MatCreate(), MatCreateIS(), MatISSetPreallocation(), MatAssemblyEnd(), MAT_FINAL_ASSEMBLY
1703: @*/
1704: PetscErrorCode  MatISFixLocalEmpty(Mat A, PetscBool fix)
1705: {

1712:   PetscTryMethod(A,"MatISFixLocalEmpty_C",(Mat,PetscBool),(A,fix));
1713:   return(0);
1714: }

1716: static PetscErrorCode MatISFixLocalEmpty_IS(Mat A, PetscBool fix)
1717: {
1718:   Mat_IS *matis = (Mat_IS*)(A->data);

1721:   matis->locempty = fix;
1722:   return(0);
1723: }

1725: /*@
1726:    MatISSetPreallocation - Preallocates memory for a MATIS parallel matrix.

1728:    Collective on MPI_Comm

1730:    Input Parameters:
1731: +  B - the matrix
1732: .  d_nz  - number of nonzeros per row in DIAGONAL portion of local submatrix
1733:            (same value is used for all local rows)
1734: .  d_nnz - array containing the number of nonzeros in the various rows of the
1735:            DIAGONAL portion of the local submatrix (possibly different for each row)
1736:            or NULL, if d_nz is used to specify the nonzero structure.
1737:            The size of this array is equal to the number of local rows, i.e 'm'.
1738:            For matrices that will be factored, you must leave room for (and set)
1739:            the diagonal entry even if it is zero.
1740: .  o_nz  - number of nonzeros per row in the OFF-DIAGONAL portion of local
1741:            submatrix (same value is used for all local rows).
1742: -  o_nnz - array containing the number of nonzeros in the various rows of the
1743:            OFF-DIAGONAL portion of the local submatrix (possibly different for
1744:            each row) or NULL, if o_nz is used to specify the nonzero
1745:            structure. The size of this array is equal to the number
1746:            of local rows, i.e 'm'.

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

1750:    Level: intermediate

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

1757: .keywords: matrix

1759: .seealso: MatCreate(), MatCreateIS(), MatMPIAIJSetPreallocation(), MatISGetLocalMat(), MATIS
1760: @*/
1761: PetscErrorCode  MatISSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
1762: {

1768:   PetscTryMethod(B,"MatISSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz));
1769:   return(0);
1770: }

1772: /* this is used by DMDA */
1773: PETSC_EXTERN PetscErrorCode  MatISSetPreallocation_IS(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
1774: {
1775:   Mat_IS         *matis = (Mat_IS*)(B->data);
1776:   PetscInt       bs,i,nlocalcols;

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

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

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

1789:   PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
1790:   MatGetSize(matis->A,NULL,&nlocalcols);
1791:   MatGetBlockSize(matis->A,&bs);
1792:   PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);

1794:   for (i=0;i<matis->sf->nleaves;i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i],nlocalcols);
1795:   MatSeqAIJSetPreallocation(matis->A,0,matis->sf_leafdata);
1796: #if defined(PETSC_HAVE_HYPRE)
1797:   MatHYPRESetPreallocation(matis->A,0,matis->sf_leafdata,0,NULL);
1798: #endif

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

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

1807:   /* for other matrix types */
1808:   MatSetUp(matis->A);
1809:   return(0);
1810: }

1812: PETSC_EXTERN PetscErrorCode MatISSetMPIXAIJPreallocation_Private(Mat A, Mat B, PetscBool maxreduce)
1813: {
1814:   Mat_IS          *matis = (Mat_IS*)(A->data);
1815:   PetscInt        *my_dnz,*my_onz,*dnz,*onz,*mat_ranges,*row_ownership;
1816:   const PetscInt  *global_indices_r,*global_indices_c;
1817:   PetscInt        i,j,bs,rows,cols;
1818:   PetscInt        lrows,lcols;
1819:   PetscInt        local_rows,local_cols;
1820:   PetscMPIInt     size;
1821:   PetscBool       isdense,issbaij;
1822:   PetscErrorCode  ierr;

1825:   MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
1826:   MatGetSize(A,&rows,&cols);
1827:   MatGetBlockSize(A,&bs);
1828:   MatGetSize(matis->A,&local_rows,&local_cols);
1829:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isdense);
1830:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);
1831:   ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&global_indices_r);
1832:   if (A->rmap->mapping != A->cmap->mapping) {
1833:     ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&global_indices_c);
1834:   } else {
1835:     global_indices_c = global_indices_r;
1836:   }

1838:   if (issbaij) {
1839:     MatGetRowUpperTriangular(matis->A);
1840:   }
1841:   /*
1842:      An SF reduce is needed to sum up properly on shared rows.
1843:      Note that generally preallocation is not exact, since it overestimates nonzeros
1844:   */
1845:   MatISSetUpSF(A);
1846:   MatGetLocalSize(A,&lrows,&lcols);
1847:   MatPreallocateInitialize(PetscObjectComm((PetscObject)A),lrows,lcols,dnz,onz);
1848:   /* All processes need to compute entire row ownership */
1849:   PetscMalloc1(rows,&row_ownership);
1850:   MatGetOwnershipRanges(A,(const PetscInt**)&mat_ranges);
1851:   for (i=0;i<size;i++) {
1852:     for (j=mat_ranges[i];j<mat_ranges[i+1];j++) {
1853:       row_ownership[j] = i;
1854:     }
1855:   }
1856:   MatGetOwnershipRangesColumn(A,(const PetscInt**)&mat_ranges);

1858:   /*
1859:      my_dnz and my_onz contains exact contribution to preallocation from each local mat
1860:      then, they will be summed up properly. This way, preallocation is always sufficient
1861:   */
1862:   PetscCalloc2(local_rows,&my_dnz,local_rows,&my_onz);
1863:   /* preallocation as a MATAIJ */
1864:   if (isdense) { /* special case for dense local matrices */
1865:     for (i=0;i<local_rows;i++) {
1866:       PetscInt owner = row_ownership[global_indices_r[i]];
1867:       for (j=0;j<local_cols;j++) {
1868:         PetscInt index_col = global_indices_c[j];
1869:         if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */
1870:           my_dnz[i] += 1;
1871:         } else { /* offdiag block */
1872:           my_onz[i] += 1;
1873:         }
1874:       }
1875:     }
1876:   } else if (matis->A->ops->getrowij) {
1877:     const PetscInt *ii,*jj,*jptr;
1878:     PetscBool      done;
1879:     MatGetRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done);
1880:     if (!done) SETERRQ(PetscObjectComm((PetscObject)(matis->A)),PETSC_ERR_PLIB,"Error in MatGetRowIJ");
1881:     jptr = jj;
1882:     for (i=0;i<local_rows;i++) {
1883:       PetscInt index_row = global_indices_r[i];
1884:       for (j=0;j<ii[i+1]-ii[i];j++,jptr++) {
1885:         PetscInt owner = row_ownership[index_row];
1886:         PetscInt index_col = global_indices_c[*jptr];
1887:         if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */
1888:           my_dnz[i] += 1;
1889:         } else { /* offdiag block */
1890:           my_onz[i] += 1;
1891:         }
1892:         /* same as before, interchanging rows and cols */
1893:         if (issbaij && index_col != index_row) {
1894:           owner = row_ownership[index_col];
1895:           if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) {
1896:             my_dnz[*jptr] += 1;
1897:           } else {
1898:             my_onz[*jptr] += 1;
1899:           }
1900:         }
1901:       }
1902:     }
1903:     MatRestoreRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done);
1904:     if (!done) SETERRQ(PetscObjectComm((PetscObject)(matis->A)),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ");
1905:   } else { /* loop over rows and use MatGetRow */
1906:     for (i=0;i<local_rows;i++) {
1907:       const PetscInt *cols;
1908:       PetscInt       ncols,index_row = global_indices_r[i];
1909:       MatGetRow(matis->A,i,&ncols,&cols,NULL);
1910:       for (j=0;j<ncols;j++) {
1911:         PetscInt owner = row_ownership[index_row];
1912:         PetscInt index_col = global_indices_c[cols[j]];
1913:         if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */
1914:           my_dnz[i] += 1;
1915:         } else { /* offdiag block */
1916:           my_onz[i] += 1;
1917:         }
1918:         /* same as before, interchanging rows and cols */
1919:         if (issbaij && index_col != index_row) {
1920:           owner = row_ownership[index_col];
1921:           if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) {
1922:             my_dnz[cols[j]] += 1;
1923:           } else {
1924:             my_onz[cols[j]] += 1;
1925:           }
1926:         }
1927:       }
1928:       MatRestoreRow(matis->A,i,&ncols,&cols,NULL);
1929:     }
1930:   }
1931:   if (global_indices_c != global_indices_r) {
1932:     ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&global_indices_c);
1933:   }
1934:   ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&global_indices_r);
1935:   PetscFree(row_ownership);

1937:   /* Reduce my_dnz and my_onz */
1938:   if (maxreduce) {
1939:     PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);
1940:     PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);
1941:     PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);
1942:     PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);
1943:   } else {
1944:     PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);
1945:     PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);
1946:     PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);
1947:     PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);
1948:   }
1949:   PetscFree2(my_dnz,my_onz);

1951:   /* Resize preallocation if overestimated */
1952:   for (i=0;i<lrows;i++) {
1953:     dnz[i] = PetscMin(dnz[i],lcols);
1954:     onz[i] = PetscMin(onz[i],cols-lcols);
1955:   }

1957:   /* Set preallocation */
1958:   MatSeqAIJSetPreallocation(B,0,dnz);
1959:   MatMPIAIJSetPreallocation(B,0,dnz,0,onz);
1960:   for (i=0;i<lrows;i+=bs) {
1961:     PetscInt b, d = dnz[i],o = onz[i];

1963:     for (b=1;b<bs;b++) {
1964:       d = PetscMax(d,dnz[i+b]);
1965:       o = PetscMax(o,onz[i+b]);
1966:     }
1967:     dnz[i/bs] = PetscMin(d/bs + d%bs,lcols/bs);
1968:     onz[i/bs] = PetscMin(o/bs + o%bs,(cols-lcols)/bs);
1969:   }
1970:   MatSeqBAIJSetPreallocation(B,bs,0,dnz);
1971:   MatMPIBAIJSetPreallocation(B,bs,0,dnz,0,onz);
1972:   MatMPISBAIJSetPreallocation(B,bs,0,dnz,0,onz);
1973:   MatPreallocateFinalize(dnz,onz);
1974:   if (issbaij) {
1975:     MatRestoreRowUpperTriangular(matis->A);
1976:   }
1977:   MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
1978:   return(0);
1979: }

1981: PETSC_INTERN PetscErrorCode MatConvert_IS_XAIJ(Mat mat, MatType mtype, MatReuse reuse, Mat *M)
1982: {
1983:   Mat_IS         *matis = (Mat_IS*)(mat->data);
1984:   Mat            local_mat,MT;
1985:   PetscInt       rbs,cbs,rows,cols,lrows,lcols;
1986:   PetscInt       local_rows,local_cols;
1987:   PetscBool      isseqdense,isseqsbaij,isseqaij,isseqbaij;
1988: #if defined (PETSC_USE_DEBUG)
1989:   PetscBool      lb[4],bb[4];
1990: #endif
1991:   PetscMPIInt    size;
1992:   PetscScalar    *array;

1996:   MPI_Comm_size(PetscObjectComm((PetscObject)mat),&size);
1997:   if (size == 1 && mat->rmap->N == matis->A->rmap->N && mat->cmap->N == matis->A->cmap->N) {
1998:     Mat      B;
1999:     IS       irows = NULL,icols = NULL;
2000:     PetscInt rbs,cbs;

2002:     ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping,&rbs);
2003:     ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping,&cbs);
2004:     if (reuse != MAT_REUSE_MATRIX) { /* check if l2g maps are one-to-one */
2005:       IS             rows,cols;
2006:       const PetscInt *ridxs,*cidxs;
2007:       PetscInt       i,nw,*work;

2009:       ISLocalToGlobalMappingGetBlockIndices(mat->rmap->mapping,&ridxs);
2010:       ISLocalToGlobalMappingGetSize(mat->rmap->mapping,&nw);
2011:       nw   = nw/rbs;
2012:       PetscCalloc1(nw,&work);
2013:       for (i=0;i<nw;i++) work[ridxs[i]] += 1;
2014:       for (i=0;i<nw;i++) if (!work[i] || work[i] > 1) break;
2015:       if (i == nw) {
2016:         ISCreateBlock(PETSC_COMM_SELF,rbs,nw,ridxs,PETSC_USE_POINTER,&rows);
2017:         ISSetPermutation(rows);
2018:         ISInvertPermutation(rows,PETSC_DECIDE,&irows);
2019:         ISDestroy(&rows);
2020:       }
2021:       ISLocalToGlobalMappingRestoreBlockIndices(mat->rmap->mapping,&ridxs);
2022:       PetscFree(work);
2023:       if (irows && mat->rmap->mapping != mat->cmap->mapping) {
2024:         ISLocalToGlobalMappingGetBlockIndices(mat->cmap->mapping,&cidxs);
2025:         ISLocalToGlobalMappingGetSize(mat->cmap->mapping,&nw);
2026:         nw   = nw/cbs;
2027:         PetscCalloc1(nw,&work);
2028:         for (i=0;i<nw;i++) work[cidxs[i]] += 1;
2029:         for (i=0;i<nw;i++) if (!work[i] || work[i] > 1) break;
2030:         if (i == nw) {
2031:           ISCreateBlock(PETSC_COMM_SELF,cbs,nw,cidxs,PETSC_USE_POINTER,&cols);
2032:           ISSetPermutation(cols);
2033:           ISInvertPermutation(cols,PETSC_DECIDE,&icols);
2034:           ISDestroy(&cols);
2035:         }
2036:         ISLocalToGlobalMappingRestoreBlockIndices(mat->cmap->mapping,&cidxs);
2037:         PetscFree(work);
2038:       } else if (irows) {
2039:         PetscObjectReference((PetscObject)irows);
2040:         icols = irows;
2041:       }
2042:     } else {
2043:       PetscObjectQuery((PetscObject)(*M),"_MatIS_IS_XAIJ_irows",(PetscObject*)&irows);
2044:       PetscObjectQuery((PetscObject)(*M),"_MatIS_IS_XAIJ_icols",(PetscObject*)&icols);
2045:       if (irows) { PetscObjectReference((PetscObject)irows); }
2046:       if (icols) { PetscObjectReference((PetscObject)icols); }
2047:     }
2048:     if (!irows || !icols) {
2049:       ISDestroy(&icols);
2050:       ISDestroy(&irows);
2051:       goto general_assembly;
2052:     }
2053:     MatConvert(matis->A,(rbs == cbs && rbs > 1) ? MATSEQBAIJ : MATSEQAIJ,MAT_INITIAL_MATRIX,&B);
2054:     if (reuse != MAT_INPLACE_MATRIX) {
2055:       MatCreateSubMatrix(B,irows,icols,reuse,M);
2056:       PetscObjectCompose((PetscObject)(*M),"_MatIS_IS_XAIJ_irows",(PetscObject)irows);
2057:       PetscObjectCompose((PetscObject)(*M),"_MatIS_IS_XAIJ_icols",(PetscObject)icols);
2058:     } else {
2059:       Mat C;

2061:       MatCreateSubMatrix(B,irows,icols,MAT_INITIAL_MATRIX,&C);
2062:       MatHeaderReplace(mat,&C);
2063:     }
2064:     MatDestroy(&B);
2065:     ISDestroy(&icols);
2066:     ISDestroy(&irows);
2067:     return(0);
2068:   }
2069: general_assembly:
2070:   MatGetSize(mat,&rows,&cols);
2071:   ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping,&rbs);
2072:   ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping,&cbs);
2073:   MatGetLocalSize(mat,&lrows,&lcols);
2074:   MatGetSize(matis->A,&local_rows,&local_cols);
2075:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isseqdense);
2076:   PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQAIJ,&isseqaij);
2077:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQBAIJ,&isseqbaij);
2078:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&isseqsbaij);
2079:   if (!isseqdense && !isseqaij && !isseqbaij && !isseqsbaij) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)(matis->A))->type_name);
2080: #if defined (PETSC_USE_DEBUG)
2081:   lb[0] = isseqdense;
2082:   lb[1] = isseqaij;
2083:   lb[2] = isseqbaij;
2084:   lb[3] = isseqsbaij;
2085:   MPIU_Allreduce(lb,bb,4,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)mat));
2086:   if (!bb[0] && !bb[1] && !bb[2] && !bb[3]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Local matrices must have the same type");
2087: #endif

2089:   if (reuse != MAT_REUSE_MATRIX) {
2090:     MatCreate(PetscObjectComm((PetscObject)mat),&MT);
2091:     MatSetSizes(MT,lrows,lcols,rows,cols);
2092:     MatSetType(MT,mtype);
2093:     MatSetBlockSizes(MT,rbs,cbs);
2094:     MatISSetMPIXAIJPreallocation_Private(mat,MT,PETSC_FALSE);
2095:   } else {
2096:     PetscInt mrbs,mcbs,mrows,mcols,mlrows,mlcols;

2098:     /* some checks */
2099:     MT   = *M;
2100:     MatGetBlockSizes(MT,&mrbs,&mcbs);
2101:     MatGetSize(MT,&mrows,&mcols);
2102:     MatGetLocalSize(MT,&mlrows,&mlcols);
2103:     if (mrows != rows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of rows (%d != %d)",rows,mrows);
2104:     if (mcols != cols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of cols (%d != %d)",cols,mcols);
2105:     if (mlrows != lrows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local rows (%d != %d)",lrows,mlrows);
2106:     if (mlcols != lcols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local cols (%d != %d)",lcols,mlcols);
2107:     if (mrbs != rbs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong row block size (%d != %d)",rbs,mrbs);
2108:     if (mcbs != cbs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong col block size (%d != %d)",cbs,mcbs);
2109:     MatZeroEntries(MT);
2110:   }

2112:   if (isseqsbaij || isseqbaij) {
2113:     MatConvert(matis->A,MATSEQAIJ,MAT_INITIAL_MATRIX,&local_mat);
2114:     isseqaij = PETSC_TRUE;
2115:   } else {
2116:     PetscObjectReference((PetscObject)matis->A);
2117:     local_mat = matis->A;
2118:   }

2120:   /* Set values */
2121:   MatSetLocalToGlobalMapping(MT,mat->rmap->mapping,mat->cmap->mapping);
2122:   if (isseqdense) { /* special case for dense local matrices */
2123:     PetscInt i,*dummy;

2125:     PetscMalloc1(PetscMax(local_rows,local_cols),&dummy);
2126:     for (i=0;i<PetscMax(local_rows,local_cols);i++) dummy[i] = i;
2127:     MatSetOption(MT,MAT_ROW_ORIENTED,PETSC_FALSE);
2128:     MatDenseGetArray(local_mat,&array);
2129:     MatSetValuesLocal(MT,local_rows,dummy,local_cols,dummy,array,ADD_VALUES);
2130:     MatDenseRestoreArray(local_mat,&array);
2131:     PetscFree(dummy);
2132:   } else if (isseqaij) {
2133:     const PetscInt *blocks;
2134:     PetscInt       i,nvtxs,*xadj,*adjncy, nb;
2135:     PetscBool      done;

2137:     MatGetRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);
2138:     if (!done) SETERRQ(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatGetRowIJ");
2139:     MatSeqAIJGetArray(local_mat,&array);
2140:     MatGetVariableBlockSizes(local_mat,&nb,&blocks);
2141:     if (nb) { /* speed up assembly for special blocked matrices (used by BDDC) */
2142:       PetscInt sum;

2144:       for (i=0,sum=0;i<nb;i++) sum += blocks[i];
2145:       if (sum == nvtxs) {
2146:         PetscInt r;

2148:         for (i=0,r=0;i<nb;i++) {
2149:           if (blocks[i] != xadj[r+1] - xadj[r]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Invalid block sizes prescribed for block %D: expected %D, got %D",i,blocks[i],xadj[r+1] - xadj[r]);
2150:           MatSetValuesLocal(MT,blocks[i],adjncy+xadj[r],blocks[i],adjncy+xadj[r],array+xadj[r],ADD_VALUES);
2151:           r   += blocks[i];
2152:         }
2153:       } else {
2154:         for (i=0;i<nvtxs;i++) {
2155:           MatSetValuesLocal(MT,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],array+xadj[i],ADD_VALUES);
2156:         }
2157:       }
2158:     } else {
2159:       for (i=0;i<nvtxs;i++) {
2160:         MatSetValuesLocal(MT,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],array+xadj[i],ADD_VALUES);
2161:       }
2162:     }
2163:     MatRestoreRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);
2164:     if (!done) SETERRQ(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ");
2165:     MatSeqAIJRestoreArray(local_mat,&array);
2166:   } else { /* very basic values insertion for all other matrix types */
2167:     PetscInt i;

2169:     for (i=0;i<local_rows;i++) {
2170:       PetscInt       j;
2171:       const PetscInt *local_indices_cols;

2173:       MatGetRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);
2174:       MatSetValuesLocal(MT,1,&i,j,local_indices_cols,array,ADD_VALUES);
2175:       MatRestoreRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);
2176:     }
2177:   }
2178:   MatAssemblyBegin(MT,MAT_FINAL_ASSEMBLY);
2179:   MatDestroy(&local_mat);
2180:   MatAssemblyEnd(MT,MAT_FINAL_ASSEMBLY);
2181:   if (isseqdense) {
2182:     MatSetOption(MT,MAT_ROW_ORIENTED,PETSC_TRUE);
2183:   }
2184:   if (reuse == MAT_INPLACE_MATRIX) {
2185:     MatHeaderReplace(mat,&MT);
2186:   } else if (reuse == MAT_INITIAL_MATRIX) {
2187:     *M = MT;
2188:   }
2189:   return(0);
2190: }

2192: /*@
2193:     MatISGetMPIXAIJ - Converts MATIS matrix into a parallel AIJ format

2195:   Input Parameter:
2196: .  mat - the matrix (should be of type MATIS)
2197: .  reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX

2199:   Output Parameter:
2200: .  newmat - the matrix in AIJ format

2202:   Level: developer

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

2207: .seealso: MATIS, MatConvert()
2208: @*/
2209: PetscErrorCode MatISGetMPIXAIJ(Mat mat, MatReuse reuse, Mat *newmat)
2210: {

2217:   if (reuse == MAT_REUSE_MATRIX) {
2220:     if (mat == *newmat) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse the same matrix");
2221:   }
2222:   PetscUseMethod(mat,"MatISGetMPIXAIJ_C",(Mat,MatType,MatReuse,Mat*),(mat,MATAIJ,reuse,newmat));
2223:   return(0);
2224: }

2226: PetscErrorCode MatDuplicate_IS(Mat mat,MatDuplicateOption op,Mat *newmat)
2227: {
2229:   Mat_IS         *matis = (Mat_IS*)(mat->data);
2230:   PetscInt       rbs,cbs,m,n,M,N;
2231:   Mat            B,localmat;

2234:   ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping,&rbs);
2235:   ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping,&cbs);
2236:   MatGetSize(mat,&M,&N);
2237:   MatGetLocalSize(mat,&m,&n);
2238:   MatCreate(PetscObjectComm((PetscObject)mat),&B);
2239:   MatSetSizes(B,m,n,M,N);
2240:   MatSetBlockSize(B,rbs == cbs ? rbs : 1);
2241:   MatSetType(B,MATIS);
2242:   MatISSetLocalMatType(B,matis->lmattype);
2243:   MatSetLocalToGlobalMapping(B,mat->rmap->mapping,mat->cmap->mapping);
2244:   MatDuplicate(matis->A,op,&localmat);
2245:   MatISSetLocalMat(B,localmat);
2246:   MatDestroy(&localmat);
2247:   if (matis->sf) {
2248:     Mat_IS *bmatis = (Mat_IS*)(B->data);

2250:     PetscObjectReference((PetscObject)matis->sf);
2251:     bmatis->sf = matis->sf;
2252:     PetscMalloc2(matis->sf->nroots,&bmatis->sf_rootdata,matis->sf->nleaves,&bmatis->sf_leafdata);
2253:     if (matis->sf != matis->csf) {
2254:       PetscObjectReference((PetscObject)matis->csf);
2255:       bmatis->csf = matis->csf;
2256:       PetscMalloc2(matis->csf->nroots,&bmatis->csf_rootdata,matis->csf->nleaves,&bmatis->csf_leafdata);
2257:     } else {
2258:       bmatis->csf          = bmatis->sf;
2259:       bmatis->csf_leafdata = bmatis->sf_leafdata;
2260:       bmatis->csf_rootdata = bmatis->sf_rootdata;
2261:     }
2262:   }
2263:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
2264:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
2265:   *newmat = B;
2266:   return(0);
2267: }

2269: static PetscErrorCode MatIsHermitian_IS(Mat A,PetscReal tol,PetscBool  *flg)
2270: {
2272:   Mat_IS         *matis = (Mat_IS*)A->data;
2273:   PetscBool      local_sym;

2276:   MatIsHermitian(matis->A,tol,&local_sym);
2277:   MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
2278:   return(0);
2279: }

2281: static PetscErrorCode MatIsSymmetric_IS(Mat A,PetscReal tol,PetscBool  *flg)
2282: {
2284:   Mat_IS         *matis = (Mat_IS*)A->data;
2285:   PetscBool      local_sym;

2288:   MatIsSymmetric(matis->A,tol,&local_sym);
2289:   MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
2290:   return(0);
2291: }

2293: static PetscErrorCode MatIsStructurallySymmetric_IS(Mat A,PetscBool  *flg)
2294: {
2296:   Mat_IS         *matis = (Mat_IS*)A->data;
2297:   PetscBool      local_sym;

2300:   if (A->rmap->mapping != A->cmap->mapping) {
2301:     *flg = PETSC_FALSE;
2302:     return(0);
2303:   }
2304:   MatIsStructurallySymmetric(matis->A,&local_sym);
2305:   MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
2306:   return(0);
2307: }

2309: static PetscErrorCode MatDestroy_IS(Mat A)
2310: {
2312:   Mat_IS         *b = (Mat_IS*)A->data;

2315:   PetscFree(b->lmattype);
2316:   MatDestroy(&b->A);
2317:   VecScatterDestroy(&b->cctx);
2318:   VecScatterDestroy(&b->rctx);
2319:   VecDestroy(&b->x);
2320:   VecDestroy(&b->y);
2321:   VecDestroy(&b->counter);
2322:   ISDestroy(&b->getsub_ris);
2323:   ISDestroy(&b->getsub_cis);
2324:   if (b->sf != b->csf) {
2325:     PetscSFDestroy(&b->csf);
2326:     PetscFree2(b->csf_rootdata,b->csf_leafdata);
2327:   } else b->csf = NULL;
2328:   PetscSFDestroy(&b->sf);
2329:   PetscFree2(b->sf_rootdata,b->sf_leafdata);
2330:   PetscFree(A->data);
2331:   PetscObjectChangeTypeName((PetscObject)A,0);
2332:   PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMatType_C",NULL);
2333:   PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL);
2334:   PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",NULL);
2335:   PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",NULL);
2336:   PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",NULL);
2337:   PetscObjectComposeFunction((PetscObject)A,"MatISSetUpSF_C",NULL);
2338:   PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",NULL);
2339:   PetscObjectComposeFunction((PetscObject)A,"MatISFixLocalEmpty_C",NULL);
2340:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpiaij_C",NULL);
2341:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpibaij_C",NULL);
2342:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpisbaij_C",NULL);
2343:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqaij_C",NULL);
2344:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqbaij_C",NULL);
2345:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqsbaij_C",NULL);
2346:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_aij_C",NULL);
2347:   return(0);
2348: }

2350: static PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y)
2351: {
2353:   Mat_IS         *is  = (Mat_IS*)A->data;
2354:   PetscScalar    zero = 0.0;

2357:   /*  scatter the global vector x into the local work vector */
2358:   VecScatterBegin(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
2359:   VecScatterEnd(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);

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

2364:   /* scatter product back into global memory */
2365:   VecSet(y,zero);
2366:   VecScatterBegin(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
2367:   VecScatterEnd(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
2368:   return(0);
2369: }

2371: static PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
2372: {
2373:   Vec            temp_vec;

2377:   if (v3 != v2) {
2378:     MatMult(A,v1,v3);
2379:     VecAXPY(v3,1.0,v2);
2380:   } else {
2381:     VecDuplicate(v2,&temp_vec);
2382:     MatMult(A,v1,temp_vec);
2383:     VecAXPY(temp_vec,1.0,v2);
2384:     VecCopy(temp_vec,v3);
2385:     VecDestroy(&temp_vec);
2386:   }
2387:   return(0);
2388: }

2390: static PetscErrorCode MatMultTranspose_IS(Mat A,Vec y,Vec x)
2391: {
2392:   Mat_IS         *is = (Mat_IS*)A->data;

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

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

2403:   /* scatter product back into global vector */
2404:   VecSet(x,0);
2405:   VecScatterBegin(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);
2406:   VecScatterEnd(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);
2407:   return(0);
2408: }

2410: static PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
2411: {
2412:   Vec            temp_vec;

2416:   if (v3 != v2) {
2417:     MatMultTranspose(A,v1,v3);
2418:     VecAXPY(v3,1.0,v2);
2419:   } else {
2420:     VecDuplicate(v2,&temp_vec);
2421:     MatMultTranspose(A,v1,temp_vec);
2422:     VecAXPY(temp_vec,1.0,v2);
2423:     VecCopy(temp_vec,v3);
2424:     VecDestroy(&temp_vec);
2425:   }
2426:   return(0);
2427: }

2429: static PetscErrorCode MatView_IS(Mat A,PetscViewer viewer)
2430: {
2431:   Mat_IS         *a = (Mat_IS*)A->data;
2433:   PetscViewer    sviewer;
2434:   PetscBool      isascii,view = PETSC_TRUE;

2437:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
2438:   if (isascii)  {
2439:     PetscViewerFormat format;

2441:     PetscViewerGetFormat(viewer,&format);
2442:     if (format == PETSC_VIEWER_ASCII_INFO) view = PETSC_FALSE;
2443:   }
2444:   if (!view) return(0);
2445:   PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
2446:   MatView(a->A,sviewer);
2447:   PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
2448:   PetscViewerFlush(viewer);
2449:   return(0);
2450: }

2452: static PetscErrorCode MatISSetUpScatters_Private(Mat A)
2453: {
2454:   Vec            cglobal,rglobal;
2455:   IS             from;
2456:   Mat_IS         *is = (Mat_IS*)A->data;
2457:   const PetscInt *garray;
2458:   PetscInt       nr,rbs,nc,cbs;
2459:   PetscBool      iscuda;

2463:   ISLocalToGlobalMappingGetSize(A->rmap->mapping,&nr);
2464:   ISLocalToGlobalMappingGetBlockSize(A->rmap->mapping,&rbs);
2465:   ISLocalToGlobalMappingGetSize(A->cmap->mapping,&nc);
2466:   ISLocalToGlobalMappingGetBlockSize(A->cmap->mapping,&cbs);
2467:   VecDestroy(&is->x);
2468:   VecDestroy(&is->y);
2469:   VecDestroy(&is->counter);
2470:   VecScatterDestroy(&is->rctx);
2471:   VecScatterDestroy(&is->cctx);
2472:   MatCreateVecs(is->A,&is->x,&is->y);
2473:   PetscObjectTypeCompare((PetscObject)is->y,VECSEQCUDA,&iscuda);
2474:   if (iscuda) {
2475:     PetscFree(A->defaultvectype);
2476:     PetscStrallocpy(VECCUDA,&A->defaultvectype);
2477:   }
2478:   MatCreateVecs(A,&cglobal,&rglobal);
2479:   ISLocalToGlobalMappingGetBlockIndices(A->rmap->mapping,&garray);
2480:   ISCreateBlock(PetscObjectComm((PetscObject)A),rbs,nr/rbs,garray,PETSC_USE_POINTER,&from);
2481:   VecScatterCreate(rglobal,from,is->y,NULL,&is->rctx);
2482:   ISLocalToGlobalMappingRestoreBlockIndices(A->rmap->mapping,&garray);
2483:   ISDestroy(&from);
2484:   if (A->rmap->mapping != A->cmap->mapping) {
2485:     ISLocalToGlobalMappingGetBlockIndices(A->cmap->mapping,&garray);
2486:     ISCreateBlock(PetscObjectComm((PetscObject)A),cbs,nc/cbs,garray,PETSC_USE_POINTER,&from);
2487:     VecScatterCreate(cglobal,from,is->x,NULL,&is->cctx);
2488:     ISLocalToGlobalMappingRestoreBlockIndices(A->cmap->mapping,&garray);
2489:     ISDestroy(&from);
2490:   } else {
2491:     PetscObjectReference((PetscObject)is->rctx);
2492:     is->cctx = is->rctx;
2493:   }
2494:   /* interface counter vector (local) */
2495:   VecDuplicate(is->y,&is->counter);
2496:   VecSet(is->y,1.);
2497:   VecScatterBegin(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);
2498:   VecScatterEnd(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);
2499:   VecScatterBegin(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);
2500:   VecScatterEnd(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);
2501:   VecDestroy(&rglobal);
2502:   VecDestroy(&cglobal);
2503:   ISDestroy(&from);
2504:   return(0);
2505: }

2507: static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping)
2508: {
2510:   PetscInt       nr,rbs,nc,cbs;
2511:   Mat_IS         *is = (Mat_IS*)A->data;

2516:   MatDestroy(&is->A);
2517:   if (is->csf != is->sf) {
2518:     PetscSFDestroy(&is->csf);
2519:     PetscFree2(is->csf_rootdata,is->csf_leafdata);
2520:   } else is->csf = NULL;
2521:   PetscSFDestroy(&is->sf);
2522:   PetscFree2(is->sf_rootdata,is->sf_leafdata);

2524:   /* Setup Layout and set local to global maps */
2525:   PetscLayoutSetUp(A->rmap);
2526:   PetscLayoutSetUp(A->cmap);
2527:   ISLocalToGlobalMappingGetSize(rmapping,&nr);
2528:   ISLocalToGlobalMappingGetBlockSize(rmapping,&rbs);
2529:   ISLocalToGlobalMappingGetSize(cmapping,&nc);
2530:   ISLocalToGlobalMappingGetBlockSize(cmapping,&cbs);
2531:   /* check if the two mappings are actually the same for square matrices (DOLFIN passes 2 different objects) */
2532:   if (rmapping != cmapping && A->rmap->N == A->cmap->N) {
2533:     PetscBool same,gsame;

2535:     same = PETSC_FALSE;
2536:     if (nr == nc && cbs == rbs) {
2537:       const PetscInt *idxs1,*idxs2;

2539:       ISLocalToGlobalMappingGetBlockIndices(rmapping,&idxs1);
2540:       ISLocalToGlobalMappingGetBlockIndices(cmapping,&idxs2);
2541:       PetscMemcmp(idxs1,idxs2,(nr/rbs)*sizeof(PetscInt),&same);
2542:       ISLocalToGlobalMappingRestoreBlockIndices(rmapping,&idxs1);
2543:       ISLocalToGlobalMappingRestoreBlockIndices(cmapping,&idxs2);
2544:     }
2545:     MPIU_Allreduce(&same,&gsame,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
2546:     if (gsame) cmapping = rmapping;
2547:   }
2548:   PetscLayoutSetBlockSize(A->rmap,rbs);
2549:   PetscLayoutSetBlockSize(A->cmap,cbs);
2550:   PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping);
2551:   PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping);

2553:   /* Create the local matrix A */
2554:   MatCreate(PETSC_COMM_SELF,&is->A);
2555:   MatSetType(is->A,is->lmattype);
2556:   MatSetSizes(is->A,nr,nc,nr,nc);
2557:   MatSetBlockSizes(is->A,rbs,cbs);
2558:   MatSetOptionsPrefix(is->A,"is_");
2559:   MatAppendOptionsPrefix(is->A,((PetscObject)A)->prefix);
2560:   PetscLayoutSetUp(is->A->rmap);
2561:   PetscLayoutSetUp(is->A->cmap);

2563:   if (!is->islocalref) { /* setup scatters and local vectors for MatMult */
2564:     MatISSetUpScatters_Private(A);
2565:   }
2566:   MatSetUp(A);
2567:   return(0);
2568: }

2570: static PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv)
2571: {
2572:   Mat_IS         *is = (Mat_IS*)mat->data;
2574: #if defined(PETSC_USE_DEBUG)
2575:   PetscInt       i,zm,zn;
2576: #endif
2577:   PetscInt       rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION];

2580: #if defined(PETSC_USE_DEBUG)
2581:   if (m > MATIS_MAX_ENTRIES_INSERTION || n > MATIS_MAX_ENTRIES_INSERTION) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of row/column indices must be <= %D: they are %D %D",MATIS_MAX_ENTRIES_INSERTION,m,n);
2582:   /* count negative indices */
2583:   for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++;
2584:   for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++;
2585: #endif
2586:   ISGlobalToLocalMappingApply(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);
2587:   ISGlobalToLocalMappingApply(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);
2588: #if defined(PETSC_USE_DEBUG)
2589:   /* count negative indices (should be the same as before) */
2590:   for (i=0;i<m;i++) if (rows_l[i] < 0) zm--;
2591:   for (i=0;i<n;i++) if (cols_l[i] < 0) zn--;
2592:   if (!is->A->rmap->mapping && zm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the row indices can not be mapped! Maybe you should not use MATIS");
2593:   if (!is->A->cmap->mapping && zn) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the column indices can not be mapped! Maybe you should not use MATIS");
2594: #endif
2595:   MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);
2596:   return(0);
2597: }

2599: static PetscErrorCode MatSetValuesBlocked_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv)
2600: {
2601:   Mat_IS         *is = (Mat_IS*)mat->data;
2603: #if defined(PETSC_USE_DEBUG)
2604:   PetscInt       i,zm,zn;
2605: #endif
2606:   PetscInt       rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION];

2609: #if defined(PETSC_USE_DEBUG)
2610:   if (m > MATIS_MAX_ENTRIES_INSERTION || n > MATIS_MAX_ENTRIES_INSERTION) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of row/column block indices must be <= %D: they are %D %D",MATIS_MAX_ENTRIES_INSERTION,m,n);
2611:   /* count negative indices */
2612:   for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++;
2613:   for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++;
2614: #endif
2615:   ISGlobalToLocalMappingApplyBlock(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);
2616:   ISGlobalToLocalMappingApplyBlock(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);
2617: #if defined(PETSC_USE_DEBUG)
2618:   /* count negative indices (should be the same as before) */
2619:   for (i=0;i<m;i++) if (rows_l[i] < 0) zm--;
2620:   for (i=0;i<n;i++) if (cols_l[i] < 0) zn--;
2621:   if (!is->A->rmap->mapping && zm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the row indices can not be mapped! Maybe you should not use MATIS");
2622:   if (!is->A->cmap->mapping && zn) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the column indices can not be mapped! Maybe you should not use MATIS");
2623: #endif
2624:   MatSetValuesBlocked(is->A,m,rows_l,n,cols_l,values,addv);
2625:   return(0);
2626: }

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

2634:   if (is->A->rmap->mapping) {
2635:     MatSetValuesLocal(is->A,m,rows,n,cols,values,addv);
2636:   } else {
2637:     MatSetValues(is->A,m,rows,n,cols,values,addv);
2638:   }
2639:   return(0);
2640: }

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

2648:   if (is->A->rmap->mapping) {
2649: #if defined(PETSC_USE_DEBUG)
2650:     PetscInt ibs,bs;

2652:     ISLocalToGlobalMappingGetBlockSize(is->A->rmap->mapping,&ibs);
2653:     MatGetBlockSize(is->A,&bs);
2654:     if (ibs != bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Different block sizes! local mat %D, local l2g map %D",bs,ibs);
2655: #endif
2656:     MatSetValuesBlockedLocal(is->A,m,rows,n,cols,values,addv);
2657:   } else {
2658:     MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv);
2659:   }
2660:   return(0);
2661: }

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

2669:   if (!n) {
2670:     is->pure_neumann = PETSC_TRUE;
2671:   } else {
2672:     PetscInt i;
2673:     is->pure_neumann = PETSC_FALSE;

2675:     if (columns) {
2676:       MatZeroRowsColumns(is->A,n,rows,diag,0,0);
2677:     } else {
2678:       MatZeroRows(is->A,n,rows,diag,0,0);
2679:     }
2680:     if (diag != 0.) {
2681:       const PetscScalar *array;
2682:       VecGetArrayRead(is->counter,&array);
2683:       for (i=0; i<n; i++) {
2684:         MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);
2685:       }
2686:       VecRestoreArrayRead(is->counter,&array);
2687:     }
2688:     MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);
2689:     MatAssemblyEnd(is->A,MAT_FINAL_ASSEMBLY);
2690:   }
2691:   return(0);
2692: }

2694: static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b,PetscBool columns)
2695: {
2696:   Mat_IS         *matis = (Mat_IS*)A->data;
2697:   PetscInt       nr,nl,len,i;
2698:   PetscInt       *lrows;

2702: #if defined(PETSC_USE_DEBUG)
2703:   if (columns || diag != 0. || (x && b)) {
2704:     PetscBool cong;

2706:     PetscLayoutCompare(A->rmap,A->cmap,&cong);
2707:     cong = (PetscBool)(cong && matis->sf == matis->csf);
2708:     if (!cong && columns) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Columns can be zeroed if and only if A->rmap and A->cmap are congruent and the l2g maps are the same for MATIS");
2709:     if (!cong && diag != 0.) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Nonzero diagonal value supported if and only if A->rmap and A->cmap are congruent and the l2g maps are the same for MATIS");
2710:     if (!cong && x && b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"A->rmap and A->cmap need to be congruent, and the l2g maps be the same");
2711:   }
2712: #endif
2713:   /* get locally owned rows */
2714:   PetscLayoutMapLocal_Private(A->rmap,n,rows,&len,&lrows,NULL);
2715:   /* fix right hand side if needed */
2716:   if (x && b) {
2717:     const PetscScalar *xx;
2718:     PetscScalar       *bb;

2720:     VecGetArrayRead(x, &xx);
2721:     VecGetArray(b, &bb);
2722:     for (i=0;i<len;++i) bb[lrows[i]] = diag*xx[lrows[i]];
2723:     VecRestoreArrayRead(x, &xx);
2724:     VecRestoreArray(b, &bb);
2725:   }
2726:   /* get rows associated to the local matrices */
2727:   MatISSetUpSF(A);
2728:   MatGetSize(matis->A,&nl,NULL);
2729:   PetscMemzero(matis->sf_leafdata,nl*sizeof(PetscInt));
2730:   PetscMemzero(matis->sf_rootdata,A->rmap->n*sizeof(PetscInt));
2731:   for (i=0;i<len;i++) matis->sf_rootdata[lrows[i]] = 1;
2732:   PetscFree(lrows);
2733:   PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
2734:   PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
2735:   PetscMalloc1(nl,&lrows);
2736:   for (i=0,nr=0;i<nl;i++) if (matis->sf_leafdata[i]) lrows[nr++] = i;
2737:   MatISZeroRowsColumnsLocal_Private(A,nr,lrows,diag,columns);
2738:   PetscFree(lrows);
2739:   return(0);
2740: }

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

2747:   MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_FALSE);
2748:   return(0);
2749: }

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

2756:   MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_TRUE);
2757:   return(0);
2758: }

2760: static PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type)
2761: {
2762:   Mat_IS         *is = (Mat_IS*)A->data;

2766:   MatAssemblyBegin(is->A,type);
2767:   return(0);
2768: }

2770: static PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type)
2771: {
2772:   Mat_IS         *is = (Mat_IS*)A->data;

2776:   MatAssemblyEnd(is->A,type);
2777:   /* fix for local empty rows/cols */
2778:   if (is->locempty && type == MAT_FINAL_ASSEMBLY) {
2779:     Mat                    newlA;
2780:     ISLocalToGlobalMapping rl2g,cl2g;
2781:     IS                     nzr,nzc;
2782:     PetscInt               nr,nc,nnzr,nnzc;
2783:     PetscBool              lnewl2g,newl2g;

2785:     MatGetSize(is->A,&nr,&nc);
2786:     MatFindNonzeroRowsOrCols_Basic(is->A,PETSC_FALSE,PETSC_SMALL,&nzr);
2787:     if (!nzr) {
2788:       ISCreateStride(PetscObjectComm((PetscObject)is->A),nr,0,1,&nzr);
2789:     }
2790:     MatFindNonzeroRowsOrCols_Basic(is->A,PETSC_TRUE,PETSC_SMALL,&nzc);
2791:     if (!nzc) {
2792:       ISCreateStride(PetscObjectComm((PetscObject)is->A),nc,0,1,&nzc);
2793:     }
2794:     ISGetSize(nzr,&nnzr);
2795:     ISGetSize(nzc,&nnzc);
2796:     if (nnzr != nr || nnzc != nc) {
2797:       ISLocalToGlobalMapping l2g;
2798:       IS                     is1,is2;

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

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

2807:       /* attach local l2g maps for successive calls of MatSetValues on the local matrix */
2808:       ISLocalToGlobalMappingCreateIS(nzr,&l2g);
2809:       ISCreateStride(PetscObjectComm((PetscObject)is->A),nr,0,1,&is1);
2810:       ISGlobalToLocalMappingApplyIS(l2g,IS_GTOLM_MASK,is1,&is2);
2811:       ISLocalToGlobalMappingDestroy(&l2g);
2812:       if (is->A->rmap->mapping) { /* the matrix has a local-to-local map already attached (probably DMDA generated) */
2813:         const PetscInt *idxs1,*idxs2;
2814:         PetscInt       j,i,nl,*tidxs;

2816:         ISLocalToGlobalMappingGetSize(is->A->rmap->mapping,&nl);
2817:         ISLocalToGlobalMappingGetIndices(is->A->rmap->mapping,&idxs1);
2818:         PetscMalloc1(nl,&tidxs);
2819:         ISGetIndices(is2,&idxs2);
2820:         for (i=0,j=0;i<nl;i++) tidxs[i] = idxs1[i] < 0 ? -1 : idxs2[j++];
2821:         if (j != nr) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected count %D != %D",j,nr);
2822:         ISLocalToGlobalMappingRestoreIndices(is->A->rmap->mapping,&idxs1);
2823:         ISRestoreIndices(is2,&idxs2);
2824:         ISDestroy(&is2);
2825:         ISCreateGeneral(PetscObjectComm((PetscObject)is->A->rmap->mapping),nl,tidxs,PETSC_OWN_POINTER,&is2);
2826:       }
2827:       ISLocalToGlobalMappingCreateIS(is2,&rl2g);
2828:       ISDestroy(&is1);
2829:       ISDestroy(&is2);

2831:       ISLocalToGlobalMappingCreateIS(nzc,&l2g);
2832:       ISCreateStride(PetscObjectComm((PetscObject)is->A),nc,0,1,&is1);
2833:       ISGlobalToLocalMappingApplyIS(l2g,IS_GTOLM_MASK,is1,&is2);
2834:       ISLocalToGlobalMappingDestroy(&l2g);
2835:       if (is->A->cmap->mapping) { /* the matrix has a local-to-local map already attached (probably DMDA generated) */
2836:         const PetscInt *idxs1,*idxs2;
2837:         PetscInt       j,i,nl,*tidxs;

2839:         ISLocalToGlobalMappingGetSize(is->A->cmap->mapping,&nl);
2840:         ISLocalToGlobalMappingGetIndices(is->A->cmap->mapping,&idxs1);
2841:         PetscMalloc1(nl,&tidxs);
2842:         ISGetIndices(is2,&idxs2);
2843:         for (i=0,j=0;i<nl;i++) tidxs[i] = idxs1[i] < 0 ? -1 : idxs2[j++];
2844:         if (j != nc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected count %D != %D",j,nc);
2845:         ISLocalToGlobalMappingRestoreIndices(is->A->cmap->mapping,&idxs1);
2846:         ISRestoreIndices(is2,&idxs2);
2847:         ISDestroy(&is2);
2848:         ISCreateGeneral(PetscObjectComm((PetscObject)is->A->cmap->mapping),nl,tidxs,PETSC_OWN_POINTER,&is2);
2849:       }
2850:       ISLocalToGlobalMappingCreateIS(is2,&cl2g);
2851:       ISDestroy(&is1);
2852:       ISDestroy(&is2);

2854:       MatSetLocalToGlobalMapping(newlA,rl2g,cl2g);

2856:       ISLocalToGlobalMappingDestroy(&rl2g);
2857:       ISLocalToGlobalMappingDestroy(&cl2g);
2858:     } else { /* local matrix fully populated */
2859:       lnewl2g = PETSC_FALSE;
2860:       MPI_Allreduce(&lnewl2g,&newl2g,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A));
2861:       PetscObjectReference((PetscObject)is->A);
2862:       newlA   = is->A;
2863:     }
2864:     /* attach new global l2g map if needed */
2865:     if (newl2g) {
2866:       IS             gnzr,gnzc;
2867:       const PetscInt *grid,*gcid;

2869:       ISLocalToGlobalMappingApplyIS(A->rmap->mapping,nzr,&gnzr);
2870:       ISLocalToGlobalMappingApplyIS(A->cmap->mapping,nzc,&gnzc);
2871:       ISGetIndices(gnzr,&grid);
2872:       ISGetIndices(gnzc,&gcid);
2873:       ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,nnzr,grid,PETSC_COPY_VALUES,&rl2g);
2874:       ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,nnzc,gcid,PETSC_COPY_VALUES,&cl2g);
2875:       ISRestoreIndices(gnzr,&grid);
2876:       ISRestoreIndices(gnzc,&gcid);
2877:       ISDestroy(&gnzr);
2878:       ISDestroy(&gnzc);
2879:       MatSetLocalToGlobalMapping(A,rl2g,cl2g);
2880:       ISLocalToGlobalMappingDestroy(&rl2g);
2881:       ISLocalToGlobalMappingDestroy(&cl2g);
2882:     }
2883:     MatISSetLocalMat(A,newlA);
2884:     MatDestroy(&newlA);
2885:     ISDestroy(&nzr);
2886:     ISDestroy(&nzc);
2887:     is->locempty = PETSC_FALSE;
2888:   }
2889:   return(0);
2890: }

2892: static PetscErrorCode MatISGetLocalMat_IS(Mat mat,Mat *local)
2893: {
2894:   Mat_IS *is = (Mat_IS*)mat->data;

2897:   *local = is->A;
2898:   return(0);
2899: }

2901: static PetscErrorCode MatISRestoreLocalMat_IS(Mat mat,Mat *local)
2902: {
2904:   *local = NULL;
2905:   return(0);
2906: }

2908: /*@
2909:     MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix.

2911:   Input Parameter:
2912: .  mat - the matrix

2914:   Output Parameter:
2915: .  local - the local matrix

2917:   Level: advanced

2919:   Notes:
2920:     This can be called if you have precomputed the nonzero structure of the
2921:   matrix and want to provide it to the inner matrix object to improve the performance
2922:   of the MatSetValues() operation.

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

2926: .seealso: MATIS
2927: @*/
2928: PetscErrorCode MatISGetLocalMat(Mat mat,Mat *local)
2929: {

2935:   PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat*),(mat,local));
2936:   return(0);
2937: }

2939: /*@
2940:     MatISRestoreLocalMat - Restores the local matrix obtained with MatISGetLocalMat()

2942:   Input Parameter:
2943: .  mat - the matrix

2945:   Output Parameter:
2946: .  local - the local matrix

2948:   Level: advanced

2950: .seealso: MATIS
2951: @*/
2952: PetscErrorCode MatISRestoreLocalMat(Mat mat,Mat *local)
2953: {

2959:   PetscUseMethod(mat,"MatISRestoreLocalMat_C",(Mat,Mat*),(mat,local));
2960:   return(0);
2961: }

2963: static PetscErrorCode MatISSetLocalMatType_IS(Mat mat,MatType mtype)
2964: {
2965:   Mat_IS         *is = (Mat_IS*)mat->data;

2969:   if (is->A) {
2970:     MatSetType(is->A,mtype);
2971:   }
2972:   PetscFree(is->lmattype);
2973:   PetscStrallocpy(mtype,&is->lmattype);
2974:   return(0);
2975: }

2977: /*@
2978:     MatISSetLocalMatType - Specifies the type of local matrix

2980:   Input Parameter:
2981: .  mat - the matrix
2982: .  mtype - the local matrix type

2984:   Output Parameter:

2986:   Level: advanced

2988: .seealso: MATIS, MatSetType(), MatType
2989: @*/
2990: PetscErrorCode MatISSetLocalMatType(Mat mat,MatType mtype)
2991: {

2996:   PetscUseMethod(mat,"MatISSetLocalMatType_C",(Mat,MatType),(mat,mtype));
2997:   return(0);
2998: }

3000: static PetscErrorCode MatISSetLocalMat_IS(Mat mat,Mat local)
3001: {
3002:   Mat_IS         *is = (Mat_IS*)mat->data;
3003:   PetscInt       nrows,ncols,orows,ocols;
3005:   MatType        mtype,otype;
3006:   PetscBool      sametype = PETSC_TRUE;

3009:   if (is->A) {
3010:     MatGetSize(is->A,&orows,&ocols);
3011:     MatGetSize(local,&nrows,&ncols);
3012:     if (orows != nrows || ocols != ncols) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local MATIS matrix should be of size %Dx%D (you passed a %Dx%D matrix)",orows,ocols,nrows,ncols);
3013:     MatGetType(local,&mtype);
3014:     MatGetType(is->A,&otype);
3015:     PetscStrcmp(mtype,otype,&sametype);
3016:   }
3017:   PetscObjectReference((PetscObject)local);
3018:   MatDestroy(&is->A);
3019:   is->A = local;
3020:   MatGetType(is->A,&mtype);
3021:   MatISSetLocalMatType(mat,mtype);
3022:   if (!sametype && !is->islocalref) {
3023:     MatISSetUpScatters_Private(mat);
3024:   }
3025:   return(0);
3026: }

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

3031:   Collective on Mat

3033:   Input Parameter:
3034: .  mat - the matrix
3035: .  local - the local matrix

3037:   Output Parameter:

3039:   Level: advanced

3041:   Notes:
3042:     This can be called if you have precomputed the local matrix and
3043:   want to provide it to the matrix object MATIS.

3045: .seealso: MATIS
3046: @*/
3047: PetscErrorCode MatISSetLocalMat(Mat mat,Mat local)
3048: {

3054:   PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local));
3055:   return(0);
3056: }

3058: static PetscErrorCode MatZeroEntries_IS(Mat A)
3059: {
3060:   Mat_IS         *a = (Mat_IS*)A->data;

3064:   MatZeroEntries(a->A);
3065:   return(0);
3066: }

3068: static PetscErrorCode MatScale_IS(Mat A,PetscScalar a)
3069: {
3070:   Mat_IS         *is = (Mat_IS*)A->data;

3074:   MatScale(is->A,a);
3075:   return(0);
3076: }

3078: static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v)
3079: {
3080:   Mat_IS         *is = (Mat_IS*)A->data;

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

3087:   /* scatter diagonal back into global vector */
3088:   VecSet(v,0);
3089:   VecScatterBegin(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);
3090:   VecScatterEnd(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);
3091:   return(0);
3092: }

3094: static PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg)
3095: {
3096:   Mat_IS         *a = (Mat_IS*)A->data;

3100:   MatSetOption(a->A,op,flg);
3101:   return(0);
3102: }

3104: static PetscErrorCode MatAXPY_IS(Mat Y,PetscScalar a,Mat X,MatStructure str)
3105: {
3106:   Mat_IS         *y = (Mat_IS*)Y->data;
3107:   Mat_IS         *x;
3108: #if defined(PETSC_USE_DEBUG)
3109:   PetscBool      ismatis;
3110: #endif

3114: #if defined(PETSC_USE_DEBUG)
3115:   PetscObjectTypeCompare((PetscObject)X,MATIS,&ismatis);
3116:   if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot call MatAXPY(Y,a,X,str) with X not of type MATIS");
3117: #endif
3118:   x = (Mat_IS*)X->data;
3119:   MatAXPY(y->A,a,x->A,str);
3120:   return(0);
3121: }

3123: static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A,IS row,IS col,Mat *submat)
3124: {
3125:   Mat                    lA;
3126:   Mat_IS                 *matis;
3127:   ISLocalToGlobalMapping rl2g,cl2g;
3128:   IS                     is;
3129:   const PetscInt         *rg,*rl;
3130:   PetscInt               nrg;
3131:   PetscInt               N,M,nrl,i,*idxs;
3132:   PetscErrorCode         ierr;

3135:   ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&rg);
3136:   ISGetLocalSize(row,&nrl);
3137:   ISGetIndices(row,&rl);
3138:   ISLocalToGlobalMappingGetSize(A->rmap->mapping,&nrg);
3139: #if defined(PETSC_USE_DEBUG)
3140:   for (i=0;i<nrl;i++) if (rl[i]>=nrg) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row index %D -> %D greater then maximum possible %D",i,rl[i],nrg);
3141: #endif
3142:   PetscMalloc1(nrg,&idxs);
3143:   /* map from [0,nrl) to row */
3144:   for (i=0;i<nrl;i++) idxs[i] = rl[i];
3145: #if defined(PETSC_USE_DEBUG)
3146:   for (i=nrl;i<nrg;i++) idxs[i] = nrg;
3147: #else
3148:   for (i=nrl;i<nrg;i++) idxs[i] = -1;
3149: #endif
3150:   ISRestoreIndices(row,&rl);
3151:   ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&rg);
3152:   ISCreateGeneral(PetscObjectComm((PetscObject)A),nrg,idxs,PETSC_OWN_POINTER,&is);
3153:   ISLocalToGlobalMappingCreateIS(is,&rl2g);
3154:   ISDestroy(&is);
3155:   /* compute new l2g map for columns */
3156:   if (col != row || A->rmap->mapping != A->cmap->mapping) {
3157:     const PetscInt *cg,*cl;
3158:     PetscInt       ncg;
3159:     PetscInt       ncl;

3161:     ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&cg);
3162:     ISGetLocalSize(col,&ncl);
3163:     ISGetIndices(col,&cl);
3164:     ISLocalToGlobalMappingGetSize(A->cmap->mapping,&ncg);
3165: #if defined(PETSC_USE_DEBUG)
3166:     for (i=0;i<ncl;i++) if (cl[i]>=ncg) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local column index %D -> %D greater then maximum possible %D",i,cl[i],ncg);
3167: #endif
3168:     PetscMalloc1(ncg,&idxs);
3169:     /* map from [0,ncl) to col */
3170:     for (i=0;i<ncl;i++) idxs[i] = cl[i];
3171: #if defined(PETSC_USE_DEBUG)
3172:     for (i=ncl;i<ncg;i++) idxs[i] = ncg;
3173: #else
3174:     for (i=ncl;i<ncg;i++) idxs[i] = -1;
3175: #endif
3176:     ISRestoreIndices(col,&cl);
3177:     ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&cg);
3178:     ISCreateGeneral(PetscObjectComm((PetscObject)A),ncg,idxs,PETSC_OWN_POINTER,&is);
3179:     ISLocalToGlobalMappingCreateIS(is,&cl2g);
3180:     ISDestroy(&is);
3181:   } else {
3182:     PetscObjectReference((PetscObject)rl2g);
3183:     cl2g = rl2g;
3184:   }
3185:   /* create the MATIS submatrix */
3186:   MatGetSize(A,&M,&N);
3187:   MatCreate(PetscObjectComm((PetscObject)A),submat);
3188:   MatSetSizes(*submat,PETSC_DECIDE,PETSC_DECIDE,M,N);
3189:   MatSetType(*submat,MATIS);
3190:   matis = (Mat_IS*)((*submat)->data);
3191:   matis->islocalref = PETSC_TRUE;
3192:   MatSetLocalToGlobalMapping(*submat,rl2g,cl2g);
3193:   MatISGetLocalMat(A,&lA);
3194:   MatISSetLocalMat(*submat,lA);
3195:   MatSetUp(*submat);
3196:   MatAssemblyBegin(*submat,MAT_FINAL_ASSEMBLY);
3197:   MatAssemblyEnd(*submat,MAT_FINAL_ASSEMBLY);
3198:   ISLocalToGlobalMappingDestroy(&rl2g);
3199:   ISLocalToGlobalMappingDestroy(&cl2g);
3200:   /* remove unsupported ops */
3201:   PetscMemzero((*submat)->ops,sizeof(struct _MatOps));
3202:   (*submat)->ops->destroy               = MatDestroy_IS;
3203:   (*submat)->ops->setvalueslocal        = MatSetValuesLocal_SubMat_IS;
3204:   (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS;
3205:   (*submat)->ops->assemblybegin         = MatAssemblyBegin_IS;
3206:   (*submat)->ops->assemblyend           = MatAssemblyEnd_IS;
3207:   return(0);
3208: }

3210: static PetscErrorCode MatSetFromOptions_IS(PetscOptionItems *PetscOptionsObject, Mat A)
3211: {
3212:   Mat_IS         *a = (Mat_IS*)A->data;
3213:   char           type[256];
3214:   PetscBool      flg;

3218:   PetscOptionsHead(PetscOptionsObject,"MATIS options");
3219:   PetscOptionsBool("-matis_fixempty","Fix local matrices in case of empty local rows/columns","MatISFixLocalEmpty",a->locempty,&a->locempty,NULL);
3220:   PetscOptionsBool("-matis_storel2l","Store local-to-local matrices generated from PtAP operations","MatISStoreL2L",a->storel2l,&a->storel2l,NULL);
3221:   PetscOptionsFList("-matis_localmat_type","Matrix type","MatISSetLocalMatType",MatList,a->lmattype,type,256,&flg);
3222:   if (flg) {
3223:     MatISSetLocalMatType(A,type);
3224:   }
3225:   if (a->A) {
3226:     MatSetFromOptions(a->A);
3227:   }
3228:   PetscOptionsTail();
3229:   return(0);
3230: }

3232: /*@
3233:     MatCreateIS - Creates a "process" unassembled matrix, assembled on each
3234:        process but not across processes.

3236:    Input Parameters:
3237: +     comm    - MPI communicator that will share the matrix
3238: .     bs      - block size of the matrix
3239: .     m,n,M,N - local and/or global sizes of the left and right vector used in matrix vector products
3240: .     rmap    - local to global map for rows
3241: -     cmap    - local to global map for cols

3243:    Output Parameter:
3244: .    A - the resulting matrix

3246:    Level: advanced

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

3254: .seealso: MATIS, MatSetLocalToGlobalMapping()
3255: @*/
3256: PetscErrorCode  MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping rmap,ISLocalToGlobalMapping cmap,Mat *A)
3257: {

3261:   if (!rmap && !cmap) SETERRQ(comm,PETSC_ERR_USER,"You need to provide at least one of the mappings");
3262:   MatCreate(comm,A);
3263:   MatSetSizes(*A,m,n,M,N);
3264:   if (bs > 0) {
3265:     MatSetBlockSize(*A,bs);
3266:   }
3267:   MatSetType(*A,MATIS);
3268:   if (rmap && cmap) {
3269:     MatSetLocalToGlobalMapping(*A,rmap,cmap);
3270:   } else if (!rmap) {
3271:     MatSetLocalToGlobalMapping(*A,cmap,cmap);
3272:   } else {
3273:     MatSetLocalToGlobalMapping(*A,rmap,rmap);
3274:   }
3275:   return(0);
3276: }

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

3284:    Operations Provided:
3285: +  MatMult()
3286: .  MatMultAdd()
3287: .  MatMultTranspose()
3288: .  MatMultTransposeAdd()
3289: .  MatZeroEntries()
3290: .  MatSetOption()
3291: .  MatZeroRows()
3292: .  MatSetValues()
3293: .  MatSetValuesBlocked()
3294: .  MatSetValuesLocal()
3295: .  MatSetValuesBlockedLocal()
3296: .  MatScale()
3297: .  MatGetDiagonal()
3298: .  MatMissingDiagonal()
3299: .  MatDuplicate()
3300: .  MatCopy()
3301: .  MatAXPY()
3302: .  MatCreateSubMatrix()
3303: .  MatGetLocalSubMatrix()
3304: .  MatTranspose()
3305: .  MatPtAP() (with P of AIJ type)
3306: -  MatSetLocalToGlobalMapping()

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

3313:    Notes:
3314:     Options prefix for the inner matrix are given by -is_mat_xxx

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

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

3321:   Level: advanced

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

3325: M*/

3327: PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A)
3328: {
3330:   Mat_IS         *b;

3333:   PetscNewLog(A,&b);
3334:   PetscStrallocpy(MATAIJ,&b->lmattype);
3335:   A->data = (void*)b;

3337:   /* matrix ops */
3338:   PetscMemzero(A->ops,sizeof(struct _MatOps));
3339:   A->ops->mult                    = MatMult_IS;
3340:   A->ops->multadd                 = MatMultAdd_IS;
3341:   A->ops->multtranspose           = MatMultTranspose_IS;
3342:   A->ops->multtransposeadd        = MatMultTransposeAdd_IS;
3343:   A->ops->destroy                 = MatDestroy_IS;
3344:   A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS;
3345:   A->ops->setvalues               = MatSetValues_IS;
3346:   A->ops->setvaluesblocked        = MatSetValuesBlocked_IS;
3347:   A->ops->setvalueslocal          = MatSetValuesLocal_IS;
3348:   A->ops->setvaluesblockedlocal   = MatSetValuesBlockedLocal_IS;
3349:   A->ops->zerorows                = MatZeroRows_IS;
3350:   A->ops->zerorowscolumns         = MatZeroRowsColumns_IS;
3351:   A->ops->assemblybegin           = MatAssemblyBegin_IS;
3352:   A->ops->assemblyend             = MatAssemblyEnd_IS;
3353:   A->ops->view                    = MatView_IS;
3354:   A->ops->zeroentries             = MatZeroEntries_IS;
3355:   A->ops->scale                   = MatScale_IS;
3356:   A->ops->getdiagonal             = MatGetDiagonal_IS;
3357:   A->ops->setoption               = MatSetOption_IS;
3358:   A->ops->ishermitian             = MatIsHermitian_IS;
3359:   A->ops->issymmetric             = MatIsSymmetric_IS;
3360:   A->ops->isstructurallysymmetric = MatIsStructurallySymmetric_IS;
3361:   A->ops->duplicate               = MatDuplicate_IS;
3362:   A->ops->missingdiagonal         = MatMissingDiagonal_IS;
3363:   A->ops->copy                    = MatCopy_IS;
3364:   A->ops->getlocalsubmatrix       = MatGetLocalSubMatrix_IS;
3365:   A->ops->createsubmatrix         = MatCreateSubMatrix_IS;
3366:   A->ops->axpy                    = MatAXPY_IS;
3367:   A->ops->diagonalset             = MatDiagonalSet_IS;
3368:   A->ops->shift                   = MatShift_IS;
3369:   A->ops->transpose               = MatTranspose_IS;
3370:   A->ops->getinfo                 = MatGetInfo_IS;
3371:   A->ops->diagonalscale           = MatDiagonalScale_IS;
3372:   A->ops->setfromoptions          = MatSetFromOptions_IS;

3374:   /* special MATIS functions */
3375:   PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMatType_C",MatISSetLocalMatType_IS);
3376:   PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS);
3377:   PetscObjectComposeFunction((PetscObject)A,"MatISRestoreLocalMat_C",MatISRestoreLocalMat_IS);
3378:   PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS);
3379:   PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",MatConvert_IS_XAIJ);
3380:   PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",MatISSetPreallocation_IS);
3381:   PetscObjectComposeFunction((PetscObject)A,"MatISSetUpSF_C",MatISSetUpSF_IS);
3382:   PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",MatISStoreL2L_IS);
3383:   PetscObjectComposeFunction((PetscObject)A,"MatISFixLocalEmpty_C",MatISFixLocalEmpty_IS);
3384:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpiaij_C",MatConvert_IS_XAIJ);
3385:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpibaij_C",MatConvert_IS_XAIJ);
3386:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpisbaij_C",MatConvert_IS_XAIJ);
3387:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqaij_C",MatConvert_IS_XAIJ);
3388:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqbaij_C",MatConvert_IS_XAIJ);
3389:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqsbaij_C",MatConvert_IS_XAIJ);
3390:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_aij_C",MatConvert_IS_XAIJ);
3391:   PetscObjectChangeTypeName((PetscObject)A,MATIS);
3392:   return(0);
3393: }