Actual source code: minsurf2.c

petsc-dev 2014-08-28
Report Typos and Errors
  1: /* Program usage: mpirun -np <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: TaoGetConvergedReason(); 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;


 47: /* -------- User-defined Routines --------- */

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

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

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

 75:   /* Specify dimension of the problem */
 76:   user.mx = 10; user.my = 10;

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

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


 86:   /* Let PETSc determine the vector distribution */
 87:   Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;

 89:   /* Create distributed array (DM) to manage parallel grid and vectors  */
 90:   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);


 93:   /* Create TAO solver and set desired solution method.*/
 94:   TaoCreate(PETSC_COMM_WORLD,&tao);
 95:   TaoSetType(tao,TAOCG);

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

107:   /*
108:      Initialize the Application context for use in function evaluations  --  application specific, see below.
109:      Set routines for function and gradient evaluation
110:   */
111:   TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&user);

113:   /*
114:      Given the command line arguments, calculate the hessian with either the user-
115:      provided function FormHessian, or the default finite-difference driven Hessian
116:      functions
117:   */
118:   PetscOptionsHasName(NULL,"-tao_fddefault",&fddefault);
119:   PetscOptionsHasName(NULL,"-tao_fdcoloring",&fdcoloring);


122:   /*
123:      Create a matrix data structure to store the Hessian and set
124:      the Hessian evalution routine.
125:      Set the matrix structure to be used for Hessian evalutions
126:   */
127:   DMCreateMatrix(user.dm,&user.H);
128:   MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE);


131:   if (fdcoloring) {
132:     DMCreateColoring(user.dm,IS_COLORING_GLOBAL,&iscoloring);
133: 

135:       MatFDColoringCreate(user.H,iscoloring,&matfdcoloring);
136: 

138:       ISColoringDestroy(&iscoloring);
139:       MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormGradient,(void*)&user);
140:       MatFDColoringSetFromOptions(matfdcoloring);

142:       TaoSetHessianRoutine(tao,user.H,user.H,TaoDefaultComputeHessianColor,(void *)matfdcoloring);

144:   } else if (fddefault){
145:       TaoSetHessianRoutine(tao,user.H,user.H,TaoDefaultComputeHessian,(void *)NULL);

147:   } else {
148:       TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void *)&user);
149:   }


152:   /*
153:      If my_monitor option is in command line, then use the user-provided
154:      monitoring function
155:   */
156:   PetscOptionsHasName(NULL,"-my_monitor",&viewmat);
157:   if (viewmat){
158:     TaoSetMonitor(tao,My_Monitor,NULL,NULL);
159:   }

161:   /* Check for any tao command line options */
162:   TaoSetFromOptions(tao);

164:   /* SOLVE THE APPLICATION */
165:   TaoSolve(tao);

167:   TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);

169:   /* Get information on termination */
170:   TaoGetConvergedReason(tao,&reason);
171:   if (reason <= 0 ){
172:       PetscPrintf(MPI_COMM_WORLD,"Try a different TAO method \n");
173:   }


176:   /* Free TAO data structures */
177:   TaoDestroy(&tao);

179:   /* Free PETSc data structures */
180:   VecDestroy(&x);
181:   MatDestroy(&user.H);
182:   if (fdcoloring) {
183:       MatFDColoringDestroy(&matfdcoloring);
184:   }
185:   PetscFree(user.bottom);
186:   PetscFree(user.top);
187:   PetscFree(user.left);
188:   PetscFree(user.right);
189:   DMDestroy(&user.dm);

191:   PetscFinalize();
192:   return 0;
193: }

197: PetscErrorCode FormGradient(Tao tao, Vec X, Vec G,void *userCtx){
199:   PetscReal fcn;
201:   FormFunctionGradient(tao,X,&fcn,G,userCtx);
202:   return(0);
203: }

205: /* -------------------------------------------------------------------- */
208: /*  FormFunctionGradient - Evaluates the function and corresponding gradient.

210:     Input Parameters:
211: .   tao     - the Tao context
212: .   XX      - input vector
213: .   userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()

215:     Output Parameters:
216: .   fcn     - the newly evaluated function
217: .   GG       - vector containing the newly evaluated gradient
218: */
219: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn,Vec G,void *userCtx){

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

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

238:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
239:   DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

241:   /* Scatter ghost points to local vector */
242:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
243:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);

