Actual source code: matis.c

petsc-master 2017-10-20
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 <../src/mat/impls/aij/mpi/mpiaij.h>
 12:  #include <petsc/private/sfimpl.h>

 14: #define MATIS_MAX_ENTRIES_INSERTION 2048

 16: static PetscErrorCode MatISContainerDestroyFields_Private(void *ptr)
 17: {
 18:   MatISLocalFields lf = (MatISLocalFields)ptr;
 19:   PetscInt         i;
 20:   PetscErrorCode   ierr;

 23:   for (i=0;i<lf->nr;i++) {
 24:     ISDestroy(&lf->rf[i]);
 25:   }
 26:   for (i=0;i<lf->nc;i++) {
 27:     ISDestroy(&lf->cf[i]);
 28:   }
 29:   PetscFree2(lf->rf,lf->cf);
 30:   PetscFree(lf);
 31:   return(0);
 32: }

 34: static PetscErrorCode MatISContainerDestroyArray_Private(void *ptr)
 35: {

 39:   PetscFree(ptr);
 40:   return(0);
 41: }

 43: PETSC_INTERN PetscErrorCode MatConvert_MPIAIJ_IS(Mat A,MatType type,MatReuse reuse,Mat *newmat)
 44: {
 45:   Mat_MPIAIJ             *aij  = (Mat_MPIAIJ*)A->data;
 46:   Mat_SeqAIJ             *diag = (Mat_SeqAIJ*)(aij->A->data);
 47:   Mat_SeqAIJ             *offd = (Mat_SeqAIJ*)(aij->B->data);
 48:   Mat                    lA;
 49:   ISLocalToGlobalMapping rl2g,cl2g;
 50:   IS                     is;
 51:   MPI_Comm               comm;
 52:   void                   *ptrs[2];
 53:   const char             *names[2] = {"_convert_csr_aux","_convert_csr_data"};
 54:   PetscScalar            *dd,*od,*aa,*data;
 55:   PetscInt               *di,*dj,*oi,*oj;
 56:   PetscInt               *aux,*ii,*jj;
 57:   PetscInt               lc,dr,dc,oc,str,stc,nnz,i,jd,jo,cum;
 58:   PetscErrorCode         ierr;

 61:   PetscObjectGetComm((PetscObject)A,&comm);
 62:   if (!aij->garray) SETERRQ(comm,PETSC_ERR_SUP,"garray not present");

 64:   /* access relevant information from MPIAIJ */
 65:   MatGetOwnershipRange(A,&str,NULL);
 66:   MatGetOwnershipRangeColumn(A,&stc,NULL);
 67:   MatGetLocalSize(A,&dr,&dc);
 68:   di   = diag->i;
 69:   dj   = diag->j;
 70:   dd   = diag->a;
 71:   oc   = aij->B->cmap->n;
 72:   oi   = offd->i;
 73:   oj   = offd->j;
 74:   od   = offd->a;
 75:   nnz  = diag->i[dr] + offd->i[dr];

 77:   /* generate l2g maps for rows and cols */
 78:   ISCreateStride(comm,dr,str,1,&is);
 79:   ISLocalToGlobalMappingCreateIS(is,&rl2g);
 80:   ISDestroy(&is);
 81:   if (dr) {
 82:     PetscMalloc1(dc+oc,&aux);
 83:     for (i=0; i<dc; i++) aux[i]    = i+stc;
 84:     for (i=0; i<oc; i++) aux[i+dc] = aij->garray[i];
 85:     ISCreateGeneral(comm,dc+oc,aux,PETSC_OWN_POINTER,&is);
 86:     lc   = dc+oc;
 87:   } else {
 88:     ISCreateGeneral(comm,0,NULL,PETSC_OWN_POINTER,&is);
 89:     lc   = 0;
 90:   }
 91:   ISLocalToGlobalMappingCreateIS(is,&cl2g);
 92:   ISDestroy(&is);

 94:   /* create MATIS object */
 95:   MatCreate(comm,newmat);
 96:   MatSetSizes(*newmat,dr,dc,PETSC_DECIDE,PETSC_DECIDE);
 97:   MatSetType(*newmat,MATIS);
 98:   MatSetLocalToGlobalMapping(*newmat,rl2g,cl2g);
 99:   ISLocalToGlobalMappingDestroy(&rl2g);
100:   ISLocalToGlobalMappingDestroy(&cl2g);

102:   /* merge local matrices */
103:   PetscMalloc1(nnz+dr+1,&aux);
104:   PetscMalloc1(nnz,&data);
105:   ii   = aux;
106:   jj   = aux+dr+1;
107:   aa   = data;
108:   *ii  = *(di++) + *(oi++);
109:   for (jd=0,jo=0,cum=0;*ii<nnz;cum++)
110:   {
111:      for (;jd<*di;jd++) { *jj++ = *dj++;      *aa++ = *dd++; }
112:      for (;jo<*oi;jo++) { *jj++ = *oj++ + dc; *aa++ = *od++; }
113:      *(++ii) = *(di++) + *(oi++);
114:   }
115:   for (;cum<dr;cum++) *(++ii) = nnz;
116:   ii   = aux;
117:   jj   = aux+dr+1;
118:   aa   = data;
119:   MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,dr,lc,ii,jj,aa,&lA);

121:   /* create containers to destroy the data */
122:   ptrs[0] = aux;
123:   ptrs[1] = data;
124:   for (i=0; i<2; i++) {
125:     PetscContainer c;

127:     PetscContainerCreate(PETSC_COMM_SELF,&c);
128:     PetscContainerSetPointer(c,ptrs[i]);
129:     PetscContainerSetUserDestroy(c,MatISContainerDestroyArray_Private);
130:     PetscObjectCompose((PetscObject)lA,names[i],(PetscObject)c);
131:     PetscContainerDestroy(&c);
132:   }

134:   /* finalize matrix */
135:   MatISSetLocalMat(*newmat,lA);
136:   MatDestroy(&lA);
137:   MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
138:   MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
139:   return(0);
140: }

142: /*@
143:    MatISSetUpSF - Setup star forest objects used by MatIS.

145:    Collective on MPI_Comm

147:    Input Parameters:
148: +  A - the matrix

150:    Level: advanced

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

154: .keywords: matrix

156: .seealso: MatCreate(), MatCreateIS(), MatISSetPreallocation(), MatISGetLocalMat()
157: @*/
158: PetscErrorCode  MatISSetUpSF(Mat A)
159: {

165:   PetscTryMethod(A,"MatISSetUpSF_C",(Mat),(A));
166:   return(0);
167: }

169: PETSC_INTERN PetscErrorCode MatConvert_Nest_IS(Mat A,MatType type,MatReuse reuse,Mat *newmat)
170: {
171:   Mat                    **nest,*snest,**rnest,lA,B;
172:   IS                     *iscol,*isrow,*islrow,*islcol;
173:   ISLocalToGlobalMapping rl2g,cl2g;
174:   MPI_Comm               comm;
175:   PetscInt               *lr,*lc,*l2gidxs;
176:   PetscInt               i,j,nr,nc,rbs,cbs;
177:   PetscBool              convert,lreuse,*istrans;
178:   PetscErrorCode         ierr;

181:   MatNestGetSubMats(A,&nr,&nc,&nest);
182:   lreuse = PETSC_FALSE;
183:   rnest  = NULL;
184:   if (reuse == MAT_REUSE_MATRIX) {
185:     PetscBool ismatis,isnest;

187:     PetscObjectTypeCompare((PetscObject)*newmat,MATIS,&ismatis);
188:     if (!ismatis) SETERRQ1(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_USER,"Cannot reuse matrix of type %s",((PetscObject)(*newmat))->type);
189:     MatISGetLocalMat(*newmat,&lA);
190:     PetscObjectTypeCompare((PetscObject)lA,MATNEST,&isnest);
191:     if (isnest) {
192:       MatNestGetSubMats(lA,&i,&j,&rnest);
193:       lreuse = (PetscBool)(i == nr && j == nc);
194:       if (!lreuse) rnest = NULL;
195:     }
196:   }
197:   PetscObjectGetComm((PetscObject)A,&comm);
198:   PetscCalloc2(nr,&lr,nc,&lc);
199:   PetscCalloc6(nr,&isrow,nc,&iscol,
200:                       nr,&islrow,nc,&islcol,
201:                       nr*nc,&snest,nr*nc,&istrans);
202:   MatNestGetISs(A,isrow,iscol);
203:   for (i=0;i<nr;i++) {
204:     for (j=0;j<nc;j++) {
205:       PetscBool ismatis;
206:       PetscInt  l1,l2,lb1,lb2,ij=i*nc+j;

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

211:       /* Nested matrices should be of type MATIS */
212:       PetscObjectTypeCompare((PetscObject)nest[i][j],MATTRANSPOSEMAT,&istrans[ij]);
213:       if (istrans[ij]) {
214:         Mat T,lT;
215:         MatTransposeGetMat(nest[i][j],&T);
216:         PetscObjectTypeCompare((PetscObject)T,MATIS,&ismatis);
217:         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);
218:         MatISGetLocalMat(T,&lT);
219:         MatCreateTranspose(lT,&snest[ij]);
220:       } else {
221:         PetscObjectTypeCompare((PetscObject)nest[i][j],MATIS,&ismatis);
222:         if (!ismatis) SETERRQ2(comm,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) is not of type MATIS",i,j);
223:         MatISGetLocalMat(nest[i][j],&snest[ij]);
224:       }

226:       /* Check compatibility of local sizes */
227:       MatGetSize(snest[ij],&l1,&l2);
228:       MatGetBlockSizes(snest[ij],&lb1,&lb2);
229:       if (!l1 || !l2) continue;
230:       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);
231:       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);
232:       lr[i] = l1;
233:       lc[j] = l2;

235:       /* check compatibilty for local matrix reusage */
236:       if (rnest && !rnest[i][j] != !snest[ij]) lreuse = PETSC_FALSE;
237:     }
238:   }

240: #if defined (PETSC_USE_DEBUG)
241:   /* Check compatibility of l2g maps for rows */
242:   for (i=0;i<nr;i++) {
243:     rl2g = NULL;
244:     for (j=0;j<nc;j++) {
245:       PetscInt n1,n2;

247:       if (!nest[i][j]) continue;
248:       if (istrans[i*nc+j]) {
249:         Mat T;

251:         MatTransposeGetMat(nest[i][j],&T);
252:         MatGetLocalToGlobalMapping(T,NULL,&cl2g);
253:       } else {
254:         MatGetLocalToGlobalMapping(nest[i][j],&cl2g,NULL);
255:       }
256:       ISLocalToGlobalMappingGetSize(cl2g,&n1);
257:       if (!n1) continue;
258:       if (!rl2g) {
259:         rl2g = cl2g;
260:       } else {
261:         const PetscInt *idxs1,*idxs2;
262:         PetscBool      same;

264:         ISLocalToGlobalMappingGetSize(rl2g,&n2);
265:         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);
266:         ISLocalToGlobalMappingGetIndices(cl2g,&idxs1);
267:         ISLocalToGlobalMappingGetIndices(rl2g,&idxs2);
268:         PetscMemcmp(idxs1,idxs2,n1*sizeof(PetscInt),&same);
269:         ISLocalToGlobalMappingRestoreIndices(cl2g,&idxs1);
270:         ISLocalToGlobalMappingRestoreIndices(rl2g,&idxs2);
271:         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);
272:       }
273:     }
274:   }
275:   /* Check compatibility of l2g maps for columns */
276:   for (i=0;i<nc;i++) {
277:     rl2g = NULL;
278:     for (j=0;j<nr;j++) {
279:       PetscInt n1,n2;

281:       if (!nest[j][i]) continue;
282:       if (istrans[j*nc+i]) {
283:         Mat T;

285:         MatTransposeGetMat(nest[j][i],&T);
286:         MatGetLocalToGlobalMapping(T,&cl2g,NULL);
287:       } else {
288:         MatGetLocalToGlobalMapping(nest[j][i],NULL,&cl2g);
289:       }
290:       ISLocalToGlobalMappingGetSize(cl2g,&n1);
291:       if (!n1) continue;
292:       if (!rl2g) {
293:         rl2g = cl2g;
294:       } else {
295:         const PetscInt *idxs1,*idxs2;
296:         PetscBool      same;

298:         ISLocalToGlobalMappingGetSize(rl2g,&n2);
299:         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);
300:         ISLocalToGlobalMappingGetIndices(cl2g,&idxs1);
301:         ISLocalToGlobalMappingGetIndices(rl2g,&idxs2);
302:         PetscMemcmp(idxs1,idxs2,n1*sizeof(PetscInt),&same);
303:         ISLocalToGlobalMappingRestoreIndices(cl2g,&idxs1);
304:         ISLocalToGlobalMappingRestoreIndices(rl2g,&idxs2);
305:         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);
306:       }
307:     }
308:   }
309: #endif

