Actual source code: minsurf2.c

tao-2.1-p0 2012-07-24
  1: /* Program usage: mpirun -np <proc> minsurf2 [-help] [all TAO options] */

  3: /*
  4:   Include "taosolver.h" so we can use TAO solvers.
  5:   petscdm.h for distributed array
  6: */
  7: #include "taosolver.h"
  8: #include "petscdm.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: TaoInitialize(); TaoFinalize();
 25:    Routines: TaoCreate(); TaoSetType();
 26:    Routines: TaoSetInitialVector(); 
 27:    Routines: TaoSetObjectiveAndGradientRoutine();
 28:    Routines: TaoSetHessianRoutine(); TaoSetFromOptions();
 29:    Routines: TaoSetMonitor();
 30:    Routines: TaoSolve(); TaoView();
 31:    Routines: TaoGetTerminationReason(); TaoDestroy();
 32:    Processors: n
 33: T*/

 35: /* 
 36:    User-defined application context - contains data needed by the 
 37:    application-provided call-back routines, FormFunctionGradient() 
 38:    and FormHessian().
 39: */
 40: typedef struct {
 41:   PetscInt      mx, my;                 /* discretization in x, y directions */
 42:   PetscReal      *bottom, *top, *left, *right;             /* boundary values */
 43:   DM          dm;                      /* distributed array data structure */
 44:   Mat         H;                       /* Hessian */
 45: } AppCtx;


 48: /* -------- User-defined Routines --------- */

 50: static PetscErrorCode MSA_BoundaryConditions(AppCtx*);
 51: static PetscErrorCode MSA_InitialPoint(AppCtx*,Vec);
 52: PetscErrorCode QuadraticH(AppCtx*,Vec,Mat);
 53: PetscErrorCode FormFunctionGradient(TaoSolver,Vec,PetscReal *,Vec,void*);
 54: PetscErrorCode FormGradient(TaoSolver,Vec,Vec,void*);
 55: PetscErrorCode FormHessian(TaoSolver,Vec,Mat*,Mat*,MatStructure *,void*);
 56: PetscErrorCode My_Monitor(TaoSolver, void *);

 60: int main( int argc, char **argv )
 61: {
 62:   PetscErrorCode    ierr;                /* used to check for functions returning nonzeros */
 63:   PetscInt          Nx, Ny;              /* number of processors in x- and y- directions */
 64:   Vec             x;                   /* solution, gradient vectors */
 65:   PetscBool      flg, viewmat;        /* flags */
 66:   PetscBool      fddefault, fdcoloring;   /* flags */
 67:   TaoSolverTerminationReason reason;           
 68:   TaoSolver       tao;                 /* TAO solver context */
 69:   AppCtx          user;                /* user-defined work context */
 70:   ISColoring     iscoloring;
 71:   MatFDColoring  matfdcoloring;

 73:   /* Initialize TAO */
 74:   PetscInitialize( &argc, &argv,(char *)0,help );
 75:   TaoInitialize( &argc, &argv,(char *)0,help );

 77:   /* Specify dimension of the problem */
 78:   user.mx = 10; user.my = 10;

 80:   /* Check for any command line arguments that override defaults */
 81:   PetscOptionsGetInt(PETSC_NULL,"-mx",&user.mx,&flg); 
 82:   PetscOptionsGetInt(PETSC_NULL,"-my",&user.my,&flg); 

 84:   PetscPrintf(MPI_COMM_WORLD,"\n---- Minimum Surface Area Problem -----\n");
 85:   PetscPrintf(MPI_COMM_WORLD,"mx: %D     my: %D   \n\n",user.mx,user.my);


 88:   /* Let PETSc determine the vector distribution */
 89:   Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;

 91:   /* Create distributed array (DM) to manage parallel grid and vectors  */
 92:   DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,
 93:                       DMDA_STENCIL_BOX,
 94:                       user.mx, user.my,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,
 95:                       &user.dm); 
 96:   

 98:   /* Create TAO solver and set desired solution method.*/
 99:   TaoCreate(PETSC_COMM_WORLD,&tao); 
100:   TaoSetType(tao,"tao_cg"); 

