Actual source code: minsurf2.c

  1: /*$Id$*/

  3: /* Program usage: mpirun -np <proc> minsurf2 [-help] [all TAO options] */

  5: /*
  6:   Include "tao.h" so we can use TAO solvers.
  7:   petscda.h for distributed array
  8: */
  9: #include "petscda.h"
 10: #include "tao.h"

 12: static  char help[] = 
 13: "This example demonstrates use of the TAO package to \n\
 14: solve an unconstrained minimization problem.  This example is based on a \n\
 15: problem from the MINPACK-2 test suite.  Given a rectangular 2-D domain and \n\
 16: boundary values along the edges of the domain, the objective is to find the\n\
 17: surface with the minimal area that satisfies the boundary conditions.\n\
 18: The command line options are:\n\
 19:   -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
 20:   -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
 21:   -start <st>, where <st> =0 for zero vector, <st> >0 for random start, and <st> <0 \n\
 22:                for an average of the boundary conditions\n\n";

 24: /*T
 25:    Concepts: TAO - Solving an unconstrained minimization problem
 26:    Routines: TaoInitialize(); TaoFinalize();
 27:    Routines: TaoCreate(); TaoDestroy();
 28:    Routines: TaoApplicationCreate(); TaoAppDestroy();
 29:    Routines: TaoAppSetInitialSolutionVec();
 30:    Routines: TaoAppSetObjectiveAndGradientRoutine();
 31:    Routines: TaoAppSetHessianMat(); TaoAppSetHessianRoutine();
 32:    Routines: TaoSetOptions();
 33:    Routines: TaoAppGetKSP(); TaoSolveApplication();
 34:    Routines: TaoAppSetMonitor(); TaoView();
 35:    Routines: TaoAppGetSolutionVec();
 36:    Processors: 1
 37: T*/

 39: /* 
 40:    User-defined application context - contains data needed by the 
 41:    application-provided call-back routines, FormFunctionGradient() 
 42:    and FormHessian().
 43: */
 44: typedef struct {
 45:   PetscInt      mx, my;                 /* discretization in x, y directions */
 46:   double      *bottom, *top, *left, *right;             /* boundary values */
 47:   DA          da;                      /* distributed array data structure */
 48:   Mat         H;                       /* Hessian */
 49:   ISColoring  iscoloring;
 50: } AppCtx;


 53: /* -------- User-defined Routines --------- */

 55: static int MSA_BoundaryConditions(AppCtx*);
 56: static int MSA_InitialPoint(AppCtx*,Vec);
 57: int QuadraticH(AppCtx*,Vec,Mat);
 58: int FormFunctionGradient(TAO_APPLICATION,Vec,double *,Vec,void*);
 59: int FormGradient(TAO_APPLICATION,Vec,Vec,void*);
 60: int FormHessian(TAO_APPLICATION,Vec,Mat*,Mat*,MatStructure *,void*);
 61: int My_Monitor(TAO_APPLICATION, void *);

 65: int main( int argc, char **argv )
 66: {
 67:   int             info;                /* used to check for functions returning nonzeros */
 68:   PetscInt          Nx, Ny;              /* number of processors in x- and y- directions */
 69:   PetscInt          iter;                /* iteration information */
 70:   double          ff,gnorm;
 71:   Vec             x;                   /* solution, gradient vectors */
 72:   PetscTruth      flg, viewmat;        /* flags */
 73:   PetscTruth      fddefault, fdcoloring;   /* flags */
 74:   KSP             ksp;                 /* Krylov subspace method */
 75:   TaoMethod       method = "tao_nls";  /* minimization method */
 76:   TaoTerminateReason reason;           
 77:   TAO_SOLVER      tao;                 /* TAO_SOLVER solver context */
 78:   TAO_APPLICATION minsurfapp;          /* The PETSc application */
 79:   AppCtx          user;                /* user-defined work context */

 81:   /* Initialize TAO */
 82:   PetscInitialize( &argc, &argv,(char *)0,help );
 83:   TaoInitialize( &argc, &argv,(char *)0,help );

 85:   /* Specify dimension of the problem */
 86:   user.mx = 10; user.my = 10;

 88:   /* Check for any command line arguments that override defaults */
 89:   info = PetscOptionsGetInt(PETSC_NULL,"-mx",&user.mx,&flg); CHKERRQ(info);
 90:   info = PetscOptionsGetInt(PETSC_NULL,"-my",&user.my,&flg); CHKERRQ(info);

 92:   PetscPrintf(MPI_COMM_WORLD,"\n---- Minimum Surface Area Problem -----\n");
 93:   PetscPrintf(MPI_COMM_WORLD,"mx: %d     my: %d   \n\n",user.mx,user.my);


 96:   /* Let PETSc determine the vector distribution */
 97:   Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;

 99:   /* Create distributed array (DA) to manage parallel grid and vectors  */
100:   info = DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_BOX,user.mx,
101:                     user.my,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,&user.da); CHKERRQ(info);
102:   

