/*
icc -g -Wall -O3 -march=native -std=c99 -qopt-report=5 -DTEST_MINMAXVAL -qopenmp -DOMP_SIMD minmax_val.c
 result on hama2 (icc 16.0.0):
     float:
minmax_val: fmin: -500000  fmax: 499999  time:   1.22s
simd      : fmin: -500000  fmax: 499999  time:   1.20s
    double:
minmax_val: fmin: -500000  fmax: 499999  time:   2.86s
orig      : fmin: -500000  fmax: 499999  time:   2.74s
simd      : fmin: -500000  fmax: 499999  time:   2.70s
avx       : fmin: -500000  fmax: 499999  time:   2.99s

gcc -g -Wall -O3 -march=native -std=c99 -DTEST_MINMAXVAL -fopenmp -DOMP_SIMD -Wa,-q minmax_val.c
 result on thunder5 (gcc 6.1.0):
float:
minmax_val: fmin: -500000  fmax: 499999  time:   8.25s
  simd    : fmin: -500000  fmax: 499999  time:   1.24s
double:
minmax_val: fmin: -500000  fmax: 499999  time:   2.73s
  orig    : fmin: -500000  fmax: 499999  time:   9.24s
  simd    : fmin: -500000  fmax: 499999  time:   2.78s
  avx     : fmin: -500000  fmax: 499999  time:   2.90s

gcc -g -Wall -O3 -march=native -std=c99 -DTEST_MINMAXVAL minmax_val.c
 result on bailung (gcc 4.8.2):
  orig    : fmin: -500000  fmax: 499999  time:   4.82s
  sse2    : fmin: -500000  fmax: 499999  time:   4.83s

gcc -g -Wall -O3 -march=native -std=c99 -DTEST_MINMAXVAL -fopenmp -DOMP_SIMD -Wa,-q minmax_val.c
 result on thunder5 (gcc 4.8.2):
  orig    : fmin: -500000  fmax: 499999  time:   3.10s
  simd    : fmin: -500000  fmax: 499999  time:   3.10s # omp simd in gcc 4.9
  avx     : fmin: -500000  fmax: 499999  time:   2.84s

icc -g -Wall -O3 -march=native -std=c99 -qopt-report=5 -DTEST_MINMAXVAL -openmp -DOMP_SIMD minmax_val.c
 result on thunder5 (icc 14.0.2):
  orig    : fmin: -500000  fmax: 499999  time:   2.83s
  simd    : fmin: -500000  fmax: 499999  time:   2.83s
  avx     : fmin: -500000  fmax: 499999  time:   2.92s

xlc_r -g -O3 -qhot -q64 -qarch=auto -qtune=auto -qreport -DTEST_MINMAXVAL minmax_val.c
 result on blizzard (xlc 12):
  orig    : fmin: -500000  fmax: 499999  time:   7.26s
  pwr6u6  : fmin: -500000  fmax: 499999  time:   5.92s
*/
#ifdef _ARCH_PWR6
#pragma options nostrict
#endif

#ifdef OMP_SIMD
#include <omp.h>
#endif

#include <stdlib.h>

//#undef _GET_X86_COUNTER
//#undef _GET_IBM_COUNTER
//#undef _GET_MACH_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

#if defined(__GNUC__) && !defined(__ICC) && !defined(__clang__)
#if (__GNUC__ >= 4) && (__GNUC_MINOR__ >= 4)
#define GNUC_PUSH_POP
#endif
#endif

#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_MINMAXVAL
#endif

#ifndef TEST_MINMAXVAL
#define DISABLE_SIMD_MINMAXVAL
#endif

#ifdef DISABLE_SIMD_MINMAXVAL
#ifdef ENABLE_AVX
#define _ENABLE_AVX
#endif
#ifdef ENABLE_SSE2
#define _ENABLE_SSE2
#endif
#endif

