mirror of
https://github.com/paboyle/Grid.git
synced 2024-11-10 07:55:35 +00:00
EvenOdd schur decomposed mpcdagmpc version of rhmc determinant.
dH is also small and plaquette looks right.
This commit is contained in:
Peter Boyle
parent
aef98b7226
commit
570150f1d3
185
lib/qcd/action/pseudofermion/OneFlavourEvenOddRational.h
Normal file
185
lib/qcd/action/pseudofermion/OneFlavourEvenOddRational.h
Normal file
@ -0,0 +1,185 @@
|
||||
#ifndef QCD_PSEUDOFERMION_ONE_FLAVOUR_EVEN_ODD_RATIONAL_H
|
||||
#define QCD_PSEUDOFERMION_ONE_FLAVOUR_EVEN_ODD_RATIONAL_H
|
||||
|
||||
namespace Grid{
|
||||
namespace QCD{
|
||||
|
||||
///////////////////////////////////////
|
||||
// One flavour rational
|
||||
///////////////////////////////////////
|
||||
|
||||
// S_f = chi^dag * N(Mpc^dag*Mpc)/D(Mpc^dag*Mpc) * chi
|
||||
//
|
||||
// Here, M is some operator
|
||||
// N and D makeup the rat. poly
|
||||
//
|
||||
|
||||
template<class Impl>
|
||||
class OneFlavourEvenOddRationalPseudoFermionAction : public Action<typename Impl::GaugeField> {
|
||||
public:
|
||||
INHERIT_IMPL_TYPES(Impl);
|
||||
|
||||
typedef OneFlavourRationalParams Params;
|
||||
Params param;
|
||||
|
||||
MultiShiftFunction PowerHalf ;
|
||||
MultiShiftFunction PowerNegHalf;
|
||||
MultiShiftFunction PowerQuarter;
|
||||
MultiShiftFunction PowerNegQuarter;
|
||||
|
||||
private:
|
||||
|
||||
FermionOperator<Impl> & FermOp;// the basic operator
|
||||
|
||||
// NOT using "Nroots"; IroIro is -- perhaps later, but this wasn't good for us historically
|
||||
// and hasenbusch works better
|
||||
|
||||
FermionField PhiEven; // the pseudo fermion field for this trajectory
|
||||
FermionField PhiOdd; // the pseudo fermion field for this trajectory
|
||||
|
||||
|
||||
public:
|
||||
|
||||
OneFlavourEvenOddRationalPseudoFermionAction(FermionOperator<Impl> &Op,
|
||||
Params & p ) : FermOp(Op),
|
||||
PhiEven(Op.FermionRedBlackGrid()),
|
||||
PhiOdd (Op.FermionRedBlackGrid()),
|
||||
param(p)
|
||||
{
|
||||
AlgRemez remez(param.lo,param.hi,param.precision);
|
||||
|
||||
// MdagM^(+- 1/2)
|
||||
std::cout<<GridLogMessage << "Generating degree "<<param.degree<<" for x^(1/2)"<<std::endl;
|
||||
remez.generateApprox(param.degree,1,2);
|
||||
PowerHalf.Init(remez,param.tolerance,false);
|
||||
PowerNegHalf.Init(remez,param.tolerance,true);
|
||||
|
||||
// MdagM^(+- 1/4)
|
||||
std::cout<<GridLogMessage << "Generating degree "<<param.degree<<" for x^(1/4)"<<std::endl;
|
||||
remez.generateApprox(param.degree,1,4);
|
||||
PowerQuarter.Init(remez,param.tolerance,false);
|
||||
PowerNegQuarter.Init(remez,param.tolerance,true);
|
||||
};
|
||||
|
||||
virtual void init(const GaugeField &U, GridParallelRNG& pRNG) {
|
||||
|
||||
// P(phi) = e^{- phi^dag (MpcdagMpc)^-1/2 phi}
|
||||
// = e^{- phi^dag (MpcdagMpc)^-1/4 (MpcdagMpc)^-1/4 phi}
|
||||
// Phi = MpcdagMpc^{1/4} eta
|
||||
//
|
||||
// P(eta) = e^{- eta^dag eta}
|
||||
//
|
||||
// e^{x^2/2 sig^2} => sig^2 = 0.5.
|
||||
//
|
||||
// So eta should be of width sig = 1/sqrt(2).
|
||||
|
||||
RealD scale = std::sqrt(0.5);
|
||||
|
||||
FermionField eta (FermOp.FermionGrid());
|
||||
FermionField etaOdd (FermOp.FermionRedBlackGrid());
|
||||
FermionField etaEven(FermOp.FermionRedBlackGrid());
|
||||
|
||||
gaussian(pRNG,eta); eta=eta*scale;
|
||||
|
||||
pickCheckerboard(Even,etaEven,eta);
|
||||
pickCheckerboard(Odd,etaOdd,eta);
|
||||
|
||||
FermOp.ImportGauge(U);
|
||||
|
||||
// mutishift CG
|
||||
SchurDifferentiableOperator<Impl> Mpc(FermOp);
|
||||
ConjugateGradientMultiShift<FermionField> msCG(param.MaxIter,PowerQuarter);
|
||||
msCG(Mpc,etaOdd,PhiOdd);
|
||||
|
||||
//////////////////////////////////////////////////////
|
||||
// FIXME : Clover term not yet..
|
||||
//////////////////////////////////////////////////////
|
||||
|
||||
assert(FermOp.ConstEE() == 1);
|
||||
PhiEven = zero;
|
||||
|
||||
};
|
||||
|
||||
//////////////////////////////////////////////////////
|
||||
// S = phi^dag (Mdag M)^-1/2 phi
|
||||
//////////////////////////////////////////////////////
|
||||
virtual RealD S(const GaugeField &U) {
|
||||
|
||||
FermOp.ImportGauge(U);
|
||||
|
||||
FermionField Y(FermOp.FermionRedBlackGrid());
|
||||
|
||||
SchurDifferentiableOperator<Impl> Mpc(FermOp);
|
||||
|
||||
ConjugateGradientMultiShift<FermionField> msCG(param.MaxIter,PowerNegQuarter);
|
||||
|
||||
msCG(Mpc,PhiOdd,Y);
|
||||
|
||||
RealD action = norm2(Y);
|
||||
std::cout << GridLogMessage << "Pseudofermion action FIXME -- is -1/4 solve or -1/2 solve faster??? "<<action<<std::endl;
|
||||
|
||||
return action;
|
||||
};
|
||||
|
||||
//////////////////////////////////////////////////////
|
||||
// Need
|
||||
// dS_f/dU = chi^dag d[N/D] chi
|
||||
//
|
||||
// N/D is expressed as partial fraction expansion:
|
||||
//
|
||||
// a0 + \sum_k ak/(M^dagM + bk)
|
||||
//
|
||||
// d[N/D] is then
|
||||
//
|
||||
// \sum_k -ak [M^dagM+bk]^{-1} [ dM^dag M + M^dag dM ] [M^dag M + bk]^{-1}
|
||||
//
|
||||
// Need
|
||||
// Mf Phi_k = [MdagM+bk]^{-1} Phi
|
||||
// Mf Phi = \sum_k ak [MdagM+bk]^{-1} Phi
|
||||
//
|
||||
// With these building blocks
|
||||
//
|
||||
// dS/dU = \sum_k -ak Mf Phi_k^dag [ dM^dag M + M^dag dM ] Mf Phi_k
|
||||
// S = innerprodReal(Phi,Mf Phi);
|
||||
//////////////////////////////////////////////////////
|
||||
virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
|
||||
|
||||
const int Npole = PowerNegHalf.poles.size();
|
||||
|
||||
std::vector<FermionField> MPhi_k (Npole,FermOp.FermionRedBlackGrid());
|
||||
|
||||
FermionField X(FermOp.FermionRedBlackGrid());
|
||||
FermionField Y(FermOp.FermionRedBlackGrid());
|
||||
|
||||
GaugeField tmp(FermOp.GaugeGrid());
|
||||
|
||||
FermOp.ImportGauge(U);
|
||||
|
||||
SchurDifferentiableOperator<Impl> Mpc(FermOp);
|
||||
|
||||
ConjugateGradientMultiShift<FermionField> msCG(param.MaxIter,PowerNegHalf);
|
||||
|
||||
msCG(Mpc,PhiOdd,MPhi_k);
|
||||
|
||||
dSdU = zero;
|
||||
for(int k=0;k<Npole;k++){
|
||||
|
||||
RealD ak = PowerNegHalf.residues[k];
|
||||
|
||||
X = MPhi_k[k];
|
||||
|
||||
Mpc.Mpc(X,Y);
|
||||
Mpc.MpcDeriv (tmp , Y, X ); dSdU=dSdU+ak*tmp;
|
||||
Mpc.MpcDagDeriv(tmp , X, Y ); dSdU=dSdU+ak*tmp;
|
||||
|
||||
}
|
||||
|
||||
dSdU = Ta(dSdU);
|
||||
|
||||
};
|
||||
};
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
#endif
|
||||
55
tests/Test_rhmc_EOWilson1p1.cc
Normal file
55
tests/Test_rhmc_EOWilson1p1.cc
Normal file
@ -0,0 +1,55 @@
|
||||
#include "Grid.h"
|
||||
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
std::vector<int> latt_size = GridDefaultLatt();
|
||||
std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
|
||||
std::vector<int> mpi_layout = GridDefaultMpi();
|
||||
|
||||
GridCartesian Fine(latt_size,simd_layout,mpi_layout);
|
||||
GridRedBlackCartesian RBFine(latt_size,simd_layout,mpi_layout);
|
||||
GridParallelRNG pRNG(&Fine);
|
||||
pRNG.SeedRandomDevice();
|
||||
LatticeLorentzColourMatrix U(&Fine);
|
||||
|
||||
SU3::HotConfiguration(pRNG, U);
|
||||
|
||||
// simplify template declaration? Strip the lorentz from the second template
|
||||
WilsonGaugeAction<LatticeLorentzColourMatrix, LatticeColourMatrix> Waction(5.6);
|
||||
|
||||
Real mass=-0.77;
|
||||
WilsonFermionR FermOp(U,Fine,RBFine,mass);
|
||||
|
||||
// 1+1 flavour
|
||||
OneFlavourRationalParams Params(1.0e-4,64.0,1000,1.0e-6);
|
||||
OneFlavourEvenOddRationalPseudoFermionAction<WilsonImplR> WilsonNf1a(FermOp,Params);
|
||||
OneFlavourEvenOddRationalPseudoFermionAction<WilsonImplR> WilsonNf1b(FermOp,Params);
|
||||
|
||||
//Collect actions
|
||||
ActionLevel Level1;
|
||||
Level1.push_back(&WilsonNf1a);
|
||||
Level1.push_back(&WilsonNf1b);
|
||||
Level1.push_back(&Waction);
|
||||
|
||||
ActionSet FullSet;
|
||||
FullSet.push_back(Level1);
|
||||
|
||||
// Create integrator
|
||||
typedef MinimumNorm2 IntegratorAlgorithm;// change here to change the algorithm
|
||||
|
||||
IntegratorParameters MDpar(12,20,1.0);
|
||||
Integrator<IntegratorAlgorithm> MDynamics(&Fine,MDpar, FullSet);
|
||||
|
||||
// Create HMC
|
||||
HMCparameters HMCpar;
|
||||
HybridMonteCarlo<IntegratorAlgorithm> HMC(HMCpar, MDynamics);
|
||||
HMC.evolve(U);
|
||||
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user