C:/programs/UNCMIN/src/include/Uncmin.h

Go to the documentation of this file.

/*! \file Uncmin.h

  \brief
  C++ implementation of UNCMIN routines for unconstrained minimization
  
  This software is available at 
  http://www.smallwaters.com/software/cpp/uncmin.html
 
  Version 0080309
  This release updates original work by Brad Hanson (http://www.b-a-h.com/) in 2001.
 
  The UNCMIN routines are based on Appendix A of:
 
  Dennis, J. E., & Schnabel, R. B. (1996). Numerical methods for
  unconstrained optimization and nonlinear equations. Philadelphia:
  Society for Industrual and Applied Mathematics.
 
  The original UNCMIN routines in FORTRAN are available from GAMS
  (http://gams.nist.gov/) as OPTIF9 and OPTIF0.
 
  FORTRAN UNCMIN routines were translated to java (Uncmin_f77.java) by
  Steve Verrill (steve@ws13.fpl.fs.fed.us).
  Uncmin_F77.java is available at http://www1.fpl.fs.fed.us/optimization.html

  The routines in Uncmin_f77.java were translated to C++ by
  Brad Hanson.
  
  (c) 2008, Werner Wothke, maintenance programming and documentation.
 
 
  Uncmin is a templated class with template parameters:
 
  V - A vector class of double precision numbers. Member functions
      that should be defined for class V are:
 
      V& operator=(const V &A); // assignment
 
      V& operator=(const double& scalar); // assign a value to all elements
 
      double& operator()(int i); // FORTRAN style one-offset subscripting
 
      int size(); // returns number of elements in vector
 
  M - A matrix class of double precision numbers. Member functions that
      should be defined for class M are:
 
      double operator()(int i, int j); // FORTRAN style one-offset subscripting
 
      int num_rows(); // returns number of rows in matrix
      int num_cols(); // returns number of columns in matrix

  FUNC - A class defining the function to minimize, its gradient and hessian.
         An example of how this class should be defined is:
  
    // V is vector class, M is matrix class
    template <class V, class M>
    class minclass
    {
       double f_to_minimize(V &x);
       // Returns value of function to minimize
 
       void gradient(V &x, V &g);
       // Returns gradient of function in g. This does not have
       // to be defined, in which case the gradient is approximated.
       // If the gradient is not defined Analytic_Gradient should 
       // return 0;
 
       void hessian(V &x, M &h);
       // Returns hessian of function in h. This does not have
       // to be defined, in which case the hessian is approximated.
       // If the hessian is not defined Analytic_Hessian should 
       // return 0;
 
       int HasAnalyticGradient();
       // Returns 1 if gradient is defined, otherwise returns 0
 
       int HasAnalyticHessian();
       // Returns 1 if hessian is defined, otherwise returns 0
 
       int ValidParameters(V &x);
       // Returns 1 if f_to_minimize(x) is defined for argument x.
 
       int dim();
       // Returns dimension of problem (number of elements in vector x
       // passed to f_to_minimize)
    }
 
    An example of how to use Uncmin to minimize a function defined in
    a class minclass with vector and matrix classes from the Simple
    C++ Numerical Toolkit (http://www.smallwaters.com/software/cpp/scppnt.html):
 
    #include "Uncmin.h"
    #include "vec.h"  // Simple C++ Numeric Toolkit (SCPPNT) vector class
    #include "cmat.h" // Simple C++ Numeric Toolkit (SCPPNT) matrix class
    #include <cstdio>
 
    using namespace SCPPNT;
 
    // Dimension of problem (number of parameters)
    const int dim = 4;
 
    // Create function object
    minclass<Vector<double>, Matrix<double> > test;
 
    // create Uncmin object
    Uncmin<Vector<double>, Matrix<double>, minclass> min(&test);
 
    // Starting values
    double start_values[4] = {1.0, 1.0, 1.0, 1.0}; // use some appropriate values here
    Vector<double> start(start_values, start_values+dim);

    // xpls will contain solution, gpls will contain
    // gradient at solution after call to min.Minimize
    Vector<double> xpls(dim), gpls(dim);

    // fpls contains the value of the function at
    // the solution given by xpls.
    double fpls;
 
    // Minimize function.
    // Minimize returns zero if successful (0, 1, 2, 3
    // returned by GetMessage), non-zero if not successful.
    if (min.Minimize(start, xpls,  fpls,  gpls))
    {
       int msg = min.GetMessage(); // find out what went wrong
       std::printf("\nMessage returned from Uncmin: %d\n", msg);
    }

    The GetMessage member function returns the status of the solution found
    by the last call to the Minimize member function.
    Possible values returned by GetMessage are:
 
    0 = Optimal solution found, terminated with gradient small
    1 = Terminated with gradient small, xpls is probably optimal.
    2 = Terminated with stepsize small, xpls is probably optimal.
    3 = Lower point cannot be found, xpls is probably optimal.
    4 = Iteration limit (default 150) exceeded.
    5 = Too many large steps, function may be unbounded.
    -1 = Analytic gradient check requested, but no analytic gradient supplied
    -2 = Analytic hessian check requested, but no analytic hessian supplied
    -3 = Illegal dimension
    -4 = Illegal tolerance
    -5 = Illegal iteration limit
    -6 = Minimization function has no good digits
    -7 = Iteration limit exceeded in lnsrch
    -20 = Function not defined at starting value
    -21 = Check of analytic gradient failed
    -22 = Check of analytic hessian failed


    The symbol BOOST_NO_STDC_NAMESPACE should be defined if the C++ compiler
    does not put the standard C library functions in the std namespace
    (this is needed for Visual C++ 6).
 
    The symbol BOOST_NO_LIMITS should be defined if the standard C++
    library header file <limits> is not available.
    (this is needed for Visual C++ 6).

 
    Version History
 
    Version 080309 (March 9, 2008)
    
    Edited doxygen-compatible documentation, second pass. Added descriptive comments,
    much as possible, incorporated any useful information from the existing FORTRAN 
    and Java documentation sets and removed these out-dated documentations afterwards. 
    
    Version 071226 (December 26, 2007)
 
    Refactored num_cols property to num_columns in UNCMIN::minimize().
    Tagged all public methods and properties with doxygen-compatible comments.
 
    Version 0.001206 (February 10, 2001)
 
    Added BOOST_NO_STDC_NAMESPACE and BOOST_NO_LIMITS symbols which
    need to be defined to compile Uncmin.h using Microsoft Visual C++ 6.

 */

