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

Merge branch 'develop' of https://github.com/paboyle/Grid into feature/block_lanczos

Browse Source
This commit is contained in:
Chulwoo Jung
2020-10-01 15:35:29 -04:00
parent 8ee26f9112 ace9cd64bb
commit 43943008bf
93 changed files with 11630 additions and 937 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
tests/IO/Test_serialisation.cc
View File

@ -47,14 +47,14 @@ public:
bool , b,
std::vector<double>, array,
std::vector<std::vector<double> >, twodimarray,
std::vector<std::vector<std::vector<Complex> > >, cmplx3darray,
std::vector<std::vector<std::vector<std::complex<double>> > >, cmplx3darray,
SpinColourMatrix, scm
);
myclass() {}
myclass(int i)
: array(4,5.1)
, twodimarray(3,std::vector<double>(5, 1.23456))
, cmplx3darray(3,std::vector<std::vector<Complex>>(5, std::vector<Complex>(7, Complex(1.2, 3.4))))
, cmplx3darray(3,std::vector<std::vector<std::complex<double>>>(5, std::vector<std::complex<double>>(7, std::complex<double>(1.2, 3.4))))
, ve(2, myenum::blue)
{
e=myenum::red;
@ -121,11 +121,11 @@ namespace Eigen {
// Perform I/O tests on a range of tensor types
// Test coverage: scalars, complex and GridVectors in single, double and default precision
class TensorIO : public Serializable {
using TestScalar = ComplexD;
using TestScalar = std::complex<double>;
using SR3 = Eigen::Sizes<9,4,2>;
using SR5 = Eigen::Sizes<5,4,3,2,1>;
using ESO = Eigen::StorageOptions;
using TensorRank3 = Eigen::Tensor<ComplexF, 3, ESO::RowMajor>;
using TensorRank3 = Eigen::Tensor<std::complex<float>, 3, ESO::RowMajor>;
using TensorR5 = Eigen::TensorFixedSize<Real, SR5>;
using TensorR5Alt = Eigen::TensorFixedSize<Real, SR5, ESO::RowMajor>;
using Tensor942 = Eigen::TensorFixedSize<TestScalar, SR3, ESO::RowMajor>;
@ -134,8 +134,8 @@ class TensorIO : public Serializable {
using LSCTensor = Eigen::TensorFixedSize<SpinColourMatrix, Eigen::Sizes<6,5>>;
static const Real FlagR;
static const Complex Flag;
static const ComplexF FlagF;
static const std::complex<double> Flag;
static const std::complex<float> FlagF;
static const TestScalar FlagTS;
static const char * const pszFilePrefix;
@ -230,8 +230,8 @@ public:
};
const Real TensorIO::FlagR {1};
const Complex TensorIO::Flag {1,-1};
const ComplexF TensorIO::FlagF {1,-1};
const std::complex<double> TensorIO::Flag {1,-1};
const std::complex<float> TensorIO::FlagF {1,-1};
const TensorIO::TestScalar TensorIO::FlagTS{1,-1};
const char * const TensorIO::pszFilePrefix = "tensor_";

66
tests/Test_simd.cc
View File

@ -101,14 +101,14 @@ public:
// FIXME still to test:
//
// innerProduct,
// norm2,
// norm2,
// Reduce,
//
// mac,mult,sub,add, vone,vzero,vcomplex_i, =Zero(),
// vset,vsplat,vstore,vstream,vload, scalar*vec, vec*scalar
// unary -,
// *= , -=, +=
// outerproduct,
// outerproduct,
// zeroit
// permute
class funcReduce {
@ -119,12 +119,12 @@ template<class reduce,class scal> void sfunc(reduce &rr,scal &i1,scal &i2) con
std::string name(void) const { return std::string("Reduce"); }
};
template<class scal, class vec,class functor >
template<class scal, class vec,class functor >
void Tester(const functor &func)
{
GridSerialRNG sRNG;
sRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
int Nsimd = vec::Nsimd();
ExtractBuffer<scal> input1(Nsimd);
@ -172,6 +172,8 @@ void Tester(const functor &func)
}
if ( ok==0 ) {
std::cout<<GridLogMessage << " OK!" <<std::endl;
} else {
std::cout<<GridLogMessage << " wrong!" <<std::endl;
}
assert(ok==0);
}
@ -229,17 +231,19 @@ void IntTester(const functor &func)
}
if ( ok==0 ) {
std::cout<<GridLogMessage << " OK!" <<std::endl;
} else {
std::cout<<GridLogMessage << " wrong!" <<std::endl;
}
assert(ok==0);
}
template<class reduced,class scal, class vec,class functor >
template<class reduced,class scal, class vec,class functor >
void ReductionTester(const functor &func)
{
GridSerialRNG sRNG;
sRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
int Nsimd = vec::Nsimd();
ExtractBuffer<scal> input1(Nsimd);
@ -278,12 +282,14 @@ void ReductionTester(const functor &func)
}
if ( ok==0 ) {
std::cout<<GridLogMessage << " OK!" <<std::endl;
} else {
std::cout<<GridLogMessage << " wrong!" <<std::endl;
}
assert(ok==0);
}
template<class reduced,class scal, class vec,class functor >
template<class reduced,class scal, class vec,class functor >
void IntReductionTester(const functor &func)
{
int Nsimd = vec::Nsimd();
@ -323,6 +329,8 @@ void IntReductionTester(const functor &func)
}
if ( ok==0 ) {
std::cout<<GridLogMessage << " OK!" <<std::endl;
} else {
std::cout<<GridLogMessage << " wrong!" <<std::endl;
}
assert(ok==0);
}
@ -333,7 +341,7 @@ public:
int n;
funcPermute(int _n) { n=_n;};
template<class vec> void operator()(vec &rr,vec &i1,vec &i2) const { permute(rr,i1,n);}
template<class scal> void apply(ExtractBuffer<scal> &rr,ExtractBuffer<scal> &in) const {
template<class scal> void apply(ExtractBuffer<scal> &rr,ExtractBuffer<scal> &in) const {
int sz=in.size();
int msk = sz>>(n+1);
for(int i=0;i<sz;i++){
@ -351,8 +359,8 @@ public:
template<class scal> void apply(ExtractBuffer<scal> &r1,
ExtractBuffer<scal> &r2,
ExtractBuffer<scal> &in1,
ExtractBuffer<scal> &in2) const
{
ExtractBuffer<scal> &in2) const
{
int sz=in1.size();
int msk = sz>>(n+1);
@ -364,7 +372,7 @@ public:
if ( (i&msk) == 0 ) { r2[i]=in1[j2];}
else { r2[i]=in2[j2];}
}
}
}
std::string name(void) const { return std::string("Exchange"); }
};
@ -374,7 +382,7 @@ public:
int n;
funcRotate(int _n) { n=_n;};
template<class vec> void operator()(vec &rr,vec &i1,vec &i2) const { rr=rotate(i1,n);}
template<class scal> void apply(ExtractBuffer<scal> &rr,ExtractBuffer<scal> &in) const {
template<class scal> void apply(ExtractBuffer<scal> &rr,ExtractBuffer<scal> &in) const {
int sz = in.size();
for(int i=0;i<sz;i++){
rr[i] = in[(i+n)%sz];
@ -384,12 +392,12 @@ public:
};
template<class scal, class vec,class functor >
template<class scal, class vec,class functor >
void PermTester(const functor &func)
{
GridSerialRNG sRNG;
sRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
int Nsimd = vec::Nsimd();
ExtractBuffer<scal> input1(Nsimd);
@ -425,37 +433,39 @@ void PermTester(const functor &func)
for(int i=0;i<Nsimd;i++){
std::cout<< input1[i]<<" ";
}
std::cout <<std::endl;
std::cout <<std::endl;
for(int i=0;i<Nsimd;i++){
std::cout<< result[i]<<" ";
}
std::cout <<std::endl;
std::cout <<std::endl;
for(int i=0;i<Nsimd;i++){
std::cout<< reference[i]<<" ";
}
std::cout <<std::endl;
std::cout <<std::endl;
std::cout<<GridLogMessage<< "*****" << std::endl;
}
for(int i=0;i<Nsimd;i++){
if ( abs(reference[i]-result[i])>1.0e-7){
std::cout<<GridLogMessage<< "*****" << std::endl;
std::cout<<GridLogMessage<< "*****" << std::endl;
std::cout<<GridLogMessage<< "["<<i<<"] "<< abs(reference[i]-result[i]) << " " <<reference[i]<< " " << result[i]<<std::endl;
ok++;
}
}
if ( ok==0 ) {
std::cout<<GridLogMessage << " OK!" <<std::endl;
} else {
std::cout<<GridLogMessage << " wrong!" <<std::endl;
}
assert(ok==0);
}
template<class scal, class vec,class functor >
template<class scal, class vec,class functor >
void ExchangeTester(const functor &func)
{
GridSerialRNG sRNG;
sRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
int Nsimd = vec::Nsimd();
ExtractBuffer<scal> input1(Nsimd);
@ -566,7 +576,7 @@ int main (int argc, char ** argv)
std::cout << " Test {1,2,3,4} " << Test <<std::endl;
}
*/
GridCartesian Grid(latt_size,simd_layout,mpi_layout);
std::vector<int> seeds({1,2,3,4});
@ -742,7 +752,7 @@ int main (int argc, char ** argv)
for(int r=0;r<vComplexD::Nsimd();r++){
PermTester<ComplexD,vComplexD>(funcRotate(r));
}
std::cout<<GridLogMessage << "==================================="<< std::endl;
std::cout<<GridLogMessage << "Testing vInteger "<< std::endl;
std::cout<<GridLogMessage << "==================================="<< std::endl;
@ -773,13 +783,13 @@ int main (int argc, char ** argv)
precisionChange(&DD[0],&F[0],Ndp);
std::cout << GridLogMessage<<"Double to single";
for(int i=0;i<Ndp;i++){
// std::cout << "DD["<<i<<"] = "<< DD[i]<<" "<<D[i]<<" "<<DD[i]-D[i] <<std::endl;
// std::cout << "DD["<<i<<"] = "<< DD[i]<<" "<<D[i]<<" "<<DD[i]-D[i] <<std::endl;
DD[i] = DD[i] - D[i];
decltype(innerProduct(DD[0],DD[0])) nrm;
nrm = innerProduct(DD[i],DD[i]);
auto tmp = Reduce(nrm);
// std::cout << tmp << std::endl;
assert( tmp < 1.0e-14 );
assert( tmp < 1.0e-14 );
}
std::cout <<" OK ! "<<std::endl;
@ -788,13 +798,13 @@ int main (int argc, char ** argv)
precisionChange(&H[0],&D[0],Ndp);
precisionChange(&DD[0],&H[0],Ndp);
for(int i=0;i<Ndp;i++){
// std::cout << "DD["<<i<<"] = "<< DD[i]<<" "<<D[i]<<" "<<DD[i]-D[i]<<std::endl;
// std::cout << "DD["<<i<<"] = "<< DD[i]<<" "<<D[i]<<" "<<DD[i]-D[i]<<std::endl;
DD[i] = DD[i] - D[i];
decltype(innerProduct(DD[0],DD[0])) nrm;
nrm = innerProduct(DD[i],DD[i]);
auto tmp = Reduce(nrm);
// std::cout << tmp << std::endl;
assert( tmp < 1.0e-3 );
assert( tmp < 1.0e-3 );
}
std::cout <<" OK ! "<<std::endl;
@ -803,13 +813,13 @@ int main (int argc, char ** argv)
precisionChange(&H[0] ,&F[0],Nsp);
precisionChange(&FF[0],&H[0],Nsp);
for(int i=0;i<Nsp;i++){
// std::cout << "FF["<<i<<"] = "<< FF[i]<<" "<<F[i]<<" "<<FF[i]-F[i]<<std::endl;
// std::cout << "FF["<<i<<"] = "<< FF[i]<<" "<<F[i]<<" "<<FF[i]-F[i]<<std::endl;
FF[i] = FF[i] - F[i];
decltype(innerProduct(FF[0],FF[0])) nrm;
nrm = innerProduct(FF[i],FF[i]);
auto tmp = Reduce(nrm);
// std::cout << tmp << std::endl;
assert( tmp < 1.0e-3 );
assert( tmp < 1.0e-3 );
}
std::cout <<" OK ! "<<std::endl;

94
tests/Test_stencil.cc
View File

@ -103,24 +103,26 @@ int main(int argc, char ** argv) {
Bar = Cshift(Foo,dir,disp);
// Implement a stencil code that should agree with cshift!
for(int i=0;i<Check.Grid()->oSites();i++){
{
autoView( check , Check, AcceleratorWrite);
autoView( foo , Foo, AcceleratorRead);
autoView(st_v ,myStencil,AcceleratorRead);
auto CBp=myStencil.CommBuf();
accelerator_for(i,Check.Grid()->oSites(), 1, {
int permute_type;
StencilEntry *SE;
SE = myStencil.GetEntry(permute_type,0,i);
autoView( check , Check, CpuWrite);
autoView( foo , Foo, CpuRead);
if ( SE->_is_local && SE->_permute )
permute(check[i],foo[SE->_offset],permute_type);
else if (SE->_is_local)
check[i] = foo[SE->_offset];
else {
check[i] = myStencil.CommBuf()[SE->_offset];
// std::cout << " receive "<<i<<" " << Check[i]<<std::endl;
// std::cout << " Foo "<<i<<" " << Foo[i]<<std::endl;
}
}
int permute_type;
StencilEntry *SE;
SE = st_v.GetEntry(permute_type,0,i);
if ( SE->_is_local && SE->_permute )
permute(check[i],foo[SE->_offset],permute_type);
else if (SE->_is_local)
check[i] = foo[SE->_offset];
else {
check[i] = CBp[SE->_offset];
}
});
}
Real nrmC = norm2(Check);
Real nrmB = norm2(Bar);
@ -204,36 +206,42 @@ int main(int argc, char ** argv) {
// Implement a stencil code that should agree with that darn cshift!
EStencil.HaloExchange(EFoo,compress);
for(int i=0;i<OCheck.Grid()->oSites();i++){
int permute_type;
StencilEntry *SE;
SE = EStencil.GetEntry(permute_type,0,i);
// std::cout << "Even source "<< i<<" -> " <<SE->_offset << " "<< SE->_is_local<<std::endl;
{
autoView( ocheck , OCheck, AcceleratorWrite);
autoView( efoo , EFoo, AcceleratorRead);
autoView( Est , EStencil,AcceleratorRead);
auto ECBp = EStencil.CommBuf();
accelerator_for(i,OCheck.Grid()->oSites(),1,{
int permute_type;
StencilEntry *SE;
SE = Est.GetEntry(permute_type,0,i);
autoView( ocheck , OCheck, CpuWrite);
autoView( efoo , EFoo, CpuRead);
if ( SE->_is_local && SE->_permute )
permute(ocheck[i],efoo[SE->_offset],permute_type);
else if (SE->_is_local)
ocheck[i] = efoo[SE->_offset];
else
ocheck[i] = EStencil.CommBuf()[SE->_offset];
if ( SE->_is_local && SE->_permute )
permute(ocheck[i],efoo[SE->_offset],permute_type);
else if (SE->_is_local)
ocheck[i] = efoo[SE->_offset];
else
ocheck[i] = ECBp[SE->_offset];
});
}
OStencil.HaloExchange(OFoo,compress);
for(int i=0;i<ECheck.Grid()->oSites();i++){
int permute_type;
StencilEntry *SE;
SE = OStencil.GetEntry(permute_type,0,i);
// std::cout << "ODD source "<< i<<" -> " <<SE->_offset << " "<< SE->_is_local<<std::endl;
{
autoView( echeck , ECheck, AcceleratorWrite);
autoView( ofoo , OFoo, AcceleratorRead);
autoView( Ost , OStencil,AcceleratorRead);
auto OCBp = OStencil.CommBuf();
accelerator_for(i,ECheck.Grid()->oSites(),1,{
int permute_type;
StencilEntry *SE;
SE = Ost.GetEntry(permute_type,0,i);
autoView( echeck , ECheck, CpuWrite);
autoView( ofoo , OFoo, CpuRead);
if ( SE->_is_local && SE->_permute )
permute(echeck[i],ofoo[SE->_offset],permute_type);
else if (SE->_is_local)
echeck[i] = ofoo[SE->_offset];
else
echeck[i] = OStencil.CommBuf()[SE->_offset];
if ( SE->_is_local && SE->_permute )
permute(echeck[i],ofoo[SE->_offset],permute_type);
else if (SE->_is_local)
echeck[i] = ofoo[SE->_offset];
else
echeck[i] = OCBp[SE->_offset];
});
}
setCheckerboard(Check,ECheck);

2
tests/core/Test_main.cc
View File

@ -137,7 +137,6 @@ int main(int argc, char **argv) {
LatticeReal iscalar(&Fine);
SpinMatrix GammaFive;
iSpinMatrix<vComplex> iGammaFive;
ColourMatrix cmat;
random(FineRNG, Foo);
@ -283,7 +282,6 @@ int main(int argc, char **argv) {
cMat = mydouble * cMat;
sMat = adj(sMat); // LatticeSpinMatrix adjoint
sMat = iGammaFive * sMat; // SpinMatrix * LatticeSpinMatrix
sMat = GammaFive * sMat; // SpinMatrix * LatticeSpinMatrix
scMat = adj(scMat);
cMat = adj(cMat);

8
tests/core/Test_staggered_naive.cc
View File

@ -261,11 +261,11 @@ int main (int argc, char ** argv)
pickCheckerboard(Odd ,phi_o,phi);
SchurDiagMooeeOperator<NaiveStaggeredFermionR,FermionField> HermOpEO(Ds);
HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2);
HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2);
HermOpEO.MpcDagMpc(chi_e,dchi_e);
HermOpEO.MpcDagMpc(chi_o,dchi_o);
HermOpEO.MpcDagMpc(phi_e,dphi_e,t1,t2);
HermOpEO.MpcDagMpc(phi_o,dphi_o,t1,t2);
HermOpEO.MpcDagMpc(phi_e,dphi_e);
HermOpEO.MpcDagMpc(phi_o,dphi_o);
pDce = innerProduct(phi_e,dchi_e);
pDco = innerProduct(phi_o,dchi_o);

80
tests/core/Test_where.cc Normal file
View File

@ -0,0 +1,80 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_poisson_fft.cc
Copyright (C) 2015
Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
Author: Peter Boyle <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 <Grid/Grid.h>
using namespace Grid;
;
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
int threads = GridThread::GetThreads();
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
int N=16;
std::vector<int> latt_size ({N,4,4});
std::vector<int> simd_layout({vComplexD::Nsimd(),1,1});
std::vector<int> mpi_layout ({1,1,1});
int vol = 1;
int nd = latt_size.size();
for(int d=0;d<nd;d++){
vol = vol * latt_size[d];
}
GridCartesian GRID(latt_size,simd_layout,mpi_layout);
LatticeComplexD zz(&GRID);
LatticeInteger coor(&GRID);
LatticeComplexD rn(&GRID);
LatticeComplexD sl(&GRID);
zz = ComplexD(0.0,0.0);
GridParallelRNG RNG(&GRID);
RNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
gaussian(RNG,rn);
RealD nn=norm2(rn);
for(int mu=0;mu<nd;mu++){
RealD ns=0.0;
for(int t=0;t<latt_size[mu];t++){
LatticeCoordinate(coor,mu);
sl=where(coor==Integer(t),rn,zz);
std::cout <<GridLogMessage<< " sl " << sl<<std::endl;
std::cout <<GridLogMessage<<" slice "<<t<<" " << norm2(sl)<<std::endl;
ns=ns+norm2(sl);
}
std::cout <<GridLogMessage <<" sliceNorm" <<mu<<" "<< nn <<" "<<ns<<" " << nn-ns<<std::endl;
}
Grid_finalize();
}

4
tests/debug/Test_cayley_coarsen_support.cc
View File

@ -122,7 +122,7 @@ int main (int argc, char ** argv)
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
Subspace Aggregates(Coarse5d,FGrid,cb);
Aggregates.CreateSubspaceRandom(RNG5);
// Aggregates.CreateSubspaceRandom(RNG5);
subspace=Aggregates.subspace;
@ -163,7 +163,7 @@ int main (int argc, char ** argv)
LittleDiracOp.M(c_src,c_res);
std::cout<<GridLogMessage<<"Testing hermiticity explicitly by inspecting matrix elements"<<std::endl;
LittleDiracOp.AssertHermitian();
// LittleDiracOp.AssertHermitian();
std::cout<<GridLogMessage << "Testing Hermiticity stochastically "<< std::endl;
CoarseVector phi(Coarse5d);

2
tests/lanczos/Test_dwf_compressed_lanczos_reorg_synthetic.cc
View File

@ -311,7 +311,7 @@ int main (int argc, char ** argv) {
std::cout << GridLogMessage << "Keep " << coarse.Nk << " total vectors" << std::endl;
assert(Nm2 >= Nm1);
const int nbasis= 70;
const int nbasis= 32;
CoarseFineIRL<vSpinColourVector,vTComplex,nbasis> IRL(FrbGrid,CoarseGrid5rb,HermOp,Odd);
std::cout << GridLogMessage << "Constructed CoarseFine IRL" << std::endl;

420
tests/solver/Test_dwf_hdcr_2level.cc Normal file
View File

@ -0,0 +1,420 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_dwf_hdcr.cc
Copyright (C) 2015
Author: Antonin Portelli <antonin.portelli@me.com>
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 <Grid/Grid.h>
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h>
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
using namespace std;
using namespace Grid;
/* Params
* Grid:
* block1(4)
* block2(4)
*
* Subspace
* * Fine : Subspace(nbasis,hi,lo,order,first,step) -- 32, 60,0.02,500,100,100
* * Coarse: Subspace(nbasis,hi,lo,order,first,step) -- 32, 18,0.02,500,100,100
* Smoother:
* * Fine: Cheby(hi, lo, order) -- 60,0.5,10
* * Coarse: Cheby(hi, lo, order) -- 12,0.1,4
* Lanczos:
* CoarseCoarse IRL( Nk, Nm, Nstop, poly(lo,hi,order)) 24,36,24,0.002,4.0,61
*/
RealD InverseApproximation(RealD x){
return 1.0/x;
}
template<class Field,class Matrix> class ChebyshevSmoother : public LinearFunction<Field>
{
public:
typedef LinearOperatorBase<Field> FineOperator;
Matrix & _SmootherMatrix;
FineOperator & _SmootherOperator;
Chebyshev<Field> Cheby;
ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator,Matrix &SmootherMatrix) :
_SmootherOperator(SmootherOperator),
_SmootherMatrix(SmootherMatrix),
Cheby(_lo,_hi,_ord,InverseApproximation)
{};
void operator() (const Field &in, Field &out)
{
Field tmp(in.Grid());
MdagMLinearOperator<Matrix,Field> MdagMOp(_SmootherMatrix);
_SmootherOperator.AdjOp(in,tmp);
Cheby(MdagMOp,tmp,out);
}
};
template<class Field,class Matrix> class MirsSmoother : public LinearFunction<Field>
{
public:
typedef LinearOperatorBase<Field> FineOperator;
Matrix & SmootherMatrix;
FineOperator & SmootherOperator;
RealD tol;
RealD shift;
int maxit;
MirsSmoother(RealD _shift,RealD _tol,int _maxit,FineOperator &_SmootherOperator,Matrix &_SmootherMatrix) :
shift(_shift),tol(_tol),maxit(_maxit),
SmootherOperator(_SmootherOperator),
SmootherMatrix(_SmootherMatrix)
{};
void operator() (const Field &in, Field &out)
{
ZeroGuesser<Field> Guess;
ConjugateGradient<Field> CG(tol,maxit,false);
Field src(in.Grid());
ShiftedMdagMLinearOperator<SparseMatrixBase<Field>,Field> MdagMOp(SmootherMatrix,shift);
SmootherOperator.AdjOp(in,src);
Guess(src,out);
CG(MdagMOp,src,out);
}
};
template<class Field,class Matrix> class RedBlackSmoother : public LinearFunction<Field>
{
public:
typedef LinearOperatorBase<Field> FineOperator;
Matrix & SmootherMatrix;
RealD tol;
RealD shift;
int maxit;
RedBlackSmoother(RealD _shift,RealD _tol,int _maxit,Matrix &_SmootherMatrix) :
shift(_shift),tol(_tol),maxit(_maxit),
SmootherMatrix(_SmootherMatrix)
{};
void operator() (const Field &in, Field &out)
{
std::cout << " Red Black Smootheer "<<norm2(in)<<" " <<norm2(out)<<std::endl;
ConjugateGradient<Field> CG(tol,maxit,false);
out =Zero();
SchurRedBlackDiagMooeeSolve<Field> RBSolver(CG);
RBSolver(SmootherMatrix,in,out);
std::cout << " Red Black Smootheer "<<norm2(in)<<" " <<norm2(out)<<std::endl;
}
};
template<class Fobj,class CComplex,int nbasis, class Matrix, class Guesser, class CoarseSolver>
class MultiGridPreconditioner : public LinearFunction< Lattice<Fobj> > {
public:
typedef Aggregation<Fobj,CComplex,nbasis> Aggregates;
typedef CoarsenedMatrix<Fobj,CComplex,nbasis> CoarseOperator;
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseVector CoarseVector;
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseMatrix CoarseMatrix;
typedef typename Aggregation<Fobj,CComplex,nbasis>::FineField FineField;
typedef LinearOperatorBase<FineField> FineOperator;
typedef LinearFunction <FineField> FineSmoother;
Aggregates & _Aggregates;
CoarseOperator & _CoarseOperator;
Matrix & _FineMatrix;
FineOperator & _FineOperator;
Guesser & _Guess;
FineSmoother & _Smoother1;
FineSmoother & _Smoother2;
CoarseSolver & _CoarseSolve;
int level; void Level(int lv) {level = lv; };
#define GridLogLevel std::cout << GridLogMessage <<std::string(level,'\t')<< " Level "<<level <<" "
MultiGridPreconditioner(Aggregates &Agg, CoarseOperator &Coarse,
FineOperator &Fine,Matrix &FineMatrix,
FineSmoother &Smoother1,
FineSmoother &Smoother2,
Guesser &Guess_,
CoarseSolver &CoarseSolve_)
: _Aggregates(Agg),
_CoarseOperator(Coarse),
_FineOperator(Fine),
_FineMatrix(FineMatrix),
_Smoother1(Smoother1),
_Smoother2(Smoother2),
_Guess(Guess_),
_CoarseSolve(CoarseSolve_),
level(1) { }
MultiGridPreconditioner(Aggregates &Agg, CoarseOperator &Coarse,
FineOperator &Fine,Matrix &FineMatrix,
FineSmoother &Smoother,
Guesser &Guess_,
CoarseSolver &CoarseSolve_)
: _Aggregates(Agg),
_CoarseOperator(Coarse),
_FineOperator(Fine),
_FineMatrix(FineMatrix),
_Smoother1(Smoother),
_Smoother2(Smoother),
_Guess(Guess_),
_CoarseSolve(CoarseSolve_),
level(1) { }
virtual void operator()(const FineField &in, FineField & out)
{
CoarseVector Csrc(_CoarseOperator.Grid());
CoarseVector Csol(_CoarseOperator.Grid());
FineField vec1(in.Grid());
FineField vec2(in.Grid());
double t;
// Fine Smoother
t=-usecond();
_Smoother1(in,out);
t+=usecond();
GridLogLevel << "Smoother took "<< t/1000.0<< "ms" <<std::endl;
// Update the residual
_FineOperator.Op(out,vec1); sub(vec1, in ,vec1);
// Fine to Coarse
t=-usecond();
_Aggregates.ProjectToSubspace (Csrc,vec1);
t+=usecond();
GridLogLevel << "Project to coarse took "<< t/1000.0<< "ms" <<std::endl;
// Coarse correction
t=-usecond();
_CoarseSolve(Csrc,Csol);
t+=usecond();
GridLogLevel << "Coarse solve took "<< t/1000.0<< "ms" <<std::endl;
// Coarse to Fine
t=-usecond();
_Aggregates.PromoteFromSubspace(Csol,vec1);
add(out,out,vec1);
t+=usecond();
GridLogLevel << "Promote to this level took "<< t/1000.0<< "ms" <<std::endl;
// Residual
_FineOperator.Op(out,vec1); sub(vec1 ,in , vec1);
// Fine Smoother
t=-usecond();
_Smoother2(vec1,vec2);
t+=usecond();
GridLogLevel << "Smoother took "<< t/1000.0<< "ms" <<std::endl;
add( out,out,vec2);
}
};
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
const int Ls=16;
const int rLs=8;
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
///////////////////////////////////////////////////
// Construct a coarsened grid; utility for this?
///////////////////////////////////////////////////
std::vector<int> block ({2,2,2,2});
const int nbasis= 32;
auto clatt = GridDefaultLatt();
for(int d=0;d<clatt.size();d++){
clatt[d] = clatt[d]/block[d];
}
GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
std::vector<int> seeds4({1,2,3,4});
std::vector<int> seeds5({5,6,7,8});
std::vector<int> cseeds({5,6,7,8});
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
GridParallelRNG CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds);
LatticeFermion src(FGrid); gaussian(RNG5,src);// src=src+g5*src;
LatticeFermion result(FGrid);
LatticeGaugeField Umu(UGrid);
FieldMetaData header;
std::string file("./ckpoint_lat");
NerscIO::readConfiguration(Umu,header,file);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Building g5R5 hermitian DWF operator" <<std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
RealD mass=0.001;
RealD M5=1.8;
DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> CoarseOperator;
typedef CoarseOperator::CoarseVector CoarseVector;
typedef CoarseOperator::siteVector siteVector;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Calling Aggregation class to build subspace" <<std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
MdagMLinearOperator<DomainWallFermionR,LatticeFermion> HermDefOp(Ddwf);
Subspace Aggregates(Coarse5d,FGrid,0);
assert ( (nbasis & 0x1)==0);
{
int nb=nbasis/2;
// Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.05,500,200,100,0.0);// 18s
// rAggregates.CreateSubspaceChebyshev(RNG5,rHermDefOp,nb,60.0,0.05,500,200,150,0.0);// 15.7 23iter
Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.05,500,200,150,0.0);//
// pad out the rAggregates.
// Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.05,500,500,150,0.0);// 19s
// Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.05,500,200,200,0.0); 15.2s
// Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.05,500,500,200,0.0); 16.3s
for(int n=0;n<nb;n++){
G5R5(Aggregates.subspace[n+nb],Aggregates.subspace[n]);
}
LatticeFermion A(FGrid);
LatticeFermion B(FGrid);
for(int n=0;n<nb;n++){
A = Aggregates.subspace[n];
B = Aggregates.subspace[n+nb];
Aggregates.subspace[n] = A+B; // 1+G5 // eigen value of G5R5 is +1
Aggregates.subspace[n+nb]= A-B; // 1-G5 // eigen value of G5R5 is -1
}
}
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Building coarse representation of Indef operator" <<std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> Level1Op;
Gamma5R5HermitianLinearOperator<DomainWallFermionR,LatticeFermion> HermIndefOp(Ddwf);
Level1Op LDOp(*Coarse5d,1); LDOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << " Running Coarse grid Lanczos "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
MdagMLinearOperator<Level1Op,CoarseVector> IRLHermOp(LDOp);
Chebyshev<CoarseVector> IRLCheby(0.002,12.,151);
FunctionHermOp<CoarseVector> IRLOpCheby(IRLCheby,IRLHermOp);
PlainHermOp<CoarseVector> IRLOp (IRLHermOp);
int Nk=48;
int Nm=64;
int Nstop=48;
int Nconv;
ImplicitlyRestartedLanczos<CoarseVector> IRL(IRLOpCheby,IRLOp,Nstop,Nk,Nm,1.0e-3,20);
std::vector<RealD> eval(Nm);
std::vector<CoarseVector> evec(Nm,Coarse5d);
CoarseVector c_src(Coarse5d);
gaussian(CRNG,c_src);
IRL.calc(eval,evec,c_src,Nconv);
// ConjugateGradient<CoarseVector> CoarseCG(0.01,1000);
ConjugateGradient<CoarseVector> CoarseCG(0.02,1000);// 14.7s
DeflatedGuesser<CoarseVector> DeflCoarseGuesser(evec,eval);
NormalEquations<CoarseVector> DeflCoarseCGNE(LDOp,CoarseCG,DeflCoarseGuesser);
c_src=1.0;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << " Fine PowerMethod "<< std::endl;
PowerMethod<LatticeFermion> PM; PM(HermDefOp,src);
std::cout<<GridLogMessage << " Coarse PowerMethod "<< std::endl;
MdagMLinearOperator<CoarseOperator,CoarseVector> PosdefLdop(LDOp);
PowerMethod<CoarseVector> cPM; cPM(PosdefLdop,c_src);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Building 2 level Multigrid "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
typedef MultiGridPreconditioner<vSpinColourVector, vTComplex,nbasis, DomainWallFermionR,DeflatedGuesser<CoarseVector> , NormalEquations<CoarseVector> > TwoLevelMG;
// MultiGrid preconditioner acting on the coarse space <-> coarsecoarse space
// ChebyshevSmoother<LatticeFermion,DomainWallFermionR> FineSmoother(0.5,60.0,14,HermIndefOp,Ddwf); // 72 iter 63s
// ChebyshevSmoother<LatticeFermion,DomainWallFermionR> FineSmoother(0.1,60.0,20,HermIndefOp,Ddwf); // 66 iter 69s
// ChebyshevSmoother<LatticeFermion,DomainWallFermionR> FineSmoother(0.5,60.0,20,HermIndefOp,Ddwf); // 63 iter 65 s
// ChebyshevSmoother<LatticeFermion,DomainWallFermionR> FineSmoother(1.0,60.0,20,HermIndefOp,Ddwf); // 69, 70
// ChebyshevSmoother<LatticeFermion,DomainWallFermionR> FineSmoother(1.0,60.0,14,HermIndefOp,Ddwf); // 77
// ChebyshevSmoother<LatticeFermion,DomainWallFermionR> FineSmoother(0.5,60.0,10,HermIndefOp,Ddwf); // 23 iter 15.9s
// ChebyshevSmoother<LatticeFermion,DomainWallFermionR> FineSmoother(0.5,60.0,14,HermIndefOp,Ddwf); // 20, 16.9s
ChebyshevSmoother<LatticeFermion,DomainWallFermionR> FineSmoother(0.5,60.0,12,HermIndefOp,Ddwf); // 21, 15.6s
// MirsSmoother<LatticeFermion,DomainWallFermionR> FineCGSmoother(0.05,0.01,20,HermIndefOp,Ddwf);
// RedBlackSmoother<LatticeFermion,DomainWallFermionR> FineRBSmoother(0.00,0.001,100,Ddwf);
// Wrap the 2nd level solver in a MultiGrid preconditioner acting on the fine space
ZeroGuesser<CoarseVector> CoarseZeroGuesser;
TwoLevelMG TwoLevelPrecon(Aggregates, LDOp,
HermIndefOp,Ddwf,
FineSmoother,
DeflCoarseGuesser,
DeflCoarseCGNE);
TwoLevelPrecon.Level(1);
// Apply the fine-coarse-coarsecoarse 2 deep MG preconditioner in an outer PGCR on the fine fgrid
PrecGeneralisedConjugateResidual<LatticeFermion> l1PGCR(1.0e-8,1000,HermIndefOp,TwoLevelPrecon,16,16);
l1PGCR.Level(1);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Calling 2 level Multigrid "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
result=Zero();
l1PGCR(src,result);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Fine CG prec DiagMooee "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
ConjugateGradient<LatticeFermion> FineCG(1.0e-8,10000);
SchurDiagMooeeOperator<DomainWallFermionR,LatticeFermion> FineDiagMooee(Ddwf); // M_ee - Meo Moo^-1 Moe
LatticeFermion f_src_e(FrbGrid); f_src_e=1.0;
LatticeFermion f_res_e(FrbGrid); f_res_e=Zero();
FineCG(FineDiagMooee,f_src_e,f_res_e);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Done "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
Grid_finalize();
}