portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-06-13 12:47:05 +01:00
Code Releases Activity

Starting reorg of Blocked lanczos

Browse Source
This commit is contained in:
paboyle
2017-10-11 10:12:07 +01:00
parent bf58557fb1
commit 9fe6ac71ea
10 changed files with 5 additions and 4 deletions
Download Patch File Download Diff File Expand all files Collapse all files

153
tests/debug/Test_synthetic_lanczos.cc
View File

@ -1,153 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_synthetic_lanczos.cc
Copyright (C) 2015
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: paboyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <fenv.h>
#include <Grid/Grid.h>
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
static int
FEenableexcept (unsigned int excepts)
{
#if 0
static fenv_t fenv;
unsigned int new_excepts = excepts & FE_ALL_EXCEPT,
old_excepts; // previous masks
if ( fegetenv (&fenv) ) return -1;
old_excepts = fenv.__control & FE_ALL_EXCEPT;
// unmask
fenv.__control &= ~new_excepts;
fenv.__mxcsr &= ~(new_excepts << 7);
return ( fesetenv (&fenv) ? -1 : old_excepts );
#else
return 0;
#endif
}
template<class Field> class DumbOperator : public LinearOperatorBase<Field> {
public:
LatticeComplex scale;
DumbOperator(GridBase *grid) : scale(grid)
{
GridParallelRNG pRNG(grid);
std::vector<int> seeds({5,6,7,8});
pRNG.SeedFixedIntegers(seeds);
random(pRNG,scale);
scale = exp(-Grid::real(scale)*3.0);
std::cout << " True matrix \n"<< scale <<std::endl;
}
// Support for coarsening to a multigrid
void OpDiag (const Field &in, Field &out) {};
void OpDir (const Field &in, Field &out,int dir,int disp){};
void Op (const Field &in, Field &out){
out = scale * in;
}
void AdjOp (const Field &in, Field &out){
out = scale * in;
}
void HermOp(const Field &in, Field &out){
double n1, n2;
HermOpAndNorm(in,out,n1,n2);
}
void HermOpAndNorm(const Field &in, Field &out,double &n1,double &n2){
ComplexD dot;
out = scale * in;
dot= innerProduct(in,out);
n1=real(dot);
dot = innerProduct(out,out);
n2=real(dot);
}
};
int main (int argc, char ** argv)
{
// FEenableexcept(FE_ALL_EXCEPT & ~FE_INEXACT);
Grid_init(&argc,&argv);
GridCartesian *grid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),
GridDefaultSimd(Nd,vComplex::Nsimd()),
GridDefaultMpi());
GridParallelRNG RNG(grid);
std::vector<int> seeds({1,2,3,4});
RNG.SeedFixedIntegers(seeds);
RealD alpha = 1.2;
RealD beta = 0.1;
RealD mu = 0.0;
int order = 11;
ChebyshevLanczos<LatticeComplex> Cheby(alpha,beta,mu,order);
std::ofstream file("cheby.dat");
Cheby.csv(file);
HermOpOperatorFunction<LatticeComplex> X;
DumbOperator<LatticeComplex> HermOp(grid);
const int Nk = 40;
const int Nm = 80;
const int Nit= 10000;
int Nconv;
RealD eresid = 1.0e-6;
ImplicitlyRestartedLanczos<LatticeComplex> IRL(HermOp,X,Nk,Nk,Nm,eresid,Nit);
ImplicitlyRestartedLanczos<LatticeComplex> ChebyIRL(HermOp,Cheby,Nk,Nk,Nm,eresid,Nit);
LatticeComplex src(grid); gaussian(RNG,src);
{
std::vector<RealD> eval(Nm);
std::vector<LatticeComplex> evec(Nm,grid);
IRL.calc(eval,evec,src, Nconv);
}
{
std::vector<RealD> eval(Nm);
std::vector<LatticeComplex> evec(Nm,grid);
ChebyIRL.calc(eval,evec,src, Nconv);
}
Grid_finalize();
}