Actual source code: matis.c

petsc-master 2015-09-02
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:   PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
631:   MatView(a->A,sviewer);
632:   PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
633:   return(0);
634: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

772:   if (x && b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support");
773:   {
774:     /*
775:        Set up is->x as a "counting vector". This is in order to MatMult_IS
776:        work properly in the interface nodes.
777:     */
778:     Vec counter;
779:     MatCreateVecs(A,NULL,&counter);
780:     VecSet(counter,0.);
781:     VecSet(is->y,1.);
782:     VecScatterBegin(is->rctx,is->y,counter,ADD_VALUES,SCATTER_REVERSE);
783:     VecScatterEnd(is->rctx,is->y,counter,ADD_VALUES,SCATTER_REVERSE);
784:     VecScatterBegin(is->rctx,counter,is->y,INSERT_VALUES,SCATTER_FORWARD);
785:     VecScatterEnd(is->rctx,counter,is->y,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->y,&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->y,&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->y);

956:   /* scatter diagonal back into global vector */
957:   VecSet(v,0);
958:   VecScatterBegin(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);
959:   VecScatterEnd(is->rctx,is->y,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      - 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: .     rmap    - local to global map for rows
986: -     cmap    - local to global map for cols

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

991:    Level: advanced

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

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

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

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

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

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

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

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

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

1052:   Level: advanced

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

1056: M*/

1060: PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A)
1061: {
1063:   Mat_IS         *b;

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

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

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