diff --git a/lib/Algorithms.h b/lib/Algorithms.h index b5a2814a..c2491cf4 100644 --- a/lib/Algorithms.h +++ b/lib/Algorithms.h @@ -5,6 +5,7 @@ #include #include +#include #include #include diff --git a/lib/Make.inc b/lib/Make.inc index 94a566fb..ba35e25f 100644 --- a/lib/Make.inc +++ b/lib/Make.inc @@ -1,4 +1,4 @@ -HFILES=./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.h ./algorithms/SparseMatrix.h ./Algorithms.h ./AlignedAllocator.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_full.h ./cartesian/Cartesian_red_black.h ./Cartesian.h ./communicator/Communicator_base.h ./Communicator.h ./Comparison.h ./cshift/Cshift_common.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./Cshift.h ./Grid.h ./GridConfig.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_comparison.h ./lattice/Lattice_conformable.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_ET.h ./lattice/Lattice_local.h ./lattice/Lattice_overload.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_reality.h ./lattice/Lattice_reduction.h ./lattice/Lattice_rng.h ./lattice/Lattice_trace.h ./lattice/Lattice_transfer.h ./lattice/Lattice_transpose.h ./lattice/Lattice_where.h ./Lattice.h ./parallelIO/NerscIO.h ./qcd/action/Actions.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/gauge/GaugeActionBase.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/Dirac.h ./qcd/LinalgUtils.h ./qcd/QCD.h ./qcd/SpaceTimeGrid.h ./qcd/TwoSpinor.h ./qcd/utils/CovariantCshift.h ./qcd/utils/WilsonLoops.h ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_qpx.h ./simd/Grid_sse4.h ./simd/Grid_vector_types.h ./simd/Old/Grid_vComplexD.h ./simd/Old/Grid_vComplexF.h ./simd/Old/Grid_vInteger.h ./simd/Old/Grid_vRealD.h ./simd/Old/Grid_vRealF.h ./Simd.h ./stencil/Lebesgue.h ./Stencil.h ./tensors/Tensor_arith.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_class.h ./tensors/Tensor_extract_merge.h ./tensors/Tensor_inner.h ./tensors/Tensor_outer.h ./tensors/Tensor_peek.h ./tensors/Tensor_poke.h ./tensors/Tensor_reality.h ./tensors/Tensor_Ta.h ./tensors/Tensor_trace.h ./tensors/Tensor_traits.h ./tensors/Tensor_transpose.h ./Tensors.h ./Threads.h +HFILES=./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.h ./algorithms/SparseMatrix.h ./Algorithms.h ./AlignedAllocator.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_full.h ./cartesian/Cartesian_red_black.h ./Cartesian.h ./communicator/Communicator_base.h ./Communicator.h ./Comparison.h ./cshift/Cshift_common.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./Cshift.h ./Grid.h ./GridConfig.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_comparison.h ./lattice/Lattice_conformable.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_ET.h ./lattice/Lattice_local.h ./lattice/Lattice_overload.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_reality.h ./lattice/Lattice_reduction.h ./lattice/Lattice_rng.h ./lattice/Lattice_trace.h ./lattice/Lattice_transfer.h ./lattice/Lattice_transpose.h ./lattice/Lattice_where.h ./Lattice.h ./parallelIO/NerscIO.h ./qcd/action/Actions.h ./qcd/action/DiffAction.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/gauge/GaugeActionBase.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/QCD.h ./qcd/spin/Dirac.h ./qcd/spin/TwoSpinor.h ./qcd/utils/CovariantCshift.h ./qcd/utils/LinalgUtils.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/WilsonLoops.h ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_qpx.h ./simd/Grid_sse4.h ./simd/Grid_vector_types.h ./simd/Old/Grid_vComplexD.h ./simd/Old/Grid_vComplexF.h ./simd/Old/Grid_vInteger.h ./simd/Old/Grid_vRealD.h ./simd/Old/Grid_vRealF.h ./Simd.h ./stencil/Lebesgue.h ./Stencil.h ./tensors/Tensor_arith.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_class.h ./tensors/Tensor_extract_merge.h ./tensors/Tensor_inner.h ./tensors/Tensor_outer.h ./tensors/Tensor_peek.h ./tensors/Tensor_poke.h ./tensors/Tensor_reality.h ./tensors/Tensor_Ta.h ./tensors/Tensor_trace.h ./tensors/Tensor_traits.h ./tensors/Tensor_transpose.h ./Tensors.h ./Threads.h -CCFILES=./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./GridInit.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/Dirac.cc ./qcd/SpaceTimeGrid.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc +CCFILES=./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./GridInit.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/spin/Dirac.cc ./qcd/utils/SpaceTimeGrid.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc diff --git a/lib/algorithms/LinearOperator.h b/lib/algorithms/LinearOperator.h index 710fcde7..492b34db 100644 --- a/lib/algorithms/LinearOperator.h +++ b/lib/algorithms/LinearOperator.h @@ -18,7 +18,8 @@ namespace Grid { public: virtual void Op (const Field &in, Field &out) = 0; // Abstract base virtual void AdjOp (const Field &in, Field &out) = 0; // Abstract base - virtual void HermOpAndNorm(const Field &in, Field &out,double &n1,double &n2)=0; + virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2)=0; + virtual void HermOp(const Field &in, Field &out)=0; }; @@ -48,9 +49,13 @@ namespace Grid { void AdjOp (const Field &in, Field &out){ _Mat.Mdag(in,out); } - void HermOpAndNorm(const Field &in, Field &out,double &n1,double &n2){ + void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ _Mat.MdagM(in,out,n1,n2); } + void HermOp(const Field &in, Field &out){ + RealD n1,n2; + HermOpAndNorm(in,out,n1,n2); + } }; //////////////////////////////////////////////////////////////////// @@ -67,7 +72,7 @@ namespace Grid { void AdjOp (const Field &in, Field &out){ _Mat.M(in,out); } - void HermOpAndNorm(const Field &in, Field &out,double &n1,double &n2){ + void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ ComplexD dot; _Mat.M(in,out); @@ -78,6 +83,9 @@ namespace Grid { dot = innerProduct(out,out); n2=real(dot); } + void HermOp(const Field &in, Field &out){ + _Mat.M(in,out); + } }; ////////////////////////////////////////////////////////// @@ -98,6 +106,10 @@ namespace Grid { void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ MpcDagMpc(in,out,n1,n2); } + void HermOp(const Field &in, Field &out){ + RealD n1,n2; + HermOpAndNorm(in,out,n1,n2); + } void Op (const Field &in, Field &out){ Mpc(in,out); } diff --git a/lib/algorithms/SparseMatrix.h b/lib/algorithms/SparseMatrix.h index 9146648f..a51657dc 100644 --- a/lib/algorithms/SparseMatrix.h +++ b/lib/algorithms/SparseMatrix.h @@ -18,6 +18,7 @@ namespace Grid { Field tmp (in._grid); ni=M(in,tmp); no=Mdag(tmp,out); + std::cout << "MdagM "<< ni<<" "< + class ConjugateResidual : public OperatorFunction { + public: + RealD Tolerance; + Integer MaxIterations; + int verbose; + + ConjugateResidual(RealD tol,Integer maxit) : Tolerance(tol), MaxIterations(maxit) { + verbose=1; + }; + + void operator() (LinearOperatorBase &Linop,const Field &src, Field &psi){ + + RealD a, b, c, d; + RealD cp, ssq,rsq; + + RealD rAr, rAAr, rArp; + RealD pAp, pAAp; + + GridBase *grid = src._grid; + psi=zero; + Field r(grid), p(grid), Ap(grid), Ar(grid); + + r=src; + p=src; + + Linop.HermOpAndNorm(p,Ap,pAp,pAAp); + Linop.HermOpAndNorm(r,Ar,rAr,rAAr); + + std::cout << "pAp, pAAp"<< pAp<<" "< -#include -#include -#include +#include +#include +#include +#include #include #include #include diff --git a/lib/qcd/action/Actions.h b/lib/qcd/action/Actions.h index edef1b53..56b453c6 100644 --- a/lib/qcd/action/Actions.h +++ b/lib/qcd/action/Actions.h @@ -65,36 +65,9 @@ #include #include - - // Chroma interface defining FermionAction - /* - template class FermAct4D : public FermionAction - virtual LinearOperator* linOp(Handle< FermState > state) const = 0; - virtual LinearOperator* lMdagM(Handle< FermState > state) const = 0; - virtual LinOpSystemSolver* invLinOp(Handle< FermState > state, - virtual MdagMSystemSolver* invMdagM(Handle< FermState > state, - virtual LinOpMultiSystemSolver* mInvLinOp(Handle< FermState > state, - virtual MdagMMultiSystemSolver* mInvMdagM(Handle< FermState > state, - virtual MdagMMultiSystemSolverAccumulate* mInvMdagMAcc(Handle< FermState > state, - virtual SystemSolver* qprop(Handle< FermState > state, - class DiffFermAct4D : public FermAct4D - virtual DiffLinearOperator* linOp(Handle< FermState > state) const = 0; - virtual DiffLinearOperator* lMdagM(Handle< FermState > state) const = 0; - */ - - - // Chroma interface defining GaugeAction - /* - template class GaugeAction - virtual const CreateGaugeState& getCreateState() const = 0; - virtual GaugeState* createState(const Q& q) const - virtual const GaugeBC& getGaugeBC() const - virtual const Set& getSet(void) const = 0; - virtual void deriv(P& result, const Handle< GaugeState >& state) const - virtual Double S(const Handle< GaugeState >& state) const = 0; - - class LinearGaugeAction : public GaugeAction< multi1d, multi1d > - typedef multi1d P; - */ +/////////////////////////////////////////////////////////////////////////////// +// G5 herm -- this has to live in QCD since dirac matrix is not in the broader sector of code +/////////////////////////////////////////////////////////////////////////////// +#include #endif diff --git a/lib/qcd/action/fermion/g5HermitianLinop.h b/lib/qcd/action/fermion/g5HermitianLinop.h new file mode 100644 index 00000000..9a8cc2eb --- /dev/null +++ b/lib/qcd/action/fermion/g5HermitianLinop.h @@ -0,0 +1,38 @@ +#ifndef G5_HERMITIAN_LINOP +#define G5_HERMITIAN_LINOP +namespace Grid { + namespace QCD { +//////////////////////////////////////////////////////////////////// +// Wrap an already herm matrix +//////////////////////////////////////////////////////////////////// +template +class Gamma5HermitianLinearOperator : public LinearOperatorBase { + Matrix &_Mat; +public: + Gamma5HermitianLinearOperator(Matrix &Mat): _Mat(Mat){}; + void Op (const Field &in, Field &out){ + _Mat.M(in,out); + } + void AdjOp (const Field &in, Field &out){ + _Mat.M(in,out); + } + void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ + + HermOp(in,out); + + ComplexD dot; + dot= innerProduct(in,out); + n1=real(dot); + + dot = innerProduct(out,out); + n2=real(dot); + } + void HermOp(const Field &in, Field &out){ + Field tmp(in._grid); + _Mat.M(in,tmp); + G5R5(out,tmp); + } +}; + +}} +#endif diff --git a/lib/qcd/Dirac.cc b/lib/qcd/spin/Dirac.cc similarity index 100% rename from lib/qcd/Dirac.cc rename to lib/qcd/spin/Dirac.cc diff --git a/lib/qcd/Dirac.h b/lib/qcd/spin/Dirac.h similarity index 100% rename from lib/qcd/Dirac.h rename to lib/qcd/spin/Dirac.h diff --git a/lib/qcd/TwoSpinor.h b/lib/qcd/spin/TwoSpinor.h similarity index 100% rename from lib/qcd/TwoSpinor.h rename to lib/qcd/spin/TwoSpinor.h diff --git a/lib/qcd/LinalgUtils.h b/lib/qcd/utils/LinalgUtils.h similarity index 88% rename from lib/qcd/LinalgUtils.h rename to lib/qcd/utils/LinalgUtils.h index 2b83d115..ef85e728 100644 --- a/lib/qcd/LinalgUtils.h +++ b/lib/qcd/utils/LinalgUtils.h @@ -109,5 +109,24 @@ PARALLEL_FOR_LOOP vstream(z._odata[ss+s],tmp); } } + +template +void G5R5(Lattice &z,const Lattice &x) +{ + GridBase *grid=x._grid; + z.checkerboard = x.checkerboard; + conformable(x,z); + int Ls = grid->_rdimensions[0]; +PARALLEL_FOR_LOOP + for(int ss=0;ssoSites();ss+=Ls){ // adds Ls + vobj tmp; + for(int s=0;s