Actual source code: da1.c
1: /*$Id: da1.c,v 1.129 2001/09/07 20:12:17 bsmith Exp $*/
3: /*
4: Code for manipulating distributed regular 1d arrays in parallel.
5: This file was created by Peter Mell 6/30/95
6: */
8: #include src/dm/da/daimpl.h
10: #if defined (PETSC_HAVE_AMS)
11: EXTERN_C_BEGIN
12: EXTERN int AMSSetFieldBlock_DA(AMS_Memory,char *,Vec);
13: EXTERN_C_END
14: #endif
16: #undef __FUNCT__
18: int DAView_1d(DA da,PetscViewer viewer)
19: {
20: int rank,ierr;
21: PetscTruth isascii,isdraw,isbinary;
24: MPI_Comm_rank(da->comm,&rank);
26: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
27: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);
28: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);
29: if (isascii) {
30: PetscViewerASCIISynchronizedPrintf(viewer,"Processor [%d] M %d m %d w %d s %dn",rank,da->M,
31: da->m,da->w,da->s);
32: PetscViewerASCIISynchronizedPrintf(viewer,"X range of indices: %d %dn",da->xs,da->xe);
33: PetscViewerFlush(viewer);
34: } else if (isdraw) {
35: PetscDraw draw;
36: double ymin = -1,ymax = 1,xmin = -1,xmax = da->M,x;
37: int base;
38: char node[10];
39: PetscTruth isnull;
41: PetscViewerDrawGetDraw(viewer,0,&draw);
42: PetscDrawIsNull(draw,&isnull); if (isnull) return(0);
44: PetscDrawSetCoordinates(draw,xmin,ymin,xmax,ymax);
45: PetscDrawSynchronizedClear(draw);
47: /* first processor draws all node lines */
48: if (!rank) {
49: int xmin_tmp;
50: ymin = 0.0; ymax = 0.3;
51:
52: /* ADIC doesn't like doubles in a for loop */
53: for (xmin_tmp =0; xmin_tmp < (int)da->M; xmin_tmp++) {
54: PetscDrawLine(draw,(double)xmin_tmp,ymin,(double)xmin_tmp,ymax,PETSC_DRAW_BLACK);
55: }
57: xmin = 0.0; xmax = da->M - 1;
58: PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_BLACK);
59: PetscDrawLine(draw,xmin,ymax,xmax,ymax,PETSC_DRAW_BLACK);
60: }
62: PetscDrawSynchronizedFlush(draw);
63: PetscDrawPause(draw);
65: /* draw my box */
66: ymin = 0; ymax = 0.3; xmin = da->xs / da->w; xmax = (da->xe / da->w) - 1;
67: PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_RED);
68: PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_RED);
69: PetscDrawLine(draw,xmin,ymax,xmax,ymax,PETSC_DRAW_RED);
70: PetscDrawLine(draw,xmax,ymin,xmax,ymax,PETSC_DRAW_RED);
72: /* Put in index numbers */
73: base = da->base / da->w;
74: for (x=xmin; x<=xmax; x++) {
75: sprintf(node,"%d",base++);
76: PetscDrawString(draw,x,ymin,PETSC_DRAW_RED,node);
77: }
79: PetscDrawSynchronizedFlush(draw);
80: PetscDrawPause(draw);
81: } else if (isbinary) {
82: DAView_Binary(da,viewer);
83: } else {
84: SETERRQ1(1,"Viewer type %s not supported for DA 1d",((PetscObject)viewer)->type_name);
85: }
86: return(0);
87: }
89: EXTERN int DAPublish_Petsc(PetscObject);
91: #undef __FUNCT__
93: /*@C
94: DACreate1d - Creates an object that will manage the communication of one-dimensional
95: regular array data that is distributed across some processors.
97: Collective on MPI_Comm
99: Input Parameters:
100: + comm - MPI communicator
101: . wrap - type of periodicity should the array have, if any. Use
102: either DA_NONPERIODIC or DA_XPERIODIC
103: . M - global dimension of the array
104: . dof - number of degrees of freedom per node
105: . lc - array containing number of nodes in the X direction on each processor,
106: or PETSC_NULL. If non-null, must be of length as m.
107: - s - stencil width
109: Output Parameter:
110: . inra - the resulting distributed array object
112: Options Database Key:
113: + -da_view - Calls DAView() at the conclusion of DACreate1d()
114: - -da_grid_x <nx> - number of grid points in x direction; can set if M < 0
116: Level: beginner
118: Notes:
119: If you are having problems with running out of memory than run with the option -da_noao
121: The array data itself is NOT stored in the DA, it is stored in Vec objects;
122: The appropriate vector objects can be obtained with calls to DACreateGlobalVector()
123: and DACreateLocalVector() and calls to VecDuplicate() if more are needed.
126: .keywords: distributed array, create, one-dimensional
128: .seealso: DADestroy(), DAView(), DACreate2d(), DACreate3d(), DAGlobalToLocalBegin(),
129: DAGlobalToLocalEnd(), DALocalToGlobal(), DALocalToLocalBegin(), DALocalToLocalEnd(),
130: DAGetInfo(), DACreateGlobalVector(), DACreateLocalVector(), DACreateNaturalVector(), DALoad(), DAView()
132: @*/
133: int DACreate1d(MPI_Comm comm,DAPeriodicType wrap,int M,int dof,int s,int *lc,DA *inra)
134: {
135: int rank,size,xs,xe,x,Xs,Xe,ierr,start,end,m;
136: int i,*idx,nn,left,refine_x = 2,tM = M;
137: PetscTruth flg1,flg2;
138: DA da;
139: Vec local,global;
140: VecScatter ltog,gtol;
141: IS to,from;
145: *inra = 0;
146: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
147: DMInitializePackage(PETSC_NULL);
148: #endif
150: if (dof < 1) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Must have 1 or more degrees of freedom per node: %d",dof);
151: if (s < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Stencil width cannot be negative: %d",s);
153: PetscOptionsBegin(comm,PETSC_NULL,"1d DA Options","DA");
154: if (M < 0) {
155: tM = -M;
156: PetscOptionsInt("-da_grid_x","Number of grid points in x direction","DACreate1d",tM,&tM,PETSC_NULL);
157: }
158: PetscOptionsInt("-da_refine_x","Refinement ratio in x direction","DACreate1d",refine_x,&refine_x,PETSC_NULL);
159: PetscOptionsEnd();
160: M = tM;
162: PetscHeaderCreate(da,_p_DA,struct _DAOps,DA_COOKIE,0,"DA",comm,DADestroy,DAView);
163: PetscLogObjectCreate(da);
164: da->bops->publish = DAPublish_Petsc;
165: da->ops->createglobalvector = DACreateGlobalVector;
166: da->ops->getinterpolation = DAGetInterpolation;
167: da->ops->getcoloring = DAGetColoring;
168: da->ops->getmatrix = DAGetMatrix;
169: da->ops->refine = DARefine;
170: PetscLogObjectMemory(da,sizeof(struct _p_DA));
171: da->dim = 1;
172: da->interptype = DA_Q1;
173: da->refine_x = refine_x;
174: PetscMalloc(dof*sizeof(char*),&da->fieldname);
175: PetscMemzero(da->fieldname,dof*sizeof(char*));
176: MPI_Comm_size(comm,&size);
177: MPI_Comm_rank(comm,&rank);
179: m = size;
181: if (M < m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"More processors than data points! %d %d",m,M);
182: if ((M-1) < s) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Array is too small for stencil! %d %d",M-1,s);
184: /*
185: Determine locally owned region
186: xs is the first local node number, x is the number of local nodes
187: */
188: if (!lc) {
189: PetscOptionsHasName(PETSC_NULL,"-da_partition_blockcomm",&flg1);
190: PetscOptionsHasName(PETSC_NULL,"-da_partition_nodes_at_end",&flg2);
191: if (flg1) { /* Block Comm type Distribution */
192: xs = rank*M/m;
193: x = (rank + 1)*M/m - xs;
194: } else if (flg2) { /* The odd nodes are evenly distributed across last nodes */
195: x = (M + rank)/m;
196: if (M/m == x) { xs = rank*x; }
197: else { xs = rank*(x-1) + (M+rank)%(x*m); }
198: } else { /* The odd nodes are evenly distributed across the first k nodes */
199: /* Regular PETSc Distribution */
200: x = M/m + ((M % m) > rank);
201: if (rank >= (M % m)) {xs = (rank * (int)(M/m) + M % m);}
202: else {xs = rank * (int)(M/m) + rank;}
203: }
204: } else {
205: x = lc[rank];
206: xs = 0;
207: for (i=0; i<rank; i++) {
208: xs += lc[i];
209: }
210: /* verify that data user provided is consistent */
211: left = xs;
212: for (i=rank; i<size; i++) {
213: left += lc[i];
214: }
215: if (left != M) {
216: SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Sum of lc across processors not equal to M %d %d",left,M);
217: }
218: }
220: /* From now on x,s,xs,xe,Xs,Xe are the exact location in the array */
221: x *= dof;
222: s *= dof; /* NOTE: here change s to be absolute stencil distance */
223: xs *= dof;
224: xe = xs + x;
226: /* determine ghost region */
227: if (wrap == DA_XPERIODIC) {
228: Xs = xs - s;
229: Xe = xe + s;
230: } else {
231: if ((xs-s) >= 0) Xs = xs-s; else Xs = 0;
232: if ((xe+s) <= M*dof) Xe = xe+s; else Xe = M*dof;
233: }
235: /* allocate the base parallel and sequential vectors */
236: da->Nlocal = x;
237: VecCreateMPIWithArray(comm,da->Nlocal,PETSC_DECIDE,0,&global);
238: VecSetBlockSize(global,dof);
239: da->nlocal = (Xe-Xs);
240: VecCreateSeqWithArray(PETSC_COMM_SELF,da->nlocal,0,&local);
241: VecSetBlockSize(local,dof);
242:
243: /* Create Local to Global Vector Scatter Context */
244: /* local to global inserts non-ghost point region into global */
245: VecGetOwnershipRange(global,&start,&end);
246: ISCreateStride(comm,x,start,1,&to);
247: ISCreateStride(comm,x,xs-Xs,1,&from);
248: VecScatterCreate(local,from,global,to,<og);
249: PetscLogObjectParent(da,to);
250: PetscLogObjectParent(da,from);
251: PetscLogObjectParent(da,ltog);
252: ISDestroy(from);
253: ISDestroy(to);
255: /* Create Global to Local Vector Scatter Context */
256: /* global to local must retrieve ghost points */
257: ISCreateStride(comm,(Xe-Xs),0,1,&to);
258:
259: PetscMalloc((x+2*s)*sizeof(int),&idx);
260: PetscLogObjectMemory(da,(x+2*s)*sizeof(int));
262: nn = 0;
263: if (wrap == DA_XPERIODIC) { /* Handle all cases with wrap first */
265: for (i=0; i<s; i++) { /* Left ghost points */
266: if ((xs-s+i)>=0) { idx[nn++] = xs-s+i;}
267: else { idx[nn++] = M*dof+(xs-s+i);}
268: }
270: for (i=0; i<x; i++) { idx [nn++] = xs + i;} /* Non-ghost points */
271:
272: for (i=0; i<s; i++) { /* Right ghost points */
273: if ((xe+i)<M*dof) { idx [nn++] = xe+i; }
274: else { idx [nn++] = (xe+i) - M*dof;}
275: }
276: } else { /* Now do all cases with no wrapping */
278: if (s <= xs) {for (i=0; i<s; i++) {idx[nn++] = xs - s + i;}}
279: else {for (i=0; i<xs; i++) {idx[nn++] = i;}}
281: for (i=0; i<x; i++) { idx [nn++] = xs + i;}
282:
283: if ((xe+s)<=M*dof) {for (i=0; i<s; i++) {idx[nn++]=xe+i;}}
284: else {for (i=xe; i<(M*dof); i++) {idx[nn++]=i; }}
285: }
287: ISCreateGeneral(comm,nn,idx,&from);
288: VecScatterCreate(global,from,local,to,>ol);
289: PetscLogObjectParent(da,to);
290: PetscLogObjectParent(da,from);
291: PetscLogObjectParent(da,gtol);
292: ISDestroy(to);
293: ISDestroy(from);
294: VecDestroy(local);
295: VecDestroy(global);
297: da->M = M; da->N = 1; da->m = m; da->n = 1;
298: da->xs = xs; da->xe = xe; da->ys = 0; da->ye = 1; da->zs = 0; da->ze = 1;
299: da->Xs = Xs; da->Xe = Xe; da->Ys = 0; da->Ye = 1; da->Zs = 0; da->Ze = 1;
300: da->P = 1; da->p = 1; da->w = dof; da->s = s/dof;
302: da->gtol = gtol;
303: da->ltog = ltog;
304: da->idx = idx;
305: da->Nl = nn;
306: da->base = xs;
307: da->ops->view = DAView_1d;
308: da->wrap = wrap;
309: da->stencil_type = DA_STENCIL_STAR;
311: /*
312: Set the local to global ordering in the global vector, this allows use
313: of VecSetValuesLocal().
314: */
315: ISLocalToGlobalMappingCreateNC(comm,nn,idx,&da->ltogmap);
316: ISLocalToGlobalMappingBlock(da->ltogmap,da->w,&da->ltogmapb);
317: PetscLogObjectParent(da,da->ltogmap);
319: da->ltol = PETSC_NULL;
320: da->ao = PETSC_NULL;
322: PetscOptionsHasName(PETSC_NULL,"-da_view",&flg1);
323: if (flg1) {DAView(da,PETSC_VIEWER_STDOUT_(da->comm));}
324: PetscOptionsHasName(PETSC_NULL,"-da_view_draw",&flg1);
325: if (flg1) {DAView(da,PETSC_VIEWER_DRAW_(da->comm));}
326: PetscOptionsHasName(PETSC_NULL,"-help",&flg1);
327: if (flg1) {DAPrintHelp(da);}
328: *inra = da;
329: PetscPublishAll(da);
330: return(0);
331: }