311:   B = NULL;
312:   if (reuse != MAT_REUSE_MATRIX) {
313:     PetscInt stl;

315:     /* Create l2g map for the rows of the new matrix and index sets for the local MATNEST */
316:     for (i=0,stl=0;i<nr;i++) stl += lr[i];
317:     PetscMalloc1(stl,&l2gidxs);
318:     for (i=0,stl=0;i<nr;i++) {
319:       Mat            usedmat;
320:       Mat_IS         *matis;
321:       const PetscInt *idxs;

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

326:       /* l2gmap */
327:       j = 0;
328:       usedmat = nest[i][j];
329:       while (!usedmat && j < nc-1) usedmat = nest[i][++j];
330:       if (!usedmat) SETERRQ(comm,PETSC_ERR_SUP,"Cannot find valid row mat");

332:       if (istrans[i*nc+j]) {
333:         Mat T;
334:         MatTransposeGetMat(usedmat,&T);
335:         usedmat = T;
336:       }
337:       MatISSetUpSF(usedmat);
338:       matis = (Mat_IS*)(usedmat->data);
339:       ISGetIndices(isrow[i],&idxs);
340:       if (istrans[i*nc+j]) {
341:         PetscSFBcastBegin(matis->csf,MPIU_INT,idxs,l2gidxs+stl);
342:         PetscSFBcastEnd(matis->csf,MPIU_INT,idxs,l2gidxs+stl);
343:       } else {
344:         PetscSFBcastBegin(matis->sf,MPIU_INT,idxs,l2gidxs+stl);
345:         PetscSFBcastEnd(matis->sf,MPIU_INT,idxs,l2gidxs+stl);
346:       }
347:       ISRestoreIndices(isrow[i],&idxs);
348:       stl += lr[i];
349:     }
350:     ISLocalToGlobalMappingCreate(comm,1,stl,l2gidxs,PETSC_OWN_POINTER,&rl2g);

352:     /* Create l2g map for columns of the new matrix and index sets for the local MATNEST */
353:     for (i=0,stl=0;i<nc;i++) stl += lc[i];
354:     PetscMalloc1(stl,&l2gidxs);
355:     for (i=0,stl=0;i<nc;i++) {
356:       Mat            usedmat;
357:       Mat_IS         *matis;
358:       const PetscInt *idxs;

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

363:       /* l2gmap */
364:       j = 0;
365:       usedmat = nest[j][i];
366:       while (!usedmat && j < nr-1) usedmat = nest[++j][i];
367:       if (!usedmat) SETERRQ(comm,PETSC_ERR_SUP,"Cannot find valid column mat");
368:       if (istrans[j*nc+i]) {
369:         Mat T;
370:         MatTransposeGetMat(usedmat,&T);
371:         usedmat = T;
372:       }
373:       MatISSetUpSF(usedmat);
374:       matis = (Mat_IS*)(usedmat->data);
375:       ISGetIndices(iscol[i],&idxs);
376:       if (istrans[j*nc+i]) {
377:         PetscSFBcastBegin(matis->sf,MPIU_INT,idxs,l2gidxs+stl);
378:         PetscSFBcastEnd(matis->sf,MPIU_INT,idxs,l2gidxs+stl);
379:       } else {
380:         PetscSFBcastBegin(matis->csf,MPIU_INT,idxs,l2gidxs+stl);
381:         PetscSFBcastEnd(matis->csf,MPIU_INT,idxs,l2gidxs+stl);
382:       }
383:       ISRestoreIndices(iscol[i],&idxs);
384:       stl += lc[i];
385:     }
386:     ISLocalToGlobalMappingCreate(comm,1,stl,l2gidxs,PETSC_OWN_POINTER,&cl2g);

388:     /* Create MATIS */
389:     MatCreate(comm,&B);
390:     MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
391:     MatGetBlockSizes(A,&rbs,&cbs);
392:     MatSetBlockSizes(B,rbs,cbs);
393:     MatSetType(B,MATIS);
394:     MatSetLocalToGlobalMapping(B,rl2g,cl2g);
395:     ISLocalToGlobalMappingDestroy(&rl2g);
396:     ISLocalToGlobalMappingDestroy(&cl2g);
397:     MatCreateNest(PETSC_COMM_SELF,nr,islrow,nc,islcol,snest,&lA);
398:     for (i=0;i<nr*nc;i++) {
399:       if (istrans[i]) {
400:         MatDestroy(&snest[i]);
401:       }
402:     }
403:     MatISSetLocalMat(B,lA);
404:     MatDestroy(&lA);
405:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
406:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
407:     if (reuse == MAT_INPLACE_MATRIX) {
408:       MatHeaderReplace(A,&B);
409:     } else {
410:       *newmat = B;
411:     }
412:   } else {
413:     if (lreuse) {
414:       MatISGetLocalMat(*newmat,&lA);
415:       for (i=0;i<nr;i++) {
416:         for (j=0;j<nc;j++) {
417:           if (snest[i*nc+j]) {
418:             MatNestSetSubMat(lA,i,j,snest[i*nc+j]);
419:             if (istrans[i*nc+j]) {
420:               MatDestroy(&snest[i*nc+j]);
421:             }
422:           }
423:         }
424:       }
425:     } else {
426:       PetscInt stl;
427:       for (i=0,stl=0;i<nr;i++) {
428:         ISCreateStride(PETSC_COMM_SELF,lr[i],stl,1,&islrow[i]);
429:         stl  += lr[i];
430:       }
431:       for (i=0,stl=0;i<nc;i++) {
432:         ISCreateStride(PETSC_COMM_SELF,lc[i],stl,1,&islcol[i]);
433:         stl  += lc[i];
434:       }
435:       MatCreateNest(PETSC_COMM_SELF,nr,islrow,nc,islcol,snest,&lA);
436:       if (istrans[i]) {
437:         MatDestroy(&snest[i]);
438:       }
439:       MatISSetLocalMat(*newmat,lA);
440:       MatDestroy(&lA);
441:     }
442:     MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
443:     MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
444:   }

446:   /* Create local matrix in MATNEST format */
447:   convert = PETSC_FALSE;
448:   PetscOptionsGetBool(NULL,((PetscObject)A)->prefix,"-matis_convert_local_nest",&convert,NULL);
449:   if (convert) {
450:     Mat              M;
451:     MatISLocalFields lf;
452:     PetscContainer   c;

454:     MatISGetLocalMat(*newmat,&lA);
455:     MatConvert(lA,MATAIJ,MAT_INITIAL_MATRIX,&M);
456:     MatISSetLocalMat(*newmat,M);
457:     MatDestroy(&M);

459:     /* attach local fields to the matrix */
460:     PetscNew(&lf);
461:     PetscCalloc2(nr,&lf->rf,nc,&lf->cf);
462:     for (i=0;i<nr;i++) {
463:       PetscInt n,st;

465:       ISGetLocalSize(islrow[i],&n);
466:       ISStrideGetInfo(islrow[i],&st,NULL);
467:       ISCreateStride(comm,n,st,1,&lf->rf[i]);
468:     }
469:     for (i=0;i<nc;i++) {
470:       PetscInt n,st;

472:       ISGetLocalSize(islcol[i],&n);
473:       ISStrideGetInfo(islcol[i],&st,NULL);
474:       ISCreateStride(comm,n,st,1,&lf->cf[i]);
475:     }
476:     lf->nr = nr;
477:     lf->nc = nc;
478:     PetscContainerCreate(PetscObjectComm((PetscObject)(*newmat)),&c);
479:     PetscContainerSetPointer(c,lf);
480:     PetscContainerSetUserDestroy(c,MatISContainerDestroyFields_Private);
481:     PetscObjectCompose((PetscObject)(*newmat),"_convert_nest_lfields",(PetscObject)c);
482:     PetscContainerDestroy(&c);
483:   }

485:   /* Free workspace */
486:   for (i=0;i<nr;i++) {
487:     ISDestroy(&islrow[i]);
488:   }
489:   for (i=0;i<nc;i++) {
490:     ISDestroy(&islcol[i]);
491:   }
492:   PetscFree6(isrow,iscol,islrow,islcol,snest,istrans);
493:   PetscFree2(lr,lc);
494:   return(0);
495: }

497: static PetscErrorCode MatDiagonalScale_IS(Mat A, Vec l, Vec r)
498: {
499:   Mat_IS            *matis = (Mat_IS*)A->data;
500:   Vec               ll,rr;
501:   const PetscScalar *Y,*X;
502:   PetscScalar       *x,*y;
503:   PetscErrorCode    ierr;

506:   MatISSetUpSF(A);
507:   if (l) {
508:     ll   = matis->y;
509:     VecGetArrayRead(l,&Y);
510:     VecGetArray(ll,&y);
511:     PetscSFBcastBegin(matis->sf,MPIU_SCALAR,Y,y);
512:   } else {
513:     ll = NULL;
514:   }
515:   if (r) {
516:     rr   = matis->x;
517:     VecGetArrayRead(r,&X);
518:     VecGetArray(rr,&x);
519:     PetscSFBcastBegin(matis->csf,MPIU_SCALAR,X,x);
520:   } else {
521:     rr = NULL;
522:   }
523:   if (ll) {
524:     PetscSFBcastEnd(matis->sf,MPIU_SCALAR,Y,y);
525:     VecRestoreArrayRead(l,&Y);
526:     VecRestoreArray(ll,&y);
527:   }
528:   if (rr) {
529:     PetscSFBcastEnd(matis->csf,MPIU_SCALAR,X,x);
530:     VecRestoreArrayRead(r,&X);
531:     VecRestoreArray(rr,&x);
532:   }
533:   MatDiagonalScale(matis->A,ll,rr);
534:   return(0);
535: }

