Actual source code: matis.c

petsc-3.4.4 2014-03-13
  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 MatDestroy_IS(Mat A)
 20: {
 22:   Mat_IS         *b = (Mat_IS*)A->data;

 25:   MatDestroy(&b->A);
 26:   VecScatterDestroy(&b->ctx);
 27:   VecDestroy(&b->x);
 28:   VecDestroy(&b->y);
 29:   ISLocalToGlobalMappingDestroy(&b->mapping);
 30:   PetscFree(A->data);
 31:   PetscObjectChangeTypeName((PetscObject)A,0);
 32:   PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL);
 33:   return(0);
 34: }

 38: PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y)
 39: {
 41:   Mat_IS         *is  = (Mat_IS*)A->data;
 42:   PetscScalar    zero = 0.0;

 45:   /*  scatter the global vector x into the local work vector */
 46:   VecScatterBegin(is->ctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
 47:   VecScatterEnd(is->ctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);

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

 52:   /* scatter product back into global memory */
 53:   VecSet(y,zero);
 54:   VecScatterBegin(is->ctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
 55:   VecScatterEnd(is->ctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
 56:   return(0);
 57: }

 61: static PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
 62: {
 63:   Vec            temp_vec;

 67:   if (v3 != v2) {
 68:     MatMult(A,v1,v3);
 69:     VecAXPY(v3,1.0,v2);
 70:   } else {
 71:     VecDuplicate(v2,&temp_vec);
 72:     MatMult(A,v1,temp_vec);
 73:     VecAXPY(temp_vec,1.0,v2);
 74:     VecCopy(temp_vec,v3);
 75:     VecDestroy(&temp_vec);
 76:   }
 77:   return(0);
 78: }

 82: PetscErrorCode MatMultTranspose_IS(Mat A,Vec x,Vec y)
 83: {
 84:   Mat_IS         *is = (Mat_IS*)A->data;

 88:   /*  scatter the global vector x into the local work vector */
 89:   VecScatterBegin(is->ctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
 90:   VecScatterEnd(is->ctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);

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

 95:   /* scatter product back into global vector */
 96:   VecSet(y,0);
 97:   VecScatterBegin(is->ctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
 98:   VecScatterEnd(is->ctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
 99:   return(0);
100: }

104: PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
105: {
106:   Vec            temp_vec;

110:   if (v3 != v2) {
111:     MatMultTranspose(A,v1,v3);
112:     VecAXPY(v3,1.0,v2);
113:   } else {
114:     VecDuplicate(v2,&temp_vec);
115:     MatMultTranspose(A,v1,temp_vec);
116:     VecAXPY(temp_vec,1.0,v2);
117:     VecCopy(temp_vec,v3);
118:     VecDestroy(&temp_vec);
119:   }
120:   return(0);
121: }

125: PetscErrorCode MatView_IS(Mat A,PetscViewer viewer)
126: {
127:   Mat_IS         *a = (Mat_IS*)A->data;
129:   PetscViewer    sviewer;

132:   PetscViewerGetSingleton(viewer,&sviewer);
133:   MatView(a->A,sviewer);
134:   PetscViewerRestoreSingleton(viewer,&sviewer);
135:   return(0);
136: }

140: PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping)
141: {
143:   PetscInt       n,bs;
144:   Mat_IS         *is = (Mat_IS*)A->data;
145:   IS             from,to;
146:   Vec            global;

149:   if (is->mapping) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Mapping already set for matrix");
151:   if (rmapping != cmapping) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"MATIS requires the row and column mappings to be identical");
152:   PetscObjectReference((PetscObject)rmapping);
153:   ISLocalToGlobalMappingDestroy(&is->mapping);
154:   is->mapping = rmapping;

156:   /* Create the local matrix A */
157:   ISLocalToGlobalMappingGetSize(rmapping,&n);
158:   MatGetBlockSize(A,&bs);
159:   MatCreate(PETSC_COMM_SELF,&is->A);
160:   MatSetSizes(is->A,n,n,n,n);
161:   MatSetBlockSize(is->A,bs);
162:   MatSetOptionsPrefix(is->A,((PetscObject)A)->prefix);
163:   MatAppendOptionsPrefix(is->A,"is_");
164:   MatSetFromOptions(is->A);

166:   /* Create the local work vectors */
167:   VecCreate(PETSC_COMM_SELF,&is->x);
168:   VecSetBlockSize(is->x,bs);
169:   VecSetSizes(is->x,n,n);
170:   VecSetOptionsPrefix(is->x,((PetscObject)A)->prefix);
171:   VecAppendOptionsPrefix(is->x,"is_");
172:   VecSetFromOptions(is->x);
173:   VecDuplicate(is->x,&is->y);

175:   /* setup the global to local scatter */
176:   ISCreateStride(PETSC_COMM_SELF,n,0,1,&to);
177:   ISLocalToGlobalMappingApplyIS(rmapping,to,&from);
178:   MatGetVecs(A,&global,NULL);
179:   VecScatterCreate(global,from,is->x,to,&is->ctx);
180:   VecDestroy(&global);
181:   ISDestroy(&to);
182:   ISDestroy(&from);
183:   return(0);
184: }

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

195: #if defined(PETSC_USE_DEBUG)
196:   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);
197: #endif
198:   ISG2LMapApply(is->mapping,m,rows,rows_l);
199:   ISG2LMapApply(is->mapping,n,cols,cols_l);
200:   MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);
201:   return(0);
202: }

204: #undef ISG2LMapSetUp
205: #undef ISG2LMapApply

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

215:   MatSetValues(is->A,m,rows,n,cols,values,addv);
216:   return(0);
217: }

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

