Actual source code: matis.c

petsc-master 2015-07-03
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) {
 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_r,*global_indices_c;
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(A->rmap->mapping,&global_indices_r);
153:   if (A->rmap->mapping != A->cmap->mapping) {
154:     ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&global_indices_c);
155:   } else {
156:     global_indices_c = global_indices_r;
157:   }

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

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

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

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

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

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

309:   if (reuse == MAT_INITIAL_MATRIX) {
310:     MatType     new_mat_type;
311:     PetscBool   issbaij_red;

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

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

354:   if (issbaij) {
355:     MatConvert(matis->A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&local_mat);
356:   } else {
357:     PetscObjectReference((PetscObject)matis->A);
358:     local_mat = matis->A;
359:   }

361:   /* Set values */
362:   MatSetLocalToGlobalMapping(*M,mat->rmap->mapping,mat->cmap->mapping);
363:   if (isdense) { /* special case for dense local matrices */
364:     PetscInt i,*dummy_rows,*dummy_cols;

366:     if (local_rows != local_cols) {
367:       PetscMalloc2(local_rows,&dummy_rows,local_cols,&dummy_cols);
368:       for (i=0;i<local_rows;i++) dummy_rows[i] = i;
369:       for (i=0;i<local_cols;i++) dummy_cols[i] = i;
370:     } else {
371:       PetscMalloc1(local_rows,&dummy_rows);
372:       for (i=0;i<local_rows;i++) dummy_rows[i] = i;
373:       dummy_cols = dummy_rows;
374:     }
375:     MatSetOption(*M,MAT_ROW_ORIENTED,PETSC_FALSE);
376:     MatDenseGetArray(local_mat,&array);
377:     MatSetValuesLocal(*M,local_rows,dummy_rows,local_cols,dummy_cols,array,ADD_VALUES);
378:     MatDenseRestoreArray(local_mat,&array);
379:     if (dummy_rows != dummy_cols) {
380:       PetscFree2(dummy_rows,dummy_cols);
381:     } else {
382:       PetscFree(dummy_rows);
383:     }
384:   } else if (isseqaij) {
385:     PetscInt  i,nvtxs,*xadj,*adjncy;
386:     PetscBool done;

388:     MatGetRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);
389:     if (!done) {
390:       SETERRQ1(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ called in %s\n",__FUNCT__);
391:     }
392:     MatSeqAIJGetArray(local_mat,&array);
393:     for (i=0;i<nvtxs;i++) {
394:       MatSetValuesLocal(*M,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],array+xadj[i],ADD_VALUES);
395:     }
396:     MatRestoreRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);
397:     if (!done) {
398:       SETERRQ1(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ called in %s\n",__FUNCT__);
399:     }
400:     MatSeqAIJRestoreArray(local_mat,&array);
401:   } else { /* very basic values insertion for all other matrix types */
402:     PetscInt i;

404:     for (i=0;i<local_rows;i++) {
405:       PetscInt       j;
406:       const PetscInt *local_indices_cols;

408:       MatGetRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);
409:       MatSetValuesLocal(*M,1,&i,j,local_indices_cols,array,ADD_VALUES);
410:       MatRestoreRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);
411:     }
412:   }
413:   MatAssemblyBegin(*M,MAT_FINAL_ASSEMBLY);
414:   MatDestroy(&local_mat);
415:   MatAssemblyEnd(*M,MAT_FINAL_ASSEMBLY);
416:   if (isdense) {
417:     MatSetOption(*M,MAT_ROW_ORIENTED,PETSC_TRUE);
418:   }
419:   return(0);
420: }

424: /*@
425:     MatISGetMPIXAIJ - Converts MATIS matrix into a parallel AIJ format

427:   Input Parameter:
428: .  mat - the matrix (should be of type MATIS)
429: .  reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX

431:   Output Parameter:
432: .  newmat - the matrix in AIJ format

434:   Level: developer

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

438: .seealso: MATIS
439: @*/
440: PetscErrorCode MatISGetMPIXAIJ(Mat mat, MatReuse reuse, Mat *newmat)
441: {

448:   if (reuse != MAT_INITIAL_MATRIX) {
451:     if (mat == *newmat) {
452:       SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse the same matrix");
453:     }
454:   }
455:   PetscUseMethod(mat,"MatISGetMPIXAIJ_C",(Mat,MatReuse,Mat*),(mat,reuse,newmat));
456:   return(0);
457: }

