/*
### new version with gribSwapByteOrder_uint16()
icc -g -Wall -O3 -march=native -std=c99 -qopt-report=5 -DTEST_ENCODE encode_array.c
 result on hama2 (icc 16.0.2):
   float:
    orig: val1: 1  val2: 1  val3: 2  valn: 66  time: 1.8731s
unrolled: val1: 1  val2: 1  val3: 2  valn: 66  time: 2.0898s
  double:
    orig: val1: 1  val2: 1  val3: 2  valn: 66  time: 3.68089s
unrolled: val1: 1  val2: 1  val3: 2  valn: 66  time: 4.30798s
     avx: val1: 1  val2: 1  val3: 2  valn: 66  time: 4.23864s

gcc -g -Wall -O3 -march=native -Wa,-q -std=c99 -DTEST_ENCODE encode_array.c
 result on hama2 (gcc 6.1.0):
float:
    orig: val1: 1  val2: 1  val3: 2  valn: 66  time: 2.22871s
unrolled: val1: 1  val2: 1  val3: 2  valn: 66  time: 2.30281s
double:
    orig: val1: 1  val2: 1  val3: 2  valn: 66  time: 4.2669s
unrolled: val1: 1  val2: 1  val3: 2  valn: 66  time: 4.81643s
     avx: val1: 1  val2: 1  val3: 2  valn: 66  time: 3.98415s

###
icc -g -Wall -O3 -march=native -std=c99 -qopt-report=5 -DTEST_ENCODE encode_array.c
 result on hama2 (icc 16.0.0):
   float:
    orig: val1: 1  val2: 1  val3: 2  valn: 66  time: 9.10691s
unrolled: val1: 1  val2: 1  val3: 2  valn: 66  time: 8.63584s
  double:
    orig: val1: 1  val2: 1  val3: 2  valn: 66  time: 13.5768s
unrolled: val1: 1  val2: 1  val3: 2  valn: 66  time: 9.17742s
     avx: val1: 1  val2: 1  val3: 2  valn: 66  time: 3.9488s

gcc -g -Wall -O3 -std=c99 -DTEST_ENCODE encode_array.c
 result on hama2 (gcc 5.2.0):
   float:
    orig: val1: 1  val2: 1  val3: 2  valn: 66  time: 5.32775s
unrolled: val1: 1  val2: 1  val3: 2  valn: 66  time: 7.87125s
  double:
    orig: val1: 1  val2: 1  val3: 2  valn: 66  time: 7.85873s
unrolled: val1: 1  val2: 1  val3: 2  valn: 66  time: 12.9979s

###
gcc -g -Wall -O3 -march=native -std=c99 -DTEST_ENCODE encode_array.c
 result on bailung (gcc 4.7):
  orig    : val1: 1  val2: 1  val3: 2  valn: 66  time: 8.4166s
  sse41   : val1: 1  val2: 1  val3: 2  valn: 66  time: 7.1522s

gcc -g -Wall -O3 -march=native -std=c99 -DTEST_ENCODE encode_array.c
 result on thunder5 (gcc 4.7):
  orig    : val1: 1  val2: 1  val3: 2  valn: 66  time: 6.21976s
  avx     : val1: 1  val2: 1  val3: 2  valn: 66  time: 4.54485s

icc -g -Wall -O3 -march=native -std=c99 -vec-report=1 -DTEST_ENCODE encode_array.c
 result on thunder5 (icc 13.2):
  orig    : val1: 1  val2: 1  val3: 2  valn: 66  time: 14.6279s
  avx     : val1: 1  val2: 1  val3: 2  valn: 66  time:  4.9776s

xlc_r -g -O3 -qhot -q64 -qarch=auto -qtune=auto -qreport -DTEST_ENCODE encode_array.c
 result on blizzard (xlc 12):
  orig    : val1: 1  val2: 1  val3: 2  valn: 66  time: 132.25s
  unrolled: val1: 1  val2: 1  val3: 2  valn: 66  time:  27.202s
  orig    : val1: 1  val2: 1  val3: 2  valn: 66  time: 106.627s  // without -qhot
  unrolled: val1: 1  val2: 1  val3: 2  valn: 66  time:  39.929s  // without -qhot
*/
#ifdef _ARCH_PWR6
#pragma options nostrict
#endif

