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: /*
23: User-defined application context - contains data needed by the
24: application-provided call-back routines, FormFunctionGradient()
25: and FormHessian().
26: */
27: typedef struct {
28: PetscInt mx, my; /* discretization in x, y directions */
29: PetscReal *bottom, *top, *left, *right; /* boundary values */
30: DM dm; /* distributed array data structure */
31: Mat H; /* Hessian */
32: } AppCtx;
34: /* -------- User-defined Routines --------- */
36: static PetscErrorCode MSA_BoundaryConditions(AppCtx*);
37: static PetscErrorCode MSA_InitialPoint(AppCtx*,Vec);
38: PetscErrorCode QuadraticH(AppCtx*,Vec,Mat);
39: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void*);
40: PetscErrorCode FormGradient(Tao,Vec,Vec,void*);
41: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat,void*);
42: PetscErrorCode My_Monitor(Tao, void *);
44: int main(int argc, char **argv)
45: {
46: PetscInt Nx, Ny; /* number of processors in x- and y- directions */
47: Vec x; /* solution, gradient vectors */
48: PetscBool flg, viewmat; /* flags */
49: PetscBool fddefault, fdcoloring; /* flags */
50: Tao tao; /* TAO solver context */
51: AppCtx user; /* user-defined work context */
52: ISColoring iscoloring;
53: MatFDColoring matfdcoloring;
55: /* Initialize TAO */
56: PetscInitialize(&argc, &argv,(char *)0,help);
58: /* Specify dimension of the problem */
59: user.mx = 10; user.my = 10;
61: /* Check for any command line arguments that override defaults */
62: PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,&flg);
63: PetscOptionsGetInt(NULL,NULL,"-my",&user.my,&flg);
65: PetscPrintf(MPI_COMM_WORLD,"\n---- Minimum Surface Area Problem -----\n");
66: PetscPrintf(MPI_COMM_WORLD,"mx: %D my: %D \n\n",user.mx,user.my);
68: /* Let PETSc determine the vector distribution */
69: Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
71: /* Create distributed array (DM) to manage parallel grid and vectors */
72: 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);
73: DMSetFromOptions(user.dm);
74: DMSetUp(user.dm);
76: /* Create TAO solver and set desired solution method.*/
77: TaoCreate(PETSC_COMM_WORLD,&tao);
78: TaoSetType(tao,TAOCG);
80: /*
81: Extract global vector from DA for the vector of variables -- PETSC routine
82: Compute the initial solution -- application specific, see below
83: Set this vector for use by TAO -- TAO routine
84: */
85: DMCreateGlobalVector(user.dm,&x);
86: MSA_BoundaryConditions(&user);
87: MSA_InitialPoint(&user,x);
88: TaoSetSolution(tao,x);
90: /*
91: Initialize the Application context for use in function evaluations -- application specific, see below.
92: Set routines for function and gradient evaluation
93: */
94: TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void *)&user);
96: /*
97: Given the command line arguments, calculate the hessian with either the user-
98: provided function FormHessian, or the default finite-difference driven Hessian
99: functions
100: */
101: PetscOptionsHasName(NULL,NULL,"-tao_fddefault",&fddefault);
102: PetscOptionsHasName(NULL,NULL,"-tao_fdcoloring",&fdcoloring);
104: /*
105: Create a matrix data structure to store the Hessian and set
106: the Hessian evalution routine.
107: Set the matrix structure to be used for Hessian evalutions
108: */
109: DMCreateMatrix(user.dm,&user.H);
110: MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE);
112: if (fdcoloring) {
113: DMCreateColoring(user.dm,IS_COLORING_GLOBAL,&iscoloring);
114: MatFDColoringCreate(user.H,iscoloring,&matfdcoloring);
115: ISColoringDestroy(&iscoloring);
116: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormGradient,(void*)&user);
117: MatFDColoringSetFromOptions(matfdcoloring);
118: TaoSetHessian(tao,user.H,user.H,TaoDefaultComputeHessianColor,(void *)matfdcoloring);
119: } else if (fddefault) {
120: TaoSetHessian(tao,user.H,user.H,TaoDefaultComputeHessian,(void *)NULL);
121: } else {
122: TaoSetHessian(tao,user.H,user.H,FormHessian,(void *)&user);
123: }
125: /*
126: If my_monitor option is in command line, then use the user-provided
127: monitoring function
128: */
129: PetscOptionsHasName(NULL,NULL,"-my_monitor",&viewmat);
130: if (viewmat) {
131: TaoSetMonitor(tao,My_Monitor,NULL,NULL);
132: }
134: /* Check for any tao command line options */
135: TaoSetFromOptions(tao);
137: /* SOLVE THE APPLICATION */
138: TaoSolve(tao);
140: TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);
142: /* Free TAO data structures */
143: TaoDestroy(&tao);
145: /* Free PETSc data structures */
146: VecDestroy(&x);
147: MatDestroy(&user.H);
148: if (fdcoloring) {
149: MatFDColoringDestroy(&matfdcoloring);
150: }
151: PetscFree(user.bottom);
152: PetscFree(user.top);
153: PetscFree(user.left);
154: PetscFree(user.right);
155: DMDestroy(&user.dm);
156: PetscFinalize();
157: return 0;
158: }
160: PetscErrorCode FormGradient(Tao tao, Vec X, Vec G,void *userCtx)
161: {
162: PetscReal fcn;
164: FormFunctionGradient(tao,X,&fcn,G,userCtx);
165: return 0;
166: }
168: /* -------------------------------------------------------------------- */
169: /* FormFunctionGradient - Evaluates the function and corresponding gradient.
171: Input Parameters:
172: . tao - the Tao context
173: . XX - input vector
174: . userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradient()
176: Output Parameters:
177: . fcn - the newly evaluated function
178: . GG - vector containing the newly evaluated gradient
179: */
180: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn,Vec G,void *userCtx)
181: {
182: AppCtx *user = (AppCtx *) userCtx;
183: PetscInt i,j;
184: PetscInt mx=user->mx, my=user->my;
185: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
186: PetscReal ft=0.0;
187: PetscReal hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy;
188: PetscReal rhx=mx+1, rhy=my+1;
189: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
190: PetscReal df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
191: PetscReal **g, **x;
192: Vec localX;
194: /* Get local mesh boundaries */
195: DMGetLocalVector(user->dm,&localX);
196: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
197: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
199: /* Scatter ghost points to local vector */
200: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
201: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
203: /* Get pointers to vector data */
204: DMDAVecGetArray(user->dm,localX,(void**)&x);
205: DMDAVecGetArray(user->dm,G,(void**)&g);
207: /* Compute function and gradient over the locally owned part of the mesh */
208: for (j=ys; j<ys+ym; j++) {
209: for (i=xs; i< xs+xm; i++) {
211: xc = x[j][i];
212: xlt=xrb=xl=xr=xb=xt=xc;
214: if (i==0) { /* left side */
215: xl= user->left[j-ys+1];
216: xlt = user->left[j-ys+2];
217: } else {
218: xl = x[j][i-1];
219: }
221: if (j==0) { /* bottom side */
222: xb=user->bottom[i-xs+1];
223: xrb =user->bottom[i-xs+2];
224: } else {
225: xb = x[j-1][i];
226: }
228: if (i+1 == gxs+gxm) { /* right side */
229: xr=user->right[j-ys+1];
230: xrb = user->right[j-ys];
231: } else {
232: xr = x[j][i+1];
233: }
235: if (j+1==gys+gym) { /* top side */
236: xt=user->top[i-xs+1];
237: xlt = user->top[i-xs];
238: }else {
239: xt = x[j+1][i];
240: }
242: if (i>gxs && j+1<gys+gym) {
243: xlt = x[j+1][i-1];
244: }
245: if (j>gys && i+1<gxs+gxm) {
246: xrb = x[j-1][i+1];
247: }
249: d1 = (xc-xl);
250: d2 = (xc-xr);
251: d3 = (xc-xt);
252: d4 = (xc-xb);
253: d5 = (xr-xrb);
254: d6 = (xrb-xb);
255: d7 = (xlt-xl);
256: d8 = (xt-xlt);
258: df1dxc = d1*hydhx;
259: df2dxc = (d1*hydhx + d4*hxdhy);
260: df3dxc = d3*hxdhy;
261: df4dxc = (d2*hydhx + d3*hxdhy);
262: df5dxc = d2*hydhx;
263: df6dxc = d4*hxdhy;
265: d1 *= rhx;
266: d2 *= rhx;
267: d3 *= rhy;
268: d4 *= rhy;
269: d5 *= rhy;
270: d6 *= rhx;
271: d7 *= rhy;
272: d8 *= rhx;
274: f1 = PetscSqrtReal(1.0 + d1*d1 + d7*d7);
275: f2 = PetscSqrtReal(1.0 + d1*d1 + d4*d4);
276: f3 = PetscSqrtReal(1.0 + d3*d3 + d8*d8);
277: f4 = PetscSqrtReal(1.0 + d3*d3 + d2*d2);
278: f5 = PetscSqrtReal(1.0 + d2*d2 + d5*d5);
279: f6 = PetscSqrtReal(1.0 + d4*d4 + d6*d6);
281: ft = ft + (f2 + f4);
283: df1dxc /= f1;
284: df2dxc /= f2;
285: df3dxc /= f3;
286: df4dxc /= f4;
287: df5dxc /= f5;
288: df6dxc /= f6;
290: g[j][i] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc) * 0.5;
292: }
293: }
295: /* Compute triangular areas along the border of the domain. */
296: if (xs==0) { /* left side */
297: for (j=ys; j<ys+ym; j++) {
298: d3=(user->left[j-ys+1] - user->left[j-ys+2])*rhy;
299: d2=(user->left[j-ys+1] - x[j][0]) *rhx;
300: ft = ft+PetscSqrtReal(1.0 + d3*d3 + d2*d2);
301: }
302: }
303: if (ys==0) { /* bottom side */
304: for (i=xs; i<xs+xm; i++) {
305: d2=(user->bottom[i+1-xs]-user->bottom[i-xs+2])*rhx;
306: d3=(user->bottom[i-xs+1]-x[0][i])*rhy;
307: ft = ft+PetscSqrtReal(1.0 + d3*d3 + d2*d2);
308: }
309: }
311: if (xs+xm==mx) { /* right side */
312: for (j=ys; j< ys+ym; j++) {
313: d1=(x[j][mx-1] - user->right[j-ys+1])*rhx;
314: d4=(user->right[j-ys]-user->right[j-ys+1])*rhy;
315: ft = ft+PetscSqrtReal(1.0 + d1*d1 + d4*d4);
316: }
317: }
318: if (ys+ym==my) { /* top side */
319: for (i=xs; i<xs+xm; i++) {
320: d1=(x[my-1][i] - user->top[i-xs+1])*rhy;
321: d4=(user->top[i-xs+1] - user->top[i-xs])*rhx;
322: ft = ft+PetscSqrtReal(1.0 + d1*d1 + d4*d4);
323: }
324: }
326: if (ys==0 && xs==0) {
327: d1=(user->left[0]-user->left[1])/hy;
328: d2=(user->bottom[0]-user->bottom[1])*rhx;
329: ft +=PetscSqrtReal(1.0 + d1*d1 + d2*d2);
330: }
331: if (ys+ym == my && xs+xm == mx) {
332: d1=(user->right[ym+1] - user->right[ym])*rhy;
333: d2=(user->top[xm+1] - user->top[xm])*rhx;
334: ft +=PetscSqrtReal(1.0 + d1*d1 + d2*d2);
335: }
337: ft=ft*area;
338: MPI_Allreduce(&ft,fcn,1,MPIU_REAL,MPIU_SUM,MPI_COMM_WORLD);
340: /* Restore vectors */
341: DMDAVecRestoreArray(user->dm,localX,(void**)&x);
342: DMDAVecRestoreArray(user->dm,G,(void**)&g);
344: /* Scatter values to global vector */
345: DMRestoreLocalVector(user->dm,&localX);
346: PetscLogFlops(67.0*xm*ym);
347: return 0;
348: }
350: /* ------------------------------------------------------------------- */
351: /*
352: FormHessian - Evaluates Hessian matrix.
354: Input Parameters:
355: . tao - the Tao context
356: . x - input vector
357: . ptr - optional user-defined context, as set by TaoSetHessian()
359: Output Parameters:
360: . H - Hessian matrix
361: . Hpre - optionally different preconditioning matrix
362: . flg - flag indicating matrix structure
364: */
365: PetscErrorCode FormHessian(Tao tao,Vec X,Mat H, Mat Hpre, void *ptr)
366: {
367: AppCtx *user = (AppCtx *) ptr;
369: /* Evaluate the Hessian entries*/
370: QuadraticH(user,X,H);
371: return 0;
372: }
374: /* ------------------------------------------------------------------- */
375: /*
376: QuadraticH - Evaluates Hessian matrix.
378: Input Parameters:
379: . user - user-defined context, as set by TaoSetHessian()
380: . X - input vector
382: Output Parameter:
383: . H - Hessian matrix
384: */
385: PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian)
386: {
387: PetscInt i,j,k;
388: PetscInt mx=user->mx, my=user->my;
389: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
390: PetscReal hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
391: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
392: PetscReal hl,hr,ht,hb,hc,htl,hbr;
393: PetscReal **x, v[7];
394: MatStencil col[7],row;
395: Vec localX;
396: PetscBool assembled;
398: /* Get local mesh boundaries */
399: DMGetLocalVector(user->dm,&localX);
401: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
402: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
404: /* Scatter ghost points to local vector */
405: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
406: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
408: /* Get pointers to vector data */
409: DMDAVecGetArray(user->dm,localX,(void**)&x);
411: /* Initialize matrix entries to zero */
412: MatAssembled(Hessian,&assembled);
413: if (assembled) MatZeroEntries(Hessian);
415: /* Set various matrix options */
416: MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
418: /* Compute Hessian over the locally owned part of the mesh */
420: for (j=ys; j<ys+ym; j++) {
422: for (i=xs; i< xs+xm; i++) {
424: xc = x[j][i];
425: xlt=xrb=xl=xr=xb=xt=xc;
427: /* Left side */
428: if (i==0) {
429: xl = user->left[j-ys+1];
430: xlt = user->left[j-ys+2];
431: } else {
432: xl = x[j][i-1];
433: }
435: if (j==0) {
436: xb = user->bottom[i-xs+1];
437: xrb = user->bottom[i-xs+2];
438: } else {
439: xb = x[j-1][i];
440: }
442: if (i+1 == mx) {
443: xr = user->right[j-ys+1];
444: xrb = user->right[j-ys];
445: } else {
446: xr = x[j][i+1];
447: }
449: if (j+1==my) {
450: xt = user->top[i-xs+1];
451: xlt = user->top[i-xs];
452: }else {
453: xt = x[j+1][i];
454: }
456: if (i>0 && j+1<my) {
457: xlt = x[j+1][i-1];
458: }
459: if (j>0 && i+1<mx) {
460: xrb = x[j-1][i+1];
461: }
463: d1 = (xc-xl)/hx;
464: d2 = (xc-xr)/hx;
465: d3 = (xc-xt)/hy;
466: d4 = (xc-xb)/hy;
467: d5 = (xrb-xr)/hy;
468: d6 = (xrb-xb)/hx;
469: d7 = (xlt-xl)/hy;
470: d8 = (xlt-xt)/hx;
472: f1 = PetscSqrtReal(1.0 + d1*d1 + d7*d7);
473: f2 = PetscSqrtReal(1.0 + d1*d1 + d4*d4);
474: f3 = PetscSqrtReal(1.0 + d3*d3 + d8*d8);
475: f4 = PetscSqrtReal(1.0 + d3*d3 + d2*d2);
476: f5 = PetscSqrtReal(1.0 + d2*d2 + d5*d5);
477: f6 = PetscSqrtReal(1.0 + d4*d4 + d6*d6);
479: hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
480: (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
481: hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
482: (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
483: ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
484: (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
485: hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
486: (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);
488: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6);
489: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3);
491: hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
492: hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
493: (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
494: (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);
496: hl/=2.0; hr/=2.0; ht/=2.0; hb/=2.0; hbr/=2.0; htl/=2.0; hc/=2.0;
498: row.j = j; row.i = i;
499: k=0;
500: if (j>0) {
501: v[k]=hb;
502: col[k].j = j - 1; col[k].i = i;
503: k++;
504: }
506: if (j>0 && i < mx -1) {
507: v[k]=hbr;
508: col[k].j = j - 1; col[k].i = i+1;
509: k++;
510: }
512: if (i>0) {
513: v[k]= hl;
514: col[k].j = j; col[k].i = i-1;
515: k++;
516: }
518: v[k]= hc;
519: col[k].j = j; col[k].i = i;
520: k++;
522: if (i < mx-1) {
523: v[k]= hr;
524: col[k].j = j; col[k].i = i+1;
525: k++;
526: }
528: if (i>0 && j < my-1) {
529: v[k]= htl;
530: col[k].j = j+1; col[k].i = i-1;
531: k++;
532: }
534: if (j < my-1) {
535: v[k]= ht;
536: col[k].j = j+1; col[k].i = i;
537: k++;
538: }
540: /*
541: Set matrix values using local numbering, which was defined
542: earlier, in the main routine.
543: */
544: MatSetValuesStencil(Hessian,1,&row,k,col,v,INSERT_VALUES);
545: }
546: }
548: DMDAVecRestoreArray(user->dm,localX,(void**)&x);
549: DMRestoreLocalVector(user->dm,&localX);
551: MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY);
552: MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY);
554: PetscLogFlops(199.0*xm*ym);
555: return 0;
556: }
558: /* ------------------------------------------------------------------- */
559: /*
560: MSA_BoundaryConditions - Calculates the boundary conditions for
561: the region.
563: Input Parameter:
564: . user - user-defined application context
566: Output Parameter:
567: . user - user-defined application context
568: */
569: static PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
570: {
571: PetscInt i,j,k,limit=0,maxits=5;
572: PetscInt xs,ys,xm,ym,gxs,gys,gxm,gym;
573: PetscInt mx=user->mx,my=user->my;
574: PetscInt bsize=0, lsize=0, tsize=0, rsize=0;
575: PetscReal one=1.0, two=2.0, three=3.0, tol=1e-10;
576: PetscReal fnorm,det,hx,hy,xt=0,yt=0;
577: PetscReal u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
578: PetscReal b=-0.5, t=0.5, l=-0.5, r=0.5;
579: PetscReal *boundary;
580: PetscBool flg;
582: /* Get local mesh boundaries */
583: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
584: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
586: bsize=xm+2;
587: lsize=ym+2;
588: rsize=ym+2;
589: tsize=xm+2;
591: PetscMalloc1(bsize,&user->bottom);
592: PetscMalloc1(tsize,&user->top);
593: PetscMalloc1(lsize,&user->left);
594: PetscMalloc1(rsize,&user->right);
596: hx= (r-l)/(mx+1); hy=(t-b)/(my+1);
598: for (j=0; j<4; j++) {
599: if (j==0) {
600: yt=b;
601: xt=l+hx*xs;
602: limit=bsize;
603: boundary=user->bottom;
604: } else if (j==1) {
605: yt=t;
606: xt=l+hx*xs;
607: limit=tsize;
608: boundary=user->top;
609: } else if (j==2) {
610: yt=b+hy*ys;
611: xt=l;
612: limit=lsize;
613: boundary=user->left;
614: } else { /* if (j==3) */
615: yt=b+hy*ys;
616: xt=r;
617: limit=rsize;
618: boundary=user->right;
619: }
621: for (i=0; i<limit; i++) {
622: u1=xt;
623: u2=-yt;
624: for (k=0; k<maxits; k++) {
625: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
626: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
627: fnorm=PetscSqrtReal(nf1*nf1+nf2*nf2);
628: if (fnorm <= tol) break;
629: njac11=one+u2*u2-u1*u1;
630: njac12=two*u1*u2;
631: njac21=-two*u1*u2;
632: njac22=-one - u1*u1 + u2*u2;
633: det = njac11*njac22-njac21*njac12;
634: u1 = u1-(njac22*nf1-njac12*nf2)/det;
635: u2 = u2-(njac11*nf2-njac21*nf1)/det;
636: }
638: boundary[i]=u1*u1-u2*u2;
639: if (j==0 || j==1) {
640: xt=xt+hx;
641: } else { /* if (j==2 || j==3) */
642: yt=yt+hy;
643: }
645: }
647: }
649: /* Scale the boundary if desired */
650: if (1==1) {
651: PetscReal scl = 1.0;
653: PetscOptionsGetReal(NULL,NULL,"-bottom",&scl,&flg);
654: if (flg) {
655: for (i=0;i<bsize;i++) user->bottom[i]*=scl;
656: }
658: PetscOptionsGetReal(NULL,NULL,"-top",&scl,&flg);
659: if (flg) {
660: for (i=0;i<tsize;i++) user->top[i]*=scl;
661: }
663: PetscOptionsGetReal(NULL,NULL,"-right",&scl,&flg);
664: if (flg) {
665: for (i=0;i<rsize;i++) user->right[i]*=scl;
666: }
668: PetscOptionsGetReal(NULL,NULL,"-left",&scl,&flg);
669: if (flg) {
670: for (i=0;i<lsize;i++) user->left[i]*=scl;
671: }
672: }
673: return 0;
674: }
676: /* ------------------------------------------------------------------- */
677: /*
678: MSA_InitialPoint - Calculates the initial guess in one of three ways.
680: Input Parameters:
681: . user - user-defined application context
682: . X - vector for initial guess
684: Output Parameters:
685: . X - newly computed initial guess
686: */
687: static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
688: {
689: PetscInt start2=-1,i,j;
690: PetscReal start1=0;
691: PetscBool flg1,flg2;
693: PetscOptionsGetReal(NULL,NULL,"-start",&start1,&flg1);
694: PetscOptionsGetInt(NULL,NULL,"-random",&start2,&flg2);
696: if (flg1) { /* The zero vector is reasonable */
698: VecSet(X, start1);
700: } else if (flg2 && start2>0) { /* Try a random start between -0.5 and 0.5 */
701: PetscRandom rctx; PetscReal np5=-0.5;
703: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
704: VecSetRandom(X, rctx);
705: PetscRandomDestroy(&rctx);
706: VecShift(X, np5);
708: } else { /* Take an average of the boundary conditions */
709: PetscInt xs,xm,ys,ym;
710: PetscInt mx=user->mx,my=user->my;
711: PetscReal **x;
713: /* Get local mesh boundaries */
714: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
716: /* Get pointers to vector data */
717: DMDAVecGetArray(user->dm,X,(void**)&x);
719: /* Perform local computations */
720: for (j=ys; j<ys+ym; j++) {
721: for (i=xs; i< xs+xm; i++) {
722: 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;
723: }
724: }
725: DMDAVecRestoreArray(user->dm,X,(void**)&x);
726: PetscLogFlops(9.0*xm*ym);
727: }
728: return 0;
729: }
731: /*-----------------------------------------------------------------------*/
732: PetscErrorCode My_Monitor(Tao tao, void *ctx)
733: {
734: Vec X;
736: TaoGetSolution(tao,&X);
737: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
738: return 0;
739: }
741: /*TEST
743: build:
744: requires: !complex
746: test:
747: args: -tao_smonitor -tao_type lmvm -mx 10 -my 8 -tao_gatol 1.e-3
748: requires: !single
750: test:
751: suffix: 2
752: nsize: 2
753: args: -tao_smonitor -tao_type nls -tao_nls_ksp_max_it 15 -tao_gatol 1.e-4
754: filter: grep -v "nls ksp"
755: requires: !single
757: test:
758: suffix: 3
759: nsize: 3
760: args: -tao_smonitor -tao_type cg -tao_cg_type fr -mx 10 -my 10 -tao_gatol 1.e-3
761: requires: !single
763: test:
764: suffix: 5
765: nsize: 2
766: args: -tao_smonitor -tao_type bmrm -mx 10 -my 8 -tao_gatol 1.e-3
767: requires: !single
769: TEST*/