Actual source code: matis.c
petsc-3.5.4 2015-05-23
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 MatDuplicate_IS(Mat mat,MatDuplicateOption op,Mat *newmat)
20: {
22: Mat_IS *matis = (Mat_IS*)(mat->data);
23: PetscInt bs,m,n,M,N;
24: Mat B,localmat;
27: MatGetBlockSize(mat,&bs);
28: MatGetSize(mat,&M,&N);
29: MatGetLocalSize(mat,&m,&n);
30: MatCreateIS(PetscObjectComm((PetscObject)mat),bs,m,n,M,N,matis->mapping,&B);
31: MatDuplicate(matis->A,op,&localmat);
32: MatISSetLocalMat(B,localmat);
33: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
34: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
35: *newmat = B;
36: return(0);
37: }
41: PetscErrorCode MatIsHermitian_IS(Mat A,PetscReal tol,PetscBool *flg)
42: {
44: Mat_IS *matis = (Mat_IS*)A->data;
45: PetscBool local_sym;
48: MatIsHermitian(matis->A,tol,&local_sym);
49: MPI_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
50: return(0);
51: }
55: PetscErrorCode MatIsSymmetric_IS(Mat A,PetscReal tol,PetscBool *flg)
56: {
58: Mat_IS *matis = (Mat_IS*)A->data;
59: PetscBool local_sym;
62: MatIsSymmetric(matis->A,tol,&local_sym);
63: MPI_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
64: return(0);
65: }
69: PetscErrorCode MatDestroy_IS(Mat A)
70: {
72: Mat_IS *b = (Mat_IS*)A->data;
75: MatDestroy(&b->A);
76: VecScatterDestroy(&b->ctx);
77: VecDestroy(&b->x);
78: VecDestroy(&b->y);
79: ISLocalToGlobalMappingDestroy(&b->mapping);
80: PetscFree(A->data);
81: PetscObjectChangeTypeName((PetscObject)A,0);
82: PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL);
83: return(0);
84: }
88: PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y)
89: {
91: Mat_IS *is = (Mat_IS*)A->data;
92: PetscScalar zero = 0.0;
95: /* scatter the global vector x into the local work vector */
96: VecScatterBegin(is->ctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
97: VecScatterEnd(is->ctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
99: /* multiply the local matrix */
100: MatMult(is->A,is->x,is->y);
102: /* scatter product back into global memory */
103: VecSet(y,zero);
104: VecScatterBegin(is->ctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
105: VecScatterEnd(is->ctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
106: return(0);
107: }
111: static PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
112: {
113: Vec temp_vec;
117: if (v3 != v2) {
118: MatMult(A,v1,v3);
119: VecAXPY(v3,1.0,v2);
120: } else {
121: VecDuplicate(v2,&temp_vec);
122: MatMult(A,v1,temp_vec);
123: VecAXPY(temp_vec,1.0,v2);
124: VecCopy(temp_vec,v3);
125: VecDestroy(&temp_vec);
126: }
127: return(0);
128: }
132: PetscErrorCode MatMultTranspose_IS(Mat A,Vec x,Vec y)
133: {
134: Mat_IS *is = (Mat_IS*)A->data;
138: /* scatter the global vector x into the local work vector */
139: VecScatterBegin(is->ctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
140: VecScatterEnd(is->ctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
142: /* multiply the local matrix */
143: MatMultTranspose(is->A,is->x,is->y);
145: /* scatter product back into global vector */
146: VecSet(y,0);
147: VecScatterBegin(is->ctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
148: VecScatterEnd(is->ctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
149: return(0);
150: }
154: PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
155: {
156: Vec temp_vec;
160: if (v3 != v2) {
161: MatMultTranspose(A,v1,v3);
162: VecAXPY(v3,1.0,v2);
163: } else {
164: VecDuplicate(v2,&temp_vec);
165: MatMultTranspose(A,v1,temp_vec);
166: VecAXPY(temp_vec,1.0,v2);
167: VecCopy(temp_vec,v3);
168: VecDestroy(&temp_vec);
169: }
170: return(0);
171: }
175: PetscErrorCode MatView_IS(Mat A,PetscViewer viewer)
176: {
177: Mat_IS *a = (Mat_IS*)A->data;
179: PetscViewer sviewer;
182: PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
183: PetscViewerGetSingleton(viewer,&sviewer);
184: MatView(a->A,sviewer);
185: PetscViewerRestoreSingleton(viewer,&sviewer);
186: return(0);
187: }
191: PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping)
192: {
194: PetscInt n,bs;
195: Mat_IS *is = (Mat_IS*)A->data;
196: IS from,to;
197: Vec global;
201: if (rmapping != cmapping) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"MATIS requires the row and column mappings to be identical");
202: if (is->mapping) { /* Currenly destroys the objects that will be created by this routine. Is there anything else that should be checked? */
203: ISLocalToGlobalMappingDestroy(&is->mapping);
204: VecDestroy(&is->x);
205: VecDestroy(&is->y);
206: VecScatterDestroy(&is->ctx);
207: MatDestroy(&is->A);
208: }
209: PetscObjectReference((PetscObject)rmapping);
210: ISLocalToGlobalMappingDestroy(&is->mapping);
211: is->mapping = rmapping;
212: /*
213: PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping);
214: PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping);
215: */
217: /* Create the local matrix A */
218: ISLocalToGlobalMappingGetSize(rmapping,&n);
219: MatGetBlockSize(A,&bs);
220: MatCreate(PETSC_COMM_SELF,&is->A);
221: MatSetSizes(is->A,n,n,n,n);
222: MatSetBlockSize(is->A,bs);
223: MatSetOptionsPrefix(is->A,((PetscObject)A)->prefix);
224: MatAppendOptionsPrefix(is->A,"is_");
225: MatSetFromOptions(is->A);
227: /* Create the local work vectors */
228: VecCreate(PETSC_COMM_SELF,&is->x);
229: VecSetBlockSize(is->x,bs);
230: VecSetSizes(is->x,n,n);
231: VecSetOptionsPrefix(is->x,((PetscObject)A)->prefix);
232: VecAppendOptionsPrefix(is->x,"is_");
233: VecSetFromOptions(is->x);
234: VecDuplicate(is->x,&is->y);
236: /* setup the global to local scatter */
237: ISCreateStride(PETSC_COMM_SELF,n,0,1,&to);
238: ISLocalToGlobalMappingApplyIS(rmapping,to,&from);
239: MatGetVecs(A,&global,NULL);
240: VecScatterCreate(global,from,is->x,to,&is->ctx);
241: VecDestroy(&global);
242: ISDestroy(&to);
243: ISDestroy(&from);
244: return(0);
245: }
249: PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv)
250: {
251: Mat_IS *is = (Mat_IS*)mat->data;
252: PetscInt rows_l[2048],cols_l[2048];
256: #if defined(PETSC_USE_DEBUG)
257: 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);
258: #endif
259: ISG2LMapApply(is->mapping,m,rows,rows_l);
260: ISG2LMapApply(is->mapping,n,cols,cols_l);
261: MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);
262: return(0);
263: }
265: #undef ISG2LMapSetUp
266: #undef ISG2LMapApply
270: PetscErrorCode MatSetValuesLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
271: {
273: Mat_IS *is = (Mat_IS*)A->data;
276: MatSetValues(is->A,m,rows,n,cols,values,addv);
277: return(0);
278: }
282: PetscErrorCode MatSetValuesBlockedLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
283: {
285: Mat_IS *is = (Mat_IS*)A->data;
288: MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv);
289: return(0);
290: }
294: PetscErrorCode MatZeroRows_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
295: {
296: Mat_IS *is = (Mat_IS*)A->data;
297: PetscInt n_l = 0, *rows_l = NULL;
301: if (x && b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support");
302: if (n) {
303: PetscMalloc1(n,&rows_l);
304: ISGlobalToLocalMappingApply(is->mapping,IS_GTOLM_DROP,n,rows,&n_l,rows_l);
305: }
306: MatZeroRowsLocal(A,n_l,rows_l,diag,x,b);
307: PetscFree(rows_l);
308: return(0);
309: }
313: PetscErrorCode MatZeroRowsLocal_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
314: {
315: Mat_IS *is = (Mat_IS*)A->data;
317: PetscInt i;
318: PetscScalar *array;
321: if (x && b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support");
322: {
323: /*
324: Set up is->x as a "counting vector". This is in order to MatMult_IS
325: work properly in the interface nodes.
326: */
327: Vec counter;
328: PetscScalar one=1.0, zero=0.0;
329: MatGetVecs(A,&counter,NULL);
330: VecSet(counter,zero);
331: VecSet(is->x,one);
332: VecScatterBegin(is->ctx,is->x,counter,ADD_VALUES,SCATTER_REVERSE);
333: VecScatterEnd (is->ctx,is->x,counter,ADD_VALUES,SCATTER_REVERSE);
334: VecScatterBegin(is->ctx,counter,is->x,INSERT_VALUES,SCATTER_FORWARD);
335: VecScatterEnd (is->ctx,counter,is->x,INSERT_VALUES,SCATTER_FORWARD);
336: VecDestroy(&counter);
337: }
338: if (!n) {
339: is->pure_neumann = PETSC_TRUE;
340: } else {
341: is->pure_neumann = PETSC_FALSE;
343: VecGetArray(is->x,&array);
344: MatZeroRows(is->A,n,rows,diag,0,0);
345: for (i=0; i<n; i++) {
346: MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);
347: }
348: MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);
349: MatAssemblyEnd (is->A,MAT_FINAL_ASSEMBLY);
350: VecRestoreArray(is->x,&array);
351: }
352: return(0);
353: }
357: PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type)
358: {
359: Mat_IS *is = (Mat_IS*)A->data;
363: MatAssemblyBegin(is->A,type);
364: return(0);
365: }
369: PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type)
370: {
371: Mat_IS *is = (Mat_IS*)A->data;
375: MatAssemblyEnd(is->A,type);
376: return(0);
377: }
381: PetscErrorCode MatISGetLocalMat_IS(Mat mat,Mat *local)
382: {
383: Mat_IS *is = (Mat_IS*)mat->data;
386: *local = is->A;
387: return(0);
388: }
392: /*@
393: MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix.
395: Input Parameter:
396: . mat - the matrix
398: Output Parameter:
399: . local - the local matrix usually MATSEQAIJ
401: Level: advanced
403: Notes:
404: This can be called if you have precomputed the nonzero structure of the
405: matrix and want to provide it to the inner matrix object to improve the performance
406: of the MatSetValues() operation.
408: .seealso: MATIS
409: @*/
410: PetscErrorCode MatISGetLocalMat(Mat mat,Mat *local)
411: {
417: PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat*),(mat,local));
418: return(0);
419: }
423: PetscErrorCode MatISSetLocalMat_IS(Mat mat,Mat local)
424: {
425: Mat_IS *is = (Mat_IS*)mat->data;
426: PetscInt nrows,ncols,orows,ocols;
430: if (is->A) {
431: MatGetSize(is->A,&orows,&ocols);
432: MatGetSize(local,&nrows,&ncols);
433: 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);
434: }
435: PetscObjectReference((PetscObject)local);
436: MatDestroy(&is->A);
437: is->A = local;
438: return(0);
439: }
443: /*@
444: MatISSetLocalMat - Set the local matrix stored inside a MATIS matrix.
446: Input Parameter:
447: . mat - the matrix
448: . local - the local matrix usually MATSEQAIJ
450: Output Parameter:
452: Level: advanced
454: Notes:
455: This can be called if you have precomputed the local matrix and
456: want to provide it to the matrix object MATIS.
458: .seealso: MATIS
459: @*/
460: PetscErrorCode MatISSetLocalMat(Mat mat,Mat local)
461: {
466: PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local));
467: return(0);
468: }
472: PetscErrorCode MatZeroEntries_IS(Mat A)
473: {
474: Mat_IS *a = (Mat_IS*)A->data;
478: MatZeroEntries(a->A);
479: return(0);
480: }
484: PetscErrorCode MatScale_IS(Mat A,PetscScalar a)
485: {
486: Mat_IS *is = (Mat_IS*)A->data;
490: MatScale(is->A,a);
491: return(0);
492: }
496: PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v)
497: {
498: Mat_IS *is = (Mat_IS*)A->data;
502: /* get diagonal of the local matrix */
503: MatGetDiagonal(is->A,is->x);
505: /* scatter diagonal back into global vector */
506: VecSet(v,0);
507: VecScatterBegin(is->ctx,is->x,v,ADD_VALUES,SCATTER_REVERSE);
508: VecScatterEnd (is->ctx,is->x,v,ADD_VALUES,SCATTER_REVERSE);
509: return(0);
510: }
514: PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg)
515: {
516: Mat_IS *a = (Mat_IS*)A->data;
520: MatSetOption(a->A,op,flg);
521: return(0);
522: }
526: /*@
527: MatCreateIS - Creates a "process" unassmembled matrix, it is assembled on each
528: process but not across processes.
530: Input Parameters:
531: + comm - MPI communicator that will share the matrix
532: . bs - local and global block size of the matrix
533: . m,n,M,N - local and/or global sizes of the left and right vector used in matrix vector products
534: - map - mapping that defines the global number for each local number
536: Output Parameter:
537: . A - the resulting matrix
539: Level: advanced
541: Notes: See MATIS for more details
542: m and n are NOT related to the size of the map, they are the size of the part of the vector owned
543: by that process. m + nghosts (or n + nghosts) is the length of map since map maps all local points
544: plus the ghost points to global indices.
546: .seealso: MATIS, MatSetLocalToGlobalMapping()
547: @*/
548: PetscErrorCode MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping map,Mat *A)
549: {
553: MatCreate(comm,A);
554: MatSetBlockSize(*A,bs);
555: MatSetSizes(*A,m,n,M,N);
556: MatSetType(*A,MATIS);
557: MatSetUp(*A);
558: MatSetLocalToGlobalMapping(*A,map,map);
559: return(0);
560: }
562: /*MC
563: MATIS - MATIS = "is" - A matrix type to be used for using the Neumann-Neumann type preconditioners, see PCBDDC.
564: This stores the matrices in globally unassembled form. Each processor
565: assembles only its local Neumann problem and the parallel matrix vector
566: product is handled "implicitly".
568: Operations Provided:
569: + MatMult()
570: . MatMultAdd()
571: . MatMultTranspose()
572: . MatMultTransposeAdd()
573: . MatZeroEntries()
574: . MatSetOption()
575: . MatZeroRows()
576: . MatZeroRowsLocal()
577: . MatSetValues()
578: . MatSetValuesLocal()
579: . MatScale()
580: . MatGetDiagonal()
581: - MatSetLocalToGlobalMapping()
583: Options Database Keys:
584: . -mat_type is - sets the matrix type to "is" during a call to MatSetFromOptions()
586: Notes: Options prefix for the inner matrix are given by -is_mat_xxx
588: You must call MatSetLocalToGlobalMapping() before using this matrix type.
590: You can do matrix preallocation on the local matrix after you obtain it with
591: MatISGetLocalMat()
593: Level: advanced
595: .seealso: PC, MatISGetLocalMat(), MatSetLocalToGlobalMapping(), PCBDDC
597: M*/
601: PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A)
602: {
604: Mat_IS *b;
607: PetscNewLog(A,&b);
608: A->data = (void*)b;
609: PetscMemzero(A->ops,sizeof(struct _MatOps));
611: A->ops->mult = MatMult_IS;
612: A->ops->multadd = MatMultAdd_IS;
613: A->ops->multtranspose = MatMultTranspose_IS;
614: A->ops->multtransposeadd = MatMultTransposeAdd_IS;
615: A->ops->destroy = MatDestroy_IS;
616: A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS;
617: A->ops->setvalues = MatSetValues_IS;
618: A->ops->setvalueslocal = MatSetValuesLocal_IS;
619: A->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_IS;
620: A->ops->zerorows = MatZeroRows_IS;
621: A->ops->zerorowslocal = MatZeroRowsLocal_IS;
622: A->ops->assemblybegin = MatAssemblyBegin_IS;
623: A->ops->assemblyend = MatAssemblyEnd_IS;
624: A->ops->view = MatView_IS;
625: A->ops->zeroentries = MatZeroEntries_IS;
626: A->ops->scale = MatScale_IS;
627: A->ops->getdiagonal = MatGetDiagonal_IS;
628: A->ops->setoption = MatSetOption_IS;
629: A->ops->ishermitian = MatIsHermitian_IS;
630: A->ops->issymmetric = MatIsSymmetric_IS;
631: A->ops->duplicate = MatDuplicate_IS;
633: PetscLayoutSetUp(A->rmap);
634: PetscLayoutSetUp(A->cmap);
636: b->A = 0;
637: b->ctx = 0;
638: b->x = 0;
639: b->y = 0;
640: b->mapping = 0;
641: PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS);
642: PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS);
643: PetscObjectChangeTypeName((PetscObject)A,MATIS);
644: return(0);
645: }