Actual source code: sliced.c
petsc-3.11.4 2019-09-28
1: #include <petscdmsliced.h>
2: #include <petscmat.h>
3: #include <petsc/private/dmimpl.h>
5: /* CSR storage of the nonzero structure of a bs*bs matrix */
6: typedef struct {
7: PetscInt bs,nz,*i,*j;
8: } DMSlicedBlockFills;
10: typedef struct {
11: PetscInt bs,n,N,Nghosts,*ghosts;
12: PetscInt d_nz,o_nz,*d_nnz,*o_nnz;
13: DMSlicedBlockFills *dfill,*ofill;
14: } DM_Sliced;
16: PetscErrorCode DMCreateMatrix_Sliced(DM dm, Mat *J)
17: {
18: PetscErrorCode ierr;
19: PetscInt *globals,*sd_nnz,*so_nnz,rstart,bs,i;
20: ISLocalToGlobalMapping lmap;
21: void (*aij)(void) = NULL;
22: DM_Sliced *slice = (DM_Sliced*)dm->data;
25: bs = slice->bs;
26: MatCreate(PetscObjectComm((PetscObject)dm),J);
27: MatSetSizes(*J,slice->n*bs,slice->n*bs,PETSC_DETERMINE,PETSC_DETERMINE);
28: MatSetBlockSize(*J,bs);
29: MatSetType(*J,dm->mattype);
30: MatSeqBAIJSetPreallocation(*J,bs,slice->d_nz,slice->d_nnz);
31: MatMPIBAIJSetPreallocation(*J,bs,slice->d_nz,slice->d_nnz,slice->o_nz,slice->o_nnz);
32: /* In general, we have to do extra work to preallocate for scalar (AIJ) matrices so we check whether it will do any
33: * good before going on with it. */
34: PetscObjectQueryFunction((PetscObject)*J,"MatMPIAIJSetPreallocation_C",&aij);
35: if (!aij) {
36: PetscObjectQueryFunction((PetscObject)*J,"MatSeqAIJSetPreallocation_C",&aij);
37: }
38: if (aij) {
39: if (bs == 1) {
40: MatSeqAIJSetPreallocation(*J,slice->d_nz,slice->d_nnz);
41: MatMPIAIJSetPreallocation(*J,slice->d_nz,slice->d_nnz,slice->o_nz,slice->o_nnz);
42: } else if (!slice->d_nnz) {
43: MatSeqAIJSetPreallocation(*J,slice->d_nz*bs,NULL);
44: MatMPIAIJSetPreallocation(*J,slice->d_nz*bs,NULL,slice->o_nz*bs,NULL);
45: } else {
46: /* The user has provided preallocation per block-row, convert it to per scalar-row respecting DMSlicedSetBlockFills() if applicable */
47: PetscMalloc2(slice->n*bs,&sd_nnz,(!!slice->o_nnz)*slice->n*bs,&so_nnz);
48: for (i=0; i<slice->n*bs; i++) {
49: sd_nnz[i] = (slice->d_nnz[i/bs]-1) * (slice->ofill ? slice->ofill->i[i%bs+1]-slice->ofill->i[i%bs] : bs)
50: + (slice->dfill ? slice->dfill->i[i%bs+1]-slice->dfill->i[i%bs] : bs);
51: if (so_nnz) {
52: so_nnz[i] = slice->o_nnz[i/bs] * (slice->ofill ? slice->ofill->i[i%bs+1]-slice->ofill->i[i%bs] : bs);
53: }
54: }
55: MatSeqAIJSetPreallocation(*J,slice->d_nz*bs,sd_nnz);
56: MatMPIAIJSetPreallocation(*J,slice->d_nz*bs,sd_nnz,slice->o_nz*bs,so_nnz);
57: PetscFree2(sd_nnz,so_nnz);
58: }
59: }
61: /* Set up the local to global map. For the scalar map, we have to translate to entry-wise indexing instead of block-wise. */
62: PetscMalloc1(slice->n+slice->Nghosts,&globals);
63: MatGetOwnershipRange(*J,&rstart,NULL);
64: for (i=0; i<slice->n; i++) globals[i] = rstart/bs + i;
66: for (i=0; i<slice->Nghosts; i++) {
67: globals[slice->n+i] = slice->ghosts[i];
68: }
69: ISLocalToGlobalMappingCreate(PETSC_COMM_SELF,bs,slice->n+slice->Nghosts,globals,PETSC_OWN_POINTER,&lmap);
70: MatSetLocalToGlobalMapping(*J,lmap,lmap);
71: ISLocalToGlobalMappingDestroy(&lmap);
72: return(0);
73: }
75: /*@C
76: DMSlicedSetGhosts - Sets the global indices of other processes elements that will
77: be ghosts on this process
79: Not Collective
81: Input Parameters:
82: + slice - the DM object
83: . bs - block size
84: . nlocal - number of local (owned, non-ghost) blocks
85: . Nghosts - number of ghost blocks on this process
86: - ghosts - global indices of each ghost block
88: Level: advanced
90: .seealso DMDestroy(), DMCreateGlobalVector()
92: @*/
93: PetscErrorCode DMSlicedSetGhosts(DM dm,PetscInt bs,PetscInt nlocal,PetscInt Nghosts,const PetscInt ghosts[])
94: {
96: DM_Sliced *slice = (DM_Sliced*)dm->data;
100: PetscFree(slice->ghosts);
101: PetscMalloc1(Nghosts,&slice->ghosts);
102: PetscMemcpy(slice->ghosts,ghosts,Nghosts*sizeof(PetscInt));
103: slice->bs = bs;
104: slice->n = nlocal;
105: slice->Nghosts = Nghosts;
106: return(0);
107: }
109: /*@C
110: DMSlicedSetPreallocation - sets the matrix memory preallocation for matrices computed by DMSliced
112: Not Collective
114: Input Parameters:
115: + slice - the DM object
116: . d_nz - number of block nonzeros per block row in diagonal portion of local
117: submatrix (same for all local rows)
118: . d_nnz - array containing the number of block nonzeros in the various block rows
119: of the in diagonal portion of the local (possibly different for each block
120: row) or NULL.
121: . o_nz - number of block nonzeros per block row in the off-diagonal portion of local
122: submatrix (same for all local rows).
123: - o_nnz - array containing the number of nonzeros in the various block rows of the
124: off-diagonal portion of the local submatrix (possibly different for
125: each block row) or NULL.
127: Notes:
128: See MatMPIBAIJSetPreallocation() for more details on preallocation. If a scalar matrix (AIJ) is
129: obtained with DMSlicedGetMatrix(), the correct preallocation will be set, respecting DMSlicedSetBlockFills().
131: Level: advanced
133: .seealso DMDestroy(), DMCreateGlobalVector(), MatMPIAIJSetPreallocation(),
134: MatMPIBAIJSetPreallocation(), DMSlicedGetMatrix(), DMSlicedSetBlockFills()
136: @*/
137: PetscErrorCode DMSlicedSetPreallocation(DM dm,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
138: {
139: DM_Sliced *slice = (DM_Sliced*)dm->data;
143: slice->d_nz = d_nz;
144: slice->d_nnz = (PetscInt*)d_nnz;
145: slice->o_nz = o_nz;
146: slice->o_nnz = (PetscInt*)o_nnz;
147: return(0);
148: }
150: static PetscErrorCode DMSlicedSetBlockFills_Private(PetscInt bs,const PetscInt *fill,DMSlicedBlockFills **inf)
151: {
152: PetscErrorCode ierr;
153: PetscInt i,j,nz,*fi,*fj;
154: DMSlicedBlockFills *f;
158: if (*inf) {PetscFree3(*inf,(*inf)->i,(*inf)->j);}
159: if (!fill) return(0);
160: for (i=0,nz=0; i<bs*bs; i++) if (fill[i]) nz++;
161: PetscMalloc3(1,&f,bs+1,&fi,nz,&fj);
162: f->bs = bs;
163: f->nz = nz;
164: f->i = fi;
165: f->j = fj;
166: for (i=0,nz=0; i<bs; i++) {
167: fi[i] = nz;
168: for (j=0; j<bs; j++) if (fill[i*bs+j]) fj[nz++] = j;
169: }
170: fi[i] = nz;
171: *inf = f;
172: return(0);
173: }
175: /*@C
176: DMSlicedSetBlockFills - Sets the fill pattern in each block for a multi-component problem
177: of the matrix returned by DMSlicedGetMatrix().
179: Logically Collective on DM
181: Input Parameter:
182: + sliced - the DM object
183: . dfill - the fill pattern in the diagonal block (may be NULL, means use dense block)
184: - ofill - the fill pattern in the off-diagonal blocks
186: Notes:
187: This only makes sense for multicomponent problems using scalar matrix formats (AIJ).
188: See DMDASetBlockFills() for example usage.
190: Level: advanced
192: .seealso DMSlicedGetMatrix(), DMDASetBlockFills()
193: @*/
194: PetscErrorCode DMSlicedSetBlockFills(DM dm,const PetscInt *dfill,const PetscInt *ofill)
195: {
196: DM_Sliced *slice = (DM_Sliced*)dm->data;
201: DMSlicedSetBlockFills_Private(slice->bs,dfill,&slice->dfill);
202: DMSlicedSetBlockFills_Private(slice->bs,ofill,&slice->ofill);
203: return(0);
204: }
206: static PetscErrorCode DMDestroy_Sliced(DM dm)
207: {
209: DM_Sliced *slice = (DM_Sliced*)dm->data;
212: PetscFree(slice->ghosts);
213: if (slice->dfill) {PetscFree3(slice->dfill,slice->dfill->i,slice->dfill->j);}
214: if (slice->ofill) {PetscFree3(slice->ofill,slice->ofill->i,slice->ofill->j);}
215: /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
216: PetscFree(slice);
217: return(0);
218: }
220: static PetscErrorCode DMCreateGlobalVector_Sliced(DM dm,Vec *gvec)
221: {
223: DM_Sliced *slice = (DM_Sliced*)dm->data;
228: *gvec = 0;
229: VecCreateGhostBlock(PetscObjectComm((PetscObject)dm),slice->bs,slice->n*slice->bs,PETSC_DETERMINE,slice->Nghosts,slice->ghosts,gvec);
230: VecSetDM(*gvec,dm);
231: return(0);
232: }
234: static PetscErrorCode DMGlobalToLocalBegin_Sliced(DM da,Vec g,InsertMode mode,Vec l)
235: {
237: PetscBool flg;
240: VecGhostIsLocalForm(g,l,&flg);
241: if (!flg) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Local vector is not local form of global vector");
242: VecGhostUpdateEnd(g,mode,SCATTER_FORWARD);
243: VecGhostUpdateBegin(g,mode,SCATTER_FORWARD);
244: return(0);
245: }
247: static PetscErrorCode DMGlobalToLocalEnd_Sliced(DM da,Vec g,InsertMode mode,Vec l)
248: {
250: PetscBool flg;
253: VecGhostIsLocalForm(g,l,&flg);
254: if (!flg) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Local vector is not local form of global vector");
255: VecGhostUpdateEnd(g,mode,SCATTER_FORWARD);
256: return(0);
257: }
259: /*MC
260: DMSLICED = "sliced" - A DM object that is used to manage data for a general graph. Uses VecCreateGhost() ghosted vectors for storing the fields
262: See DMCreateSliced() for details.
264: Level: intermediate
266: .seealso: DMType, DMCOMPOSITE, DMCreateSliced(), DMCreate()
267: M*/
269: PETSC_EXTERN PetscErrorCode DMCreate_Sliced(DM p)
270: {
272: DM_Sliced *slice;
275: PetscNewLog(p,&slice);
276: p->data = slice;
278: PetscObjectChangeTypeName((PetscObject)p,DMSLICED);
280: p->ops->createglobalvector = DMCreateGlobalVector_Sliced;
281: p->ops->creatematrix = DMCreateMatrix_Sliced;
282: p->ops->globaltolocalbegin = DMGlobalToLocalBegin_Sliced;
283: p->ops->globaltolocalend = DMGlobalToLocalEnd_Sliced;
284: p->ops->destroy = DMDestroy_Sliced;
285: return(0);
286: }
288: /*@C
289: DMSlicedCreate - Creates a DM object, used to manage data for a unstructured problem
291: Collective on MPI_Comm
293: Input Parameter:
294: + comm - the processors that will share the global vector
295: . bs - the block size
296: . nlocal - number of vector entries on this process
297: . Nghosts - number of ghost points needed on this process
298: . ghosts - global indices of all ghost points for this process
299: . d_nnz - matrix preallocation information representing coupling within this process
300: - o_nnz - matrix preallocation information representing coupling between this process and other processes
302: Output Parameters:
303: . slice - the slice object
305: Notes:
306: This DM does not support DMCreateLocalVector(), DMGlobalToLocalBegin(), and DMGlobalToLocalEnd() instead one directly uses
307: VecGhostGetLocalForm() and VecGhostRestoreLocalForm() to access the local representation and VecGhostUpdateBegin() and VecGhostUpdateEnd() to update
308: the ghost points.
310: One can use DMGlobalToLocalBegin(), and DMGlobalToLocalEnd() instead of VecGhostUpdateBegin() and VecGhostUpdateEnd().
312: Level: advanced
314: .seealso DMDestroy(), DMCreateGlobalVector(), DMSetType(), DMSLICED, DMSlicedSetGhosts(), DMSlicedSetPreallocation(), VecGhostUpdateBegin(), VecGhostUpdateEnd(),
315: VecGhostGetLocalForm(), VecGhostRestoreLocalForm()
317: @*/
318: PetscErrorCode DMSlicedCreate(MPI_Comm comm,PetscInt bs,PetscInt nlocal,PetscInt Nghosts,const PetscInt ghosts[], const PetscInt d_nnz[],const PetscInt o_nnz[],DM *dm)
319: {
324: DMCreate(comm,dm);
325: DMSetType(*dm,DMSLICED);
326: DMSlicedSetGhosts(*dm,bs,nlocal,Nghosts,ghosts);
327: if (d_nnz) {
328: DMSlicedSetPreallocation(*dm,0, d_nnz,0,o_nnz);
329: }
330: return(0);
331: }