104:   /* Create TAO solver and set desired solution method. Create an TAO application structure */
105:   info = TaoCreate(PETSC_COMM_WORLD,method,&tao); CHKERRQ(info);
106:   info = TaoApplicationCreate(PETSC_COMM_WORLD,&minsurfapp); CHKERRQ(info);

108:   /*
109:      Extract global vector from DA for the vector of variables --  PETSC routine
110:      Compute the initial solution                              --  application specific, see below
111:      Set this vector for use by TAO                            --  TAO routine
112:   */
113:   info = DACreateGlobalVector(user.da,&x); CHKERRQ(info);
114:   info = MSA_BoundaryConditions(&user); CHKERRQ(info);         
115:   info = MSA_InitialPoint(&user,x); CHKERRQ(info);
116:   info = TaoAppSetInitialSolutionVec(minsurfapp,x); CHKERRQ(info);

118:   /* 
119:      Initialize the Application context for use in function evaluations  --  application specific, see below.
120:      Set routines for function and gradient evaluation 
121:   */
122:   info = TaoAppSetObjectiveAndGradientRoutine(minsurfapp,FormFunctionGradient,(void *)&user); CHKERRQ(info);

124:   /* 
125:      Given the command line arguments, calculate the hessian with either the user-
126:      provided function FormHessian, or the default finite-difference driven Hessian
127:      functions 
128:   */
129:   info = PetscOptionsHasName(PETSC_NULL,"-tao_fddefault",&fddefault);CHKERRQ(info);
130:   info = PetscOptionsHasName(PETSC_NULL,"-tao_fdcoloring",&fdcoloring);CHKERRQ(info);


133:   /* 
134:      Create a matrix data structure to store the Hessian and set 
135:      the Hessian evalution routine.
136:      Set the matrix structure to be used for Hessian evalutions
137:   */
138:   info = DAGetMatrix(user.da,MATAIJ,&user.H);CHKERRQ(info);
139:   info = MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE); CHKERRQ(info);

141:   info = TaoAppSetHessianMat(minsurfapp,user.H,user.H); CHKERRQ(info);

143:   if (fdcoloring) {
144:     info = DAGetColoring(user.da,IS_COLORING_GLOBAL,MATAIJ,&(user.iscoloring));
145:     CHKERRQ(info);
146:     info = TaoAppSetColoring(minsurfapp, user.iscoloring); CHKERRQ(info);
147:     info = TaoAppSetHessianRoutine(minsurfapp,TaoAppDefaultComputeHessianColor,(void *)&user); CHKERRQ(info);
148:   } else if (fddefault){
149:     info = TaoAppSetHessianRoutine(minsurfapp,TaoAppDefaultComputeHessian,(void *)&user); CHKERRQ(info);
150:   } else {
151:     info = TaoAppSetHessianRoutine(minsurfapp,FormHessian,(void *)&user); CHKERRQ(info);
152:   }


155:   /* 
156:      If my_monitor option is in command line, then use the user-provided
157:      monitoring function
158:   */
159:   info = PetscOptionsHasName(PETSC_NULL,"-my_monitor",&viewmat); CHKERRQ(info);
160:   if (viewmat){
161:     info = TaoAppSetMonitor(minsurfapp,My_Monitor,TAO_NULL); CHKERRQ(info);
162:   }

