Actual source code: matis.c

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

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

144:   MPI_Comm_size(PetscObjectComm((PetscObject)A),&nsubdomains);
145:   MatGetSize(A,&rows,&cols);
146:   MatGetBlockSize(A,&bs);
147:   MatGetSize(matis->A,&local_rows,&local_cols);
148:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isdense);
149:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);
150:   ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&global_indices_r);
151:   if (A->rmap->mapping != A->cmap->mapping) {
152:     ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&global_indices_c);
153:   } else {
154:     global_indices_c = global_indices_r;
155:   }

157:   if (issbaij) {
158:     MatGetRowUpperTriangular(matis->A);
159:   }
160:   /*
161:      An SF reduce is needed to sum up properly on shared rows.
162:      Note that generally preallocation is not exact, since it overestimates nonzeros
163:   */
164:   if (!matis->sf) { /* setup SF if not yet created and allocate rootdata and leafdata */
165:     MatISComputeSF_Private(A);
166:   }
167:   MatGetLocalSize(A,&lrows,&lcols);
168:   MatPreallocateInitialize(PetscObjectComm((PetscObject)A),lrows,lcols,dnz,onz);
169:   /* All processes need to compute entire row ownership */
170:   PetscMalloc1(rows,&row_ownership);
171:   MatGetOwnershipRanges(A,(const PetscInt**)&mat_ranges);
172:   for (i=0;i<nsubdomains;i++) {
173:     for (j=mat_ranges[i];j<mat_ranges[i+1];j++) {
174:       row_ownership[j] = i;
175:     }
176:   }

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

238:   /* Reduce my_dnz and my_onz */
239:   if (maxreduce) {
240:     PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);
241:     PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);
242:     PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);
243:     PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);
244:   } else {
245:     PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);
246:     PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);
247:     PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);
248:     PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);
249:   }
250:   PetscFree2(my_dnz,my_onz);

252:   /* Resize preallocation if overestimated */
253:   for (i=0;i<lrows;i++) {
254:     dnz[i] = PetscMin(dnz[i],lcols);
255:     onz[i] = PetscMin(onz[i],cols-lcols);
256:   }
257:   /* set preallocation */
258:   MatMPIAIJSetPreallocation(B,0,dnz,0,onz);
259:   for (i=0;i<lrows/bs;i++) {
260:     dnz[i] = dnz[i*bs]/bs;
261:     onz[i] = onz[i*bs]/bs;
262:   }
263:   MatMPIBAIJSetPreallocation(B,bs,0,dnz,0,onz);
264:   MatMPISBAIJSetPreallocation(B,bs,0,dnz,0,onz);
265:   MatPreallocateFinalize(dnz,onz);
266:   if (issbaij) {
267:     MatRestoreRowUpperTriangular(matis->A);
268:   }
269:   return(0);
270: }

274: PetscErrorCode MatISGetMPIXAIJ_IS(Mat mat, MatReuse reuse, Mat *M)
275: {
276:   Mat_IS         *matis = (Mat_IS*)(mat->data);
277:   Mat            local_mat;
278:   /* info on mat */
279:   PetscInt       bs,rows,cols,lrows,lcols;
280:   PetscInt       local_rows,local_cols;
281:   PetscBool      isdense,issbaij,isseqaij;
282:   PetscMPIInt    nsubdomains;
283:   /* values insertion */
284:   PetscScalar    *array;
285:   /* work */

289:   /* get info from mat */
290:   MPI_Comm_size(PetscObjectComm((PetscObject)mat),&nsubdomains);
291:   if (nsubdomains == 1) {
292:     if (reuse == MAT_INITIAL_MATRIX) {
293:       MatDuplicate(matis->A,MAT_COPY_VALUES,&(*M));
294:     } else {
295:       MatCopy(matis->A,*M,SAME_NONZERO_PATTERN);
296:     }
297:     return(0);
298:   }
299:   MatGetSize(mat,&rows,&cols);
300:   MatGetBlockSize(mat,&bs);
301:   MatGetLocalSize(mat,&lrows,&lcols);
302:   MatGetSize(matis->A,&local_rows,&local_cols);
303:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isdense);
304:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQAIJ,&isseqaij);
305:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);

