Dirichlet

Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatio.h

@@ -59,7 +59,7 @@ NAMESPACE_BEGIN(Grid);
 
       typedef RationalActionParams Params;
       Params param;
-
+      RealD  RefreshAction;
       //For action evaluation
       MultiShiftFunction ApproxPowerAction   ;  //rational approx for X^{1/inv_pow}
       MultiShiftFunction ApproxNegPowerAction;  //rational approx for X^{-1/inv_pow}
@@ -270,12 +270,18 @@ NAMESPACE_BEGIN(Grid);
 	assert(NumOp.ConstEE() == 1);
 	assert(DenOp.ConstEE() == 1);
 	PhiEven = Zero();
-	std::cout<<GridLogMessage << action_name() << " refresh: starting" << std::endl;
+
+	RefreshAction = norm2( etaOdd );
+        std::cout<<GridLogMessage << action_name() << " refresh: action is " << RefreshAction << std::endl;
       };
 
       //////////////////////////////////////////////////////
       // S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi       
       //////////////////////////////////////////////////////
+      virtual RealD Sinitial(const GaugeField &U) {
+	return RefreshAction;
+      }
+
       virtual RealD S(const GaugeField &U) {
 	std::cout<<GridLogMessage << action_name() << " compute action: starting" << std::endl;
 	ImportGauge(U);

Grid/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h

@@ -50,6 +50,8 @@ NAMESPACE_BEGIN(Grid);
       FermionField PhiOdd;   // the pseudo fermion field for this trajectory
       FermionField PhiEven;  // the pseudo fermion field for this trajectory
 
+      RealD RefreshAction;
+      
     public:
       TwoFlavourEvenOddRatioPseudoFermionAction(FermionOperator<Impl>  &_NumOp, 
                                                 FermionOperator<Impl>  &_DenOp, 
@@ -132,6 +134,9 @@ NAMESPACE_BEGIN(Grid);
         // Even det factors
         DenOp.MooeeDag(etaEven,tmp);
         NumOp.MooeeInvDag(tmp,PhiEven);
+
+	RefreshAction = norm2(etaEven)+norm2(etaOdd);
+	std::cout << " refresh " <<action_name()<< " action "<<RefreshAction<<std::endl;
       };
 
       //////////////////////////////////////////////////////

Grid/qcd/hmc/HMC.h

@@ -189,7 +189,7 @@ private:
     //////////////////////////////////////////////////////////////////////////////////////////////////////
     std::cout << GridLogMessage << "--------------------------------------------------\n";
     std::cout << GridLogMessage << "Compute initial action";
-    RealD H0 = TheIntegrator.S(U);  
+    RealD H0 = TheIntegrator.Sinitial(U);  
     std::cout << GridLogMessage << "--------------------------------------------------\n";
 
     std::streamsize current_precision = std::cout.precision();

Grid/qcd/hmc/integrators/Integrator.h

@@ -419,6 +419,48 @@ public:
     return H;
   }
 
