Actual source code: matis.c

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

  8:      We provide:
  9:          MatMult()

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

 13: */

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

 19: PetscErrorCode MatDuplicate_IS(Mat mat,MatDuplicateOption op,Mat *newmat)
 20: {
 22:   Mat_IS         *matis = (Mat_IS*)(mat->data);
 23:   PetscInt       bs,m,n,M,N;
 24:   Mat            B,localmat;

 27:   MatGetBlockSize(mat,&bs);
 28:   MatGetSize(mat,&M,&N);
 29:   MatGetLocalSize(mat,&m,&n);
 30:   MatCreateIS(PetscObjectComm((PetscObject)mat),bs,m,n,M,N,matis->mapping,&B);
 31:   MatDuplicate(matis->A,op,&localmat);
 32:   MatISSetLocalMat(B,localmat);
 33:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 34:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
 35:   *newmat = B;
 36:   return(0);
 37: }

 41: PetscErrorCode MatIsHermitian_IS(Mat A,PetscReal tol,PetscBool  *flg)
 42: {
 44:   Mat_IS         *matis = (Mat_IS*)A->data;
 45:   PetscBool      local_sym;

 48:   MatIsHermitian(matis->A,tol,&local_sym);
 49:   MPI_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
 50:   return(0);
 51: }

 55: PetscErrorCode MatIsSymmetric_IS(Mat A,PetscReal tol,PetscBool  *flg)
 56: {
 58:   Mat_IS         *matis = (Mat_IS*)A->data;
 59:   PetscBool      local_sym;

 62:   MatIsSymmetric(matis->A,tol,&local_sym);
 63:   MPI_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
 64:   return(0);
 65: }

 69: PetscErrorCode MatDestroy_IS(Mat A)
 70: {
 72:   Mat_IS         *b = (Mat_IS*)A->data;

 75:   MatDestroy(&b->A);
 76:   VecScatterDestroy(&b->ctx);
 77:   VecDestroy(&b->x);
 78:   VecDestroy(&b->y);
 79:   ISLocalToGlobalMappingDestroy(&b->mapping);
 80:   PetscFree(A->data);
 81:   PetscObjectChangeTypeName((PetscObject)A,0);
 82:   PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL);
 83:   return(0);
 84: }

 88: PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y)
 89: {
 91:   Mat_IS         *is  = (Mat_IS*)A->data;
 92:   PetscScalar    zero = 0.0;

 95:   /*  scatter the global vector x into the local work vector */
 96:   VecScatterBegin(is->ctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
 97:   VecScatterEnd(is->ctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);

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

102:   /* scatter product back into global memory */
103:   VecSet(y,zero);
104:   VecScatterBegin(is->ctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
105:   VecScatterEnd(is->ctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
106:   return(0);
107: }

111: static PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
112: {
113:   Vec            temp_vec;

117:   if (v3 != v2) {
118:     MatMult(A,v1,v3);
119:     VecAXPY(v3,1.0,v2);
120:   } else {
121:     VecDuplicate(v2,&temp_vec);
122:     MatMult(A,v1,temp_vec);
123:     VecAXPY(temp_vec,1.0,v2);
124:     VecCopy(temp_vec,v3);
125:     VecDestroy(&temp_vec);
126:   }
127:   return(0);
128: }

132: PetscErrorCode MatMultTranspose_IS(Mat A,Vec x,Vec y)
133: {
134:   Mat_IS         *is = (Mat_IS*)A->data;

138:   /*  scatter the global vector x into the local work vector */
139:   VecScatterBegin(is->ctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
140:   VecScatterEnd(is->ctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);

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

145:   /* scatter product back into global vector */
146:   VecSet(y,0);
147:   VecScatterBegin(is->ctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
148:   VecScatterEnd(is->ctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
149:   return(0);
150: }

154: PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
155: {
156:   Vec            temp_vec;

160:   if (v3 != v2) {
161:     MatMultTranspose(A,v1,v3);
162:     VecAXPY(v3,1.0,v2);
163:   } else {
164:     VecDuplicate(v2,&temp_vec);
165:     MatMultTranspose(A,v1,temp_vec);
166:     VecAXPY(temp_vec,1.0,v2);
167:     VecCopy(temp_vec,v3);
168:     VecDestroy(&temp_vec);
169:   }
170:   return(0);
171: }

175: PetscErrorCode MatView_IS(Mat A,PetscViewer viewer)
176: {
177:   Mat_IS         *a = (Mat_IS*)A->data;
179:   PetscViewer    sviewer;

182:   PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
183:   PetscViewerGetSingleton(viewer,&sviewer);
184:   MatView(a->A,sviewer);
185:   PetscViewerRestoreSingleton(viewer,&sviewer);
186:   return(0);
187: }

191: PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping)
192: {
194:   PetscInt       n,bs;
195:   Mat_IS         *is = (Mat_IS*)A->data;
196:   IS             from,to;
197:   Vec            global;

201:   if (rmapping != cmapping) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"MATIS requires the row and column mappings to be identical");
202:   if (is->mapping) { /* Currenly destroys the objects that will be created by this routine. Is there anything else that should be checked? */
203:     ISLocalToGlobalMappingDestroy(&is->mapping);
204:     VecDestroy(&is->x);
205:     VecDestroy(&is->y);
206:     VecScatterDestroy(&is->ctx);
207:     MatDestroy(&is->A);
208:   }
209:   PetscObjectReference((PetscObject)rmapping);
210:   ISLocalToGlobalMappingDestroy(&is->mapping);
211:   is->mapping = rmapping;
212: /*
213:   PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping);
214:   PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping);
215: */

217:   /* Create the local matrix A */
218:   ISLocalToGlobalMappingGetSize(rmapping,&n);
219:   MatGetBlockSize(A,&bs);
220:   MatCreate(PETSC_COMM_SELF,&is->A);
221:   MatSetSizes(is->A,n,n,n,n);
222:   MatSetBlockSize(is->A,bs);
223:   MatSetOptionsPrefix(is->A,((PetscObject)A)->prefix);
224:   MatAppendOptionsPrefix(is->A,"is_");
225:   MatSetFromOptions(is->A);

227:   /* Create the local work vectors */
228:   VecCreate(PETSC_COMM_SELF,&is->x);
229:   VecSetBlockSize(is->x,bs);
230:   VecSetSizes(is->x,n,n);
231:   VecSetOptionsPrefix(is->x,((PetscObject)A)->prefix);
232:   VecAppendOptionsPrefix(is->x,"is_");
233:   VecSetFromOptions(is->x);
234:   VecDuplicate(is->x,&is->y);

236:   /* setup the global to local scatter */
237:   ISCreateStride(PETSC_COMM_SELF,n,0,1,&to);
238:   ISLocalToGlobalMappingApplyIS(rmapping,to,&from);
239:   MatGetVecs(A,&global,NULL);
240:   VecScatterCreate(global,from,is->x,to,&is->ctx);
241:   VecDestroy(&global);
242:   ISDestroy(&to);
243:   ISDestroy(&from);
244:   return(0);
245: }

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

256: #if defined(PETSC_USE_DEBUG)
257:   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);
258: #endif
259:   ISG2LMapApply(is->mapping,m,rows,rows_l);
260:   ISG2LMapApply(is->mapping,n,cols,cols_l);
261:   MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);
262:   return(0);
263: }

265: #undef ISG2LMapSetUp
266: #undef ISG2LMapApply

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

276:   MatSetValues(is->A,m,rows,n,cols,values,addv);
277:   return(0);
278: }

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

