Actual source code: matis.c

petsc-3.6.1 2015-08-06
Report Typos and Errors
  2: /*
  3:     Creates a matrix class for using the Neumann-Neumann type preconditioners.
  4:    This stores the matrices in globally unassembled form. Each processor
  5:    assembles only its local Neumann problem and the parallel matrix vector
  6:    product is handled "implicitly".

  8:      We provide:
  9:          MatMult()

 11:     Currently this allows for only one subdomain per processor.

 13: */

 15: #include <../src/mat/impls/is/matis.h>      /*I "petscmat.h" I*/

 17: static PetscErrorCode MatISComputeSF_Private(Mat);

 21: static PetscErrorCode MatISComputeSF_Private(Mat B)
 22: {
 23:   Mat_IS         *matis = (Mat_IS*)(B->data);
 24:   const PetscInt *gidxs;

 28:   MatGetSize(matis->A,&matis->sf_nleaves,NULL);
 29:   MatGetLocalSize(B,&matis->sf_nroots,NULL);
 30:   PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->sf);
 31:   ISLocalToGlobalMappingGetIndices(matis->mapping,&gidxs);
 32:   /* PETSC_OWN_POINTER refers to ilocal which is NULL */
 33:   PetscSFSetGraphLayout(matis->sf,B->rmap,matis->sf_nleaves,NULL,PETSC_OWN_POINTER,gidxs);
 34:   ISLocalToGlobalMappingRestoreIndices(matis->mapping,&gidxs);
 35:   PetscMalloc2(matis->sf_nroots,&matis->sf_rootdata,matis->sf_nleaves,&matis->sf_leafdata);
 36:   return(0);
 37: }

 41: /*@
 42:    MatISSetPreallocation - Preallocates memory for a MATIS parallel matrix.

 44:    Collective on MPI_Comm

 46:    Input Parameters:
 47: +  B - the matrix
 48: .  d_nz  - number of nonzeros per row in DIAGONAL portion of local submatrix
 49:            (same value is used for all local rows)
 50: .  d_nnz - array containing the number of nonzeros in the various rows of the
 51:            DIAGONAL portion of the local submatrix (possibly different for each row)
 52:            or NULL, if d_nz is used to specify the nonzero structure.
 53:            The size of this array is equal to the number of local rows, i.e 'm'.
 54:            For matrices that will be factored, you must leave room for (and set)
 55:            the diagonal entry even if it is zero.
 56: .  o_nz  - number of nonzeros per row in the OFF-DIAGONAL portion of local
 57:            submatrix (same value is used for all local rows).
 58: -  o_nnz - array containing the number of nonzeros in the various rows of the
 59:            OFF-DIAGONAL portion of the local submatrix (possibly different for
 60:            each row) or NULL, if o_nz is used to specify the nonzero
 61:            structure. The size of this array is equal to the number
 62:            of local rows, i.e 'm'.

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

 66:    Level: intermediate

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

 72: .keywords: matrix

 74: .seealso: MatCreate(), MatCreateIS(), MatMPIAIJSetPreallocation(), MatISGetLocalMat()
 75: @*/
 76: PetscErrorCode  MatISSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
 77: {

 83:   PetscTryMethod(B,"MatISSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz));
 84:   return(0);
 85: }

 89: PetscErrorCode  MatISSetPreallocation_IS(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
 90: {
 91:   Mat_IS         *matis = (Mat_IS*)(B->data);
 92:   PetscInt       bs,i,nlocalcols;

 96:   if (!matis->A) {
 97:     SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"You should first call MatSetLocalToGlobalMapping");
 98:   }
 99:   if (!matis->sf) { /* setup SF if not yet created and allocate rootdata and leafdata */
100:     MatISComputeSF_Private(B);
101:   }
102:   if (!d_nnz) {
103:     for (i=0;i<matis->sf_nroots;i++) matis->sf_rootdata[i] = d_nz;
104:   } else {
105:     for (i=0;i<matis->sf_nroots;i++) matis->sf_rootdata[i] = d_nnz[i];
106:   }
107:   if (!o_nnz) {
108:     for (i=0;i<matis->sf_nroots;i++) matis->sf_rootdata[i] += o_nz;
109:   } else {
110:     for (i=0;i<matis->sf_nroots;i++) matis->sf_rootdata[i] += o_nnz[i];
111:   }
112:   PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
113:   MatGetSize(matis->A,NULL,&nlocalcols);
114:   MatGetBlockSize(matis->A,&bs);
115:   PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
116:   for (i=0;i<matis->sf_nleaves;i++) {
117:     matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i],nlocalcols);
118:   }
119:   MatSeqAIJSetPreallocation(matis->A,0,matis->sf_leafdata);
120:   for (i=0;i<matis->sf_nleaves/bs;i++) {
121:     matis->sf_leafdata[i] = matis->sf_leafdata[i*bs]/bs;
122:   }
123:   MatSeqBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata);
124:   for (i=0;i<matis->sf_nleaves/bs;i++) {
125:     matis->sf_leafdata[i] = matis->sf_leafdata[i]-i;
126:   }
127:   MatSeqSBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata);
128:   return(0);
129: }