245:   /* Get pointers to vector data */
246:   DMDAVecGetArray(user->dm,localX,(void**)&x);
247:   DMDAVecGetArray(user->dm,G,(void**)&g);

249:   /* Compute function and gradient over the locally owned part of the mesh */
250:   for (j=ys; j<ys+ym; j++){
251:     for (i=xs; i< xs+xm; i++){

253:       xc = x[j][i];
254:       xlt=xrb=xl=xr=xb=xt=xc;

256:       if (i==0){ /* left side */
257:         xl= user->left[j-ys+1];
258:         xlt = user->left[j-ys+2];
259:       } else {
260:         xl = x[j][i-1];
261:       }

263:       if (j==0){ /* bottom side */
264:         xb=user->bottom[i-xs+1];
265:         xrb =user->bottom[i-xs+2];
266:       } else {
267:         xb = x[j-1][i];
268:       }

270:       if (i+1 == gxs+gxm){ /* right side */
271:         xr=user->right[j-ys+1];
272:         xrb = user->right[j-ys];
273:       } else {
274:         xr = x[j][i+1];
275:       }

277:       if (j+1==gys+gym){ /* top side */
278:         xt=user->top[i-xs+1];
279:         xlt = user->top[i-xs];
280:       }else {
281:         xt = x[j+1][i];
282:       }

284:       if (i>gxs && j+1<gys+gym){
285:         xlt = x[j+1][i-1];
286:       }
287:       if (j>gys && i+1<gxs+gxm){
288:         xrb = x[j-1][i+1];
289:       }

291:       d1 = (xc-xl);
292:       d2 = (xc-xr);
293:       d3 = (xc-xt);
294:       d4 = (xc-xb);
295:       d5 = (xr-xrb);
296:       d6 = (xrb-xb);
297:       d7 = (xlt-xl);
298:       d8 = (xt-xlt);

300:       df1dxc = d1*hydhx;
301:       df2dxc = ( d1*hydhx + d4*hxdhy );
302:       df3dxc = d3*hxdhy;
303:       df4dxc = ( d2*hydhx + d3*hxdhy );
304:       df5dxc = d2*hydhx;
305:       df6dxc = d4*hxdhy;

307:       d1 *= rhx;
308:       d2 *= rhx;
309:       d3 *= rhy;
310:       d4 *= rhy;
311:       d5 *= rhy;
312:       d6 *= rhx;
313:       d7 *= rhy;
314:       d8 *= rhx;

316:       f1 = PetscSqrtReal( 1.0 + d1*d1 + d7*d7);
317:       f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
318:       f3 = PetscSqrtReal( 1.0 + d3*d3 + d8*d8);
319:       f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
320:       f5 = PetscSqrtReal( 1.0 + d2*d2 + d5*d5);
321:       f6 = PetscSqrtReal( 1.0 + d4*d4 + d6*d6);

323:       f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
324:       f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);

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

328:       df1dxc /= f1;
329:       df2dxc /= f2;
330:       df3dxc /= f3;
331:       df4dxc /= f4;
332:       df5dxc /= f5;
333:       df6dxc /= f6;

335:       g[j][i] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc ) * 0.5;

337:     }
338:   }

340:   /* Compute triangular areas along the border of the domain. */
341:   if (xs==0){ /* left side */
342:     for (j=ys; j<ys+ym; j++){
343:       d3=(user->left[j-ys+1] - user->left[j-ys+2])*rhy;
344:       d2=(user->left[j-ys+1] - x[j][0]) *rhx;
345:       ft = ft+PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
346:     }
347:   }
348:   if (ys==0){ /* bottom side */
349:     for (i=xs; i<xs+xm; i++){
350:       d2=(user->bottom[i+1-xs]-user->bottom[i-xs+2])*rhx;
351:       d3=(user->bottom[i-xs+1]-x[0][i])*rhy;
352:       ft = ft+PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
353:     }
354:   }

356:   if (xs+xm==mx){ /* right side */
357:     for (j=ys; j< ys+ym; j++){
358:       d1=(x[j][mx-1] - user->right[j-ys+1])*rhx;
359:       d4=(user->right[j-ys]-user->right[j-ys+1])*rhy;
360:       ft = ft+PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
361:     }
362:   }
363:   if (ys+ym==my){ /* top side */
364:     for (i=xs; i<xs+xm; i++){
365:       d1=(x[my-1][i] - user->top[i-xs+1])*rhy;
366:       d4=(user->top[i-xs+1] - user->top[i-xs])*rhx;
367:       ft = ft+PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
368:     }
369:   }

