Actual source code: matis.c

petsc-dev 2014-08-28
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*/

 19: PetscErrorCode MatISGetMPIXAIJ_IS(Mat mat, MatReuse reuse, Mat *M)
 20: {
 21:   Mat_IS                 *matis = (Mat_IS*)(mat->data);
 22:   /* info on mat */
 23:   /* ISLocalToGlobalMapping rmapping,cmapping; */
 24:   PetscInt               bs,rows,cols;
 25:   PetscInt               lrows,lcols;
 26:   PetscInt               local_rows,local_cols;
 27:   PetscBool              isdense,issbaij,issbaij_red;
 28:   /* values insertion */
 29:   PetscScalar            *array;
 30:   PetscInt               *local_indices,*global_indices;
 31:   /* work */
 32:   PetscInt               i,j,index_row;
 33:   PetscErrorCode         ierr;

 36:   /* MISSING CHECKS
 37:     - rectangular case not covered (it is not allowed by MATIS)
 38:   */
 39:   /* get info from mat */
 40:   /* MatGetLocalToGlobalMapping(mat,&rmapping,&cmapping); */
 41:   MatGetSize(mat,&rows,&cols);
 42:   MatGetBlockSize(mat,&bs);
 43:   MatGetSize(matis->A,&local_rows,&local_cols);
 44:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isdense);
 45:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);

 47:   /* work */
 48:   PetscMalloc1(local_rows,&local_indices);
 49:   for (i=0;i<local_rows;i++) local_indices[i]=i;
 50:   /* map indices of local mat to global values */
 51:   PetscMalloc(PetscMax(local_cols,local_rows)*sizeof(*global_indices),&global_indices);
 52:   /* ISLocalToGlobalMappingApply(rmapping,local_rows,local_indices,global_indices); */
 53:   ISLocalToGlobalMappingApply(matis->mapping,local_rows,local_indices,global_indices);

 55:   if (issbaij) {
 56:     MatGetRowUpperTriangular(matis->A);
 57:   }

 59:   if (reuse == MAT_INITIAL_MATRIX) {
 60:     Mat         new_mat;
 61:     MatType     new_mat_type;
 62:     Vec         vec_dnz,vec_onz;
 63:     PetscScalar *my_dnz,*my_onz;
 64:     PetscInt    *dnz,*onz,*mat_ranges,*row_ownership;
 65:     PetscInt    index_col,owner;
 66:     PetscMPIInt nsubdomains;

 68:     /* determining new matrix type */
 69:     MPI_Allreduce(&issbaij,&issbaij_red,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)mat));
 70:     if (issbaij_red) {
 71:       new_mat_type = MATSBAIJ;
 72:     } else {
 73:       if (bs>1) {
 74:         new_mat_type = MATBAIJ;
 75:       } else {
 76:         new_mat_type = MATAIJ;
 77:       }
 78:     }

 80:     MPI_Comm_size(PetscObjectComm((PetscObject)mat),&nsubdomains);
 81:     MatCreate(PetscObjectComm((PetscObject)mat),&new_mat);
 82:     MatSetSizes(new_mat,PETSC_DECIDE,PETSC_DECIDE,rows,cols);
 83:     MatSetBlockSize(new_mat,bs);
 84:     MatSetType(new_mat,new_mat_type);
 85:     MatSetUp(new_mat);
 86:     MatGetLocalSize(new_mat,&lrows,&lcols);

 88:     /*
 89:       preallocation
 90:     */

 92:     MatPreallocateInitialize(PetscObjectComm((PetscObject)new_mat),lrows,lcols,dnz,onz);
 93:     /*
 94:        Some vectors are needed to sum up properly on shared interface dofs.
 95:        Preallocation macros cannot do the job.
 96:        Note that preallocation is not exact, since it overestimates nonzeros
 97:     */
 98:     MatGetVecs(new_mat,NULL,&vec_dnz);
 99:     /* VecSetLocalToGlobalMapping(vec_dnz,rmapping); */