#ifndef UNCMIN_H_
#define UNCMIN_H_

#include <cstdio>
#include <cmath>

#ifdef BOOST_NO_STDC_NAMESPACE
// If the standard C library functions are not in the std namespace (like Microsoft Visual C++ 6)
namespace std
{ using ::sqrt; using ::pow; using ::fabs; using ::log; using ::FILE; using ::fprintf;}
#endif

#ifndef BOOST_NO_LIMITS
// include limits for epsilon function used in SetTolerances.
#include <limits>
#else
// If there is no limits header.
#include <float.h>
#endif

/*! Documentation of member functions is given in their definitions which follow
 the Uncmin class definition */
template<class V, class M, class FUNC> class Uncmin
{

public:

  typedef FUNC function_type;   //!< Type definition of function object to minimize.

  /* Constructor. */
  Uncmin(FUNC *f);
  
  /* Constructor. */
  Uncmin(int d = 1);

  /* Destructor */
  ~Uncmin();

  /* Calculate function minimum. */
  int Minimize(V &start, V &xpls, double &fpls, V &gpls, M &hpls);

  /* Version of Minimize where hessian matrix is allocated locally rather
   than passed as an argument to the function. */
  int Minimize(V &start, V &xpls, double &fpls, V &gpls);

  /* Get object containing function to minimize, gradient, and hessian. */
  FUNC *GetFunction();

  /* Set object containing function to minimize, gradient, and hessian. */
  void SetFunction(FUNC *f);

  /* Methods to set options. */
  void SetFuncExpensive(int iexp);

  /* Set typical magnitude of each argument to function. */
  int SetScaleArg(V typsiz);

  /* Set typical magnitude of function near minimum. */
  void SetScaleFunc(double fscale);

  /* Set the number of reliable digits returned by f_to_minimize. */
  void SetNumDigits(double ndigit);

  /* Set maximum number of iterations. */
  void SetMaxIter(int maxiter);

  /* Set step tolerance, gradient tolerance, and machine epsilon. */
  void SetTolerances(double step, double grad, double macheps = -1.0);

  /* Set maximum allowable scaled step length at any iteration. */
  void SetMaxStep(double stepmx);

  /* Set initial trust region radius. */
  void SetTrustRegionRadius(double dlt);

  /* Set flags to check analytic gradient and hessian. */
  void SetChecks(int check_gradient, int check_hessian = 0);

  /* Set flags to indicate whether results are printed. */
  void SetPrint(std::FILE *file, int print_results, int print_iterations = 0);

  /* Set method to use to solve minimization problem.
   1 = line search, 2 = double dogleg, 3 = More-Hebdon.  */
  void SetMethod(int method);

  /* Return stopping criteria computed at the end of the last call to Minimize. */
  void GetStoppingCriteria(double &gradtol, double &steptol, int &iterations);

  /*! Returns message generated in the last call to Minimize. */
  int GetMessage();

private:

  FUNC *minclass; //!<  Class that contains function to miminize, gradient, and hessian. 


  int algorithm; //!<  Indicates algorithm to use to solve minimization problem: 1 = line search, 2 = double dogleg, 3 = More-Hebdon.

  std::FILE *mfile; //!< File to print messages to.

  int n; //!< Dimension of problem, provided in minclass by minclass->dim().

  double epsm; //!< Machine epsilon.


  double mLastGradCrit;  //!< Stopping criterion for gradient.
  

  double mLastStepCrit; //!<  Stopping criterion for argument change.

  
  int mLastIteration; //!<  Stopping criterion by number of iterations.


  int mLastMessage; //!<  Message generated by last call to Minimize.

  /* Options */

  int iexp;  //!< Flag indicating computational complexity. Set iexp = 1, for expensive calculations, iexp = 0 otherwise. If iexp = 1, the Hessian will be evaluated by secant (BFGS) update.
  
  /* These two flags are set from minclass using minclass.Analytic_Gradient() and
   minclass.Analytic_Methods() */
  int iagflg; //!< Flag: iagflag = 0 if an analytic gradient is not supplied, provided by minclass.
  int iahflg; //!< Flag: iahflag = 0 if an analytic Hessian is not supplied, provided by minclass.

  int fCheckGradient; //!< fCheckGradient = 0 if analytic gradient is not checked.

  int fCheckHessian; //!< fCheckHessian = 0 if analytic hessian is not checked.

  int fPrintResults; //!< fPrintResults = 0 if results are not printed to mfile.

  int fPrintIterationResults; //!< fPrintIterationResults = 0 if results of each iteration are not printed to mfile.

  /* Tolerances */

  V typsiz; //!< Typical size of of each argument of function to minimize.

  double fscale; //!< Estimate of scale of objective function.

  int ndigit; //!< Number of good digits in the minimization function. Special case, set ndigit = -1, for function providing within 1 or 2 of full number of significant digits.

  int itnlim; //!< Maximum number of allowable iterations.

  double mdlt; //!< Trust region radius. Special case,  mdlt = -1, uses length of initial scaled gradient instead.

  double gradtl; //!< Tolerance at which the gradient is considered close enough to zero to terminate the algorithm.

  double steptl; //!< Tolerance at which scaled distance between two successive iterates is considered close enough to zero to terminate algorithm. 
  
  double stepmx; //!< Maximum allowable step size.

  /* Private functions used in minimization. */

  void dfault();

  int optchk(V &x, V &sx, int &msg, int &method);

  void fstofd(V &xpls, double fpls, V &g, V &sx, double rnoise);

  void fstofd(V &xpls, V &fpls, M &a, V &sx, double rnoise, V &fhat);

  int grdchk(V &x, double f, V &g, V &typsiz, V &sx, double fscale, double rnf, double analtl,
      V &gest, int &msg);

  void optstp(V &xpls, double fpls, V &gpls, V &x, int &itncnt, int &icscmx, int &itrmcd, V &sx,
      double &fscale, int &itnlim, int &iretcd, bool &mxtake, int &msg);

  void hsnint(M &a, V &sx, int method);

  void sndofd(V &xpls, double &fpls, M &a, V &sx, double &rnoise, V &stepsz, V &anbr);

  int heschk(V &x, double &f, V &g, M &a, V &typsiz, V &sx, double rnf, double analtl, int &iagflg,
      V &udiag, V &wrk1, V &wrk2, int &msg);

  void result(V &x, double &f, V &g, M &a, V &p, int &itncnt, int iflg);

  void bakslv(M &a, V &x, V &b);

  void chlhsn(M &a, V &sx, V &udiag);

  void choldc(M &a, double diagmx, double tol, double &addmax);

  void dogdrv(V &x, double &f, V &g, M &a, V &p, V &xpls, double &fpls, V &sx, double &dlt,
      int &iretcd, bool &mxtake, V &sc, V &wrk1, V &wrk2, V &wrk3);

  void dogstp(V &g, M &a, V &p, V &sx, double rnwtln, double &dlt, bool &nwtake, bool &fstdog,
      V &ssd, V &v, double &cln, double &eta, V &sc);

  void forslv(M &a, V &x, V &b);

  void fstocd(V &x, V &sx, double rnoise, V &g);

  void hookdr(V &x, double &f, V &g, M &a, V &udiag, V &p, V &xpls, double &fpls, V &sx,
      double &dlt, int &iretcd, bool &mxtake, double &amu, double &dltp, double &phi,
      double &phip0, V &sc, V &xplsp, V &wrk0, int &itncnt);

  void hookst(V &g, M &a, V &udiag, V &p, V &sx, double rnwtln, double &dlt, double &amu,
      double &dltp, double &phi, double &phip0, bool &fstime, V &sc, bool &nwtake, V &wrk0);

  void lltslv(M &a, V &x, V &b);

  void lnsrch(V &x, double &f, V &g, V &p, V &xpls, double &fpls, bool &mxtake, int &iretcd, V &sx);

  void mvmltl(M &a, V &x, V &y);

  void mvmlts(M &a, V &x, V &y);

  void mvmltu(M &a, V &x, V &y);

  void qraux1(M &r, int i);

  void qraux2(M &r, int i, double a, double b);

  void qrupdt(M &a, V &u, V &v);

  void secfac(V &x, V &g, M &a, V &xpls, V &gpls, int &itncnt, double rnf, int &iagflg,
      bool &noupdt, V &s, V &y, V &u, V &w);

  void secunf(V &x, V &g, M &a, V &udiag, V &xpls, V &gpls, int &itncnt, double rnf, int &iagflg,
      bool &noupdt, V &s, V &y, V &t);

  void tregup(V &x, double &f, V &g, M &a, V &sc, V &sx, bool &nwtake, double &dlt, int &iretcd,
      V &xplsp, double &fplsp, V &xpls, double &fpls, bool &mxtake, int method, V &udiag);

  double ddot(V &x, V &y);

  double dnrm2(V &x);

  double max(double a, double b);

  double min(double a, double b);

};

