src/mdlib/nbnxn_kernels/nbnxn_kernel_simd_utils_x86_128d.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_128d_h_
  38 #define _nbnxn_kernel_simd_utils_x86_128d_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 #define gmx_exclfilter gmx_epi32
  49 static const int filter_stride = GMX_SIMD_EPI32_WIDTH/GMX_SIMD_WIDTH_HERE;
  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 __m128d
  62 gmx_mm_transpose_sum2_pr(__m128d in0, __m128d in1)
  63 {
  64     __m128d tr0, tr1;
  65
  66     gmx_mm_transpose2_op_pd(in0, in1, &tr0, &tr1);
  67
  68     return _mm_add_pd(tr0, tr1);
  69 }
  70
  71 static inline __m128
  72 gmx_mm128_invsqrt_ps_single(__m128 x)
  73 {
  74     const __m128 half  = _mm_set_ps(0.5, 0.5, 0.5, 0.5);
  75     const __m128 three = _mm_set_ps(3.0, 3.0, 3.0, 3.0);
  76
  77     __m128       lu = _mm_rsqrt_ps(x);
  78
  79     return _mm_mul_ps(half, _mm_mul_ps(_mm_sub_ps(three, _mm_mul_ps(_mm_mul_ps(lu, lu), x)), lu));
  80 }
  81
  82 /* Do 2 double precision invsqrt operations.
  83  * Doing the SIMD rsqrt and the first Newton Raphson iteration
  84  * in single precision gives full double precision accuracy.
  85  */
  86 static gmx_inline void
  87 gmx_mm_invsqrt2_pd(__m128d in0, __m128d in1,
  88                    __m128d *out0, __m128d *out1)
  89 {
  90     const __m128d half  = _mm_set1_pd(0.5);
  91     const __m128d three = _mm_set1_pd(3.0);
  92     __m128        s, ir;
  93     __m128d       lu0, lu1;
  94
  95     s     = _mm_movelh_ps(_mm_cvtpd_ps(in0), _mm_cvtpd_ps(in1));
  96     ir    = gmx_mm128_invsqrt_ps_single(s);
  97     lu0   = _mm_cvtps_pd(ir);
  98     lu1   = _mm_cvtps_pd(_mm_movehl_ps(ir, ir));
  99     *out0 = _mm_mul_pd(half, _mm_mul_pd(_mm_sub_pd(three, _mm_mul_pd(_mm_mul_pd(lu0, lu0), in0)), lu0));
 100     *out1 = _mm_mul_pd(half, _mm_mul_pd(_mm_sub_pd(three, _mm_mul_pd(_mm_mul_pd(lu1, lu1), in1)), lu1));
 101 }
 102
 103 static gmx_inline void
 104 load_lj_pair_params(const real *nbfp, const int *type, int aj,
 105                     __m128d *c6_S, __m128d *c12_S)
 106 {
 107     __m128d clj_S[UNROLLJ];
 108     int     p;
 109
 110     for (p = 0; p < UNROLLJ; p++)
 111     {
 112         clj_S[p] = _mm_load_pd(nbfp+type[aj+p]*nbfp_stride);
 113     }
 114     gmx_mm_transpose2_op_pd(clj_S[0], clj_S[1], c6_S, c12_S);
 115 }
 116
 117 /* The load_table functions below are performance critical.
 118  * The routines issue UNROLLI*UNROLLJ _mm_load_ps calls.
 119  * As these all have latencies, scheduling is crucial.
 120  * The Intel compilers and CPUs seem to do a good job at this.
 121  * But AMD CPUs perform significantly worse with gcc than with icc.
 122  * Performance is improved a bit by using the extract function UNROLLJ times,
 123  * instead of doing an _mm_store_si128 for every i-particle.
 124  * This is only faster when we use FDV0 formatted tables, where we also need
 125  * to multiple the index by 4, which can be done by a SIMD bit shift.
 126  * With single precision AVX, 8 extracts are much slower than 1 store.
 127  * Because of this, the load_table_f function always takes the ti
 128  * parameter, which should contain a buffer that is aligned with
 129  * prepare_table_load_buffer(), but it is only used with full-width
 130  * AVX_256. */
 131
 132 static gmx_inline void
 133 load_table_f(const real *tab_coul_F, gmx_epi32 ti_S, int *ti,
 134              __m128d *ctab0_S, __m128d *ctab1_S)
 135 {
 136     int     idx[2];
 137     __m128d ctab_S[2];
 138
 139     /* Without SSE4.1 the extract macro needs an immediate: unroll */
 140     idx[0]    = gmx_mm_extract_epi32(ti_S, 0);
 141     ctab_S[0] = _mm_loadu_pd(tab_coul_F+idx[0]);
 142     idx[1]    = gmx_mm_extract_epi32(ti_S, 1);
 143     ctab_S[1] = _mm_loadu_pd(tab_coul_F+idx[1]);
 144
 145     /* Shuffle the force table entries to a convenient order */
 146     gmx_mm_transpose2_op_pd(ctab_S[0], ctab_S[1], ctab0_S, ctab1_S);
 147     /* The second force table entry should contain the difference */
 148     *ctab1_S = _mm_sub_pd(*ctab1_S, *ctab0_S);
 149 }
 150
 151 static gmx_inline void
 152 load_table_f_v(const real *tab_coul_F, const real *tab_coul_V,
 153                gmx_epi32 ti_S, int *ti,
 154                __m128d *ctab0_S, __m128d *ctab1_S, __m128d *ctabv_S)
 155 {
 156     int     idx[2];
 157     __m128d ctab_S[4];
 158
 159     /* Without SSE4.1 the extract macro needs an immediate: unroll */
 160     idx[0]    = gmx_mm_extract_epi32(ti_S, 0);
 161     ctab_S[0] = _mm_loadu_pd(tab_coul_F+idx[0]);
 162     idx[1]    = gmx_mm_extract_epi32(ti_S, 1);
 163     ctab_S[1] = _mm_loadu_pd(tab_coul_F+idx[1]);
 164
 165     /* Shuffle the force table entries to a convenient order */
 166     gmx_mm_transpose2_op_pd(ctab_S[0], ctab_S[1], ctab0_S, ctab1_S);
 167     /* The second force table entry should contain the difference */
 168     *ctab1_S = _mm_sub_pd(*ctab1_S, *ctab0_S);
 169
 170     ctab_S[2] = _mm_loadu_pd(tab_coul_V+idx[0]);
 171     ctab_S[3] = _mm_loadu_pd(tab_coul_V+idx[1]);
 172
 173     /* Shuffle the energy table entries to a single register */
 174     *ctabv_S = _mm_shuffle_pd(ctab_S[2], ctab_S[3], _MM_SHUFFLE2(0, 0));
 175 }
 176
 177 static gmx_inline gmx_exclfilter
 178 gmx_load1_exclfilter(int e)
 179 {
 180     return _mm_set1_epi32(e);
 181 }
 182
 183 static gmx_inline gmx_exclfilter
 184 gmx_load_exclusion_filter(const unsigned *i)
 185 {
 186     return _mm_load_si128((__m128i *) i);
 187 }
 188
 189 static gmx_inline gmx_mm_pb
 190 gmx_checkbitmask_pb(gmx_exclfilter m0, gmx_exclfilter m1)
 191 {
 192     return gmx_mm_castsi128_pd(_mm_cmpeq_epi32(_mm_andnot_si128(m0, m1), _mm_setzero_si128()));
 193 }
 194
 195 #endif /* _nbnxn_kernel_simd_utils_x86_s128d_h_ */