164:   /* Check for any tao command line options */
165:   info = TaoSetOptions(minsurfapp,tao); CHKERRQ(info);

167:   /* Limit the number of iterations in the KSP linear solver */
168:   info = TaoAppGetKSP(minsurfapp,&ksp); CHKERRQ(info);
169:   if (ksp) {                                              /* Modify the PETSc KSP structure */
170:     info = KSPSetTolerances(ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,user.mx*user.my);
171:     CHKERRQ(info);
172:   }

174:   /* SOLVE THE APPLICATION */
175:   info = TaoSolveApplication(minsurfapp,tao);  CHKERRQ(info);

177:   /* Get information on termination */
178:   info = TaoGetSolutionStatus(tao,&iter,&ff,&gnorm,0,0,&reason); CHKERRQ(info);
179:   if (reason <= 0 ){
180:     PetscPrintf(MPI_COMM_WORLD,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n");
181:     PetscPrintf(MPI_COMM_WORLD," Iterations: %d,  Function Value: %4.2e, Residual: %4.2e \n",iter,ff,gnorm);
182:   }

184:   /* 
185:      To View TAO solver information use
186:       info = TaoView(tao); CHKERRQ(info); 
187:   */

189:   /* Free TAO data structures */
190:   info = TaoDestroy(tao); CHKERRQ(info);
191:   info = TaoAppDestroy(minsurfapp); CHKERRQ(info);

193:   /* Free PETSc data structures */
194:   info = VecDestroy(x); CHKERRQ(info);
195:   info = MatDestroy(user.H); CHKERRQ(info);
196:   if (fdcoloring) {
197:     info = ISColoringDestroy(user.iscoloring);CHKERRQ(info);
198:   }
199:   info = PetscFree(user.bottom); CHKERRQ(info);
200:   info = PetscFree(user.top); CHKERRQ(info);
201:   info = PetscFree(user.left); CHKERRQ(info);
202:   info = PetscFree(user.right); CHKERRQ(info);
203:   info = DADestroy(user.da); CHKERRQ(info);

205:   /* Finalize TAO */
206:   TaoFinalize();
207:   PetscFinalize();
208:   
209:   return 0;
210: }

214: int FormGradient(TAO_APPLICATION taoapp, Vec X, Vec G,void *userCtx){
215:   int info;
216:   double fcn;
217:   TaoFunctionBegin;
218:   info = FormFunctionGradient(taoapp,X,&fcn,G,userCtx);CHKERRQ(info);
219:   TaoFunctionReturn(0);
220: }

222: /* -------------------------------------------------------------------- */
225: /*  FormFunctionGradient - Evaluates the function and corresponding gradient.

227:     Input Parameters:
228: .   taoapp     - the TAO_APPLICATION context
229: .   XX      - input vector
230: .   userCtx - optional user-defined context, as set by TaoSetFunctionGradient()
231:     
232:     Output Parameters:
233: .   fcn     - the newly evaluated function
234: .   GG       - vector containing the newly evaluated gradient
235: */
236: int FormFunctionGradient(TAO_APPLICATION taoapp, Vec X, double *fcn,Vec G,void *userCtx){

238:   AppCtx * user = (AppCtx *) userCtx;
239:   int    info;
240:   PetscInt i,j;
241:   PetscInt mx=user->mx, my=user->my;
242:   PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
243:   double ft=0;
244:   double hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy;
245:   double rhx=mx+1, rhy=my+1;
246:   double f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
247:   double df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
248:   PetscScalar **g, **x;
249:   Vec    localX;

251:   /* Get local mesh boundaries */
252:   info = DAGetLocalVector(user->da,&localX);CHKERRQ(info);

254:   info = DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);
255:   info = DAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); CHKERRQ(info);

257:   /* Scatter ghost points to local vector */
258:   info = DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX); CHKERRQ(info);
259:   info = DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX); CHKERRQ(info);

261:   /* Get pointers to vector data */
262:   info = DAVecGetArray(user->da,localX,(void**)&x);
263:   info = DAVecGetArray(user->da,G,(void**)&g);