288:   MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv);
289:   return(0);
290: }

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

301:   if (x && b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support");
302:   if (n) {
303:     PetscMalloc1(n,&rows_l);
304:     ISGlobalToLocalMappingApply(is->mapping,IS_GTOLM_DROP,n,rows,&n_l,rows_l);
305:   }
306:   MatZeroRowsLocal(A,n_l,rows_l,diag,x,b);
307:   PetscFree(rows_l);
308:   return(0);
309: }

313: PetscErrorCode MatZeroRowsLocal_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
314: {
315:   Mat_IS         *is = (Mat_IS*)A->data;
317:   PetscInt       i;
318:   PetscScalar    *array;

321:   if (x && b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support");
322:   {
323:     /*
324:        Set up is->x as a "counting vector". This is in order to MatMult_IS
325:        work properly in the interface nodes.
326:     */
327:     Vec         counter;
328:     PetscScalar one=1.0, zero=0.0;
329:     MatGetVecs(A,&counter,NULL);
330:     VecSet(counter,zero);
331:     VecSet(is->x,one);
332:     VecScatterBegin(is->ctx,is->x,counter,ADD_VALUES,SCATTER_REVERSE);
333:     VecScatterEnd  (is->ctx,is->x,counter,ADD_VALUES,SCATTER_REVERSE);
334:     VecScatterBegin(is->ctx,counter,is->x,INSERT_VALUES,SCATTER_FORWARD);
335:     VecScatterEnd  (is->ctx,counter,is->x,INSERT_VALUES,SCATTER_FORWARD);
336:     VecDestroy(&counter);
337:   }
338:   if (!n) {
339:     is->pure_neumann = PETSC_TRUE;
340:   } else {
341:     is->pure_neumann = PETSC_FALSE;

343:     VecGetArray(is->x,&array);
344:     MatZeroRows(is->A,n,rows,diag,0,0);
345:     for (i=0; i<n; i++) {
346:       MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);
347:     }
348:     MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);
349:     MatAssemblyEnd  (is->A,MAT_FINAL_ASSEMBLY);
350:     VecRestoreArray(is->x,&array);
351:   }
352:   return(0);
353: }

