Actual source code: matis.c
petsc-3.3-p7 2013-05-11
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","",PETSC_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);
57: return(0);
58: }
62: static PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
63: {
65: Mat_IS *is = (Mat_IS*)A->data;
68: /* scatter the global vector v1 into the local work vector */
69: VecScatterBegin(is->ctx,v1,is->x,INSERT_VALUES,SCATTER_FORWARD);
70: VecScatterEnd (is->ctx,v1,is->x,INSERT_VALUES,SCATTER_FORWARD);
71: VecScatterBegin(is->ctx,v2,is->y,INSERT_VALUES,SCATTER_FORWARD);
72: VecScatterEnd (is->ctx,v2,is->y,INSERT_VALUES,SCATTER_FORWARD);
74: /* multiply the local matrix */
75: MatMultAdd(is->A,is->x,is->y,is->y);
77: /* scatter result back into global vector */
78: VecSet(v3,0);
79: VecScatterBegin(is->ctx,is->y,v3,ADD_VALUES,SCATTER_REVERSE);
80: VecScatterEnd (is->ctx,is->y,v3,ADD_VALUES,SCATTER_REVERSE);
81: return(0);
82: }
86: PetscErrorCode MatMultTranspose_IS(Mat A,Vec x,Vec y)
87: {
88: Mat_IS *is = (Mat_IS*)A->data;
92: /* scatter the global vector x into the local work vector */
93: VecScatterBegin(is->ctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
94: VecScatterEnd(is->ctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
96: /* multiply the local matrix */
97: MatMultTranspose(is->A,is->x,is->y);
99: /* scatter product back into global vector */
100: VecSet(y,0);
101: VecScatterBegin(is->ctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
102: VecScatterEnd(is->ctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
103: return(0);
104: }
108: PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
109: {
110: Mat_IS *is = (Mat_IS*)A->data;
114: /* scatter the global vectors v1 and v2 into the local work vectors */
115: VecScatterBegin(is->ctx,v1,is->x,INSERT_VALUES,SCATTER_FORWARD);
116: VecScatterEnd (is->ctx,v1,is->x,INSERT_VALUES,SCATTER_FORWARD);
117: VecScatterBegin(is->ctx,v2,is->y,INSERT_VALUES,SCATTER_FORWARD);
118: VecScatterEnd (is->ctx,v2,is->y,INSERT_VALUES,SCATTER_FORWARD);
120: /* multiply the local matrix */
121: MatMultTransposeAdd(is->A,is->x,is->y,is->y);
123: /* scatter result back into global vector */
124: VecSet(v3,0);
125: VecScatterBegin(is->ctx,is->y,v3,ADD_VALUES,SCATTER_REVERSE);
126: VecScatterEnd (is->ctx,is->y,v3,ADD_VALUES,SCATTER_REVERSE);
127: return(0);
128: }
132: PetscErrorCode MatView_IS(Mat A,PetscViewer viewer)
133: {
134: Mat_IS *a = (Mat_IS*)A->data;
136: PetscViewer sviewer;
139: PetscViewerGetSingleton(viewer,&sviewer);
140: MatView(a->A,sviewer);
141: PetscViewerRestoreSingleton(viewer,&sviewer);
142: return(0);
143: }
147: PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping)
148: {
150: PetscInt n,bs;
151: Mat_IS *is = (Mat_IS*)A->data;
152: IS from,to;
153: Vec global;
156: if (is->mapping) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Mapping already set for matrix");
158: if (rmapping != cmapping) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_INCOMP,"MATIS requires the row and column mappings to be identical");
159: PetscObjectReference((PetscObject)rmapping);
160: ISLocalToGlobalMappingDestroy(&is->mapping);
161: is->mapping = rmapping;
163: /* Create the local matrix A */
164: ISLocalToGlobalMappingGetSize(rmapping,&n);
165: MatGetBlockSize(A,&bs);
166: MatCreate(PETSC_COMM_SELF,&is->A);
167: MatSetSizes(is->A,n,n,n,n);
168: MatSetBlockSize(is->A,bs);
169: MatSetOptionsPrefix(is->A,((PetscObject)A)->prefix);
170: MatAppendOptionsPrefix(is->A,"is_");
171: MatSetFromOptions(is->A);
173: /* Create the local work vectors */
174: VecCreate(PETSC_COMM_SELF,&is->x);
175: VecSetBlockSize(is->x,bs);
176: VecSetSizes(is->x,n,n);
177: VecSetOptionsPrefix(is->x,((PetscObject)A)->prefix);
178: VecAppendOptionsPrefix(is->x,"is_");
179: VecSetFromOptions(is->x);
180: VecDuplicate(is->x,&is->y);
182: /* setup the global to local scatter */
183: ISCreateStride(PETSC_COMM_SELF,n,0,1,&to);
184: ISLocalToGlobalMappingApplyIS(rmapping,to,&from);
185: MatGetVecs(A,&global,PETSC_NULL);
186: VecScatterCreate(global,from,is->x,to,&is->ctx);
187: VecDestroy(&global);
188: ISDestroy(&to);
189: ISDestroy(&from);
190: return(0);
191: }
193: #define ISG2LMapApply(mapping,n,in,out) 0;\
194: if (!(mapping)->globals) {\
195: PetscErrorCode _ISGlobalToLocalMappingApply((mapping),IS_GTOLM_MASK,0,0,0,0);CHKERRQ(_ierr);\
196: }\
197: {\
198: PetscInt _i,*_globals = (mapping)->globals,_start = (mapping)->globalstart,_end = (mapping)->globalend;\
199: for (_i=0; _i<n; _i++) {\
200: if (in[_i] < 0) out[_i] = in[_i];\
201: else if (in[_i] < _start) out[_i] = -1;\
202: else if (in[_i] > _end) out[_i] = -1;\
203: else out[_i] = _globals[in[_i] - _start];\
204: }\
205: }
209: PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv)
210: {
211: Mat_IS *is = (Mat_IS*)mat->data;
212: PetscInt rows_l[2048],cols_l[2048];
216: #if defined(PETSC_USE_DEBUG)
217: if (m > 2048 || n > 2048) {
218: SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of row/column indices must be <= 2048: they are %D %D",m,n);
219: }
220: #endif
221: ISG2LMapApply(is->mapping,m,rows,rows_l);
222: ISG2LMapApply(is->mapping,n,cols,cols_l);
223: MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);
224: return(0);
225: }
227: #undef ISG2LMapSetUp
228: #undef ISG2LMapApply
232: PetscErrorCode MatSetValuesLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
233: {
235: Mat_IS *is = (Mat_IS*)A->data;
238: MatSetValues(is->A,m,rows,n,cols,values,addv);
239: return(0);
240: }
244: PetscErrorCode MatSetValuesBlockedLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
245: {
247: Mat_IS *is = (Mat_IS*)A->data;
250: MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv);
251: return(0);
252: }
256: PetscErrorCode MatZeroRows_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
257: {
258: Mat_IS *is = (Mat_IS*)A->data;
259: PetscInt n_l=0, *rows_l = PETSC_NULL;
263: if (x && b) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_SUP,"No support");
264: if (n) {
265: PetscMalloc(n*sizeof(PetscInt),&rows_l);
266: ISGlobalToLocalMappingApply(is->mapping,IS_GTOLM_DROP,n,rows,&n_l,rows_l);
267: }
268: MatZeroRowsLocal(A,n_l,rows_l,diag,x,b);
269: PetscFree(rows_l);
270: return(0);
271: }
275: PetscErrorCode MatZeroRowsLocal_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
276: {
277: Mat_IS *is = (Mat_IS*)A->data;
279: PetscInt i;
280: PetscScalar *array;
283: if (x && b) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_SUP,"No support");
284: {
285: /*
286: Set up is->x as a "counting vector". This is in order to MatMult_IS
287: work properly in the interface nodes.
288: */
289: Vec counter;
290: PetscScalar one=1.0, zero=0.0;
291: MatGetVecs(A,&counter,PETSC_NULL);
292: VecSet(counter,zero);
293: VecSet(is->x,one);
294: VecScatterBegin(is->ctx,is->x,counter,ADD_VALUES,SCATTER_REVERSE);
295: VecScatterEnd (is->ctx,is->x,counter,ADD_VALUES,SCATTER_REVERSE);
296: VecScatterBegin(is->ctx,counter,is->x,INSERT_VALUES,SCATTER_FORWARD);
297: VecScatterEnd (is->ctx,counter,is->x,INSERT_VALUES,SCATTER_FORWARD);
298: VecDestroy(&counter);
299: }
300: if (!n) {
301: is->pure_neumann = PETSC_TRUE;
302: } else {
303: is->pure_neumann = PETSC_FALSE;
304: VecGetArray(is->x,&array);
305: MatZeroRows(is->A,n,rows,diag,0,0);
306: for (i=0; i<n; i++) {
307: MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);
308: }
309: MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);
310: MatAssemblyEnd (is->A,MAT_FINAL_ASSEMBLY);
311: VecRestoreArray(is->x,&array);
312: }
314: return(0);
315: }
319: PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type)
320: {
321: Mat_IS *is = (Mat_IS*)A->data;
325: MatAssemblyBegin(is->A,type);
326: return(0);
327: }
331: PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type)
332: {
333: Mat_IS *is = (Mat_IS*)A->data;
337: MatAssemblyEnd(is->A,type);
338: return(0);
339: }
341: EXTERN_C_BEGIN
344: PetscErrorCode MatISGetLocalMat_IS(Mat mat,Mat *local)
345: {
346: Mat_IS *is = (Mat_IS *)mat->data;
347:
349: *local = is->A;
350: return(0);
351: }
352: EXTERN_C_END
356: /*@
357: MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix.
359: Input Parameter:
360: . mat - the matrix
362: Output Parameter:
363: . local - the local matrix usually MATSEQAIJ
365: Level: advanced
367: Notes:
368: This can be called if you have precomputed the nonzero structure of the
369: matrix and want to provide it to the inner matrix object to improve the performance
370: of the MatSetValues() operation.
372: .seealso: MATIS
373: @*/
374: PetscErrorCode MatISGetLocalMat(Mat mat,Mat *local)
375: {
381: PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat *),(mat,local));
382: return(0);
383: }
385: EXTERN_C_BEGIN
388: PetscErrorCode MatISSetLocalMat_IS(Mat mat,Mat local)
389: {
390: Mat_IS *is = (Mat_IS *)mat->data;
391: PetscInt nrows,ncols,orows,ocols;
395: if(is->A) {
396: MatGetSize(is->A,&orows,&ocols);
397: MatGetSize(local,&nrows,&ncols);
398: if(orows != nrows || ocols != ncols ) {
399: 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);
400: }
401: }
402: PetscObjectReference((PetscObject)local);
403: MatDestroy(&is->A);
404: is->A = local;
406: return(0);
407: }
408: EXTERN_C_END
412: /*@
413: MatISSetLocalMat - Set the local matrix stored inside a MATIS matrix.
415: Input Parameter:
416: . mat - the matrix
417: . local - the local matrix usually MATSEQAIJ
419: Output Parameter:
421: Level: advanced
423: Notes:
424: This can be called if you have precomputed the local matrix and
425: want to provide it to the matrix object MATIS.
427: .seealso: MATIS
428: @*/
429: PetscErrorCode MatISSetLocalMat(Mat mat,Mat local)
430: {
435: PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local));
436: return(0);
437: }
441: PetscErrorCode MatZeroEntries_IS(Mat A)
442: {
443: Mat_IS *a = (Mat_IS*)A->data;
447: MatZeroEntries(a->A);
448: return(0);
449: }
453: PetscErrorCode MatScale_IS(Mat A,PetscScalar a)
454: {
455: Mat_IS *is = (Mat_IS*)A->data;
459: MatScale(is->A,a);
460: return(0);
461: }
465: PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v)
466: {
467: Mat_IS *is = (Mat_IS*)A->data;
471: /* get diagonal of the local matrix */
472: MatGetDiagonal(is->A,is->x);
474: /* scatter diagonal back into global vector */
475: VecSet(v,0);
476: VecScatterBegin(is->ctx,is->x,v,ADD_VALUES,SCATTER_REVERSE);
477: VecScatterEnd (is->ctx,is->x,v,ADD_VALUES,SCATTER_REVERSE);
478: return(0);
479: }
483: PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg)
484: {
485: Mat_IS *a = (Mat_IS*)A->data;
489: MatSetOption(a->A,op,flg);
490: return(0);
491: }
495: /*@
496: MatCreateIS - Creates a "process" unassmembled matrix, it is assembled on each
497: process but not across processes.
499: Input Parameters:
500: + comm - MPI communicator that will share the matrix
501: . bs - local and global block size of the matrix
502: . m,n,M,N - local and/or global sizes of the the left and right vector used in matrix vector products
503: - map - mapping that defines the global number for each local number
505: Output Parameter:
506: . A - the resulting matrix
508: Level: advanced
510: Notes: See MATIS for more details
511: m and n are NOT related to the size of the map, they are the size of the part of the vector owned
512: by that process. m + nghosts (or n + nghosts) is the length of map since map maps all local points
513: plus the ghost points to global indices.
515: .seealso: MATIS, MatSetLocalToGlobalMapping()
516: @*/
517: PetscErrorCode MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping map,Mat *A)
518: {
522: MatCreate(comm,A);
523: MatSetBlockSize(*A,bs);
524: MatSetSizes(*A,m,n,M,N);
525: MatSetType(*A,MATIS);
526: MatSetUp(*A);
527: MatSetLocalToGlobalMapping(*A,map,map);
528: return(0);
529: }
531: /*MC
532: MATIS - MATIS = "is" - A matrix type to be used for using the Neumann-Neumann type preconditioners.
533: This stores the matrices in globally unassembled form. Each processor
534: assembles only its local Neumann problem and the parallel matrix vector
535: product is handled "implicitly".
537: Operations Provided:
538: + MatMult()
539: . MatMultAdd()
540: . MatMultTranspose()
541: . MatMultTransposeAdd()
542: . MatZeroEntries()
543: . MatSetOption()
544: . MatZeroRows()
545: . MatZeroRowsLocal()
546: . MatSetValues()
547: . MatSetValuesLocal()
548: . MatScale()
549: . MatGetDiagonal()
550: - MatSetLocalToGlobalMapping()
552: Options Database Keys:
553: . -mat_type is - sets the matrix type to "is" during a call to MatSetFromOptions()
555: Notes: Options prefix for the inner matrix are given by -is_mat_xxx
556:
557: You must call MatSetLocalToGlobalMapping() before using this matrix type.
559: You can do matrix preallocation on the local matrix after you obtain it with
560: MatISGetLocalMat()
562: Level: advanced
564: .seealso: PC, MatISGetLocalMat(), MatSetLocalToGlobalMapping()
566: M*/
568: EXTERN_C_BEGIN
571: PetscErrorCode MatCreate_IS(Mat A)
572: {
574: Mat_IS *b;
577: PetscNewLog(A,Mat_IS,&b);
578: A->data = (void*)b;
579: PetscMemzero(A->ops,sizeof(struct _MatOps));
581: A->ops->mult = MatMult_IS;
582: A->ops->multadd = MatMultAdd_IS;
583: A->ops->multtranspose = MatMultTranspose_IS;
584: A->ops->multtransposeadd = MatMultTransposeAdd_IS;
585: A->ops->destroy = MatDestroy_IS;
586: A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS;
587: A->ops->setvalues = MatSetValues_IS;
588: A->ops->setvalueslocal = MatSetValuesLocal_IS;
589: A->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_IS;
590: A->ops->zerorows = MatZeroRows_IS;
591: A->ops->zerorowslocal = MatZeroRowsLocal_IS;
592: A->ops->assemblybegin = MatAssemblyBegin_IS;
593: A->ops->assemblyend = MatAssemblyEnd_IS;
594: A->ops->view = MatView_IS;
595: A->ops->zeroentries = MatZeroEntries_IS;
596: A->ops->scale = MatScale_IS;
597: A->ops->getdiagonal = MatGetDiagonal_IS;
598: A->ops->setoption = MatSetOption_IS;
600: PetscLayoutSetUp(A->rmap);
601: PetscLayoutSetUp(A->cmap);
603: b->A = 0;
604: b->ctx = 0;
605: b->x = 0;
606: b->y = 0;
607: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatISGetLocalMat_C","MatISGetLocalMat_IS",MatISGetLocalMat_IS);
608: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatISSetLocalMat_C","MatISSetLocalMat_IS",MatISSetLocalMat_IS);
609: PetscObjectChangeTypeName((PetscObject)A,MATIS);
611: return(0);
612: }
613: EXTERN_C_END