265:   /* Compute function and gradient over the locally owned part of the mesh */
266:   for (j=ys; j<ys+ym; j++){
267:     for (i=xs; i< xs+xm; i++){
268:       
269:       xc = x[j][i];
270:       xlt=xrb=xl=xr=xb=xt=xc;
271:       
272:       if (i==0){ /* left side */
273:         xl= user->left[j-ys+1];
274:         xlt = user->left[j-ys+2];
275:       } else {
276:         xl = x[j][i-1];
277:       }

279:       if (j==0){ /* bottom side */
280:         xb=user->bottom[i-xs+1];
281:         xrb =user->bottom[i-xs+2];
282:       } else {
283:         xb = x[j-1][i];
284:       }
285:       
286:       if (i+1 == gxs+gxm){ /* right side */
287:         xr=user->right[j-ys+1];
288:         xrb = user->right[j-ys];
289:       } else {
290:         xr = x[j][i+1];
291:       }

293:       if (j+1==gys+gym){ /* top side */
294:         xt=user->top[i-xs+1];
295:         xlt = user->top[i-xs];
296:       }else {
297:         xt = x[j+1][i];
298:       }

300:       if (i>gxs && j+1<gys+gym){
301:         xlt = x[j+1][i-1];
302:       }
303:       if (j>gys && i+1<gxs+gxm){
304:         xrb = x[j-1][i+1];
305:       }

307:       d1 = (xc-xl);
308:       d2 = (xc-xr);
309:       d3 = (xc-xt);
310:       d4 = (xc-xb);
311:       d5 = (xr-xrb);
312:       d6 = (xrb-xb);
313:       d7 = (xlt-xl);
314:       d8 = (xt-xlt);
315:       
316:       df1dxc = d1*hydhx;
317:       df2dxc = ( d1*hydhx + d4*hxdhy );
318:       df3dxc = d3*hxdhy;
319:       df4dxc = ( d2*hydhx + d3*hxdhy );
320:       df5dxc = d2*hydhx;
321:       df6dxc = d4*hxdhy;

323:       d1 *= rhx;
324:       d2 *= rhx;
325:       d3 *= rhy;
326:       d4 *= rhy;
327:       d5 *= rhy;
328:       d6 *= rhx;
329:       d7 *= rhy;
330:       d8 *= rhx;

332:       f1 = sqrt( 1.0 + d1*d1 + d7*d7);
333:       f2 = sqrt( 1.0 + d1*d1 + d4*d4);
334:       f3 = sqrt( 1.0 + d3*d3 + d8*d8);
335:       f4 = sqrt( 1.0 + d3*d3 + d2*d2);
336:       f5 = sqrt( 1.0 + d2*d2 + d5*d5);
337:       f6 = sqrt( 1.0 + d4*d4 + d6*d6);
338:       
339:       f2 = sqrt( 1.0 + d1*d1 + d4*d4);
340:       f4 = sqrt( 1.0 + d3*d3 + d2*d2);

342:       ft = ft + (f2 + f4);

344:       df1dxc /= f1;
345:       df2dxc /= f2;
346:       df3dxc /= f3;
347:       df4dxc /= f4;
348:       df5dxc /= f5;
349:       df6dxc /= f6;

351:       g[j][i] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc ) * 0.5;
352:       
353:     }
354:   }

356:   /* Compute triangular areas along the border of the domain. */
357:   if (xs==0){ /* left side */
358:     for (j=ys; j<ys+ym; j++){
359:       d3=(user->left[j-ys+1] - user->left[j-ys+2])*rhy;
360:       d2=(user->left[j-ys+1] - x[j][0]) *rhx;
361:       ft = ft+sqrt( 1.0 + d3*d3 + d2*d2);
362:     }
363:   }
364:   if (ys==0){ /* bottom side */
365:     for (i=xs; i<xs+xm; i++){
366:       d2=(user->bottom[i+1-xs]-user->bottom[i-xs+2])*rhx;
367:       d3=(user->bottom[i-xs+1]-x[0][i])*rhy;
368:       ft = ft+sqrt( 1.0 + d3*d3 + d2*d2);
369:     }
370:   }