102:   /*
103:      Extract global vector from DA for the vector of variables --  PETSC routine
104:      Compute the initial solution                              --  application specific, see below
105:      Set this vector for use by TAO                            --  TAO routine
106:   */
107:   DMCreateGlobalVector(user.dm,&x); 
108:   MSA_BoundaryConditions(&user);          
109:   MSA_InitialPoint(&user,x); 
110:   TaoSetInitialVector(tao,x); 

112:   /* 
113:      Initialize the Application context for use in function evaluations  --  application specific, see below.
114:      Set routines for function and gradient evaluation 
115:   */
116:   TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&user); 

118:   /* 
119:      Given the command line arguments, calculate the hessian with either the user-
120:      provided function FormHessian, or the default finite-difference driven Hessian
121:      functions 
122:   */
123:   PetscOptionsHasName(PETSC_NULL,"-tao_fddefault",&fddefault);
124:   PetscOptionsHasName(PETSC_NULL,"-tao_fdcoloring",&fdcoloring);


127:   /* 
128:      Create a matrix data structure to store the Hessian and set 
129:      the Hessian evalution routine.
130:      Set the matrix structure to be used for Hessian evalutions
131:   */
132:   DMCreateMatrix(user.dm,MATAIJ,&user.H);
133:   MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE); 


136:   if (fdcoloring) {
137:       DMCreateColoring(user.dm,IS_COLORING_GLOBAL,MATAIJ,&iscoloring); 
138:       

140:       MatFDColoringCreate(user.H,iscoloring,&matfdcoloring);
141:       

143:       ISColoringDestroy(&iscoloring);
144:       MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormGradient,(void*)&user); 
145:       MatFDColoringSetFromOptions(matfdcoloring); 
146:       
147:       TaoSetHessianRoutine(tao,user.H,user.H,TaoDefaultComputeHessianColor,(void *)matfdcoloring); 

149:   } else if (fddefault){
150:       TaoSetHessianRoutine(tao,user.H,user.H,TaoDefaultComputeHessian,(void *)PETSC_NULL); 

152:   } else { 
153:       TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void *)&user); 
154:   }


157:   /* 
158:      If my_monitor option is in command line, then use the user-provided
159:      monitoring function
160:   */
161:   PetscOptionsHasName(PETSC_NULL,"-my_monitor",&viewmat); 
162:   if (viewmat){
163:     TaoSetMonitor(tao,My_Monitor,PETSC_NULL,PETSC_NULL); 
164:   }

166:   /* Check for any tao command line options */
167:   TaoSetFromOptions(tao); 

169:   /* SOLVE THE APPLICATION */
170:   TaoSolve(tao); 

172:   TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);

174:   /* Get information on termination */
175:   TaoGetTerminationReason(tao,&reason); 
176:   if (reason <= 0 ){
177:       PetscPrintf(MPI_COMM_WORLD,"Try a different TAO method \n");
178:   }


181:   /* Free TAO data structures */
182:   TaoDestroy(&tao); 

184:   /* Free PETSc data structures */
185:   VecDestroy(&x); 
186:   MatDestroy(&user.H); 
187:   if (fdcoloring) {
188:       MatFDColoringDestroy(&matfdcoloring); 
189:   }
190:   PetscFree(user.bottom); 
191:   PetscFree(user.top); 
192:   PetscFree(user.left); 
193:   PetscFree(user.right); 
194:   DMDestroy(&user.dm); 

196:   /* Finalize TAO */
197:   TaoFinalize();
198:   PetscFinalize();
199:   
200:   return 0;
201: }

205: PetscErrorCode FormGradient(TaoSolver tao, Vec X, Vec G,void *userCtx){
207:   PetscReal fcn;
209:   FormFunctionGradient(tao,X,&fcn,G,userCtx);
210:   return(0);
211: }

213: /* -------------------------------------------------------------------- */
216: /*  FormFunctionGradient - Evaluates the function and corresponding gradient.

218:     Input Parameters:
219: .   tao     - the TaoSolver context
220: .   XX      - input vector
221: .   userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()
222:     
223:     Output Parameters:
224: .   fcn     - the newly evaluated function
225: .   GG       - vector containing the newly evaluated gradient
226: */
227: PetscErrorCode FormFunctionGradient(TaoSolver tao, Vec X, PetscReal *fcn,Vec G,void *userCtx){

229:   AppCtx * user = (AppCtx *) userCtx;
230:   PetscErrorCode    ierr;
231:   PetscInt i,j;
232:   PetscInt mx=user->mx, my=user->my;
233:   PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
234:   PetscReal ft=0.0;
235:   PetscReal hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy;
236:   PetscReal rhx=mx+1, rhy=my+1;
237:   PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
238:   PetscReal df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
239:   PetscReal **g, **x;
240:   Vec    localX;

243:   /* Get local mesh boundaries */
244:   DMGetLocalVector(user->dm,&localX);

246:   DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); 
247:   DMDAGetGhostCorners(user->dm,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); 

