Checked the hermiticity of the op in derivative, ok

Still CG fails to converge

lib/Make.inc

@@ -1,4 +1,4 @@
 
-HFILES=./tensors/Tensor_extract_merge.h ./tensors/Tensor_exp.h ./tensors/Tensor_index.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_determinant.h ./tensors/Tensor_outer.h ./tensors/Tensor_trace.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_logical.h ./tensors/Tensor_reality.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_Ta.h ./tensors/Tensor_class.h ./tensors/Tensor_inner.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_transpose.h ./tensors/Tensor_arith.h ./tensors/Tensor_unary.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_traits.h ./stencil/Lebesgue.h ./PerfCount.h ./Communicator.h ./simd/Grid_vector_types.h ./simd/Grid_qpx.h ./simd/Grid_imci.h ./simd/Grid_avx.h ./simd/Grid_empty.h ./simd/Intel512single.h ./simd/Grid_vector_unops.h ./simd/Grid_neon.h ./simd/Intel512common.h ./simd/Intel512wilson.h ./simd/Intel512imci.h ./simd/Intel512double.h ./simd/Grid_avx512.h ./simd/Grid_sse4.h ./simd/Intel512avx.h ./Threads.h ./Timer.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/bigfloat.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/SparseMatrix.h ./algorithms/LinearOperator.h ./algorithms/Preconditioner.h ./algorithms/CoarsenedMatrix.h ./algorithms/iterative/AdefGeneric.h ./algorithms/iterative/ImplicitlyRestartedLanczos.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/DenseMatrix.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/EigenSort.h ./algorithms/iterative/MatrixUtils.h ./algorithms/iterative/ConjugateGradientMixedPrec.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/Francis.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/Householder.h ./algorithms/iterative/PrecConjugateResidual.h ./algorithms/iterative/Matrix.h ./Algorithms.h ./lattice/Lattice_base.h ./lattice/Lattice_overload.h ./lattice/Lattice_ET.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_trace.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_arith.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_local.h ./lattice/Lattice_rng.h ./lattice/Lattice_comparison.h ./lattice/Lattice_reality.h ./lattice/Lattice_transpose.h ./lattice/Lattice_transfer.h ./lattice/Lattice_conformable.h ./lattice/Lattice_unary.h ./lattice/Lattice_where.h ./lattice/Lattice_reduction.h ./Cartesian.h ./Simd.h ./Init.h ./serialisation/TextIO.h ./serialisation/Serialisation.h ./serialisation/MacroMagic.h ./serialisation/BaseIO.h ./serialisation/XmlIO.h ./serialisation/BinaryIO.h ./communicator/Communicator_base.h ./Log.h ./Lattice.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./cshift/Cshift_common.h ./Lexicographic.h ./Tensors.h ./Cshift.h ./parallelIO/NerscIO.h ./parallelIO/BinaryIO.h ./qcd/smearing/BaseSmearing.h ./qcd/smearing/GaugeConfiguration.h ./qcd/smearing/StoutSmearing.h ./qcd/smearing/APEsmearing.h ./qcd/smearing/Smearing.h ./qcd/hmc/HmcRunner.h ./qcd/hmc/NerscCheckpointer.h ./qcd/hmc/integrators/Integrator.h ./qcd/hmc/integrators/Integrator_algorithm.h ./qcd/hmc/HMC.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/FermionOperatorImpl.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/WilsonTMFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonKernelsAsmBody.h ./qcd/action/ActionParams.h ./qcd/action/Actions.h ./qcd/action/pseudofermion/OneFlavourRationalRatio.h ./qcd/action/pseudofermion/OneFlavourEvenOddRational.h ./qcd/action/pseudofermion/TwoFlavourEvenOdd.h ./qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h ./qcd/action/pseudofermion/TwoFlavourRatio.h ./qcd/action/pseudofermion/EvenOddSchurDifferentiable.h ./qcd/action/pseudofermion/TwoFlavour.h ./qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h ./qcd/action/pseudofermion/OneFlavourRational.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/action/gauge/PlaqPlusRectangleAction.h ./qcd/action/gauge/GaugeImpl.h ./qcd/action/ActionBase.h ./qcd/QCD.h ./qcd/spin/TwoSpinor.h ./qcd/spin/Dirac.h ./qcd/representations/hmc_types.h ./qcd/representatiLine truncated
+HFILES=./tensors/Tensor_extract_merge.h ./tensors/Tensor_exp.h ./tensors/Tensor_index.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_determinant.h ./tensors/Tensor_outer.h ./tensors/Tensor_trace.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_logical.h ./tensors/Tensor_reality.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_Ta.h ./tensors/Tensor_class.h ./tensors/Tensor_inner.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_transpose.h ./tensors/Tensor_arith.h ./tensors/Tensor_unary.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_traits.h ./stencil/Lebesgue.h ./PerfCount.h ./Communicator.h ./simd/Grid_vector_types.h ./simd/Grid_qpx.h ./simd/Grid_imci.h ./simd/Grid_avx.h ./simd/Grid_empty.h ./simd/Intel512single.h ./simd/Grid_vector_unops.h ./simd/Grid_neon.h ./simd/Intel512common.h ./simd/Intel512wilson.h ./simd/Intel512imci.h ./simd/Intel512double.h ./simd/Grid_avx512.h ./simd/Grid_sse4.h ./simd/Intel512avx.h ./Threads.h ./Timer.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/bigfloat.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/SparseMatrix.h ./algorithms/LinearOperator.h ./algorithms/Preconditioner.h ./algorithms/CoarsenedMatrix.h ./algorithms/iterative/AdefGeneric.h ./algorithms/iterative/ImplicitlyRestartedLanczos.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/DenseMatrix.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/EigenSort.h ./algorithms/iterative/MatrixUtils.h ./algorithms/iterative/ConjugateGradientMixedPrec.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/Francis.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/Householder.h ./algorithms/iterative/PrecConjugateResidual.h ./algorithms/iterative/Matrix.h ./Algorithms.h ./lattice/Lattice_base.h ./lattice/Lattice_overload.h ./lattice/Lattice_ET.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_trace.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_arith.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_local.h ./lattice/Lattice_rng.h ./lattice/Lattice_comparison.h ./lattice/Lattice_reality.h ./lattice/Lattice_transpose.h ./lattice/Lattice_transfer.h ./lattice/Lattice_conformable.h ./lattice/Lattice_unary.h ./lattice/Lattice_where.h ./lattice/Lattice_reduction.h ./Cartesian.h ./Simd.h ./Init.h ./serialisation/TextIO.h ./serialisation/Serialisation.h ./serialisation/MacroMagic.h ./serialisation/BaseIO.h ./serialisation/XmlIO.h ./serialisation/BinaryIO.h ./communicator/Communicator_base.h ./Log.h ./Lattice.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./cshift/Cshift_common.h ./Lexicographic.h ./Tensors.h ./Cshift.h ./parallelIO/NerscIO.h ./parallelIO/BinaryIO.h ./qcd/smearing/BaseSmearing.h ./qcd/smearing/GaugeConfiguration.h ./qcd/smearing/StoutSmearing.h ./qcd/smearing/APEsmearing.h ./qcd/smearing/Smearing.h ./qcd/hmc/HmcRunner.h ./qcd/hmc/NerscCheckpointer.h ./qcd/hmc/integrators/Integrator.h ./qcd/hmc/integrators/Integrator_algorithm.h ./qcd/hmc/HMC.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/FermionOperatorImpl.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/WilsonTMFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonKernelsAsmBody.h ./qcd/action/ActionParams.h ./qcd/action/Actions.h ./qcd/action/pseudofermion/OneFlavourRationalRatio.h ./qcd/action/pseudofermion/OneFlavourEvenOddRational.h ./qcd/action/pseudofermion/TwoFlavourEvenOdd.h ./qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h ./qcd/action/pseudofermion/TwoFlavourRatio.h ./qcd/action/pseudofermion/EvenOddSchurDifferentiable.h ./qcd/action/pseudofermion/TwoFlavour.h ./qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h ./qcd/action/pseudofermion/OneFlavourRational.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/action/gauge/PlaqPlusRectangleAction.h ./qcd/action/gauge/GaugeImpl.h ./qcd/action/ActionBase.h ./qcd/QCD.h ./qcd/spin/TwoSpinor.h ./qcd/spin/Dirac.h ./qcd/representations/fundamental.h ./qcd/representaLine truncated
 
 CCFILES=./stencil/Stencil_common.cc ./stencil/Lebesgue.cc ./algorithms/approx/MultiShiftFunction.cc ./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./PerfCount.cc ./serialisation/BinaryIO.cc ./serialisation/XmlIO.cc ./serialisation/TextIO.cc ./qcd/hmc/HMC.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/action/fermion/WilsonTMFermion.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/WilsonKernelsAsm.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/spin/Dirac.cc ./qcd/utils/SpaceTimeGrid.cc ./Init.cc ./pugixml/pugixml.cc ./Log.cc

