Actual source code: matis.c

petsc-master 2016-12-07
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>

 12: #define MATIS_MAX_ENTRIES_INSERTION 2048

 14: static PetscErrorCode MatISComputeSF_Private(Mat);

 18: static PetscErrorCode MatGetInfo_IS(Mat A,MatInfoType flag,MatInfo *ginfo)
 19: {
 20:   Mat_IS         *matis = (Mat_IS*)A->data;
 21:   MatInfo        info;
 22:   PetscReal      isend[6],irecv[6];
 23:   PetscInt       bs;

 27:   MatGetBlockSize(A,&bs);
 28:   if (matis->A->ops->getinfo) {
 29:     MatGetInfo(matis->A,MAT_LOCAL,&info);
 30:     isend[0] = info.nz_used;
 31:     isend[1] = info.nz_allocated;
 32:     isend[2] = info.nz_unneeded;
 33:     isend[3] = info.memory;
 34:     isend[4] = info.mallocs;
 35:   } else {
 36:     isend[0] = 0.;
 37:     isend[1] = 0.;
 38:     isend[2] = 0.;
 39:     isend[3] = 0.;
 40:     isend[4] = 0.;
 41:   }
 42:   isend[5] = matis->A->num_ass;
 43:   if (flag == MAT_LOCAL) {
 44:     ginfo->nz_used      = isend[0];
 45:     ginfo->nz_allocated = isend[1];
 46:     ginfo->nz_unneeded  = isend[2];
 47:     ginfo->memory       = isend[3];
 48:     ginfo->mallocs      = isend[4];
 49:     ginfo->assemblies   = isend[5];
 50:   } else if (flag == MAT_GLOBAL_MAX) {
 51:     MPIU_Allreduce(isend,irecv,6,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)A));

 53:     ginfo->nz_used      = irecv[0];
 54:     ginfo->nz_allocated = irecv[1];
 55:     ginfo->nz_unneeded  = irecv[2];
 56:     ginfo->memory       = irecv[3];
 57:     ginfo->mallocs      = irecv[4];
 58:     ginfo->assemblies   = irecv[5];
 59:   } else if (flag == MAT_GLOBAL_SUM) {
 60:     MPIU_Allreduce(isend,irecv,5,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)A));

 62:     ginfo->nz_used      = irecv[0];
 63:     ginfo->nz_allocated = irecv[1];
 64:     ginfo->nz_unneeded  = irecv[2];
 65:     ginfo->memory       = irecv[3];
 66:     ginfo->mallocs      = irecv[4];
 67:     ginfo->assemblies   = A->num_ass;
 68:   }
 69:   ginfo->block_size        = bs;
 70:   ginfo->fill_ratio_given  = 0;
 71:   ginfo->fill_ratio_needed = 0;
 72:   ginfo->factor_mallocs    = 0;
 73:   return(0);
 74: }

 78: PETSC_INTERN PetscErrorCode MatConvert_Nest_IS(Mat A,MatType type,MatReuse reuse,Mat *newmat)
 79: {
 80:   Mat                    **nest,*snest,**rnest,lA,B;
 81:   IS                     *iscol,*isrow,*islrow,*islcol;
 82:   ISLocalToGlobalMapping rl2g,cl2g;
 83:   MPI_Comm               comm;
 84:   PetscInt               *lr,*lc,*lrb,*lcb,*l2gidxs;
 85:   PetscInt               sti,stl,i,j,nr,nc,rbs,cbs;
 86:   PetscBool              convert,lreuse,*istrans;
 87:   PetscErrorCode         ierr;

 90:   MatNestGetSubMats(A,&nr,&nc,&nest);
 91:   lreuse = PETSC_FALSE;
 92:   rnest  = NULL;
 93:   if (reuse == MAT_REUSE_MATRIX) {
 94:     PetscBool ismatis,isnest;

 96:     PetscObjectTypeCompare((PetscObject)*newmat,MATIS,&ismatis);
 97:     if (!ismatis) SETERRQ1(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_USER,"Cannot reuse matrix of type %s",((PetscObject)(*newmat))->type);
 98:     MatISGetLocalMat(*newmat,&lA);
 99:     PetscObjectTypeCompare((PetscObject)lA,MATNEST,&isnest);
100:     if (isnest) {
101:       MatNestGetSubMats(lA,&i,&j,&rnest);
102:       lreuse = (PetscBool)(i == nr && j == nc);
103:       if (!lreuse) rnest = NULL;
104:     }
105:   }
106:   PetscObjectGetComm((PetscObject)A,&comm);
107:   PetscCalloc4(nr,&lr,nc,&lc,nr,&lrb,nc,&lcb);
108:   PetscCalloc6(nr,&isrow,nc,&iscol,
109:                       nr,&islrow,nc,&islcol,
110:                       nr*nc,&snest,nr*nc,&istrans);
111:   MatNestGetISs(A,isrow,iscol);
112:   for (i=0;i<nr;i++) {
113:     for (j=0;j<nc;j++) {
114:       PetscBool ismatis;
115:       PetscInt  l1,l2,lb1,lb2,ij=i*nc+j;

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

120:       /* Nested matrices should be of type MATIS */
121:       PetscObjectTypeCompare((PetscObject)nest[i][j],MATTRANSPOSEMAT,&istrans[ij]);
122:       if (istrans[ij]) {
123:         Mat T,lT;
124:         MatTransposeGetMat(nest[i][j],&T);
125:         PetscObjectTypeCompare((PetscObject)T,MATIS,&ismatis);
126:         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);
127:         MatISGetLocalMat(T,&lT);
128:         MatCreateTranspose(lT,&snest[ij]);
129:       } else {
130:         PetscObjectTypeCompare((PetscObject)nest[i][j],MATIS,&ismatis);
131:         if (!ismatis) SETERRQ2(comm,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) is not of type MATIS",i,j);
132:         MatISGetLocalMat(nest[i][j],&snest[ij]);
133:       }

135:       /* Check compatibility of local sizes */
136:       MatGetSize(snest[ij],&l1,&l2);
137:       MatGetBlockSizes(snest[ij],&lb1,&lb2);
138:       if (!l1 || !l2) continue;
139:       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);
140:       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);
141:       lr[i] = l1;
142:       lc[j] = l2;
143:       if (lrb[i] && lb1 != lrb[i]) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid block size %D != %D",i,j,lrb[i],lb1);
144:       if (lcb[j] && lb2 != lcb[j]) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid block size %D != %D",i,j,lcb[j],lb2);
145:       lrb[i] = lb1;
146:       lcb[j] = lb2;

148:       /* check compatibilty for local matrix reusage */
149:       if (rnest && !rnest[i][j] != !snest[ij]) lreuse = PETSC_FALSE;
150:     }
151:   }

153: #if defined (PETSC_USE_DEBUG)
154:   /* Check compatibility of l2g maps for rows */
155:   for (i=0;i<nr;i++) {
156:     rl2g = NULL;
157:     for (j=0;j<nc;j++) {
158:       PetscInt n1,n2;

160:       if (!nest[i][j]) continue;
161:       if (istrans[i*nc+j]) {
162:         Mat T;

164:         MatTransposeGetMat(nest[i][j],&T);
165:         MatGetLocalToGlobalMapping(T,NULL,&cl2g);
166:       } else {
167:         MatGetLocalToGlobalMapping(nest[i][j],&cl2g,NULL);
168:       }
169:       ISLocalToGlobalMappingGetSize(cl2g,&n1);
170:       if (!n1) continue;
171:       if (!rl2g) {
172:         rl2g = cl2g;
173:       } else {
174:         const PetscInt *idxs1,*idxs2;
175:         PetscBool      same;

177:         ISLocalToGlobalMappingGetSize(rl2g,&n2);
178:         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);
179:         ISLocalToGlobalMappingGetIndices(cl2g,&idxs1);
180:         ISLocalToGlobalMappingGetIndices(rl2g,&idxs2);
181:         PetscMemcmp(idxs1,idxs2,n1*sizeof(PetscInt),&same);
182:         ISLocalToGlobalMappingRestoreIndices(cl2g,&idxs1);
183:         ISLocalToGlobalMappingRestoreIndices(rl2g,&idxs2);
184:         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);
185:       }
186:     }
187:   }
188:   /* Check compatibility of l2g maps for columns */
189:   for (i=0;i<nc;i++) {
190:     rl2g = NULL;
191:     for (j=0;j<nr;j++) {
192:       PetscInt n1,n2;

194:       if (!nest[j][i]) continue;
195:       if (istrans[j*nc+i]) {
196:         Mat T;

198:         MatTransposeGetMat(nest[j][i],&T);
199:         MatGetLocalToGlobalMapping(T,&cl2g,NULL);
200:       } else {
201:         MatGetLocalToGlobalMapping(nest[j][i],NULL,&cl2g);
202:       }
203:       ISLocalToGlobalMappingGetSize(cl2g,&n1);
204:       if (!n1) continue;
205:       if (!rl2g) {
206:         rl2g = cl2g;
207:       } else {
208:         const PetscInt *idxs1,*idxs2;
209:         PetscBool      same;

211:         ISLocalToGlobalMappingGetSize(rl2g,&n2);
212:         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);
213:         ISLocalToGlobalMappingGetIndices(cl2g,&idxs1);
214:         ISLocalToGlobalMappingGetIndices(rl2g,&idxs2);
215:         PetscMemcmp(idxs1,idxs2,n1*sizeof(PetscInt),&same);
216:         ISLocalToGlobalMappingRestoreIndices(cl2g,&idxs1);
217:         ISLocalToGlobalMappingRestoreIndices(rl2g,&idxs2);
218:         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);
219:       }
220:     }
221:   }
222: #endif