/*!
  \brief 
  Constructor that sets dimension and function to minimize.
  
  The minclass function object must define:
                    
  1. a method, f_to_minimize, to minimize. f_to_minimize must have the form 
         public static double f_to_minimize(double x[])
     where x is the vector of arguments to the function and the return value is the 
     value of the function evaluated at x.  

  2. a method, gradient, that has the form
         public static void gradient(double x[], double g[])
     where g is the gradient of f evaluated at x.  This method will have an empty body if 
     the user does not wish to provide an analytic estimate of the gradient.
     
  3. a method, hessian, that has the form
         public static void hessian(double x[],double h[][])
     where h is the Hessian of f evaluated at x.  This method will have an empty body if the 
     user does not wish to provide an analytic estimate of the Hessian. If the user wants 
     Uncmin++ to check the Hessian, then the hessian method should only fill the lower 
     triangle (and diagonal)of h.
  
   Minimization parameters are set at default values.
  
  \section template_args Template Parameters
  
  \param M  SCPPNT Matrix type.
  \param V  SCPPNT Vector type.
  \param FUNC Type of function to minimize.
  
  \section function_args Function Parameters
   
  \param[in]  *f  Pointer to function object f.
 */
template<class V, class M, class FUNC> Uncmin<V, M, FUNC>::Uncmin(FUNC *f) :
  minclass(f)
{
  if (minclass)
    n = minclass->dim();
  else
    n = 0;

  mfile = 0;

  mLastGradCrit = 0.0;
  mLastStepCrit = 0.0;
  mLastIteration = 0;
  mLastMessage = 0;

  /* Set default parameter values */
  dfault();
}

/*!
  \brief 
  Constructor that sets dimension, but not function to minimize.
  
  Note: Function to minimize must be set by calling SetFunction. 
  Minimization parameters are set to default values.
  
  \section template_args Template Parameters
  
  \param M  SCPPNT Matrix type.
  \param V  SCPPNT Vector type.
  \param FUNC Type of function to minimize.
  
  \section function_args Function Parameters
   
  \param[in]  d  Dimension d (number of function arguments).
 */
template<class V, class M, class FUNC> Uncmin<V, M, FUNC>::Uncmin(int d) :
  minclass(0), typsiz(d)
{

  n = d;

  mfile = 0;

  mLastGradCrit = 0.0;
  mLastStepCrit = 0.0;
  mLastIteration = 0;
  mLastMessage = 0;

  /* Set default parameter values */
  dfault();

}

/*!
  \brief Destructor of uncmin.
 */
template<class V, class M, class FUNC> Uncmin<V, M, FUNC>::~Uncmin()
{

}

