diff --git a/INSTALL b/INSTALL index 80a61507..f812f5a2 120000 --- a/INSTALL +++ b/INSTALL @@ -1 +1 @@ -/opt/local/share/automake-1.15/INSTALL \ No newline at end of file +/usr/share/automake-1.14/INSTALL \ No newline at end of file diff --git a/configure b/configure index 06ca42d5..935f1018 100755 --- a/configure +++ b/configure @@ -2577,7 +2577,7 @@ test -n "$target_alias" && NONENONEs,x,x, && program_prefix=${target_alias}- -am__api_version='1.15' +am__api_version='1.14' # Find a good install program. We prefer a C program (faster), # so one script is as good as another. But avoid the broken or @@ -2749,8 +2749,8 @@ test "$program_suffix" != NONE && ac_script='s/[\\$]/&&/g;s/;s,x,x,$//' program_transform_name=`$as_echo "$program_transform_name" | sed "$ac_script"` -# Expand $ac_aux_dir to an absolute path. -am_aux_dir=`cd "$ac_aux_dir" && pwd` +# expand $ac_aux_dir to an absolute path +am_aux_dir=`cd $ac_aux_dir && pwd` if test x"${MISSING+set}" != xset; then case $am_aux_dir in @@ -2769,7 +2769,7 @@ else $as_echo "$as_me: WARNING: 'missing' script is too old or missing" >&2;} fi -if test x"${install_sh+set}" != xset; then +if test x"${install_sh}" != xset; then case $am_aux_dir in *\ * | *\ *) install_sh="\${SHELL} '$am_aux_dir/install-sh'" ;; @@ -3097,8 +3097,8 @@ MAKEINFO=${MAKEINFO-"${am_missing_run}makeinfo"} # mkdir_p='$(MKDIR_P)' -# We need awk for the "check" target (and possibly the TAP driver). The -# system "awk" is bad on some platforms. +# We need awk for the "check" target. The system "awk" is bad on +# some platforms. # Always define AMTAR for backward compatibility. Yes, it's still used # in the wild :-( We should find a proper way to deprecate it ... AMTAR='$${TAR-tar}' @@ -3157,7 +3157,6 @@ END fi - ac_config_headers="$ac_config_headers lib/GridConfig.h" # Check whether --enable-silent-rules was given. diff --git a/lib/GridConfig.h b/lib/GridConfig.h index c9dea1cb..f3a1be6a 100644 --- a/lib/GridConfig.h +++ b/lib/GridConfig.h @@ -19,6 +19,12 @@ /* GRID_COMMS_NONE */ #define GRID_COMMS_NONE 1 +/* GRID_DEFAULT_PRECISION is DOUBLE */ +/* #undef GRID_DEFAULT_PRECISION_DOUBLE */ + +/* GRID_DEFAULT_PRECISION is SINGLE */ +#define GRID_DEFAULT_PRECISION_SINGLE 1 + /* Support Altivec instructions */ /* #undef HAVE_ALTIVEC */ diff --git a/lib/Make.inc b/lib/Make.inc index a1fe777e..18a29228 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/MultiShiftFunction.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/CoarsenedMatrix.h ./algorithms/iterative/AdefGeneric.h ./algorithms/iterative/BfmHDCG.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/PrecConjugateResidual.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.h ./algorithms/Preconditioner.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 ./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_comparison_utils.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_unary.h ./lattice/Lattice_where.h ./Lattice.h ./Old/Tensor_peek.h ./Old/Tensor_poke.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/g5HermitianLinop.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/SUn.h ./qcd/utils/WilsonLoops.h ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./simd/Grid_sse4.h ./simd/Grid_vector_types.h ./simd/Grid_vector_unops.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_determinant.h ./tensors/Tensor_exp.h ./tensors/Tensor_extract_merge.h ./tensors/Tensor_index.h ./tensors/Tensor_inner.h ./tensors/Tensor_logical.h ./tensors/Tensor_outer.h ./tensors/Tensor_reality.h ./tensors/Tensor_Ta.h ./tensors/Tensor_trace.h ./tensors/Tensor_traits.h ./tensors/Tensor_transpose.h ./tensors/Tensor_unary.h ./Tensors.h ./Threads.h +HFILES=./Cshift.h ./simd/Grid_avx.h ./simd/Grid_vector_types.h ./simd/Grid_sse4.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_vector_unops.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./Tensors.h ./Algorithms.h ./communicator/Communicator_base.h ./lattice/Lattice_rng.h ./lattice/Lattice_reduction.h ./lattice/Lattice_transfer.h ./lattice/Lattice_unary.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_comparison.h ./lattice/Lattice_overload.h ./lattice/Lattice_reality.h ./lattice/Lattice_local.h ./lattice/Lattice_conformable.h ./lattice/Lattice_where.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_ET.h ./lattice/Lattice_transpose.h ./lattice/Lattice_trace.h ./Stencil.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_exp.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_class.h ./tensors/Tensor_logical.h ./tensors/Tensor_transpose.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_reality.h ./tensors/Tensor_trace.h ./tensors/Tensor_index.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_outer.h ./tensors/Tensor_inner.h ./tensors/Tensor_traits.h ./tensors/Tensor_Ta.h ./tensors/Tensor_unary.h ./tensors/Tensor_determinant.h ./tensors/Tensor_arith.h ./tensors/Tensor_extract_merge.h ./Communicator.h ./Cartesian.h ./parallelIO/NerscIO.h ./qcd/QCD.h ./qcd/hmc/integrators/Integrator_base.h ./qcd/hmc/HMC.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/LinalgUtils.h ./qcd/utils/CovariantCshift.h ./qcd/utils/WilsonLoops.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/action/gauge/GaugeActionBase.h ./qcd/action/Actions.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/spin/TwoSpinor.h ./qcd/spin/Dirac.h ./cshift/Cshift_common.h ./cshift/Cshift_none.h ./cshift/Cshift_mpi.h ./Simd.h ./GridConfig.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_red_black.h ./cartesian/Cartesian_full.h ./AlignedAllocator.h ./Lattice.h ./Old/Tensor_poke.h ./Old/Tensor_peek.h ./Threads.h ./Grid.h ./algorithms/Preconditioner.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/AdefGeneric.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/Zolotarev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Remez.h ./algorithms/LinearOperator.h ./algorithms/SparseMatrix.h ./algorithms/CoarsenedMatrix.h ./stencil/Lebesgue.h -CCFILES=./algorithms/approx/MultiShiftFunction.cc ./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 +CCFILES=./qcd/hmc/integrators/Integrator_base.cc ./qcd/hmc/HMC.cc ./qcd/utils/SpaceTimeGrid.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/spin/Dirac.cc ./GridInit.cc ./algorithms/approx/MultiShiftFunction.cc ./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc diff --git a/lib/qcd/QCD.h b/lib/qcd/QCD.h index 0e841d8e..cf8ff6d8 100644 --- a/lib/qcd/QCD.h +++ b/lib/qcd/QCD.h @@ -413,5 +413,8 @@ namespace QCD { #include #include #include +#include +#include + #endif diff --git a/lib/qcd/hmc/.dirstamp b/lib/qcd/hmc/.dirstamp new file mode 100644 index 00000000..e69de29b diff --git a/lib/qcd/hmc/HMC.cc b/lib/qcd/hmc/HMC.cc new file mode 100644 index 00000000..32ef8ca0 --- /dev/null +++ b/lib/qcd/hmc/HMC.cc @@ -0,0 +1,98 @@ +#include + +namespace Grid{ + namespace QCD{ + + HMCparameters::HMCparameters(){ + // FIXME fill this constructor now just default values + + ////////////////////////////// Default values + Nsweeps = 10; + TotalSweeps = 10; + ThermalizationSteps = 0; + StartingConfig = 0; + SaveInterval = 1; + Filename_prefix = "Conf_"; + ///////////////////////////////// + + } + + + ///////////////////////////////////////////////////////////////// + + HybridMonteCarlo::HybridMonteCarlo(GridSerialRNG& R):RNG(R){ + //FIXME + } + + + void HybridMonteCarlo::evolve(LatticeColourMatrix& Uin)const{ + Real DeltaH; + Real timer; + + // Thermalizations + for(int iter=1; iter <= Params.ThermalizationSteps; ++iter){ + std::cout << "-- # Thermalization step = "<< iter << "\n"; + + DeltaH = evolve_step(Uin); + + std::cout<< "[Timing] Trajectory time (s) : "<< timer/1000.0 << "\n"; + std::cout<< "dH = "<< DeltaH << "\n"; + + // Update matrix + //Uin = md_->get_U(); //accept every time + } + + // Actual updates + for(int iter=Params.StartingConfig; + iter < Params.Nsweeps+Params.StartingConfig; ++iter){ + std::cout << "-- # Sweep = "<< iter << "\n"; + + DeltaH = evolve_step(Uin); + + if(metropolis_test(DeltaH)) {};//Uin = md_->get_U(); + // need sync? + } + + + } + + + RealD HybridMonteCarlo::evolve_step(LatticeColourMatrix& Uin)const{ + /* + md_->init(Uin,RNG); // set U and initialize P and phi's + RealD H0 = md_->S(); // current state + std::cout<<"Total H_before = "<< H0 << "\n"; + + md_->integrator(); + + RealD H1 = md_->calc_H(); // updated state + std::cout<<"Total H_after = "<< H1 << "\n"; + + return (H1-H0); + */ + return 0; + } + + + + + bool HybridMonteCarlo::metropolis_test(const RealD DeltaH)const{ + RealD rn_test; + RealD prob = std::exp(-DeltaH); + random(RNG,rn_test); + + std::cout<< "--------------------------------------------\n"; + std::cout<< "dH = "<1.0) || (rn_test <= prob)){ // accepted + std::cout <<"-- ACCEPTED\n"; + return true; + } else { // rejected + std::cout <<"-- REJECTED\n"; + return false; + } + } + + } +} diff --git a/lib/qcd/hmc/HMC.h b/lib/qcd/hmc/HMC.h new file mode 100644 index 00000000..e25bb3ed --- /dev/null +++ b/lib/qcd/hmc/HMC.h @@ -0,0 +1,50 @@ +//-------------------------------------------------------------------- +/*! @file HMC.h + * @brief Declaration of classes for HybridMonteCarlo update + * + * @author Guido Cossu + */ +//-------------------------------------------------------------------- +#ifndef HMC_INCLUDED +#define HMC_INCLUDED + +#include +#include + +namespace Grid{ + namespace QCD{ + + struct HMCparameters{ + Integer Nsweeps; /* @brief Number of sweeps in this run */ + Integer TotalSweeps; /* @brief If provided, the total number of sweeps */ + Integer ThermalizationSteps; + Integer StartingConfig; + Integer SaveInterval; //Setting to 0 does not save configurations + std::string Filename_prefix; // To save configurations + + HMCparameters(); + }; + + + class HybridMonteCarlo{ + const HMCparameters Params; + GridSerialRNG& RNG; + // FIXME need the integrator + + bool metropolis_test(const RealD DeltaH)const; + RealD evolve_step(LatticeColourMatrix&)const; + + + public: + HybridMonteCarlo(GridSerialRNG&); + ~HybridMonteCarlo(){}; + + void evolve(LatticeColourMatrix& Uin)const; + + }; + + }// QCD +}// Grid + + +#endif diff --git a/lib/qcd/hmc/integrators/.dirstamp b/lib/qcd/hmc/integrators/.dirstamp new file mode 100644 index 00000000..e69de29b diff --git a/lib/qcd/hmc/integrators/Integrator_base.cc b/lib/qcd/hmc/integrators/Integrator_base.cc new file mode 100644 index 00000000..5fc89482 --- /dev/null +++ b/lib/qcd/hmc/integrators/Integrator_base.cc @@ -0,0 +1,22 @@ +/*! + @file Integrator_base.cc + @brief utilities for MD including funcs to generate initial HMC momentum + */ + + +#include + +static const double sq3i = 1.0/sqrt(3.0); + +namespace Grid{ + namespace QCD{ + + + + void MDutils::generate_momenta_su3(LatticeColourMatrix& P,GridParallelRNG& pRNG){ + SU3::GaussianLieAlgebraMatrix(pRNG, P); + } + + + } +} diff --git a/lib/qcd/hmc/integrators/Integrator_base.h b/lib/qcd/hmc/integrators/Integrator_base.h new file mode 100644 index 00000000..f1d6f711 --- /dev/null +++ b/lib/qcd/hmc/integrators/Integrator_base.h @@ -0,0 +1,50 @@ +//-------------------------------------------------------------------- +/*! @file Integrator_base.h + * @brief Declaration of classes for the abstract Molecular Dynamics integrator + * + * @author Guido Cossu + */ +//-------------------------------------------------------------------- + +#ifndef INTEGRATOR_INCLUDED +#define INTEGRATOR_INCLUDED + +class Action; +class RandNum; +class Observer; + +typedef std::vector ActionLevel; +typedef std::vector ActionSet; +typedef std::vector ObserverList; + +/*! @brief Abstract base class for Molecular Dynamics management */ + +namespace Grid{ + namespace QCD{ + + class Integrator{ + private: + virtual void update_P(int lv,double ep) = 0; + virtual void update_U(double ep) = 0; + + virtual void register_observers() = 0; + virtual void notify_observers() = 0; + + public: + virtual ~Integrator(){} + virtual void init(const LatticeColourMatrix&, + const GridParallelRNG& RNG)=0; + virtual double S()const =0; + virtual void integrate(int level) =0; + virtual const LatticeColourMatrix get_U() const =0; + + void generate_momenta(LatticeColourMatrix& P,const RandNum& rand); + }; + + namespace MDutils{ + void generate_momenta_su3(LatticeColourMatrix& P,GridParallelRNG& RNG); + } + + } +} +#endif//INTEGRATOR_INCLUDED diff --git a/lib/qcd/utils/SUn.h b/lib/qcd/utils/SUn.h index 69de3f9f..9a34f3e3 100644 --- a/lib/qcd/utils/SUn.h +++ b/lib/qcd/utils/SUn.h @@ -543,6 +543,26 @@ Note that in step D setting B ~ X - A and using B in place of A in step E will g taExp(lie,out); } + static void GaussianLieAlgebraMatrix(GridParallelRNG &pRNG,LatticeMatrix &out,double scale=1.0){ + GridBase *grid = out._grid; + LatticeComplex ca (grid); + LatticeMatrix lie(grid); + LatticeMatrix la (grid); + Complex ci(0.0,scale); + Matrix ta; + + out=zero; + for(int a=0;a