224:   B = NULL;
225:   if (reuse != MAT_REUSE_MATRIX) {
226:     /* Create l2g map for the rows of the new matrix and index sets for the local MATNEST */
227:     for (i=0,stl=0;i<nr;i++) stl += lr[i];
228:     PetscMalloc1(stl,&l2gidxs);
229:     for (i=0,sti=0,stl=0;i<nr;i++) {
230:       Mat            usedmat;
231:       Mat_IS         *matis;
232:       const PetscInt *idxs;
233:       PetscInt       n = lr[i]/lrb[i];

235:       /* local IS for local NEST */
236:       ISCreateStride(PETSC_COMM_SELF,n,sti,lrb[i],&islrow[i]);
237:       sti  += n;

239:       /* l2gmap */
240:       j = 0;
241:       usedmat = nest[i][j];
242:       while (!usedmat && j < nc-1) usedmat = nest[i][++j];
243:       if (!usedmat) SETERRQ(comm,PETSC_ERR_SUP,"Cannot find valid row mat");

245:       if (istrans[i*nc+j]) {
246:         Mat T;
247:         MatTransposeGetMat(usedmat,&T);
248:         usedmat = T;
249:       }
250:       matis = (Mat_IS*)(usedmat->data);
251:       if (!matis->sf) {
252:         MatISComputeSF_Private(usedmat);
253:       }
254:       ISGetIndices(isrow[i],&idxs);
255:       if (istrans[i*nc+j]) {
256:         PetscSFBcastBegin(matis->csf,MPIU_INT,idxs,l2gidxs+stl);
257:         PetscSFBcastEnd(matis->csf,MPIU_INT,idxs,l2gidxs+stl);
258:       } else {
259:         PetscSFBcastBegin(matis->sf,MPIU_INT,idxs,l2gidxs+stl);
260:         PetscSFBcastEnd(matis->sf,MPIU_INT,idxs,l2gidxs+stl);
261:       }
262:       ISRestoreIndices(isrow[i],&idxs);
263:       stl += lr[i];
264:     }
265:     ISLocalToGlobalMappingCreate(comm,1,stl,l2gidxs,PETSC_OWN_POINTER,&rl2g);

267:     /* Create l2g map for columns of the new matrix and index sets for the local MATNEST */
268:     for (i=0,stl=0;i<nc;i++) stl += lc[i];
269:     PetscMalloc1(stl,&l2gidxs);
270:     for (i=0,sti=0,stl=0;i<nc;i++) {
271:       Mat            usedmat;
272:       Mat_IS         *matis;
273:       const PetscInt *idxs;
274:       PetscInt       n = lc[i]/lcb[i];

276:       /* local IS for local NEST */
277:       ISCreateStride(PETSC_COMM_SELF,n,sti,lcb[i],&islcol[i]);
278:       sti  += n;

280:       /* l2gmap */
281:       j = 0;
282:       usedmat = nest[j][i];
283:       while (!usedmat && j < nr-1) usedmat = nest[++j][i];
284:       if (!usedmat) SETERRQ(comm,PETSC_ERR_SUP,"Cannot find valid column mat");
285:       if (istrans[j*nc+i]) {
286:         Mat T;
287:         MatTransposeGetMat(usedmat,&T);
288:         usedmat = T;
289:       }
290:       matis = (Mat_IS*)(usedmat->data);
291:       if (!matis->sf) {
292:         MatISComputeSF_Private(usedmat);
293:       }
294:       ISGetIndices(iscol[i],&idxs);
295:       if (istrans[j*nc+i]) {
296:         PetscSFBcastBegin(matis->sf,MPIU_INT,idxs,l2gidxs+stl);
297:         PetscSFBcastEnd(matis->sf,MPIU_INT,idxs,l2gidxs+stl);
298:       } else {
299:         PetscSFBcastBegin(matis->csf,MPIU_INT,idxs,l2gidxs+stl);
300:         PetscSFBcastEnd(matis->csf,MPIU_INT,idxs,l2gidxs+stl);
301:       }
302:       ISRestoreIndices(iscol[i],&idxs);
303:       stl += lc[i];
304:     }
305:     ISLocalToGlobalMappingCreate(comm,1,stl,l2gidxs,PETSC_OWN_POINTER,&cl2g);

307:     /* Create MATIS */
308:     MatCreate(comm,&B);
309:     MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
310:     MatGetBlockSizes(A,&rbs,&cbs);
311:     MatSetBlockSizes(B,rbs,cbs);
312:     MatSetType(B,MATIS);
313:     MatSetLocalToGlobalMapping(B,rl2g,cl2g);
314:     ISLocalToGlobalMappingDestroy(&rl2g);
315:     ISLocalToGlobalMappingDestroy(&cl2g);
316:     MatCreateNest(PETSC_COMM_SELF,nr,islrow,nc,islcol,snest,&lA);
317:     for (i=0;i<nr*nc;i++) {
318:       if (istrans[i]) {
319:         MatDestroy(&snest[i]);
320:       }
321:     }
322:     MatISSetLocalMat(B,lA);
323:     MatDestroy(&lA);
324:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
325:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
326:     if (reuse == MAT_INPLACE_MATRIX) {
327:       MatHeaderReplace(A,&B);
328:     } else {
329:       *newmat = B;
330:     }
331:   } else {
332:     if (lreuse) {
333:       MatISGetLocalMat(*newmat,&lA);
334:       for (i=0;i<nr;i++) {
335:         for (j=0;j<nc;j++) {
336:           if (snest[i*nc+j]) {
337:             MatNestSetSubMat(lA,i,j,snest[i*nc+j]);
338:             if (istrans[i*nc+j]) {
339:               MatDestroy(&snest[i*nc+j]);
340:             }
341:           }
342:         }
343:       }
344:     } else {
345:       for (i=0,sti=0;i<nr;i++) {
346:         PetscInt n = lr[i]/lrb[i];

348:         ISCreateStride(PETSC_COMM_SELF,n,sti,lrb[i],&islrow[i]);
349:         sti  += n;
350:       }
351:       for (i=0,sti=0;i<nc;i++) {
352:         PetscInt n = lc[i]/lcb[i];

354:         ISCreateStride(PETSC_COMM_SELF,n,sti,lcb[i],&islcol[i]);
355:         sti  += n;
356:       }
357:       MatCreateNest(PETSC_COMM_SELF,nr,islrow,nc,islcol,snest,&lA);
358:       if (istrans[i]) {
359:         MatDestroy(&snest[i]);
360:       }
361:       MatISSetLocalMat(*newmat,lA);
362:       MatDestroy(&lA);
363:     }
364:     MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
365:     MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
366:   }

368:   /* Free workspace */
369:   for (i=0;i<nr;i++) {
370:     ISDestroy(&islrow[i]);
371:   }
372:   for (i=0;i<nc;i++) {
373:     ISDestroy(&islcol[i]);
374:   }
375:   PetscFree4(lr,lc,lrb,lcb);
376:   PetscFree6(isrow,iscol,islrow,islcol,snest,istrans);

378:   /* Create local matrix in MATNEST format */
379:   convert = PETSC_FALSE;
380:   PetscOptionsGetBool(NULL,((PetscObject)A)->prefix,"-matis_convert_local_nest",&convert,NULL);
381:   if (convert) {
382:     Mat M;

384:     MatISGetLocalMat(*newmat,&lA);
385:     MatConvert(lA,MATAIJ,MAT_INITIAL_MATRIX,&M);
386:     MatISSetLocalMat(*newmat,M);
387:     MatDestroy(&M);
388:   }
389:   return(0);
390: }

394: PetscErrorCode MatTranspose_IS(Mat A,MatReuse reuse,Mat *B)
395: {
396:   Mat                    C,lC,lA;
397:   ISLocalToGlobalMapping rl2g,cl2g;
398:   PetscBool              notset = PETSC_FALSE,cong = PETSC_TRUE;
399:   PetscErrorCode         ierr;

402:   if (reuse == MAT_REUSE_MATRIX) {
403:     PetscBool rcong,ccong;

405:     PetscLayoutCompare((*B)->cmap,A->rmap,&rcong);
406:     PetscLayoutCompare((*B)->rmap,A->cmap,&ccong);
407:     cong = (PetscBool)(rcong || ccong);
408:   }
409:   if (reuse == MAT_INITIAL_MATRIX || *B == A || !cong) {
410:     MatCreate(PetscObjectComm((PetscObject)A),&C);
411:   } else {
412:     PetscObjectTypeCompare((PetscObject)*B,MATIS,&notset);
413:     C = *B;
414:   }

416:   if (!notset) {
417:     MatSetSizes(C,A->cmap->n,A->rmap->n,A->cmap->N,A->rmap->N);
418:     MatSetBlockSizes(C,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));
419:     MatSetType(C,MATIS);
420:     MatGetLocalToGlobalMapping(A,&rl2g,&cl2g);
421:     MatSetLocalToGlobalMapping(C,cl2g,rl2g);
422:   }

424:   /* perform local transposition */
425:   MatISGetLocalMat(A,&lA);
426:   MatTranspose(lA,MAT_INITIAL_MATRIX,&lC);
427:   MatISSetLocalMat(C,lC);
428:   MatDestroy(&lC);

430:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
431:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

433:   if (reuse == MAT_INITIAL_MATRIX || *B != A) {
434:     if (!cong) {
435:       MatDestroy(B);
436:     }
437:     *B = C;
438:   } else {
439:     MatHeaderMerge(A,&C);
440:   }
441:   return(0);
442: }

446: PetscErrorCode  MatDiagonalSet_IS(Mat A,Vec D,InsertMode insmode)
447: {
448:   Mat_IS         *is = (Mat_IS*)A->data;

452:   if (!D) { /* this code branch is used by MatShift_IS */
453:     VecSet(is->y,1.);
454:   } else {
455:     VecScatterBegin(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD);
456:     VecScatterEnd(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD);
457:   }
458:   VecPointwiseDivide(is->y,is->y,is->counter);
459:   MatDiagonalSet(is->A,is->y,insmode);
460:   return(0);
461: }

465: PetscErrorCode  MatShift_IS(Mat A,PetscScalar a)
466: {

470:   MatDiagonalSet_IS(A,NULL,ADD_VALUES);
471:   return(0);
472: }

476: static PetscErrorCode MatSetValuesLocal_SubMat_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
477: {
479:   Mat_IS         *is = (Mat_IS*)A->data;
480:   PetscInt       rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION];

483: #if defined(PETSC_USE_DEBUG)
484:   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);
485: #endif
486:   ISLocalToGlobalMappingApply(A->rmap->mapping,m,rows,rows_l);
487:   ISLocalToGlobalMappingApply(A->cmap->mapping,n,cols,cols_l);
488:   MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);
489:   return(0);
490: }

494: static PetscErrorCode MatSetValuesBlockedLocal_SubMat_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
495: {
497:   Mat_IS         *is = (Mat_IS*)A->data;
498:   PetscInt       rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION];

501: #if defined(PETSC_USE_DEBUG)
502:   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);
503: #endif
504:   ISLocalToGlobalMappingApplyBlock(A->rmap->mapping,m,rows,rows_l);
505:   ISLocalToGlobalMappingApplyBlock(A->cmap->mapping,n,cols,cols_l);
506:   MatSetValuesBlocked(is->A,m,rows_l,n,cols_l,values,addv);
507:   return(0);
508: }

512: static PetscErrorCode PetscLayoutMapLocal_Private(PetscLayout map,PetscInt N,const PetscInt idxs[], PetscInt *on,PetscInt **oidxs,PetscInt **ogidxs)
513: {
514:   PetscInt      *owners = map->range;
515:   PetscInt       n      = map->n;
516:   PetscSF        sf;
517:   PetscInt      *lidxs,*work = NULL;
518:   PetscSFNode   *ridxs;
519:   PetscMPIInt    rank;
520:   PetscInt       r, p = 0, len = 0;

524:   if (on) *on = 0;              /* squelch -Wmaybe-uninitialized */
525:   /* Create SF where leaves are input idxs and roots are owned idxs (code adapted from MatZeroRowsMapLocal_Private) */
526:   MPI_Comm_rank(map->comm,&rank);
527:   PetscMalloc1(n,&lidxs);
528:   for (r = 0; r < n; ++r) lidxs[r] = -1;
529:   PetscMalloc1(N,&ridxs);
530:   for (r = 0; r < N; ++r) {
531:     const PetscInt idx = idxs[r];
532:     if (idx < 0 || map->N <= idx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Index %D out of range [0,%D)",idx,map->N);
533:     if (idx < owners[p] || owners[p+1] <= idx) { /* short-circuit the search if the last p owns this idx too */
534:       PetscLayoutFindOwner(map,idx,&p);
535:     }
536:     ridxs[r].rank = p;
537:     ridxs[r].index = idxs[r] - owners[p];
538:   }
539:   PetscSFCreate(map->comm,&sf);
540:   PetscSFSetGraph(sf,n,N,NULL,PETSC_OWN_POINTER,ridxs,PETSC_OWN_POINTER);
541:   PetscSFReduceBegin(sf,MPIU_INT,(PetscInt*)idxs,lidxs,MPI_LOR);
542:   PetscSFReduceEnd(sf,MPIU_INT,(PetscInt*)idxs,lidxs,MPI_LOR);
543:   if (ogidxs) { /* communicate global idxs */
544:     PetscInt cum = 0,start,*work2;

546:     PetscMalloc1(n,&work);
547:     PetscCalloc1(N,&work2);
548:     for (r = 0; r < N; ++r) if (idxs[r] >=0) cum++;
549:     MPI_Scan(&cum,&start,1,MPIU_INT,MPI_SUM,map->comm);
550:     start -= cum;
551:     cum = 0;
552:     for (r = 0; r < N; ++r) if (idxs[r] >=0) work2[r] = start+cum++;
553:     PetscSFReduceBegin(sf,MPIU_INT,work2,work,MPIU_REPLACE);
554:     PetscSFReduceEnd(sf,MPIU_INT,work2,work,MPIU_REPLACE);
555:     PetscFree(work2);
556:   }
557:   PetscSFDestroy(&sf);
558:   /* Compress and put in indices */
559:   for (r = 0; r < n; ++r)
560:     if (lidxs[r] >= 0) {
561:       if (work) work[len] = work[r];
562:       lidxs[len++] = r;
563:     }
564:   if (on) *on = len;
565:   if (oidxs) *oidxs = lidxs;
566:   if (ogidxs) *ogidxs = work;
567:   return(0);
568: }

572: static PetscErrorCode MatGetSubMatrix_IS(Mat mat,IS irow,IS icol,MatReuse scall,Mat *newmat)
573: {
574:   Mat               locmat,newlocmat;
575:   Mat_IS            *newmatis;
576: #if defined(PETSC_USE_DEBUG)
577:   Vec               rtest,ltest;
578:   const PetscScalar *array;
579: #endif
580:   const PetscInt    *idxs;
581:   PetscInt          i,m,n;
582:   PetscErrorCode    ierr;

585:   if (scall == MAT_REUSE_MATRIX) {
586:     PetscBool ismatis;

588:     PetscObjectTypeCompare((PetscObject)*newmat,MATIS,&ismatis);
589:     if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Not of MATIS type");
590:     newmatis = (Mat_IS*)(*newmat)->data;
591:     if (!newmatis->getsub_ris) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Misses local row IS");
592:     if (!newmatis->getsub_cis) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Misses local col IS");
593:   }
594:   /* irow and icol may not have duplicate entries */
595: #if defined(PETSC_USE_DEBUG)
596:   MatCreateVecs(mat,&ltest,&rtest);
597:   ISGetLocalSize(irow,&n);
598:   ISGetIndices(irow,&idxs);
599:   for (i=0;i<n;i++) {
600:     VecSetValue(rtest,idxs[i],1.0,ADD_VALUES);
601:   }
602:   VecAssemblyBegin(rtest);
603:   VecAssemblyEnd(rtest);
604:   VecGetLocalSize(rtest,&n);
605:   VecGetOwnershipRange(rtest,&m,NULL);
606:   VecGetArrayRead(rtest,&array);
607:   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]));
608:   VecRestoreArrayRead(rtest,&array);
609:   ISRestoreIndices(irow,&idxs);
610:   ISGetLocalSize(icol,&n);
611:   ISGetIndices(icol,&idxs);
612:   for (i=0;i<n;i++) {
613:     VecSetValue(ltest,idxs[i],1.0,ADD_VALUES);
614:   }
615:   VecAssemblyBegin(ltest);
616:   VecAssemblyEnd(ltest);
617:   VecGetLocalSize(ltest,&n);
618:   VecGetOwnershipRange(ltest,&m,NULL);
619:   VecGetArrayRead(ltest,&array);
620:   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]));
621:   VecRestoreArrayRead(ltest,&array);
622:   ISRestoreIndices(icol,&idxs);
623:   VecDestroy(&rtest);
624:   VecDestroy(&ltest);
625: #endif
626:   if (scall == MAT_INITIAL_MATRIX) {
627:     Mat_IS                 *matis = (Mat_IS*)mat->data;
628:     ISLocalToGlobalMapping rl2g;
629:     IS                     is;
630:     PetscInt               *lidxs,*lgidxs,*newgidxs;
631:     PetscInt               ll,newloc;
632:     MPI_Comm               comm;

634:     PetscObjectGetComm((PetscObject)mat,&comm);
635:     ISGetLocalSize(irow,&m);
636:     ISGetLocalSize(icol,&n);
637:     MatCreate(comm,newmat);
638:     MatSetType(*newmat,MATIS);
639:     MatSetSizes(*newmat,m,n,PETSC_DECIDE,PETSC_DECIDE);
640:     /* communicate irow to their owners in the layout */
641:     ISGetIndices(irow,&idxs);
642:     PetscLayoutMapLocal_Private(mat->rmap,m,idxs,&ll,&lidxs,&lgidxs);
643:     ISRestoreIndices(irow,&idxs);
644:     if (!matis->sf) { /* setup SF if not yet created and allocate rootdata and leafdata */
645:       MatISComputeSF_Private(mat);
646:     }
647:     PetscMemzero(matis->sf_rootdata,matis->sf_nroots*sizeof(PetscInt));
648:     for (i=0;i<ll;i++) matis->sf_rootdata[lidxs[i]] = lgidxs[i]+1;
649:     PetscFree(lidxs);
650:     PetscFree(lgidxs);
651:     PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
652:     PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
653:     for (i=0,newloc=0;i<matis->sf_nleaves;i++) if (matis->sf_leafdata[i]) newloc++;
654:     PetscMalloc1(newloc,&newgidxs);
655:     PetscMalloc1(newloc,&lidxs);
656:     for (i=0,newloc=0;i<matis->sf_nleaves;i++)
657:       if (matis->sf_leafdata[i]) {
658:         lidxs[newloc] = i;
659:         newgidxs[newloc++] = matis->sf_leafdata[i]-1;
660:       }
661:     ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is);
662:     ISLocalToGlobalMappingCreateIS(is,&rl2g);
663:     ISDestroy(&is);
664:     /* local is to extract local submatrix */
665:     newmatis = (Mat_IS*)(*newmat)->data;
666:     ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_ris);
667:     if (mat->congruentlayouts == -1) { /* first time we compare rows and cols layouts */
668:       PetscBool cong;
669:       PetscLayoutCompare(mat->rmap,mat->cmap,&cong);
670:       if (cong) mat->congruentlayouts = 1;
671:       else      mat->congruentlayouts = 0;
672:     }
673:     if (mat->congruentlayouts && irow == icol) {
674:       MatSetLocalToGlobalMapping(*newmat,rl2g,rl2g);
675:       PetscObjectReference((PetscObject)newmatis->getsub_ris);
676:       newmatis->getsub_cis = newmatis->getsub_ris;
677:     } else {
678:       ISLocalToGlobalMapping cl2g;

680:       /* communicate icol to their owners in the layout */
681:       ISGetIndices(icol,&idxs);
682:       PetscLayoutMapLocal_Private(mat->cmap,n,idxs,&ll,&lidxs,&lgidxs);
683:       ISRestoreIndices(icol,&idxs);
684:       PetscMemzero(matis->csf_rootdata,matis->csf_nroots*sizeof(PetscInt));
685:       for (i=0;i<ll;i++) matis->csf_rootdata[lidxs[i]] = lgidxs[i]+1;
686:       PetscFree(lidxs);
687:       PetscFree(lgidxs);
688:       PetscSFBcastBegin(matis->csf,MPIU_INT,matis->csf_rootdata,matis->csf_leafdata);
689:       PetscSFBcastEnd(matis->csf,MPIU_INT,matis->csf_rootdata,matis->csf_leafdata);
690:       for (i=0,newloc=0;i<matis->csf_nleaves;i++) if (matis->csf_leafdata[i]) newloc++;
691:       PetscMalloc1(newloc,&newgidxs);
692:       PetscMalloc1(newloc,&lidxs);
693:       for (i=0,newloc=0;i<matis->csf_nleaves;i++)
694:         if (matis->csf_leafdata[i]) {
695:           lidxs[newloc] = i;
696:           newgidxs[newloc++] = matis->csf_leafdata[i]-1;
697:         }
698:       ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is);
699:       ISLocalToGlobalMappingCreateIS(is,&cl2g);
700:       ISDestroy(&is);
701:       /* local is to extract local submatrix */
702:       ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_cis);
703:       MatSetLocalToGlobalMapping(*newmat,rl2g,cl2g);
704:       ISLocalToGlobalMappingDestroy(&cl2g);
705:     }
706:     ISLocalToGlobalMappingDestroy(&rl2g);
707:   } else {
708:     MatISGetLocalMat(*newmat,&newlocmat);
709:   }
710:   MatISGetLocalMat(mat,&locmat);
711:   newmatis = (Mat_IS*)(*newmat)->data;
712:   MatGetSubMatrix(locmat,newmatis->getsub_ris,newmatis->getsub_cis,scall,&newlocmat);
713:   if (scall == MAT_INITIAL_MATRIX) {
714:     MatISSetLocalMat(*newmat,newlocmat);
715:     MatDestroy(&newlocmat);
716:   }
717:   MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
718:   MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
719:   return(0);
720: }