#ifndef DISABLE_SIMD_MINMAXVAL
#ifdef __AVX__
#define _ENABLE_AVX
#endif
#ifdef __SSE2__
#define _ENABLE_SSE2
#endif
#endif

#include <float.h>
#include <stdint.h>
#include <inttypes.h>

#if defined(_ENABLE_AVX)
#include <immintrin.h>
#elif defined(_ENABLE_SSE2)
#include <emmintrin.h>
#endif

#if defined(_ENABLE_AVX)

static void
avx_minmax_val_double(const double *restrict buf, size_t nframes, double *min, double *max)
{
  double fmin[4], fmax[4];
  __m256d current_max, current_min, work;

  // load max and min values into all four slots of the YMM registers
  current_min = _mm256_set1_pd(*min);
  current_max = _mm256_set1_pd(*max);

  // Work input until "buf" reaches 32 byte alignment
  while (((unsigned long) buf) % 32 != 0 && nframes > 0)
    {

      // Load the next double into the work buffer
      work = _mm256_set1_pd(*buf);
      current_min = _mm256_min_pd(current_min, work);
      current_max = _mm256_max_pd(current_max, work);
      buf++;
      nframes--;
    }

  while (nframes >= 16)
    {

      (void) _mm_prefetch((const char *) (buf + 8), _MM_HINT_NTA);

      work = _mm256_load_pd(buf);
      current_min = _mm256_min_pd(current_min, work);
      current_max = _mm256_max_pd(current_max, work);
      buf += 4;

      work = _mm256_load_pd(buf);
      current_min = _mm256_min_pd(current_min, work);
      current_max = _mm256_max_pd(current_max, work);
      buf += 4;

      (void) _mm_prefetch((const char *) (buf + 8), _MM_HINT_NTA);

      work = _mm256_load_pd(buf);
      current_min = _mm256_min_pd(current_min, work);
      current_max = _mm256_max_pd(current_max, work);
      buf += 4;

      work = _mm256_load_pd(buf);
      current_min = _mm256_min_pd(current_min, work);
      current_max = _mm256_max_pd(current_max, work);
      buf += 4;
      nframes -= 16;
    }

  // work through aligned buffers
  while (nframes >= 4)
    {
      work = _mm256_load_pd(buf);
      current_min = _mm256_min_pd(current_min, work);
      current_max = _mm256_max_pd(current_max, work);
      buf += 4;
      nframes -= 4;
    }

  // work through the remainung values
  while (nframes > 0)
    {
      work = _mm256_set1_pd(*buf);
      current_min = _mm256_min_pd(current_min, work);
      current_max = _mm256_max_pd(current_max, work);
      buf++;
      nframes--;
    }

  // find min & max value through shuffle tricks

  work = current_min;
  work = _mm256_shuffle_pd(work, work, 5);
  work = _mm256_min_pd(work, current_min);
  current_min = work;
  work = _mm256_permute2f128_pd(work, work, 1);
  work = _mm256_min_pd(work, current_min);
  _mm256_storeu_pd(fmin, work);

  work = current_max;
  work = current_max;
  work = _mm256_shuffle_pd(work, work, 5);
  work = _mm256_max_pd(work, current_max);
  current_max = work;
  work = _mm256_permute2f128_pd(work, work, 1);
  work = _mm256_max_pd(work, current_max);
  _mm256_storeu_pd(fmax, work);

  *min = fmin[0];
  *max = fmax[0];

  return;
}

#elif defined(_ENABLE_SSE2)