lib/qcd/action/pseudofermion/TwoFlavour.h

@@ -134,9 +134,9 @@ class TwoFlavourPseudoFermionAction : public Action<typename Impl::GaugeField> {
     MdagMOp.Op(X, Y);
 
     // Check hermiticity
-    /*
+    
     std::vector<int> seeds({1,2,3,4});
-    GridParallelRNG          RNG(U._grid);   RNG.SeedFixedIntegers(seeds);
+    GridParallelRNG RNG(U._grid);   RNG.SeedFixedIntegers(seeds);
     FermionField RNGphi(FermOp.FermionGrid());
     FermionField RNGchi(FermOp.FermionGrid());
     FermionField Achi(FermOp.FermionGrid());
@@ -144,8 +144,8 @@ class TwoFlavourPseudoFermionAction : public Action<typename Impl::GaugeField> {
 
     random(RNG, RNGphi);
     random(RNG, RNGchi);
-    MdagMOp.Op(RNGchi, Achi);
-    MdagMOp.Op(RNGphi, Aphi);
+    MdagMOp.HermOp(RNGchi, Achi);
+    MdagMOp.HermOp(RNGphi, Aphi);
     ComplexD pAc = innerProduct(RNGphi, Achi);
     ComplexD cAp = innerProduct(RNGchi, Aphi);
     //these should be real
@@ -153,8 +153,9 @@ class TwoFlavourPseudoFermionAction : public Action<typename Impl::GaugeField> {
     ComplexD cAc = innerProduct(RNGchi, Achi);
 
     std::cout<<GridLogMessage<< "pAc "<<pAc<<" cAp "<< cAp<< " diff "<<pAc-adj(cAp)<<std::endl;
+    // These ones should be real
     std::cout << GridLogMessage << "pAp " << pAp << " cAc " << cAc << std::endl;
-*/
+
     // Our conventions really make this UdSdU; We do not differentiate wrt Udag
     // here.
     // So must take dSdU - adj(dSdU) and left multiply by mom to get dS/dt.

lib/qcd/hmc/HmcRunner.h

@@ -109,12 +109,12 @@ class NerscHmcRunnerTemplate {
     int Nsmear = 1;    // number of smearing levels
     Smear_Stout<Gimpl> Stout(rho);
     std::cout << GridLogDebug << " Creating the SmearedConfiguration class\n";
-    SmearedConfiguration<Gimpl> SmearingPolicy(UGrid, Nsmear, Stout);
+    //SmearedConfiguration<Gimpl> SmearingPolicy(UGrid, Nsmear, Stout);
     std::cout << GridLogDebug << " done\n";
     
     //////////////
-    //NoSmearing<Gimpl> SmearingPolicy;
-    typedef MinimumNorm2<GaugeField, SmearedConfiguration<Gimpl>, RepresentationsPolicy >
+    NoSmearing<Gimpl> SmearingPolicy;
+    typedef MinimumNorm2<GaugeField, NoSmearing<Gimpl>, RepresentationsPolicy >
         IntegratorType;  // change here to change the algorithm
     IntegratorParameters MDpar(20, 1.0);
     IntegratorType MDynamics(UGrid, MDpar, TheAction, SmearingPolicy);

lib/qcd/hmc/integrators/Integrator.h

@@ -149,7 +149,7 @@ class Integrator {
     }
 
     // Force from the other representations
-    //as[level].apply(update_P_hireps, Representations, Mom, U, ep);
+    as[level].apply(update_P_hireps, Representations, Mom, U, ep);
   }
 
   void update_U(GaugeField& U, double ep) {
@@ -167,13 +167,12 @@ class Integrator {
       auto Umu = PeekIndex<LorentzIndex>(U, mu);
       auto Pmu = PeekIndex<LorentzIndex>(Mom, mu);
       Umu = expMat(Pmu, ep, Params.Nexp) * Umu;
-      ProjectOnGroup(Umu);
-      PokeIndex<LorentzIndex>(U, Umu, mu);
+      PokeIndex<LorentzIndex>(U, ProjectOnGroup(Umu), mu);
     }
     // Update the smeared fields, can be implemented as observer
     Smearer.set_GaugeField(U);
     // Update the higher representations fields
-    //Representations.update(U);  // void functions if fundamental representation
+    Representations.update(U);  // void functions if fundamental representation
   }
 
   virtual void step(GaugeField& U, int level, int first, int last) = 0;
@@ -217,7 +216,7 @@ class Integrator {
     // of the Metropolis
     Smearer.set_GaugeField(U);
     // Set the (eventual) representations gauge fields
-    //Representations.update(U);
+    Representations.update(U);
 
     // The Smearer is attached to a pointer of the gauge field
     // automatically gets the correct field
@@ -232,7 +231,7 @@ class Integrator {
       }
 
       // Refresh the higher representation actions
-      //as[level].apply(refresh_hireps, Representations, pRNG);
+      as[level].apply(refresh_hireps, Representations, pRNG);
     }
   }
 
