Commit 0a0b7908 by Francois Gygi

Adaptive preconditioner

git-svn-id: http://qboxcode.org/svn/qb/trunk@1212 cba15fb0-1239-40c8-b417-11db7ca47a34
parent d3ea27f4
......@@ -32,7 +32,6 @@
#include "BMDIonicStepper.h"
#include "SDCellStepper.h"
#include "CGCellStepper.h"
#include "Preconditioner.h"
#include "AndersonMixer.h"
#include "MLWFTransform.h"
......@@ -201,17 +200,6 @@ void BOSampleStepper::step(int niter)
Timer tm_iter;
const bool use_preconditioner = wf_dyn == "PSD" ||
wf_dyn == "PSDA"||
wf_dyn == "JD";
Preconditioner *preconditioner = 0;
if ( use_preconditioner )
{
// create a preconditioner using the information about wf in s_.wf
// and the information about the hessian in df
preconditioner = new Preconditioner(s_,ef_);
}
WavefunctionStepper* wf_stepper = 0;
if ( wf_dyn == "SD" )
{
......@@ -227,11 +215,11 @@ void BOSampleStepper::step(int niter)
wf_stepper = new SDWavefunctionStepper(wf,dt2bye,tmap);
}
else if ( wf_dyn == "PSD" )
wf_stepper = new PSDWavefunctionStepper(wf,*preconditioner,tmap);
wf_stepper = new PSDWavefunctionStepper(wf,s_.ctrl.ecutprec,tmap);
else if ( wf_dyn == "PSDA" )
wf_stepper = new PSDAWavefunctionStepper(wf,*preconditioner,tmap);
wf_stepper = new PSDAWavefunctionStepper(wf,s_.ctrl.ecutprec,tmap);
else if ( wf_dyn == "JD" )
wf_stepper = new JDWavefunctionStepper(wf,*preconditioner,ef_,tmap);
wf_stepper = new JDWavefunctionStepper(wf,s_.ctrl.ecutprec,ef_,tmap);
// wf_stepper == 0 indicates that wf_dyn == LOCKED
......@@ -508,9 +496,6 @@ void BOSampleStepper::step(int niter)
ef_.cell_moved();
ef_.atoms_moved(); // modifications of the cell also move ions
if ( use_preconditioner )
preconditioner->update();
}
}
} // if !gs_only
......@@ -865,8 +850,8 @@ void BOSampleStepper::step(int niter)
ehart_m = ehart;
ehart = ef_.ehart();
const double delta_ehart = fabs(ehart - ehart_m);
if ( onpe0 && nite_ > 1 )
cout << " delta_ehart = " << delta_ehart << endl;
// if ( onpe0 && nite_ > 1 )
// cout << " delta_ehart = " << delta_ehart << endl;
int ite = 0;
double etotal_int;
double eigenvalue_sum;
......@@ -1232,8 +1217,6 @@ void BOSampleStepper::step(int niter)
delete ionic_stepper;
delete cell_stepper;
// delete preconditioner
if ( use_preconditioner ) delete preconditioner;
if ( ntc_extrapolation || asp_extrapolation ) delete wfmm;
initial_atomic_density = false;
......
......@@ -15,7 +15,6 @@
// JDWavefunctionStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: JDWavefunctionStepper.C,v 1.2 2009-09-08 05:36:49 fgygi Exp $
#include "JDWavefunctionStepper.h"
#include "Wavefunction.h"
......@@ -27,9 +26,17 @@ using namespace std;
////////////////////////////////////////////////////////////////////////////////
JDWavefunctionStepper::JDWavefunctionStepper(Wavefunction& wf,
Preconditioner& p, EnergyFunctional& ef, TimerMap& tmap) :
WavefunctionStepper(wf,tmap), prec_(p), wft_(wf), dwft_(wf), ef_(ef)
{}
double ecutprec, EnergyFunctional& ef, TimerMap& tmap) :
WavefunctionStepper(wf,tmap), prec_(0), wft_(wf), dwft_(wf), ef_(ef)
{
prec_ = new Preconditioner(wf,ecutprec);
}
////////////////////////////////////////////////////////////////////////////////
JDWavefunctionStepper::~JDWavefunctionStepper(void)
{
delete prec_;
}
////////////////////////////////////////////////////////////////////////////////
void JDWavefunctionStepper::update(Wavefunction& dwf)
......@@ -37,6 +44,7 @@ void JDWavefunctionStepper::update(Wavefunction& dwf)
// save copy of wf in wft_ and dwf in dwft_
wft_ = wf_;
dwft_ = dwf;
tmap_["jd_residual"].start();
for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
{
for ( int ikp = 0; ikp < wf_.nkp(); ikp++ )
......@@ -46,9 +54,7 @@ void JDWavefunctionStepper::update(Wavefunction& dwf)
// compute A = Y^T H Y and descent direction HY - YA
// proxy real matrices c, cp
DoubleMatrix c_proxy(wf_.sd(ispin,ikp)->c());
DoubleMatrix c_proxy_t(wft_.sd(ispin,ikp)->c());
DoubleMatrix cp_proxy(dwf.sd(ispin,ikp)->c());
DoubleMatrix a(c_proxy.context(),c_proxy.n(),c_proxy.n(),
c_proxy.nb(),c_proxy.nb());
......@@ -59,13 +65,43 @@ void JDWavefunctionStepper::update(Wavefunction& dwf)
// cp = cp - c * a
cp_proxy.gemm('n','n',-1.0,c_proxy,a,1.0);
}
else // real
{
// compute A = Y^T H Y and descent direction HY - YA
ComplexMatrix& c = wf_.sd(ispin,ikp)->c();
ComplexMatrix& cp = dwf.sd(ispin,ikp)->c();
ComplexMatrix a(c.context(),c.n(),c.n(),c.nb(),c.nb());
// (Y,HY)
a.gemm('c','n',1.0,c,cp,0.0);
// cp = cp - c * a
cp.gemm('n','n',-1.0,c,a,1.0);
// dwf.sd->c() now contains the descent direction (HV-VA)
}
} // ikp
} // ispin
tmap_["jd_residual"].stop();
// dwf.sd->c() now contains the descent direction (HV-VA)
prec_->update(wf_);
tmap_["jd_compute_z"].start();
for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
{
for ( int ikp = 0; ikp < wf_.nkp(); ikp++ )
{
if ( wf_.sd(ispin,ikp)->basis().real() )
{
// Apply preconditioner K and store -K(HV-VA) in wf
const valarray<double>& diag = prec_.diag(ispin,ikp);
double* c = (double*) wf_.sd(ispin,ikp)->c().valptr();
double* dc = (double*) dwf.sd(ispin,ikp)->c().valptr();
DoubleMatrix c_proxy(wf_.sd(ispin,ikp)->c());
DoubleMatrix a(c_proxy.context(),c_proxy.n(),c_proxy.n(),
c_proxy.nb(),c_proxy.nb());
DoubleMatrix c_proxy_t(wft_.sd(ispin,ikp)->c());
const int mloc = wf_.sd(ispin,ikp)->c().mloc();
const int ngwl = wf_.sd(ispin,ikp)->basis().localsize();
const int nloc = wf_.sd(ispin,ikp)->c().nloc();
......@@ -75,10 +111,10 @@ void JDWavefunctionStepper::update(Wavefunction& dwf)
// note: double mloc length for complex<double> indices
double* dcn = &dc[2*mloc*n];
double* cn = &c[2*mloc*n];
// loop to ngwl only since diag[i] is defined on [0:mloc-1]
for ( int i = 0; i < ngwl; i++ )
{
const double fac = diag[i];
const double fac = prec_->diag(ispin,ikp,n,i);
const double f0 = -fac * dcn[2*i];
const double f1 = -fac * dcn[2*i+1];
cn[2*i] = f0;
......@@ -106,25 +142,13 @@ void JDWavefunctionStepper::update(Wavefunction& dwf)
} // if real
else
{
// compute A = Y^T H Y and descent direction HY - YA
ComplexMatrix& c = wf_.sd(ispin,ikp)->c();
ComplexMatrix& ct = wft_.sd(ispin,ikp)->c();
ComplexMatrix& cp = dwf.sd(ispin,ikp)->c();
ComplexMatrix a(c.context(),c.n(),c.n(),c.nb(),c.nb());
// (Y,HY)
a.gemm('c','n',1.0,c,cp,0.0);
// cp = cp - c * a
cp.gemm('n','n',-1.0,c,a,1.0);
// dwf.sd->c() now contains the descent direction (HV-VA)
// Apply preconditioner K and store -K(HV-VA) in wf
const valarray<double>& diag = prec_.diag(ispin,ikp);
ComplexMatrix& c = wf_.sd(ispin,ikp)->c();
complex<double>* cv = c.valptr();
ComplexMatrix& cp = dwf.sd(ispin,ikp)->c();
complex<double>* cpv = cp.valptr();
ComplexMatrix& ct = wft_.sd(ispin,ikp)->c();
ComplexMatrix a(c.context(),c.n(),c.n(),c.nb(),c.nb());
const int ngwl = wf_.sd(ispin,ikp)->basis().localsize();
const int mloc = c.mloc();
const int nloc = c.nloc();
......@@ -133,10 +157,10 @@ void JDWavefunctionStepper::update(Wavefunction& dwf)
{
complex<double>* cpn = &cpv[mloc*n];
complex<double>* cn = &cv[mloc*n];
// loop to ngwl only since diag[i] is defined on [0:mloc-1]
for ( int i = 0; i < ngwl; i++ )
{
const double fac = diag[i];
const double fac = prec_->diag(ispin,ikp,n,i);
cn[i] = -fac * cpn[i];
}
}
......@@ -158,7 +182,9 @@ void JDWavefunctionStepper::update(Wavefunction& dwf)
}
} // ikp
} // ispin
tmap_["jd_compute_z"].stop();
tmap_["jd_hz"].start();
// compute HZ
const bool compute_hpsi = true;
const bool compute_forces = false;
......@@ -167,7 +193,9 @@ void JDWavefunctionStepper::update(Wavefunction& dwf)
valarray<double> sigma;
ef_.energy(compute_hpsi,dwf,compute_forces,fion,
compute_stress,sigma);
tmap_["jd_hz"].stop();
tmap_["jd_blocks"].start();
for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
{
for ( int ikp = 0; ikp < wf_.nkp(); ikp++ )
......@@ -291,4 +319,5 @@ void JDWavefunctionStepper::update(Wavefunction& dwf)
} // if real
} // ikp
} // ispin
tmap_["jd_blocks"].stop();
}
......@@ -15,7 +15,6 @@
// JDWavefunctionStepper.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: JDWavefunctionStepper.h,v 1.2 2009-09-08 05:36:49 fgygi Exp $
#ifndef JDWAVEFUNCTIONSTEPPER_H
#define JDWAVEFUNCTIONSTEPPER_H
......@@ -29,7 +28,7 @@ class JDWavefunctionStepper : public WavefunctionStepper
{
private:
Preconditioner& prec_;
Preconditioner *prec_;
EnergyFunctional& ef_;
Wavefunction wft_, dwft_;
......@@ -38,8 +37,8 @@ class JDWavefunctionStepper : public WavefunctionStepper
void update(Wavefunction& dwf);
virtual void preprocess(void) {}
JDWavefunctionStepper(Wavefunction& wf, Preconditioner& p,
JDWavefunctionStepper(Wavefunction& wf, double ecutprec,
EnergyFunctional& ef, TimerMap& tmap);
~JDWavefunctionStepper() {};
~JDWavefunctionStepper();
};
#endif
......@@ -287,12 +287,26 @@ BOSampleStepper.o: JDWavefunctionStepper.h SDIonicStepper.h IonicStepper.h
BOSampleStepper.o: Species.h SDAIonicStepper.h LineMinimizer.h
BOSampleStepper.o: CGIonicStepper.h CGOptimizer.h MDIonicStepper.h
BOSampleStepper.o: BMDIonicStepper.h SDCellStepper.h CellStepper.h
BOSampleStepper.o: CGCellStepper.h Preconditioner.h AndersonMixer.h
BOSampleStepper.o: MLWFTransform.h BasisMapping.h
BOSampleStepper.o: CGCellStepper.h AndersonMixer.h MLWFTransform.h
BOSampleStepper.o: BasisMapping.h
BOSampleStepper.o: SampleStepper.h Timer.h EnergyFunctional.h ChargeDensity.h
BOSampleStepper.o: Context.h StructureFactor.h Sample.h AtomSet.h Atom.h
BOSampleStepper.o: D3vector.h UnitCell.h ConstraintSet.h ExtForceSet.h
BOSampleStepper.o: Wavefunction.h Control.h
BOSampleStepper-merge.o: BOSampleStepper.h SampleStepper.h Timer.h
BOSampleStepper-merge.o: EnergyFunctional.h ChargeDensity.h Context.h
BOSampleStepper-merge.o: StructureFactor.h Sample.h AtomSet.h Atom.h
BOSampleStepper-merge.o: D3vector.h UnitCell.h ConstraintSet.h ExtForceSet.h
BOSampleStepper-merge.o: Wavefunction.h Control.h SlaterDet.h Basis.h
BOSampleStepper-merge.o: Matrix.h WavefunctionStepper.h
BOSampleStepper-merge.o: SDWavefunctionStepper.h PSDWavefunctionStepper.h
BOSampleStepper-merge.o: PSDAWavefunctionStepper.h JDWavefunctionStepper.h
BOSampleStepper-merge.o: SDIonicStepper.h IonicStepper.h Species.h
BOSampleStepper-merge.o: SDAIonicStepper.h LineMinimizer.h CGIonicStepper.h
BOSampleStepper-merge.o: CGOptimizer.h MDIonicStepper.h BMDIonicStepper.h
BOSampleStepper-merge.o: SDCellStepper.h CellStepper.h CGCellStepper.h
BOSampleStepper-merge.o: Preconditioner.h AndersonMixer.h MLWFTransform.h
BOSampleStepper-merge.o: BasisMapping.h
BtHF.o: Sample.h AtomSet.h Context.h Atom.h D3vector.h UnitCell.h
BtHF.o: ConstraintSet.h ExtForceSet.h Wavefunction.h Control.h
CellDyn.o: Sample.h AtomSet.h Context.h Atom.h D3vector.h UnitCell.h
......@@ -511,11 +525,8 @@ PlotCmd.o: UnitCell.h ConstraintSet.h ExtForceSet.h Wavefunction.h Control.h
PositionConstraint.o: PositionConstraint.h Constraint.h AtomSet.h Context.h
PositionConstraint.o: Atom.h D3vector.h UnitCell.h Species.h
PositionConstraint.o: Constraint.h
Preconditioner.o: Preconditioner.h EnergyFunctional.h ChargeDensity.h Timer.h
Preconditioner.o: Context.h StructureFactor.h Sample.h AtomSet.h Atom.h
Preconditioner.o: D3vector.h UnitCell.h ConstraintSet.h ExtForceSet.h
Preconditioner.o: Wavefunction.h Control.h Basis.h SlaterDet.h Matrix.h
Preconditioner.o: ConfinementPotential.h
Preconditioner.o: Preconditioner.h Basis.h D3vector.h UnitCell.h Context.h
Preconditioner.o: Wavefunction.h SlaterDet.h Matrix.h Timer.h
PrintCmd.o: UserInterface.h Sample.h AtomSet.h Context.h Atom.h D3vector.h
PrintCmd.o: UnitCell.h ConstraintSet.h ExtForceSet.h Wavefunction.h Control.h
PSDAWavefunctionStepper.o: PSDAWavefunctionStepper.h WavefunctionStepper.h
......
......@@ -15,7 +15,6 @@
// PSDAWavefunctionStepper.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: PSDAWavefunctionStepper.h,v 1.7 2008-09-08 15:56:18 fgygi Exp $
#ifndef PSDAWAVEFUNCTIONSTEPPER_H
#define PSDAWAVEFUNCTIONSTEPPER_H
......@@ -28,7 +27,7 @@ class PSDAWavefunctionStepper : public WavefunctionStepper
{
private:
Preconditioner& prec_;
Preconditioner *prec_;
Wavefunction wf_last_, dwf_last_;
// Anderson acceleration flag
......@@ -39,7 +38,7 @@ class PSDAWavefunctionStepper : public WavefunctionStepper
void update(Wavefunction& dwf);
virtual void preprocess(void) { extrapolate_ = false; }
PSDAWavefunctionStepper(Wavefunction& wf, Preconditioner& p, TimerMap& tmap);
~PSDAWavefunctionStepper() {};
PSDAWavefunctionStepper(Wavefunction& wf, double ecutprec, TimerMap& tmap);
~PSDAWavefunctionStepper();
};
#endif
......@@ -15,7 +15,6 @@
// PSDWavefunctionStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: PSDWavefunctionStepper.C,v 1.12 2008-09-08 15:56:18 fgygi Exp $
#include "PSDWavefunctionStepper.h"
#include "Wavefunction.h"
......@@ -26,9 +25,16 @@ using namespace std;
////////////////////////////////////////////////////////////////////////////////
PSDWavefunctionStepper::PSDWavefunctionStepper(Wavefunction& wf,
Preconditioner& p, TimerMap& tmap) :
WavefunctionStepper(wf,tmap), prec_(p)
{}
double ecutprec, TimerMap& tmap) : WavefunctionStepper(wf,tmap)
{
prec_ = new Preconditioner(wf,ecutprec);
}
////////////////////////////////////////////////////////////////////////////////
PSDWavefunctionStepper::~PSDWavefunctionStepper(void)
{
delete prec_;
}
////////////////////////////////////////////////////////////////////////////////
void PSDWavefunctionStepper::update(Wavefunction& dwf)
......@@ -38,7 +44,6 @@ void PSDWavefunctionStepper::update(Wavefunction& dwf)
for ( int ikp = 0; ikp < wf_.nkp(); ikp++ )
{
// compute A = V^T H V and descent direction HV - VA
tmap_["psd_residual"].start();
if ( wf_.sd(ispin,ikp)->basis().real() )
{
......@@ -66,12 +71,19 @@ void PSDWavefunctionStepper::update(Wavefunction& dwf)
cp.gemm('n','n',-1.0,c,a,1.0);
}
tmap_["psd_residual"].stop();
}
}
// dwf.sd->c() now contains the descent direction (HV-VA)
// dwf.sd->c() now contains the descent direction (HV-VA)
tmap_["psd_update_wf"].start();
const valarray<double>& diag = prec_.diag(ispin,ikp);
// update preconditioner using the residual
prec_->update(dwf);
for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
{
for ( int ikp = 0; ikp < wf_.nkp(); ikp++ )
{
tmap_["psd_update_wf"].start();
double* coeff = (double*) wf_.sd(ispin,ikp)->c().valptr();
const double* dcoeff =
(const double*) dwf.sd(ispin,ikp)->c().cvalptr();
......@@ -83,10 +95,10 @@ void PSDWavefunctionStepper::update(Wavefunction& dwf)
// note: double mloc length for complex<double> indices
double* c = &coeff[2*mloc*n];
const double* dc = &dcoeff[2*mloc*n];
// loop to ngwl only since diag[i] is defined on [0:mloc-1]
for ( int i = 0; i < ngwl; i++ )
{
const double fac = diag[i];
const double fac = prec_->diag(ispin,ikp,n,i);
const double delta_re = fac * dc[2*i];
const double delta_im = fac * dc[2*i+1];
c[2*i] -= delta_re;
......
......@@ -15,7 +15,6 @@
// PSDWavefunctionStepper.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: PSDWavefunctionStepper.h,v 1.7 2008-09-08 15:56:19 fgygi Exp $
#ifndef PSDWAVEFUNCTIONSTEPPER_H
#define PSDWAVEFUNCTIONSTEPPER_H
......@@ -27,13 +26,13 @@ class PSDWavefunctionStepper : public WavefunctionStepper
{
private:
Preconditioner& prec_;
Preconditioner *prec_;
public:
void update(Wavefunction& dwf);
PSDWavefunctionStepper(Wavefunction& wf, Preconditioner& p, TimerMap& tmap);
~PSDWavefunctionStepper() {};
PSDWavefunctionStepper(Wavefunction& wf, double ecutprec, TimerMap& tmap);
~PSDWavefunctionStepper();
};
#endif
......@@ -15,61 +15,122 @@
// Preconditioner.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: Preconditioner.C,v 1.8 2008-09-08 15:56:19 fgygi Exp $
#include "Preconditioner.h"
#include "EnergyFunctional.h"
#include "Sample.h"
#include "Basis.h"
#include "Wavefunction.h"
#include "SlaterDet.h"
#include "ConfinementPotential.h"
using namespace std;
////////////////////////////////////////////////////////////////////////////////
Preconditioner::Preconditioner(const Sample& s, const EnergyFunctional& ef) :
s_(s), ef_(ef)
Preconditioner::Preconditioner(const Wavefunction& wf, double ecutprec)
: ecutprec_(ecutprec)
{
update();
kpg2_.resize(wf.nspin());
ekin_.resize(wf.nspin());
for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
{
kpg2_[ispin].resize(wf.nkp());
ekin_[ispin].resize(wf.nkp());
for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
{
const SlaterDet& sd = *(wf.sd(ispin,ikp));
const Basis& wfbasis = sd.basis();
kpg2_[ispin][ikp] = wfbasis.kpg2_ptr();
ekin_[ispin][ikp].resize(sd.nstloc());
}
}
update(wf);
}
////////////////////////////////////////////////////////////////////////////////
void Preconditioner::update(void)
double Preconditioner::diag(int ispin, int ikp, int n, int ig) const
{
// reinitialize preconditioner
bool use_confinement = s_.ctrl.ecuts > 0.0;
const Wavefunction& wf = s_.wf;
// If ecutprec is zero, use ecut
const double ecutpr = s_.ctrl.ecutprec > 0.0 ? s_.ctrl.ecutprec : wf.ecut();
if ( ecutprec_ == 0.0 )
{
const double ekin_n = ekin_[ispin][ikp][n];
if ( ekin_n == 0.0 )
return 0.0;
const double q2 = kpg2_[ispin][ikp][ig];
#if 0
// Payne Teter Allan adaptive preconditioner
const double x = 0.5 * q2 / ekin_n;
const double num = 27.0 + x * ( 18.0 + x * ( 12 + x * 8 ));
const double den = num + 16.0 * x*x*x*x;
return (num/den);
#else
// modified Payne Teter Allan adaptive preconditioner (Kresse)
const double x = 0.5 * q2 / ( 1.5 * ekin_n );
const double num = 27.0 + x * ( 18.0 + x * ( 12 + x * 8 ));
const double den = num + 16.0 * x*x*x*x;
return ( 2.0 / (1.5 * ekin_n )) * (num/den);
#endif
}
else
{
// next lines: inclusion of fstress if confinement potential is used
// const valarray<double>& fstress = ef_.confpot(ikp)->fstress();
// double e = 0.5 * ( kpg2_ptr[ig] + fstress[ig] );
// diag_[ispin][ikp][ig] = ( e < ecutpr ) ? 0.5 / ecutpr : 0.5 / e;
const double e = 0.5 * kpg2_[ispin][ikp][ig];
return ( e < ecutprec_ ) ? 0.5 / ecutprec_ : 0.5 / e;
}
}
diag_.resize(wf.nspin());
////////////////////////////////////////////////////////////////////////////////
void Preconditioner::update(const Wavefunction& wf)
{
// update the kinetic energy ekin_[ispin][ikp][n] of states in wf
for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
{
diag_[ispin].resize(wf.nkp());
for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
{
// Only resize and initialize diag_ if ikp is active on this task
const Basis& basis = wf.sd(ispin,ikp)->basis();
const int ngwloc = basis.localsize();
diag_[ispin][ikp].resize(ngwloc);
const double *kpg2_ptr = basis.kpg2_ptr();
const SlaterDet& sd = *(wf.sd(ispin,ikp));
const Basis& wfbasis = sd.basis();
// factor fac in next lines: 2.0 for G and -G (if basis is real) and
// 0.5 from 1/(2m)
const double fac = wfbasis.real() ? 1.0 : 0.5;
const ComplexMatrix& c = sd.c();
const Context& sdctxt = sd.context();
if ( use_confinement )
const int ngwloc = wfbasis.localsize();
const complex<double>* p = c.cvalptr();
const int mloc = c.mloc();
const int nloc = c.nloc();
valarray<double> buf(2*nloc);
const double* pkpg2 = kpg2_[ispin][ikp];
for ( int n = 0; n < nloc; n++ )
{
const valarray<double>& fstress = ef_.confpot(ikp)->fstress();
double sum_norm = 0.0;
double sum_ekin = 0.0;
for ( int ig = 0; ig < ngwloc; ig++ )
{
double e = 0.5 * ( kpg2_ptr[ig] + fstress[ig] );
diag_[ispin][ikp][ig] = ( e < ecutpr ) ? 0.5 / ecutpr : 0.5 / e;
const double psi2 = norm(p[ig+n*mloc]);
sum_norm += psi2;
sum_ekin += psi2 * pkpg2[ig];
}
// correct for double counting of G=0 in norm
if ( sdctxt.myrow() == 0 )
sum_norm -= 0.5 * norm(p[n*mloc]);
// store norm in buf[n] and ekin in buf[n+nloc]
buf[n] = fac * sum_norm;
buf[n+nloc] = fac * sum_ekin;
}
else
sdctxt.dsum('C',2*nloc,1,&buf[0],2*nloc);
for ( int n = 0; n < nloc; n++ )
ekin_[ispin][ikp][n] = buf[n] > 0.0 ? buf[n+nloc]/buf[n] : 0;
#ifdef DEBUG
if ( sdctxt.onpe0() )
{
for ( int ig = 0; ig < ngwloc; ig++ )
{
double e = 0.5 * kpg2_ptr[ig];
diag_[ispin][ikp][ig] = ( e < ecutpr ) ? 0.5 / ecutpr : 0.5 / e;
}
for ( int n = 0; n < nloc; n++ )
cout << "Preconditioner::update ekin[" << n << "] = "
<< ekin_[ispin][ikp][n] << endl;
}
#endif
}
}
}
......@@ -15,13 +15,11 @@
// Preconditioner.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: Preconditioner.h,v 1.5 2008-09-08 15:56:19 fgygi Exp $
#ifndef PRECONDITIONER_H
#define PRECONDITIONER_H
class Sample;
class EnergyFunctional;
class Wavefunction;
#include <vector>
#include <valarray>
......@@ -30,19 +28,21 @@ class Preconditioner
{
private:
const Sample& s_;
const EnergyFunctional& ef_;
// diag_[ispin][ikp][ig]
std::vector<std::vector<std::valarray<double> > > diag_;
// kpg2_[ispin][ikp][ig]
std::vector<std::vector<const double *> > kpg2_;
// ekin_[ispin][ikp][n]
std::vector<std::vector<std::valarray<double> > > ekin_;
double ecutprec_;
public:
void update(void);
// update values of ekin_
void update(const Wavefunction& wf);
const std::valarray<double>& diag(int ispin, int ikp) const
{ return diag_[ispin][ikp]; }
double diag(int ispin, int ikp, int n, int ig) const;
Preconditioner(const Sample& s, const EnergyFunctional& ef);
Preconditioner(const Wavefunction& wf, double ecutprec);
//~Preconditioner();
};
#endif
--------------------------------------------------------------------------------
wrk
Adaptive preconditioner.
r1211: tolerance for non-self consistent iterations
--------------------------------------------------------------------------------
rel1_55_6
Redesign of Bisection and ExchangeOperator separates bisection transformation
and load-balancing permutation.
Fixed KPGridConnectivity
......
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