249:   /* Scatter ghost points to local vector */
250:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX); 
251:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX); 

253:   /* Get pointers to vector data */
254:   DMDAVecGetArray(user->dm,localX,(void**)&x);
255:   DMDAVecGetArray(user->dm,G,(void**)&g);

257:   /* Compute function and gradient over the locally owned part of the mesh */
258:   for (j=ys; j<ys+ym; j++){
259:     for (i=xs; i< xs+xm; i++){
260:       
261:       xc = x[j][i];
262:       xlt=xrb=xl=xr=xb=xt=xc;
263:       
264:       if (i==0){ /* left side */
265:         xl= user->left[j-ys+1];
266:         xlt = user->left[j-ys+2];
267:       } else {
268:         xl = x[j][i-1];
269:       }

271:       if (j==0){ /* bottom side */
272:         xb=user->bottom[i-xs+1];
273:         xrb =user->bottom[i-xs+2];
274:       } else {
275:         xb = x[j-1][i];
276:       }
277:       
278:       if (i+1 == gxs+gxm){ /* right side */
279:         xr=user->right[j-ys+1];
280:         xrb = user->right[j-ys];
281:       } else {
282:         xr = x[j][i+1];
283:       }

285:       if (j+1==gys+gym){ /* top side */
286:         xt=user->top[i-xs+1];
287:         xlt = user->top[i-xs];
288:       }else {
289:         xt = x[j+1][i];
290:       }

292:       if (i>gxs && j+1<gys+gym){
293:         xlt = x[j+1][i-1];
294:       }
295:       if (j>gys && i+1<gxs+gxm){
296:         xrb = x[j-1][i+1];
297:       }

299:       d1 = (xc-xl);
300:       d2 = (xc-xr);
301:       d3 = (xc-xt);
302:       d4 = (xc-xb);
303:       d5 = (xr-xrb);
304:       d6 = (xrb-xb);
305:       d7 = (xlt-xl);
306:       d8 = (xt-xlt);
307:       
308:       df1dxc = d1*hydhx;
309:       df2dxc = ( d1*hydhx + d4*hxdhy );
310:       df3dxc = d3*hxdhy;
311:       df4dxc = ( d2*hydhx + d3*hxdhy );
312:       df5dxc = d2*hydhx;
313:       df6dxc = d4*hxdhy;

315:       d1 *= rhx;
316:       d2 *= rhx;
317:       d3 *= rhy;
318:       d4 *= rhy;
319:       d5 *= rhy;
320:       d6 *= rhx;
321:       d7 *= rhy;
322:       d8 *= rhx;

324:       f1 = PetscSqrtReal( 1.0 + d1*d1 + d7*d7);
325:       f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
326:       f3 = PetscSqrtReal( 1.0 + d3*d3 + d8*d8);
327:       f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
328:       f5 = PetscSqrtReal( 1.0 + d2*d2 + d5*d5);
329:       f6 = PetscSqrtReal( 1.0 + d4*d4 + d6*d6);
330:       
331:       f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
332:       f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);

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

336:       df1dxc /= f1;
337:       df2dxc /= f2;
338:       df3dxc /= f3;
339:       df4dxc /= f4;
340:       df5dxc /= f5;
341:       df6dxc /= f6;

343:       g[j][i] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc ) * 0.5;
344:       
345:     }
346:   }

348:   /* Compute triangular areas along the border of the domain. */
349:   if (xs==0){ /* left side */
350:     for (j=ys; j<ys+ym; j++){
351:       d3=(user->left[j-ys+1] - user->left[j-ys+2])*rhy;
352:       d2=(user->left[j-ys+1] - x[j][0]) *rhx;
353:       ft = ft+PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
354:     }
355:   }
356:   if (ys==0){ /* bottom side */
357:     for (i=xs; i<xs+xm; i++){
358:       d2=(user->bottom[i+1-xs]-user->bottom[i-xs+2])*rhx;
359:       d3=(user->bottom[i-xs+1]-x[0][i])*rhy;
360:       ft = ft+PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
361:     }
362:   }