100:     VecSetLocalToGlobalMapping(vec_dnz,matis->mapping);
101:     VecDuplicate(vec_dnz,&vec_onz);
102:     /* All processes need to compute entire row ownership */
103:     PetscMalloc1(rows,&row_ownership);
104:     MatGetOwnershipRanges(new_mat,(const PetscInt**)&mat_ranges);
105:     for (i=0;i<nsubdomains;i++) {
106:       for (j=mat_ranges[i];j<mat_ranges[i+1];j++) {
107:         row_ownership[j]=i;
108:       }
109:     }

111:     /*
112:        my_dnz and my_onz contains exact contribution to preallocation from each local mat
113:        then, they will be summed up properly. This way, preallocation is always sufficient
114:     */
115:     PetscMalloc1(local_rows,&my_dnz);
116:     PetscMalloc1(local_rows,&my_onz);
117:     PetscMemzero(my_dnz,local_rows*sizeof(*my_dnz));
118:     PetscMemzero(my_onz,local_rows*sizeof(*my_onz));
119:     /* preallocation as a MATAIJ */
120:     if (isdense) { /* special case for dense local matrices */
121:       for (i=0;i<local_rows;i++) {
122:         index_row = global_indices[i];
123:         for (j=i;j<local_rows;j++) {
124:           owner = row_ownership[index_row];
125:           index_col = global_indices[j];
126:           if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */
127:             my_dnz[i] += 1.0;
128:           } else { /* offdiag block */
129:             my_onz[i] += 1.0;
130:           }
131:           /* same as before, interchanging rows and cols */
132:           if (i != j) {
133:             owner = row_ownership[index_col];
134:             if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) {
135:               my_dnz[j] += 1.0;
136:             } else {
137:               my_onz[j] += 1.0;
138:             }
139:           }
140:         }
141:       }
142:     } else {
143:       for (i=0;i<local_rows;i++) {
144:         PetscInt ncols;
145:         const PetscInt *cols;
146:         index_row = global_indices[i];
147:         MatGetRow(matis->A,i,&ncols,&cols,NULL);
148:         for (j=0;j<ncols;j++) {
149:           owner = row_ownership[index_row];
150:           index_col = global_indices[cols[j]];
151:           if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */
152:             my_dnz[i] += 1.0;
153:           } else { /* offdiag block */
154:             my_onz[i] += 1.0;
155:           }
156:           /* same as before, interchanging rows and cols */
157:           if (issbaij) {
158:             owner = row_ownership[index_col];
159:             if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) {
160:               my_dnz[j] += 1.0;
161:             } else {
162:               my_onz[j] += 1.0;
163:             }
164:           }
165:         }
166:         MatRestoreRow(matis->A,i,&ncols,&cols,NULL);
167:       }
168:     }
169:     VecSet(vec_dnz,0.0);
170:     VecSet(vec_onz,0.0);
171:     if (local_rows) { /* multilevel guard */
172:       VecSetValuesLocal(vec_dnz,local_rows,local_indices,my_dnz,ADD_VALUES);
173:       VecSetValuesLocal(vec_onz,local_rows,local_indices,my_onz,ADD_VALUES);
174:     }
175:     VecAssemblyBegin(vec_dnz);
176:     VecAssemblyBegin(vec_onz);
177:     VecAssemblyEnd(vec_dnz);
178:     VecAssemblyEnd(vec_onz);
179:     PetscFree(my_dnz);
180:     PetscFree(my_onz);
181:     PetscFree(row_ownership);

183:     /* set computed preallocation in dnz and onz */
184:     VecGetArray(vec_dnz,&array);
185:     for (i=0; i<lrows; i++) dnz[i] = (PetscInt)PetscRealPart(array[i]);
186:     VecRestoreArray(vec_dnz,&array);
187:     VecGetArray(vec_onz,&array);
188:     for (i=0;i<lrows;i++) onz[i] = (PetscInt)PetscRealPart(array[i]);
189:     VecRestoreArray(vec_onz,&array);
190:     VecDestroy(&vec_dnz);
191:     VecDestroy(&vec_onz);

