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

improevd lanczos

Browse Source
This commit is contained in:
Jiri Pittner
2026-02-10 15:18:56 +01:00
parent 07b923379d
commit 36983222e8
2 changed files with 20 additions and 5 deletions
Download Patch File Download Diff File Expand all files Collapse all files

14
lanczos.h
View File

@@ -114,12 +114,22 @@ for(j=1; j<maxkrylov;++j)
{ {
//generate an arbitrary vector and orthonormalize it to all previous v_j //generate an arbitrary vector and orthonormalize it to all previous v_j
vec2.randomize(1.); vec2.randomize(1.);
NRVec<T> vec3(n);
for(int k=0; k<j; ++k) for(int k=0; k<j; ++k)
{ {
T f = vec3.dot(vec2); T f;
if(incore)
{
f = v0[k].dot(vec2);
vec2.axpy(-f,v0[k]);
}
else
{
NRVec<T> vec3(n);
s0->get(vec3,k);
f = vec3.dot(vec2);
vec2.axpy(-f,vec3); vec2.axpy(-f,vec3);
} }
}
vec2.normalize(); vec2.normalize();
} }
//vec2 now stores v_j //vec2 now stores v_j

9
t.cc
View File

@@ -61,6 +61,11 @@ inline int randind(const int n)
complex<double> mycident (const complex<double>&x) {return x;} complex<double> mycident (const complex<double>&x) {return x;}
void randomguess(NRVec<double> &v)
{
v.randomize(1);
}
void printme(const NRPerm<int> &p) void printme(const NRPerm<int> &p)
{ {
PERM_RANK_TYPE rank=p.rank(); PERM_RANK_TYPE rank=p.rank();
@@ -4678,8 +4683,8 @@ NRVec<double> r(m);
NRVec<double> *eivecs = new NRVec<double>[m]; NRVec<double> *eivecs = new NRVec<double>[m];
cout <<"Exact energies " <<rr<<endl; cout <<"Exact energies " <<rr<<endl;
if(which) lanczos(a,r,eivecs,NULL,m,1,1e-6,1,200,300); if(which) lanczos(a,r,eivecs,NULL,m,1,1e-6,1,500,500,randomguess);
else davidson(a,r,eivecs,NULL,m,1,1e-6,1,200,300); else davidson(a,r,eivecs,NULL,m,1,1e-6,1,500,500,randomguess);
cout <<"Iterative energies " <<r; cout <<"Iterative energies " <<r;
cout <<"Eigenvectors compare:\n"; cout <<"Eigenvectors compare:\n";
for(int i=0; i<m; ++i) for(int i=0; i<m; ++i)