/*!
  \brief 
  Assigns flag indicating whether function is expensive to evaluate.
  
  Set iexp =1 if function is expensive to evaluate, = 0 otherwise.  
  If iexp = 1, then the Hessian will be evaluated by secant (BFGS) update 
  rather than analytically or by finite differences.

  Default value is 0 if this function is not called.
  
  \section template_args Template Parameters
  
  \param M  SCPPNT Matrix type.
  \param V  SCPPNT Vector type.
  \param FUNC Type of function to minimize.
  
  \section function_args Function Parameters
   
  \param[in]  iexp  Set iexp =1 if function is expensive to evaluate, iexp = 0 otherwise.
 */
template<class V, class M, class FUNC> void Uncmin<V, M, FUNC>::SetFuncExpensive(int iexp)
{
  this->iexp = iexp;
}

/*!
  \brief
  Assigns typical magnitude of each argument to function.

  Returns zero if dimension matches, otherwise returns 1.
 
  The default is to set all elements of typsiz to 1 if this function is not called.
  
  \section template_args Template Parameters
  
  \param M  SCPPNT Matrix type.
  \param V  SCPPNT Vector type.
  \param FUNC Type of function to minimize.
  
  \section function_args Function Parameters
   
  \param[in]  typsiz  Vector describing the magnitude of all function arguments.
 */
template<class V, class M, class FUNC> int Uncmin<V, M, FUNC>::SetScaleArg(V typsiz)
{
  if (typsiz.size() != n)
    return 1;

  this->typsiz = typsiz;

  return 0;
}

/*!
  \brief
  Assigns typical magnitude of function near minimum.
  
  \section template_args Template Parameters
  
  \param M  SCPPNT Matrix type.
  \param V  SCPPNT Vector type.
  \param FUNC Type of function to minimize.
  
  \section function_args Function Parameters
   
  \param[in]  fscale  Magnitude of function values near minimum.
                      A default value of 1 is used for fscale if this
                      function is not called.
 */
template<class V, class M, class FUNC> void Uncmin<V, M, FUNC>::SetScaleFunc(double fscale)
{
  this->fscale = fscale;
}

/*!
  \brief
  Assigns the number of reliable digits returned by f_to_minimize.
  
  If the argument to the function is -1 this
  means that f_to_minimize is expected to
  provide within one or two of the full number of
  significant digits available.
 
  The default value of -1 is used for ndigit when
  this function is not called.

  \section template_args Template Parameters
  
  \param M  SCPPNT Matrix type.
  \param V  SCPPNT Vector type.
  \param FUNC Type of function to minimize.
  
  \section function_args Function Parameters
   
  \param[in]  ndigit  Number of digits precision required for the solution.
 */
template<class V, class M, class FUNC> void Uncmin<V, M, FUNC>::SetNumDigits(double ndigit)
{
  this->ndigit = ndigit;
}

/*!
  \brief
  Assigns the maximum number of iterations.
  
  The default value of 150 is used if this function is not called.

  \section template_args Template Parameters
  
  \param M  SCPPNT Matrix type.
  \param V  SCPPNT Vector type.
  \param FUNC Type of function to minimize.
  
  \section function_args Function Parameters
   
  \param[in]  maxiter  Maximum number of iterations.
 */
template<class V, class M, class FUNC> void Uncmin<V, M, FUNC>::SetMaxIter(int maxiter)
{
  this->itnlim = maxiter;
}

/*!
  \brief
  Assigns step tolerance, gradient tolerance, and the machine epsilon.
  
  Note: For any of these values <= zero, default values are used.
  
  \section template_args Template Parameters
  
  \param M  SCPPNT Matrix type.
  \param V  SCPPNT Vector type.
  \param FUNC Type of function to minimize.
  
  \section function_args Function Parameters
   
  \param[in]  step  Tolerance at which scaled distance between two successive iterates
    is considered close enough to zero to terminate algorithm. If the value of
    step is <= zero a default value is used.
  \param[in]  grad  Tolerance at which gradient is considered close enough
    to zero to terminate algorithm. If the value of
    grad is <= zero a default value is used.
  \param[in]  macheps Largest relative spacing. Machine epsilon == 2**(-T+1) where
    (1 + 2**-T) == 1. This argument has a default value of -1 if not included
    in function call.
 */
template<class V, class M, class FUNC> void Uncmin<V, M, FUNC>::SetTolerances(double step,
    double grad, double macheps)
{
  if (macheps <= 0.0) // use default value
  {

#ifdef BOOST_NO_LIMITS
    epsm = DBL_EPSILON; // from float.h
#else
    // use epsilon from numeric_limits in standard C++ library
    epsm = std::numeric_limits<double>::epsilon();
#endif

  }
  else
  {
    epsm = macheps;
  }

  if (grad <= 0.0)
  {
    gradtl = std::pow(epsm, 1.0/3.0); // default value
  }
  else
  {
    gradtl = grad;
  }

  if (step <= 0.0)
  {
    steptl = std::sqrt(epsm); // default value
  }
  else
  {
    steptl = step;
  }

}

/*!
  \brief
  Assigns maximum allowable scaled step length at any iteration.
  
  The maximum step length is used to prevent steps that would cause the 
  optimization algorithm to overflow or leave the domain of iterest,
  as well as to detect divergence. It should be chosen small enough
  to prevent the first two of these occurrences, but larger
  than any anticipated reasonable step size. The algorithm will
  halt if it takes steps larger than the maximum allowable step
  length on 5 consecutive iterations.
 
  If the argument to the function is <= 0 then a default value of
  stepmx will be computed in optchk.
 
  A default value of 0 is used for stepmx when this function is
  not called.

  \section template_args Template Parameters
  
  \param M  SCPPNT Matrix type.
  \param V  SCPPNT Vector type.
  \param FUNC Type of function to minimize.
  
  \section function_args Function Parameters
   
  \param[in]  stepmx  Maximum (scaled) step size.
 */
template<class V, class M, class FUNC> void Uncmin<V, M, FUNC>::SetMaxStep(double stepmx)
{
  this->stepmx = stepmx;
}