193:     /* Resize preallocation if overestimated */
194:     for (i=0;i<lrows;i++) {
195:       dnz[i] = PetscMin(dnz[i],lcols);
196:       onz[i] = PetscMin(onz[i],cols-lcols);
197:     }
198:     /* set preallocation */
199:     MatMPIAIJSetPreallocation(new_mat,0,dnz,0,onz);
200:     for (i=0;i<lrows/bs;i++) {
201:       dnz[i] = dnz[i*bs]/bs;
202:       onz[i] = onz[i*bs]/bs;
203:     }
204:     MatMPIBAIJSetPreallocation(new_mat,bs,0,dnz,0,onz);
205:     for (i=0;i<lrows/bs;i++) {
206:       dnz[i] = dnz[i]-i;
207:     }
208:     MatMPISBAIJSetPreallocation(new_mat,bs,0,dnz,0,onz);
209:     MatPreallocateFinalize(dnz,onz);
210:     *M = new_mat;
211:   } else {
212:     PetscInt mbs,mrows,mcols;
213:     /* some checks */
214:     MatGetBlockSize(*M,&mbs);
215:     MatGetSize(*M,&mrows,&mcols);
216:     if (mrows != rows) {
217:       SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of rows (%d != %d)",rows,mrows);
218:     }
219:     if (mrows != rows) {
220:       SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of cols (%d != %d)",cols,mcols);
221:     }
222:     if (mbs != bs) {
223:       SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong block size (%d != %d)",bs,mbs);
224:     }
225:     MatZeroEntries(*M);
226:   }
227:   /* set local to global mappings */
228:   /* MatSetLocalToGlobalMapping(*M,rmapping,cmapping); */
229:   /* Set values */
230:   if (isdense) { /* special case for dense local matrices */
231:     MatSetOption(*M,MAT_ROW_ORIENTED,PETSC_FALSE);
232:     MatDenseGetArray(matis->A,&array);
233:     MatSetValues(*M,local_rows,global_indices,local_cols,global_indices,array,ADD_VALUES);
234:     MatDenseRestoreArray(matis->A,&array);
235:     PetscFree(local_indices);
236:     PetscFree(global_indices);
237:   } else { /* very basic values insertion for all other matrix types */
238:     PetscFree(local_indices);
239:     for (i=0;i<local_rows;i++) {
240:       MatGetRow(matis->A,i,&j,(const PetscInt**)&local_indices,(const PetscScalar**)&array);
241:       /* MatSetValuesLocal(*M,1,&i,j,local_indices,array,ADD_VALUES); */
242:       ISLocalToGlobalMappingApply(matis->mapping,j,local_indices,global_indices);
243:       ISLocalToGlobalMappingApply(matis->mapping,1,&i,&index_row);
244:       MatSetValues(*M,1,&index_row,j,global_indices,array,ADD_VALUES);
245:       MatRestoreRow(matis->A,i,&j,(const PetscInt**)&local_indices,(const PetscScalar**)&array);
246:     }
247:     PetscFree(global_indices);
248:   }
249:   MatAssemblyBegin(*M,MAT_FINAL_ASSEMBLY);
250:   MatAssemblyEnd(*M,MAT_FINAL_ASSEMBLY);
251:   if (isdense) {
252:     MatSetOption(*M,MAT_ROW_ORIENTED,PETSC_TRUE);
253:   }
254:   if (issbaij) {
255:     MatRestoreRowUpperTriangular(matis->A);
256:   }
257:   return(0);
258: }

262: /*@
263:     MatISGetMPIXAIJ - Converts MATIS matrix into a parallel AIJ format

265:   Input Parameter:
266: .  mat - the matrix (should be of type MATIS)
267: .  reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX

269:   Output Parameter:
270: .  newmat - the matrix in AIJ format

272:   Level: developer

274:   Notes:

276: .seealso: MATIS
277: @*/
278: PetscErrorCode MatISGetMPIXAIJ(Mat mat, MatReuse reuse, Mat *newmat)
279: {

286:   if (reuse != MAT_INITIAL_MATRIX) {
289:   }
290:   PetscUseMethod(mat,"MatISGetMPIXAIJ_C",(Mat,MatReuse,Mat*),(mat,reuse,newmat));
291:   return(0);
292: }

