Added two spinor functionality required to support the Wilson hopping term.

TODO

@@ -13,7 +13,26 @@
 * Consider switch std::vector to boost arrays or something lighter weight
   boost::multi_array<type, 3> A()...    to replace multi1d, multi2d etc..
 
-* How to define simple matrix operations, such as flavour matrices?
+* TensorRemove is a hack, come up with a long term rationalised approach to Complex vs. Scalar<Scalar<Scalar<Complex > > >
+  QDP forces use of "toDouble" to get back to non tensor scalar. This role is presently taken TensorRemove, but I
+  want to introduce a syntax that does not require this.
+
+* norm2l is a hack. figure out syntax error and make this norm2. c.f. tests/Grid_gamma.cc
+
+*  std::vector replacement;
+
+  Had to change "reserve" back to "resize" on std::vector in Lattice class.
+
+  This forces the constructor call on EVERY element of the array with negative
+  performance effects on temporaries.
+  
+  The reversion was required because copy constructur has to work.
+
+  CONCLUSION: I must implement a similar to vector without construction/fill on
+  resize. Find out if valarray or alternative works differently prior to 
+  doing this since there may still be something I can use..
+
+* Flavour matrices?
 
 * Make the Tensor types and Complex etc... play more nicely.
 
@@ -36,7 +55,6 @@
 
 * Conformable test in Cshift routines.
 
-* QDP++ regression suite and comparative benchmark
 
 AUDITS:
 

lib/math/Grid_math_inner.h

@@ -6,7 +6,7 @@ namespace Grid {
     // innerProduct Vector x Vector -> Scalar
     // innerProduct Matrix x Matrix -> Scalar
     ///////////////////////////////////////////////////////////////////////////////////////
-    template<class sobj> inline RealD norm2l(sobj &arg){
+    template<class sobj> inline RealD norm2l(const sobj &arg){
       typedef typename sobj::scalar_type scalar;
       decltype(innerProduct(arg,arg)) nrm;
       nrm = innerProduct(arg,arg);

lib/qcd/Grid_qcd.h

@@ -7,6 +7,7 @@ namespace QCD {
     static const int Nc=3;
     static const int Ns=4;
     static const int Nd=4;
+    static const int Nhs=2; // half spinor
 
     static const int CbRed  =0;
     static const int CbBlack=1;
@@ -38,6 +39,9 @@ namespace QCD {
     template<typename vtype> using iColourVector     = iScalar<iScalar<iVector<vtype, Nc> > >;
     template<typename vtype> using iSpinColourVector = iScalar<iVector<iVector<vtype, Nc>, Ns> >;
 
+    template<typename vtype> using iHalfSpinVector       = iScalar<iVector<iScalar<vtype>, Nhs> >;
+    template<typename vtype> using iHalfSpinColourVector = iScalar<iVector<iVector<vtype, Nc>, Nhs> >;
+
     typedef iSpinMatrix<Complex >           SpinMatrix;
     typedef iColourMatrix<Complex >         ColourMatrix;
     typedef iSpinColourMatrix<Complex >     SpinColourMatrix;
@@ -48,6 +52,8 @@ namespace QCD {
     typedef iSpinVector<Complex >           SpinVector;
     typedef iColourVector<Complex >         ColourVector;
     typedef iSpinColourVector<Complex >     SpinColourVector;
+    typedef iHalfSpinVector<Complex >       HalfSpinVector;
+    typedef iHalfSpinColourVector<Complex > HalfSpinColourVector;
 
     
     typedef iSpinMatrix<vComplex >          vSpinMatrix;
@@ -58,8 +64,10 @@ namespace QCD {
     typedef iSpinVector<vComplex >           vSpinVector;
     typedef iColourVector<vComplex >         vColourVector;
     typedef iSpinColourVector<vComplex >     vSpinColourVector;
+    typedef iHalfSpinVector<vComplex >       vHalfSpinVector;
+    typedef iHalfSpinColourVector<vComplex > vHalfSpinColourVector;
     
-    typedef iSinglet<Complex >          TComplex;    // This is painful. Tensor singlet complex type.
+    typedef iSinglet<Complex >         TComplex;    // FIXME This is painful. Tensor singlet complex type.
     typedef iSinglet<vComplex >        vTComplex;   // what if we don't know the tensor structure
     typedef iSinglet<Real >            TReal;       // Shouldn't need these; can I make it work without?
     typedef iSinglet<vReal >           vTReal;      
@@ -74,17 +82,25 @@ namespace QCD {
     typedef Lattice<vSpinMatrix>       LatticeSpinMatrix;
     typedef Lattice<vSpinColourMatrix> LatticeSpinColourMatrix;
 
-    typedef Lattice<vSpinColourVector> LatticeSpinColourVector;
     typedef Lattice<vSpinVector>           LatticeSpinVector;
     typedef Lattice<vColourVector>         LatticeColourVector;
+    typedef Lattice<vSpinColourVector>     LatticeSpinColourVector;
+    typedef Lattice<vHalfSpinVector>       LatticeHalfSpinVector;
+    typedef Lattice<vHalfSpinColourVector> LatticeHalfSpinColourVector;
 
     ///////////////////////////////////////////
     // Physical names for things
     ///////////////////////////////////////////
+    typedef Lattice<vHalfSpinColourVector> LatticeHalfFermion;
     typedef Lattice<vSpinColourVector>     LatticeFermion;
+
     typedef Lattice<vSpinColourMatrix>     LatticePropagator;
     typedef Lattice<vLorentzColourMatrix>  LatticeGaugeField;
 
+    // Uhgg... typing this hurt  ;)
+    // (my keyboard got burning hot when I typed this, must be the anti-Fermion)
+    typedef Lattice<vColourVector>         LatticeStaggeredFermion;    
+    typedef Lattice<vColourMatrix>         LatticeStaggeredPropagator; 
 
     //////////////////////////////////////////////////////////////////////////////
     // Peek and Poke named after physics attributes
@@ -145,5 +161,7 @@ namespace QCD {
 } // Grid
 
 #include <qcd/Grid_qcd_dirac.h>
+#include <qcd/Grid_qcd_2spinor.h>
+//#include <qcd/Grid_qcd_pauli.h>
 
 #endif

tests/Grid_gamma.cc

@@ -29,42 +29,16 @@ int main (int argc, char ** argv)
   SpinMatrix rr=zero;
   SpinMatrix result;
 
-  SpinVector lv=zero;
-  SpinVector rv=zero;
+  SpinVector lv; random(sRNG,lv);
+  SpinVector rv; random(sRNG,rv);
 
   for(int a=0;a<Ns;a++){
     ident()(a,a) = 1.0;
   }
 
-  Gamma::GammaMatrix g [] = {
-    Gamma::Identity,
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT,
-    Gamma::Gamma5,
-    Gamma::MinusIdentity,
-    Gamma::MinusGammaX,
-    Gamma::MinusGammaY,
-    Gamma::MinusGammaZ,
-    Gamma::MinusGammaT,
-    Gamma::MinusGamma5
-  };
-  const char *list[] = { 
-    "Identity ",
-    "GammaX   ",
-    "GammaY   ",
-    "GammaZ   ",
-    "GammaT   ",
-    "Gamma5   ",
-    "-Identity",
-    "-GammaX  ",
-    "-GammaY  ",
-    "-GammaZ  ",
-    "-GammaT  ",
-    "-Gamma5  ",
-    "         "
-  };
+  const Gamma::GammaMatrix *g = Gamma::GammaMatrices;
+  const char **list          = Gamma::GammaMatrixNames;
+
   result =ll*Gamma(g[0])*rr;
   result =ll*Gamma(g[0]);
   rv = Gamma(g[0])*lv;
@@ -124,5 +98,81 @@ int main (int argc, char ** argv)
     std::cout << list[mu]<<" " << mag<<std::endl;
   }
 
+  // Testing spins and reconstructs
+  SpinVector     recon; random(sRNG,rv);
+  SpinVector      full;
+  HalfSpinVector hsp,hsm;
+
+  // Xp
+  double mag;
+  spProjXp(hsm,rv);
+  spReconXp(recon,hsm);
+  full = rv + Gamma(Gamma::GammaX) *rv;
+  mag = TensorRemove(norm2l(full-recon));
+  std::cout << "Xp "<< mag<<std::endl;
+
+  // Xm
+  spProjXm(hsm,rv);
+  spReconXm(recon,hsm);
+  full = rv - Gamma(Gamma::GammaX) *rv;
+  mag = TensorRemove(norm2l(full-recon));
+  std::cout << "Xm "<< mag<<std::endl;
+
+  // Yp
+  spProjYp(hsm,rv);
+  spReconYp(recon,hsm);
+  full = rv + Gamma(Gamma::GammaY) *rv;
+  mag = TensorRemove(norm2l(full-recon));
+  std::cout << "Yp "<< mag<<std::endl;
+
+  // Ym
+  spProjYm(hsm,rv);
+  spReconYm(recon,hsm);
+  full = rv - Gamma(Gamma::GammaY) *rv;
+  mag = TensorRemove(norm2l(full-recon));
+  std::cout << "Ym "<< mag<<std::endl;
+
+  // Zp
+  spProjZp(hsm,rv);
+  spReconZp(recon,hsm);
+  full = rv + Gamma(Gamma::GammaZ) *rv;
+  mag = TensorRemove(norm2l(full-recon));
+  std::cout << "Zp "<< mag<<std::endl;
+
+  // Zm
+  spProjZm(hsm,rv);
+  spReconZm(recon,hsm);
+  full = rv - Gamma(Gamma::GammaZ) *rv;
+  mag = TensorRemove(norm2l(full-recon));
+  std::cout << "Zm "<< mag<<std::endl;
+
+  // Tp
+  spProjTp(hsm,rv);
+  spReconTp(recon,hsm);
+  full = rv + Gamma(Gamma::GammaT) *rv;
+  mag = TensorRemove(norm2l(full-recon));
+  std::cout << "Tp "<< mag<<std::endl;
+
+  // Tm
+  spProjTm(hsm,rv);
+  spReconTm(recon,hsm);
+  full = rv - Gamma(Gamma::GammaT) *rv;
+  mag = TensorRemove(norm2l(full-recon));
+  std::cout << "Tm "<< mag<<std::endl;
+
+  // 5p
+  spProj5p(hsm,rv);
+  spRecon5p(recon,hsm);
+  full = rv + Gamma(Gamma::Gamma5) *rv;
+  mag = TensorRemove(norm2l(full-recon));
+  std::cout << "5p "<< mag<<std::endl;
+
+  // 5m
+  spProj5m(hsm,rv);
+  spRecon5m(recon,hsm);
+  full = rv - Gamma(Gamma::Gamma5) *rv;
+  mag = TensorRemove(norm2l(full-recon));
+  std::cout << "5m "<< mag<<std::endl;
+
   Grid_finalize();
 }