Actual source code: sliced.c
petsc-3.4.5 2014-06-29
1: #include <petscdmsliced.h> /*I "petscdmsliced.h" I*/
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;
18: PetscErrorCode DMCreateMatrix_Sliced(DM dm, MatType mtype,Mat *J)
19: {
20: PetscErrorCode ierr;
21: PetscInt *globals,*sd_nnz,*so_nnz,rstart,bs,i;
22: ISLocalToGlobalMapping lmap,blmap;
23: void (*aij)(void) = NULL;
24: DM_Sliced *slice = (DM_Sliced*)dm->data;
27: bs = slice->bs;
28: MatCreate(PetscObjectComm((PetscObject)dm),J);
29: MatSetSizes(*J,slice->n*bs,slice->n*bs,PETSC_DETERMINE,PETSC_DETERMINE);
30: MatSetBlockSize(*J,bs);
31: MatSetType(*J,mtype);
32: MatSeqBAIJSetPreallocation(*J,bs,slice->d_nz,slice->d_nnz);
33: MatMPIBAIJSetPreallocation(*J,bs,slice->d_nz,slice->d_nnz,slice->o_nz,slice->o_nnz);
34: /* In general, we have to do extra work to preallocate for scalar (AIJ) matrices so we check whether it will do any
35: * good before going on with it. */
36: PetscObjectQueryFunction((PetscObject)*J,"MatMPIAIJSetPreallocation_C",&aij);
37: if (!aij) {
38: PetscObjectQueryFunction((PetscObject)*J,"MatSeqAIJSetPreallocation_C",&aij);
39: }
40: if (aij) {
41: if (bs == 1) {
42: MatSeqAIJSetPreallocation(*J,slice->d_nz,slice->d_nnz);
43: MatMPIAIJSetPreallocation(*J,slice->d_nz,slice->d_nnz,slice->o_nz,slice->o_nnz);
44: } else if (!slice->d_nnz) {
45: MatSeqAIJSetPreallocation(*J,slice->d_nz*bs,NULL);
46: MatMPIAIJSetPreallocation(*J,slice->d_nz*bs,NULL,slice->o_nz*bs,NULL);
47: } else {
48: /* The user has provided preallocation per block-row, convert it to per scalar-row respecting DMSlicedSetBlockFills() if applicable */
49: PetscMalloc2(slice->n*bs,PetscInt,&sd_nnz,(!!slice->o_nnz)*slice->n*bs,PetscInt,&so_nnz);
50: for (i=0; i<slice->n*bs; i++) {
51: sd_nnz[i] = (slice->d_nnz[i/bs]-1) * (slice->ofill ? slice->ofill->i[i%bs+1]-slice->ofill->i[i%bs] : bs)
52: + (slice->dfill ? slice->dfill->i[i%bs+1]-slice->dfill->i[i%bs] : bs);
53: if (so_nnz) {
54: so_nnz[i] = slice->o_nnz[i/bs] * (slice->ofill ? slice->ofill->i[i%bs+1]-slice->ofill->i[i%bs] : bs);
55: }
56: }
57: MatSeqAIJSetPreallocation(*J,slice->d_nz*bs,sd_nnz);
58: MatMPIAIJSetPreallocation(*J,slice->d_nz*bs,sd_nnz,slice->o_nz*bs,so_nnz);
59: PetscFree2(sd_nnz,so_nnz);
60: }
61: }
63: /* Set up the local to global map. For the scalar map, we have to translate to entry-wise indexing instead of block-wise. */
64: PetscMalloc((slice->n+slice->Nghosts)*bs*sizeof(PetscInt),&globals);
65: MatGetOwnershipRange(*J,&rstart,NULL);
66: for (i=0; i<slice->n*bs; i++) globals[i] = rstart + i;
68: for (i=0; i<slice->Nghosts*bs; i++) {
69: globals[slice->n*bs+i] = slice->ghosts[i/bs]*bs + i%bs;
70: }
71: ISLocalToGlobalMappingCreate(PETSC_COMM_SELF,(slice->n+slice->Nghosts)*bs,globals,PETSC_OWN_POINTER,&lmap);
72: ISLocalToGlobalMappingBlock(lmap,bs,&blmap);
73: MatSetLocalToGlobalMapping(*J,lmap,lmap);
74: MatSetLocalToGlobalMappingBlock(*J,blmap,blmap);
75: ISLocalToGlobalMappingDestroy(&lmap);
76: ISLocalToGlobalMappingDestroy(&blmap);
77: return(0);
78: }
82: /*@C
83: DMSlicedSetGhosts - Sets the global indices of other processes elements that will
84: be ghosts on this process
86: Not Collective
88: Input Parameters:
89: + slice - the DM object
90: . bs - block size
91: . nlocal - number of local (owned, non-ghost) blocks
92: . Nghosts - number of ghost blocks on this process
93: - ghosts - global indices of each ghost block
95: Level: advanced
97: .seealso DMDestroy(), DMCreateGlobalVector()
99: @*/
100: PetscErrorCode DMSlicedSetGhosts(DM dm,PetscInt bs,PetscInt nlocal,PetscInt Nghosts,const PetscInt ghosts[])
101: {
103: DM_Sliced *slice = (DM_Sliced*)dm->data;
107: PetscFree(slice->ghosts);
108: PetscMalloc(Nghosts*sizeof(PetscInt),&slice->ghosts);
109: PetscMemcpy(slice->ghosts,ghosts,Nghosts*sizeof(PetscInt));
110: slice->bs = bs;
111: slice->n = nlocal;
112: slice->Nghosts = Nghosts;
113: return(0);
114: }
118: /*@C
119: DMSlicedSetPreallocation - sets the matrix memory preallocation for matrices computed by DMSliced
121: Not Collective
123: Input Parameters:
124: + slice - the DM object
125: . d_nz - number of block nonzeros per block row in diagonal portion of local
126: submatrix (same for all local rows)
127: . d_nnz - array containing the number of block nonzeros in the various block rows
128: of the in diagonal portion of the local (possibly different for each block
129: row) or NULL.
130: . o_nz - number of block nonzeros per block row in the off-diagonal portion of local
131: submatrix (same for all local rows).
132: - o_nnz - array containing the number of nonzeros in the various block rows of the
133: off-diagonal portion of the local submatrix (possibly different for
134: each block row) or NULL.
136: Notes:
137: See MatMPIBAIJSetPreallocation() for more details on preallocation. If a scalar matrix (AIJ) is
138: obtained with DMSlicedGetMatrix(), the correct preallocation will be set, respecting DMSlicedSetBlockFills().
140: Level: advanced
142: .seealso DMDestroy(), DMCreateGlobalVector(), MatMPIAIJSetPreallocation(),
143: MatMPIBAIJSetPreallocation(), DMSlicedGetMatrix(), DMSlicedSetBlockFills()
145: @*/
146: PetscErrorCode DMSlicedSetPreallocation(DM dm,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
147: {
148: DM_Sliced *slice = (DM_Sliced*)dm->data;
152: slice->d_nz = d_nz;
153: slice->d_nnz = (PetscInt*)d_nnz;
154: slice->o_nz = o_nz;
155: slice->o_nnz = (PetscInt*)o_nnz;
156: return(0);
157: }
161: static PetscErrorCode DMSlicedSetBlockFills_Private(PetscInt bs,const PetscInt *fill,DMSlicedBlockFills **inf)
162: {
163: PetscErrorCode ierr;
164: PetscInt i,j,nz,*fi,*fj;
165: DMSlicedBlockFills *f;
169: if (*inf) {PetscFree3((*inf)->i,(*inf)->j,*inf);}
170: if (!fill) return(0);
171: for (i=0,nz=0; i<bs*bs; i++) if (fill[i]) nz++;
172: PetscMalloc3(1,DMSlicedBlockFills,&f,bs+1,PetscInt,&fi,nz,PetscInt,&fj);
173: f->bs = bs;
174: f->nz = nz;
175: f->i = fi;
176: f->j = fj;
177: for (i=0,nz=0; i<bs; i++) {
178: fi[i] = nz;
179: for (j=0; j<bs; j++) if (fill[i*bs+j]) fj[nz++] = j;
180: }
181: fi[i] = nz;
182: *inf = f;
183: return(0);
184: }
188: /*@C
189: DMSlicedSetBlockFills - Sets the fill pattern in each block for a multi-component problem
190: of the matrix returned by DMSlicedGetMatrix().
192: Logically Collective on DM
194: Input Parameter:
195: + sliced - the DM object
196: . dfill - the fill pattern in the diagonal block (may be NULL, means use dense block)
197: - ofill - the fill pattern in the off-diagonal blocks
199: Notes:
200: This only makes sense for multicomponent problems using scalar matrix formats (AIJ).
201: See DMDASetBlockFills() for example usage.
203: Level: advanced
205: .seealso DMSlicedGetMatrix(), DMDASetBlockFills()
206: @*/
207: PetscErrorCode DMSlicedSetBlockFills(DM dm,const PetscInt *dfill,const PetscInt *ofill)
208: {
209: DM_Sliced *slice = (DM_Sliced*)dm->data;
214: DMSlicedSetBlockFills_Private(slice->bs,dfill,&slice->dfill);
215: DMSlicedSetBlockFills_Private(slice->bs,ofill,&slice->ofill);
216: return(0);
217: }
221: static PetscErrorCode DMDestroy_Sliced(DM dm)
222: {
224: DM_Sliced *slice = (DM_Sliced*)dm->data;
227: PetscFree(slice->ghosts);
228: if (slice->dfill) {PetscFree3(slice->dfill,slice->dfill->i,slice->dfill->j);}
229: if (slice->ofill) {PetscFree3(slice->ofill,slice->ofill->i,slice->ofill->j);}
230: /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
231: PetscFree(slice);
232: return(0);
233: }
237: static PetscErrorCode DMCreateGlobalVector_Sliced(DM dm,Vec *gvec)
238: {
240: DM_Sliced *slice = (DM_Sliced*)dm->data;
245: *gvec = 0;
246: VecCreateGhostBlock(PetscObjectComm((PetscObject)dm),slice->bs,slice->n*slice->bs,PETSC_DETERMINE,slice->Nghosts,slice->ghosts,gvec);
247: VecSetDM(*gvec,dm);
248: return(0);
249: }
253: static PetscErrorCode DMGlobalToLocalBegin_Sliced(DM da,Vec g,InsertMode mode,Vec l)
254: {
256: PetscBool flg;
259: VecGhostIsLocalForm(g,l,&flg);
260: if (!flg) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Local vector is not local form of global vector");
261: VecGhostUpdateEnd(g,mode,SCATTER_FORWARD);
262: VecGhostUpdateBegin(g,mode,SCATTER_FORWARD);
263: return(0);
264: }
268: static PetscErrorCode DMGlobalToLocalEnd_Sliced(DM da,Vec g,InsertMode mode,Vec l)
269: {
271: PetscBool flg;
274: VecGhostIsLocalForm(g,l,&flg);
275: if (!flg) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Local vector is not local form of global vector");
276: VecGhostUpdateEnd(g,mode,SCATTER_FORWARD);
277: return(0);
278: }
280: /*MC
281: DMSLICED = "sliced" - A DM object that is used to manage data for a general graph. Uses VecCreateGhost() ghosted vectors for storing the fields
283: See DMCreateSliced() for details.
285: Level: intermediate
287: .seealso: DMType, DMCOMPOSITE, DMCreateSliced(), DMCreate()
288: M*/
292: PETSC_EXTERN PetscErrorCode DMCreate_Sliced(DM p)
293: {
295: DM_Sliced *slice;
298: PetscNewLog(p,DM_Sliced,&slice);
299: p->data = slice;
301: PetscObjectChangeTypeName((PetscObject)p,DMSLICED);
303: p->ops->createglobalvector = DMCreateGlobalVector_Sliced;
304: p->ops->creatematrix = DMCreateMatrix_Sliced;
305: p->ops->globaltolocalbegin = DMGlobalToLocalBegin_Sliced;
306: p->ops->globaltolocalend = DMGlobalToLocalEnd_Sliced;
307: p->ops->destroy = DMDestroy_Sliced;
308: return(0);
309: }
313: /*@C
314: DMSlicedCreate - Creates a DM object, used to manage data for a unstructured problem
316: Collective on MPI_Comm
318: Input Parameter:
319: + comm - the processors that will share the global vector
320: . bs - the block size
321: . nlocal - number of vector entries on this process
322: . Nghosts - number of ghost points needed on this process
323: . ghosts - global indices of all ghost points for this process
324: . d_nnz - matrix preallocation information representing coupling within this process
325: - o_nnz - matrix preallocation information representing coupling between this process and other processes
327: Output Parameters:
328: . slice - the slice object
330: Notes:
331: This DM does not support DMCreateLocalVector(), DMGlobalToLocalBegin(), and DMGlobalToLocalEnd() instead one directly uses
332: VecGhostGetLocalForm() and VecGhostRestoreLocalForm() to access the local representation and VecGhostUpdateBegin() and VecGhostUpdateEnd() to update
333: the ghost points.
335: One can use DMGlobalToLocalBegin(), and DMGlobalToLocalEnd() instead of VecGhostUpdateBegin() and VecGhostUpdateEnd().
337: Level: advanced
339: .seealso DMDestroy(), DMCreateGlobalVector(), DMSetType(), DMSLICED, DMSlicedSetGhosts(), DMSlicedSetPreallocation(), VecGhostUpdateBegin(), VecGhostUpdateEnd(),
340: VecGhostGetLocalForm(), VecGhostRestoreLocalForm()
342: @*/
343: PetscErrorCode DMSlicedCreate(MPI_Comm comm,PetscInt bs,PetscInt nlocal,PetscInt Nghosts,const PetscInt ghosts[], const PetscInt d_nnz[],const PetscInt o_nnz[],DM *dm)
344: {
349: DMCreate(comm,dm);
350: DMSetType(*dm,DMSLICED);
351: DMSlicedSetGhosts(*dm,bs,nlocal,Nghosts,ghosts);
352: if (d_nnz) {
353: DMSlicedSetPreallocation(*dm,0, d_nnz,0,o_nnz);
354: }
355: return(0);
356: }