#ifdef TEST_ENCODE
#include <stdio.h>
#include <stdlib.h>
#define GRIBPACK unsigned char

#ifndef HOST_ENDIANNESS
#ifdef __cplusplus
static const uint32_t HOST_ENDIANNESS_temp[1] = { UINT32_C(0x00030201) };
#define HOST_ENDIANNESS (((const unsigned char *) HOST_ENDIANNESS_temp)[0])
#else
#define HOST_ENDIANNESS (((const unsigned char *) &(const uint32_t[1]){ UINT32_C(0x00030201) })[0])
#endif
#endif

#define IS_BIGENDIAN() (HOST_ENDIANNESS == 0)
#define Error(x, y)
#endif

//#undef _GET_X86_COUNTER
//#undef _GET_MACH_COUNTER
//#undef _GET_IBM_COUNTER
//#undef _ARCH_PWR6

#if defined _GET_IBM_COUNTER
#include <libhpc.h>
#elif defined _GET_X86_COUNTER
#include <x86intrin.h>
#elif defined _GET_MACH_COUNTER
#include <mach/mach_time.h>
#endif

#include <stdint.h>
#include <math.h>

#ifndef DISABLE_SIMD
#if defined(__GNUC__) && (__GNUC__ >= 4)
#elif defined(__ICC) && (__ICC >= 1100)
#elif defined(__clang__)
#else
#define DISABLE_SIMD
#endif
#endif

#ifdef DISABLE_SIMD
#define DISABLE_SIMD_ENCODE
#endif

//#define DISABLE_SIMD_ENCODE

#ifdef DISABLE_SIMD_ENCODE
#ifdef ENABLE_AVX
#define _ENABLE_AVX
#endif
#ifdef ENABLE_SSE4_1
#define _ENABLE_SSE4_1
#endif
#endif

#ifndef DISABLE_SIMD_ENCODE
#ifdef __AVX__
#define _ENABLE_AVX
#endif
#ifdef __SSE4_1__
#define _ENABLE_SSE4_1
#endif
#endif

#if defined _ENABLE_AVX
#include <immintrin.h>
#elif defined _ENABLE_SSE4_1
#include <smmintrin.h>
#endif

#if defined _ENABLE_AVX

static void avx_encode_array_2byte_double(size_t datasize, unsigned char *restrict lGrib, const double *restrict data, double zref,
                                          double factor, size_t *gz) __attribute__((optimize(2)));
