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: }