Actual source code: matis.c

  1: #define PETSCMAT_DLL

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

  9:      We provide:
 10:          MatMult()

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

 14: */

 16:  #include src/mat/impls/is/matis.h

 20: PetscErrorCode MatDestroy_IS(Mat A)
 21: {
 23:   Mat_IS         *b = (Mat_IS*)A->data;

 26:   if (b->A) {
 27:     MatDestroy(b->A);
 28:   }
 29:   if (b->ctx) {
 30:     VecScatterDestroy(b->ctx);
 31:   }
 32:   if (b->x) {
 33:     VecDestroy(b->x);
 34:   }
 35:   if (b->y) {
 36:     VecDestroy(b->y);
 37:   }
 38:   if (b->mapping) {
 39:     ISLocalToGlobalMappingDestroy(b->mapping);
 40:   }
 41:   PetscFree(b);
 42:   PetscObjectChangeTypeName((PetscObject)A,0);
 43:   PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C","",PETSC_NULL);
 44:   return(0);
 45: }

 49: PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y)
 50: {
 52:   Mat_IS         *is = (Mat_IS*)A->data;
 53:   PetscScalar    zero = 0.0;

 56:   /*  scatter the global vector x into the local work vector */
 57:   VecScatterBegin(x,is->x,INSERT_VALUES,SCATTER_FORWARD,is->ctx);
 58:   VecScatterEnd(x,is->x,INSERT_VALUES,SCATTER_FORWARD,is->ctx);

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

 63:   /* scatter product back into global memory */
 64:   VecSet(y,zero);
 65:   VecScatterBegin(is->y,y,ADD_VALUES,SCATTER_REVERSE,is->ctx);
 66:   VecScatterEnd(is->y,y,ADD_VALUES,SCATTER_REVERSE,is->ctx);

 68:   return(0);
 69: }

 73: PetscErrorCode MatView_IS(Mat A,PetscViewer viewer)
 74: {
 75:   Mat_IS         *a = (Mat_IS*)A->data;
 77:   PetscViewer    sviewer;

 80:   PetscViewerGetSingleton(viewer,&sviewer);
 81:   MatView(a->A,sviewer);
 82:   PetscViewerRestoreSingleton(viewer,&sviewer);
 83:   return(0);
 84: }

 88: PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping mapping)
 89: {
 91:   PetscInt       n;
 92:   Mat_IS         *is = (Mat_IS*)A->data;
 93:   IS             from,to;
 94:   Vec            global;

 97:   is->mapping = mapping;
 98:   PetscObjectReference((PetscObject)mapping);

100:   /* Create the local matrix A */
101:   ISLocalToGlobalMappingGetSize(mapping,&n);
102:   MatCreate(PETSC_COMM_SELF,&is->A);
103:   MatSetSizes(is->A,n,n,n,n);
104:   PetscObjectSetOptionsPrefix((PetscObject)is->A,"is");
105:   MatSetFromOptions(is->A);

107:   /* Create the local work vectors */
108:   VecCreateSeq(PETSC_COMM_SELF,n,&is->x);
109:   VecDuplicate(is->x,&is->y);

111:   /* setup the global to local scatter */
112:   ISCreateStride(PETSC_COMM_SELF,n,0,1,&to);
113:   ISLocalToGlobalMappingApplyIS(mapping,to,&from);
114:   VecCreateMPI(A->comm,A->cmap.n,A->cmap.N,&global);
115:   VecScatterCreate(global,from,is->x,to,&is->ctx);
116:   VecDestroy(global);
117:   ISDestroy(to);
118:   ISDestroy(from);
119:   return(0);
120: }


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

131:   MatSetValues(is->A,m,rows,n,cols,values,addv);
132:   return(0);
133: }

137: PetscErrorCode MatZeroRowsLocal_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag)
138: {
139:   Mat_IS         *is = (Mat_IS*)A->data;
141:   PetscInt       i;
142:   PetscScalar    *array;

145:   {
146:     /*
147:        Set up is->x as a "counting vector". This is in order to MatMult_IS
148:        work properly in the interface nodes.
149:     */
150:     Vec         counter;
151:     PetscScalar one=1.0, zero=0.0;
152:     VecCreateMPI(A->comm,A->cmap.n,A->cmap.N,&counter);
153:     VecSet(counter,zero);
154:     VecSet(is->x,one);
155:     VecScatterBegin(is->x,counter,ADD_VALUES,SCATTER_REVERSE,is->ctx);
156:     VecScatterEnd  (is->x,counter,ADD_VALUES,SCATTER_REVERSE,is->ctx);
157:     VecScatterBegin(counter,is->x,INSERT_VALUES,SCATTER_FORWARD,is->ctx);
158:     VecScatterEnd  (counter,is->x,INSERT_VALUES,SCATTER_FORWARD,is->ctx);
159:     VecDestroy(counter);
160:   }
161:   if (!n) {
162:     is->pure_neumann = PETSC_TRUE;
163:   } else {
164:     is->pure_neumann = PETSC_FALSE;
165:     VecGetArray(is->x,&array);
166:     MatZeroRows(is->A,n,rows,diag);
167:     for (i=0; i<n; i++) {
168:       MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);
169:     }
170:     MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);
171:     MatAssemblyEnd  (is->A,MAT_FINAL_ASSEMBLY);
172:     VecRestoreArray(is->x,&array);
173:   }

