mirror of
				https://github.com/paboyle/Grid.git
				synced 2025-11-04 05:54:32 +00:00 
			
		
		
		
	NAMESPACE and formatting
This commit is contained in:
		@@ -1,4 +1,4 @@
 | 
			
		||||
    /*************************************************************************************
 | 
			
		||||
/*************************************************************************************
 | 
			
		||||
 | 
			
		||||
    Grid physics library, www.github.com/paboyle/Grid 
 | 
			
		||||
 | 
			
		||||
@@ -25,413 +25,413 @@ Author: Antonin Portelli <antonin.portelli@me.com>
 | 
			
		||||
    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 */
 | 
			
		||||
*************************************************************************************/
 | 
			
		||||
/*  END LEGAL */
 | 
			
		||||
 | 
			
		||||
#include "Grid_generic_types.h"
 | 
			
		||||
 | 
			
		||||
namespace Grid {
 | 
			
		||||
namespace Optimization {
 | 
			
		||||
NAMESPACE_BEGIN(Grid);
 | 
			
		||||
NAMESPACE_BEGIN(Optimization);
 | 
			
		||||
  
 | 
			
		||||
  struct Vsplat{
 | 
			
		||||
    // Complex
 | 
			
		||||
    template <typename T>
 | 
			
		||||
    inline vec<T> operator()(T a, T b){
 | 
			
		||||
      vec<T> out;
 | 
			
		||||
struct Vsplat{
 | 
			
		||||
  // Complex
 | 
			
		||||
  template <typename T>
 | 
			
		||||
  inline vec<T> operator()(T a, T b){
 | 
			
		||||
    vec<T> out;
 | 
			
		||||
      
 | 
			
		||||
      VECTOR_FOR(i, W<T>::r, 2)
 | 
			
		||||
    VECTOR_FOR(i, W<T>::r, 2)
 | 
			
		||||
      {
 | 
			
		||||
        out.v[i]   = a;
 | 
			
		||||
        out.v[i+1] = b;
 | 
			
		||||
      }
 | 
			
		||||
 | 
			
		||||
      return out;
 | 
			
		||||
    }
 | 
			
		||||
    return out;
 | 
			
		||||
  }
 | 
			
		||||
    
 | 
			
		||||
    // Real
 | 
			
		||||
    template <typename T>
 | 
			
		||||
    inline vec<T> operator()(T a){
 | 
			
		||||
      vec<T> out;
 | 
			
		||||
  // Real
 | 
			
		||||
  template <typename T>
 | 
			
		||||
  inline vec<T> operator()(T a){
 | 
			
		||||
    vec<T> out;
 | 
			
		||||
      
 | 
			
		||||
      VECTOR_FOR(i, W<T>::r, 1)
 | 
			
		||||
    VECTOR_FOR(i, W<T>::r, 1)
 | 
			
		||||
      {
 | 
			
		||||
        out.v[i] = a;
 | 
			
		||||
      }
 | 
			
		||||
      
 | 
			
		||||
      return out;
 | 
			
		||||
    }
 | 
			
		||||
  };
 | 
			
		||||
    return out;
 | 
			
		||||
  }
 | 
			
		||||
};
 | 
			
		||||
 | 
			
		||||
  struct Vstore{
 | 
			
		||||
    // Real
 | 
			
		||||
    template <typename T>
 | 
			
		||||
    inline void operator()(vec<T> a, T *D){
 | 
			
		||||
      *((vec<T> *)D) = a;
 | 
			
		||||
    }
 | 
			
		||||
  };
 | 
			
		||||
struct Vstore{
 | 
			
		||||
  // Real
 | 
			
		||||
  template <typename T>
 | 
			
		||||
  inline void operator()(vec<T> a, T *D){
 | 
			
		||||
    *((vec<T> *)D) = a;
 | 
			
		||||
  }
 | 
			
		||||
};
 | 
			
		||||
 | 
			
		||||
  struct Vstream{
 | 
			
		||||
    // Real
 | 
			
		||||
    template <typename T>
 | 
			
		||||
    inline void operator()(T * a, vec<T> b){
 | 
			
		||||
      *((vec<T> *)a) = b;
 | 
			
		||||
    }
 | 
			
		||||
  };
 | 
			
		||||
struct Vstream{
 | 
			
		||||
  // Real
 | 
			
		||||
  template <typename T>
 | 
			
		||||
  inline void operator()(T * a, vec<T> b){
 | 
			
		||||
    *((vec<T> *)a) = b;
 | 
			
		||||
  }
 | 
			
		||||
};
 | 
			
		||||
 | 
			
		||||
  struct Vset{
 | 
			
		||||
    // Complex
 | 
			
		||||
    template <typename T>
 | 
			
		||||
    inline vec<T> operator()(std::complex<T> *a){
 | 
			
		||||
      vec<T> out;
 | 
			
		||||
struct Vset{
 | 
			
		||||
  // Complex
 | 
			
		||||
  template <typename T>
 | 
			
		||||
  inline vec<T> operator()(std::complex<T> *a){
 | 
			
		||||
    vec<T> out;
 | 
			
		||||
      
 | 
			
		||||
      VECTOR_FOR(i, W<T>::c, 1)
 | 
			
		||||
    VECTOR_FOR(i, W<T>::c, 1)
 | 
			
		||||
      {
 | 
			
		||||
        out.v[2*i]   = a[i].real();
 | 
			
		||||
        out.v[2*i+1] = a[i].imag();
 | 
			
		||||
      }
 | 
			
		||||
      
 | 
			
		||||
      return out;
 | 
			
		||||
    }
 | 
			
		||||
    return out;
 | 
			
		||||
  }
 | 
			
		||||
    
 | 
			
		||||
    // Real
 | 
			
		||||
    template <typename T>
 | 
			
		||||
    inline vec<T> operator()(T *a){
 | 
			
		||||
      vec<T> out;
 | 
			
		||||
  // Real
 | 
			
		||||
  template <typename T>
 | 
			
		||||
  inline vec<T> operator()(T *a){
 | 
			
		||||
    vec<T> out;
 | 
			
		||||
      
 | 
			
		||||
      out = *((vec<T> *)a);
 | 
			
		||||
    out = *((vec<T> *)a);
 | 
			
		||||
      
 | 
			
		||||
      return out;
 | 
			
		||||
    }
 | 
			
		||||
  };
 | 
			
		||||
    return out;
 | 
			
		||||
  }
 | 
			