static void
avx_encode_array_2byte_double(size_t datasize, unsigned char *restrict lGrib, const double *restrict data, double zref,
                              double factor, size_t *gz)
{
  const double *dval = data;
  __m128i *sgrib = (__m128i *) (lGrib + (*gz));

  const __m128i swap = _mm_set_epi8(14, 15, 12, 13, 10, 11, 8, 9, 6, 7, 4, 5, 2, 3, 0, 1);

  const __m256d c0 = _mm256_set1_pd(zref);
  const __m256d c1 = _mm256_set1_pd(factor);
  const __m256d c2 = _mm256_set1_pd(0.5);

  __m256d d0, d3, d2, d1;
  __m128i i0, i1, i2, i3;
  __m128i s0, s1;

  size_t residual = datasize % 16;

  size_t i;
  for (i = 0; i < (datasize - residual); i += 16)
    {
      (void) _mm_prefetch((const char *) (dval + 8), _MM_HINT_NTA);
      //_____________________________________________________________________________

      d0 = _mm256_loadu_pd(dval);
      d0 = _mm256_sub_pd(d0, c0);
      d0 = _mm256_mul_pd(d0, c1);
      d0 = _mm256_add_pd(d0, c2);

      i0 = _mm256_cvttpd_epi32(d0);

      //_____________________________________________________________________________

      d1 = _mm256_loadu_pd(dval + 4);
      d1 = _mm256_sub_pd(d1, c0);
      d1 = _mm256_mul_pd(d1, c1);
      d1 = _mm256_add_pd(d1, c2);

      i1 = _mm256_cvttpd_epi32(d1);

      //_____________________________________________________________________________

      s0 = _mm_packus_epi32(i0, i1);
      s0 = _mm_shuffle_epi8(s0, swap);
      (void) _mm_storeu_si128(sgrib, s0);

      //_____________________________________________________________________________

      (void) _mm_prefetch((const char *) (dval + 16), _MM_HINT_NTA);

      //_____________________________________________________________________________

      d2 = _mm256_loadu_pd(dval + 8);
      d2 = _mm256_sub_pd(d2, c0);
      d2 = _mm256_mul_pd(d2, c1);
      d2 = _mm256_add_pd(d2, c2);

      i2 = _mm256_cvttpd_epi32(d2);

      //_____________________________________________________________________________

      d3 = _mm256_loadu_pd(dval + 12);
      d3 = _mm256_sub_pd(d3, c0);
      d3 = _mm256_mul_pd(d3, c1);
      d3 = _mm256_add_pd(d3, c2);

      i3 = _mm256_cvttpd_epi32(d3);

      //_____________________________________________________________________________

      s1 = _mm_packus_epi32(i2, i3);
      s1 = _mm_shuffle_epi8(s1, swap);
      (void) _mm_storeu_si128(sgrib + 1, s1);

      //_____________________________________________________________________________

      dval += 16;
      sgrib += 2;
    }

  if (i != datasize)
    {
      uint16_t ui16;
      for (size_t j = i; j < datasize; ++j)
        {
          ui16 = (uint16_t) ((data[j] - zref) * factor + 0.5);
          lGrib[*gz + 2 * j] = ui16 >> 8;
          lGrib[*gz + 2 * j + 1] = ui16;
        }
    }

  *gz += 2 * datasize;

  return;
}

#define grib_encode_array_2byte_double avx_encode_array_2byte_double

#elif defined _ENABLE_SSE4_1