133: PETSC_EXTERN PetscErrorCode MatISSetMPIXAIJPreallocation_Private(Mat A, Mat B, PetscBool maxreduce)
134: {
135:   Mat_IS          *matis = (Mat_IS*)(A->data);
136:   PetscInt        *my_dnz,*my_onz,*dnz,*onz,*mat_ranges,*row_ownership;
137:   const PetscInt* global_indices;
138:   PetscInt        i,j,bs,rows,cols;
139:   PetscInt        lrows,lcols;
140:   PetscInt        local_rows,local_cols;
141:   PetscMPIInt     nsubdomains;
142:   PetscBool       isdense,issbaij;
143:   PetscErrorCode  ierr;

146:   MPI_Comm_size(PetscObjectComm((PetscObject)A),&nsubdomains);
147:   MatGetSize(A,&rows,&cols);
148:   MatGetBlockSize(A,&bs);
149:   MatGetSize(matis->A,&local_rows,&local_cols);
150:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isdense);
151:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);
152:   ISLocalToGlobalMappingGetIndices(matis->mapping,&global_indices);
153:   if (issbaij) {
154:     MatGetRowUpperTriangular(matis->A);
155:   }
156:   /*
157:      An SF reduce is needed to sum up properly on shared interface dofs.
158:      Note that generally preallocation is not exact, since it overestimates nonzeros
159:   */
160:   if (!matis->sf) { /* setup SF if not yet created and allocate rootdata and leafdata */
161:     MatISComputeSF_Private(A);
162:   }
163:   MatGetLocalSize(A,&lrows,&lcols);
164:   MatPreallocateInitialize(PetscObjectComm((PetscObject)A),lrows,lcols,dnz,onz);
165:   /* All processes need to compute entire row ownership */
166:   PetscMalloc1(rows,&row_ownership);
167:   MatGetOwnershipRanges(A,(const PetscInt**)&mat_ranges);
168:   for (i=0;i<nsubdomains;i++) {
169:     for (j=mat_ranges[i];j<mat_ranges[i+1];j++) {
170:       row_ownership[j]=i;
171:     }
172:   }

174:   /*
175:      my_dnz and my_onz contains exact contribution to preallocation from each local mat
176:      then, they will be summed up properly. This way, preallocation is always sufficient
177:   */
178:   PetscCalloc2(local_rows,&my_dnz,local_rows,&my_onz);
179:   /* preallocation as a MATAIJ */
180:   if (isdense) { /* special case for dense local matrices */
181:     for (i=0;i<local_rows;i++) {
182:       PetscInt index_row = global_indices[i];
183:       for (j=i;j<local_rows;j++) {
184:         PetscInt owner = row_ownership[index_row];
185:         PetscInt index_col = global_indices[j];
186:         if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */
187:           my_dnz[i] += 1;
188:         } else { /* offdiag block */
189:           my_onz[i] += 1;
190:         }
191:         /* same as before, interchanging rows and cols */
192:         if (i != j) {
193:           owner = row_ownership[index_col];
194:           if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) {
195:             my_dnz[j] += 1;
196:           } else {
197:             my_onz[j] += 1;
198:           }
199:         }
200:       }
201:     }
202:   } else {
203:     for (i=0;i<local_rows;i++) {
204:       const PetscInt *cols;
205:       PetscInt       ncols,index_row = global_indices[i];
206:       MatGetRow(matis->A,i,&ncols,&cols,NULL);
207:       for (j=0;j<ncols;j++) {
208:         PetscInt owner = row_ownership[index_row];
209:         PetscInt index_col = global_indices[cols[j]];
210:         if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */
211:           my_dnz[i] += 1;
212:         } else { /* offdiag block */
213:           my_onz[i] += 1;
214:         }
215:         /* same as before, interchanging rows and cols */
216:         if (issbaij && index_col != index_row) {
217:           owner = row_ownership[index_col];
218:           if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) {
219:             my_dnz[cols[j]] += 1;
220:           } else {
221:             my_onz[cols[j]] += 1;
222:           }
223:         }
224:       }
225:       MatRestoreRow(matis->A,i,&ncols,&cols,NULL);
226:     }
227:   }
228:   ISLocalToGlobalMappingRestoreIndices(matis->mapping,&global_indices);
229:   PetscFree(row_ownership);
230:   if (maxreduce) {
231:     PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);
232:     PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);
233:     PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);
234:     PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);
235:   } else {
236:     PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);
237:     PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);
238:     PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);
239:     PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);
240:   }
241:   PetscFree2(my_dnz,my_onz);