724: static PetscErrorCode MatCopy_IS(Mat A,Mat B,MatStructure str)
725: {
726:   Mat_IS         *a = (Mat_IS*)A->data,*b;
727:   PetscBool      ismatis;

731:   PetscObjectTypeCompare((PetscObject)B,MATIS,&ismatis);
732:   if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Need to be implemented");
733:   b = (Mat_IS*)B->data;
734:   MatCopy(a->A,b->A,str);
735:   return(0);
736: }

740: static PetscErrorCode MatMissingDiagonal_IS(Mat A,PetscBool  *missing,PetscInt *d)
741: {
742:   Vec               v;
743:   const PetscScalar *array;
744:   PetscInt          i,n;
745:   PetscErrorCode    ierr;

748:   *missing = PETSC_FALSE;
749:   MatCreateVecs(A,NULL,&v);
750:   MatGetDiagonal(A,v);
751:   VecGetLocalSize(v,&n);
752:   VecGetArrayRead(v,&array);
753:   for (i=0;i<n;i++) if (array[i] == 0.) break;
754:   VecRestoreArrayRead(v,&array);
755:   VecDestroy(&v);
756:   if (i != n) *missing = PETSC_TRUE;
757:   if (d) {
758:     *d = -1;
759:     if (*missing) {
760:       PetscInt rstart;
761:       MatGetOwnershipRange(A,&rstart,NULL);
762:       *d = i+rstart;
763:     }
764:   }
765:   return(0);
766: }

770: static PetscErrorCode MatISComputeSF_Private(Mat B)
771: {
772:   Mat_IS         *matis = (Mat_IS*)(B->data);
773:   const PetscInt *gidxs;

777:   MatGetSize(matis->A,&matis->sf_nleaves,NULL);
778:   MatGetLocalSize(B,&matis->sf_nroots,NULL);
779:   PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->sf);
780:   ISLocalToGlobalMappingGetIndices(B->rmap->mapping,&gidxs);
781:   PetscSFSetGraphLayout(matis->sf,B->rmap,matis->sf_nleaves,NULL,PETSC_OWN_POINTER,gidxs);
782:   ISLocalToGlobalMappingRestoreIndices(B->rmap->mapping,&gidxs);
783:   PetscMalloc2(matis->sf_nroots,&matis->sf_rootdata,matis->sf_nleaves,&matis->sf_leafdata);
784:   if (B->rmap->mapping != B->cmap->mapping) { /* setup SF for columns */
785:     MatGetSize(matis->A,NULL,&matis->csf_nleaves);
786:     MatGetLocalSize(B,NULL,&matis->csf_nroots);
787:     PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->csf);
788:     ISLocalToGlobalMappingGetIndices(B->cmap->mapping,&gidxs);
789:     PetscSFSetGraphLayout(matis->csf,B->cmap,matis->csf_nleaves,NULL,PETSC_OWN_POINTER,gidxs);
790:     ISLocalToGlobalMappingRestoreIndices(B->cmap->mapping,&gidxs);
791:     PetscMalloc2(matis->csf_nroots,&matis->csf_rootdata,matis->csf_nleaves,&matis->csf_leafdata);
792:   } else {
793:     matis->csf = matis->sf;
794:     matis->csf_nleaves = matis->sf_nleaves;
795:     matis->csf_nroots = matis->sf_nroots;
796:     matis->csf_leafdata = matis->sf_leafdata;
797:     matis->csf_rootdata = matis->sf_rootdata;
798:   }
799:   return(0);
800: }

804: /*@
805:    MatISSetPreallocation - Preallocates memory for a MATIS parallel matrix.

807:    Collective on MPI_Comm

809:    Input Parameters:
810: +  B - the matrix
811: .  d_nz  - number of nonzeros per row in DIAGONAL portion of local submatrix
812:            (same value is used for all local rows)
813: .  d_nnz - array containing the number of nonzeros in the various rows of the
814:            DIAGONAL portion of the local submatrix (possibly different for each row)
815:            or NULL, if d_nz is used to specify the nonzero structure.
816:            The size of this array is equal to the number of local rows, i.e 'm'.
817:            For matrices that will be factored, you must leave room for (and set)
818:            the diagonal entry even if it is zero.
819: .  o_nz  - number of nonzeros per row in the OFF-DIAGONAL portion of local
820:            submatrix (same value is used for all local rows).
821: -  o_nnz - array containing the number of nonzeros in the various rows of the
822:            OFF-DIAGONAL portion of the local submatrix (possibly different for
823:            each row) or NULL, if o_nz is used to specify the nonzero
824:            structure. The size of this array is equal to the number
825:            of local rows, i.e 'm'.

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

829:    Level: intermediate

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

835: .keywords: matrix

837: .seealso: MatCreate(), MatCreateIS(), MatMPIAIJSetPreallocation(), MatISGetLocalMat(), MATIS
838: @*/
839: PetscErrorCode  MatISSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
840: {

846:   PetscTryMethod(B,"MatISSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz));
847:   return(0);
848: }

852: static PetscErrorCode  MatISSetPreallocation_IS(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
853: {
854:   Mat_IS         *matis = (Mat_IS*)(B->data);
855:   PetscInt       bs,i,nlocalcols;

859:   if (!matis->A) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"You should first call MatSetLocalToGlobalMapping");
860:   if (!matis->sf) { /* setup SF if not yet created and allocate rootdata and leafdata */
861:     MatISComputeSF_Private(B);
862:   }
863:   if (!d_nnz) {
864:     for (i=0;i<matis->sf_nroots;i++) matis->sf_rootdata[i] = d_nz;
865:   } else {
866:     for (i=0;i<matis->sf_nroots;i++) matis->sf_rootdata[i] = d_nnz[i];
867:   }
868:   if (!o_nnz) {
869:     for (i=0;i<matis->sf_nroots;i++) matis->sf_rootdata[i] += o_nz;
870:   } else {
871:     for (i=0;i<matis->sf_nroots;i++) matis->sf_rootdata[i] += o_nnz[i];
872:   }
873:   PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
874:   MatGetSize(matis->A,NULL,&nlocalcols);
875:   MatGetBlockSize(matis->A,&bs);
876:   PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
877:   for (i=0;i<matis->sf_nleaves;i++) {
878:     matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i],nlocalcols);
879:   }
880:   MatSeqAIJSetPreallocation(matis->A,0,matis->sf_leafdata);
881:   for (i=0;i<matis->sf_nleaves/bs;i++) {
882:     matis->sf_leafdata[i] = matis->sf_leafdata[i*bs]/bs;
883:   }
884:   MatSeqBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata);
885:   for (i=0;i<matis->sf_nleaves/bs;i++) {
886:     matis->sf_leafdata[i] = matis->sf_leafdata[i]-i;
887:   }
888:   MatSeqSBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata);
889:   return(0);
890: }

894: PETSC_EXTERN PetscErrorCode MatISSetMPIXAIJPreallocation_Private(Mat A, Mat B, PetscBool maxreduce)
895: {
896:   Mat_IS          *matis = (Mat_IS*)(A->data);
897:   PetscInt        *my_dnz,*my_onz,*dnz,*onz,*mat_ranges,*row_ownership;
898:   const PetscInt  *global_indices_r,*global_indices_c;
899:   PetscInt        i,j,bs,rows,cols;
900:   PetscInt        lrows,lcols;
901:   PetscInt        local_rows,local_cols;
902:   PetscMPIInt     nsubdomains;
903:   PetscBool       isdense,issbaij;
904:   PetscErrorCode  ierr;

907:   MPI_Comm_size(PetscObjectComm((PetscObject)A),&nsubdomains);
908:   MatGetSize(A,&rows,&cols);
909:   MatGetBlockSize(A,&bs);
910:   MatGetSize(matis->A,&local_rows,&local_cols);
911:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isdense);
912:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);
913:   ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&global_indices_r);
914:   if (A->rmap->mapping != A->cmap->mapping) {
915:     ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&global_indices_c);
916:   } else {
917:     global_indices_c = global_indices_r;
918:   }

920:   if (issbaij) {
921:     MatGetRowUpperTriangular(matis->A);
922:   }
923:   /*
924:      An SF reduce is needed to sum up properly on shared rows.
925:      Note that generally preallocation is not exact, since it overestimates nonzeros
926:   */
927:   if (!matis->sf) { /* setup SF if not yet created and allocate rootdata and leafdata */
928:     MatISComputeSF_Private(A);
929:   }
930:   MatGetLocalSize(A,&lrows,&lcols);
931:   MatPreallocateInitialize(PetscObjectComm((PetscObject)A),lrows,lcols,dnz,onz);
932:   /* All processes need to compute entire row ownership */
933:   PetscMalloc1(rows,&row_ownership);
934:   MatGetOwnershipRanges(A,(const PetscInt**)&mat_ranges);
935:   for (i=0;i<nsubdomains;i++) {
936:     for (j=mat_ranges[i];j<mat_ranges[i+1];j++) {
937:       row_ownership[j] = i;
938:     }
939:   }
940:   MatGetOwnershipRangesColumn(A,(const PetscInt**)&mat_ranges);