461: PetscErrorCode MatDuplicate_IS(Mat mat,MatDuplicateOption op,Mat *newmat)
462: {
464:   Mat_IS         *matis = (Mat_IS*)(mat->data);
465:   PetscInt       bs,m,n,M,N;
466:   Mat            B,localmat;

469:   MatGetBlockSize(mat,&bs);
470:   MatGetSize(mat,&M,&N);
471:   MatGetLocalSize(mat,&m,&n);
472:   MatCreateIS(PetscObjectComm((PetscObject)mat),bs,m,n,M,N,mat->rmap->mapping,mat->cmap->mapping,&B);
473:   MatDuplicate(matis->A,op,&localmat);
474:   MatISSetLocalMat(B,localmat);
475:   MatDestroy(&localmat);
476:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
477:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
478:   *newmat = B;
479:   return(0);
480: }

484: PetscErrorCode MatIsHermitian_IS(Mat A,PetscReal tol,PetscBool  *flg)
485: {
487:   Mat_IS         *matis = (Mat_IS*)A->data;
488:   PetscBool      local_sym;

491:   MatIsHermitian(matis->A,tol,&local_sym);
492:   MPI_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
493:   return(0);
494: }

498: PetscErrorCode MatIsSymmetric_IS(Mat A,PetscReal tol,PetscBool  *flg)
499: {
501:   Mat_IS         *matis = (Mat_IS*)A->data;
502:   PetscBool      local_sym;

505:   MatIsSymmetric(matis->A,tol,&local_sym);
506:   MPI_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
507:   return(0);
508: }

512: PetscErrorCode MatDestroy_IS(Mat A)
513: {
515:   Mat_IS         *b = (Mat_IS*)A->data;

518:   MatDestroy(&b->A);
519:   VecScatterDestroy(&b->cctx);
520:   VecScatterDestroy(&b->rctx);
521:   VecDestroy(&b->x);
522:   VecDestroy(&b->y);
523:   PetscSFDestroy(&b->sf);
524:   PetscFree2(b->sf_rootdata,b->sf_leafdata);
525:   PetscFree(A->data);
526:   PetscObjectChangeTypeName((PetscObject)A,0);
527:   PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL);
528:   PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",NULL);
529:   PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",NULL);
530:   PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",NULL);
531:   return(0);
532: }

