src/mdlib/nbnxn_kernels/nbnxn_kernel_simd_utils_x86_256d.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2012, The GROMACS Development Team
   6  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   7  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   8  * others, as listed in the AUTHORS file in the top-level source
   9  * directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #ifndef _nbnxn_kernel_simd_utils_x86_256d_h_
  38 #define _nbnxn_kernel_simd_utils_x86_256d_h_
  39
  40 /* This files contains all functions/macros for the SIMD kernels
  41  * which have explicit dependencies on the j-cluster size and/or SIMD-width.
  42  * The functionality which depends on the j-cluster size is:
  43  *   LJ-parameter lookup
  44  *   force table lookup
  45  *   energy group pair energy storage
  46  */
  47
  48 typedef gmx_mm_pr gmx_exclfilter;
  49 static const int filter_stride = 2;
  50
  51 /* Transpose 2 double precision registers */
  52 static gmx_inline void
  53 gmx_mm_transpose2_op_pd(__m128d in0, __m128d in1,
  54                         __m128d *out0, __m128d *out1)
  55 {
  56     *out0 = _mm_unpacklo_pd(in0, in1);
  57     *out1 = _mm_unpackhi_pd(in0, in1);
  58 }
  59
  60 /* Sum the elements within each input register and store the sums in out */
  61 static gmx_inline __m256d
  62 gmx_mm_transpose_sum4_pr(__m256d in0, __m256d in1,
  63                          __m256d in2, __m256d in3)
  64 {
  65     in0 = _mm256_hadd_pd(in0, in1);
  66     in2 = _mm256_hadd_pd(in2, in3);
  67
  68     return _mm256_add_pd(_mm256_permute2f128_pd(in0, in2, 0x20), _mm256_permute2f128_pd(in0, in2, 0x31));
  69 }
  70
  71 static gmx_inline __m256
  72 gmx_mm256_invsqrt_ps_single(__m256 x)
  73 {
  74     const __m256 half  = _mm256_set_ps(0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5);
  75     const __m256 three = _mm256_set_ps(3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0);
  76
  77     __m256       lu = _mm256_rsqrt_ps(x);
  78
  79     return _mm256_mul_ps(half, _mm256_mul_ps(_mm256_sub_ps(three, _mm256_mul_ps(_mm256_mul_ps(lu, lu), x)), lu));
  80 }
  81
  82 /* Put two 128-bit 4-float registers into one 256-bit 8-float register */
  83 static gmx_inline __m256
  84 gmx_2_m128_to_m256(__m128 in0, __m128 in1)
  85 {
  86     return _mm256_insertf128_ps(_mm256_castps128_ps256(in0), in1, 1);
  87 }
  88
  89 /* Put two 128-bit 2-double registers into one 256-bit 4-double register */
  90 static gmx_inline __m256d
  91 gmx_2_m128d_to_m256d(__m128d in0, __m128d in1)
  92 {
  93     return _mm256_insertf128_pd(_mm256_castpd128_pd256(in0), in1, 1);
  94 }
  95
  96 /* Do 2 double precision invsqrt operations.
  97  * Doing the SIMD rsqrt and the first Newton Raphson iteration
  98  * in single precision gives full double precision accuracy.
  99  */
 100 static gmx_inline void
 101 gmx_mm_invsqrt2_pd(__m256d in0, __m256d in1,
 102                    __m256d *out0, __m256d *out1)
 103 {
 104     const __m256d half  = _mm256_set1_pd(0.5);
 105     const __m256d three = _mm256_set1_pd(3.0);
 106     __m256        s, ir;
 107     __m256d       lu0, lu1;
 108
 109     s     =  gmx_2_m128_to_m256(_mm256_cvtpd_ps(in0), _mm256_cvtpd_ps(in1));
 110     ir    = gmx_mm256_invsqrt_ps_single(s);
 111     lu0   = _mm256_cvtps_pd(_mm256_castps256_ps128(ir));
 112     lu1   = _mm256_cvtps_pd(_mm256_extractf128_ps(ir, 1));
 113     *out0 = _mm256_mul_pd(half, _mm256_mul_pd(_mm256_sub_pd(three, _mm256_mul_pd(_mm256_mul_pd(lu0, lu0), in0)), lu0));
 114     *out1 = _mm256_mul_pd(half, _mm256_mul_pd(_mm256_sub_pd(three, _mm256_mul_pd(_mm256_mul_pd(lu1, lu1), in1)), lu1));
 115 }
 116
 117 static gmx_inline void
 118 load_lj_pair_params(const real *nbfp, const int *type, int aj,
 119                     __m256d *c6_S, __m256d *c12_S)
 120 {
 121     __m128d clj_S[UNROLLJ], c6t_S[2], c12t_S[2];
 122     int     p;
 123
 124     for (p = 0; p < UNROLLJ; p++)
 125     {
 126         clj_S[p] = _mm_load_pd(nbfp+type[aj+p]*nbfp_stride);
 127     }
 128     gmx_mm_transpose2_op_pd(clj_S[0], clj_S[1], &c6t_S[0], &c12t_S[0]);
 129     gmx_mm_transpose2_op_pd(clj_S[2], clj_S[3], &c6t_S[1], &c12t_S[1]);
 130     *c6_S  = gmx_2_m128d_to_m256d(c6t_S[0], c6t_S[1]);
 131     *c12_S = gmx_2_m128d_to_m256d(c12t_S[0], c12t_S[1]);
 132 }
 133
 134 /* The load_table functions below are performance critical. They
 135  * always take the ti parameter, which should contain a buffer that
 136  * is aligned with prepare_table_load_buffer(), but it is only used
 137  * with full-width AVX_256. */
 138
 139 static gmx_inline void
 140 load_table_f(const real *tab_coul_F, __m128i ti_S, int *ti,
 141              __m256d *ctab0_S, __m256d *ctab1_S)
 142 {
 143     __m128d ctab_S[4], tr_S[4];
 144     int     j;
 145
 146     _mm_store_si128((__m128i *)ti, ti_S);
 147     for (j = 0; j < 4; j++)
 148     {
 149         ctab_S[j] = _mm_loadu_pd(tab_coul_F+ti[j]);
 150     }
 151     /* Shuffle the force table entries to a convenient order */
 152     gmx_mm_transpose2_op_pd(ctab_S[0], ctab_S[1], &tr_S[0], &tr_S[1]);
 153     gmx_mm_transpose2_op_pd(ctab_S[2], ctab_S[3], &tr_S[2], &tr_S[3]);
 154     *ctab0_S = gmx_2_m128d_to_m256d(tr_S[0], tr_S[2]);
 155     *ctab1_S = gmx_2_m128d_to_m256d(tr_S[1], tr_S[3]);
 156     /* The second force table entry should contain the difference */
 157     *ctab1_S = _mm256_sub_pd(*ctab1_S, *ctab0_S);
 158 }
 159
 160 static gmx_inline void
 161 load_table_f_v(const real *tab_coul_F, const real *tab_coul_V,
 162                __m128i ti_S, int *ti,
 163                __m256d *ctab0_S, __m256d *ctab1_S, __m256d *ctabv_S)
 164 {
 165     __m128d ctab_S[8], tr_S[4];
 166     int     j;
 167
 168     _mm_store_si128((__m128i *)ti, ti_S);
 169     for (j = 0; j < 4; j++)
 170     {
 171         ctab_S[j] = _mm_loadu_pd(tab_coul_F+ti[j]);
 172     }
 173     /* Shuffle the force table entries to a convenient order */
 174     gmx_mm_transpose2_op_pd(ctab_S[0], ctab_S[1], &tr_S[0], &tr_S[1]);
 175     gmx_mm_transpose2_op_pd(ctab_S[2], ctab_S[3], &tr_S[2], &tr_S[3]);
 176     *ctab0_S = gmx_2_m128d_to_m256d(tr_S[0], tr_S[2]);
 177     *ctab1_S = gmx_2_m128d_to_m256d(tr_S[1], tr_S[3]);
 178     /* The second force table entry should contain the difference */
 179     *ctab1_S = _mm256_sub_pd(*ctab1_S, *ctab0_S);
 180
 181     for (j = 0; j < 4; j++)
 182     {
 183         ctab_S[4+j] = _mm_loadu_pd(tab_coul_V+ti[j]);
 184     }
 185     /* Shuffle the energy table entries to a single register */
 186     *ctabv_S = gmx_2_m128d_to_m256d(_mm_shuffle_pd(ctab_S[4], ctab_S[5], _MM_SHUFFLE2(0, 0)), _mm_shuffle_pd(ctab_S[6], ctab_S[7], _MM_SHUFFLE2(0, 0)));
 187 }
 188
 189 static gmx_inline gmx_exclfilter
 190 gmx_load1_exclfilter(int e)
 191 {
 192     return _mm256_castsi256_pd(_mm256_set1_epi32(e));
 193 }
 194
 195 static gmx_inline gmx_exclfilter
 196 gmx_load_exclusion_filter(const unsigned *i)
 197 {
 198     return gmx_load_pr((real *) (i));
 199 }
 200
 201 static gmx_inline gmx_mm_pb
 202 gmx_checkbitmask_pb(gmx_exclfilter m0, gmx_exclfilter m1)
 203 {
 204     /* With <= 16 bits used the cast and conversion should not be
 205      * required, since only mantissa bits are set and that would give
 206      * a non-zero float, but with the Intel compiler this does not
 207      * work correctly. Because AVX does not have int->double
 208      * conversion, we convert via float. */
 209     return _mm256_cmp_pd(_mm256_castps_pd(_mm256_cvtepi32_ps(_mm256_castpd_si256(_mm256_and_pd(m0, m1)))), _mm256_setzero_pd(), 0x0c);
 210 }
 211
 212 #endif /* _nbnxn_kernel_simd_utils_x86_s256d_h_ */