		||||
};
 | 
			
		||||
 | 
			
		||||
  /////////////////////////////////////////////////////
 | 
			
		||||
  // Arithmetic operations
 | 
			
		||||
  /////////////////////////////////////////////////////
 | 
			
		||||
  struct Sum{
 | 
			
		||||
    // Complex/Real
 | 
			
		||||
    template <typename T>
 | 
			
		||||
    inline vec<T> operator()(vec<T> a, vec<T> b){
 | 
			
		||||
      vec<T> out;
 | 
			
		||||
/////////////////////////////////////////////////////
 | 
			
		||||
// Arithmetic operations
 | 
			
		||||
/////////////////////////////////////////////////////
 | 
			
		||||
struct Sum{
 | 
			
		||||
  // Complex/Real
 | 
			
		||||
  template <typename T>
 | 
			
		||||
  inline vec<T> operator()(vec<T> a, vec<T> b){
 | 
			
		||||
    vec<T> out;
 | 
			
		||||
      
 | 
			
		||||
      VECTOR_FOR(i, W<T>::r, 1)
 | 
			
		||||
    VECTOR_FOR(i, W<T>::r, 1)
 | 
			
		||||
      {
 | 
			
		||||
        out.v[i] = a.v[i] + b.v[i];
 | 
			
		||||
      }
 | 
			
		||||
      
 | 
			
		||||
      return out;
 | 
			
		||||
    }
 | 
			
		||||
  };
 | 
			
		||||
    return out;
 | 
			
		||||
  }
 | 
			
		||||
};
 | 
			
		||||
 | 
			
		||||
  struct Sub{
 | 
			
		||||
    // Complex/Real
 | 
			
		||||
    template <typename T>
 | 
			
		||||
    inline vec<T> operator()(vec<T> a, vec<T> b){
 | 
			
		||||
      vec<T> out;
 | 
			
		||||
struct Sub{
 | 
			
		||||
  // Complex/Real
 | 
			
		||||
  template <typename T>
 | 
			
		||||
  inline vec<T> operator()(vec<T> a, vec<T> b){
 | 
			
		||||
    vec<T> out;
 | 
			
		||||
      
 | 
			
		||||
      VECTOR_FOR(i, W<T>::r, 1)
 | 
			
		||||
    VECTOR_FOR(i, W<T>::r, 1)
 | 
			
		||||
      {
 | 
			
		||||
        out.v[i] = a.v[i] - b.v[i];
 | 
			
		||||
      }
 | 
			
		||||
      
 | 
			
		||||
      return out;
 | 
			
		||||
    }
 | 
			
		||||
  };
 | 
			
		||||
    return out;
 | 
			
		||||
  }
 | 
			
		||||
};
 | 
			
		||||
 | 
			
		||||
  struct Mult{
 | 
			
		||||
    // Real
 | 
			
		||||
    template <typename T>
 | 
			
		||||
    inline vec<T> operator()(vec<T> a, vec<T> b){
 | 
			
		||||
      vec<T> out;
 | 
			
		||||
struct Mult{
 | 
			
		||||
  // Real
 | 
			
		||||
  template <typename T>
 | 
			
		||||
  inline vec<T> operator()(vec<T> a, vec<T> b){
 | 
			
		||||
    vec<T> out;
 | 
			
		||||
      
 | 
			
		||||
      VECTOR_FOR(i, W<T>::r, 1)
 | 
			
		||||
    VECTOR_FOR(i, W<T>::r, 1)
 | 
			
		||||
      {
 | 
			
		||||
        out.v[i] = a.v[i]*b.v[i];
 | 
			
		||||
      }
 | 
			
		||||
      
 | 
			
		||||
      return out;
 | 
			
		||||
    }
 | 
			
		||||
  };
 | 
			
		||||
    return out;
 | 
			
		||||
  }
 | 
			
		||||
};
 | 
			
		||||
  
 | 
			
		||||
  #define cmul(a, b, c, i)\
 | 
			
		||||
  c[i]   = a[i]*b[i]   - a[i+1]*b[i+1];\
 | 
			
		||||
#define cmul(a, b, c, i)			\
 | 
			
		||||
  c[i]   = a[i]*b[i]   - a[i+1]*b[i+1];		\
 | 
			
		||||
  c[i+1] = a[i]*b[i+1] + a[i+1]*b[i];
 | 
			
		||||
 | 
			
		||||
  struct MultRealPart{
 | 
			
		||||
    template <typename T>
 | 
			
		||||
    inline vec<T> operator()(vec<T> a, vec<T> b){
 | 
			
		||||
      vec<T> out;
 | 
			
		||||
struct MultRealPart{
 | 
			
		||||
  template <typename T>
 | 
			
		||||
  inline vec<T> operator()(vec<T> a, vec<T> b){
 | 
			
		||||
    vec<T> out;
 | 
			
		||||
      
 | 
			
		||||
      VECTOR_FOR(i, W<T>::c, 1)
 | 
			
		||||
    VECTOR_FOR(i, W<T>::c, 1)
 | 
			
		||||
      {
 | 
			
		||||
         out.v[2*i]   = a.v[2*i]*b.v[2*i];
 | 
			
		||||
         out.v[2*i+1] = a.v[2*i]*b.v[2*i+1];
 | 
			
		||||
	out.v[2*i]   = a.v[2*i]*b.v[2*i];
 | 
			
		||||
	out.v[2*i+1] = a.v[2*i]*b.v[2*i+1];
 | 
			
		||||
      }      
 | 
			
		||||
      return out;
 | 
			
		||||
    }
 | 
			
		||||
  };
 | 
			
		||||
    return out;
 | 
			
		||||
  }
 | 
			
		||||
};
 | 
			
		||||
 | 
			
		||||
  struct MaddRealPart{
 | 
			
		||||
    template <typename T>
 | 
			
