/*$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;
}