307:   if (reuse == MAT_INITIAL_MATRIX) {
308:     MatType     new_mat_type;
309:     PetscBool   issbaij_red;

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

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

342:   if (issbaij) {
343:     MatConvert(matis->A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&local_mat);
344:   } else {
345:     PetscObjectReference((PetscObject)matis->A);
346:     local_mat = matis->A;
347:   }

349:   /* Set values */
350:   MatSetLocalToGlobalMapping(*M,mat->rmap->mapping,mat->cmap->mapping);
351:   if (isdense) { /* special case for dense local matrices */
352:     PetscInt i,*dummy_rows,*dummy_cols;

354:     if (local_rows != local_cols) {
355:       PetscMalloc2(local_rows,&dummy_rows,local_cols,&dummy_cols);
356:       for (i=0;i<local_rows;i++) dummy_rows[i] = i;
357:       for (i=0;i<local_cols;i++) dummy_cols[i] = i;
358:     } else {
359:       PetscMalloc1(local_rows,&dummy_rows);
360:       for (i=0;i<local_rows;i++) dummy_rows[i] = i;
361:       dummy_cols = dummy_rows;
362:     }
363:     MatSetOption(*M,MAT_ROW_ORIENTED,PETSC_FALSE);
364:     MatDenseGetArray(local_mat,&array);
365:     MatSetValuesLocal(*M,local_rows,dummy_rows,local_cols,dummy_cols,array,ADD_VALUES);
366:     MatDenseRestoreArray(local_mat,&array);
367:     if (dummy_rows != dummy_cols) {
368:       PetscFree2(dummy_rows,dummy_cols);
369:     } else {
370:       PetscFree(dummy_rows);
371:     }
372:   } else if (isseqaij) {
373:     PetscInt  i,nvtxs,*xadj,*adjncy;
374:     PetscBool done;

376:     MatGetRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);
377:     if (!done) SETERRQ1(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ called in %s\n",__FUNCT__);
378:     MatSeqAIJGetArray(local_mat,&array);
379:     for (i=0;i<nvtxs;i++) {
380:       MatSetValuesLocal(*M,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],array+xadj[i],ADD_VALUES);
381:     }
382:     MatRestoreRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);
383:     if (!done) SETERRQ1(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ called in %s\n",__FUNCT__);
384:     MatSeqAIJRestoreArray(local_mat,&array);
385:   } else { /* very basic values insertion for all other matrix types */
386:     PetscInt i;

388:     for (i=0;i<local_rows;i++) {
389:       PetscInt       j;
390:       const PetscInt *local_indices_cols;

392:       MatGetRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);
393:       MatSetValuesLocal(*M,1,&i,j,local_indices_cols,array,ADD_VALUES);
394:       MatRestoreRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);
395:     }
396:   }
397:   MatAssemblyBegin(*M,MAT_FINAL_ASSEMBLY);
398:   MatDestroy(&local_mat);
399:   MatAssemblyEnd(*M,MAT_FINAL_ASSEMBLY);
400:   if (isdense) {
401:     MatSetOption(*M,MAT_ROW_ORIENTED,PETSC_TRUE);
402:   }
403:   return(0);
404: }

408: /*@
409:     MatISGetMPIXAIJ - Converts MATIS matrix into a parallel AIJ format

411:   Input Parameter:
412: .  mat - the matrix (should be of type MATIS)
413: .  reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX

415:   Output Parameter:
416: .  newmat - the matrix in AIJ format

418:   Level: developer

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

422: .seealso: MATIS
423: @*/
424: PetscErrorCode MatISGetMPIXAIJ(Mat mat, MatReuse reuse, Mat *newmat)
425: {

432:   if (reuse != MAT_INITIAL_MATRIX) {
435:     if (mat == *newmat) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse the same matrix");
436:   }
437:   PetscUseMethod(mat,"MatISGetMPIXAIJ_C",(Mat,MatReuse,Mat*),(mat,reuse,newmat));
438:   return(0);
439: }