		||||
    inline vec<T> operator()(vec<T> a, vec<T> b, vec<T> c){
 | 
			
		||||
      vec<T> out;
 | 
			
		||||
struct MaddRealPart{
 | 
			
		||||
  template <typename T>
 | 
			
		||||
  inline vec<T> operator()(vec<T> a, vec<T> b, vec<T> c){
 | 
			
		||||
    vec<T> out;
 | 
			
		||||
      
 | 
			
		||||
      VECTOR_FOR(i, W<T>::c, 1)
 | 
			
		||||
    VECTOR_FOR(i, W<T>::c, 1)
 | 
			
		||||
      {
 | 
			
		||||
         out.v[2*i]   = a.v[2*i]*b.v[2*i] + c.v[2*i];
 | 
			
		||||
         out.v[2*i+1] = a.v[2*i]*b.v[2*i+1] + c.v[2*i+1];
 | 
			
		||||
	out.v[2*i]   = a.v[2*i]*b.v[2*i] + c.v[2*i];
 | 
			
		||||
	out.v[2*i+1] = a.v[2*i]*b.v[2*i+1] + c.v[2*i+1];
 | 
			
		||||
      }      
 | 
			
		||||
      return out;
 | 
			
		||||
    }
 | 
			
		||||
  };
 | 
			
		||||
    return out;
 | 
			
		||||
  }
 | 
			
		||||
};
 | 
			
		||||
  
 | 
			
		||||
  struct MultComplex{
 | 
			
		||||
    // Complex
 | 
			
		||||
    template <typename T>
 | 
			
		||||
    inline vec<T> operator()(vec<T> a, vec<T> b){
 | 
			
		||||
      vec<T> out;
 | 
			
		||||
struct MultComplex{
 | 
			
		||||
  // Complex
 | 
			
		||||
  template <typename T>
 | 
			
		||||
  inline vec<T> operator()(vec<T> a, vec<T> b){
 | 
			
		||||
    vec<T> out;
 | 
			
		||||
      
 | 
			
		||||
      VECTOR_FOR(i, W<T>::c, 1)
 | 
			
		||||
    VECTOR_FOR(i, W<T>::c, 1)
 | 
			
		||||
      {
 | 
			
		||||
        cmul(a.v, b.v, out.v, 2*i);
 | 
			
		||||
      }      
 | 
			
		||||
      
 | 
			
		||||
      return out;
 | 
			
		||||
    }
 | 
			
		||||
  };
 | 
			
		||||
    return out;
 | 
			
		||||
  }
 | 
			
		||||
};
 | 
			
		||||
  
 | 
			
		||||
  #undef cmul
 | 
			
		||||
#undef cmul
 | 
			
		||||
 | 
			
		||||
  struct Div{
 | 
			
		||||
    // Real
 | 
			
		||||
    template <typename T>
 | 
			
		||||
    inline vec<T> operator()(vec<T> a, vec<T> b){
 | 
			
		||||
      vec<T> out;
 | 
			
		||||
struct Div{
 | 
			
		||||
  // Real
 | 
			
		||||
  template <typename T>
 | 
			
		||||
  inline vec<T> operator()(vec<T> a, vec<T> b){
 | 
			
		||||
    vec<T> out;
 | 
			
		||||
      
 | 
			
		||||
      VECTOR_FOR(i, W<T>::r, 1)
 | 
			
		||||
    VECTOR_FOR(i, W<T>::r, 1)
 | 
			
		||||
      {
 | 
			
		||||
        out.v[i] = a.v[i]/b.v[i];
 | 
			
		||||
      }
 | 
			
		||||
      
 | 
			
		||||
      return out;
 | 
			
		||||
    }
 | 
			
		||||
  };
 | 
			
		||||
    return out;
 | 
			
		||||
  }
 | 
			
		||||
};
 | 
			
		||||
  
 | 
			
		||||
  #define conj(a, b, i)\
 | 
			
		||||
  b[i]   = a[i];\
 | 
			
		||||
#define conj(a, b, i)				\
 | 
			
		||||
  b[i]   = a[i];				\
 | 
			
		||||
  b[i+1] = -a[i+1];
 | 
			
		||||
  
 | 
			
		||||
  struct Conj{
 | 
			
		||||
    // Complex
 | 
			
		||||
    template <typename T>
 | 
			
		||||
    inline vec<T> operator()(vec<T> a){
 | 
			
		||||
      vec<T> out;
 | 
			
		||||
struct Conj{
 | 
			
		||||
  // Complex
 | 
			
		||||
  template <typename T>
 | 
			
		||||
  inline vec<T> operator()(vec<T> a){
 | 
			
		||||
    vec<T> out;
 | 
			
		||||
      
 | 
			
		||||
      VECTOR_FOR(i, W<T>::c, 1)
 | 
			
		||||
    VECTOR_FOR(i, W<T>::c, 1)
 | 
			
		||||
      {
 | 
			
		||||
        conj(a.v, out.v, 2*i);
 | 
			
		||||
      }
 | 
			
		||||
      
 | 
			
		||||
      return out;
 | 
			
		||||
    }
 | 
			
		||||
  };
 | 
			
		||||
    return out;
 | 
			
		||||
  }
 | 
			
		||||
};
 | 
			
		||||
  
 | 
			
		||||
  #undef conj
 | 
			
		||||
#undef conj
 | 
			
		||||
 | 
			
		||||
  #define timesmi(a, b, i)\
 | 
			
		||||
  b[i]   = a[i+1];\
 | 
			
		||||
#define timesmi(a, b, i)			\
 | 
			
		||||
  b[i]   = a[i+1];				\
 | 
			
		||||
  b[i+1] = -a[i];
 | 
			
		||||
  
 | 
			
		||||
  struct TimesMinusI{
 | 
			
		||||
    // Complex
 | 
			
		||||
    template <typename T>
 | 
			
		||||
    inline vec<T> operator()(vec<T> a, vec<T> b){
 | 
			
		||||
      vec<T> out;
 | 
			
		||||
struct TimesMinusI{
 | 
			
		||||
  // Complex
 | 
			
		||||
  template <typename T>
 | 
			
		||||
  inline vec<T> operator()(vec<T> a, vec<T> b){
 | 
			
		||||
    vec<T> out;
 | 
			
		||||
      
 | 
			
		||||
      VECTOR_FOR(i, W<T>::c, 1)
 | 
			
		||||
    VECTOR_FOR(i, W<T>::c, 1)
 | 
			
		||||
      {
 | 
			
		||||
        timesmi(a.v, out.v, 2*i);
 | 
			
		||||
      }
 | 
			
		||||
      
 | 
			
		||||
      return out;
 | 
			
		||||
    }
 | 
			
		||||
  };
 | 
			
		||||
    return out;
 | 
			
		||||
  }
 | 
			
