Actual source code: matis.c
petsc-3.7.5 2017-01-01
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*/
17: static PetscErrorCode MatISComputeSF_Private(Mat);
21: static PetscErrorCode MatISComputeSF_Private(Mat B)
22: {
23: Mat_IS *matis = (Mat_IS*)(B->data);
24: const PetscInt *gidxs;
28: MatGetSize(matis->A,&matis->sf_nleaves,NULL);
29: MatGetLocalSize(B,&matis->sf_nroots,NULL);
30: PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->sf);
31: ISLocalToGlobalMappingGetIndices(B->rmap->mapping,&gidxs);
32: /* PETSC_OWN_POINTER refers to ilocal which is NULL */
33: PetscSFSetGraphLayout(matis->sf,B->rmap,matis->sf_nleaves,NULL,PETSC_OWN_POINTER,gidxs);
34: ISLocalToGlobalMappingRestoreIndices(B->rmap->mapping,&gidxs);
35: PetscMalloc2(matis->sf_nroots,&matis->sf_rootdata,matis->sf_nleaves,&matis->sf_leafdata);
36: return(0);
37: }
41: /*@
42: MatISSetPreallocation - Preallocates memory for a MATIS parallel matrix.
44: Collective on MPI_Comm
46: Input Parameters:
47: + B - the matrix
48: . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix
49: (same value is used for all local rows)
50: . d_nnz - array containing the number of nonzeros in the various rows of the
51: DIAGONAL portion of the local submatrix (possibly different for each row)
52: or NULL, if d_nz is used to specify the nonzero structure.
53: The size of this array is equal to the number of local rows, i.e 'm'.
54: For matrices that will be factored, you must leave room for (and set)
55: the diagonal entry even if it is zero.
56: . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local
57: submatrix (same value is used for all local rows).
58: - o_nnz - array containing the number of nonzeros in the various rows of the
59: OFF-DIAGONAL portion of the local submatrix (possibly different for
60: each row) or NULL, if o_nz is used to specify the nonzero
61: structure. The size of this array is equal to the number
62: of local rows, i.e 'm'.
64: If the *_nnz parameter is given then the *_nz parameter is ignored
66: Level: intermediate
68: Notes: This function has the same interface as the MPIAIJ preallocation routine in order to simplify the transition
69: from the asssembled format to the unassembled one. It overestimates the preallocation of MATIS local
70: matrices; for exact preallocation, the user should set the preallocation directly on local matrix objects.
72: .keywords: matrix
74: .seealso: MatCreate(), MatCreateIS(), MatMPIAIJSetPreallocation(), MatISGetLocalMat()
75: @*/
76: PetscErrorCode MatISSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
77: {
83: PetscTryMethod(B,"MatISSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz));
84: return(0);
85: }
89: PetscErrorCode MatISSetPreallocation_IS(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
90: {
91: Mat_IS *matis = (Mat_IS*)(B->data);
92: PetscInt bs,i,nlocalcols;
96: if (!matis->A) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"You should first call MatSetLocalToGlobalMapping");
97: if (!matis->sf) { /* setup SF if not yet created and allocate rootdata and leafdata */
98: MatISComputeSF_Private(B);
99: }
100: if (!d_nnz) {
101: for (i=0;i<matis->sf_nroots;i++) matis->sf_rootdata[i] = d_nz;
102: } else {
103: for (i=0;i<matis->sf_nroots;i++) matis->sf_rootdata[i] = d_nnz[i];
104: }
105: if (!o_nnz) {
106: for (i=0;i<matis->sf_nroots;i++) matis->sf_rootdata[i] += o_nz;
107: } else {
108: for (i=0;i<matis->sf_nroots;i++) matis->sf_rootdata[i] += o_nnz[i];
109: }
110: PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
111: MatGetSize(matis->A,NULL,&nlocalcols);
112: MatGetBlockSize(matis->A,&bs);
113: PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);
114: for (i=0;i<matis->sf_nleaves;i++) {
115: matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i],nlocalcols);
116: }
117: MatSeqAIJSetPreallocation(matis->A,0,matis->sf_leafdata);
118: for (i=0;i<matis->sf_nleaves/bs;i++) {
119: matis->sf_leafdata[i] = matis->sf_leafdata[i*bs]/bs;
120: }
121: MatSeqBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata);
122: for (i=0;i<matis->sf_nleaves/bs;i++) {
123: matis->sf_leafdata[i] = matis->sf_leafdata[i]-i;
124: }
125: MatSeqSBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata);
126: return(0);
127: }
131: PETSC_EXTERN PetscErrorCode MatISSetMPIXAIJPreallocation_Private(Mat A, Mat B, PetscBool maxreduce)
132: {
133: Mat_IS *matis = (Mat_IS*)(A->data);
134: PetscInt *my_dnz,*my_onz,*dnz,*onz,*mat_ranges,*row_ownership;
135: const PetscInt *global_indices_r,*global_indices_c;
136: PetscInt i,j,bs,rows,cols;
137: PetscInt lrows,lcols;
138: PetscInt local_rows,local_cols;
139: PetscMPIInt nsubdomains;
140: PetscBool isdense,issbaij;
141: PetscErrorCode ierr;
144: MPI_Comm_size(PetscObjectComm((PetscObject)A),&nsubdomains);
145: MatGetSize(A,&rows,&cols);
146: MatGetBlockSize(A,&bs);
147: MatGetSize(matis->A,&local_rows,&local_cols);
148: PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isdense);
149: PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);
150: ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&global_indices_r);
151: if (A->rmap->mapping != A->cmap->mapping) {
152: ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&global_indices_c);
153: } else {
154: global_indices_c = global_indices_r;
155: }
157: if (issbaij) {
158: MatGetRowUpperTriangular(matis->A);
159: }
160: /*
161: An SF reduce is needed to sum up properly on shared rows.
162: Note that generally preallocation is not exact, since it overestimates nonzeros
163: */
164: if (!matis->sf) { /* setup SF if not yet created and allocate rootdata and leafdata */
165: MatISComputeSF_Private(A);
166: }
167: MatGetLocalSize(A,&lrows,&lcols);
168: MatPreallocateInitialize(PetscObjectComm((PetscObject)A),lrows,lcols,dnz,onz);
169: /* All processes need to compute entire row ownership */
170: PetscMalloc1(rows,&row_ownership);
171: MatGetOwnershipRanges(A,(const PetscInt**)&mat_ranges);
172: for (i=0;i<nsubdomains;i++) {
173: for (j=mat_ranges[i];j<mat_ranges[i+1];j++) {
174: row_ownership[j] = i;
175: }
176: }
178: /*
179: my_dnz and my_onz contains exact contribution to preallocation from each local mat
180: then, they will be summed up properly. This way, preallocation is always sufficient
181: */
182: PetscCalloc2(local_rows,&my_dnz,local_rows,&my_onz);
183: /* preallocation as a MATAIJ */
184: if (isdense) { /* special case for dense local matrices */
185: for (i=0;i<local_rows;i++) {
186: PetscInt index_row = global_indices_r[i];
187: for (j=i;j<local_rows;j++) {
188: PetscInt owner = row_ownership[index_row];
189: PetscInt index_col = global_indices_c[j];
190: if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */
191: my_dnz[i] += 1;
192: } else { /* offdiag block */
193: my_onz[i] += 1;
194: }
195: /* same as before, interchanging rows and cols */
196: if (i != j) {
197: owner = row_ownership[index_col];
198: if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) {
199: my_dnz[j] += 1;
200: } else {
201: my_onz[j] += 1;
202: }
203: }
204: }
205: }
206: } else { /* TODO: this could be optimized using MatGetRowIJ */
207: for (i=0;i<local_rows;i++) {
208: const PetscInt *cols;
209: PetscInt ncols,index_row = global_indices_r[i];
210: MatGetRow(matis->A,i,&ncols,&cols,NULL);
211: for (j=0;j<ncols;j++) {
212: PetscInt owner = row_ownership[index_row];
213: PetscInt index_col = global_indices_c[cols[j]];
214: if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */
215: my_dnz[i] += 1;
216: } else { /* offdiag block */
217: my_onz[i] += 1;
218: }
219: /* same as before, interchanging rows and cols */
220: if (issbaij && index_col != index_row) {
221: owner = row_ownership[index_col];
222: if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) {
223: my_dnz[cols[j]] += 1;
224: } else {
225: my_onz[cols[j]] += 1;
226: }
227: }
228: }
229: MatRestoreRow(matis->A,i,&ncols,&cols,NULL);
230: }
231: }
232: ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&global_indices_r);
233: if (global_indices_c != global_indices_r) {
234: ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&global_indices_c);
235: }
236: PetscFree(row_ownership);
238: /* Reduce my_dnz and my_onz */
239: if (maxreduce) {
240: PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);
241: PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);
242: PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);
243: PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);
244: } else {
245: PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);
246: PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);
247: PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);
248: PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);
249: }
250: PetscFree2(my_dnz,my_onz);
252: /* Resize preallocation if overestimated */
253: for (i=0;i<lrows;i++) {
254: dnz[i] = PetscMin(dnz[i],lcols);
255: onz[i] = PetscMin(onz[i],cols-lcols);
256: }
257: /* set preallocation */
258: MatMPIAIJSetPreallocation(B,0,dnz,0,onz);
259: for (i=0;i<lrows/bs;i++) {
260: dnz[i] = dnz[i*bs]/bs;
261: onz[i] = onz[i*bs]/bs;
262: }
263: MatMPIBAIJSetPreallocation(B,bs,0,dnz,0,onz);
264: MatMPISBAIJSetPreallocation(B,bs,0,dnz,0,onz);
265: MatPreallocateFinalize(dnz,onz);
266: if (issbaij) {
267: MatRestoreRowUpperTriangular(matis->A);
268: }
269: return(0);
270: }
274: PetscErrorCode MatISGetMPIXAIJ_IS(Mat mat, MatReuse reuse, Mat *M)
275: {
276: Mat_IS *matis = (Mat_IS*)(mat->data);
277: Mat local_mat;
278: /* info on mat */
279: PetscInt bs,rows,cols,lrows,lcols;
280: PetscInt local_rows,local_cols;
281: PetscBool isdense,issbaij,isseqaij;
282: PetscMPIInt nsubdomains;
283: /* values insertion */
284: PetscScalar *array;
285: /* work */
289: /* get info from mat */
290: MPI_Comm_size(PetscObjectComm((PetscObject)mat),&nsubdomains);
291: if (nsubdomains == 1) {
292: if (reuse == MAT_INITIAL_MATRIX) {
293: MatDuplicate(matis->A,MAT_COPY_VALUES,&(*M));
294: } else {
295: MatCopy(matis->A,*M,SAME_NONZERO_PATTERN);
296: }
297: return(0);
298: }
299: MatGetSize(mat,&rows,&cols);
300: MatGetBlockSize(mat,&bs);
301: MatGetLocalSize(mat,&lrows,&lcols);
302: MatGetSize(matis->A,&local_rows,&local_cols);
303: PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isdense);
304: PetscObjectTypeCompare((PetscObject)matis->A,MATSEQAIJ,&isseqaij);
305: PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);
307: if (reuse == MAT_INITIAL_MATRIX) {
308: MatType new_mat_type;
309: PetscBool issbaij_red;
311: /* determining new matrix type */
312: MPIU_Allreduce(&issbaij,&issbaij_red,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)mat));
313: if (issbaij_red) {
314: new_mat_type = MATSBAIJ;
315: } else {
316: if (bs>1) {
317: new_mat_type = MATBAIJ;
318: } else {
319: new_mat_type = MATAIJ;
320: }
321: }
323: MatCreate(PetscObjectComm((PetscObject)mat),M);
324: MatSetSizes(*M,lrows,lcols,rows,cols);
325: MatSetBlockSize(*M,bs);
326: MatSetType(*M,new_mat_type);
327: MatISSetMPIXAIJPreallocation_Private(mat,*M,PETSC_FALSE);
328: } else {
329: PetscInt mbs,mrows,mcols,mlrows,mlcols;
330: /* some checks */
331: MatGetBlockSize(*M,&mbs);
332: MatGetSize(*M,&mrows,&mcols);
333: MatGetLocalSize(*M,&mlrows,&mlcols);
334: if (mrows != rows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of rows (%d != %d)",rows,mrows);
335: if (mcols != cols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of cols (%d != %d)",cols,mcols);
336: if (mlrows != lrows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local rows (%d != %d)",lrows,mlrows);
337: if (mlcols != lcols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local cols (%d != %d)",lcols,mlcols);
338: if (mbs != bs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong block size (%d != %d)",bs,mbs);
339: MatZeroEntries(*M);
340: }
342: if (issbaij) {
343: MatConvert(matis->A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&local_mat);
344: } else {
345: PetscObjectReference((PetscObject)matis->A);
346: local_mat = matis->A;
347: }
349: /* Set values */
350: MatSetLocalToGlobalMapping(*M,mat->rmap->mapping,mat->cmap->mapping);
351: if (isdense) { /* special case for dense local matrices */
352: PetscInt i,*dummy_rows,*dummy_cols;
354: if (local_rows != local_cols) {
355: PetscMalloc2(local_rows,&dummy_rows,local_cols,&dummy_cols);
356: for (i=0;i<local_rows;i++) dummy_rows[i] = i;
357: for (i=0;i<local_cols;i++) dummy_cols[i] = i;
358: } else {
359: PetscMalloc1(local_rows,&dummy_rows);
360: for (i=0;i<local_rows;i++) dummy_rows[i] = i;
361: dummy_cols = dummy_rows;
362: }
363: MatSetOption(*M,MAT_ROW_ORIENTED,PETSC_FALSE);
364: MatDenseGetArray(local_mat,&array);
365: MatSetValuesLocal(*M,local_rows,dummy_rows,local_cols,dummy_cols,array,ADD_VALUES);
366: MatDenseRestoreArray(local_mat,&array);
367: if (dummy_rows != dummy_cols) {
368: PetscFree2(dummy_rows,dummy_cols);
369: } else {
370: PetscFree(dummy_rows);
371: }
372: } else if (isseqaij) {
373: PetscInt i,nvtxs,*xadj,*adjncy;
374: PetscBool done;
376: MatGetRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);
377: if (!done) SETERRQ1(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ called in %s\n",__FUNCT__);
378: MatSeqAIJGetArray(local_mat,&array);
379: for (i=0;i<nvtxs;i++) {
380: MatSetValuesLocal(*M,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],array+xadj[i],ADD_VALUES);
381: }
382: MatRestoreRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);
383: if (!done) SETERRQ1(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ called in %s\n",__FUNCT__);
384: MatSeqAIJRestoreArray(local_mat,&array);
385: } else { /* very basic values insertion for all other matrix types */
386: PetscInt i;
388: for (i=0;i<local_rows;i++) {
389: PetscInt j;
390: const PetscInt *local_indices_cols;
392: MatGetRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);
393: MatSetValuesLocal(*M,1,&i,j,local_indices_cols,array,ADD_VALUES);
394: MatRestoreRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);
395: }
396: }
397: MatAssemblyBegin(*M,MAT_FINAL_ASSEMBLY);
398: MatDestroy(&local_mat);
399: MatAssemblyEnd(*M,MAT_FINAL_ASSEMBLY);
400: if (isdense) {
401: MatSetOption(*M,MAT_ROW_ORIENTED,PETSC_TRUE);
402: }
403: return(0);
404: }
408: /*@
409: MatISGetMPIXAIJ - Converts MATIS matrix into a parallel AIJ format
411: Input Parameter:
412: . mat - the matrix (should be of type MATIS)
413: . reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
415: Output Parameter:
416: . newmat - the matrix in AIJ format
418: Level: developer
420: Notes: mat and *newmat cannot be the same object when MAT_REUSE_MATRIX is requested.
422: .seealso: MATIS
423: @*/
424: PetscErrorCode MatISGetMPIXAIJ(Mat mat, MatReuse reuse, Mat *newmat)
425: {
432: if (reuse != MAT_INITIAL_MATRIX) {
435: if (mat == *newmat) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse the same matrix");
436: }
437: PetscUseMethod(mat,"MatISGetMPIXAIJ_C",(Mat,MatReuse,Mat*),(mat,reuse,newmat));
438: return(0);
439: }
443: PetscErrorCode MatDuplicate_IS(Mat mat,MatDuplicateOption op,Mat *newmat)
444: {
446: Mat_IS *matis = (Mat_IS*)(mat->data);
447: PetscInt bs,m,n,M,N;
448: Mat B,localmat;
451: MatGetBlockSize(mat,&bs);
452: MatGetSize(mat,&M,&N);
453: MatGetLocalSize(mat,&m,&n);
454: MatCreateIS(PetscObjectComm((PetscObject)mat),bs,m,n,M,N,mat->rmap->mapping,mat->cmap->mapping,&B);
455: MatDuplicate(matis->A,op,&localmat);
456: MatISSetLocalMat(B,localmat);
457: MatDestroy(&localmat);
458: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
459: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
460: *newmat = B;
461: return(0);
462: }
466: PetscErrorCode MatIsHermitian_IS(Mat A,PetscReal tol,PetscBool *flg)
467: {
469: Mat_IS *matis = (Mat_IS*)A->data;
470: PetscBool local_sym;
473: MatIsHermitian(matis->A,tol,&local_sym);
474: MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
475: return(0);
476: }
480: PetscErrorCode MatIsSymmetric_IS(Mat A,PetscReal tol,PetscBool *flg)
481: {
483: Mat_IS *matis = (Mat_IS*)A->data;
484: PetscBool local_sym;
487: MatIsSymmetric(matis->A,tol,&local_sym);
488: MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
489: return(0);
490: }
494: PetscErrorCode MatDestroy_IS(Mat A)
495: {
497: Mat_IS *b = (Mat_IS*)A->data;
500: MatDestroy(&b->A);
501: VecScatterDestroy(&b->cctx);
502: VecScatterDestroy(&b->rctx);
503: VecDestroy(&b->x);
504: VecDestroy(&b->y);
505: PetscSFDestroy(&b->sf);
506: PetscFree2(b->sf_rootdata,b->sf_leafdata);
507: PetscFree(A->data);
508: PetscObjectChangeTypeName((PetscObject)A,0);
509: PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL);
510: PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",NULL);
511: PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",NULL);
512: PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",NULL);
513: return(0);
514: }
518: PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y)
519: {
521: Mat_IS *is = (Mat_IS*)A->data;
522: PetscScalar zero = 0.0;
525: /* scatter the global vector x into the local work vector */
526: VecScatterBegin(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
527: VecScatterEnd(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);
529: /* multiply the local matrix */
530: MatMult(is->A,is->x,is->y);
532: /* scatter product back into global memory */
533: VecSet(y,zero);
534: VecScatterBegin(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
535: VecScatterEnd(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);
536: return(0);
537: }
541: PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
542: {
543: Vec temp_vec;
547: if (v3 != v2) {
548: MatMult(A,v1,v3);
549: VecAXPY(v3,1.0,v2);
550: } else {
551: VecDuplicate(v2,&temp_vec);
552: MatMult(A,v1,temp_vec);
553: VecAXPY(temp_vec,1.0,v2);
554: VecCopy(temp_vec,v3);
555: VecDestroy(&temp_vec);
556: }
557: return(0);
558: }
562: PetscErrorCode MatMultTranspose_IS(Mat A,Vec y,Vec x)
563: {
564: Mat_IS *is = (Mat_IS*)A->data;
568: /* scatter the global vector x into the local work vector */
569: VecScatterBegin(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);
570: VecScatterEnd(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);
572: /* multiply the local matrix */
573: MatMultTranspose(is->A,is->y,is->x);
575: /* scatter product back into global vector */
576: VecSet(x,0);
577: VecScatterBegin(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);
578: VecScatterEnd(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);
579: return(0);
580: }
584: PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3)
585: {
586: Vec temp_vec;
590: if (v3 != v2) {
591: MatMultTranspose(A,v1,v3);
592: VecAXPY(v3,1.0,v2);
593: } else {
594: VecDuplicate(v2,&temp_vec);
595: MatMultTranspose(A,v1,temp_vec);
596: VecAXPY(temp_vec,1.0,v2);
597: VecCopy(temp_vec,v3);
598: VecDestroy(&temp_vec);
599: }
600: return(0);
601: }
605: PetscErrorCode MatView_IS(Mat A,PetscViewer viewer)
606: {
607: Mat_IS *a = (Mat_IS*)A->data;
609: PetscViewer sviewer;
612: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
613: MatView(a->A,sviewer);
614: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
615: return(0);
616: }
620: PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping)
621: {
623: PetscInt nr,rbs,nc,cbs;
624: Mat_IS *is = (Mat_IS*)A->data;
625: IS from,to;
626: Vec cglobal,rglobal;
631: /* Destroy any previous data */
632: VecDestroy(&is->x);
633: VecDestroy(&is->y);
634: VecScatterDestroy(&is->rctx);
635: VecScatterDestroy(&is->cctx);
636: MatDestroy(&is->A);
637: PetscSFDestroy(&is->sf);
638: PetscFree2(is->sf_rootdata,is->sf_leafdata);
640: /* Setup Layout and set local to global maps */
641: PetscLayoutSetUp(A->rmap);
642: PetscLayoutSetUp(A->cmap);
643: PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping);
644: PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping);
646: /* Create the local matrix A */
647: ISLocalToGlobalMappingGetSize(rmapping,&nr);
648: ISLocalToGlobalMappingGetBlockSize(rmapping,&rbs);
649: ISLocalToGlobalMappingGetSize(cmapping,&nc);
650: ISLocalToGlobalMappingGetBlockSize(cmapping,&cbs);
651: MatCreate(PETSC_COMM_SELF,&is->A);
652: MatSetType(is->A,MATAIJ);
653: MatSetSizes(is->A,nr,nc,nr,nc);
654: MatSetBlockSizes(is->A,rbs,cbs);
655: MatSetOptionsPrefix(is->A,((PetscObject)A)->prefix);
656: MatAppendOptionsPrefix(is->A,"is_");
657: MatSetFromOptions(is->A);
659: /* Create the local work vectors */
660: MatCreateVecs(is->A,&is->x,&is->y);
662: /* setup the global to local scatters */
663: MatCreateVecs(A,&cglobal,&rglobal);
664: ISCreateStride(PETSC_COMM_SELF,nr,0,1,&to);
665: ISLocalToGlobalMappingApplyIS(rmapping,to,&from);
666: VecScatterCreate(rglobal,from,is->y,to,&is->rctx);
667: if (rmapping != cmapping) {
668: ISDestroy(&to);
669: ISDestroy(&from);
670: ISCreateStride(PETSC_COMM_SELF,nc,0,1,&to);
671: ISLocalToGlobalMappingApplyIS(cmapping,to,&from);
672: VecScatterCreate(cglobal,from,is->x,to,&is->cctx);
673: } else {
674: PetscObjectReference((PetscObject)is->rctx);
675: is->cctx = is->rctx;
676: }
677: VecDestroy(&rglobal);
678: VecDestroy(&cglobal);
679: ISDestroy(&to);
680: ISDestroy(&from);
681: return(0);
682: }
686: PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv)
687: {
688: Mat_IS *is = (Mat_IS*)mat->data;
689: PetscInt rows_l[2048],cols_l[2048];
693: #if defined(PETSC_USE_DEBUG)
694: 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);
695: #endif
696: ISG2LMapApply(mat->rmap->mapping,m,rows,rows_l);
697: ISG2LMapApply(mat->cmap->mapping,n,cols,cols_l);
698: MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);
699: return(0);
700: }
704: PetscErrorCode MatSetValuesLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
705: {
707: Mat_IS *is = (Mat_IS*)A->data;
710: MatSetValues(is->A,m,rows,n,cols,values,addv);
711: return(0);
712: }
716: PetscErrorCode MatSetValuesBlockedLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
717: {
719: Mat_IS *is = (Mat_IS*)A->data;
722: MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv);
723: return(0);
724: }
728: PetscErrorCode MatZeroRows_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
729: {
730: PetscInt n_l = 0, *rows_l = NULL;
734: if (x && b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support");
735: if (n) {
736: PetscMalloc1(n,&rows_l);
737: ISGlobalToLocalMappingApply(A->rmap->mapping,IS_GTOLM_DROP,n,rows,&n_l,rows_l);
738: }
739: MatZeroRowsLocal(A,n_l,rows_l,diag,x,b);
740: PetscFree(rows_l);
741: return(0);
742: }
746: PetscErrorCode MatZeroRowsLocal_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
747: {
748: Mat_IS *is = (Mat_IS*)A->data;
750: PetscInt i;
751: PetscScalar *array;
754: if (x && b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support");
755: {
756: /*
757: Set up is->x as a "counting vector". This is in order to MatMult_IS
758: work properly in the interface nodes.
759: */
760: Vec counter;
761: MatCreateVecs(A,NULL,&counter);
762: VecSet(counter,0.);
763: VecSet(is->y,1.);
764: VecScatterBegin(is->rctx,is->y,counter,ADD_VALUES,SCATTER_REVERSE);
765: VecScatterEnd(is->rctx,is->y,counter,ADD_VALUES,SCATTER_REVERSE);
766: VecScatterBegin(is->rctx,counter,is->y,INSERT_VALUES,SCATTER_FORWARD);
767: VecScatterEnd(is->rctx,counter,is->y,INSERT_VALUES,SCATTER_FORWARD);
768: VecDestroy(&counter);
769: }
770: if (!n) {
771: is->pure_neumann = PETSC_TRUE;
772: } else {
773: is->pure_neumann = PETSC_FALSE;
775: VecGetArray(is->y,&array);
776: MatZeroRows(is->A,n,rows,diag,0,0);
777: for (i=0; i<n; i++) {
778: MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);
779: }
780: MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);
781: MatAssemblyEnd(is->A,MAT_FINAL_ASSEMBLY);
782: VecRestoreArray(is->y,&array);
783: }
784: return(0);
785: }
789: PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type)
790: {
791: Mat_IS *is = (Mat_IS*)A->data;
795: MatAssemblyBegin(is->A,type);
796: return(0);
797: }
801: PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type)
802: {
803: Mat_IS *is = (Mat_IS*)A->data;
807: MatAssemblyEnd(is->A,type);
808: return(0);
809: }
813: PetscErrorCode MatISGetLocalMat_IS(Mat mat,Mat *local)
814: {
815: Mat_IS *is = (Mat_IS*)mat->data;
818: *local = is->A;
819: return(0);
820: }
824: /*@
825: MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix.
827: Input Parameter:
828: . mat - the matrix
830: Output Parameter:
831: . local - the local matrix
833: Level: advanced
835: Notes:
836: This can be called if you have precomputed the nonzero structure of the
837: matrix and want to provide it to the inner matrix object to improve the performance
838: of the MatSetValues() operation.
840: .seealso: MATIS
841: @*/
842: PetscErrorCode MatISGetLocalMat(Mat mat,Mat *local)
843: {
849: PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat*),(mat,local));
850: return(0);
851: }
855: PetscErrorCode MatISSetLocalMat_IS(Mat mat,Mat local)
856: {
857: Mat_IS *is = (Mat_IS*)mat->data;
858: PetscInt nrows,ncols,orows,ocols;
862: if (is->A) {
863: MatGetSize(is->A,&orows,&ocols);
864: MatGetSize(local,&nrows,&ncols);
865: 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);
866: }
867: PetscObjectReference((PetscObject)local);
868: MatDestroy(&is->A);
869: is->A = local;
870: return(0);
871: }
875: /*@
876: MatISSetLocalMat - Replace the local matrix stored inside a MATIS object.
878: Input Parameter:
879: . mat - the matrix
880: . local - the local matrix
882: Output Parameter:
884: Level: advanced
886: Notes:
887: This can be called if you have precomputed the local matrix and
888: want to provide it to the matrix object MATIS.
890: .seealso: MATIS
891: @*/
892: PetscErrorCode MatISSetLocalMat(Mat mat,Mat local)
893: {
899: PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local));
900: return(0);
901: }
905: PetscErrorCode MatZeroEntries_IS(Mat A)
906: {
907: Mat_IS *a = (Mat_IS*)A->data;
911: MatZeroEntries(a->A);
912: return(0);
913: }
917: PetscErrorCode MatScale_IS(Mat A,PetscScalar a)
918: {
919: Mat_IS *is = (Mat_IS*)A->data;
923: MatScale(is->A,a);
924: return(0);
925: }
929: PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v)
930: {
931: Mat_IS *is = (Mat_IS*)A->data;
935: /* get diagonal of the local matrix */
936: MatGetDiagonal(is->A,is->y);
938: /* scatter diagonal back into global vector */
939: VecSet(v,0);
940: VecScatterBegin(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);
941: VecScatterEnd(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);
942: return(0);
943: }
947: PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg)
948: {
949: Mat_IS *a = (Mat_IS*)A->data;
953: MatSetOption(a->A,op,flg);
954: return(0);
955: }
959: /*@
960: MatCreateIS - Creates a "process" unassmembled matrix, it is assembled on each
961: process but not across processes.
963: Input Parameters:
964: + comm - MPI communicator that will share the matrix
965: . bs - block size of the matrix
966: . m,n,M,N - local and/or global sizes of the left and right vector used in matrix vector products
967: . rmap - local to global map for rows
968: - cmap - local to global map for cols
970: Output Parameter:
971: . A - the resulting matrix
973: Level: advanced
975: Notes: See MATIS for more details
976: m and n are NOT related to the size of the map, they are the size of the part of the vector owned
977: by that process. The sizes of rmap and cmap define the size of the local matrices.
978: If either rmap or cmap are NULL, than the matrix is assumed to be square
980: .seealso: MATIS, MatSetLocalToGlobalMapping()
981: @*/
982: PetscErrorCode MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping rmap,ISLocalToGlobalMapping cmap,Mat *A)
983: {
987: if (!rmap && !cmap) SETERRQ(comm,PETSC_ERR_USER,"You need to provide at least one of the mapping");
988: MatCreate(comm,A);
989: MatSetBlockSize(*A,bs);
990: MatSetSizes(*A,m,n,M,N);
991: MatSetType(*A,MATIS);
992: MatSetUp(*A);
993: if (rmap && cmap) {
994: MatSetLocalToGlobalMapping(*A,rmap,cmap);
995: } else if (!rmap) {
996: MatSetLocalToGlobalMapping(*A,cmap,cmap);
997: } else {
998: MatSetLocalToGlobalMapping(*A,rmap,rmap);
999: }
1000: return(0);
1001: }
1003: /*MC
1004: MATIS - MATIS = "is" - A matrix type to be used for using the non-overlapping domain decomposition type preconditioners (e.g. PCBDDC).
1005: This stores the matrices in globally unassembled form. Each processor
1006: assembles only its local Neumann problem and the parallel matrix vector
1007: product is handled "implicitly".
1009: Operations Provided:
1010: + MatMult()
1011: . MatMultAdd()
1012: . MatMultTranspose()
1013: . MatMultTransposeAdd()
1014: . MatZeroEntries()
1015: . MatSetOption()
1016: . MatZeroRows()
1017: . MatZeroRowsLocal()
1018: . MatSetValues()
1019: . MatSetValuesLocal()
1020: . MatScale()
1021: . MatGetDiagonal()
1022: - MatSetLocalToGlobalMapping()
1024: Options Database Keys:
1025: . -mat_type is - sets the matrix type to "is" during a call to MatSetFromOptions()
1027: Notes: Options prefix for the inner matrix are given by -is_mat_xxx
1029: You must call MatSetLocalToGlobalMapping() before using this matrix type.
1031: You can do matrix preallocation on the local matrix after you obtain it with
1032: MatISGetLocalMat(); otherwise, you could use MatISSetPreallocation()
1034: Level: advanced
1036: .seealso: Mat, MatISGetLocalMat(), MatSetLocalToGlobalMapping(), MatISSetPreallocation(), PCBDDC
1038: M*/
1042: PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A)
1043: {
1045: Mat_IS *b;
1048: PetscNewLog(A,&b);
1049: A->data = (void*)b;
1051: /* matrix ops */
1052: PetscMemzero(A->ops,sizeof(struct _MatOps));
1053: A->ops->mult = MatMult_IS;
1054: A->ops->multadd = MatMultAdd_IS;
1055: A->ops->multtranspose = MatMultTranspose_IS;
1056: A->ops->multtransposeadd = MatMultTransposeAdd_IS;
1057: A->ops->destroy = MatDestroy_IS;
1058: A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS;
1059: A->ops->setvalues = MatSetValues_IS;
1060: A->ops->setvalueslocal = MatSetValuesLocal_IS;
1061: A->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_IS;
1062: A->ops->zerorows = MatZeroRows_IS;
1063: A->ops->zerorowslocal = MatZeroRowsLocal_IS;
1064: A->ops->assemblybegin = MatAssemblyBegin_IS;
1065: A->ops->assemblyend = MatAssemblyEnd_IS;
1066: A->ops->view = MatView_IS;
1067: A->ops->zeroentries = MatZeroEntries_IS;
1068: A->ops->scale = MatScale_IS;
1069: A->ops->getdiagonal = MatGetDiagonal_IS;
1070: A->ops->setoption = MatSetOption_IS;
1071: A->ops->ishermitian = MatIsHermitian_IS;
1072: A->ops->issymmetric = MatIsSymmetric_IS;
1073: A->ops->duplicate = MatDuplicate_IS;
1075: /* special MATIS functions */
1076: PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS);
1077: PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS);
1078: PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",MatISGetMPIXAIJ_IS);
1079: PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",MatISSetPreallocation_IS);
1080: PetscObjectChangeTypeName((PetscObject)A,MATIS);
1081: return(0);
1082: }