243:   /* Resize preallocation if overestimated */
244:   for (i=0;i<lrows;i++) {
245:     dnz[i] = PetscMin(dnz[i],lcols);
246:     onz[i] = PetscMin(onz[i],cols-lcols);
247:   }
248:   /* set preallocation */
249:   MatMPIAIJSetPreallocation(B,0,dnz,0,onz);
250:   for (i=0;i<lrows/bs;i++) {
251:     dnz[i] = dnz[i*bs]/bs;
252:     onz[i] = onz[i*bs]/bs;
253:   }
254:   MatMPIBAIJSetPreallocation(B,bs,0,dnz,0,onz);
255:   MatMPISBAIJSetPreallocation(B,bs,0,dnz,0,onz);
256:   MatPreallocateFinalize(dnz,onz);
257:   if (issbaij) {
258:     MatRestoreRowUpperTriangular(matis->A);
259:   }
260:   return(0);
261: }

265: PetscErrorCode MatISGetMPIXAIJ_IS(Mat mat, MatReuse reuse, Mat *M)
266: {
267:   Mat_IS         *matis = (Mat_IS*)(mat->data);
268:   Mat            local_mat;
269:   /* info on mat */
270:   PetscInt       bs,rows,cols,lrows,lcols;
271:   PetscInt       local_rows,local_cols;
272:   PetscBool      isdense,issbaij,isseqaij;
273:   const PetscInt *global_indices_rows;
274:   PetscMPIInt    nsubdomains;
275:   /* values insertion */
276:   PetscScalar    *array;
277:   /* work */

281:   /* MISSING CHECKS
282:     - rectangular case not covered (it is not allowed by MATIS)
283:   */
284:   /* get info from mat */
285:   MPI_Comm_size(PetscObjectComm((PetscObject)mat),&nsubdomains);
286:   if (nsubdomains == 1) {
287:     if (reuse == MAT_INITIAL_MATRIX) {
288:       MatDuplicate(matis->A,MAT_COPY_VALUES,&(*M));
289:     } else {
290:       MatCopy(matis->A,*M,SAME_NONZERO_PATTERN);
291:     }
292:     return(0);
293:   }
294:   MatGetSize(mat,&rows,&cols);
295:   MatGetBlockSize(mat,&bs);
296:   MatGetLocalSize(mat,&lrows,&lcols);
297:   MatGetSize(matis->A,&local_rows,&local_cols);
298:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isdense);
299:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQAIJ,&isseqaij);
300:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);

302:   /* map indices of local mat rows to global values */
303:   ISLocalToGlobalMappingGetIndices(matis->mapping,&global_indices_rows);

305:   if (reuse == MAT_INITIAL_MATRIX) {
306:     MatType     new_mat_type;
307:     PetscBool   issbaij_red;

309:     /* determining new matrix type */
310:     MPI_Allreduce(&issbaij,&issbaij_red,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)mat));
311:     if (issbaij_red) {
312:       new_mat_type = MATSBAIJ;
313:     } else {
314:       if (bs>1) {
315:         new_mat_type = MATBAIJ;
316:       } else {
317:         new_mat_type = MATAIJ;
318:       }
319:     }

321:     MatCreate(PetscObjectComm((PetscObject)mat),M);
322:     MatSetSizes(*M,lrows,lcols,rows,cols);
323:     MatSetBlockSize(*M,bs);
324:     MatSetType(*M,new_mat_type);
325:     MatISSetMPIXAIJPreallocation_Private(mat,*M,PETSC_FALSE);
326:   } else {
327:     PetscInt mbs,mrows,mcols,mlrows,mlcols;
328:     /* some checks */
329:     MatGetBlockSize(*M,&mbs);
330:     MatGetSize(*M,&mrows,&mcols);
331:     MatGetLocalSize(*M,&mlrows,&mlcols);
332:     if (mrows != rows) {
333:       SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of rows (%d != %d)",rows,mrows);
334:     }
335:     if (mcols != cols) {
336:       SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of cols (%d != %d)",cols,mcols);
337:     }
338:     if (mlrows != lrows) {
339:       SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local rows (%d != %d)",lrows,mlrows);
340:     }
341:     if (mlcols != lcols) {
342:       SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local cols (%d != %d)",lcols,mlcols);
343:     }
344:     if (mbs != bs) {
345:       SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong block size (%d != %d)",bs,mbs);
346:     }
347:     MatZeroEntries(*M);
348:   }