		||||
};
 | 
			
		||||
 | 
			
		||||
  #undef timesmi
 | 
			
		||||
#undef timesmi
 | 
			
		||||
  
 | 
			
		||||
  #define timesi(a, b, i)\
 | 
			
		||||
  b[i]   = -a[i+1];\
 | 
			
		||||
#define timesi(a, b, i)				\
 | 
			
		||||
  b[i]   = -a[i+1];				\
 | 
			
		||||
  b[i+1] = a[i];
 | 
			
		||||
  
 | 
			
		||||
  struct TimesI{
 | 
			
		||||
    // Complex
 | 
			
		||||
    template <typename T>
 | 
			
		||||
    inline vec<T> operator()(vec<T> a, vec<T> b){
 | 
			
		||||
      vec<T> out;
 | 
			
		||||
struct TimesI{
 | 
			
		||||
  // Complex
 | 
			
		||||
  template <typename T>
 | 
			
		||||
  inline vec<T> operator()(vec<T> a, vec<T> b){
 | 
			
		||||
    vec<T> out;
 | 
			
		||||
      
 | 
			
		||||
      VECTOR_FOR(i, W<T>::c, 1)
 | 
			
		||||
    VECTOR_FOR(i, W<T>::c, 1)
 | 
			
		||||
      {
 | 
			
		||||
        timesi(a.v, out.v, 2*i);
 | 
			
		||||
      }
 | 
			
		||||
      
 | 
			
		||||
      return out;
 | 
			
		||||
    }
 | 
			
		||||
  };
 | 
			
		||||
    return out;
 | 
			
		||||
  }
 | 
			
		||||
};
 | 
			
		||||
  
 | 
			
		||||
  #undef timesi
 | 
			
		||||
#undef timesi
 | 
			
		||||
 | 
			
		||||
  struct PrecisionChange {
 | 
			
		||||
    static inline vech StoH (const vecf &a,const vecf &b) {
 | 
			
		||||
      vech ret;
 | 
			
		||||
struct PrecisionChange {
 | 
			
		||||
  static inline vech StoH (const vecf &a,const vecf &b) {
 | 
			
		||||
    vech ret;
 | 
			
		||||
#ifdef USE_FP16
 | 
			
		||||
      vech *ha = (vech *)&a;
 | 
			
		||||
      vech *hb = (vech *)&b;
 | 
			
		||||
      const int nf = W<float>::r;
 | 
			
		||||
      //      VECTOR_FOR(i, nf,1){ ret.v[i]    = ( (uint16_t *) &a.v[i])[1] ; }
 | 
			
		||||
      //      VECTOR_FOR(i, nf,1){ ret.v[i+nf] = ( (uint16_t *) &b.v[i])[1] ; }
 | 
			
		||||
      VECTOR_FOR(i, nf,1){ ret.v[i]    = ha->v[2*i+1]; }
 | 
			
		||||
      VECTOR_FOR(i, nf,1){ ret.v[i+nf] = hb->v[2*i+1]; }
 | 
			
		||||
    vech *ha = (vech *)&a;
 | 
			
		||||
    vech *hb = (vech *)&b;
 | 
			
		||||
    const int nf = W<float>::r;
 | 
			
		||||
    //      VECTOR_FOR(i, nf,1){ ret.v[i]    = ( (uint16_t *) &a.v[i])[1] ; }
 | 
			
		||||
    //      VECTOR_FOR(i, nf,1){ ret.v[i+nf] = ( (uint16_t *) &b.v[i])[1] ; }
 | 
			
		||||
    VECTOR_FOR(i, nf,1){ ret.v[i]    = ha->v[2*i+1]; }
 | 
			
		||||
    VECTOR_FOR(i, nf,1){ ret.v[i+nf] = hb->v[2*i+1]; }
 | 
			
		||||
#else
 | 
			
		||||
      assert(0);
 | 
			
		||||
    assert(0);
 | 
			
		||||
#endif
 | 
			
		||||
      return ret;
 | 
			
		||||
    }
 | 
			
		||||
    static inline void  HtoS (vech h,vecf &sa,vecf &sb) {
 | 
			
		||||
    return ret;
 | 
			
		||||
  }
 | 
			
		||||
  static inline void  HtoS (vech h,vecf &sa,vecf &sb) {
 | 
			
		||||
#ifdef USE_FP16
 | 
			
		||||
      const int nf = W<float>::r;
 | 
			
		||||
      const int nh = W<uint16_t>::r;
 | 
			
		||||
      vech *ha = (vech *)&sa;
 | 
			
		||||
      vech *hb = (vech *)&sb;
 | 
			
		||||
      VECTOR_FOR(i, nf, 1){ sb.v[i]= sa.v[i] = 0; }
 | 
			
		||||
      //      VECTOR_FOR(i, nf, 1){ ( (uint16_t *) (&sa.v[i]))[1] = h.v[i];}
 | 
			
		||||
      //      VECTOR_FOR(i, nf, 1){ ( (uint16_t *) (&sb.v[i]))[1] = h.v[i+nf];}
 | 
			
		||||
      VECTOR_FOR(i, nf, 1){ ha->v[2*i+1]=h.v[i]; }
 | 
			
		||||
      VECTOR_FOR(i, nf, 1){ hb->v[2*i+1]=h.v[i+nf]; }
 | 
			
		||||
    const int nf = W<float>::r;
 | 
			
		||||
    const int nh = W<uint16_t>::r;
 | 
			
		||||
    vech *ha = (vech *)&sa;
 | 
			
		||||
    vech *hb = (vech *)&sb;
 | 
			
		||||
    VECTOR_FOR(i, nf, 1){ sb.v[i]= sa.v[i] = 0; }
 | 
			
		||||
    //      VECTOR_FOR(i, nf, 1){ ( (uint16_t *) (&sa.v[i]))[1] = h.v[i];}
 | 
			
		||||
    //      VECTOR_FOR(i, nf, 1){ ( (uint16_t *) (&sb.v[i]))[1] = h.v[i+nf];}
 | 
			
		||||
    VECTOR_FOR(i, nf, 1){ ha->v[2*i+1]=h.v[i]; }
 | 
			
		||||
    VECTOR_FOR(i, nf, 1){ hb->v[2*i+1]=h.v[i+nf]; }
 | 
			
		||||
#else
 | 
			
		||||
      assert(0);
 | 
			
		||||
    assert(0);
 | 
			
		||||
#endif
 | 
			
		||||
    }
 | 
			
		||||
    static inline vecf DtoS (vecd a,vecd b) {
 | 
			
		||||
      const int nd = W<double>::r;
 | 
			
		||||
      const int nf = W<float>::r;
 | 
			
		||||
      vecf ret;
 | 
			
		||||
      VECTOR_FOR(i, nd,1){ ret.v[i]    = a.v[i] ; }
 | 
			
		||||
      VECTOR_FOR(i, nd,1){ ret.v[i+nd] = b.v[i] ; }
 | 
			
		||||
      return ret;
 | 
			
		||||
    }
 | 
			
		||||
    static inline void StoD (vecf s,vecd &a,vecd &b) {
 | 
			
		||||
      const int nd = W<double>::r;
 | 
			
		||||
      VECTOR_FOR(i, nd,1){ a.v[i] = s.v[i] ; }
 | 
			
		||||
      VECTOR_FOR(i, nd,1){ b.v[i] = s.v[i+nd] ; }
 | 
			
		||||
    }
 | 
			
		||||
    static inline vech DtoH (vecd a,vecd b,vecd c,vecd d) {
 | 
			
		||||
      vecf sa,sb;
 | 
			
		||||
      sa = DtoS(a,b);
 | 
			
		||||
      sb = DtoS(c,d);
 | 
			
		||||
      return StoH(sa,sb);
 | 
			
		||||
    }
 | 
			
		||||
    static inline void HtoD (vech h,vecd &a,vecd &b,vecd &c,vecd &d) {
 | 
			
		||||
      vecf sa,sb;
 | 
			
		||||
      HtoS(h,sa,sb);
 | 
			
		||||
      StoD(sa,a,b);
 | 
			
		||||
      StoD(sb,c,d);
 | 
			
		||||
    }
 | 
			
		||||
  };
 | 
			
		||||
  }
 | 
			
