Commit f45b5044 by Francois Gygi

Implemented multidimensional Anderson acceleration.


git-svn-id: http://qboxcode.org/svn/qb/trunk@680 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 602e7022
......@@ -15,67 +15,176 @@
// AndersonMixer.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: AndersonMixer.C,v 1.7 2008-09-08 15:56:17 fgygi Exp $
// $Id: AndersonMixer.C,v 1.8 2009-03-08 01:10:30 fgygi Exp $
#include "AndersonMixer.h"
#include "blas.h"
#include <iostream>
using namespace std;
////////////////////////////////////////////////////////////////////////////////
void AndersonMixer::restart(void)
{
extrapolate_ = false;
n_ = -1;
k_ = -1;
}
////////////////////////////////////////////////////////////////////////////////
void AndersonMixer::update(const double* f, double* theta, double* fbar)
void AndersonMixer::update(double* x, double* f, double* xbar, double* fbar)
{
*theta = 0.0;
if ( extrapolate_ )
// update:
// input: x, f
// output: xbar, fbar
//
// Computes the pair (xbar,fbar) using pairs (x,f) used
// in previous updates, according to the Anderson algorithm.
//
// increment index of current vector
k_ = ( k_ + 1 ) % ( nmax_ + 1 );
// increment current number of vectors
if ( n_ < nmax_ ) n_++;
// save vectors
for ( int i = 0; i < m_; i++ )
{
x_[k_][i] = x[i];
f_[k_][i] = f[i];
}
valarray<double> a;
valarray<double> b;
valarray<double> theta;
if ( n_ > 0 )
{
// compute matrix A = F^T F and rhs b = F^T f
// compute the lower part of A only (i>=j)
a.resize(n_*n_);
b.resize(n_);
theta.resize(n_);
for ( int i = 0; i < n_; i++ )
{
int ione=1;
valarray<double> tmp0(n_);
const int kmi = ( k_ - i + nmax_ ) % ( nmax_ + 1 );
assert(kmi>=0);
assert(kmi<nmax_+1);
// cout << "k=" << k_ << " i=" << i << " kmi=" << kmi << endl;
for ( int j = 0; j <= i; j++ )
{
const int kmj = ( k_ - j + nmax_ ) % ( nmax_ + 1 );
assert(kmj>=0);
assert(kmj<nmax_+1);
// cout << "k=" << k_ << " j=" << j << " kmj=" << kmj << endl;
double sum = 0.0;
for ( int l = 0; l < m_; l++ )
sum += (f_[k_][l] - f_[kmi][l]) * (f_[k_][l] - f_[kmj][l]);
a[i+j*n_] = sum;
}
// compute theta = - a / b
// tmp0 = delta_F = f - flast
for ( int i = 0; i < n_; i++ )
tmp0[i] = f[i] - flast_[i];
double bsum = 0.0;
for ( int l = 0; l < m_; l++ )
bsum += ( f_[k_][l] - f_[kmi][l] ) * f_[k_][l];
b[i] = bsum;
}
// a = delta_F * F
double a = ddot(&n_, &tmp0[0], &ione, f, &ione);
#if 0
// print matrix a and rhs b
for ( int i = 0; i < n_; i++ )
for ( int j = 0; j <=i; j++ )
cout << "a("<<i<<","<<j<<")=" << a[i+j*n_] << endl;
for ( int i = 0; i < n_; i++ )
cout << "b("<<i<<")=" << b[i] << endl;
#endif
// b = delta_F * delta_F
double b = ddot(&n_, &tmp0[0], &ione, &tmp0[0], &ione);
if ( pctxt_ != 0 )
{
pctxt_->dsum(n_*n_,1,&a[0],n_*n_);
pctxt_->dsum(n_,1,&b[0],n_);
}
if ( pctxt_ != 0 )
// solve the linear system a * theta = b
const bool diag = false;
if ( diag )
{
// solve the linear system using eigenvalues and eigenvectors
// compute eigenvalues of a
char jobz = 'V';
char uplo = 'L';
valarray<double> w(n_);
int lwork = 3*n_;
valarray<double> work(lwork);
int info;
dsyev(&jobz,&uplo,&n_,&a[0],&n_,&w[0],&work[0],&lwork,&info);
for ( int i = 0; i < n_; i++ )
cout << w[i] << " ";
cout << endl;
if ( info != 0 )
{
double work[2] = { a, b };
pctxt_->dsum(2,1,work,2);
a = work[0];
b = work[1];
cerr << " AndersonMixer: Error in dsyev" << endl;
cerr << " info = " << info << endl;
exit(1);
}
if ( b != 0.0 )
*theta = - a / b;
else
*theta = 0.0;
// test if residual has increased
if ( *theta <= -1.0 )
// solve for theta
// theta_i = sum_j
for ( int k = 0; k < n_; k++ )
{
*theta = theta_nc_;
// correct only if eigenvalue w[k] is large enough compared to the
// largest eigenvalue
if ( w[n_-1] < 5.0 * w[k] )
{
const double fac = 1.0 / w[k];
for ( int i = 0; i < n_; i++ )
for ( int j = 0; j < n_; j++ )
theta[i] += fac * a[i+k*n_] * a[j+k*n_] * b[j];
}
}
}
else
{
// solve the linear system directly
char uplo = 'L';
int nrhs = 1;
valarray<int> ipiv(n_);
valarray<double> work(n_);
int info;
dsysv(&uplo,&n_,&nrhs,&a[0],&n_,&ipiv[0],&b[0],&n_,&work[0],&n_,&info);
if ( info != 0 )
{
cerr << " AndersonMixer: Error in dsysv" << endl;
cerr << " info = " << info << endl;
exit(1);
}
// the vector b now contains the solution
theta = b;
if ( pctxt_ != 0 )
{
if ( pctxt_->onpe0() )
{
cout << " AndersonMixer: theta = ";
for ( int i = 0; i < theta.size(); i++ )
cout << theta[i] << " ";
cout << endl;
}
}
*theta = min(theta_max_,*theta);
}
} // if n_ > 0
// fbar = f + theta * ( f - flast )
// flast_ = f_
for ( int i = 0; i < n_; i++ )
// fbar = f[k] + sum_i theta_i * ( f[k] - f[kmi] )
// xbar = x[k] + sum_i theta_i * ( x[k] - x[kmi] )
for ( int l = 0; l < m_; l++ )
{
fbar[l] = f_[k_][l];
xbar[l] = x_[k_][l];
}
for ( int i = 0; i < n_; i++ )
{
const int kmi = ( k_ - i + nmax_ ) % ( nmax_ + 1 );
assert(kmi>=0);
assert(kmi<nmax_+1);
for ( int l = 0; l < m_; l++ )
{
const double ftmp = f[i];
fbar[i] = ftmp + *theta * ( ftmp - flast_[i] );
flast_[i] = ftmp;
fbar[l] -= theta[i] * ( f_[k_][l] - f_[kmi][l] );
xbar[l] -= theta[i] * ( x_[k_][l] - x_[kmi][l] );
}
extrapolate_ = true;
}
}
......@@ -15,39 +15,46 @@
// AndersonMixer.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: AndersonMixer.h,v 1.7 2008-09-08 15:56:17 fgygi Exp $
// $Id: AndersonMixer.h,v 1.8 2009-03-08 01:10:30 fgygi Exp $
#ifndef ANDERSONMIXER_H
#define ANDERSONMIXER_H
#include <vector>
#include <valarray>
#include <cassert>
#include "Context.h"
class AndersonMixer
{
int n_; // size of vectors
// nmax is the dimension of the subspace of previous search directions
// nmax=0: use simple mixing (no acceleration)
// nmax=1: use one previous direction
int m_; // dimension of vectors
int nmax_; // maximum number of vectors (without current)
int n_; // number of vectors
int k_; // index of current vector
const Context* const pctxt_; // pointer to relevant Context, null if local
double theta_max_; // maximum extrapolation
double theta_nc_; // negative curvature value
std::valarray<double> flast_; // last residual
bool extrapolate_; // state variable
std::vector<std::valarray<double> > x_,f_;
public:
AndersonMixer(const int n, const Context* const pctxt) :
n_(n), pctxt_(pctxt), extrapolate_(false), theta_max_(2.0), theta_nc_(0.0)
AndersonMixer(const int m, const int nmax, const Context* const pctxt) :
m_(m), nmax_(nmax), pctxt_(pctxt)
{
assert( n > 0 );
flast_.resize(n);
assert( nmax >= 0 );
x_.resize(nmax_+1);
f_.resize(nmax_+1);
for ( int n = 0; n < nmax_+1; n++ )
{
x_[n].resize(m_);
f_[n].resize(m_);
}
restart();
}
void update(const double* f, double* theta, double* fbar);
void update(double* x, double* f, double* xbar, double* fbar);
void restart(void);
void set_theta_max(double theta_max) { theta_max_ = theta_max; }
void set_theta_nc(double theta_nc) { theta_nc_ = theta_nc; }
double theta_max(void) const { return theta_max_; }
double theta_nc(void) const { return theta_nc_; }
};
#endif
......@@ -15,7 +15,7 @@
// BOSampleStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: BOSampleStepper.C,v 1.47 2008-09-15 14:57:11 fgygi Exp $
// $Id: BOSampleStepper.C,v 1.48 2009-03-08 01:10:59 fgygi Exp $
#include "BOSampleStepper.h"
#include "EnergyFunctional.h"
......@@ -73,8 +73,7 @@ void BOSampleStepper::step(int niter)
{
const bool onpe0 = s_.ctxt_.onpe0();
const bool anderson_charge_mixing =
( s_.ctrl.debug.find("AND_CHMIX") != string::npos );
const bool anderson_charge_mixing = ( s_.ctrl.charge_mix_ndim > 0 );
// determine whether eigenvectors must be computed
// eigenvectors are computed if explicitly requested with wf_diag==T
......@@ -196,6 +195,40 @@ void BOSampleStepper::step(int niter)
mlwft = new MLWFTransform(*wf.sd(0,0));
}
// Charge mixing variables
vector<complex<double> > rhog_in(cd_.rhog[0]);
vector<complex<double> > drhog(rhog_in.size());
vector<complex<double> > rhobar(rhog_in.size());
vector<complex<double> > drhobar(rhog_in.size());
vector<double> wkerker(rhog_in.size());
vector<double> wkerker_inv(rhog_in.size());
const int anderson_ndim = s_.ctrl.charge_mix_ndim;
AndersonMixer mixer(2*rhog_in.size(),anderson_ndim,&cd_.vcontext());
// compute Kerker preconditioning
// real space Kerker cutoff in a.u.
const double rc_Kerker = s_.ctrl.charge_mix_rcut;
const double *const g2 = cd_.vbasis()->g2_ptr();
const double *const g2i = cd_.vbasis()->g2i_ptr();
if ( rc_Kerker > 0.0 )
{
const double q0_kerker = 2 * M_PI / rc_Kerker;
const double q0_kerker2 = q0_kerker * q0_kerker;
for ( int i=0; i < wkerker.size(); i++ )
{
wkerker[i] = g2[i] / ( g2[i] + q0_kerker2 );
wkerker_inv[i] = g2i[i] * ( g2[i] + q0_kerker2 );
}
}
else
{
for ( int i=0; i < wkerker.size(); i++ )
{
wkerker[i] = 1.0;
wkerker_inv[i] = 1.0;
}
}
// Next line: special case of niter=0: compute GS only
for ( int iter = 0; iter < max(niter,1); iter++ )
{
......@@ -541,12 +574,6 @@ void BOSampleStepper::step(int niter)
if ( wf_stepper != 0 )
{
assert(cd_.rhog.size()==1); // works for nspin=1 only
vector<complex<double> > rhog_current(cd_.rhog[0]);
vector<complex<double> > rhog_last(rhog_current);
vector<complex<double> > drhog(rhog_current.size());
vector<complex<double> > drhog_bar(rhog_current.size());
AndersonMixer mixer(2*rhog_current.size(),&cd_.vcontext());
mixer.set_theta_max(2.0);
wf_stepper->preprocess();
for ( int itscf = 0; itscf < nitscf_; itscf++ )
......@@ -564,74 +591,42 @@ void BOSampleStepper::step(int niter)
{
if ( itscf == 0 )
{
rhog_current = cd_.rhog[0];
rhog_last = cd_.rhog[0];
}
// compute correction drhog
for ( int i=0; i < rhog_current.size(); i++ )
{
drhog[i] = (cd_.rhog[0][i] - rhog_current[i]);
}
// Apply Kerker preconditioner to drhog
// Use Kerker preconditioning if rc_Kerker > 0.0,
// no preconditioning otherwise
const double alpha = s_.ctrl.charge_mix_coeff;
const double *const g2 = cd_.vbasis()->g2_ptr();
// real space Kerker cutoff in a.u.
const double rc_Kerker = s_.ctrl.charge_mix_rcut;
if ( rc_Kerker > 0.0 )
{
const double q0_kerker = 2 * M_PI / rc_Kerker;
const double q0_kerker2 = q0_kerker * q0_kerker;
for ( int i=0; i < rhog_current.size(); i++ )
{
drhog[i] *= alpha * g2[i] / ( g2[i] + q0_kerker2 );
}
// at first scf iteration, nothing to mix
// memorize rhog in rhog_in
rhog_in = cd_.rhog[0];
}
else
{
for ( int i=0; i < rhog_current.size(); i++ )
// itscf > 0
// compute unscreened correction drhog
for ( int i=0; i < rhog_in.size(); i++ )
drhog[i] = (cd_.rhog[0][i] - rhog_in[i]);
const double alpha = s_.ctrl.charge_mix_coeff;
// Anderson acceleration
if ( anderson_charge_mixing )
{
drhog[i] *= alpha;
}
}
// row weighting of LS calculation (preconditioning)
for ( int i=0; i < drhog.size(); i++ )
drhog[i] *= wkerker_inv[i];
// Anderson acceleration
double theta = 0.0;
for ( int i=0; i < rhog_current.size(); i++ )
{
drhog_bar[i] = drhog[i];
}
mixer.update((double*)&rhog_in[0],(double*)&drhog[0],
(double*)&rhobar[0],(double*)&drhobar[0]);
for ( int i=0; i < drhog.size(); i++ )
drhobar[i] *= wkerker[i];
if ( anderson_charge_mixing )
{
mixer.update((double*)&drhog[0],&theta,(double*)&drhog_bar[0]);
if ( onpe0 )
for ( int i=0; i < rhog_in.size(); i++ )
rhog_in[i] = rhobar[i] + alpha * drhobar[i];
}
else
{
cout << " Charge mixing: Anderson theta="
<< theta << endl;
for ( int i=0; i < rhog_in.size(); i++ )
rhog_in[i] += alpha * drhog[i] * wkerker[i];
}
}
// update rhog_current
// rhog_current = rhog_current + theta*(rhog_current-rhog_last)
for ( int i=0; i < rhog_current.size(); i++ )
{
complex<double> rhotmp = rhog_current[i];
rhog_current[i] += theta * (rhog_current[i] - rhog_last[i]);
rhog_last[i] = rhotmp;
}
// Apply correction
for ( int i=0; i < rhog_current.size(); i++ )
{
cd_.rhog[0][i] = rhog_current[i] + drhog_bar[i];
cd_.rhog[0] = rhog_in;
cd_.update_rhor();
}
rhog_current = cd_.rhog[0];
cd_.update_rhor();
} // if nite > 1
ef_.update_vhxc();
......
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