350:   if (issbaij) {
351:     MatConvert(matis->A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&local_mat);
352:   } else {
353:     PetscObjectReference((PetscObject)matis->A);
354:     local_mat = matis->A;
355:   }

357:   /* Set values */
358:   if (isdense) { /* special case for dense local matrices */
359:     MatSetOption(*M,MAT_ROW_ORIENTED,PETSC_FALSE);
360:     MatDenseGetArray(local_mat,&array);
361:     MatSetValues(*M,local_rows,global_indices_rows,local_cols,global_indices_rows,array,ADD_VALUES);
362:     MatDenseRestoreArray(local_mat,&array);
363:   } else if (isseqaij) {
364:     PetscInt  i,nvtxs,*xadj,*adjncy,*global_indices_cols;
365:     PetscBool done;

367:     MatGetRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);
368:     if (!done) {
369:       SETERRQ1(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ called in %s\n",__FUNCT__);
370:     }
371:     MatSeqAIJGetArray(local_mat,&array);
372:     PetscMalloc1(xadj[nvtxs],&global_indices_cols);
373:     ISLocalToGlobalMappingApply(matis->mapping,xadj[nvtxs],adjncy,global_indices_cols);
374:     for (i=0;i<nvtxs;i++) {
375:       PetscInt    row,*cols,ncols;
376:       PetscScalar *mat_vals;

378:       row = global_indices_rows[i];
379:       ncols = xadj[i+1]-xadj[i];
380:       cols = global_indices_cols+xadj[i];
381:       mat_vals = array+xadj[i];
382:       MatSetValues(*M,1,&row,ncols,cols,mat_vals,ADD_VALUES);
383:     }
384:     MatRestoreRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);
385:     if (!done) {
386:       SETERRQ1(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ called in %s\n",__FUNCT__);
387:     }
388:     MatSeqAIJRestoreArray(local_mat,&array);
389:     PetscFree(global_indices_cols);
390:   } else { /* very basic values insertion for all other matrix types */
391:     PetscInt i,*global_indices_cols;

393:     PetscMalloc1(local_cols,&global_indices_cols);
394:     for (i=0;i<local_rows;i++) {
395:       PetscInt       j;
396:       const PetscInt *local_indices;

398:       MatGetRow(local_mat,i,&j,&local_indices,(const PetscScalar**)&array);
399:       ISLocalToGlobalMappingApply(matis->mapping,j,local_indices,global_indices_cols);
400:       MatSetValues(*M,1,&global_indices_rows[i],j,global_indices_cols,array,ADD_VALUES);
401:       MatRestoreRow(local_mat,i,&j,&local_indices,(const PetscScalar**)&array);
402:     }
403:     PetscFree(global_indices_cols);
404:   }
405:   MatAssemblyBegin(*M,MAT_FINAL_ASSEMBLY);
406:   ISLocalToGlobalMappingRestoreIndices(matis->mapping,&global_indices_rows);
407:   MatDestroy(&local_mat);
408:   MatAssemblyEnd(*M,MAT_FINAL_ASSEMBLY);
409:   if (isdense) {
410:     MatSetOption(*M,MAT_ROW_ORIENTED,PETSC_TRUE);
411:   }
412:   return(0);
413: }

417: /*@
418:     MatISGetMPIXAIJ - Converts MATIS matrix into a parallel AIJ format

420:   Input Parameter:
421: .  mat - the matrix (should be of type MATIS)
422: .  reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX

424:   Output Parameter:
425: .  newmat - the matrix in AIJ format

427:   Level: developer

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

431: .seealso: MATIS
432: @*/
433: PetscErrorCode MatISGetMPIXAIJ(Mat mat, MatReuse reuse, Mat *newmat)
434: {

441:   if (reuse != MAT_INITIAL_MATRIX) {
444:     if (mat == *newmat) {
445:       SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse the same matrix");
446:     }
447:   }
448:   PetscUseMethod(mat,"MatISGetMPIXAIJ_C",(Mat,MatReuse,Mat*),(mat,reuse,newmat));
449:   return(0);
450: }