296: PetscErrorCode MatDuplicate_IS(Mat mat,MatDuplicateOption op,Mat *newmat)
297: {
299:   Mat_IS         *matis = (Mat_IS*)(mat->data);
300:   PetscInt       bs,m,n,M,N;
301:   Mat            B,localmat;

304:   MatGetBlockSize(mat,&bs);
305:   MatGetSize(mat,&M,&N);
306:   MatGetLocalSize(mat,&m,&n);
307:   MatCreateIS(PetscObjectComm((PetscObject)mat),bs,m,n,M,N,matis->mapping,&B);
308:   MatDuplicate(matis->A,op,&localmat);
309:   MatISSetLocalMat(B,localmat);
310:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
311:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
312:   *newmat = B;
313:   return(0);
314: }

318: PetscErrorCode MatIsHermitian_IS(Mat A,PetscReal tol,PetscBool  *flg)
319: {
321:   Mat_IS         *matis = (Mat_IS*)A->data;
322:   PetscBool      local_sym;

325:   MatIsHermitian(matis->A,tol,&local_sym);
326:   MPI_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
327:   return(0);
328: }

332: PetscErrorCode MatIsSymmetric_IS(Mat A,PetscReal tol,PetscBool  *flg)
333: {
335:   Mat_IS         *matis = (Mat_IS*)A->data;
336:   PetscBool      local_sym;

339:   MatIsSymmetric(matis->A,tol,&local_sym);
340:   MPI_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
341:   return(0);
342: }

346: PetscErrorCode MatDestroy_IS(Mat A)
347: {
349:   Mat_IS         *b = (Mat_IS*)A->data;

352:   MatDestroy(&b->A);
353:   VecScatterDestroy(&b->ctx);
354:   VecDestroy(&b->x);
355:   VecDestroy(&b->y);
356:   ISLocalToGlobalMappingDestroy(&b->mapping);
357:   PetscFree(A->data);
358:   PetscObjectChangeTypeName((PetscObject)A,0);
359:   PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL);
360:   PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",NULL);
361:   PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",NULL);
362:   return(0);
363: }

