hot start for sp2n

Grid/qcd/QCD.h

@@ -47,6 +47,7 @@ static constexpr int Ym = 5;
 static constexpr int Zm = 6;
 static constexpr int Tm = 7;
 
+static constexpr int isSymplectic=Sp2n_config;
 static constexpr int Nc=Config_Nc;
 static constexpr int Ns=4;
 static constexpr int Nd=4;

Grid/qcd/action/gauge/GaugeImplTypes.h

@@ -191,11 +191,26 @@ public:
 
   static inline void HotConfiguration(GridParallelRNG &pRNG, Field &U)
   {
-      Group::HotConfiguration(pRNG, U);
+      if (isSp2n == true)
+      {
+          const int nSp = Nrepresentation/2;
+          Sp<nSp>::HotConfiguration(pRNG, U);
+      } else
+      {
+          Group::HotConfiguration(pRNG, U);
+      }
   }
 
-  static inline void TepidConfiguration(GridParallelRNG &pRNG, Field &U) {
+  static inline void TepidConfiguration(GridParallelRNG &pRNG, Field &U)
-    Group::TepidConfiguration(pRNG, U);
+  {
+      if (isSp2n == true)
+    {
+        const int nSp = Nrepresentation/2;
+        Sp<nSp>::TepidConfiguration(pRNG, U);
+    } else
+    {
+        Group::TepidConfiguration(pRNG, U);
+    }
   }
 
 
@@ -212,47 +227,6 @@ public:
 
   }
 
-    
-  //sp2n... see sp2n.h
-/*
-    static inline void generate_sp2n_momenta(Field &P, GridSerialRNG & sRNG, GridParallelRNG &pRNG)
-    {
-      
-      const int nSp = Nrepresentation/2;
-      LinkField Pmu(P.Grid());
-      Pmu = Zero();
-      for (int mu = 0; mu < Nd; mu++) {
-        
-        Sp<nSp>::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu);
-        RealD scale = ::sqrt(HMC_MOMENTUM_DENOMINATOR) ;
-        Pmu = Pmu*scale;
-        PokeIndex<LorentzIndex>(P, Pmu, mu);
-      }
-    }
-    
-    static inline void update_sp2n_field(Field& P, Field& U, double ep){
- 
-      autoView(U_v,U,AcceleratorWrite);
-      autoView(P_v,P,AcceleratorRead);
-      accelerator_for(ss, P.Grid()->oSites(),1,{
-        for (int mu = 0; mu < Nd; mu++) {
-          U_v[ss](mu) = ProjectOnSpGroup(Exponentiate(P_v[ss](mu), ep, Nexp) * U_v[ss](mu));
-        }
-      });
-    }
-    
-    static inline void SpProject(Field &U) {
-      ProjectSp2n(U);
-    }
-
-
-    static inline void ColdSpConfiguration(GridParallelRNG &pRNG, Field &U) {
-      const int nSp = Nrepresentation/2;
-      Sp<nSp>::ColdConfiguration(pRNG, U);
-    }
-*/
-
-
 };
 
 

Grid/qcd/utils/Sp2n.h

@@ -310,6 +310,45 @@ public:
     }
     
     
