Replace all mdrun rngs with cycle based rng

[alexxy/gromacs.git] / src / gromacs / simd / four_wide_macros.h

/*
 * This file is part of the GROMACS molecular simulation package.
 *
 * Copyright (c) 2012,2013,2014, by the GROMACS development team, led by
 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
 * and including many others, as listed in the AUTHORS file in the
 * top-level source directory and at http://www.gromacs.org.
 *
 * GROMACS is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License
 * as published by the Free Software Foundation; either version 2.1
 * of the License, or (at your option) any later version.
 *
 * GROMACS is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with GROMACS; if not, see
 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
 *
 * If you want to redistribute modifications to GROMACS, please
 * consider that scientific software is very special. Version
 * control is crucial - bugs must be traceable. We will be happy to
 * consider code for inclusion in the official distribution, but
 * derived work must not be called official GROMACS. Details are found
 * in the README & COPYING files - if they are missing, get the
 * official version at http://www.gromacs.org.
 *
 * To help us fund GROMACS development, we humbly ask that you cite
 * the research papers on the package. Check out http://www.gromacs.org.
 */

/* The macros in this file are intended to be used for writing
 * architecture-independent SIMD intrinsics code with a SIMD width of 4.
 * To support a new architecture, adding macros here should be all
 * that is needed.
 *
 * Note that this file is intended only for SIMD operations that require
 * a SIMD width of 4. In general gmx_simd_macros.h provides wider hardware
 * support, more functionality and higher performance, but the SIMD width is
 * not necessarily equal to 4.
 */

#ifdef GMX_SIMD_FOUR_WIDE_MACROS_H
#error "four_wide_macros.h included twice"
#else
#define GMX_SIMD_FOUR_WIDE_MACROS_H


/* The SIMD width here is always 4, since that is the whole point */
#define GMX_SIMD4_WIDTH  4


#if defined GMX_SIMD4_SINGLE || defined GMX_SIMD4_DOUBLE
/* Precision set before inclusion, honour that request */
#else
/* Match precision to the Gromacs real precision */
#ifdef GMX_DOUBLE
#define GMX_SIMD4_DOUBLE
#else
#define GMX_SIMD4_SINGLE
#endif
#endif

#ifdef GMX_SIMD4_DOUBLE
typedef double  gmx_simd4_real;
#endif
#ifdef GMX_SIMD4_SINGLE
typedef float   gmx_simd4_real;
#endif

/* Uncomment the next line, without other SIMD active, for testing plain-C */
/* #define GMX_SIMD4_REFERENCE */
#ifdef GMX_SIMD4_REFERENCE
/* Plain C SIMD reference implementation, also serves as documentation */
#define GMX_HAVE_SIMD4_MACROS

/* Include plain-C reference implementation, also serves as documentation */
#include "four_wide_macros_ref.h"

/* float/double SIMD register type */
#define gmx_simd4_real_t  gmx_simd4_ref_pr

/* boolean SIMD register type */
#define gmx_simd4_bool_t  gmx_simd4_ref_pb

#define gmx_simd4_load_r       gmx_simd4_ref_load_pr
#define gmx_simd4_load_bb_pr    gmx_simd4_ref_load_pr
#define gmx_simd4_set1_r       gmx_simd4_ref_set1_pr
#define gmx_simd4_setzero_r    gmx_simd4_ref_setzero_pr
#define gmx_simd4_store_r      gmx_simd4_ref_store_pr

/* Unaligned load+store are not required,
 * but they can speed up the PME spread+gather operations.
 */
#define GMX_SIMD4_HAVE_UNALIGNED
#ifdef GMX_SIMD4_HAVE_UNALIGNED
#define gmx_simd4_loadu_r      gmx_simd4_ref_load_pr
#define gmx_simd4_storeu_r     gmx_simd4_ref_store_pr
#endif