227:   MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv);
228:   return(0);
229: }

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

240:   if (x && b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support");
241:   if (n) {
242:     PetscMalloc(n*sizeof(PetscInt),&rows_l);
243:     ISGlobalToLocalMappingApply(is->mapping,IS_GTOLM_DROP,n,rows,&n_l,rows_l);
244:   }
245:   MatZeroRowsLocal(A,n_l,rows_l,diag,x,b);
246:   PetscFree(rows_l);
247:   return(0);
248: }

252: PetscErrorCode MatZeroRowsLocal_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
253: {
254:   Mat_IS         *is = (Mat_IS*)A->data;
256:   PetscInt       i;
257:   PetscScalar    *array;

260:   if (x && b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support");
261:   {
262:     /*
263:        Set up is->x as a "counting vector". This is in order to MatMult_IS
264:        work properly in the interface nodes.
265:     */
266:     Vec         counter;
267:     PetscScalar one=1.0, zero=0.0;
268:     MatGetVecs(A,&counter,NULL);
269:     VecSet(counter,zero);
270:     VecSet(is->x,one);
271:     VecScatterBegin(is->ctx,is->x,counter,ADD_VALUES,SCATTER_REVERSE);
272:     VecScatterEnd  (is->ctx,is->x,counter,ADD_VALUES,SCATTER_REVERSE);
273:     VecScatterBegin(is->ctx,counter,is->x,INSERT_VALUES,SCATTER_FORWARD);
274:     VecScatterEnd  (is->ctx,counter,is->x,INSERT_VALUES,SCATTER_FORWARD);
275:     VecDestroy(&counter);
276:   }
277:   if (!n) {
278:     is->pure_neumann = PETSC_TRUE;
279:   } else {
280:     is->pure_neumann = PETSC_FALSE;

282:     VecGetArray(is->x,&array);
283:     MatZeroRows(is->A,n,rows,diag,0,0);
284:     for (i=0; i<n; i++) {
285:       MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);
286:     }
287:     MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);
288:     MatAssemblyEnd  (is->A,MAT_FINAL_ASSEMBLY);
289:     VecRestoreArray(is->x,&array);
290:   }
291:   return(0);
292: }

296: PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type)
297: {
298:   Mat_IS         *is = (Mat_IS*)A->data;

302:   MatAssemblyBegin(is->A,type);
303:   return(0);
304: }

308: PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type)
309: {
310:   Mat_IS         *is = (Mat_IS*)A->data;

314:   MatAssemblyEnd(is->A,type);
315:   return(0);
316: }

320: PetscErrorCode  MatISGetLocalMat_IS(Mat mat,Mat *local)
321: {
322:   Mat_IS *is = (Mat_IS*)mat->data;

325:   *local = is->A;
326:   return(0);
327: }

331: /*@
332:     MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix.

334:   Input Parameter:
335: .  mat - the matrix

337:   Output Parameter:
338: .  local - the local matrix usually MATSEQAIJ

340:   Level: advanced

342:   Notes:
343:     This can be called if you have precomputed the nonzero structure of the
344:   matrix and want to provide it to the inner matrix object to improve the performance
345:   of the MatSetValues() operation.

347: .seealso: MATIS
348: @*/
349: PetscErrorCode  MatISGetLocalMat(Mat mat,Mat *local)
350: {

356:   PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat*),(mat,local));
357:   return(0);
358: }

362: PetscErrorCode  MatISSetLocalMat_IS(Mat mat,Mat local)
363: {
364:   Mat_IS         *is = (Mat_IS*)mat->data;
365:   PetscInt       nrows,ncols,orows,ocols;

369:   if (is->A) {
370:     MatGetSize(is->A,&orows,&ocols);
371:     MatGetSize(local,&nrows,&ncols);
372:     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);
373:   }
374:   PetscObjectReference((PetscObject)local);
375:   MatDestroy(&is->A);
376:   is->A = local;
377:   return(0);
378: }

