// Instantiante some functions to be able to inspect the generated
// machine code

#define VML_NODEBUG

#define restrict __restrict__

#include "vecmathlib.h"



namespace vecmathlib {
  
  template<typename realvec_t, int n>
  typename realvec_t::real_t get_elt(realvec_t x)
  {
    return x[n];
  }
  template<typename realvec_t, int n>
  realvec_t set_elt(realvec_t x, typename realvec_t::real_t a)
  {
    return x.set_elt(n, a);
  }
  
  // template realbuiltinvec<float,1> fabs(realbuiltinvec<float,1> x);
  // template realbuiltinvec<float,1> fmin(realbuiltinvec<float,1> x, realbuiltinvec<float,1> y);
  // template intbuiltinvec<float,1> lsr(intbuiltinvec<float,1> x, intbuiltinvec<float,1>::int_t n);
  // template intbuiltinvec<double,1> lsr(intbuiltinvec<double,1> x, intbuiltinvec<double,1>::int_t n);
  // template intbuiltinvec<double,2> lsr(intbuiltinvec<double,2> x, intbuiltinvec<double,2>::int_t n);
  // template intbuiltinvec<double,2> lsr(intbuiltinvec<double,2> x, intbuiltinvec<double,2> n);
  // template realbuiltinvec<float,1> ifthen(realbuiltinvec<float,1>::boolvec_t c, realbuiltinvec<float,1> x, realbuiltinvec<float,1> y);
  // template realbuiltinvec<double,1> ifthen(realbuiltinvec<double,1>::boolvec_t c, realbuiltinvec<double,1> x, realbuiltinvec<double,1> y);
  // template realbuiltinvec<float,4> ifthen(realbuiltinvec<float,4>::boolvec_t c, realbuiltinvec<float,4> x, realbuiltinvec<float,4> y);
  // template realbuiltinvec<double,2> ifthen(realbuiltinvec<double,2>::boolvec_t c, realbuiltinvec<double,2> x, realbuiltinvec<double,2> y);
  
#ifdef VECMATHLIB_HAVE_VEC_FLOAT_1
  template realvec<float,1> round(realvec<float,1> x);
#endif
  
#ifdef VECMATHLIB_HAVE_VEC_FLOAT_8
  template intvec<float,8> popcount(intvec<float,8>);
#endif
  
#ifdef VECMATHLIB_HAVE_VEC_DOUBLE_1
  template realvec<double,1> exp(realvec<double,1> x);
  template realvec<double,1> log(realvec<double,1> x);
  template realvec<double,1> sin(realvec<double,1> x);
  template realvec<double,1> sqrt(realvec<double,1> x);
  template realvec<double,1>::real_t get_elt<realvec<double,1>,0>(realvec<double,1> x);
  template realvec<double,1> set_elt<realvec<double,1>,0>(realvec<double,1> x, realvec<double,1>::real_t a);
#endif
  
#ifdef VECMATHLIB_HAVE_VEC_DOUBLE_2
  template realvec<double,2> exp(realvec<double,2> x);
  template realvec<double,2> log(realvec<double,2> x);
  template realvec<double,2> sin(realvec<double,2> x);
  template realvec<double,2> sqrt(realvec<double,2> x);
  template realvec<double,2>::real_t get_elt<realvec<double,2>,0>(realvec<double,2>);
  template realvec<double,2>::real_t get_elt<realvec<double,2>,1>(realvec<double,2>);
  template realvec<double,2> set_elt<realvec<double,2>,0>(realvec<double,2> x, realvec<double,2>::real_t a);
  template realvec<double,2> set_elt<realvec<double,2>,1>(realvec<double,2> x, realvec<double,2>::real_t a);
#endif
  
#ifdef VECMATHLIB_HAVE_VEC_DOUBLE_4
  template realvec<double,4> exp(realvec<double,4> x);
  template realvec<double,4> log(realvec<double,4> x);
  template realvec<double,4> sin(realvec<double,4> x);
  template realvec<double,4> sqrt(realvec<double,4> x);
  template realvec<double,4>::real_t get_elt<realvec<double,4>,0>(realvec<double,4>);
  template realvec<double,4>::real_t get_elt<realvec<double,4>,1>(realvec<double,4>);
  template realvec<double,4>::real_t get_elt<realvec<double,4>,2>(realvec<double,4>);
  template realvec<double,4>::real_t get_elt<realvec<double,4>,3>(realvec<double,4>);
  template realvec<double,4> set_elt<realvec<double,4>,0>(realvec<double,4> x, realvec<double,4>::real_t a);
  template realvec<double,4> set_elt<realvec<double,4>,1>(realvec<double,4> x, realvec<double,4>::real_t a);
  template realvec<double,4> set_elt<realvec<double,4>,2>(realvec<double,4> x, realvec<double,4>::real_t a);
  template realvec<double,4> set_elt<realvec<double,4>,3>(realvec<double,4> x, realvec<double,4>::real_t a);
  template intvec<double,4> popcount(intvec<double,4>);
#endif
  
}



