Actual source code: da2.c
1: /*$Id: da2.c,v 1.180 2001/09/07 20:12:17 bsmith Exp $*/
2:
3: #include src/dm/da/daimpl.h
5: #undef __FUNCT__
7: int DAGetOwnershipRange(DA da,int **lx,int **ly,int **lz)
8: {
11: if (lx) *lx = da->lx;
12: if (ly) *ly = da->ly;
13: if (lz) *lz = da->lz;
14: return(0);
15: }
17: #undef __FUNCT__
19: int DAView_2d(DA da,PetscViewer viewer)
20: {
21: int rank,ierr;
22: PetscTruth isascii,isdraw,isbinary;
25: MPI_Comm_rank(da->comm,&rank);
27: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
28: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);
29: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);
30: if (isascii) {
31: PetscViewerASCIISynchronizedPrintf(viewer,"Processor [%d] M %d N %d m %d n %d w %d s %dn",rank,da->M,
32: da->N,da->m,da->n,da->w,da->s);
33: PetscViewerASCIISynchronizedPrintf(viewer,"X range of indices: %d %d, Y range of indices: %d %dn",da->xs,da->xe,da->ys,da->ye);
34: PetscViewerFlush(viewer);
35: } else if (isdraw) {
36: PetscDraw draw;
37: double ymin = -1*da->s-1,ymax = da->N+da->s;
38: double xmin = -1*da->s-1,xmax = da->M+da->s;
39: double x,y;
40: int base,*idx;
41: char node[10];
42: PetscTruth isnull;
43:
44: PetscViewerDrawGetDraw(viewer,0,&draw);
45: PetscDrawIsNull(draw,&isnull); if (isnull) return(0);
46: PetscDrawSetCoordinates(draw,xmin,ymin,xmax,ymax);
47: PetscDrawSynchronizedClear(draw);
49: /* first processor draw all node lines */
50: if (!rank) {
51: ymin = 0.0; ymax = da->N - 1;
52: for (xmin=0; xmin<da->M; xmin++) {
53: PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_BLACK);
54: }
55: xmin = 0.0; xmax = da->M - 1;
56: for (ymin=0; ymin<da->N; ymin++) {
57: PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_BLACK);
58: }
59: }
60: PetscDrawSynchronizedFlush(draw);
61: PetscDrawPause(draw);
63: /* draw my box */
64: ymin = da->ys; ymax = da->ye - 1; xmin = da->xs/da->w;
65: xmax =(da->xe-1)/da->w;
66: PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_RED);
67: PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_RED);
68: PetscDrawLine(draw,xmin,ymax,xmax,ymax,PETSC_DRAW_RED);
69: PetscDrawLine(draw,xmax,ymin,xmax,ymax,PETSC_DRAW_RED);
71: /* put in numbers */
72: base = (da->base)/da->w;
73: for (y=ymin; y<=ymax; y++) {
74: for (x=xmin; x<=xmax; x++) {
75: sprintf(node,"%d",base++);
76: PetscDrawString(draw,x,y,PETSC_DRAW_BLACK,node);
77: }
78: }
80: PetscDrawSynchronizedFlush(draw);
81: PetscDrawPause(draw);
82: /* overlay ghost numbers, useful for error checking */
83: /* put in numbers */
85: base = 0; idx = da->idx;
86: ymin = da->Ys; ymax = da->Ye; xmin = da->Xs; xmax = da->Xe;
87: for (y=ymin; y<ymax; y++) {
88: for (x=xmin; x<xmax; x++) {
89: if ((base % da->w) == 0) {
90: sprintf(node,"%d",idx[base]/da->w);
91: PetscDrawString(draw,x/da->w,y,PETSC_DRAW_BLUE,node);
92: }
93: base++;
94: }
95: }
96: PetscDrawSynchronizedFlush(draw);
97: PetscDrawPause(draw);
98: } else if (isbinary) {
99: DAView_Binary(da,viewer);
100: } else {
101: SETERRQ1(1,"Viewer type %s not supported for DA2d",((PetscObject)viewer)->type_name);
102: }
103: return(0);
104: }
106: #if defined(PETSC_HAVE_AMS)
107: /*
108: This function tells the AMS the layout of the vectors, it is called
109: in the VecPublish_xx routines.
110: */
111: EXTERN_C_BEGIN
112: #undef __FUNCT__
114: int AMSSetFieldBlock_DA(AMS_Memory amem,char *name,Vec vec)
115: {
116: int ierr,dof,dim,ends[4],shift = 0,starts[] = {0,0,0,0};
117: DA da = 0;
118: PetscTruth isseq,ismpi;
121: if (((PetscObject)vec)->amem < 0) return(0); /* return if not published */
123: PetscObjectQuery((PetscObject)vec,"DA",(PetscObject*)&da);
124: if (!da) return(0);
125: DAGetInfo(da,&dim,0,0,0,0,0,0,&dof,0,0,0);
126: if (dof > 1) {dim++; shift = 1; ends[0] = dof;}
128: PetscTypeCompare((PetscObject)vec,VECSEQ,&isseq);
129: PetscTypeCompare((PetscObject)vec,VECMPI,&ismpi);
130: if (isseq) {
131: DAGetGhostCorners(da,0,0,0,ends+shift,ends+shift+1,ends+shift+2);
132: ends[shift] += starts[shift]-1;
133: ends[shift+1] += starts[shift+1]-1;
134: ends[shift+2] += starts[shift+2]-1;
135: AMS_Memory_set_field_block(amem,name,dim,starts,ends);
136: if (ierr) {
137: char *message;
138: AMS_Explain_error(ierr,&message);
139: SETERRQ(ierr,message);
140: }
141: } else if (ismpi) {
142: DAGetCorners(da,starts+shift,starts+shift+1,starts+shift+2,
143: ends+shift,ends+shift+1,ends+shift+2);
144: ends[shift] += starts[shift]-1;
145: ends[shift+1] += starts[shift+1]-1;
146: ends[shift+2] += starts[shift+2]-1;
147: AMS_Memory_set_field_block(amem,name,dim,starts,ends);
148: if (ierr) {
149: char *message;
150: AMS_Explain_error(ierr,&message);
151: SETERRQ(ierr,message);
152: }
153: } else {
154: SETERRQ1(1,"Wrong vector type %s for this call",((PetscObject)vec)->type_name);
155: }
157: return(0);
158: }
159: EXTERN_C_END
160: #endif
162: #undef __FUNCT__
164: int DAPublish_Petsc(PetscObject obj)
165: {
166: #if defined(PETSC_HAVE_AMS)
167: DA v = (DA) obj;
168: int ierr;
169: #endif
173: #if defined(PETSC_HAVE_AMS)
174: /* if it is already published then return */
175: if (v->amem >=0) return(0);
177: PetscObjectPublishBaseBegin(obj);
178: PetscObjectPublishBaseEnd(obj);
179: #endif
181: return(0);
182: }
185: #undef __FUNCT__
187: /*@C
188: DACreate2d - Creates an object that will manage the communication of two-dimensional
189: regular array data that is distributed across some processors.
191: Collective on MPI_Comm
193: Input Parameters:
194: + comm - MPI communicator
195: . wrap - type of periodicity should the array have.
196: Use one of DA_NONPERIODIC, DA_XPERIODIC, DA_YPERIODIC, or DA_XYPERIODIC.
197: . stencil_type - stencil type. Use either DA_STENCIL_BOX or DA_STENCIL_STAR.
198: . M,N - global dimension in each direction of the array
199: . m,n - corresponding number of processors in each dimension
200: (or PETSC_DECIDE to have calculated)
201: . dof - number of degrees of freedom per node
202: . s - stencil width
203: - lx, ly - arrays containing the number of nodes in each cell along
204: the x and y coordinates, or PETSC_NULL. If non-null, these
205: must be of length as m and n, and the corresponding
206: m and n cannot be PETSC_DECIDE. The sum of the lx[] entries
207: must be M, and the sum of the ly[] entries must be N.
209: Output Parameter:
210: . inra - the resulting distributed array object
212: Options Database Key:
213: + -da_view - Calls DAView() at the conclusion of DACreate2d()
214: . -da_grid_x <nx> - number of grid points in x direction, if M < 0
215: . -da_grid_y <ny> - number of grid points in y direction, if N < 0
216: . -da_processors_x <nx> - number of processors in x direction
217: - -da_processors_y <ny> - number of processors in y direction
219: Level: beginner
221: Notes:
222: The stencil type DA_STENCIL_STAR with width 1 corresponds to the
223: standard 5-pt stencil, while DA_STENCIL_BOX with width 1 denotes
224: the standard 9-pt stencil.
226: The array data itself is NOT stored in the DA, it is stored in Vec objects;
227: The appropriate vector objects can be obtained with calls to DACreateGlobalVector()
228: and DACreateLocalVector() and calls to VecDuplicate() if more are needed.
230: .keywords: distributed array, create, two-dimensional
232: .seealso: DADestroy(), DAView(), DACreate1d(), DACreate3d(), DAGlobalToLocalBegin(),
233: DAGlobalToLocalEnd(), DALocalToGlobal(), DALocalToLocalBegin(), DALocalToLocalEnd(),
234: DAGetInfo(), DACreateGlobalVector(), DACreateLocalVector(), DACreateNaturalVector(), DALoad(), DAView()
236: @*/
237: int DACreate2d(MPI_Comm comm,DAPeriodicType wrap,DAStencilType stencil_type,
238: int M,int N,int m,int n,int dof,int s,int *lx,int *ly,DA *inra)
239: {
240: int rank,size,xs,xe,ys,ye,x,y,Xs,Xe,Ys,Ye,ierr,start,end;
241: int up,down,left,i,n0,n1,n2,n3,n5,n6,n7,n8,*idx,nn;
242: int xbase,*bases,*ldims,j,x_t,y_t,s_t,base,count;
243: int s_x,s_y; /* s proportionalized to w */
244: int *flx = 0,*fly = 0;
245: int sn0 = 0,sn2 = 0,sn6 = 0,sn8 = 0,refine_x = 2, refine_y = 2,tM = M,tN = N;
246: PetscTruth flg1,flg2;
247: DA da;
248: Vec local,global;
249: VecScatter ltog,gtol;
250: IS to,from;
254: *inra = 0;
255: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
256: DMInitializePackage(PETSC_NULL);
257: #endif
259: if (dof < 1) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Must have 1 or more degrees of freedom per node: %d",dof);
260: if (s < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Stencil width cannot be negative: %d",s);
262: PetscOptionsBegin(comm,PETSC_NULL,"2d DA Options","DA");
263: if (M < 0){
264: tM = -M;
265: PetscOptionsInt("-da_grid_x","Number of grid points in x direction","DACreate2d",tM,&tM,PETSC_NULL);
266: }
267: if (N < 0){
268: tN = -N;
269: PetscOptionsInt("-da_grid_y","Number of grid points in y direction","DACreate2d",tN,&tN,PETSC_NULL);
270: }
271: PetscOptionsInt("-da_processors_x","Number of processors in x direction","DACreate2d",m,&m,PETSC_NULL);
272: PetscOptionsInt("-da_processors_y","Number of processors in y direction","DACreate2d",n,&n,PETSC_NULL);
273: PetscOptionsInt("-da_refine_x","Refinement ratio in x direction","DACreate2d",refine_x,&refine_x,PETSC_NULL);
274: PetscOptionsInt("-da_refine_y","Refinement ratio in y direction","DACreate2d",refine_y,&refine_y,PETSC_NULL);
275: PetscOptionsEnd();
276: M = tM; N = tN;
278: PetscHeaderCreate(da,_p_DA,struct _DAOps,DA_COOKIE,0,"DA",comm,DADestroy,DAView);
279: PetscLogObjectCreate(da);
280: da->bops->publish = DAPublish_Petsc;
281: da->ops->createglobalvector = DACreateGlobalVector;
282: da->ops->getinterpolation = DAGetInterpolation;
283: da->ops->getcoloring = DAGetColoring;
284: da->ops->getmatrix = DAGetMatrix;
285: da->ops->refine = DARefine;
286: PetscLogObjectMemory(da,sizeof(struct _p_DA));
287: da->dim = 2;
288: da->interptype = DA_Q1;
289: da->refine_x = refine_x;
290: da->refine_y = refine_y;
291: PetscMalloc(dof*sizeof(char*),&da->fieldname);
292: PetscMemzero(da->fieldname,dof*sizeof(char*));
294: MPI_Comm_size(comm,&size);
295: MPI_Comm_rank(comm,&rank);
297: if (m != PETSC_DECIDE) {
298: if (m < 1) {SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in X direction: %d",m);}
299: else if (m > size) {SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in X direction: %d %d",m,size);}
300: }
301: if (n != PETSC_DECIDE) {
302: if (n < 1) {SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in Y direction: %d",n);}
303: else if (n > size) {SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in Y direction: %d %d",n,size);}
304: }
306: if (m == PETSC_DECIDE || n == PETSC_DECIDE) {
307: /* try for squarish distribution */
308: /* This should use MPI_Dims_create instead */
309: m = (int)(0.5 + sqrt(((double)M)*((double)size)/((double)N)));
310: if (!m) m = 1;
311: while (m > 0) {
312: n = size/m;
313: if (m*n == size) break;
314: m--;
315: }
316: if (M > N && m < n) {int _m = m; m = n; n = _m;}
317: if (m*n != size) SETERRQ(PETSC_ERR_PLIB,"Internally Created Bad Partition");
318: } else if (m*n != size) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Given Bad partition");
320: if (M < m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Partition in x direction is too fine! %d %d",M,m);
321: if (N < n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Partition in y direction is too fine! %d %d",N,n);
323: /*
324: We should create an MPI Cartesian topology here, with reorder
325: set to true. That would create a NEW communicator that we would
326: need to use for operations on this distributed array
327: */
328: PetscOptionsHasName(PETSC_NULL,"-da_partition_nodes_at_end",&flg2);
330: /*
331: Determine locally owned region
332: xs is the first local node number, x is the number of local nodes
333: */
334: if (lx) { /* user sets distribution */
335: x = lx[rank % m];
336: xs = 0;
337: for (i=0; i<(rank % m); i++) {
338: xs += lx[i];
339: }
340: left = xs;
341: for (i=(rank % m); i<m; i++) {
342: left += lx[i];
343: }
344: if (left != M) {
345: SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Sum of lx across processors not equal to M: %d %d",left,M);
346: }
347: } else if (flg2) {
348: x = (M + rank%m)/m;
349: if (x < s) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column width is too thin for stencil! %d %d",x,s);
350: if (M/m == x) { xs = (rank % m)*x; }
351: else { xs = (rank % m)*(x-1) + (M+(rank % m))%(x*m); }
352: SETERRQ(PETSC_ERR_SUP,"-da_partition_nodes_at_end not supported");
353: } else { /* Normal PETSc distribution */
354: x = M/m + ((M % m) > (rank % m));
355: if (x < s) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column width is too thin for stencil! %d %d",x,s);
356: if ((M % m) > (rank % m)) { xs = (rank % m)*x; }
357: else { xs = (M % m)*(x+1) + ((rank % m)-(M % m))*x; }
358: PetscMalloc(m*sizeof(int),&lx);
359: flx = lx;
360: for (i=0; i<m; i++) {
361: lx[i] = M/m + ((M % m) > i);
362: }
363: }
365: /*
366: Determine locally owned region
367: ys is the first local node number, y is the number of local nodes
368: */
369: if (ly) { /* user sets distribution */
370: y = ly[rank/m];
371: ys = 0;
372: for (i=0; i<(rank/m); i++) {
373: ys += ly[i];
374: }
375: left = ys;
376: for (i=(rank/m); i<n; i++) {
377: left += ly[i];
378: }
379: if (left != N) {
380: SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Sum of ly across processors not equal to N: %d %d",left,N);
381: }
382: } else if (flg2) {
383: y = (N + rank/m)/n;
384: if (y < s) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row width is too thin for stencil! %d %d",y,s);
385: if (N/n == y) { ys = (rank/m)*y; }
386: else { ys = (rank/m)*(y-1) + (N+(rank/m))%(y*n); }
387: SETERRQ(PETSC_ERR_SUP,"-da_partition_nodes_at_end not supported");
388: } else { /* Normal PETSc distribution */
389: y = N/n + ((N % n) > (rank/m));
390: if (y < s) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row width is too thin for stencil! %d %d",y,s);
391: if ((N % n) > (rank/m)) { ys = (rank/m)*y; }
392: else { ys = (N % n)*(y+1) + ((rank/m)-(N % n))*y; }
393: PetscMalloc(n*sizeof(int),&ly);
394: fly = ly;
395: for (i=0; i<n; i++) {
396: ly[i] = N/n + ((N % n) > i);
397: }
398: }
400: xe = xs + x;
401: ye = ys + y;
403: /* determine ghost region */
404: /* Assume No Periodicity */
405: if (xs-s > 0) Xs = xs - s; else Xs = 0;
406: if (ys-s > 0) Ys = ys - s; else Ys = 0;
407: if (xe+s <= M) Xe = xe + s; else Xe = M;
408: if (ye+s <= N) Ye = ye + s; else Ye = N;
410: /* X Periodic */
411: if (DAXPeriodic(wrap)){
412: Xs = xs - s;
413: Xe = xe + s;
414: }
416: /* Y Periodic */
417: if (DAYPeriodic(wrap)){
418: Ys = ys - s;
419: Ye = ye + s;
420: }
422: /* Resize all X parameters to reflect w */
423: x *= dof;
424: xs *= dof;
425: xe *= dof;
426: Xs *= dof;
427: Xe *= dof;
428: s_x = s*dof;
429: s_y = s;
431: /* determine starting point of each processor */
432: nn = x*y;
433: PetscMalloc((2*size+1)*sizeof(int),&bases);
434: ldims = (int*)(bases+size+1);
435: MPI_Allgather(&nn,1,MPI_INT,ldims,1,MPI_INT,comm);
436: bases[0] = 0;
437: for (i=1; i<=size; i++) {
438: bases[i] = ldims[i-1];
439: }
440: for (i=1; i<=size; i++) {
441: bases[i] += bases[i-1];
442: }
444: /* allocate the base parallel and sequential vectors */
445: da->Nlocal = x*y;
446: VecCreateMPIWithArray(comm,da->Nlocal,PETSC_DECIDE,0,&global);
447: VecSetBlockSize(global,dof);
448: da->nlocal = (Xe-Xs)*(Ye-Ys) ;
449: VecCreateSeqWithArray(PETSC_COMM_SELF,da->nlocal,0,&local);
450: VecSetBlockSize(local,dof);
453: /* generate appropriate vector scatters */
454: /* local to global inserts non-ghost point region into global */
455: VecGetOwnershipRange(global,&start,&end);
456: ISCreateStride(comm,x*y,start,1,&to);
458: left = xs - Xs; down = ys - Ys; up = down + y;
459: PetscMalloc(x*(up - down)*sizeof(int),&idx);
460: count = 0;
461: for (i=down; i<up; i++) {
462: for (j=0; j<x; j++) {
463: idx[count++] = left + i*(Xe-Xs) + j;
464: }
465: }
466: ISCreateGeneral(comm,count,idx,&from);
467: PetscFree(idx);
469: VecScatterCreate(local,from,global,to,<og);
470: PetscLogObjectParent(da,to);
471: PetscLogObjectParent(da,from);
472: PetscLogObjectParent(da,ltog);
473: ISDestroy(from);
474: ISDestroy(to);
476: /* global to local must include ghost points */
477: if (stencil_type == DA_STENCIL_BOX) {
478: ISCreateStride(comm,(Xe-Xs)*(Ye-Ys),0,1,&to);
479: } else {
480: /* must drop into cross shape region */
481: /* ---------|
482: | top |
483: |--- ---|
484: | middle |
485: | |
486: ---- ----
487: | bottom |
488: -----------
489: Xs xs xe Xe */
490: /* bottom */
491: left = xs - Xs; down = ys - Ys; up = down + y;
492: count = down*(xe-xs) + (up-down)*(Xe-Xs) + (Ye-Ys-up)*(xe-xs);
493: ierr = PetscMalloc(count*sizeof(int),&idx);
494: count = 0;
495: for (i=0; i<down; i++) {
496: for (j=0; j<xe-xs; j++) {
497: idx[count++] = left + i*(Xe-Xs) + j;
498: }
499: }
500: /* middle */
501: for (i=down; i<up; i++) {
502: for (j=0; j<Xe-Xs; j++) {
503: idx[count++] = i*(Xe-Xs) + j;
504: }
505: }
506: /* top */
507: for (i=up; i<Ye-Ys; i++) {
508: for (j=0; j<xe-xs; j++) {
509: idx[count++] = left + i*(Xe-Xs) + j;
510: }
511: }
512: ISCreateGeneral(comm,count,idx,&to);
513: PetscFree(idx);
514: }
517: /* determine who lies on each side of us stored in n6 n7 n8
518: n3 n5
519: n0 n1 n2
520: */
522: /* Assume the Non-Periodic Case */
523: n1 = rank - m;
524: if (rank % m) {
525: n0 = n1 - 1;
526: } else {
527: n0 = -1;
528: }
529: if ((rank+1) % m) {
530: n2 = n1 + 1;
531: n5 = rank + 1;
532: n8 = rank + m + 1; if (n8 >= m*n) n8 = -1;
533: } else {
534: n2 = -1; n5 = -1; n8 = -1;
535: }
536: if (rank % m) {
537: n3 = rank - 1;
538: n6 = n3 + m; if (n6 >= m*n) n6 = -1;
539: } else {
540: n3 = -1; n6 = -1;
541: }
542: n7 = rank + m; if (n7 >= m*n) n7 = -1;
545: /* Modify for Periodic Cases */
546: if (wrap == DA_YPERIODIC) { /* Handle Top and Bottom Sides */
547: if (n1 < 0) n1 = rank + m * (n-1);
548: if (n7 < 0) n7 = rank - m * (n-1);
549: if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1;
550: if ((n3 >= 0) && (n6 < 0)) n6 = (rank%m)-1;
551: if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1;
552: if ((n5 >= 0) && (n8 < 0)) n8 = (rank%m)+1;
553: } else if (wrap == DA_XPERIODIC) { /* Handle Left and Right Sides */
554: if (n3 < 0) n3 = rank + (m-1);
555: if (n5 < 0) n5 = rank - (m-1);
556: if ((n1 >= 0) && (n0 < 0)) n0 = rank-1;
557: if ((n1 >= 0) && (n2 < 0)) n2 = rank-2*m+1;
558: if ((n7 >= 0) && (n6 < 0)) n6 = rank+2*m-1;
559: if ((n7 >= 0) && (n8 < 0)) n8 = rank+1;
560: } else if (wrap == DA_XYPERIODIC) {
562: /* Handle all four corners */
563: if ((n6 < 0) && (n7 < 0) && (n3 < 0)) n6 = m-1;
564: if ((n8 < 0) && (n7 < 0) && (n5 < 0)) n8 = 0;
565: if ((n2 < 0) && (n5 < 0) && (n1 < 0)) n2 = size-m;
566: if ((n0 < 0) && (n3 < 0) && (n1 < 0)) n0 = size-1;
568: /* Handle Top and Bottom Sides */
569: if (n1 < 0) n1 = rank + m * (n-1);
570: if (n7 < 0) n7 = rank - m * (n-1);
571: if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1;
572: if ((n3 >= 0) && (n6 < 0)) n6 = (rank%m)-1;
573: if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1;
574: if ((n5 >= 0) && (n8 < 0)) n8 = (rank%m)+1;
576: /* Handle Left and Right Sides */
577: if (n3 < 0) n3 = rank + (m-1);
578: if (n5 < 0) n5 = rank - (m-1);
579: if ((n1 >= 0) && (n0 < 0)) n0 = rank-1;
580: if ((n1 >= 0) && (n2 < 0)) n2 = rank-2*m+1;
581: if ((n7 >= 0) && (n6 < 0)) n6 = rank+2*m-1;
582: if ((n7 >= 0) && (n8 < 0)) n8 = rank+1;
583: }
585: if (stencil_type == DA_STENCIL_STAR) {
586: /* save corner processor numbers */
587: sn0 = n0; sn2 = n2; sn6 = n6; sn8 = n8;
588: n0 = n2 = n6 = n8 = -1;
589: }
591: PetscMalloc((x+2*s_x)*(y+2*s_y)*sizeof(int),&idx);
592: PetscLogObjectMemory(da,(x+2*s_x)*(y+2*s_y)*sizeof(int));
593: nn = 0;
595: xbase = bases[rank];
596: for (i=1; i<=s_y; i++) {
597: if (n0 >= 0) { /* left below */
598: x_t = lx[n0 % m]*dof;
599: y_t = ly[(n0/m)];
600: s_t = bases[n0] + x_t*y_t - (s_y-i)*x_t - s_x;
601: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
602: }
603: if (n1 >= 0) { /* directly below */
604: x_t = x;
605: y_t = ly[(n1/m)];
606: s_t = bases[n1] + x_t*y_t - (s_y+1-i)*x_t;
607: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
608: }
609: if (n2 >= 0) { /* right below */
610: x_t = lx[n2 % m]*dof;
611: y_t = ly[(n2/m)];
612: s_t = bases[n2] + x_t*y_t - (s_y+1-i)*x_t;
613: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
614: }
615: }
617: for (i=0; i<y; i++) {
618: if (n3 >= 0) { /* directly left */
619: x_t = lx[n3 % m]*dof;
620: /* y_t = y; */
621: s_t = bases[n3] + (i+1)*x_t - s_x;
622: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
623: }
625: for (j=0; j<x; j++) { idx[nn++] = xbase++; } /* interior */
627: if (n5 >= 0) { /* directly right */
628: x_t = lx[n5 % m]*dof;
629: /* y_t = y; */
630: s_t = bases[n5] + (i)*x_t;
631: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
632: }
633: }
635: for (i=1; i<=s_y; i++) {
636: if (n6 >= 0) { /* left above */
637: x_t = lx[n6 % m]*dof;
638: /* y_t = ly[(n6/m)]; */
639: s_t = bases[n6] + (i)*x_t - s_x;
640: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
641: }
642: if (n7 >= 0) { /* directly above */
643: x_t = x;
644: /* y_t = ly[(n7/m)]; */
645: s_t = bases[n7] + (i-1)*x_t;
646: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
647: }
648: if (n8 >= 0) { /* right above */
649: x_t = lx[n8 % m]*dof;
650: /* y_t = ly[(n8/m)]; */
651: s_t = bases[n8] + (i-1)*x_t;
652: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
653: }
654: }
656: base = bases[rank];
657: ISCreateGeneral(comm,nn,idx,&from);
658: VecScatterCreate(global,from,local,to,>ol);
659: PetscLogObjectParent(da,to);
660: PetscLogObjectParent(da,from);
661: PetscLogObjectParent(da,gtol);
662: ISDestroy(to);
663: ISDestroy(from);
665: if (stencil_type == DA_STENCIL_STAR) {
666: /*
667: Recompute the local to global mappings, this time keeping the
668: information about the cross corner processor numbers.
669: */
670: n0 = sn0; n2 = sn2; n6 = sn6; n8 = sn8;
671: nn = 0;
672: xbase = bases[rank];
673: for (i=1; i<=s_y; i++) {
674: if (n0 >= 0) { /* left below */
675: x_t = lx[n0 % m]*dof;
676: y_t = ly[(n0/m)];
677: s_t = bases[n0] + x_t*y_t - (s_y-i)*x_t - s_x;
678: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
679: }
680: if (n1 >= 0) { /* directly below */
681: x_t = x;
682: y_t = ly[(n1/m)];
683: s_t = bases[n1] + x_t*y_t - (s_y+1-i)*x_t;
684: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
685: }
686: if (n2 >= 0) { /* right below */
687: x_t = lx[n2 % m]*dof;
688: y_t = ly[(n2/m)];
689: s_t = bases[n2] + x_t*y_t - (s_y+1-i)*x_t;
690: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
691: }
692: }
694: for (i=0; i<y; i++) {
695: if (n3 >= 0) { /* directly left */
696: x_t = lx[n3 % m]*dof;
697: /* y_t = y; */
698: s_t = bases[n3] + (i+1)*x_t - s_x;
699: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
700: }
702: for (j=0; j<x; j++) { idx[nn++] = xbase++; } /* interior */
704: if (n5 >= 0) { /* directly right */
705: x_t = lx[n5 % m]*dof;
706: /* y_t = y; */
707: s_t = bases[n5] + (i)*x_t;
708: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
709: }
710: }
712: for (i=1; i<=s_y; i++) {
713: if (n6 >= 0) { /* left above */
714: x_t = lx[n6 % m]*dof;
715: /* y_t = ly[(n6/m)]; */
716: s_t = bases[n6] + (i)*x_t - s_x;
717: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
718: }
719: if (n7 >= 0) { /* directly above */
720: x_t = x;
721: /* y_t = ly[(n7/m)]; */
722: s_t = bases[n7] + (i-1)*x_t;
723: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
724: }
725: if (n8 >= 0) { /* right above */
726: x_t = lx[n8 % m]*dof;
727: /* y_t = ly[(n8/m)]; */
728: s_t = bases[n8] + (i-1)*x_t;
729: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
730: }
731: }
732: }
733: PetscFree(bases);
735: da->M = M; da->N = N; da->m = m; da->n = n; da->w = dof; da->s = s;
736: da->xs = xs; da->xe = xe; da->ys = ys; da->ye = ye; da->zs = 0; da->ze = 1;
737: da->Xs = Xs; da->Xe = Xe; da->Ys = Ys; da->Ye = Ye; da->Zs = 0; da->Ze = 1;
738: da->P = 1; da->p = 1;
740: VecDestroy(local);
741: VecDestroy(global);
743: da->gtol = gtol;
744: da->ltog = ltog;
745: da->idx = idx;
746: da->Nl = nn;
747: da->base = base;
748: da->wrap = wrap;
749: da->ops->view = DAView_2d;
750: da->stencil_type = stencil_type;
752: /*
753: Set the local to global ordering in the global vector, this allows use
754: of VecSetValuesLocal().
755: */
756: ISLocalToGlobalMappingCreateNC(comm,nn,idx,&da->ltogmap);
757: ISLocalToGlobalMappingBlock(da->ltogmap,da->w,&da->ltogmapb);
758: PetscLogObjectParent(da,da->ltogmap);
760: *inra = da;
762: da->ltol = PETSC_NULL;
763: da->ao = PETSC_NULL;
766: if (!flx) {
767: PetscMalloc(m*sizeof(int),&flx);
768: PetscMemcpy(flx,lx,m*sizeof(int));
769: }
770: if (!fly) {
771: PetscMalloc(n*sizeof(int),&fly);
772: PetscMemcpy(fly,ly,n*sizeof(int));
773: }
774: da->lx = flx;
775: da->ly = fly;
777: PetscOptionsHasName(PETSC_NULL,"-da_view",&flg1);
778: if (flg1) {DAView(da,PETSC_VIEWER_STDOUT_(da->comm));}
779: PetscOptionsHasName(PETSC_NULL,"-da_view_draw",&flg1);
780: if (flg1) {DAView(da,PETSC_VIEWER_DRAW_(da->comm));}
781: PetscOptionsHasName(PETSC_NULL,"-help",&flg1);
782: if (flg1) {DAPrintHelp(da);}
784: PetscPublishAll(da);
785: return(0);
786: }
788: #undef __FUNCT__
790: /*@
791: DAPrintHelp - Prints command line options for DA.
793: Collective on DA
795: Input Parameters:
796: . da - the distributed array
798: Level: intermediate
800: .seealso: DACreate1d(), DACreate2d(), DACreate3d()
802: .keywords: DA, help
804: @*/
805: int DAPrintHelp(DA da)
806: {
807: static PetscTruth called = PETSC_FALSE;
808: MPI_Comm comm;
809: int ierr;
814: comm = da->comm;
815: if (!called) {
816: (*PetscHelpPrintf)(comm,"General Distributed Array (DA) options:n");
817: (*PetscHelpPrintf)(comm," -da_view: print DA distribution to screenn");
818: (*PetscHelpPrintf)(comm," -da_view_draw: display DA in windown");
819: called = PETSC_TRUE;
820: }
821: return(0);
822: }
824: #undef __FUNCT__
826: /*@C
827: DARefine - Creates a new distributed array that is a refinement of a given
828: distributed array.
830: Collective on DA
832: Input Parameter:
833: + da - initial distributed array
834: - comm - communicator to contain refined DA, must be either same as the da communicator or include the
835: da communicator and be 2, 4, or 8 times larger. Currently ignored
837: Output Parameter:
838: . daref - refined distributed array
840: Level: advanced
842: Note:
843: Currently, refinement consists of just doubling the number of grid spaces
844: in each dimension of the DA.
846: .keywords: distributed array, refine
848: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy()
849: @*/
850: int DARefine(DA da,MPI_Comm comm,DA *daref)
851: {
852: int M,N,P,ierr;
853: DA da2;
858: if (DAXPeriodic(da->wrap) || da->interptype == DA_Q0){
859: M = da->refine_x*da->M;
860: } else {
861: M = 1 + da->refine_x*(da->M - 1);
862: }
863: if (DAYPeriodic(da->wrap) || da->interptype == DA_Q0){
864: N = da->refine_y*da->N;
865: } else {
866: N = 1 + da->refine_y*(da->N - 1);
867: }
868: if (DAZPeriodic(da->wrap) || da->interptype == DA_Q0){
869: P = da->refine_z*da->P;
870: } else {
871: P = 1 + da->refine_z*(da->P - 1);
872: }
873: if (da->dim == 1) {
874: DACreate1d(da->comm,da->wrap,M,da->w,da->s,PETSC_NULL,&da2);
875: } else if (da->dim == 2) {
876: DACreate2d(da->comm,da->wrap,da->stencil_type,M,N,da->m,da->n,da->w,da->s,PETSC_NULL,PETSC_NULL,&da2);
877: } else if (da->dim == 3) {
878: DACreate3d(da->comm,da->wrap,da->stencil_type,M,N,P,da->m,da->n,da->p,da->w,da->s,0,0,0,&da2);
879: }
880: *daref = da2;
881: return(0);
882: }
884: /*
885: M is number of grid points
886: m is number of processors
888: */
889: #undef __FUNCT__
891: int DASplitComm2d(MPI_Comm comm,int M,int N,int sw,MPI_Comm *outcomm)
892: {
893: int ierr,m,n = 0,csize,size,rank,x = 0,y = 0;
896: MPI_Comm_size(comm,&size);
897: MPI_Comm_rank(comm,&rank);
899: csize = 4*size;
900: do {
901: if (csize % 4) SETERRQ4(1,"Cannot split communicator of size %d tried %d %d %d",size,csize,x,y);
902: csize = csize/4;
903:
904: m = (int)(0.5 + sqrt(((double)M)*((double)csize)/((double)N)));
905: if (!m) m = 1;
906: while (m > 0) {
907: n = csize/m;
908: if (m*n == csize) break;
909: m--;
910: }
911: if (M > N && m < n) {int _m = m; m = n; n = _m;}
913: x = M/m + ((M % m) > ((csize-1) % m));
914: y = (N + (csize-1)/m)/n;
915: } while ((x < 4 || y < 4) && csize > 1);
916: if (size != csize) {
917: MPI_Group entire_group,sub_group;
918: int i,*groupies;
920: ierr = MPI_Comm_group(comm,&entire_group);
921: PetscMalloc(csize*sizeof(int),&groupies);
922: for (i=0; i<csize; i++) {
923: groupies[i] = (rank/csize)*csize + i;
924: }
925: ierr = MPI_Group_incl(entire_group,csize,groupies,&sub_group);
926: ierr = PetscFree(groupies);
927: ierr = MPI_Comm_create(comm,sub_group,outcomm);
928: ierr = MPI_Group_free(&entire_group);
929: ierr = MPI_Group_free(&sub_group);
930: PetscLogInfo(0,"Creating redundant coarse problems of size %dn",csize);
931: } else {
932: *outcomm = comm;
933: }
934: return(0);
935: }
937: #undef __FUNCT__
939: /*@C
940: DASetLocalFunction - Caches in a DA a local function.
942: Collective on DA
944: Input Parameter:
945: + da - initial distributed array
946: - lf - the local function
948: Level: intermediate
950: Notes: The routine SNESDAFormFunction() uses this the cached function to evaluate the user provided function.
952: .keywords: distributed array, refine
954: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy(), DAGetLocalFunction(), DASetLocalFunctioni()
955: @*/
956: int DASetLocalFunction(DA da,DALocalFunction1 lf)
957: {
960: da->lf = lf;
961: return(0);
962: }
964: #undef __FUNCT__
966: /*@C
967: DASetLocalFunctioni - Caches in a DA a local function that evaluates a single component
969: Collective on DA
971: Input Parameter:
972: + da - initial distributed array
973: - lfi - the local function
975: Level: intermediate
977: .keywords: distributed array, refine
979: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy(), DAGetLocalFunction(), DASetLocalFunction()
980: @*/
981: int DASetLocalFunctioni(DA da,int (*lfi)(DALocalInfo*,MatStencil*,void*,PetscScalar*,void*))
982: {
985: da->lfi = lfi;
986: return(0);
987: }
989: /*MC
990: DASetLocalAdicFunction - Caches in a DA a local function computed by ADIC/ADIFOR
992: Collective on DA
994: Synopsis:
995: int int DASetLocalAdicFunction(DA da,DALocalFunction1 ad_lf)
996:
997: Input Parameter:
998: + da - initial distributed array
999: - ad_lf - the local function as computed by ADIC/ADIFOR
1001: Level: intermediate
1003: .keywords: distributed array, refine
1005: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy(), DAGetLocalFunction(), DASetLocalFunction(),
1006: DASetLocalJacobian()
1007: M*/
1009: #undef __FUNCT__
1011: int DASetLocalAdicFunction_Private(DA da,DALocalFunction1 ad_lf)
1012: {
1015: da->adic_lf = ad_lf;
1016: return(0);
1017: }
1019: /*MC
1020: DASetLocalAdicFunctioni - Caches in a DA a local functioni computed by ADIC/ADIFOR
1022: Collective on DA
1024: Synopsis:
1025: int int DASetLocalAdicFunctioni(DA da,int (ad_lf*)(DALocalInfo*,MatStencil*,void*,void*,void*)
1026:
1027: Input Parameter:
1028: + da - initial distributed array
1029: - ad_lfi - the local function as computed by ADIC/ADIFOR
1031: Level: intermediate
1033: .keywords: distributed array, refine
1035: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy(), DAGetLocalFunction(), DASetLocalFunction(),
1036: DASetLocalJacobian(), DASetLocalFunctioni()
1037: M*/
1039: #undef __FUNCT__
1041: int DASetLocalAdicFunctioni_Private(DA da,int (*ad_lfi)(DALocalInfo*,MatStencil*,void*,void*,void*))
1042: {
1045: da->adic_lfi = ad_lfi;
1046: return(0);
1047: }
1049: /*MC
1050: DASetLocalAdicMFFunctioni - Caches in a DA a local functioni computed by ADIC/ADIFOR
1052: Collective on DA
1054: Synopsis:
1055: int int DASetLocalAdicFunctioni(DA da,int (ad_lf*)(DALocalInfo*,MatStencil*,void*,void*,void*)
1056:
1057: Input Parameter:
1058: + da - initial distributed array
1059: - admf_lfi - the local matrix-free function as computed by ADIC/ADIFOR
1061: Level: intermediate
1063: .keywords: distributed array, refine
1065: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy(), DAGetLocalFunction(), DASetLocalFunction(),
1066: DASetLocalJacobian(), DASetLocalFunctioni()
1067: M*/
1069: #undef __FUNCT__
1071: int DASetLocalAdicMFFunctioni_Private(DA da,int (*admf_lfi)(DALocalInfo*,MatStencil*,void*,void*,void*))
1072: {
1075: da->adicmf_lfi = admf_lfi;
1076: return(0);
1077: }
1079: /*MC
1080: DASetLocalAdicMFFunction - Caches in a DA a local function computed by ADIC/ADIFOR
1082: Collective on DA
1084: Synopsis:
1085: int int DASetLocalAdicMFFunction(DA da,DALocalFunction1 ad_lf)
1086:
1087: Input Parameter:
1088: + da - initial distributed array
1089: - ad_lf - the local function as computed by ADIC/ADIFOR
1091: Level: intermediate
1093: .keywords: distributed array, refine
1095: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy(), DAGetLocalFunction(), DASetLocalFunction(),
1096: DASetLocalJacobian()
1097: M*/
1099: #undef __FUNCT__
1101: int DASetLocalAdicMFFunction_Private(DA da,DALocalFunction1 ad_lf)
1102: {
1105: da->adicmf_lf = ad_lf;
1106: return(0);
1107: }
1109: /*@C
1110: DASetLocalJacobian - Caches in a DA a local Jacobian
1112: Collective on DA
1114:
1115: Input Parameter:
1116: + da - initial distributed array
1117: - lj - the local Jacobian
1119: Level: intermediate
1121: Notes: The routine SNESDAFormFunction() uses this the cached function to evaluate the user provided function.
1123: .keywords: distributed array, refine
1125: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy(), DAGetLocalFunction(), DASetLocalFunction()
1126: @*/
1127: #undef __FUNCT__
1129: int DASetLocalJacobian(DA da,DALocalFunction1 lj)
1130: {
1133: da->lj = lj;
1134: return(0);
1135: }
1137: #undef __FUNCT__
1139: /*@C
1140: DAGetLocalFunction - Gets from a DA a local function and its ADIC/ADIFOR Jacobian
1142: Collective on DA
1144: Input Parameter:
1145: . da - initial distributed array
1147: Output Parameters:
1148: . lf - the local function
1150: Level: intermediate
1152: .keywords: distributed array, refine
1154: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy(), DASetLocalFunction()
1155: @*/
1156: int DAGetLocalFunction(DA da,DALocalFunction1 *lf)
1157: {
1160: if (lf) *lf = da->lf;
1161: return(0);
1162: }
1164: #undef __FUNCT__
1166: /*@
1167: DAFormFunction1 - Evaluates a user provided function on each processor that
1168: share a DA
1170: Input Parameters:
1171: + da - the DA that defines the grid
1172: . vu - input vector
1173: . vfu - output vector
1174: - w - any user data
1176: Notes: Does NOT do ghost updates on vu upon entry
1178: Level: advanced
1180: .seealso: DAComputeJacobian1WithAdic()
1182: @*/
1183: int DAFormFunction1(DA da,Vec vu,Vec vfu,void *w)
1184: {
1185: int ierr;
1186: void *u,*fu;
1187: DALocalInfo info;
1188:
1191: DAGetLocalInfo(da,&info);
1192: DAVecGetArray(da,vu,(void**)&u);
1193: DAVecGetArray(da,vfu,(void**)&fu);
1195: (*da->lf)(&info,u,fu,w);
1197: DAVecRestoreArray(da,vu,(void**)&u);
1198: DAVecRestoreArray(da,vfu,(void**)&fu);
1199: return(0);
1200: }
1202: #undef __FUNCT__
1204: int DAFormFunctioniTest1(DA da,void *w)
1205: {
1206: Vec vu,fu,fui;
1207: int ierr,i,n;
1208: PetscScalar *ui,mone = -1.0;
1209: PetscRandom rnd;
1210: PetscReal norm;
1213: DAGetLocalVector(da,&vu);
1214: PetscRandomCreate(PETSC_COMM_SELF,RANDOM_DEFAULT,&rnd);
1215: VecSetRandom(rnd,vu);
1216: PetscRandomDestroy(rnd);
1218: DAGetGlobalVector(da,&fu);
1219: DAGetGlobalVector(da,&fui);
1220:
1221: DAFormFunction1(da,vu,fu,w);
1223: VecGetArray(fui,&ui);
1224: VecGetLocalSize(fui,&n);
1225: for (i=0; i<n; i++) {
1226: DAFormFunctioni1(da,i,vu,ui+i,w);
1227: }
1228: VecRestoreArray(fui,&ui);
1230: VecAXPY(&mone,fu,fui);
1231: VecNorm(fui,NORM_2,&norm);
1232: PetscPrintf(da->comm,"Norm of difference in vectors %gn",norm);
1233: VecView(fu,0);
1234: VecView(fui,0);
1236: DARestoreLocalVector(da,&vu);
1237: DARestoreGlobalVector(da,&fu);
1238: DARestoreGlobalVector(da,&fui);
1239: return(0);
1240: }
1242: #undef __FUNCT__
1244: /*@
1245: DAFormFunctioni1 - Evaluates a user provided function
1247: Input Parameters:
1248: + da - the DA that defines the grid
1249: . i - the component of the function we wish to compute (must be local)
1250: . vu - input vector
1251: . vfu - output value
1252: - w - any user data
1254: Notes: Does NOT do ghost updates on vu upon entry
1256: Level: advanced
1258: .seealso: DAComputeJacobian1WithAdic()
1260: @*/
1261: int DAFormFunctioni1(DA da,int i,Vec vu,PetscScalar *vfu,void *w)
1262: {
1263: int ierr;
1264: void *u;
1265: DALocalInfo info;
1266: MatStencil stencil;
1267:
1270: DAGetLocalInfo(da,&info);
1271: DAVecGetArray(da,vu,(void**)&u);
1273: /* figure out stencil value from i */
1274: stencil.c = i % info.dof;
1275: stencil.i = (i % (info.xm*info.dof))/info.dof;
1276: stencil.j = (i % (info.xm*info.ym*info.dof))/(info.xm*info.dof);
1277: stencil.k = i/(info.xm*info.ym*info.dof);
1279: (*da->lfi)(&info,&stencil,u,vfu,w);
1281: DAVecRestoreArray(da,vu,(void**)&u);
1282: return(0);
1283: }
1285: #if defined(new)
1286: #undef __FUNCT__
1288: /*
1289: DAGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix where local
1290: function lives on a DA
1292: y ~= (F(u + ha) - F(u))/h,
1293: where F = nonlinear function, as set by SNESSetFunction()
1294: u = current iterate
1295: h = difference interval
1296: */
1297: int DAGetDiagonal_MFFD(DA da,Vec U,Vec a)
1298: {
1299: PetscScalar h,*aa,*ww,v;
1300: PetscReal epsilon = PETSC_SQRT_MACHINE_EPSILON,umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
1301: int ierr,gI,nI;
1302: MatStencil stencil;
1303: DALocalInfo info;
1304:
1306: (*ctx->func)(0,U,a,ctx->funcctx);
1307: (*ctx->funcisetbase)(U,ctx->funcctx);
1309: VecGetArray(U,&ww);
1310: VecGetArray(a,&aa);
1311:
1312: nI = 0;
1313: h = ww[gI];
1314: if (h == 0.0) h = 1.0;
1315: #if !defined(PETSC_USE_COMPLEX)
1316: if (h < umin && h >= 0.0) h = umin;
1317: else if (h < 0.0 && h > -umin) h = -umin;
1318: #else
1319: if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0) h = umin;
1320: else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
1321: #endif
1322: h *= epsilon;
1323:
1324: ww[gI += h;
1325: ierr = (*ctx->funci)(i,w,&v,ctx->funcctx);
1326: aa[nI] = (v - aa[nI])/h;
1327: ww[gI] -= h;
1328: nI++;
1329: }
1330: VecRestoreArray(U,&ww);
1331: VecRestoreArray(a,&aa);
1332: return(0);
1333: }
1334: #endif
1336: #if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE)
1337: EXTERN_C_BEGIN
1338: #include "adic/ad_utils.h"
1339: EXTERN_C_END
1341: #undef __FUNCT__
1343: /*@
1344: DAComputeJacobian1WithAdic - Evaluates a adiC provided Jacobian function on each processor that
1345: share a DA
1347: Input Parameters:
1348: + da - the DA that defines the grid
1349: . vu - input vector (ghosted)
1350: . J - output matrix
1351: - w - any user data
1353: Level: advanced
1355: Notes: Does NOT do ghost updates on vu upon entry
1357: .seealso: DAFormFunction1()
1359: @*/
1360: int DAComputeJacobian1WithAdic(DA da,Vec vu,Mat J,void *w)
1361: {
1362: int ierr,gtdof,tdof;
1363: PetscScalar *ustart;
1364: DALocalInfo info;
1365: void *ad_u,*ad_f,*ad_ustart,*ad_fstart;
1366: ISColoring iscoloring;
1369: DAGetLocalInfo(da,&info);
1371: /* get space for derivative objects. */
1372: DAGetAdicArray(da,PETSC_TRUE,(void **)&ad_u,&ad_ustart,>dof);
1373: DAGetAdicArray(da,PETSC_FALSE,(void **)&ad_f,&ad_fstart,&tdof);
1374: VecGetArray(vu,&ustart);
1375: PetscADSetValArray(((DERIV_TYPE*)ad_ustart),gtdof,ustart);
1376: VecRestoreArray(vu,&ustart);
1378: PetscADResetIndep();
1379: DAGetColoring(da,IS_COLORING_GHOSTED,&iscoloring);
1380: PetscADSetIndepArrayColored(ad_ustart,gtdof,iscoloring->colors);
1381: PetscADIncrementTotalGradSize(iscoloring->n);
1382: ISColoringDestroy(iscoloring);
1383: PetscADSetIndepDone();
1385: (*da->adic_lf)(&info,ad_u,ad_f,w);
1387: /* stick the values into the matrix */
1388: MatSetValuesAdic(J,(PetscScalar**)ad_fstart);
1390: /* return space for derivative objects. */
1391: DARestoreAdicArray(da,PETSC_TRUE,(void **)&ad_u,&ad_ustart,>dof);
1392: DARestoreAdicArray(da,PETSC_FALSE,(void **)&ad_f,&ad_fstart,&tdof);
1393: return(0);
1394: }
1396: #undef __FUNCT__
1398: /*@C
1399: DAMultiplyByJacobian1WithAdic - Applies an ADIC-provided Jacobian function to a vector on
1400: each processor that shares a DA.
1402: Input Parameters:
1403: + da - the DA that defines the grid
1404: . vu - Jacobian is computed at this point (ghosted)
1405: . v - product is done on this vector (ghosted)
1406: . fu - output vector = J(vu)*v (not ghosted)
1407: - w - any user data
1409: Notes:
1410: This routine does NOT do ghost updates on vu upon entry.
1412: Level: advanced
1414: .seealso: DAFormFunction1()
1416: @*/
1417: int DAMultiplyByJacobian1WithAdic(DA da,Vec vu,Vec v,Vec f,void *w)
1418: {
1419: int ierr,i,gtdof,tdof;
1420: PetscScalar *avu,*av,*af,*ad_vustart,*ad_fstart;
1421: DALocalInfo info;
1422: void *ad_vu,*ad_f;
1425: DAGetLocalInfo(da,&info);
1427: /* get space for derivative objects. */
1428: DAGetAdicMFArray(da,PETSC_TRUE,(void **)&ad_vu,(void**)&ad_vustart,>dof);
1429: DAGetAdicMFArray(da,PETSC_FALSE,(void **)&ad_f,(void**)&ad_fstart,&tdof);
1431: /* copy input vector into derivative object */
1432: VecGetArray(vu,&avu);
1433: VecGetArray(v,&av);
1434: for (i=0; i<gtdof; i++) {
1435: ad_vustart[2*i] = avu[i];
1436: ad_vustart[2*i+1] = av[i];
1437: }
1438: VecRestoreArray(vu,&avu);
1439: VecRestoreArray(v,&av);
1441: PetscADResetIndep();
1442: PetscADIncrementTotalGradSize(1);
1443: PetscADSetIndepDone();
1445: (*da->adicmf_lf)(&info,ad_vu,ad_f,w);
1447: /* stick the values into the vector */
1448: VecGetArray(f,&af);
1449: for (i=0; i<tdof; i++) {
1450: af[i] = ad_fstart[2*i+1];
1451: }
1452: VecRestoreArray(f,&af);
1454: /* return space for derivative objects. */
1455: DARestoreAdicMFArray(da,PETSC_TRUE,(void **)&ad_vu,(void**)&ad_vustart,>dof);
1456: DARestoreAdicMFArray(da,PETSC_FALSE,(void **)&ad_f,(void**)&ad_fstart,&tdof);
1457: return(0);
1458: }
1461: #else
1463: #undef __FUNCT__
1465: int DAComputeJacobian1WithAdic(DA da,Vec vu,Mat J,void *w)
1466: {
1468: SETERRQ(1,"Must compile with bmake/PETSC_ARCH/packages flag PETSC_HAVE_ADIC for this routine");
1469: }
1471: #undef __FUNCT__
1473: int DAMultiplyByJacobian1WithAdic(DA da,Vec vu,Vec v,Vec f,void *w)
1474: {
1476: SETERRQ(1,"Must compile with bmake/PETSC_ARCH/packages flag PETSC_HAVE_ADIC for this routine");
1477: }
1479: #endif
1481: #undef __FUNCT__
1483: /*@
1484: DAComputeJacobian1 - Evaluates a local Jacobian function on each processor that
1485: share a DA
1487: Input Parameters:
1488: + da - the DA that defines the grid
1489: . vu - input vector (ghosted)
1490: . J - output matrix
1491: - w - any user data
1493: Notes: Does NOT do ghost updates on vu upon entry
1495: Level: advanced
1497: .seealso: DAFormFunction1()
1499: @*/
1500: int DAComputeJacobian1(DA da,Vec vu,Mat J,void *w)
1501: {
1502: int ierr;
1503: void *u;
1504: DALocalInfo info;
1507: DAGetLocalInfo(da,&info);
1508: DAVecGetArray(da,vu,&u);
1509: (*da->lj)(&info,u,J,w);
1510: DAVecRestoreArray(da,vu,&u);
1511: return(0);
1512: }
1515: #undef __FUNCT__
1517: /*
1518: DAComputeJacobian1WithAdifor - Evaluates a ADIFOR provided Jacobian local function on each processor that
1519: share a DA
1521: Input Parameters:
1522: + da - the DA that defines the grid
1523: . vu - input vector (ghosted)
1524: . J - output matrix
1525: - w - any user data
1527: Notes: Does NOT do ghost updates on vu upon entry
1529: .seealso: DAFormFunction1()
1531: */
1532: int DAComputeJacobian1WithAdifor(DA da,Vec vu,Mat J,void *w)
1533: {
1534: int i,ierr,Nc,*color,N;
1535: DALocalInfo info;
1536: PetscScalar *u,*g_u,*g_f,*f,*p_u;
1537: ISColoring iscoloring;
1538: void (*lf)(int *,DALocalInfo*,PetscScalar*,PetscScalar*,int*,PetscScalar*,PetscScalar*,int*,void*,int*) =
1539: (void (*)(int *,DALocalInfo*,PetscScalar*,PetscScalar*,int*,PetscScalar*,PetscScalar*,int*,void*,int*))*da->adifor_lf;
1542: DAGetColoring(da,IS_COLORING_GHOSTED,&iscoloring);
1543: Nc = iscoloring->n;
1544: DAGetLocalInfo(da,&info);
1545: N = info.gxm*info.gym*info.gzm*info.dof;
1547: /* get space for derivative objects. */
1548: ierr = PetscMalloc(Nc*info.gxm*info.gym*info.gzm*info.dof*sizeof(PetscScalar),&g_u);
1549: ierr = PetscMemzero(g_u,Nc*info.gxm*info.gym*info.gzm*info.dof*sizeof(PetscScalar));
1550: p_u = g_u;
1551: color = iscoloring->colors;
1552: for (i=0; i<N; i++) {
1553: p_u[*color++] = 1.0;
1554: p_u += Nc;
1555: }
1556: ISColoringDestroy(iscoloring);
1557: PetscMalloc(Nc*info.xm*info.ym*info.zm*info.dof*sizeof(PetscScalar),&g_f);
1558: PetscMalloc(info.xm*info.ym*info.zm*info.dof*sizeof(PetscScalar),&f);
1560: /* Seed the input array g_u with coloring information */
1561:
1562: VecGetArray(vu,&u);
1563: (lf)(&Nc,&info,u,g_u,&Nc,f,g_f,&Nc,w,&ierr);
1564: VecRestoreArray(vu,&u);
1566: /* stick the values into the matrix */
1567: /* PetscScalarView(Nc*info.xm*info.ym,g_f,0); */
1568: MatSetValuesAdifor(J,Nc,g_f);
1570: /* return space for derivative objects. */
1571: PetscFree(g_u);
1572: PetscFree(g_f);
1573: PetscFree(f);
1574: return(0);
1575: }
1577: #undef __FUNCT__
1579: /*@C
1580: DAMultiplyByJacobian1WithAD - Applies a Jacobian function supplied by ADIFOR or ADIC
1581: to a vector on each processor that shares a DA.
1583: Input Parameters:
1584: + da - the DA that defines the grid
1585: . vu - Jacobian is computed at this point (ghosted)
1586: . v - product is done on this vector (ghosted)
1587: . fu - output vector = J(vu)*v (not ghosted)
1588: - w - any user data
1590: Notes:
1591: This routine does NOT do ghost updates on vu and v upon entry.
1592:
1593: Automatically calls DAMultiplyByJacobian1WithAdifor() or DAMultiplyByJacobian1WithAdic()
1594: depending on whether DASetLocalAdicMFFunction() or DASetLocalAdiforMFFunction() was called.
1596: Level: advanced
1598: .seealso: DAFormFunction1(), DAMultiplyByJacobian1WithAdifor(), DAMultiplyByJacobian1WithAdic()
1600: @*/
1601: int DAMultiplyByJacobian1WithAD(DA da,Vec u,Vec v,Vec f,void *w)
1602: {
1603: int ierr;
1606: if (da->adicmf_lf) {
1607: #if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE)
1608: DAMultiplyByJacobian1WithAdic(da,u,v,f,w);
1609: #else
1610: SETERRQ(1,"Requires ADIC to be installed and cannot use complex numbers");
1611: #endif
1612: } else if (da->adiformf_lf) {
1613: DAMultiplyByJacobian1WithAdifor(da,u,v,f,w);
1614: } else {
1615: SETERRQ(1,"Must call DASetLocalAdiforMFFunction() or DASetLocalAdicMFFunction() before using");
1616: }
1617: return(0);
1618: }
1621: #undef __FUNCT__
1623: /*@C
1624: DAMultiplyByJacobian1WithAdifor - Applies a ADIFOR provided Jacobian function on each processor that
1625: share a DA to a vector
1627: Input Parameters:
1628: + da - the DA that defines the grid
1629: . vu - Jacobian is computed at this point (ghosted)
1630: . v - product is done on this vector (ghosted)
1631: . fu - output vector = J(vu)*v (not ghosted)
1632: - w - any user data
1634: Notes: Does NOT do ghost updates on vu and v upon entry
1636: Level: advanced
1638: .seealso: DAFormFunction1()
1640: @*/
1641: int DAMultiplyByJacobian1WithAdifor(DA da,Vec u,Vec v,Vec f,void *w)
1642: {
1643: int ierr;
1644: PetscScalar *au,*av,*af,*awork;
1645: Vec work;
1646: DALocalInfo info;
1647: void (*lf)(DALocalInfo*,PetscScalar*,PetscScalar*,PetscScalar*,PetscScalar*,void*,int*) =
1648: (void (*)(DALocalInfo*,PetscScalar*,PetscScalar*,PetscScalar*,PetscScalar*,void*,int*))*da->adiformf_lf;
1651: DAGetLocalInfo(da,&info);
1653: DAGetGlobalVector(da,&work);
1654: VecGetArray(u,&au);
1655: VecGetArray(v,&av);
1656: VecGetArray(f,&af);
1657: VecGetArray(work,&awork);
1658: (lf)(&info,au,av,awork,af,w,&ierr);
1659: VecRestoreArray(u,&au);
1660: VecRestoreArray(v,&av);
1661: VecRestoreArray(f,&af);
1662: VecRestoreArray(work,&awork);
1663: DARestoreGlobalVector(da,&work);
1665: return(0);
1666: }
1668: #undef __FUNCT__
1670: /*@C
1671: DASetInterpolationType - Sets the type of interpolation that will be
1672: returned by DAGetInterpolation()
1674: Collective on DA
1676: Input Parameter:
1677: + da - initial distributed array
1678: . ctype - DA_Q1 and DA_Q0 are currently the only supported forms
1680: Level: intermediate
1682: .keywords: distributed array, interpolation
1684: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy(), DA, DAInterpolationType
1685: @*/
1686: int DASetInterpolationType(DA da,DAInterpolationType ctype)
1687: {
1690: da->interptype = ctype;
1691: return(0);
1692: }