static void
sse2_minmax_val_double(const double *restrict buf, size_t nframes, double *min, double *max)
{
  __m128d current_max, current_min, work;

  // load starting max and min values into all slots of the XMM registers
  current_min = _mm_set1_pd(*min);
  current_max = _mm_set1_pd(*max);

  // work on input until buf reaches 16 byte alignment
  while (((unsigned long) buf) % 16 != 0 && nframes > 0)
    {

      // load one double and replicate
      work = _mm_set1_pd(*buf);
      current_min = _mm_min_pd(current_min, work);
      current_max = _mm_max_pd(current_max, work);
      buf++;
      nframes--;
    }

  while (nframes >= 8)
    {
      // use 64 byte prefetch for double octetts
      // __builtin_prefetch(buf+64,0,0); // for GCC 4.3.2 +

      work = _mm_load_pd(buf);
      current_min = _mm_min_pd(current_min, work);
      current_max = _mm_max_pd(current_max, work);
      buf += 2;
      work = _mm_load_pd(buf);
      current_min = _mm_min_pd(current_min, work);
      current_max = _mm_max_pd(current_max, work);
      buf += 2;
      work = _mm_load_pd(buf);
      current_min = _mm_min_pd(current_min, work);
      current_max = _mm_max_pd(current_max, work);
      buf += 2;
      work = _mm_load_pd(buf);
      current_min = _mm_min_pd(current_min, work);
      current_max = _mm_max_pd(current_max, work);
      buf += 2;
      nframes -= 8;
    }

  // work through smaller chunks of aligned buffers without prefetching
  while (nframes >= 2)
    {
      work = _mm_load_pd(buf);
      current_min = _mm_min_pd(current_min, work);
      current_max = _mm_max_pd(current_max, work);
      buf += 2;
      nframes -= 2;
    }

  // work through the remaining value
  while (nframes > 0)
    {
      // load the last double and replicate
      work = _mm_set1_pd(*buf);
      current_min = _mm_min_pd(current_min, work);
      current_max = _mm_max_pd(current_max, work);
      buf++;
      nframes--;
    }

  // find final min and max value through shuffle tricks
  work = current_min;
  work = _mm_shuffle_pd(work, work, _MM_SHUFFLE2(0, 1));
  work = _mm_min_pd(work, current_min);
  _mm_store_sd(min, work);
  work = current_max;
  work = _mm_shuffle_pd(work, work, _MM_SHUFFLE2(0, 1));
  work = _mm_max_pd(work, current_max);
  _mm_store_sd(max, work);

  return;
}

#endif  // SIMD

#if defined(_ARCH_PWR6)
static void
pwr6_minmax_val_double_unrolled6(const double *restrict data, size_t datasize, double *fmin, double *fmax)
{
#define __UNROLL_DEPTH_1 6

  // to allow pipelining we have to unroll

  {
    size_t residual = datasize % __UNROLL_DEPTH_1;
    size_t ofs = datasize - residual;
    double dmin[__UNROLL_DEPTH_1];
    double dmax[__UNROLL_DEPTH_1];

    for (size_t j = 0; j < __UNROLL_DEPTH_1; ++j)
      {
        dmin[j] = data[0];
        dmax[j] = data[0];
      }

    for (size_t i = 0; i < datasize - residual; i += __UNROLL_DEPTH_1)
      {
        for (size_t j = 0; j < __UNROLL_DEPTH_1; ++j)
          {
            dmin[j] = __fsel(dmin[j] - data[i + j], data[i + j], dmin[j]);
            dmax[j] = __fsel(data[i + j] - dmax[j], data[i + j], dmax[j]);
          }
      }

    for (size_t j = 0; j < residual; ++j)
      {
        dmin[j] = __fsel(dmin[j] - data[ofs + j], data[ofs + j], dmin[j]);
        dmax[j] = __fsel(data[ofs + j] - dmax[j], data[ofs + j], dmax[j]);
      }

    for (size_t j = 0; j < __UNROLL_DEPTH_1; ++j)
      {
        *fmin = __fsel(*fmin - dmin[j], dmin[j], *fmin);
        *fmax = __fsel(dmax[j] - *fmax, dmax[j], *fmax);
      }
  }
#undef __UNROLL_DEPTH_1
}
#endif