+    template <typename LatticeMatrixType>
+    static void LieRandomize(GridParallelRNG &pRNG, LatticeMatrixType &out, double scale = 1.0)
+    {
+      GridBase *grid = out.Grid();
+
+      typedef typename LatticeMatrixType::vector_type vector_type;
+      typedef typename LatticeMatrixType::scalar_type scalar_type;
+
+      typedef iSinglet<vector_type> vTComplexType;
+
+      typedef Lattice<vTComplexType> LatticeComplexType;
+      typedef typename GridTypeMapper<typename LatticeMatrixType::vector_object>::scalar_object MatrixType;
+
+      LatticeComplexType ca(grid);
+      LatticeMatrixType lie(grid);
+      LatticeMatrixType la(grid);
+      ComplexD ci(0.0, scale);
+      //    ComplexD cone(1.0, 0.0);
+      MatrixType ta;
+
+      lie = Zero();
+
+      for (int a = 0; a < AlgebraDimension; a++) {
+        random(pRNG, ca);
+
+        ca = (ca + conjugate(ca)) * 0.5;
+        ca = ca - 0.5;
+
+        generator(a, ta);
+
+        la = ci * ca * ta;
+
+        lie = lie + la;  // e^{i la ta}
+
+      }
+      taExp(lie, out);
+    }
+    
+    
     static void GaussianFundamentalLieAlgebraMatrix(GridParallelRNG &pRNG, //same as sun
                                                     LatticeMatrix &out,
                                                     Real scale = 1.0) {
@@ -352,8 +391,8 @@ public:
  
     
     
-/*    template <typename GaugeField>
+    template <typename GaugeField>
-    static void HotSpConfiguration(GridParallelRNG &pRNG, GaugeField &out) {
+    static void HotConfiguration(GridParallelRNG &pRNG, GaugeField &out) {
       typedef typename GaugeField::vector_type vector_type;
       typedef iSp2nMatrix<vector_type> vMatrixType;
       typedef Lattice<vMatrixType> LatticeMatrixType;
@@ -365,7 +404,7 @@ public:
       }
     }
     template<typename GaugeField>
-    static void TepidSpConfiguration(GridParallelRNG &pRNG,GaugeField &out){
+    static void TepidConfiguration(GridParallelRNG &pRNG,GaugeField &out){
       typedef typename GaugeField::vector_type vector_type;
       typedef iSp2nMatrix<vector_type> vMatrixType;
       typedef Lattice<vMatrixType> LatticeMatrixType;
@@ -375,7 +414,7 @@ public:
         LieRandomize(pRNG,Umu,0.01);    //def
         PokeIndex<LorentzIndex>(out,Umu,mu);
       }
-    } */
+    }
     
     template<typename GaugeField>
     static void ColdConfiguration(GaugeField &out){
@@ -398,39 +437,12 @@ public:
     
 }; // end of class Sp
     
