portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-06-11 03:46:55 +01:00
Code Releases Activity

Eigen fixes and HDCR work

Browse Source
This commit is contained in:
paboyle
2016-03-30 00:16:02 -07:00
parent f54e0ec9bd
commit 340428a1fe
6 changed files with 54 additions and 10106 deletions
Download Patch File Download Diff File Expand all files Collapse all files

30
lib/algorithms/CoarsenedMatrix.h
View File

@ -152,16 +152,18 @@ namespace Grid {
{
// Run a Lanczos with sloppy convergence
const int Nstop = nn;
const int Nk = nn+10;
const int Np = nn+10;
const int Nk = nn+20;
const int Np = nn+20;
const int Nm = Nk+Np;
const int MaxIt= 10000;
RealD resid = 1.0e-5;
RealD resid = 1.0e-3;
Chebyshev<FineField> Cheb(0.2,5.,11);
Chebyshev<FineField> Cheb(0.5,64.0,21);
ImplicitlyRestartedLanczos<FineField> IRL(hermop,Cheb,Nstop,Nk,Nm,resid,MaxIt);
// IRL.lock = 1;
FineField noise(FineGrid); gaussian(RNG,noise);
FineField tmp(FineGrid);
std::vector<RealD> eval(Nm);
std::vector<FineField> evec(Nm,FineGrid);
@ -172,16 +174,34 @@ namespace Grid {
// pull back nn vectors
for(int b=0;b<nn;b++){
subspace[b] = evec[b];
std::cout << GridLogMessage <<"subspace["<<b<<"] = "<<norm2(subspace[b])<<std::endl;
hermop.Op(subspace[b],tmp);
std::cout<<GridLogMessage << "filtered["<<b<<"] <f|MdagM|f> "<<norm2(tmp)<<std::endl;
noise = tmp - sqrt(eval[b])*subspace[b] ;
std::cout<<GridLogMessage << " lambda_"<<b<<" = "<< eval[b] <<" ; [ M - Lambda ]_"<<b<<" vec_"<<b<<" = " <<norm2(noise)<<std::endl;
noise = tmp + eval[b]*subspace[b] ;
std::cout<<GridLogMessage << " lambda_"<<b<<" = "<< eval[b] <<" ; [ M - Lambda ]_"<<b<<" vec_"<<b<<" = " <<norm2(noise)<<std::endl;
}
Orthogonalise();
for(int b=0;b<nn;b++){
std::cout << GridLogMessage <<"subspace["<<b<<"] = "<<norm2(subspace[b])<<std::endl;
}
}
virtual void CreateSubspace(GridParallelRNG &RNG,LinearOperatorBase<FineField> &hermop,int nn=nbasis) {
RealD scale;
ConjugateGradient<FineField> CG(5.0e-3,10000);
ConjugateGradient<FineField> CG(1.0e-2,10000);
FineField noise(FineGrid);
FineField Mn(FineGrid);

33
lib/algorithms/iterative/EigenSort.h
View File

@ -39,42 +39,33 @@ class SortEigen {
private:
//hacking for testing for now
#if 0
static bool less_lmd(RealD left,RealD right){
return fabs(left) < fabs(right);
}
static bool less_pair(std::pair<RealD,Field>& left,
std::pair<RealD,Field>& right){
return fabs(left.first) < fabs(right.first);
}
#else
private:
static bool less_lmd(RealD left,RealD right){
return left > right;
}
static bool less_pair(std::pair<RealD,Field>& left,
std::pair<RealD,Field>& right){
static bool less_pair(std::pair<RealD,Field const*>& left,
std::pair<RealD,Field const*>& right){
return left.first > (right.first);
}
#endif
public:
void push(DenseVector<RealD>& lmd,
DenseVector<Field>& evec,int N) {
DenseVector<std::pair<RealD, Field> > emod;
typename DenseVector<std::pair<RealD, Field> >::iterator it;
DenseVector<Field>& evec,int N) {
DenseVector<Field> cpy(lmd.size(),evec[0]._grid);
for(int i=0;i<lmd.size();i++) cpy[i] = evec[i];
for(int i=0;i<lmd.size();++i){
emod.push_back(std::pair<RealD,Field>(lmd[i],evec[i]));
}
DenseVector<std::pair<RealD, Field const*> > emod(lmd.size());
for(int i=0;i<lmd.size();++i)
emod[i] = std::pair<RealD,Field const*>(lmd[i],&cpy[i]);
partial_sort(emod.begin(),emod.begin()+N,emod.end(),less_pair);
it=emod.begin();
typename DenseVector<std::pair<RealD, Field const*> >::iterator it = emod.begin();
for(int i=0;i<N;++i){
lmd[i]=it->first;
evec[i]=it->second;
evec[i]=*(it->second);
++it;
}
}

29
lib/algorithms/iterative/ImplicitlyRestartedLanczos.h
View File

@ -637,21 +637,20 @@ until convergence
abort();
converged:
// Sorting
eval.clear();
evec.clear();
for(int i=0; i<Nconv; ++i){
eval.push_back(eval2[Iconv[i]]);
evec.push_back(B[Iconv[i]]);
}
_sort.push(eval,evec,Nconv);
std::cout << "\n Converged\n Summary :\n";
std::cout << " -- Iterations = "<< Nconv << "\n";
std::cout << " -- beta(k) = "<< beta_k << "\n";
std::cout << " -- Nconv = "<< Nconv << "\n";
}
// Sorting
eval.resize(Nconv);
evec.resize(Nconv,grid);
for(int i=0; i<Nconv; ++i){
eval[i] = eval2[Iconv[i]];
evec[i] = B[Iconv[i]];
}
_sort.push(eval,evec,Nconv);
std::cout << "\n Converged\n Summary :\n";
std::cout << " -- Iterations = "<< Nconv << "\n";
std::cout << " -- beta(k) = "<< beta_k << "\n";
std::cout << " -- Nconv = "<< Nconv << "\n";
}
/////////////////////////////////////////////////
// Adapted from Rudy's lanczos factor routine