// clang-format off
#if defined(TEST_MINMAXVAL) && defined(__GNUC__)
static void minmax_val_double_orig(const double *restrict data, size_t datasize, double *fmin, double *fmax) __attribute__((noinline));
static void minmax_val_double_simd(const double *restrict data, size_t datasize, double *fmin, double *fmax) __attribute__((noinline));
static void minmax_val_double_omp(const double *restrict data, size_t datasize, double *fmin, double *fmax) __attribute__((noinline));
static void minmax_val_float(const float *restrict data, long datasize, float *fmin, float *fmax) __attribute__((noinline));
static void minmax_val_float_simd(const float *restrict data, size_t datasize, float *fmin, float *fmax) __attribute__((noinline));
#endif
// clang-format on

#if defined(GNUC_PUSH_POP) && defined __OPTIMIZE__
#pragma GCC push_options
#pragma GCC optimize("O3", "fast-math")
#endif
static void
minmax_val_double_orig(const double *restrict data, size_t datasize, double *fmin, double *fmax)
{
  double dmin = *fmin, dmax = *fmax;

#if defined(CRAY)
#pragma _CRI ivdep
#elif defined(SX)
#pragma vdir nodep
#elif defined(__uxp__)
#pragma loop novrec
#elif defined(__ICC)
#pragma ivdep
#endif
  for (size_t i = 0; i < datasize; ++i)
    {
      dmin = (dmin < data[i]) ? dmin : data[i];
      dmax = (dmax > data[i]) ? dmax : data[i];
    }

  *fmin = dmin;
  *fmax = dmax;
}

static void
minmax_val_float(const float *restrict data, long idatasize, float *fmin, float *fmax)
{
  size_t datasize = (size_t) idatasize;
  float dmin = *fmin, dmax = *fmax;

#if defined(CRAY)
#pragma _CRI ivdep
#elif defined(SX)
#pragma vdir nodep
#elif defined(__uxp__)
#pragma loop novrec
#elif defined(__ICC)
#pragma ivdep
#endif
  for (size_t i = 0; i < datasize; ++i)
    {
      dmin = (dmin < data[i]) ? dmin : data[i];
      dmax = (dmax > data[i]) ? dmax : data[i];
    }

  *fmin = dmin;
  *fmax = dmax;
}
#if defined(GNUC_PUSH_POP) && defined __OPTIMIZE__
#pragma GCC pop_options
#endif

// TEST
#if defined(OMP_SIMD)

#if defined(GNUC_PUSH_POP) && defined __OPTIMIZE__
#pragma GCC push_options
#pragma GCC optimize("O3", "fast-math")
#endif
static void
minmax_val_double_omp(const double *restrict data, size_t datasize, double *fmin, double *fmax)
{
  double dmin = *fmin, dmax = *fmax;

#if defined(_OPENMP)
#pragma omp parallel for simd reduction(min : dmin) reduction(max : dmax)
#endif
  for (size_t i = 0; i < datasize; ++i)
    {
      dmin = (dmin < data[i]) ? dmin : data[i];
      dmax = (dmax > data[i]) ? dmax : data[i];
    }

  *fmin = dmin;
  *fmax = dmax;
}

static void
minmax_val_double_simd(const double *restrict data, size_t datasize, double *fmin, double *fmax)
{
  double dmin = *fmin, dmax = *fmax;

#ifdef _OPENMP
#pragma omp simd reduction(min : dmin) reduction(max : dmax)
#endif
  for (size_t i = 0; i < datasize; ++i)
    {
      dmin = (dmin < data[i]) ? dmin : data[i];
      dmax = (dmax > data[i]) ? dmax : data[i];
    }

  *fmin = dmin;
  *fmax = dmax;
}

static void
minmax_val_float_simd(const float *restrict data, size_t datasize, float *fmin, float *fmax)
{
  float dmin = *fmin, dmax = *fmax;

#if defined(_OPENMP)
#pragma omp simd reduction(min : dmin) reduction(max : dmax)
#endif
  for (size_t i = 0; i < datasize; ++i)
    {
      dmin = (dmin < data[i]) ? dmin : data[i];
      dmax = (dmax > data[i]) ? dmax : data[i];
    }

  *fmin = dmin;
  *fmax = dmax;
}
#if defined(GNUC_PUSH_POP) && defined __OPTIMIZE__
#pragma GCC pop_options
#endif
#endif