372:   if (xs+xm==mx){ /* right side */
373:     for (j=ys; j< ys+ym; j++){
374:       d1=(x[j][mx-1] - user->right[j-ys+1])*rhx;
375:       d4=(user->right[j-ys]-user->right[j-ys+1])*rhy;
376:       ft = ft+sqrt( 1.0 + d1*d1 + d4*d4);
377:     }
378:   }
379:   if (ys+ym==my){ /* top side */
380:     for (i=xs; i<xs+xm; i++){
381:       d1=(x[my-1][i] - user->top[i-xs+1])*rhy;
382:       d4=(user->top[i-xs+1] - user->top[i-xs])*rhx;
383:       ft = ft+sqrt( 1.0 + d1*d1 + d4*d4);
384:     }
385:   }

387:   if (ys==0 && xs==0){
388:     d1=(user->left[0]-user->left[1])/hy;
389:     d2=(user->bottom[0]-user->bottom[1])*rhx;
390:     ft +=sqrt( 1.0 + d1*d1 + d2*d2);
391:   }
392:   if (ys+ym == my && xs+xm == mx){
393:     d1=(user->right[ym+1] - user->right[ym])*rhy;
394:     d2=(user->top[xm+1] - user->top[xm])*rhx;
395:     ft +=sqrt( 1.0 + d1*d1 + d2*d2);
396:   }

398:   ft=ft*area;
399:   info = MPI_Allreduce(&ft,fcn,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);CHKERRQ(info);

401:   /* Restore vectors */
402:   info = DAVecRestoreArray(user->da,localX,(void**)&x);
403:   info = DAVecRestoreArray(user->da,G,(void**)&g);

405:   /* Scatter values to global vector */
406:   info = DARestoreLocalVector(user->da,&localX); CHKERRQ(info);

408:   info = PetscLogFlops(67*xm*ym); CHKERRQ(info);

410:   return 0;
411: }

413: /* ------------------------------------------------------------------- */
416: /*
417:    FormHessian - Evaluates Hessian matrix.

419:    Input Parameters:
420: .  taoapp  - the TAO_APPLICATION context
421: .  x    - input vector
422: .  ptr  - optional user-defined context, as set by TaoSetHessian()

424:    Output Parameters:
425: .  H    - Hessian matrix
426: .  Hpre - optionally different preconditioning matrix
427: .  flg  - flag indicating matrix structure

429: */
430: int FormHessian(TAO_APPLICATION taoapp,Vec X,Mat *H, Mat *Hpre, MatStructure *flg, void *ptr)
431: { 
432:   int    info;
433:   AppCtx *user = (AppCtx *) ptr;

435:   /* Evaluate the Hessian entries*/
436:   info = QuadraticH(user,X,*H); CHKERRQ(info);


439:   /* Indicate that this matrix has the same sparsity pattern during
440:      successive iterations; setting this flag can save significant work
441:      in computing the preconditioner for some methods. */
442:   *flg=SAME_NONZERO_PATTERN;

444:   return 0;
445: }

447: /* ------------------------------------------------------------------- */
450: /*
451:    QuadraticH - Evaluates Hessian matrix.

453:    Input Parameters:
454: .  user - user-defined context, as set by TaoSetHessian()
455: .  X    - input vector

457:    Output Parameter:
458: .  H    - Hessian matrix
459: */
460: int QuadraticH(AppCtx *user, Vec X, Mat Hessian)
461: {
462:   int    info;
463:   PetscInt i,j,k;
464:   PetscInt mx=user->mx, my=user->my;
465:   PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
466:   double hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
467:   double f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
468:   double hl,hr,ht,hb,hc,htl,hbr;
469:   PetscScalar **x, v[7];
470:   MatStencil col[7],row;
471:   Vec    localX;
472:   PetscTruth assembled;

474:   /* Get local mesh boundaries */
475:   info = DAGetLocalVector(user->da,&localX);CHKERRQ(info);

477:   info = DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);
478:   info = DAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); CHKERRQ(info);

480:   /* Scatter ghost points to local vector */
481:   info = DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX); CHKERRQ(info);
482:   info = DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX); CHKERRQ(info);

