Actual source code: matis.c

petsc-dev 2014-04-15
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 MatIsHermitian_IS(Mat A,PetscReal tol,PetscBool  *flg)
 20: {
 22:   Mat_IS         *matis = (Mat_IS*)A->data;
 23:   PetscBool      local_sym;

 26:   MatIsHermitian(matis->A,tol,&local_sym);
 27:   MPI_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
 28:   return(0);
 29: }

 33: PetscErrorCode MatIsSymmetric_IS(Mat A,PetscReal tol,PetscBool  *flg)
 34: {
 36:   Mat_IS         *matis = (Mat_IS*)A->data;
 37:   PetscBool      local_sym;

 40:   MatIsSymmetric(matis->A,tol,&local_sym);
 41:   MPI_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
 42:   return(0);
 43: }

 47: PetscErrorCode MatDestroy_IS(Mat A)
 48: {
 50:   Mat_IS         *b = (Mat_IS*)A->data;

 53:   MatDestroy(&b->A);
 54:   VecScatterDestroy(&b->ctx);
 55:   VecDestroy(&b->x);
 56:   VecDestroy(&b->y);
 57:   ISLocalToGlobalMappingDestroy(&b->mapping);
 58:   PetscFree(A->data);
 59:   PetscObjectChangeTypeName((PetscObject)A,0);
 60:   PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL);
 61:   return(0);
 62: }

 66: PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y)
 67: {
 69:   Mat_IS         *is  = (Mat_IS*)A->data;
 70:   PetscScalar    zero = 0.0;

 73:   /*  scatter the global vector x into the local work vector */
 74:   VecScatterBegin(is->ctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
 75:   VecScatterEnd(is->ctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);

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

 80:   /* scatter product back into global memory */
 81:   VecSet(y,zero);
 82:   VecScatterBegin(is->ctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
 83:   VecScatterEnd(is->ctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
 84:   return(0);
 85: }

 89: static PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
 90: {
 91:   Vec            temp_vec;

 95:   if (v3 != v2) {
 96:     MatMult(A,v1,v3);
 97:     VecAXPY(v3,1.0,v2);
 98:   } else {
 99:     VecDuplicate(v2,&temp_vec);
100:     MatMult(A,v1,temp_vec);
101:     VecAXPY(temp_vec,1.0,v2);
102:     VecCopy(temp_vec,v3);
103:     VecDestroy(&temp_vec);
104:   }
105:   return(0);
106: }

110: PetscErrorCode MatMultTranspose_IS(Mat A,Vec x,Vec y)
111: {
112:   Mat_IS         *is = (Mat_IS*)A->data;

116:   /*  scatter the global vector x into the local work vector */
117:   VecScatterBegin(is->ctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
118:   VecScatterEnd(is->ctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);

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

123:   /* scatter product back into global vector */
124:   VecSet(y,0);
125:   VecScatterBegin(is->ctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
126:   VecScatterEnd(is->ctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
127:   return(0);
128: }

132: PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
133: {
134:   Vec            temp_vec;

138:   if (v3 != v2) {
139:     MatMultTranspose(A,v1,v3);
140:     VecAXPY(v3,1.0,v2);
141:   } else {
142:     VecDuplicate(v2,&temp_vec);
143:     MatMultTranspose(A,v1,temp_vec);
144:     VecAXPY(temp_vec,1.0,v2);
145:     VecCopy(temp_vec,v3);
146:     VecDestroy(&temp_vec);
147:   }
148:   return(0);
149: }

153: PetscErrorCode MatView_IS(Mat A,PetscViewer viewer)
154: {
155:   Mat_IS         *a = (Mat_IS*)A->data;
157:   PetscViewer    sviewer;

160:   PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
161:   PetscViewerGetSingleton(viewer,&sviewer);
162:   MatView(a->A,sviewer);
163:   PetscViewerRestoreSingleton(viewer,&sviewer);
164:   return(0);
165: }

169: PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping)
170: {
172:   PetscInt       n,bs;
173:   Mat_IS         *is = (Mat_IS*)A->data;
174:   IS             from,to;
175:   Vec            global;

178:   if (is->mapping) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Mapping already set for matrix");
180:   if (rmapping != cmapping) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"MATIS requires the row and column mappings to be identical");
181:   PetscObjectReference((PetscObject)rmapping);
182:   ISLocalToGlobalMappingDestroy(&is->mapping);
183:   is->mapping = rmapping;
184: /*
185:   PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping);
186:   PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping);
187: */

189:   /* Create the local matrix A */
190:   ISLocalToGlobalMappingGetSize(rmapping,&n);
191:   MatGetBlockSize(A,&bs);
192:   MatCreate(PETSC_COMM_SELF,&is->A);
193:   MatSetSizes(is->A,n,n,n,n);
194:   MatSetBlockSize(is->A,bs);
195:   MatSetOptionsPrefix(is->A,((PetscObject)A)->prefix);
196:   MatAppendOptionsPrefix(is->A,"is_");
197:   MatSetFromOptions(is->A);

199:   /* Create the local work vectors */
200:   VecCreate(PETSC_COMM_SELF,&is->x);
201:   VecSetBlockSize(is->x,bs);
202:   VecSetSizes(is->x,n,n);
203:   VecSetOptionsPrefix(is->x,((PetscObject)A)->prefix);
204:   VecAppendOptionsPrefix(is->x,"is_");
205:   VecSetFromOptions(is->x);
206:   VecDuplicate(is->x,&is->y);

208:   /* setup the global to local scatter */
209:   ISCreateStride(PETSC_COMM_SELF,n,0,1,&to);
210:   ISLocalToGlobalMappingApplyIS(rmapping,to,&from);
211:   MatGetVecs(A,&global,NULL);
212:   VecScatterCreate(global,from,is->x,to,&is->ctx);
213:   VecDestroy(&global);
214:   ISDestroy(&to);
215:   ISDestroy(&from);
216:   return(0);
217: }

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

228: #if defined(PETSC_USE_DEBUG)
229:   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);
230: #endif
231:   ISG2LMapApply(is->mapping,m,rows,rows_l);
232:   ISG2LMapApply(is->mapping,n,cols,cols_l);
233:   MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);
234:   return(0);
235: }

237: #undef ISG2LMapSetUp
238: #undef ISG2LMapApply

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

248:   MatSetValues(is->A,m,rows,n,cols,values,addv);
249:   return(0);
250: }

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