static void
sse41_encode_array_2byte_double(size_t datasize, unsigned char *restrict lGrib, const double *restrict data, double zref,
                                double factor, size_t *gz)
{
  const double *dval = data;
  __m128i *sgrib = (__m128i *) (lGrib + (*gz));

  const __m128i swap = _mm_set_epi8(14, 15, 12, 13, 10, 11, 8, 9, 6, 7, 4, 5, 2, 3, 0, 1);

  const __m128d c0 = _mm_set1_pd(zref);
  const __m128d c1 = _mm_set1_pd(factor);
  const __m128d c2 = _mm_set1_pd(0.5);

  __m128d d0, d4, d3, d2, d1;
  __m128i i0, i1, i2, i3, i4;
  __m128i s0, s1;

  size_t residual = datasize % 16;

  size_t i;
  for (i = 0; i < (datasize - residual); i += 16)
    {
      (void) _mm_prefetch((const char *) (dval + 8), _MM_HINT_NTA);
      //_____________________________________________________________________________

      d0 = _mm_loadu_pd(dval);
      d0 = _mm_sub_pd(d0, c0);
      d0 = _mm_mul_pd(d0, c1);
      d0 = _mm_add_pd(d0, c2);

      d4 = _mm_loadu_pd(dval + 2);
      d4 = _mm_sub_pd(d4, c0);
      d4 = _mm_mul_pd(d4, c1);
      d4 = _mm_add_pd(d4, c2);

      i0 = _mm_cvttpd_epi32(d0);
      i4 = _mm_cvttpd_epi32(d4);
      i0 = _mm_unpacklo_epi64(i0, i4);

      //_____________________________________________________________________________

      d1 = _mm_loadu_pd(dval + 4);
      d1 = _mm_sub_pd(d1, c0);
      d1 = _mm_mul_pd(d1, c1);
      d1 = _mm_add_pd(d1, c2);

      d4 = _mm_loadu_pd(dval + 6);
      d4 = _mm_sub_pd(d4, c0);
      d4 = _mm_mul_pd(d4, c1);
      d4 = _mm_add_pd(d4, c2);

      i1 = _mm_cvttpd_epi32(d1);
      i4 = _mm_cvttpd_epi32(d4);
      i1 = _mm_unpacklo_epi64(i1, i4);

      //_____________________________________________________________________________

      s0 = _mm_packus_epi32(i0, i1);
      s0 = _mm_shuffle_epi8(s0, swap);
      (void) _mm_storeu_si128(sgrib, s0);

      //_____________________________________________________________________________

      (void) _mm_prefetch((const char *) (dval + 16), _MM_HINT_NTA);

      //_____________________________________________________________________________

      d2 = _mm_loadu_pd(dval + 8);
      d2 = _mm_sub_pd(d2, c0);
      d2 = _mm_mul_pd(d2, c1);
      d2 = _mm_add_pd(d2, c2);

      d4 = _mm_loadu_pd(dval + 10);
      d4 = _mm_sub_pd(d4, c0);
      d4 = _mm_mul_pd(d4, c1);
      d4 = _mm_add_pd(d4, c2);

      i2 = _mm_cvttpd_epi32(d2);
      i4 = _mm_cvttpd_epi32(d4);
      i2 = _mm_unpacklo_epi64(i2, i4);

      //_____________________________________________________________________________

      d3 = _mm_loadu_pd(dval + 12);
      d3 = _mm_sub_pd(d3, c0);
      d3 = _mm_mul_pd(d3, c1);
      d3 = _mm_add_pd(d3, c2);

      d4 = _mm_loadu_pd(dval + 14);
      d4 = _mm_sub_pd(d4, c0);
      d4 = _mm_mul_pd(d4, c1);
      d4 = _mm_add_pd(d4, c2);

      i3 = _mm_cvttpd_epi32(d3);
      i4 = _mm_cvttpd_epi32(d4);
      i3 = _mm_unpacklo_epi64(i3, i4);

      //_____________________________________________________________________________

      s1 = _mm_packus_epi32(i2, i3);
      s1 = _mm_shuffle_epi8(s1, swap);
      (void) _mm_storeu_si128(sgrib + 1, s1);

      //_____________________________________________________________________________

      dval += 16;
      sgrib += 2;
    }

  if (i != datasize)
    {
      uint16_t ui16;
      for (size_t j = i; j < datasize; ++j)
        {
          ui16 = (uint16_t) ((data[j] - zref) * factor + 0.5);
          lGrib[*gz + 2 * j] = ui16 >> 8;
          lGrib[*gz + 2 * j + 1] = ui16;
        }
    }

  *gz += 2 * datasize;

  return;
}

#define grib_encode_array_2byte_double sse41_encode_array_2byte_double

#else

#define grib_encode_array_2byte_double encode_array_2byte_double

#endif  // SIMD variants

#ifdef TEST_ENCODE

// clang-format off
#define CAT(X,Y)      X##_##Y
#define TEMPLATE(X,Y) CAT(X,Y)

#ifdef T
#undef T
#endif
#define T double
#include "encode_array.tmpl"

#ifdef T
#undef T
#endif
#define T float
#include "encode_array.tmpl"
// clang-format on

#include <sys/time.h>