/*!
  \brief
  Assigns initial trust region radius.
  
  Used when method = 2 or 3, ignored when method = 1. 
  The value should be what the user considers a reasonable 
  scaled step length for the first iteration, and should be 
  less than the maximum step length.
 
  If dlt = -1, then the length of the initial scaled gradient is used instead.
 
  A default value of -1 is used if this function is not called.
 
  \section template_args Template Parameters
  
  \param M  SCPPNT Matrix type.
  \param V  SCPPNT Vector type.
  \param FUNC Type of function to minimize.
  
  \section function_args Function Parameters
   
  \param[in]  dlt  Maximum (scaled) step size of first iteration.
 */
template<class V, class M, class FUNC> void Uncmin<V, M, FUNC>::SetTrustRegionRadius(double dlt)
{
  this->mdlt = dlt;
}

/*!
  \brief
  Returns pointer to object containing function to minimize, gradient, and hessian.
 
  \section template_args Template Parameters
  
  \param M  SCPPNT Matrix type.
  \param V  SCPPNT Vector type.
  \param FUNC Type of function to minimize.
 */
template<class V, class M, class FUNC> FUNC * Uncmin<V, M, FUNC>::GetFunction()
{
  return minclass;
}


/*!
  \brief
  Assigns pointer to object containing function to minimize, gradient, and hessian.

  \section template_args Template Parameters
  
  \param M  SCPPNT Matrix type.
  \param V  SCPPNT Vector type.
  \param FUNC Type of function to minimize.
  
  \section function_args Function Parameters
   
  \param[in]  *f  Pointer to function object f.
 */
template<class V, class M, class FUNC> void Uncmin<V, M, FUNC>::SetFunction(FUNC *f)
{
  minclass = f;

  /* If dimension does not match change n and re-intialize typsiz */
  if (n != f->dim())
  {
    n = f->dim();
    V defaultTypSiz(n, 1.0);
    typsiz = defaultTypSiz;
  }
}

/*!
  \brief
  Assigns flags to check analytic gradient and hessian.

  The default is to not check the gradient and hessian, if this function is not called
  (fCheckGradient = 0 and fCheckHessian = 0).
 
  Note that check_hessian can be left off function call in which case default value of
  0 is used.
  
  \section template_args Template Parameters
  
  \param M  SCPPNT Matrix type.
  \param V  SCPPNT Vector type.
  \param FUNC Type of function to minimize.
  
  \section function_args Function Parameters
   
  \param[in]  check_gradient  Flag to check analytic gradient of function: 1 = yes, 0 = no.
  \param[in]  check_hessian  Flag to check analytic hessian of function: 1 = yes, 0 = no.
 */
template<class V, class M, class FUNC> void Uncmin<V, M, FUNC>::SetChecks(int check_gradient,
    int check_hessian)
{
  fCheckGradient = check_gradient;
  fCheckHessian = check_hessian;
}

/*!
  \brief
  Assigns flags to indicate whether results are printed and specify which file is 
  receive the print output.

  Printing is done to open file indicated by first argument.

  The default action is not to print any output.
  
  \section template_args Template Parameters
  
  \param M  SCPPNT Matrix type.
  \param V  SCPPNT Vector type.
  \param FUNC Type of function to minimize.
  \param FILE Type of file for writing print output to.
  
  \section function_args Function Parameters
   
  \param[in]  *file  Pointer to the open file that will receive the printed results.
  \param[in]  print_results  Flag (if non-zero) to print standard results.
  \param[in]  print_iterations  Flag (if non-zero) to print intermediate results at each 
      iteration (this argument can be left off function call in which case the
      default value of 0 is used).
  */
template<class V, class M, class FUNC> void Uncmin<V, M, FUNC>::SetPrint(std::FILE *file,
    int print_results, int print_iterations)
{

  mfile = file;

  if (!file)
  {
    fPrintResults = 0;
    fPrintIterationResults = 0;
  }
  else
  {
    fPrintResults = print_results;
    fPrintIterationResults = print_iterations;
  }

}

/*!
  \brief
  Assigns method to use to solve minimization problem.
  
  A default value is 1 (line search) if SetMethod is not called.

  \section template_args Template Parameters
  
  \param M  SCPPNT Matrix type.
  \param V  SCPPNT Vector type.
  \param FUNC Type of function to minimize.
  
  \section function_args Function Parameters
   
  \param[in]  method  Method indicator: 1 = line search, 2 = double dogleg, 3 = More-Hebdon.
 */
template<class V, class M, class FUNC> void Uncmin<V, M, FUNC>::SetMethod(int method)
{
  if (method == 2 || method == 3)
    algorithm = method;
  else
    algorithm = 1; // use default value of 1 for any values of argument other than 2 or 3
}

/*!
  \brief
  Returns stopping criteria computed at the end of the last call to Minimize.

  \section template_args Template Parameters
  
  \param M  SCPPNT Matrix type.
  \param V  SCPPNT Vector type.
  \param FUNC Type of function to minimize.
  
  \section function_args Function Parameters
   
  \param[out]  &gradtol  Address of gradient tolerance, i.e., how close the gradient is to zero.
  \param[out]  &steptol  Address of maximum scaled difference between function parameters in 
      consecutive iterations.
  \param[out]  &iterations  Address of number of iterations used.
 */
template<class V, class M, class FUNC> void Uncmin<V, M, FUNC>::GetStoppingCriteria(
    double &gradtol, double &steptol, int &iterations)
{
  gradtol = mLastGradCrit;
  steptol = mLastStepCrit;
  iterations = mLastIteration;
}

