diff --git a/Grid/algorithms/LinearOperator.h b/Grid/algorithms/LinearOperator.h index 86d3f030..b86863b8 100644 --- a/Grid/algorithms/LinearOperator.h +++ b/Grid/algorithms/LinearOperator.h @@ -380,6 +380,12 @@ namespace Grid { template class OperatorFunction { public: virtual void operator() (LinearOperatorBase &Linop, const Field &in, Field &out) = 0; + virtual void operator() (LinearOperatorBase &Linop, const std::vector &in,std::vector &out) { + assert(in.size()==out.size()); + for(int k=0;k class LinearFunction { @@ -421,7 +427,7 @@ namespace Grid { // Hermitian operator Linear function and operator function //////////////////////////////////////////////////////////////////////////////////////////// template - class HermOpOperatorFunction : public OperatorFunction { + class HermOpOperatorFunction : public OperatorFunction { void operator() (LinearOperatorBase &Linop, const Field &in, Field &out) { Linop.HermOp(in,out); }; diff --git a/Grid/algorithms/SparseMatrix.h b/Grid/algorithms/SparseMatrix.h index 1611a6f4..9b0f7659 100644 --- a/Grid/algorithms/SparseMatrix.h +++ b/Grid/algorithms/SparseMatrix.h @@ -55,6 +55,14 @@ namespace Grid { template class CheckerBoardedSparseMatrixBase : public SparseMatrixBase { public: virtual GridBase *RedBlackGrid(void)=0; + + ////////////////////////////////////////////////////////////////////// + // Query the even even properties to make algorithmic decisions + ////////////////////////////////////////////////////////////////////// + virtual RealD Mass(void) { return 0.0; }; + virtual int ConstEE(void) { return 0; }; // Disable assumptions unless overridden + virtual int isTrivialEE(void) { return 0; }; // by a derived class that knows better + // half checkerboard operaions virtual void Meooe (const Field &in, Field &out)=0; virtual void Mooee (const Field &in, Field &out)=0; diff --git a/Grid/algorithms/iterative/BlockConjugateGradient.h b/Grid/algorithms/iterative/BlockConjugateGradient.h index e0eeddcb..972aa608 100644 --- a/Grid/algorithms/iterative/BlockConjugateGradient.h +++ b/Grid/algorithms/iterative/BlockConjugateGradient.h @@ -33,7 +33,7 @@ directory namespace Grid { -enum BlockCGtype { BlockCG, BlockCGrQ, CGmultiRHS }; +enum BlockCGtype { BlockCG, BlockCGrQ, CGmultiRHS, BlockCGVec, BlockCGrQVec }; ////////////////////////////////////////////////////////////////////////// // Block conjugate gradient. Dimension zero should be the block direction @@ -42,7 +42,6 @@ template class BlockConjugateGradient : public OperatorFunction { public: - typedef typename Field::scalar_type scomplex; int blockDim ; @@ -54,21 +53,15 @@ class BlockConjugateGradient : public OperatorFunction { RealD Tolerance; Integer MaxIterations; Integer IterationsToComplete; //Number of iterations the CG took to finish. Filled in upon completion + Integer PrintInterval; //GridLogMessages or Iterative BlockConjugateGradient(BlockCGtype cgtype,int _Orthog,RealD tol, Integer maxit, bool err_on_no_conv = true) - : Tolerance(tol), CGtype(cgtype), blockDim(_Orthog), MaxIterations(maxit), ErrorOnNoConverge(err_on_no_conv) + : Tolerance(tol), CGtype(cgtype), blockDim(_Orthog), MaxIterations(maxit), ErrorOnNoConverge(err_on_no_conv),PrintInterval(100) {}; //////////////////////////////////////////////////////////////////////////////////////////////////// // Thin QR factorisation (google it) //////////////////////////////////////////////////////////////////////////////////////////////////// -void ThinQRfact (Eigen::MatrixXcd &m_rr, - Eigen::MatrixXcd &C, - Eigen::MatrixXcd &Cinv, - Field & Q, - const Field & R) -{ - int Orthog = blockDim; // First dimension is block dim; this is an assumption //////////////////////////////////////////////////////////////////////////////////////////////////// //Dimensions // R_{ferm x Nblock} = Q_{ferm x Nblock} x C_{Nblock x Nblock} -> ferm x Nblock @@ -85,22 +78,20 @@ void ThinQRfact (Eigen::MatrixXcd &m_rr, // Cdag C = Rdag R ; passes. // QdagQ = 1 ; passes //////////////////////////////////////////////////////////////////////////////////////////////////// +void ThinQRfact (Eigen::MatrixXcd &m_rr, + Eigen::MatrixXcd &C, + Eigen::MatrixXcd &Cinv, + Field & Q, + const Field & R) +{ + int Orthog = blockDim; // First dimension is block dim; this is an assumption sliceInnerProductMatrix(m_rr,R,R,Orthog); // Force manifest hermitian to avoid rounding related m_rr = 0.5*(m_rr+m_rr.adjoint()); -#if 0 - std::cout << " Calling Cholesky ldlt on m_rr " << m_rr < & Q, + const std::vector & R) +{ + InnerProductMatrix(m_rr,R,R); + + m_rr = 0.5*(m_rr+m_rr.adjoint()); + + Eigen::MatrixXcd L = m_rr.llt().matrixL(); + + C = L.adjoint(); + Cinv = C.inverse(); + + MulMatrix(Q,Cinv,R); +} + //////////////////////////////////////////////////////////////////////////////////////////////////// // Call one of several implementations //////////////////////////////////////////////////////////////////////////////////////////////////// @@ -119,14 +129,20 @@ void operator()(LinearOperatorBase &Linop, const Field &Src, Field &Psi) { if ( CGtype == BlockCGrQ ) { BlockCGrQsolve(Linop,Src,Psi); - } else if (CGtype == BlockCG ) { - BlockCGsolve(Linop,Src,Psi); } else if (CGtype == CGmultiRHS ) { CGmultiRHSsolve(Linop,Src,Psi); } else { assert(0); } } +virtual void operator()(LinearOperatorBase &Linop, const std::vector &Src, std::vector &Psi) +{ + if ( CGtype == BlockCGrQVec ) { + BlockCGrQsolveVec(Linop,Src,Psi); + } else { + assert(0); + } +} //////////////////////////////////////////////////////////////////////////// // BlockCGrQ implementation: @@ -139,7 +155,8 @@ void BlockCGrQsolve(LinearOperatorBase &Linop, const Field &B, Field &X) { int Orthog = blockDim; // First dimension is block dim; this is an assumption Nblock = B._grid->_fdimensions[Orthog]; - +/* FAKE */ + Nblock=8; std::cout< &Linop, const Field &B, Field &X) std::cout << GridLogMessage<<"BlockCGrQ algorithm initialisation " < Thin QR factorisation (google it) - Linop.HermOp(X, AD); tmp = B - AD; - //std::cout << GridLogMessage << " initial tmp " << norm2(tmp)<< std::endl; + ThinQRfact (m_rr, m_C, m_Cinv, Q, tmp); - //std::cout << GridLogMessage << " initial Q " << norm2(Q)<< std::endl; - //std::cout << GridLogMessage << " m_rr " << m_rr< &Linop, const Field &B, Field &X) MatrixTimer.Start(); Linop.HermOp(D, Z); MatrixTimer.Stop(); - //std::cout << GridLogMessage << " norm2 Z " < &Linop, const Field &B, Field &X) //7. D = Q + D S^dag m_tmp = m_S.adjoint(); + sliceMaddTimer.Start(); sliceMaddMatrix(D,m_tmp,D,Q,Orthog); sliceMaddTimer.Stop(); @@ -317,152 +328,6 @@ void BlockCGrQsolve(LinearOperatorBase &Linop, const Field &B, Field &X) IterationsToComplete = k; } ////////////////////////////////////////////////////////////////////////// -// Block conjugate gradient; Original O'Leary Dimension zero should be the block direction -////////////////////////////////////////////////////////////////////////// -void BlockCGsolve(LinearOperatorBase &Linop, const Field &Src, Field &Psi) -{ - int Orthog = blockDim; // First dimension is block dim; this is an assumption - Nblock = Src._grid->_fdimensions[Orthog]; - - std::cout< residuals(Nblock); - std::vector ssq(Nblock); - - sliceNorm(ssq,Src,Orthog); - RealD sssum=0; - for(int b=0;b max_resid ) max_resid = rr; - } - - if ( max_resid < Tolerance*Tolerance ) { - - SolverTimer.Stop(); - - std::cout << GridLogMessage<<"BlockCG converged in "< &Linop, const Field &Src, Field & IterationsToComplete = k; } +void InnerProductMatrix(Eigen::MatrixXcd &m , const std::vector &X, const std::vector &Y){ + for(int b=0;b &AP, Eigen::MatrixXcd &m , const std::vector &X,const std::vector &Y,RealD scale=1.0){ + // Should make this cache friendly with site outermost, parallel_for + // Deal with case AP aliases with either Y or X + std::vector tmp(Nblock,X[0]); + for(int b=0;b &AP, Eigen::MatrixXcd &m , const std::vector &X){ + // Should make this cache friendly with site outermost, parallel_for + for(int b=0;b &P){ + double nn = 0.0; + for(int b=0;b &Linop, const std::vector &B, std::vector &X) +{ + Nblock = B.size(); + assert(Nblock == X.size()); + + std::cout< tmp(Nblock,Fake); + std::vector Q(Nblock,Fake); + std::vector D(Nblock,Fake); + std::vector Z(Nblock,Fake); + std::vector AD(Nblock,Fake); + + Eigen::MatrixXcd m_DZ = Eigen::MatrixXcd::Identity(Nblock,Nblock); + Eigen::MatrixXcd m_M = Eigen::MatrixXcd::Identity(Nblock,Nblock); + Eigen::MatrixXcd m_rr = Eigen::MatrixXcd::Zero(Nblock,Nblock); + + Eigen::MatrixXcd m_C = Eigen::MatrixXcd::Zero(Nblock,Nblock); + Eigen::MatrixXcd m_Cinv = Eigen::MatrixXcd::Zero(Nblock,Nblock); + Eigen::MatrixXcd m_S = Eigen::MatrixXcd::Zero(Nblock,Nblock); + Eigen::MatrixXcd m_Sinv = Eigen::MatrixXcd::Zero(Nblock,Nblock); + + Eigen::MatrixXcd m_tmp = Eigen::MatrixXcd::Identity(Nblock,Nblock); + Eigen::MatrixXcd m_tmp1 = Eigen::MatrixXcd::Identity(Nblock,Nblock); + + // Initial residual computation & set up + std::vector residuals(Nblock); + std::vector ssq(Nblock); + + RealD sssum=0; + for(int b=0;b Thin QR factorisation (google it) + * for k: + * Z = AD + * M = [D^dag Z]^{-1} + * X = X + D MC + * QS = Q - ZM + * D = Q + D S^dag + * C = S C + */ + /////////////////////////////////////// + // Initial block: initial search dir is guess + /////////////////////////////////////// + std::cout << GridLogMessage<<"BlockCGrQvec algorithm initialisation " < Thin QR factorisation (google it) + for(int b=0;b max_resid ) max_resid = rr; + } + + std::cout << GridLogIterative << "\t Block Iteration "< { LinalgTimer.Stop(); std::cout << GridLogIterative << "ConjugateGradient: Iteration " << k - << " residual " << cp << " target " << rsq << std::endl; + << " residual^2 " << sqrt(cp/ssq) << " target " << Tolerance << std::endl; // Stopping condition if (cp <= rsq) { diff --git a/Grid/algorithms/iterative/SchurRedBlack.h b/Grid/algorithms/iterative/SchurRedBlack.h index 5abc4d9f..fdb17a98 100644 --- a/Grid/algorithms/iterative/SchurRedBlack.h +++ b/Grid/algorithms/iterative/SchurRedBlack.h @@ -86,229 +86,23 @@ Author: Peter Boyle */ namespace Grid { + /////////////////////////////////////////////////////////////////////////////////////////////////////// + // Use base class to share code + /////////////////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////////////////// // Take a matrix and form a Red Black solver calling a Herm solver // Use of RB info prevents making SchurRedBlackSolve conform to standard interface /////////////////////////////////////////////////////////////////////////////////////////////////////// - // Now make the norm reflect extra factor of Mee - template class SchurRedBlackStaggeredSolve { - private: + template class SchurRedBlackBase { + protected: + typedef CheckerBoardedSparseMatrixBase Matrix; OperatorFunction & _HermitianRBSolver; int CBfactorise; bool subGuess; public: - ///////////////////////////////////////////////////// - // Wrap the usual normal equations Schur trick - ///////////////////////////////////////////////////// - SchurRedBlackStaggeredSolve(OperatorFunction &HermitianRBSolver, const bool initSubGuess = false) : - _HermitianRBSolver(HermitianRBSolver) - { - CBfactorise=0; - subtractGuess(initSubGuess); - }; - void subtractGuess(const bool initSubGuess) - { - subGuess = initSubGuess; - } - bool isSubtractGuess(void) - { - return subGuess; - } - - template - void operator() (Matrix & _Matrix,const Field &in, Field &out){ - ZeroGuesser guess; - (*this)(_Matrix,in,out,guess); - } - template - void operator() (Matrix & _Matrix,const Field &in, Field &out, Guesser &guess){ - - // FIXME CGdiagonalMee not implemented virtual function - // FIXME use CBfactorise to control schur decomp - GridBase *grid = _Matrix.RedBlackGrid(); - GridBase *fgrid= _Matrix.Grid(); - - SchurStaggeredOperator _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); - - std::cout << GridLogMessage << " SchurRedBlackStaggeredSolve " < using SchurRedBlackStagSolve = SchurRedBlackStaggeredSolve; - - /////////////////////////////////////////////////////////////////////////////////////////////////////// - // Take a matrix and form a Red Black solver calling a Herm solver - // Use of RB info prevents making SchurRedBlackSolve conform to standard interface - /////////////////////////////////////////////////////////////////////////////////////////////////////// - template class SchurRedBlackDiagMooeeSolve { - private: - OperatorFunction & _HermitianRBSolver; - int CBfactorise; - bool subGuess; - public: - - ///////////////////////////////////////////////////// - // Wrap the usual normal equations Schur trick - ///////////////////////////////////////////////////// - SchurRedBlackDiagMooeeSolve(OperatorFunction &HermitianRBSolver,int cb=0, const bool initSubGuess = false) : _HermitianRBSolver(HermitianRBSolver) - { - CBfactorise=cb; - subtractGuess(initSubGuess); - }; - void subtractGuess(const bool initSubGuess) - { - subGuess = initSubGuess; - } - bool isSubtractGuess(void) - { - return subGuess; - } - template - void operator() (Matrix & _Matrix,const Field &in, Field &out){ - ZeroGuesser guess; - (*this)(_Matrix,in,out,guess); - } - template - void operator() (Matrix & _Matrix,const Field &in, Field &out,Guesser &guess){ - - // FIXME CGdiagonalMee not implemented virtual function - // FIXME use CBfactorise to control schur decomp - GridBase *grid = _Matrix.RedBlackGrid(); - GridBase *fgrid= _Matrix.Grid(); - - SchurDiagMooeeOperator _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 SchurRedBlackDiagTwoSolve { - private: - OperatorFunction & _HermitianRBSolver; - int CBfactorise; - bool subGuess; - public: - - ///////////////////////////////////////////////////// - // Wrap the usual normal equations Schur trick - ///////////////////////////////////////////////////// - SchurRedBlackDiagTwoSolve(OperatorFunction &HermitianRBSolver, const bool initSubGuess = false) : - _HermitianRBSolver(HermitianRBSolver) + SchurRedBlackBase(OperatorFunction &HermitianRBSolver, const bool initSubGuess = false) : + _HermitianRBSolver(HermitianRBSolver) { CBfactorise = 0; subtractGuess(initSubGuess); @@ -322,12 +116,86 @@ namespace Grid { return subGuess; } - template + ///////////////////////////////////////////////////////////// + // Shared code + ///////////////////////////////////////////////////////////// void operator() (Matrix & _Matrix,const Field &in, Field &out){ ZeroGuesser guess; (*this)(_Matrix,in,out,guess); } - template + void operator()(Matrix &_Matrix, const std::vector &in, std::vector &out) + { + ZeroGuesser guess; + (*this)(_Matrix,in,out,guess); + } + + template + void operator()(Matrix &_Matrix, const std::vector &in, std::vector &out,Guesser &guess) + { + GridBase *grid = _Matrix.RedBlackGrid(); + GridBase *fgrid= _Matrix.Grid(); + int nblock = in.size(); + + std::vector src_o(nblock,grid); + std::vector sol_o(nblock,grid); + + std::vector guess_save; + + Field resid(fgrid); + Field tmp(grid); + + //////////////////////////////////////////////// + // Prepare RedBlack source + //////////////////////////////////////////////// + for(int b=0;b void operator() (Matrix & _Matrix,const Field &in, Field &out,Guesser &guess){ // FIXME CGdiagonalMee not implemented virtual function @@ -335,52 +203,39 @@ namespace Grid { 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); + Field src_o(grid); + Field src_e(grid); + Field sol_o(grid); - 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); + //////////////////////////////////////////////// + // RedBlack source + //////////////////////////////////////////////// + RedBlackSource(_Matrix,in,src_e,src_o); - // get the right MpcDag - _HermOpEO.MpcDag(tmp,src_o); assert(src_o.checkerboard ==Odd); + //////////////////////////////// + // Construct the guess + //////////////////////////////// + Field tmp(grid); + guess(src_o,sol_o); + + Field guess_save(grid); + guess_save = sol_o; ////////////////////////////////////////////////////////////// // Call the red-black solver ////////////////////////////////////////////////////////////// - std::cout< &src_o, std::vector &sol_o)=0; + }; - /////////////////////////////////////////////////////////////////////////////////////////////////////// - // Take a matrix and form a Red Black solver calling a Herm solver - // Use of RB info prevents making SchurRedBlackSolve conform to standard interface - /////////////////////////////////////////////////////////////////////////////////////////////////////// - template class SchurRedBlackDiagTwoMixed { - private: - LinearFunction & _HermitianRBSolver; - int CBfactorise; - bool subGuess; + + template class SchurRedBlackStaggeredSolve : public SchurRedBlackBase { public: + typedef CheckerBoardedSparseMatrixBase Matrix; + + SchurRedBlackStaggeredSolve(OperatorFunction &HermitianRBSolver, const bool initSubGuess = false) + : SchurRedBlackBase (HermitianRBSolver,initSubGuess) + { + } + + ////////////////////////////////////////////////////// + // Override RedBlack specialisation + ////////////////////////////////////////////////////// + virtual void RedBlackSource(Matrix & _Matrix,const Field &src, Field &src_e,Field &src_o) + { + GridBase *grid = _Matrix.RedBlackGrid(); + GridBase *fgrid= _Matrix.Grid(); + + Field tmp(grid); + Field Mtmp(grid); + + pickCheckerboard(Even,src_e,src); + pickCheckerboard(Odd ,src_o,src); + + ///////////////////////////////////////////////////// + // src_o = (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); + + _Matrix.Mooee(tmp,src_o); // Extra factor of "m" in source from dumb choice of matrix norm. + } + virtual void RedBlackSolution(Matrix & _Matrix,const Field &sol_o, const Field &src_e_c,Field &sol) + { + GridBase *grid = _Matrix.RedBlackGrid(); + GridBase *fgrid= _Matrix.Grid(); + + Field tmp(grid); + Field sol_e(grid); + Field src_e(grid); + + src_e = src_e_c; // Const correctness + + /////////////////////////////////////////////////// + // sol_e = M_ee^-1 * ( src_e - Meo sol_o )... + /////////////////////////////////////////////////// + _Matrix.Meooe(sol_o,tmp); assert( tmp.checkerboard ==Even); + src_e = src_e-tmp; assert( src_e.checkerboard ==Even); + _Matrix.MooeeInv(src_e,sol_e); assert( sol_e.checkerboard ==Even); + + setCheckerboard(sol,sol_e); assert( sol_e.checkerboard ==Even); + setCheckerboard(sol,sol_o); assert( sol_o.checkerboard ==Odd ); + } + virtual void RedBlackSolve (Matrix & _Matrix,const Field &src_o, Field &sol_o) + { + SchurStaggeredOperator _HermOpEO(_Matrix); + this->_HermitianRBSolver(_HermOpEO,src_o,sol_o); assert(sol_o.checkerboard==Odd); + }; + virtual void RedBlackSolve (Matrix & _Matrix,const std::vector &src_o, std::vector &sol_o) + { + SchurStaggeredOperator _HermOpEO(_Matrix); + this->_HermitianRBSolver(_HermOpEO,src_o,sol_o); + } + }; + template using SchurRedBlackStagSolve = SchurRedBlackStaggeredSolve; + + /////////////////////////////////////////////////////////////////////////////////////////////////////// + // Site diagonal has Mooee on it. + /////////////////////////////////////////////////////////////////////////////////////////////////////// + template class SchurRedBlackDiagMooeeSolve : public SchurRedBlackBase { + public: + typedef CheckerBoardedSparseMatrixBase Matrix; + + SchurRedBlackDiagMooeeSolve(OperatorFunction &HermitianRBSolver, const bool initSubGuess = false) + : SchurRedBlackBase (HermitianRBSolver,initSubGuess) {}; + + + ////////////////////////////////////////////////////// + // Override RedBlack specialisation + ////////////////////////////////////////////////////// + virtual void RedBlackSource(Matrix & _Matrix,const Field &src, Field &src_e,Field &src_o) + { + GridBase *grid = _Matrix.RedBlackGrid(); + GridBase *fgrid= _Matrix.Grid(); + + Field tmp(grid); + Field Mtmp(grid); + + pickCheckerboard(Even,src_e,src); + pickCheckerboard(Odd ,src_o,src); + + ///////////////////////////////////////////////////// + // 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 + SchurDiagMooeeOperator _HermOpEO(_Matrix); + _HermOpEO.MpcDag(tmp,src_o); assert(src_o.checkerboard ==Odd); + + } + virtual void RedBlackSolution(Matrix & _Matrix,const Field &sol_o, const Field &src_e,Field &sol) + { + GridBase *grid = _Matrix.RedBlackGrid(); + GridBase *fgrid= _Matrix.Grid(); + + Field tmp(grid); + Field sol_e(grid); + Field src_e_i(grid); + /////////////////////////////////////////////////// + // sol_e = M_ee^-1 * ( src_e - Meo sol_o )... + /////////////////////////////////////////////////// + _Matrix.Meooe(sol_o,tmp); assert( tmp.checkerboard ==Even); + src_e_i = src_e-tmp; assert( src_e_i.checkerboard ==Even); + _Matrix.MooeeInv(src_e_i,sol_e); assert( sol_e.checkerboard ==Even); + + setCheckerboard(sol,sol_e); assert( sol_e.checkerboard ==Even); + setCheckerboard(sol,sol_o); assert( sol_o.checkerboard ==Odd ); + } + virtual void RedBlackSolve (Matrix & _Matrix,const Field &src_o, Field &sol_o) + { + SchurDiagMooeeOperator _HermOpEO(_Matrix); + this->_HermitianRBSolver(_HermOpEO,src_o,sol_o); assert(sol_o.checkerboard==Odd); + }; + virtual void RedBlackSolve (Matrix & _Matrix,const std::vector &src_o, std::vector &sol_o) + { + SchurDiagMooeeOperator _HermOpEO(_Matrix); + this->_HermitianRBSolver(_HermOpEO,src_o,sol_o); + } + }; + + /////////////////////////////////////////////////////////////////////////////////////////////////////// + // Site diagonal is identity, right preconditioned by Mee^inv + // ( 1 - Meo Moo^inv Moe Mee^inv ) phi =( 1 - Meo Moo^inv Moe Mee^inv ) Mee psi = = eta = eta + //=> psi = MeeInv phi + /////////////////////////////////////////////////////////////////////////////////////////////////////// + template class SchurRedBlackDiagTwoSolve : public SchurRedBlackBase { + public: + typedef CheckerBoardedSparseMatrixBase Matrix; ///////////////////////////////////////////////////// // Wrap the usual normal equations Schur trick ///////////////////////////////////////////////////// - SchurRedBlackDiagTwoMixed(LinearFunction &HermitianRBSolver, const bool initSubGuess = false) : - _HermitianRBSolver(HermitianRBSolver) - { - CBfactorise=0; - subtractGuess(initSubGuess); - }; - void subtractGuess(const bool initSubGuess) - { - subGuess = initSubGuess; - } - bool isSubtractGuess(void) - { - return subGuess; - } + SchurRedBlackDiagTwoSolve(OperatorFunction &HermitianRBSolver, const bool initSubGuess = false) + : SchurRedBlackBase(HermitianRBSolver,initSubGuess) {}; - template - void operator() (Matrix & _Matrix,const Field &in, Field &out){ - ZeroGuesser guess; - (*this)(_Matrix,in,out,guess); - } - template - void operator() (Matrix & _Matrix,const Field &in, Field &out,Guesser &guess){ - - // FIXME CGdiagonalMee not implemented virtual function - // FIXME use CBfactorise to control schur decomp + virtual void RedBlackSource(Matrix & _Matrix,const Field &src, Field &src_e,Field &src_o) + { 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); + pickCheckerboard(Even,src_e,src); + pickCheckerboard(Odd ,src_o,src); ///////////////////////////////////////////////////// // src_o = Mdag * (source_o - Moe MeeInv source_e) @@ -461,43 +430,44 @@ namespace Grid { // get the right MpcDag _HermOpEO.MpcDag(tmp,src_o); assert(src_o.checkerboard ==Odd); + } - ////////////////////////////////////////////////////////////// - // Call the red-black solver - ////////////////////////////////////////////////////////////// - std::cout< _HermOpEO(_Matrix); + this->_HermitianRBSolver(_HermOpEO,src_o,sol_o); + }; + virtual void RedBlackSolve (Matrix & _Matrix,const std::vector &src_o, std::vector &sol_o) + { + SchurDiagTwoOperator _HermOpEO(_Matrix); + this->_HermitianRBSolver(_HermOpEO,src_o,sol_o); + } }; - } #endif diff --git a/Grid/communicator/Communicator_mpi3.cc b/Grid/communicator/Communicator_mpi3.cc index 424b7973..e6800afe 100644 --- a/Grid/communicator/Communicator_mpi3.cc +++ b/Grid/communicator/Communicator_mpi3.cc @@ -50,8 +50,6 @@ void CartesianCommunicator::Init(int *argc, char ***argv) assert(0); } - Grid_quiesce_nodes(); - // Never clean up as done once. MPI_Comm_dup (MPI_COMM_WORLD,&communicator_world); @@ -124,10 +122,8 @@ CartesianCommunicator::CartesianCommunicator(const std::vector &processors, // split the communicator ////////////////////////////////////////////////////////////////////////////////////////////////////// // int Nparent = parent._processors ; - // std::cout << " splitting from communicator "< &processors, int Nchild = Nparent/childsize; assert (childsize * Nchild == Nparent); - // std::cout << " child size "< ccoor(_ndimension); // coor within subcommunicator std::vector scoor(_ndimension); // coor of split within parent std::vector ssize(_ndimension); // coor of split within parent diff --git a/Grid/communicator/SharedMemoryMPI.cc b/Grid/communicator/SharedMemoryMPI.cc index b425d97b..d3b094ad 100644 --- a/Grid/communicator/SharedMemoryMPI.cc +++ b/Grid/communicator/SharedMemoryMPI.cc @@ -413,7 +413,7 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags) assert(((uint64_t)ptr&0x3F)==0); close(fd); WorldShmCommBufs[r] =ptr; - std::cout << "Set WorldShmCommBufs["< seeds; - std::stringstream sha; seeds = GridChecksum::sha256_seeds(s); - for(int i=0;i &seeds){ diff --git a/Grid/lattice/Lattice_transfer.h b/Grid/lattice/Lattice_transfer.h index f988f310..4eca50d6 100644 --- a/Grid/lattice/Lattice_transfer.h +++ b/Grid/lattice/Lattice_transfer.h @@ -464,8 +464,10 @@ void InsertSliceLocal(const Lattice &lowDim, Lattice & higherDim,int assert(orthog>=0); for(int d=0;d_processors[d] == hg->_processors[d]); - assert(lg->_ldimensions[d] == hg->_ldimensions[d]); + if ( d!=orthog ) { + assert(lg->_processors[d] == hg->_processors[d]); + assert(lg->_ldimensions[d] == hg->_ldimensions[d]); + } } // the above should guarantee that the operations are local @@ -499,8 +501,10 @@ void ExtractSliceLocal(Lattice &lowDim, Lattice & higherDim,int slic assert(orthog>=0); for(int d=0;d_processors[d] == hg->_processors[d]); - assert(lg->_ldimensions[d] == hg->_ldimensions[d]); + if ( d!=orthog ) { + assert(lg->_processors[d] == hg->_processors[d]); + assert(lg->_ldimensions[d] == hg->_ldimensions[d]); + } } // the above should guarantee that the operations are local diff --git a/Grid/log/Log.h b/Grid/log/Log.h index b58c5d16..5d97ee5a 100644 --- a/Grid/log/Log.h +++ b/Grid/log/Log.h @@ -146,9 +146,11 @@ public: if ( log.timestamp ) { log.StopWatch->Stop(); GridTime now = log.StopWatch->Elapsed(); + if ( log.timing_mode==1 ) log.StopWatch->Reset(); log.StopWatch->Start(); - stream << log.evidence()<< std::setw(6)< munge; BinaryIO::readLatticeObject< vobj, sobj >(field, filename, munge, offset, format,nersc_csum,scidac_csuma,scidac_csumb); - + std::cout << GridLogMessage << "SciDAC checksum A " << std::hex << scidac_csuma << std::dec << std::endl; + std::cout << GridLogMessage << "SciDAC checksum B " << std::hex << scidac_csumb << std::dec << std::endl; ///////////////////////////////////////////// // Insist checksum is next record ///////////////////////////////////////////// diff --git a/Grid/perfmon/Timer.h b/Grid/perfmon/Timer.h index 6fc421b8..07c5febd 100644 --- a/Grid/perfmon/Timer.h +++ b/Grid/perfmon/Timer.h @@ -49,21 +49,35 @@ inline double usecond(void) { typedef std::chrono::system_clock GridClock; typedef std::chrono::time_point GridTimePoint; -typedef std::chrono::milliseconds GridMillisecs; -typedef std::chrono::microseconds GridTime; -typedef std::chrono::microseconds GridUsecs; -inline std::ostream& operator<< (std::ostream & stream, const std::chrono::milliseconds & time) +typedef std::chrono::seconds GridSecs; +typedef std::chrono::milliseconds GridMillisecs; +typedef std::chrono::microseconds GridUsecs; +typedef std::chrono::microseconds GridTime; + +inline std::ostream& operator<< (std::ostream & stream, const GridSecs & time) { - stream << time.count()<<" ms"; + stream << time.count()<<" s"; return stream; } -inline std::ostream& operator<< (std::ostream & stream, const std::chrono::microseconds & time) +inline std::ostream& operator<< (std::ostream & stream, const GridMillisecs & now) { - stream << time.count()<<" usec"; + GridSecs second(1); + auto secs = now/second ; + auto subseconds = now%second ; + stream << secs<<"."< using iSinglet = iScalar > >; - template using iSpinMatrix = iScalar, Ns> >; - template using iColourMatrix = iScalar > > ; - template using iSpinColourMatrix = iScalar, Ns> >; - template using iLorentzColourMatrix = iVector >, Nd > ; - template using iDoubleStoredColourMatrix = iVector >, Nds > ; - template using iSpinVector = iScalar, Ns> >; - template using iColourVector = iScalar > >; - template using iSpinColourVector = iScalar, Ns> >; - template using iHalfSpinVector = iScalar, Nhs> >; - template using iHalfSpinColourVector = iScalar, Nhs> >; + + template using iSinglet = iScalar > >; + template using iSpinMatrix = iScalar, Ns> >; + template using iColourMatrix = iScalar > > ; + template using iSpinColourMatrix = iScalar, Ns> >; + template using iLorentzColourMatrix = iVector >, Nd > ; + template using iDoubleStoredColourMatrix = iVector >, Nds > ; + template using iSpinVector = iScalar, Ns> >; + template using iColourVector = iScalar > >; + template using iSpinColourVector = iScalar, Ns> >; + template using iHalfSpinVector = iScalar, Nhs> >; + template using iHalfSpinColourVector = iScalar, Nhs> >; + template using iSpinColourSpinColourMatrix = iScalar, Ns>, Nc>, Ns> >; + template using iGparitySpinColourVector = iVector, Ns>, Ngp >; template using iGparityHalfSpinColourVector = iVector, Nhs>, Ngp >; @@ -127,10 +130,28 @@ namespace QCD { typedef iSpinColourMatrix SpinColourMatrix; typedef iSpinColourMatrix SpinColourMatrixF; typedef iSpinColourMatrix SpinColourMatrixD; - + typedef iSpinColourMatrix vSpinColourMatrix; typedef iSpinColourMatrix vSpinColourMatrixF; typedef iSpinColourMatrix vSpinColourMatrixD; + + // SpinColourSpinColour matrix + typedef iSpinColourSpinColourMatrix SpinColourSpinColourMatrix; + typedef iSpinColourSpinColourMatrix SpinColourSpinColourMatrixF; + typedef iSpinColourSpinColourMatrix SpinColourSpinColourMatrixD; + + typedef iSpinColourSpinColourMatrix vSpinColourSpinColourMatrix; + typedef iSpinColourSpinColourMatrix vSpinColourSpinColourMatrixF; + typedef iSpinColourSpinColourMatrix vSpinColourSpinColourMatrixD; + + // SpinColourSpinColour matrix + typedef iSpinColourSpinColourMatrix SpinColourSpinColourMatrix; + typedef iSpinColourSpinColourMatrix SpinColourSpinColourMatrixF; + typedef iSpinColourSpinColourMatrix SpinColourSpinColourMatrixD; + + typedef iSpinColourSpinColourMatrix vSpinColourSpinColourMatrix; + typedef iSpinColourSpinColourMatrix vSpinColourSpinColourMatrixF; + typedef iSpinColourSpinColourMatrix vSpinColourSpinColourMatrixD; // LorentzColour typedef iLorentzColourMatrix LorentzColourMatrix; @@ -229,6 +250,9 @@ namespace QCD { typedef Lattice LatticeSpinColourMatrixF; typedef Lattice LatticeSpinColourMatrixD; + typedef Lattice LatticeSpinColourSpinColourMatrix; + typedef Lattice LatticeSpinColourSpinColourMatrixF; + typedef Lattice LatticeSpinColourSpinColourMatrixD; typedef Lattice LatticeLorentzColourMatrix; typedef Lattice LatticeLorentzColourMatrixF; diff --git a/Grid/qcd/action/ActionParams.h b/Grid/qcd/action/ActionParams.h index d25b60a9..76868abc 100644 --- a/Grid/qcd/action/ActionParams.h +++ b/Grid/qcd/action/ActionParams.h @@ -44,12 +44,15 @@ namespace QCD { struct WilsonImplParams { bool overlapCommsCompute; + std::vector twist_n_2pi_L; std::vector boundary_phases; WilsonImplParams() : overlapCommsCompute(false) { boundary_phases.resize(Nd, 1.0); + twist_n_2pi_L.resize(Nd, 0.0); }; - WilsonImplParams(const std::vector phi) - : boundary_phases(phi), overlapCommsCompute(false) {} + WilsonImplParams(const std::vector phi) : boundary_phases(phi), overlapCommsCompute(false) { + twist_n_2pi_L.resize(Nd, 0.0); + } }; struct StaggeredImplParams { diff --git a/Grid/qcd/action/fermion/CayleyFermion5D.cc b/Grid/qcd/action/fermion/CayleyFermion5D.cc index d4ed5a7c..a95ea4a0 100644 --- a/Grid/qcd/action/fermion/CayleyFermion5D.cc +++ b/Grid/qcd/action/fermion/CayleyFermion5D.cc @@ -68,6 +68,26 @@ void CayleyFermion5D::ExportPhysicalFermionSolution(const FermionField &so ExtractSlice(exported4d, tmp, 0, 0); } template +void CayleyFermion5D::P(const FermionField &psi, FermionField &chi) +{ + int Ls= this->Ls; + chi=zero; + for(int s=0;s +void CayleyFermion5D::Pdag(const FermionField &psi, FermionField &chi) +{ + int Ls= this->Ls; + chi=zero; + for(int s=0;s void CayleyFermion5D::ExportPhysicalFermionSource(const FermionField &solution5d,FermionField &exported4d) { int Ls = this->Ls; @@ -465,9 +485,13 @@ void CayleyFermion5D::SetCoefficientsInternal(RealD zolo_hi,std::vector _gamma; + RealD _zolo_hi; + RealD _b; + RealD _c; + // Cayley form Moebius (tanh and zolotarev) std::vector omega; std::vector bs; // S dependent coeffs diff --git a/Grid/qcd/action/fermion/Fermion.h b/Grid/qcd/action/fermion/Fermion.h index 2a008cb7..77a4681f 100644 --- a/Grid/qcd/action/fermion/Fermion.h +++ b/Grid/qcd/action/fermion/Fermion.h @@ -80,12 +80,24 @@ Author: Peter Boyle /////////////////////////////////////////////////////////////////////////////// #include +/////////////////////////////////////////////////////////////////////////////// +// Fourier accelerated Pauli Villars inverse support +/////////////////////////////////////////////////////////////////////////////// +#include + +//////////////////////////////////////////////////////////////////////////////// +// Move this group to a DWF specific tools/algorithms subdir? +//////////////////////////////////////////////////////////////////////////////// +#include +#include +#include +#include + //////////////////////////////////////////////////////////////////////////////////////////////////// // More maintainable to maintain the following typedef list centrally, as more "impl" targets // are added, (e.g. extension for gparity, half precision project in comms etc..) //////////////////////////////////////////////////////////////////////////////////////////////////// - // Cayley 5d namespace Grid { namespace QCD { diff --git a/Grid/qcd/action/fermion/FermionOperator.h b/Grid/qcd/action/fermion/FermionOperator.h index dc1ab924..281c077e 100644 --- a/Grid/qcd/action/fermion/FermionOperator.h +++ b/Grid/qcd/action/fermion/FermionOperator.h @@ -64,11 +64,6 @@ namespace Grid { virtual RealD M (const FermionField &in, FermionField &out)=0; virtual RealD Mdag (const FermionField &in, FermionField &out)=0; - // Query the even even properties to make algorithmic decisions - virtual int ConstEE(void) { return 1; }; // clover returns zero as EE depends on gauge field - virtual int isTrivialEE(void) { return 0; }; - virtual RealD Mass(void) {return 0.0;}; - // half checkerboard operaions virtual void Meooe (const FermionField &in, FermionField &out)=0; virtual void MeooeDag (const FermionField &in, FermionField &out)=0; diff --git a/Grid/qcd/action/fermion/FermionOperatorImpl.h b/Grid/qcd/action/fermion/FermionOperatorImpl.h index e6f6e136..721004e1 100644 --- a/Grid/qcd/action/fermion/FermionOperatorImpl.h +++ b/Grid/qcd/action/fermion/FermionOperatorImpl.h @@ -141,6 +141,7 @@ namespace QCD { //////////////////////////////////////////////////////////////////////// #define INHERIT_FIMPL_TYPES(Impl)\ + typedef Impl Impl_t; \ typedef typename Impl::FermionField FermionField; \ typedef typename Impl::PropagatorField PropagatorField; \ typedef typename Impl::DoubledGaugeField DoubledGaugeField; \ @@ -239,16 +240,30 @@ namespace QCD { GaugeLinkField tmp(GaugeGrid); Lattice > coor(GaugeGrid); + //////////////////////////////////////////////////// + // apply any boundary phase or twists + //////////////////////////////////////////////////// for (int mu = 0; mu < Nd; mu++) { - auto pha = Params.boundary_phases[mu]; - scalar_type phase( real(pha),imag(pha) ); + ////////// boundary phase ///////////// + auto pha = Params.boundary_phases[mu]; + scalar_type phase( real(pha),imag(pha) ); - int Lmu = GaugeGrid->GlobalDimensions()[mu] - 1; + int L = GaugeGrid->GlobalDimensions()[mu]; + int Lmu = L - 1; LatticeCoordinate(coor, mu); U = PeekIndex(Umu, mu); + + // apply any twists + RealD theta = Params.twist_n_2pi_L[mu] * 2*M_PI / L; + if ( theta != 0.0) { + scalar_type twphase(::cos(theta),::sin(theta)); + U = twphase*U; + std::cout << GridLogMessage << " Twist ["<(Uds, tmp, mu); diff --git a/Grid/qcd/action/fermion/FourierAcceleratedPV.h b/Grid/qcd/action/fermion/FourierAcceleratedPV.h new file mode 100644 index 00000000..d6196eee --- /dev/null +++ b/Grid/qcd/action/fermion/FourierAcceleratedPV.h @@ -0,0 +1,237 @@ + + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/qcd/action/fermion/FourierAcceleratedPV.h + + Copyright (C) 2015 + +Author: Christoph Lehner (lifted with permission by Peter Boyle, brought back to Grid) +Author: Peter Boyle + + 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 */ +#pragma once +namespace Grid { +namespace QCD { + + template + void get_real_const_bc(M& m, RealD& _b, RealD& _c) { + ComplexD b,c; + b=m.bs[0]; + c=m.cs[0]; + std::cout << GridLogMessage << "b=" << b << ", c=" << c << std::endl; + for (size_t i=1;i +class FourierAcceleratedPV { + public: + + ConjugateGradient &cg; + M& dwfPV; + G& Umu; + GridCartesian* grid5D; + GridRedBlackCartesian* gridRB5D; + int group_in_s; + + FourierAcceleratedPV(M& _dwfPV, G& _Umu, ConjugateGradient &_cg, int _group_in_s = 2) + : dwfPV(_dwfPV), Umu(_Umu), cg(_cg), group_in_s(_group_in_s) + { + assert( dwfPV.FermionGrid()->_fdimensions[0] % (2*group_in_s) == 0); + grid5D = QCD::SpaceTimeGrid::makeFiveDimGrid(2*group_in_s, (GridCartesian*)Umu._grid); + gridRB5D = QCD::SpaceTimeGrid::makeFiveDimRedBlackGrid(2*group_in_s, (GridCartesian*)Umu._grid); + } + + void rotatePV(const Vi& _src, Vi& dst, bool forward) const { + + GridStopWatch gsw1, gsw2; + + typedef typename Vi::scalar_type Coeff_t; + int Ls = dst._grid->_fdimensions[0]; + + Vi _tmp(dst._grid); + double phase = M_PI / (double)Ls; + Coeff_t bzero(0.0,0.0); + + FFT theFFT((GridCartesian*)dst._grid); + + if (!forward) { + gsw1.Start(); + for (int s=0;s_fdimensions[0]; + + GridStopWatch gswT; + gswT.Start(); + + RealD b,c; + get_real_const_bc(dwfPV,b,c); + RealD M5 = dwfPV.M5; + + // U(true) Rightinv TMinv U(false) = Minv + + Vi _src_diag(_dst._grid); + Vi _src_diag_slice(dwfPV.GaugeGrid()); + Vi _dst_diag_slice(dwfPV.GaugeGrid()); + Vi _src_diag_slices(grid5D); + Vi _dst_diag_slices(grid5D); + Vi _dst_diag(_dst._grid); + + rotatePV(_src,_src_diag,false); + + // now do TM solves + Gamma G5(Gamma::Algebra::Gamma5); + + GridStopWatch gswA, gswB; + + gswA.Start(); + + typedef typename M::Impl_t Impl; + //WilsonTMFermion tm(x.Umu,*x.UGridF,*x.UrbGridF,0.0,0.0,solver_outer.parent.par.wparams_f); + std::vector vmass(grid5D->_fdimensions[0],0.0); + std::vector vmu(grid5D->_fdimensions[0],0.0); + + WilsonTMFermion5D tm(Umu,*grid5D,*gridRB5D, + *(GridCartesian*)dwfPV.GaugeGrid(), + *(GridRedBlackCartesian*)dwfPV.GaugeRedBlackGrid(), + vmass,vmu); + + //SchurRedBlackDiagTwoSolve sol(cg); + SchurRedBlackDiagMooeeSolve sol(cg); // same performance as DiagTwo + gswA.Stop(); + + gswB.Start(); + + for (int sgroup=0;sgroup + + 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 */ +#pragma once + +namespace Grid { +namespace QCD { + +template X=0> +inline void convert(const Fieldi &from,Fieldo &to) +{ + precisionChange(to,from); +} +template X=0> +inline void convert(const Fieldi &from,Fieldo &to) +{ + to=from; +} + +template +class MADWF +{ + private: + typedef typename Matrixo::FermionField FermionFieldo; + typedef typename Matrixi::FermionField FermionFieldi; + + PVinverter & PauliVillarsSolvero;// For the outer field + SchurSolver & SchurSolveri; // For the inner approx field + Guesser & Guesseri; // To deflate the inner approx solves + + Matrixo & Mato; // Action object for outer + Matrixi & Mati; // Action object for inner + + RealD target_resid; + int maxiter; + public: + + MADWF(Matrixo &_Mato, + Matrixi &_Mati, + PVinverter &_PauliVillarsSolvero, + SchurSolver &_SchurSolveri, + Guesser & _Guesseri, + RealD resid, + int _maxiter) : + + Mato(_Mato),Mati(_Mati), + SchurSolveri(_SchurSolveri), + PauliVillarsSolvero(_PauliVillarsSolvero),Guesseri(_Guesseri) + { + target_resid=resid; + maxiter =_maxiter; + }; + + void operator() (const FermionFieldo &src4,FermionFieldo &sol5) + { + std::cout << GridLogMessage<< " ************************************************" << std::endl; + std::cout << GridLogMessage<< " MADWF-like algorithm " << std::endl; + std::cout << GridLogMessage<< " ************************************************" << std::endl; + + FermionFieldi c0i(Mati.GaugeGrid()); // 4d + FermionFieldi y0i(Mati.GaugeGrid()); // 4d + FermionFieldo c0 (Mato.GaugeGrid()); // 4d + FermionFieldo y0 (Mato.GaugeGrid()); // 4d + + FermionFieldo A(Mato.FermionGrid()); // Temporary outer + FermionFieldo B(Mato.FermionGrid()); // Temporary outer + FermionFieldo b(Mato.FermionGrid()); // 5d source + + FermionFieldo c(Mato.FermionGrid()); // PVinv source; reused so store + FermionFieldo defect(Mato.FermionGrid()); // 5d source + + FermionFieldi ci(Mati.FermionGrid()); + FermionFieldi yi(Mati.FermionGrid()); + FermionFieldi xi(Mati.FermionGrid()); + FermionFieldi srci(Mati.FermionGrid()); + FermionFieldi Ai(Mati.FermionGrid()); + + RealD m=Mati.Mass(); + + /////////////////////////////////////// + //Import source, include Dminus factors + /////////////////////////////////////// + Mato.ImportPhysicalFermionSource(src4,b); + std::cout << GridLogMessage << " src4 " < + + 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 */ +#pragma once + +namespace Grid { +namespace QCD { + +template +class PauliVillarsSolverUnprec +{ + public: + ConjugateGradient & CG; + PauliVillarsSolverUnprec( ConjugateGradient &_CG) : CG(_CG){}; + + template + void operator() (Matrix &_Matrix,const Field &src,Field &sol) + { + RealD m = _Matrix.Mass(); + Field A (_Matrix.FermionGrid()); + + MdagMLinearOperator HermOp(_Matrix); + + _Matrix.SetMass(1.0); + _Matrix.Mdag(src,A); + CG(HermOp,A,sol); + _Matrix.SetMass(m); + }; +}; + +template +class PauliVillarsSolverRBprec +{ + public: + SchurSolverType & SchurSolver; + PauliVillarsSolverRBprec( SchurSolverType &_SchurSolver) : SchurSolver(_SchurSolver){}; + + template + void operator() (Matrix &_Matrix,const Field &src,Field &sol) + { + RealD m = _Matrix.Mass(); + Field A (_Matrix.FermionGrid()); + + _Matrix.SetMass(1.0); + SchurSolver(_Matrix,src,sol); + _Matrix.SetMass(m); + }; +}; + +template +class PauliVillarsSolverFourierAccel +{ + public: + GaugeField & Umu; + ConjugateGradient & CG; + + PauliVillarsSolverFourierAccel(GaugeField &_Umu,ConjugateGradient &_CG) : Umu(_Umu), CG(_CG) + { + }; + + template + void operator() (Matrix &_Matrix,const Field &src,Field &sol) + { + FourierAcceleratedPV faPV(_Matrix,Umu,CG) ; + faPV.pvInv(src,sol); + }; +}; + + +} +} diff --git a/Grid/qcd/action/fermion/Reconstruct5Dprop.h b/Grid/qcd/action/fermion/Reconstruct5Dprop.h new file mode 100644 index 00000000..6862c5ee --- /dev/null +++ b/Grid/qcd/action/fermion/Reconstruct5Dprop.h @@ -0,0 +1,135 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/algorithms/iterative/SchurRedBlack.h + + Copyright (C) 2015 + +Author: Peter Boyle + + 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 */ +#pragma once + +namespace Grid { +namespace QCD { + +template class Reconstruct5DfromPhysical { + private: + PVinverter & PauliVillarsSolver; + public: + + ///////////////////////////////////////////////////// + // First cut works, 10 Oct 2018. + // + // Must form a plan to get this into production for Zmobius acceleration + // of the Mobius exact AMA corrections. + // + // TODO : understand absence of contact term in eqns in Hantao's thesis + // sol4 is contact term subtracted, but thesis & Brower's paper suggests not. + // + // Step 1: Localise PV inverse in a routine. [DONE] + // Step 2: Schur based PV inverse [DONE] + // Step 3: Fourier accelerated PV inverse [DONE] + // + ///////////////////////////////////////////////////// + + Reconstruct5DfromPhysical(PVinverter &_PauliVillarsSolver) + : PauliVillarsSolver(_PauliVillarsSolver) + { + }; + + + template + void PV(Matrix &_Matrix,const Field &src,Field &sol) + { + RealD m = _Matrix.Mass(); + _Matrix.SetMass(1.0); + _Matrix.M(src,sol); + _Matrix.SetMass(m); + } + template + void PVdag(Matrix &_Matrix,const Field &src,Field &sol) + { + RealD m = _Matrix.Mass(); + _Matrix.SetMass(1.0); + _Matrix.Mdag(src,sol); + _Matrix.SetMass(m); + } + template + void operator() (Matrix & _Matrix,const Field &sol4,const Field &src4, Field &sol5){ + + int Ls = _Matrix.Ls; + + Field psi4(_Matrix.GaugeGrid()); + Field psi(_Matrix.FermionGrid()); + Field A (_Matrix.FermionGrid()); + Field B (_Matrix.FermionGrid()); + Field c (_Matrix.FermionGrid()); + + typedef typename Matrix::Coeff_t Coeff_t; + + std::cout << GridLogMessage<< " ************************************************" << std::endl; + std::cout << GridLogMessage<< " Reconstruct5Dprop: c.f. MADWF algorithm " << std::endl; + std::cout << GridLogMessage<< " ************************************************" << std::endl; + + /////////////////////////////////////// + //Import source, include Dminus factors + /////////////////////////////////////// + _Matrix.ImportPhysicalFermionSource(src4,B); + + /////////////////////////////////////// + // Set up c from src4 + /////////////////////////////////////// + PauliVillarsSolver(_Matrix,B,A); + _Matrix.Pdag(A,c); + + ////////////////////////////////////// + // Build Pdag PV^-1 Dm P [-sol4,c2,c3... cL] + ////////////////////////////////////// + psi4 = - sol4; + InsertSlice(psi4, psi, 0 , 0); + for (int s=1;s ; NB Christoph did similar in GPT + + 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 */ +#pragma once + +#include +#include + + +namespace Grid { + + namespace QCD { + + template + class WilsonTMFermion5D : public WilsonFermion5D + { + public: + INHERIT_IMPL_TYPES(Impl); + public: + + virtual void Instantiatable(void) {}; + + // Constructors + WilsonTMFermion5D(GaugeField &_Umu, + GridCartesian &Fgrid, + GridRedBlackCartesian &Frbgrid, + GridCartesian &Ugrid, + GridRedBlackCartesian &Urbgrid, + const std::vector _mass, + const std::vector _mu, + const ImplParams &p= ImplParams() + ) : + WilsonFermion5D(_Umu, + Fgrid, + Frbgrid, + Ugrid, + Urbgrid, + 4.0,p) + + { + update(_mass,_mu); + } + + virtual void Meooe(const FermionField &in, FermionField &out) { + if (in.checkerboard == Odd) { + this->DhopEO(in, out, DaggerNo); + } else { + this->DhopOE(in, out, DaggerNo); + } + } + + virtual void MeooeDag(const FermionField &in, FermionField &out) { + if (in.checkerboard == Odd) { + this->DhopEO(in, out, DaggerYes); + } else { + this->DhopOE(in, out, DaggerYes); + } + } + + // allow override for twisted mass and clover + virtual void Mooee(const FermionField &in, FermionField &out) { + out.checkerboard = in.checkerboard; + //axpibg5x(out,in,a,b); // out = a*in + b*i*G5*in + for (int s=0;s<(int)this->mass.size();s++) { + ComplexD a = 4.0+this->mass[s]; + ComplexD b(0.0,this->mu[s]); + axpbg5y_ssp(out,a,in,b,in,s,s); + } + } + + virtual void MooeeDag(const FermionField &in, FermionField &out) { + out.checkerboard = in.checkerboard; + for (int s=0;s<(int)this->mass.size();s++) { + ComplexD a = 4.0+this->mass[s]; + ComplexD b(0.0,-this->mu[s]); + axpbg5y_ssp(out,a,in,b,in,s,s); + } + } + virtual void MooeeInv(const FermionField &in, FermionField &out) { + for (int s=0;s<(int)this->mass.size();s++) { + RealD m = this->mass[s]; + RealD tm = this->mu[s]; + RealD mtil = 4.0+this->mass[s]; + RealD sq = mtil*mtil+tm*tm; + ComplexD a = mtil/sq; + ComplexD b(0.0, -tm /sq); + axpbg5y_ssp(out,a,in,b,in,s,s); + } + } + virtual void MooeeInvDag(const FermionField &in, FermionField &out) { + for (int s=0;s<(int)this->mass.size();s++) { + RealD m = this->mass[s]; + RealD tm = this->mu[s]; + RealD mtil = 4.0+this->mass[s]; + RealD sq = mtil*mtil+tm*tm; + ComplexD a = mtil/sq; + ComplexD b(0.0,tm /sq); + axpbg5y_ssp(out,a,in,b,in,s,s); + } + } + + virtual RealD M(const FermionField &in, FermionField &out) { + out.checkerboard = in.checkerboard; + this->Dhop(in, out, DaggerNo); + FermionField tmp(out._grid); + for (int s=0;s<(int)this->mass.size();s++) { + ComplexD a = 4.0+this->mass[s]; + ComplexD b(0.0,this->mu[s]); + axpbg5y_ssp(tmp,a,in,b,in,s,s); + } + return axpy_norm(out, 1.0, tmp, out); + } + + // needed for fast PV + void update(const std::vector& _mass, const std::vector& _mu) { + assert(_mass.size() == _mu.size()); + assert(_mass.size() == this->FermionGrid()->_fdimensions[0]); + this->mass = _mass; + this->mu = _mu; + } + + private: + std::vector mu; + std::vector mass; + + }; + + typedef WilsonTMFermion5D WilsonTMFermion5DF; + typedef WilsonTMFermion5D WilsonTMFermion5DD; + +}} diff --git a/Grid/util/Sha.h b/Grid/util/Sha.h index 6cfbe0bd..b0a8cc10 100644 --- a/Grid/util/Sha.h +++ b/Grid/util/Sha.h @@ -38,7 +38,21 @@ public: { return ::crc32(0L,(unsigned char *)data,bytes); } - static inline std::vector sha256(void *data,size_t bytes) + template + static inline std::string sha256_string(const std::vector &hash) + { + std::stringstream sha; + std::string s; + + for(unsigned int i = 0; i < hash.size(); i++) + { + sha << std::hex << static_cast(hash[i]); + } + s = sha.str(); + + return s; + } + static inline std::vector sha256(const void *data,size_t bytes) { std::vector hash(SHA256_DIGEST_LENGTH); SHA256_CTX sha256; diff --git a/Hadrons/A2AVectors.hpp b/Hadrons/A2AVectors.hpp index 423a5f08..f55eb6d7 100644 --- a/Hadrons/A2AVectors.hpp +++ b/Hadrons/A2AVectors.hpp @@ -36,7 +36,7 @@ See the full license in the file "LICENSE" in the top level distribution directo BEGIN_HADRONS_NAMESPACE /****************************************************************************** - * Classes to generate V & W all-to-all vectors * + * Class to generate V & W all-to-all vectors * ******************************************************************************/ template class A2AVectorsSchurDiagTwo @@ -70,6 +70,42 @@ private: SchurDiagTwoOperator op_; }; +/****************************************************************************** + * Methods for V & W all-to-all vectors I/O * + ******************************************************************************/ +class A2AVectorsIo +{ +public: + struct Record: Serializable + { + GRID_SERIALIZABLE_CLASS_MEMBERS(Record, + unsigned int, index); + Record(void): index(0) {} + }; +public: + template + static void write(const std::string fileStem, std::vector &vec, + const bool multiFile, const int trajectory = -1); + template + static void read(std::vector &vec, const std::string fileStem, + const bool multiFile, const int trajectory = -1); +private: + static inline std::string vecFilename(const std::string stem, const int traj, + const bool multiFile) + { + std::string t = (traj < 0) ? "" : ("." + std::to_string(traj)); + + if (multiFile) + { + return stem + t; + } + else + { + return stem + t + ".bin"; + } + } +}; + /****************************************************************************** * A2AVectorsSchurDiagTwo template implementation * ******************************************************************************/ @@ -217,6 +253,90 @@ void A2AVectorsSchurDiagTwo::makeHighModeW5D(FermionField &wout_4d, } } +/****************************************************************************** + * all-to-all vectors I/O template implementation * + ******************************************************************************/ +template +void A2AVectorsIo::write(const std::string fileStem, std::vector &vec, + const bool multiFile, const int trajectory) +{ + Record record; + GridBase *grid = vec[0]._grid; + ScidacWriter binWriter(grid->IsBoss()); + std::string filename = vecFilename(fileStem, trajectory, multiFile); + + if (multiFile) + { + std::string fullFilename; + + for (unsigned int i = 0; i < vec.size(); ++i) + { + fullFilename = filename + "/elem" + std::to_string(i) + ".bin"; + + LOG(Message) << "Writing vector " << i << std::endl; + makeFileDir(fullFilename, grid); + binWriter.open(fullFilename); + record.index = i; + binWriter.writeScidacFieldRecord(vec[i], record); + binWriter.close(); + } + } + else + { + makeFileDir(filename, grid); + binWriter.open(filename); + for (unsigned int i = 0; i < vec.size(); ++i) + { + LOG(Message) << "Writing vector " << i << std::endl; + record.index = i; + binWriter.writeScidacFieldRecord(vec[i], record); + } + binWriter.close(); + } +} + +template +void A2AVectorsIo::read(std::vector &vec, const std::string fileStem, + const bool multiFile, const int trajectory) +{ + Record record; + ScidacReader binReader; + std::string filename = vecFilename(fileStem, trajectory, multiFile); + + if (multiFile) + { + std::string fullFilename; + + for (unsigned int i = 0; i < vec.size(); ++i) + { + fullFilename = filename + "/elem" + std::to_string(i) + ".bin"; + + LOG(Message) << "Reading vector " << i << std::endl; + binReader.open(fullFilename); + binReader.readScidacFieldRecord(vec[i], record); + binReader.close(); + if (record.index != i) + { + HADRONS_ERROR(Io, "vector index mismatch"); + } + } + } + else + { + binReader.open(filename); + for (unsigned int i = 0; i < vec.size(); ++i) + { + LOG(Message) << "Reading vector " << i << std::endl; + binReader.readScidacFieldRecord(vec[i], record); + if (record.index != i) + { + HADRONS_ERROR(Io, "vector index mismatch"); + } + } + binReader.close(); + } +} + END_HADRONS_NAMESPACE #endif // A2A_Vectors_hpp_ diff --git a/Hadrons/DilutedNoise.hpp b/Hadrons/DilutedNoise.hpp index 1f8186ae..1eb7ff40 100644 --- a/Hadrons/DilutedNoise.hpp +++ b/Hadrons/DilutedNoise.hpp @@ -7,6 +7,7 @@ Source file: Hadrons/DilutedNoise.hpp Copyright (C) 2015-2018 Author: Antonin Portelli +Author: Vera Guelpers 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 @@ -76,6 +77,22 @@ private: unsigned int nt_; }; +template +class FullVolumeSpinColorDiagonalNoise: public DilutedNoise +{ +public: + typedef typename FImpl::FermionField FermionField; +public: + // constructor/destructor + FullVolumeSpinColorDiagonalNoise(GridCartesian *g, unsigned int n_src); + virtual ~FullVolumeSpinColorDiagonalNoise(void) = default; + // generate noise + virtual void generateNoise(GridParallelRNG &rng); +private: + unsigned int nSrc_; +}; + + /****************************************************************************** * DilutedNoise template implementation * ******************************************************************************/ @@ -186,6 +203,47 @@ void TimeDilutedSpinColorDiagonalNoise::generateNoise(GridParallelRNG &rn } } +/****************************************************************************** + * FullVolumeSpinColorDiagonalNoise template implementation * + ******************************************************************************/ +template +FullVolumeSpinColorDiagonalNoise:: +FullVolumeSpinColorDiagonalNoise(GridCartesian *g, unsigned int nSrc) +: DilutedNoise(g, nSrc*Ns*FImpl::Dimension), nSrc_(nSrc) +{} + +template +void FullVolumeSpinColorDiagonalNoise::generateNoise(GridParallelRNG &rng) +{ + typedef decltype(peekColour((*this)[0], 0)) SpinField; + + auto &noise = *this; + auto g = this->getGrid(); + auto nd = g->GlobalDimensions().size(); + auto nc = FImpl::Dimension; + Complex shift(1., 1.); + LatticeComplex eta(g); + SpinField etas(g); + unsigned int i = 0; + + bernoulli(rng, eta); + eta = (2.*eta - shift)*(1./::sqrt(2.)); + for (unsigned int n = 0; n < nSrc_; ++n) + { + for (unsigned int s = 0; s < Ns; ++s) + { + etas = zero; + pokeSpin(etas, eta, s); + for (unsigned int c = 0; c < nc; ++c) + { + noise[i] = zero; + pokeColour(noise[i], etas, c); + i++; + } + } + } +} + END_HADRONS_NAMESPACE #endif // Hadrons_DilutedNoise_hpp_ diff --git a/Hadrons/DiskVector.hpp b/Hadrons/DiskVector.hpp index 231854f3..0b4bffe1 100644 --- a/Hadrons/DiskVector.hpp +++ b/Hadrons/DiskVector.hpp @@ -34,6 +34,12 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include +#ifdef DV_DEBUG +#define DV_DEBUG_MSG(dv, stream) LOG(Debug) << "diskvector " << (dv) << ": " << stream << std::endl +#else +#define DV_DEBUG_MSG(dv, stream) +#endif + BEGIN_HADRONS_NAMESPACE /****************************************************************************** @@ -56,9 +62,7 @@ public: // write to disk and cache T &operator=(const T &obj) const { -#ifdef DV_DEBUG - LOG(Debug) << "diskvector " << &master_ << ": writing to " << i_ << std::endl; -#endif + DV_DEBUG_MSG(&master_, "writing to " << i_); master_.cacheInsert(i_, obj); master_.save(master_.filename(i_), obj); @@ -82,6 +86,8 @@ public: virtual ~DiskVectorBase(void); const T & operator[](const unsigned int i) const; RwAccessHelper operator[](const unsigned int i); + double hitRatio(void) const; + void resetStat(void); private: virtual void load(T &obj, const std::string filename) const = 0; virtual void save(const std::string filename, const T &obj) const = 0; @@ -93,6 +99,7 @@ private: private: std::string dirname_; unsigned int size_, cacheSize_; + double access_{0.}, hit_{0.}; bool clean_; // using pointers to allow modifications when class is const // semantic: const means data unmodified, but cache modification allowed @@ -115,6 +122,7 @@ private: read(reader, basename(filename), obj); } + virtual void save(const std::string filename, const T &obj) const { Writer writer(filename); @@ -123,13 +131,76 @@ private: } }; +/****************************************************************************** + * Specialisation for Eigen matrices * + ******************************************************************************/ +template +using EigenDiskVectorMat = Eigen::Matrix; + +template +class EigenDiskVector: public DiskVectorBase> +{ +public: + using DiskVectorBase>::DiskVectorBase; + typedef EigenDiskVectorMat Matrix; +public: + T operator()(const unsigned int i, const Eigen::Index j, + const Eigen::Index k) const + { + return (*this)[i](j, k); + } +private: + virtual void load(EigenDiskVectorMat &obj, const std::string filename) const + { + std::ifstream f(filename, std::ios::binary); + std::vector hash(SHA256_DIGEST_LENGTH); + Eigen::Index nRow, nCol; + size_t matSize; + double t; + + f.read(reinterpret_cast(hash.data()), hash.size()*sizeof(unsigned char)); + f.read(reinterpret_cast(&nRow), sizeof(Eigen::Index)); + f.read(reinterpret_cast(&nCol), sizeof(Eigen::Index)); + obj.resize(nRow, nCol); + matSize = nRow*nCol*sizeof(T); + t = -usecond(); + f.read(reinterpret_cast(obj.data()), matSize); + t += usecond(); + DV_DEBUG_MSG(this, "Eigen read " << matSize/t*1.0e6/1024/1024 << " MB/s"); + auto check = GridChecksum::sha256(obj.data(), matSize); + DV_DEBUG_MSG(this, "Eigen sha256 " << GridChecksum::sha256_string(check)); + if (hash != check) + { + HADRONS_ERROR(Io, "checksum failed") + } + } + + virtual void save(const std::string filename, const EigenDiskVectorMat &obj) const + { + std::ofstream f(filename, std::ios::binary); + std::vector hash(SHA256_DIGEST_LENGTH); + Eigen::Index nRow, nCol; + size_t matSize; + double t; + + nRow = obj.rows(); + nCol = obj.cols(); + matSize = nRow*nCol*sizeof(T); + hash = GridChecksum::sha256(obj.data(), matSize); + DV_DEBUG_MSG(this, "Eigen sha256 " << GridChecksum::sha256_string(hash)); + f.write(reinterpret_cast(hash.data()), hash.size()*sizeof(unsigned char)); + f.write(reinterpret_cast(&nRow), sizeof(Eigen::Index)); + f.write(reinterpret_cast(&nCol), sizeof(Eigen::Index)); + t = -usecond(); + f.write(reinterpret_cast(obj.data()), matSize); + t += usecond(); + DV_DEBUG_MSG(this, "Eigen write " << matSize/t*1.0e6/1024/1024 << " MB/s"); + } +}; + /****************************************************************************** * DiskVectorBase implementation * ******************************************************************************/ -#ifdef DV_DEBUG -#define DV_DEBUG_MSG(stream) LOG(Debug) << "diskvector " << this << ": " << stream << std::endl -#endif - template DiskVectorBase::DiskVectorBase(const std::string dirname, const unsigned int size, @@ -160,28 +231,29 @@ DiskVectorBase::~DiskVectorBase(void) template const T & DiskVectorBase::operator[](const unsigned int i) const { - auto &cache = *cachePtr_; - auto &loads = *loadsPtr_; + auto &cache = *cachePtr_; + auto &loads = *loadsPtr_; - DV_DEBUG_MSG("accessing " << i << " (RO)"); + DV_DEBUG_MSG(this, "accessing " << i << " (RO)"); if (i >= size_) { HADRONS_ERROR(Size, "index out of range"); } - + const_cast(access_)++; if (cache.find(i) == cache.end()) { // cache miss - DV_DEBUG_MSG("cache miss"); + DV_DEBUG_MSG(this, "cache miss"); fetch(i); } else { - DV_DEBUG_MSG("cache hit"); + DV_DEBUG_MSG(this, "cache hit"); auto pos = std::find(loads.begin(), loads.end(), i); + const_cast(hit_)++; loads.erase(pos); loads.push_back(i); } @@ -193,7 +265,7 @@ const T & DiskVectorBase::operator[](const unsigned int i) const { msg += std::to_string(p) + " "; } - DV_DEBUG_MSG("in cache: " << msg); + DV_DEBUG_MSG(this, "in cache: " << msg); #endif return cache.at(i); @@ -202,7 +274,7 @@ const T & DiskVectorBase::operator[](const unsigned int i) const template typename DiskVectorBase::RwAccessHelper DiskVectorBase::operator[](const unsigned int i) { - DV_DEBUG_MSG("accessing " << i << " (RW)"); + DV_DEBUG_MSG(this, "accessing " << i << " (RW)"); if (i >= size_) { @@ -212,6 +284,19 @@ typename DiskVectorBase::RwAccessHelper DiskVectorBase::operator[](const u return RwAccessHelper(*this, i); } +template +double DiskVectorBase::hitRatio(void) const +{ + return hit_/access_; +} + +template +void DiskVectorBase::resetStat(void) +{ + access_ = 0.; + hit_ = 0.; +} + template std::string DiskVectorBase::filename(const unsigned int i) const { @@ -226,7 +311,7 @@ void DiskVectorBase::evict(void) const if (cache.size() >= cacheSize_) { - DV_DEBUG_MSG("evicting " << loads.front()); + DV_DEBUG_MSG(this, "evicting " << loads.front()); cache.erase(loads.front()); loads.pop_front(); } @@ -239,7 +324,7 @@ void DiskVectorBase::fetch(const unsigned int i) const auto &loads = *loadsPtr_; struct stat s; - DV_DEBUG_MSG("loading " << i << " from disk"); + DV_DEBUG_MSG(this, "loading " << i << " from disk"); evict(); if(stat(filename(i).c_str(), &s) != 0) @@ -267,7 +352,7 @@ void DiskVectorBase::cacheInsert(const unsigned int i, const T &obj) const { msg += std::to_string(p) + " "; } - DV_DEBUG_MSG("in cache: " << msg); + DV_DEBUG_MSG(this, "in cache: " << msg); #endif } diff --git a/Hadrons/Modules.hpp b/Hadrons/Modules.hpp index 2859f0fc..da289307 100644 --- a/Hadrons/Modules.hpp +++ b/Hadrons/Modules.hpp @@ -16,6 +16,7 @@ #include #include #include +#include #include #include #include @@ -23,13 +24,17 @@ #include #include #include +#include #include #include #include +#include #include +#include #include #include #include +#include #include #include #include @@ -40,6 +45,9 @@ #include #include #include +#include +#include +#include #include #include #include @@ -50,7 +58,6 @@ #include #include #include -#include #include #include #include @@ -61,6 +68,7 @@ #include #include #include +#include #include #include #include diff --git a/Hadrons/Modules/MAction/DWF.cc b/Hadrons/Modules/MAction/DWF.cc index 6cda84ac..38d25cb9 100644 --- a/Hadrons/Modules/MAction/DWF.cc +++ b/Hadrons/Modules/MAction/DWF.cc @@ -32,4 +32,6 @@ using namespace Hadrons; using namespace MAction; template class Grid::Hadrons::MAction::TDWF; +#ifdef GRID_DEFAULT_PRECISION_DOUBLE template class Grid::Hadrons::MAction::TDWF; +#endif diff --git a/Hadrons/Modules/MAction/DWF.hpp b/Hadrons/Modules/MAction/DWF.hpp index 67cfeb1b..a7104b42 100644 --- a/Hadrons/Modules/MAction/DWF.hpp +++ b/Hadrons/Modules/MAction/DWF.hpp @@ -49,7 +49,8 @@ public: unsigned int, Ls, double , mass, double , M5, - std::string , boundary); + std::string , boundary, + std::string , twist); }; template @@ -73,7 +74,9 @@ protected: }; MODULE_REGISTER_TMP(DWF, TDWF, MAction); +#ifdef GRID_DEFAULT_PRECISION_DOUBLE MODULE_REGISTER_TMP(DWFF, TDWF, MAction); +#endif /****************************************************************************** * DWF template implementation * @@ -117,8 +120,9 @@ void TDWF::setup(void) auto &grb4 = *envGetRbGrid(FermionField); auto &g5 = *envGetGrid(FermionField, par().Ls); auto &grb5 = *envGetRbGrid(FermionField, par().Ls); - std::vector boundary = strToVec(par().boundary); - typename DomainWallFermion::ImplParams implParams(boundary); + typename DomainWallFermion::ImplParams implParams; + implParams.boundary_phases = strToVec(par().boundary); + implParams.twist_n_2pi_L = strToVec(par().twist); envCreateDerived(FMat, DomainWallFermion, getName(), par().Ls, U, g5, grb5, g4, grb4, par().mass, par().M5, implParams); } diff --git a/Hadrons/Modules/MAction/MobiusDWF.cc b/Hadrons/Modules/MAction/MobiusDWF.cc index 41683920..879452d8 100644 --- a/Hadrons/Modules/MAction/MobiusDWF.cc +++ b/Hadrons/Modules/MAction/MobiusDWF.cc @@ -32,4 +32,6 @@ using namespace Hadrons; using namespace MAction; template class Grid::Hadrons::MAction::TMobiusDWF; +#ifdef GRID_DEFAULT_PRECISION_DOUBLE template class Grid::Hadrons::MAction::TMobiusDWF; +#endif diff --git a/Hadrons/Modules/MAction/MobiusDWF.hpp b/Hadrons/Modules/MAction/MobiusDWF.hpp index 22964c9a..0ba9c4c3 100644 --- a/Hadrons/Modules/MAction/MobiusDWF.hpp +++ b/Hadrons/Modules/MAction/MobiusDWF.hpp @@ -49,7 +49,8 @@ public: double , M5, double , b, double , c, - std::string , boundary); + std::string , boundary, + std::string , twist); }; template @@ -72,7 +73,9 @@ public: }; MODULE_REGISTER_TMP(MobiusDWF, TMobiusDWF, MAction); +#ifdef GRID_DEFAULT_PRECISION_DOUBLE MODULE_REGISTER_TMP(MobiusDWFF, TMobiusDWF, MAction); +#endif /****************************************************************************** * TMobiusDWF implementation * @@ -117,8 +120,9 @@ void TMobiusDWF::setup(void) auto &grb4 = *envGetRbGrid(FermionField); auto &g5 = *envGetGrid(FermionField, par().Ls); auto &grb5 = *envGetRbGrid(FermionField, par().Ls); - std::vector boundary = strToVec(par().boundary); - typename MobiusFermion::ImplParams implParams(boundary); + typename MobiusFermion::ImplParams implParams; + implParams.boundary_phases = strToVec(par().boundary); + implParams.twist_n_2pi_L = strToVec(par().twist); envCreateDerived(FMat, MobiusFermion, getName(), par().Ls, U, g5, grb5, g4, grb4, par().mass, par().M5, par().b, par().c, implParams); diff --git a/Hadrons/Modules/MAction/ScaledDWF.cc b/Hadrons/Modules/MAction/ScaledDWF.cc index 70eafac5..7008bf5d 100644 --- a/Hadrons/Modules/MAction/ScaledDWF.cc +++ b/Hadrons/Modules/MAction/ScaledDWF.cc @@ -32,4 +32,6 @@ using namespace Hadrons; using namespace MAction; template class Grid::Hadrons::MAction::TScaledDWF; +#ifdef GRID_DEFAULT_PRECISION_DOUBLE template class Grid::Hadrons::MAction::TScaledDWF; +#endif diff --git a/Hadrons/Modules/MAction/ScaledDWF.hpp b/Hadrons/Modules/MAction/ScaledDWF.hpp index c890f0e9..3b8066e6 100644 --- a/Hadrons/Modules/MAction/ScaledDWF.hpp +++ b/Hadrons/Modules/MAction/ScaledDWF.hpp @@ -48,7 +48,8 @@ public: double , mass, double , M5, double , scale, - std::string , boundary); + std::string , boundary, + std::string , twist); }; template @@ -71,7 +72,9 @@ public: }; MODULE_REGISTER_TMP(ScaledDWF, TScaledDWF, MAction); +#ifdef GRID_DEFAULT_PRECISION_DOUBLE MODULE_REGISTER_TMP(ScaledDWFF, TScaledDWF, MAction); +#endif /****************************************************************************** * TScaledDWF implementation * @@ -116,8 +119,9 @@ void TScaledDWF::setup(void) auto &grb4 = *envGetRbGrid(FermionField); auto &g5 = *envGetGrid(FermionField, par().Ls); auto &grb5 = *envGetRbGrid(FermionField, par().Ls); - std::vector boundary = strToVec(par().boundary); - typename MobiusFermion::ImplParams implParams(boundary); + typename ScaledShamirFermion::ImplParams implParams; + implParams.boundary_phases = strToVec(par().boundary); + implParams.twist_n_2pi_L = strToVec(par().twist); envCreateDerived(FMat, ScaledShamirFermion, getName(), par().Ls, U, g5, grb5, g4, grb4, par().mass, par().M5, par().scale, implParams); diff --git a/Hadrons/Modules/MAction/Wilson.cc b/Hadrons/Modules/MAction/Wilson.cc index 4ce3e7ef..1e801ed6 100644 --- a/Hadrons/Modules/MAction/Wilson.cc +++ b/Hadrons/Modules/MAction/Wilson.cc @@ -32,4 +32,6 @@ using namespace Hadrons; using namespace MAction; template class Grid::Hadrons::MAction::TWilson; +#ifdef GRID_DEFAULT_PRECISION_DOUBLE template class Grid::Hadrons::MAction::TWilson; +#endif diff --git a/Hadrons/Modules/MAction/Wilson.hpp b/Hadrons/Modules/MAction/Wilson.hpp index a9327c2f..b5e53837 100644 --- a/Hadrons/Modules/MAction/Wilson.hpp +++ b/Hadrons/Modules/MAction/Wilson.hpp @@ -47,7 +47,9 @@ public: GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonPar, std::string, gauge, double , mass, - std::string, boundary); + std::string, boundary, + std::string, string, + std::string, twist); }; template @@ -71,7 +73,9 @@ protected: }; MODULE_REGISTER_TMP(Wilson, TWilson, MAction); +#ifdef GRID_DEFAULT_PRECISION_DOUBLE MODULE_REGISTER_TMP(WilsonF, TWilson, MAction); +#endif /****************************************************************************** * TWilson template implementation * @@ -111,8 +115,9 @@ void TWilson::setup(void) auto &U = envGet(GaugeField, par().gauge); auto &grid = *envGetGrid(FermionField); auto &gridRb = *envGetRbGrid(FermionField); - std::vector boundary = strToVec(par().boundary); - typename WilsonFermion::ImplParams implParams(boundary); + typename WilsonFermion::ImplParams implParams; + implParams.boundary_phases = strToVec(par().boundary); + implParams.twist_n_2pi_L = strToVec(par().twist); envCreateDerived(FMat, WilsonFermion, getName(), 1, U, grid, gridRb, par().mass, implParams); } diff --git a/Hadrons/Modules/MAction/WilsonClover.cc b/Hadrons/Modules/MAction/WilsonClover.cc index 2c5c0e66..eed1582c 100644 --- a/Hadrons/Modules/MAction/WilsonClover.cc +++ b/Hadrons/Modules/MAction/WilsonClover.cc @@ -32,4 +32,6 @@ using namespace Hadrons; using namespace MAction; template class Grid::Hadrons::MAction::TWilsonClover; +#ifdef GRID_DEFAULT_PRECISION_DOUBLE template class Grid::Hadrons::MAction::TWilsonClover; +#endif diff --git a/Hadrons/Modules/MAction/WilsonClover.hpp b/Hadrons/Modules/MAction/WilsonClover.hpp index 349abe84..ad301380 100644 --- a/Hadrons/Modules/MAction/WilsonClover.hpp +++ b/Hadrons/Modules/MAction/WilsonClover.hpp @@ -51,7 +51,8 @@ public: double , csw_r, double , csw_t, WilsonAnisotropyCoefficients ,clover_anisotropy, - std::string, boundary + std::string, boundary, + std::string, twist ); }; @@ -75,7 +76,9 @@ public: }; MODULE_REGISTER_TMP(WilsonClover, TWilsonClover, MAction); +#ifdef GRID_DEFAULT_PRECISION_DOUBLE MODULE_REGISTER_TMP(WilsonCloverF, TWilsonClover, MAction); +#endif /****************************************************************************** * TWilsonClover template implementation * @@ -117,8 +120,9 @@ void TWilsonClover::setup(void) auto &U = envGet(GaugeField, par().gauge); auto &grid = *envGetGrid(FermionField); auto &gridRb = *envGetRbGrid(FermionField); - std::vector boundary = strToVec(par().boundary); - typename WilsonCloverFermion::ImplParams implParams(boundary); + typename WilsonCloverFermion::ImplParams implParams; + implParams.boundary_phases = strToVec(par().boundary); + implParams.twist_n_2pi_L = strToVec(par().twist); envCreateDerived(FMat, WilsonCloverFermion, getName(), 1, U, grid, gridRb, par().mass, par().csw_r, par().csw_t, par().clover_anisotropy, implParams); diff --git a/Hadrons/Modules/MAction/ZMobiusDWF.cc b/Hadrons/Modules/MAction/ZMobiusDWF.cc index ef8e4799..609b76cc 100644 --- a/Hadrons/Modules/MAction/ZMobiusDWF.cc +++ b/Hadrons/Modules/MAction/ZMobiusDWF.cc @@ -32,4 +32,6 @@ using namespace Hadrons; using namespace MAction; template class Grid::Hadrons::MAction::TZMobiusDWF; +#ifdef GRID_DEFAULT_PRECISION_DOUBLE template class Grid::Hadrons::MAction::TZMobiusDWF; +#endif diff --git a/Hadrons/Modules/MAction/ZMobiusDWF.hpp b/Hadrons/Modules/MAction/ZMobiusDWF.hpp index f7959127..12ad82ea 100644 --- a/Hadrons/Modules/MAction/ZMobiusDWF.hpp +++ b/Hadrons/Modules/MAction/ZMobiusDWF.hpp @@ -50,7 +50,8 @@ public: double , b, double , c, std::vector>, omega, - std::string , boundary); + std::string , boundary, + std::string , twist); }; template @@ -73,7 +74,9 @@ public: }; MODULE_REGISTER_TMP(ZMobiusDWF, TZMobiusDWF, MAction); +#ifdef GRID_DEFAULT_PRECISION_DOUBLE MODULE_REGISTER_TMP(ZMobiusDWFF, TZMobiusDWF, MAction); +#endif /****************************************************************************** * TZMobiusDWF implementation * @@ -125,8 +128,9 @@ void TZMobiusDWF::setup(void) auto &g5 = *envGetGrid(FermionField, par().Ls); auto &grb5 = *envGetRbGrid(FermionField, par().Ls); auto omega = par().omega; - std::vector boundary = strToVec(par().boundary); - typename ZMobiusFermion::ImplParams implParams(boundary); + typename ZMobiusFermion::ImplParams implParams; + implParams.boundary_phases = strToVec(par().boundary); + implParams.twist_n_2pi_L = strToVec(par().twist); envCreateDerived(FMat, ZMobiusFermion, getName(), par().Ls, U, g5, grb5, g4, grb4, par().mass, par().M5, omega, par().b, par().c, implParams); diff --git a/Hadrons/Modules/MContraction/A2AAslashField.cc b/Hadrons/Modules/MContraction/A2AAslashField.cc index eb76932e..65c49198 100644 --- a/Hadrons/Modules/MContraction/A2AAslashField.cc +++ b/Hadrons/Modules/MContraction/A2AAslashField.cc @@ -1,3 +1,30 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MContraction/A2AAslashField.cc + +Copyright (C) 2015-2018 + +Author: Antonin Portelli + +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 using namespace Grid; diff --git a/Hadrons/Modules/MContraction/A2AAslashField.hpp b/Hadrons/Modules/MContraction/A2AAslashField.hpp index 792ff6ee..cc8b1559 100644 --- a/Hadrons/Modules/MContraction/A2AAslashField.hpp +++ b/Hadrons/Modules/MContraction/A2AAslashField.hpp @@ -1,3 +1,30 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MContraction/A2AAslashField.hpp + +Copyright (C) 2015-2018 + +Author: Antonin Portelli + +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 */ #ifndef Hadrons_MContraction_A2AAslashField_hpp_ #define Hadrons_MContraction_A2AAslashField_hpp_ diff --git a/Hadrons/Modules/MGauge/Electrify.cc b/Hadrons/Modules/MGauge/Electrify.cc new file mode 100644 index 00000000..1feea9ec --- /dev/null +++ b/Hadrons/Modules/MGauge/Electrify.cc @@ -0,0 +1,34 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MGauge/Electrify.cc + +Copyright (C) 2015-2018 + +Author: Vera Guelpers + +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 + +using namespace Grid; +using namespace Hadrons; +using namespace MGauge; + +template class Grid::Hadrons::MGauge::TElectrify; diff --git a/Hadrons/Modules/MGauge/Electrify.hpp b/Hadrons/Modules/MGauge/Electrify.hpp new file mode 100644 index 00000000..58d65eba --- /dev/null +++ b/Hadrons/Modules/MGauge/Electrify.hpp @@ -0,0 +1,151 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MGauge/Electrify.hpp + +Copyright (C) 2015-2018 + +Author: Vera Guelpers + +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 */ + +#ifndef Hadrons_MGauge_Electrify_hpp_ +#define Hadrons_MGauge_Electrify_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Electrify gauge * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MGauge) + +/**************************************************************************** +* Electrify a gauge field: +* +* Ue_mu(x) = U_mu(x)*exp(ieqA_mu(x)) +* +* with +* +* - gauge: U_mu(x): gauge field +* - emField: A_mu(x): electromagnetic photon field +* - e: value for the elementary charge +* - q: charge in units of e +* +*****************************************************************************/ + + +class ElectrifyPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(ElectrifyPar, + std::string, gauge, + std::string, emField, + double, e, + double, charge); +}; + +template +class TElectrify: public Module +{ +public: + GAUGE_TYPE_ALIASES(GImpl,); +public: + typedef PhotonR::GaugeField EmField; +public: + // constructor + TElectrify(const std::string name); + // destructor + virtual ~TElectrify(void) {}; + // dependencies/products + virtual std::vector getInput(void); + virtual std::vector getOutput(void); +protected: + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_TMP(Electrify, TElectrify, MGauge); + +/****************************************************************************** +* TElectrify implementation * +******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TElectrify::TElectrify(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TElectrify::getInput(void) +{ + std::vector in = {par().gauge, par().emField}; + + return in; +} + +template +std::vector TElectrify::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TElectrify::setup(void) +{ + envCreateLat(GaugeField, getName()); + envTmpLat(LatticeComplex, "eiAmu"); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TElectrify::execute(void) +{ + LOG(Message) << "Electrify the gauge field " << par().gauge << " using the photon field " + << par().emField << " with charge e*q= " << par().e << "*" << par().charge << std::endl; + + auto &Ue = envGet(GaugeField, getName()); + auto &U = envGet(GaugeField, par().gauge); + auto &A = envGet(EmField, par().emField); + envGetTmp(LatticeComplex, eiAmu); + + Complex i(0.0,1.0); + + for(unsigned int mu = 0; mu < env().getNd(); mu++) + { + eiAmu = exp(i * (Real)(par().e * par().charge) * PeekIndex(A, mu)); + PokeIndex(Ue, PeekIndex(U, mu) * eiAmu, mu); + } +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MGauge_Electrify_hpp_ diff --git a/Hadrons/Modules/MScalarSUN/TimeMomProbe.cc b/Hadrons/Modules/MGauge/GaugeFix.cc similarity index 68% rename from Hadrons/Modules/MScalarSUN/TimeMomProbe.cc rename to Hadrons/Modules/MGauge/GaugeFix.cc index 5bfa203a..53aa16da 100644 --- a/Hadrons/Modules/MScalarSUN/TimeMomProbe.cc +++ b/Hadrons/Modules/MGauge/GaugeFix.cc @@ -2,11 +2,12 @@ Grid physics library, www.github.com/paboyle/Grid -Source file: Hadrons/Modules/MScalarSUN/TimeMomProbe.cc +Source file: Hadrons/Modules/MGauge/GaugeFix.cc Copyright (C) 2015-2018 Author: Antonin Portelli +Author: Peter Boyle 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 @@ -25,14 +26,11 @@ with this program; if not, write to the Free Software Foundation, Inc., See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ -#include + +#include using namespace Grid; using namespace Hadrons; -using namespace MScalarSUN; +using namespace MGauge; -template class Grid::Hadrons::MScalarSUN::TTimeMomProbe>; -template class Grid::Hadrons::MScalarSUN::TTimeMomProbe>; -template class Grid::Hadrons::MScalarSUN::TTimeMomProbe>; -template class Grid::Hadrons::MScalarSUN::TTimeMomProbe>; -template class Grid::Hadrons::MScalarSUN::TTimeMomProbe>; +template class Grid::Hadrons::MGauge::TGaugeFix; diff --git a/Hadrons/Modules/MGauge/GaugeFix.hpp b/Hadrons/Modules/MGauge/GaugeFix.hpp new file mode 100644 index 00000000..ece8c19d --- /dev/null +++ b/Hadrons/Modules/MGauge/GaugeFix.hpp @@ -0,0 +1,135 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MGauge/GaugeFix.hpp + +Copyright (C) 2015-2018 + +Author: Antonin Portelli +Author: Peter Boyle + +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 */ + +#ifndef Hadrons_MGaugeFix_hpp_ +#define Hadrons_MGaugeFix_hpp_ + +#include +#include +#include +#include +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Fix gauge * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MGauge) + +class GaugeFixPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(GaugeFixPar, + std::string, gauge, + Real, alpha, + int, maxiter, + Real, Omega_tol, + Real, Phi_tol, + bool, Fourier); +}; + +template +class TGaugeFix: public Module +{ +public: + GAUGE_TYPE_ALIASES(GImpl,); +public: + // constructor + TGaugeFix(const std::string name); + // destructor + virtual ~TGaugeFix(void) {}; + // dependencies/products + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_TMP(GaugeFix, TGaugeFix, MGauge); + +/****************************************************************************** +* TGaugeFix implementation * +******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TGaugeFix::TGaugeFix(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TGaugeFix::getInput(void) +{ + std::vector in = {par().gauge}; + return in; +} + +template +std::vector TGaugeFix::getOutput(void) +{ + std::vector out = {getName()}; + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TGaugeFix::setup(void) +{ + envCreateLat(GaugeField, getName()); +} + + +// execution /////////////////////////////////////////////////////////////////// +template +void TGaugeFix::execute(void) +//Loads the gauge and fixes it +{ + std::cout << "executing" << std::endl; + LOG(Message) << "Fixing the Gauge" << std::endl; + LOG(Message) << par().gauge << std::endl; + auto &U = envGet(GaugeField, par().gauge); + auto &Umu = envGet(GaugeField, getName()); + LOG(Message) << "Gauge Field fetched" << std::endl; + //do we allow maxiter etc to be user set? + Real alpha = par().alpha; + int maxiter = par().maxiter; + Real Omega_tol = par().Omega_tol; + Real Phi_tol = par().Phi_tol; + bool Fourier = par().Fourier; + FourierAcceleratedGaugeFixer::SteepestDescentGaugeFix(U,alpha,maxiter,Omega_tol,Phi_tol,Fourier); + Umu = U; + LOG(Message) << "Gauge Fixed" << std::endl; +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MGaugeFix_hpp_ diff --git a/Hadrons/Modules/MIO/LoadA2AVectors.cc b/Hadrons/Modules/MIO/LoadA2AVectors.cc new file mode 100644 index 00000000..7a40a6f5 --- /dev/null +++ b/Hadrons/Modules/MIO/LoadA2AVectors.cc @@ -0,0 +1,34 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MIO/LoadA2AVectors.cc + +Copyright (C) 2015-2018 + +Author: Antonin Portelli + +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 + +using namespace Grid; +using namespace Hadrons; +using namespace MIO; + +template class Grid::Hadrons::MIO::TLoadA2AVectors; diff --git a/Hadrons/Modules/MIO/LoadA2AVectors.hpp b/Hadrons/Modules/MIO/LoadA2AVectors.hpp new file mode 100644 index 00000000..5b194c16 --- /dev/null +++ b/Hadrons/Modules/MIO/LoadA2AVectors.hpp @@ -0,0 +1,120 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MIO/LoadA2AVectors.hpp + +Copyright (C) 2015-2018 + +Author: Antonin Portelli + +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 */ +#ifndef Hadrons_MIO_LoadA2AVectors_hpp_ +#define Hadrons_MIO_LoadA2AVectors_hpp_ + +#include +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Module to load all-to-all vectors * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MIO) + +class LoadA2AVectorsPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(LoadA2AVectorsPar, + std::string, filestem, + bool, multiFile, + unsigned int, size); +}; + +template +class TLoadA2AVectors: public Module +{ +public: + FERM_TYPE_ALIASES(FImpl,); +public: + // constructor + TLoadA2AVectors(const std::string name); + // destructor + virtual ~TLoadA2AVectors(void) {}; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_TMP(LoadA2AVectors, TLoadA2AVectors, MIO); + +/****************************************************************************** + * TLoadA2AVectors implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TLoadA2AVectors::TLoadA2AVectors(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TLoadA2AVectors::getInput(void) +{ + std::vector in; + + return in; +} + +template +std::vector TLoadA2AVectors::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TLoadA2AVectors::setup(void) +{ + envCreate(std::vector, getName(), 1, par().size, + envGetGrid(FermionField)); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TLoadA2AVectors::execute(void) +{ + auto &vec = envGet(std::vector, getName()); + + A2AVectorsIo::read(vec, par().filestem, par().multiFile, vm().getTrajectory()); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MIO_LoadA2AVectors_hpp_ diff --git a/Hadrons/Modules/MIO/LoadEigenPack.cc b/Hadrons/Modules/MIO/LoadEigenPack.cc index f0ae63c3..28fdeb01 100644 --- a/Hadrons/Modules/MIO/LoadEigenPack.cc +++ b/Hadrons/Modules/MIO/LoadEigenPack.cc @@ -32,4 +32,6 @@ using namespace Hadrons; using namespace MIO; template class Grid::Hadrons::MIO::TLoadEigenPack>; +#ifdef GRID_DEFAULT_PRECISION_DOUBLE template class Grid::Hadrons::MIO::TLoadEigenPack>; +#endif diff --git a/Hadrons/Modules/MIO/LoadEigenPack.hpp b/Hadrons/Modules/MIO/LoadEigenPack.hpp index 2302b15d..016675c9 100644 --- a/Hadrons/Modules/MIO/LoadEigenPack.hpp +++ b/Hadrons/Modules/MIO/LoadEigenPack.hpp @@ -72,7 +72,9 @@ public: }; MODULE_REGISTER_TMP(LoadFermionEigenPack, TLoadEigenPack>, MIO); +#ifdef GRID_DEFAULT_PRECISION_DOUBLE MODULE_REGISTER_TMP(LoadFermionEigenPackIo32, ARG(TLoadEigenPack>), MIO); +#endif /****************************************************************************** * TLoadEigenPack implementation * diff --git a/Hadrons/Modules/MNPR/Amputate.cc b/Hadrons/Modules/MNPR/Amputate.cc new file mode 100644 index 00000000..ec7c5940 --- /dev/null +++ b/Hadrons/Modules/MNPR/Amputate.cc @@ -0,0 +1,36 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MNPR/Amputate.cc + +Copyright (C) 2015-2018 + +Author: Antonin Portelli +Author: Peter Boyle + +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 + +using namespace Grid; +using namespace Hadrons; +using namespace MNPR; + +template class Grid::Hadrons::MNPR::TAmputate; + diff --git a/Hadrons/Modules/MNPR/Amputate.hpp b/Hadrons/Modules/MNPR/Amputate.hpp new file mode 100644 index 00000000..953bb354 --- /dev/null +++ b/Hadrons/Modules/MNPR/Amputate.hpp @@ -0,0 +1,200 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MNPR/Amputate.hpp + +Copyright (C) 2015-2018 + +Author: Antonin Portelli +Author: Julia Kettle J.R.Kettle-2@sms.ed.ac.uk +Author: Peter Boyle + +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 */ + +#ifndef Hadrons_Amputate_hpp_ +#define Hadrons_Amputate_hpp_ + +#include +#include +#include +#include +//#include +//#include +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * TAmputate * + Performs bilinear contractions of the type tr[g5*adj(Sout)*g5*G*Sin] + Suitable for non exceptional momenta +******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MNPR) + +class AmputatePar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(AmputatePar, + std::string, Sin, //need to make this a propogator type? + std::string, Sout, //same + std::string, vertex, + std::string, pin, + std::string, pout, + std::string, output, + std::string, input); +}; + +template +class TAmputate: public Module +{ +public: + FERM_TYPE_ALIASES(FImpl1, 1); + FERM_TYPE_ALIASES(FImpl2, 2); + class Result: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Result, + std::vector, Vamp, + ); + }; +public: + // constructor + TAmputate(const std::string name); + // destructor + virtual ~TAmputate(void) {}; + // dependencies/products + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + virtual SpinColourMatrix invertspincolmat(SpinColourMatrix &scmat); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_TMP(Amputate, ARG(TAmputate), MNPR); + +/****************************************************************************** + * TAmputate implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TAmputate::TAmputate(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TAmputate::getInput(void) +{ + std::vector input = {par().Sin, par().Sout, par().vertex}; + + return input; +} + +template +std::vector TAmputate::getOutput(void) +{ + std::vector output = {getName()}; + + + return output; +} + +// Invert spin colour matrix using Eigen +template +SpinColourMatrix TAmputate::invertspincolmat(SpinColourMatrix &scmat) +{ + Eigen::MatrixXcf scmat_2d(Ns*Nc,Ns*Nc); + for(int ic=0; ic +void TAmputate::execute(void) +{ + LOG(Message) << "Computing bilinear amputations '" << getName() << "' using" + << " momentum '" << par().Sin << "' and '" << par().Sout << "'" + << std::endl; + BinaryWriter writer(par().output); + PropagatorField1 &Sin = *env().template getObject(par().Sin); //Do these have the phases taken into account?? Don't think so. FIX + PropagatorField2 &Sout = *env().template getObject(par().Sout); + std::vector pin = strToVec(par().pin), pout = strToVec(par().pout); + std::vector latt_size(pin.begin(), pin.end()); + LatticeComplex pdotxin(env().getGrid()), pdotxout(env().getGrid()), coor(env().getGrid()); + LOG(Message) << "Propagators set up " << std::endl; + std::vector vertex; // Let's read from file here + Gamma g5(Gamma::Algebra::Gamma5); + Result result; + LOG(Message) << "reading file - " << par().input << std::endl; + BinaryReader reader(par().input); + Complex Ci(0.0,1.0); + + std::string svertex; + read(reader,"vertex", vertex); + LOG(Message) << "vertex read" << std::endl; + + pdotxin=zero; + pdotxout=zero; + for (unsigned int mu = 0; mu < 4; ++mu) + { + Real TwoPiL = M_PI * 2.0/ latt_size[mu]; + LatticeCoordinate(coor,mu); + pdotxin = pdotxin +(TwoPiL * pin[mu]) * coor; + pdotxout= pdotxout +(TwoPiL * pout[mu]) * coor; + } + Sin = Sin*exp(-Ci*pdotxin); //phase corrections + Sout = Sout*exp(-Ci*pdotxout); + + SpinColourMatrix Sin_mom = sum(Sin); + SpinColourMatrix Sout_mom = sum(Sout); + LOG(Message) << "summed over lattice" << std::endl; + + LOG(Message) << "Lattice -> spincolourmatrix conversion" << std::endl; + + SpinColourMatrix Sin_inv = invertspincolmat(Sin_mom); + SpinColourMatrix Sout_inv = invertspincolmat(Sout_mom); + LOG(Message) << "Inversions done" << std::endl; + + result.Vamp.resize(Gamma::nGamma/2); + for( int mu=0; mu < Gamma::nGamma/2; mu++){ + Gamma::Algebra gam = mu; + result.Vamp[mu] = 1/12.0*trace(adj(Gamma(mu*2+1))*g5*Sout_inv*g5*vertex[mu]*Sin_inv); + LOG(Message) << "Vamp[" << mu << "] - " << result.Vamp[mu] << std::endl; + } +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_Amputate_hpp_ diff --git a/Hadrons/Modules/MNPR/Bilinear.cc b/Hadrons/Modules/MNPR/Bilinear.cc new file mode 100644 index 00000000..c5b38d37 --- /dev/null +++ b/Hadrons/Modules/MNPR/Bilinear.cc @@ -0,0 +1,36 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MNPR/Bilinear.cc + +Copyright (C) 2015-2018 + +Author: Antonin Portelli +Author: Peter Boyle + +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 + +using namespace Grid; +using namespace Hadrons; +using namespace MNPR; + +template class Grid::Hadrons::MNPR::TBilinear; + diff --git a/Hadrons/Modules/MNPR/Bilinear.hpp b/Hadrons/Modules/MNPR/Bilinear.hpp new file mode 100644 index 00000000..66cd22a6 --- /dev/null +++ b/Hadrons/Modules/MNPR/Bilinear.hpp @@ -0,0 +1,225 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MNPR/Bilinear.hpp + +Copyright (C) 2015-2018 + +Author: Antonin Portelli +Author: Julia Kettle J.R.Kettle-2@sms.ed.ac.uk +Author: Peter Boyle + +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 */ + +#ifndef Hadrons_Bilinear_hpp_ +#define Hadrons_Bilinear_hpp_ + +#include +#include +#include +#include +//#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * TBilinear * + Performs bilinear contractions of the type tr[g5*adj(Sout)*g5*G*Sin] + Suitable for non exceptional momenta in Rome-Southampton NPR +******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MNPR) + +class BilinearPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(BilinearPar, + std::string, Sin, + std::string, Sout, + std::string, pin, + std::string, pout, + std::string, output); +}; + +template +class TBilinear: public Module +{ +public: + FERM_TYPE_ALIASES(FImpl1, 1); + FERM_TYPE_ALIASES(FImpl2, 2); + class Result: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Result, + std::vector, bilinear); + }; +public: + // constructor + TBilinear(const std::string name); + // destructor + virtual ~TBilinear(void) {}; + // dependencies/products + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + //LatticeSpinColourMatrix PhaseProps(LatticeSpinColourMatrix S, std::vector p); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_TMP(Bilinear, ARG(TBilinear), MNPR); + +/****************************************************************************** + * TBilinear implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TBilinear::TBilinear(const std::string name) +: Module(name) +{} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TBilinear::setup(void) +{ + //env().template registerLattice(getName()); + //env().template registerObject(getName()); +} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TBilinear::getInput(void) +{ + std::vector input = {par().Sin, par().Sout}; + + return input; +} + +template +std::vector TBilinear::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +/* +/////Phase propagators////////////////////////// +template +LatticeSpinColourMatrix TBilinear::PhaseProps(LatticeSpinColourMatrix S, std::vector p) +{ + GridBase *grid = S._grid; + LatticeComplex pdotx(grid), coor(grid); + std::vector latt_size = grid->_fdimensions; + Complex Ci(0.0,1.0); + pdotx=zero; + for (unsigned int mu = 0; mu < 4; ++mu) + { + Real TwoPiL = M_PI * 2.0/ latt_size[mu]; + LatticeCoordinate(coor,mu); + pdotx = pdotx +(TwoPiL * p[mu]) * coor; + } + S = S*exp(-Ci*pdotx); + return S; +} +*/ +// execution /////////////////////////////////////////////////////////////////// +template +void TBilinear::execute(void) +{ +/************************************************************************** + +Compute the bilinear vertex needed for the NPR. +V(G) = sum_x [ g5 * adj(S'(x,p2)) * g5 * G * S'(x,p1) ]_{si,sj,ci,cj} +G is one of the 16 gamma vertices [I,gmu,g5,g5gmu,sig(mu,nu)] + + * G + / \ + p1/ \p2 + / \ + / \ + +Returns a spin-colour matrix, with indices si,sj, ci,cj + +Conventions: +p1 - incoming momenta +p2 - outgoing momenta +q = (p1-p2) +**************************************************************************/ + + LOG(Message) << "Computing bilinear contractions '" << getName() << "' using" + << " momentum '" << par().Sin << "' and '" << par().Sout << "'" + << std::endl; + + BinaryWriter writer(par().output); + + + // Propogators + LatticeSpinColourMatrix &Sin = *env().template getObject(par().Sin); + LatticeSpinColourMatrix &Sout = *env().template getObject(par().Sout); + LatticeComplex pdotxin(env().getGrid()), pdotxout(env().getGrid()), coor(env().getGrid()); + // momentum on legs + std::vector pin = strToVec(par().pin), pout = strToVec(par().pout); + std::vector latt_size(pin.begin(), pin.end()); + //bilinears + LatticeSpinColourMatrix bilinear_x(env().getGrid()); + SpinColourMatrix bilinear; + Gamma g5(Gamma::Algebra::Gamma5); + Result result; + Complex Ci(0.0,1.0); + + // + + pdotxin=zero; + pdotxout=zero; + for (unsigned int mu = 0; mu < 4; ++mu) + { + Real TwoPiL = M_PI * 2.0/ latt_size[mu]; + LatticeCoordinate(coor,mu); + pdotxin = pdotxin +(TwoPiL * pin[mu]) * coor; + pdotxout= pdotxout +(TwoPiL * pout[mu]) * coor; + } + Sin = Sin*exp(-Ci*pdotxin); //phase corrections + Sout = Sout*exp(-Ci*pdotxout); + + ////Set up gamma vector////////////////////////// + std::vector gammavector; + for( int i=0; i +Author: Peter Boyle + +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 + +using namespace Grid; +using namespace Hadrons; +using namespace MNPR; + +template class Grid::Hadrons::MNPR::TFourQuark; + diff --git a/Hadrons/Modules/MNPR/FourQuark.hpp b/Hadrons/Modules/MNPR/FourQuark.hpp new file mode 100644 index 00000000..f961a366 --- /dev/null +++ b/Hadrons/Modules/MNPR/FourQuark.hpp @@ -0,0 +1,274 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MNPR/FourQuark.hpp + +Copyright (C) 2015-2018 + +Author: Antonin Portelli +Author: Julia Kettle J.R.Kettle-2@sms.ed.ac.uk +Author: Peter Boyle + +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 */ + +#ifndef Hadrons_FourQuark_hpp_ +#define Hadrons_FourQuark_hpp_ + +#include +#include +#include +#include +#include +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * TFourQuark * + Performs fourquark contractions of the type tr[g5*adj(Sout)*g5*G*Sin] + Suitable for non exceptional momenta +******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MNPR) + +class FourQuarkPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(FourQuarkPar, + std::string, Sin, //need to make this a propogator type? + std::string, Sout, //same + std::string, pin, + std::string, pout, + bool, fullbasis, + std::string, output); +}; + +template +class TFourQuark: public Module +{ +public: + FERM_TYPE_ALIASES(FImpl1, 1); + FERM_TYPE_ALIASES(FImpl2, 2); + class Result: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Result, + std::vector, fourquark); + }; +public: + // constructor + TFourQuark(const std::string name); + // destructor + virtual ~TFourQuark(void) {}; + // dependencies/products + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void tensorprod(LatticeSpinColourSpinColourMatrix &lret, LatticeSpinColourMatrix a, LatticeSpinColourMatrix b); + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_TMP(FourQuark, ARG(TFourQuark), MNPR); + +/****************************************************************************** + * TFourQuark implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TFourQuark::TFourQuark(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TFourQuark::getInput(void) +{ + std::vector input = {par().Sin, par().Sout}; + + return input; +} + +template +std::vector TFourQuark::getOutput(void) +{ + std::vector output = {getName()}; + + return output; +} + + +template +void TFourQuark::tensorprod(LatticeSpinColourSpinColourMatrix &lret, LatticeSpinColourMatrix a, LatticeSpinColourMatrix b) +{ +#if 0 + parallel_for(auto site=lret.begin();site +void TFourQuark::setup(void) +{ + envCreateLat(LatticeSpinColourMatrix, getName()); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TFourQuark::execute(void) +{ + +/********************************************************************************* + +TFourQuark : Creates the four quark vertex required for the NPR of four-quark ops + +V_{Gamma_1,Gamma_2} = sum_x [ ( g5 * adj(S'(x,p2)) * g5 * G1 * S'(x,p1) )_ci,cj;si,sj x ( g5 * adj(S'(x,p2)) * g5 * G2 S'(x,p1) )_ck,cl;sk,cl ] + +Create a bilinear vertex for G1 and G2 the spin and colour indices are kept free. Where there are 16 potential Gs. +We then find the outer product of V1 and V2, keeping the spin and colour indices uncontracted +Then this is summed over the lattice coordinate +Result is a SpinColourSpinColourMatrix - with 4 colour and 4 spin indices. +We have up to 256 of these including the offdiag (G1 != G2). + + \ / + \p1 p1/ + \ / + \ / + G1 * * G2 + / \ + / \ + /p2 p2\ + / \ + +*********************************************************************************/ + + + + + LOG(Message) << "Computing fourquark contractions '" << getName() << "' using" + << " momentum '" << par().Sin << "' and '" << par().Sout << "'" + << std::endl; + + BinaryWriter writer(par().output); + + PropagatorField1 &Sin = *env().template getObject(par().Sin); + PropagatorField2 &Sout = *env().template getObject(par().Sout); + std::vector pin = strToVec(par().pin), pout = strToVec(par().pout); + bool fullbasis = par().fullbasis; + Gamma g5(Gamma::Algebra::Gamma5); + Result result; + std::vector latt_size(pin.begin(), pin.end()); + LatticeComplex pdotxin(env().getGrid()), pdotxout(env().getGrid()), coor(env().getGrid()); + LatticeSpinColourMatrix bilinear_mu(env().getGrid()), bilinear_nu(env().getGrid()); + LatticeSpinColourSpinColourMatrix lret(env().getGrid()); + Complex Ci(0.0,1.0); + + //Phase propagators + //Sin = Grid::QCD::PropUtils::PhaseProps(Sin,pin); + //Sout = Grid::QCD::PropUtils::PhaseProps(Sout,pout); + + //find p.x for in and out so phase can be accounted for in propagators + pdotxin=zero; + pdotxout=zero; + for (unsigned int mu = 0; mu < 4; ++mu) + { + Real TwoPiL = M_PI * 2.0/ latt_size[mu]; + LatticeCoordinate(coor,mu); + pdotxin = pdotxin +(TwoPiL * pin[mu]) * coor; + pdotxout= pdotxout +(TwoPiL * pout[mu]) * coor; + } + Sin = Sin*exp(-Ci*pdotxin); //phase corrections + Sout = Sout*exp(-Ci*pdotxout); + + + //Set up Gammas + std::vector gammavector; + for( int i=1; i +Author: Vera Guelpers + +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 + +using namespace Grid; +using namespace Hadrons; +using namespace MNoise; + +template class Grid::Hadrons::MNoise::TFullVolumeSpinColorDiagonal; +template class Grid::Hadrons::MNoise::TFullVolumeSpinColorDiagonal; diff --git a/Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.hpp b/Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.hpp new file mode 100644 index 00000000..93990882 --- /dev/null +++ b/Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.hpp @@ -0,0 +1,121 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.hpp + +Copyright (C) 2015-2018 + +Author: Antonin Portelli +Author: Vera Guelpers + +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 */ +#ifndef Hadrons_MNoise_FullVolumeSpinColorDiagonal_hpp_ +#define Hadrons_MNoise_FullVolumeSpinColorDiagonal_hpp_ + +#include +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Generate full volume spin-color diagonal noise * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MNoise) + +class FullVolumeSpinColorDiagonalPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(FullVolumeSpinColorDiagonalPar, + unsigned int, nsrc); +}; + +template +class TFullVolumeSpinColorDiagonal: public Module +{ +public: + FERM_TYPE_ALIASES(FImpl,); +public: + // constructor + TFullVolumeSpinColorDiagonal(const std::string name); + // destructor + virtual ~TFullVolumeSpinColorDiagonal(void) {}; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_TMP(FullVolumeSpinColorDiagonal, TFullVolumeSpinColorDiagonal, MNoise); +MODULE_REGISTER_TMP(ZFullVolumeSpinColorDiagonal, TFullVolumeSpinColorDiagonal, MNoise); + +/****************************************************************************** + * TFullVolumeSpinColorDiagonal implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TFullVolumeSpinColorDiagonal::TFullVolumeSpinColorDiagonal(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TFullVolumeSpinColorDiagonal::getInput(void) +{ + std::vector in; + + return in; +} + +template +std::vector TFullVolumeSpinColorDiagonal::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TFullVolumeSpinColorDiagonal::setup(void) +{ + envCreateDerived(DilutedNoise, + FullVolumeSpinColorDiagonalNoise, + getName(), 1, envGetGrid(FermionField), par().nsrc); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TFullVolumeSpinColorDiagonal::execute(void) +{ + auto &noise = envGet(DilutedNoise, getName()); + LOG(Message) << "Generating full volume, spin-color diagonal noise" << std::endl; + noise.generateNoise(rng4d()); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MNoise_FullVolumeSpinColorDiagonal_hpp_ diff --git a/Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp b/Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp deleted file mode 100644 index a73da8de..00000000 --- a/Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp +++ /dev/null @@ -1,268 +0,0 @@ -/************************************************************************************* - -Grid physics library, www.github.com/paboyle/Grid - -Source file: Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp - -Copyright (C) 2015-2018 - -Author: Antonin Portelli - -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 */ -#ifndef Hadrons_MScalarSUN_TimeMomProbe_hpp_ -#define Hadrons_MScalarSUN_TimeMomProbe_hpp_ - -#include -#include -#include -#include - -BEGIN_HADRONS_NAMESPACE - -/****************************************************************************** - * n-point functions O(t,p)*tr(phi(t_1,p_1)*...*phi(t_n,p_n)) * - ******************************************************************************/ -BEGIN_MODULE_NAMESPACE(MScalarSUN) - -class TimeMomProbePar: Serializable -{ -public: - GRID_SERIALIZABLE_CLASS_MEMBERS(TimeMomProbePar, - std::string, field, - std::vector, op, - std::vector>, timeMom, - std::string, output); -}; - -class TimeMomProbeResult: Serializable -{ -public: - GRID_SERIALIZABLE_CLASS_MEMBERS(TimeMomProbeResult, - std::string, op, - std::vector>, timeMom, - std::vector, data); -}; - -template -class TTimeMomProbe: public Module -{ -public: - typedef typename SImpl::Field Field; - typedef typename SImpl::SiteField::scalar_object Site; - typedef typename SImpl::ComplexField ComplexField; - typedef std::vector SlicedOp; -public: - // constructor - TTimeMomProbe(const std::string name); - // destructor - virtual ~TTimeMomProbe(void) {}; - // dependency relation - virtual std::vector getInput(void); - virtual std::vector getOutput(void); - // setup - virtual void setup(void); - // execution - virtual void execute(void); -private: - void vectorModulo(std::vector &v); -}; - -MODULE_REGISTER_TMP(TimeMomProbeSU2, TTimeMomProbe>, MScalarSUN); -MODULE_REGISTER_TMP(TimeMomProbeSU3, TTimeMomProbe>, MScalarSUN); -MODULE_REGISTER_TMP(TimeMomProbeSU4, TTimeMomProbe>, MScalarSUN); -MODULE_REGISTER_TMP(TimeMomProbeSU5, TTimeMomProbe>, MScalarSUN); -MODULE_REGISTER_TMP(TimeMomProbeSU6, TTimeMomProbe>, MScalarSUN); - -/****************************************************************************** - * TTimeMomProbe implementation * - ******************************************************************************/ -// constructor ///////////////////////////////////////////////////////////////// -template -TTimeMomProbe::TTimeMomProbe(const std::string name) -: Module(name) -{} - -// dependencies/products /////////////////////////////////////////////////////// -template -std::vector TTimeMomProbe::getInput(void) -{ - std::vector in = par().op; - - in.push_back(par().field); - - return in; -} - -template -std::vector TTimeMomProbe::getOutput(void) -{ - std::vector out; - - return out; -} - -// setup /////////////////////////////////////////////////////////////////////// -template -void TTimeMomProbe::setup(void) -{ - envTmpLat(ComplexField, "ftBuf"); - envTmpLat(Field, "ftMatBuf"); -} - -// execution /////////////////////////////////////////////////////////////////// -// NB: time is direction 0 -template -void TTimeMomProbe::vectorModulo(std::vector &v) -{ - for (unsigned int mu = 0; mu < env().getNd(); ++mu) - { - auto d = env().getDim(mu); - v[mu] = ((v[mu] % d) + d) % d; - } -} - -template -void TTimeMomProbe::execute(void) -{ - const unsigned int nd = env().getNd(); - const unsigned int nt = env().getDim(0); - double partVol = 1.; - std::set> timeMomSet; - std::vector>> timeMom; - std::vector> transferMom; - FFT fft(envGetGrid(Field)); - std::vector dMask(nd, 1); - std::vector result; - std::map> slicedOp; - std::vector slicedProbe; - auto &phi = envGet(Field, par().field); - - envGetTmp(ComplexField, ftBuf); - envGetTmp(Field, ftMatBuf); - dMask[0] = 0; - for (unsigned int mu = 1; mu < nd; ++mu) - { - partVol *= env().getDim(mu); - } - timeMom.resize(par().timeMom.size()); - for (unsigned int p = 0; p < timeMom.size(); ++p) - { - for (auto &tms: par().timeMom[p]) - { - std::vector tm = strToVec(tms); - - timeMom[p].push_back(tm); - timeMomSet.insert(tm); - } - transferMom.push_back(std::vector(nd - 1, 0)); - for (auto &tm: timeMom[p]) - { - for (unsigned int j = 1; j < nd; ++j) - { - transferMom[p][j - 1] -= tm[j]; - } - } - LOG(Message) << "Probe " << p << " (" << timeMom[p].size() << " points) : " << std::endl; - LOG(Message) << " phi(t_i, p_i) for (t_i, p_i) in " << timeMom[p] << std::endl; - LOG(Message) << " operator with momentum " << transferMom[p] << std::endl; - } - LOG(Message) << "FFT: field '" << par().field << "'" << std::endl; - fft.FFT_dim_mask(ftMatBuf, phi, dMask, FFT::forward); - slicedProbe.resize(timeMom.size()); - for (unsigned int p = 0; p < timeMom.size(); ++p) - { - std::vector qt; - - LOG(Message) << "Making probe " << p << std::endl; - slicedProbe[p].resize(nt); - for (unsigned int t = 0; t < nt; ++t) - { - Site acc; - - for (unsigned int i = 0; i < timeMom[p].size(); ++i) - { - Site buf; - - qt = timeMom[p][i]; - qt[0] += t; - vectorModulo(qt); - peekSite(buf, ftMatBuf, qt); - if (i == 0) - { - acc = buf; - } - else - { - acc *= buf; - } - } - slicedProbe[p][t] = TensorRemove(trace(acc)); - } - //std::cout << slicedProbe[p]<< std::endl; - } - for (auto &o: par().op) - { - auto &op = envGet(ComplexField, o); - - slicedOp[o].resize(transferMom.size()); - LOG(Message) << "FFT: operator '" << o << "'" << std::endl; - fft.FFT_dim_mask(ftBuf, op, dMask, FFT::forward); - //std::cout << ftBuf << std::endl; - for (unsigned int p = 0; p < transferMom.size(); ++p) - { - std::vector qt(nd, 0); - - for (unsigned int j = 1; j < nd; ++j) - { - qt[j] = transferMom[p][j - 1]; - } - slicedOp[o][p].resize(nt); - for (unsigned int t = 0; t < nt; ++t) - { - TComplex buf; - - qt[0] = t; - vectorModulo(qt); - peekSite(buf, ftBuf, qt); - slicedOp[o][p][t] = TensorRemove(buf); - } - //std::cout << ftBuf << std::endl; - //std::cout << slicedOp[o][p] << std::endl; - } - } - LOG(Message) << "Making correlators" << std::endl; - for (auto &o: par().op) - for (unsigned int p = 0; p < timeMom.size(); ++p) - { - TimeMomProbeResult r; - - LOG(Message) << " <" << o << " probe_" << p << ">" << std::endl; - r.op = o; - r.timeMom = timeMom[p]; - r.data = makeTwoPoint(slicedOp[o][p], slicedProbe[p], 1./partVol); - result.push_back(r); - } - saveResult(par().output, "timemomprobe", result); -} - -END_MODULE_NAMESPACE - -END_HADRONS_NAMESPACE - -#endif // Hadrons_MScalarSUN_TimeMomProbe_hpp_ diff --git a/Hadrons/Modules/MScalarSUN/TrMag.hpp b/Hadrons/Modules/MScalarSUN/TrMag.hpp index d7d40e88..b9602be3 100644 --- a/Hadrons/Modules/MScalarSUN/TrMag.hpp +++ b/Hadrons/Modules/MScalarSUN/TrMag.hpp @@ -124,7 +124,8 @@ void TTrMag::execute(void) std::vector result; auto &phi = envGet(Field, par().field); - auto m2 = sum(phi), mn = m2; + auto m2 = sum(phi); + auto mn = m2; m2 = -m2*m2; mn = 1.; diff --git a/Hadrons/Modules/MScalarSUN/Utils.hpp b/Hadrons/Modules/MScalarSUN/Utils.hpp index fdc42c9a..7eba5900 100644 --- a/Hadrons/Modules/MScalarSUN/Utils.hpp +++ b/Hadrons/Modules/MScalarSUN/Utils.hpp @@ -103,7 +103,7 @@ std::vector makeTwoPoint(const std::vector &sink, { for (unsigned int t = 0; t < nt; ++t) { - res[dt] += trace(sink[(t+dt)%nt]*source[t]); + res[dt] += trace(sink[(t+dt)%nt]*adj(source[t])); } res[dt] *= factor/static_cast(nt); } diff --git a/Hadrons/Modules/MSolver/A2AAslashVector.cc b/Hadrons/Modules/MSolver/A2AAslashVector.cc new file mode 100644 index 00000000..27344397 --- /dev/null +++ b/Hadrons/Modules/MSolver/A2AAslashVector.cc @@ -0,0 +1,35 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MSolver/A2AAslashVector.cc + +Copyright (C) 2015-2018 + +Author: Vera Guelpers + +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 + +using namespace Grid; +using namespace Hadrons; +using namespace MSolver; + +template class Grid::Hadrons::MSolver::TA2AAslashVector; +template class Grid::Hadrons::MSolver::TA2AAslashVector; diff --git a/Hadrons/Modules/MSolver/A2AAslashVector.hpp b/Hadrons/Modules/MSolver/A2AAslashVector.hpp new file mode 100644 index 00000000..3142d1e7 --- /dev/null +++ b/Hadrons/Modules/MSolver/A2AAslashVector.hpp @@ -0,0 +1,189 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MSolver/A2AAslashVector.hpp + +Copyright (C) 2015-2018 + +Author: Vera Guelpers + +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 */ +#ifndef Hadrons_MSolver_A2AAslashVector_hpp_ +#define Hadrons_MSolver_A2AAslashVector_hpp_ + +#include +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Create all-to-all V & W vectors * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MSolver) + +/**************************************************************************** +* Calculate a sequential propagator on an insertion of i*g_mu*A_mu +* on an A2A vector +* +* vv_i(y) = S(y,x) * i * g_mu*A_mu(x) * v_i(x) +* +* with +* +* - vector: A2A vector v_i(x) +* - emField: A_mu(x): electromagnetic photon field +* - solver: the solver for calculating the sequential propagator +* +*****************************************************************************/ + + +class A2AAslashVectorPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(A2AAslashVectorPar, + std::string, vector, + std::string, emField, + std::string, solver); +}; + +template +class TA2AAslashVector : public Module +{ +public: + FERM_TYPE_ALIASES(FImpl,); + SOLVER_TYPE_ALIASES(FImpl,); +public: + typedef PhotonR::GaugeField EmField; +public: + // constructor + TA2AAslashVector(const std::string name); + // destructor + virtual ~TA2AAslashVector(void) {}; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +private: + unsigned int Ls_; +}; + +MODULE_REGISTER_TMP(A2AAslashVector,TA2AAslashVector, MSolver); +MODULE_REGISTER_TMP(ZA2AAslashVector,TA2AAslashVector, MSolver); + +/****************************************************************************** + * TA2AAslashVector implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TA2AAslashVector::TA2AAslashVector(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TA2AAslashVector::getInput(void) +{ + std::vector in = {par().vector, par().emField, par().solver}; + + return in; +} + +template +std::vector TA2AAslashVector::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TA2AAslashVector::setup(void) +{ + Ls_ = env().getObjectLs(par().solver); + auto &vvector = envGet(std::vector, par().vector); + unsigned int Nmodes = vvector.size(); + envCreate(std::vector, getName(), 1, + Nmodes, envGetGrid(FermionField)); + + envTmpLat(FermionField, "v4dtmp"); + envTmpLat(FermionField, "v5dtmp", Ls_); + envTmpLat(FermionField, "v5dtmp_sol", Ls_); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TA2AAslashVector::execute(void) +{ + auto &solver = envGet(Solver, par().solver); + auto &stoch_photon = envGet(EmField, par().emField); + auto &vvector = envGet(std::vector, par().vector); + auto &Aslashv = envGet(std::vector, getName()); + unsigned int Nmodes = vvector.size(); + auto &mat = solver.getFMat(); + envGetTmp(FermionField, v4dtmp); + envGetTmp(FermionField, v5dtmp); + envGetTmp(FermionField, v5dtmp_sol); + + Complex ci(0.0,1.0); + + + startTimer("Seq Aslash"); + + LOG(Message) << "Calculate Sequential propagator on Aslash * v with the A2A vector " << par().vector + << " and the photon field " << par().emField << std::endl; + + + for(unsigned int i=0; i(stoch_photon, mu) * (gmu * vvector[i]); + } + stopTimer("Multiply Aslash"); + + if (Ls_ == 1) + { + solver(Aslashv[i], v4dtmp); + } + else + { + mat.ImportPhysicalFermionSource(v4dtmp, v5dtmp); + solver(v5dtmp_sol, v5dtmp); + mat.ExportPhysicalFermionSolution(v5dtmp_sol, v4dtmp); + Aslashv[i] = v4dtmp; + } + } + + stopTimer("Seq Aslash"); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MSolver_A2AAslashVector_hpp_ diff --git a/Hadrons/Modules/MSolver/A2AVectors.hpp b/Hadrons/Modules/MSolver/A2AVectors.hpp index ee7cfa30..f9980ee3 100644 --- a/Hadrons/Modules/MSolver/A2AVectors.hpp +++ b/Hadrons/Modules/MSolver/A2AVectors.hpp @@ -51,7 +51,9 @@ public: std::string, noise, std::string, action, std::string, eigenPack, - std::string, solver); + std::string, solver, + std::string, output, + bool, multiFile); }; template @@ -236,6 +238,17 @@ void TA2AVectors::execute(void) } stopTimer("W high mode"); } + + // I/O if necessary + if (!par().output.empty()) + { + startTimer("V I/O"); + A2AVectorsIo::write(par().output + "_v", v, par().multiFile, vm().getTrajectory()); + stopTimer("V I/O"); + startTimer("W I/O"); + A2AVectorsIo::write(par().output + "_w", w, par().multiFile, vm().getTrajectory()); + stopTimer("W I/O"); + } } END_MODULE_NAMESPACE diff --git a/Hadrons/Modules/MSolver/LocalCoherenceLanczos.cc b/Hadrons/Modules/MSolver/LocalCoherenceLanczos.cc index fa73e5b6..dacc871f 100644 --- a/Hadrons/Modules/MSolver/LocalCoherenceLanczos.cc +++ b/Hadrons/Modules/MSolver/LocalCoherenceLanczos.cc @@ -33,4 +33,7 @@ using namespace MSolver; template class Grid::Hadrons::MSolver::TLocalCoherenceLanczos; template class Grid::Hadrons::MSolver::TLocalCoherenceLanczos; - +#ifdef GRID_DEFAULT_PRECISION_DOUBLE +template class Grid::Hadrons::MSolver::TLocalCoherenceLanczos; +template class Grid::Hadrons::MSolver::TLocalCoherenceLanczos; +#endif diff --git a/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp b/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp index 7242eaf9..492ff39e 100644 --- a/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp +++ b/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp @@ -55,17 +55,17 @@ public: bool, multiFile); }; -template +template class TLocalCoherenceLanczos: public Module { public: FERM_TYPE_ALIASES(FImpl,); typedef LocalCoherenceLanczos LCL; - typedef BaseFermionEigenPack BasePack; - typedef CoarseFermionEigenPack CoarsePack; - typedef HADRONS_DEFAULT_SCHUR_OP SchurFMat; + nBasis> LCL; + typedef BaseFermionEigenPack BasePack; + typedef CoarseFermionEigenPack CoarsePack; + typedef HADRONS_DEFAULT_SCHUR_OP SchurFMat; public: // constructor TLocalCoherenceLanczos(const std::string name); @@ -82,27 +82,31 @@ public: MODULE_REGISTER_TMP(LocalCoherenceLanczos, ARG(TLocalCoherenceLanczos), MSolver); MODULE_REGISTER_TMP(ZLocalCoherenceLanczos, ARG(TLocalCoherenceLanczos), MSolver); +#ifdef GRID_DEFAULT_PRECISION_DOUBLE +MODULE_REGISTER_TMP(LocalCoherenceLanczosIo32, ARG(TLocalCoherenceLanczos), MSolver); +MODULE_REGISTER_TMP(ZLocalCoherenceLanczosIo32, ARG(TLocalCoherenceLanczos), MSolver); +#endif /****************************************************************************** * TLocalCoherenceLanczos implementation * ******************************************************************************/ // constructor ///////////////////////////////////////////////////////////////// -template -TLocalCoherenceLanczos::TLocalCoherenceLanczos(const std::string name) +template +TLocalCoherenceLanczos::TLocalCoherenceLanczos(const std::string name) : Module(name) {} // dependencies/products /////////////////////////////////////////////////////// -template -std::vector TLocalCoherenceLanczos::getInput(void) +template +std::vector TLocalCoherenceLanczos::getInput(void) { std::vector in = {par().action}; return in; } -template -std::vector TLocalCoherenceLanczos::getOutput(void) +template +std::vector TLocalCoherenceLanczos::getOutput(void) { std::vector out = {getName()}; @@ -110,8 +114,8 @@ std::vector TLocalCoherenceLanczos::getOutput(void) } // setup /////////////////////////////////////////////////////////////////////// -template -void TLocalCoherenceLanczos::setup(void) +template +void TLocalCoherenceLanczos::setup(void) { LOG(Message) << "Setting up local coherence Lanczos eigensolver for" << " action '" << par().action << "' (" << nBasis @@ -138,8 +142,8 @@ void TLocalCoherenceLanczos::setup(void) } // execution /////////////////////////////////////////////////////////////////// -template -void TLocalCoherenceLanczos::execute(void) +template +void TLocalCoherenceLanczos::execute(void) { auto &finePar = par().fineParams; auto &coarsePar = par().coarseParams; diff --git a/Hadrons/Modules/MSource/Momentum.cc b/Hadrons/Modules/MSource/Momentum.cc new file mode 100644 index 00000000..9bcf65ae --- /dev/null +++ b/Hadrons/Modules/MSource/Momentum.cc @@ -0,0 +1,36 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MSource/Momentum.cc + +Copyright (C) 2015-2018 + +Author: Antonin Portelli +Author: Peter Boyle + +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 + +using namespace Grid; +using namespace Hadrons; +using namespace MSource; + +template class Grid::Hadrons::MSource::TMomentum; + diff --git a/Hadrons/Modules/MSource/Momentum.hpp b/Hadrons/Modules/MSource/Momentum.hpp new file mode 100644 index 00000000..8336ea5c --- /dev/null +++ b/Hadrons/Modules/MSource/Momentum.hpp @@ -0,0 +1,149 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MSource/Momentum.hpp + +Copyright (C) 2015-2018 + +Author: Antonin Portelli +Author: Peter Boyle + +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 */ +#ifndef Hadrons_Momentum_hpp_ +#define Hadrons_Momentum_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/* +Plane Wave source +----------------- +src_x = e^i2pi/L * p *position +*/ + +/****************************************************************************** + * Plane Wave source * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MSource) + +class MomentumPar: Serializable +{ +public: +//What is meant by serializable in this context +GRID_SERIALIZABLE_CLASS_MEMBERS(MomentumPar, +std::string, mom); +}; + + +template +class TMomentum: public Module +{ +public: +FERM_TYPE_ALIASES(FImpl,); +public: +// constructor +TMomentum(const std::string name); +// destructor +virtual ~TMomentum(void) {}; +// dependency relation +virtual std::vector getInput(void); +virtual std::vector getOutput(void); +// setup +virtual void setup(void); +// execution +virtual void execute(void); +}; + +MODULE_REGISTER_TMP(Momentum, TMomentum, MSource); +//MODULE_REGISTER_NS(Momentum, TMomentum, MSource); + +/****************************************************************************** +* TMomentum template implementation * +******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TMomentum::TMomentum(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TMomentum::getInput(void) +{ + std::vector in; + return in; +} + +template +std::vector TMomentum::getOutput(void) +{ + std::vector out = {getName()}; + return out; +} + + +// setup /////////////////////////////////////////////////////////////////////// +template +void TMomentum::setup(void) +{ + envCreateLat(PropagatorField, getName()); +} + + +//execution////////////////////////////////////////////////////////////////// +template +void TMomentum::execute(void) +{ + LOG(Message) << "Generating planewave momentum source with momentum " << par().mom << std::endl; + //what does this env do? + PropagatorField &src = envGet(PropagatorField, getName()); + Lattice> t(env().getGrid()); + LatticeComplex C(env().getGrid()), coor(env().getGrid()); + std::vector p; + std::vector latt_size(GridDefaultLatt().begin(), GridDefaultLatt().end()); + Complex i(0.0,1.0); + + LOG(Message) << " " << std::endl; + //get the momentum from parameters + p = strToVec(par().mom); + C = zero; + LOG(Message) << "momentum converted from string - " << std::to_string(p[0]) < + +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 using namespace Grid; diff --git a/Hadrons/TimerArray.hpp b/Hadrons/TimerArray.hpp index 477f11cd..77cc2b8c 100644 --- a/Hadrons/TimerArray.hpp +++ b/Hadrons/TimerArray.hpp @@ -1,3 +1,30 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/TimerArray.hpp + +Copyright (C) 2015-2018 + +Author: Antonin Portelli + +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 */ #ifndef Hadrons_TimerArray_hpp_ #define Hadrons_TimerArray_hpp_ diff --git a/Hadrons/Utilities/EigenPackCast.cc b/Hadrons/Utilities/EigenPackCast.cc index 7247c92d..8c3a2b3f 100644 --- a/Hadrons/Utilities/EigenPackCast.cc +++ b/Hadrons/Utilities/EigenPackCast.cc @@ -1,3 +1,30 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Utilities/EigenPackCast.cc + +Copyright (C) 2015-2018 + +Author: Antonin Portelli + +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 #include @@ -8,7 +35,7 @@ using namespace Hadrons; template void convert(const std::string outFilename, const std::string inFilename, const unsigned int Ls, const bool rb, const unsigned int size, - const bool multiFile) + const bool multiFile, const bool testRead) { assert(outFilename != inFilename); @@ -75,6 +102,7 @@ void convert(const std::string outFilename, const std::string inFilename, LOG(Message) << "Out type : " << typeName() << std::endl; LOG(Message) << "#vectors : " << size << std::endl; LOG(Message) << "Multifile : " << (multiFile ? "yes" : "no") << std::endl; + LOG(Message) << "Test read : " << (testRead ? "yes" : "no") << std::endl; if (multiFile) { for(unsigned int k = 0; k < size; ++k) @@ -85,6 +113,8 @@ void convert(const std::string outFilename, const std::string inFilename, LOG(Message) << "==== Converting vector " << k << std::endl; LOG(Message) << "In : " << inV << std::endl; LOG(Message) << "Out: " << outV << std::endl; + // conversion + LOG(Message) << "-- Doing conversion" << std::endl; makeFileDir(outV, gOut); binWriter.open(outV); binReader.open(inV); @@ -94,10 +124,20 @@ void convert(const std::string outFilename, const std::string inFilename, EigenPackIo::writeElement(binWriter, bufIn, eval, k, &bufOut, &testIn); binWriter.close(); binReader.close(); + // read test + if (testRead) + { + LOG(Message) << "-- Test read" << std::endl; + binReader.open(outV); + EigenPackIo::readElement(bufOut, eval, k, binReader); + binReader.close(); + } } } else { + // conversion + LOG(Message) << "-- Doing conversion" << std::endl; makeFileDir(outFilename, gOut); binWriter.open(outFilename); binReader.open(inFilename); @@ -110,6 +150,18 @@ void convert(const std::string outFilename, const std::string inFilename, } binWriter.close(); binReader.close(); + // read test + if (testRead) + { + LOG(Message) << "-- Test read" << std::endl; + binReader.open(outFilename); + EigenPackIo::readHeader(record, binReader); + for(unsigned int k = 0; k < size; ++k) + { + EigenPackIo::readElement(bufOut, eval, k, binReader); + } + binReader.close(); + } } } @@ -127,11 +179,11 @@ int main(int argc, char *argv[]) // parse command line std::string outFilename, inFilename; unsigned int size, Ls; - bool rb, multiFile; + bool rb, multiFile, testRead; - if (argc < 7) + if (argc < 8) { - std::cerr << "usage: " << argv[0] << " <#vector> [Grid options]"; + std::cerr << "usage: " << argv[0] << " <#vector> [Grid options]"; std::cerr << std::endl; std::exit(EXIT_FAILURE); } @@ -141,6 +193,7 @@ int main(int argc, char *argv[]) rb = (std::string(argv[4]) != "0"); size = std::stoi(std::string(argv[5])); multiFile = (std::string(argv[6]) != "0"); + testRead = (std::string(argv[7]) != "0"); // initialization Grid_init(&argc, &argv); @@ -149,7 +202,7 @@ int main(int argc, char *argv[]) // execution try { - convert(outFilename, inFilename, Ls, rb, size, multiFile); + convert(outFilename, inFilename, Ls, rb, size, multiFile, testRead); } catch (const std::exception& e) { diff --git a/Hadrons/Utilities/HadronsXmlRun.cc b/Hadrons/Utilities/HadronsXmlRun.cc index c6ea7dc5..9b8fc901 100644 --- a/Hadrons/Utilities/HadronsXmlRun.cc +++ b/Hadrons/Utilities/HadronsXmlRun.cc @@ -2,7 +2,7 @@ Grid physics library, www.github.com/paboyle/Grid -Source file: Hadrons/HadronsXmlRun.cc +Source file: Hadrons/Utilities/HadronsXmlRun.cc Copyright (C) 2015-2018 diff --git a/Hadrons/modules.inc b/Hadrons/modules.inc index 866dfd91..2add645c 100644 --- a/Hadrons/modules.inc +++ b/Hadrons/modules.inc @@ -11,6 +11,7 @@ modules_cc =\ Modules/MContraction/Gamma3pt.cc \ Modules/MFermion/FreeProp.cc \ Modules/MFermion/GaugeProp.cc \ + Modules/MSource/Momentum.cc \ Modules/MSource/Point.cc \ Modules/MSource/Wall.cc \ Modules/MSource/SeqConserved.cc \ @@ -19,16 +20,20 @@ modules_cc =\ Modules/MSink/Point.cc \ Modules/MSink/Smear.cc \ Modules/MSolver/A2AVectors.cc \ + Modules/MSolver/A2AAslashVector.cc \ Modules/MSolver/RBPrecCG.cc \ Modules/MSolver/MixedPrecisionRBPrecCG.cc \ Modules/MSolver/LocalCoherenceLanczos.cc \ Modules/MGauge/StoutSmearing.cc \ Modules/MGauge/Unit.cc \ + Modules/MGauge/Electrify.cc \ Modules/MGauge/UnitEm.cc \ Modules/MGauge/StochEm.cc \ Modules/MGauge/Random.cc \ Modules/MGauge/FundtoHirep.cc \ + Modules/MGauge/GaugeFix.cc \ Modules/MNoise/TimeDilutedSpinColorDiagonal.cc \ + Modules/MNoise/FullVolumeSpinColorDiagonal.cc \ Modules/MUtilities/RandomVectors.cc \ Modules/MUtilities/TestSeqGamma.cc \ Modules/MUtilities/PrecisionCast.cc \ @@ -38,6 +43,9 @@ modules_cc =\ Modules/MScalar/VPCounterTerms.cc \ Modules/MScalar/ChargedProp.cc \ Modules/MScalar/ScalarVP.cc \ + Modules/MNPR/Amputate.cc \ + Modules/MNPR/Bilinear.cc \ + Modules/MNPR/FourQuark.cc \ Modules/MAction/Wilson.cc \ Modules/MAction/MobiusDWF.cc \ Modules/MAction/ZMobiusDWF.cc \ @@ -46,7 +54,6 @@ modules_cc =\ Modules/MAction/ScaledDWF.cc \ Modules/MScalarSUN/TrPhi.cc \ Modules/MScalarSUN/Grad.cc \ - Modules/MScalarSUN/TimeMomProbe.cc \ Modules/MScalarSUN/TrMag.cc \ Modules/MScalarSUN/TrKinetic.cc \ Modules/MScalarSUN/EMT.cc \ @@ -60,7 +67,8 @@ modules_cc =\ Modules/MIO/LoadBinary.cc \ Modules/MIO/LoadNersc.cc \ Modules/MIO/LoadCoarseEigenPack.cc \ - Modules/MIO/LoadCosmHol.cc + Modules/MIO/LoadCosmHol.cc \ + Modules/MIO/LoadA2AVectors.cc modules_hpp =\ Modules/MContraction/Baryon.hpp \ @@ -81,6 +89,7 @@ modules_hpp =\ Modules/MSource/Wall.hpp \ Modules/MSource/Z2.hpp \ Modules/MSource/SeqConserved.hpp \ + Modules/MSource/Momentum.hpp \ Modules/MSink/Smear.hpp \ Modules/MSink/Point.hpp \ Modules/MSolver/MixedPrecisionRBPrecCG.hpp \ @@ -88,13 +97,17 @@ modules_hpp =\ Modules/MSolver/Guesser.hpp \ Modules/MSolver/RBPrecCG.hpp \ Modules/MSolver/A2AVectors.hpp \ + Modules/MSolver/A2AAslashVector.hpp \ Modules/MGauge/UnitEm.hpp \ Modules/MGauge/StoutSmearing.hpp \ Modules/MGauge/Unit.hpp \ + Modules/MGauge/Electrify.hpp \ Modules/MGauge/Random.hpp \ + Modules/MGauge/GaugeFix.hpp \ Modules/MGauge/FundtoHirep.hpp \ Modules/MGauge/StochEm.hpp \ Modules/MNoise/TimeDilutedSpinColorDiagonal.hpp \ + Modules/MNoise/FullVolumeSpinColorDiagonal.hpp \ Modules/MUtilities/PrecisionCast.hpp \ Modules/MUtilities/RandomVectors.hpp \ Modules/MUtilities/TestSeqGamma.hpp \ @@ -105,6 +118,9 @@ modules_hpp =\ Modules/MScalar/ScalarVP.hpp \ Modules/MScalar/Scalar.hpp \ Modules/MScalar/ChargedProp.hpp \ + Modules/MNPR/Bilinear.hpp \ + Modules/MNPR/Amputate.hpp \ + Modules/MNPR/FourQuark.hpp \ Modules/MAction/DWF.hpp \ Modules/MAction/MobiusDWF.hpp \ Modules/MAction/Wilson.hpp \ @@ -115,7 +131,6 @@ modules_hpp =\ Modules/MScalarSUN/TwoPointNPR.hpp \ Modules/MScalarSUN/ShiftProbe.hpp \ Modules/MScalarSUN/Div.hpp \ - Modules/MScalarSUN/TimeMomProbe.hpp \ Modules/MScalarSUN/TrMag.hpp \ Modules/MScalarSUN/EMT.hpp \ Modules/MScalarSUN/TwoPoint.hpp \ @@ -126,6 +141,7 @@ modules_hpp =\ Modules/MScalarSUN/TrKinetic.hpp \ Modules/MIO/LoadEigenPack.hpp \ Modules/MIO/LoadNersc.hpp \ + Modules/MIO/LoadA2AVectors.hpp \ Modules/MIO/LoadCosmHol.hpp \ Modules/MIO/LoadCoarseEigenPack.hpp \ Modules/MIO/LoadBinary.hpp diff --git a/benchmarks/Benchmark_IO.cc b/benchmarks/Benchmark_IO.cc index 479ae037..58e0340d 100644 --- a/benchmarks/Benchmark_IO.cc +++ b/benchmarks/Benchmark_IO.cc @@ -1,108 +1,48 @@ -#include -#ifdef HAVE_LIME -using namespace std; -using namespace Grid; -using namespace Grid::QCD; +#include "Benchmark_IO.hpp" -#define MSG cout << GridLogMessage -#define SEP \ -"=============================================================================" #ifndef BENCH_IO_LMAX #define BENCH_IO_LMAX 40 #endif -typedef function WriterFn; -typedef function ReaderFn; +using namespace Grid; +using namespace QCD; -string filestem(const int l) +std::string filestem(const int l) { - return "iobench_l" + to_string(l); -} - -void limeWrite(const string filestem, LatticeFermion &vec) -{ - emptyUserRecord record; - ScidacWriter binWriter(vec._grid->IsBoss()); - - binWriter.open(filestem + ".bin"); - binWriter.writeScidacFieldRecord(vec, record); - binWriter.close(); -} - -void limeRead(LatticeFermion &vec, const string filestem) -{ - emptyUserRecord record; - ScidacReader binReader; - - binReader.open(filestem + ".bin"); - binReader.readScidacFieldRecord(vec, record); - binReader.close(); -} - -void writeBenchmark(const int l, const WriterFn &write) -{ - auto mpi = GridDefaultMpi(); - auto simd = GridDefaultSimd(Nd, vComplex::Nsimd()); - vector latt = {l*mpi[0], l*mpi[1], l*mpi[2], l*mpi[3]}; - unique_ptr gPt(SpaceTimeGrid::makeFourDimGrid(latt, simd, mpi)); - GridCartesian *g = gPt.get(); - GridParallelRNG rng(g); - LatticeFermion vec(g); - emptyUserRecord record; - ScidacWriter binWriter(g->IsBoss()); - - cout << "-- Local volume " << l << "^4" << endl; - random(rng, vec); - write(filestem(l), vec); -} - -void readBenchmark(const int l, const ReaderFn &read) -{ - auto mpi = GridDefaultMpi(); - auto simd = GridDefaultSimd(Nd, vComplex::Nsimd()); - vector latt = {l*mpi[0], l*mpi[1], l*mpi[2], l*mpi[3]}; - unique_ptr gPt(SpaceTimeGrid::makeFourDimGrid(latt, simd, mpi)); - GridCartesian *g = gPt.get(); - LatticeFermion vec(g); - emptyUserRecord record; - ScidacReader binReader; - - cout << "-- Local volume " << l << "^4" << endl; - read(vec, filestem(l)); + return "iobench_l" + std::to_string(l); } int main (int argc, char ** argv) { Grid_init(&argc,&argv); - auto simd = GridDefaultSimd(Nd,vComplex::Nsimd()); - auto mpi = GridDefaultMpi(); - int64_t threads = GridThread::GetThreads(); - MSG << "Grid is setup to use " << threads << " threads" << endl; - MSG << SEP << endl; - MSG << "Benchmark Lime write" << endl; - MSG << SEP << endl; + MSG << "Grid is setup to use " << threads << " threads" << std::endl; + MSG << SEP << std::endl; + MSG << "Benchmark Lime write" << std::endl; + MSG << SEP << std::endl; for (int l = 4; l <= BENCH_IO_LMAX; l += 2) { - writeBenchmark(l, limeWrite); + auto mpi = GridDefaultMpi(); + std::vector latt = {l*mpi[0], l*mpi[1], l*mpi[2], l*mpi[3]}; + + std::cout << "-- Local volume " << l << "^4" << std::endl; + writeBenchmark(latt, filestem(l), limeWrite); } - MSG << "Benchmark Lime read" << endl; - MSG << SEP << endl; + MSG << "Benchmark Lime read" << std::endl; + MSG << SEP << std::endl; for (int l = 4; l <= BENCH_IO_LMAX; l += 2) { - readBenchmark(l, limeRead); + auto mpi = GridDefaultMpi(); + std::vector latt = {l*mpi[0], l*mpi[1], l*mpi[2], l*mpi[3]}; + + std::cout << "-- Local volume " << l << "^4" << std::endl; + readBenchmark(latt, filestem(l), limeRead); } Grid_finalize(); return EXIT_SUCCESS; } -#else -int main (int argc, char ** argv) -{ - return EXIT_SUCCESS; -} -#endif diff --git a/benchmarks/Benchmark_IO.hpp b/benchmarks/Benchmark_IO.hpp new file mode 100644 index 00000000..9565c928 --- /dev/null +++ b/benchmarks/Benchmark_IO.hpp @@ -0,0 +1,107 @@ +#ifndef Benchmark_IO_hpp_ +#define Benchmark_IO_hpp_ + +#include + +#define MSG std::cout << GridLogMessage +#define SEP \ +"=============================================================================" + +namespace Grid { + +template +using WriterFn = std::function ; +template +using ReaderFn = std::function; + +template +void limeWrite(const std::string filestem, Field &vec) +{ + emptyUserRecord record; + QCD::ScidacWriter binWriter(vec._grid->IsBoss()); + + binWriter.open(filestem + ".bin"); + binWriter.writeScidacFieldRecord(vec, record); + binWriter.close(); +} + +template +void limeRead(Field &vec, const std::string filestem) +{ + emptyUserRecord record; + QCD::ScidacReader binReader; + + binReader.open(filestem + ".bin"); + binReader.readScidacFieldRecord(vec, record); + binReader.close(); +} + +inline void makeGrid(std::shared_ptr &gPt, + const std::shared_ptr &gBasePt, + const unsigned int Ls = 1, const bool rb = false) +{ + if (rb) + { + if (Ls > 1) + { + gPt.reset(QCD::SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, gBasePt.get())); + } + else + { + gPt.reset(QCD::SpaceTimeGrid::makeFourDimRedBlackGrid(gBasePt.get())); + } + } + else + { + if (Ls > 1) + { + gPt.reset(QCD::SpaceTimeGrid::makeFiveDimGrid(Ls, gBasePt.get())); + } + else + { + gPt = gBasePt; + } + } +} + +template +void writeBenchmark(const std::vector &latt, const std::string filename, + const WriterFn &write, + const unsigned int Ls = 1, const bool rb = false) +{ + auto mpi = GridDefaultMpi(); + auto simd = GridDefaultSimd(latt.size(), Field::vector_type::Nsimd()); + std::shared_ptr gBasePt(QCD::SpaceTimeGrid::makeFourDimGrid(latt, simd, mpi)); + std::shared_ptr gPt; + + makeGrid(gPt, gBasePt, Ls, rb); + + GridBase *g = gPt.get(); + GridParallelRNG rng(g); + Field vec(g); + + random(rng, vec); + write(filename, vec); +} + +template +void readBenchmark(const std::vector &latt, const std::string filename, + const ReaderFn &read, + const unsigned int Ls = 1, const bool rb = false) +{ + auto mpi = GridDefaultMpi(); + auto simd = GridDefaultSimd(latt.size(), Field::vector_type::Nsimd()); + std::shared_ptr gBasePt(QCD::SpaceTimeGrid::makeFourDimGrid(latt, simd, mpi)); + std::shared_ptr gPt; + + makeGrid(gPt, gBasePt, Ls, rb); + + GridBase *g = gPt.get(); + Field vec(g); + + read(vec, filename); +} + +} + +#endif // Benchmark_IO_hpp_ diff --git a/benchmarks/Benchmark_IO_vs_dir.cc b/benchmarks/Benchmark_IO_vs_dir.cc new file mode 100644 index 00000000..15f6194f --- /dev/null +++ b/benchmarks/Benchmark_IO_vs_dir.cc @@ -0,0 +1,79 @@ +#include "Benchmark_IO.hpp" + +#define MSG std::cout << GridLogMessage +#define SEP \ +"=============================================================================" + +using namespace Grid; +using namespace QCD; + +int main (int argc, char ** argv) +{ + std::vector dir; + unsigned int Ls; + bool rb; + if (argc < 4) + { + std::cerr << "usage: " << argv[0] << " [ ... ] [Grid options]"; + std::cerr << std::endl; + } + Ls = std::stoi(argv[1]); + rb = (std::string(argv[2]) == "1"); + for (unsigned int i = 3; i < argc; ++i) + { + std::string a = argv[i]; + + if (a[0] != '-') + { + dir.push_back(std::string(argv[i])); + } + else + { + break; + } + } + Grid_init(&argc,&argv); + + + int64_t threads = GridThread::GetThreads(); + MSG << "Grid is setup to use " << threads << " threads" << std::endl; + MSG << SEP << std::endl; + MSG << "Benchmark double precision Lime write" << std::endl; + MSG << SEP << std::endl; + for (auto &d: dir) + { + MSG << "-- Directory " << d << std::endl; + writeBenchmark(GridDefaultLatt(), d + "/ioBench", limeWrite, Ls, rb); + } + + MSG << SEP << std::endl; + MSG << "Benchmark double precision Lime read" << std::endl; + MSG << SEP << std::endl; + for (auto &d: dir) + { + MSG << "-- Directory " << d << std::endl; + readBenchmark(GridDefaultLatt(), d + "/ioBench", limeRead, Ls, rb); + } + + MSG << SEP << std::endl; + MSG << "Benchmark single precision Lime write" << std::endl; + MSG << SEP << std::endl; + for (auto &d: dir) + { + MSG << "-- Directory " << d << std::endl; + writeBenchmark(GridDefaultLatt(), d + "/ioBench", limeWrite, Ls, rb); + } + + MSG << SEP << std::endl; + MSG << "Benchmark single precision Lime read" << std::endl; + MSG << SEP << std::endl; + for (auto &d: dir) + { + MSG << "-- Directory " << d << std::endl; + readBenchmark(GridDefaultLatt(), d + "/ioBench", limeRead, Ls, rb); + } + + Grid_finalize(); + + return EXIT_SUCCESS; +} diff --git a/configure.ac b/configure.ac index d362ed5a..ad2b8bba 100644 --- a/configure.ac +++ b/configure.ac @@ -485,6 +485,7 @@ DX_INIT_DOXYGEN([$PACKAGE_NAME], [doxygen.cfg]) ############### Ouput cwd=`pwd -P`; cd ${srcdir}; abs_srcdir=`pwd -P`; cd ${cwd} +GRID_CXX="$CXX" GRID_CXXFLAGS="$AM_CXXFLAGS $CXXFLAGS" GRID_LDFLAGS="$AM_LDFLAGS $LDFLAGS" GRID_LIBS=$LIBS @@ -497,6 +498,7 @@ AM_LDFLAGS="-L${cwd}/Grid $AM_LDFLAGS" AC_SUBST([AM_CFLAGS]) AC_SUBST([AM_CXXFLAGS]) AC_SUBST([AM_LDFLAGS]) +AC_SUBST([GRID_CXX]) AC_SUBST([GRID_CXXFLAGS]) AC_SUBST([GRID_LDFLAGS]) AC_SUBST([GRID_LIBS]) diff --git a/grid-config.in b/grid-config.in index bd340846..f39b01bd 100755 --- a/grid-config.in +++ b/grid-config.in @@ -61,6 +61,10 @@ while test $# -gt 0; do echo @GRID_CXXFLAGS@ ;; + --cxx) + echo @GRID_CXX@ + ;; + --ldflags) echo @GRID_LDFLAGS@ ;; diff --git a/tests/debug/Test_cayley_cg.cc b/tests/debug/Test_cayley_cg.cc index 57fe7717..5150268f 100644 --- a/tests/debug/Test_cayley_cg.cc +++ b/tests/debug/Test_cayley_cg.cc @@ -1,5 +1,4 @@ /************************************************************************************* - Grid physics library, www.github.com/paboyle/Grid Source file: ./tests/Test_cayley_cg.cc @@ -27,6 +26,7 @@ Author: paboyle *************************************************************************************/ /* END LEGAL */ #include +#include using namespace std; using namespace Grid; @@ -46,6 +46,7 @@ struct scal { template void TestCGinversions(What & Ddwf, + LatticeGaugeField &Umu, GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, RealD mass, RealD M5, @@ -75,6 +76,25 @@ void TestCGprec(What & Ddwf, GridParallelRNG *RNG4, GridParallelRNG *RNG5); +template +void TestReconstruct5D(What & Ddwf, + LatticeGaugeField &Umu, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5); + +template +void TestReconstruct5DFA(What & Ddwf, + WhatF & DdwfF, + LatticeGaugeField &Umu, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5); + int main (int argc, char ** argv) { Grid_init(&argc,&argv); @@ -83,63 +103,104 @@ int main (int argc, char ** argv) std::cout< seeds4({1,2,3,4}); std::vector seeds5({5,6,7,8}); GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); LatticeGaugeField Umu(UGrid); + LatticeGaugeFieldF UmuF(UGridF); SU3::HotConfiguration(RNG4,Umu); + precisionChange(UmuF,Umu); std::vector U(4,UGrid); RealD mass=0.1; RealD M5 =1.8; + std::cout<(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + DomainWallFermionF DdwfF(UmuF,*FGridF,*FrbGridF,*UGridF,*UrbGridF,mass,M5); + TestCGinversions(Ddwf,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestReconstruct5DFA(Ddwf,DdwfF,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); RealD b=1.5;// Scale factor b+c=2, b-c=1 RealD c=0.5; std::vector gamma(Ls,ComplexD(1.0,0.0)); + std::cout<(Dmob,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + MobiusFermionF DmobF(UmuF,*FGridF,*FrbGridF,*UGridF,*UrbGridF,mass,M5,b,c); + TestCGinversions(Dmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestReconstruct5DFA(Dmob,DmobF,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + std::cout<(ZDmob,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestCGinversions(ZDmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestReconstruct5D(ZDmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + std::cout<(Dzolo,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestCGinversions(Dzolo,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestReconstruct5D(Dzolo,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + std::cout<(Dsham,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + ScaledShamirFermionF DshamF(UmuF,*FGridF,*FrbGridF,*UGridF,*UrbGridF,mass,M5,2.0); + TestCGinversions(Dsham,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestReconstruct5DFA(Dsham,DshamF,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + std::cout<(Dshamz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestCGinversions(Dshamz,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestReconstruct5D(Dshamz,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + std::cout<(Dov,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + std::cout<(Dov,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestReconstruct5DFA(Dov,DovF,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + std::cout<(Dovz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestCGinversions(Dovz,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestReconstruct5D(Dovz,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); Grid_finalize(); } template void TestCGinversions(What & Ddwf, + LatticeGaugeField &Umu, GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, RealD mass, RealD M5, @@ -154,6 +215,7 @@ void TestCGinversions(What & Ddwf, TestCGschur(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5); } + template void TestCGunprec(What & Ddwf, GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, @@ -189,6 +251,147 @@ void TestCGprec(What & Ddwf, CG(HermOpEO,src_o,result_o); } +template +void TestReconstruct5D(What & Ddwf, + LatticeGaugeField & Umu, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5) +{ + LatticeFermion src4 (UGrid); random(*RNG4,src4); + LatticeFermion res4 (UGrid); res4 = zero; + + LatticeFermion src (FGrid); + LatticeFermion src_NE(FGrid); + LatticeFermion result(FGrid); + LatticeFermion result_rec(FGrid); + LatticeFermion result_madwf(FGrid); + + MdagMLinearOperator HermOp(Ddwf); + double Resid = 1.0e-12; + double Residi = 1.0e-6; + ConjugateGradient CG(Resid,10000); + ConjugateGradient CGi(Residi,10000); + + Ddwf.ImportPhysicalFermionSource(src4,src); + Ddwf.Mdag(src,src_NE); + CG(HermOp,src_NE,result); + + Ddwf.ExportPhysicalFermionSolution(result, res4); + + Ddwf.M(result,src_NE); + src_NE = src_NE - src; + std::cout < SchurSolverType; + typedef SchurRedBlackDiagTwoSolve SchurSolverTypei; + typedef PauliVillarsSolverRBprec PVinverter; + SchurSolverType SchurSolver(CG); + PVinverter PVinverse(SchurSolver); + + Reconstruct5DfromPhysical reconstructor(PVinverse); + + reconstructor(Ddwf,res4,src4,result_rec); + + std::cout < Guess; + MADWF > + madwf(Ddwf,Ddwf,PVinverse,SchurSolveri,Guess,Resid,10); + + madwf(src4,result_madwf); + result_madwf = result_madwf - result; + std::cout < +void TestReconstruct5DFA(What & Ddwf, + WhatF & DdwfF, + LatticeGaugeField & Umu, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5) +{ + LatticeFermion src4 (UGrid); random(*RNG4,src4); + LatticeFermion res4 (UGrid); res4 = zero; + + LatticeFermion src (FGrid); + LatticeFermion src_NE(FGrid); + LatticeFermion result(FGrid); + LatticeFermion result_rec(FGrid); + LatticeFermion result_madwf(FGrid); + + MdagMLinearOperator HermOp(Ddwf); + double Resid = 1.0e-12; + double Residi = 1.0e-5; + ConjugateGradient CG(Resid,10000); + ConjugateGradient CGi(Residi,10000); + + Ddwf.ImportPhysicalFermionSource(src4,src); + Ddwf.Mdag(src,src_NE); + CG(HermOp,src_NE,result); + + Ddwf.ExportPhysicalFermionSolution(result, res4); + + Ddwf.M(result,src_NE); + src_NE = src_NE - src; + std::cout < SchurSolverTypei; + typedef PauliVillarsSolverFourierAccel PVinverter; + PVinverter PVinverse(Umu,CG); + + Reconstruct5DfromPhysical reconstructor(PVinverse); + + reconstructor(Ddwf,res4,src4,result_rec); + + std::cout < Guess; + MADWF > + madwf(Ddwf,DdwfF,PVinverse,SchurSolver,Guess,Resid,10); + + madwf(src4,result_madwf); + result_madwf = result_madwf - result; + std::cout < void TestCGschur(What & Ddwf, diff --git a/tests/debug/Test_split_laplacian.cc b/tests/debug/Test_split_laplacian.cc new file mode 100644 index 00000000..174dc1b1 --- /dev/null +++ b/tests/debug/Test_split_laplacian.cc @@ -0,0 +1,104 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_dwf_mrhs_cg.cc + + Copyright (C) 2015 + +Author: Peter Boyle + + 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 + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + + +int main (int argc, char ** argv) +{ + typedef LatticeComplex ComplexField; + + Grid_init(&argc,&argv); + + std::vector latt_size = GridDefaultLatt(); + int nd = latt_size.size(); + int ndm1 = nd-1; + + std::vector simd_layout = GridDefaultSimd(nd,vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + std::vector mpi_split (mpi_layout.size(),1); + + std::cout << " Full " << GridCmdVectorIntToString(latt_size) << " subgrid" < latt_m = latt_size; latt_m[nd-1] = 1; + std::vector mpi_m = mpi_layout; mpi_m [nd-1] = 1; + std::vector simd_m = GridDefaultSimd(ndm1,vComplex::Nsimd()); simd_m.push_back(1); + + + std::cout << " Requesting " << GridCmdVectorIntToString(latt_m)<< " subgrid" <GlobalSum(tmp); + std::cout << GridLogMessage<< " Full nodes "<< tmp <GlobalSum(tmp); + std::cout << GridLogMessage<< " Split nodes "<< tmp <Barrier(); + + auto local_latt = GridN->LocalDimensions(); + + Full_cpy = zero; + std::vector seeds({1,2,3,4}); + GridParallelRNG RNG(GridN); RNG.SeedFixedIntegers(seeds); + + random(RNG,Full); + for(int t=0;t("DWF_" + flavour[i], actionPar); diff --git a/tests/hadrons/Test_diskvector.cc b/tests/hadrons/Test_diskvector.cc index 363ae2ce..10bc4db1 100644 --- a/tests/hadrons/Test_diskvector.cc +++ b/tests/hadrons/Test_diskvector.cc @@ -91,6 +91,22 @@ int main(int argc, char *argv[]) v13r = v[13]; LOG(Message) << "v[13] correct? " << ((v13r == v13w) ? "yes" : "no" ) << std::endl; + LOG(Message) << "hit ratio " << v.hitRatio() << std::endl; + + EigenDiskVector w("eigendiskvector_test", 1000, 4); + EigenDiskVector::Matrix m,n; + + w[2] = EigenDiskVectorMat::Random(2000, 2000); + m = w[2]; + w[3] = EigenDiskVectorMat::Random(2000, 2000); + w[4] = EigenDiskVectorMat::Random(2000, 2000); + w[5] = EigenDiskVectorMat::Random(2000, 2000); + w[6] = EigenDiskVectorMat::Random(2000, 2000); + w[7] = EigenDiskVectorMat::Random(2000, 2000); + n = w[2]; + LOG(Message) << "w[2] correct? " + << ((m == n) ? "yes" : "no" ) << std::endl; + LOG(Message) << "hit ratio " << w.hitRatio() << std::endl; Grid_finalize(); diff --git a/tests/hadrons/Test_hadrons.hpp b/tests/hadrons/Test_hadrons.hpp index 1ee72356..7f07bea5 100644 --- a/tests/hadrons/Test_hadrons.hpp +++ b/tests/hadrons/Test_hadrons.hpp @@ -126,6 +126,7 @@ inline void makeWilsonAction(Application &application, std::string actionName, actionPar.gauge = gaugeField; actionPar.mass = mass; actionPar.boundary = boundary; + actionPar.twist = "0. 0. 0. 0."; application.createModule(actionName, actionPar); } } @@ -154,6 +155,7 @@ inline void makeDWFAction(Application &application, std::string actionName, actionPar.M5 = M5; actionPar.mass = mass; actionPar.boundary = boundary; + actionPar.twist = "0. 0. 0. 0."; application.createModule(actionName, actionPar); } } diff --git a/tests/hadrons/Test_hadrons_meson_3pt.cc b/tests/hadrons/Test_hadrons_meson_3pt.cc index 5b7cda22..741e2a7c 100644 --- a/tests/hadrons/Test_hadrons_meson_3pt.cc +++ b/tests/hadrons/Test_hadrons_meson_3pt.cc @@ -66,6 +66,7 @@ int main(int argc, char *argv[]) // set fermion boundary conditions to be periodic space, antiperiodic time. std::string boundary = "1 1 1 -1"; + std::string twist = "0. 0. 0. 0."; // sink MSink::Point::Par sinkPar; @@ -80,6 +81,7 @@ int main(int argc, char *argv[]) actionPar.M5 = 1.8; actionPar.mass = mass[i]; actionPar.boundary = boundary; + actionPar.twist = twist; application.createModule("DWF_" + flavour[i], actionPar); // solvers diff --git a/tests/hadrons/Test_hadrons_spectrum.cc b/tests/hadrons/Test_hadrons_spectrum.cc index b57deb38..af1dccd7 100644 --- a/tests/hadrons/Test_hadrons_spectrum.cc +++ b/tests/hadrons/Test_hadrons_spectrum.cc @@ -72,6 +72,7 @@ int main(int argc, char *argv[]) // set fermion boundary conditions to be periodic space, antiperiodic time. std::string boundary = "1 1 1 -1"; + std::string twist = "0. 0. 0. 0."; for (unsigned int i = 0; i < flavour.size(); ++i) { @@ -82,6 +83,7 @@ int main(int argc, char *argv[]) actionPar.M5 = 1.8; actionPar.mass = mass[i]; actionPar.boundary = boundary; + actionPar.twist = twist; application.createModule("DWF_" + flavour[i], actionPar); // solvers diff --git a/tests/solver/Test_dwf_mrhs_cg_mpi.cc b/tests/solver/Test_dwf_mrhs_cg_mpi.cc index aa36ebbc..c1d8ce62 100644 --- a/tests/solver/Test_dwf_mrhs_cg_mpi.cc +++ b/tests/solver/Test_dwf_mrhs_cg_mpi.cc @@ -38,6 +38,7 @@ int main (int argc, char ** argv) typedef typename DomainWallFermionR::ComplexField ComplexField; typename DomainWallFermionR::ImplParams params; + double stp=1.0e-5; const int Ls=4; Grid_init(&argc,&argv); @@ -197,7 +198,7 @@ int main (int argc, char ** argv) MdagMLinearOperator HermOp(Ddwf); MdagMLinearOperator HermOpCk(Dchk); - ConjugateGradient CG((1.0e-2),10000); + ConjugateGradient CG((stp),10000); s_res = zero; CG(HermOp,s_src,s_res); @@ -227,5 +228,11 @@ int main (int argc, char ** argv) std::cout << GridLogMessage<<" resid["< BCGV (BlockCGVec,blockDim,stp,10000); + BCGV.PrintInterval=10; + BCGV(HermOpCk,src,result); + Grid_finalize(); } diff --git a/tests/solver/Test_mobius_bcg.cc b/tests/solver/Test_mobius_bcg.cc new file mode 100644 index 00000000..e59cb7e0 --- /dev/null +++ b/tests/solver/Test_mobius_bcg.cc @@ -0,0 +1,220 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_dwf_mrhs_cg.cc + + Copyright (C) 2015 + +Author: Peter Boyle + + 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 +#include + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +int main (int argc, char ** argv) +{ + typedef typename MobiusFermionR::FermionField FermionField; + typedef typename MobiusFermionR::ComplexField ComplexField; + typename MobiusFermionR::ImplParams params; + + const int Ls=12; + + Grid_init(&argc,&argv); + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + std::vector mpi_split (mpi_layout.size(),1); + std::vector split_coor (mpi_layout.size(),1); + std::vector split_dim (mpi_layout.size(),1); + + std::vector boundary_phases(Nd,1.); + boundary_phases[Nd-1]=-1.; + params.boundary_phases = boundary_phases; + + GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), + GridDefaultSimd(Nd,vComplex::Nsimd()), + GridDefaultMpi()); + GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); + GridRedBlackCartesian * rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); + + ///////////////////////////////////////////// + // Split into 1^4 mpi communicators + ///////////////////////////////////////////// + + for(int i=0;i> mpi_split[k]; + } + break; + } + } + + + double stp = 1.e-8; + int nrhs = 1; + int me; + for(int i=0;i seeds({1,2,3,4}); + + std::vector src(nrhs,FGrid); + std::vector src_chk(nrhs,FGrid); + std::vector result(nrhs,FGrid); + FermionField tmp(FGrid); + std::cout << GridLogMessage << "Made the Fermion Fields"< HermOp(Ddwf); + MdagMLinearOperator HermOpCk(Dchk); + ConjugateGradient CG((stp),100000); + s_res = zero; + + CG(HermOp,s_src,s_res); + + std::cout << GridLogMessage << " split residual norm "< iterations(nrhs,0); + iterations[me] = CG.IterationsToComplete; + + for(int n=0;nGlobalSum(iterations[n]); + std::cout << GridLogMessage<<" Rank "< BCGV (BlockCGrQVec,blockDim,stp,100000); + { + BCGV(HermOpCk,src,result); + } + + + + Grid_finalize(); +} diff --git a/tests/solver/Test_mobius_bcg_nosplit.cc b/tests/solver/Test_mobius_bcg_nosplit.cc new file mode 100644 index 00000000..f3ed621f --- /dev/null +++ b/tests/solver/Test_mobius_bcg_nosplit.cc @@ -0,0 +1,144 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_dwf_mrhs_cg.cc + + Copyright (C) 2015 + +Author: Peter Boyle + + 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 + +#include +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +int main (int argc, char ** argv) +{ + typedef typename DomainWallFermionR::FermionField FermionField; + typedef typename DomainWallFermionR::ComplexField ComplexField; + typename DomainWallFermionR::ImplParams params; + + const int Ls=16; + + Grid_init(&argc,&argv); + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + + std::vector boundary_phases(Nd,1.); + boundary_phases[Nd-1]=-1.; + params.boundary_phases = boundary_phases; + + GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), + GridDefaultSimd(Nd,vComplex::Nsimd()), + GridDefaultMpi()); + GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); + GridRedBlackCartesian * rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); + + double stp = 1.e-8; + int nrhs = 2; + + /////////////////////////////////////////////// + // Set up the problem as a 4d spreadout job + /////////////////////////////////////////////// + std::vector seeds({1,2,3,4}); + + std::vector src(nrhs,FGrid); + std::vector src_chk(nrhs,FGrid); + std::vector result(nrhs,FGrid); + FermionField tmp(FGrid); + std::cout << GridLogMessage << "Made the Fermion Fields"< HermOp(Ddwf); + ConjugateGradient CG((stp),100000); + + for(int rhs=0;rhs<1;rhs++){ + result[rhs] = zero; + CG(HermOp,src[rhs],result[rhs]); + } + + for(int rhs=0;rhs<1;rhs++){ + std::cout << " Result["< + +#include +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +int main (int argc, char ** argv) +{ + typedef typename DomainWallFermionR::FermionField FermionField; + typedef typename DomainWallFermionR::ComplexField ComplexField; + typename DomainWallFermionR::ImplParams params; + + const int Ls=16; + + Grid_init(&argc,&argv); + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + + std::vector boundary_phases(Nd,1.); + boundary_phases[Nd-1]=-1.; + params.boundary_phases = boundary_phases; + + GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), + GridDefaultSimd(Nd,vComplex::Nsimd()), + GridDefaultMpi()); + GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); + GridRedBlackCartesian * rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); + + double stp = 1.e-8; + int nrhs = 2; + + /////////////////////////////////////////////// + // Set up the problem as a 4d spreadout job + /////////////////////////////////////////////// + std::vector seeds({1,2,3,4}); + + std::vector src4(nrhs,UGrid); + std::vector src(nrhs,FGrid); + std::vector src_chk(nrhs,FGrid); + std::vector result(nrhs,FGrid); + FermionField tmp(FGrid); + std::cout << GridLogMessage << "Made the Fermion Fields"< HermOp(Ddwf); + ConjugateGradient CG((stp),100000); + + for(int rhs=0;rhs<1;rhs++){ + result[rhs] = zero; + // CG(HermOp,src[rhs],result[rhs]); + } + + for(int rhs=0;rhs<1;rhs++){ + std::cout << " Result["< + +#include +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +int main (int argc, char ** argv) +{ + typedef typename DomainWallFermionR::FermionField FermionField; + typedef typename DomainWallFermionR::ComplexField ComplexField; + typename DomainWallFermionR::ImplParams params; + + const int Ls=16; + + Grid_init(&argc,&argv); + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + + std::vector boundary_phases(Nd,1.); + boundary_phases[Nd-1]=-1.; + params.boundary_phases = boundary_phases; + + GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), + GridDefaultSimd(Nd,vComplex::Nsimd()), + GridDefaultMpi()); + GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); + GridRedBlackCartesian * rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); + + double stp = 1.e-8; + int nrhs = 4; + + /////////////////////////////////////////////// + // Set up the problem as a 4d spreadout job + /////////////////////////////////////////////// + std::vector seeds({1,2,3,4}); + + std::vector src(nrhs,FGrid); + std::vector src_chk(nrhs,FGrid); + std::vector result(nrhs,FGrid); + FermionField tmp(FGrid); + std::cout << GridLogMessage << "Made the Fermion Fields"< HermOp(Ddwf); + ConjugateGradient CG((stp),100000); + + for(int rhs=0;rhs<1;rhs++){ + result[rhs] = zero; + CG(HermOp,src[rhs],result[rhs]); + } + + for(int rhs=0;rhs<1;rhs++){ + std::cout << " Result["< src_v (Ls,UrbGrid); + std::vector result_v(Ls,UrbGrid); + for(int s=0;s