367: PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y)
368: {
370:   Mat_IS         *is  = (Mat_IS*)A->data;
371:   PetscScalar    zero = 0.0;

374:   /*  scatter the global vector x into the local work vector */
375:   VecScatterBegin(is->ctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
376:   VecScatterEnd(is->ctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);

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

381:   /* scatter product back into global memory */
382:   VecSet(y,zero);
383:   VecScatterBegin(is->ctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
384:   VecScatterEnd(is->ctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
385:   return(0);
386: }

390: PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
391: {
392:   Vec            temp_vec;

396:   if (v3 != v2) {
397:     MatMult(A,v1,v3);
398:     VecAXPY(v3,1.0,v2);
399:   } else {
400:     VecDuplicate(v2,&temp_vec);
401:     MatMult(A,v1,temp_vec);
402:     VecAXPY(temp_vec,1.0,v2);
403:     VecCopy(temp_vec,v3);
404:     VecDestroy(&temp_vec);
405:   }
406:   return(0);
407: }

411: PetscErrorCode MatMultTranspose_IS(Mat A,Vec x,Vec y)
412: {
413:   Mat_IS         *is = (Mat_IS*)A->data;

417:   /*  scatter the global vector x into the local work vector */
418:   VecScatterBegin(is->ctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
419:   VecScatterEnd(is->ctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);

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

424:   /* scatter product back into global vector */
425:   VecSet(y,0);
426:   VecScatterBegin(is->ctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
427:   VecScatterEnd(is->ctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
428:   return(0);
429: }

433: PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
434: {
435:   Vec            temp_vec;

439:   if (v3 != v2) {
440:     MatMultTranspose(A,v1,v3);
441:     VecAXPY(v3,1.0,v2);
442:   } else {
443:     VecDuplicate(v2,&temp_vec);
444:     MatMultTranspose(A,v1,temp_vec);
445:     VecAXPY(temp_vec,1.0,v2);
446:     VecCopy(temp_vec,v3);
447:     VecDestroy(&temp_vec);
448:   }
449:   return(0);
450: }

454: PetscErrorCode MatView_IS(Mat A,PetscViewer viewer)
455: {
456:   Mat_IS         *a = (Mat_IS*)A->data;
458:   PetscViewer    sviewer;

461:   PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
462:   PetscViewerGetSingleton(viewer,&sviewer);
463:   MatView(a->A,sviewer);
464:   PetscViewerRestoreSingleton(viewer,&sviewer);
465:   return(0);
466: }

470: PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping)
471: {
473:   PetscInt       n,bs;
474:   Mat_IS         *is = (Mat_IS*)A->data;
475:   IS             from,to;
476:   Vec            global;

480:   if (rmapping != cmapping) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"MATIS requires the row and column mappings to be identical");
481:   if (is->mapping) { /* Currenly destroys the objects that will be created by this routine. Is there anything else that should be checked? */
482:     ISLocalToGlobalMappingDestroy(&is->mapping);
483:     VecDestroy(&is->x);
484:     VecDestroy(&is->y);
485:     VecScatterDestroy(&is->ctx);
486:     MatDestroy(&is->A);
487:   }
488:   PetscObjectReference((PetscObject)rmapping);
489:   ISLocalToGlobalMappingDestroy(&is->mapping);
490:   is->mapping = rmapping;
491: /*
492:   PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping);
493:   PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping);
494: */

496:   /* Create the local matrix A */
497:   ISLocalToGlobalMappingGetSize(rmapping,&n);
498:   MatGetBlockSize(A,&bs);
499:   MatCreate(PETSC_COMM_SELF,&is->A);
500:   MatSetSizes(is->A,n,n,n,n);
501:   MatSetBlockSize(is->A,bs);
502:   MatSetOptionsPrefix(is->A,((PetscObject)A)->prefix);
503:   MatAppendOptionsPrefix(is->A,"is_");
504:   MatSetFromOptions(is->A);

506:   /* Create the local work vectors */
507:   VecCreate(PETSC_COMM_SELF,&is->x);
508:   VecSetBlockSize(is->x,bs);
509:   VecSetSizes(is->x,n,n);
510:   VecSetOptionsPrefix(is->x,((PetscObject)A)->prefix);
511:   VecAppendOptionsPrefix(is->x,"is_");
512:   VecSetFromOptions(is->x);
513:   VecDuplicate(is->x,&is->y);

515:   /* setup the global to local scatter */
516:   ISCreateStride(PETSC_COMM_SELF,n,0,1,&to);
517:   ISLocalToGlobalMappingApplyIS(rmapping,to,&from);
518:   MatGetVecs(A,&global,NULL);
519:   VecScatterCreate(global,from,is->x,to,&is->ctx);
520:   VecDestroy(&global);
521:   ISDestroy(&to);
522:   ISDestroy(&from);
523:   return(0);
524: }

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

535: #if defined(PETSC_USE_DEBUG)
536:   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);
537: #endif
538:   ISG2LMapApply(is->mapping,m,rows,rows_l);
539:   ISG2LMapApply(is->mapping,n,cols,cols_l);
540:   MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);
541:   return(0);
542: }

544: #undef ISG2LMapSetUp
545: #undef ISG2LMapApply

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

555:   MatSetValues(is->A,m,rows,n,cols,values,addv);
556:   return(0);
557: }

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

567:   MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv);
568:   return(0);
569: }

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

580:   if (x && b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support");
581:   if (n) {
582:     PetscMalloc1(n,&rows_l);
583:     ISGlobalToLocalMappingApply(is->mapping,IS_GTOLM_DROP,n,rows,&n_l,rows_l);
584:   }
585:   MatZeroRowsLocal(A,n_l,rows_l,diag,x,b);
586:   PetscFree(rows_l);
587:   return(0);
588: }