443: PetscErrorCode MatDuplicate_IS(Mat mat,MatDuplicateOption op,Mat *newmat)
444: {
446:   Mat_IS         *matis = (Mat_IS*)(mat->data);
447:   PetscInt       bs,m,n,M,N;
448:   Mat            B,localmat;

451:   MatGetBlockSize(mat,&bs);
452:   MatGetSize(mat,&M,&N);
453:   MatGetLocalSize(mat,&m,&n);
454:   MatCreateIS(PetscObjectComm((PetscObject)mat),bs,m,n,M,N,mat->rmap->mapping,mat->cmap->mapping,&B);
455:   MatDuplicate(matis->A,op,&localmat);
456:   MatISSetLocalMat(B,localmat);
457:   MatDestroy(&localmat);
458:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
459:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
460:   *newmat = B;
461:   return(0);
462: }

466: PetscErrorCode MatIsHermitian_IS(Mat A,PetscReal tol,PetscBool  *flg)
467: {
469:   Mat_IS         *matis = (Mat_IS*)A->data;
470:   PetscBool      local_sym;

473:   MatIsHermitian(matis->A,tol,&local_sym);
474:   MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
475:   return(0);
476: }

480: PetscErrorCode MatIsSymmetric_IS(Mat A,PetscReal tol,PetscBool  *flg)
481: {
483:   Mat_IS         *matis = (Mat_IS*)A->data;
484:   PetscBool      local_sym;

487:   MatIsSymmetric(matis->A,tol,&local_sym);
488:   MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
489:   return(0);
490: }

494: PetscErrorCode MatDestroy_IS(Mat A)
495: {
497:   Mat_IS         *b = (Mat_IS*)A->data;

500:   MatDestroy(&b->A);
501:   VecScatterDestroy(&b->cctx);
502:   VecScatterDestroy(&b->rctx);
503:   VecDestroy(&b->x);
504:   VecDestroy(&b->y);
505:   PetscSFDestroy(&b->sf);
506:   PetscFree2(b->sf_rootdata,b->sf_leafdata);
507:   PetscFree(A->data);
508:   PetscObjectChangeTypeName((PetscObject)A,0);
509:   PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL);
510:   PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",NULL);
511:   PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",NULL);
512:   PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",NULL);
513:   return(0);
514: }