#define gmx_simd4_add_r        gmx_simd4_ref_add_pr
#define gmx_simd4_sub_r        gmx_simd4_ref_sub_pr
#define gmx_simd4_mul_r        gmx_simd4_ref_mul_pr
/* For the FMA macros below, aim for c=d in code, so FMA3 uses 1 instruction */
#define gmx_simd4_fmadd_r       gmx_simd4_ref_madd_pr
#define gmx_simd4_fnmadd_r      gmx_simd4_ref_nmsub_pr

#define gmx_simd4_dotproduct3_r   gmx_simd4_ref_dotproduct3

#define gmx_simd4_min_r        gmx_simd4_ref_min_pr
#define gmx_simd4_max_r        gmx_simd4_ref_max_pr

#define gmx_simd4_blendzero_r  gmx_simd4_ref_blendzero_pr

/* Comparison */
#define gmx_simd4_cmplt_r      gmx_simd4_ref_cmplt_pr

/* Logical operations on SIMD booleans */
#define gmx_simd4_and_b        gmx_simd4_ref_and_pb
#define gmx_simd4_or_b         gmx_simd4_ref_or_pb

/* Returns a single int (0/1) which tells if any of the 4 booleans is True */
#define gmx_simd4_anytrue_b    gmx_simd4_ref_anytrue_pb

#endif /* GMX_SIMD4_REFERENCE */


/* The same SIMD macros can be translated to SIMD intrinsics (and compiled
 * to instructions for) different SIMD width and float precision.
 *
 * On x86: The gmx_simd4 prefix is replaced by _mm_ or _mm256_ (SSE or AVX).
 * The _pr suffix is replaced by _ps or _pd (for single or double precision).
 * Compiler settings will decide if 128-bit intrinsics will
 * be translated into SSE or AVX instructions.
 */


#ifdef GMX_SIMD_X86_SSE2_OR_HIGHER
/* This is for general x86 SIMD instruction sets that also support SSE2 */

#ifdef GMX_SIMD4_SINGLE
#define GMX_HAVE_SIMD4_MACROS
#endif

#ifdef GMX_SIMD4_DOUBLE
/* Note that here we will use 256-bit SIMD with GMX_SIMD_X86_AVX_128_FMA_OR_HIGHER.
 * This is inconsistent naming wise, but should give the best performance.
 */
#if defined GMX_SIMD_X86_AVX_128_FMA_OR_HIGHER || defined GMX_SIMD_X86_AVX_256_OR_HIGHER
#define GMX_HAVE_SIMD4_MACROS
#endif
#endif

#ifdef GMX_HAVE_SIMD4_MACROS

#if defined GMX_SIMD_X86_AVX_128_FMA_OR_HIGHER || defined GMX_SIMD_X86_AVX_256_OR_HIGHER

#include <immintrin.h>
#ifdef HAVE_X86INTRIN_H
#include <x86intrin.h> /* FMA */
#endif
#ifdef HAVE_INTRIN_H
#include <intrin.h> /* FMA MSVC */
#endif

#else
#ifdef GMX_SIMD_X86_SSE4_1_OR_HIGHER
#include <smmintrin.h>
#else
/* We only have SSE2 */
#include <emmintrin.h>
#endif
#endif

#ifdef GMX_SIMD4_SINGLE

#define gmx_simd4_real_t  __m128

#define gmx_simd4_bool_t  __m128

#define gmx_simd4_load_r       _mm_load_ps
#define gmx_simd4_load_bb_pr    _mm_load_ps
#define gmx_simd4_set1_r       _mm_set1_ps
#define gmx_simd4_setzero_r    _mm_setzero_ps
#define gmx_simd4_store_r      _mm_store_ps

/* Some old AMD processors could have problems with unaligned loads+stores */
#ifndef GMX_FAHCORE
#define GMX_SIMD4_HAVE_UNALIGNED
#endif
#ifdef GMX_SIMD4_HAVE_UNALIGNED
#define gmx_simd4_loadu_r      _mm_loadu_ps
#define gmx_simd4_storeu_r     _mm_storeu_ps
#endif