484:   /* Get pointers to vector data */
485:   info = DAVecGetArray(user->da,localX,(void**)&x);

487:   /* Initialize matrix entries to zero */
488:   info = MatAssembled(Hessian,&assembled); CHKERRQ(info);
489:   if (assembled){info = MatZeroEntries(Hessian);  CHKERRQ(info);}


492:   /* Set various matrix options */
493:   info = MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE); CHKERRQ(info);

495:   /* Compute Hessian over the locally owned part of the mesh */

497:   for (j=ys; j<ys+ym; j++){
498:       
499:     for (i=xs; i< xs+xm; i++){

501:       xc = x[j][i];
502:       xlt=xrb=xl=xr=xb=xt=xc;

504:       /* Left side */
505:       if (i==0){
506:         xl  = user->left[j-ys+1];
507:         xlt = user->left[j-ys+2];
508:       } else {
509:         xl  = x[j][i-1];
510:       }
511:       
512:       if (j==0){
513:         xb  = user->bottom[i-xs+1];
514:         xrb = user->bottom[i-xs+2];
515:       } else {
516:         xb  = x[j-1][i];
517:       }
518:       
519:       if (i+1 == mx){
520:         xr  = user->right[j-ys+1];
521:         xrb = user->right[j-ys];
522:       } else {
523:         xr  = x[j][i+1];
524:       }

526:       if (j+1==my){
527:         xt  = user->top[i-xs+1];
528:         xlt = user->top[i-xs];
529:       }else {
530:         xt  = x[j+1][i];
531:       }

533:       if (i>0 && j+1<my){
534:         xlt = x[j+1][i-1];
535:       }
536:       if (j>0 && i+1<mx){
537:         xrb = x[j-1][i+1];
538:       }


541:       d1 = (xc-xl)/hx;
542:       d2 = (xc-xr)/hx;
543:       d3 = (xc-xt)/hy;
544:       d4 = (xc-xb)/hy;
545:       d5 = (xrb-xr)/hy;
546:       d6 = (xrb-xb)/hx;
547:       d7 = (xlt-xl)/hy;
548:       d8 = (xlt-xt)/hx;
549:       
550:       f1 = sqrt( 1.0 + d1*d1 + d7*d7);
551:       f2 = sqrt( 1.0 + d1*d1 + d4*d4);
552:       f3 = sqrt( 1.0 + d3*d3 + d8*d8);
553:       f4 = sqrt( 1.0 + d3*d3 + d2*d2);
554:       f5 = sqrt( 1.0 + d2*d2 + d5*d5);
555:       f6 = sqrt( 1.0 + d4*d4 + d6*d6);


558:       hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
559:         (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
560:       hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
561:         (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
562:       ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
563:         (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
564:       hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
565:         (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);

567:       hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6);
568:       htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3);

570:       hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
571:         hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
572:         (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
573:         (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);

575:       hl/=2.0; hr/=2.0; ht/=2.0; hb/=2.0; hbr/=2.0; htl/=2.0;  hc/=2.0; 

577:       row.j = j; row.i = i;
578:       k=0;
579:       if (j>0){ 
580:         v[k]=hb;
581:         col[k].j = j - 1; col[k].i = i;
582:         k++;
583:       }
584:       
585:       if (j>0 && i < mx -1){
586:         v[k]=hbr;
587:         col[k].j = j - 1; col[k].i = i+1;
588:         k++;
589:       }
590:       
591:       if (i>0){
592:         v[k]= hl;
593:         col[k].j = j; col[k].i = i-1;
594:         k++;
595:       }
596:       
597:       v[k]= hc;
598:       col[k].j = j; col[k].i = i;
599:       k++;
600:       
601:       if (i < mx-1 ){
602:         v[k]= hr; 
603:         col[k].j = j; col[k].i = i+1;
604:         k++;
605:       }
606:       
607:       if (i>0 && j < my-1 ){
608:         v[k]= htl;
609:         col[k].j = j+1; col[k].i = i-1;
610:         k++;
611:       }
612:       
613:       if (j < my-1 ){
614:         v[k]= ht; 
615:         col[k].j = j+1; col[k].i = i;
616:         k++;
617:       }
618:       
619:       /* 
620:          Set matrix values using local numbering, which was defined
621:          earlier, in the main routine.
622:       */
623:       info = MatSetValuesStencil(Hessian,1,&row,k,col,v,INSERT_VALUES);
624:       CHKERRQ(info);
625:       
626:     }
627:   }
628:   
629:   /* Restore vectors */
630:   info = DAVecRestoreArray(user->da,localX,(void**)&x);

632:   info = DARestoreLocalVector(user->da,&localX); CHKERRQ(info);

634:   /* Assemble the matrix */
635:   info = MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
636:   info = MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY); CHKERRQ(info);

638:   info = PetscLogFlops(199*xm*ym); CHKERRQ(info);
639:   return 0;
640: }

