Adding Binrary IO, untested

lib/algorithms/approx/Chebyshev.h

@@ -197,8 +197,9 @@ namespace Grid {
     void operator() (LinearOperatorBase<Field> &Linop, const Field &in, Field &out) {
 
       GridBase *grid=in._grid;
-//std::cout << "Chevyshef(): in._grid="<<in._grid<<std::endl;
+
-//<<" Linop.Grid()="<<Linop.Grid()<<"Linop.RedBlackGrid()="<<Linop.RedBlackGrid()<<std::endl;
+      // std::cout << "Chevyshef(): in._grid="<<in._grid<<std::endl;
+      //std::cout <<" Linop.Grid()="<<Linop.Grid()<<"Linop.RedBlackGrid()="<<Linop.RedBlackGrid()<<std::endl;
 
       int vol=grid->gSites();
 

lib/parallelIO/NerscIO.h

@@ -397,8 +397,8 @@ static inline void writeConfiguration(Lattice<iLorentzColourMatrix<vsimd> > &Umu
 
   typedef LorentzColourMatrixD fobj3D;
   typedef LorentzColour2x3D    fobj2D;
-  typedef LorentzColourMatrixF fobj3f;
+  //typedef LorentzColourMatrixF fobj3f;
-  typedef LorentzColour2x3F    fobj2f;
+  //typedef LorentzColour2x3F    fobj2f;
 
   GridBase *grid = Umu._grid;
 

lib/qcd/hmc/BinaryCheckpointer.h

@@ -1,106 +1,146 @@
-    /*************************************************************************************
+/*************************************************************************************
 
-    Grid physics library, www.github.com/paboyle/Grid 
+Grid physics library, www.github.com/paboyle/Grid
 
-    Source file: ./lib/qcd/hmc/NerscCheckpointer.h
+Source file: ./lib/qcd/hmc/NerscCheckpointer.h
 
-    Copyright (C) 2015
+Copyright (C) 2015
 
 Author: paboyle <paboyle@ph.ed.ac.uk>
 
-    This program is free software; you can redistribute it and/or modify
+This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
+it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
+the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
+(at your option) any later version.
 
-    This program is distributed in the hope that it will be useful,
+This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
+but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
+GNU General Public License for more details.
 
-    You should have received a copy of the GNU General Public License along
+You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
+with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 
-    See the full license in the file "LICENSE" in the top level distribution directory
+See the full license in the file "LICENSE" in the top level distribution
-    *************************************************************************************/
+directory
-    /*  END LEGAL */
+*************************************************************************************/
-#ifndef NERSC_CHECKPOINTER
+/*  END LEGAL */
-#define NERSC_CHECKPOINTER
+#ifndef BINARY_CHECKPOINTER
+#define BINARY_CHECKPOINTER
 
-#include <string>
 #include <iostream>
 #include <sstream>
+#include <string>
 
+namespace Grid {
+namespace QCD {
 
-namespace Grid{
+template <class fobj, class sobj>
-  namespace QCD{
+struct BinarySimpleUnmunger {
+  typedef typename getPrecision<fobj>::real_scalar_type fobj_stype;
+  typedef typename getPrecision<sobj>::real_scalar_type sobj_stype;
+
+  void operator()(sobj &in, fobj &out, uint32_t &csum) {
+    // take word by word and transform accoding to the status
+  	fobj_stype* out_buffer = (fobj_stype*)&out;
+  	sobj_stype* in_buffer = (sobj_stype*)&in;
+  	size_t fobj_words = sizeof(out)/sizeof(fobj_stype);
+  	size_t sobj_words = sizeof(in)/sizeof(sobj_stype);
+  	assert(fobj_words == sobj_words);
+
+  	for (unsigned int word = 0; word < sobj_words; word++)
+			out_buffer[word] = in_buffer[word];  // type conversion on the fly		
+
+			BinaryIO::Uint32Checksum((uint32_t*)&out,sizeof(out),csum);
     
-    template<class fobj,class sobj>
+  };
-    struct BinarySimpleUnmunger{
+
-      void operator() (sobj &in,fobj &out,uint32_t &csum){
+template <class fobj, class sobj>
-        fobj = sobj;
+struct BinarySimpleMunger {
-        csum =0;
+  typedef typename getPrecision<fobj>::real_scalar_type fobj_stype;
-    };
+  typedef typename getPrecision<sobj>::real_scalar_type sobj_stype;
+
+  void operator()(sobj &out, fobj &in, uint32_t &csum) {
+    // take word by word and transform accoding to the status
+  	fobj_stype* in_buffer = (fobj_stype*)&in;
+  	sobj_stype* out_buffer = (sobj_stype*)&out;
+  	size_t fobj_words = sizeof(in)/sizeof(fobj_stype);
+  	size_t sobj_words = sizeof(out)/sizeof(sobj_stype);
+  	assert(fobj_words == sobj_words);
+
+  	for (unsigned int word = 0; word < sobj_words; word++)
+			out_buffer[word] = in_buffer[word];  // type conversion on the fly		
+
+			BinaryIO::Uint32Checksum((uint32_t*)&in,sizeof(in),csum);
+    
+  };
 
 
-    // Only for gauge fields
+  // Only for the main field in the hmc
-    template<class Impl> 
+  template <class Impl>
-    class BinaryHmcCheckpointer : public HmcObservable<typename Impl::Field> {
+  class BinaryHmcCheckpointer : public HmcObservable<typename Impl::Field> {
-    private:
+   private:
-      std::string configStem;
+    std::string configStem;
-      std::string rngStem;
+    std::string rngStem;
-      int SaveInterval;
+    int SaveInterval;
-    public:
-      INHERIT_FIELD_TYPES(Impl); // The Field is a Lattice object
 
-      typedef typename Field::vector_object vobj;
+   public:
-      typedef typename vobj::scalar_object sobj;
+    INHERIT_FIELD_TYPES(Impl);  // The Field is a Lattice object
- 
 
-      BinaryHmcCheckpointer(std::string cf, std::string rn,int savemodulo) {
+    typedef typename Field::vector_object vobj;
-        configStem  = cf;
+    typedef typename vobj::scalar_object sobj;
-        rngStem     = rn;
+		typedef typename getPrecision<sobj>::real_scalar_type sobj_stype;
-        SaveInterval= savemodulo;
+		typedef typename sobj::DoublePrecision sobj_double;
-      };
 
-      void TrajectoryComplete(int traj, Field &U, GridSerialRNG &sRNG,
+    BinaryHmcCheckpointer(std::string cf, std::string rn, int savemodulo, const std::string &format)
-                              GridParallelRNG &pRNG) {
+        : configStem(cf),
-        if ((traj % SaveInterval) == 0) {
+          rngStem(rn),
-          std::string rng;
+          SaveInterval(savemodulo){};
-          {
-            std::ostringstream os;
-            os << rngStem << "." << traj;
-            rng = os.str();
-          }
-          std::string config;
-          {
-            std::ostringstream os;
-            os << configStem << "." << traj;
-            config = os.str();
-          }
 
-          int prec = getPrecision<Field>::value; //get precision of oject
+    void TrajectoryComplete(int traj, Field &U, GridSerialRNG &sRNG,
-          std::string floating_point = "IEEE64BIG"; //default precision
+                            GridParallelRNG &pRNG) {
-          if (prec == 1)
+      if ((traj % SaveInterval) == 0) {
-            floating_point = "IEEE32BIG"
+        std::string rng;
-
+        {
-          BinarySimpleUnmunger<sobj, sobj> munge;
+          std::ostringstream os;
-          BinaryIO::writeRNGSerial(sRNG, pRNG, rng, 0);
+          os << rngStem << "." << traj;
-          BinaryIO::writeObjectParallel<vobj, sobj>(U, config, munge, 0 ,floating_point);
+          rng = os.str();
+        }
+        std::string config;
+        {
+          std::ostringstream os;
+          os << configStem << "." << traj;
+          config = os.str();
         }
-      };
-
-      void CheckpointRestore(int traj, GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG & pRNG ){
-
-  std::string rng;   { std::ostringstream os; os << rngStem     <<"."<< traj; rng = os.str(); }
-  std::string config;{ std::ostringstream os; os << configStem  <<"."<< traj; config = os.str();}
-
-  NerscField header;
-  NerscIO::readRNGState(sRNG,pRNG,header,rng);
-  NerscIO::readConfiguration(U,header,config);
-      };
 
+        // Save always in double precision
+        BinarySimpleUnmunger<sobj_double, sobj> munge;
+        BinaryIO::writeRNGSerial(sRNG, pRNG, rng, 0);
+        BinaryIO::writeObjectParallel<vobj, sobj_double>(U, config, munge, 0, format);
+      }
     };
-}}
+
+    void CheckpointRestore(int traj, Field &U, GridSerialRNG &sRNG,
+                           GridParallelRNG &pRNG) {
+      std::string rng;
+      {
+        std::ostringstream os;
+        os << rngStem << "." << traj;
+        rng = os.str();
+      }
+      std::string config;
+      {
+        std::ostringstream os;
+        os << configStem << "." << traj;
+        config = os.str();
+      }
+
+      BinarySimpleMunger<sobj_double, sobj> munge;
+      BinaryIO::readRNGSerial(sRNG, pRNG, rng, header);
+      BinaryIO::readObjectParallel<vobj, sobj_double>(U, config, munge, 0, format);
+    };
+  };
+}
+}
 #endif

lib/tensors/Tensor_class.h

@@ -65,6 +65,9 @@ class iScalar {
   typedef iScalar<typename GridTypeMapper<vtype>::Complexified> Complexified;
   typedef iScalar<typename GridTypeMapper<vtype>::Realified> Realified;
 
+  // get double precision version
+  typedef iScalar<typename GridTypeMapper<vtype>::DoublePrecision> DoublePrecision;
+
   enum { TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1 };
 
   // Scalar no action
@@ -197,6 +200,10 @@ class iVector {
   typedef iVector<typename GridTypeMapper<vtype>::Complexified, N> Complexified;
   typedef iVector<typename GridTypeMapper<vtype>::Realified, N> Realified;
 
+  // get double precision version
+  typedef iVector<typename GridTypeMapper<vtype>::DoublePrecision, N> DoublePrecision;
+
+
   template <class T, typename std::enable_if<!isGridTensor<T>::value, T>::type
                          * = nullptr>
   strong_inline auto operator=(T arg) -> iVector<vtype, N> {
@@ -300,7 +307,11 @@ class iMatrix {
   typedef iMatrix<typename GridTypeMapper<vtype>::Complexified, N> Complexified;
   typedef iMatrix<typename GridTypeMapper<vtype>::Realified, N> Realified;
 
-  // Tensure removal
+  // get double precision version
+  typedef iMatrix<typename GridTypeMapper<vtype>::DoublePrecision, N> DoublePrecision;
+
+
+  // Tensor removal
   typedef iScalar<tensor_reduced_v> tensor_reduced;
   typedef iMatrix<recurse_scalar_object, N> scalar_object;
 

lib/tensors/Tensor_traits.h

@@ -57,6 +57,7 @@ namespace Grid {
     typedef typename T::scalar_object scalar_object;
     typedef typename T::Complexified Complexified;
     typedef typename T::Realified Realified;
+    typedef typename T::DoublePrecision DoublePrecision;
     enum { TensorLevel = T::TensorLevel };
   };
 
@@ -71,6 +72,7 @@ namespace Grid {
     typedef RealF scalar_object;
     typedef ComplexF Complexified;
     typedef RealF Realified;
+    typedef RealD DoublePrecision;
     enum { TensorLevel = 0 };
   };
   template<> class GridTypeMapper<RealD> {
@@ -81,6 +83,7 @@ namespace Grid {
     typedef RealD scalar_object;
     typedef ComplexD Complexified;
     typedef RealD Realified;
+    typedef RealD DoublePrecision;
     enum { TensorLevel = 0 };
   };
   template<> class GridTypeMapper<ComplexF> {
@@ -91,6 +94,7 @@ namespace Grid {
     typedef ComplexF scalar_object;
     typedef ComplexF Complexified;
     typedef RealF Realified;
+    typedef ComplexD DoublePrecision;
     enum { TensorLevel = 0 };
   };
   template<> class GridTypeMapper<ComplexD> {
@@ -101,6 +105,7 @@ namespace Grid {
     typedef ComplexD scalar_object;
     typedef ComplexD Complexified;
     typedef RealD Realified;
+    typedef ComplexD DoublePrecision;
     enum { TensorLevel = 0 };
   };
   template<> class GridTypeMapper<Integer> {
@@ -111,6 +116,7 @@ namespace Grid {
     typedef Integer scalar_object;
     typedef void Complexified;
     typedef void Realified;
+    typedef void DoublePrecision;
     enum { TensorLevel = 0 };
   };
 
@@ -122,6 +128,7 @@ namespace Grid {
     typedef RealF  scalar_object;
     typedef vComplexF Complexified;
     typedef vRealF Realified;
+    typedef vRealD DoublePrecision;
     enum { TensorLevel = 0 };
   };
   template<> class GridTypeMapper<vRealD> {
@@ -132,6 +139,7 @@ namespace Grid {
     typedef RealD  scalar_object;
     typedef vComplexD Complexified;
     typedef vRealD Realified;
+    typedef vRealD DoublePrecision;
     enum { TensorLevel = 0 };
   };
   template<> class GridTypeMapper<vComplexF> {
@@ -142,6 +150,7 @@ namespace Grid {
     typedef ComplexF  scalar_object;
     typedef vComplexF Complexified;
     typedef vRealF Realified;
+    typedef vComplexD DoublePrecision;
     enum { TensorLevel = 0 };
   };
   template<> class GridTypeMapper<vComplexD> {
@@ -152,6 +161,7 @@ namespace Grid {
     typedef ComplexD  scalar_object;
     typedef vComplexD Complexified;
     typedef vRealD Realified;
+    typedef vComplexD DoublePrecision;
     enum { TensorLevel = 0 };
   };
   template<> class GridTypeMapper<vInteger> {
@@ -162,6 +172,7 @@ namespace Grid {
     typedef  Integer scalar_object;
     typedef void Complexified;
     typedef void Realified;
+    typedef void DoublePrecision;
     enum { TensorLevel = 0 };
   };
 
@@ -256,8 +267,8 @@ namespace Grid {
     typedef typename getVectorType<T>::type vector_obj; //get the vector_obj (i.e. a grid Tensor) if its a Lattice<vobj>, do nothing otherwise (i.e. if fundamental or grid Tensor)
   
     typedef typename GridTypeMapper<vector_obj>::scalar_type scalar_type; //get the associated scalar type. Works on fundamental and tensor types
-    typedef typename GridTypeMapper<scalar_type>::Realified real_scalar_type; //remove any std::complex wrapper, should get us to the fundamental type
   public:
+    typedef typename GridTypeMapper<scalar_type>::Realified real_scalar_type; //remove any std::complex wrapper, should get us to the fundamental type
     enum { value = sizeof(real_scalar_type)/sizeof(float) };
   };
 }

tests/solver/Test_wilson_lanczos.cc

@@ -0,0 +1,108 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./tests/Test_dwf_lanczos.cc
+
+Copyright (C) 2015
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Grid.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+typedef WilsonFermionR FermionOp;
+typedef typename WilsonFermionR::FermionField FermionField;
+
+
+RealD AllZero(RealD x) { return 0.; }
+
+int main(int argc, char** argv) {
+  Grid_init(&argc, &argv);
+
+  GridCartesian* UGrid = SpaceTimeGrid::makeFourDimGrid(
+      GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()),
+      GridDefaultMpi());
+  GridRedBlackCartesian* UrbGrid =
+      SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  GridCartesian* FGrid = UGrid;
+  GridRedBlackCartesian* FrbGrid = UrbGrid;
+  printf("UGrid=%p UrbGrid=%p FGrid=%p FrbGrid=%p\n", UGrid, UrbGrid, FGrid,
+         FrbGrid);
+
+  std::vector<int> seeds4({1, 2, 3, 4});
+  std::vector<int> seeds5({5, 6, 7, 8});
+  GridParallelRNG RNG5(FGrid);
+  RNG5.SeedFixedIntegers(seeds5);
+  GridParallelRNG RNG4(UGrid);
+  RNG4.SeedFixedIntegers(seeds4);
+  GridParallelRNG RNG5rb(FrbGrid);
+  RNG5.SeedFixedIntegers(seeds5);
+
+  LatticeGaugeField Umu(UGrid);
+  SU3::HotConfiguration(RNG4, Umu);
+
+/*
+  std::vector<LatticeColourMatrix> U(4, UGrid);
+  for (int mu = 0; mu < Nd; mu++) {
+    U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
+  }
+*/
+
+  RealD mass = -0.1;
+  RealD M5 = 1.8;
+  RealD mob_b = 1.5;
+  FermionOp WilsonOperator(Umu,*FGrid,*FrbGrid,mass);
+  MdagMLinearOperator<FermionOp,LatticeFermion> HermOp(WilsonOperator); /// <-----
+  //SchurDiagTwoOperator<FermionOp,FermionField> HermOp(WilsonOperator);
+
+  const int Nstop = 20;
+  const int Nk = 60;
+  const int Np = 60;
+  const int Nm = Nk + Np;
+  const int MaxIt = 10000;
+  RealD resid = 1.0e-6;
+
+  std::vector<double> Coeffs{0, 1.};
+  Polynomial<FermionField> PolyX(Coeffs);
+  Chebyshev<FermionField> Cheb(0.0, 10., 12);
+  ImplicitlyRestartedLanczos<FermionField> IRL(HermOp, PolyX, Nstop, Nk, Nm,
+                                               resid, MaxIt);
+
+  std::vector<RealD> eval(Nm);
+  FermionField src(FGrid);
+  gaussian(RNG5, src);
+  std::vector<FermionField> evec(Nm, FGrid);
+  for (int i = 0; i < 1; i++) {
+    std::cout << i << " / " << Nm << " grid pointer " << evec[i]._grid
+              << std::endl;
+  };
+
+  int Nconv;
+  IRL.calc(eval, evec, src, Nconv);
+
+  std::cout << eval << std::endl;
+
+  Grid_finalize();
+}