592: PetscErrorCode MatZeroRowsLocal_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
593: {
594:   Mat_IS         *is = (Mat_IS*)A->data;
596:   PetscInt       i;
597:   PetscScalar    *array;

600:   if (x && b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support");
601:   {
602:     /*
603:        Set up is->x as a "counting vector". This is in order to MatMult_IS
604:        work properly in the interface nodes.
605:     */
606:     Vec         counter;
607:     PetscScalar one=1.0, zero=0.0;
608:     MatGetVecs(A,&counter,NULL);
609:     VecSet(counter,zero);
610:     VecSet(is->x,one);
611:     VecScatterBegin(is->ctx,is->x,counter,ADD_VALUES,SCATTER_REVERSE);
612:     VecScatterEnd  (is->ctx,is->x,counter,ADD_VALUES,SCATTER_REVERSE);
613:     VecScatterBegin(is->ctx,counter,is->x,INSERT_VALUES,SCATTER_FORWARD);
614:     VecScatterEnd  (is->ctx,counter,is->x,INSERT_VALUES,SCATTER_FORWARD);
615:     VecDestroy(&counter);
616:   }
617:   if (!n) {
618:     is->pure_neumann = PETSC_TRUE;
619:   } else {
620:     is->pure_neumann = PETSC_FALSE;

622:     VecGetArray(is->x,&array);
623:     MatZeroRows(is->A,n,rows,diag,0,0);
624:     for (i=0; i<n; i++) {
625:       MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);
626:     }
627:     MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);
628:     MatAssemblyEnd  (is->A,MAT_FINAL_ASSEMBLY);
629:     VecRestoreArray(is->x,&array);
630:   }
631:   return(0);
632: }

636: PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type)
637: {
638:   Mat_IS         *is = (Mat_IS*)A->data;

642:   MatAssemblyBegin(is->A,type);
643:   return(0);
644: }

648: PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type)
649: {
650:   Mat_IS         *is = (Mat_IS*)A->data;

654:   MatAssemblyEnd(is->A,type);
655:   return(0);
656: }

660: PetscErrorCode MatISGetLocalMat_IS(Mat mat,Mat *local)
661: {
662:   Mat_IS *is = (Mat_IS*)mat->data;

665:   *local = is->A;
666:   return(0);
667: }

671: /*@
672:     MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix.

674:   Input Parameter:
675: .  mat - the matrix

677:   Output Parameter:
678: .  local - the local matrix usually MATSEQAIJ

680:   Level: advanced

682:   Notes:
683:     This can be called if you have precomputed the nonzero structure of the
684:   matrix and want to provide it to the inner matrix object to improve the performance
685:   of the MatSetValues() operation.

687: .seealso: MATIS
688: @*/
689: PetscErrorCode MatISGetLocalMat(Mat mat,Mat *local)
690: {

696:   PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat*),(mat,local));
697:   return(0);
698: }

702: PetscErrorCode MatISSetLocalMat_IS(Mat mat,Mat local)
703: {
704:   Mat_IS         *is = (Mat_IS*)mat->data;
705:   PetscInt       nrows,ncols,orows,ocols;

709:   if (is->A) {
710:     MatGetSize(is->A,&orows,&ocols);
711:     MatGetSize(local,&nrows,&ncols);
712:     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);
713:   }
714:   PetscObjectReference((PetscObject)local);
715:   MatDestroy(&is->A);
716:   is->A = local;
717:   return(0);
718: }

722: /*@
723:     MatISSetLocalMat - Set the local matrix stored inside a MATIS matrix.

725:   Input Parameter:
726: .  mat - the matrix
727: .  local - the local matrix usually MATSEQAIJ

729:   Output Parameter:

731:   Level: advanced

733:   Notes:
734:     This can be called if you have precomputed the local matrix and
735:   want to provide it to the matrix object MATIS.

737: .seealso: MATIS
738: @*/
739: PetscErrorCode MatISSetLocalMat(Mat mat,Mat local)
740: {

746:   PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local));
747:   return(0);
748: }

752: PetscErrorCode MatZeroEntries_IS(Mat A)
753: {
754:   Mat_IS         *a = (Mat_IS*)A->data;

758:   MatZeroEntries(a->A);
759:   return(0);
760: }

764: PetscErrorCode MatScale_IS(Mat A,PetscScalar a)
765: {
766:   Mat_IS         *is = (Mat_IS*)A->data;

770:   MatScale(is->A,a);
771:   return(0);
772: }

776: PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v)
777: {
778:   Mat_IS         *is = (Mat_IS*)A->data;

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

785:   /* scatter diagonal back into global vector */
786:   VecSet(v,0);
787:   VecScatterBegin(is->ctx,is->x,v,ADD_VALUES,SCATTER_REVERSE);
788:   VecScatterEnd  (is->ctx,is->x,v,ADD_VALUES,SCATTER_REVERSE);
789:   return(0);
790: }