371:   if (ys==0 && xs==0){
372:     d1=(user->left[0]-user->left[1])/hy;
373:     d2=(user->bottom[0]-user->bottom[1])*rhx;
374:     ft +=PetscSqrtReal( 1.0 + d1*d1 + d2*d2);
375:   }
376:   if (ys+ym == my && xs+xm == mx){
377:     d1=(user->right[ym+1] - user->right[ym])*rhy;
378:     d2=(user->top[xm+1] - user->top[xm])*rhx;
379:     ft +=PetscSqrtReal( 1.0 + d1*d1 + d2*d2);
380:   }

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

385:   /* Restore vectors */
386:   DMDAVecRestoreArray(user->dm,localX,(void**)&x);
387:   DMDAVecRestoreArray(user->dm,G,(void**)&g);

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

392:   PetscLogFlops(67*xm*ym);

394:   return(0);
395: }

397: /* ------------------------------------------------------------------- */
400: /*
401:    FormHessian - Evaluates Hessian matrix.

403:    Input Parameters:
404: .  tao  - the Tao context
405: .  x    - input vector
406: .  ptr  - optional user-defined context, as set by TaoSetHessianRoutine()

408:    Output Parameters:
409: .  H    - Hessian matrix
410: .  Hpre - optionally different preconditioning matrix
411: .  flg  - flag indicating matrix structure

413: */
414: PetscErrorCode FormHessian(Tao tao,Vec X,Mat H, Mat Hpre, void *ptr)
415: {
417:   AppCtx         *user = (AppCtx *) ptr;

420:   /* Evaluate the Hessian entries*/
421:   QuadraticH(user,X,H);
422:   return(0);
423: }

425: /* ------------------------------------------------------------------- */
428: /*
429:    QuadraticH - Evaluates Hessian matrix.

431:    Input Parameters:
432: .  user - user-defined context, as set by TaoSetHessian()
433: .  X    - input vector

435:    Output Parameter:
436: .  H    - Hessian matrix
437: */
438: PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian)
439: {
440:   PetscErrorCode    ierr;
441:   PetscInt i,j,k;
442:   PetscInt mx=user->mx, my=user->my;
443:   PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
444:   PetscReal hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
445:   PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
446:   PetscReal hl,hr,ht,hb,hc,htl,hbr;
447:   PetscReal **x, v[7];
448:   MatStencil col[7],row;
449:   Vec    localX;
450:   PetscBool assembled;

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

456:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
457:   DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

459:   /* Scatter ghost points to local vector */
460:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
461:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);

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

466:   /* Initialize matrix entries to zero */
467:   MatAssembled(Hessian,&assembled);
468:   if (assembled){MatZeroEntries(Hessian); }


471:   /* Set various matrix options */
472:   MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);

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