942:   /*
943:      my_dnz and my_onz contains exact contribution to preallocation from each local mat
944:      then, they will be summed up properly. This way, preallocation is always sufficient
945:   */
946:   PetscCalloc2(local_rows,&my_dnz,local_rows,&my_onz);
947:   /* preallocation as a MATAIJ */
948:   if (isdense) { /* special case for dense local matrices */
949:     for (i=0;i<local_rows;i++) {
950:       PetscInt index_row = global_indices_r[i];
951:       for (j=i;j<local_rows;j++) {
952:         PetscInt owner = row_ownership[index_row];
953:         PetscInt index_col = global_indices_c[j];
954:         if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */
955:           my_dnz[i] += 1;
956:         } else { /* offdiag block */
957:           my_onz[i] += 1;
958:         }
959:         /* same as before, interchanging rows and cols */
960:         if (i != j) {
961:           owner = row_ownership[index_col];
962:           if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) {
963:             my_dnz[j] += 1;
964:           } else {
965:             my_onz[j] += 1;
966:           }
967:         }
968:       }
969:     }
970:   } else if (matis->A->ops->getrowij) {
971:     const PetscInt *ii,*jj,*jptr;
972:     PetscBool      done;
973:     MatGetRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done);
974:     if (!done) SETERRQ1(PetscObjectComm((PetscObject)(matis->A)),PETSC_ERR_PLIB,"Error in MatGetRowIJ called from %s\n",__FUNCT__);
975:     jptr = jj;
976:     for (i=0;i<local_rows;i++) {
977:       PetscInt index_row = global_indices_r[i];
978:       for (j=0;j<ii[i+1]-ii[i];j++,jptr++) {
979:         PetscInt owner = row_ownership[index_row];
980:         PetscInt index_col = global_indices_c[*jptr];
981:         if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */
982:           my_dnz[i] += 1;
983:         } else { /* offdiag block */
984:           my_onz[i] += 1;
985:         }
986:         /* same as before, interchanging rows and cols */
987:         if (issbaij && index_col != index_row) {
988:           owner = row_ownership[index_col];
989:           if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) {
990:             my_dnz[*jptr] += 1;
991:           } else {
992:             my_onz[*jptr] += 1;
993:           }
994:         }
995:       }
996:     }
997:     MatRestoreRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done);
998:     if (!done) SETERRQ1(PetscObjectComm((PetscObject)(matis->A)),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ called from %s\n",__FUNCT__);
999:   } else { /* loop over rows and use MatGetRow */
1000:     for (i=0;i<local_rows;i++) {
1001:       const PetscInt *cols;
1002:       PetscInt       ncols,index_row = global_indices_r[i];
1003:       MatGetRow(matis->A,i,&ncols,&cols,NULL);
1004:       for (j=0;j<ncols;j++) {
1005:         PetscInt owner = row_ownership[index_row];
1006:         PetscInt index_col = global_indices_c[cols[j]];
1007:         if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */
1008:           my_dnz[i] += 1;
1009:         } else { /* offdiag block */
1010:           my_onz[i] += 1;
1011:         }
1012:         /* same as before, interchanging rows and cols */
1013:         if (issbaij && index_col != index_row) {
1014:           owner = row_ownership[index_col];
1015:           if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) {
1016:             my_dnz[cols[j]] += 1;
1017:           } else {
1018:             my_onz[cols[j]] += 1;
1019:           }
1020:         }
1021:       }
1022:       MatRestoreRow(matis->A,i,&ncols,&cols,NULL);
1023:     }
1024:   }
1025:   ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&global_indices_r);
1026:   if (global_indices_c != global_indices_r) {
1027:     ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&global_indices_c);
1028:   }
1029:   PetscFree(row_ownership);

1031:   /* Reduce my_dnz and my_onz */
1032:   if (maxreduce) {
1033:     PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);
1034:     PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);
1035:     PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);
1036:     PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);
1037:   } else {
1038:     PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);
1039:     PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);
1040:     PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);
1041:     PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);
1042:   }
1043:   PetscFree2(my_dnz,my_onz);

1045:   /* Resize preallocation if overestimated */
1046:   for (i=0;i<lrows;i++) {
1047:     dnz[i] = PetscMin(dnz[i],lcols);
1048:     onz[i] = PetscMin(onz[i],cols-lcols);
1049:   }

1051:   /* Set preallocation */
1052:   MatMPIAIJSetPreallocation(B,0,dnz,0,onz);
1053:   for (i=0;i<lrows/bs;i++) {
1054:     dnz[i] = dnz[i*bs]/bs;
1055:     onz[i] = onz[i*bs]/bs;
1056:   }
1057:   MatMPIBAIJSetPreallocation(B,bs,0,dnz,0,onz);
1058:   MatMPISBAIJSetPreallocation(B,bs,0,dnz,0,onz);
1059:   MatPreallocateFinalize(dnz,onz);
1060:   if (issbaij) {
1061:     MatRestoreRowUpperTriangular(matis->A);
1062:   }
1063:   return(0);
1064: }

1068: static PetscErrorCode MatISGetMPIXAIJ_IS(Mat mat, MatReuse reuse, Mat *M)
1069: {
1070:   Mat_IS         *matis = (Mat_IS*)(mat->data);
1071:   Mat            local_mat;
1072:   /* info on mat */
1073:   PetscInt       bs,rows,cols,lrows,lcols;
1074:   PetscInt       local_rows,local_cols;
1075:   PetscBool      isdense,issbaij,isseqaij;
1076:   PetscMPIInt    nsubdomains;
1077:   /* values insertion */
1078:   PetscScalar    *array;
1079:   /* work */

1083:   /* get info from mat */
1084:   MPI_Comm_size(PetscObjectComm((PetscObject)mat),&nsubdomains);
1085:   if (nsubdomains == 1) {
1086:     Mat            B;
1087:     IS             rows,cols;
1088:     const PetscInt *ridxs,*cidxs;

1090:     ISLocalToGlobalMappingGetIndices(mat->rmap->mapping,&ridxs);
1091:     ISLocalToGlobalMappingGetIndices(mat->cmap->mapping,&cidxs);
1092:     ISCreateGeneral(PETSC_COMM_SELF,mat->rmap->n,ridxs,PETSC_USE_POINTER,&rows);
1093:     ISCreateGeneral(PETSC_COMM_SELF,mat->cmap->n,cidxs,PETSC_USE_POINTER,&cols);
1094:     MatConvert(matis->A,MATAIJ,MAT_INITIAL_MATRIX,&B);
1095:     MatGetSubMatrix(B,rows,cols,reuse,M);
1096:     MatDestroy(&B);
1097:     ISDestroy(&cols);
1098:     ISDestroy(&rows);
1099:     ISLocalToGlobalMappingRestoreIndices(mat->rmap->mapping,&ridxs);
1100:     ISLocalToGlobalMappingRestoreIndices(mat->cmap->mapping,&cidxs);
1101:     return(0);
1102:   }
1103:   MatGetSize(mat,&rows,&cols);
1104:   MatGetBlockSize(mat,&bs);
1105:   MatGetLocalSize(mat,&lrows,&lcols);
1106:   MatGetSize(matis->A,&local_rows,&local_cols);
1107:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isdense);
1108:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQAIJ,&isseqaij);
1109:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);

1111:   if (reuse == MAT_INITIAL_MATRIX) {
1112:     MatType     new_mat_type;
1113:     PetscBool   issbaij_red;

1115:     /* determining new matrix type */
1116:     MPIU_Allreduce(&issbaij,&issbaij_red,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)mat));
1117:     if (issbaij_red) {
1118:       new_mat_type = MATSBAIJ;
1119:     } else {
1120:       if (bs>1) {
1121:         new_mat_type = MATBAIJ;
1122:       } else {
1123:         new_mat_type = MATAIJ;
1124:       }
1125:     }