364:   if (xs+xm==mx){ /* right side */
365:     for (j=ys; j< ys+ym; j++){
366:       d1=(x[j][mx-1] - user->right[j-ys+1])*rhx;
367:       d4=(user->right[j-ys]-user->right[j-ys+1])*rhy;
368:       ft = ft+PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
369:     }
370:   }
371:   if (ys+ym==my){ /* top side */
372:     for (i=xs; i<xs+xm; i++){
373:       d1=(x[my-1][i] - user->top[i-xs+1])*rhy;
374:       d4=(user->top[i-xs+1] - user->top[i-xs])*rhx;
375:       ft = ft+PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
376:     }
377:   }

379:   if (ys==0 && xs==0){
380:     d1=(user->left[0]-user->left[1])/hy;
381:     d2=(user->bottom[0]-user->bottom[1])*rhx;
382:     ft +=PetscSqrtReal( 1.0 + d1*d1 + d2*d2);
383:   }
384:   if (ys+ym == my && xs+xm == mx){
385:     d1=(user->right[ym+1] - user->right[ym])*rhy;
386:     d2=(user->top[xm+1] - user->top[xm])*rhx;
387:     ft +=PetscSqrtReal( 1.0 + d1*d1 + d2*d2);
388:   }

390:   ft=ft*area;
391:   MPI_Allreduce(&ft,fcn,1,MPIU_REAL,MPIU_SUM,MPI_COMM_WORLD);

393:   /* Restore vectors */
394:   DMDAVecRestoreArray(user->dm,localX,(void**)&x);
395:   DMDAVecRestoreArray(user->dm,G,(void**)&g);

397:   /* Scatter values to global vector */
398:   DMRestoreLocalVector(user->dm,&localX); 

400:   PetscLogFlops(67*xm*ym); 

402:   return(0);
403: }

405: /* ------------------------------------------------------------------- */
408: /*
409:    FormHessian - Evaluates Hessian matrix.

411:    Input Parameters:
412: .  tao  - the TaoSolver context
413: .  x    - input vector
414: .  ptr  - optional user-defined context, as set by TaoSetHessianRoutine()

416:    Output Parameters:
417: .  H    - Hessian matrix
418: .  Hpre - optionally different preconditioning matrix
419: .  flg  - flag indicating matrix structure

421: */
422: PetscErrorCode FormHessian(TaoSolver tao,Vec X,Mat *H, Mat *Hpre, MatStructure *flg, void *ptr)
423: { 
424:   PetscErrorCode    ierr;
425:   AppCtx *user = (AppCtx *) ptr;

428:   /* Evaluate the Hessian entries*/
429:   QuadraticH(user,X,*H); 


432:   /* Indicate that this matrix has the same sparsity pattern during
433:      successive iterations; setting this flag can save significant work
434:      in computing the preconditioner for some methods. */
435:   *flg=SAME_NONZERO_PATTERN;
436:   return(0);
437: }

439: /* ------------------------------------------------------------------- */
442: /*
443:    QuadraticH - Evaluates Hessian matrix.

445:    Input Parameters:
446: .  user - user-defined context, as set by TaoSetHessian()
447: .  X    - input vector

449:    Output Parameter:
450: .  H    - Hessian matrix
451: */
452: PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian)
453: {
454:   PetscErrorCode    ierr;
455:   PetscInt i,j,k;
456:   PetscInt mx=user->mx, my=user->my;
457:   PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
458:   PetscReal hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
459:   PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
460:   PetscReal hl,hr,ht,hb,hc,htl,hbr;
461:   PetscReal **x, v[7];
462:   MatStencil col[7],row;
463:   Vec    localX;
464:   PetscBool assembled;

467:   /* Get local mesh boundaries */
468:   DMGetLocalVector(user->dm,&localX);

470:   DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); 
471:   DMDAGetGhostCorners(user->dm,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); 

473:   /* Scatter ghost points to local vector */
474:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX); 
475:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX); 

477:   /* Get pointers to vector data */
478:   DMDAVecGetArray(user->dm,localX,(void**)&x);