		||||
  static inline vecf DtoS (vecd a,vecd b) {
 | 
			
		||||
    const int nd = W<double>::r;
 | 
			
		||||
    const int nf = W<float>::r;
 | 
			
		||||
    vecf ret;
 | 
			
		||||
    VECTOR_FOR(i, nd,1){ ret.v[i]    = a.v[i] ; }
 | 
			
		||||
    VECTOR_FOR(i, nd,1){ ret.v[i+nd] = b.v[i] ; }
 | 
			
		||||
    return ret;
 | 
			
		||||
  }
 | 
			
		||||
  static inline void StoD (vecf s,vecd &a,vecd &b) {
 | 
			
		||||
    const int nd = W<double>::r;
 | 
			
		||||
    VECTOR_FOR(i, nd,1){ a.v[i] = s.v[i] ; }
 | 
			
		||||
    VECTOR_FOR(i, nd,1){ b.v[i] = s.v[i+nd] ; }
 | 
			
		||||
  }
 | 
			
		||||
  static inline vech DtoH (vecd a,vecd b,vecd c,vecd d) {
 | 
			
		||||
    vecf sa,sb;
 | 
			
		||||
    sa = DtoS(a,b);
 | 
			
		||||
    sb = DtoS(c,d);
 | 
			
		||||
    return StoH(sa,sb);
 | 
			
		||||
  }
 | 
			
		||||
  static inline void HtoD (vech h,vecd &a,vecd &b,vecd &c,vecd &d) {
 | 
			
		||||
    vecf sa,sb;
 | 
			
		||||
    HtoS(h,sa,sb);
 | 
			
		||||
    StoD(sa,a,b);
 | 
			
		||||
    StoD(sb,c,d);
 | 
			
		||||
  }
 | 
			
		||||
};
 | 
			
		||||
 | 
			
		||||
  //////////////////////////////////////////////
 | 
			
		||||
  // Exchange support
 | 
			
		||||
  struct Exchange{
 | 
			
		||||
//////////////////////////////////////////////
 | 
			
		||||
// Exchange support
 | 
			
		||||
struct Exchange{
 | 
			
		||||
 | 
			
		||||
    template <typename T,int n>
 | 
			
		||||
  template <typename T,int n>
 | 
			
		||||
    static inline void ExchangeN(vec<T> &out1,vec<T> &out2,vec<T> &in1,vec<T> &in2){
 | 
			
		||||
      const int w = W<T>::r;
 | 
			
		||||
      unsigned int mask = w >> (n + 1);
 | 
			
		||||
      //      std::cout << " Exchange "<<n<<" nsimd "<<w<<" mask 0x" <<std::hex<<mask<<std::dec<<std::endl;
 | 
			
		||||
      VECTOR_FOR(i, w, 1) {	
 | 
			
		||||
	int j1 = i&(~mask);
 | 
			
		||||
	if  ( (i&mask) == 0 ) { out1.v[i]=in1.v[j1];}
 | 
			
		||||
	else                  { out1.v[i]=in2.v[j1];}
 | 
			
		||||
	int j2 = i|mask;
 | 
			
		||||
	if  ( (i&mask) == 0 ) { out2.v[i]=in1.v[j2];}
 | 
			
		||||
	else                  { out2.v[i]=in2.v[j2];}
 | 
			
		||||
      }      
 | 
			
		||||
    }
 | 
			
		||||
    template <typename T>
 | 
			
		||||
    static inline void Exchange0(vec<T> &out1,vec<T> &out2,vec<T> &in1,vec<T> &in2){
 | 
			
		||||
      ExchangeN<T,0>(out1,out2,in1,in2);
 | 
			
		||||
    };
 | 
			
		||||
    template <typename T>
 | 
			
		||||
    static inline void Exchange1(vec<T> &out1,vec<T> &out2,vec<T> &in1,vec<T> &in2){
 | 
			
		||||
      ExchangeN<T,1>(out1,out2,in1,in2);
 | 
			
		||||
    };
 | 
			