static void
minmax_val_double(const double *restrict data, long idatasize, double *fmin, double *fmax)
{
#if defined(_GET_X86_COUNTER) || defined(_GET_MACH_COUNTER)
  uint64_t start_minmax, end_minmax;
#endif
  size_t datasize = (size_t) idatasize;

  if (idatasize >= 1)
    ;
  else
    return;

#if defined(_GET_X86_COUNTER)
  start_minmax = _rdtsc();
#endif
#if defined(_GET_MACH_COUNTER)
  start_minmax = mach_absolute_time();
#endif

#if defined(_ENABLE_AVX)

  avx_minmax_val_double(data, datasize, fmin, fmax);

#elif defined(_ENABLE_SSE2)

  sse2_minmax_val_double(data, datasize, fmin, fmax);

#else

#if defined(_ARCH_PWR6)
#define __UNROLL_DEPTH_1 6

    // to allow pipelining we have to unroll

#if defined(_GET_IBM_COUNTER)
  hpmStart(1, "minmax fsel");
#endif

  pwr6_minmax_val_double_unrolled6(data, datasize, fmin, fmax);

#if defined(_GET_IBM_COUNTER)
  hpmStop(1);
#endif

#undef __UNROLL_DEPTH_1

#else  // original loop

#if defined(_GET_IBM_COUNTER)
  hpmStart(1, "minmax base");
#endif

  minmax_val_double_orig(data, datasize, fmin, fmax);

#if defined(_GET_IBM_COUNTER)
  hpmStop(1);
#endif

#endif  // _ARCH_PWR6 && original loop
#endif  // SIMD

#if defined(_GET_X86_COUNTER) || defined(_GET_MACH_COUNTER)
#if defined(_GET_X86_COUNTER)
  end_minmax = _rdtsc();
#endif
#if defined(_GET_MACH_COUNTER)
  end_minmax = mach_absolute_time();
#endif
#if defined(_ENABLE_AVX)
  printf("AVX minmax cycles:: %" PRIu64 "\n", end_minmax - start_minmax);
  fprintf(stderr, "AVX min: %lf max: %lf\n", *fmin, *fmax);
#elif defined(_ENABLE_SSE2)
  printf("SSE2 minmax cycles:: %" PRIu64 "\n", end_minmax - start_minmax);
  fprintf(stderr, "SSE2 min: %lf max: %lf\n", *fmin, *fmax);
#else
  printf("loop minmax cycles:: %" PRIu64 "\n", end_minmax - start_minmax);
  fprintf(stderr, "loop min: %lf max: %lf\n", *fmin, *fmax);
#endif
#endif

  return;
}

#if defined(TEST_MINMAXVAL)

#include <stdio.h>
#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);
}

#define NRUN 10000