382: /*@
383:     MatISSetLocalMat - Set the local matrix stored inside a MATIS matrix.

385:   Input Parameter:
386: .  mat - the matrix
387: .  local - the local matrix usually MATSEQAIJ

389:   Output Parameter:

391:   Level: advanced

393:   Notes:
394:     This can be called if you have precomputed the local matrix and
395:   want to provide it to the matrix object MATIS.

397: .seealso: MATIS
398: @*/
399: PetscErrorCode  MatISSetLocalMat(Mat mat,Mat local)
400: {

405:   PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local));
406:   return(0);
407: }

411: PetscErrorCode MatZeroEntries_IS(Mat A)
412: {
413:   Mat_IS         *a = (Mat_IS*)A->data;

417:   MatZeroEntries(a->A);
418:   return(0);
419: }

423: PetscErrorCode MatScale_IS(Mat A,PetscScalar a)
424: {
425:   Mat_IS         *is = (Mat_IS*)A->data;

429:   MatScale(is->A,a);
430:   return(0);
431: }

435: PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v)
436: {
437:   Mat_IS         *is = (Mat_IS*)A->data;

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

444:   /* scatter diagonal back into global vector */
445:   VecSet(v,0);
446:   VecScatterBegin(is->ctx,is->x,v,ADD_VALUES,SCATTER_REVERSE);
447:   VecScatterEnd  (is->ctx,is->x,v,ADD_VALUES,SCATTER_REVERSE);
448:   return(0);
449: }

453: PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg)
454: {
455:   Mat_IS         *a = (Mat_IS*)A->data;

459:   MatSetOption(a->A,op,flg);
460:   return(0);
461: }

465: /*@
466:     MatCreateIS - Creates a "process" unassmembled matrix, it is assembled on each
467:        process but not across processes.

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

475:    Output Parameter:
476: .    A - the resulting matrix

478:    Level: advanced

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

485: .seealso: MATIS, MatSetLocalToGlobalMapping()
486: @*/
487: PetscErrorCode  MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping map,Mat *A)
488: {

492:   MatCreate(comm,A);
493:   MatSetBlockSize(*A,bs);
494:   MatSetSizes(*A,m,n,M,N);
495:   MatSetType(*A,MATIS);
496:   MatSetUp(*A);
497:   MatSetLocalToGlobalMapping(*A,map,map);
498:   return(0);
499: }

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

507:    Operations Provided:
508: +  MatMult()
509: .  MatMultAdd()
510: .  MatMultTranspose()
511: .  MatMultTransposeAdd()
512: .  MatZeroEntries()
513: .  MatSetOption()
514: .  MatZeroRows()
515: .  MatZeroRowsLocal()
516: .  MatSetValues()
517: .  MatSetValuesLocal()
518: .  MatScale()
519: .  MatGetDiagonal()
520: -  MatSetLocalToGlobalMapping()

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

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

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

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

532:   Level: advanced

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

536: M*/

540: PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A)
541: {
543:   Mat_IS         *b;

546:   PetscNewLog(A,Mat_IS,&b);
547:   A->data = (void*)b;
548:   PetscMemzero(A->ops,sizeof(struct _MatOps));

550:   A->ops->mult                    = MatMult_IS;
551:   A->ops->multadd                 = MatMultAdd_IS;
552:   A->ops->multtranspose           = MatMultTranspose_IS;
553:   A->ops->multtransposeadd        = MatMultTransposeAdd_IS;
554:   A->ops->destroy                 = MatDestroy_IS;
555:   A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS;
556:   A->ops->setvalues               = MatSetValues_IS;
557:   A->ops->setvalueslocal          = MatSetValuesLocal_IS;
558:   A->ops->setvaluesblockedlocal   = MatSetValuesBlockedLocal_IS;
559:   A->ops->zerorows                = MatZeroRows_IS;
560:   A->ops->zerorowslocal           = MatZeroRowsLocal_IS;
561:   A->ops->assemblybegin           = MatAssemblyBegin_IS;
562:   A->ops->assemblyend             = MatAssemblyEnd_IS;
563:   A->ops->view                    = MatView_IS;
564:   A->ops->zeroentries             = MatZeroEntries_IS;
565:   A->ops->scale                   = MatScale_IS;
566:   A->ops->getdiagonal             = MatGetDiagonal_IS;
567:   A->ops->setoption               = MatSetOption_IS;

569:   PetscLayoutSetUp(A->rmap);
570:   PetscLayoutSetUp(A->cmap);

572:   b->A   = 0;
573:   b->ctx = 0;
574:   b->x   = 0;
575:   b->y   = 0;
576:   PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS);
577:   PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS);
578:   PetscObjectChangeTypeName((PetscObject)A,MATIS);
579:   return(0);
580: }