Actual source code: rosenbrock1.c

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

      3: /* Program usage: mpirun -np 1 rosenbrock1 [-help] [all TAO options] */

      5: /*  Include "tao.h" so we can use TAO solvers.  */

     7:  #include tao.h

      9: static  char help[] = "This example demonstrates use of the TAO package to \n\
     10: solve an unconstrained minimization problem on a single processor.  We \n\
     11: minimize the extended Rosenbrock function: \n\
     12:    sum_{i=0}^{n/2-1} ( alpha*(x_{2i+1}-x_{2i}^2)^2 + (1-x_{2i})^2 ) \n";

     14: /*T
     15:    Concepts: TAO - Solving an unconstrained minimization problem
     16:    Routines: TaoInitialize(); TaoFinalize();
     17:    Routines: TaoApplicationCreate(); TaoAppDestroy();
     18:    Routines: TaoCreate(); TaoDestroy();
     19:    Routines: TaoAppSetObjectiveAndGradientRoutine();
     20:    Routines: TaoAppSetHessianMat(); TaoAppSetHessianRoutine();
     21:    Routines: TaoSetOptions();
     22:    Routines: TaoAppSetInitialSolutionVec();
     23:    Routines: TaoSolveApplication();
     24:    Routines: TaoGetTerminationReason();
     25:    Processors: 1
     26: T*/

     29: /*
     30:    User-defined application context - contains data needed by the
     31:    application-provided call-back routines that evaluate the function,
     32:    gradient, and hessian.
     33: */
     34: typedef struct {
     35:   int n;          /* dimension */
     36:   PetscReal alpha;   /* condition parameter */
     37: } AppCtx;

     39: /* -------------- User-defined routines ---------- */
     40: int FormFunctionGradient(TAO_APPLICATION,Vec,double*,Vec,void*);
     41: int FormHessian(TAO_APPLICATION,Vec,Mat*,Mat*,MatStructure*,void*);

     45: int main(int argc,char **argv)
     46: {
     47:   int        info;                  /* used to check for functions returning nonzeros */
     48:   PetscScalar zero=0.0;
     49:   Vec        x;                     /* solution vector */
     50:   Mat        H;                     /* Hessian matrix */
     51:   TAO_SOLVER tao;                   /* TAO_SOLVER solver context */
     52:   TAO_APPLICATION taoapp;           /* TAO application context */
     53:   PetscTruth  flg;
     54:   int        size,rank;                  /* number of processes running */
     55:   TaoTerminateReason reason;
     56:   AppCtx     user;                  /* user-defined application context */

     58:   /* Initialize TAO and PETSc */
     59:   PetscInitialize(&argc,&argv,(char *)0,help);
     60:   TaoInitialize(&argc,&argv,(char *)0,help);

     62:   info = MPI_Comm_size(PETSC_COMM_WORLD,&size); CHKERRQ(info);
     63:   info = MPI_Comm_rank(PETSC_COMM_WORLD,&rank); CHKERRQ(info);

     65:   if (size >1) {
     66:     if (rank == 0)
     67:       PetscPrintf(PETSC_COMM_SELF,"This example is intended for single processor use!\n");
     68:     SETERRQ(1,"Incorrect number of processors");
     69:   }

     72:   /* Initialize problem parameters */
     73:   user.n = 2; user.alpha = 99.0;

     75:   /* Check for command line arguments to override defaults */
     76:   info = PetscOptionsGetInt(PETSC_NULL,"-n",&user.n,&flg); CHKERRQ(info);
     77:   info = PetscOptionsGetReal(PETSC_NULL,"-alpha",&user.alpha,&flg); CHKERRQ(info);

     79:   /* Allocate vectors for the solution and gradient */
     80:   info = VecCreateSeq(PETSC_COMM_SELF,user.n,&x); CHKERRQ(info);

     82:   /*
     83:      Allocate storage space for Hessian matrix;
     84:      Hessian information is optional -- unless a Newton method is selected
     85:   */
     86:   info = MatCreateSeqBDiag(PETSC_COMM_SELF,user.n,user.n,0,2,0,0,&H);CHKERRQ(info);
     87:   info = MatSetOption(H,MAT_SYMMETRIC); CHKERRQ(info);

     89:   /* The TAO code begins here */

     91:   /* Create TAO solver with desired solution method */
     92:   info = TaoCreate(PETSC_COMM_SELF,"tao_lmvm",&tao); CHKERRQ(info);
     93:   info = TaoApplicationCreate(PETSC_COMM_SELF,&taoapp); CHKERRQ(info);

     95:   /* Set solution vec and an initial guess */
     96:   info = VecSet(x, zero); CHKERRQ(info);
     97:   info = TaoAppSetInitialSolutionVec(taoapp,x); CHKERRQ(info);

     99:   /* Set routines for function, gradient, hessian evaluation */
    100:   info = TaoAppSetObjectiveAndGradientRoutine(taoapp,FormFunctionGradient,(void *)&user);
    101:   CHKERRQ(info);
    102:   info = TaoAppSetHessianMat(taoapp,H,H); CHKERRQ(info);
    103:   info = TaoAppSetHessianRoutine(taoapp,FormHessian,(void *)&user); CHKERRQ(info);

    105:   /* Check for TAO command line options */
    106:   info = TaoSetOptions(taoapp,tao); CHKERRQ(info);

    108:   /* SOLVE THE APPLICATION */
    109:   info = TaoSolveApplication(taoapp,tao); CHKERRQ(info);

    111:   /* Get termination information */
    112:   info = TaoGetTerminationReason(tao,&reason); CHKERRQ(info);
    113:   if (reason <= 0)
    114:     PetscPrintf(MPI_COMM_WORLD,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n");

    117:   /* Free TAO data structures */
    118:   info = TaoDestroy(tao); CHKERRQ(info);
    119:   info = TaoAppDestroy(taoapp); CHKERRQ(info);

    121:   /* Free PETSc data structures */
    122:   info = VecDestroy(x); CHKERRQ(info);
    123:   info = MatDestroy(H); CHKERRQ(info);

    125:   /* Finalize TAO */
    126:   TaoFinalize();
    127:   PetscFinalize();

    129:   return 0;
    130: }

    132: /* -------------------------------------------------------------------- */
    135: /*
    136:     FormFunctionGradient - Evaluates the function, f(X), and gradient, G(X).

    138:     Input Parameters:
    139: .   taoapp  - the TAO_APPLICATION context
    140: .   X    - input vector
    141: .   ptr  - optional user-defined context, as set by TaoSetFunctionGradient()
    142:
    143:     Output Parameters:
    144: .   G - vector containing the newly evaluated gradient
    145: .   f - function value

    147:     Note:
    148:     Some optimization methods ask for the function and the gradient evaluation
    149:     at the same time.  Evaluating both at once may be more efficient that
    150:     evaluating each separately.
    151: */
    152: int FormFunctionGradient(TAO_APPLICATION taoapp,Vec X,double *f, Vec G,void *ptr)
    153: {
    154:   AppCtx *user = (AppCtx *) ptr;
    155:   int    i,info,nn=user->n/2;
    156:   double ff=0,t1,t2,alpha=user->alpha;
    157:   PetscScalar *x,*g;

    159:   /* Get pointers to vector data */
    160:   info = VecGetArray(X,&x); CHKERRQ(info);
    161:   info = VecGetArray(G,&g); CHKERRQ(info);

    163:   /* Compute G(X) */
    164:   for (i=0; i<nn; i++){
    165:     t1 = x[2*i+1]-x[2*i]*x[2*i]; t2= 1-x[2*i];
    166:     ff += alpha*t1*t1 + t2*t2;
    167:     g[2*i] = -4*alpha*t1*x[2*i]-2.0*t2;
    168:     g[2*i+1] = 2*alpha*t1;
    169:   }

    171:   /* Restore vectors */
    172:   info = VecRestoreArray(X,&x); CHKERRQ(info);
    173:   info = VecRestoreArray(G,&g); CHKERRQ(info);
    174:   *f=ff;

    176:   info = PetscLogFlops(nn*15); CHKERRQ(info);
    177:   return 0;
    178: }

    180: /* ------------------------------------------------------------------- */
    183: /*
    184:    FormHessian - Evaluates Hessian matrix.

    186:    Input Parameters:
    187: .  taoapp   - the TAO_APPLICATION context
    188: .  x     - input vector
    189: .  ptr   - optional user-defined context, as set by TaoSetHessian()

    191:    Output Parameters:
    192: .  H     - Hessian matrix

    194:    Note:  Providing the Hessian may not be necessary.  Only some solvers
    195:    require this matrix.
    196: */
    197: int FormHessian(TAO_APPLICATION taoapp,Vec X,Mat *HH, Mat *Hpre, MatStructure *flag,void *ptr)
    198: {
    199:   AppCtx  *user = (AppCtx*)ptr;
    200:   int     i, nn=user->n/2, info, ind[2];
    201:   double  alpha=user->alpha;
    202:   PetscScalar  v[2][2],*x;
    203:   Mat H=*HH;
    204:   PetscTruth assembled;

    206:   /* Zero existing matrix entries */
    207:   info = MatAssembled(H,&assembled); CHKERRQ(info);
    208:   if (assembled){info = MatZeroEntries(H);  CHKERRQ(info);}

    211:   /* Get a pointer to vector data */
    212:   info = VecGetArray(X,&x); CHKERRQ(info);

    214:   /* Compute H(X) entries */
    215:   for (i=0; i<user->n/2; i++){
    216:     v[1][1] = 2*alpha;
    217:     v[0][0] = -4*alpha*(x[2*i+1]-3*x[2*i]*x[2*i]) + 2;
    218:     v[1][0] = v[0][1] = -4.0*alpha*x[2*i];
    219:     ind[0]=2*i; ind[1]=2*i+1;
    220:     info = MatSetValues(H,2,ind,2,ind,v[0],INSERT_VALUES); CHKERRQ(info);
    221:   }
    222:   info = VecRestoreArray(X,&x); CHKERRQ(info);

    224:   /* Assemble matrix */
    225:   info = MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
    226:   info = MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
    227:   *flag=SAME_NONZERO_PATTERN;
    228:
    229:   info = PetscLogFlops(nn*9); CHKERRQ(info);
    230:   return 0;
    231: }