Actual source code: minsurf2.c

      1: /*$Id: minsurf2.c 1.80 05/05/11 08:44:55-05:00 sarich@zorak.(none) $*/

      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: TaoGetKSP(); 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:   int         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:   int             Nx, Ny;              /* number of processors in x- and y- directions */
     69:   int             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);

    132:   if (fdcoloring) {
    133:     info = TaoAppSetColoring(minsurfapp, user.iscoloring); CHKERRQ(info);
    134:     info = TaoAppSetHessianRoutine(minsurfapp,TaoAppDefaultComputeHessianColor,(void *)&user); CHKERRQ(info);
    135:   } else if (fddefault){
    136:     info = TaoAppSetHessianRoutine(minsurfapp,TaoAppDefaultComputeHessian,(void *)&user); CHKERRQ(info);
    137:   } else {
    138:     info = TaoAppSetHessianRoutine(minsurfapp,FormHessian,(void *)&user); CHKERRQ(info);
    139:   }

    141:   /*
    142:      Create a matrix data structure to store the Hessian and set
    143:      the Hessian evalution routine.
    144:      Set the matrix structure to be used for Hessian evalutions
    145:   */
    146:   info = DAGetMatrix(user.da,MATAIJ,&user.H);CHKERRQ(info);
    147:   info = DAGetColoring(user.da,IS_COLORING_GHOSTED,&user.iscoloring);
    148:   CHKERRQ(info);
    149:   info = MatSetOption(user.H,MAT_SYMMETRIC); CHKERRQ(info);

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

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

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

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

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

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

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

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

    191:   /* Free PETSc data structures */
    192:   info = VecDestroy(x); CHKERRQ(info);
    193:   info = MatDestroy(user.H); CHKERRQ(info);
    194:   info = ISColoringDestroy(user.iscoloring);CHKERRQ(info);
    195:   PetscFree(user.bottom);
    196:   PetscFree(user.top);
    197:   PetscFree(user.left);
    198:   PetscFree(user.right);
    199:   info = DADestroy(user.da); CHKERRQ(info);

    201:   /* Finalize TAO */
    202:   TaoFinalize();
    203:   PetscFinalize();
    204:
    205:   return 0;
    206: }

    210: int FormGradient(TAO_APPLICATION taoapp, Vec X, Vec G,void *userCtx){
    211:   int info;
    212:   double fcn;
    213:   TaoFunctionBegin;
    214:   info = FormFunctionGradient(taoapp,X,&fcn,G,userCtx);CHKERRQ(info);
    215:   TaoFunctionReturn(0);
    216: }

    218: /* -------------------------------------------------------------------- */
    221: /*  FormFunctionGradient - Evaluates the function and corresponding gradient.

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

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

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

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

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

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

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

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

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

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

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

    318:       d1 *= rhx;
    319:       d2 *= rhx;
    320:       d3 *= rhy;
    321:       d4 *= rhy;
    322:       d5 *= rhy;
    323:       d6 *= rhx;
    324:       d7 *= rhy;
    325:       d8 *= rhx;

    327:       f1 = sqrt( 1.0 + d1*d1 + d7*d7);
    328:       f2 = sqrt( 1.0 + d1*d1 + d4*d4);
    329:       f3 = sqrt( 1.0 + d3*d3 + d8*d8);
    330:       f4 = sqrt( 1.0 + d3*d3 + d2*d2);
    331:       f5 = sqrt( 1.0 + d2*d2 + d5*d5);
    332:       f6 = sqrt( 1.0 + d4*d4 + d6*d6);
    333:
    334:       f2 = sqrt( 1.0 + d1*d1 + d4*d4);
    335:       f4 = sqrt( 1.0 + d3*d3 + d2*d2);

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

    339:       df1dxc /= f1;
    340:       df2dxc /= f2;
    341:       df3dxc /= f3;
    342:       df4dxc /= f4;
    343:       df5dxc /= f5;
    344:       df6dxc /= f6;

    346:       g[j][i] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc ) * 0.5;
    347:
    348:     }
    349:   }

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

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

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

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

    396:   /* Restore vectors */
    397:   info = DAVecRestoreArray(user->da,localX,(void**)&x);
    398:   info = DAVecRestoreArray(user->da,G,(void**)&g);

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

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

    405:   return 0;
    406: }

    408: /* ------------------------------------------------------------------- */
    411: /*
    412:    FormHessian - Evaluates Hessian matrix.

    414:    Input Parameters:
    415: .  taoapp  - the TAO_APPLICATION context
    416: .  x    - input vector
    417: .  ptr  - optional user-defined context, as set by TaoSetHessian()

    419:    Output Parameters:
    420: .  H    - Hessian matrix
    421: .  Hpre - optionally different preconditioning matrix
    422: .  flg  - flag indicating matrix structure

    424: */
    425: int FormHessian(TAO_APPLICATION taoapp,Vec X,Mat *H, Mat *Hpre, MatStructure *flg, void *ptr)
    426: {
    427:   int    info;
    428:   AppCtx *user = (AppCtx *) ptr;

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

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

    439:   return 0;
    440: }

    442: /* ------------------------------------------------------------------- */
    445: /*
    446:    QuadraticH - Evaluates Hessian matrix.

    448:    Input Parameters:
    449: .  user - user-defined context, as set by TaoSetHessian()
    450: .  X    - input vector

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

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

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

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

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

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

    486:   /* Set various matrix options */
    487:   info = MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES); CHKERRQ(info);
    488:   info = MatSetOption(Hessian,MAT_COLUMNS_SORTED); CHKERRQ(info);
    489:   info = MatSetOption(Hessian,MAT_ROWS_SORTED); CHKERRQ(info);

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

    493:   for (j=ys; j<ys+ym; j++){
    494:
    495:     for (i=xs; i< xs+xm; i++){

    497:       xc = x[j][i];
    498:       xlt=xrb=xl=xr=xb=xt=xc;

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

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

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

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

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

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

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

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

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

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

    630:   /* Assemble the matrix */
    631:   info = MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
    632:   info = MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY); CHKERRQ(info);

    634:   info = PetscLogFlops(199*xm*ym); CHKERRQ(info);
    635:   return 0;
    636: }

    638: /* ------------------------------------------------------------------- */
    641: /*
    642:    MSA_BoundaryConditions -  Calculates the boundary conditions for
    643:    the region.

    645:    Input Parameter:
    646: .  user - user-defined application context

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

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

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

    673:   info = PetscMalloc(bsize*sizeof(double),&user->bottom); CHKERRQ(info);
    674:   info = PetscMalloc(tsize*sizeof(double),&user->top); CHKERRQ(info);
    675:   info = PetscMalloc(lsize*sizeof(double),&user->left); CHKERRQ(info);
    676:   info = PetscMalloc(rsize*sizeof(double),&user->right); CHKERRQ(info);

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

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

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

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

    729:   }

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

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

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

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

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

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

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

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

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

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

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

    791:     PetscRandom rctx;  PetscScalar np5=-0.5;

    793:     info = PetscRandomCreate(PETSC_COMM_WORLD,RANDOM_DEFAULT,&rctx);
    794:     CHKERRQ(info);
    795:     for (i=0; i<start2; i++){
    796:       info = VecSetRandom(X, rctx); CHKERRQ(info);
    797:     }
    798:     info = PetscRandomDestroy(rctx); CHKERRQ(info);
    799:     info = VecShift(X, np5); CHKERRQ(info);

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

    803:     int    xs,xm,ys,ym;
    804:     int    mx=user->mx,my=user->my;
    805:     PetscScalar **x;
    806:
    807:     /* Get local mesh boundaries */
    808:     info = DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);
    809:
    810:     /* Get pointers to vector data */
    811:     info = DAVecGetArray(user->da,X,(void**)&x);

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

    824:     info = PetscLogFlops(9*xm*ym); CHKERRQ(info);
    825:
    826:   }
    827:   return 0;
    828: }

    830: /*-----------------------------------------------------------------------*/
    833: int My_Monitor(TAO_APPLICATION minsurfapp, void *ctx){
    834:   int info;
    835:   Vec X;

    837:   info = TaoAppGetSolutionVec(minsurfapp,&X); CHKERRQ(info);
    838:   info = VecView(X,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(info);
    839:   return 0;
    840: }