537: static PetscErrorCode MatGetInfo_IS(Mat A,MatInfoType flag,MatInfo *ginfo)
538: {
539:   Mat_IS         *matis = (Mat_IS*)A->data;
540:   MatInfo        info;
541:   PetscReal      isend[6],irecv[6];
542:   PetscInt       bs;

546:   MatGetBlockSize(A,&bs);
547:   if (matis->A->ops->getinfo) {
548:     MatGetInfo(matis->A,MAT_LOCAL,&info);
549:     isend[0] = info.nz_used;
550:     isend[1] = info.nz_allocated;
551:     isend[2] = info.nz_unneeded;
552:     isend[3] = info.memory;
553:     isend[4] = info.mallocs;
554:   } else {
555:     isend[0] = 0.;
556:     isend[1] = 0.;
557:     isend[2] = 0.;
558:     isend[3] = 0.;
559:     isend[4] = 0.;
560:   }
561:   isend[5] = matis->A->num_ass;
562:   if (flag == MAT_LOCAL) {
563:     ginfo->nz_used      = isend[0];
564:     ginfo->nz_allocated = isend[1];
565:     ginfo->nz_unneeded  = isend[2];
566:     ginfo->memory       = isend[3];
567:     ginfo->mallocs      = isend[4];
568:     ginfo->assemblies   = isend[5];
569:   } else if (flag == MAT_GLOBAL_MAX) {
570:     MPIU_Allreduce(isend,irecv,6,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)A));

572:     ginfo->nz_used      = irecv[0];
573:     ginfo->nz_allocated = irecv[1];
574:     ginfo->nz_unneeded  = irecv[2];
575:     ginfo->memory       = irecv[3];
576:     ginfo->mallocs      = irecv[4];
577:     ginfo->assemblies   = irecv[5];
578:   } else if (flag == MAT_GLOBAL_SUM) {
579:     MPIU_Allreduce(isend,irecv,5,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)A));

581:     ginfo->nz_used      = irecv[0];
582:     ginfo->nz_allocated = irecv[1];
583:     ginfo->nz_unneeded  = irecv[2];
584:     ginfo->memory       = irecv[3];
585:     ginfo->mallocs      = irecv[4];
586:     ginfo->assemblies   = A->num_ass;
587:   }
588:   ginfo->block_size        = bs;
589:   ginfo->fill_ratio_given  = 0;
590:   ginfo->fill_ratio_needed = 0;
591:   ginfo->factor_mallocs    = 0;
592:   return(0);
593: }

595: PetscErrorCode MatTranspose_IS(Mat A,MatReuse reuse,Mat *B)
596: {
597:   Mat                    C,lC,lA;
598:   PetscErrorCode         ierr;

601:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
602:     ISLocalToGlobalMapping rl2g,cl2g;
603:     MatCreate(PetscObjectComm((PetscObject)A),&C);
604:     MatSetSizes(C,A->cmap->n,A->rmap->n,A->cmap->N,A->rmap->N);
605:     MatSetBlockSizes(C,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));
606:     MatSetType(C,MATIS);
607:     MatGetLocalToGlobalMapping(A,&rl2g,&cl2g);
608:     MatSetLocalToGlobalMapping(C,cl2g,rl2g);
609:   } else {
610:     C = *B;
611:   }

613:   /* perform local transposition */
614:   MatISGetLocalMat(A,&lA);
615:   MatTranspose(lA,MAT_INITIAL_MATRIX,&lC);
616:   MatISSetLocalMat(C,lC);
617:   MatDestroy(&lC);

619:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
620:     *B = C;
621:   } else {
622:     MatHeaderMerge(A,&C);
623:   }
624:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
625:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
626:   return(0);
627: }

629: PetscErrorCode  MatDiagonalSet_IS(Mat A,Vec D,InsertMode insmode)
630: {
631:   Mat_IS         *is = (Mat_IS*)A->data;

635:   if (!D) { /* this code branch is used by MatShift_IS */
636:     VecSet(is->y,1.);
637:   } else {
638:     VecScatterBegin(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD);
639:     VecScatterEnd(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD);
640:   }
641:   VecPointwiseDivide(is->y,is->y,is->counter);
642:   MatDiagonalSet(is->A,is->y,insmode);
643:   return(0);
644: }

646: PetscErrorCode  MatShift_IS(Mat A,PetscScalar a)
647: {

651:   MatDiagonalSet_IS(A,NULL,ADD_VALUES);
652:   return(0);
653: }

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

662: #if defined(PETSC_USE_DEBUG)
663:   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);
664: #endif
665:   ISLocalToGlobalMappingApply(A->rmap->mapping,m,rows,rows_l);
666:   ISLocalToGlobalMappingApply(A->cmap->mapping,n,cols,cols_l);
667:   MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);
668:   return(0);
669: }

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

678: #if defined(PETSC_USE_DEBUG)
679:   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);
680: #endif
681:   ISLocalToGlobalMappingApplyBlock(A->rmap->mapping,m,rows,rows_l);
682:   ISLocalToGlobalMappingApplyBlock(A->cmap->mapping,n,cols,cols_l);
683:   MatSetValuesBlocked(is->A,m,rows_l,n,cols_l,values,addv);
684:   return(0);
685: }

687: static PetscErrorCode PetscLayoutMapLocal_Private(PetscLayout map,PetscInt N,const PetscInt idxs[], PetscInt *on,PetscInt **oidxs,PetscInt **ogidxs)
688: {
689:   PetscInt      *owners = map->range;
690:   PetscInt       n      = map->n;
691:   PetscSF        sf;
692:   PetscInt      *lidxs,*work = NULL;
693:   PetscSFNode   *ridxs;
694:   PetscMPIInt    rank;
695:   PetscInt       r, p = 0, len = 0;

699:   if (on) *on = 0;              /* squelch -Wmaybe-uninitialized */
700:   /* Create SF where leaves are input idxs and roots are owned idxs (code adapted from MatZeroRowsMapLocal_Private) */
701:   MPI_Comm_rank(map->comm,&rank);
702:   PetscMalloc1(n,&lidxs);
703:   for (r = 0; r < n; ++r) lidxs[r] = -1;
704:   PetscMalloc1(N,&ridxs);
705:   for (r = 0; r < N; ++r) {
706:     const PetscInt idx = idxs[r];
707:     if (idx < 0 || map->N <= idx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Index %D out of range [0,%D)",idx,map->N);
708:     if (idx < owners[p] || owners[p+1] <= idx) { /* short-circuit the search if the last p owns this idx too */
709:       PetscLayoutFindOwner(map,idx,&p);
710:     }
711:     ridxs[r].rank = p;
712:     ridxs[r].index = idxs[r] - owners[p];
713:   }
714:   PetscSFCreate(map->comm,&sf);
715:   PetscSFSetGraph(sf,n,N,NULL,PETSC_OWN_POINTER,ridxs,PETSC_OWN_POINTER);
716:   PetscSFReduceBegin(sf,MPIU_INT,(PetscInt*)idxs,lidxs,MPI_LOR);
717:   PetscSFReduceEnd(sf,MPIU_INT,(PetscInt*)idxs,lidxs,MPI_LOR);
718:   if (ogidxs) { /* communicate global idxs */
719:     PetscInt cum = 0,start,*work2;

721:     PetscMalloc1(n,&work);
722:     PetscCalloc1(N,&work2);
723:     for (r = 0; r < N; ++r) if (idxs[r] >=0) cum++;
724:     MPI_Scan(&cum,&start,1,MPIU_INT,MPI_SUM,map->comm);
725:     start -= cum;
726:     cum = 0;
727:     for (r = 0; r < N; ++r) if (idxs[r] >=0) work2[r] = start+cum++;
728:     PetscSFReduceBegin(sf,MPIU_INT,work2,work,MPIU_REPLACE);
729:     PetscSFReduceEnd(sf,MPIU_INT,work2,work,MPIU_REPLACE);
730:     PetscFree(work2);
731:   }
732:   PetscSFDestroy(&sf);
733:   /* Compress and put in indices */
734:   for (r = 0; r < n; ++r)
735:     if (lidxs[r] >= 0) {
736:       if (work) work[len] = work[r];
737:       lidxs[len++] = r;
738:     }
739:   if (on) *on = len;
740:   if (oidxs) *oidxs = lidxs;
741:   if (ogidxs) *ogidxs = work;
742:   return(0);
743: }

745: static PetscErrorCode MatCreateSubMatrix_IS(Mat mat,IS irow,IS icol,MatReuse scall,Mat *newmat)
746: {
747:   Mat               locmat,newlocmat;
748:   Mat_IS            *newmatis;
749: #if defined(PETSC_USE_DEBUG)
750:   Vec               rtest,ltest;
751:   const PetscScalar *array;
752: #endif
753:   const PetscInt    *idxs;
754:   PetscInt          i,m,n;
755:   PetscErrorCode    ierr;

758:   if (scall == MAT_REUSE_MATRIX) {
759:     PetscBool ismatis;

761:     PetscObjectTypeCompare((PetscObject)*newmat,MATIS,&ismatis);
762:     if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Not of MATIS type");
763:     newmatis = (Mat_IS*)(*newmat)->data;
764:     if (!newmatis->getsub_ris) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Misses local row IS");
765:     if (!newmatis->getsub_cis) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Misses local col IS");
766:   }
767:   /* irow and icol may not have duplicate entries */
768: #if defined(PETSC_USE_DEBUG)
769:   MatCreateVecs(mat,&ltest,&rtest);
770:   ISGetLocalSize(irow,&n);
771:   ISGetIndices(irow,&idxs);
772:   for (i=0;i<n;i++) {
773:     VecSetValue(rtest,idxs[i],1.0,ADD_VALUES);
774:   }
775:   VecAssemblyBegin(rtest);
776:   VecAssemblyEnd(rtest);
777:   VecGetLocalSize(rtest,&n);
778:   VecGetOwnershipRange(rtest,&m,NULL);
779:   VecGetArrayRead(rtest,&array);
780:   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]));
781:   VecRestoreArrayRead(rtest,&array);
782:   ISRestoreIndices(irow,&idxs);
783:   ISGetLocalSize(icol,&n);
784:   ISGetIndices(icol,&idxs);
785:   for (i=0;i<n;i++) {
786:     VecSetValue(ltest,idxs[i],1.0,ADD_VALUES);
787:   }
788:   VecAssemblyBegin(ltest);
789:   VecAssemblyEnd(ltest);
790:   VecGetLocalSize(ltest,&n);
791:   VecGetOwnershipRange(ltest,&m,NULL);
792:   VecGetArrayRead(ltest,&array);
793:   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]));
794:   VecRestoreArrayRead(ltest,&array);
795:   ISRestoreIndices(icol,&idxs);
796:   VecDestroy(&rtest);
797:   VecDestroy(&ltest);
798: #endif
799:   if (scall == MAT_INITIAL_MATRIX) {
800:     Mat_IS                 *matis = (Mat_IS*)mat->data;
801:     ISLocalToGlobalMapping rl2g;
802:     IS                     is;
803:     PetscInt               *lidxs,*lgidxs,*newgidxs;
804:     PetscInt               ll,newloc;
805:     MPI_Comm               comm;

807:     PetscObjectGetComm((PetscObject)mat,&comm);
808:     ISGetLocalSize(irow,&m);
809:     ISGetLocalSize(icol,&n);
810:     MatCreate(comm,newmat);
811:     MatSetType(*newmat,MATIS);
812:     MatSetSizes(*newmat,m,n,PETSC_DECIDE,PETSC_DECIDE);
813:     /* communicate irow to their owners in the layout */
814:     ISGetIndices(irow,&idxs);
815:     PetscLayoutMapLocal_Private(mat->rmap,m,idxs,&ll,&lidxs,&lgidxs);
816:     ISRestoreIndices(irow,&idxs);
817:     MatISSetUpSF(mat);
818:     PetscMemzero(matis->sf_rootdata,matis->sf->nroots*sizeof(PetscInt));
819:     for (i=0;i<ll;i++) matis->sf_rootdata[lidxs[i]] = lgidxs[i]+1;
820:     PetscFree(lidxs);
821:     PetscFree(lgidxs);
822:     PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
823:     PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
824:     for (i=0,newloc=0;i<matis->sf->nleaves;i++) if (matis->sf_leafdata[i]) newloc++;
825:     PetscMalloc1(newloc,&newgidxs);
826:     PetscMalloc1(newloc,&lidxs);
827:     for (i=0,newloc=0;i<matis->sf->nleaves;i++)
828:       if (matis->sf_leafdata[i]) {
829:         lidxs[newloc] = i;
830:         newgidxs[newloc++] = matis->sf_leafdata[i]-1;
831:       }
832:     ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is);
833:     ISLocalToGlobalMappingCreateIS(is,&rl2g);
834:     ISDestroy(&is);
835:     /* local is to extract local submatrix */
836:     newmatis = (Mat_IS*)(*newmat)->data;
837:     ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_ris);
838:     if (mat->congruentlayouts == -1) { /* first time we compare rows and cols layouts */
839:       PetscBool cong;
840:       PetscLayoutCompare(mat->rmap,mat->cmap,&cong);
841:       if (cong) mat->congruentlayouts = 1;
842:       else      mat->congruentlayouts = 0;
843:     }
844:     if (mat->congruentlayouts && irow == icol) {
845:       MatSetLocalToGlobalMapping(*newmat,rl2g,rl2g);
846:       PetscObjectReference((PetscObject)newmatis->getsub_ris);
847:       newmatis->getsub_cis = newmatis->getsub_ris;
848:     } else {
849:       ISLocalToGlobalMapping cl2g;

851:       /* communicate icol to their owners in the layout */
852:       ISGetIndices(icol,&idxs);
853:       PetscLayoutMapLocal_Private(mat->cmap,n,idxs,&ll,&lidxs,&lgidxs);
854:       ISRestoreIndices(icol,&idxs);
855:       PetscMemzero(matis->csf_rootdata,matis->csf->nroots*sizeof(PetscInt));
856:       for (i=0;i<ll;i++) matis->csf_rootdata[lidxs[i]] = lgidxs[i]+1;
857:       PetscFree(lidxs);
858:       PetscFree(lgidxs);
859:       PetscSFBcastBegin(matis->csf,MPIU_INT,matis->csf_rootdata,matis->csf_leafdata);
860:       PetscSFBcastEnd(matis->csf,MPIU_INT,matis->csf_rootdata,matis->csf_leafdata);
861:       for (i=0,newloc=0;i<matis->csf->nleaves;i++) if (matis->csf_leafdata[i]) newloc++;
862:       PetscMalloc1(newloc,&newgidxs);
863:       PetscMalloc1(newloc,&lidxs);
864:       for (i=0,newloc=0;i<matis->csf->nleaves;i++)
865:         if (matis->csf_leafdata[i]) {
866:           lidxs[newloc] = i;
867:           newgidxs[newloc++] = matis->csf_leafdata[i]-1;
868:         }
869:       ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is);
870:       ISLocalToGlobalMappingCreateIS(is,&cl2g);
871:       ISDestroy(&is);
872:       /* local is to extract local submatrix */
873:       ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_cis);
874:       MatSetLocalToGlobalMapping(*newmat,rl2g,cl2g);
875:       ISLocalToGlobalMappingDestroy(&cl2g);
876:     }
877:     ISLocalToGlobalMappingDestroy(&rl2g);
878:   } else {
879:     MatISGetLocalMat(*newmat,&newlocmat);
880:   }
881:   MatISGetLocalMat(mat,&locmat);
882:   newmatis = (Mat_IS*)(*newmat)->data;
883:   MatCreateSubMatrix(locmat,newmatis->getsub_ris,newmatis->getsub_cis,scall,&newlocmat);
884:   if (scall == MAT_INITIAL_MATRIX) {
885:     MatISSetLocalMat(*newmat,newlocmat);
886:     MatDestroy(&newlocmat);
887:   }
888:   MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
889:   MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
890:   return(0);
891: }

