refactoring of gaugeconfig and adjointfield, remove adjointfield next

CMakeLists.txt

@@ -6,7 +6,7 @@ cmake_minimum_required(VERSION 2.8.11)
 project(su2 C CXX)
 project(su2hmc C CXX)
 
-add_definitions(--std=c++11)
+add_definitions(--std=c++14)
 
 ## here we generate version.hh
 ## to track the git version of the executables
@@ -56,7 +56,6 @@ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS} -D_USE_OMP_")
 # generated objects files or create a simple static library. The latter is done
 # here.
 add_library(su2
-    exp_gauge.cc
     print_program_options.cc
     parse_commandline.cc
 )

adjointfield.hh

@@ -120,9 +120,10 @@ template<typename Float> struct adjoint_type<Float, _u1> {
 };
 
 
-template<typename Float, class Group> class adjointfield {
+template<typename Adjoint> class adjointfield {
 public:
-  using value_type = typename adjoint_type<Float, Group>::type;
+  //using value_type = typename adjoint_type<Float, Group>::type;
+  using value_type = Adjoint;
   adjointfield(const size_t Lx, const size_t Ly, const size_t Lz, const size_t Lt, const size_t ndims = 4) : 
     Lx(Lx), Ly(Ly), Lz(Lz), Lt(Lt), volume(Lx*Ly*Lz*Lt), ndims(ndims) {
     data.resize(volume*4);
@@ -161,7 +162,7 @@ public:
   size_t getSize() const {
     return(volume*ndims);
   }
-  void operator=(const adjointfield<Float, Group> &U) {
+  void operator=(const adjointfield<Adjoint> &U) {
     Lx = U.getLx();
     Ly = U.getLz();
     Lz = U.getLz();
@@ -222,15 +223,15 @@ template<typename Float> inline adjointu1<Float> operator*(const Float &x, const
 }
 
 
-template<typename Float, class Group>  adjointfield<Float, su2> operator*(const Float &x, const adjointfield<Float, Group> &A) {
-  adjointfield<Float, Group> res(A.getLx(), A.getLy(), A.getLz(), A.getLt());
+template<typename Float, typename Adjoint>  adjointfield<Adjoint> operator*(const Float &x, const adjointfield<Adjoint> &A) {
+  adjointfield<Adjoint> res(A.getLx(), A.getLy(), A.getLz(), A.getLt());
   for(size_t i = 0; i < A.getSize(); i++) {
     res[i] = x * A[i];
   }
   return res;
 }
 
-template<typename Float> Float operator*(const adjointfield<Float, su2> &A, const adjointfield<Float, su2> &B) {
+template<typename Float> Float operator*(const adjointfield<adjointsu2<Float>> &A, const adjointfield<adjointsu2<Float>> &B) {
   Float res = 0.;
   assert(A.getSize() == B.getSize());
   for(size_t i = 0; i < A.getSize(); i++) {
@@ -239,7 +240,7 @@ template<typename Float> Float operator*(const adjointfield<Float, su2> &A, cons
   return res;
 }
 
-template<typename Float> Float operator*(const adjointfield<Float, _u1> &A, const adjointfield<Float, _u1> &B) {
+template<typename Float> Float operator*(const adjointfield<adjointu1<Float>> &A, const adjointfield<adjointu1<Float>> &B) {
   Float res = 0.;
   assert(A.getSize() == B.getSize());
   for(size_t i = 0; i < A.getSize(); i++) {
@@ -249,7 +250,7 @@ template<typename Float> Float operator*(const adjointfield<Float, _u1> &A, cons
 }
 
 
-template<class URNG, typename Float> void initnormal(URNG &engine, adjointfield<Float, su2> &A) {
+template<class URNG, typename Float> void initnormal(URNG &engine, adjointfield<adjointsu2<Float>> &A) {
   std::normal_distribution<double> normal(0., 1.);
   for(size_t i = 0; i < A.getSize(); i++) {
     A[i].seta(Float(normal(engine)));
@@ -259,15 +260,15 @@ template<class URNG, typename Float> void initnormal(URNG &engine, adjointfield<
   return;
 }
 
-template<class URNG, typename Float> void initnormal(URNG &engine, adjointfield<Float, _u1> &A) {
+template<class URNG, typename Float> void initnormal(URNG &engine, adjointfield<adjointu1<Float>> &A) {
   std::normal_distribution<double> normal(0., 1.);
   for(size_t i = 0; i < A.getSize(); i++) {
     A[i].seta(Float(normal(engine)));
   }
   return;
 }
 
-template<typename Float, class Group> inline void zeroadjointfield(adjointfield<Float, Group> &A) {
+template<typename Adjoint> inline void zeroadjointfield(adjointfield<Adjoint> &A) {
   for(size_t i = 0; i < A.getSize(); i++) {
     A[i].setzero();
   }

exp_gauge.cc

@@ -5,16 +5,5 @@
 #include<vector>
 #include<complex>
 
-_su2 exp(adjointsu2<double> const & x) {
-  double a = x.geta(), b = x.getb(), c = x.getc();
-  // normalise the vector
-  const double alpha = sqrt(a*a+b*b+c*c);
-  std::vector<double> n = {a/alpha, b/alpha, c/alpha};
-
-  const double salpha = sin(alpha);
-  _su2 res(Complex(cos(alpha), salpha*n[2]), 
-           salpha*Complex(n[1], n[0]));
-  return res;
-}
 
 

exp_gauge.hh

@@ -4,7 +4,18 @@
 #include "u1.hh"
 #include "adjointfield.hh"
 
-_su2 exp(adjointsu2<double> const & x);
+template<typename Float=double>
+_su2<Float> exp(adjointsu2<Float> const & x) {
+  double a = x.geta(), b = x.getb(), c = x.getc();
+  // normalise the vector
+  const double alpha = sqrt(a*a+b*b+c*c);
+  std::vector<double> n = {a/alpha, b/alpha, c/alpha};
+
+  const double salpha = sin(alpha);
+  _su2<Float> res(typename _su2<Float>::Complex(cos(alpha), salpha*n[2]), 
+                  salpha*(typename _su2<Float>::Complex(n[1], n[0])));
+  return res;
+}
 
 inline _u1 exp(adjointu1<double> const & x) {
   return _u1(x.geta());

gaugeconfig.hh

@@ -1,7 +1,6 @@
 #pragma once
 
-#include"su2.hh"
-#include"u1.hh"
+#include"site_types.hh"
 #include"random_element.hh"
 #include<random>
 #include<vector>
@@ -10,7 +9,7 @@
 #include<complex>
 #include<iostream>
 #include<cassert>
-
+#include<typeinfo>
 using std::vector;
 
 template<class T> class gaugeconfig {
@@ -68,6 +67,14 @@ public:
     }
   }
 
+  //  template<class Q=T>
+  template<class Q = T, typename = std::enable_if_t< (typeid(Q) == typeid(su2)) > >
+  void flipsign() {
+    for(size_t i = 0; i < getSize(); i++) {
+      data[i].flipsign();
+    }
+  }
+
   void operator=(const gaugeconfig &U) {
     volume = U.getVolume();
     data.resize(U.getSize());

gaugemonomial.hh

@@ -24,7 +24,7 @@ public:
     monomial<Float, Group>::Hnew = h.U->getBeta()*(h.U->getVolume()*6 - gauge_energy(*(h.U))/double(h.U->getNc()));
     return;
   }
-  void derivative(adjointfield<Float, Group> &deriv, hamiltonian_field<Float, Group> const &h, const Float fac = 1.) const override {
+  void derivative(adjointfield<typename adjoint_type<Float, Group>::type> &deriv, hamiltonian_field<Float, Group> const &h, const Float fac = 1.) const override {
     std::vector<size_t> x = {0, 0, 0, 0};
     typedef typename accum_type<Group>::type accum;
 #pragma omp parallel for

gradient_flow.hh

@@ -63,7 +63,7 @@ template<class Group> void gradient_flow(gaugeconfig<Group> &U, std::string cons
   E[2] = density;
 
   gaugeconfig<Group> Vt(U);
-  adjointfield<double, Group> deriv(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims());
+  adjointfield<typename adjoint_type<double, Group>::type> deriv(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims());
   hamiltonian_field<double, Group> h(deriv, Vt);
   gaugemonomial<double, Group> SW(0);
 

hamiltonian_field.hh

@@ -5,8 +5,8 @@
 #include<vector>
 
 template<typename Float, class Group> struct hamiltonian_field {
-  adjointfield<Float, Group> * momenta;
+  adjointfield<typename adjoint_type<Float, Group>::type> * momenta;
   gaugeconfig<Group> * U;
-  hamiltonian_field(adjointfield<Float, Group> &momenta, gaugeconfig<Group> &U) :
+  hamiltonian_field(adjointfield<typename adjoint_type<Float, Group>::type> &momenta, gaugeconfig<Group> &U) :
     momenta(&momenta), U(&U) {}
 };

integrator.hh

@@ -28,7 +28,7 @@ public:
   void integrate(std::list<monomial<Float, Group>*> &monomial_list, hamiltonian_field<Float, Group> &h, 
                  md_params const &params, const bool restore=true) {
 
-    adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
+    adjointfield<typename adjoint_type<Float, Group>::type> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
 
     Float dtau = params.gettau()/Float(params.getnsteps());
     // initial half-step for the  momenta
@@ -55,7 +55,7 @@ public:
   lp_leapfrog(size_t n) : n_prec(n) {}
   void integrate(std::list<monomial<Float, Group>*> &monomial_list, hamiltonian_field<Float, Group> &h, 
                  md_params const &params, const bool restore = true) {
-    adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
+    adjointfield<typename adjoint_type<Float, Group>::type> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
     const size_t N = pow(10, n_prec);
 
     Float dtau = params.gettau()/Float(params.getnsteps());
@@ -85,7 +85,7 @@ public:
   void integrate(std::list<monomial<Float, Group>*> &monomial_list, hamiltonian_field<Float, Group> &h, 
                  md_params const &params, const bool restore = true) {
 
-    adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
+    adjointfield<typename adjoint_type<Float, Group>::type> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
     Float dtau = params.gettau()/Float(params.getnsteps());
     Float oneminus2lambda = (1.-2.*lambda);
 
@@ -121,7 +121,7 @@ public:
   void integrate(std::list<monomial<Float, Group>*> &monomial_list, hamiltonian_field<Float, Group> &h, 
                  md_params const &params, const bool restore = true) {
 
-    adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
+    adjointfield<typename adjoint_type<Float, Group>::type> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
     Float dtau = params.gettau()/Float(params.getnsteps());
     Float eps[10] = {rho*dtau, lambda*dtau, 
                  theta*dtau, 0.5*(1-2.*(lambda+vartheta))*dtau, 
@@ -161,7 +161,7 @@ public:
   void integrate(std::list<monomial<Float, Group>*> &monomial_list, hamiltonian_field<Float, Group> &h, 
                  md_params const &params, const bool restore = true) {
 
-    adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
+    adjointfield<typename adjoint_type<Float, Group>::type> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
     const size_t N = pow(10, n_prec);
 
     Float dtau = params.gettau()/Float(params.getnsteps());
@@ -203,7 +203,7 @@ public:
   void integrate(std::list<monomial<Float, Group>*> &monomial_list, hamiltonian_field<Float, Group> &h, 
                  md_params const &params, const bool restore = true) {
 
-    adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
+    adjointfield<typename adjoint_type<Float, Group>::type> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
     Float dtau = params.gettau()/Float(params.getnsteps());
     // nsteps full steps
     for(size_t i = 0; i < params.getnsteps(); i++) {
@@ -222,7 +222,7 @@ public:
   void integrate(std::list<monomial<Float, Group>*> &monomial_list, hamiltonian_field<Float, Group> &h, 
                  md_params const &params, const bool restore = true) {
 
-    adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
+    adjointfield<typename adjoint_type<Float, Group>::type> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
 
     Float dtau = params.gettau()/Float(params.getnsteps());
     // nsteps full steps

kramers_md_update.hh

@@ -22,7 +22,7 @@ template<class URNG, typename Float, class Group> void kramers_md_update(gaugeco
                                                                          md_params &params,
                                                                          std::list<monomial<Float, Group>*> &monomial_list, 
                                                                          integrator<Float, Group> &md_integ) {
-  adjointfield<Float, Group> momenta(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims());
+  adjointfield<typename adjoint_type<Float, Group>::type> momenta(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims());
   // generate standard normal distributed random momenta
   // normal distribution checked!
   initnormal(engine, momenta);
@@ -31,8 +31,8 @@ template<class URNG, typename Float, class Group> void kramers_md_update(gaugeco
   std::uniform_real_distribution<Float> uniform(0., 1.);
 
   const double gamma = params.getgamma();
-  adjointfield<Float, Group> eta(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims());
-  adjointfield<Float, Group> momenta_old(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims());
+  adjointfield<typename adjoint_type<Float, Group>::type> eta(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims());
+  adjointfield<typename adjoint_type<Float, Group>::type> momenta_old(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims());
   gaugeconfig<su2> U_old(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims(), U.getBeta());
 
   for(size_t k = 0; k < params.getkmax(); k++) {

lyapunov.hh

@@ -30,10 +30,10 @@ template<typename Float, class URNG, class Group> void compute_lyapunov(gaugecon
 
   const size_t n = pow(10, d);
 
-  adjointfield<Float, Group> momenta(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims());
+  adjointfield<typename adjoint_type<Float, Group>::type> momenta(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims());
   // generate standard normal distributed random momenta
   initnormal(engine, momenta);
-  adjointfield<Float, Group> momenta2(momenta);
+  adjointfield<typename adjoint_type<Float, Group>::type> momenta2(momenta);
 
   // generate copy of U, but round to d decimal digits
   gaugeconfig<Group> U2(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getBeta());

md_update.hh

@@ -20,7 +20,7 @@ template<class URNG, typename Float, class Group> void md_update(gaugeconfig<Gro
                                                                  md_params &params,
                                                                  std::list<monomial<Float, Group>*> &monomial_list, 
                                                                  integrator<Float, Group> &md_integ) {
-  adjointfield<Float, Group> momenta(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims());
+  adjointfield<typename adjoint_type<Float, Group>::type> momenta(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims());
   // generate standard normal distributed random momenta
   // normal distribution checked!
   initnormal(engine, momenta);

monomial.hh

@@ -10,7 +10,7 @@ public:
   monomial() : Hold(0.), Hnew(0.), timescale(0), mdactive(true) {}
   virtual void heatbath(hamiltonian_field<Float, Group> const &h) = 0;
   virtual void accept(hamiltonian_field<Float, Group> const &h) = 0;
-  virtual void derivative(adjointfield<Float, Group> &deriv, hamiltonian_field<Float, Group> const &h, const Float fac) const = 0;
+  virtual void derivative(adjointfield<typename adjoint_type<Float, Group>::type> &deriv, hamiltonian_field<Float, Group> const &h, const Float fac) const = 0;
   void reset() {
     Hold = 0.;
     Hnew = 0.;
@@ -47,7 +47,7 @@ public:
   void accept(hamiltonian_field<Float, Group> const &h) override {
     monomial<Float, Group>::Hnew = 0.5 * ((*h.momenta)*(*h.momenta));
   }
-  virtual void derivative(adjointfield<Float, Group> &deriv, hamiltonian_field<Float, Group> const &h, const Float fac) const {}
+  virtual void derivative(adjointfield<typename adjoint_type<Float, Group>::type> &deriv, hamiltonian_field<Float, Group> const &h, const Float fac) const {}
 };
 
 #include"gaugemonomial.hh"

site_types.hh

@@ -0,0 +1,5 @@
+#pragma once
+
+#include"su2.hh"
+#include"u1.hh"
+#include"adjointfield.hh"