int
main(void)
{
  long datasize = 1000000;
  double t_begin, t_end;

  printf("datasize %ld\n", datasize);
#if defined(_OPENMP)
  printf("_OPENMP=%d\n", _OPENMP);
#endif

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

  {
    float fmin, fmax;
    float *data_sp = (float *) malloc(datasize * sizeof(float));

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

    printf("float:\n");

    t_begin = dtime();
    for (int i = 0; i < NRUN; ++i)
      {
        fmin = fmax = data_sp[0];
        minmax_val_float(data_sp, datasize, &fmin, &fmax);
      }
    t_end = dtime();
    printf("minmax_val: fmin: %ld  fmax: %ld  time: %6.2fs\n", (long) fmin, (long) fmax, t_end - t_begin);

#if defined(OMP_SIMD)
    t_begin = dtime();
    for (int i = 0; i < NRUN; ++i)
      {
        fmin = fmax = data_sp[0];
        minmax_val_float_simd(data_sp, datasize, &fmin, &fmax);
      }
    t_end = dtime();
    printf("simd      : fmin: %ld  fmax: %ld  time: %6.2fs\n", (long) fmin, (long) fmax, t_end - t_begin);
#endif

    free(data_sp);
  }

  {
    double fmin, fmax;
    double *data_dp = (double *) malloc(datasize * sizeof(double));

    // for (long i = datasize-1; i >= 0; --i) data[i] = (double) (-datasize/2 + i);
    for (long i = 0; i < datasize / 2; ++i) data_dp[i] = (double) (i);
    for (long i = datasize / 2; i < datasize; ++i) data_dp[i] = (double) (-datasize + i);

    printf("double:\n");

    t_begin = dtime();
    for (int i = 0; i < NRUN; ++i)
      {
        fmin = fmax = data_dp[0];
        minmax_val_double(data_dp, datasize, &fmin, &fmax);
      }
    t_end = dtime();
    printf("minmax_val: fmin: %ld  fmax: %ld  time: %6.2fs\n", (long) fmin, (long) fmax, t_end - t_begin);

    t_begin = dtime();
    for (int i = 0; i < NRUN; ++i)
      {
        fmin = fmax = data_dp[0];
        minmax_val_double_orig(data_dp, datasize, &fmin, &fmax);
      }
    t_end = dtime();
    printf("orig      : fmin: %ld  fmax: %ld  time: %6.2fs\n", (long) fmin, (long) fmax, t_end - t_begin);

#if defined(OMP_SIMD)
    t_begin = dtime();
    for (int i = 0; i < NRUN; ++i)
      {
        fmin = fmax = data_dp[0];
        minmax_val_double_simd(data_dp, datasize, &fmin, &fmax);
      }
    t_end = dtime();
    printf("simd      : fmin: %ld  fmax: %ld  time: %6.2fs\n", (long) fmin, (long) fmax, t_end - t_begin);

    t_begin = dtime();
    for (int i = 0; i < NRUN; ++i)
      {
        fmin = fmax = data_dp[0];
        minmax_val_double_omp(data_dp, datasize, &fmin, &fmax);
      }
    t_end = dtime();
    printf("openmp %d  : fmin: %ld  fmax: %ld  time: %6.2fs\n", omp_get_max_threads(), (long) fmin, (long) fmax, t_end - t_begin);
#endif

#if defined(_ENABLE_AVX)
    t_begin = dtime();
    for (int i = 0; i < NRUN; ++i)
      {
        fmin = fmax = data_dp[0];
        avx_minmax_val_double(data_dp, datasize, &fmin, &fmax);
      }
    t_end = dtime();
    printf("avx       : fmin: %ld  fmax: %ld  time: %6.2fs\n", (long) fmin, (long) fmax, t_end - t_begin);
#elif defined(_ENABLE_SSE2)
    t_begin = dtime();
    for (int i = 0; i < NRUN; ++i)
      {
        fmin = fmax = data_dp[0];
        sse2_minmax_val_double(data_dp, datasize, &fmin, &fmax);
      }
    t_end = dtime();
    printf("sse2      : fmin: %ld  fmax: %ld  time: %6.2fs\n", (long) fmin, (long) fmax, t_end - t_begin);
#endif
#if defined(_ARCH_PWR6)
    t_begin = dtime();
    for (int i = 0; i < NRUN; ++i)
      {
        fmin = fmax = data_dp[0];
        pwr6_minmax_val_double_unrolled6(data_dp, datasize, &fmin, &fmax);
      }
    t_end = dtime();
    printf("pwr6u6  : fmin: %ld  fmax: %ld  time: %6.2fs\n", (long) fmin, (long) fmax, t_end - t_begin);
#endif
    free(data_dp);
  }

  return 0;
}
#endif  // TEST_MINMAXVAL

#undef DISABLE_SIMD_MINMAXVAL
#undef _ENABLE_AVX
#undef _ENABLE_SSE2
#undef GNUC_PUSH_POP
