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: TaoSetInitialVector();
26: Routines: TaoSetObjectiveAndGradientRoutine();
27: Routines: TaoSetHessianRoutine(); 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: PetscErrorCode ierr; /* used to check for functions returning nonzeros */
59: PetscInt Nx, Ny; /* number of processors in x- and y- directions */
60: Vec x; /* solution, gradient vectors */
61: PetscBool flg, viewmat; /* flags */
62: PetscBool fddefault, fdcoloring; /* flags */
63: Tao tao; /* TAO solver context */
64: AppCtx user; /* user-defined work context */
65: ISColoring iscoloring;
66: MatFDColoring matfdcoloring;
68: /* Initialize TAO */
69: PetscInitialize(&argc, &argv,(char *)0,help);if (ierr) return ierr;
71: /* Specify dimension of the problem */
72: user.mx = 10; user.my = 10;
74: /* Check for any command line arguments that override defaults */
75: PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,&flg);
76: PetscOptionsGetInt(NULL,NULL,"-my",&user.my,&flg);
78: PetscPrintf(MPI_COMM_WORLD,"\n---- Minimum Surface Area Problem -----\n");
79: PetscPrintf(MPI_COMM_WORLD,"mx: %D my: %D \n\n",user.mx,user.my);
81: /* Let PETSc determine the vector distribution */
82: Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
84: /* Create distributed array (DM) to manage parallel grid and vectors */
85: 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);
86: DMSetFromOptions(user.dm);
87: DMSetUp(user.dm);
89: /* Create TAO solver and set desired solution method.*/
90: TaoCreate(PETSC_COMM_WORLD,&tao);
91: TaoSetType(tao,TAOCG);
93: /*
94: Extract global vector from DA for the vector of variables -- PETSC routine
95: Compute the initial solution -- application specific, see below
96: Set this vector for use by TAO -- TAO routine
97: */
98: DMCreateGlobalVector(user.dm,&x);
99: MSA_BoundaryConditions(&user);
100: MSA_InitialPoint(&user,x);
101: TaoSetInitialVector(tao,x);
103: /*
104: Initialize the Application context for use in function evaluations -- application specific, see below.
105: Set routines for function and gradient evaluation
106: */
107: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&user);
109: /*
110: Given the command line arguments, calculate the hessian with either the user-
111: provided function FormHessian, or the default finite-difference driven Hessian
112: functions
113: */
114: PetscOptionsHasName(NULL,NULL,"-tao_fddefault",&fddefault);
115: PetscOptionsHasName(NULL,NULL,"-tao_fdcoloring",&fdcoloring);
117: /*
118: Create a matrix data structure to store the Hessian and set
119: the Hessian evalution routine.
120: Set the matrix structure to be used for Hessian evalutions
121: */
122: DMCreateMatrix(user.dm,&user.H);
123: MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE);
125: if (fdcoloring) {
126: DMCreateColoring(user.dm,IS_COLORING_GLOBAL,&iscoloring);
127: MatFDColoringCreate(user.H,iscoloring,&matfdcoloring);
128: ISColoringDestroy(&iscoloring);
129: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormGradient,(void*)&user);
130: MatFDColoringSetFromOptions(matfdcoloring);
131: TaoSetHessianRoutine(tao,user.H,user.H,TaoDefaultComputeHessianColor,(void *)matfdcoloring);
132: } else if (fddefault) {
133: TaoSetHessianRoutine(tao,user.H,user.H,TaoDefaultComputeHessian,(void *)NULL);
134: } else {
135: TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void *)&user);
136: }
138: /*
139: If my_monitor option is in command line, then use the user-provided
140: monitoring function
141: */
142: PetscOptionsHasName(NULL,NULL,"-my_monitor",&viewmat);
143: if (viewmat) {
144: TaoSetMonitor(tao,My_Monitor,NULL,NULL);
145: }
147: /* Check for any tao command line options */
148: TaoSetFromOptions(tao);
150: /* SOLVE THE APPLICATION */
151: TaoSolve(tao);
153: TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);
155: /* Free TAO data structures */
156: TaoDestroy(&tao);
158: /* Free PETSc data structures */
159: VecDestroy(&x);
160: MatDestroy(&user.H);
161: if (fdcoloring) {
162: MatFDColoringDestroy(&matfdcoloring);
163: }
164: PetscFree(user.bottom);
165: PetscFree(user.top);
166: PetscFree(user.left);
167: PetscFree(user.right);
168: DMDestroy(&user.dm);
169: PetscFinalize();
170: return ierr;
171: }
173: PetscErrorCode FormGradient(Tao tao, Vec X, Vec G,void *userCtx)
174: {
176: PetscReal fcn;
179: FormFunctionGradient(tao,X,&fcn,G,userCtx);
180: return(0);
181: }
183: /* -------------------------------------------------------------------- */
184: /* FormFunctionGradient - Evaluates the function and corresponding gradient.
186: Input Parameters:
187: . tao - the Tao context
188: . XX - input vector
189: . userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()
191: Output Parameters:
192: . fcn - the newly evaluated function
193: . GG - vector containing the newly evaluated gradient
194: */
195: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn,Vec G,void *userCtx)
196: {
197: AppCtx *user = (AppCtx *) userCtx;
199: PetscInt i,j;
200: PetscInt mx=user->mx, my=user->my;
201: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
202: PetscReal ft=0.0;
203: PetscReal hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy;
204: PetscReal rhx=mx+1, rhy=my+1;
205: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
206: PetscReal df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
207: PetscReal **g, **x;
208: Vec localX;
211: /* Get local mesh boundaries */
212: DMGetLocalVector(user->dm,&localX);
213: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
214: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
216: /* Scatter ghost points to local vector */
217: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
218: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
220: /* Get pointers to vector data */
221: DMDAVecGetArray(user->dm,localX,(void**)&x);
222: DMDAVecGetArray(user->dm,G,(void**)&g);
224: /* Compute function and gradient over the locally owned part of the mesh */
225: for (j=ys; j<ys+ym; j++) {
226: for (i=xs; i< xs+xm; i++) {
228: xc = x[j][i];
229: xlt=xrb=xl=xr=xb=xt=xc;
231: if (i==0) { /* left side */
232: xl= user->left[j-ys+1];
233: xlt = user->left[j-ys+2];
234: } else {
235: xl = x[j][i-1];
236: }
238: if (j==0) { /* bottom side */
239: xb=user->bottom[i-xs+1];
240: xrb =user->bottom[i-xs+2];
241: } else {
242: xb = x[j-1][i];
243: }
245: if (i+1 == gxs+gxm) { /* right side */
246: xr=user->right[j-ys+1];
247: xrb = user->right[j-ys];
248: } else {
249: xr = x[j][i+1];
250: }
252: if (j+1==gys+gym) { /* top side */
253: xt=user->top[i-xs+1];
254: xlt = user->top[i-xs];
255: }else {
256: xt = x[j+1][i];
257: }
259: if (i>gxs && j+1<gys+gym) {
260: xlt = x[j+1][i-1];
261: }
262: if (j>gys && i+1<gxs+gxm) {
263: xrb = x[j-1][i+1];
264: }
266: d1 = (xc-xl);
267: d2 = (xc-xr);
268: d3 = (xc-xt);
269: d4 = (xc-xb);
270: d5 = (xr-xrb);
271: d6 = (xrb-xb);
272: d7 = (xlt-xl);
273: d8 = (xt-xlt);
275: df1dxc = d1*hydhx;
276: df2dxc = (d1*hydhx + d4*hxdhy);
277: df3dxc = d3*hxdhy;
278: df4dxc = (d2*hydhx + d3*hxdhy);
279: df5dxc = d2*hydhx;
280: df6dxc = d4*hxdhy;
282: d1 *= rhx;
283: d2 *= rhx;
284: d3 *= rhy;
285: d4 *= rhy;
286: d5 *= rhy;
287: d6 *= rhx;
288: d7 *= rhy;
289: d8 *= rhx;
291: f1 = PetscSqrtReal(1.0 + d1*d1 + d7*d7);
292: f2 = PetscSqrtReal(1.0 + d1*d1 + d4*d4);
293: f3 = PetscSqrtReal(1.0 + d3*d3 + d8*d8);
294: f4 = PetscSqrtReal(1.0 + d3*d3 + d2*d2);
295: f5 = PetscSqrtReal(1.0 + d2*d2 + d5*d5);
296: f6 = PetscSqrtReal(1.0 + d4*d4 + d6*d6);
298: ft = ft + (f2 + f4);
300: df1dxc /= f1;
301: df2dxc /= f2;
302: df3dxc /= f3;
303: df4dxc /= f4;
304: df5dxc /= f5;
305: df6dxc /= f6;
307: g[j][i] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc) * 0.5;
309: }
310: }
312: /* Compute triangular areas along the border of the domain. */
313: if (xs==0) { /* left side */
314: for (j=ys; j<ys+ym; j++) {
315: d3=(user->left[j-ys+1] - user->left[j-ys+2])*rhy;
316: d2=(user->left[j-ys+1] - x[j][0]) *rhx;
317: ft = ft+PetscSqrtReal(1.0 + d3*d3 + d2*d2);
318: }
319: }
320: if (ys==0) { /* bottom side */
321: for (i=xs; i<xs+xm; i++) {
322: d2=(user->bottom[i+1-xs]-user->bottom[i-xs+2])*rhx;
323: d3=(user->bottom[i-xs+1]-x[0][i])*rhy;
324: ft = ft+PetscSqrtReal(1.0 + d3*d3 + d2*d2);
325: }
326: }
328: if (xs+xm==mx) { /* right side */
329: for (j=ys; j< ys+ym; j++) {
330: d1=(x[j][mx-1] - user->right[j-ys+1])*rhx;
331: d4=(user->right[j-ys]-user->right[j-ys+1])*rhy;
332: ft = ft+PetscSqrtReal(1.0 + d1*d1 + d4*d4);
333: }
334: }
335: if (ys+ym==my) { /* top side */
336: for (i=xs; i<xs+xm; i++) {
337: d1=(x[my-1][i] - user->top[i-xs+1])*rhy;
338: d4=(user->top[i-xs+1] - user->top[i-xs])*rhx;
339: ft = ft+PetscSqrtReal(1.0 + d1*d1 + d4*d4);
340: }
341: }
343: if (ys==0 && xs==0) {
344: d1=(user->left[0]-user->left[1])/hy;
345: d2=(user->bottom[0]-user->bottom[1])*rhx;
346: ft +=PetscSqrtReal(1.0 + d1*d1 + d2*d2);
347: }
348: if (ys+ym == my && xs+xm == mx) {
349: d1=(user->right[ym+1] - user->right[ym])*rhy;
350: d2=(user->top[xm+1] - user->top[xm])*rhx;
351: ft +=PetscSqrtReal(1.0 + d1*d1 + d2*d2);
352: }
354: ft=ft*area;
355: MPI_Allreduce(&ft,fcn,1,MPIU_REAL,MPIU_SUM,MPI_COMM_WORLD);
357: /* Restore vectors */
358: DMDAVecRestoreArray(user->dm,localX,(void**)&x);
359: DMDAVecRestoreArray(user->dm,G,(void**)&g);
361: /* Scatter values to global vector */
362: DMRestoreLocalVector(user->dm,&localX);
363: PetscLogFlops(67.0*xm*ym);
364: return(0);
365: }
367: /* ------------------------------------------------------------------- */
368: /*
369: FormHessian - Evaluates Hessian matrix.
371: Input Parameters:
372: . tao - the Tao context
373: . x - input vector
374: . ptr - optional user-defined context, as set by TaoSetHessianRoutine()
376: Output Parameters:
377: . H - Hessian matrix
378: . Hpre - optionally different preconditioning matrix
379: . flg - flag indicating matrix structure
381: */
382: PetscErrorCode FormHessian(Tao tao,Vec X,Mat H, Mat Hpre, void *ptr)
383: {
385: AppCtx *user = (AppCtx *) ptr;
388: /* Evaluate the Hessian entries*/
389: QuadraticH(user,X,H);
390: return(0);
391: }
393: /* ------------------------------------------------------------------- */
394: /*
395: QuadraticH - Evaluates Hessian matrix.
397: Input Parameters:
398: . user - user-defined context, as set by TaoSetHessian()
399: . X - input vector
401: Output Parameter:
402: . H - Hessian matrix
403: */
404: PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian)
405: {
406: PetscErrorCode ierr;
407: PetscInt i,j,k;
408: PetscInt mx=user->mx, my=user->my;
409: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
410: PetscReal hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
411: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
412: PetscReal hl,hr,ht,hb,hc,htl,hbr;
413: PetscReal **x, v[7];
414: MatStencil col[7],row;
415: Vec localX;
416: PetscBool assembled;
419: /* Get local mesh boundaries */
420: DMGetLocalVector(user->dm,&localX);
422: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
423: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
425: /* Scatter ghost points to local vector */
426: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
427: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
429: /* Get pointers to vector data */
430: DMDAVecGetArray(user->dm,localX,(void**)&x);
432: /* Initialize matrix entries to zero */
433: MatAssembled(Hessian,&assembled);
434: if (assembled) {MatZeroEntries(Hessian);}
436: /* Set various matrix options */
437: MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
439: /* Compute Hessian over the locally owned part of the mesh */
441: for (j=ys; j<ys+ym; j++) {
443: for (i=xs; i< xs+xm; i++) {
445: xc = x[j][i];
446: xlt=xrb=xl=xr=xb=xt=xc;
448: /* Left side */
449: if (i==0) {
450: xl = user->left[j-ys+1];
451: xlt = user->left[j-ys+2];
452: } else {
453: xl = x[j][i-1];
454: }
456: if (j==0) {
457: xb = user->bottom[i-xs+1];
458: xrb = user->bottom[i-xs+2];
459: } else {
460: xb = x[j-1][i];
461: }
463: if (i+1 == mx) {
464: xr = user->right[j-ys+1];
465: xrb = user->right[j-ys];
466: } else {
467: xr = x[j][i+1];
468: }
470: if (j+1==my) {
471: xt = user->top[i-xs+1];
472: xlt = user->top[i-xs];
473: }else {
474: xt = x[j+1][i];
475: }
477: if (i>0 && j+1<my) {
478: xlt = x[j+1][i-1];
479: }
480: if (j>0 && i+1<mx) {
481: xrb = x[j-1][i+1];
482: }
484: d1 = (xc-xl)/hx;
485: d2 = (xc-xr)/hx;
486: d3 = (xc-xt)/hy;
487: d4 = (xc-xb)/hy;
488: d5 = (xrb-xr)/hy;
489: d6 = (xrb-xb)/hx;
490: d7 = (xlt-xl)/hy;
491: d8 = (xlt-xt)/hx;
493: f1 = PetscSqrtReal(1.0 + d1*d1 + d7*d7);
494: f2 = PetscSqrtReal(1.0 + d1*d1 + d4*d4);
495: f3 = PetscSqrtReal(1.0 + d3*d3 + d8*d8);
496: f4 = PetscSqrtReal(1.0 + d3*d3 + d2*d2);
497: f5 = PetscSqrtReal(1.0 + d2*d2 + d5*d5);
498: f6 = PetscSqrtReal(1.0 + d4*d4 + d6*d6);
500: hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
501: (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
502: hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
503: (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
504: ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
505: (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
506: hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
507: (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);
509: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6);
510: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3);
512: hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
513: hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
514: (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
515: (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);
517: hl/=2.0; hr/=2.0; ht/=2.0; hb/=2.0; hbr/=2.0; htl/=2.0; hc/=2.0;
519: row.j = j; row.i = i;
520: k=0;
521: if (j>0) {
522: v[k]=hb;
523: col[k].j = j - 1; col[k].i = i;
524: k++;
525: }
527: if (j>0 && i < mx -1) {
528: v[k]=hbr;
529: col[k].j = j - 1; col[k].i = i+1;
530: k++;
531: }
533: if (i>0) {
534: v[k]= hl;
535: col[k].j = j; col[k].i = i-1;
536: k++;
537: }
539: v[k]= hc;
540: col[k].j = j; col[k].i = i;
541: k++;
543: if (i < mx-1) {
544: v[k]= hr;
545: col[k].j = j; col[k].i = i+1;
546: k++;
547: }
549: if (i>0 && j < my-1) {
550: v[k]= htl;
551: col[k].j = j+1; col[k].i = i-1;
552: k++;
553: }
555: if (j < my-1) {
556: v[k]= ht;
557: col[k].j = j+1; col[k].i = i;
558: k++;
559: }
561: /*
562: Set matrix values using local numbering, which was defined
563: earlier, in the main routine.
564: */
565: MatSetValuesStencil(Hessian,1,&row,k,col,v,INSERT_VALUES);
566: }
567: }
569: DMDAVecRestoreArray(user->dm,localX,(void**)&x);
570: DMRestoreLocalVector(user->dm,&localX);
572: MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY);
573: MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY);
575: PetscLogFlops(199.0*xm*ym);
576: return(0);
577: }
579: /* ------------------------------------------------------------------- */
580: /*
581: MSA_BoundaryConditions - Calculates the boundary conditions for
582: the region.
584: Input Parameter:
585: . user - user-defined application context
587: Output Parameter:
588: . user - user-defined application context
589: */
590: static PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
591: {
593: PetscInt i,j,k,limit=0,maxits=5;
594: PetscInt xs,ys,xm,ym,gxs,gys,gxm,gym;
595: PetscInt mx=user->mx,my=user->my;
596: PetscInt bsize=0, lsize=0, tsize=0, rsize=0;
597: PetscReal one=1.0, two=2.0, three=3.0, tol=1e-10;
598: PetscReal fnorm,det,hx,hy,xt=0,yt=0;
599: PetscReal u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
600: PetscReal b=-0.5, t=0.5, l=-0.5, r=0.5;
601: PetscReal *boundary;
602: PetscBool flg;
605: /* Get local mesh boundaries */
606: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
607: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
609: bsize=xm+2;
610: lsize=ym+2;
611: rsize=ym+2;
612: tsize=xm+2;
614: PetscMalloc1(bsize,&user->bottom);
615: PetscMalloc1(tsize,&user->top);
616: PetscMalloc1(lsize,&user->left);
617: PetscMalloc1(rsize,&user->right);
619: hx= (r-l)/(mx+1); hy=(t-b)/(my+1);
621: for (j=0; j<4; j++) {
622: if (j==0) {
623: yt=b;
624: xt=l+hx*xs;
625: limit=bsize;
626: boundary=user->bottom;
627: } else if (j==1) {
628: yt=t;
629: xt=l+hx*xs;
630: limit=tsize;
631: boundary=user->top;
632: } else if (j==2) {
633: yt=b+hy*ys;
634: xt=l;
635: limit=lsize;
636: boundary=user->left;
637: } else { /* if (j==3) */
638: yt=b+hy*ys;
639: xt=r;
640: limit=rsize;
641: boundary=user->right;
642: }
644: for (i=0; i<limit; i++) {
645: u1=xt;
646: u2=-yt;
647: for (k=0; k<maxits; k++) {
648: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
649: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
650: fnorm=PetscSqrtReal(nf1*nf1+nf2*nf2);
651: if (fnorm <= tol) break;
652: njac11=one+u2*u2-u1*u1;
653: njac12=two*u1*u2;
654: njac21=-two*u1*u2;
655: njac22=-one - u1*u1 + u2*u2;
656: det = njac11*njac22-njac21*njac12;
657: u1 = u1-(njac22*nf1-njac12*nf2)/det;
658: u2 = u2-(njac11*nf2-njac21*nf1)/det;
659: }
661: boundary[i]=u1*u1-u2*u2;
662: if (j==0 || j==1) {
663: xt=xt+hx;
664: } else { /* if (j==2 || j==3) */
665: yt=yt+hy;
666: }
668: }
670: }
672: /* Scale the boundary if desired */
673: if (1==1) {
674: PetscReal scl = 1.0;
676: PetscOptionsGetReal(NULL,NULL,"-bottom",&scl,&flg);
677: if (flg) {
678: for (i=0;i<bsize;i++) user->bottom[i]*=scl;
679: }
681: PetscOptionsGetReal(NULL,NULL,"-top",&scl,&flg);
682: if (flg) {
683: for (i=0;i<tsize;i++) user->top[i]*=scl;
684: }
686: PetscOptionsGetReal(NULL,NULL,"-right",&scl,&flg);
687: if (flg) {
688: for (i=0;i<rsize;i++) user->right[i]*=scl;
689: }
691: PetscOptionsGetReal(NULL,NULL,"-left",&scl,&flg);
692: if (flg) {
693: for (i=0;i<lsize;i++) user->left[i]*=scl;
694: }
695: }
696: return(0);
697: }
699: /* ------------------------------------------------------------------- */
700: /*
701: MSA_InitialPoint - Calculates the initial guess in one of three ways.
703: Input Parameters:
704: . user - user-defined application context
705: . X - vector for initial guess
707: Output Parameters:
708: . X - newly computed initial guess
709: */
710: static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
711: {
713: PetscInt start2=-1,i,j;
714: PetscReal start1=0;
715: PetscBool flg1,flg2;
718: PetscOptionsGetReal(NULL,NULL,"-start",&start1,&flg1);
719: PetscOptionsGetInt(NULL,NULL,"-random",&start2,&flg2);
721: if (flg1) { /* The zero vector is reasonable */
723: VecSet(X, start1);
725: } else if (flg2 && start2>0) { /* Try a random start between -0.5 and 0.5 */
726: PetscRandom rctx; PetscReal np5=-0.5;
728: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
729: VecSetRandom(X, rctx);
730: PetscRandomDestroy(&rctx);
731: VecShift(X, np5);
733: } else { /* Take an average of the boundary conditions */
734: PetscInt xs,xm,ys,ym;
735: PetscInt mx=user->mx,my=user->my;
736: PetscReal **x;
738: /* Get local mesh boundaries */
739: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
741: /* Get pointers to vector data */
742: DMDAVecGetArray(user->dm,X,(void**)&x);
744: /* Perform local computations */
745: for (j=ys; j<ys+ym; j++) {
746: for (i=xs; i< xs+xm; i++) {
747: 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;
748: }
749: }
750: DMDAVecRestoreArray(user->dm,X,(void**)&x);
751: PetscLogFlops(9.0*xm*ym);
752: }
753: return(0);
754: }
756: /*-----------------------------------------------------------------------*/
757: PetscErrorCode My_Monitor(Tao tao, void *ctx)
758: {
760: Vec X;
763: TaoGetSolutionVector(tao,&X);
764: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
765: return(0);
766: }
768: /*TEST
770: build:
771: requires: !complex
773: test:
774: args: -tao_smonitor -tao_type lmvm -mx 10 -my 8 -tao_gatol 1.e-3
775: requires: !single
777: test:
778: suffix: 2
779: nsize: 2
780: args: -tao_smonitor -tao_type nls -tao_nls_ksp_max_it 15 -tao_gatol 1.e-4
781: filter: grep -v "nls ksp"
782: requires: !single
784: test:
785: suffix: 3
786: nsize: 3
787: args: -tao_smonitor -tao_type cg -tao_cg_type fr -mx 10 -my 10 -tao_gatol 1.e-3
788: requires: !single
790: test:
791: suffix: 5
792: nsize: 2
793: args: -tao_smonitor -tao_type bmrm -mx 10 -my 8 -tao_gatol 1.e-3
794: requires: !single
796: TEST*/