893: static PetscErrorCode MatCopy_IS(Mat A,Mat B,MatStructure str)
894: {
895:   Mat_IS         *a = (Mat_IS*)A->data,*b;
896:   PetscBool      ismatis;

900:   PetscObjectTypeCompare((PetscObject)B,MATIS,&ismatis);
901:   if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Need to be implemented");
902:   b = (Mat_IS*)B->data;
903:   MatCopy(a->A,b->A,str);
904:   PetscObjectStateIncrease((PetscObject)B);
905:   return(0);
906: }

908: static PetscErrorCode MatMissingDiagonal_IS(Mat A,PetscBool  *missing,PetscInt *d)
909: {
910:   Vec               v;
911:   const PetscScalar *array;
912:   PetscInt          i,n;
913:   PetscErrorCode    ierr;

916:   *missing = PETSC_FALSE;
917:   MatCreateVecs(A,NULL,&v);
918:   MatGetDiagonal(A,v);
919:   VecGetLocalSize(v,&n);
920:   VecGetArrayRead(v,&array);
921:   for (i=0;i<n;i++) if (array[i] == 0.) break;
922:   VecRestoreArrayRead(v,&array);
923:   VecDestroy(&v);
924:   if (i != n) *missing = PETSC_TRUE;
925:   if (d) {
926:     *d = -1;
927:     if (*missing) {
928:       PetscInt rstart;
929:       MatGetOwnershipRange(A,&rstart,NULL);
930:       *d = i+rstart;
931:     }
932:   }
933:   return(0);
934: }

936: static PetscErrorCode MatISSetUpSF_IS(Mat B)
937: {
938:   Mat_IS         *matis = (Mat_IS*)(B->data);
939:   const PetscInt *gidxs;
940:   PetscInt       nleaves;

944:   if (matis->sf) return(0);
945:   PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->sf);
946:   ISLocalToGlobalMappingGetIndices(B->rmap->mapping,&gidxs);
947:   ISLocalToGlobalMappingGetSize(B->rmap->mapping,&nleaves);
948:   PetscSFSetGraphLayout(matis->sf,B->rmap,nleaves,NULL,PETSC_OWN_POINTER,gidxs);
949:   ISLocalToGlobalMappingRestoreIndices(B->rmap->mapping,&gidxs);
950:   PetscMalloc2(matis->sf->nroots,&matis->sf_rootdata,matis->sf->nleaves,&matis->sf_leafdata);
951:   if (B->rmap->mapping != B->cmap->mapping) { /* setup SF for columns */
952:     ISLocalToGlobalMappingGetSize(B->cmap->mapping,&nleaves);
953:     PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->csf);
954:     ISLocalToGlobalMappingGetIndices(B->cmap->mapping,&gidxs);
955:     PetscSFSetGraphLayout(matis->csf,B->cmap,nleaves,NULL,PETSC_OWN_POINTER,gidxs);
956:     ISLocalToGlobalMappingRestoreIndices(B->cmap->mapping,&gidxs);
957:     PetscMalloc2(matis->csf->nroots,&matis->csf_rootdata,matis->csf->nleaves,&matis->csf_leafdata);
958:   } else {
959:     matis->csf = matis->sf;
960:     matis->csf_leafdata = matis->sf_leafdata;
961:     matis->csf_rootdata = matis->sf_rootdata;
962:   }
963:   return(0);
964: }

966: /*@
967:    MatISSetPreallocation - Preallocates memory for a MATIS parallel matrix.

969:    Collective on MPI_Comm

971:    Input Parameters:
972: +  B - the matrix
973: .  d_nz  - number of nonzeros per row in DIAGONAL portion of local submatrix
974:            (same value is used for all local rows)
975: .  d_nnz - array containing the number of nonzeros in the various rows of the
976:            DIAGONAL portion of the local submatrix (possibly different for each row)
977:            or NULL, if d_nz is used to specify the nonzero structure.
978:            The size of this array is equal to the number of local rows, i.e 'm'.
979:            For matrices that will be factored, you must leave room for (and set)
980:            the diagonal entry even if it is zero.
981: .  o_nz  - number of nonzeros per row in the OFF-DIAGONAL portion of local
982:            submatrix (same value is used for all local rows).
983: -  o_nnz - array containing the number of nonzeros in the various rows of the
984:            OFF-DIAGONAL portion of the local submatrix (possibly different for
985:            each row) or NULL, if o_nz is used to specify the nonzero
986:            structure. The size of this array is equal to the number
987:            of local rows, i.e 'm'.

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

991:    Level: intermediate

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

997: .keywords: matrix

999: .seealso: MatCreate(), MatCreateIS(), MatMPIAIJSetPreallocation(), MatISGetLocalMat(), MATIS
1000: @*/
1001: PetscErrorCode  MatISSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
1002: {

1008:   PetscTryMethod(B,"MatISSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz));
1009:   return(0);
1010: }

1012: static PetscErrorCode  MatISSetPreallocation_IS(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
1013: {
1014:   Mat_IS         *matis = (Mat_IS*)(B->data);
1015:   PetscInt       bs,i,nlocalcols;

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

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

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

1028:   PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
1029:   MatGetSize(matis->A,NULL,&nlocalcols);
1030:   MatGetBlockSize(matis->A,&bs);
1031:   PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);

1033:   for (i=0;i<matis->sf->nleaves;i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i],nlocalcols);
1034:   MatSeqAIJSetPreallocation(matis->A,0,matis->sf_leafdata);
1035: #if defined(PETSC_HAVE_HYPRE)
1036:   MatHYPRESetPreallocation(matis->A,0,matis->sf_leafdata,0,NULL);
1037: #endif

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

1042:   for (i=0;i<matis->sf->nleaves/bs;i++) matis->sf_leafdata[i] = matis->sf_leafdata[i]-i;
1043:   MatSeqSBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata);

1045:   /* for other matrix types */
1046:   MatSetUp(matis->A);
1047:   return(0);
1048: }

1050: PETSC_EXTERN PetscErrorCode MatISSetMPIXAIJPreallocation_Private(Mat A, Mat B, PetscBool maxreduce)
1051: {
1052:   Mat_IS          *matis = (Mat_IS*)(A->data);
1053:   PetscInt        *my_dnz,*my_onz,*dnz,*onz,*mat_ranges,*row_ownership;
1054:   const PetscInt  *global_indices_r,*global_indices_c;
1055:   PetscInt        i,j,bs,rows,cols;
1056:   PetscInt        lrows,lcols;
1057:   PetscInt        local_rows,local_cols;
1058:   PetscMPIInt     nsubdomains;
1059:   PetscBool       isdense,issbaij;
1060:   PetscErrorCode  ierr;

1063:   MPI_Comm_size(PetscObjectComm((PetscObject)A),&nsubdomains);
1064:   MatGetSize(A,&rows,&cols);
1065:   MatGetBlockSize(A,&bs);
1066:   MatGetSize(matis->A,&local_rows,&local_cols);
1067:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isdense);
1068:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);
1069:   ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&global_indices_r);
1070:   if (A->rmap->mapping != A->cmap->mapping) {
1071:     ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&global_indices_c);
1072:   } else {
1073:     global_indices_c = global_indices_r;
1074:   }

1076:   if (issbaij) {
1077:     MatGetRowUpperTriangular(matis->A);
1078:   }
1079:   /*
1080:      An SF reduce is needed to sum up properly on shared rows.
1081:      Note that generally preallocation is not exact, since it overestimates nonzeros
1082:   */
1083:   MatISSetUpSF(A);
1084:   MatGetLocalSize(A,&lrows,&lcols);
1085:   MatPreallocateInitialize(PetscObjectComm((PetscObject)A),lrows,lcols,dnz,onz);
1086:   /* All processes need to compute entire row ownership */
1087:   PetscMalloc1(rows,&row_ownership);
1088:   MatGetOwnershipRanges(A,(const PetscInt**)&mat_ranges);
1089:   for (i=0;i<nsubdomains;i++) {
1090:     for (j=mat_ranges[i];j<mat_ranges[i+1];j++) {
1091:       row_ownership[j] = i;
1092:     }
1093:   }
1094:   MatGetOwnershipRangesColumn(A,(const PetscInt**)&mat_ranges);