#define gmx_simd4_add_r        _mm_add_ps
#define gmx_simd4_sub_r        _mm_sub_ps
#define gmx_simd4_mul_r        _mm_mul_ps

#ifdef GMX_SIMD_X86_AVX_128_FMA_OR_HIGHER
#define gmx_simd4_fmadd_r(a, b, c)   _mm_macc_ps(a, b, c)
#define gmx_simd4_fnmadd_r(a, b, c)  _mm_nmacc_ps(a, b, c)
#else
#define gmx_simd4_fmadd_r(a, b, c)   _mm_add_ps(c, _mm_mul_ps(a, b))
#define gmx_simd4_fnmadd_r(a, b, c)  _mm_sub_ps(c, _mm_mul_ps(a, b))
#endif

static inline float gmx_simd4_dotproduct3_r(__m128 a, __m128 b)
#ifdef GMX_SIMD_X86_SSE4_1_OR_HIGHER
{
    float dp;

    /* SSE4.1 dot product of components 0,1,2, stored in component 0 */
    _mm_store_ss(&dp, _mm_dp_ps(a, b, 0x71));

    return dp;
}
#else
{
    float        dp_array[7], *dp;

    /* Generate an aligned pointer */
    dp = (float *)(((size_t)(dp_array+3)) & (~((size_t)15)));

    _mm_store_ps(dp, _mm_mul_ps(a, b));

    return dp[0] + dp[1] + dp[2];
}
#endif

#define gmx_simd4_min_r        _mm_min_ps
#define gmx_simd4_max_r        _mm_max_ps

#define gmx_simd4_blendzero_r  _mm_and_ps

#define gmx_simd4_cmplt_r      _mm_cmplt_ps
#define gmx_simd4_and_b        _mm_and_ps
#define gmx_simd4_or_b         _mm_or_ps

#define gmx_simd4_anytrue_b    _mm_movemask_ps

#endif /* GMX_SIMD4_SINGLE */


#ifdef GMX_SIMD4_DOUBLE

#define gmx_simd4_real_t  __m256d

#define gmx_simd4_bool_t  __m256d

#define gmx_simd4_load_r       _mm256_load_pd
#define gmx_simd4_load_bb_pr    _mm256_load_pd
#define gmx_simd4_set1_r       _mm256_set1_pd
#define gmx_simd4_setzero_r    _mm256_setzero_pd
#define gmx_simd4_store_r      _mm256_store_pd

#define GMX_SIMD4_HAVE_UNALIGNED
#define gmx_simd4_loadu_r      _mm256_loadu_pd
#define gmx_simd4_storeu_r     _mm256_storeu_pd

#define gmx_simd4_add_r        _mm256_add_pd
#define gmx_simd4_sub_r        _mm256_sub_pd
#define gmx_simd4_mul_r        _mm256_mul_pd
#ifdef GMX_SIMD_X86_AVX_128_FMA_OR_HIGHER
#define gmx_simd4_fmadd_r(a, b, c)   _mm256_macc_pd(a, b, c)
#define gmx_simd4_fnmadd_r(a, b, c)  _mm256_nmacc_pd(a, b, c)
#else
#define gmx_simd4_fmadd_r(a, b, c)   _mm256_add_pd(c, _mm256_mul_pd(a, b))
#define gmx_simd4_fnmadd_r(a, b, c)  _mm256_sub_pd(c, _mm256_mul_pd(a, b))
#endif
#define gmx_simd4_min_r        _mm256_min_pd
#define gmx_simd4_max_r        _mm256_max_pd

#define gmx_simd4_blendzero_r  _mm256_and_pd