480:   /* Initialize matrix entries to zero */
481:   MatAssembled(Hessian,&assembled); 
482:   if (assembled){MatZeroEntries(Hessian);  }


485:   /* Set various matrix options */
486:   MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE); 

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

490:   for (j=ys; j<ys+ym; j++){
491:       
492:     for (i=xs; i< xs+xm; i++){

494:       xc = x[j][i];
495:       xlt=xrb=xl=xr=xb=xt=xc;

497:       /* Left side */
498:       if (i==0){
499:         xl  = user->left[j-ys+1];
500:         xlt = user->left[j-ys+2];
501:       } else {
502:         xl  = x[j][i-1];
503:       }
504:       
505:       if (j==0){
506:         xb  = user->bottom[i-xs+1];
507:         xrb = user->bottom[i-xs+2];
508:       } else {
509:         xb  = x[j-1][i];
510:       }
511:       
512:       if (i+1 == mx){
513:         xr  = user->right[j-ys+1];
514:         xrb = user->right[j-ys];
515:       } else {
516:         xr  = x[j][i+1];
517:       }

519:       if (j+1==my){
520:         xt  = user->top[i-xs+1];
521:         xlt = user->top[i-xs];
522:       }else {
523:         xt  = x[j+1][i];
524:       }

526:       if (i>0 && j+1<my){
527:         xlt = x[j+1][i-1];
528:       }
529:       if (j>0 && i+1<mx){
530:         xrb = x[j-1][i+1];
531:       }


534:       d1 = (xc-xl)/hx;
535:       d2 = (xc-xr)/hx;
536:       d3 = (xc-xt)/hy;
537:       d4 = (xc-xb)/hy;
538:       d5 = (xrb-xr)/hy;
539:       d6 = (xrb-xb)/hx;
540:       d7 = (xlt-xl)/hy;
541:       d8 = (xlt-xt)/hx;
542:       
543:       f1 = PetscSqrtReal( 1.0 + d1*d1 + d7*d7);
544:       f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
545:       f3 = PetscSqrtReal( 1.0 + d3*d3 + d8*d8);
546:       f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
547:       f5 = PetscSqrtReal( 1.0 + d2*d2 + d5*d5);
548:       f6 = PetscSqrtReal( 1.0 + d4*d4 + d6*d6);


551:       hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
552:         (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
553:       hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
554:         (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
555:       ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
556:         (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
557:       hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
558:         (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);

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

563:       hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
564:         hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
565:         (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
566:         (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);

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

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

625:   DMRestoreLocalVector(user->dm,&localX); 

627:   /* Assemble the matrix */
628:   MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY); 
629:   MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY); 

631:   PetscLogFlops(199*xm*ym); 
632:   return(0);
633: }

635: /* ------------------------------------------------------------------- */
638: /* 
639:    MSA_BoundaryConditions -  Calculates the boundary conditions for
640:    the region.

642:    Input Parameter:
643: .  user - user-defined application context

645:    Output Parameter:
646: .  user - user-defined application context
647: */
648: static PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
649: {
650:   PetscErrorCode        ierr;
651:   PetscInt     i,j,k,limit=0,maxits=5;
652:   PetscInt     xs,ys,xm,ym,gxs,gys,gxm,gym;
653:   PetscInt     mx=user->mx,my=user->my;
654:   PetscInt     bsize=0, lsize=0, tsize=0, rsize=0;
655:   PetscReal     one=1.0, two=2.0, three=3.0, tol=1e-10;
656:   PetscReal     fnorm,det,hx,hy,xt=0,yt=0;
657:   PetscReal     u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
658:   PetscReal     b=-0.5, t=0.5, l=-0.5, r=0.5;
659:   PetscReal     *boundary;
660:   PetscBool   flg;

663:   /* Get local mesh boundaries */
664:   DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); 
665:   DMDAGetGhostCorners(user->dm,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); 

667:   bsize=xm+2;
668:   lsize=ym+2;
669:   rsize=ym+2;
670:   tsize=xm+2;

672:   PetscMalloc(bsize*sizeof(PetscReal),&user->bottom); 
673:   PetscMalloc(tsize*sizeof(PetscReal),&user->top); 
674:   PetscMalloc(lsize*sizeof(PetscReal),&user->left); 
675:   PetscMalloc(rsize*sizeof(PetscReal),&user->right); 

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

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

702:     for (i=0; i<limit; i++){
703:       u1=xt;
704:       u2=-yt;
705:       for (k=0; k<maxits; k++){
706:         nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
707:         nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
708:         fnorm=PetscSqrtReal(nf1*nf1+nf2*nf2);
709:         if (fnorm <= tol) break;
710:         njac11=one+u2*u2-u1*u1;
711:         njac12=two*u1*u2;
712:         njac21=-two*u1*u2;
713:         njac22=-one - u1*u1 + u2*u2;
714:         det = njac11*njac22-njac21*njac12;
715:         u1 = u1-(njac22*nf1-njac12*nf2)/det;
716:         u2 = u2-(njac11*nf2-njac21*nf1)/det;
717:       }

719:       boundary[i]=u1*u1-u2*u2;
720:       if (j==0 || j==1) {
721:         xt=xt+hx;
722:       } else { /*  if (j==2 || j==3) */
723:         yt=yt+hy;
724:       }
725:       
726:     }

728:   }

730:   /* Scale the boundary if desired */
731:   if (1==1){ 
732:     PetscReal scl = 1.0;

734:     PetscOptionsGetReal(PETSC_NULL,"-bottom",&scl,&flg); 
735:     
736:     if (flg){
737:       for (i=0;i<bsize;i++) user->bottom[i]*=scl;
738:     }

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

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

752:     PetscOptionsGetReal(PETSC_NULL,"-left",&scl,&flg); 
753:     
754:     if (flg){
755:       for (i=0;i<lsize;i++) user->left[i]*=scl;
756:     }
757:   }
758:   
759:   return(0);
760: }

762: /* ------------------------------------------------------------------- */
765: /*
766:    MSA_InitialPoint - Calculates the initial guess in one of three ways. 

768:    Input Parameters:
769: .  user - user-defined application context
770: .  X - vector for initial guess

772:    Output Parameters:
773: .  X - newly computed initial guess
774: */
775: static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
776: {
777:   PetscErrorCode      ierr;
778:   PetscInt   start2=-1,i,j;
779:   PetscReal   start1=0;
780:   PetscBool flg1,flg2;

783:   PetscOptionsGetReal(PETSC_NULL,"-start",&start1,&flg1); 
784:   PetscOptionsGetInt(PETSC_NULL,"-random",&start2,&flg2); 

786:   if (flg1){ /* The zero vector is reasonable */
787:  
788:     VecSet(X, start1); 

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

792:     PetscRandom rctx;  PetscReal np5=-0.5;

794:     PetscRandomCreate(PETSC_COMM_WORLD,&rctx); 
795:     PetscRandomSetType(rctx,PETSCRAND); 
796:     
797:     for (i=0; i<start2; i++){
798:       VecSetRandom(X, rctx); 
799:     }
800:     PetscRandomDestroy(&rctx); 
801:     VecShift(X, np5); 

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

805:     PetscInt xs,xm,ys,ym;
806:     PetscInt mx=user->mx,my=user->my;
807:     PetscReal **x;
808:     
809:     /* Get local mesh boundaries */
810:     DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); 
811:     
812:     /* Get pointers to vector data */
813:     DMDAVecGetArray(user->dm,X,(void**)&x);

815:     /* Perform local computations */    
816:     for (j=ys; j<ys+ym; j++){
817:       for (i=xs; i< xs+xm; i++){
818:         x[j][i] = ( ((j+1)*user->bottom[i-xs+1]+(my-j+1)*user->top[i-xs+1])/(my+2)+
819:                    ((i+1)*user->left[j-ys+1]+(mx-i+1)*user->right[j-ys+1])/(mx+2))/2.0; 
820:       }
821:     }
822:     
823:     /* Restore vectors */
824:     DMDAVecRestoreArray(user->dm,X,(void**)&x);  

826:     PetscLogFlops(9*xm*ym); 
827:     
828:   }
829:   return(0);
830: }

832: /*-----------------------------------------------------------------------*/
835: PetscErrorCode My_Monitor(TaoSolver tao, void *ctx){
837:   Vec X;
839:   TaoGetSolutionVector(tao,&X); 
840:   VecView(X,PETSC_VIEWER_STDOUT_WORLD); 
841:   return(0);
842: }