Actual source code: matis.c
petsc-3.4.5 2014-06-29
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: }