/* Less-than (we use ordered, non-signaling, but that's not required) */
#define gmx_simd4_cmplt_r(x, y) _mm256_cmp_pd(x, y, 0x11)
#define gmx_simd4_and_b        _mm256_and_pd
#define gmx_simd4_or_b         _mm256_or_pd

#define gmx_simd4_anytrue_b    _mm256_movemask_pd

#endif /* GMX_SIMD4_DOUBLE */


#endif /* GMX_HAVE_SIMD4_MACROS */

#endif /* GMX_SIMD_X86_SSE2_OR_HIGHER */

#ifdef GMX_SIMD_IBM_QPX
/* i.e. BlueGene/Q */

/* This hack works on the compilers that can reach this code. A real
   solution with broader scope will be proposed in master branch. */
#define gmx_always_inline __attribute__((always_inline))

#ifdef GMX_SIMD4_SINGLE
#define GMX_HAVE_SIMD4_MACROS
#endif

typedef vector4double gmx_simd4_real_t;
typedef vector4double gmx_simd4_bool_t;

/* The declarations of vec_ld* use non-const pointers, and IBM
   can't/won't fix this any time soon. So GROMACS has to cast away the
   const-ness of its pointers before loads. Four-wide SIMD loads
   sometimes occur from variables of type real, and sometimes from
   variables of type float (even at double precison), so the correct
   cast cannot be done easily. The correct cast is necessary because
   the resulting type determines the alignment assumption of vec_ld*,
   which is different for float and double. So the loads of
   always-float variables have to be done with a function that does
   the correct cast. Since functions cannot be overloaded by type in
   C, they have to have different names. Thus we have
   gmx_simd4_load_r and gmx_simd4_load_bb_pr.
 */

static gmx_inline gmx_simd4_real_t gmx_always_inline gmx_simd4_load_r(const real *a)
{
#ifdef NDEBUG
    return vec_ld(0, (real *) a);
#else
    return vec_lda(0, (real *) a);
#endif
}

static gmx_inline gmx_simd4_real_t gmx_always_inline gmx_simd4_load_bb_pr(const float *a)
{
#ifdef NDEBUG
    return vec_ld(0, (float *) a);
#else
    return vec_lda(0, (float *) a);
#endif
}

static gmx_inline gmx_simd4_real_t gmx_always_inline gmx_simd4_set1_r(const real a)
{
    return vec_splats(a);
}

static gmx_inline gmx_simd4_real_t gmx_always_inline gmx_simd4_setzero_r()
{
    return vec_splats(0.0);
}

/* TODO this will not yet work, because the function might be passed a
   pointer to a float when running in double precision.
 */
static gmx_inline void gmx_always_inline gmx_simd4_store_r(real *a, gmx_simd4_real_t b)
{
#ifdef NDEBUG
    vec_st(b, 0, a);
#else
    vec_sta(b, 0, a);
#endif
}

static gmx_inline gmx_simd4_real_t gmx_always_inline gmx_simd4_add_r(gmx_simd4_real_t a, gmx_simd4_real_t b)
{
    return vec_add(a, b);
}

static gmx_inline gmx_simd4_real_t gmx_always_inline gmx_simd4_sub_r(gmx_simd4_real_t a, gmx_simd4_real_t b)
{
    return vec_sub(a, b);
}

static gmx_inline gmx_simd4_real_t gmx_always_inline gmx_simd4_mul_r(gmx_simd4_real_t a, gmx_simd4_real_t b)
{
    return vec_mul(a, b);
}

static gmx_inline gmx_simd4_real_t gmx_always_inline gmx_simd4_fmadd_r(gmx_simd4_real_t a, gmx_simd4_real_t b, gmx_simd4_real_t c)
{
    return vec_madd(a, b, c);
}

static gmx_inline gmx_simd4_real_t gmx_always_inline gmx_simd4_fnmadd_r(gmx_simd4_real_t a, gmx_simd4_real_t b, gmx_simd4_real_t c)
{
    return vec_nmsub(a, b, c);
}