794: PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg)
795: {
796:   Mat_IS         *a = (Mat_IS*)A->data;

800:   MatSetOption(a->A,op,flg);
801:   return(0);
802: }

806: /*@
807:     MatCreateIS - Creates a "process" unassmembled matrix, it is assembled on each
808:        process but not across processes.

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

816:    Output Parameter:
817: .    A - the resulting matrix

819:    Level: advanced

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

826: .seealso: MATIS, MatSetLocalToGlobalMapping()
827: @*/
828: PetscErrorCode  MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping map,Mat *A)
829: {

833:   MatCreate(comm,A);
834:   MatSetBlockSize(*A,bs);
835:   MatSetSizes(*A,m,n,M,N);
836:   MatSetType(*A,MATIS);
837:   MatSetUp(*A);
838:   MatSetLocalToGlobalMapping(*A,map,map);
839:   return(0);
840: }

842: /*MC
843:    MATIS - MATIS = "is" - A matrix type to be used for using the Neumann-Neumann type preconditioners, see PCBDDC.
844:    This stores the matrices in globally unassembled form. Each processor
845:    assembles only its local Neumann problem and the parallel matrix vector
846:    product is handled "implicitly".

848:    Operations Provided:
849: +  MatMult()
850: .  MatMultAdd()
851: .  MatMultTranspose()
852: .  MatMultTransposeAdd()
853: .  MatZeroEntries()
854: .  MatSetOption()
855: .  MatZeroRows()
856: .  MatZeroRowsLocal()
857: .  MatSetValues()
858: .  MatSetValuesLocal()
859: .  MatScale()
860: .  MatGetDiagonal()
861: -  MatSetLocalToGlobalMapping()

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

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

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

870:           You can do matrix preallocation on the local matrix after you obtain it with
871:           MatISGetLocalMat()

873:   Level: advanced

875: .seealso: PC, MatISGetLocalMat(), MatSetLocalToGlobalMapping(), PCBDDC

877: M*/

881: PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A)
882: {
884:   Mat_IS         *b;

887:   PetscNewLog(A,&b);
888:   A->data = (void*)b;
889:   PetscMemzero(A->ops,sizeof(struct _MatOps));

891:   A->ops->mult                    = MatMult_IS;
892:   A->ops->multadd                 = MatMultAdd_IS;
893:   A->ops->multtranspose           = MatMultTranspose_IS;
894:   A->ops->multtransposeadd        = MatMultTransposeAdd_IS;
895:   A->ops->destroy                 = MatDestroy_IS;
896:   A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS;
897:   A->ops->setvalues               = MatSetValues_IS;
898:   A->ops->setvalueslocal          = MatSetValuesLocal_IS;
899:   A->ops->setvaluesblockedlocal   = MatSetValuesBlockedLocal_IS;
900:   A->ops->zerorows                = MatZeroRows_IS;
901:   A->ops->zerorowslocal           = MatZeroRowsLocal_IS;
902:   A->ops->assemblybegin           = MatAssemblyBegin_IS;
903:   A->ops->assemblyend             = MatAssemblyEnd_IS;
904:   A->ops->view                    = MatView_IS;
905:   A->ops->zeroentries             = MatZeroEntries_IS;
906:   A->ops->scale                   = MatScale_IS;
907:   A->ops->getdiagonal             = MatGetDiagonal_IS;
908:   A->ops->setoption               = MatSetOption_IS;
909:   A->ops->ishermitian             = MatIsHermitian_IS;
910:   A->ops->issymmetric             = MatIsSymmetric_IS;
911:   A->ops->duplicate               = MatDuplicate_IS;

913:   PetscLayoutSetUp(A->rmap);
914:   PetscLayoutSetUp(A->cmap);

916:   /* zeros pointer */
917:   b->A       = 0;
918:   b->ctx     = 0;
919:   b->x       = 0;
920:   b->y       = 0;
921:   b->mapping = 0;

923:   /* special MATIS functions */
924:   PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS);
925:   PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS);
926:   PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",MatISGetMPIXAIJ_IS);
927:   PetscObjectChangeTypeName((PetscObject)A,MATIS);
928:   return(0);
929: }