454: PetscErrorCode MatDuplicate_IS(Mat mat,MatDuplicateOption op,Mat *newmat)
455: {
457:   Mat_IS         *matis = (Mat_IS*)(mat->data);
458:   PetscInt       bs,m,n,M,N;
459:   Mat            B,localmat;

462:   MatGetBlockSize(mat,&bs);
463:   MatGetSize(mat,&M,&N);
464:   MatGetLocalSize(mat,&m,&n);
465:   MatCreateIS(PetscObjectComm((PetscObject)mat),bs,m,n,M,N,matis->mapping,&B);
466:   MatDuplicate(matis->A,op,&localmat);
467:   MatISSetLocalMat(B,localmat);
468:   MatDestroy(&localmat);
469:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
470:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
471:   *newmat = B;
472:   return(0);
473: }

477: PetscErrorCode MatIsHermitian_IS(Mat A,PetscReal tol,PetscBool  *flg)
478: {
480:   Mat_IS         *matis = (Mat_IS*)A->data;
481:   PetscBool      local_sym;

484:   MatIsHermitian(matis->A,tol,&local_sym);
485:   MPI_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
486:   return(0);
487: }

491: PetscErrorCode MatIsSymmetric_IS(Mat A,PetscReal tol,PetscBool  *flg)
492: {
494:   Mat_IS         *matis = (Mat_IS*)A->data;
495:   PetscBool      local_sym;

498:   MatIsSymmetric(matis->A,tol,&local_sym);
499:   MPI_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
500:   return(0);
501: }

505: PetscErrorCode MatDestroy_IS(Mat A)
506: {
508:   Mat_IS         *b = (Mat_IS*)A->data;

511:   MatDestroy(&b->A);
512:   VecScatterDestroy(&b->ctx);
513:   VecDestroy(&b->x);
514:   VecDestroy(&b->y);
515:   ISLocalToGlobalMappingDestroy(&b->mapping);
516:   PetscSFDestroy(&b->sf);
517:   PetscFree2(b->sf_rootdata,b->sf_leafdata);
518:   PetscFree(A->data);
519:   PetscObjectChangeTypeName((PetscObject)A,0);
520:   PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL);
521:   PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",NULL);
522:   PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",NULL);
523:   PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",NULL);
524:   return(0);
525: }

