Commit 24f7991c by Francois Gygi

removed use of inner_product

git-svn-id: http://qboxcode.org/svn/qb/trunk@966 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 7e897a56
......@@ -18,14 +18,18 @@
#include "CGOptimizer.h"
#include <iostream>
#include <cassert>
#include <numeric>
using namespace std;
////////////////////////////////////////////////////////////////////////////////
double CGOptimizer::norm2(const std::valarray<double>& v)
double norm2(const std::valarray<double>&v)
{
return inner_product(&v[0],&v[v.size()],&v[0],0.0);
double sum = 0.0;
for ( int i = 0; i < v.size(); i++ )
sum += v[i]*v[i];
return sum;
}
////////////////////////////////////////////////////////////////////////////////
void CGOptimizer::compute_xp(const valarray<double>& x, const double f,
void CGOptimizer::compute_xp(const valarray<double>& x, const double f,
const valarray<double>& g, valarray<double>& xp)
{
// Use the function value f and the gradient g at x to generate a new point xp
......@@ -55,9 +59,9 @@ void CGOptimizer::compute_xp(const valarray<double>& x, const double f,
linmin_.reset();
alpha_ = linmin_.next_alpha(alpha_,f,fp);
if ( debug_print )
cout << " CGOptimizer: first_step: alpha=" << alpha_
cout << " CGOptimizer: first_step: alpha=" << alpha_
<< " f=" << f << " fp=" << fp << endl;
xp = x0_ + alpha_ * p_;
first_step_ = false;
}
......@@ -65,10 +69,12 @@ void CGOptimizer::compute_xp(const valarray<double>& x, const double f,
{
// fp: derivative along the current descent direction p_
// fp = df(x0+alpha*p)/dalpha at x
fp = inner_product(&g[0],&g[g.size()],&p_[0],0.0);
fp = 0.0;
for ( int i = 0; i < g.size(); i++ )
fp += g[i] * p_[i];
double new_alpha = linmin_.next_alpha(alpha_,f,fp);
if ( debug_print )
cout << " CGOptimizer: alpha=" << alpha_
cout << " CGOptimizer: alpha=" << alpha_
<< " f=" << f << " fp=" << fp << endl;
if ( linmin_.fail() )
......@@ -94,9 +100,9 @@ void CGOptimizer::compute_xp(const valarray<double>& x, const double f,
linmin_.reset();
alpha_ = linmin_.next_alpha(alpha_,f,fp);
if ( debug_print )
cout << " CGOptimizer: restart after fail: alpha=" << alpha_
cout << " CGOptimizer: restart after fail: alpha=" << alpha_
<< " f=" << f << " fp=" << fp << endl;
xp = x0_ + alpha_ * p_;
first_step_ = false;
}
......@@ -123,7 +129,9 @@ void CGOptimizer::compute_xp(const valarray<double>& x, const double f,
f0_ = f;
// recalculate f0, fp0
// fp0 = d_e / d_alpha in direction pc_
fp0_ = inner_product(&g[0],&g[g.size()],&p_[0],0.0);
fp0_ = 0.0;
for ( int i = 0; i < g.size(); i++ )
fp0_ += g[i] * p_[i];
g0norm2_ = norm2(g);
fp = fp0_;
assert(fp<0.0 && "CGOptimizer: p_ not a descent direction");
......@@ -132,7 +140,7 @@ void CGOptimizer::compute_xp(const valarray<double>& x, const double f,
linmin_.reset();
alpha_ = linmin_.next_alpha(alpha_,f,fp);
if ( debug_print )
cout << " CGOptimizer: restart: alpha=" << alpha_
cout << " CGOptimizer: restart: alpha=" << alpha_
<< " f=" << f << " fp=" << fp << endl;
}
else
......
......@@ -28,17 +28,16 @@ class CGOptimizer
int n_;
bool first_step_, debug_print;
std::valarray<double> x0_, p_;
std::valarray<double> x0_, p_;
double f0_, fp0_, g0norm2_, alpha_, beta_max_;
LineMinimizer linmin_;
double norm2(const std::valarray<double>& v);
public:
CGOptimizer(int n): n_(n), first_step_(true), alpha_(0.0), beta_max_(0.0),
debug_print(false)
{
x0_.resize(n);
{
x0_.resize(n);
p_.resize(n);
}
......@@ -55,7 +54,7 @@ class CGOptimizer
double alpha(void) const { return alpha_; }
double alpha_start(void) const { return linmin_.alpha_start(); }
double beta_max(void) const { return beta_max_; }
void compute_xp(const std::valarray<double>& x, const double f,
void compute_xp(const std::valarray<double>& x, const double f,
const std::valarray<double>& g, std::valarray<double>& xp);
};
#endif
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