357: PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type)
358: {
359:   Mat_IS         *is = (Mat_IS*)A->data;

363:   MatAssemblyBegin(is->A,type);
364:   return(0);
365: }

369: PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type)
370: {
371:   Mat_IS         *is = (Mat_IS*)A->data;

375:   MatAssemblyEnd(is->A,type);
376:   return(0);
377: }

381: PetscErrorCode  MatISGetLocalMat_IS(Mat mat,Mat *local)
382: {
383:   Mat_IS *is = (Mat_IS*)mat->data;

386:   *local = is->A;
387:   return(0);
388: }

392: /*@
393:     MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix.

395:   Input Parameter:
396: .  mat - the matrix

398:   Output Parameter:
399: .  local - the local matrix usually MATSEQAIJ

401:   Level: advanced

403:   Notes:
404:     This can be called if you have precomputed the nonzero structure of the
405:   matrix and want to provide it to the inner matrix object to improve the performance
406:   of the MatSetValues() operation.

408: .seealso: MATIS
409: @*/
410: PetscErrorCode  MatISGetLocalMat(Mat mat,Mat *local)
411: {

417:   PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat*),(mat,local));
418:   return(0);
419: }

423: PetscErrorCode  MatISSetLocalMat_IS(Mat mat,Mat local)
424: {
425:   Mat_IS         *is = (Mat_IS*)mat->data;
426:   PetscInt       nrows,ncols,orows,ocols;

430:   if (is->A) {
431:     MatGetSize(is->A,&orows,&ocols);
432:     MatGetSize(local,&nrows,&ncols);
433:     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);
434:   }
435:   PetscObjectReference((PetscObject)local);
436:   MatDestroy(&is->A);
437:   is->A = local;
438:   return(0);
439: }

443: /*@
444:     MatISSetLocalMat - Set the local matrix stored inside a MATIS matrix.

446:   Input Parameter:
447: .  mat - the matrix
448: .  local - the local matrix usually MATSEQAIJ

450:   Output Parameter:

452:   Level: advanced

454:   Notes:
455:     This can be called if you have precomputed the local matrix and
456:   want to provide it to the matrix object MATIS.

458: .seealso: MATIS
459: @*/
460: PetscErrorCode  MatISSetLocalMat(Mat mat,Mat local)
461: {

466:   PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local));
467:   return(0);
468: }

472: PetscErrorCode MatZeroEntries_IS(Mat A)
473: {
474:   Mat_IS         *a = (Mat_IS*)A->data;

478:   MatZeroEntries(a->A);
479:   return(0);
480: }

484: PetscErrorCode MatScale_IS(Mat A,PetscScalar a)
485: {
486:   Mat_IS         *is = (Mat_IS*)A->data;

490:   MatScale(is->A,a);
491:   return(0);
492: }

496: PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v)
497: {
498:   Mat_IS         *is = (Mat_IS*)A->data;

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

505:   /* scatter diagonal back into global vector */
506:   VecSet(v,0);
507:   VecScatterBegin(is->ctx,is->x,v,ADD_VALUES,SCATTER_REVERSE);
508:   VecScatterEnd  (is->ctx,is->x,v,ADD_VALUES,SCATTER_REVERSE);
509:   return(0);
510: }

514: PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg)
515: {
516:   Mat_IS         *a = (Mat_IS*)A->data;

520:   MatSetOption(a->A,op,flg);
521:   return(0);
522: }

526: /*@
527:     MatCreateIS - Creates a "process" unassmembled matrix, it is assembled on each
528:        process but not across processes.

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

536:    Output Parameter:
537: .    A - the resulting matrix

539:    Level: advanced

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

546: .seealso: MATIS, MatSetLocalToGlobalMapping()
547: @*/
548: PetscErrorCode  MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping map,Mat *A)
549: {

553:   MatCreate(comm,A);
554:   MatSetBlockSize(*A,bs);
555:   MatSetSizes(*A,m,n,M,N);
556:   MatSetType(*A,MATIS);
557:   MatSetUp(*A);
558:   MatSetLocalToGlobalMapping(*A,map,map);
559:   return(0);
560: }

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

568:    Operations Provided:
569: +  MatMult()
570: .  MatMultAdd()
571: .  MatMultTranspose()
572: .  MatMultTransposeAdd()
573: .  MatZeroEntries()
574: .  MatSetOption()
575: .  MatZeroRows()
576: .  MatZeroRowsLocal()
577: .  MatSetValues()
578: .  MatSetValuesLocal()
579: .  MatScale()
580: .  MatGetDiagonal()
581: -  MatSetLocalToGlobalMapping()

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

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

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

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

593:   Level: advanced

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

597: M*/

601: PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A)
602: {
604:   Mat_IS         *b;

607:   PetscNewLog(A,&b);
608:   A->data = (void*)b;
609:   PetscMemzero(A->ops,sizeof(struct _MatOps));

611:   A->ops->mult                    = MatMult_IS;
612:   A->ops->multadd                 = MatMultAdd_IS;
613:   A->ops->multtranspose           = MatMultTranspose_IS;
614:   A->ops->multtransposeadd        = MatMultTransposeAdd_IS;
615:   A->ops->destroy                 = MatDestroy_IS;
616:   A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS;
617:   A->ops->setvalues               = MatSetValues_IS;
618:   A->ops->setvalueslocal          = MatSetValuesLocal_IS;
619:   A->ops->setvaluesblockedlocal   = MatSetValuesBlockedLocal_IS;
620:   A->ops->zerorows                = MatZeroRows_IS;
621:   A->ops->zerorowslocal           = MatZeroRowsLocal_IS;
622:   A->ops->assemblybegin           = MatAssemblyBegin_IS;
623:   A->ops->assemblyend             = MatAssemblyEnd_IS;
624:   A->ops->view                    = MatView_IS;
625:   A->ops->zeroentries             = MatZeroEntries_IS;
626:   A->ops->scale                   = MatScale_IS;
627:   A->ops->getdiagonal             = MatGetDiagonal_IS;
628:   A->ops->setoption               = MatSetOption_IS;
629:   A->ops->ishermitian             = MatIsHermitian_IS;
630:   A->ops->issymmetric             = MatIsSymmetric_IS;
631:   A->ops->duplicate               = MatDuplicate_IS;

633:   PetscLayoutSetUp(A->rmap);
634:   PetscLayoutSetUp(A->cmap);

636:   b->A       = 0;
637:   b->ctx     = 0;
638:   b->x       = 0;
639:   b->y       = 0;
640:   b->mapping = 0;
641:   PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS);
642:   PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS);
643:   PetscObjectChangeTypeName((PetscObject)A,MATIS);
644:   return(0);
645: }