@@ -281,7 +280,7 @@ class Integrator {
                   << actionID << " H = " << Hterm << std::endl;
         H += Hterm;
       }
-      //as[level].apply(S_hireps, Representations, level, H);
+      as[level].apply(S_hireps, Representations, level, H);
     }
 
     return H;

lib/qcd/representations/adjoint.h

@@ -1,6 +1,6 @@
 /*
- *	Policy classes for the HMC
- *	Author: Guido Cossu
+ *  Policy classes for the HMC
+ *  Author: Guido Cossu
 */
 
 #ifndef ADJOINT_H
@@ -27,8 +27,9 @@ class AdjointRep {
   LatticeField U;
 
   explicit AdjointRep(GridBase *grid) : U(grid) {}
+
   void update_representation(const LatticeGaugeField &Uin) {
-  	std::cout << GridLogDebug << "Updating adjoint representation\n" ;
+    std::cout << GridLogDebug << "Updating adjoint representation\n" ;
     // Uin is in the fundamental representation
     // get the U in AdjointRep
     // (U_adj)_B = tr[e^a U e^b U^dag]
@@ -41,28 +42,36 @@ class AdjointRep {
 
     Vector<typename SU<ncolour>::Matrix> ta(Dimension);
 
+    // Debug lines
+    //LatticeMatrix uno(Uin._grid);
+    //uno = 1.0;
+    ////////////////
+
     // FIXME probably not very efficient to get all the generators
     // everytime
     for (int a = 0; a < Dimension; a++) SU<ncolour>::generator(a, ta[a]);
 
     for (int mu = 0; mu < Nd; mu++) {
       auto Uin_mu = peekLorentz(Uin, mu);
-      auto U_mu = peekLorentz(U, mu);
+      auto U_mu   = peekLorentz(U, mu);
       for (int a = 0; a < Dimension; a++) {
         tmp = 2.0 * adj(Uin_mu) * ta[a] * Uin_mu;
         for (int b = 0; b < Dimension; b++)
-          pokeColour(U_mu, trace(tmp * ta[b]), a, b);
+          pokeColour(U_mu, trace(tmp * ta[b]), b, a);
       }
-      // Check matrix U_mu, must be real orthogonal
-      //reality
-      LatticeMatrix Ucheck = U_mu - conjugate(U_mu);
-      std::cout << GridLogMessage << "Reality check: " << norm2(Ucheck) << std::endl;
-      LatticeMatrix uno(Uin._grid);
-      uno = 1.0;
-      Ucheck = U_mu * adj(U_mu) - uno;
-      std::cout << GridLogMessage << "orthogonality check: " << norm2(Ucheck) << std::endl;
-
       pokeLorentz(U, U_mu, mu);
+      // Check matrix U_mu, must be real orthogonal
+      // reality
+      /*
+      LatticeMatrix Ucheck = U_mu - conjugate(U_mu);
+      std::cout << GridLogMessage << "Reality check: " << norm2(Ucheck) <<
+      std::endl;
+
+      Ucheck = U_mu * adj(U_mu) - uno;
+      std::cout << GridLogMessage << "orthogonality check: " << norm2(Ucheck) <<
+      std::endl;
+      */
+
     }
   }
 
@@ -81,8 +90,8 @@ class AdjointRep {
       projectOnAlgebra(h, in_mu, scale);
       FundamentalLieAlgebraMatrix(h, out_mu, 1.0);  // apply scale only once
       pokeLorentz(out, out_mu, mu);
-	  // Returns traceless antihermitian matrix Nc * Nc.
-	  // Confirmed            
+    // Returns traceless antihermitian matrix Nc * Nc.
+    // Confirmed            
     }
     return out;
   }