518: PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y)
519: {
521:   Mat_IS         *is  = (Mat_IS*)A->data;
522:   PetscScalar    zero = 0.0;

525:   /*  scatter the global vector x into the local work vector */
526:   VecScatterBegin(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
527:   VecScatterEnd(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);

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

532:   /* scatter product back into global memory */
533:   VecSet(y,zero);
534:   VecScatterBegin(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
535:   VecScatterEnd(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
536:   return(0);
537: }

541: PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
542: {
543:   Vec            temp_vec;

547:   if (v3 != v2) {
548:     MatMult(A,v1,v3);
549:     VecAXPY(v3,1.0,v2);
550:   } else {
551:     VecDuplicate(v2,&temp_vec);
552:     MatMult(A,v1,temp_vec);
553:     VecAXPY(temp_vec,1.0,v2);
554:     VecCopy(temp_vec,v3);
555:     VecDestroy(&temp_vec);
556:   }
557:   return(0);
558: }

562: PetscErrorCode MatMultTranspose_IS(Mat A,Vec y,Vec x)
563: {
564:   Mat_IS         *is = (Mat_IS*)A->data;

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

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

575:   /* scatter product back into global vector */
576:   VecSet(x,0);
577:   VecScatterBegin(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);
578:   VecScatterEnd(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);
579:   return(0);
580: }

584: PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
585: {
586:   Vec            temp_vec;

590:   if (v3 != v2) {
591:     MatMultTranspose(A,v1,v3);
592:     VecAXPY(v3,1.0,v2);
593:   } else {
594:     VecDuplicate(v2,&temp_vec);
595:     MatMultTranspose(A,v1,temp_vec);
596:     VecAXPY(temp_vec,1.0,v2);
597:     VecCopy(temp_vec,v3);
598:     VecDestroy(&temp_vec);
599:   }
600:   return(0);
601: }

605: PetscErrorCode MatView_IS(Mat A,PetscViewer viewer)
606: {
607:   Mat_IS         *a = (Mat_IS*)A->data;
609:   PetscViewer    sviewer;

612:   PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
613:   MatView(a->A,sviewer);
614:   PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
615:   return(0);
616: }

620: PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping)
621: {
623:   PetscInt       nr,rbs,nc,cbs;
624:   Mat_IS         *is = (Mat_IS*)A->data;
625:   IS             from,to;
626:   Vec            cglobal,rglobal;

631:   /* Destroy any previous data */
632:   VecDestroy(&is->x);
633:   VecDestroy(&is->y);
634:   VecScatterDestroy(&is->rctx);
635:   VecScatterDestroy(&is->cctx);
636:   MatDestroy(&is->A);
637:   PetscSFDestroy(&is->sf);
638:   PetscFree2(is->sf_rootdata,is->sf_leafdata);

640:   /* Setup Layout and set local to global maps */
641:   PetscLayoutSetUp(A->rmap);
642:   PetscLayoutSetUp(A->cmap);
643:   PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping);
644:   PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping);

646:   /* Create the local matrix A */
647:   ISLocalToGlobalMappingGetSize(rmapping,&nr);
648:   ISLocalToGlobalMappingGetBlockSize(rmapping,&rbs);
649:   ISLocalToGlobalMappingGetSize(cmapping,&nc);
650:   ISLocalToGlobalMappingGetBlockSize(cmapping,&cbs);
651:   MatCreate(PETSC_COMM_SELF,&is->A);
652:   MatSetType(is->A,MATAIJ);
653:   MatSetSizes(is->A,nr,nc,nr,nc);
654:   MatSetBlockSizes(is->A,rbs,cbs);
655:   MatSetOptionsPrefix(is->A,((PetscObject)A)->prefix);
656:   MatAppendOptionsPrefix(is->A,"is_");
657:   MatSetFromOptions(is->A);

659:   /* Create the local work vectors */
660:   MatCreateVecs(is->A,&is->x,&is->y);

662:   /* setup the global to local scatters */
663:   MatCreateVecs(A,&cglobal,&rglobal);
664:   ISCreateStride(PETSC_COMM_SELF,nr,0,1,&to);
665:   ISLocalToGlobalMappingApplyIS(rmapping,to,&from);
666:   VecScatterCreate(rglobal,from,is->y,to,&is->rctx);
667:   if (rmapping != cmapping) {
668:     ISDestroy(&to);
669:     ISDestroy(&from);
670:     ISCreateStride(PETSC_COMM_SELF,nc,0,1,&to);
671:     ISLocalToGlobalMappingApplyIS(cmapping,to,&from);
672:     VecScatterCreate(cglobal,from,is->x,to,&is->cctx);
673:   } else {
674:     PetscObjectReference((PetscObject)is->rctx);
675:     is->cctx = is->rctx;
676:   }
677:   VecDestroy(&rglobal);
678:   VecDestroy(&cglobal);
679:   ISDestroy(&to);
680:   ISDestroy(&from);
681:   return(0);
682: }

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

693: #if defined(PETSC_USE_DEBUG)
694:   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);
695: #endif
696:   ISG2LMapApply(mat->rmap->mapping,m,rows,rows_l);
697:   ISG2LMapApply(mat->cmap->mapping,n,cols,cols_l);
698:   MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);
699:   return(0);
700: }

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

710:   MatSetValues(is->A,m,rows,n,cols,values,addv);
711:   return(0);
712: }

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

722:   MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv);
723:   return(0);
724: }

728: PetscErrorCode MatZeroRows_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
729: {
730:   PetscInt       n_l = 0, *rows_l = NULL;

734:   if (x && b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support");
735:   if (n) {
736:     PetscMalloc1(n,&rows_l);
737:     ISGlobalToLocalMappingApply(A->rmap->mapping,IS_GTOLM_DROP,n,rows,&n_l,rows_l);
738:   }
739:   MatZeroRowsLocal(A,n_l,rows_l,diag,x,b);
740:   PetscFree(rows_l);
741:   return(0);
742: }

746: PetscErrorCode MatZeroRowsLocal_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
747: {
748:   Mat_IS         *is = (Mat_IS*)A->data;
750:   PetscInt       i;
751:   PetscScalar    *array;

754:   if (x && b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support");
755:   {
756:     /*
757:        Set up is->x as a "counting vector". This is in order to MatMult_IS
758:        work properly in the interface nodes.
759:     */
760:     Vec counter;
761:     MatCreateVecs(A,NULL,&counter);
762:     VecSet(counter,0.);
763:     VecSet(is->y,1.);
764:     VecScatterBegin(is->rctx,is->y,counter,ADD_VALUES,SCATTER_REVERSE);
765:     VecScatterEnd(is->rctx,is->y,counter,ADD_VALUES,SCATTER_REVERSE);
766:     VecScatterBegin(is->rctx,counter,is->y,INSERT_VALUES,SCATTER_FORWARD);
767:     VecScatterEnd(is->rctx,counter,is->y,INSERT_VALUES,SCATTER_FORWARD);
768:     VecDestroy(&counter);
769:   }
770:   if (!n) {
771:     is->pure_neumann = PETSC_TRUE;
772:   } else {
773:     is->pure_neumann = PETSC_FALSE;

775:     VecGetArray(is->y,&array);
776:     MatZeroRows(is->A,n,rows,diag,0,0);
777:     for (i=0; i<n; i++) {
778:       MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);
779:     }
780:     MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);
781:     MatAssemblyEnd(is->A,MAT_FINAL_ASSEMBLY);
782:     VecRestoreArray(is->y,&array);
783:   }
784:   return(0);
785: }