		||||
    template <typename T>
 | 
			
		||||
    static inline void Exchange2(vec<T> &out1,vec<T> &out2,vec<T> &in1,vec<T> &in2){
 | 
			
		||||
      ExchangeN<T,2>(out1,out2,in1,in2);
 | 
			
		||||
    };
 | 
			
		||||
    template <typename T>
 | 
			
		||||
    static inline void Exchange3(vec<T> &out1,vec<T> &out2,vec<T> &in1,vec<T> &in2){
 | 
			
		||||
      ExchangeN<T,3>(out1,out2,in1,in2);
 | 
			
		||||
    };
 | 
			
		||||
    const int w = W<T>::r;
 | 
			
		||||
    unsigned int mask = w >> (n + 1);
 | 
			
		||||
    //      std::cout << " Exchange "<<n<<" nsimd "<<w<<" mask 0x" <<std::hex<<mask<<std::dec<<std::endl;
 | 
			
		||||
    VECTOR_FOR(i, w, 1) {	
 | 
			
		||||
      int j1 = i&(~mask);
 | 
			
		||||
      if  ( (i&mask) == 0 ) { out1.v[i]=in1.v[j1];}
 | 
			
		||||
      else                  { out1.v[i]=in2.v[j1];}
 | 
			
		||||
      int j2 = i|mask;
 | 
			
		||||
      if  ( (i&mask) == 0 ) { out2.v[i]=in1.v[j2];}
 | 
			
		||||
      else                  { out2.v[i]=in2.v[j2];}
 | 
			
		||||
    }      
 | 
			
		||||
  }
 | 
			
		||||
  template <typename T>
 | 
			
		||||
  static inline void Exchange0(vec<T> &out1,vec<T> &out2,vec<T> &in1,vec<T> &in2){
 | 
			
		||||
    ExchangeN<T,0>(out1,out2,in1,in2);
 | 
			
		||||
  };
 | 
			
		||||
  template <typename T>
 | 
			
		||||
  static inline void Exchange1(vec<T> &out1,vec<T> &out2,vec<T> &in1,vec<T> &in2){
 | 
			
		||||
    ExchangeN<T,1>(out1,out2,in1,in2);
 | 
			
		||||
  };
 | 
			
		||||
  template <typename T>
 | 
			
		||||
  static inline void Exchange2(vec<T> &out1,vec<T> &out2,vec<T> &in1,vec<T> &in2){
 | 
			
		||||
    ExchangeN<T,2>(out1,out2,in1,in2);
 | 
			
		||||
  };
 | 
			
		||||
  template <typename T>
 | 
			
		||||
  static inline void Exchange3(vec<T> &out1,vec<T> &out2,vec<T> &in1,vec<T> &in2){
 | 
			
		||||
    ExchangeN<T,3>(out1,out2,in1,in2);
 | 
			
		||||
  };
 | 
			
		||||
};
 | 
			
		||||
 | 
			
		||||
 | 
			
		||||
  //////////////////////////////////////////////
 | 
			
		||||
  // Some Template specialization
 | 
			
		||||
  #define perm(a, b, n, w)\
 | 
			
		||||
  unsigned int _mask = w >> (n + 1);\
 | 
			
		||||
  VECTOR_FOR(i, w, 1)\
 | 
			
		||||
  {\
 | 
			
		||||
    b[i] = a[i^_mask];\
 | 
			
		||||
//////////////////////////////////////////////
 | 
			
		||||
// Some Template specialization
 | 
			
		||||
#define perm(a, b, n, w)			\
 | 
			
		||||
  unsigned int _mask = w >> (n + 1);		\
 | 
			
		||||
  VECTOR_FOR(i, w, 1)				\
 | 
			
		||||
  {						\
 | 
			
		||||
    b[i] = a[i^_mask];				\
 | 
			
		||||
  }
 | 
			
		||||
  
 | 
			
		||||
  #define DECL_PERMUTE_N(n)\
 | 
			
		||||
  template <typename T>\
 | 
			
		||||
  static inline vec<T> Permute##n(vec<T> in) {\
 | 
			
		||||
    vec<T> out;\
 | 
			
		||||
    perm(in.v, out.v, n, W<T>::r);\
 | 
			
		||||
    return out;\
 | 
			
		||||
#define DECL_PERMUTE_N(n)			\
 | 
			
		||||
  template <typename T>				\
 | 
			
		||||
  static inline vec<T> Permute##n(vec<T> in) {	\
 | 
			
		||||
    vec<T> out;					\
 | 
			
		||||
    perm(in.v, out.v, n, W<T>::r);		\
 | 
			
		||||
    return out;					\
 | 
			
		||||
  }
 | 
			
		||||
  
 | 
			
		||||
  struct Permute{
 | 
			
		||||
    DECL_PERMUTE_N(0);
 | 
			
		||||
    DECL_PERMUTE_N(1);
 | 
			
		||||
    DECL_PERMUTE_N(2);
 | 
			
		||||
    DECL_PERMUTE_N(3);
 | 
			
		||||
  };
 | 
			
		||||
struct Permute{
 | 
			
		||||
  DECL_PERMUTE_N(0);
 | 
			
		||||
  DECL_PERMUTE_N(1);
 | 
			
		||||
  DECL_PERMUTE_N(2);
 | 
			
		||||
  DECL_PERMUTE_N(3);
 | 
			
		||||
};
 | 
			
		||||
  
 | 
			
		||||
  #undef perm
 | 
			
		||||
  #undef DECL_PERMUTE_N
 | 
			
		||||
#undef perm
 | 
			
		||||
#undef DECL_PERMUTE_N
 | 
			
		||||
  
 | 
			
		||||
  #define rot(a, b, n, w)\
 | 
			
		||||
  VECTOR_FOR(i, w, 1)\
 | 
			