536: PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y)
537: {
539:   Mat_IS         *is  = (Mat_IS*)A->data;
540:   PetscScalar    zero = 0.0;

543:   /*  scatter the global vector x into the local work vector */
544:   VecScatterBegin(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
545:   VecScatterEnd(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);

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

550:   /* scatter product back into global memory */
551:   VecSet(y,zero);
552:   VecScatterBegin(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
553:   VecScatterEnd(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
554:   return(0);
555: }

559: PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
560: {
561:   Vec            temp_vec;

565:   if (v3 != v2) {
566:     MatMult(A,v1,v3);
567:     VecAXPY(v3,1.0,v2);
568:   } else {
569:     VecDuplicate(v2,&temp_vec);
570:     MatMult(A,v1,temp_vec);
571:     VecAXPY(temp_vec,1.0,v2);
572:     VecCopy(temp_vec,v3);
573:     VecDestroy(&temp_vec);
574:   }
575:   return(0);
576: }

580: PetscErrorCode MatMultTranspose_IS(Mat A,Vec y,Vec x)
581: {
582:   Mat_IS         *is = (Mat_IS*)A->data;

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

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

593:   /* scatter product back into global vector */
594:   VecSet(x,0);
595:   VecScatterBegin(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);
596:   VecScatterEnd(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);
597:   return(0);
598: }

602: PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
603: {
604:   Vec            temp_vec;

608:   if (v3 != v2) {
609:     MatMultTranspose(A,v1,v3);
610:     VecAXPY(v3,1.0,v2);
611:   } else {
612:     VecDuplicate(v2,&temp_vec);
613:     MatMultTranspose(A,v1,temp_vec);
614:     VecAXPY(temp_vec,1.0,v2);
615:     VecCopy(temp_vec,v3);
616:     VecDestroy(&temp_vec);
617:   }
618:   return(0);
619: }

623: PetscErrorCode MatView_IS(Mat A,PetscViewer viewer)
624: {
625:   Mat_IS         *a = (Mat_IS*)A->data;
627:   PetscViewer    sviewer;

630:   PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
631:   PetscViewerGetSingleton(viewer,&sviewer);
632:   MatView(a->A,sviewer);
633:   PetscViewerRestoreSingleton(viewer,&sviewer);
634:   return(0);
635: }

639: PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping)
640: {
642:   PetscInt       nr,rbs,nc,cbs;
643:   Mat_IS         *is = (Mat_IS*)A->data;
644:   IS             from,to;
645:   Vec            cglobal,rglobal;

650:   /* Destroy any previous data */
651:   VecDestroy(&is->x);
652:   VecDestroy(&is->y);
653:   VecScatterDestroy(&is->rctx);
654:   VecScatterDestroy(&is->cctx);
655:   MatDestroy(&is->A);
656:   PetscSFDestroy(&is->sf);
657:   PetscFree2(is->sf_rootdata,is->sf_leafdata);

659:   /* Setup Layout and set local to global maps */
660:   PetscLayoutSetUp(A->rmap);
661:   PetscLayoutSetUp(A->cmap);
662:   PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping);
663:   PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping);

665:   /* Create the local matrix A */
666:   ISLocalToGlobalMappingGetSize(rmapping,&nr);
667:   ISLocalToGlobalMappingGetBlockSize(rmapping,&rbs);
668:   ISLocalToGlobalMappingGetSize(cmapping,&nc);
669:   ISLocalToGlobalMappingGetBlockSize(cmapping,&cbs);
670:   MatCreate(PETSC_COMM_SELF,&is->A);
671:   MatSetType(is->A,MATAIJ);
672:   MatSetSizes(is->A,nr,nc,nr,nc);
673:   MatSetBlockSizes(is->A,rbs,cbs);
674:   MatSetOptionsPrefix(is->A,((PetscObject)A)->prefix);
675:   MatAppendOptionsPrefix(is->A,"is_");
676:   MatSetFromOptions(is->A);

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

681:   /* setup the global to local scatters */
682:   MatCreateVecs(A,&cglobal,&rglobal);
683:   ISCreateStride(PETSC_COMM_SELF,nr,0,1,&to);
684:   ISLocalToGlobalMappingApplyIS(rmapping,to,&from);
685:   VecScatterCreate(rglobal,from,is->y,to,&is->rctx);
686:   if (rmapping != cmapping) {
687:     ISDestroy(&to);
688:     ISDestroy(&from);
689:     ISCreateStride(PETSC_COMM_SELF,nc,0,1,&to);
690:     ISLocalToGlobalMappingApplyIS(cmapping,to,&from);
691:     VecScatterCreate(cglobal,from,is->x,to,&is->cctx);
692:   } else {
693:     PetscObjectReference((PetscObject)is->rctx);
694:     is->cctx = is->rctx;
695:   }
696:   VecDestroy(&rglobal);
697:   VecDestroy(&cglobal);
698:   ISDestroy(&to);
699:   ISDestroy(&from);
700:   return(0);
701: }

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

712: #if defined(PETSC_USE_DEBUG)
713:   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);
714: #endif
715:   ISG2LMapApply(mat->rmap->mapping,m,rows,rows_l);
716:   ISG2LMapApply(mat->cmap->mapping,n,cols,cols_l);
717:   MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);
718:   return(0);
719: }

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