1096:   /*
1097:      my_dnz and my_onz contains exact contribution to preallocation from each local mat
1098:      then, they will be summed up properly. This way, preallocation is always sufficient
1099:   */
1100:   PetscCalloc2(local_rows,&my_dnz,local_rows,&my_onz);
1101:   /* preallocation as a MATAIJ */
1102:   if (isdense) { /* special case for dense local matrices */
1103:     for (i=0;i<local_rows;i++) {
1104:       PetscInt owner = row_ownership[global_indices_r[i]];
1105:       for (j=0;j<local_cols;j++) {
1106:         PetscInt index_col = global_indices_c[j];
1107:         if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */
1108:           my_dnz[i] += 1;
1109:         } else { /* offdiag block */
1110:           my_onz[i] += 1;
1111:         }
1112:       }
1113:     }
1114:   } else if (matis->A->ops->getrowij) {
1115:     const PetscInt *ii,*jj,*jptr;
1116:     PetscBool      done;
1117:     MatGetRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done);
1118:     if (!done) SETERRQ(PetscObjectComm((PetscObject)(matis->A)),PETSC_ERR_PLIB,"Error in MatGetRowIJ");
1119:     jptr = jj;
1120:     for (i=0;i<local_rows;i++) {
1121:       PetscInt index_row = global_indices_r[i];
1122:       for (j=0;j<ii[i+1]-ii[i];j++,jptr++) {
1123:         PetscInt owner = row_ownership[index_row];
1124:         PetscInt index_col = global_indices_c[*jptr];
1125:         if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */
1126:           my_dnz[i] += 1;
1127:         } else { /* offdiag block */
1128:           my_onz[i] += 1;
1129:         }
1130:         /* same as before, interchanging rows and cols */
1131:         if (issbaij && index_col != index_row) {
1132:           owner = row_ownership[index_col];
1133:           if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) {
1134:             my_dnz[*jptr] += 1;
1135:           } else {
1136:             my_onz[*jptr] += 1;
1137:           }
1138:         }
1139:       }
1140:     }
1141:     MatRestoreRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done);
1142:     if (!done) SETERRQ(PetscObjectComm((PetscObject)(matis->A)),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ");
1143:   } else { /* loop over rows and use MatGetRow */
1144:     for (i=0;i<local_rows;i++) {
1145:       const PetscInt *cols;
1146:       PetscInt       ncols,index_row = global_indices_r[i];
1147:       MatGetRow(matis->A,i,&ncols,&cols,NULL);
1148:       for (j=0;j<ncols;j++) {
1149:         PetscInt owner = row_ownership[index_row];
1150:         PetscInt index_col = global_indices_c[cols[j]];
1151:         if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */
1152:           my_dnz[i] += 1;
1153:         } else { /* offdiag block */
1154:           my_onz[i] += 1;
1155:         }
1156:         /* same as before, interchanging rows and cols */
1157:         if (issbaij && index_col != index_row) {
1158:           owner = row_ownership[index_col];
1159:           if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) {
1160:             my_dnz[cols[j]] += 1;
1161:           } else {
1162:             my_onz[cols[j]] += 1;
1163:           }
1164:         }
1165:       }
1166:       MatRestoreRow(matis->A,i,&ncols,&cols,NULL);
1167:     }
1168:   }
1169:   if (global_indices_c != global_indices_r) {
1170:     ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&global_indices_c);
1171:   }
1172:   ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&global_indices_r);
1173:   PetscFree(row_ownership);

1175:   /* Reduce my_dnz and my_onz */
1176:   if (maxreduce) {
1177:     PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);
1178:     PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);
1179:     PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);
1180:     PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);
1181:   } else {
1182:     PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);
1183:     PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);
1184:     PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);
1185:     PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);
1186:   }
1187:   PetscFree2(my_dnz,my_onz);

1189:   /* Resize preallocation if overestimated */
1190:   for (i=0;i<lrows;i++) {
1191:     dnz[i] = PetscMin(dnz[i],lcols);
1192:     onz[i] = PetscMin(onz[i],cols-lcols);
1193:   }

1195:   /* Set preallocation */
1196:   MatMPIAIJSetPreallocation(B,0,dnz,0,onz);
1197:   for (i=0;i<lrows/bs;i++) {
1198:     dnz[i] = dnz[i*bs]/bs;
1199:     onz[i] = onz[i*bs]/bs;
1200:   }
1201:   MatMPIBAIJSetPreallocation(B,bs,0,dnz,0,onz);
1202:   MatMPISBAIJSetPreallocation(B,bs,0,dnz,0,onz);
1203:   MatPreallocateFinalize(dnz,onz);
1204:   if (issbaij) {
1205:     MatRestoreRowUpperTriangular(matis->A);
1206:   }
1207:   return(0);
1208: }