		||||
  {\
 | 
			
		||||
    b[i] = a[(i + n)%w];\
 | 
			
		||||
#define rot(a, b, n, w)				\
 | 
			
		||||
  VECTOR_FOR(i, w, 1)				\
 | 
			
		||||
  {						\
 | 
			
		||||
    b[i] = a[(i + n)%w];			\
 | 
			
		||||
  }
 | 
			
		||||
  
 | 
			
		||||
  struct Rotate{
 | 
			
		||||
struct Rotate{
 | 
			
		||||
      
 | 
			
		||||
    template <int n, typename T> static inline vec<T> tRotate(vec<T> in){
 | 
			
		||||
      return rotate(in, n);
 | 
			
		||||
    }
 | 
			
		||||
  template <int n, typename T> static inline vec<T> tRotate(vec<T> in){
 | 
			
		||||
    return rotate(in, n);
 | 
			
		||||
  }
 | 
			
		||||
    
 | 
			
		||||
    template <typename T>
 | 
			
		||||
    static inline vec<T> rotate(vec<T> in, int n){
 | 
			
		||||
      vec<T> out;
 | 
			
		||||
  template <typename T>
 | 
			
		||||
  static inline vec<T> rotate(vec<T> in, int n){
 | 
			
		||||
    vec<T> out;
 | 
			
		||||
      
 | 
			
		||||
      rot(in.v, out.v, n, W<T>::r);
 | 
			
		||||
    rot(in.v, out.v, n, W<T>::r);
 | 
			
		||||
      
 | 
			
		||||
      return out;
 | 
			
		||||
    }
 | 
			
		||||
  };
 | 
			
		||||
 | 
			
		||||
  #undef rot
 | 
			
		||||
  
 | 
			
		||||
  #define acc(v, a, off, step, n)\
 | 
			
		||||
  for (unsigned int i = off; i < n; i += step)\
 | 
			
		||||
  {\
 | 
			
		||||
    a += v[i];\
 | 
			
		||||
    return out;
 | 
			
		||||
  }
 | 
			
		||||
};
 | 
			
		||||
 | 
			
		||||
#undef rot
 | 
			
		||||
  
 | 
			
		||||
  template <typename Out_type, typename In_type>
 | 
			
		||||
#define acc(v, a, off, step, n)			\
 | 
			
		||||
  for (unsigned int i = off; i < n; i += step)	\
 | 
			
		||||
    {						\
 | 
			
		||||
      a += v[i];				\
 | 
			
		||||
    }
 | 
			
		||||
  
 | 
			
		||||
template <typename Out_type, typename In_type>
 | 
			
		||||
  struct Reduce{
 | 
			
		||||
    //Need templated class to overload output type
 | 
			
		||||
    //General form must generate error if compiled
 | 
			
		||||
@@ -442,89 +442,91 @@ namespace Optimization {
 | 
			
		||||
    }
 | 
			
		||||
  };
 | 
			
		||||
  
 | 
			
		||||
  //Complex float Reduce
 | 
			
		||||
  template <>
 | 
			
		||||
  inline Grid::ComplexF Reduce<Grid::ComplexF, vecf>::operator()(vecf in){
 | 
			
		||||
    float a = 0.f, b = 0.f;
 | 
			
		||||
//Complex float Reduce
 | 
			
		||||
template <>
 | 
			
		||||
inline Grid::ComplexF Reduce<Grid::ComplexF, vecf>::operator()(vecf in){
 | 
			
		||||
  float a = 0.f, b = 0.f;
 | 
			
		||||
    
 | 
			
		||||
    acc(in.v, a, 0, 2, W<float>::r);
 | 
			
		||||
    acc(in.v, b, 1, 2, W<float>::r);
 | 
			
		||||
  acc(in.v, a, 0, 2, W<float>::r);
 | 
			
		||||
  acc(in.v, b, 1, 2, W<float>::r);
 | 
			
		||||
    
 | 
			
		||||
    return Grid::ComplexF(a, b);
 | 
			
		||||
  }
 | 
			
		||||
  
 | 
			
		||||
  //Real float Reduce
 | 
			
		||||
  template<>
 | 
			
		||||
  inline Grid::RealF Reduce<Grid::RealF, vecf>::operator()(vecf in){
 | 
			
		||||
    float a = 0.;
 | 
			
		||||
    
 | 
			
		||||
    acc(in.v, a, 0, 1, W<float>::r);
 | 
			
		||||
    
 | 
			
		||||
    return a;
 | 
			
		||||
  }
 | 
			
		||||
  
 | 
			
		||||
  //Complex double Reduce
 | 
			
		||||
  template<>
 | 
			
		||||
  inline Grid::ComplexD Reduce<Grid::ComplexD, vecd>::operator()(vecd in){
 | 
			
		||||
    double a = 0., b = 0.;
 | 
			
		||||
    
 | 
			
		||||
    acc(in.v, a, 0, 2, W<double>::r);
 | 
			
		||||
    acc(in.v, b, 1, 2, W<double>::r);
 | 
			
		||||
    
 | 
			
		||||
    return Grid::ComplexD(a, b);
 | 
			
		||||
  }
 | 
			
		||||
  
 | 
			
		||||
  //Real double Reduce
 | 
			
		||||
  template<>
 | 
			
		||||
  inline Grid::RealD Reduce<Grid::RealD, vecd>::operator()(vecd in){
 | 
			
		||||
    double a = 0.f;
 | 
			
		||||
    
 | 
			
		||||
    acc(in.v, a, 0, 1, W<double>::r);
 | 
			
		||||
    
 | 
			
		||||
    return a;
 | 
			
		||||
  }
 | 
			
		||||
 | 
			
		||||
  //Integer Reduce
 | 
			
		||||
  template<>
 | 
			
		||||
  inline Integer Reduce<Integer, veci>::operator()(veci in){
 | 
			
		||||
    Integer a = 0;
 | 
			
		||||
    
 | 
			
		||||
    acc(in.v, a, 0, 1, W<Integer>::r);
 | 
			
		||||
    
 | 
			
		||||
    return a;
 | 
			
		||||
  }
 | 
			
		||||
 | 
			
		||||
  #undef acc  // EIGEN compatibility
 | 
			
		||||
  return Grid::ComplexF(a, b);
 | 
			
		||||
}
 | 
			
		||||
  
 | 
			
		||||
//Real float Reduce
 | 
			
		||||
template<>
 | 
			
		||||
inline Grid::RealF Reduce<Grid::RealF, vecf>::operator()(vecf in){
 | 
			
		||||
  float a = 0.;
 | 
			
		||||
    
 | 
			
		||||
  acc(in.v, a, 0, 1, W<float>::r);
 | 
			
		||||
    
 | 
			
		||||
  return a;
 | 
			
		||||
}
 | 
			
		||||
  
 | 
			
		||||
//Complex double Reduce
 | 
			
		||||
template<>
 | 
			
		||||
inline Grid::ComplexD Reduce<Grid::ComplexD, vecd>::operator()(vecd in){
 | 
			
		||||
  double a = 0., b = 0.;
 | 
			
		||||
    
 | 
			
		||||
  acc(in.v, a, 0, 2, W<double>::r);
 | 
			
		||||
  acc(in.v, b, 1, 2, W<double>::r);
 | 
			
		||||
    
 | 
			
		||||
  return Grid::ComplexD(a, b);
 | 
			
		||||
}
 | 
			
		||||
  
 | 
			
		||||
//Real double Reduce
 | 
			
		||||
template<>
 | 
			
		||||
inline Grid::RealD Reduce<Grid::RealD, vecd>::operator()(vecd in){
 | 
			
		||||
  double a = 0.f;
 | 
			
		||||
    
 | 
			
		||||
  acc(in.v, a, 0, 1, W<double>::r);
 | 
			
		||||
    
 | 
			
		||||
  return a;
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
//Integer Reduce
 | 
			
		||||
template<>
 | 
			
		||||
inline Integer Reduce<Integer, veci>::operator()(veci in){
 | 
			
		||||
  Integer a = 0;
 | 
			
		||||
    
 | 
			
		||||
  acc(in.v, a, 0, 1, W<Integer>::r);
 | 
			
		||||
    
 | 
			
		||||
  return a;
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
#undef acc  // EIGEN compatibility
 | 
			
		||||
NAMESPACE_END(Optimization)
 | 
			
		||||
 | 
			
		||||
//////////////////////////////////////////////////////////////////////////////////////
 | 
			
		||||
// Here assign types 
 | 
			
		||||
 | 
			
		||||
  typedef Optimization::vech SIMD_Htype; // Reduced precision type
 | 
			
		||||
  typedef Optimization::vecf SIMD_Ftype; // Single precision type
 | 
			
		||||
  typedef Optimization::vecd SIMD_Dtype; // Double precision type
 | 
			
		||||
  typedef Optimization::veci SIMD_Itype; // Integer type
 | 
			
		||||
typedef Optimization::vech SIMD_Htype; // Reduced precision type
 | 
			
		||||
typedef Optimization::vecf SIMD_Ftype; // Single precision type
 | 
			
		||||
typedef Optimization::vecd SIMD_Dtype; // Double precision type
 | 
			
		||||
typedef Optimization::veci SIMD_Itype; // Integer type
 | 
			
		||||
 | 
			
		||||
  // prefetch utilities
 | 
			
		||||
  inline void v_prefetch0(int size, const char *ptr){};
 | 
			
		||||
  inline void prefetch_HINT_T0(const char *ptr){};
 | 
			
		||||
// prefetch utilities
 | 
			
		||||
inline void v_prefetch0(int size, const char *ptr){};
 | 
			
		||||
inline void prefetch_HINT_T0(const char *ptr){};
 | 
			
		||||
 | 
			
		||||
  // Function name aliases
 | 
			
		||||
  typedef Optimization::Vsplat   VsplatSIMD;
 | 
			
		||||
  typedef Optimization::Vstore   VstoreSIMD;
 | 
			
		||||
  typedef Optimization::Vset     VsetSIMD;
 | 
			
		||||
  typedef Optimization::Vstream  VstreamSIMD;
 | 
			
		||||
  template <typename S, typename T> using ReduceSIMD = Optimization::Reduce<S,T>;
 | 
			
		||||
// Function name aliases
 | 
			
		||||
typedef Optimization::Vsplat   VsplatSIMD;
 | 
			
		||||
typedef Optimization::Vstore   VstoreSIMD;
 | 
			
		||||
typedef Optimization::Vset     VsetSIMD;
 | 
			
		||||
typedef Optimization::Vstream  VstreamSIMD;
 | 
			
		||||
template <typename S, typename T> using ReduceSIMD = Optimization::Reduce<S,T>;
 | 
			
		||||
 | 
			
		||||
// Arithmetic operations
 | 
			
		||||
typedef Optimization::Sum         SumSIMD;
 | 
			
		||||
typedef Optimization::Sub         SubSIMD;
 | 
			
		||||
typedef Optimization::Div         DivSIMD;
 | 
			
		||||
typedef Optimization::Mult        MultSIMD;
 | 
			
		||||
typedef Optimization::MultComplex MultComplexSIMD;
 | 
			
		||||
typedef Optimization::MultRealPart MultRealPartSIMD;
 | 
			
		||||
typedef Optimization::MaddRealPart MaddRealPartSIMD;
 | 
			
		||||
typedef Optimization::Conj        ConjSIMD;
 | 
			
		||||
typedef Optimization::TimesMinusI TimesMinusISIMD;
 | 
			
		||||
typedef Optimization::TimesI      TimesISIMD;
 | 
			
		||||
 | 
			
		||||
NAMESPACE_END(Grid)
 | 
			
		||||
 | 
			
		||||
  // Arithmetic operations
 | 
			
		||||
  typedef Optimization::Sum         SumSIMD;
 | 
			
		||||
  typedef Optimization::Sub         SubSIMD;
 | 
			
		||||
  typedef Optimization::Div         DivSIMD;
 | 
			
		||||
  typedef Optimization::Mult        MultSIMD;
 | 
			
		||||
  typedef Optimization::MultComplex MultComplexSIMD;
 | 
			
		||||
  typedef Optimization::MultRealPart MultRealPartSIMD;
 | 
			
		||||
  typedef Optimization::MaddRealPart MaddRealPartSIMD;
 | 
			
		||||
  typedef Optimization::Conj        ConjSIMD;
 | 
			
		||||
  typedef Optimization::TimesMinusI TimesMinusISIMD;
 | 
			
		||||
  typedef Optimization::TimesI      TimesISIMD;
 | 
			
		||||
}
 | 
			
		||||
 
 | 
			
		||||
		Reference in New Issue
	
	Block a user