Commit ca1fe06c by Francois Gygi

Implemented complex case


git-svn-id: http://qboxcode.org/svn/qb/trunk@723 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 4b25eb1a
......@@ -15,7 +15,7 @@
// JDWavefunctionStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: JDWavefunctionStepper.C,v 1.1 2009-09-07 19:24:41 fgygi Exp $
// $Id: JDWavefunctionStepper.C,v 1.2 2009-09-08 05:36:49 fgygi Exp $
#include "JDWavefunctionStepper.h"
#include "Wavefunction.h"
......@@ -48,7 +48,6 @@ void JDWavefunctionStepper::update(Wavefunction& dwf)
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 cp_proxy_t(dwft_.sd(ispin,ikp)->c());
DoubleMatrix a(c_proxy.context(),c_proxy.n(),c_proxy.n(),
c_proxy.nb(),c_proxy.nb());
......@@ -107,8 +106,55 @@ void JDWavefunctionStepper::update(Wavefunction& dwf)
} // if real
else
{
// not implemented in the complex case
assert(false);
// 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);
complex<double>* cv = c.valptr();
complex<double>* cpv = cp.valptr();
const int ngwl = wf_.sd(ispin,ikp)->basis().localsize();
const int mloc = c.mloc();
const int nloc = c.nloc();
for ( int n = 0; n < nloc; n++ )
{
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];
cn[i] = -fac * cpn[i];
}
}
// wf_ now contains the preconditioned descent
// direction Z = -K(HY-YA)
// orthogonalize Z to Y
// Z = Z - YY^T Z
// A = Y^T * Z
a.gemm('c','n',1.0,ct,c,0.0);
// Z = Z - Y * A
c.gemm('n','n',-1.0,ct,a,1.0);
// orthogonalize Z: gram(Z)
wf_.sd(ispin,ikp)->gram();
// wf now contains Z, orthonormal
}
} // ikp
} // ispin
......@@ -143,7 +189,7 @@ void JDWavefunctionStepper::update(Wavefunction& dwf)
DoubleMatrix h(c_proxy.context(),2*c_proxy.n(),2*c_proxy.n(),
c_proxy.nb(),c_proxy.nb());
// (Y,HY)
// (Y,HY)
// factor 2.0 in next line: G and -G
a.gemm('t','n',2.0,c_proxy_t,cp_proxy_t,0.0);
// rank-1 update correction
......@@ -151,22 +197,22 @@ void JDWavefunctionStepper::update(Wavefunction& dwf)
// a contains (Y,HY), copy to h11 block
h.getsub(a,a.m(),a.n(),0,0,0,0);
// (Z,HY)
// (Z,HY)
// factor 2.0 in next line: G and -G
a.gemm('t','n',2.0,c_proxy,cp_proxy_t,0.0);
// rank-1 update correction
a.ger(-1.0,c_proxy,0,cp_proxy_t,0);
// a contains (Z,HY), copy to h21 block
h.getsub(a,a.m(),a.n(),0,0,a.m(),0);
// (Z,HZ)
// (Z,HZ)
// factor 2.0 in next line: G and -G
a.gemm('t','n',2.0,c_proxy,cp_proxy,0.0);
// rank-1 update correction
a.ger(-1.0,c_proxy,0,cp_proxy,0);
// a contains (Z,HZ), copy to h22 block
h.getsub(a,a.m(),a.n(),0,0,a.m(),a.n());
// diagonalize h
// Note: we only need the first n eigenvectors of the (2n x 2n) matrix
valarray<double> w(h.m());
......@@ -192,9 +238,57 @@ void JDWavefunctionStepper::update(Wavefunction& dwf)
} // if real
else
{
// not implemented in the complex case
assert(false);
}
// wf now contains Z
// dwf now contains HZ
// compute blocks (Y,HY), (Y,HZ), (Z,HZ)
ComplexMatrix& c = wf_.sd(ispin,ikp)->c();
ComplexMatrix& ct = wft_.sd(ispin,ikp)->c();
ComplexMatrix& cp = dwf.sd(ispin,ikp)->c();
ComplexMatrix& cpt = dwft_.sd(ispin,ikp)->c();
ComplexMatrix a(c.context(),c.n(),c.n(),c.nb(),c.nb());
ComplexMatrix h(c.context(),2*c.n(),2*c.n(),c.nb(),c.nb());
// (Y,HY)
// factor 2.0 in next line: G and -G
a.gemm('c','n',1.0,ct,cpt,0.0);
// a contains (Y,HY), copy to h11 block
h.getsub(a,a.m(),a.n(),0,0,0,0);
// (Z,HY)
a.gemm('c','n',1.0,c,cpt,0.0);
// a contains (Z,HY), copy to h21 block
h.getsub(a,a.m(),a.n(),0,0,a.m(),0);
// (Z,HZ)
a.gemm('c','n',1.0,c,cp,0.0);
// a contains (Z,HZ), copy to h22 block
h.getsub(a,a.m(),a.n(),0,0,a.m(),a.n());
// diagonalize h
// Note: we only need the first n eigenvectors of the (2n x 2n) matrix
valarray<double> w(h.m());
// q is (2n,2n)
ComplexMatrix q(h.context(),h.n(),h.n(),h.nb(),h.nb());
h.heev('l',w,q);
// compute the first n eigenvectors and store in wf
// Y = Z Q21 (store result in dwf)
// get Q21 in a
a.getsub(q,a.n(),a.n(),a.n(),0);
cp.gemm('n','n',1.0,c,a,0.0);
// Y = Y Q11 (store result in wf)
// get Q11 in a
a.getsub(q,a.n(),a.n(),0,0);
c.gemm('n','n',1.0,ct,a,0.0);
// add two contributions
c += cp;
// wf now contains the corrected eigenvectors
} // if real
} // ikp
} // ispin
}
......@@ -15,7 +15,7 @@
// JDWavefunctionStepper.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: JDWavefunctionStepper.h,v 1.1 2009-09-07 19:24:41 fgygi Exp $
// $Id: JDWavefunctionStepper.h,v 1.2 2009-09-08 05:36:49 fgygi Exp $
#ifndef JDWAVEFUNCTIONSTEPPER_H
#define JDWAVEFUNCTIONSTEPPER_H
......@@ -38,7 +38,7 @@ class JDWavefunctionStepper : public WavefunctionStepper
void update(Wavefunction& dwf);
virtual void preprocess(void) {}
JDWavefunctionStepper(Wavefunction& wf, Preconditioner& p,
JDWavefunctionStepper(Wavefunction& wf, Preconditioner& p,
EnergyFunctional& ef, TimerMap& tmap);
~JDWavefunctionStepper() {};
};
......
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