1210: static PetscErrorCode MatISGetMPIXAIJ_IS(Mat mat, MatReuse reuse, Mat *M)
1211: {
1212:   Mat_IS         *matis = (Mat_IS*)(mat->data);
1213:   Mat            local_mat;
1214:   /* info on mat */
1215:   PetscInt       bs,rows,cols,lrows,lcols;
1216:   PetscInt       local_rows,local_cols;
1217:   PetscBool      isseqdense,isseqsbaij,isseqaij,isseqbaij;
1218: #if defined (PETSC_USE_DEBUG)
1219:   PetscBool      lb[4],bb[4];
1220: #endif
1221:   PetscMPIInt    nsubdomains;
1222:   /* values insertion */
1223:   PetscScalar    *array;
1224:   /* work */

1228:   /* get info from mat */
1229:   MPI_Comm_size(PetscObjectComm((PetscObject)mat),&nsubdomains);
1230:   if (nsubdomains == 1) {
1231:     Mat            B;
1232:     IS             rows,cols;
1233:     IS             irows,icols;
1234:     const PetscInt *ridxs,*cidxs;

1236:     ISLocalToGlobalMappingGetIndices(mat->rmap->mapping,&ridxs);
1237:     ISLocalToGlobalMappingGetIndices(mat->cmap->mapping,&cidxs);
1238:     ISCreateGeneral(PETSC_COMM_SELF,mat->rmap->n,ridxs,PETSC_USE_POINTER,&rows);
1239:     ISCreateGeneral(PETSC_COMM_SELF,mat->cmap->n,cidxs,PETSC_USE_POINTER,&cols);
1240:     ISLocalToGlobalMappingRestoreIndices(mat->rmap->mapping,&ridxs);
1241:     ISLocalToGlobalMappingRestoreIndices(mat->cmap->mapping,&cidxs);
1242:     ISSetPermutation(rows);
1243:     ISSetPermutation(cols);
1244:     ISInvertPermutation(rows,mat->rmap->n,&irows);
1245:     ISInvertPermutation(cols,mat->cmap->n,&icols);
1246:     ISDestroy(&cols);
1247:     ISDestroy(&rows);
1248:     MatConvert(matis->A,MATSEQAIJ,MAT_INITIAL_MATRIX,&B);
1249:     MatCreateSubMatrix(B,irows,icols,reuse,M);
1250:     MatDestroy(&B);
1251:     ISDestroy(&icols);
1252:     ISDestroy(&irows);
1253:     return(0);
1254:   }
1255:   MatGetSize(mat,&rows,&cols);
1256:   MatGetBlockSize(mat,&bs);
1257:   MatGetLocalSize(mat,&lrows,&lcols);
1258:   MatGetSize(matis->A,&local_rows,&local_cols);
1259:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isseqdense);
1260:   PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQAIJ,&isseqaij);
1261:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQBAIJ,&isseqbaij);
1262:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&isseqsbaij);
1263:   if (!isseqdense && !isseqaij && !isseqbaij && !isseqsbaij) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)(matis->A))->type_name);
1264: #if defined (PETSC_USE_DEBUG)
1265:   lb[0] = isseqdense;
1266:   lb[1] = isseqaij;
1267:   lb[2] = isseqbaij;
1268:   lb[3] = isseqsbaij;
1269:   MPIU_Allreduce(lb,bb,4,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)mat));
1270:   if (!bb[0] && !bb[1] && !bb[2] && !bb[3]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Local matrices must have the same type");
1271: #endif

1273:   if (reuse == MAT_INITIAL_MATRIX) {
1274:     MatCreate(PetscObjectComm((PetscObject)mat),M);
1275:     MatSetSizes(*M,lrows,lcols,rows,cols);
1276:     MatSetBlockSize(*M,bs);
1277:     if (!isseqsbaij) {
1278:       MatSetType(*M,MATAIJ);
1279:     } else {
1280:       MatSetType(*M,MATSBAIJ);
1281:     }
1282:     MatISSetMPIXAIJPreallocation_Private(mat,*M,PETSC_FALSE);
1283:   } else {
1284:     PetscInt mbs,mrows,mcols,mlrows,mlcols;
1285:     /* some checks */
1286:     MatGetBlockSize(*M,&mbs);
1287:     MatGetSize(*M,&mrows,&mcols);
1288:     MatGetLocalSize(*M,&mlrows,&mlcols);
1289:     if (mrows != rows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of rows (%d != %d)",rows,mrows);
1290:     if (mcols != cols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of cols (%d != %d)",cols,mcols);
1291:     if (mlrows != lrows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local rows (%d != %d)",lrows,mlrows);
1292:     if (mlcols != lcols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local cols (%d != %d)",lcols,mlcols);
1293:     if (mbs != bs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong block size (%d != %d)",bs,mbs);
1294:     MatZeroEntries(*M);
1295:   }

1297:   if (isseqsbaij) {
1298:     MatConvert(matis->A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&local_mat);
1299:   } else {
1300:     PetscObjectReference((PetscObject)matis->A);
1301:     local_mat = matis->A;
1302:   }

1304:   /* Set values */
1305:   MatSetLocalToGlobalMapping(*M,mat->rmap->mapping,mat->cmap->mapping);
1306:   if (isseqdense) { /* special case for dense local matrices */
1307:     PetscInt i,*dummy;

1309:     PetscMalloc1(PetscMax(local_rows,local_cols),&dummy);
1310:     for (i=0;i<PetscMax(local_rows,local_cols);i++) dummy[i] = i;
1311:     MatSetOption(*M,MAT_ROW_ORIENTED,PETSC_FALSE);
1312:     MatDenseGetArray(local_mat,&array);
1313:     MatSetValuesLocal(*M,local_rows,dummy,local_cols,dummy,array,ADD_VALUES);
1314:     MatDenseRestoreArray(local_mat,&array);
1315:     PetscFree(dummy);
1316:   } else if (isseqaij) {
1317:     PetscInt  i,nvtxs,*xadj,*adjncy;
1318:     PetscBool done;

1320:     MatGetRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);
1321:     if (!done) SETERRQ(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatGetRowIJ");
1322:     MatSeqAIJGetArray(local_mat,&array);
1323:     for (i=0;i<nvtxs;i++) {
1324:       MatSetValuesLocal(*M,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],array+xadj[i],ADD_VALUES);
1325:     }
1326:     MatRestoreRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);
1327:     if (!done) SETERRQ(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ");
1328:     MatSeqAIJRestoreArray(local_mat,&array);
1329:   } else { /* very basic values insertion for all other matrix types */
1330:     PetscInt i;

1332:     for (i=0;i<local_rows;i++) {
1333:       PetscInt       j;
1334:       const PetscInt *local_indices_cols;

1336:       MatGetRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);
1337:       MatSetValuesLocal(*M,1,&i,j,local_indices_cols,array,ADD_VALUES);
1338:       MatRestoreRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);
1339:     }
1340:   }
1341:   MatAssemblyBegin(*M,MAT_FINAL_ASSEMBLY);
1342:   MatDestroy(&local_mat);
1343:   MatAssemblyEnd(*M,MAT_FINAL_ASSEMBLY);
1344:   if (isseqdense) {
1345:     MatSetOption(*M,MAT_ROW_ORIENTED,PETSC_TRUE);
1346:   }
1347:   return(0);
1348: }

1350: /*@
1351:     MatISGetMPIXAIJ - Converts MATIS matrix into a parallel AIJ format

1353:   Input Parameter:
1354: .  mat - the matrix (should be of type MATIS)
1355: .  reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX

1357:   Output Parameter:
1358: .  newmat - the matrix in AIJ format

1360:   Level: developer

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

1364: .seealso: MATIS
1365: @*/
1366: PetscErrorCode MatISGetMPIXAIJ(Mat mat, MatReuse reuse, Mat *newmat)
1367: {

1374:   if (reuse != MAT_INITIAL_MATRIX) {
1377:     if (mat == *newmat) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse the same matrix");
1378:   }
1379:   PetscUseMethod(mat,"MatISGetMPIXAIJ_C",(Mat,MatReuse,Mat*),(mat,reuse,newmat));
1380:   return(0);
1381: }

1383: PetscErrorCode MatDuplicate_IS(Mat mat,MatDuplicateOption op,Mat *newmat)
1384: {
1386:   Mat_IS         *matis = (Mat_IS*)(mat->data);
1387:   PetscInt       bs,m,n,M,N;
1388:   Mat            B,localmat;

1391:   MatGetBlockSize(mat,&bs);
1392:   MatGetSize(mat,&M,&N);
1393:   MatGetLocalSize(mat,&m,&n);
1394:   MatCreateIS(PetscObjectComm((PetscObject)mat),bs,m,n,M,N,mat->rmap->mapping,mat->cmap->mapping,&B);
1395:   MatDuplicate(matis->A,op,&localmat);
1396:   MatISSetLocalMat(B,localmat);
1397:   MatDestroy(&localmat);
1398:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1399:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1400:   *newmat = B;
1401:   return(0);
1402: }

1404: static PetscErrorCode MatIsHermitian_IS(Mat A,PetscReal tol,PetscBool  *flg)
1405: {
1407:   Mat_IS         *matis = (Mat_IS*)A->data;
1408:   PetscBool      local_sym;

1411:   MatIsHermitian(matis->A,tol,&local_sym);
1412:   MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
1413:   return(0);
1414: }

1416: static PetscErrorCode MatIsSymmetric_IS(Mat A,PetscReal tol,PetscBool  *flg)
1417: {
1419:   Mat_IS         *matis = (Mat_IS*)A->data;
1420:   PetscBool      local_sym;

1423:   MatIsSymmetric(matis->A,tol,&local_sym);
1424:   MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
1425:   return(0);
1426: }

1428: static PetscErrorCode MatIsStructurallySymmetric_IS(Mat A,PetscBool  *flg)
1429: {
1431:   Mat_IS         *matis = (Mat_IS*)A->data;
1432:   PetscBool      local_sym;

1435:   if (A->rmap->mapping != A->cmap->mapping) {
1436:     *flg = PETSC_FALSE;
1437:     return(0);
1438:   }
1439:   MatIsStructurallySymmetric(matis->A,&local_sym);
1440:   MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
1441:   return(0);
1442: }

1444: static PetscErrorCode MatDestroy_IS(Mat A)
1445: {
1447:   Mat_IS         *b = (Mat_IS*)A->data;

1450:   MatDestroy(&b->A);
1451:   VecScatterDestroy(&b->cctx);
1452:   VecScatterDestroy(&b->rctx);
1453:   VecDestroy(&b->x);
1454:   VecDestroy(&b->y);
1455:   VecDestroy(&b->counter);
1456:   ISDestroy(&b->getsub_ris);
1457:   ISDestroy(&b->getsub_cis);
1458:   if (b->sf != b->csf) {
1459:     PetscSFDestroy(&b->csf);
1460:     PetscFree2(b->csf_rootdata,b->csf_leafdata);
1461:   }
1462:   PetscSFDestroy(&b->sf);
1463:   PetscFree2(b->sf_rootdata,b->sf_leafdata);
1464:   PetscFree(A->data);
1465:   PetscObjectChangeTypeName((PetscObject)A,0);
1466:   PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL);
1467:   PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",NULL);
1468:   PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",NULL);
1469:   PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",NULL);
1470:   PetscObjectComposeFunction((PetscObject)A,"MatISSetUpSF_C",NULL);
1471:   return(0);
1472: }

1474: static PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y)
1475: {
1477:   Mat_IS         *is  = (Mat_IS*)A->data;
1478:   PetscScalar    zero = 0.0;

1481:   /*  scatter the global vector x into the local work vector */
1482:   VecScatterBegin(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
1483:   VecScatterEnd(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);

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

1488:   /* scatter product back into global memory */
1489:   VecSet(y,zero);
1490:   VecScatterBegin(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
1491:   VecScatterEnd(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
1492:   return(0);
1493: }

1495: static PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
1496: {
1497:   Vec            temp_vec;

1501:   if (v3 != v2) {
1502:     MatMult(A,v1,v3);
1503:     VecAXPY(v3,1.0,v2);
1504:   } else {
1505:     VecDuplicate(v2,&temp_vec);
1506:     MatMult(A,v1,temp_vec);
1507:     VecAXPY(temp_vec,1.0,v2);
1508:     VecCopy(temp_vec,v3);
1509:     VecDestroy(&temp_vec);
1510:   }
1511:   return(0);
1512: }

1514: static PetscErrorCode MatMultTranspose_IS(Mat A,Vec y,Vec x)
1515: {
1516:   Mat_IS         *is = (Mat_IS*)A->data;

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

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

1527:   /* scatter product back into global vector */
1528:   VecSet(x,0);
1529:   VecScatterBegin(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);
1530:   VecScatterEnd(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);
1531:   return(0);
1532: }

1534: static PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
1535: {
1536:   Vec            temp_vec;

1540:   if (v3 != v2) {
1541:     MatMultTranspose(A,v1,v3);
1542:     VecAXPY(v3,1.0,v2);
1543:   } else {
1544:     VecDuplicate(v2,&temp_vec);
1545:     MatMultTranspose(A,v1,temp_vec);
1546:     VecAXPY(temp_vec,1.0,v2);
1547:     VecCopy(temp_vec,v3);
1548:     VecDestroy(&temp_vec);
1549:   }
1550:   return(0);
1551: }

1553: static PetscErrorCode MatView_IS(Mat A,PetscViewer viewer)
1554: {
1555:   Mat_IS         *a = (Mat_IS*)A->data;
1557:   PetscViewer    sviewer;
1558:   PetscBool      isascii,view = PETSC_TRUE;

1561:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1562:   if (isascii)  {
1563:     PetscViewerFormat format;

1565:     PetscViewerGetFormat(viewer,&format);
1566:     if (format == PETSC_VIEWER_ASCII_INFO) view = PETSC_FALSE;
1567:   }
1568:   if (!view) return(0);
1569:   PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
1570:   MatView(a->A,sviewer);
1571:   PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
1572:   PetscViewerFlush(viewer);
1573:   return(0);
1574: }

1576: static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping)
1577: {
1579:   PetscInt       nr,rbs,nc,cbs;
1580:   Mat_IS         *is = (Mat_IS*)A->data;
1581:   IS             from,to;
1582:   Vec            cglobal,rglobal;

1587:   /* Destroy any previous data */
1588:   VecDestroy(&is->x);
1589:   VecDestroy(&is->y);
1590:   VecDestroy(&is->counter);
1591:   VecScatterDestroy(&is->rctx);
1592:   VecScatterDestroy(&is->cctx);
1593:   MatDestroy(&is->A);
1594:   if (is->csf != is->sf) {
1595:     PetscSFDestroy(&is->csf);
1596:     PetscFree2(is->csf_rootdata,is->csf_leafdata);
1597:   }
1598:   PetscSFDestroy(&is->sf);
1599:   PetscFree2(is->sf_rootdata,is->sf_leafdata);

1601:   /* Setup Layout and set local to global maps */
1602:   PetscLayoutSetUp(A->rmap);
1603:   PetscLayoutSetUp(A->cmap);
1604:   PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping);
1605:   PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping);

1607:   /* Create the local matrix A */
1608:   ISLocalToGlobalMappingGetSize(rmapping,&nr);
1609:   ISLocalToGlobalMappingGetBlockSize(rmapping,&rbs);
1610:   ISLocalToGlobalMappingGetSize(cmapping,&nc);
1611:   ISLocalToGlobalMappingGetBlockSize(cmapping,&cbs);
1612:   MatCreate(PETSC_COMM_SELF,&is->A);
1613:   MatSetType(is->A,MATAIJ);
1614:   MatSetSizes(is->A,nr,nc,nr,nc);
1615:   MatSetBlockSizes(is->A,rbs,cbs);
1616:   MatSetOptionsPrefix(is->A,((PetscObject)A)->prefix);
1617:   MatAppendOptionsPrefix(is->A,"is_");
1618:   MatSetFromOptions(is->A);
1619:   PetscLayoutSetUp(is->A->rmap);
1620:   PetscLayoutSetUp(is->A->cmap);

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

1626:     /* setup the global to local scatters */
1627:     MatCreateVecs(A,&cglobal,&rglobal);
1628:     ISCreateStride(PETSC_COMM_SELF,nr,0,1,&to);
1629:     ISLocalToGlobalMappingApplyIS(rmapping,to,&from);
1630:     VecScatterCreate(rglobal,from,is->y,to,&is->rctx);
1631:     if (rmapping != cmapping) {
1632:       ISDestroy(&to);
1633:       ISDestroy(&from);
1634:       ISCreateStride(PETSC_COMM_SELF,nc,0,1,&to);
1635:       ISLocalToGlobalMappingApplyIS(cmapping,to,&from);
1636:       VecScatterCreate(cglobal,from,is->x,to,&is->cctx);
1637:     } else {
1638:       PetscObjectReference((PetscObject)is->rctx);
1639:       is->cctx = is->rctx;
1640:     }

1642:     /* interface counter vector (local) */
1643:     VecDuplicate(is->y,&is->counter);
1644:     VecSet(is->y,1.);
1645:     VecScatterBegin(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);
1646:     VecScatterEnd(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);
1647:     VecScatterBegin(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);
1648:     VecScatterEnd(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);

1650:     /* free workspace */
1651:     VecDestroy(&rglobal);
1652:     VecDestroy(&cglobal);
1653:     ISDestroy(&to);
1654:     ISDestroy(&from);
1655:   }
1656:   MatSetUp(A);
1657:   return(0);
1658: }

1660: static PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv)
1661: {
1662:   Mat_IS         *is = (Mat_IS*)mat->data;
1664: #if defined(PETSC_USE_DEBUG)
1665:   PetscInt       i,zm,zn;
1666: #endif
1667:   PetscInt       rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION];