-    
-
-
-template<int N> //same of su
- LatticeComplexD SpDeterminant(const Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu)
- {
-   GridBase *grid=Umu.Grid();
-   auto lvol = grid->lSites();
-   LatticeComplexD ret(grid);
-
-   autoView(Umu_v,Umu,CpuRead);
-   autoView(ret_v,ret,CpuWrite);
-   thread_for(site,lvol,{
-     Eigen::MatrixXcd EigenU = Eigen::MatrixXcd::Zero(N,N);
-     Coordinate lcoor;
-     grid->LocalIndexToLocalCoor(site, lcoor);
-     iScalar<iScalar<iMatrix<ComplexD, N> > > Us;
-     peekLocalSite(Us, Umu_v, lcoor);
-     for(int i=0;i<N;i++){
-       for(int j=0;j<N;j++){
-     EigenU(i,j) = Us()()(i,j);
-       }}
-     ComplexD det = EigenU.determinant();
-     pokeLocalSite(det,ret_v,lcoor);
-   });
-   return ret;
- }
  
  template<int N>
  static void ProjectSp2n(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu)
  {
    Umu      = ProjectOnSpGroup(Umu);
-   auto det = SpDeterminant(Umu); // ok ?
+   auto det = Determinant(Umu); // ok ?
 
    det = conjugate(det);
 

configure.ac

@@ -787,7 +787,7 @@ os (target)                 : $target_os
 compiler vendor             : ${ax_cv_cxx_compiler_vendor}
 compiler version            : ${ax_cv_gxx_version}
 ----- BUILD OPTIONS -----------------------------------
-gauge group                 : `if test "${ac_ENABLE_SP}x" == "yesx"; then echo Sp2n; else echo SUn; fi`
+#gauge group                 : `if test "${ac_ENABLE_SP}x" == "yesx"; then echo Sp2n; else echo SUn; fi`
 Nc                          : ${ac_Nc}
 SIMD                        : ${ac_SIMD}${SIMD_GEN_WIDTH_MSG}
 Threading                   : ${ac_openmp}

tests/sp2n/Test_Sp_start.cc

@@ -0,0 +1,158 @@
+#include <Grid/Grid.h>
+
+using namespace Grid;
+
+int main (int argc, char **argv)
+{
+    Grid_init(&argc,&argv);
+    
+    Coordinate latt_size   = GridDefaultLatt();
+    Coordinate simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
+    Coordinate mpi_layout  = GridDefaultMpi();
+
+    GridCartesian             Grid(latt_size,simd_layout,mpi_layout);
+    GridRedBlackCartesian     RBGrid(&Grid);
+    
+    LatticeGaugeField Umu(&Grid);
+    LatticeGaugeField identity(&Grid);
+    LatticeColourMatrixD U(&Grid);
+    LatticeColourMatrixD Cidentity(&Grid);
+    LatticeColourMatrixD aux(&Grid);
+    
+    identity = 1.0;
+    Cidentity = 1.0;
+    Complex i(0., 1.);
+    
+    double vol = Umu.Grid()->gSites();
+    
+    const int nsp = Nc / 2;
+    std::vector<int> pseeds({1,2,3,4,5});
+    std::vector<int> sseeds({6,7,8,9,10});
+    GridParallelRNG  pRNG(&Grid); pRNG.SeedFixedIntegers(pseeds);
+    GridSerialRNG    sRNG;       sRNG.SeedFixedIntegers(sseeds);
+    
+    std::cout << GridLogMessage << std::endl;
+    
+    
+    Sp<nsp>::ColdConfiguration(pRNG,Umu);
+    Umu = Umu - identity;
+    U = PeekIndex<LorentzIndex>(Umu,1);
+    
+    std::cout << GridLogMessage << "U - 1 = " << norm2(sum(U)) << std::endl;
+    assert ( norm2(sum(U)) < 1e-6);
+    std::cout << GridLogMessage << std::endl;
+    
+    std::cout << GridLogMessage << "Hot Start " << std::endl;
+    Sp<nsp>::HotConfiguration(pRNG,Umu);
+    U = PeekIndex<LorentzIndex>(Umu,1);
+    
+    std::cout << GridLogMessage << "Checking unitarity " << std::endl;
+    
+    aux = U*adj(U) - Cidentity;
+
+    assert( norm2(aux) < 1e-6 );
+    
+    std::cout << GridLogMessage << "ok" << std::endl;
+    std::cout << GridLogMessage << std::endl;
+    
+    std::cout << GridLogMessage << "Checking determinant " << std::endl;
+    std::cout << GridLogMessage << "Det = " <<  norm2( Determinant(U) )/vol << std::endl;
+    std::cout << GridLogMessage << std::endl;
+    
+    std::cout << GridLogMessage << "Checking the matrix is in Sp2n " << std::endl;
+    std::cout << GridLogMessage << std::endl;
+    
+    for (int c1 = 0; c1 < nsp; c1++) //check on W
+    {
+        for (int c2 = 0; c2 < nsp; c2++)
+        {
+            auto W = PeekIndex<ColourIndex>(U,c1,c2);
+            auto Wstar =  PeekIndex<ColourIndex>(U,c1+nsp,c2+nsp);
+            auto Ww = conjugate( Wstar );
+            auto amizero = sum(W - Ww);
+            auto amizeroo = TensorRemove(amizero);
+            std::cout << GridLogMessage << "diff(real,im) = " << amizeroo << std::endl;
+            assert(  amizeroo.real() < 10e-6 );
+            amizeroo *= i;
+            assert(  amizeroo.real() < 10e-6 );
+            std::cout << GridLogMessage << "ok " << std::endl;
+        }
+        
+    }
+    std::cout <<GridLogMessage << std::endl;
+    std::cout << GridLogMessage << "Checking that U [ i , j+N/2 ] + conj U [ i+N/2, j+N/2 ] = 0 for i = 0, ... , N/2 " << std::endl;
+    std::cout <<GridLogMessage << std::endl;
+    for (int c1 = 0; c1 < nsp ; c1++)
+    {
+        for (int c2 = 0; c2 < nsp; c2++)
+        {
+            auto X = PeekIndex<ColourIndex>(U,c1,c2+nsp);
+            auto minusXstar = PeekIndex<ColourIndex>(U,c1+nsp,c2);
+            auto minusXx = conjugate(minusXstar);
+            auto amizero = sum (X + minusXx);
+            auto amizeroo = TensorRemove(amizero);
+            std::cout << GridLogMessage << "diff(real,im) = " << amizeroo << std::endl;
+            assert(  amizeroo.real() < 10e-6 );
+            amizeroo *= i;
+            assert(  amizeroo.real() < 10e-6 );
+            std::cout << GridLogMessage << "ok " << std::endl;
+        }
+    }
+    
+    
+    std::vector<int> ppseeds({11,12,13,14,15});
+    std::vector<int> ssseeds({16,17,18,19,20});
+    GridParallelRNG  ppRNG(&Grid); pRNG.SeedFixedIntegers(ppseeds);
+    GridSerialRNG    ssRNG;       sRNG.SeedFixedIntegers(ssseeds);
+    
+    std::cout << GridLogMessage << std::endl;
+    std::cout << GridLogMessage << "Testing gauge implementation " << std::endl;
+    
+    SymplPeriodicGimplR::HotConfiguration(ppRNG, Umu);
+    
+    U = PeekIndex<LorentzIndex>(Umu,1);
+    
+    for (int c1 = 0; c1 < nsp; c1++) //check on W
+    {
+        for (int c2 = 0; c2 < nsp; c2++)
+        {
+            auto W = PeekIndex<ColourIndex>(U,c1,c2);
+            auto Wstar =  PeekIndex<ColourIndex>(U,c1+nsp,c2+nsp);
+            auto Ww = conjugate( Wstar );
+            auto amizero = sum(W - Ww);
+            auto amizeroo = TensorRemove(amizero);
+            std::cout << GridLogMessage << "diff(real,im) = " << amizeroo << std::endl;
+            assert(  amizeroo.real() < 10e-6 );
+            amizeroo *= i;
+            assert(  amizeroo.real() < 10e-6 );
+            std::cout << GridLogMessage << "ok " << std::endl;
+        }
+        
+    }
+    std::cout <<GridLogMessage << std::endl;
+    std::cout << GridLogMessage << "Checking that U [ i , j+N/2 ] + conj U [ i+N/2, j+N/2 ] = 0 for i = 0, ... , N/2 " << std::endl;
+    std::cout <<GridLogMessage << std::endl;
+    for (int c1 = 0; c1 < nsp ; c1++)
+    {
+        for (int c2 = 0; c2 < nsp; c2++)
+        {
+            auto X = PeekIndex<ColourIndex>(U,c1,c2+nsp);
+            auto minusXstar = PeekIndex<ColourIndex>(U,c1+nsp,c2);
+            auto minusXx = conjugate(minusXstar);
+            auto amizero = sum (X + minusXx);
+            auto amizeroo = TensorRemove(amizero);
+            std::cout << GridLogMessage << "diff(real,im) = " << amizeroo << std::endl;
+            assert(  amizeroo.real() < 10e-6 );
+            amizeroo *= i;
+            assert(  amizeroo.real() < 10e-6 );
+            std::cout << GridLogMessage << "ok " << std::endl;
+        }
+    }
+    
+    
+    
+    Grid_finalize();
+
+
+    }
+