// Various tests to detect auto-vectorization features



#include <cassert>
#include <cstdlib>
using namespace std;

using namespace vecmathlib;

#if defined VECMATHLIB_HAVE_VEC_DOUBLE_4
typedef realvec<double,4> realV;
#elif defined VECMATHLIB_HAVE_VEC_DOUBLE_2
typedef realvec<double,2> realV;
#elif defined VECMATHLIB_HAVE_VEC_FLOAT_8
typedef realvec<float,8> realV;
#elif defined VECMATHLIB_HAVE_VEC_FLOAT_4
typedef realvec<float,4> realV;
#elif defined VECMATHLIB_HAVE_VEC_FLOAT_2
typedef realvec<float,2> realV;
#else
#  error "There are no vector types"
#endif

typedef realV::scalar_t real;
const int vecsize = realV::size;



// Simple, naive loop adding two arrays
extern "C"
void loop_add(real* a,
              real* b,
              real* c,
              ptrdiff_t n)
{
  for (ptrdiff_t i=0; i<n; i+=vecsize) {
    realV tmpb = realV::loadu(&b[i]);
    realV tmpc = realV::loadu(&c[i]);
    realV tmpa = tmpb + tmpc;
    storeu(tmpa, &a[i]);
  }
}



// Declare pointers as restrict
extern "C"
void loop_add_restrict(real *restrict a,
                       real *restrict b,
                       real *restrict c,
                       ptrdiff_t n)
{
  for (ptrdiff_t i=0; i<n; i+=vecsize) {
    realV tmpb = realV::loadu(&b[i]);
    realV tmpc = realV::loadu(&c[i]);
    realV tmpa = tmpb + tmpc;
    storeu(tmpa, &a[i]);
  }
}



// Declare pointers as restrict and aligned
extern "C"
void loop_add_aligned(real *restrict a,
                      real *restrict b,
                      real *restrict c,
                      ptrdiff_t n)
{
  for (ptrdiff_t i=0; i<n; i+=vecsize) {
    realV tmpb = realV::loada(&b[i]);
    realV tmpc = realV::loada(&c[i]);
    realV tmpa = tmpb + tmpc;
    storea(tmpa, &a[i]);
  }
}



// Reduction loop
extern "C"
real loop_dot_reduce(real *restrict a,
                     real *restrict b,
                     ptrdiff_t n)
{
  realV sumV = 0.0;
  for (ptrdiff_t i=0; i<n; i+=vecsize) {
    realV tmpa = realV::loada(&a[i]);
    realV tmpb = realV::loada(&b[i]);
    sumV += tmpa * tmpb;
  }
  return sum(sumV);
}



// Loop with a simple if condition (fmax)
extern "C"
void loop_if_simple(real *restrict a,
                    real *restrict b,
                    real *restrict c,
                    ptrdiff_t n)
{
  for (ptrdiff_t i=0; i<n; i+=vecsize) {
    realV tmpb = realV::loada(&b[i]);
    realV tmpc = realV::loada(&c[i]);
    realV tmpa = ifthen(tmpb > tmpc, tmpb, tmpc);
    storea(tmpa, &a[i]);
  }
}



// Loop with a complex if condition (select)
extern "C"
void loop_if(real *restrict a,
             real *restrict b,
             real *restrict c,
             ptrdiff_t n)
{
  for (ptrdiff_t i=0; i<n; i+=vecsize) {
    realV tmpb = realV::loada(&b[i]);
    realV tmpc = realV::loada(&c[i]);
    realV tmpa = ifthen(tmpb > realV(0.0), tmpb * tmpc, realV(1.0));
    storea(tmpa, &a[i]);
  }
}



// Skip ghost points
extern "C"
void loop_add_masked(real *restrict a,
                     real *restrict b,
                     real *restrict c,
                     ptrdiff_t n)
{
  for (realV::mask_t mask(1, n-1, 0); mask; ++mask) {
    ptrdiff_t i = mask.index();
    realV tmpb = realV::loada(&b[i]);
    realV tmpc = realV::loada(&c[i]);
    realV tmpa = tmpb + tmpc;
    storea(tmpa, &a[i], mask);
  }
}