Actual source code: minsurf1.c

      1: #include "tao.h"

      3: /*
      4:    Include "tao.h" so we can use the TAO solvers with PETSc support.
      5: */
      6: static char  help[] =
      7: "This example demonstrates use of the TAO package to\n\
      8: solve an unconstrained system of equations.  This example is based on a\n\
      9: problem from the MINPACK-2 test suite.  Given a rectangular 2-D domain and\n\
     10: boundary values along the edges of the domain, the objective is to find the\n\
     11: surface with the minimal area that satisfies the boundary conditions.\n\
     12: This application solves this problem using complimentarity -- We are actually\n\
     13: solving the system  (grad f)_i >= 0, if x_i == l_i \n\
     14:                     (grad f)_i = 0, if l_i < x_i < u_i \n\
     15:                     (grad f)_i <= 0, if x_i == u_i  \n\
     16: where f is the function to be minimized. \n\
     17: \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, and an average of the boundary conditions otherwise \n\n";

     23: /*T
     24:    Concepts: TAO - Solving a complementarity problem
     25:    Routines: TaoInitialize(); TaoFinalize();
     26:    Routines: TaoCreate(); TaoDestroy();
     27:    Routines: TaoApplicationCreate(); TaoAppDestroy();
     28:    Routines: TaoSetConstraintRoutine(); TaoSetJacobianRoutine();
     29:    Routines: TaoGetVariableBounds(); TaoAppSetInitialSolutionVec();
     30:    Routines: TaoSetOptions();
     31:    Routines: TaoSolve();
     32:    Processors: 1
     33: T*/

     36: /*
     37:    User-defined application context - contains data needed by the
     38:    application-provided call-back routines, FormFunctionGradient(),
     39:    FormHessian().
     40: */
     41: typedef struct {
     42:   int mx, my;
     43:   double *bottom, *top, *left, *right;
     44: } AppCtx;

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

     49: static int MSA_BoundaryConditions(AppCtx *);
     50: static int MSA_InitialPoint(AppCtx *, Vec);
     51: int FormConstraints(TAO_APPLICATION, Vec, Vec, void *);
     52: int FormJacobian(TAO_APPLICATION, Vec, Mat *, Mat*, MatStructure*,void *);

     56: int main(int argc, char **argv)
     57: {
     58:   int          info;              /* used to check for functions returning nonzeros */
     59:   Vec          x;                 /* solution vector */
     60:   Vec          c;                 /* Constraints function vector */
     61:   Vec          xl,xu;             /* Bounds on the variables */
     62:   PetscTruth   flg;               /* A return variable when checking for user options */
     63:   TAO_SOLVER   tao;               /* TAO_SOLVER solver context */
     64:   TAO_APPLICATION my_app;         /* The PETSc application */
     65:   Mat          J;                 /* Jacobian matrix */
     66:   int          N;                 /* Number of elements in vector */
     67:   PetscScalar lb = -TAO_INFINITY;      /* lower bound constant */
     68:   PetscScalar ub =  TAO_INFINITY;      /* upper bound constant */
     69:   AppCtx user;                    /* user-defined work context */

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

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

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

     82:   /* Calculate any derived values from parameters */
     83:   N = user.mx*user.my;

     85:
     86:   PetscPrintf(PETSC_COMM_SELF,"\n---- Minimum Surface Area Problem -----\n");
     87:   PetscPrintf(PETSC_COMM_SELF,"mx:%d, my:%d\n", user.mx,user.my);

     90:   /* Create appropriate vectors and matrices */
     91:   info = VecCreateSeq(MPI_COMM_SELF, N, &x);
     92:   info = VecDuplicate(x, &c); CHKERRQ(info);
     93:   info = MatCreateSeqAIJ(MPI_COMM_SELF, N, N, 7, PETSC_NULL, &J); CHKERRQ(info);

     95:   /* The TAO code begins here */

     97:   /* Create TAO solver and set desired solution method */
     98:   info = TaoCreate(PETSC_COMM_SELF, "tao_ssils", &tao); CHKERRQ(info);
     99:   info = TaoApplicationCreate(PETSC_COMM_SELF, &my_app); CHKERRQ(info);

    101:   /* Set data structure */
    102:   info = TaoAppSetInitialSolutionVec(my_app, x); CHKERRQ(info);

    104:   /*  Set routines for constraints function and Jacobian evaluation */
    105:   info = TaoAppSetConstraintRoutine(my_app, FormConstraints, (void *)&user);
    106:   CHKERRQ(info);
    107:   info = TaoAppSetJacobianRoutine(my_app, FormJacobian, (void *)&user); CHKERRQ(info);
    108:   info = TaoAppSetJacobianMat(my_app, J, J); CHKERRQ(info);

    110:   /* Set the variable bounds */
    111:   info = MSA_BoundaryConditions(&user); CHKERRQ(info);

    113:   /* Set initial solution guess */
    114:   info = MSA_InitialPoint(&user, x); CHKERRQ(info);

    116:   /* Now that the PETSc application is all set, attach to tao context */
    117:   info = TaoSetupApplicationSolver(my_app,tao); CHKERRQ(info);

    119:   /* Set Bounds on variables */
    120:   info = VecDuplicate(x, &xl); CHKERRQ(info);
    121:   info = VecDuplicate(x, &xu); CHKERRQ(info);
    122:   info = VecSet(xl, lb); CHKERRQ(info);
    123:   info = VecSet(xu, ub); CHKERRQ(info);
    124:   info = TaoAppSetVariableBounds(my_app,xl,xu); CHKERRQ(info);

    126:   /* Check for any tao command line options */
    127:   info = TaoSetOptions(my_app,tao); CHKERRQ(info);

    129:   /* Solve the application */
    130:   info = TaoSolveApplication(my_app,tao);  CHKERRQ(info);

    133:   /* Free Tao data structures */
    134:   info = TaoDestroy(tao); CHKERRQ(info);
    135:   info = TaoAppDestroy(my_app);  CHKERRQ(info);

    137:   /* Free PETSc data structures */
    138:   info = VecDestroy(x); CHKERRQ(info);
    139:   info = VecDestroy(xl); CHKERRQ(info);
    140:   info = VecDestroy(xu); CHKERRQ(info);
    141:   info = VecDestroy(c); CHKERRQ(info);
    142:   info = MatDestroy(J); CHKERRQ(info);

    144:   /* Free user-created data structures */
    145:   PetscFree(user.bottom);
    146:   PetscFree(user.top);
    147:   PetscFree(user.left);
    148:   PetscFree(user.right);

    150:   /* Finalize TAO and PETSc */
    151:   PetscFinalize();
    152:   TaoFinalize();

    154:   return 0;
    155: }

    157: /* -------------------------------------------------------------------- */

    161: /*  FormConstraints - Evaluates gradient of f.

    163:     Input Parameters:
    164: .   tao  - the TAO_APPLICATION context
    165: .   X    - input vector
    166: .   ptr  - optional user-defined context, as set by TaoAppSetConstraintRoutine()
    167:
    168:     Output Parameters:
    169: .   G - vector containing the newly evaluated gradient
    170: */
    171: int FormConstraints(TAO_APPLICATION tao, Vec X, Vec G, void *ptr){
    172:   AppCtx *user = (AppCtx *) ptr;
    173:   int    info,i,j,row;
    174:   int    mx=user->mx, my=user->my;
    175:   double hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
    176:   double f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
    177:   double df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
    178:   PetscScalar zero=0.0;
    179:   PetscScalar *g, *x;

    181:   /* Initialize vector to zero */
    182:   info = VecSet(G, zero); CHKERRQ(info);

    184:   /* Get pointers to vector data */
    185:   info = VecGetArray(X, &x); CHKERRQ(info);
    186:   info = VecGetArray(G, &g); CHKERRQ(info);

    188:   /* Compute function over the locally owned part of the mesh */
    189:   for (j=0; j<my; j++){
    190:     for (i=0; i< mx; i++){
    191:       row= j*mx + i;
    192:
    193:       xc = x[row];
    194:       xlt=xrb=xl=xr=xb=xt=xc;
    195:
    196:       if (i==0){ /* left side */
    197:         xl= user->left[j+1];
    198:         xlt = user->left[j+2];
    199:       } else {
    200:         xl = x[row-1];
    201:       }

    203:       if (j==0){ /* bottom side */
    204:         xb=user->bottom[i+1];
    205:         xrb = user->bottom[i+2];
    206:       } else {
    207:         xb = x[row-mx];
    208:       }
    209:
    210:       if (i+1 == mx){ /* right side */
    211:         xr=user->right[j+1];
    212:         xrb = user->right[j];
    213:       } else {
    214:         xr = x[row+1];
    215:       }

    217:       if (j+1==0+my){ /* top side */
    218:         xt=user->top[i+1];
    219:         xlt = user->top[i];
    220:       }else {
    221:         xt = x[row+mx];
    222:       }

    224:       if (i>0 && j+1<my){
    225:         xlt = x[row-1+mx];
    226:       }
    227:       if (j>0 && i+1<mx){
    228:         xrb = x[row+1-mx];
    229:       }

    231:       d1 = (xc-xl);
    232:       d2 = (xc-xr);
    233:       d3 = (xc-xt);
    234:       d4 = (xc-xb);
    235:       d5 = (xr-xrb);
    236:       d6 = (xrb-xb);
    237:       d7 = (xlt-xl);
    238:       d8 = (xt-xlt);
    239:
    240:       df1dxc = d1*hydhx;
    241:       df2dxc = ( d1*hydhx + d4*hxdhy );
    242:       df3dxc = d3*hxdhy;
    243:       df4dxc = ( d2*hydhx + d3*hxdhy );
    244:       df5dxc = d2*hydhx;
    245:       df6dxc = d4*hxdhy;

    247:       d1 /= hx;
    248:       d2 /= hx;
    249:       d3 /= hy;
    250:       d4 /= hy;
    251:       d5 /= hy;
    252:       d6 /= hx;
    253:       d7 /= hy;
    254:       d8 /= hx;

    256:       f1 = sqrt( 1.0 + d1*d1 + d7*d7);
    257:       f2 = sqrt( 1.0 + d1*d1 + d4*d4);
    258:       f3 = sqrt( 1.0 + d3*d3 + d8*d8);
    259:       f4 = sqrt( 1.0 + d3*d3 + d2*d2);
    260:       f5 = sqrt( 1.0 + d2*d2 + d5*d5);
    261:       f6 = sqrt( 1.0 + d4*d4 + d6*d6);
    262:
    263:       df1dxc /= f1;
    264:       df2dxc /= f2;
    265:       df3dxc /= f3;
    266:       df4dxc /= f4;
    267:       df5dxc /= f5;
    268:       df6dxc /= f6;

    270:       g[row] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc )/2.0;
    271:
    272:     }
    273:   }
    274:
    275:   /* Restore vectors */
    276:   info = VecRestoreArray(X, &x); CHKERRQ(info);
    277:   info = VecRestoreArray(G, &g); CHKERRQ(info);
    278:   info = PetscLogFlops(67*mx*my); CHKERRQ(info);
    279:   return 0;
    280: }

    282: /* ------------------------------------------------------------------- */
    285: /*
    286:    FormJacobian - Evaluates Jacobian matrix.

    288:    Input Parameters:
    289: .  tao  - the TAO_APPLICATION context
    290: .  X    - input vector
    291: .  ptr  - optional user-defined context, as set by TaoSetJacobian()

    293:    Output Parameters:
    294: .  tH    - Jacobian matrix

    296: */
    297: int FormJacobian(TAO_APPLICATION tao, Vec X, Mat *tH, Mat* tHPre, MatStructure* flag, void *ptr)
    298: {
    299:   AppCtx *user = (AppCtx *) ptr;
    300:   Mat H = *tH;
    301:   int    i,j,k,row,info;
    302:   int    mx=user->mx, my=user->my;
    303:   int    col[7];
    304:   double hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
    305:   double f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
    306:   double hl,hr,ht,hb,hc,htl,hbr;
    307:   PetscScalar *x, v[7];
    308:   PetscTruth assembled;

    310:   /* Set various matrix options */
    311:   info = MatSetOption(H,MAT_IGNORE_OFF_PROC_ENTRIES); CHKERRQ(info);
    312:   info = MatSetOption(H,MAT_COLUMNS_SORTED); CHKERRQ(info);
    313:   info = MatSetOption(H,MAT_ROWS_SORTED); CHKERRQ(info);
    314:   info = MatAssembled(H,&assembled); CHKERRQ(info);
    315:   if (assembled){info = MatZeroEntries(H);  CHKERRQ(info);}
    316:   *flag=SAME_NONZERO_PATTERN;

    318:   /* Get pointers to vector data */
    319:   info = VecGetArray(X, &x); CHKERRQ(info);

    321:   /* Compute Jacobian over the locally owned part of the mesh */
    322:   for (i=0; i< mx; i++){
    323:     for (j=0; j<my; j++){
    324:       row= j*mx + i;
    325:
    326:       xc = x[row];
    327:       xlt=xrb=xl=xr=xb=xt=xc;

    329:       /* Left side */
    330:       if (i==0){
    331:         xl= user->left[j+1];
    332:         xlt = user->left[j+2];
    333:       } else {
    334:         xl = x[row-1];
    335:       }
    336:
    337:       if (j==0){
    338:         xb=user->bottom[i+1];
    339:         xrb = user->bottom[i+2];
    340:       } else {
    341:         xb = x[row-mx];
    342:       }
    343:
    344:       if (i+1 == mx){
    345:         xr=user->right[j+1];
    346:         xrb = user->right[j];
    347:       } else {
    348:         xr = x[row+1];
    349:       }

    351:       if (j+1==my){
    352:         xt=user->top[i+1];
    353:         xlt = user->top[i];
    354:       }else {
    355:         xt = x[row+mx];
    356:       }

    358:       if (i>0 && j+1<my){
    359:         xlt = x[row-1+mx];
    360:       }
    361:       if (j>0 && i+1<mx){
    362:         xrb = x[row+1-mx];
    363:       }

    366:       d1 = (xc-xl)/hx;
    367:       d2 = (xc-xr)/hx;
    368:       d3 = (xc-xt)/hy;
    369:       d4 = (xc-xb)/hy;
    370:       d5 = (xrb-xr)/hy;
    371:       d6 = (xrb-xb)/hx;
    372:       d7 = (xlt-xl)/hy;
    373:       d8 = (xlt-xt)/hx;
    374:
    375:       f1 = sqrt( 1.0 + d1*d1 + d7*d7);
    376:       f2 = sqrt( 1.0 + d1*d1 + d4*d4);
    377:       f3 = sqrt( 1.0 + d3*d3 + d8*d8);
    378:       f4 = sqrt( 1.0 + d3*d3 + d2*d2);
    379:       f5 = sqrt( 1.0 + d2*d2 + d5*d5);
    380:       f6 = sqrt( 1.0 + d4*d4 + d6*d6);

    383:       hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
    384:         (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
    385:       hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
    386:         (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
    387:       ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
    388:         (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
    389:       hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
    390:         (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);

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

    395:       hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
    396:         hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
    397:         (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
    398:         (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);

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

    402:       k=0;
    403:       if (j>0){
    404:         v[k]=hb; col[k]=row - mx; k++;
    405:       }
    406:
    407:       if (j>0 && i < mx -1){
    408:         v[k]=hbr; col[k]=row - mx+1; k++;
    409:       }
    410:
    411:       if (i>0){
    412:         v[k]= hl; col[k]=row - 1; k++;
    413:       }
    414:
    415:       v[k]= hc; col[k]=row; k++;
    416:
    417:       if (i < mx-1 ){
    418:         v[k]= hr; col[k]=row+1; k++;
    419:       }
    420:
    421:       if (i>0 && j < my-1 ){
    422:         v[k]= htl; col[k] = row+mx-1; k++;
    423:       }
    424:
    425:       if (j < my-1 ){
    426:         v[k]= ht; col[k] = row+mx; k++;
    427:       }
    428:
    429:       /*
    430:          Set matrix values using local numbering, which was defined
    431:          earlier, in the main routine.
    432:       */
    433:       info = MatSetValues(H,1,&row,k,col,v,INSERT_VALUES);
    434:       CHKERRQ(info);
    435:     }
    436:   }
    437:
    438:   /* Restore vectors */
    439:   info = VecRestoreArray(X,&x); CHKERRQ(info);

    441:   /* Assemble the matrix */
    442:   info = MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
    443:   info = MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
    444:   info = PetscLogFlops(199*mx*my); CHKERRQ(info);
    445:   return 0;
    446: }

    448: /* ------------------------------------------------------------------- */
    451: /*
    452:    MSA_BoundaryConditions -  Calculates the boundary conditions for
    453:    the region.

    455:    Input Parameter:
    456: .  user - user-defined application context

    458:    Output Parameter:
    459: .  user - user-defined application context
    460: */
    461: static int MSA_BoundaryConditions(AppCtx * user)
    462: {
    463:   int        i,j,k,limit=0,info,maxits=5;
    464:   int        mx=user->mx,my=user->my;
    465:   int        bsize=0, lsize=0, tsize=0, rsize=0;
    466:   double     one=1.0, two=2.0, three=3.0, tol=1e-10;
    467:   double     fnorm,det,hx,hy,xt=0,yt=0;
    468:   double     u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
    469:   double     b=-0.5, t=0.5, l=-0.5, r=0.5;
    470:   double     *boundary;

    472:   bsize=mx+2; lsize=my+2; rsize=my+2; tsize=mx+2;

    474:   info = PetscMalloc(bsize*sizeof(double), &user->bottom);CHKERRQ(info);
    475:   info = PetscMalloc(tsize*sizeof(double), &user->top);CHKERRQ(info);
    476:   info = PetscMalloc(lsize*sizeof(double), &user->left);CHKERRQ(info);
    477:   info = PetscMalloc(rsize*sizeof(double), &user->right);CHKERRQ(info);

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

    481:   for (j=0; j<4; j++){
    482:     if (j==0){
    483:       yt=b;
    484:       xt=l;
    485:       limit=bsize;
    486:       boundary=user->bottom;
    487:     } else if (j==1){
    488:       yt=t;
    489:       xt=l;
    490:       limit=tsize;
    491:       boundary=user->top;
    492:     } else if (j==2){
    493:       yt=b;
    494:       xt=l;
    495:       limit=lsize;
    496:       boundary=user->left;
    497:     } else { // if  (j==3)
    498:       yt=b;
    499:       xt=r;
    500:       limit=rsize;
    501:       boundary=user->right;
    502:     }

    504:     for (i=0; i<limit; i++){
    505:       u1=xt;
    506:       u2=-yt;
    507:       for (k=0; k<maxits; k++){
    508:         nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
    509:         nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
    510:         fnorm=sqrt(nf1*nf1+nf2*nf2);
    511:         if (fnorm <= tol) break;
    512:         njac11=one+u2*u2-u1*u1;
    513:         njac12=two*u1*u2;
    514:         njac21=-two*u1*u2;
    515:         njac22=-one - u1*u1 + u2*u2;
    516:         det = njac11*njac22-njac21*njac12;
    517:         u1 = u1-(njac22*nf1-njac12*nf2)/det;
    518:         u2 = u2-(njac11*nf2-njac21*nf1)/det;
    519:       }

    521:       boundary[i]=u1*u1-u2*u2;
    522:       if (j==0 || j==1) {
    523:         xt=xt+hx;
    524:       } else { // if (j==2 || j==3)
    525:         yt=yt+hy;
    526:       }
    527:     }
    528:   }

    530:   return 0;
    531: }

    533: /* ------------------------------------------------------------------- */
    536: /*
    537:    MSA_InitialPoint - Calculates the initial guess in one of three ways.

    539:    Input Parameters:
    540: .  user - user-defined application context
    541: .  X - vector for initial guess

    543:    Output Parameters:
    544: .  X - newly computed initial guess
    545: */
    546: static int MSA_InitialPoint(AppCtx * user, Vec X)
    547: {
    548:   int           start=-1,i,j,info;
    549:   PetscScalar   zero=0.0;
    550:   PetscTruth    flg;

    552:   info = PetscOptionsGetInt(PETSC_NULL,"-start",&start,&flg); CHKERRQ(info);

    554:   if (flg && start==0){ /* The zero vector is reasonable */
    555:
    556:     info = VecSet(X, zero); CHKERRQ(info);
    557:     /* PLogInfo(user,"Min. Surface Area Problem: Start with 0 vector \n"); */

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

    562:     int    row;
    563:     int    mx=user->mx,my=user->my;
    564:     PetscScalar *x;
    565:
    566:     /* Get pointers to vector data */
    567:     info = VecGetArray(X,&x); CHKERRQ(info);

    569:     /* Perform local computations */
    570:     for (j=0; j<my; j++){
    571:       for (i=0; i< mx; i++){
    572:         row=(j)*mx + (i);
    573:         x[row] = ( ((j+1)*user->bottom[i+1]+(my-j+1)*user->top[i+1])/(my+2)+
    574:                    ((i+1)*user->left[j+1]+(mx-i+1)*user->right[j+1])/(mx+2))/2.0;
    575:       }
    576:     }
    577:
    578:     /* Restore vectors */
    579:     info = VecRestoreArray(X,&x); CHKERRQ(info);
    580:
    581:   }
    582:   return 0;
    583: }