Commit a9986f22 by Francois Gygi

Added Polak-Ribiere option. simplified dot product calculations using BLAS.


git-svn-id: http://qboxcode.org/svn/qb/trunk@1006 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 11e562ef
......@@ -26,7 +26,7 @@ void CGOptimizer::compute_xp(const valarray<double>& x, const double f,
valarray<double>& g, valarray<double>& xp)
{
// Use the function value f and the gradient g at x to generate a new point xp
// using the Fletcher-Powell CG algorithm
// using the Fletcher-Reeves or Polak-Ribiere CG algorithm
// return xp=x if the 2-norm of g is smaller than tol
const double tol = 1.0e-18;
const int one = 1;
......@@ -38,6 +38,7 @@ void CGOptimizer::compute_xp(const valarray<double>& x, const double f,
if ( first_step_ )
{
p_ = -g;
gm_ = g;
x0_ = x;
f0_ = f;
......@@ -77,6 +78,7 @@ void CGOptimizer::compute_xp(const valarray<double>& x, const double f,
// restart from current point
p_ = -g;
gm_ = g;
x0_ = x;
f0_ = f;
......@@ -108,15 +110,21 @@ void CGOptimizer::compute_xp(const valarray<double>& x, const double f,
// define a new descent direction p_ using the Fletcher-Reeves formula
assert(g0norm2_ > 0.0);
double beta = ddot_(&n_,&g[0],&one,&g[0],&one) / g0norm2_;
#if 0
// Fletcher-Reeves
double beta = ddot_(&n_,&g[0],&one,&g[0],&one) / g0norm2_;
#else
// Polak-Ribiere
double beta = (ddot_(&n_,&g[0],&one,&g[0],&one)-
ddot_(&n_,&gm_[0],&one,&g[0],&one)) / g0norm2_;
#endif
if ( beta_max_ > 0.0 && beta > beta_max_ )
{
if ( debug_print )
cout << " CGOptimizer: beta exceeds beta_max " << endl;
beta = min(beta, beta_max_);
if ( beta > beta_max_ ) beta = 0.0;
}
#endif
if ( debug_print )
cout << " CGOptimizer: beta = " << beta << endl;
p_ = beta * p_ - g;
......@@ -127,6 +135,8 @@ void CGOptimizer::compute_xp(const valarray<double>& x, const double f,
// fp0 = d_e / d_alpha in direction pc_
fp = ddot_(&n_,&g[0],&one,&p_[0],&one);
g0norm2_ = ddot_(&n_,&g[0],&one,&g[0],&one);
gm_ = g;
fp0_ = fp;
if ( fp > 0.0 )
......@@ -136,6 +146,7 @@ void CGOptimizer::compute_xp(const valarray<double>& x, const double f,
if ( debug_print )
cout << " CGOptimizer: p_ not a descent direction" << endl;
p_ = -g;
gm_ = g;
x0_ = x;
f0_ = f;
......
......@@ -28,7 +28,7 @@ class CGOptimizer
int n_;
bool first_step_, debug_print;
std::valarray<double> x0_, p_;
std::valarray<double> x0_, p_, gm_;
double f0_, fp0_, g0norm2_, alpha_, beta_max_;
LineMinimizer linmin_;
double norm2(std::valarray<double>& v);
......@@ -40,6 +40,7 @@ class CGOptimizer
{
x0_.resize(n);
p_.resize(n);
gm_.resize(n);
}
void reset(void) { first_step_ = true; }
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment