Commit 47b3418a by Francois Gygi

skip diagonalization if nst==0

git-svn-id: http://qboxcode.org/svn/qb/trunk@1214 cba15fb0-1239-40c8-b417-11db7ca47a34
parent bd5a57c3
......@@ -585,96 +585,99 @@ void Wavefunction::diag(Wavefunction& dwf, bool eigvec)
// overwritten
for ( int ispin = 0; ispin < nspin(); ispin++ )
{
for ( int ikp = 0; ikp < nkp(); ikp++ )
if ( nst_[ispin] > 0 )
{
// compute eigenvalues
if ( sd(ispin,ikp)->basis().real() )
for ( int ikp = 0; ikp < nkp(); ikp++ )
{
// proxy real matrices c, cp
DoubleMatrix c(sd(ispin,ikp)->c());
DoubleMatrix cp(dwf.sd(ispin,ikp)->c());
// compute eigenvalues
if ( sd(ispin,ikp)->basis().real() )
{
// proxy real matrices c, cp
DoubleMatrix c(sd(ispin,ikp)->c());
DoubleMatrix cp(dwf.sd(ispin,ikp)->c());
DoubleMatrix h(c.context(),c.n(),c.n(),c.nb(),c.nb());
DoubleMatrix h(c.context(),c.n(),c.n(),c.nb(),c.nb());
// factor 2.0 in next line: G and -G
h.gemm('t','n',2.0,c,cp,0.0);
// rank-1 update correction
h.ger(-1.0,c,0,cp,0);
// factor 2.0 in next line: G and -G
h.gemm('t','n',2.0,c,cp,0.0);
// rank-1 update correction
h.ger(-1.0,c,0,cp,0);
// cout << " Hamiltonian at k = " << sd(ispin,ikp)->kpoint()
// << endl;
// cout << h;
// cout << " Hamiltonian at k = " << sd(ispin,ikp)->kpoint()
// << endl;
// cout << h;
#if 1
valarray<double> w(h.m());
if ( eigvec )
{
valarray<double> w(h.m());
if ( eigvec )
{
DoubleMatrix z(c.context(),c.n(),c.n(),c.nb(),c.nb());
h.syev('l',w,z);
//h.syevx('l',w,z,1.e-6);
cp = c;
c.gemm('n','n',1.0,cp,z,0.0);
}
else
{
h.syev('l',w);
}
#else
vector<double> w(h.m());
DoubleMatrix z(c.context(),c.n(),c.n(),c.nb(),c.nb());
h.syev('l',w,z);
//h.syevx('l',w,z,1.e-6);
cp = c;
c.gemm('n','n',1.0,cp,z,0.0);
const int maxsweep = 30;
int nsweep = jacobi(maxsweep,1.e-6,h,z,w);
if ( eigvec )
{
cp = c;
c.gemm('n','n',1.0,cp,z,0.0);
}
#endif
// set eigenvalues in SlaterDet
sd(ispin,ikp)->set_eig(w);
}
else
{
h.syev('l',w);
}
#else
vector<double> w(h.m());
DoubleMatrix z(c.context(),c.n(),c.n(),c.nb(),c.nb());
const int maxsweep = 30;
int nsweep = jacobi(maxsweep,1.e-6,h,z,w);
if ( eigvec )
{
cp = c;
c.gemm('n','n',1.0,cp,z,0.0);
}
#endif
// set eigenvalues in SlaterDet
sd(ispin,ikp)->set_eig(w);
}
else
{
ComplexMatrix& c(sd(ispin,ikp)->c());
ComplexMatrix& cp(dwf.sd(ispin,ikp)->c());
ComplexMatrix h(c.context(),c.n(),c.n(),c.nb(),c.nb());
h.gemm('c','n',1.0,c,cp,0.0);
// cout << " Hamiltonian at k = "
// << sd(ispin,ikp)->kpoint() << endl;
// cout << h;
valarray<double> w(h.m());
if ( eigvec )
{
ComplexMatrix& c(sd(ispin,ikp)->c());
ComplexMatrix& cp(dwf.sd(ispin,ikp)->c());
ComplexMatrix h(c.context(),c.n(),c.n(),c.nb(),c.nb());
h.gemm('c','n',1.0,c,cp,0.0);
// cout << " Hamiltonian at k = "
// << sd(ispin,ikp)->kpoint() << endl;
// cout << h;
valarray<double> w(h.m());
if ( eigvec )
{
#if DEBUG
ComplexMatrix hcopy(h);
ComplexMatrix hcopy(h);
#endif
ComplexMatrix z(c.context(),c.n(),c.n(),c.nb(),c.nb());
h.heev('l',w,z);
cp = c;
c.gemm('n','n',1.0,cp,z,0.0);
ComplexMatrix z(c.context(),c.n(),c.n(),c.nb(),c.nb());
h.heev('l',w,z);
cp = c;
c.gemm('n','n',1.0,cp,z,0.0);
#if DEBUG
// check that z contains eigenvectors of h
// diagonal matrix with eigenvalues on diagonal
ComplexMatrix d(c.context(),c.n(),c.n(),c.nb(),c.nb());
// the following test works only on one task
assert(ctxt_.size()==1);
for ( int i = 0; i < d.m(); i++ )
d[i+d.n()*i] = w[i];
ComplexMatrix dz(c.context(),c.n(),c.n(),c.nb(),c.nb());
dz.gemm('n','c',1.0,d,z,0.0);
ComplexMatrix zdz(c.context(),c.n(),c.n(),c.nb(),c.nb());
zdz.gemm('n','n',1.0,z,dz,0.0);
// zdz should be equal to hcopy
zdz -= hcopy;
cout << " heev: norm of error: " << zdz.nrm2() << endl;
// check that z contains eigenvectors of h
// diagonal matrix with eigenvalues on diagonal
ComplexMatrix d(c.context(),c.n(),c.n(),c.nb(),c.nb());
// the following test works only on one task
assert(ctxt_.size()==1);
for ( int i = 0; i < d.m(); i++ )
d[i+d.n()*i] = w[i];
ComplexMatrix dz(c.context(),c.n(),c.n(),c.nb(),c.nb());
dz.gemm('n','c',1.0,d,z,0.0);
ComplexMatrix zdz(c.context(),c.n(),c.n(),c.nb(),c.nb());
zdz.gemm('n','n',1.0,z,dz,0.0);
// zdz should be equal to hcopy
zdz -= hcopy;
cout << " heev: norm of error: " << zdz.nrm2() << endl;
#endif
}
else
{
h.heev('l',w);
}
// set eigenvalues in SlaterDet
sd(ispin,ikp)->set_eig(w);
}
else
{
h.heev('l',w);
}
// set eigenvalues in SlaterDet
sd(ispin,ikp)->set_eig(w);
}
}
}
......
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