+  struct _Sinitial {
+    template <class FieldType, class Repr>
+    void operator()(std::vector<Action<FieldType>*> repr_set, Repr& Rep, int level, RealD& H) {
+      
+      for (int a = 0; a < repr_set.size(); ++a) {
+        RealD Hterm = repr_set.at(a)->Sinitial(Rep.U);
+        std::cout << GridLogMessage << "Sinitial Level " << level << " term " << a << " H Hirep = " << Hterm << std::endl;
+        H += Hterm;
+
+      }
+    }
+  } Sinitial_hireps{};
+
+  RealD Sinitial(Field& U) 
+  {  // here also U not used
+
+    std::cout << GridLogIntegrator << "Integrator initial action\n";
+
+    RealD H = - FieldImplementation::FieldSquareNorm(P)/HMC_MOMENTUM_DENOMINATOR; // - trace (P*P)/denom
+
+    RealD Hterm;
+
+    // Actions
+    for (int level = 0; level < as.size(); ++level) {
+      for (int actionID = 0; actionID < as[level].actions.size(); ++actionID) {
+        // get gauge field from the SmearingPolicy and
+        // based on the boolean is_smeared in actionID
+        Field& Us = Smearer.get_U(as[level].actions.at(actionID)->is_smeared);
+        std::cout << GridLogMessage << "S [" << level << "][" << actionID << "] action eval " << std::endl;
+	        as[level].actions.at(actionID)->S_timer_start();
+        Hterm = as[level].actions.at(actionID)->Sinitial(Us);
+   	        as[level].actions.at(actionID)->S_timer_stop();
+        std::cout << GridLogMessage << "S [" << level << "][" << actionID << "] H = " << Hterm << std::endl;
+        H += Hterm;
+      }
+      as[level].apply(Sinitial_hireps, Representations, level, H);
+    }
+
+    return H;
+  }
+
+  
   void integrate(Field& U) 
   {
     // reset the clocks

Grid/simd/Grid_vector_types.h

@@ -148,23 +148,6 @@ accelerator_inline Grid_half sfw_float_to_half(float ff) {
   #endif
 #endif
 
-/*
-#ifdef A64FXVLA
-#pragma message("building A64FX VLA")
-#if defined(ARMCLANGCOMPAT)
-  #pragma message("applying data types patch")
-#endif
-#include <arm_sve.h>
-#include "Grid_a64fx-2.h"
-#endif
-
-#ifdef A64FXVLS
-#pragma message("building A64FX VLS")
-#include <arm_sve.h>
-#include "Grid_a64fx-fixedsize.h"
-#endif
-*/
-
 #ifdef SSE4
 #include "Grid_sse4.h"
 #endif
@@ -1147,42 +1130,6 @@ static_assert(sizeof(SIMD_Ftype) == sizeof(SIMD_Itype), "SIMD vector lengths inc
 #endif
 #endif
 
-/////////////////////////////////////////
-// Some traits to recognise the types
-/////////////////////////////////////////
-template <typename T>
-struct is_simd : public std::false_type {};
-template <> struct is_simd<vRealF>     : public std::true_type {};
-template <> struct is_simd<vRealD>     : public std::true_type {};
-template <> struct is_simd<vRealH>     : public std::true_type {};
-template <> struct is_simd<vComplexF>  : public std::true_type {};
-template <> struct is_simd<vComplexD>  : public std::true_type {};
-template <> struct is_simd<vComplexH>  : public std::true_type {};
-template <> struct is_simd<vInteger>   : public std::true_type {};
-
-template <typename T> using IfSimd    = Invoke<std::enable_if<is_simd<T>::value, int> >;
-template <typename T> using IfNotSimd = Invoke<std::enable_if<!is_simd<T>::value, unsigned> >;
-
-///////////////////////////////////////////////
-// Convenience insert / extract with complex support
-///////////////////////////////////////////////
-template <class S, class V>
-accelerator_inline S getlane(const Grid_simd<S, V> &in,int lane) {
-  return in.getlane(lane);
-}
-template <class S, class V>
-accelerator_inline void putlane(Grid_simd<S, V> &vec,const S &_S, int lane){
-  vec.putlane(_S,lane);
-}
-template <class S,IfNotSimd<S> = 0 >
-accelerator_inline S getlane(const S &in,int lane) {
-  return in;
-}
-template <class S,IfNotSimd<S> = 0 >
-accelerator_inline void putlane(S &vec,const S &_S, int lane){
-  vec = _S;
-}
-
 
 NAMESPACE_END(Grid);
 

configure.ac

@@ -1,4 +1,4 @@
-AC_PREREQ([2.71])
+AC_PREREQ([2.69])
 AC_INIT([Grid],[0.7.0],[https://github.com/paboyle/Grid],[Grid])
 AC_CANONICAL_BUILD
 AC_CANONICAL_HOST

systems/Crusher/config-command

@@ -14,4 +14,6 @@ CXXFLAGS="-fPIC -I{$ROCM_PATH}/include/ -std=c++14 -I${MPICH_DIR}/include -L/lib
  LDFLAGS="-L/lib64 -L/opt/rocm-5.2.0/lib/ -L${MPICH_DIR}/lib -lmpi -L${CRAY_MPICH_ROOTDIR}/gtl/lib -lmpi_gtl_hsa -lamdhip64 "
 
 
+#--enable-simd=GPU-RRII \
+