1127:     MatCreate(PetscObjectComm((PetscObject)mat),M);
1128:     MatSetSizes(*M,lrows,lcols,rows,cols);
1129:     MatSetBlockSize(*M,bs);
1130:     MatSetType(*M,new_mat_type);
1131:     MatISSetMPIXAIJPreallocation_Private(mat,*M,PETSC_FALSE);
1132:   } else {
1133:     PetscInt mbs,mrows,mcols,mlrows,mlcols;
1134:     /* some checks */
1135:     MatGetBlockSize(*M,&mbs);
1136:     MatGetSize(*M,&mrows,&mcols);
1137:     MatGetLocalSize(*M,&mlrows,&mlcols);
1138:     if (mrows != rows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of rows (%d != %d)",rows,mrows);
1139:     if (mcols != cols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of cols (%d != %d)",cols,mcols);
1140:     if (mlrows != lrows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local rows (%d != %d)",lrows,mlrows);
1141:     if (mlcols != lcols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local cols (%d != %d)",lcols,mlcols);
1142:     if (mbs != bs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong block size (%d != %d)",bs,mbs);
1143:     MatZeroEntries(*M);
1144:   }

1146:   if (issbaij) {
1147:     MatConvert(matis->A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&local_mat);
1148:   } else {
1149:     PetscObjectReference((PetscObject)matis->A);
1150:     local_mat = matis->A;
1151:   }

1153:   /* Set values */
1154:   MatSetLocalToGlobalMapping(*M,mat->rmap->mapping,mat->cmap->mapping);
1155:   if (isdense) { /* special case for dense local matrices */
1156:     PetscInt i,*dummy_rows,*dummy_cols;

1158:     if (local_rows != local_cols) {
1159:       PetscMalloc2(local_rows,&dummy_rows,local_cols,&dummy_cols);
1160:       for (i=0;i<local_rows;i++) dummy_rows[i] = i;
1161:       for (i=0;i<local_cols;i++) dummy_cols[i] = i;
1162:     } else {
1163:       PetscMalloc1(local_rows,&dummy_rows);
1164:       for (i=0;i<local_rows;i++) dummy_rows[i] = i;
1165:       dummy_cols = dummy_rows;
1166:     }
1167:     MatSetOption(*M,MAT_ROW_ORIENTED,PETSC_FALSE);
1168:     MatDenseGetArray(local_mat,&array);
1169:     MatSetValuesLocal(*M,local_rows,dummy_rows,local_cols,dummy_cols,array,ADD_VALUES);
1170:     MatDenseRestoreArray(local_mat,&array);
1171:     if (dummy_rows != dummy_cols) {
1172:       PetscFree2(dummy_rows,dummy_cols);
1173:     } else {
1174:       PetscFree(dummy_rows);
1175:     }
1176:   } else if (isseqaij) {
1177:     PetscInt  i,nvtxs,*xadj,*adjncy;
1178:     PetscBool done;

1180:     MatGetRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);
1181:     if (!done) SETERRQ1(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatGetRowIJ called from %s\n",__FUNCT__);
1182:     MatSeqAIJGetArray(local_mat,&array);
1183:     for (i=0;i<nvtxs;i++) {
1184:       MatSetValuesLocal(*M,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],array+xadj[i],ADD_VALUES);
1185:     }
1186:     MatRestoreRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);
1187:     if (!done) SETERRQ1(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ called from %s\n",__FUNCT__);
1188:     MatSeqAIJRestoreArray(local_mat,&array);
1189:   } else { /* very basic values insertion for all other matrix types */
1190:     PetscInt i;

1192:     for (i=0;i<local_rows;i++) {
1193:       PetscInt       j;
1194:       const PetscInt *local_indices_cols;

1196:       MatGetRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);
1197:       MatSetValuesLocal(*M,1,&i,j,local_indices_cols,array,ADD_VALUES);
1198:       MatRestoreRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);
1199:     }
1200:   }
1201:   MatAssemblyBegin(*M,MAT_FINAL_ASSEMBLY);
1202:   MatDestroy(&local_mat);
1203:   MatAssemblyEnd(*M,MAT_FINAL_ASSEMBLY);
1204:   if (isdense) {
1205:     MatSetOption(*M,MAT_ROW_ORIENTED,PETSC_TRUE);
1206:   }
1207:   return(0);
1208: }

1212: /*@
1213:     MatISGetMPIXAIJ - Converts MATIS matrix into a parallel AIJ format

1215:   Input Parameter:
1216: .  mat - the matrix (should be of type MATIS)
1217: .  reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX

1219:   Output Parameter:
1220: .  newmat - the matrix in AIJ format

1222:   Level: developer

1224:   Notes: mat and *newmat cannot be the same object when MAT_REUSE_MATRIX is requested.

1226: .seealso: MATIS
1227: @*/
1228: PetscErrorCode MatISGetMPIXAIJ(Mat mat, MatReuse reuse, Mat *newmat)
1229: {

1236:   if (reuse != MAT_INITIAL_MATRIX) {
1239:     if (mat == *newmat) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse the same matrix");
1240:   }
1241:   PetscUseMethod(mat,"MatISGetMPIXAIJ_C",(Mat,MatReuse,Mat*),(mat,reuse,newmat));
1242:   return(0);
1243: }

1247: PetscErrorCode MatDuplicate_IS(Mat mat,MatDuplicateOption op,Mat *newmat)
1248: {
1250:   Mat_IS         *matis = (Mat_IS*)(mat->data);
1251:   PetscInt       bs,m,n,M,N;
1252:   Mat            B,localmat;

1255:   MatGetBlockSize(mat,&bs);
1256:   MatGetSize(mat,&M,&N);
1257:   MatGetLocalSize(mat,&m,&n);
1258:   MatCreateIS(PetscObjectComm((PetscObject)mat),bs,m,n,M,N,mat->rmap->mapping,mat->cmap->mapping,&B);
1259:   MatDuplicate(matis->A,op,&localmat);
1260:   MatISSetLocalMat(B,localmat);
1261:   MatDestroy(&localmat);
1262:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1263:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1264:   *newmat = B;
1265:   return(0);
1266: }

1270: static PetscErrorCode MatIsHermitian_IS(Mat A,PetscReal tol,PetscBool  *flg)
1271: {
1273:   Mat_IS         *matis = (Mat_IS*)A->data;
1274:   PetscBool      local_sym;

1277:   MatIsHermitian(matis->A,tol,&local_sym);
1278:   MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
1279:   return(0);
1280: }

1284: static PetscErrorCode MatIsSymmetric_IS(Mat A,PetscReal tol,PetscBool  *flg)
1285: {
1287:   Mat_IS         *matis = (Mat_IS*)A->data;
1288:   PetscBool      local_sym;

1291:   MatIsSymmetric(matis->A,tol,&local_sym);
1292:   MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
1293:   return(0);
1294: }

1298: static PetscErrorCode MatDestroy_IS(Mat A)
1299: {
1301:   Mat_IS         *b = (Mat_IS*)A->data;

1304:   MatDestroy(&b->A);
1305:   VecScatterDestroy(&b->cctx);
1306:   VecScatterDestroy(&b->rctx);
1307:   VecDestroy(&b->x);
1308:   VecDestroy(&b->y);
1309:   VecDestroy(&b->counter);
1310:   ISDestroy(&b->getsub_ris);
1311:   ISDestroy(&b->getsub_cis);
1312:   if (b->sf != b->csf) {
1313:     PetscSFDestroy(&b->csf);
1314:     PetscFree2(b->csf_rootdata,b->csf_leafdata);
1315:   }
1316:   PetscSFDestroy(&b->sf);
1317:   PetscFree2(b->sf_rootdata,b->sf_leafdata);
1318:   PetscFree(A->data);
1319:   PetscObjectChangeTypeName((PetscObject)A,0);
1320:   PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL);
1321:   PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",NULL);
1322:   PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",NULL);
1323:   PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",NULL);
1324:   return(0);
1325: }

1329: static PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y)
1330: {
1332:   Mat_IS         *is  = (Mat_IS*)A->data;
1333:   PetscScalar    zero = 0.0;

1336:   /*  scatter the global vector x into the local work vector */
1337:   VecScatterBegin(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
1338:   VecScatterEnd(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);

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

1343:   /* scatter product back into global memory */
1344:   VecSet(y,zero);
1345:   VecScatterBegin(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
1346:   VecScatterEnd(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
1347:   return(0);
1348: }

1352: static PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
1353: {
1354:   Vec            temp_vec;

1358:   if (v3 != v2) {
1359:     MatMult(A,v1,v3);
1360:     VecAXPY(v3,1.0,v2);
1361:   } else {
1362:     VecDuplicate(v2,&temp_vec);
1363:     MatMult(A,v1,temp_vec);
1364:     VecAXPY(temp_vec,1.0,v2);
1365:     VecCopy(temp_vec,v3);
1366:     VecDestroy(&temp_vec);
1367:   }
1368:   return(0);
1369: }

1373: static PetscErrorCode MatMultTranspose_IS(Mat A,Vec y,Vec x)
1374: {
1375:   Mat_IS         *is = (Mat_IS*)A->data;

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

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

1386:   /* scatter product back into global vector */
1387:   VecSet(x,0);
1388:   VecScatterBegin(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);
1389:   VecScatterEnd(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);
1390:   return(0);
1391: }

1395: static PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
1396: {
1397:   Vec            temp_vec;

1401:   if (v3 != v2) {
1402:     MatMultTranspose(A,v1,v3);
1403:     VecAXPY(v3,1.0,v2);
1404:   } else {
1405:     VecDuplicate(v2,&temp_vec);
1406:     MatMultTranspose(A,v1,temp_vec);
1407:     VecAXPY(temp_vec,1.0,v2);
1408:     VecCopy(temp_vec,v3);
1409:     VecDestroy(&temp_vec);
1410:   }
1411:   return(0);
1412: }

1416: static PetscErrorCode MatView_IS(Mat A,PetscViewer viewer)
1417: {
1418:   Mat_IS         *a = (Mat_IS*)A->data;
1420:   PetscViewer    sviewer;

1423:   PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
1424:   MatView(a->A,sviewer);
1425:   PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
1426:   PetscViewerFlush(viewer);
1427:   return(0);
1428: }

1432: static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping)
1433: {
1435:   PetscInt       nr,rbs,nc,cbs;
1436:   Mat_IS         *is = (Mat_IS*)A->data;
1437:   IS             from,to;
1438:   Vec            cglobal,rglobal;

1443:   /* Destroy any previous data */
1444:   VecDestroy(&is->x);
1445:   VecDestroy(&is->y);
1446:   VecDestroy(&is->counter);
1447:   VecScatterDestroy(&is->rctx);
1448:   VecScatterDestroy(&is->cctx);
1449:   MatDestroy(&is->A);
1450:   PetscSFDestroy(&is->sf);
1451:   PetscFree2(is->sf_rootdata,is->sf_leafdata);

1453:   /* Setup Layout and set local to global maps */
1454:   PetscLayoutSetUp(A->rmap);
1455:   PetscLayoutSetUp(A->cmap);
1456:   PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping);
1457:   PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping);

1459:   /* Create the local matrix A */
1460:   ISLocalToGlobalMappingGetSize(rmapping,&nr);
1461:   ISLocalToGlobalMappingGetBlockSize(rmapping,&rbs);
1462:   ISLocalToGlobalMappingGetSize(cmapping,&nc);
1463:   ISLocalToGlobalMappingGetBlockSize(cmapping,&cbs);
1464:   MatCreate(PETSC_COMM_SELF,&is->A);
1465:   MatSetType(is->A,MATAIJ);
1466:   MatSetSizes(is->A,nr,nc,nr,nc);
1467:   MatSetBlockSizes(is->A,rbs,cbs);
1468:   MatSetOptionsPrefix(is->A,((PetscObject)A)->prefix);
1469:   MatAppendOptionsPrefix(is->A,"is_");
1470:   MatSetFromOptions(is->A);

1472:   if (!is->islocalref) { /* setup scatters and local vectors for MatMult */
1473:     /* Create the local work vectors */
1474:     MatCreateVecs(is->A,&is->x,&is->y);

1476:     /* setup the global to local scatters */
1477:     MatCreateVecs(A,&cglobal,&rglobal);
1478:     ISCreateStride(PETSC_COMM_SELF,nr,0,1,&to);
1479:     ISLocalToGlobalMappingApplyIS(rmapping,to,&from);
1480:     VecScatterCreate(rglobal,from,is->y,to,&is->rctx);
1481:     if (rmapping != cmapping) {
1482:       ISDestroy(&to);
1483:       ISDestroy(&from);
1484:       ISCreateStride(PETSC_COMM_SELF,nc,0,1,&to);
1485:       ISLocalToGlobalMappingApplyIS(cmapping,to,&from);
1486:       VecScatterCreate(cglobal,from,is->x,to,&is->cctx);
1487:     } else {
1488:       PetscObjectReference((PetscObject)is->rctx);
1489:       is->cctx = is->rctx;
1490:     }

1492:     /* interface counter vector (local) */
1493:     VecDuplicate(is->y,&is->counter);
1494:     VecSet(is->y,1.);
1495:     VecScatterBegin(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);
1496:     VecScatterEnd(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);
1497:     VecScatterBegin(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);
1498:     VecScatterEnd(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);

1500:     /* free workspace */
1501:     VecDestroy(&rglobal);
1502:     VecDestroy(&cglobal);
1503:     ISDestroy(&to);
1504:     ISDestroy(&from);
1505:   }
1506:   MatSetUp(A);
1507:   return(0);
1508: }