529: PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y)
530: {
532:   Mat_IS         *is  = (Mat_IS*)A->data;
533:   PetscScalar    zero = 0.0;

536:   /*  scatter the global vector x into the local work vector */
537:   VecScatterBegin(is->ctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
538:   VecScatterEnd(is->ctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);

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

543:   /* scatter product back into global memory */
544:   VecSet(y,zero);
545:   VecScatterBegin(is->ctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
546:   VecScatterEnd(is->ctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
547:   return(0);
548: }

552: PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
553: {
554:   Vec            temp_vec;

558:   if (v3 != v2) {
559:     MatMult(A,v1,v3);
560:     VecAXPY(v3,1.0,v2);
561:   } else {
562:     VecDuplicate(v2,&temp_vec);
563:     MatMult(A,v1,temp_vec);
564:     VecAXPY(temp_vec,1.0,v2);
565:     VecCopy(temp_vec,v3);
566:     VecDestroy(&temp_vec);
567:   }
568:   return(0);
569: }

573: PetscErrorCode MatMultTranspose_IS(Mat A,Vec x,Vec y)
574: {
575:   Mat_IS         *is = (Mat_IS*)A->data;

579:   /*  scatter the global vector x into the local work vector */
580:   VecScatterBegin(is->ctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
581:   VecScatterEnd(is->ctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);

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

586:   /* scatter product back into global vector */
587:   VecSet(y,0);
588:   VecScatterBegin(is->ctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
589:   VecScatterEnd(is->ctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
590:   return(0);
591: }

595: PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
596: {
597:   Vec            temp_vec;

601:   if (v3 != v2) {
602:     MatMultTranspose(A,v1,v3);
603:     VecAXPY(v3,1.0,v2);
604:   } else {
605:     VecDuplicate(v2,&temp_vec);
606:     MatMultTranspose(A,v1,temp_vec);
607:     VecAXPY(temp_vec,1.0,v2);
608:     VecCopy(temp_vec,v3);
609:     VecDestroy(&temp_vec);
610:   }
611:   return(0);
612: }

616: PetscErrorCode MatView_IS(Mat A,PetscViewer viewer)
617: {
618:   Mat_IS         *a = (Mat_IS*)A->data;
620:   PetscViewer    sviewer;

623:   PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
624:   PetscViewerGetSingleton(viewer,&sviewer);
625:   MatView(a->A,sviewer);
626:   PetscViewerRestoreSingleton(viewer,&sviewer);
627:   return(0);
628: }

632: PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping)
633: {
635:   PetscInt       n,bs;
636:   Mat_IS         *is = (Mat_IS*)A->data;
637:   IS             from,to;
638:   Vec            global;

642:   if (rmapping != cmapping) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"MATIS requires the row and column mappings to be identical");
643:   PetscLayoutSetUp(A->rmap);
644:   PetscLayoutSetUp(A->cmap);
645:   if (is->mapping) { /* Currenly destroys the objects that will be created by this routine. Is there anything else that should be checked? */
646:     ISLocalToGlobalMappingDestroy(&is->mapping);
647:     VecDestroy(&is->x);
648:     VecDestroy(&is->y);
649:     VecScatterDestroy(&is->ctx);
650:     MatDestroy(&is->A);
651:     PetscSFDestroy(&is->sf);
652:     PetscFree2(is->sf_rootdata,is->sf_leafdata);
653:   }
654:   PetscObjectReference((PetscObject)rmapping);
655:   ISLocalToGlobalMappingDestroy(&is->mapping);
656:   is->mapping = rmapping;
657: /*
658:   PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping);
659:   PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping);
660: */

662:   /* Create the local matrix A */
663:   ISLocalToGlobalMappingGetSize(rmapping,&n);
664:   ISLocalToGlobalMappingGetBlockSize(rmapping,&bs);
665:   MatCreate(PETSC_COMM_SELF,&is->A);
666:   if (bs > 1) {
667:     MatSetType(is->A,MATSEQBAIJ);
668:   } else {
669:     MatSetType(is->A,MATSEQAIJ);
670:   }
671:   MatSetSizes(is->A,n,n,n,n);
672:   MatSetBlockSize(is->A,bs);
673:   MatSetOptionsPrefix(is->A,((PetscObject)A)->prefix);
674:   MatAppendOptionsPrefix(is->A,"is_");
675:   MatSetFromOptions(is->A);

677:   /* Create the local work vectors */
678:   VecCreate(PETSC_COMM_SELF,&is->x);
679:   VecSetBlockSize(is->x,bs);
680:   VecSetSizes(is->x,n,n);
681:   VecSetOptionsPrefix(is->x,((PetscObject)A)->prefix);
682:   VecAppendOptionsPrefix(is->x,"is_");
683:   VecSetFromOptions(is->x);
684:   VecDuplicate(is->x,&is->y);

686:   /* setup the global to local scatter */
687:   ISCreateStride(PETSC_COMM_SELF,n,0,1,&to);
688:   ISLocalToGlobalMappingApplyIS(rmapping,to,&from);
689:   MatCreateVecs(A,&global,NULL);
690:   VecScatterCreate(global,from,is->x,to,&is->ctx);
691:   VecDestroy(&global);
692:   ISDestroy(&to);
693:   ISDestroy(&from);
694:   return(0);
695: }

699: PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv)
700: {
701:   Mat_IS         *is = (Mat_IS*)mat->data;
702:   PetscInt       rows_l[2048],cols_l[2048];

706: #if defined(PETSC_USE_DEBUG)
707:   if (m > 2048 || n > 2048) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of row/column indices must be <= 2048: they are %D %D",m,n);
708: #endif
709:   ISG2LMapApply(is->mapping,m,rows,rows_l);
710:   ISG2LMapApply(is->mapping,n,cols,cols_l);
711:   MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);
712:   return(0);
713: }

715: #undef ISG2LMapSetUp
716: #undef ISG2LMapApply

720: PetscErrorCode MatSetValuesLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
721: {
723:   Mat_IS         *is = (Mat_IS*)A->data;

726:   MatSetValues(is->A,m,rows,n,cols,values,addv);
727:   return(0);
728: }

732: PetscErrorCode MatSetValuesBlockedLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
733: {
735:   Mat_IS         *is = (Mat_IS*)A->data;

738:   MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv);
739:   return(0);
740: }

744: PetscErrorCode MatZeroRows_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
745: {
746:   Mat_IS         *is = (Mat_IS*)A->data;
747:   PetscInt       n_l = 0, *rows_l = NULL;

751:   if (x && b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support");
752:   if (n) {
753:     PetscMalloc1(n,&rows_l);
754:     ISGlobalToLocalMappingApply(is->mapping,IS_GTOLM_DROP,n,rows,&n_l,rows_l);
755:   }
756:   MatZeroRowsLocal(A,n_l,rows_l,diag,x,b);
757:   PetscFree(rows_l);
758:   return(0);
759: }

763: PetscErrorCode MatZeroRowsLocal_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
764: {
765:   Mat_IS         *is = (Mat_IS*)A->data;
767:   PetscInt       i;
768:   PetscScalar    *array;

771:   if (x && b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support");
772:   {
773:     /*
774:        Set up is->x as a "counting vector". This is in order to MatMult_IS
775:        work properly in the interface nodes.
776:     */
777:     Vec         counter;
778:     PetscScalar one=1.0, zero=0.0;
779:     MatCreateVecs(A,&counter,NULL);
780:     VecSet(counter,zero);
781:     VecSet(is->x,one);
782:     VecScatterBegin(is->ctx,is->x,counter,ADD_VALUES,SCATTER_REVERSE);
783:     VecScatterEnd  (is->ctx,is->x,counter,ADD_VALUES,SCATTER_REVERSE);
784:     VecScatterBegin(is->ctx,counter,is->x,INSERT_VALUES,SCATTER_FORWARD);
785:     VecScatterEnd  (is->ctx,counter,is->x,INSERT_VALUES,SCATTER_FORWARD);
786:     VecDestroy(&counter);
787:   }
788:   if (!n) {
789:     is->pure_neumann = PETSC_TRUE;
790:   } else {
791:     is->pure_neumann = PETSC_FALSE;

793:     VecGetArray(is->x,&array);
794:     MatZeroRows(is->A,n,rows,diag,0,0);
795:     for (i=0; i<n; i++) {
796:       MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);
797:     }
798:     MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);
799:     MatAssemblyEnd  (is->A,MAT_FINAL_ASSEMBLY);
800:     VecRestoreArray(is->x,&array);
801:   }
802:   return(0);
803: }

807: PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type)
808: {
809:   Mat_IS         *is = (Mat_IS*)A->data;

813:   MatAssemblyBegin(is->A,type);
814:   return(0);
815: }

