xref: /petsc/src/tao/unconstrained/tutorials/minsurf1.c (revision 40badf4fbc550ac1f60bd080eaff6de6d55b946d)
1 /* Program usage: mpiexec -n 1 minsurf1 [-help] [all TAO options] */
2 
3 /*  Include "petsctao.h" so we can use TAO solvers.  */
4 #include <petsctao.h>
5 
6 static char  help[] =
7 "This example demonstrates use of the TAO package to \n\
8 solve an unconstrained minimization problem.  This example is based on a \n\
9 problem from the MINPACK-2 test suite.  Given a rectangular 2-D domain and \n\
10 boundary values along the edges of the domain, the objective is to find the\n\
11 surface with the minimal area that satisfies the boundary conditions.\n\
12 The command line options are:\n\
13   -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
14   -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
15   -start <st>, where <st> =0 for zero vector, and an average of the boundary conditions otherwise \n\n";
16 
17 /*T
18    Concepts: TAO^Solving an unconstrained minimization problem
19    Routines: TaoCreate(); TaoSetType();
20    Routines: TaoSetSolution();
21    Routines: TaoSetObjectiveAndGradient();
22    Routines: TaoSetHessian(); TaoSetFromOptions();
23    Routines: TaoGetKSP(); TaoSolve();
24    Routines: TaoDestroy();
25    Processors: 1
26 T*/
27 
28 /*
29    User-defined application context - contains data needed by the
30    application-provided call-back routines, FormFunctionGradient()
31    and FormHessian().
32 */
33 typedef struct {
34   PetscInt    mx, my;                 /* discretization in x, y directions */
35   PetscReal   *bottom, *top, *left, *right;             /* boundary values */
36   Mat         H;
37 } AppCtx;
38 
39 /* -------- User-defined Routines --------- */
40 
41 static PetscErrorCode MSA_BoundaryConditions(AppCtx*);
42 static PetscErrorCode MSA_InitialPoint(AppCtx*,Vec);
43 static PetscErrorCode QuadraticH(AppCtx*,Vec,Mat);
44 PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*);
45 PetscErrorCode FormHessian(Tao,Vec,Mat,Mat,void*);
46 
47 int main(int argc, char **argv)
48 {
49   PetscErrorCode     ierr;              /* used to check for functions returning nonzeros */
50   PetscInt           N;                 /* Size of vector */
51   PetscMPIInt        size;              /* Number of processors */
52   Vec                x;                 /* solution, gradient vectors */
53   KSP                ksp;               /*  PETSc Krylov subspace method */
54   PetscBool          flg;               /* A return value when checking for user options */
55   Tao                tao;               /* Tao solver context */
56   AppCtx             user;              /* user-defined work context */
57 
58   /* Initialize TAO,PETSc */
59   ierr = PetscInitialize(&argc, &argv,(char *)0,help);if (ierr) return ierr;
60 
61   CHKERRMPI(MPI_Comm_size(MPI_COMM_WORLD,&size));
62   PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Incorrect number of processors");
63 
64   /* Specify default dimension of the problem */
65   user.mx = 4; user.my = 4;
66 
67   /* Check for any command line arguments that override defaults */
68   CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,&flg));
69   CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-my",&user.my,&flg));
70 
71   CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"\n---- Minimum Surface Area Problem -----\n"));
72   CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"mx: %D     my: %D   \n\n",user.mx,user.my));
73 
74   /* Calculate any derived values from parameters */
75   N    = user.mx*user.my;
76 
77   /* Create TAO solver and set desired solution method  */
78   CHKERRQ(TaoCreate(PETSC_COMM_SELF,&tao));
79   CHKERRQ(TaoSetType(tao,TAOLMVM));
80 
81   /* Initialize minsurf application data structure for use in the function evaluations  */
82   CHKERRQ(MSA_BoundaryConditions(&user));            /* Application specific routine */
83 
84   /*
85      Create a vector to hold the variables.  Compute an initial solution.
86      Set this vector, which will be used by TAO.
87   */
88   CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,N,&x));
89   CHKERRQ(MSA_InitialPoint(&user,x));                /* Application specific routine */
90   CHKERRQ(TaoSetSolution(tao,x));   /* A TAO routine                */
91 
92   /* Provide TAO routines for function, gradient, and Hessian evaluation */
93   CHKERRQ(TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void *)&user));
94 
95   /* Create a matrix data structure to store the Hessian.  This structure will be used by TAO */
96   CHKERRQ(MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,7,NULL,&(user.H)));
97   CHKERRQ(MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE));
98   CHKERRQ(TaoSetHessian(tao,user.H,user.H,FormHessian,(void *)&user));
99 
100   /* Check for any TAO command line options */
101   CHKERRQ(TaoSetFromOptions(tao));
102 
103   /* Limit the number of iterations in the KSP linear solver */
104   CHKERRQ(TaoGetKSP(tao,&ksp));
105   if (ksp) {
106     CHKERRQ(KSPSetTolerances(ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,user.mx*user.my));
107   }
108 
109   /* SOLVE THE APPLICATION */
110   CHKERRQ(TaoSolve(tao));
111 
112   CHKERRQ(TaoDestroy(&tao));
113   CHKERRQ(VecDestroy(&x));
114   CHKERRQ(MatDestroy(&user.H));
115   CHKERRQ(PetscFree(user.bottom));
116   CHKERRQ(PetscFree(user.top));
117   CHKERRQ(PetscFree(user.left));
118   CHKERRQ(PetscFree(user.right));
119 
120   ierr = PetscFinalize();
121   return ierr;
122 }
123 
124 /* -------------------------------------------------------------------- */
125 
126 /*  FormFunctionGradient - Evaluates function and corresponding gradient.
127 
128     Input Parameters:
129 .   tao     - the Tao context
130 .   X       - input vector
131 .   userCtx - optional user-defined context, as set by TaoSetFunctionGradient()
132 
133     Output Parameters:
134 .   fcn     - the newly evaluated function
135 .   G       - vector containing the newly evaluated gradient
136 */
137 PetscErrorCode FormFunctionGradient(Tao tao,Vec X,PetscReal *fcn,Vec G,void *userCtx)
138 {
139   AppCtx            *user = (AppCtx *) userCtx;
140   PetscInt          i,j,row;
141   PetscInt          mx=user->mx, my=user->my;
142   PetscReal         rhx=mx+1, rhy=my+1;
143   PetscReal         hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy, ft=0;
144   PetscReal         f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
145   PetscReal         df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
146   PetscReal         zero=0.0;
147   PetscScalar       *g;
148   const PetscScalar *x;
149 
150   PetscFunctionBeginUser;
151   CHKERRQ(VecSet(G, zero));
152 
153   CHKERRQ(VecGetArrayRead(X,&x));
154   CHKERRQ(VecGetArray(G,&g));
155 
156   /* Compute function over the locally owned part of the mesh */
157   for (j=0; j<my; j++) {
158     for (i=0; i< mx; i++) {
159       row=(j)*mx + (i);
160       xc = x[row];
161       xlt=xrb=xl=xr=xb=xt=xc;
162       if (i==0) { /* left side */
163         xl  = user->left[j+1];
164         xlt = user->left[j+2];
165       } else {
166         xl = x[row-1];
167       }
168 
169       if (j==0) { /* bottom side */
170         xb  = user->bottom[i+1];
171         xrb = user->bottom[i+2];
172       } else {
173         xb = x[row-mx];
174       }
175 
176       if (i+1 == mx) { /* right side */
177         xr  = user->right[j+1];
178         xrb = user->right[j];
179       } else {
180         xr = x[row+1];
181       }
182 
183       if (j+1==0+my) { /* top side */
184         xt  = user->top[i+1];
185         xlt = user->top[i];
186       }else {
187         xt = x[row+mx];
188       }
189 
190       if (i>0 && j+1<my) {
191         xlt = x[row-1+mx];
192       }
193       if (j>0 && i+1<mx) {
194         xrb = x[row+1-mx];
195       }
196 
197       d1 = (xc-xl);
198       d2 = (xc-xr);
199       d3 = (xc-xt);
200       d4 = (xc-xb);
201       d5 = (xr-xrb);
202       d6 = (xrb-xb);
203       d7 = (xlt-xl);
204       d8 = (xt-xlt);
205 
206       df1dxc = d1*hydhx;
207       df2dxc = (d1*hydhx + d4*hxdhy);
208       df3dxc = d3*hxdhy;
209       df4dxc = (d2*hydhx + d3*hxdhy);
210       df5dxc = d2*hydhx;
211       df6dxc = d4*hxdhy;
212 
213       d1 *= rhx;
214       d2 *= rhx;
215       d3 *= rhy;
216       d4 *= rhy;
217       d5 *= rhy;
218       d6 *= rhx;
219       d7 *= rhy;
220       d8 *= rhx;
221 
222       f1 = PetscSqrtScalar(1.0 + d1*d1 + d7*d7);
223       f2 = PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
224       f3 = PetscSqrtScalar(1.0 + d3*d3 + d8*d8);
225       f4 = PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
226       f5 = PetscSqrtScalar(1.0 + d2*d2 + d5*d5);
227       f6 = PetscSqrtScalar(1.0 + d4*d4 + d6*d6);
228 
229       ft = ft + (f2 + f4);
230 
231       df1dxc /= f1;
232       df2dxc /= f2;
233       df3dxc /= f3;
234       df4dxc /= f4;
235       df5dxc /= f5;
236       df6dxc /= f6;
237 
238       g[row] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc)/2.0;
239     }
240   }
241 
242   for (j=0; j<my; j++) {   /* left side */
243     d3 = (user->left[j+1] - user->left[j+2])*rhy;
244     d2 = (user->left[j+1] - x[j*mx])*rhx;
245     ft = ft+PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
246   }
247 
248   for (i=0; i<mx; i++) { /* bottom */
249     d2 = (user->bottom[i+1]-user->bottom[i+2])*rhx;
250     d3 = (user->bottom[i+1]-x[i])*rhy;
251     ft = ft+PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
252   }
253 
254   for (j=0; j< my; j++) { /* right side */
255     d1 = (x[(j+1)*mx-1]-user->right[j+1])*rhx;
256     d4 = (user->right[j]-user->right[j+1])*rhy;
257     ft = ft+PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
258   }
259 
260   for (i=0; i<mx; i++) { /* top side */
261     d1 = (x[(my-1)*mx + i] - user->top[i+1])*rhy;
262     d4 = (user->top[i+1] - user->top[i])*rhx;
263     ft = ft+PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
264   }
265 
266   /* Bottom left corner */
267   d1  = (user->left[0]-user->left[1])*rhy;
268   d2  = (user->bottom[0]-user->bottom[1])*rhx;
269   ft += PetscSqrtScalar(1.0 + d1*d1 + d2*d2);
270 
271   /* Top right corner */
272   d1  = (user->right[my+1] - user->right[my])*rhy;
273   d2  = (user->top[mx+1] - user->top[mx])*rhx;
274   ft += PetscSqrtScalar(1.0 + d1*d1 + d2*d2);
275 
276   (*fcn)=ft*area;
277 
278   /* Restore vectors */
279   CHKERRQ(VecRestoreArrayRead(X,&x));
280   CHKERRQ(VecRestoreArray(G,&g));
281   CHKERRQ(PetscLogFlops(67.0*mx*my));
282   PetscFunctionReturn(0);
283 }
284 
285 /* ------------------------------------------------------------------- */
286 /*
287    FormHessian - Evaluates the Hessian matrix.
288 
289    Input Parameters:
290 .  tao  - the Tao context
291 .  x    - input vector
292 .  ptr  - optional user-defined context, as set by TaoSetHessian()
293 
294    Output Parameters:
295 .  H    - Hessian matrix
296 .  Hpre - optionally different preconditioning matrix
297 .  flg  - flag indicating matrix structure
298 
299 */
300 PetscErrorCode FormHessian(Tao tao,Vec X,Mat H, Mat Hpre, void *ptr)
301 {
302   AppCtx         *user = (AppCtx *) ptr;
303 
304   PetscFunctionBeginUser;
305   /* Evaluate the Hessian entries*/
306   CHKERRQ(QuadraticH(user,X,H));
307   PetscFunctionReturn(0);
308 }
309 
310 /* ------------------------------------------------------------------- */
311 /*
312    QuadraticH - Evaluates the Hessian matrix.
313 
314    Input Parameters:
315 .  user - user-defined context, as set by TaoSetHessian()
316 .  X    - input vector
317 
318    Output Parameter:
319 .  H    - Hessian matrix
320 */
321 PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian)
322 {
323   PetscInt          i,j,k,row;
324   PetscInt          mx=user->mx, my=user->my;
325   PetscInt          col[7];
326   PetscReal         hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
327   PetscReal         rhx=mx+1, rhy=my+1;
328   PetscReal         f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
329   PetscReal         hl,hr,ht,hb,hc,htl,hbr;
330   const PetscScalar *x;
331   PetscReal         v[7];
332 
333   PetscFunctionBeginUser;
334   /* Get pointers to vector data */
335   CHKERRQ(VecGetArrayRead(X,&x));
336 
337   /* Initialize matrix entries to zero */
338   CHKERRQ(MatZeroEntries(Hessian));
339 
340   /* Set various matrix options */
341   CHKERRQ(MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE));
342 
343   /* Compute Hessian over the locally owned part of the mesh */
344   for (i=0; i< mx; i++) {
345     for (j=0; j<my; j++) {
346 
347       row=(j)*mx + (i);
348 
349       xc = x[row];
350       xlt=xrb=xl=xr=xb=xt=xc;
351 
352       /* Left side */
353       if (i==0) {
354         xl= user->left[j+1];
355         xlt = user->left[j+2];
356       } else {
357         xl = x[row-1];
358       }
359 
360       if (j==0) {
361         xb=user->bottom[i+1];
362         xrb = user->bottom[i+2];
363       } else {
364         xb = x[row-mx];
365       }
366 
367       if (i+1 == mx) {
368         xr=user->right[j+1];
369         xrb = user->right[j];
370       } else {
371         xr = x[row+1];
372       }
373 
374       if (j+1==my) {
375         xt=user->top[i+1];
376         xlt = user->top[i];
377       }else {
378         xt = x[row+mx];
379       }
380 
381       if (i>0 && j+1<my) {
382         xlt = x[row-1+mx];
383       }
384       if (j>0 && i+1<mx) {
385         xrb = x[row+1-mx];
386       }
387 
388       d1 = (xc-xl)*rhx;
389       d2 = (xc-xr)*rhx;
390       d3 = (xc-xt)*rhy;
391       d4 = (xc-xb)*rhy;
392       d5 = (xrb-xr)*rhy;
393       d6 = (xrb-xb)*rhx;
394       d7 = (xlt-xl)*rhy;
395       d8 = (xlt-xt)*rhx;
396 
397       f1 = PetscSqrtScalar(1.0 + d1*d1 + d7*d7);
398       f2 = PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
399       f3 = PetscSqrtScalar(1.0 + d3*d3 + d8*d8);
400       f4 = PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
401       f5 = PetscSqrtScalar(1.0 + d2*d2 + d5*d5);
402       f6 = PetscSqrtScalar(1.0 + d4*d4 + d6*d6);
403 
404       hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+(-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
405       hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+(-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
406       ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+(-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
407       hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+(-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);
408 
409       hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6);
410       htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3);
411 
412       hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) + hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
413            (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +  (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);
414 
415       hl*=0.5; hr*=0.5; ht*=0.5; hb*=0.5; hbr*=0.5; htl*=0.5;  hc*=0.5;
416 
417       k=0;
418       if (j>0) {
419         v[k]=hb; col[k]=row - mx; k++;
420       }
421 
422       if (j>0 && i < mx -1) {
423         v[k]=hbr; col[k]=row - mx+1; k++;
424       }
425 
426       if (i>0) {
427         v[k]= hl; col[k]=row - 1; k++;
428       }
429 
430       v[k]= hc; col[k]=row; k++;
431 
432       if (i < mx-1) {
433         v[k]= hr; col[k]=row+1; k++;
434       }
435 
436       if (i>0 && j < my-1) {
437         v[k]= htl; col[k] = row+mx-1; k++;
438       }
439 
440       if (j < my-1) {
441         v[k]= ht; col[k] = row+mx; k++;
442       }
443 
444       /*
445          Set matrix values using local numbering, which was defined
446          earlier, in the main routine.
447       */
448       CHKERRQ(MatSetValues(Hessian,1,&row,k,col,v,INSERT_VALUES));
449     }
450   }
451 
452   /* Restore vectors */
453   CHKERRQ(VecRestoreArrayRead(X,&x));
454 
455   /* Assemble the matrix */
456   CHKERRQ(MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY));
457   CHKERRQ(MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY));
458 
459   CHKERRQ(PetscLogFlops(199.0*mx*my));
460   PetscFunctionReturn(0);
461 }
462 
463 /* ------------------------------------------------------------------- */
464 /*
465    MSA_BoundaryConditions -  Calculates the boundary conditions for
466    the region.
467 
468    Input Parameter:
469 .  user - user-defined application context
470 
471    Output Parameter:
472 .  user - user-defined application context
473 */
474 static PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
475 {
476   PetscInt       i,j,k,limit=0;
477   PetscInt       maxits=5;
478   PetscInt       mx=user->mx,my=user->my;
479   PetscInt       bsize=0, lsize=0, tsize=0, rsize=0;
480   PetscReal      one=1.0, two=2.0, three=3.0, tol=1e-10;
481   PetscReal      fnorm,det,hx,hy,xt=0,yt=0;
482   PetscReal      u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
483   PetscReal      b=-0.5, t=0.5, l=-0.5, r=0.5;
484   PetscReal      *boundary;
485 
486   PetscFunctionBeginUser;
487   bsize=mx+2; lsize=my+2; rsize=my+2; tsize=mx+2;
488 
489   CHKERRQ(PetscMalloc1(bsize,&user->bottom));
490   CHKERRQ(PetscMalloc1(tsize,&user->top));
491   CHKERRQ(PetscMalloc1(lsize,&user->left));
492   CHKERRQ(PetscMalloc1(rsize,&user->right));
493 
494   hx= (r-l)/(mx+1); hy=(t-b)/(my+1);
495 
496   for (j=0; j<4; j++) {
497     if (j==0) {
498       yt=b;
499       xt=l;
500       limit=bsize;
501       boundary=user->bottom;
502     } else if (j==1) {
503       yt=t;
504       xt=l;
505       limit=tsize;
506       boundary=user->top;
507     } else if (j==2) {
508       yt=b;
509       xt=l;
510       limit=lsize;
511       boundary=user->left;
512     } else {  /*  if (j==3) */
513       yt=b;
514       xt=r;
515       limit=rsize;
516       boundary=user->right;
517     }
518 
519     for (i=0; i<limit; i++) {
520       u1=xt;
521       u2=-yt;
522       for (k=0; k<maxits; k++) {
523         nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
524         nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
525         fnorm=PetscSqrtScalar(nf1*nf1+nf2*nf2);
526         if (fnorm <= tol) break;
527         njac11=one+u2*u2-u1*u1;
528         njac12=two*u1*u2;
529         njac21=-two*u1*u2;
530         njac22=-one - u1*u1 + u2*u2;
531         det = njac11*njac22-njac21*njac12;
532         u1 = u1-(njac22*nf1-njac12*nf2)/det;
533         u2 = u2-(njac11*nf2-njac21*nf1)/det;
534       }
535 
536       boundary[i]=u1*u1-u2*u2;
537       if (j==0 || j==1) {
538         xt=xt+hx;
539       } else { /*  if (j==2 || j==3) */
540         yt=yt+hy;
541       }
542     }
543   }
544   PetscFunctionReturn(0);
545 }
546 
547 /* ------------------------------------------------------------------- */
548 /*
549    MSA_InitialPoint - Calculates the initial guess in one of three ways.
550 
551    Input Parameters:
552 .  user - user-defined application context
553 .  X - vector for initial guess
554 
555    Output Parameters:
556 .  X - newly computed initial guess
557 */
558 static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
559 {
560   PetscInt       start=-1,i,j;
561   PetscReal      zero=0.0;
562   PetscBool      flg;
563 
564   PetscFunctionBeginUser;
565   CHKERRQ(VecSet(X, zero));
566   CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-start",&start,&flg));
567 
568   if (flg && start==0) { /* The zero vector is reasonable */
569      CHKERRQ(VecSet(X, zero));
570    } else { /* Take an average of the boundary conditions */
571     PetscInt    row;
572     PetscInt    mx=user->mx,my=user->my;
573     PetscReal *x;
574 
575     /* Get pointers to vector data */
576     CHKERRQ(VecGetArray(X,&x));
577     /* Perform local computations */
578     for (j=0; j<my; j++) {
579       for (i=0; i< mx; i++) {
580         row=(j)*mx + (i);
581         x[row] = (((j+1)*user->bottom[i+1]+(my-j+1)*user->top[i+1])/(my+2)+ ((i+1)*user->left[j+1]+(mx-i+1)*user->right[j+1])/(mx+2))/2.0;
582       }
583     }
584     /* Restore vectors */
585     CHKERRQ(VecRestoreArray(X,&x));
586   }
587   PetscFunctionReturn(0);
588 }
589 
590 /*TEST
591 
592    build:
593       requires: !complex
594 
595    test:
596       args: -tao_smonitor -tao_type nls -mx 10 -my 8 -tao_gatol 1.e-4
597       requires: !single
598 
599    test:
600       suffix: 2
601       args: -tao_smonitor -tao_type bmrm -mx 10 -my 8 -tao_gatol 1.e-3
602       requires: !single
603 
604    test:
605       suffix: 3
606       args: -tao_smonitor -tao_type lmvm -mx 10 -my 8 -tao_gatol 1.e-3
607       requires: !single
608 
609    test:
610       suffix: 4
611       args: -tao_smonitor -tao_type bntr -mx 10 -my 8 -tao_gatol 1.e-4
612       requires: !single
613 
614    test:
615       suffix: 5
616       args: -tao_smonitor -tao_type bntl -mx 10 -my 8 -tao_gatol 1.e-4
617       requires: !single
618 
619    test:
620       suffix: 6
621       args: -tao_smonitor -tao_type bnls -mx 10 -my 8 -tao_gatol 1.e-4
622       requires: !single
623 
624    test:
625       suffix: 7
626       args: -tao_smonitor -tao_type bntr -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4
627       requires: !single
628 
629    test:
630       suffix: 8
631       args: -tao_smonitor -tao_type bntl -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4
632       requires: !single
633 
634    test:
635       suffix: 9
636       args: -tao_smonitor -tao_type bnls -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4
637       requires: !single
638 
639    test:
640       suffix: 10
641       args: -tao_smonitor -tao_type bnls -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 -tao_mf_hessian
642 
643    test:
644       suffix: 11
645       args: -tao_smonitor -tao_type bntr -mx 10 -my 8 -tao_gatol 1.e-4 -tao_mf_hessian
646       requires: !single
647 
648    test:
649       suffix: 12
650       args: -tao_smonitor -tao_type bntl -mx 10 -my 8 -tao_gatol 1.e-4 -tao_mf_hessian
651       requires: !single
652 
653 TEST*/
654