1670: #if defined(PETSC_USE_DEBUG)
1671:   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);
1672:   /* count negative indices */
1673:   for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++;
1674:   for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++;
1675: #endif
1676:   ISGlobalToLocalMappingApply(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);
1677:   ISGlobalToLocalMappingApply(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);
1678: #if defined(PETSC_USE_DEBUG)
1679:   /* count negative indices (should be the same as before) */
1680:   for (i=0;i<m;i++) if (rows_l[i] < 0) zm--;
1681:   for (i=0;i<n;i++) if (cols_l[i] < 0) zn--;
1682:   /* disable check when usesetlocal is true */
1683:   if (!is->usesetlocal && zm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the row indices can not be mapped! Maybe you should not use MATIS");
1684:   if (!is->usesetlocal && zn) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the column indices can not be mapped! Maybe you should not use MATIS");
1685: #endif
1686:   if (is->usesetlocal) {
1687:     MatSetValuesLocal(is->A,m,rows_l,n,cols_l,values,addv);
1688:   } else {
1689:     MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);
1690:   }
1691:   return(0);
1692: }

1694: static PetscErrorCode MatSetValuesBlocked_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv)
1695: {
1696:   Mat_IS         *is = (Mat_IS*)mat->data;
1698: #if defined(PETSC_USE_DEBUG)
1699:   PetscInt       i,zm,zn;
1700: #endif
1701:   PetscInt       rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION];

1704: #if defined(PETSC_USE_DEBUG)
1705:   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);
1706:   /* count negative indices */
1707:   for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++;
1708:   for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++;
1709: #endif
1710:   ISGlobalToLocalMappingApplyBlock(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);
1711:   ISGlobalToLocalMappingApplyBlock(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);
1712: #if defined(PETSC_USE_DEBUG)
1713:   /* count negative indices (should be the same as before) */
1714:   for (i=0;i<m;i++) if (rows_l[i] < 0) zm--;
1715:   for (i=0;i<n;i++) if (cols_l[i] < 0) zn--;
1716:   if (zm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the row indices can not be mapped! Maybe you should not use MATIS");
1717:   if (zn) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the column indices can not be mapped! Maybe you should not use MATIS");
1718: #endif
1719:   MatSetValuesBlocked(is->A,m,rows_l,n,cols_l,values,addv);
1720:   return(0);
1721: }

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

1729:   if (is->usesetlocal) {
1730:     MatSetValuesLocal(is->A,m,rows,n,cols,values,addv);
1731:   } else {
1732:     MatSetValues(is->A,m,rows,n,cols,values,addv);
1733:   }
1734:   return(0);
1735: }

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

1743:   MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv);
1744:   return(0);
1745: }

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

1753:   if (!n) {
1754:     is->pure_neumann = PETSC_TRUE;
1755:   } else {
1756:     PetscInt i;
1757:     is->pure_neumann = PETSC_FALSE;

1759:     if (columns) {
1760:       MatZeroRowsColumns(is->A,n,rows,diag,0,0);
1761:     } else {
1762:       MatZeroRows(is->A,n,rows,diag,0,0);
1763:     }
1764:     if (diag != 0.) {
1765:       const PetscScalar *array;
1766:       VecGetArrayRead(is->counter,&array);
1767:       for (i=0; i<n; i++) {
1768:         MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);
1769:       }
1770:       VecRestoreArrayRead(is->counter,&array);
1771:     }
1772:     MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);
1773:     MatAssemblyEnd(is->A,MAT_FINAL_ASSEMBLY);
1774:   }
1775:   return(0);
1776: }

1778: static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b,PetscBool columns)
1779: {
1780:   Mat_IS         *matis = (Mat_IS*)A->data;
1781:   PetscInt       nr,nl,len,i;
1782:   PetscInt       *lrows;

1786: #if defined(PETSC_USE_DEBUG)
1787:   if (columns || diag != 0. || (x && b)) {
1788:     PetscBool cong;
1789:     PetscLayoutCompare(A->rmap,A->cmap,&cong);
1790:     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");
1791:     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");
1792:     if (!cong && x && b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"A->rmap and A->cmap need to be congruent");
1793:   }
1794: #endif
1795:   /* get locally owned rows */
1796:   PetscLayoutMapLocal_Private(A->rmap,n,rows,&len,&lrows,NULL);
1797:   /* fix right hand side if needed */
1798:   if (x && b) {
1799:     const PetscScalar *xx;
1800:     PetscScalar       *bb;

1802:     VecGetArrayRead(x, &xx);
1803:     VecGetArray(b, &bb);
1804:     for (i=0;i<len;++i) bb[lrows[i]] = diag*xx[lrows[i]];
1805:     VecRestoreArrayRead(x, &xx);
1806:     VecRestoreArray(b, &bb);
1807:   }
1808:   /* get rows associated to the local matrices */
1809:   MatISSetUpSF(A);
1810:   MatGetSize(matis->A,&nl,NULL);
1811:   PetscMemzero(matis->sf_leafdata,nl*sizeof(PetscInt));
1812:   PetscMemzero(matis->sf_rootdata,A->rmap->n*sizeof(PetscInt));
1813:   for (i=0;i<len;i++) matis->sf_rootdata[lrows[i]] = 1;
1814:   PetscFree(lrows);
1815:   PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
1816:   PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
1817:   PetscMalloc1(nl,&lrows);
1818:   for (i=0,nr=0;i<nl;i++) if (matis->sf_leafdata[i]) lrows[nr++] = i;
1819:   MatISZeroRowsColumnsLocal_Private(A,nr,lrows,diag,columns);
1820:   PetscFree(lrows);
1821:   return(0);
1822: }

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

1829:   MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_FALSE);
1830:   return(0);
1831: }

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

1838:   MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_TRUE);
1839:   return(0);
1840: }

1842: static PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type)
1843: {
1844:   Mat_IS         *is = (Mat_IS*)A->data;

1848:   MatAssemblyBegin(is->A,type);
1849:   return(0);
1850: }

1852: static PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type)
1853: {
1854:   Mat_IS         *is = (Mat_IS*)A->data;

1858:   MatAssemblyEnd(is->A,type);
1859:   /* fix for local empty rows/cols */
1860:   if (is->locempty && type == MAT_FINAL_ASSEMBLY) {
1861:     Mat                    newlA;
1862:     ISLocalToGlobalMapping l2g;
1863:     IS                     tis;
1864:     const PetscScalar      *v;
1865:     PetscInt               i,n,cf,*loce,*locf;
1866:     PetscBool              sym;

1868:     MatGetRowMaxAbs(is->A,is->y,NULL);
1869:     MatIsSymmetric(is->A,PETSC_SMALL,&sym);
1870:     if (!sym) SETERRQ(PetscObjectComm((PetscObject)is->A),PETSC_ERR_SUP,"Not yet implemented for unsymmetric case");
1871:     VecGetLocalSize(is->y,&n);
1872:     PetscMalloc1(n,&loce);
1873:     PetscMalloc1(n,&locf);
1874:     VecGetArrayRead(is->y,&v);
1875:     for (i=0,cf=0;i<n;i++) {
1876:       if (v[i] == 0.0) {
1877:         loce[i] = -1;
1878:       } else {
1879:         loce[i]    = cf;
1880:         locf[cf++] = i;
1881:       }
1882:     }
1883:     VecRestoreArrayRead(is->y,&v);
1884:     /* extract valid submatrix */
1885:     ISCreateGeneral(PetscObjectComm((PetscObject)is->A),cf,locf,PETSC_USE_POINTER,&tis);
1886:     MatCreateSubMatrix(is->A,tis,tis,MAT_INITIAL_MATRIX,&newlA);
1887:     ISDestroy(&tis);
1888:     /* attach local l2g map for successive calls of MatSetValues */
1889:     ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)is->A),1,n,loce,PETSC_OWN_POINTER,&l2g);
1890:     MatSetLocalToGlobalMapping(newlA,l2g,l2g);
1891:     ISLocalToGlobalMappingDestroy(&l2g);
1892:     /* flag MatSetValues */
1893:     is->usesetlocal = PETSC_TRUE;
1894:     /* attach new global l2g map */
1895:     ISLocalToGlobalMappingApply(A->rmap->mapping,cf,locf,locf);
1896:     ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,cf,locf,PETSC_OWN_POINTER,&l2g);
1897:     MatSetLocalToGlobalMapping(A,l2g,l2g);
1898:     MatISSetLocalMat(A,newlA);
1899:     MatDestroy(&newlA);
1900:     ISLocalToGlobalMappingDestroy(&l2g);
1901:   }
1902:   is->locempty = PETSC_FALSE;
1903:   return(0);
1904: }

1906: static PetscErrorCode MatISGetLocalMat_IS(Mat mat,Mat *local)
1907: {
1908:   Mat_IS *is = (Mat_IS*)mat->data;

1911:   *local = is->A;
1912:   return(0);
1913: }

1915: static PetscErrorCode MatISRestoreLocalMat_IS(Mat mat,Mat *local)
1916: {
1918:   *local = NULL;
1919:   return(0);
1920: }

1922: /*@
1923:     MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix.

1925:   Input Parameter:
1926: .  mat - the matrix

1928:   Output Parameter:
1929: .  local - the local matrix

1931:   Level: advanced

1933:   Notes:
1934:     This can be called if you have precomputed the nonzero structure of the
1935:   matrix and want to provide it to the inner matrix object to improve the performance
1936:   of the MatSetValues() operation.

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

1940: .seealso: MATIS
1941: @*/
1942: PetscErrorCode MatISGetLocalMat(Mat mat,Mat *local)
1943: {

1949:   PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat*),(mat,local));
1950:   return(0);
1951: }

1953: /*@
1954:     MatISRestoreLocalMat - Restores the local matrix obtained with MatISGetLocalMat()

1956:   Input Parameter:
1957: .  mat - the matrix

1959:   Output Parameter:
1960: .  local - the local matrix

1962:   Level: advanced

1964: .seealso: MATIS
1965: @*/
1966: PetscErrorCode MatISRestoreLocalMat(Mat mat,Mat *local)
1967: {

1973:   PetscUseMethod(mat,"MatISRestoreLocalMat_C",(Mat,Mat*),(mat,local));
1974:   return(0);
1975: }

1977: static PetscErrorCode MatISSetLocalMat_IS(Mat mat,Mat local)
1978: {
1979:   Mat_IS         *is = (Mat_IS*)mat->data;
1980:   PetscInt       nrows,ncols,orows,ocols;

1984:   if (is->A) {
1985:     MatGetSize(is->A,&orows,&ocols);
1986:     MatGetSize(local,&nrows,&ncols);
1987:     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);
1988:   }
1989:   PetscObjectReference((PetscObject)local);
1990:   MatDestroy(&is->A);
1991:   is->A = local;
1992:   return(0);
1993: }

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

1998:   Input Parameter:
1999: .  mat - the matrix
2000: .  local - the local matrix

2002:   Output Parameter:

2004:   Level: advanced

2006:   Notes:
2007:     This can be called if you have precomputed the local matrix and
2008:   want to provide it to the matrix object MATIS.

2010: .seealso: MATIS
2011: @*/
2012: PetscErrorCode MatISSetLocalMat(Mat mat,Mat local)
2013: {

2019:   PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local));
2020:   return(0);
2021: }

2023: static PetscErrorCode MatZeroEntries_IS(Mat A)
2024: {
2025:   Mat_IS         *a = (Mat_IS*)A->data;

2029:   MatZeroEntries(a->A);
2030:   return(0);
2031: }

2033: static PetscErrorCode MatScale_IS(Mat A,PetscScalar a)
2034: {
2035:   Mat_IS         *is = (Mat_IS*)A->data;

2039:   MatScale(is->A,a);
2040:   return(0);
2041: }

2043: static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v)
2044: {
2045:   Mat_IS         *is = (Mat_IS*)A->data;

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

2052:   /* scatter diagonal back into global vector */
2053:   VecSet(v,0);
2054:   VecScatterBegin(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);
2055:   VecScatterEnd(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);
2056:   return(0);
2057: }