static double
dtime()
{
  double tseconds = 0.0;
  struct timeval mytime;
  gettimeofday(&mytime, NULL);
  tseconds = (double) (mytime.tv_sec + (double) mytime.tv_usec * 1.0e-6);
  return (tseconds);
}

static void
pout(char *name, int s, unsigned char *lgrib, long datasize, double tt)
{
  printf("%8s: val1: %d  val2: %d  val3: %d  valn: %d  time: %gs\n", name, (int) lgrib[s * 1 + 1], (int) lgrib[s * 2 + 1],
         (int) lgrib[s * 3 + 1], (int) lgrib[2 * datasize - 1], tt);
}

int
main(void)
{
  enum
  {
    datasize = 1000000,
    NRUN = 10000,
  };

  double t_begin, t_end;

  float *dataf = (float *) malloc(datasize * sizeof(float));
  double *data = (double *) malloc(datasize * sizeof(double));
  unsigned char *lgrib = (unsigned char *) malloc(2 * datasize * sizeof(unsigned char));

  for (long i = 0; i < datasize; ++i) dataf[i] = (float) (-datasize / 2 + i);
  for (long i = 0; i < datasize; ++i) data[i] = (double) (-datasize / 2 + i);

  int PackStart = 0;
  int nbpv = 16;
  double zref = data[0];
  size_t z;
  double factor = 0.00390625;
  int s = 256;

  if (0)
    {
      encode_array_float(0, 0, 0, NULL, NULL, 0, 0, NULL);
      encode_array_double(0, 0, 0, NULL, NULL, 0, 0, NULL);
    }

#if defined(__ICC)
  printf("icc\n");
#elif defined(__clang__)
  printf("clang\n");
#elif defined(__GNUC__)
  printf("gcc\n");
#endif

  printf("float:\n");

  t_begin = dtime();
  for (int i = 0; i < NRUN; ++i)
    {
      z = 0;
      encode_array_2byte_float(datasize, lgrib, dataf, (float) zref, (float) factor, &z);
    }
  t_end = dtime();
  pout("orig", s, lgrib, datasize, t_end - t_begin);

  t_begin = dtime();
  for (int i = 0; i < NRUN; ++i)
    {
      z = 0;
      encode_array_unrolled_float(nbpv, PackStart, datasize, lgrib, dataf, (float) zref, (float) factor, &z);
    }
  t_end = dtime();
  pout("unrolled", s, lgrib, datasize, t_end - t_begin);

  printf("double:\n");

  t_begin = dtime();
  for (int i = 0; i < NRUN; ++i)
    {
      z = 0;
      encode_array_2byte_double(datasize, lgrib, data, zref, factor, &z);
    }
  t_end = dtime();
  pout("orig", s, lgrib, datasize, t_end - t_begin);

  t_begin = dtime();
  for (int i = 0; i < NRUN; ++i)
    {
      z = 0;
      encode_array_unrolled_double(nbpv, PackStart, datasize, lgrib, data, zref, factor, &z);
    }
  t_end = dtime();
  pout("unrolled", s, lgrib, datasize, t_end - t_begin);

#if defined _ENABLE_AVX
  t_begin = dtime();
  for (int i = 0; i < NRUN; ++i)
    {
      z = 0;
      avx_encode_array_2byte_double(datasize, lgrib, data, zref, factor, &z);
    }
  t_end = dtime();
  pout("avx", s, lgrib, datasize, t_end - t_begin);
#elif defined _ENABLE_SSE4_1
  t_begin = dtime();
  for (int i = 0; i < NRUN; ++i)
    {
      z = 0;
      sse41_encode_array_2byte_double(datasize, lgrib, data, zref, factor, &z);
    }
  t_end = dtime();
  pout("sse41", s, lgrib, datasize, t_end - t_begin);
#endif

  return 0;
}
#endif  // TEST_ENCODE

#undef DISABLE_SIMD_ENCODE
#undef _ENABLE_AVX
#undef _ENABLE_SSE4_1