175:   return(0);
176: }

180: PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type)
181: {
182:   Mat_IS         *is = (Mat_IS*)A->data;

186:   MatAssemblyBegin(is->A,type);
187:   return(0);
188: }

192: PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type)
193: {
194:   Mat_IS         *is = (Mat_IS*)A->data;

198:   MatAssemblyEnd(is->A,type);
199:   return(0);
200: }

205: PetscErrorCode  MatISGetLocalMat_IS(Mat mat,Mat *local)
206: {
207:   Mat_IS *is = (Mat_IS *)mat->data;
208: 
210:   *local = is->A;
211:   return(0);
212: }

217: /*@
218:     MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix.

220:   Input Parameter:
221: .  mat - the matrix

223:   Output Parameter:
224: .  local - the local matrix usually MATSEQAIJ

226:   Level: advanced

228:   Notes:
229:     This can be called if you have precomputed the nonzero structure of the 
230:   matrix and want to provide it to the inner matrix object to improve the performance
231:   of the MatSetValues() operation.

233: .seealso: MATIS
234: @*/
235: PetscErrorCode  MatISGetLocalMat(Mat mat,Mat *local)
236: {
237:   PetscErrorCode ierr,(*f)(Mat,Mat *);

242:   PetscObjectQueryFunction((PetscObject)mat,"MatISGetLocalMat_C",(void (**)(void))&f);
243:   if (f) {
244:     (*f)(mat,local);
245:   } else {
246:     local = 0;
247:   }
248:   return(0);
249: }

253: PetscErrorCode MatZeroEntries_IS(Mat A)
254: {
255:   Mat_IS         *a = (Mat_IS*)A->data;

259:   MatZeroEntries(a->A);
260:   return(0);
261: }

265: PetscErrorCode MatSetOption_IS(Mat A,MatOption op)
266: {
267:   Mat_IS         *a = (Mat_IS*)A->data;

271:   MatSetOption(a->A,op);
272:   return(0);
273: }

277: /*@
278:     MatCreateIS - Creates a "process" unassmembled matrix, it is assembled on each 
279:        process but not across processes.

281:    Input Parameters:
282: +     comm - MPI communicator that will share the matrix
283: .     m,n,M,N - local and/or global sizes of the the left and right vector used in matrix vector products
284: -     map - mapping that defines the global number for each local number

286:    Output Parameter:
287: .    A - the resulting matrix

289:    Level: advanced

291:    Notes: See MATIS for more details
292:           m and n are NOT related to the size of the map

294: .seealso: MATIS, MatSetLocalToGlobalMapping()
295: @*/
296: PetscErrorCode  MatCreateIS(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping map,Mat *A)
297: {

301:   MatCreate(comm,A);
302:   MatSetSizes(*A,m,n,M,N);
303:   MatSetType(*A,MATIS);
304:   MatSetLocalToGlobalMapping(*A,map);
305:   return(0);
306: }

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

314:    Operations Provided:
315: +  MatMult()
316: .  MatZeroEntries()
317: .  MatSetOption()
318: .  MatZeroRowsLocal()
319: .  MatSetValuesLocal()
320: -  MatSetLocalToGlobalMapping()

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

325:    Notes: Options prefix for the inner matrix are given by -is_mat_xxx
326:     
327:           You must call MatSetLocalToGlobalMapping() before using this matrix type.

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

332:   Level: advanced

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

336: M*/

341: PetscErrorCode  MatCreate_IS(Mat A)
342: {
344:   Mat_IS         *b;

347:   PetscNew(Mat_IS,&b);
348:   A->data             = (void*)b;
349:   PetscMemzero(A->ops,sizeof(struct _MatOps));
350:   A->factor           = 0;
351:   A->mapping          = 0;

353:   A->ops->mult                    = MatMult_IS;
354:   A->ops->destroy                 = MatDestroy_IS;
355:   A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS;
356:   A->ops->setvalueslocal          = MatSetValuesLocal_IS;
357:   A->ops->zerorowslocal           = MatZeroRowsLocal_IS;
358:   A->ops->assemblybegin           = MatAssemblyBegin_IS;
359:   A->ops->assemblyend             = MatAssemblyEnd_IS;
360:   A->ops->view                    = MatView_IS;
361:   A->ops->zeroentries             = MatZeroEntries_IS;
362:   A->ops->setoption               = MatSetOption_IS;

364:   PetscMapInitialize(A->comm,&A->rmap);
365:   PetscMapInitialize(A->comm,&A->cmap);

367:   b->A          = 0;
368:   b->ctx        = 0;
369:   b->x          = 0;
370:   b->y          = 0;
371:   PetscObjectComposeFunctionDynamic((PetscObject)A,"MatISGetLocalMat_C","MatISGetLocalMat_IS",MatISGetLocalMat_IS);
372:   PetscObjectChangeTypeName((PetscObject)A,MATIS);

374:   return(0);
375: }