Actual source code: minsurf2.c
1: /* Program usage: mpiexec -n <proc> minsurf2 [-help] [all TAO options] */
3: /*
4: Include "petsctao.h" so we can use TAO solvers.
5: petscdmda.h for distributed array
6: */
7: #include <petsctao.h>
8: #include <petscdmda.h>
10: static char help[] =
11: "This example demonstrates use of the TAO package to \n\
12: solve an unconstrained minimization problem. This example is based on a \n\
13: problem from the MINPACK-2 test suite. Given a rectangular 2-D domain and \n\
14: boundary values along the edges of the domain, the objective is to find the\n\
15: surface with the minimal area that satisfies the boundary conditions.\n\
16: The command line options are:\n\
17: -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
18: -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
19: -start <st>, where <st> =0 for zero vector, <st> >0 for random start, and <st> <0 \n\
20: for an average of the boundary conditions\n\n";
22: /*T
23: Concepts: TAO^Solving an unconstrained minimization problem
24: Routines: TaoCreate(); TaoSetType();
25: Routines: TaoSetSolution();
26: Routines: TaoSetObjectiveAndGradient();
27: Routines: TaoSetHessian(); TaoSetFromOptions();
28: Routines: TaoSetMonitor();
29: Routines: TaoSolve(); TaoView();
30: Routines: TaoDestroy();
31: Processors: n
32: T*/
34: /*
35: User-defined application context - contains data needed by the
36: application-provided call-back routines, FormFunctionGradient()
37: and FormHessian().
38: */
39: typedef struct {
40: PetscInt mx, my; /* discretization in x, y directions */
41: PetscReal *bottom, *top, *left, *right; /* boundary values */
42: DM dm; /* distributed array data structure */
43: Mat H; /* Hessian */
44: } AppCtx;
46: /* -------- User-defined Routines --------- */
48: static PetscErrorCode MSA_BoundaryConditions(AppCtx*);
49: static PetscErrorCode MSA_InitialPoint(AppCtx*,Vec);
50: PetscErrorCode QuadraticH(AppCtx*,Vec,Mat);
51: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void*);
52: PetscErrorCode FormGradient(Tao,Vec,Vec,void*);
53: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat,void*);
54: PetscErrorCode My_Monitor(Tao, void *);
56: int main(int argc, char **argv)
57: {
58: PetscInt Nx, Ny; /* number of processors in x- and y- directions */
59: Vec x; /* solution, gradient vectors */
60: PetscBool flg, viewmat; /* flags */
61: PetscBool fddefault, fdcoloring; /* flags */
62: Tao tao; /* TAO solver context */
63: AppCtx user; /* user-defined work context */
64: ISColoring iscoloring;
65: MatFDColoring matfdcoloring;
67: /* Initialize TAO */
68: PetscInitialize(&argc, &argv,(char *)0,help);
70: /* Specify dimension of the problem */
71: user.mx = 10; user.my = 10;
73: /* Check for any command line arguments that override defaults */
74: PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,&flg);
75: PetscOptionsGetInt(NULL,NULL,"-my",&user.my,&flg);
77: PetscPrintf(MPI_COMM_WORLD,"\n---- Minimum Surface Area Problem -----\n");
78: PetscPrintf(MPI_COMM_WORLD,"mx: %D my: %D \n\n",user.mx,user.my);
80: /* Let PETSc determine the vector distribution */
81: Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
83: /* Create distributed array (DM) to manage parallel grid and vectors */
84: DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,user.mx, user.my,Nx,Ny,1,1,NULL,NULL,&user.dm);
85: DMSetFromOptions(user.dm);
86: DMSetUp(user.dm);
88: /* Create TAO solver and set desired solution method.*/
89: TaoCreate(PETSC_COMM_WORLD,&tao);
90: TaoSetType(tao,TAOCG);
92: /*
93: Extract global vector from DA for the vector of variables -- PETSC routine
94: Compute the initial solution -- application specific, see below
95: Set this vector for use by TAO -- TAO routine
96: */
97: DMCreateGlobalVector(user.dm,&x);
98: MSA_BoundaryConditions(&user);
99: MSA_InitialPoint(&user,x);
100: TaoSetSolution(tao,x);
102: /*
103: Initialize the Application context for use in function evaluations -- application specific, see below.
104: Set routines for function and gradient evaluation
105: */
106: TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void *)&user);
108: /*
109: Given the command line arguments, calculate the hessian with either the user-
110: provided function FormHessian, or the default finite-difference driven Hessian
111: functions
112: */
113: PetscOptionsHasName(NULL,NULL,"-tao_fddefault",&fddefault);
114: PetscOptionsHasName(NULL,NULL,"-tao_fdcoloring",&fdcoloring);
116: /*
117: Create a matrix data structure to store the Hessian and set
118: the Hessian evalution routine.
119: Set the matrix structure to be used for Hessian evalutions
120: */
121: DMCreateMatrix(user.dm,&user.H);
122: MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE);
124: if (fdcoloring) {
125: DMCreateColoring(user.dm,IS_COLORING_GLOBAL,&iscoloring);
126: MatFDColoringCreate(user.H,iscoloring,&matfdcoloring);
127: ISColoringDestroy(&iscoloring);
128: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormGradient,(void*)&user);
129: MatFDColoringSetFromOptions(matfdcoloring);
130: TaoSetHessian(tao,user.H,user.H,TaoDefaultComputeHessianColor,(void *)matfdcoloring);
131: } else if (fddefault) {
132: TaoSetHessian(tao,user.H,user.H,TaoDefaultComputeHessian,(void *)NULL);
133: } else {
134: TaoSetHessian(tao,user.H,user.H,FormHessian,(void *)&user);
135: }
137: /*
138: If my_monitor option is in command line, then use the user-provided
139: monitoring function
140: */
141: PetscOptionsHasName(NULL,NULL,"-my_monitor",&viewmat);
142: if (viewmat) {
143: TaoSetMonitor(tao,My_Monitor,NULL,NULL);
144: }
146: /* Check for any tao command line options */
147: TaoSetFromOptions(tao);
149: /* SOLVE THE APPLICATION */
150: TaoSolve(tao);
152: TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);
154: /* Free TAO data structures */
155: TaoDestroy(&tao);
157: /* Free PETSc data structures */
158: VecDestroy(&x);
159: MatDestroy(&user.H);
160: if (fdcoloring) {
161: MatFDColoringDestroy(&matfdcoloring);
162: }
163: PetscFree(user.bottom);
164: PetscFree(user.top);
165: PetscFree(user.left);
166: PetscFree(user.right);
167: DMDestroy(&user.dm);
168: PetscFinalize();
169: return 0;
170: }
172: PetscErrorCode FormGradient(Tao tao, Vec X, Vec G,void *userCtx)
173: {
174: PetscReal fcn;
176: FormFunctionGradient(tao,X,&fcn,G,userCtx);
177: return 0;
178: }
180: /* -------------------------------------------------------------------- */
181: /* FormFunctionGradient - Evaluates the function and corresponding gradient.
183: Input Parameters:
184: . tao - the Tao context
185: . XX - input vector
186: . userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradient()
188: Output Parameters:
189: . fcn - the newly evaluated function
190: . GG - vector containing the newly evaluated gradient
191: */
192: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn,Vec G,void *userCtx)
193: {
194: AppCtx *user = (AppCtx *) userCtx;
195: PetscInt i,j;
196: PetscInt mx=user->mx, my=user->my;
197: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
198: PetscReal ft=0.0;
199: PetscReal hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy;
200: PetscReal rhx=mx+1, rhy=my+1;
201: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
202: PetscReal df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
203: PetscReal **g, **x;
204: Vec localX;
206: /* Get local mesh boundaries */
207: DMGetLocalVector(user->dm,&localX);
208: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
209: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
211: /* Scatter ghost points to local vector */
212: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
213: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
215: /* Get pointers to vector data */
216: DMDAVecGetArray(user->dm,localX,(void**)&x);
217: DMDAVecGetArray(user->dm,G,(void**)&g);
219: /* Compute function and gradient over the locally owned part of the mesh */
220: for (j=ys; j<ys+ym; j++) {
221: for (i=xs; i< xs+xm; i++) {
223: xc = x[j][i];
224: xlt=xrb=xl=xr=xb=xt=xc;
226: if (i==0) { /* left side */
227: xl= user->left[j-ys+1];
228: xlt = user->left[j-ys+2];
229: } else {
230: xl = x[j][i-1];
231: }
233: if (j==0) { /* bottom side */
234: xb=user->bottom[i-xs+1];
235: xrb =user->bottom[i-xs+2];
236: } else {
237: xb = x[j-1][i];
238: }
240: if (i+1 == gxs+gxm) { /* right side */
241: xr=user->right[j-ys+1];
242: xrb = user->right[j-ys];
243: } else {
244: xr = x[j][i+1];
245: }
247: if (j+1==gys+gym) { /* top side */
248: xt=user->top[i-xs+1];
249: xlt = user->top[i-xs];
250: }else {
251: xt = x[j+1][i];
252: }
254: if (i>gxs && j+1<gys+gym) {
255: xlt = x[j+1][i-1];
256: }
257: if (j>gys && i+1<gxs+gxm) {
258: xrb = x[j-1][i+1];
259: }
261: d1 = (xc-xl);
262: d2 = (xc-xr);
263: d3 = (xc-xt);
264: d4 = (xc-xb);
265: d5 = (xr-xrb);
266: d6 = (xrb-xb);
267: d7 = (xlt-xl);
268: d8 = (xt-xlt);
270: df1dxc = d1*hydhx;
271: df2dxc = (d1*hydhx + d4*hxdhy);
272: df3dxc = d3*hxdhy;
273: df4dxc = (d2*hydhx + d3*hxdhy);
274: df5dxc = d2*hydhx;
275: df6dxc = d4*hxdhy;
277: d1 *= rhx;
278: d2 *= rhx;
279: d3 *= rhy;
280: d4 *= rhy;
281: d5 *= rhy;
282: d6 *= rhx;
283: d7 *= rhy;
284: d8 *= rhx;
286: f1 = PetscSqrtReal(1.0 + d1*d1 + d7*d7);
287: f2 = PetscSqrtReal(1.0 + d1*d1 + d4*d4);
288: f3 = PetscSqrtReal(1.0 + d3*d3 + d8*d8);
289: f4 = PetscSqrtReal(1.0 + d3*d3 + d2*d2);
290: f5 = PetscSqrtReal(1.0 + d2*d2 + d5*d5);
291: f6 = PetscSqrtReal(1.0 + d4*d4 + d6*d6);
293: ft = ft + (f2 + f4);
295: df1dxc /= f1;
296: df2dxc /= f2;
297: df3dxc /= f3;
298: df4dxc /= f4;
299: df5dxc /= f5;
300: df6dxc /= f6;
302: g[j][i] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc) * 0.5;
304: }
305: }
307: /* Compute triangular areas along the border of the domain. */
308: if (xs==0) { /* left side */
309: for (j=ys; j<ys+ym; j++) {
310: d3=(user->left[j-ys+1] - user->left[j-ys+2])*rhy;
311: d2=(user->left[j-ys+1] - x[j][0]) *rhx;
312: ft = ft+PetscSqrtReal(1.0 + d3*d3 + d2*d2);
313: }
314: }
315: if (ys==0) { /* bottom side */
316: for (i=xs; i<xs+xm; i++) {
317: d2=(user->bottom[i+1-xs]-user->bottom[i-xs+2])*rhx;
318: d3=(user->bottom[i-xs+1]-x[0][i])*rhy;
319: ft = ft+PetscSqrtReal(1.0 + d3*d3 + d2*d2);
320: }
321: }
323: if (xs+xm==mx) { /* right side */
324: for (j=ys; j< ys+ym; j++) {
325: d1=(x[j][mx-1] - user->right[j-ys+1])*rhx;
326: d4=(user->right[j-ys]-user->right[j-ys+1])*rhy;
327: ft = ft+PetscSqrtReal(1.0 + d1*d1 + d4*d4);
328: }
329: }
330: if (ys+ym==my) { /* top side */
331: for (i=xs; i<xs+xm; i++) {
332: d1=(x[my-1][i] - user->top[i-xs+1])*rhy;
333: d4=(user->top[i-xs+1] - user->top[i-xs])*rhx;
334: ft = ft+PetscSqrtReal(1.0 + d1*d1 + d4*d4);
335: }
336: }
338: if (ys==0 && xs==0) {
339: d1=(user->left[0]-user->left[1])/hy;
340: d2=(user->bottom[0]-user->bottom[1])*rhx;
341: ft +=PetscSqrtReal(1.0 + d1*d1 + d2*d2);
342: }
343: if (ys+ym == my && xs+xm == mx) {
344: d1=(user->right[ym+1] - user->right[ym])*rhy;
345: d2=(user->top[xm+1] - user->top[xm])*rhx;
346: ft +=PetscSqrtReal(1.0 + d1*d1 + d2*d2);
347: }
349: ft=ft*area;
350: MPI_Allreduce(&ft,fcn,1,MPIU_REAL,MPIU_SUM,MPI_COMM_WORLD);
352: /* Restore vectors */
353: DMDAVecRestoreArray(user->dm,localX,(void**)&x);
354: DMDAVecRestoreArray(user->dm,G,(void**)&g);
356: /* Scatter values to global vector */
357: DMRestoreLocalVector(user->dm,&localX);
358: PetscLogFlops(67.0*xm*ym);
359: return 0;
360: }
362: /* ------------------------------------------------------------------- */
363: /*
364: FormHessian - Evaluates Hessian matrix.
366: Input Parameters:
367: . tao - the Tao context
368: . x - input vector
369: . ptr - optional user-defined context, as set by TaoSetHessian()
371: Output Parameters:
372: . H - Hessian matrix
373: . Hpre - optionally different preconditioning matrix
374: . flg - flag indicating matrix structure
376: */
377: PetscErrorCode FormHessian(Tao tao,Vec X,Mat H, Mat Hpre, void *ptr)
378: {
379: AppCtx *user = (AppCtx *) ptr;
381: /* Evaluate the Hessian entries*/
382: QuadraticH(user,X,H);
383: return 0;
384: }
386: /* ------------------------------------------------------------------- */
387: /*
388: QuadraticH - Evaluates Hessian matrix.
390: Input Parameters:
391: . user - user-defined context, as set by TaoSetHessian()
392: . X - input vector
394: Output Parameter:
395: . H - Hessian matrix
396: */
397: PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian)
398: {
399: PetscInt i,j,k;
400: PetscInt mx=user->mx, my=user->my;
401: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
402: PetscReal hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
403: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
404: PetscReal hl,hr,ht,hb,hc,htl,hbr;
405: PetscReal **x, v[7];
406: MatStencil col[7],row;
407: Vec localX;
408: PetscBool assembled;
410: /* Get local mesh boundaries */
411: DMGetLocalVector(user->dm,&localX);
413: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
414: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
416: /* Scatter ghost points to local vector */
417: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
418: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
420: /* Get pointers to vector data */
421: DMDAVecGetArray(user->dm,localX,(void**)&x);
423: /* Initialize matrix entries to zero */
424: MatAssembled(Hessian,&assembled);
425: if (assembled) MatZeroEntries(Hessian);
427: /* Set various matrix options */
428: MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
430: /* Compute Hessian over the locally owned part of the mesh */
432: for (j=ys; j<ys+ym; j++) {
434: for (i=xs; i< xs+xm; i++) {
436: xc = x[j][i];
437: xlt=xrb=xl=xr=xb=xt=xc;
439: /* Left side */
440: if (i==0) {
441: xl = user->left[j-ys+1];
442: xlt = user->left[j-ys+2];
443: } else {
444: xl = x[j][i-1];
445: }
447: if (j==0) {
448: xb = user->bottom[i-xs+1];
449: xrb = user->bottom[i-xs+2];
450: } else {
451: xb = x[j-1][i];
452: }
454: if (i+1 == mx) {
455: xr = user->right[j-ys+1];
456: xrb = user->right[j-ys];
457: } else {
458: xr = x[j][i+1];
459: }
461: if (j+1==my) {
462: xt = user->top[i-xs+1];
463: xlt = user->top[i-xs];
464: }else {
465: xt = x[j+1][i];
466: }
468: if (i>0 && j+1<my) {
469: xlt = x[j+1][i-1];
470: }
471: if (j>0 && i+1<mx) {
472: xrb = x[j-1][i+1];
473: }
475: d1 = (xc-xl)/hx;
476: d2 = (xc-xr)/hx;
477: d3 = (xc-xt)/hy;
478: d4 = (xc-xb)/hy;
479: d5 = (xrb-xr)/hy;
480: d6 = (xrb-xb)/hx;
481: d7 = (xlt-xl)/hy;
482: d8 = (xlt-xt)/hx;
484: f1 = PetscSqrtReal(1.0 + d1*d1 + d7*d7);
485: f2 = PetscSqrtReal(1.0 + d1*d1 + d4*d4);
486: f3 = PetscSqrtReal(1.0 + d3*d3 + d8*d8);
487: f4 = PetscSqrtReal(1.0 + d3*d3 + d2*d2);
488: f5 = PetscSqrtReal(1.0 + d2*d2 + d5*d5);
489: f6 = PetscSqrtReal(1.0 + d4*d4 + d6*d6);
491: hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
492: (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
493: hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
494: (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
495: ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
496: (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
497: hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
498: (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);
500: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6);
501: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3);
503: hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
504: hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
505: (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
506: (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);
508: hl/=2.0; hr/=2.0; ht/=2.0; hb/=2.0; hbr/=2.0; htl/=2.0; hc/=2.0;
510: row.j = j; row.i = i;
511: k=0;
512: if (j>0) {
513: v[k]=hb;
514: col[k].j = j - 1; col[k].i = i;
515: k++;
516: }
518: if (j>0 && i < mx -1) {
519: v[k]=hbr;
520: col[k].j = j - 1; col[k].i = i+1;
521: k++;
522: }
524: if (i>0) {
525: v[k]= hl;
526: col[k].j = j; col[k].i = i-1;
527: k++;
528: }
530: v[k]= hc;
531: col[k].j = j; col[k].i = i;
532: k++;
534: if (i < mx-1) {
535: v[k]= hr;
536: col[k].j = j; col[k].i = i+1;
537: k++;
538: }
540: if (i>0 && j < my-1) {
541: v[k]= htl;
542: col[k].j = j+1; col[k].i = i-1;
543: k++;
544: }
546: if (j < my-1) {
547: v[k]= ht;
548: col[k].j = j+1; col[k].i = i;
549: k++;
550: }
552: /*
553: Set matrix values using local numbering, which was defined
554: earlier, in the main routine.
555: */
556: MatSetValuesStencil(Hessian,1,&row,k,col,v,INSERT_VALUES);
557: }
558: }
560: DMDAVecRestoreArray(user->dm,localX,(void**)&x);
561: DMRestoreLocalVector(user->dm,&localX);
563: MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY);
564: MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY);
566: PetscLogFlops(199.0*xm*ym);
567: return 0;
568: }
570: /* ------------------------------------------------------------------- */
571: /*
572: MSA_BoundaryConditions - Calculates the boundary conditions for
573: the region.
575: Input Parameter:
576: . user - user-defined application context
578: Output Parameter:
579: . user - user-defined application context
580: */
581: static PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
582: {
583: PetscInt i,j,k,limit=0,maxits=5;
584: PetscInt xs,ys,xm,ym,gxs,gys,gxm,gym;
585: PetscInt mx=user->mx,my=user->my;
586: PetscInt bsize=0, lsize=0, tsize=0, rsize=0;
587: PetscReal one=1.0, two=2.0, three=3.0, tol=1e-10;
588: PetscReal fnorm,det,hx,hy,xt=0,yt=0;
589: PetscReal u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
590: PetscReal b=-0.5, t=0.5, l=-0.5, r=0.5;
591: PetscReal *boundary;
592: PetscBool flg;
594: /* Get local mesh boundaries */
595: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
596: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
598: bsize=xm+2;
599: lsize=ym+2;
600: rsize=ym+2;
601: tsize=xm+2;
603: PetscMalloc1(bsize,&user->bottom);
604: PetscMalloc1(tsize,&user->top);
605: PetscMalloc1(lsize,&user->left);
606: PetscMalloc1(rsize,&user->right);
608: hx= (r-l)/(mx+1); hy=(t-b)/(my+1);
610: for (j=0; j<4; j++) {
611: if (j==0) {
612: yt=b;
613: xt=l+hx*xs;
614: limit=bsize;
615: boundary=user->bottom;
616: } else if (j==1) {
617: yt=t;
618: xt=l+hx*xs;
619: limit=tsize;
620: boundary=user->top;
621: } else if (j==2) {
622: yt=b+hy*ys;
623: xt=l;
624: limit=lsize;
625: boundary=user->left;
626: } else { /* if (j==3) */
627: yt=b+hy*ys;
628: xt=r;
629: limit=rsize;
630: boundary=user->right;
631: }
633: for (i=0; i<limit; i++) {
634: u1=xt;
635: u2=-yt;
636: for (k=0; k<maxits; k++) {
637: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
638: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
639: fnorm=PetscSqrtReal(nf1*nf1+nf2*nf2);
640: if (fnorm <= tol) break;
641: njac11=one+u2*u2-u1*u1;
642: njac12=two*u1*u2;
643: njac21=-two*u1*u2;
644: njac22=-one - u1*u1 + u2*u2;
645: det = njac11*njac22-njac21*njac12;
646: u1 = u1-(njac22*nf1-njac12*nf2)/det;
647: u2 = u2-(njac11*nf2-njac21*nf1)/det;
648: }
650: boundary[i]=u1*u1-u2*u2;
651: if (j==0 || j==1) {
652: xt=xt+hx;
653: } else { /* if (j==2 || j==3) */
654: yt=yt+hy;
655: }
657: }
659: }
661: /* Scale the boundary if desired */
662: if (1==1) {
663: PetscReal scl = 1.0;
665: PetscOptionsGetReal(NULL,NULL,"-bottom",&scl,&flg);
666: if (flg) {
667: for (i=0;i<bsize;i++) user->bottom[i]*=scl;
668: }
670: PetscOptionsGetReal(NULL,NULL,"-top",&scl,&flg);
671: if (flg) {
672: for (i=0;i<tsize;i++) user->top[i]*=scl;
673: }
675: PetscOptionsGetReal(NULL,NULL,"-right",&scl,&flg);
676: if (flg) {
677: for (i=0;i<rsize;i++) user->right[i]*=scl;
678: }
680: PetscOptionsGetReal(NULL,NULL,"-left",&scl,&flg);
681: if (flg) {
682: for (i=0;i<lsize;i++) user->left[i]*=scl;
683: }
684: }
685: return 0;
686: }
688: /* ------------------------------------------------------------------- */
689: /*
690: MSA_InitialPoint - Calculates the initial guess in one of three ways.
692: Input Parameters:
693: . user - user-defined application context
694: . X - vector for initial guess
696: Output Parameters:
697: . X - newly computed initial guess
698: */
699: static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
700: {
701: PetscInt start2=-1,i,j;
702: PetscReal start1=0;
703: PetscBool flg1,flg2;
705: PetscOptionsGetReal(NULL,NULL,"-start",&start1,&flg1);
706: PetscOptionsGetInt(NULL,NULL,"-random",&start2,&flg2);
708: if (flg1) { /* The zero vector is reasonable */
710: VecSet(X, start1);
712: } else if (flg2 && start2>0) { /* Try a random start between -0.5 and 0.5 */
713: PetscRandom rctx; PetscReal np5=-0.5;
715: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
716: VecSetRandom(X, rctx);
717: PetscRandomDestroy(&rctx);
718: VecShift(X, np5);
720: } else { /* Take an average of the boundary conditions */
721: PetscInt xs,xm,ys,ym;
722: PetscInt mx=user->mx,my=user->my;
723: PetscReal **x;
725: /* Get local mesh boundaries */
726: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
728: /* Get pointers to vector data */
729: DMDAVecGetArray(user->dm,X,(void**)&x);
731: /* Perform local computations */
732: for (j=ys; j<ys+ym; j++) {
733: for (i=xs; i< xs+xm; i++) {
734: x[j][i] = (((j+1)*user->bottom[i-xs+1]+(my-j+1)*user->top[i-xs+1])/(my+2)+((i+1)*user->left[j-ys+1]+(mx-i+1)*user->right[j-ys+1])/(mx+2))/2.0;
735: }
736: }
737: DMDAVecRestoreArray(user->dm,X,(void**)&x);
738: PetscLogFlops(9.0*xm*ym);
739: }
740: return 0;
741: }
743: /*-----------------------------------------------------------------------*/
744: PetscErrorCode My_Monitor(Tao tao, void *ctx)
745: {
746: Vec X;
748: TaoGetSolution(tao,&X);
749: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
750: return 0;
751: }
753: /*TEST
755: build:
756: requires: !complex
758: test:
759: args: -tao_smonitor -tao_type lmvm -mx 10 -my 8 -tao_gatol 1.e-3
760: requires: !single
762: test:
763: suffix: 2
764: nsize: 2
765: args: -tao_smonitor -tao_type nls -tao_nls_ksp_max_it 15 -tao_gatol 1.e-4
766: filter: grep -v "nls ksp"
767: requires: !single
769: test:
770: suffix: 3
771: nsize: 3
772: args: -tao_smonitor -tao_type cg -tao_cg_type fr -mx 10 -my 10 -tao_gatol 1.e-3
773: requires: !single
775: test:
776: suffix: 5
777: nsize: 2
778: args: -tao_smonitor -tao_type bmrm -mx 10 -my 8 -tao_gatol 1.e-3
779: requires: !single
781: TEST*/