2059: static PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg)
2060: {
2061:   Mat_IS         *a = (Mat_IS*)A->data;

2065:   MatSetOption(a->A,op,flg);
2066:   return(0);
2067: }

2069: static PetscErrorCode MatAXPY_IS(Mat Y,PetscScalar a,Mat X,MatStructure str)
2070: {
2071:   Mat_IS         *y = (Mat_IS*)Y->data;
2072:   Mat_IS         *x;
2073: #if defined(PETSC_USE_DEBUG)
2074:   PetscBool      ismatis;
2075: #endif

2079: #if defined(PETSC_USE_DEBUG)
2080:   PetscObjectTypeCompare((PetscObject)X,MATIS,&ismatis);
2081:   if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot call MatAXPY(Y,a,X,str) with X not of type MATIS");
2082: #endif
2083:   x = (Mat_IS*)X->data;
2084:   MatAXPY(y->A,a,x->A,str);
2085:   return(0);
2086: }

2088: static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A,IS row,IS col,Mat *submat)
2089: {
2090:   Mat                    lA;
2091:   Mat_IS                 *matis;
2092:   ISLocalToGlobalMapping rl2g,cl2g;
2093:   IS                     is;
2094:   const PetscInt         *rg,*rl;
2095:   PetscInt               nrg;
2096:   PetscInt               N,M,nrl,i,*idxs;
2097:   PetscErrorCode         ierr;

2100:   ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&rg);
2101:   ISGetLocalSize(row,&nrl);
2102:   ISGetIndices(row,&rl);
2103:   ISLocalToGlobalMappingGetSize(A->rmap->mapping,&nrg);
2104: #if defined(PETSC_USE_DEBUG)
2105:   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);
2106: #endif
2107:   PetscMalloc1(nrg,&idxs);
2108:   /* map from [0,nrl) to row */
2109:   for (i=0;i<nrl;i++) idxs[i] = rl[i];
2110: #if defined(PETSC_USE_DEBUG)
2111:   for (i=nrl;i<nrg;i++) idxs[i] = nrg;
2112: #else
2113:   for (i=nrl;i<nrg;i++) idxs[i] = -1;
2114: #endif
2115:   ISRestoreIndices(row,&rl);
2116:   ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&rg);
2117:   ISCreateGeneral(PetscObjectComm((PetscObject)A),nrg,idxs,PETSC_OWN_POINTER,&is);
2118:   ISLocalToGlobalMappingCreateIS(is,&rl2g);
2119:   ISDestroy(&is);
2120:   /* compute new l2g map for columns */
2121:   if (col != row || A->rmap->mapping != A->cmap->mapping) {
2122:     const PetscInt *cg,*cl;
2123:     PetscInt       ncg;
2124:     PetscInt       ncl;

2126:     ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&cg);
2127:     ISGetLocalSize(col,&ncl);
2128:     ISGetIndices(col,&cl);
2129:     ISLocalToGlobalMappingGetSize(A->cmap->mapping,&ncg);
2130: #if defined(PETSC_USE_DEBUG)
2131:     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);
2132: #endif
2133:     PetscMalloc1(ncg,&idxs);
2134:     /* map from [0,ncl) to col */
2135:     for (i=0;i<ncl;i++) idxs[i] = cl[i];
2136: #if defined(PETSC_USE_DEBUG)
2137:     for (i=ncl;i<ncg;i++) idxs[i] = ncg;
2138: #else
2139:     for (i=ncl;i<ncg;i++) idxs[i] = -1;
2140: #endif
2141:     ISRestoreIndices(col,&cl);
2142:     ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&cg);
2143:     ISCreateGeneral(PetscObjectComm((PetscObject)A),ncg,idxs,PETSC_OWN_POINTER,&is);
2144:     ISLocalToGlobalMappingCreateIS(is,&cl2g);
2145:     ISDestroy(&is);
2146:   } else {
2147:     PetscObjectReference((PetscObject)rl2g);
2148:     cl2g = rl2g;
2149:   }
2150:   /* create the MATIS submatrix */
2151:   MatGetSize(A,&M,&N);
2152:   MatCreate(PetscObjectComm((PetscObject)A),submat);
2153:   MatSetSizes(*submat,PETSC_DECIDE,PETSC_DECIDE,M,N);
2154:   MatSetType(*submat,MATIS);
2155:   matis = (Mat_IS*)((*submat)->data);
2156:   matis->islocalref = PETSC_TRUE;
2157:   MatSetLocalToGlobalMapping(*submat,rl2g,cl2g);
2158:   MatISGetLocalMat(A,&lA);
2159:   MatISSetLocalMat(*submat,lA);
2160:   MatSetUp(*submat);
2161:   MatAssemblyBegin(*submat,MAT_FINAL_ASSEMBLY);
2162:   MatAssemblyEnd(*submat,MAT_FINAL_ASSEMBLY);
2163:   ISLocalToGlobalMappingDestroy(&rl2g);
2164:   ISLocalToGlobalMappingDestroy(&cl2g);
2165:   /* remove unsupported ops */
2166:   PetscMemzero((*submat)->ops,sizeof(struct _MatOps));
2167:   (*submat)->ops->destroy               = MatDestroy_IS;
2168:   (*submat)->ops->setvalueslocal        = MatSetValuesLocal_SubMat_IS;
2169:   (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS;
2170:   (*submat)->ops->assemblybegin         = MatAssemblyBegin_IS;
2171:   (*submat)->ops->assemblyend           = MatAssemblyEnd_IS;
2172:   return(0);
2173: }

2175: static PetscErrorCode MatSetFromOptions_IS(PetscOptionItems *PetscOptionsObject, Mat A)
2176: {
2177:   Mat_IS         *a = (Mat_IS*)A->data;

2181:   PetscOptionsHead(PetscOptionsObject,"MATIS options");
2182:   PetscObjectOptionsBegin((PetscObject)A);
2183:   PetscOptionsBool("-matis_fixempty","Fix local matrices in case of empty local rows/columns",NULL,a->locempty,&a->locempty,NULL);
2184:   PetscOptionsEnd();
2185:   return(0);
2186: }

2188: /*@
2189:     MatCreateIS - Creates a "process" unassembled matrix, assembled on each
2190:        process but not across processes.

2192:    Input Parameters:
2193: +     comm    - MPI communicator that will share the matrix
2194: .     bs      - block size of the matrix
2195: .     m,n,M,N - local and/or global sizes of the left and right vector used in matrix vector products
2196: .     rmap    - local to global map for rows
2197: -     cmap    - local to global map for cols

2199:    Output Parameter:
2200: .    A - the resulting matrix

2202:    Level: advanced

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

2209: .seealso: MATIS, MatSetLocalToGlobalMapping()
2210: @*/
2211: PetscErrorCode  MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping rmap,ISLocalToGlobalMapping cmap,Mat *A)
2212: {

2216:   if (!rmap && !cmap) SETERRQ(comm,PETSC_ERR_USER,"You need to provide at least one of the mappings");
2217:   MatCreate(comm,A);
2218:   MatSetSizes(*A,m,n,M,N);
2219:   if (bs > 0) {
2220:     MatSetBlockSize(*A,bs);
2221:   }
2222:   MatSetType(*A,MATIS);
2223:   if (rmap && cmap) {
2224:     MatSetLocalToGlobalMapping(*A,rmap,cmap);
2225:   } else if (!rmap) {
2226:     MatSetLocalToGlobalMapping(*A,cmap,cmap);
2227:   } else {
2228:     MatSetLocalToGlobalMapping(*A,rmap,rmap);
2229:   }
2230:   return(0);
2231: }

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

2239:    Operations Provided:
2240: +  MatMult()
2241: .  MatMultAdd()
2242: .  MatMultTranspose()
2243: .  MatMultTransposeAdd()
2244: .  MatZeroEntries()
2245: .  MatSetOption()
2246: .  MatZeroRows()
2247: .  MatSetValues()
2248: .  MatSetValuesBlocked()
2249: .  MatSetValuesLocal()
2250: .  MatSetValuesBlockedLocal()
2251: .  MatScale()
2252: .  MatGetDiagonal()
2253: .  MatMissingDiagonal()
2254: .  MatDuplicate()
2255: .  MatCopy()
2256: .  MatAXPY()
2257: .  MatCreateSubMatrix()
2258: .  MatGetLocalSubMatrix()
2259: .  MatTranspose()
2260: -  MatSetLocalToGlobalMapping()

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

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

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

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

2272:   Level: advanced

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

2276: M*/

2278: PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A)
2279: {
2281:   Mat_IS         *b;

2284:   PetscNewLog(A,&b);
2285:   A->data = (void*)b;

2287:   /* matrix ops */
2288:   PetscMemzero(A->ops,sizeof(struct _MatOps));
2289:   A->ops->mult                    = MatMult_IS;
2290:   A->ops->multadd                 = MatMultAdd_IS;
2291:   A->ops->multtranspose           = MatMultTranspose_IS;
2292:   A->ops->multtransposeadd        = MatMultTransposeAdd_IS;
2293:   A->ops->destroy                 = MatDestroy_IS;
2294:   A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS;
2295:   A->ops->setvalues               = MatSetValues_IS;
2296:   A->ops->setvaluesblocked        = MatSetValuesBlocked_IS;
2297:   A->ops->setvalueslocal          = MatSetValuesLocal_IS;
2298:   A->ops->setvaluesblockedlocal   = MatSetValuesBlockedLocal_IS;
2299:   A->ops->zerorows                = MatZeroRows_IS;
2300:   A->ops->zerorowscolumns         = MatZeroRowsColumns_IS;
2301:   A->ops->assemblybegin           = MatAssemblyBegin_IS;
2302:   A->ops->assemblyend             = MatAssemblyEnd_IS;
2303:   A->ops->view                    = MatView_IS;
2304:   A->ops->zeroentries             = MatZeroEntries_IS;
2305:   A->ops->scale                   = MatScale_IS;
2306:   A->ops->getdiagonal             = MatGetDiagonal_IS;
2307:   A->ops->setoption               = MatSetOption_IS;
2308:   A->ops->ishermitian             = MatIsHermitian_IS;
2309:   A->ops->issymmetric             = MatIsSymmetric_IS;
2310:   A->ops->isstructurallysymmetric = MatIsStructurallySymmetric_IS;
2311:   A->ops->duplicate               = MatDuplicate_IS;
2312:   A->ops->missingdiagonal         = MatMissingDiagonal_IS;
2313:   A->ops->copy                    = MatCopy_IS;
2314:   A->ops->getlocalsubmatrix       = MatGetLocalSubMatrix_IS;
2315:   A->ops->createsubmatrix         = MatCreateSubMatrix_IS;
2316:   A->ops->axpy                    = MatAXPY_IS;
2317:   A->ops->diagonalset             = MatDiagonalSet_IS;
2318:   A->ops->shift                   = MatShift_IS;
2319:   A->ops->transpose               = MatTranspose_IS;
2320:   A->ops->getinfo                 = MatGetInfo_IS;
2321:   A->ops->diagonalscale           = MatDiagonalScale_IS;
2322:   A->ops->setfromoptions          = MatSetFromOptions_IS;

2324:   /* special MATIS functions */
2325:   PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS);
2326:   PetscObjectComposeFunction((PetscObject)A,"MatISRestoreLocalMat_C",MatISRestoreLocalMat_IS);
2327:   PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS);
2328:   PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",MatISGetMPIXAIJ_IS);
2329:   PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",MatISSetPreallocation_IS);
2330:   PetscObjectComposeFunction((PetscObject)A,"MatISSetUpSF_C",MatISSetUpSF_IS);
2331:   PetscObjectChangeTypeName((PetscObject)A,MATIS);
2332:   return(0);
2333: }