642: /* ------------------------------------------------------------------- */
645: /* 
646:    MSA_BoundaryConditions -  Calculates the boundary conditions for
647:    the region.

649:    Input Parameter:
650: .  user - user-defined application context

652:    Output Parameter:
653: .  user - user-defined application context
654: */
655: static int MSA_BoundaryConditions(AppCtx * user)
656: {
657:   int        info;
658:   PetscInt     i,j,k,limit=0,maxits=5;
659:   PetscInt     xs,ys,xm,ym,gxs,gys,gxm,gym;
660:   PetscInt     mx=user->mx,my=user->my;
661:   PetscInt     bsize=0, lsize=0, tsize=0, rsize=0;
662:   double     one=1.0, two=2.0, three=3.0, tol=1e-10;
663:   double     fnorm,det,hx,hy,xt=0,yt=0;
664:   double     u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
665:   double     b=-0.5, t=0.5, l=-0.5, r=0.5;
666:   double     *boundary;
667:   PetscTruth   flg;

669:   /* Get local mesh boundaries */
670:   info = DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);
671:   info = DAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); CHKERRQ(info);

673:   bsize=xm+2;
674:   lsize=ym+2;
675:   rsize=ym+2;
676:   tsize=xm+2;

678:   info = PetscMalloc(bsize*sizeof(double),&user->bottom); CHKERRQ(info);
679:   info = PetscMalloc(tsize*sizeof(double),&user->top); CHKERRQ(info);
680:   info = PetscMalloc(lsize*sizeof(double),&user->left); CHKERRQ(info);
681:   info = PetscMalloc(rsize*sizeof(double),&user->right); CHKERRQ(info);

683:   hx= (r-l)/(mx+1); hy=(t-b)/(my+1);

685:   for (j=0; j<4; j++){
686:     if (j==0){
687:       yt=b;
688:       xt=l+hx*xs;
689:       limit=bsize;
690:       boundary=user->bottom;
691:     } else if (j==1){
692:       yt=t;
693:       xt=l+hx*xs;
694:       limit=tsize;
695:       boundary=user->top;
696:     } else if (j==2){
697:       yt=b+hy*ys;
698:       xt=l;
699:       limit=lsize;
700:       boundary=user->left;
701:     } else { //if (j==3)
702:       yt=b+hy*ys;
703:       xt=r;
704:       limit=rsize;
705:       boundary=user->right;
706:     }

708:     for (i=0; i<limit; i++){
709:       u1=xt;
710:       u2=-yt;
711:       for (k=0; k<maxits; k++){
712:         nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
713:         nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
714:         fnorm=sqrt(nf1*nf1+nf2*nf2);
715:         if (fnorm <= tol) break;
716:         njac11=one+u2*u2-u1*u1;
717:         njac12=two*u1*u2;
718:         njac21=-two*u1*u2;
719:         njac22=-one - u1*u1 + u2*u2;
720:         det = njac11*njac22-njac21*njac12;
721:         u1 = u1-(njac22*nf1-njac12*nf2)/det;
722:         u2 = u2-(njac11*nf2-njac21*nf1)/det;
723:       }

725:       boundary[i]=u1*u1-u2*u2;
726:       if (j==0 || j==1) {
727:         xt=xt+hx;
728:       } else { // if (j==2 || j==3)
729:         yt=yt+hy;
730:       }
731:       
732:     }

734:   }