789: PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type)
790: {
791:   Mat_IS         *is = (Mat_IS*)A->data;

795:   MatAssemblyBegin(is->A,type);
796:   return(0);
797: }

801: PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type)
802: {
803:   Mat_IS         *is = (Mat_IS*)A->data;

807:   MatAssemblyEnd(is->A,type);
808:   return(0);
809: }

813: PetscErrorCode MatISGetLocalMat_IS(Mat mat,Mat *local)
814: {
815:   Mat_IS *is = (Mat_IS*)mat->data;

818:   *local = is->A;
819:   return(0);
820: }

824: /*@
825:     MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix.

827:   Input Parameter:
828: .  mat - the matrix

830:   Output Parameter:
831: .  local - the local matrix

833:   Level: advanced

835:   Notes:
836:     This can be called if you have precomputed the nonzero structure of the
837:   matrix and want to provide it to the inner matrix object to improve the performance
838:   of the MatSetValues() operation.

840: .seealso: MATIS
841: @*/
842: PetscErrorCode MatISGetLocalMat(Mat mat,Mat *local)
843: {

849:   PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat*),(mat,local));
850:   return(0);
851: }

855: PetscErrorCode MatISSetLocalMat_IS(Mat mat,Mat local)
856: {
857:   Mat_IS         *is = (Mat_IS*)mat->data;
858:   PetscInt       nrows,ncols,orows,ocols;

862:   if (is->A) {
863:     MatGetSize(is->A,&orows,&ocols);
864:     MatGetSize(local,&nrows,&ncols);
865:     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);
866:   }
867:   PetscObjectReference((PetscObject)local);
868:   MatDestroy(&is->A);
869:   is->A = local;
870:   return(0);
871: }

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

878:   Input Parameter:
879: .  mat - the matrix
880: .  local - the local matrix

882:   Output Parameter:

884:   Level: advanced

886:   Notes:
887:     This can be called if you have precomputed the local matrix and
888:   want to provide it to the matrix object MATIS.

890: .seealso: MATIS
891: @*/
892: PetscErrorCode MatISSetLocalMat(Mat mat,Mat local)
893: {

899:   PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local));
900:   return(0);
901: }

905: PetscErrorCode MatZeroEntries_IS(Mat A)
906: {
907:   Mat_IS         *a = (Mat_IS*)A->data;

911:   MatZeroEntries(a->A);
912:   return(0);
913: }

917: PetscErrorCode MatScale_IS(Mat A,PetscScalar a)
918: {
919:   Mat_IS         *is = (Mat_IS*)A->data;

923:   MatScale(is->A,a);
924:   return(0);
925: }

929: PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v)
930: {
931:   Mat_IS         *is = (Mat_IS*)A->data;

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

938:   /* scatter diagonal back into global vector */
939:   VecSet(v,0);
940:   VecScatterBegin(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);
941:   VecScatterEnd(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);
942:   return(0);
943: }

947: PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg)
948: {
949:   Mat_IS         *a = (Mat_IS*)A->data;

953:   MatSetOption(a->A,op,flg);
954:   return(0);
955: }

