diff --git a/tests/solver/Test_hw_multigrid.cc b/tests/solver/Test_hw_multigrid.cc index 6fd7811a..b531decc 100644 --- a/tests/solver/Test_hw_multigrid.cc +++ b/tests/solver/Test_hw_multigrid.cc @@ -39,12 +39,13 @@ using namespace Grid; // Coarse Grid axpby_ssp_pplus template -class CayleyBase +class CayleyBase : public SparseMatrixBase { public: int Ls; // protected: RealD mass; + RealD M5; // Save arguments to SetCoefficientsInternal Vector _gamma; RealD _zolo_hi; @@ -70,6 +71,18 @@ public: Vector ueem; Vector dee; public: + CayleyBase(RealD _M5, RealD _mass, int _Ls, RealD b_, RealD c_) : + M5(_M5), + mass(_mass), + Ls(_Ls), + _b(b_), + _c(c_) + { + RealD eps = 1.0; + Approx::zolotarev_data *zdata = Approx::higham(eps,this->Ls);// eps is ignored for higham + this->SetCoefficientsTanh(zdata,1.0,0.0); + Approx::zolotarev_free(zdata); + } ///////////////////////////////////////////////////////// // Replicates functionality // Use a common base class approach @@ -192,6 +205,12 @@ public: ////////////////////////////// // M and Mdag ////////////////////////////// + virtual void Mdiag (const Field &in, Field &out) {} + virtual void Mdir (const Field &in, Field &out,int dir, int disp){}; + virtual void MdirAll (const Field &in, std::vector &out){}; + virtual void DW (const Field &psi, Field &chi)=0; + virtual void DWDag (const Field &psi, Field &chi)=0; + void M (const Field &psi, Field &chi) { Field Din(psi.Grid()); @@ -359,9 +378,7 @@ public: }; template -class CoarseCayleyFermion : public - CayleyBase< Lattice > , ComplexD >, - SparseMatrixBase > > +class CoarseCayleyFermion : public CayleyBase< Lattice > , ComplexD > { public: typedef iVector siteVector; @@ -383,18 +400,15 @@ public: CoarseCayleyFermion(GridCartesian &CoarseGrid4, GridCartesian &CoarseGrid5, - CoarsenedMatrix &_Dw ): + CoarsenedMatrix &_Dw, + RealD M5, RealD mass, int Ls, RealD b, RealD c) : + CayleyBase(M5,mass,Ls,b,c), Coarse4D(&CoarseGrid4), Coarse5D(&CoarseGrid5), Dw(_Dw), geom(CoarseGrid4._ndimension), Stencil( &CoarseGrid4,geom.npoint,Even,geom.directions,geom.displacements,0) { - this->Ls=Coarse5D->_fdimensions[0]; - RealD eps = 1.0; - Approx::zolotarev_data *zdata = Approx::higham(eps,this->Ls);// eps is ignored for higham - this->SetCoefficientsTanh(zdata,1.0,0.0); - Approx::zolotarev_free(zdata); }; public: @@ -462,7 +476,7 @@ public: // Inefficient G5 hermitian use CoarseVector tmp(Grid()); G5C(tmp, in); //There has to be a better way - M(tmp, out); + DW(tmp, out); G5C(out, out); }; @@ -741,13 +755,11 @@ int main (int argc, char ** argv) std::cout< Level1Op; + typedef CoarsenedMatrix Level1Op4; + typedef CoarseCayleyFermion Level1Op5; - // NonHermitianLinearOperator LinOpDwf(Ddwf); - - Level1Op c_Dw (*Coarse4d,0); - - + Level1Op4 c_Dw (*Coarse4d,0); + Level1Op5 c_Dwf (*Coarse4d,*Coarse5d,c_Dw,M5, mass, Ls, 1.0,0.0); std::cout< LinOpDw(Dw); c_Dw.CoarsenOperator(UGrid,LinOpDw,Aggregates4D);