819: PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type)
820: {
821:   Mat_IS         *is = (Mat_IS*)A->data;

825:   MatAssemblyEnd(is->A,type);
826:   return(0);
827: }

831: PetscErrorCode MatISGetLocalMat_IS(Mat mat,Mat *local)
832: {
833:   Mat_IS *is = (Mat_IS*)mat->data;

836:   *local = is->A;
837:   return(0);
838: }

842: /*@
843:     MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix.

845:   Input Parameter:
846: .  mat - the matrix

848:   Output Parameter:
849: .  local - the local matrix

851:   Level: advanced

853:   Notes:
854:     This can be called if you have precomputed the nonzero structure of the
855:   matrix and want to provide it to the inner matrix object to improve the performance
856:   of the MatSetValues() operation.

858: .seealso: MATIS
859: @*/
860: PetscErrorCode MatISGetLocalMat(Mat mat,Mat *local)
861: {

867:   PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat*),(mat,local));
868:   return(0);
869: }

873: PetscErrorCode MatISSetLocalMat_IS(Mat mat,Mat local)
874: {
875:   Mat_IS         *is = (Mat_IS*)mat->data;
876:   PetscInt       nrows,ncols,orows,ocols;

880:   if (is->A) {
881:     MatGetSize(is->A,&orows,&ocols);
882:     MatGetSize(local,&nrows,&ncols);
883:     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)\n",orows,ocols,nrows,ncols);
884:   }
885:   PetscObjectReference((PetscObject)local);
886:   MatDestroy(&is->A);
887:   is->A = local;
888:   return(0);
889: }

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

896:   Input Parameter:
897: .  mat - the matrix
898: .  local - the local matrix

900:   Output Parameter:

902:   Level: advanced

904:   Notes:
905:     This can be called if you have precomputed the local matrix and
906:   want to provide it to the matrix object MATIS.

908: .seealso: MATIS
909: @*/
910: PetscErrorCode MatISSetLocalMat(Mat mat,Mat local)
911: {

917:   PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local));
918:   return(0);
919: }

923: PetscErrorCode MatZeroEntries_IS(Mat A)
924: {
925:   Mat_IS         *a = (Mat_IS*)A->data;

929:   MatZeroEntries(a->A);
930:   return(0);
931: }

935: PetscErrorCode MatScale_IS(Mat A,PetscScalar a)
936: {
937:   Mat_IS         *is = (Mat_IS*)A->data;

941:   MatScale(is->A,a);
942:   return(0);
943: }

947: PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v)
948: {
949:   Mat_IS         *is = (Mat_IS*)A->data;

953:   /* get diagonal of the local matrix */
954:   MatGetDiagonal(is->A,is->x);

956:   /* scatter diagonal back into global vector */
957:   VecSet(v,0);
958:   VecScatterBegin(is->ctx,is->x,v,ADD_VALUES,SCATTER_REVERSE);
959:   VecScatterEnd  (is->ctx,is->x,v,ADD_VALUES,SCATTER_REVERSE);
960:   return(0);
961: }

965: PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg)
966: {
967:   Mat_IS         *a = (Mat_IS*)A->data;

971:   MatSetOption(a->A,op,flg);
972:   return(0);
973: }

