Actual source code: matis.c

petsc-3.14.6 2021-03-30
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>
 12: #include <petsc/private/vecimpl.h>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1567:    Collective

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

1573:    Level: advanced

1575:    Notes:

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

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

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

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

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

1607:    Collective

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

1613:    Level: advanced

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

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

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

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

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

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

1643:    Collective

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

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

1665:    Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2113:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2784:       MatSetLocalToGlobalMapping(newlA,rl2g,cl2g);

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

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

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

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

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

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

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

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

2847:   Level: advanced

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

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

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

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

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

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

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

2878:   Level: advanced

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

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

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

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

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

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

2914:   Output Parameter:

2916:   Level: advanced

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

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

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

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

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

2961:   Collective on Mat

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

2967:   Output Parameter:

2969:   Level: advanced

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3166:    Level: advanced

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

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

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

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

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

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

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

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

3216:   Level: advanced

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

3220: M*/

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

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

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

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