260:   MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv);
261:   return(0);
262: }

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

273:   if (x && b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support");
274:   if (n) {
275:     PetscMalloc1(n,&rows_l);
276:     ISGlobalToLocalMappingApply(is->mapping,IS_GTOLM_DROP,n,rows,&n_l,rows_l);
277:   }
278:   MatZeroRowsLocal(A,n_l,rows_l,diag,x,b);
279:   PetscFree(rows_l);
280:   return(0);
281: }

285: PetscErrorCode MatZeroRowsLocal_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
286: {
287:   Mat_IS         *is = (Mat_IS*)A->data;
289:   PetscInt       i;
290:   PetscScalar    *array;

293:   if (x && b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support");
294:   {
295:     /*
296:        Set up is->x as a "counting vector". This is in order to MatMult_IS
297:        work properly in the interface nodes.
298:     */
299:     Vec         counter;
300:     PetscScalar one=1.0, zero=0.0;
301:     MatGetVecs(A,&counter,NULL);
302:     VecSet(counter,zero);
303:     VecSet(is->x,one);
304:     VecScatterBegin(is->ctx,is->x,counter,ADD_VALUES,SCATTER_REVERSE);
305:     VecScatterEnd  (is->ctx,is->x,counter,ADD_VALUES,SCATTER_REVERSE);
306:     VecScatterBegin(is->ctx,counter,is->x,INSERT_VALUES,SCATTER_FORWARD);
307:     VecScatterEnd  (is->ctx,counter,is->x,INSERT_VALUES,SCATTER_FORWARD);
308:     VecDestroy(&counter);
309:   }
310:   if (!n) {
311:     is->pure_neumann = PETSC_TRUE;
312:   } else {
313:     is->pure_neumann = PETSC_FALSE;

315:     VecGetArray(is->x,&array);
316:     MatZeroRows(is->A,n,rows,diag,0,0);
317:     for (i=0; i<n; i++) {
318:       MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);
319:     }
320:     MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);
321:     MatAssemblyEnd  (is->A,MAT_FINAL_ASSEMBLY);
322:     VecRestoreArray(is->x,&array);
323:   }
324:   return(0);
325: }

