portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-10 07:55:35 +00:00
Code Releases Activity

Merge remote-tracking branch 'upstream/master'

Browse Source
This commit is contained in:
neo 2015-07-21 17:17:50 +09:00
commit 5fc6af1c77
13 changed files with 494 additions and 45 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
lib/algorithms/CoarsenedMatrix.h
View File

@ -117,15 +117,15 @@ namespace Grid {
}
Orthogonalise();
}
virtual void CreateSubspace(GridParallelRNG &RNG,LinearOperatorBase<FineField> &hermop) {
virtual void CreateSubspace(GridParallelRNG &RNG,LinearOperatorBase<FineField> &hermop,int nn=nbasis) {
RealD scale;
ConjugateGradient<FineField> CG(1.0e-4,10000);
ConjugateGradient<FineField> CG(2.0e-3,10000);
FineField noise(FineGrid);
FineField Mn(FineGrid);
for(int b=0;b<nbasis;b++){
for(int b=0;b<nn;b++){
gaussian(RNG,noise);
scale = std::pow(norm2(noise),-0.5);
@ -133,7 +133,7 @@ namespace Grid {
hermop.Op(noise,Mn); std::cout << "noise ["<<b<<"] <n|MdagM|n> "<<norm2(Mn)<<std::endl;
for(int i=0;i<2;i++){
for(int i=0;i<1;i++){
CG(hermop,noise,subspace[b]);
@ -144,7 +144,8 @@ namespace Grid {
}
hermop.Op(noise,Mn); std::cout << "filtered["<<b<<"] <f|MdagM|f> "<<norm2(Mn)<<std::endl;
subspace[b] = noise;
subspace[b] = noise;
}
Orthogonalise();

41
lib/algorithms/LinearOperator.h
View File

@ -71,6 +71,47 @@ namespace Grid {
}
};
////////////////////////////////////////////////////////////////////
// Construct herm op and shift it for mgrid smoother
////////////////////////////////////////////////////////////////////
template<class Matrix,class Field>
class ShiftedMdagMLinearOperator : public LinearOperatorBase<Field> {
Matrix &_Mat;
RealD _shift;
public:
ShiftedMdagMLinearOperator(Matrix &Mat,RealD shift): _Mat(Mat), _shift(shift){};
// Support for coarsening to a multigrid
void OpDiag (const Field &in, Field &out) {
_Mat.Mdiag(in,out);
assert(0);
}
void OpDir (const Field &in, Field &out,int dir,int disp) {
_Mat.Mdir(in,out,dir,disp);
assert(0);
}
void Op (const Field &in, Field &out){
_Mat.M(in,out);
assert(0);
}
void AdjOp (const Field &in, Field &out){
_Mat.Mdag(in,out);
assert(0);
}
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
_Mat.MdagM(in,out,n1,n2);
out = out + _shift*in;
ComplexD dot;
dot= innerProduct(in,out);
n1=real(dot);
n2=norm2(out);
}
void HermOp(const Field &in, Field &out){
RealD n1,n2;
HermOpAndNorm(in,out,n1,n2);
}
};
////////////////////////////////////////////////////////////////////
// Wrap an already herm matrix
////////////////////////////////////////////////////////////////////

46
lib/algorithms/approx/Chebyshev.h
View File

@ -50,6 +50,15 @@ namespace Grid {
return;
}
// Convenience for plotting the approximation
void PlotApprox(std::ostream &out) {
out<<"Polynomial approx ["<<lo<<","<<hi<<"]"<<std::endl;
for(double x=lo;x<hi;x+=(hi-lo)/50.0){
out <<x<<"\t"<<approx(x)<<std::endl;
}
};
Chebyshev(double _lo,double _hi,int _order, double (* func)(double) ){
lo=_lo;
hi=_hi;
@ -95,46 +104,39 @@ namespace Grid {
return sum;
};
// Convenience for plotting the approximation
void PlotApprox(std::ostream &out) {
out<<"Polynomial approx ["<<lo<<","<<hi<<"]"<<std::endl;
for(double x=lo;x<hi;x+=(hi-lo)/50.0){
out <<x<<"\t"<<approx(x)<<std::endl;
}
};
// Implement the required interface; could require Lattice base class
// Implement the required interface
void operator() (LinearOperatorBase<Field> &Linop, const Field &in, Field &out) {
Field T0 = in;
Field T1 = T0; // Field T1(T0._grid); more efficient but hardwires Lattice class
Field T2 = T1;
GridBase *grid=in._grid;
int vol=grid->gSites();
Field T0(grid); T0 = in;
Field T1(grid);
Field T2(grid);
Field y(grid);
// use a pointer trick to eliminate copies
Field *Tnm = &T0;
Field *Tn = &T1;
Field *Tnp = &T2;
Field y = in;
std::cout << "Chebyshev ["<<lo<<","<<hi<<"]"<< " order "<<order <<std::endl;
// Tn=T1 = (xscale M + mscale)in
double xscale = 2.0/(hi-lo);
double mscale = -(hi+lo)/(hi-lo);
// Tn=T1 = (xscale M + mscale)in
Linop.Op(T0,y);
Linop.HermOp(T0,y);
T1=y*xscale+in*mscale;
// sum = .5 c[0] T0 + c[1] T1
out = (0.5*Coeffs[0])*T0 + Coeffs[1]*T1;
for(int n=2;n<order;n++){
Linop.Op(*Tn,y);
Linop.HermOp(*Tn,y);
y=xscale*y+mscale*(*Tn);
*Tnp=2.0*y-(*Tnm);
out=out+Coeffs[n]* (*Tnp);
// Cycle pointers to avoid copies

2
lib/algorithms/approx/Remez.h
View File

@ -15,7 +15,7 @@
#ifndef INCLUDED_ALG_REMEZ_H
#define INCLUDED_ALG_REMEZ_H
#include <algorithms/approx/bigfloat_double.h>
#include <algorithms/approx/bigfloat.h>
#define JMAX 10000 //Maximum number of iterations of Newton's approximation
#define SUM_MAX 10 // Maximum number of terms in exponential

2
lib/algorithms/iterative/ConjugateGradient.h
View File

@ -15,7 +15,7 @@ public:
Integer MaxIterations;
int verbose;
ConjugateGradient(RealD tol,Integer maxit) : Tolerance(tol), MaxIterations(maxit) {
verbose=1;
verbose=0;
};

2
lib/algorithms/iterative/PrecGeneralisedConjugateResidual.h
View File

@ -137,7 +137,7 @@ namespace Grid {
cp = axpy_norm(r,-a,q[peri_k],r);
std::cout<< " VPCG_step resid" <<sqrt(cp/rsq)<<std::endl;
std::cout<< " VPGCR_step resid" <<sqrt(cp/rsq)<<std::endl;
if((k==nstep-1)||(cp<rsq)){
return cp;
}

2
lib/qcd/action/fermion/WilsonKernels.cc
View File

@ -1,5 +1,4 @@
#include <Grid.h>
namespace Grid {
namespace QCD {
@ -13,7 +12,6 @@ void DiracOptDhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U,
vHalfSpinColourVector Uchi;
int offset,local,perm, ptype;
// Xp
int ss = sF;
offset = st._offsets [Xp][ss];

2
lib/qcd/action/fermion/g5HermitianLinop.h
View File

@ -1,7 +1,9 @@
#ifndef G5_HERMITIAN_LINOP
#define G5_HERMITIAN_LINOP
namespace Grid {
namespace QCD {
////////////////////////////////////////////////////////////////////
// Wrap an already herm matrix
////////////////////////////////////////////////////////////////////

2
lib/qcd/hmc/integrators/Integrator.cc
View File

@ -18,7 +18,7 @@ namespace Grid{
Pmu = zero;
for(int mu=0;mu<Nd;mu++){
SU3::GaussianLieAlgebraMatrix(pRNG, Pmu);
pokeLorentz(P, Pmu, mu);
PokeIndex<LorentzIndex>(P, Pmu, mu);
}
}

6
lib/qcd/hmc/integrators/Integrator.h
View File

@ -77,10 +77,10 @@ namespace Grid{
LatticeColourMatrix Umu(U._grid);
LatticeColourMatrix Pmu(U._grid);
for (int mu = 0; mu < Nd; mu++){
Umu=peekLorentz(U, mu);
Pmu=peekLorentz(*P, mu);
Umu=PeekIndex<LorentzIndex>(U, mu);
Pmu=PeekIndex<LorentzIndex>(*P, mu);
Umu = expMat(Pmu, ep, Params.Nexp)*Umu;
pokeLorentz(U, Umu, mu);
PokeIndex<LorentzIndex>(U, Umu, mu);
}
}

212
lib/qcd/hmc/integrators/Integrator_base.h Normal file
View File

@ -0,0 +1,212 @@
//--------------------------------------------------------------------
/*! @file Integrator_base.h
* @brief Declaration of classes for the abstract Molecular Dynamics integrator
*
* @author Guido Cossu
*/
//--------------------------------------------------------------------
#ifndef INTEGRATOR_INCLUDED
#define INTEGRATOR_INCLUDED
#include <memory>
class Observer;
/*! @brief Abstract base class for Molecular Dynamics management */
namespace Grid{
namespace QCD{
typedef Action<LatticeLorentzColourMatrix>* ActPtr; // now force the same size as the rest of the code
typedef std::vector<ActPtr> ActionLevel;
typedef std::vector<ActionLevel> ActionSet;
typedef std::vector<Observer*> ObserverList;
class LeapFrog;
struct IntegratorParameters{
int Nexp;
int MDsteps; // number of outer steps
RealD trajL; // trajectory length
RealD stepsize;
IntegratorParameters(int Nexp_,
int MDsteps_,
RealD trajL_):
Nexp(Nexp_),MDsteps(MDsteps_),trajL(trajL_),stepsize(trajL/MDsteps){};
};
namespace MDutils{
void generate_momenta(LatticeLorentzColourMatrix&,GridParallelRNG&);
void generate_momenta_su3(LatticeLorentzColourMatrix&,GridParallelRNG&);
}
template< class IntegratorPolicy >
class Integrator{
private:
IntegratorParameters Params;
const ActionSet as;
const std::vector<int> Nrel; //relative step size per level
//ObserverList observers; // not yet
std::unique_ptr<LatticeLorentzColourMatrix> P;
IntegratorPolicy TheIntegrator;// contains parameters too
void update_P(LatticeLorentzColourMatrix&U, int level,double ep){
for(int a=0; a<as[level].size(); ++a){
LatticeLorentzColourMatrix force(U._grid);
as[level].at(a)->deriv(U,force);
*P -= force*ep;
}
}
void update_U(LatticeLorentzColourMatrix&U, double ep){
//rewrite exponential to deal with the lorentz index?
LatticeColourMatrix Umu(U._grid);
LatticeColourMatrix Pmu(U._grid);
for (int mu = 0; mu < Nd; mu++){
Umu=PeekIndex<LorentzIndex>(U, mu);
Pmu=PeekIndex<LorentzIndex>(*P, mu);
Umu = expMat(Pmu, Complex(ep, 0.0))*Umu;
}
}
void register_observers();
void notify_observers();
friend void IntegratorPolicy::step (LatticeLorentzColourMatrix& U,
int level, std::vector<int>& clock,
Integrator<LeapFrog>* Integ);
public:
Integrator(IntegratorParameters Par,
ActionSet& Aset, std::vector<int> Nrel_):
Params(Par),as(Aset),Nrel(Nrel_){
assert(as.size() == Nrel.size());
};
~Integrator(){}
//Initialization of momenta and actions
void init(LatticeLorentzColourMatrix& U,
GridParallelRNG& pRNG){
std::cout<< "Integrator init\n";
if (!P)
P = new LatticeLorentzColourMatrix(U._grid);
MDutils::generate_momenta(*P,pRNG);
for(int level=0; level< as.size(); ++level){
for(int actionID=0; actionID<as.at(level).size(); ++actionID){
as[level].at(actionID)->init(U, pRNG);
}
}
}
RealD S(LatticeLorentzColourMatrix& U){
// Momenta
LatticeComplex Hloc = - trace((*P)*adj(*P));
Complex Hsum = sum(Hloc);
RealD H = Hsum.real();
// Actions
for(int level=0; level<as.size(); ++level)
for(int actionID=0; actionID<as.at(level).size(); ++actionID)
H += as[level].at(actionID)->S(U);
return H;
}
void integrate(LatticeLorentzColourMatrix& U, int level){
std::vector<int> clock;
clock.resize(as.size(),0);
for(int step=0; step< Params.MDsteps; ++step) // MD step
TheIntegrator.step(U,0,clock, *(this));
}
};
class MinimumNorm2{
const double lambda = 0.1931833275037836;
public:
void step (LatticeLorentzColourMatrix& U, int level, std::vector<int>& clock);
};
class LeapFrog{
public:
void step (LatticeLorentzColourMatrix& U,
int level, std::vector<int>& clock,
Integrator<LeapFrog>* Integ){
// cl : current level
// fl : final level
// eps : current step size
int fl = Integ->as.size() -1;
double eps = Integ->Params.stepsize;
// Get current level step size
for(int l=0; l<=level; ++l) eps/= Integ->Nrel[l];
int fin = 1;
for(int l=0; l<=level; ++l) fin*= Integ->Nrel[l];
fin = 2*Integ->Params.MDsteps*fin - 1;
for(int e=0; e<Integ->Nrel[level]; ++e){
if(clock[level] == 0){ // initial half step
Integ->update_P(U, level,eps/2);
++clock[level];
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"P "<< 0.5*clock[level] <<std::endl;
}
if(level == fl){ // lowest level
Integ->update_U(U, eps);
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"U "<< 0.5*(clock[level]+1) <<std::endl;
}else{ // recursive function call
step(U, level+1,clock, Integ);
}
if(clock[level] == fin){ // final half step
Integ->update_P(U, level,eps/2);
++clock[level];
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"P "<< 0.5*clock[level] <<std::endl;
}else{ // bulk step
Integ->update_P(U, level,eps);
clock[level]+=2;
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"P "<< 0.5*clock[level] <<std::endl;
}
}
}
};
}
}
#endif//INTEGRATOR_INCLUDED

176
tests/Test_dwf_hdcr.cc
View File

@ -6,6 +6,10 @@ using namespace std;
using namespace Grid;
using namespace Grid::QCD;
RealD InverseApproximation(RealD x){
return 1.0/x;
}
template<class Fobj,class CComplex,int nbasis, class Matrix>
class MultiGridPreconditioner : public LinearFunction< Lattice<Fobj> > {
public:
@ -33,6 +37,27 @@ public:
_Matrix(FineMatrix)
{
}
void PowerMethod(const FineField &in) {
FineField p1(in._grid);
FineField p2(in._grid);
MdagMLinearOperator<Matrix,FineField> fMdagMOp(_Matrix);
p1=in;
RealD absp2;
for(int i=0;i<20;i++){
RealD absp1=std::sqrt(norm2(p1));
fMdagMOp.HermOp(p1,p2);// this is the G5 herm bit
// _FineOperator.Op(p1,p2);// this is the G5 herm bit
RealD absp2=std::sqrt(norm2(p2));
if(i%10==9)
std::cout << "Power method on mdagm "<<i<<" " << absp2/absp1<<std::endl;
p1=p2*(1.0/std::sqrt(absp2));
}
}
#if 0
void operator()(const FineField &in, FineField & out) {
@ -49,7 +74,7 @@ public:
std::cout<<"Completeness: "<<std::sqrt(norm2(out)/norm2(in))<<std::endl;
// Build some solvers
ConjugateGradient<FineField> fCG(1.0e-1,1000);
ConjugateGradient<FineField> fCG(1.0e-3,1000);
ConjugateGradient<CoarseVector> CG(1.0e-8,100000);
////////////////////////////////////////////////////////////////////////
@ -164,6 +189,7 @@ public:
}
#endif
// ADEF1: [MP+Q ] in =M [1 - A Q] in + Q in
#if 0
void operator()(const FineField &in, FineField & out) {
CoarseVector Csrc(_CoarseOperator.Grid());
@ -171,11 +197,11 @@ public:
CoarseVector Csol(_CoarseOperator.Grid()); Csol=zero;
ConjugateGradient<CoarseVector> CG(1.0e-10,100000);
ConjugateGradient<FineField> fCG(1.0e-3,1000);
ConjugateGradient<FineField> fCG(3.0e-2,1000);
HermitianLinearOperator<CoarseOperator,CoarseVector> HermOp(_CoarseOperator);
MdagMLinearOperator<CoarseOperator,CoarseVector> MdagMOp(_CoarseOperator);
MdagMLinearOperator<Matrix,FineField> fMdagMOp(_Matrix);
ShiftedMdagMLinearOperator<Matrix,FineField> fMdagMOp(_Matrix,0.1);
FineField tmp(in._grid);
FineField res(in._grid);
@ -191,7 +217,7 @@ public:
CG(MdagMOp,Ctmp,Csol);
_Aggregates.PromoteFromSubspace(Csol,Qin);
// Qin=0;
_FineOperator.Op(Qin,tmp);// A Q in
tmp = in - tmp; // in - A Q in
@ -206,6 +232,116 @@ public:
std::cout<<"preconditioner thinks residual is "<<std::sqrt(norm2(tmp)/norm2(in))<<std::endl;
}
#endif
void SmootherTest (const FineField & in){
FineField vec1(in._grid);
FineField vec2(in._grid);
RealD lo[3] = { 0.5, 1.0, 2.0};
// MdagMLinearOperator<Matrix,FineField> fMdagMOp(_Matrix);
ShiftedMdagMLinearOperator<Matrix,FineField> fMdagMOp(_Matrix,0.5);
RealD Ni,r;
Ni = norm2(in);
for(int ilo=0;ilo<4;ilo++){
for(int ord=10;ord<60;ord+=10){
_FineOperator.AdjOp(in,vec1);
Chebyshev<FineField> Cheby (lo[ilo],70.0,ord,InverseApproximation);
Cheby(fMdagMOp,vec1,vec2); // solves MdagM = g5 M g5M
_FineOperator.Op(vec2,vec1);// this is the G5 herm bit
vec1 = in - vec1; // tmp = in - A Min
r=norm2(vec1);
std::cout << "Smoother resid "<<std::sqrt(r/Ni)<<std::endl;
}
}
}
void operator()(const FineField &in, FineField & out) {
CoarseVector Csrc(_CoarseOperator.Grid());
CoarseVector Ctmp(_CoarseOperator.Grid());
CoarseVector Csol(_CoarseOperator.Grid()); Csol=zero;
ConjugateGradient<CoarseVector> CG(1.0e-5,100000);
ConjugateGradient<FineField> fCG(3.0e-2,1000);
HermitianLinearOperator<CoarseOperator,CoarseVector> HermOp(_CoarseOperator);
MdagMLinearOperator<CoarseOperator,CoarseVector> MdagMOp(_CoarseOperator);
// MdagMLinearOperator<Matrix,FineField> fMdagMOp(_Matrix);
ShiftedMdagMLinearOperator<Matrix,FineField> fMdagMOp(_Matrix,0.5);
FineField vec1(in._grid);
FineField vec2(in._grid);
// Chebyshev<FineField> Cheby (0.5,70.0,30,InverseApproximation);
// Chebyshev<FineField> ChebyAccu(0.5,70.0,30,InverseApproximation);
Chebyshev<FineField> Cheby (1.0,70.0,20,InverseApproximation);
Chebyshev<FineField> ChebyAccu(1.0,70.0,20,InverseApproximation);
_Aggregates.ProjectToSubspace (Csrc,in);
_Aggregates.PromoteFromSubspace(Csrc,out);
std::cout<<"Completeness: "<<std::sqrt(norm2(out)/norm2(in))<<std::endl;
// ofstream fout("smoother");
// Cheby.csv(fout);
// V11 multigrid.
// Use a fixed chebyshev and hope hermiticity helps.
// To make a working smoother for indefinite operator
// must multiply by "Mdag" (ouch loses all low mode content)
// and apply to poly approx of (mdagm)^-1.
// so that we end up with an odd polynomial.
RealD Ni = norm2(in);
_FineOperator.AdjOp(in,vec1);// this is the G5 herm bit
ChebyAccu(fMdagMOp,vec1,out); // solves MdagM = g5 M g5M
std::cout << "Smoother norm "<<norm2(out)<<std::endl;
// Update with residual for out
_FineOperator.Op(out,vec1);// this is the G5 herm bit
vec1 = in - vec1; // tmp = in - A Min
RealD r = norm2(vec1);
std::cout << "Smoother resid "<<std::sqrt(r/Ni)<< " " << r << " " << Ni <<std::endl;
_Aggregates.ProjectToSubspace (Csrc,vec1);
HermOp.AdjOp(Csrc,Ctmp);// Normal equations
CG(MdagMOp,Ctmp,Csol);
_Aggregates.PromoteFromSubspace(Csol,vec1); // Ass^{-1} [in - A Min]_s
// Q = Q[in - A Min]
out = out+vec1;
// Three preconditioner smoothing -- hermitian if C3 = C1
// Recompute error
_FineOperator.Op(out,vec1);// this is the G5 herm bit
vec1 = in - vec1; // tmp = in - A Min
r=norm2(vec1);
std::cout << "Coarse resid "<<std::sqrt(r/Ni)<<std::endl;
// Reapply smoother
_FineOperator.Op(vec1,vec2); // this is the G5 herm bit
ChebyAccu(fMdagMOp,vec2,vec1); // solves MdagM = g5 M g5M
out =out+vec1;
_FineOperator.Op(out,vec1);// this is the G5 herm bit
vec1 = in - vec1; // tmp = in - A Min
r=norm2(vec1);
std::cout << "Smoother resid "<<std::sqrt(r/Ni)<<std::endl;
}
};
@ -224,7 +360,7 @@ int main (int argc, char ** argv)
///////////////////////////////////////////////////
// Construct a coarsened grid; utility for this?
///////////////////////////////////////////////////
const int block=4;
const int block=2;
std::vector<int> clatt = GridDefaultLatt();
for(int d=0;d<clatt.size();d++){
clatt[d] = clatt[d]/block;
@ -265,7 +401,8 @@ int main (int argc, char ** argv)
std::cout << "**************************************************"<< std::endl;
DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
const int nbasis = 6;
const int nbasis = 32;
// const int nbasis = 4;
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> CoarseOperator;
@ -276,7 +413,19 @@ int main (int argc, char ** argv)
std::cout << "**************************************************"<< std::endl;
MdagMLinearOperator<DomainWallFermion,LatticeFermion> HermDefOp(Ddwf);
Subspace Aggregates(Coarse5d,FGrid);
Aggregates.CreateSubspace(RNG5,HermDefOp);
// Aggregates.CreateSubspace(RNG5,HermDefOp,nbasis);
assert ( (nbasis & 0x1)==0);
int nb=nbasis/2;
std::cout << " nbasis/2 = "<<nb<<std::endl;
Aggregates.CreateSubspace(RNG5,HermDefOp,nb);
for(int n=0;n<nb;n++){
G5R5(Aggregates.subspace[n+nb],Aggregates.subspace[n]);
std::cout<<n<<" subspace "<<norm2(Aggregates.subspace[n+nb])<<" "<<norm2(Aggregates.subspace[n]) <<std::endl;
}
for(int n=0;n<nbasis;n++){
std::cout << "vec["<<n<<"] = "<<norm2(Aggregates.subspace[n]) <<std::endl;
}
// for(int i=0;i<nbasis;i++){
// result = Aggregates.subspace[i];
// Aggregates.subspace[i]=result+g5*result;
@ -319,12 +468,18 @@ int main (int argc, char ** argv)
MultiGridPreconditioner <vSpinColourVector,vTComplex,nbasis,DomainWallFermion> Precon(Aggregates, LDOp,HermIndefOp,Ddwf);
TrivialPrecon<LatticeFermion> simple;
std::cout << "**************************************************"<< std::endl;
std::cout << "Testing smoother efficacy"<< std::endl;
std::cout << "**************************************************"<< std::endl;
Precon.SmootherTest(src);
std::cout << "**************************************************"<< std::endl;
std::cout << "Unprec CG "<< std::endl;
std::cout << "**************************************************"<< std::endl;
// TrivialPrecon<LatticeFermion> simple;
ConjugateGradient<LatticeFermion> fCG(1.0e-8,100000);
fCG(HermDefOp,src,result);
// exit(0);
std::cout << "**************************************************"<< std::endl;
std::cout << "Testing GCR on indef matrix "<< std::endl;
@ -332,6 +487,13 @@ int main (int argc, char ** argv)
// PrecGeneralisedConjugateResidual<LatticeFermion> UPGCR(1.0e-8,100000,simple,8,128);
// UPGCR(HermIndefOp,src,result);
/// Get themax eval
std::cout << "**************************************************"<< std::endl;
std::cout <<" Applying power method to find spectral range "<<std::endl;
std::cout << "**************************************************"<< std::endl;
Precon.PowerMethod(src);
std::cout << "**************************************************"<< std::endl;
std::cout << "Building a two level PGCR "<< std::endl;
std::cout << "**************************************************"<< std::endl;

35
tests/Test_multishift_sqrt.cc
View File

@ -17,8 +17,12 @@ public:
pRNG.SeedFixedIntegers(seeds);
random(pRNG,sqrtscale);
sqrtscale = sqrtscale * adj(sqrtscale);// force real pos def
scale = sqrtscale * sqrtscale;
sqrtscale = real(sqrtscale)*3.0+0.5;// force real pos def
scale = sqrtscale *sqrtscale; //scale should be bounded by 12.25
//
// sqrtscale = 2.0;
// scale = 4.0;
}
// Support for coarsening to a multigrid
void OpDiag (const Field &in, Field &out) {};
@ -48,8 +52,18 @@ public:
void ApplySqrt(const Field &in, Field &out){
out = sqrtscale * in;
}
void ApplyInverse(const Field &in, Field &out){
out = pow(scale,-1.0) * in;
}
};
RealD InverseApproximation(RealD x){
return 1.0/x;
}
RealD SqrtApproximation(RealD x){
return std::sqrt(x);
}
int main (int argc, char ** argv)
{
@ -114,5 +128,22 @@ int main (int argc, char ** argv)
error = summed - combined;
std::cout << " summed-combined "<<norm2(error) <<std::endl;
src=1.0;
Chebyshev<LatticeFermion> Cheby(0.1,40.0,200,InverseApproximation);
std::cout<<"Chebuy approx vector "<<std::endl;
Cheby(Diagonal,src,combined);
std::ofstream of("cheby");
Cheby.csv(of);
Diagonal.ApplyInverse(src,reference);
error = reference - combined;
std::cout << "Chebyshev inverse test "<<std::endl;
std::cout << " Reference "<<norm2(reference)<<std::endl;
std::cout << " combined "<<norm2(combined) <<std::endl;
std::cout << " error "<<norm2(error) <<std::endl;
Grid_finalize();
}