From 09651c33262e0ddf9c9a44e6d9e8fa260124d753 Mon Sep 17 00:00:00 2001 From: Chulwoo Jung Date: Tue, 2 May 2017 00:47:18 -0400 Subject: [PATCH] Checking in before rearranging Lanczos --- lib/algorithms/LinearOperator.h | 69 +++++++++++ .../iterative/ImplicitlyRestartedLanczos.h | 4 +- lib/algorithms/iterative/SchurRedBlack.h | 1 + lib/qcd/action/fermion/Fermion.h | 2 +- lib/qcd/action/fermion/SchurDiagTwoKappa.h | 112 ++++++++++++++++++ .../EvenOddSchurDifferentiable.h | 1 + 6 files changed, 186 insertions(+), 3 deletions(-) diff --git a/lib/algorithms/LinearOperator.h b/lib/algorithms/LinearOperator.h index ea47d43b..9859d032 100644 --- a/lib/algorithms/LinearOperator.h +++ b/lib/algorithms/LinearOperator.h @@ -300,6 +300,75 @@ namespace Grid { } }; +#if 0 + // This is specific to (Z)mobius fermions + template + class KappaSimilarityTransform { + public: +// INHERIT_IMPL_TYPES(Matrix); + typedef typename Matrix::Coeff_t Coeff_t; + std::vector kappa, kappaDag, kappaInv, kappaInvDag; + + KappaSimilarityTransform (Matrix &zmob) { + for (int i=0;i<(int)zmob.bs.size();i++) { + Coeff_t k = 1.0 / ( 2.0 * (zmob.bs[i] *(4 - zmob.M5) + 1.0) ); + kappa.push_back( k ); + kappaDag.push_back( conj(k) ); + kappaInv.push_back( 1.0 / k ); + kappaInvDag.push_back( 1.0 / conj(k) ); + } + } + + template + void sscale(const Lattice& in, Lattice& out, Coeff_t* s) { + GridBase *grid=out._grid; + out.checkerboard = in.checkerboard; + assert(grid->_simd_layout[0] == 1); // should be fine for ZMobius for now + int Ls = grid->_rdimensions[0]; + parallel_for(int ss=0;ssoSites();ss++){ + vobj tmp = s[ss % Ls]*in._odata[ss]; + vstream(out._odata[ss],tmp); + } + } + + RealD sscale_norm(const Field& in, Field& out, Coeff_t* s) { + sscale(in,out,s); + return norm2(out); + } + + virtual RealD M (const Field& in, Field& out) { return sscale_norm(in,out,&kappa[0]); } + virtual RealD MDag (const Field& in, Field& out) { return sscale_norm(in,out,&kappaDag[0]);} + virtual RealD MInv (const Field& in, Field& out) { return sscale_norm(in,out,&kappaInv[0]);} + virtual RealD MInvDag (const Field& in, Field& out) { return sscale_norm(in,out,&kappaInvDag[0]);} + + }; + + template + class SchurDiagTwoKappaOperator : public SchurOperatorBase { + public: + KappaSimilarityTransform _S; + SchurDiagTwoOperator _Mat; + + SchurDiagTwoKappaOperator (Matrix &Mat): _S(Mat), _Mat(Mat) {}; + + virtual RealD Mpc (const Field &in, Field &out) { + Field tmp(in._grid); + + _S.MInv(in,out); + _Mat.Mpc(out,tmp); + return _S.M(tmp,out); + + } + virtual RealD MpcDag (const Field &in, Field &out){ + Field tmp(in._grid); + + _S.MDag(in,out); + _Mat.MpcDag(out,tmp); + return _S.MInvDag(tmp,out); + } + }; +#endif + ///////////////////////////////////////////////////////////// // Base classes for functions of operators diff --git a/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h b/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h index a16a7ebf..33248b4f 100644 --- a/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h +++ b/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h @@ -654,8 +654,8 @@ PARALLEL_FOR_LOOP } for(int j=k1-1; j #include #include #include -#include #include +//#include #include #include #include diff --git a/lib/qcd/action/fermion/SchurDiagTwoKappa.h b/lib/qcd/action/fermion/SchurDiagTwoKappa.h index 8305f98a..e5406778 100644 --- a/lib/qcd/action/fermion/SchurDiagTwoKappa.h +++ b/lib/qcd/action/fermion/SchurDiagTwoKappa.h @@ -1,3 +1,4 @@ +#if 1 /************************************************************************************* Grid physics library, www.github.com/paboyle/Grid @@ -97,6 +98,117 @@ namespace Grid { } }; +#if 0 + /////////////////////////////////////////////////////////////////////////////////////////////////////// + // Copied from DiagTwoSolve + /////////////////////////////////////////////////////////////////////////////////////////////////////// + template class SchurRedBlackDiagTwoSolve { + private: + OperatorFunction & _HermitianRBSolver; + int CBfactorise; + public: + + ///////////////////////////////////////////////////// + // Wrap the usual normal equations Schur trick + ///////////////////////////////////////////////////// + SchurRedBlackDiagTwoSolve(OperatorFunction &HermitianRBSolver) : + _HermitianRBSolver(HermitianRBSolver) + { + CBfactorise=0; + }; + + template + void operator() (Matrix & _Matrix,const Field &in, Field &out){ + + // FIXME CGdiagonalMee not implemented virtual function + // FIXME use CBfactorise to control schur decomp + GridBase *grid = _Matrix.RedBlackGrid(); + GridBase *fgrid= _Matrix.Grid(); + + SchurDiagTwoOperator _HermOpEO(_Matrix); + + Field src_e(grid); + Field src_o(grid); + Field sol_e(grid); + Field sol_o(grid); + Field tmp(grid); + Field Mtmp(grid); + Field resid(fgrid); + + pickCheckerboard(Even,src_e,in); + pickCheckerboard(Odd ,src_o,in); + pickCheckerboard(Even,sol_e,out); + pickCheckerboard(Odd ,sol_o,out); + + ///////////////////////////////////////////////////// + // src_o = Mdag * (source_o - Moe MeeInv source_e) + ///////////////////////////////////////////////////// + _Matrix.MooeeInv(src_e,tmp); assert( tmp.checkerboard ==Even); + _Matrix.Meooe (tmp,Mtmp); assert( Mtmp.checkerboard ==Odd); + tmp=src_o-Mtmp; assert( tmp.checkerboard ==Odd); + + // get the right MpcDag + _HermOpEO.MpcDag(tmp,src_o); assert(src_o.checkerboard ==Odd); + + ////////////////////////////////////////////////////////////// + // Call the red-black solver + ////////////////////////////////////////////////////////////// + std::cout< + class SchurDifferentiableDiagTwo: public SchurDiagTwoOperator,typename Impl::FermionField> + { + public: + INHERIT_IMPL_TYPES(Impl); + + typedef FermionOperator Matrix; + + SchurDifferentiableDiagTwo (Matrix &Mat) : SchurDiagTwoOperator(Mat) {}; + }; +#if 0 + template + class SchurDifferentiableDiagTwoKappa : public SchurDiagTwoKappaOperator,typename Impl::FermionField> + { + public: + INHERIT_IMPL_TYPES(Impl); + + typedef FermionOperator Matrix; + + SchurDifferentiableDiagTwoKappa (Matrix &Mat) : SchurDiagTwoKappaOperator(Mat) {}; + }; +#endif +} + } #endif +#endif diff --git a/lib/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h b/lib/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h index 6837bb19..b47540c7 100644 --- a/lib/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h +++ b/lib/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h @@ -134,6 +134,7 @@ namespace Grid{ }; + } } #endif