1512: static PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv)
1513: {
1514:   Mat_IS         *is = (Mat_IS*)mat->data;
1516: #if defined(PETSC_USE_DEBUG)
1517:   PetscInt       i,zm,zn;
1518: #endif
1519:   PetscInt       rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION];

1522: #if defined(PETSC_USE_DEBUG)
1523:   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);
1524:   /* count negative indices */
1525:   for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++;
1526:   for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++;
1527: #endif
1528:   ISGlobalToLocalMappingApply(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);
1529:   ISGlobalToLocalMappingApply(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);
1530: #if defined(PETSC_USE_DEBUG)
1531:   /* count negative indices (should be the same as before) */
1532:   for (i=0;i<m;i++) if (rows_l[i] < 0) zm--;
1533:   for (i=0;i<n;i++) if (cols_l[i] < 0) zn--;
1534:   if (zm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the row indices can not be mapped! Maybe you should not use MATIS");
1535:   if (zn) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the column indices can not be mapped! Maybe you should not use MATIS");
1536: #endif
1537:   MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);
1538:   return(0);
1539: }

1543: static PetscErrorCode MatSetValuesBlocked_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv)
1544: {
1545:   Mat_IS         *is = (Mat_IS*)mat->data;
1547: #if defined(PETSC_USE_DEBUG)
1548:   PetscInt       i,zm,zn;
1549: #endif
1550:   PetscInt       rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION];

1553: #if defined(PETSC_USE_DEBUG)
1554:   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);
1555:   /* count negative indices */
1556:   for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++;
1557:   for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++;
1558: #endif
1559:   ISGlobalToLocalMappingApplyBlock(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);
1560:   ISGlobalToLocalMappingApplyBlock(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);
1561: #if defined(PETSC_USE_DEBUG)
1562:   /* count negative indices (should be the same as before) */
1563:   for (i=0;i<m;i++) if (rows_l[i] < 0) zm--;
1564:   for (i=0;i<n;i++) if (cols_l[i] < 0) zn--;
1565:   if (zm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the row indices can not be mapped! Maybe you should not use MATIS");
1566:   if (zn) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the column indices can not be mapped! Maybe you should not use MATIS");
1567: #endif
1568:   MatSetValuesBlocked(is->A,m,rows_l,n,cols_l,values,addv);
1569:   return(0);
1570: }

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

1580:   MatSetValues(is->A,m,rows,n,cols,values,addv);
1581:   return(0);
1582: }

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

1592:   MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv);
1593:   return(0);
1594: }

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

1604:   if (!n) {
1605:     is->pure_neumann = PETSC_TRUE;
1606:   } else {
1607:     PetscInt i;
1608:     is->pure_neumann = PETSC_FALSE;

1610:     if (columns) {
1611:       MatZeroRowsColumns(is->A,n,rows,diag,0,0);
1612:     } else {
1613:       MatZeroRows(is->A,n,rows,diag,0,0);
1614:     }
1615:     if (diag != 0.) {
1616:       const PetscScalar *array;
1617:       VecGetArrayRead(is->counter,&array);
1618:       for (i=0; i<n; i++) {
1619:         MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);
1620:       }
1621:       VecRestoreArrayRead(is->counter,&array);
1622:     }
1623:     MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);
1624:     MatAssemblyEnd(is->A,MAT_FINAL_ASSEMBLY);
1625:   }
1626:   return(0);
1627: }

1631: static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b,PetscBool columns)
1632: {
1633:   Mat_IS         *matis = (Mat_IS*)A->data;
1634:   PetscInt       nr,nl,len,i;
1635:   PetscInt       *lrows;

1639: #if defined(PETSC_USE_DEBUG)
1640:   if (columns || diag != 0. || (x && b)) {
1641:     PetscBool cong;
1642:     PetscLayoutCompare(A->rmap,A->cmap,&cong);
1643:     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 for MATIS");
1644:     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 for MATIS");
1645:     if (!cong && x && b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"A->rmap and A->cmap need to be congruent");
1646:   }
1647: #endif
1648:   /* get locally owned rows */
1649:   PetscLayoutMapLocal_Private(A->rmap,n,rows,&len,&lrows,NULL);
1650:   /* fix right hand side if needed */
1651:   if (x && b) {
1652:     const PetscScalar *xx;
1653:     PetscScalar       *bb;

1655:     VecGetArrayRead(x, &xx);
1656:     VecGetArray(b, &bb);
1657:     for (i=0;i<len;++i) bb[lrows[i]] = diag*xx[lrows[i]];
1658:     VecRestoreArrayRead(x, &xx);
1659:     VecRestoreArray(b, &bb);
1660:   }
1661:   /* get rows associated to the local matrices */
1662:   if (!matis->sf) { /* setup SF if not yet created and allocate rootdata and leafdata */
1663:     MatISComputeSF_Private(A);
1664:   }
1665:   MatGetSize(matis->A,&nl,NULL);
1666:   PetscMemzero(matis->sf_leafdata,nl*sizeof(PetscInt));
1667:   PetscMemzero(matis->sf_rootdata,A->rmap->n*sizeof(PetscInt));
1668:   for (i=0;i<len;i++) matis->sf_rootdata[lrows[i]] = 1;
1669:   PetscFree(lrows);
1670:   PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
1671:   PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
1672:   PetscMalloc1(nl,&lrows);
1673:   for (i=0,nr=0;i<nl;i++) if (matis->sf_leafdata[i]) lrows[nr++] = i;
1674:   MatISZeroRowsColumnsLocal_Private(A,nr,lrows,diag,columns);
1675:   PetscFree(lrows);
1676:   return(0);
1677: }

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

1686:   MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_FALSE);
1687:   return(0);
1688: }

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

1697:   MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_TRUE);
1698:   return(0);
1699: }

1703: static PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type)
1704: {
1705:   Mat_IS         *is = (Mat_IS*)A->data;

1709:   MatAssemblyBegin(is->A,type);
1710:   return(0);
1711: }

1715: static PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type)
1716: {
1717:   Mat_IS         *is = (Mat_IS*)A->data;

1721:   MatAssemblyEnd(is->A,type);
1722:   return(0);
1723: }

1727: static PetscErrorCode MatISGetLocalMat_IS(Mat mat,Mat *local)
1728: {
1729:   Mat_IS *is = (Mat_IS*)mat->data;

1732:   *local = is->A;
1733:   return(0);
1734: }

1738: /*@
1739:     MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix.

1741:   Input Parameter:
1742: .  mat - the matrix

1744:   Output Parameter:
1745: .  local - the local matrix

1747:   Level: advanced

1749:   Notes:
1750:     This can be called if you have precomputed the nonzero structure of the
1751:   matrix and want to provide it to the inner matrix object to improve the performance
1752:   of the MatSetValues() operation.

1754: .seealso: MATIS
1755: @*/
1756: PetscErrorCode MatISGetLocalMat(Mat mat,Mat *local)
1757: {

1763:   PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat*),(mat,local));
1764:   return(0);
1765: }

1769: static PetscErrorCode MatISSetLocalMat_IS(Mat mat,Mat local)
1770: {
1771:   Mat_IS         *is = (Mat_IS*)mat->data;
1772:   PetscInt       nrows,ncols,orows,ocols;

1776:   if (is->A) {
1777:     MatGetSize(is->A,&orows,&ocols);
1778:     MatGetSize(local,&nrows,&ncols);
1779:     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);
1780:   }
1781:   PetscObjectReference((PetscObject)local);
1782:   MatDestroy(&is->A);
1783:   is->A = local;
1784:   return(0);
1785: }

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

1792:   Input Parameter:
1793: .  mat - the matrix
1794: .  local - the local matrix

1796:   Output Parameter:

1798:   Level: advanced

1800:   Notes:
1801:     This can be called if you have precomputed the local matrix and
1802:   want to provide it to the matrix object MATIS.

1804: .seealso: MATIS
1805: @*/
1806: PetscErrorCode MatISSetLocalMat(Mat mat,Mat local)
1807: {

1813:   PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local));
1814:   return(0);
1815: }

1819: static PetscErrorCode MatZeroEntries_IS(Mat A)
1820: {
1821:   Mat_IS         *a = (Mat_IS*)A->data;

1825:   MatZeroEntries(a->A);
1826:   return(0);
1827: }

1831: static PetscErrorCode MatScale_IS(Mat A,PetscScalar a)
1832: {
1833:   Mat_IS         *is = (Mat_IS*)A->data;

1837:   MatScale(is->A,a);
1838:   return(0);
1839: }

1843: static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v)
1844: {
1845:   Mat_IS         *is = (Mat_IS*)A->data;

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

1852:   /* scatter diagonal back into global vector */
1853:   VecSet(v,0);
1854:   VecScatterBegin(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);
1855:   VecScatterEnd(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);
1856:   return(0);
1857: }

1861: static PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg)
1862: {
1863:   Mat_IS         *a = (Mat_IS*)A->data;

1867:   MatSetOption(a->A,op,flg);
1868:   return(0);
1869: }

1873: static PetscErrorCode MatAXPY_IS(Mat Y,PetscScalar a,Mat X,MatStructure str)
1874: {
1875:   Mat_IS         *y = (Mat_IS*)Y->data;
1876:   Mat_IS         *x;
1877: #if defined(PETSC_USE_DEBUG)
1878:   PetscBool      ismatis;
1879: #endif

1883: #if defined(PETSC_USE_DEBUG)
1884:   PetscObjectTypeCompare((PetscObject)X,MATIS,&ismatis);
1885:   if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot call MatAXPY(Y,a,X,str) with X not of type MATIS");
1886: #endif
1887:   x = (Mat_IS*)X->data;
1888:   MatAXPY(y->A,a,x->A,str);
1889:   return(0);
1890: }