959: /*@
960:     MatCreateIS - Creates a "process" unassmembled matrix, it is assembled on each
961:        process but not across processes.

963:    Input Parameters:
964: +     comm    - MPI communicator that will share the matrix
965: .     bs      - block size of the matrix
966: .     m,n,M,N - local and/or global sizes of the left and right vector used in matrix vector products
967: .     rmap    - local to global map for rows
968: -     cmap    - local to global map for cols

970:    Output Parameter:
971: .    A - the resulting matrix

973:    Level: advanced

975:    Notes: See MATIS for more details
976:           m and n are NOT related to the size of the map, they are the size of the part of the vector owned
977:           by that process. The sizes of rmap and cmap define the size of the local matrices.
978:           If either rmap or cmap are NULL, than the matrix is assumed to be square

980: .seealso: MATIS, MatSetLocalToGlobalMapping()
981: @*/
982: PetscErrorCode  MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping rmap,ISLocalToGlobalMapping cmap,Mat *A)
983: {

987:   if (!rmap && !cmap) SETERRQ(comm,PETSC_ERR_USER,"You need to provide at least one of the mapping");
988:   MatCreate(comm,A);
989:   MatSetBlockSize(*A,bs);
990:   MatSetSizes(*A,m,n,M,N);
991:   MatSetType(*A,MATIS);
992:   MatSetUp(*A);
993:   if (rmap && cmap) {
994:     MatSetLocalToGlobalMapping(*A,rmap,cmap);
995:   } else if (!rmap) {
996:     MatSetLocalToGlobalMapping(*A,cmap,cmap);
997:   } else {
998:     MatSetLocalToGlobalMapping(*A,rmap,rmap);
999:   }
1000:   return(0);
1001: }

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

1009:    Operations Provided:
1010: +  MatMult()
1011: .  MatMultAdd()
1012: .  MatMultTranspose()
1013: .  MatMultTransposeAdd()
1014: .  MatZeroEntries()
1015: .  MatSetOption()
1016: .  MatZeroRows()
1017: .  MatZeroRowsLocal()
1018: .  MatSetValues()
1019: .  MatSetValuesLocal()
1020: .  MatScale()
1021: .  MatGetDiagonal()
1022: -  MatSetLocalToGlobalMapping()

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

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

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

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

1034:   Level: advanced

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

1038: M*/

1042: PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A)
1043: {
1045:   Mat_IS         *b;

1048:   PetscNewLog(A,&b);
1049:   A->data = (void*)b;

1051:   /* matrix ops */
1052:   PetscMemzero(A->ops,sizeof(struct _MatOps));
1053:   A->ops->mult                    = MatMult_IS;
1054:   A->ops->multadd                 = MatMultAdd_IS;
1055:   A->ops->multtranspose           = MatMultTranspose_IS;
1056:   A->ops->multtransposeadd        = MatMultTransposeAdd_IS;
1057:   A->ops->destroy                 = MatDestroy_IS;
1058:   A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS;
1059:   A->ops->setvalues               = MatSetValues_IS;
1060:   A->ops->setvalueslocal          = MatSetValuesLocal_IS;
1061:   A->ops->setvaluesblockedlocal   = MatSetValuesBlockedLocal_IS;
1062:   A->ops->zerorows                = MatZeroRows_IS;
1063:   A->ops->zerorowslocal           = MatZeroRowsLocal_IS;
1064:   A->ops->assemblybegin           = MatAssemblyBegin_IS;
1065:   A->ops->assemblyend             = MatAssemblyEnd_IS;
1066:   A->ops->view                    = MatView_IS;
1067:   A->ops->zeroentries             = MatZeroEntries_IS;
1068:   A->ops->scale                   = MatScale_IS;
1069:   A->ops->getdiagonal             = MatGetDiagonal_IS;
1070:   A->ops->setoption               = MatSetOption_IS;
1071:   A->ops->ishermitian             = MatIsHermitian_IS;
1072:   A->ops->issymmetric             = MatIsSymmetric_IS;
1073:   A->ops->duplicate               = MatDuplicate_IS;

1075:   /* special MATIS functions */
1076:   PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS);
1077:   PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS);
1078:   PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",MatISGetMPIXAIJ_IS);
1079:   PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",MatISSetPreallocation_IS);
1080:   PetscObjectChangeTypeName((PetscObject)A,MATIS);
1081:   return(0);
1082: }