977: /*@
978:     MatCreateIS - Creates a "process" unassmembled matrix, it is assembled on each
979:        process but not across processes.

981:    Input Parameters:
982: +     comm - MPI communicator that will share the matrix
983: .     bs - local and global block size of the matrix
984: .     m,n,M,N - local and/or global sizes of the left and right vector used in matrix vector products
985: -     map - mapping that defines the global number for each local number

987:    Output Parameter:
988: .    A - the resulting matrix

990:    Level: advanced

992:    Notes: See MATIS for more details
993:           m and n are NOT related to the size of the map, they are the size of the part of the vector owned
994:           by that process. m + nghosts (or n + nghosts) is the length of map since map maps all local points
995:           plus the ghost points to global indices.

997: .seealso: MATIS, MatSetLocalToGlobalMapping()
998: @*/
999: PetscErrorCode  MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping map,Mat *A)
1000: {

1004:   MatCreate(comm,A);
1005:   MatSetBlockSize(*A,bs);
1006:   MatSetSizes(*A,m,n,M,N);
1007:   MatSetType(*A,MATIS);
1008:   MatSetUp(*A);
1009:   MatSetLocalToGlobalMapping(*A,map,map);
1010:   return(0);
1011: }

1013: /*MC
1014:    MATIS - MATIS = "is" - A matrix type to be used for using the non-overlapping domain decomposition type preconditioners (e.g. PCBDDC).
1015:    This stores the matrices in globally unassembled form. Each processor
1016:    assembles only its local Neumann problem and the parallel matrix vector
1017:    product is handled "implicitly".

1019:    Operations Provided:
1020: +  MatMult()
1021: .  MatMultAdd()
1022: .  MatMultTranspose()
1023: .  MatMultTransposeAdd()
1024: .  MatZeroEntries()
1025: .  MatSetOption()
1026: .  MatZeroRows()
1027: .  MatZeroRowsLocal()
1028: .  MatSetValues()
1029: .  MatSetValuesLocal()
1030: .  MatScale()
1031: .  MatGetDiagonal()
1032: -  MatSetLocalToGlobalMapping()

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

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

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

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

1044:   Level: advanced

1046: .seealso: Mat, MatISGetLocalMat(), MatSetLocalToGlobalMapping(), MatISSetPreallocation(), PCBDDC

1048: M*/

1052: PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A)
1053: {
1055:   Mat_IS         *b;

1058:   PetscNewLog(A,&b);
1059:   A->data = (void*)b;
1060:   PetscMemzero(A->ops,sizeof(struct _MatOps));

1062:   A->ops->mult                    = MatMult_IS;
1063:   A->ops->multadd                 = MatMultAdd_IS;
1064:   A->ops->multtranspose           = MatMultTranspose_IS;
1065:   A->ops->multtransposeadd        = MatMultTransposeAdd_IS;
1066:   A->ops->destroy                 = MatDestroy_IS;
1067:   A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS;
1068:   A->ops->setvalues               = MatSetValues_IS;
1069:   A->ops->setvalueslocal          = MatSetValuesLocal_IS;
1070:   A->ops->setvaluesblockedlocal   = MatSetValuesBlockedLocal_IS;
1071:   A->ops->zerorows                = MatZeroRows_IS;
1072:   A->ops->zerorowslocal           = MatZeroRowsLocal_IS;
1073:   A->ops->assemblybegin           = MatAssemblyBegin_IS;
1074:   A->ops->assemblyend             = MatAssemblyEnd_IS;
1075:   A->ops->view                    = MatView_IS;
1076:   A->ops->zeroentries             = MatZeroEntries_IS;
1077:   A->ops->scale                   = MatScale_IS;
1078:   A->ops->getdiagonal             = MatGetDiagonal_IS;
1079:   A->ops->setoption               = MatSetOption_IS;
1080:   A->ops->ishermitian             = MatIsHermitian_IS;
1081:   A->ops->issymmetric             = MatIsSymmetric_IS;
1082:   A->ops->duplicate               = MatDuplicate_IS;

1084:   /* zeros pointer */
1085:   b->A           = 0;
1086:   b->ctx         = 0;
1087:   b->x           = 0;
1088:   b->y           = 0;
1089:   b->mapping     = 0;
1090:   b->sf          = 0;
1091:   b->sf_rootdata = 0;
1092:   b->sf_leafdata = 0;

1094:   /* special MATIS functions */
1095:   PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS);
1096:   PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS);
1097:   PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",MatISGetMPIXAIJ_IS);
1098:   PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",MatISSetPreallocation_IS);
1099:   PetscObjectChangeTypeName((PetscObject)A,MATIS);
1100:   return(0);
1101: }