1894: static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A,IS row,IS col,Mat *submat)
1895: {
1896:   Mat                    lA;
1897:   Mat_IS                 *matis;
1898:   ISLocalToGlobalMapping rl2g,cl2g;
1899:   IS                     is;
1900:   const PetscInt         *rg,*rl;
1901:   PetscInt               nrg;
1902:   PetscInt               N,M,nrl,i,*idxs;
1903:   PetscErrorCode         ierr;

1906:   ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&rg);
1907:   ISGetLocalSize(row,&nrl);
1908:   ISGetIndices(row,&rl);
1909:   ISLocalToGlobalMappingGetSize(A->rmap->mapping,&nrg);
1910: #if defined(PETSC_USE_DEBUG)
1911:   for (i=0;i<nrl;i++) if (rl[i]>=nrg) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row index %D -> %D greater than maximum possible %D",i,rl[i],nrg);
1912: #endif
1913:   PetscMalloc1(nrg,&idxs);
1914:   /* map from [0,nrl) to row */
1915:   for (i=0;i<nrl;i++) idxs[i] = rl[i];
1916: #if defined(PETSC_USE_DEBUG)
1917:   for (i=nrl;i<nrg;i++) idxs[i] = nrg;
1918: #else
1919:   for (i=nrl;i<nrg;i++) idxs[i] = -1;
1920: #endif
1921:   ISRestoreIndices(row,&rl);
1922:   ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&rg);
1923:   ISCreateGeneral(PetscObjectComm((PetscObject)A),nrg,idxs,PETSC_OWN_POINTER,&is);
1924:   ISLocalToGlobalMappingCreateIS(is,&rl2g);
1925:   ISDestroy(&is);
1926:   /* compute new l2g map for columns */
1927:   if (col != row || A->rmap->mapping != A->cmap->mapping) {
1928:     const PetscInt *cg,*cl;
1929:     PetscInt       ncg;
1930:     PetscInt       ncl;

1932:     ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&cg);
1933:     ISGetLocalSize(col,&ncl);
1934:     ISGetIndices(col,&cl);
1935:     ISLocalToGlobalMappingGetSize(A->cmap->mapping,&ncg);
1936: #if defined(PETSC_USE_DEBUG)
1937:     for (i=0;i<ncl;i++) if (cl[i]>=ncg) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local column index %D -> %D greater than maximum possible %D",i,cl[i],ncg);
1938: #endif
1939:     PetscMalloc1(ncg,&idxs);
1940:     /* map from [0,ncl) to col */
1941:     for (i=0;i<ncl;i++) idxs[i] = cl[i];
1942: #if defined(PETSC_USE_DEBUG)
1943:     for (i=ncl;i<ncg;i++) idxs[i] = ncg;
1944: #else
1945:     for (i=ncl;i<ncg;i++) idxs[i] = -1;
1946: #endif
1947:     ISRestoreIndices(col,&cl);
1948:     ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&cg);
1949:     ISCreateGeneral(PetscObjectComm((PetscObject)A),ncg,idxs,PETSC_OWN_POINTER,&is);
1950:     ISLocalToGlobalMappingCreateIS(is,&cl2g);
1951:     ISDestroy(&is);
1952:   } else {
1953:     PetscObjectReference((PetscObject)rl2g);
1954:     cl2g = rl2g;
1955:   }
1956:   /* create the MATIS submatrix */
1957:   MatGetSize(A,&M,&N);
1958:   MatCreate(PetscObjectComm((PetscObject)A),submat);
1959:   MatSetSizes(*submat,PETSC_DECIDE,PETSC_DECIDE,M,N);
1960:   MatSetType(*submat,MATIS);
1961:   matis = (Mat_IS*)((*submat)->data);
1962:   matis->islocalref = PETSC_TRUE;
1963:   MatSetLocalToGlobalMapping(*submat,rl2g,cl2g);
1964:   MatISGetLocalMat(A,&lA);
1965:   MatISSetLocalMat(*submat,lA);
1966:   MatSetUp(*submat);
1967:   MatAssemblyBegin(*submat,MAT_FINAL_ASSEMBLY);
1968:   MatAssemblyEnd(*submat,MAT_FINAL_ASSEMBLY);
1969:   ISLocalToGlobalMappingDestroy(&rl2g);
1970:   ISLocalToGlobalMappingDestroy(&cl2g);
1971:   /* remove unsupported ops */
1972:   PetscMemzero((*submat)->ops,sizeof(struct _MatOps));
1973:   (*submat)->ops->destroy               = MatDestroy_IS;
1974:   (*submat)->ops->setvalueslocal        = MatSetValuesLocal_SubMat_IS;
1975:   (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS;
1976:   (*submat)->ops->assemblybegin         = MatAssemblyBegin_IS;
1977:   (*submat)->ops->assemblyend           = MatAssemblyEnd_IS;
1978:   return(0);
1979: }

1983: /*@
1984:     MatCreateIS - Creates a "process" unassembled matrix, assembled on each
1985:        process but not across processes.

1987:    Input Parameters:
1988: +     comm    - MPI communicator that will share the matrix
1989: .     bs      - block size of the matrix
1990: .     m,n,M,N - local and/or global sizes of the left and right vector used in matrix vector products
1991: .     rmap    - local to global map for rows
1992: -     cmap    - local to global map for cols

1994:    Output Parameter:
1995: .    A - the resulting matrix

1997:    Level: advanced

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

2004: .seealso: MATIS, MatSetLocalToGlobalMapping()
2005: @*/
2006: PetscErrorCode  MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping rmap,ISLocalToGlobalMapping cmap,Mat *A)
2007: {

2011:   if (!rmap && !cmap) SETERRQ(comm,PETSC_ERR_USER,"You need to provide at least one of the mappings");
2012:   MatCreate(comm,A);
2013:   if (bs > 0) {
2014:     MatSetBlockSize(*A,bs);
2015:   }
2016:   MatSetSizes(*A,m,n,M,N);
2017:   MatSetType(*A,MATIS);
2018:   if (rmap && cmap) {
2019:     MatSetLocalToGlobalMapping(*A,rmap,cmap);
2020:   } else if (!rmap) {
2021:     MatSetLocalToGlobalMapping(*A,cmap,cmap);
2022:   } else {
2023:     MatSetLocalToGlobalMapping(*A,rmap,rmap);
2024:   }
2025:   return(0);
2026: }

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

2034:    Operations Provided:
2035: +  MatMult()
2036: .  MatMultAdd()
2037: .  MatMultTranspose()
2038: .  MatMultTransposeAdd()
2039: .  MatZeroEntries()
2040: .  MatSetOption()
2041: .  MatZeroRows()
2042: .  MatSetValues()
2043: .  MatSetValuesBlocked()
2044: .  MatSetValuesLocal()
2045: .  MatSetValuesBlockedLocal()
2046: .  MatScale()
2047: .  MatGetDiagonal()
2048: .  MatMissingDiagonal()
2049: .  MatDuplicate()
2050: .  MatCopy()
2051: .  MatAXPY()
2052: .  MatGetSubMatrix()
2053: .  MatGetLocalSubMatrix()
2054: .  MatTranspose()
2055: -  MatSetLocalToGlobalMapping()

2057:    Options Database Keys:
2058: . -mat_type is - sets the matrix type to "is" during a call to MatSetFromOptions()

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

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

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

2067:   Level: advanced

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

2071: M*/

2075: PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A)
2076: {
2078:   Mat_IS         *b;

2081:   PetscNewLog(A,&b);
2082:   A->data = (void*)b;

2084:   /* matrix ops */
2085:   PetscMemzero(A->ops,sizeof(struct _MatOps));
2086:   A->ops->mult                    = MatMult_IS;
2087:   A->ops->multadd                 = MatMultAdd_IS;
2088:   A->ops->multtranspose           = MatMultTranspose_IS;
2089:   A->ops->multtransposeadd        = MatMultTransposeAdd_IS;
2090:   A->ops->destroy                 = MatDestroy_IS;
2091:   A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS;
2092:   A->ops->setvalues               = MatSetValues_IS;
2093:   A->ops->setvaluesblocked        = MatSetValuesBlocked_IS;
2094:   A->ops->setvalueslocal          = MatSetValuesLocal_IS;
2095:   A->ops->setvaluesblockedlocal   = MatSetValuesBlockedLocal_IS;
2096:   A->ops->zerorows                = MatZeroRows_IS;
2097:   A->ops->zerorowscolumns         = MatZeroRowsColumns_IS;
2098:   A->ops->assemblybegin           = MatAssemblyBegin_IS;
2099:   A->ops->assemblyend             = MatAssemblyEnd_IS;
2100:   A->ops->view                    = MatView_IS;
2101:   A->ops->zeroentries             = MatZeroEntries_IS;
2102:   A->ops->scale                   = MatScale_IS;
2103:   A->ops->getdiagonal             = MatGetDiagonal_IS;
2104:   A->ops->setoption               = MatSetOption_IS;
2105:   A->ops->ishermitian             = MatIsHermitian_IS;
2106:   A->ops->issymmetric             = MatIsSymmetric_IS;
2107:   A->ops->duplicate               = MatDuplicate_IS;
2108:   A->ops->missingdiagonal         = MatMissingDiagonal_IS;
2109:   A->ops->copy                    = MatCopy_IS;
2110:   A->ops->getlocalsubmatrix       = MatGetLocalSubMatrix_IS;
2111:   A->ops->getsubmatrix            = MatGetSubMatrix_IS;
2112:   A->ops->axpy                    = MatAXPY_IS;
2113:   A->ops->diagonalset             = MatDiagonalSet_IS;
2114:   A->ops->shift                   = MatShift_IS;
2115:   A->ops->transpose               = MatTranspose_IS;
2116:   A->ops->getinfo                 = MatGetInfo_IS;

2118:   /* special MATIS functions */
2119:   PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS);
2120:   PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS);
2121:   PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",MatISGetMPIXAIJ_IS);
2122:   PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",MatISSetPreallocation_IS);
2123:   PetscObjectChangeTypeName((PetscObject)A,MATIS);
2124:   return(0);
2125: }