static gmx_inline gmx_simd4_real_t gmx_always_inline gmx_simd4_min_r(gmx_simd4_real_t a, gmx_simd4_real_t b)
{
    /* Implemented the same way as max, but with the subtraction
       operands swapped. */
    return vec_sel(b, a, vec_sub(b, a));
}

static gmx_inline gmx_simd4_real_t gmx_always_inline gmx_simd4_max_r(gmx_simd4_real_t a, gmx_simd4_real_t b)
{
    return vec_sel(b, a, vec_sub(a, b));
}

static gmx_inline gmx_simd4_real_t gmx_always_inline gmx_simd4_blendzero_r(gmx_simd4_real_t a, gmx_simd4_real_t b)
{
    return vec_sel(gmx_simd_setzero_r(), a, b);
}

static gmx_inline gmx_simd4_bool_t gmx_always_inline gmx_simd4_cmplt_r(gmx_simd4_real_t a, gmx_simd4_real_t b)
{
    return vec_cmplt(a, b);
}

static gmx_inline gmx_simd4_bool_t gmx_always_inline gmx_simd4_and_b(gmx_simd4_bool_t a, gmx_simd4_bool_t b)
{
    return vec_and(a, b);
}

static gmx_inline gmx_simd4_bool_t gmx_always_inline gmx_simd4_or_b(gmx_simd4_bool_t a, gmx_simd4_bool_t b)
{
    return vec_or(a, b);
}

static gmx_inline float gmx_always_inline gmx_simd4_dotproduct3_r(gmx_simd4_real_t a, gmx_simd4_real_t b)
{
    /* The dot product is done solely on the QPX AXU (which is the
       only available FPU). This is awkward, because pretty much no
       "horizontal" SIMD-vector operations exist, unlike x86 where
       SSE4.1 added various kinds of horizontal operations. So we have
       to make do with shifting vector elements and operating on the
       results. This makes for lots of data dependency, but the main
       alternative of storing to memory and reloading is not going to
       help, either. OpenMP over 2 or 4 hardware threads per core will
       hide much of the latency from the data dependency. The
       vec_extract() lets the compiler correctly use a floating-point
       comparison on the zeroth vector element, which avoids needing
       memory at all.
     */

    gmx_simd4_real_t dp_shifted_left_0 = vec_mul(a, b);
    gmx_simd4_real_t dp_shifted_left_1 = vec_sldw(dp_shifted_left_0, dp_shifted_left_0, 1);
    gmx_simd4_real_t dp_shifted_left_2 = vec_sldw(dp_shifted_left_0, dp_shifted_left_0, 2);
    gmx_simd4_real_t dp                = vec_add(dp_shifted_left_2,
                                                 vec_add(dp_shifted_left_0, dp_shifted_left_1));

    /* See comment in nbnxn_make_pairlist_part() about how this should
       be able to return a double on PowerPC. */
    return (float) vec_extract(dp, 0);
}

static gmx_inline int gmx_always_inline gmx_simd4_anytrue_b(gmx_simd4_bool_t a)
{
    return gmx_simd_anytrue_b(a);
}

#undef gmx_always_inline

#endif /* GMX_SIMD_IBM_QPX */

#ifdef GMX_HAVE_SIMD4_MACROS
/* Generic functions to extract a SIMD4 aligned pointer from a pointer x.
 * x should have at least GMX_SIMD4_WIDTH=4 elements extra compared
 * to how many you want to use, to avoid indexing outside the aligned region.
 */

static gmx_inline gmx_simd4_real *
gmx_simd4_align_real(const gmx_simd4_real *x)
{
    return (gmx_simd4_real *)(((size_t)((x)+GMX_SIMD4_WIDTH)) & (~((size_t)(GMX_SIMD4_WIDTH*sizeof(gmx_simd4_real)-1))));
}
#endif


#endif