729:   MatSetValues(is->A,m,rows,n,cols,values,addv);
730:   return(0);
731: }

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

741:   MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv);
742:   return(0);
743: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

852:   Level: advanced

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

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

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

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

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

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

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

901:   Output Parameter:

903:   Level: advanced

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

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

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

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

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

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

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

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

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

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

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

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

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

982:    Input Parameters:
983: +     comm    - MPI communicator that will share the matrix
984: .     bs      - block size of the matrix
985: .     m,n,M,N - local and/or global sizes of the left and right vector used in matrix vector products
986: .     rmap    - local to global map for rows
987: -     cmap    - local to global map for cols

989:    Output Parameter:
990: .    A - the resulting matrix

992:    Level: advanced

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

999: .seealso: MATIS, MatSetLocalToGlobalMapping()
1000: @*/
1001: PetscErrorCode  MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping rmap,ISLocalToGlobalMapping cmap,Mat *A)
1002: {

1006:   if (!rmap && !cmap) SETERRQ(comm,PETSC_ERR_USER,"You need to provide at least one of the mapping");
1007:   MatCreate(comm,A);
1008:   MatSetBlockSize(*A,bs);
1009:   MatSetSizes(*A,m,n,M,N);
1010:   MatSetType(*A,MATIS);
1011:   MatSetUp(*A);
1012:   if (rmap && cmap) {
1013:     MatSetLocalToGlobalMapping(*A,rmap,cmap);
1014:   } else if (!rmap) {
1015:     MatSetLocalToGlobalMapping(*A,cmap,cmap);
1016:   } else {
1017:     MatSetLocalToGlobalMapping(*A,rmap,rmap);
1018:   }
1019:   return(0);
1020: }

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

1028:    Operations Provided:
1029: +  MatMult()
1030: .  MatMultAdd()
1031: .  MatMultTranspose()
1032: .  MatMultTransposeAdd()
1033: .  MatZeroEntries()
1034: .  MatSetOption()
1035: .  MatZeroRows()
1036: .  MatZeroRowsLocal()
1037: .  MatSetValues()
1038: .  MatSetValuesLocal()
1039: .  MatScale()
1040: .  MatGetDiagonal()
1041: -  MatSetLocalToGlobalMapping()

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

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

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

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

1053:   Level: advanced

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

1057: M*/

1061: PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A)
1062: {
1064:   Mat_IS         *b;

1067:   PetscNewLog(A,&b);
1068:   A->data = (void*)b;

1070:   /* matrix ops */
1071:   PetscMemzero(A->ops,sizeof(struct _MatOps));
1072:   A->ops->mult                    = MatMult_IS;
1073:   A->ops->multadd                 = MatMultAdd_IS;
1074:   A->ops->multtranspose           = MatMultTranspose_IS;
1075:   A->ops->multtransposeadd        = MatMultTransposeAdd_IS;
1076:   A->ops->destroy                 = MatDestroy_IS;
1077:   A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS;
1078:   A->ops->setvalues               = MatSetValues_IS;
1079:   A->ops->setvalueslocal          = MatSetValuesLocal_IS;
1080:   A->ops->setvaluesblockedlocal   = MatSetValuesBlockedLocal_IS;
1081:   A->ops->zerorows                = MatZeroRows_IS;
1082:   A->ops->zerorowslocal           = MatZeroRowsLocal_IS;
1083:   A->ops->assemblybegin           = MatAssemblyBegin_IS;
1084:   A->ops->assemblyend             = MatAssemblyEnd_IS;
1085:   A->ops->view                    = MatView_IS;
1086:   A->ops->zeroentries             = MatZeroEntries_IS;
1087:   A->ops->scale                   = MatScale_IS;
1088:   A->ops->getdiagonal             = MatGetDiagonal_IS;
1089:   A->ops->setoption               = MatSetOption_IS;
1090:   A->ops->ishermitian             = MatIsHermitian_IS;
1091:   A->ops->issymmetric             = MatIsSymmetric_IS;
1092:   A->ops->duplicate               = MatDuplicate_IS;

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: }