First cut at higher precision reduction

lib/lattice/Lattice_reduction.h

@@ -47,7 +47,7 @@ namespace Grid {
 inline ComplexD innerProduct(const Lattice<vobj> &left,const Lattice<vobj> &right) 
 {
   typedef typename vobj::scalar_type scalar_type;
-      typedef typename vobj::vector_type vector_type;
+  typedef typename vobj::vector_typeD vector_type;
   scalar_type  nrm;
   
   GridBase *grid = left._grid;
@@ -61,9 +61,9 @@ namespace Grid {
     int nwork, mywork, myoff;
     GridThread::GetWork(left._grid->oSites(),thr,mywork,myoff);
     
-	decltype(innerProduct(left._odata[0],right._odata[0])) vnrm=zero; // private to thread; sub summation
+    decltype(innerProductD(left._odata[0],right._odata[0])) vnrm=zero; // private to thread; sub summation
     for(int ss=myoff;ss<mywork+myoff; ss++){
-	  vnrm = vnrm + innerProduct(left._odata[ss],right._odata[ss]);
+      vnrm = vnrm + innerProductD(left._odata[ss],right._odata[ss]);
     }
     sumarray[thr]=TensorRemove(vnrm) ;
   }
@@ -103,8 +103,8 @@ namespace Grid {
 }
 
 template<class vobj>
-    inline typename vobj::scalar_object sum(const Lattice<vobj> &arg){
+inline typename vobj::scalar_object sum(const Lattice<vobj> &arg)
-
+{
   GridBase *grid=arg._grid;
   int Nsimd = grid->Nsimd();
   
@@ -141,8 +141,6 @@ namespace Grid {
   return ssum;
 }
 
-
-
 template<class vobj> inline void sliceSum(const Lattice<vobj> &Data,std::vector<typename vobj::scalar_object> &result,int orthogdim)
 {
   typedef typename vobj::scalar_object sobj;
@@ -214,7 +212,6 @@ template<class vobj> inline void sliceSum(const Lattice<vobj> &Data,std::vector<
 
     result[t]=gsum;
   }
-
 }
 
 

lib/simd/Grid_vector_types.h

@@ -694,7 +694,6 @@ inline Grid_simd<S, V> innerProduct(const Grid_simd<S, V> &l,
                                     const Grid_simd<S, V> &r) {
   return conjugate(l) * r;
 }
-
 template <class S, class V>
 inline Grid_simd<S, V> outerProduct(const Grid_simd<S, V> &l,
                                     const Grid_simd<S, V> &r) {

lib/tensors/Tensor_class.h

@@ -56,6 +56,7 @@ class iScalar {
   typedef vtype element;
   typedef typename GridTypeMapper<vtype>::scalar_type scalar_type;
   typedef typename GridTypeMapper<vtype>::vector_type vector_type;
+  typedef typename GridTypeMapper<vtype>::vector_typeD vector_typeD;
   typedef typename GridTypeMapper<vtype>::tensor_reduced tensor_reduced_v;
   typedef iScalar<tensor_reduced_v> tensor_reduced;
   typedef typename GridTypeMapper<vtype>::scalar_object recurse_scalar_object;
@@ -193,6 +194,7 @@ class iVector {
   typedef vtype element;
   typedef typename GridTypeMapper<vtype>::scalar_type scalar_type;
   typedef typename GridTypeMapper<vtype>::vector_type vector_type;
+  typedef typename GridTypeMapper<vtype>::vector_typeD vector_typeD;
   typedef typename GridTypeMapper<vtype>::tensor_reduced tensor_reduced_v;
   typedef typename GridTypeMapper<vtype>::scalar_object recurse_scalar_object;
   typedef iScalar<tensor_reduced_v> tensor_reduced;
@@ -305,6 +307,7 @@ class iMatrix {
   typedef vtype element;
   typedef typename GridTypeMapper<vtype>::scalar_type scalar_type;
   typedef typename GridTypeMapper<vtype>::vector_type vector_type;
+  typedef typename GridTypeMapper<vtype>::vector_typeD vector_typeD;
   typedef typename GridTypeMapper<vtype>::tensor_reduced tensor_reduced_v;
   typedef typename GridTypeMapper<vtype>::scalar_object recurse_scalar_object;
 

lib/tensors/Tensor_inner.h

@@ -35,13 +35,71 @@ namespace Grid {
   // innerProduct Matrix x Matrix -> Scalar
   ///////////////////////////////////////////////////////////////////////////////////////
   template<class sobj> inline RealD norm2(const sobj &arg){
-      typedef typename sobj::scalar_type scalar;
+    auto nrm = innerProductD(arg,arg);
-      decltype(innerProduct(arg,arg)) nrm;
-      nrm = innerProduct(arg,arg);
     RealD ret = real(nrm);
     return ret;
   }
+  //////////////////////////////////////
+  // If single promote to double and sum 2x
+  //////////////////////////////////////
 
+inline ComplexD innerProductD(const ComplexF &l,const ComplexF &r){  return innerProduct(l,r); }
+inline ComplexD innerProductD(const ComplexD &l,const ComplexD &r){  return innerProduct(l,r); }
+inline RealD    innerProductD(const RealD    &l,const RealD    &r){  return innerProduct(l,r); }
+inline RealD    innerProductD(const RealF    &l,const RealF    &r){  return innerProduct(l,r); }
+
+inline vComplexD innerProductD(const vComplexD &l,const vComplexD &r){  return innerProduct(l,r); }
+inline vRealD    innerProductD(const vRealD    &l,const vRealD    &r){  return innerProduct(l,r); }
+inline vComplexD innerProductD(const vComplexF &l,const vComplexF &r){  
+  vComplexD la,lb;
+  vComplexD ra,rb;
+  Optimization::PrecisionChange::StoD(l.v,la.v,lb.v);
+  Optimization::PrecisionChange::StoD(r.v,ra.v,rb.v);
+  return innerProduct(la,ra) + innerProduct(lb,rb); 
+}
+inline vRealD innerProductD(const vRealF &l,const vRealF &r){  
+  vRealD la,lb;
+  vRealD ra,rb;
+  Optimization::PrecisionChange::StoD(l.v,la.v,lb.v);
+  Optimization::PrecisionChange::StoD(r.v,ra.v,rb.v);
+  return innerProduct(la,ra) + innerProduct(lb,rb); 
+}
+
+  template<class l,class r,int N> inline
+  auto innerProductD (const iVector<l,N>& lhs,const iVector<r,N>& rhs) -> iScalar<decltype(innerProductD(lhs._internal[0],rhs._internal[0]))>
+  {
+    typedef decltype(innerProductD(lhs._internal[0],rhs._internal[0])) ret_t;
+    iScalar<ret_t> ret;
+    ret=zero;
+    for(int c1=0;c1<N;c1++){
+      ret._internal += innerProductD(lhs._internal[c1],rhs._internal[c1]);
+    }
+    return ret;
+  }
+  template<class l,class r,int N> inline
+  auto innerProductD (const iMatrix<l,N>& lhs,const iMatrix<r,N>& rhs) -> iScalar<decltype(innerProductD(lhs._internal[0][0],rhs._internal[0][0]))>
+  {
+    typedef decltype(innerProductD(lhs._internal[0][0],rhs._internal[0][0])) ret_t;
+    iScalar<ret_t> ret;
+    iScalar<ret_t> tmp;
+    ret=zero;
+    for(int c1=0;c1<N;c1++){
+    for(int c2=0;c2<N;c2++){
+      ret._internal+=innerProductD(lhs._internal[c1][c2],rhs._internal[c1][c2]);
+    }}
+    return ret;
+  }
+  template<class l,class r> inline
+  auto innerProductD (const iScalar<l>& lhs,const iScalar<r>& rhs) -> iScalar<decltype(innerProductD(lhs._internal,rhs._internal))>
+  {
+    typedef decltype(innerProductD(lhs._internal,rhs._internal)) ret_t;
+    iScalar<ret_t> ret;
+    ret._internal = innerProductD(lhs._internal,rhs._internal);
+    return ret;
+  }
+  //////////////////////
+  // Keep same precison
+  //////////////////////
   template<class l,class r,int N> inline
   auto innerProduct (const iVector<l,N>& lhs,const iVector<r,N>& rhs) -> iScalar<decltype(innerProduct(lhs._internal[0],rhs._internal[0]))>
   {

lib/tensors/Tensor_traits.h

@@ -53,6 +53,7 @@ namespace Grid {
   public:
     typedef typename T::scalar_type scalar_type;
     typedef typename T::vector_type vector_type;
+    typedef typename T::vector_typeD vector_typeD;
     typedef typename T::tensor_reduced tensor_reduced;
     typedef typename T::scalar_object scalar_object;
     typedef typename T::Complexified Complexified;
@@ -67,6 +68,7 @@ namespace Grid {
   public:
     typedef RealF scalar_type;
     typedef RealF vector_type;
+    typedef RealD vector_typeD;
     typedef RealF tensor_reduced ;
     typedef RealF scalar_object;
     typedef ComplexF Complexified;
@@ -77,6 +79,7 @@ namespace Grid {
   public:
     typedef RealD scalar_type;
     typedef RealD vector_type;
+    typedef RealD vector_typeD;
     typedef RealD tensor_reduced;
     typedef RealD scalar_object;
     typedef ComplexD Complexified;
@@ -87,6 +90,7 @@ namespace Grid {
   public:
     typedef ComplexF scalar_type;
     typedef ComplexF vector_type;
+    typedef ComplexD vector_typeD;
     typedef ComplexF tensor_reduced;
     typedef ComplexF scalar_object;
     typedef ComplexF Complexified;
@@ -97,6 +101,7 @@ namespace Grid {
   public:
     typedef ComplexD scalar_type;
     typedef ComplexD vector_type;
+    typedef ComplexD vector_typeD;
     typedef ComplexD tensor_reduced;
     typedef ComplexD scalar_object;
     typedef ComplexD Complexified;
@@ -118,6 +123,7 @@ namespace Grid {
   public:
     typedef RealF  scalar_type;
     typedef vRealF vector_type;
+    typedef vRealD vector_typeD;
     typedef vRealF tensor_reduced;
     typedef RealF  scalar_object;
     typedef vComplexF Complexified;
@@ -128,6 +134,7 @@ namespace Grid {
   public:
     typedef RealD  scalar_type;
     typedef vRealD vector_type;
+    typedef vRealD vector_typeD;
     typedef vRealD tensor_reduced;
     typedef RealD  scalar_object;
     typedef vComplexD Complexified;
@@ -138,6 +145,7 @@ namespace Grid {
   public:
     typedef ComplexF  scalar_type;
     typedef vComplexF vector_type;
+    typedef vComplexD vector_typeD;
     typedef vComplexF tensor_reduced;
     typedef ComplexF  scalar_object;
     typedef vComplexF Complexified;
@@ -148,6 +156,7 @@ namespace Grid {
   public:
     typedef ComplexD  scalar_type;
     typedef vComplexD vector_type;
+    typedef vComplexD vector_typeD;
     typedef vComplexD tensor_reduced;
     typedef ComplexD  scalar_object;
     typedef vComplexD Complexified;
@@ -158,6 +167,7 @@ namespace Grid {
   public:
     typedef  Integer scalar_type;
     typedef vInteger vector_type;
+    typedef vInteger vector_typeD;
     typedef vInteger tensor_reduced;
     typedef  Integer scalar_object;
     typedef void Complexified;