Commit 0797b6b3 by Francois Gygi

steepest descent with line minimization


git-svn-id: http://qboxcode.org/svn/qb/trunk@593 cba15fb0-1239-40c8-b417-11db7ca47a34
parent df62fad8
......@@ -3,67 +3,119 @@
// SDAIonicStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: SDAIonicStepper.C,v 1.6 2008-03-26 04:57:54 fgygi Exp $
// $Id: SDAIonicStepper.C,v 1.7 2008-04-06 17:47:12 fgygi Exp $
#include "SDAIonicStepper.h"
#include <iostream>
using namespace std;
////////////////////////////////////////////////////////////////////////////////
void SDAIonicStepper::compute_r(double e0, const vector<vector< double> >& f0)
{
// Steepest descent step
atoms_.get_positions(r0_);
double fp0;
bool wolfe1, wolfe2;
if ( first_step_ )
// check if the largest component of the force f0 is larger than max_force
const double max_force = 0.1;
double largest_force = 0.0;
for ( int is = 0; is < r0_.size(); is++ )
{
rm_ = r0_;
for ( int i = 0; i < r0_[is].size(); i++ )
{
largest_force = max(largest_force,fabs(f0[is][i]));
}
}
// copy forces f0 in linear array f_
int k = 0;
double sum = 0.0;
for ( int is = 0; is < f0.size(); is++ )
if ( largest_force > max_force )
{
for ( int i = 0; i < f0[is].size(); i++ )
if ( s_.ctxt_.onpe0() )
cout << " SDAIonicStepper: force exceeds limit, taking SD step " << endl;
// take a steepest descent step with limited displacement and exit
const double alpha_sd = max_force/largest_force;
// SD step
for ( int is = 0; is < r0_.size(); is++ )
{
f_[k++] = f0[is][i];
sum += f0[is][i]*f0[is][i];
for ( int i = 0; i < r0_[is].size(); i++ )
{
rp_[is][i] = r0_[is][i] + alpha_sd * f0[is][i];
}
}
constraints_.enforce_r(r0_,rp_);
rm_ = r0_;
r0_ = rp_;
atoms_.set_positions(r0_);
// reset the SDA algorithm
first_step_ = true;
return;
}
if ( s_.ctxt_.onpe0() )
cout << " <sda_residual> " << sum
<< "</sda_residual>" << endl;
mixer_.update(&f_[0],&theta_,&fbar_[0]);
if ( e0 > em_ )
// SDA algorithm
if ( !first_step_ )
{
theta_ = 0.0;
mixer_.restart();
wolfe1 = e0 < ec_ + fpc_ * sigma1_ * alpha_;
// fp0 = -proj(f0,pc)
fp0 = 0.0;
for ( int is = 0; is < r0_.size(); is++ )
{
for ( int i = 0; i < r0_[is].size(); i++ )
{
fp0 -= f0[is][i] * pc_[is][i];
}
}
wolfe2 = fabs(fp0) < sigma2_ * fabs(fpc_);
if ( s_.ctxt_.onpe0() )
cout << " SDAIonicStepper: sda_theta = " << theta_
<< " reset to 0" << endl;
{
cout << " SDAIonicStepper: fpc = " << fpc_ << endl;
cout << " SDAIonicStepper: fp0 = " << fp0 << endl;
cout << " SDAIonicStepper: ec = " << ec_ << " e0 = " << e0 << endl;
cout << " SDAIonicStepper: ec_ + fpc_ * sigma1_ * alpha_ ="
<< ec_ + fpc_ * sigma1_ * alpha_ << endl;
cout << " SDAIonicStepper: wolfe1/wolfe2 = "
<< wolfe1 << "/" << wolfe2 << endl;
}
}
if ( s_.ctxt_.onpe0() )
cout << " <sda_theta> " << theta_
<< "</sda_theta>" << endl;
if ( first_step_ || (wolfe1 && wolfe2) )
{
// set new descent direction
// pc = f0
fc_ = f0;
pc_ = fc_;
// fpc = d_e / d_alpha in direction pc
fpc_ = 0.0;
for ( int is = 0; is < r0_.size(); is++ )
{
for ( int i = 0; i < r0_[is].size(); i++ )
{
fpc_ -= fc_[is][i] * pc_[is][i];
}
}
ec_ = e0;
rc_ = r0_;
fp0 = fpc_;
// reset line minimizer
linmin_.reset();
}
alpha_ = linmin_.newalpha(alpha_,e0,fp0);
k = 0;
if ( s_.ctxt_.onpe0() )
cout << " SDAIonicStepper: alpha = " << alpha_ << endl;
// rp = rc + alpha * pc
for ( int is = 0; is < r0_.size(); is++ )
{
const double dt2bym = dt_ * dt_ / pmass_[is];
for ( int i = 0; i < r0_[is].size(); i++ )
{
rp_[is][i] = r0_[is][i] +
theta_ * ( r0_[is][i] - rm_[is][i] ) +
dt2bym * fbar_[k++];
rp_[is][i] = rc_[is][i] + alpha_ * pc_[is][i];
}
}
first_step_ = false;
constraints_.enforce_r(r0_,rp_);
rm_ = r0_;
r0_ = rp_;
atoms_.set_positions(r0_);
em_ = e0;
first_step_ = false;
}
......@@ -3,41 +3,34 @@
// SDAIonicStepper.h:
//
////////////////////////////////////////////////////////////////////////////////
// $Id: SDAIonicStepper.h,v 1.9 2008-03-26 04:57:54 fgygi Exp $
// $Id: SDAIonicStepper.h,v 1.10 2008-04-06 17:47:12 fgygi Exp $
#ifndef SDAIONICSTEPPER_H
#define SDAIONICSTEPPER_H
#include "IonicStepper.h"
#include "AndersonMixer.h"
#include "LineMinimizer.h"
#include <vector>
class SDAIonicStepper : public IonicStepper
{
private:
std::vector<double> f_;
std::vector<double> fbar_;
double theta_;
bool first_step_;
AndersonMixer mixer_;
double em_;
std::vector<std::vector< double> > rc_;
std::vector<std::vector< double> > pc_;
std::vector<std::vector< double> > fc_;
double ec_, fpc_;
double alpha_, sigma1_, sigma2_;
LineMinimizer linmin_;
public:
SDAIonicStepper(Sample& s) : IonicStepper(s), first_step_(true), theta_(0),
mixer_(3*atoms_.size(), 0)
{
f_.resize(3*atoms_.size());
fbar_.resize(3*atoms_.size());
rm_ = r0_;
mixer_.set_theta_max(1.0);
mixer_.set_theta_nc(1.0);
em_ = 1.0e50;
}
SDAIonicStepper(Sample& s) : IonicStepper(s), first_step_(true),
sigma1_(0.1), sigma2_(0.5) { linmin_.set_sigma1(sigma1_); }
void compute_r(double e0, const std::vector<std::vector< double> >& f0);
void compute_v(double e0, const std::vector<std::vector< double> >& f0) {}
void reset(void) { first_step_ = true; }
};
#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