329: PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type)
330: {
331:   Mat_IS         *is = (Mat_IS*)A->data;

335:   MatAssemblyBegin(is->A,type);
336:   return(0);
337: }

341: PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type)
342: {
343:   Mat_IS         *is = (Mat_IS*)A->data;

347:   MatAssemblyEnd(is->A,type);
348:   return(0);
349: }

353: PetscErrorCode  MatISGetLocalMat_IS(Mat mat,Mat *local)
354: {
355:   Mat_IS *is = (Mat_IS*)mat->data;

358:   *local = is->A;
359:   return(0);
360: }

364: /*@
365:     MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix.

367:   Input Parameter:
368: .  mat - the matrix

370:   Output Parameter:
371: .  local - the local matrix usually MATSEQAIJ

373:   Level: advanced

375:   Notes:
376:     This can be called if you have precomputed the nonzero structure of the
377:   matrix and want to provide it to the inner matrix object to improve the performance
378:   of the MatSetValues() operation.

380: .seealso: MATIS
381: @*/
382: PetscErrorCode  MatISGetLocalMat(Mat mat,Mat *local)
383: {

389:   PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat*),(mat,local));
390:   return(0);
391: }

395: PetscErrorCode  MatISSetLocalMat_IS(Mat mat,Mat local)
396: {
397:   Mat_IS         *is = (Mat_IS*)mat->data;
398:   PetscInt       nrows,ncols,orows,ocols;

402:   if (is->A) {
403:     MatGetSize(is->A,&orows,&ocols);
404:     MatGetSize(local,&nrows,&ncols);
405:     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);
406:   }
407:   PetscObjectReference((PetscObject)local);
408:   MatDestroy(&is->A);
409:   is->A = local;
410:   return(0);
411: }

415: /*@
416:     MatISSetLocalMat - Set the local matrix stored inside a MATIS matrix.

418:   Input Parameter:
419: .  mat - the matrix
420: .  local - the local matrix usually MATSEQAIJ

422:   Output Parameter:

424:   Level: advanced

426:   Notes:
427:     This can be called if you have precomputed the local matrix and
428:   want to provide it to the matrix object MATIS.

430: .seealso: MATIS
431: @*/
432: PetscErrorCode  MatISSetLocalMat(Mat mat,Mat local)
433: {

438:   PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local));
439:   return(0);
440: }

444: PetscErrorCode MatZeroEntries_IS(Mat A)
445: {
446:   Mat_IS         *a = (Mat_IS*)A->data;

450:   MatZeroEntries(a->A);
451:   return(0);
452: }

456: PetscErrorCode MatScale_IS(Mat A,PetscScalar a)
457: {
458:   Mat_IS         *is = (Mat_IS*)A->data;

462:   MatScale(is->A,a);
463:   return(0);
464: }

468: PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v)
469: {
470:   Mat_IS         *is = (Mat_IS*)A->data;

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

477:   /* scatter diagonal back into global vector */
478:   VecSet(v,0);
479:   VecScatterBegin(is->ctx,is->x,v,ADD_VALUES,SCATTER_REVERSE);
480:   VecScatterEnd  (is->ctx,is->x,v,ADD_VALUES,SCATTER_REVERSE);
481:   return(0);
482: }

486: PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg)
487: {
488:   Mat_IS         *a = (Mat_IS*)A->data;

492:   MatSetOption(a->A,op,flg);
493:   return(0);
494: }

498: /*@
499:     MatCreateIS - Creates a "process" unassmembled matrix, it is assembled on each
500:        process but not across processes.

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

508:    Output Parameter:
509: .    A - the resulting matrix

511:    Level: advanced

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

518: .seealso: MATIS, MatSetLocalToGlobalMapping()
519: @*/
520: PetscErrorCode  MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping map,Mat *A)
521: {

525:   MatCreate(comm,A);
526:   MatSetBlockSize(*A,bs);
527:   MatSetSizes(*A,m,n,M,N);
528:   MatSetType(*A,MATIS);
529:   MatSetUp(*A);
530:   MatSetLocalToGlobalMapping(*A,map,map);
531:   return(0);
532: }

534: /*MC
535:    MATIS - MATIS = "is" - A matrix type to be used for using the Neumann-Neumann type preconditioners.
536:    This stores the matrices in globally unassembled form. Each processor
537:    assembles only its local Neumann problem and the parallel matrix vector
538:    product is handled "implicitly".

540:    Operations Provided:
541: +  MatMult()
542: .  MatMultAdd()
543: .  MatMultTranspose()
544: .  MatMultTransposeAdd()
545: .  MatZeroEntries()
546: .  MatSetOption()
547: .  MatZeroRows()
548: .  MatZeroRowsLocal()
549: .  MatSetValues()
550: .  MatSetValuesLocal()
551: .  MatScale()
552: .  MatGetDiagonal()
553: -  MatSetLocalToGlobalMapping()

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

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

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

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

565:   Level: advanced

567: .seealso: PC, MatISGetLocalMat(), MatSetLocalToGlobalMapping()

569: M*/

573: PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A)
574: {
576:   Mat_IS         *b;

579:   PetscNewLog(A,&b);
580:   A->data = (void*)b;
581:   PetscMemzero(A->ops,sizeof(struct _MatOps));

583:   A->ops->mult                    = MatMult_IS;
584:   A->ops->multadd                 = MatMultAdd_IS;
585:   A->ops->multtranspose           = MatMultTranspose_IS;
586:   A->ops->multtransposeadd        = MatMultTransposeAdd_IS;
587:   A->ops->destroy                 = MatDestroy_IS;
588:   A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS;
589:   A->ops->setvalues               = MatSetValues_IS;
590:   A->ops->setvalueslocal          = MatSetValuesLocal_IS;
591:   A->ops->setvaluesblockedlocal   = MatSetValuesBlockedLocal_IS;
592:   A->ops->zerorows                = MatZeroRows_IS;
593:   A->ops->zerorowslocal           = MatZeroRowsLocal_IS;
594:   A->ops->assemblybegin           = MatAssemblyBegin_IS;
595:   A->ops->assemblyend             = MatAssemblyEnd_IS;
596:   A->ops->view                    = MatView_IS;
597:   A->ops->zeroentries             = MatZeroEntries_IS;
598:   A->ops->scale                   = MatScale_IS;
599:   A->ops->getdiagonal             = MatGetDiagonal_IS;
600:   A->ops->setoption               = MatSetOption_IS;
601:   A->ops->ishermitian             = MatIsHermitian_IS;
602:   A->ops->issymmetric             = MatIsSymmetric_IS;

604:   PetscLayoutSetUp(A->rmap);
605:   PetscLayoutSetUp(A->cmap);

607:   b->A   = 0;
608:   b->ctx = 0;
609:   b->x   = 0;
610:   b->y   = 0;
611:   PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS);
612:   PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS);
613:   PetscObjectChangeTypeName((PetscObject)A,MATIS);
614:   return(0);
615: }