/*!
  \brief
  Returns message index generated in the last call to Minimize.
 
  Possible values of message are:
 
    0 = Optimal solution found, terminated with gradient small.
    1 = Terminated with gradient small, xpls is probably optimal.
    2 = Terminated with stepsize small, xpls is probably optimal.
    3 = Lower point cannot be found, xpls is probably optimal.
    4 = Iteration limit (default 150) exceeded.
    5 = Too many large steps, function may be unbounded.
   -1 = Analytic gradient check requested, but no analytic gradient supplied
   -2 = Analytic hessian check requested, but no analytic hessian supplied
   -3 = Illegal dimension
   -4 = Illegal tolerance
   -5 = Illegal iteration limit
   -6 = Minimization function has no good digits
   -7 = Iteration limit exceeded in lnsrch
  -20 = Function not defined at starting value
  -21 = Check of analytic gradient failed
  -22 = Check of analytic hessian failed

  \section template_args Template Parameters
  
  \param M  SCPPNT Matrix type.
  \param V  SCPPNT Vector type.
  \param FUNC Type of function to minimize.
 */
template<class V, class M, class FUNC> inline int Uncmin<V, M, FUNC>::GetMessage()
{
  return mLastMessage;
}

/*!
  \brief
  Method to compute values that minimize the function.

  Translated by Brad Hanson from optdrv_f77 by Steve Verrill.
  
  Returns 0 if mLastMessage is 0, 1, 2, or 3 (indicating optimal solution probably found), 
  otherwise returns non-zero.
  
  On exit data member mLastMessage is set to indicate status of solution.

  Possible values of mLastMessage are:
    0 = Optimal solution found, terminated with gradient small.
    1 = Terminated with gradient small, xpls is probably optimal.
    2 = Terminated with stepsize small, xpls is probably optimal.
    3 = Lower point cannot be found, xpls is probably optimal.
    4 = Iteration limit (default 150) exceeded.
    5 = Too many large steps, function may be unbounded.
   -1 = Analytic gradient check requested, but no analytic gradient supplied.
   -2 = Analytic hessian check requested, but no analytic hessian supplied.
   -3 = Illegal dimension.
   -4 = Illegal tolerance.
   -5 = Illegal iteration limit.
   -6 = Minimization function has no good digits. 
   -7 = Iteration limit exceeded in lnsrch.
  -20 = Function not defined at starting value.
  -21 = Check of analytic gradient failed.
  -22 = Check of analytic hessian failed.

  \section template_args Template Parameters
  
  \param M  SCPPNT Matrix type.
  \param V  SCPPNT Vector type.
  \param FUNC Type of function to minimize.
  
  \section function_args Function Parameters
   
  \param[in]   &start  Address of vector containing initial estimate of location of minimum.
  \param[out]  &xpls  Address of vector containing local minimum (size must be n on input).
  \param[out]  &fpls  Address of function value at local minimum.
  \param[out]  &gpls  Address of vector containing gradient at local minimum (size must be 
      n on input). 
  \param[out]  &hpls  Address of matrix containing hessian at local minimum (size must be 
      n X n on input). Only the lower half should be input; the upper off-diagonal entries 
      can be set to zero. On return, the lower half of hpls contains the Cholesky factor 
      of the Hessian, evaluated at the solution (i.e., H = LL').
 */
