jiri/LA_library
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

davidson: added root targeting

Browse Source
This commit is contained in:
Jiri Pittner
2025-12-09 16:03:04 +01:00
parent 3d284d544a
commit 1cb536421f
2 changed files with 46 additions and 3 deletions
Download Patch File Download Diff File Expand all files Collapse all files

22
davidson.h
View File

@@ -33,7 +33,7 @@ namespace LA {
//Note that for efficiency in a direct CI case the diagonalof() should cache its result
//works for complex hermitian-only too, but does not converge sowell for more than 1 root
//works for complex hermitian-only too, but does not converge so well for more than 1 root
//@@@ for large krylov spaces >200 it can occur 'convergence problem in sygv/syev in diagonalize()'
//@@@options: left eigenvectors by matrix transpose, overridesymmetric (for nrmat)
//@@@small matrix gdiagonalize - shift complex roots up (option to gdiagonalize?)
@@ -48,7 +48,7 @@ template <typename T, typename Matrix>
extern void davidson(const Matrix &bigmat, NRVec<T> &eivals, NRVec<T> *eivecs, const char *eivecsfile,
int nroots=1, const bool verbose=0, const double eps=1e-6,
const bool incore=1, int maxit=100, const int maxkrylov = 500,
void (*initguess)(NRVec<T> &)=NULL)
void (*initguess)(NRVec<T> &)=NULL, const typename LA_traits<T>::normtype *target=NULL)
{
bool flag=0;
int n=bigmat.nrows();
@@ -155,6 +155,24 @@ else
for(int i=0; i<=krylovsize; ++i) {r[i]/=LA_traits<T>::realpart(beta[i]); ri[i]/=LA_traits<T>::realpart(beta[i]);}
}
if(target) //resort eigenpairs by distance from the target
{
NRVec<typename LA_traits<T>::normtype> key(krylovsize+1);
for(int i=0; i<=krylovsize; ++i) key[i] = abs(r[i] - *target);
NRPerm<int> p(krylovsize+1);
key.sort(0,p);
NRVec<typename LA_traits<T>::normtype> rp(maxkrylov);
NRMat<T> smallVp(maxkrylov,maxkrylov);
for(int i=0; i<=krylovsize; ++i)
{
rp[i]= r[p[i+1]-1];
for(int j=0; j<=krylovsize; ++j) smallVp(j,i) = smallV(j,p[i+1]-1);
}
r = rp;
smallV = smallVp;
}
typename LA_traits<T>::normtype eival_n=r[nroot];
oldnroot=nroot;
typename LA_traits<T>::normtype test=std::abs(smallV(krylovsize,nroot));

27
t.cc
View File

@@ -4500,7 +4500,7 @@ cout <<"Error = "<<(xt-y).norm()<<endl;
}
if(1)
if(0)
{
//n must not be too small
int n,m;
@@ -4537,4 +4537,29 @@ for(int i=0; i<m; ++i)
}
if(1)
{
int n,m;
cin>>n >>m;
NRSMat<double> a(n);
NRVec<double> rr(n);
for(int i=0;i<n;++i) for(int j=0;j<=i;++j)
{
a(i,j)= RANDDOUBLE();
if(i==j) a(i,i)+= .5*(i-n*.5);
}
NRSMat<double> aa;
NRMat<double> vv(n,n);
aa=a; diagonalize(aa,rr,&vv);
NRVec<double> r(m);
NRVec<double> *eivecs = new NRVec<double>[m];
double target= 0;
cout <<"Exact energies " <<rr<<endl;
davidson(a,r,eivecs,NULL,m,1,1e-6,1,200,300,(void (*)(LA::NRVec<double>&))NULL,&target);
cout <<"Davidson energies " <<r;
}
}//main