736:   /* Scale the boundary if desired */
737:   if (1==1){ 
738:     PetscReal scl = 1.0;

740:     info = PetscOptionsGetReal(PETSC_NULL,"-bottom",&scl,&flg); 
741:     CHKERRQ(info);
742:     if (flg){
743:       for (i=0;i<bsize;i++) user->bottom[i]*=scl;
744:     }

746:     info = PetscOptionsGetReal(PETSC_NULL,"-top",&scl,&flg); 
747:     CHKERRQ(info);
748:     if (flg){
749:       for (i=0;i<tsize;i++) user->top[i]*=scl;
750:     }

752:     info = PetscOptionsGetReal(PETSC_NULL,"-right",&scl,&flg); 
753:     CHKERRQ(info);
754:     if (flg){
755:       for (i=0;i<rsize;i++) user->right[i]*=scl;
756:     }

758:     info = PetscOptionsGetReal(PETSC_NULL,"-left",&scl,&flg); 
759:     CHKERRQ(info);
760:     if (flg){
761:       for (i=0;i<lsize;i++) user->left[i]*=scl;
762:     }
763:   }
764:   
765:   return 0;
766: }

768: /* ------------------------------------------------------------------- */
771: /*
772:    MSA_InitialPoint - Calculates the initial guess in one of three ways. 

774:    Input Parameters:
775: .  user - user-defined application context
776: .  X - vector for initial guess

778:    Output Parameters:
779: .  X - newly computed initial guess
780: */
781: static int MSA_InitialPoint(AppCtx * user, Vec X)
782: {
783:   int      info;
784:   PetscInt   start2=-1,i,j;
785:   PetscReal   start1=0;
786:   PetscTruth flg1,flg2;

788:   info = PetscOptionsGetReal(PETSC_NULL,"-start",&start1,&flg1); CHKERRQ(info);
789:   info = PetscOptionsGetInt(PETSC_NULL,"-random",&start2,&flg2); CHKERRQ(info);

791:   if (flg1){ /* The zero vector is reasonable */
792:  
793:     info = VecSet(X, start1); CHKERRQ(info);

795:   } else if (flg2 && start2>0){ /* Try a random start between -0.5 and 0.5 */

797:     PetscRandom rctx;  PetscScalar np5=-0.5;

799:     info = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
800:     CHKERRQ(info);
801:     for (i=0; i<start2; i++){
802:       info = VecSetRandom(X, rctx); CHKERRQ(info);
803:     }
804:     info = PetscRandomDestroy(rctx); CHKERRQ(info);
805:     info = VecShift(X, np5); CHKERRQ(info);

807:   } else { /* Take an average of the boundary conditions */

809:     PetscInt xs,xm,ys,ym;
810:     PetscInt mx=user->mx,my=user->my;
811:     PetscScalar **x;
812:     
813:     /* Get local mesh boundaries */
814:     info = DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);
815:     
816:     /* Get pointers to vector data */
817:     info = DAVecGetArray(user->da,X,(void**)&x);

819:     /* Perform local computations */    
820:     for (j=ys; j<ys+ym; j++){
821:       for (i=xs; i< xs+xm; i++){
822:         x[j][i] = ( ((j+1)*user->bottom[i-xs+1]+(my-j+1)*user->top[i-xs+1])/(my+2)+
823:                    ((i+1)*user->left[j-ys+1]+(mx-i+1)*user->right[j-ys+1])/(mx+2))/2.0; 
824:       }
825:     }
826:     
827:     /* Restore vectors */
828:     info = DAVecRestoreArray(user->da,X,(void**)&x);  CHKERRQ(info);

830:     info = PetscLogFlops(9*xm*ym); CHKERRQ(info);
831:     
832:   }
833:   return 0;
834: }

836: /*-----------------------------------------------------------------------*/
839: int My_Monitor(TAO_APPLICATION minsurfapp, void *ctx){
840:   int info;
841:   Vec X;

843:   info = TaoAppGetSolutionVec(minsurfapp,&X); CHKERRQ(info);
844:   info = VecView(X,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(info);
845:   return 0;
846: }