Actual source code: minsurf2.c

      1: /* Program usage: mpiexec -n <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: TaoDestroy();
     31:    Processors: n
     32: T*/

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

     49: /* -------- User-defined Routines --------- */

     51: static PetscErrorCode MSA_BoundaryConditions(AppCtx*);
     52: static PetscErrorCode MSA_InitialPoint(AppCtx*,Vec);
     53: PetscErrorCode QuadraticH(AppCtx*,Vec,Mat);
     54: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void*);
     55: PetscErrorCode FormGradient(Tao,Vec,Vec,void*);
     56: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat,void*);
     57: 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:   Tao                tao;                 /* TAO solver context */
     67:   AppCtx             user;                /* user-defined work context */
     68:   ISColoring         iscoloring;
     69:   MatFDColoring      matfdcoloring;

     71:   /* Initialize TAO */
     72:   PetscInitialize(&argc, &argv,(char *)0,help);if (ierr) return ierr;

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

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

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

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

     88:   /* Create distributed array (DM) to manage parallel grid and vectors  */
     89:   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);
     90:   DMSetFromOptions(user.dm);
     91:   DMSetUp(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,NULL,"-tao_fddefault",&fddefault);
    119:   PetscOptionsHasName(NULL,NULL,"-tao_fdcoloring",&fdcoloring);

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

    130:   if (fdcoloring) {
    131:     DMCreateColoring(user.dm,IS_COLORING_GLOBAL,&iscoloring);
    132:     MatFDColoringCreate(user.H,iscoloring,&matfdcoloring);
    133:     ISColoringDestroy(&iscoloring);
    134:     MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormGradient,(void*)&user);
    135:     MatFDColoringSetFromOptions(matfdcoloring);
    136:     TaoSetHessianRoutine(tao,user.H,user.H,TaoDefaultComputeHessianColor,(void *)matfdcoloring);
    137:   } else if (fddefault){
    138:     TaoSetHessianRoutine(tao,user.H,user.H,TaoDefaultComputeHessian,(void *)NULL);
    139:   } else {
    140:     TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void *)&user);
    141:   }

    143:   /*
    144:      If my_monitor option is in command line, then use the user-provided
    145:      monitoring function
    146:   */
    147:   PetscOptionsHasName(NULL,NULL,"-my_monitor",&viewmat);
    148:   if (viewmat){
    149:     TaoSetMonitor(tao,My_Monitor,NULL,NULL);
    150:   }

    152:   /* Check for any tao command line options */
    153:   TaoSetFromOptions(tao);

    155:   /* SOLVE THE APPLICATION */
    156:   TaoSolve(tao);

    158:   TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);

    160:   /* Free TAO data structures */
    161:   TaoDestroy(&tao);

    163:   /* Free PETSc data structures */
    164:   VecDestroy(&x);
    165:   MatDestroy(&user.H);
    166:   if (fdcoloring) {
    167:     MatFDColoringDestroy(&matfdcoloring);
    168:   }
    169:   PetscFree(user.bottom);
    170:   PetscFree(user.top);
    171:   PetscFree(user.left);
    172:   PetscFree(user.right);
    173:   DMDestroy(&user.dm);
    174:   PetscFinalize();
    175:   return ierr;
    176: }

    178: PetscErrorCode FormGradient(Tao tao, Vec X, Vec G,void *userCtx)
    179: {
    181:   PetscReal      fcn;

    184:   FormFunctionGradient(tao,X,&fcn,G,userCtx);
    185:   return(0);
    186: }

    188: /* -------------------------------------------------------------------- */
    189: /*  FormFunctionGradient - Evaluates the function and corresponding gradient.

    191:     Input Parameters:
    192: .   tao     - the Tao context
    193: .   XX      - input vector
    194: .   userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()

    196:     Output Parameters:
    197: .   fcn     - the newly evaluated function
    198: .   GG       - vector containing the newly evaluated gradient
    199: */
    200: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn,Vec G,void *userCtx){

    202:   AppCtx         *user = (AppCtx *) userCtx;
    204:   PetscInt       i,j;
    205:   PetscInt       mx=user->mx, my=user->my;
    206:   PetscInt       xs,xm,gxs,gxm,ys,ym,gys,gym;
    207:   PetscReal      ft=0.0;
    208:   PetscReal      hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy;
    209:   PetscReal      rhx=mx+1, rhy=my+1;
    210:   PetscReal      f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
    211:   PetscReal      df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
    212:   PetscReal      **g, **x;
    213:   Vec            localX;

    216:   /* Get local mesh boundaries */
    217:   DMGetLocalVector(user->dm,&localX);
    218:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
    219:   DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

    221:   /* Scatter ghost points to local vector */
    222:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
    223:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);

    225:   /* Get pointers to vector data */
    226:   DMDAVecGetArray(user->dm,localX,(void**)&x);
    227:   DMDAVecGetArray(user->dm,G,(void**)&g);

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

    233:       xc = x[j][i];
    234:       xlt=xrb=xl=xr=xb=xt=xc;

    236:       if (i==0){ /* left side */
    237:         xl= user->left[j-ys+1];
    238:         xlt = user->left[j-ys+2];
    239:       } else {
    240:         xl = x[j][i-1];
    241:       }

    243:       if (j==0){ /* bottom side */
    244:         xb=user->bottom[i-xs+1];
    245:         xrb =user->bottom[i-xs+2];
    246:       } else {
    247:         xb = x[j-1][i];
    248:       }

    250:       if (i+1 == gxs+gxm){ /* right side */
    251:         xr=user->right[j-ys+1];
    252:         xrb = user->right[j-ys];
    253:       } else {
    254:         xr = x[j][i+1];
    255:       }

    257:       if (j+1==gys+gym){ /* top side */
    258:         xt=user->top[i-xs+1];
    259:         xlt = user->top[i-xs];
    260:       }else {
    261:         xt = x[j+1][i];
    262:       }

    264:       if (i>gxs && j+1<gys+gym){
    265:         xlt = x[j+1][i-1];
    266:       }
    267:       if (j>gys && i+1<gxs+gxm){
    268:         xrb = x[j-1][i+1];
    269:       }

    271:       d1 = (xc-xl);
    272:       d2 = (xc-xr);
    273:       d3 = (xc-xt);
    274:       d4 = (xc-xb);
    275:       d5 = (xr-xrb);
    276:       d6 = (xrb-xb);
    277:       d7 = (xlt-xl);
    278:       d8 = (xt-xlt);

    280:       df1dxc = d1*hydhx;
    281:       df2dxc = (d1*hydhx + d4*hxdhy);
    282:       df3dxc = d3*hxdhy;
    283:       df4dxc = (d2*hydhx + d3*hxdhy);
    284:       df5dxc = d2*hydhx;
    285:       df6dxc = d4*hxdhy;

    287:       d1 *= rhx;
    288:       d2 *= rhx;
    289:       d3 *= rhy;
    290:       d4 *= rhy;
    291:       d5 *= rhy;
    292:       d6 *= rhx;
    293:       d7 *= rhy;
    294:       d8 *= rhx;

    296:       f1 = PetscSqrtReal(1.0 + d1*d1 + d7*d7);
    297:       f2 = PetscSqrtReal(1.0 + d1*d1 + d4*d4);
    298:       f3 = PetscSqrtReal(1.0 + d3*d3 + d8*d8);
    299:       f4 = PetscSqrtReal(1.0 + d3*d3 + d2*d2);
    300:       f5 = PetscSqrtReal(1.0 + d2*d2 + d5*d5);
    301:       f6 = PetscSqrtReal(1.0 + d4*d4 + d6*d6);

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

    305:       df1dxc /= f1;
    306:       df2dxc /= f2;
    307:       df3dxc /= f3;
    308:       df4dxc /= f4;
    309:       df5dxc /= f5;
    310:       df6dxc /= f6;

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

    314:     }
    315:   }

    317:   /* Compute triangular areas along the border of the domain. */
    318:   if (xs==0){ /* left side */
    319:     for (j=ys; j<ys+ym; j++){
    320:       d3=(user->left[j-ys+1] - user->left[j-ys+2])*rhy;
    321:       d2=(user->left[j-ys+1] - x[j][0]) *rhx;
    322:       ft = ft+PetscSqrtReal(1.0 + d3*d3 + d2*d2);
    323:     }
    324:   }
    325:   if (ys==0){ /* bottom side */
    326:     for (i=xs; i<xs+xm; i++){
    327:       d2=(user->bottom[i+1-xs]-user->bottom[i-xs+2])*rhx;
    328:       d3=(user->bottom[i-xs+1]-x[0][i])*rhy;
    329:       ft = ft+PetscSqrtReal(1.0 + d3*d3 + d2*d2);
    330:     }
    331:   }

    333:   if (xs+xm==mx){ /* right side */
    334:     for (j=ys; j< ys+ym; j++){
    335:       d1=(x[j][mx-1] - user->right[j-ys+1])*rhx;
    336:       d4=(user->right[j-ys]-user->right[j-ys+1])*rhy;
    337:       ft = ft+PetscSqrtReal(1.0 + d1*d1 + d4*d4);
    338:     }
    339:   }
    340:   if (ys+ym==my){ /* top side */
    341:     for (i=xs; i<xs+xm; i++){
    342:       d1=(x[my-1][i] - user->top[i-xs+1])*rhy;
    343:       d4=(user->top[i-xs+1] - user->top[i-xs])*rhx;
    344:       ft = ft+PetscSqrtReal(1.0 + d1*d1 + d4*d4);
    345:     }
    346:   }

    348:   if (ys==0 && xs==0){
    349:     d1=(user->left[0]-user->left[1])/hy;
    350:     d2=(user->bottom[0]-user->bottom[1])*rhx;
    351:     ft +=PetscSqrtReal(1.0 + d1*d1 + d2*d2);
    352:   }
    353:   if (ys+ym == my && xs+xm == mx){
    354:     d1=(user->right[ym+1] - user->right[ym])*rhy;
    355:     d2=(user->top[xm+1] - user->top[xm])*rhx;
    356:     ft +=PetscSqrtReal(1.0 + d1*d1 + d2*d2);
    357:   }

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

    362:   /* Restore vectors */
    363:   DMDAVecRestoreArray(user->dm,localX,(void**)&x);
    364:   DMDAVecRestoreArray(user->dm,G,(void**)&g);

    366:   /* Scatter values to global vector */
    367:   DMRestoreLocalVector(user->dm,&localX);
    368:   PetscLogFlops(67.0*xm*ym);
    369:   return(0);
    370: }

    372: /* ------------------------------------------------------------------- */
    373: /*
    374:    FormHessian - Evaluates Hessian matrix.

    376:    Input Parameters:
    377: .  tao  - the Tao context
    378: .  x    - input vector
    379: .  ptr  - optional user-defined context, as set by TaoSetHessianRoutine()

    381:    Output Parameters:
    382: .  H    - Hessian matrix
    383: .  Hpre - optionally different preconditioning matrix
    384: .  flg  - flag indicating matrix structure

    386: */
    387: PetscErrorCode FormHessian(Tao tao,Vec X,Mat H, Mat Hpre, void *ptr)
    388: {
    390:   AppCtx         *user = (AppCtx *) ptr;

    393:   /* Evaluate the Hessian entries*/
    394:   QuadraticH(user,X,H);
    395:   return(0);
    396: }

    398: /* ------------------------------------------------------------------- */
    399: /*
    400:    QuadraticH - Evaluates Hessian matrix.

    402:    Input Parameters:
    403: .  user - user-defined context, as set by TaoSetHessian()
    404: .  X    - input vector

    406:    Output Parameter:
    407: .  H    - Hessian matrix
    408: */
    409: PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian)
    410: {
    411:   PetscErrorCode    ierr;
    412:   PetscInt i,j,k;
    413:   PetscInt mx=user->mx, my=user->my;
    414:   PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
    415:   PetscReal hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
    416:   PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
    417:   PetscReal hl,hr,ht,hb,hc,htl,hbr;
    418:   PetscReal **x, v[7];
    419:   MatStencil col[7],row;
    420:   Vec    localX;
    421:   PetscBool assembled;

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

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

    430:   /* Scatter ghost points to local vector */
    431:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
    432:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);

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

    437:   /* Initialize matrix entries to zero */
    438:   MatAssembled(Hessian,&assembled);
    439:   if (assembled){MatZeroEntries(Hessian);}

    442:   /* Set various matrix options */
    443:   MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);

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

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

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

    451:       xc = x[j][i];
    452:       xlt=xrb=xl=xr=xb=xt=xc;

    454:       /* Left side */
    455:       if (i==0){
    456:         xl  = user->left[j-ys+1];
    457:         xlt = user->left[j-ys+2];
    458:       } else {
    459:         xl  = x[j][i-1];
    460:       }

    462:       if (j==0){
    463:         xb  = user->bottom[i-xs+1];
    464:         xrb = user->bottom[i-xs+2];
    465:       } else {
    466:         xb  = x[j-1][i];
    467:       }

    469:       if (i+1 == mx){
    470:         xr  = user->right[j-ys+1];
    471:         xrb = user->right[j-ys];
    472:       } else {
    473:         xr  = x[j][i+1];
    474:       }

    476:       if (j+1==my){
    477:         xt  = user->top[i-xs+1];
    478:         xlt = user->top[i-xs];
    479:       }else {
    480:         xt  = x[j+1][i];
    481:       }

    483:       if (i>0 && j+1<my){
    484:         xlt = x[j+1][i-1];
    485:       }
    486:       if (j>0 && i+1<mx){
    487:         xrb = x[j-1][i+1];
    488:       }

    491:       d1 = (xc-xl)/hx;
    492:       d2 = (xc-xr)/hx;
    493:       d3 = (xc-xt)/hy;
    494:       d4 = (xc-xb)/hy;
    495:       d5 = (xrb-xr)/hy;
    496:       d6 = (xrb-xb)/hx;
    497:       d7 = (xlt-xl)/hy;
    498:       d8 = (xlt-xt)/hx;

    500:       f1 = PetscSqrtReal(1.0 + d1*d1 + d7*d7);
    501:       f2 = PetscSqrtReal(1.0 + d1*d1 + d4*d4);
    502:       f3 = PetscSqrtReal(1.0 + d3*d3 + d8*d8);
    503:       f4 = PetscSqrtReal(1.0 + d3*d3 + d2*d2);
    504:       f5 = PetscSqrtReal(1.0 + d2*d2 + d5*d5);
    505:       f6 = PetscSqrtReal(1.0 + d4*d4 + d6*d6);

    508:       hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
    509:         (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
    510:       hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
    511:         (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
    512:       ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
    513:         (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
    514:       hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
    515:         (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);

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

    520:       hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
    521:         hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
    522:         (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
    523:         (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);

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

    527:       row.j = j; row.i = i;
    528:       k=0;
    529:       if (j>0){
    530:         v[k]=hb;
    531:         col[k].j = j - 1; col[k].i = i;
    532:         k++;
    533:       }

    535:       if (j>0 && i < mx -1){
    536:         v[k]=hbr;
    537:         col[k].j = j - 1; col[k].i = i+1;
    538:         k++;
    539:       }

    541:       if (i>0){
    542:         v[k]= hl;
    543:         col[k].j = j; col[k].i = i-1;
    544:         k++;
    545:       }

    547:       v[k]= hc;
    548:       col[k].j = j; col[k].i = i;
    549:       k++;

    551:       if (i < mx-1){
    552:         v[k]= hr;
    553:         col[k].j = j; col[k].i = i+1;
    554:         k++;
    555:       }

    557:       if (i>0 && j < my-1){
    558:         v[k]= htl;
    559:         col[k].j = j+1; col[k].i = i-1;
    560:         k++;
    561:       }

    563:       if (j < my-1){
    564:         v[k]= ht;
    565:         col[k].j = j+1; col[k].i = i;
    566:         k++;
    567:       }

    569:       /*
    570:          Set matrix values using local numbering, which was defined
    571:          earlier, in the main routine.
    572:       */
    573:       MatSetValuesStencil(Hessian,1,&row,k,col,v,INSERT_VALUES);
    574:     }
    575:   }

    577:   DMDAVecRestoreArray(user->dm,localX,(void**)&x);
    578:   DMRestoreLocalVector(user->dm,&localX);

    580:   MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY);
    581:   MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY);

    583:   PetscLogFlops(199.0*xm*ym);
    584:   return(0);
    585: }

    587: /* ------------------------------------------------------------------- */
    588: /*
    589:    MSA_BoundaryConditions -  Calculates the boundary conditions for
    590:    the region.

    592:    Input Parameter:
    593: .  user - user-defined application context

    595:    Output Parameter:
    596: .  user - user-defined application context
    597: */
    598: static PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
    599: {
    601:   PetscInt       i,j,k,limit=0,maxits=5;
    602:   PetscInt       xs,ys,xm,ym,gxs,gys,gxm,gym;
    603:   PetscInt       mx=user->mx,my=user->my;
    604:   PetscInt       bsize=0, lsize=0, tsize=0, rsize=0;
    605:   PetscReal      one=1.0, two=2.0, three=3.0, tol=1e-10;
    606:   PetscReal      fnorm,det,hx,hy,xt=0,yt=0;
    607:   PetscReal      u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
    608:   PetscReal      b=-0.5, t=0.5, l=-0.5, r=0.5;
    609:   PetscReal      *boundary;
    610:   PetscBool      flg;

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

    617:   bsize=xm+2;
    618:   lsize=ym+2;
    619:   rsize=ym+2;
    620:   tsize=xm+2;

    622:   PetscMalloc1(bsize,&user->bottom);
    623:   PetscMalloc1(tsize,&user->top);
    624:   PetscMalloc1(lsize,&user->left);
    625:   PetscMalloc1(rsize,&user->right);

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

    629:   for (j=0; j<4; j++){
    630:     if (j==0){
    631:       yt=b;
    632:       xt=l+hx*xs;
    633:       limit=bsize;
    634:       boundary=user->bottom;
    635:     } else if (j==1){
    636:       yt=t;
    637:       xt=l+hx*xs;
    638:       limit=tsize;
    639:       boundary=user->top;
    640:     } else if (j==2){
    641:       yt=b+hy*ys;
    642:       xt=l;
    643:       limit=lsize;
    644:       boundary=user->left;
    645:     } else { /* if (j==3) */
    646:       yt=b+hy*ys;
    647:       xt=r;
    648:       limit=rsize;
    649:       boundary=user->right;
    650:     }

    652:     for (i=0; i<limit; i++){
    653:       u1=xt;
    654:       u2=-yt;
    655:       for (k=0; k<maxits; k++){
    656:         nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
    657:         nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
    658:         fnorm=PetscSqrtReal(nf1*nf1+nf2*nf2);
    659:         if (fnorm <= tol) break;
    660:         njac11=one+u2*u2-u1*u1;
    661:         njac12=two*u1*u2;
    662:         njac21=-two*u1*u2;
    663:         njac22=-one - u1*u1 + u2*u2;
    664:         det = njac11*njac22-njac21*njac12;
    665:         u1 = u1-(njac22*nf1-njac12*nf2)/det;
    666:         u2 = u2-(njac11*nf2-njac21*nf1)/det;
    667:       }

    669:       boundary[i]=u1*u1-u2*u2;
    670:       if (j==0 || j==1) {
    671:         xt=xt+hx;
    672:       } else { /*  if (j==2 || j==3) */
    673:         yt=yt+hy;
    674:       }

    676:     }

    678:   }

    680:   /* Scale the boundary if desired */
    681:   if (1==1){
    682:     PetscReal scl = 1.0;

    684:     PetscOptionsGetReal(NULL,NULL,"-bottom",&scl,&flg);
    685:     if (flg){
    686:       for (i=0;i<bsize;i++) user->bottom[i]*=scl;
    687:     }

    689:     PetscOptionsGetReal(NULL,NULL,"-top",&scl,&flg);
    690:     if (flg){
    691:       for (i=0;i<tsize;i++) user->top[i]*=scl;
    692:     }

    694:     PetscOptionsGetReal(NULL,NULL,"-right",&scl,&flg);
    695:     if (flg){
    696:       for (i=0;i<rsize;i++) user->right[i]*=scl;
    697:     }

    699:     PetscOptionsGetReal(NULL,NULL,"-left",&scl,&flg);
    700:     if (flg){
    701:       for (i=0;i<lsize;i++) user->left[i]*=scl;
    702:     }
    703:   }
    704:   return(0);
    705: }

    707: /* ------------------------------------------------------------------- */
    708: /*
    709:    MSA_InitialPoint - Calculates the initial guess in one of three ways.

    711:    Input Parameters:
    712: .  user - user-defined application context
    713: .  X - vector for initial guess

    715:    Output Parameters:
    716: .  X - newly computed initial guess
    717: */
    718: static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
    719: {
    721:   PetscInt       start2=-1,i,j;
    722:   PetscReal      start1=0;
    723:   PetscBool      flg1,flg2;

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

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

    731:     VecSet(X, start1);

    733:   } else if (flg2 && start2>0){ /* Try a random start between -0.5 and 0.5 */
    734:     PetscRandom rctx;  PetscReal np5=-0.5;

    736:     PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
    737:     VecSetRandom(X, rctx);
    738:     PetscRandomDestroy(&rctx);
    739:     VecShift(X, np5);

    741:   } else { /* Take an average of the boundary conditions */
    742:     PetscInt  xs,xm,ys,ym;
    743:     PetscInt  mx=user->mx,my=user->my;
    744:     PetscReal **x;

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

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

    752:     /* Perform local computations */
    753:     for (j=ys; j<ys+ym; j++){
    754:       for (i=xs; i< xs+xm; i++){
    755:         x[j][i] = (((j+1)*user->bottom[i-xs+1]+(my-j+1)*user->top[i-xs+1])/(my+2)+((i+1)*user->left[j-ys+1]+(mx-i+1)*user->right[j-ys+1])/(mx+2))/2.0;
    756:       }
    757:     }
    758:     DMDAVecRestoreArray(user->dm,X,(void**)&x);
    759:     PetscLogFlops(9.0*xm*ym);
    760:   }
    761:   return(0);
    762: }

    764: /*-----------------------------------------------------------------------*/
    765: PetscErrorCode My_Monitor(Tao tao, void *ctx)
    766: {
    768:   Vec            X;

    771:   TaoGetSolutionVector(tao,&X);
    772:   VecView(X,PETSC_VIEWER_STDOUT_WORLD);
    773:   return(0);
    774: }

    777: /*TEST

    779:    build:
    780:       requires: !complex

    782:    test:
    783:       args: -tao_smonitor -tao_type lmvm -mx 10 -my 8 -tao_gatol 1.e-3
    784:       requires: !single

    786:    test:
    787:       suffix: 2
    788:       nsize: 2
    789:       args: -tao_smonitor -tao_type nls -tao_nls_ksp_max_it 15 -tao_gatol 1.e-4
    790:       filter: grep -v "nls ksp"
    791:       requires: !single

    793:    test:
    794:       suffix: 3
    795:       nsize: 3
    796:       args: -tao_smonitor -tao_type cg -tao_cg_type fr -mx 10 -my 10 -tao_gatol 1.e-3
    797:       requires: !single

    799:    test:
    800:       suffix: 5
    801:       nsize: 2
    802:       args: -tao_smonitor -tao_type bmrm -mx 10 -my 8 -tao_gatol 1.e-3
    803:       requires: !single

    805: TEST*/