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

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

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

    #include "tao.h"

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

    /*T
       Concepts: TAO - Solving an unconstrained minimization problem
       Routines: TaoInitialize(); TaoFinalize();
       Routines: TaoApplicationCreate(); TaoAppDestroy();
       Routines: TaoCreate(); TaoDestroy();
       Routines: TaoAppSetObjectiveAndGradientRoutine();
       Routines: TaoAppSetHessianMat(); TaoAppSetHessianRoutine();
       Routines: TaoSetOptions();
       Routines: TaoAppSetInitialSolutionVec();
       Routines: TaoSolveApplication();
       Routines: TaoGetTerminationReason();
       Processors: 1
    T*/

    /*
       User-defined application context - contains data needed by the
       application-provided call-back routines that evaluate the function,
       gradient, and hessian.
    */
    typedef struct {
      int n;          /* dimension */
      PetscReal alpha;   /* condition parameter */
    } AppCtx;

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

    #undef __FUNCT__
    #define __FUNCT__ "main"
    int main(int argc,char **argv)
    {
      int        info;                  /* used to check for functions returning nonzeros */
      PetscScalar zero=0.0;
      Vec        x;                     /* solution vector */
      Mat        H;                     /* Hessian matrix */
      TAO_SOLVER tao;                   /* TAO_SOLVER solver context */
      TAO_APPLICATION taoapp;           /* TAO application context */
      PetscTruth  flg;
      int        size,rank;                  /* number of processes running */
      TaoTerminateReason reason;
      AppCtx     user;                  /* user-defined application context */

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

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

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

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

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

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

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

      /* The TAO code begins here */

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

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

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

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

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

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

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

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

      /* Finalize TAO */
      TaoFinalize();
      PetscFinalize();

      return 0;
    }

    /* -------------------------------------------------------------------- */
    #undef __FUNCT__
    #define __FUNCT__ "FormFunctionGradient"
    /*
        FormFunctionGradient - Evaluates the function, f(X), and gradient, G(X).

        Input Parameters:
    .   taoapp  - the TAO_APPLICATION context
    .   X    - input vector
    .   ptr  - optional user-defined context, as set by TaoSetFunctionGradient()

        Output Parameters:
    .   G - vector containing the newly evaluated gradient
    .   f - function value

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

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

      /* Compute G(X) */
      for (i=0; i<nn; i++){
        t1 = x[2*i+1]-x[2*i]*x[2*i]; t2= 1-x[2*i];
        ff += alpha*t1*t1 + t2*t2;
        g[2*i] = -4*alpha*t1*x[2*i]-2.0*t2;
        g[2*i+1] = 2*alpha*t1;
      }

      /* Restore vectors */
      info = VecRestoreArray(X,&x); CHKERRQ(info);
      info = VecRestoreArray(G,&g); CHKERRQ(info);
      *f=ff;

      info = PetscLogFlops(nn*15); CHKERRQ(info);
      return 0;
    }

    /* ------------------------------------------------------------------- */
    #undef __FUNCT__
    #define __FUNCT__ "FormHessian"
    /*
       FormHessian - Evaluates Hessian matrix.

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

       Output Parameters:
    .  H     - Hessian matrix

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

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

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

      /* Compute H(X) entries */
      for (i=0; i<user->n/2; i++){
        v[1][1] = 2*alpha;
        v[0][0] = -4*alpha*(x[2*i+1]-3*x[2*i]*x[2*i]) + 2;
        v[1][0] = v[0][1] = -4.0*alpha*x[2*i];
        ind[0]=2*i; ind[1]=2*i+1;
        info = MatSetValues(H,2,ind,2,ind,v[0],INSERT_VALUES); CHKERRQ(info);
      }
      info = VecRestoreArray(X,&x); CHKERRQ(info);

      /* Assemble matrix */
      info = MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
      info = MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
      *flag=SAME_NONZERO_PATTERN;

      info = PetscLogFlops(nn*9); CHKERRQ(info);
      return 0;
    }