template<class V, class M, class FUNC> int Uncmin<V, M, FUNC>::Minimize(V &start, V &xpls,
    double &fpls, V &gpls, M &hpls)
{

  bool noupdt;
  bool mxtake;

  int i, j, num5, remain, ilow, ihigh;
  int icscmx;
  int iretcd;
  int itncnt;

  /*
    Return code.
      itrmcd =  0:    Optimal solution found
      itrmcd =  1:    Terminated with gradient small, X is probably optimal
      itrmcd =  2:    Terminated with stepsize small, X is probably optimal
      itrmcd =  3:    Lower point cannot be found, X is probably optimal
      itrmcd =  4:    Iteration limit (150) exceeded
      itrmcd =  5:    Too many large steps, function may be unbounded.
   */
  int itrmcd;

  double rnf, analtl, dltsav, amusav, dlpsav, phisav, phpsav;

  double f;

  double amu = 0.0;
  double dltp = 0.0;
  double phi = 0.0;
  double phip0 = 0.0;

  int method = algorithm;

  mLastGradCrit = 0.0;
  mLastStepCrit = 0.0;
  mLastIteration = 0;
  mLastMessage = -9999;

  iagflg = minclass->HasAnalyticGradient();
  iahflg = minclass->HasAnalyticHessian();

  /*
    msg
    Flag to set various options.
    Bits which indicate flags which control checks of gradient and hessian and what output is printed.
    msg should be zero or a sum of the following values which indicate which flags to set or not set.

      1 = Turn off check of analytic gradient.
      2 = Turn off check of analytic hessian.
      4 = Turn off printing of results.
      8 = Turn on printing of results at each iteration.
   */
  int msg = (fCheckGradient == 0) * 1 + (fCheckHessian == 0) * 2 + (fPrintResults == 0) * 4
      + (fPrintIterationResults != 0) * 8;

  /* Check that dimensions match */
  if (start.size() != n || xpls.size() != n || gpls.size() != n || hpls.num_rows() != n
      || hpls.num_columns() != n) // Refactored function "num_columns", ww, 12-22-2007
  {
    mLastMessage = -3;
    return -1;
  }

  /* Check that starting values are valid */
  if (!(minclass->ValidParameters(start)))
  {
    mLastMessage = -20;
    return -1;
  }

  /* Workspace */
  M &a = hpls;
  V udiag(n); //workspace (for diagonal of Hessian)
  V g(n); // workspace (for gradient at current iterate)
  V p(n); // workspace for step
  V sx(n); // workspace (for scaling vector)
  V wrk0(n), wrk1(n), wrk2(n), wrk3(n);

  V x(start);

  double dlt = mdlt;

  dltsav = amusav = dlpsav = phisav = phpsav = 0.0;

  // INITIALIZATION

  for (i = 1; i <= n; i++)
  {
    p(i) = 0.0;
  }

  itncnt = 0;
  iretcd = -1;

  /* Check for valid options */
  if (optchk(x, sx, msg, method))
  {
    mLastMessage = msg;
    return -1;
  }

  rnf = max(std::pow(10.0, -ndigit), epsm);
  analtl = max(.01, std::sqrt(rnf));

  if (!(msg & 4))
  {
    num5 = n/5;
    remain = n%5;

    std::fprintf(mfile, "\n\nOPTDRV          Typical x\n\n");

    ilow = -4;
    ihigh = 0;

    for (i = 1; i <= num5; i++)
    {
      ilow += 5;
      ihigh += 5;

      std::fprintf(mfile, "%d--%d     ", ilow, ihigh);

      for (j = 1; j <= 5; j++)
      {
        std::fprintf(mfile, "%lf  ", (double) typsiz(ilow+j-1));
      }
      std::fprintf(mfile, "\n");
    }

    ilow += 5;
    ihigh = ilow + remain - 1;

    std::fprintf(mfile, "%d--%d     ", ilow, ihigh);

    for (j = 1; j <= remain; j++)
    {
      std::fprintf(mfile, "%lf  ", (double) typsiz(ilow+j-1));
    }

    std::fprintf(mfile, "\n");

    std::fprintf(mfile, "\n\nOPTDRV      Scaling vector for x\n\n");

    ilow = -4;
    ihigh = 0;

    for (i = 1; i <= num5; i++)
    {
      ilow += 5;
      ihigh += 5;

      std::fprintf(mfile, "%d--%d     ", ilow, ihigh);

      for (j = 1; j <= 5; j++)
      {
        std::fprintf(mfile, "%lf  ", (double) sx(ilow+j-1));
      }
      std::fprintf(mfile, "\n");
    }

    ilow += 5;
    ihigh = ilow + remain - 1;

    std::fprintf(mfile, "%d--%d     ", ilow, ihigh);

    for (j = 1; j <= remain; j++)
    {
      std::fprintf(mfile, "%lf  ", (double) sx(ilow+j-1));
    }

    std::fprintf(mfile, "\n");

    std::fprintf(mfile, "\n\nOPTDRV      Typical f = %f\n", fscale);
    std::fprintf(mfile, "OPTDRV      Number of good digits in");
    std::fprintf(mfile, " f_to_minimize = %d\n", ndigit);
    std::fprintf(mfile, "OPTDRV      Gradient flag");
    std::fprintf(mfile, " = %d\n", iagflg);
    std::fprintf(mfile, "OPTDRV      Hessian flag");
    std::fprintf(mfile, " = %d\n", iahflg);
    std::fprintf(mfile, "OPTDRV      Expensive function calculation flag");
    std::fprintf(mfile, " = %d\n", iexp);
    std::fprintf(mfile, "OPTDRV      Method to use");
    std::fprintf(mfile, " = %d\n", method);
    std::fprintf(mfile, "OPTDRV      Iteration limit");
    std::fprintf(mfile, " = %d\n", itnlim);
    std::fprintf(mfile, "OPTDRV      Machine epsilon");
    std::fprintf(mfile, " = %e\n", epsm);
    std::fprintf(mfile, "OPTDRV      Maximum step size");
    std::fprintf(mfile, " = %f\n", stepmx);
    std::fprintf(mfile, "OPTDRV      Step tolerance");
    std::fprintf(mfile, " = %e\n", steptl);
    std::fprintf(mfile, "OPTDRV      Gradient tolerance");
    std::fprintf(mfile, " = %e\n", gradtl);
    std::fprintf(mfile, "OPTDRV      Trust region radius");
    std::fprintf(mfile, " = %f\n", dlt);
    std::fprintf(mfile, "OPTDRV      Relative noise in");
    std::fprintf(mfile, " f_to_minimize = %e\n", rnf);
    std::fprintf(mfile, "OPTDRV      Analytical fd tolerance");
    std::fprintf(mfile, " = %e\n", analtl);
  }

  // EVALUATE FCN(X)
  f = minclass->f_to_minimize(x);

  // EVALUATE ANALYTIC OR FINITE DIFFERENCE GRADIENT AND CHECK ANALYTIC
  // GRADIENT, IF REQUESTED.
  if (iagflg == 0)
  {
    fstofd(x, f, g, sx, rnf);
  }
  else
  {
    minclass->gradient(x, g);

    if (!(msg & 1))
    {
      if (grdchk(x, f, g, typsiz, sx, fscale, rnf, analtl, wrk1, msg))
      {
        mLastMessage = msg;
        return -1; // msg = -21
      }
    }
  }

  optstp(x, f, g, wrk1, itncnt, icscmx, itrmcd, sx, fscale, itnlim, iretcd, mxtake, msg);

  if (itrmcd != 0)
  {
    fpls = f;

    for (i = 1; i <= n; i++)
    {
      xpls(i) = x(i);
      gpls(i) = g(i);
    }
  }
  else
  {
    if (iexp == 1)
    {
      // IF OPTIMIZATION FUNCTION EXPENSIVE TO EVALUATE (IEXP=1), THEN
      // HESSIAN WILL BE OBTAINED BY SECANT (BFGS) UPDATES.  GET INITIAL HESSIAN.

      hsnint(a, sx, method);
    }
    else
    {
      // EVALUATE ANALYTIC OR FINITE DIFFERENCE HESSIAN AND CHECK ANALYTIC
      // HESSIAN IF REQUESTED (ONLY IF USER-SUPPLIED ANALYTIC HESSIAN
      // ROUTINE minclass->hessian FILLS ONLY LOWER TRIANGULAR PART AND DIAGONAL OF A).

      if (iahflg == 0)
      {
        if (iagflg == 1)
        {
          fstofd(x, g, a, sx, rnf, wrk1);
        }
        else
        {
          sndofd(x, f, a, sx, rnf, wrk1, wrk2);
        }
      }
      else
      {
        if (msg & 2)
        {
          minclass->hessian(x, a);
        }
        else
        {
          if (heschk(x, f, g, a, typsiz, sx, rnf, analtl, iagflg, udiag, wrk1, wrk2, msg))
          {
            mLastMessage = msg;
            return -1; // msg = -22
          }
          // HESCHK EVALUATES minclass->hessian AND CHECKS IT AGAINST THE FINITE
          // DIFFERENCE HESSIAN WHICH IT CALCULATES BY CALLING FSTOFD
          // (IF IAGFLG .EQ. 1) OR SNDOFD (OTHERWISE).
        }
      }
    }

    if (!(msg & 4))
      result(x, f, g, a, p, itncnt, 1);

    // ITERATION
    while (itrmcd == 0)
    {
      itncnt++;
      mLastIteration = itncnt; // BAH

      // FIND PERTURBED LOCAL MODEL HESSIAN AND ITS LL+ DECOMPOSITION
      // (SKIP THIS STEP IF LINE SEARCH OR DOGSTEP TECHNIQUES BEING USED WITH
      // SECANT UPDATES.  CHOLESKY DECOMPOSITION L ALREADY OBTAINED FROM
      // SECFAC.)
      if (iexp != 1 || method == 3)
      {
        chlhsn(a, sx, udiag);
      }

      // SOLVE FOR NEWTON STEP:  AP=-G
      for (i = 1; i <= n; i++)
      {
        wrk1(i) = -g(i);
      }

      lltslv(a, p, wrk1);

      // DECIDE WHETHER TO ACCEPT NEWTON STEP  XPLS=X + P
      // OR TO CHOOSE XPLS BY A GLOBAL STRATEGY.
      if (iagflg == 0 && method != 1)
      {
        dltsav = dlt;

        if (method != 2)
        {
          amusav = amu;
          dlpsav = dltp;
          phisav = phi;
          phpsav = phip0;
        }
      }

      if (method == 1)
      {
        lnsrch(x, f, g, p, xpls, fpls, mxtake, iretcd, sx);
        if (iretcd == -7) // BAH
        {
          mLastMessage = iretcd;
          return -1;
        }
      }
      else if (method == 2)
      {
        dogdrv(x, f, g, a, p, xpls, fpls, sx, dlt, iretcd, mxtake, wrk0, wrk1, wrk2, wrk3);
      }
      else
      {
        hookdr(x, f, g, a, udiag, p, xpls, fpls, sx, dlt, iretcd, mxtake, amu, dltp, phi, phip0,
            wrk0, wrk1, wrk2, itncnt);
      }

      // IF COULD NOT FIND SATISFACTORY STEP AND FORWARD DIFFERENCE
      // GRADIENT WAS USED, RETRY USING CENTRAL DIFFERENCE GRADIENT.
      if (iretcd == 1 && iagflg == 0)
      {
        // SET IAGFLG FOR CENTRAL DIFFERENCES
        iagflg = -1;

        if (!(msg & 4))
        {
          std::fprintf(mfile, "\nOPTDRV      Shift from forward to central");
          std::fprintf(mfile, " differences");
          std::fprintf(mfile, "\nOPTDRV      in iteration %d", itncnt);
          std::fprintf(mfile, "\n");
        }

        fstocd(x, sx, rnf, g);

        if (method == 1)
        {
          // SOLVE FOR NEWTON STEP:  AP=-G
          for (i = 1; i <= n; i++)
          {
            wrk1(i) = -g(i);
          }

          lltslv(a, p, wrk1);

          lnsrch(x, f, g, p, xpls, fpls, mxtake, iretcd, sx);
          if (iretcd == -7) // BAH
          {
            mLastMessage = iretcd;
            return -1;
          }
        }
        else
        {
          dlt = dltsav;
          if (method == 2)
          {
            // SOLVE FOR NEWTON STEP:  AP=-G
            for (i = 1; i <= n; i++)
            {
              wrk1(i) = -g(i);
            }

            lltslv(a, p, wrk1);

            dogdrv(x, f, g, a, p, xpls, fpls, sx, dlt, iretcd, mxtake, wrk0, wrk1, wrk2, wrk3);
          }
          else
          {
            amu = amusav;
            dltp = dlpsav;
            phi = phisav;
            phip0 = phpsav;

            chlhsn(a, sx, udiag);

            // SOLVE FOR NEWTON STEP:  AP=-G
            for (i = 1; i <= n; i++)
            {
              wrk1(i) = -g(i);
            }

            lltslv(a, p, wrk1);

            hookdr(x, f, g, a, udiag, p, xpls, fpls, sx, dlt, iretcd, mxtake, amu, dltp, phi,
                phip0, wrk0, wrk1, wrk2, itncnt);
          }
        }
      }
      // CALCULATE STEP FOR OUTPUT
      for (i = 1; i <= n; i++)
      {
        p(i) = xpls(i) - x(i);
      }

      // CALCULATE GRADIENT AT XPLS
      if (iagflg == -1)
      {
        // CENTRAL DIFFERENCE GRADIENT
        fstocd(xpls, sx, rnf, gpls);
      }
      else if (iagflg == 0)
      {
        // FORWARD DIFFERENCE GRADIENT
        fstofd(xpls, fpls, gpls, sx, rnf);
      }
      else
      {
        // ANALYTIC GRADIENT
        minclass->gradient(xpls, gpls);
      }

      // CHECK WHETHER STOPPING CRITERIA SATISFIED
      optstp(xpls, fpls, gpls, x, itncnt, icscmx, itrmcd, sx, fscale, itnlim, iretcd, mxtake, msg);

      if (itrmcd == 0)
      {
        // EVALUATE HESSIAN AT XPLS
        if (iexp != 0)
        {
          if (method == 3)
          {
            secunf(x, g, a, udiag, xpls, gpls, itncnt, rnf, iagflg, noupdt, wrk1, wrk2, wrk3);
          }
          else
          {
            secfac(x, g, a, xpls, gpls, itncnt, rnf, iagflg, noupdt, wrk0, wrk1, wrk2, wrk3);
          }
        }
        else
        {
          if (iahflg == 1)
          {
            minclass->hessian(xpls, a);
          }
          else
          {
            if (iagflg == 1)
            {
              fstofd(xpls, gpls, a, sx, rnf, wrk1);
            }
            else
            {
              sndofd(xpls, fpls, a, sx, rnf, wrk1, wrk2);
            }
          }
        }

        if