476:   for (j=ys; j<ys+ym; j++){

478:     for (i=xs; i< xs+xm; i++){

480:       xc = x[j][i];
481:       xlt=xrb=xl=xr=xb=xt=xc;

483:       /* Left side */
484:       if (i==0){
485:         xl  = user->left[j-ys+1];
486:         xlt = user->left[j-ys+2];
487:       } else {
488:         xl  = x[j][i-1];
489:       }

491:       if (j==0){
492:         xb  = user->bottom[i-xs+1];
493:         xrb = user->bottom[i-xs+2];
494:       } else {
495:         xb  = x[j-1][i];
496:       }

498:       if (i+1 == mx){
499:         xr  = user->right[j-ys+1];
500:         xrb = user->right[j-ys];
501:       } else {
502:         xr  = x[j][i+1];
503:       }

505:       if (j+1==my){
506:         xt  = user->top[i-xs+1];
507:         xlt = user->top[i-xs];
508:       }else {
509:         xt  = x[j+1][i];
510:       }

512:       if (i>0 && j+1<my){
513:         xlt = x[j+1][i-1];
514:       }
515:       if (j>0 && i+1<mx){
516:         xrb = x[j-1][i+1];
517:       }


520:       d1 = (xc-xl)/hx;
521:       d2 = (xc-xr)/hx;
522:       d3 = (xc-xt)/hy;
523:       d4 = (xc-xb)/hy;
524:       d5 = (xrb-xr)/hy;
525:       d6 = (xrb-xb)/hx;
526:       d7 = (xlt-xl)/hy;
527:       d8 = (xlt-xt)/hx;

529:       f1 = PetscSqrtReal( 1.0 + d1*d1 + d7*d7);
530:       f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
531:       f3 = PetscSqrtReal( 1.0 + d3*d3 + d8*d8);
532:       f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
533:       f5 = PetscSqrtReal( 1.0 + d2*d2 + d5*d5);
534:       f6 = PetscSqrtReal( 1.0 + d4*d4 + d6*d6);


537:       hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
538:         (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
539:       hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
540:         (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
541:       ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
542:         (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
543:       hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
544:         (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);

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

549:       hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
550:         hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
551:         (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
552:         (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);

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

556:       row.j = j; row.i = i;
557:       k=0;
558:       if (j>0){
559:         v[k]=hb;
560:         col[k].j = j - 1; col[k].i = i;
561:         k++;
562:       }

564:       if (j>0 && i < mx -1){
565:         v[k]=hbr;
566:         col[k].j = j - 1; col[k].i = i+1;
567:         k++;
568:       }

570:       if (i>0){
571:         v[k]= hl;
572:         col[k].j = j; col[k].i = i-1;
573:         k++;
574:       }

576:       v[k]= hc;
577:       col[k].j = j; col[k].i = i;
578:       k++;

580:       if (i < mx-1 ){
581:         v[k]= hr;
582:         col[k].j = j; col[k].i = i+1;
583:         k++;
584:       }

586:       if (i>0 && j < my-1 ){
587:         v[k]= htl;
588:         col[k].j = j+1; col[k].i = i-1;
589:         k++;
590:       }

592:       if (j < my-1 ){
593:         v[k]= ht;
594:         col[k].j = j+1; col[k].i = i;
595:         k++;
596:       }

598:       /*
599:          Set matrix values using local numbering, which was defined
600:          earlier, in the main routine.
601:       */
602:       MatSetValuesStencil(Hessian,1,&row,k,col,v,INSERT_VALUES);
603: 

605:     }
606:   }

608:   /* Restore vectors */
609:   DMDAVecRestoreArray(user->dm,localX,(void**)&x);

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

613:   /* Assemble the matrix */
614:   MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY);
615:   MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY);

617:   PetscLogFlops(199*xm*ym);
618:   return(0);
619: }

621: /* ------------------------------------------------------------------- */
624: /*
625:    MSA_BoundaryConditions -  Calculates the boundary conditions for
626:    the region.

628:    Input Parameter:
629: .  user - user-defined application context

631:    Output Parameter:
632: .  user - user-defined application context
633: */
634: static PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
635: {
636:   PetscErrorCode        ierr;
637:   PetscInt     i,j,k,limit=0,maxits=5;
638:   PetscInt     xs,ys,xm,ym,gxs,gys,gxm,gym;
639:   PetscInt     mx=user->mx,my=user->my;
640:   PetscInt     bsize=0, lsize=0, tsize=0, rsize=0;
641:   PetscReal     one=1.0, two=2.0, three=3.0, tol=1e-10;
642:   PetscReal     fnorm,det,hx,hy,xt=0,yt=0;
643:   PetscReal     u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
644:   PetscReal     b=-0.5, t=0.5, l=-0.5, r=0.5;
645:   PetscReal     *boundary;
646:   PetscBool   flg;

649:   /* Get local mesh boundaries */
650:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
651:   DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

653:   bsize=xm+2;
654:   lsize=ym+2;
655:   rsize=ym+2;
656:   tsize=xm+2;

658:   PetscMalloc1(bsize,&user->bottom);
659:   PetscMalloc1(tsize,&user->top);
660:   PetscMalloc1(lsize,&user->left);
661:   PetscMalloc1(rsize,&user->right);

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

665:   for (j=0; j<4; j++){
666:     if (j==0){
667:       yt=b;
668:       xt=l+hx*xs;
669:       limit=bsize;
670:       boundary=user->bottom;
671:     } else if (j==1){
672:       yt=t;
673:       xt=l+hx*xs;
674:       limit=tsize;
675:       boundary=user->top;
676:     } else if (j==2){
677:       yt=b+hy*ys;
678:       xt=l;
679:       limit=lsize;
680:       boundary=user->left;
681:     } else { /* if (j==3) */
682:       yt=b+hy*ys;
683:       xt=r;
684:       limit=rsize;
685:       boundary=user->right;
686:     }

688:     for (i=0; i<limit; i++){
689:       u1=xt;
690:       u2=-yt;
691:       for (k=0; k<maxits; k++){
692:         nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
693:         nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
694:         fnorm=PetscSqrtReal(nf1*nf1+nf2*nf2);
695:         if (fnorm <= tol) break;
696:         njac11=one+u2*u2-u1*u1;
697:         njac12=two*u1*u2;
698:         njac21=-two*u1*u2;
699:         njac22=-one - u1*u1 + u2*u2;
700:         det = njac11*njac22-njac21*njac12;
701:         u1 = u1-(njac22*nf1-njac12*nf2)/det;
702:         u2 = u2-(njac11*nf2-njac21*nf1)/det;
703:       }

705:       boundary[i]=u1*u1-u2*u2;
706:       if (j==0 || j==1) {
707:         xt=xt+hx;
708:       } else { /*  if (j==2 || j==3) */
709:         yt=yt+hy;
710:       }

712:     }

714:   }

716:   /* Scale the boundary if desired */
717:   if (1==1){
718:     PetscReal scl = 1.0;

720:     PetscOptionsGetReal(NULL,"-bottom",&scl,&flg);
721: 
722:     if (flg){
723:       for (i=0;i<bsize;i++) user->bottom[i]*=scl;
724:     }

726:     PetscOptionsGetReal(NULL,"-top",&scl,&flg);
727: 
728:     if (flg){
729:       for (i=0;i<tsize;i++) user->top[i]*=scl;
730:     }

732:     PetscOptionsGetReal(NULL,"-right",&scl,&flg);
733: 
734:     if (flg){
735:       for (i=0;i<rsize;i++) user->right[i]*=scl;
736:     }

738:     PetscOptionsGetReal(NULL,"-left",&scl,&flg);
739: 
740:     if (flg){
741:       for (i=0;i<lsize;i++) user->left[i]*=scl;
742:     }
743:   }

745:   return(0);
746: }

748: /* ------------------------------------------------------------------- */
751: /*
752:    MSA_InitialPoint - Calculates the initial guess in one of three ways.

754:    Input Parameters:
755: .  user - user-defined application context
756: .  X - vector for initial guess

758:    Output Parameters:
759: .  X - newly computed initial guess
760: */
761: static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
762: {
763:   PetscErrorCode      ierr;
764:   PetscInt   start2=-1,i,j;
765:   PetscReal   start1=0;
766:   PetscBool flg1,flg2;

769:   PetscOptionsGetReal(NULL,"-start",&start1,&flg1);
770:   PetscOptionsGetInt(NULL,"-random",&start2,&flg2);

772:   if (flg1){ /* The zero vector is reasonable */

774:     VecSet(X, start1);

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

778:     PetscRandom rctx;  PetscReal np5=-0.5;

780:     PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
781:     PetscRandomSetType(rctx,PETSCRAND);
782: 
783:     for (i=0; i<start2; i++){
784:       VecSetRandom(X, rctx);
785:     }
786:     PetscRandomDestroy(&rctx);
787:     VecShift(X, np5);

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

791:     PetscInt xs,xm,ys,ym;
792:     PetscInt mx=user->mx,my=user->my;
793:     PetscReal **x;

795:     /* Get local mesh boundaries */
796:     DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);

798:     /* Get pointers to vector data */
799:     DMDAVecGetArray(user->dm,X,(void**)&x);

801:     /* Perform local computations */
802:     for (j=ys; j<ys+ym; j++){
803:       for (i=xs; i< xs+xm; i++){
804:         x[j][i] = ( ((j+1)*user->bottom[i-xs+1]+(my-j+1)*user->top[i-xs+1])/(my+2)+
805:                    ((i+1)*user->left[j-ys+1]+(mx-i+1)*user->right[j-ys+1])/(mx+2))/2.0;
806:       }
807:     }

809:     /* Restore vectors */
810:     DMDAVecRestoreArray(user->dm,X,(void**)&x);

812:     PetscLogFlops(9*xm*ym);

814:   }
815:   return(0);
816: }

818: /*-----------------------------------------------------------------------*/
821: PetscErrorCode My_Monitor(Tao tao, void *ctx){
823:   Vec X;
825:   TaoGetSolutionVector(tao,&X);
826:   VecView(X,PETSC_VIEWER_STDOUT_WORLD);
827:   return(0);
828: }