src/gromacs/gmxlib/nonbonded/nb_kernel_sse4_1_double/kernelutil_x86_sse4_1_double.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013,2014, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 #ifndef _kernelutil_x86_sse4_1_double_h_
  36 #define _kernelutil_x86_sse4_1_double_h_
  37
  38 #include <math.h>
  39 #include <stdio.h>
  40
  41 #include "config.h"
  42
  43 #define gmx_mm_extract_epi32(x, imm) _mm_extract_epi32((x), (imm))
  44 #define gmx_mm_castsi128_pd(a) _mm_castsi128_pd(a)
  45
  46 #define GMX_MM_TRANSPOSE2_PD(row0, row1) {           \
  47         __m128d __gmx_t1 = row0;                         \
  48         row0           = _mm_unpacklo_pd(row0, row1);     \
  49         row1           = _mm_unpackhi_pd(__gmx_t1, row1); \
  50 }
  51
  52 /* Normal sum of four ymm registers */
  53 #define gmx_mm_sum4_pd(t0, t1, t2, t3)  _mm_add_pd(_mm_add_pd(t0, t1), _mm_add_pd(t2, t3))
  54
  55 static gmx_inline int gmx_simdcall
  56 gmx_mm_any_lt(__m128d a, __m128d b)
  57 {
  58     return _mm_movemask_pd(_mm_cmplt_pd(a, b));
  59 }
  60
  61 static gmx_inline __m128d gmx_simdcall
  62 gmx_mm_calc_rsq_pd(__m128d dx, __m128d dy, __m128d dz)
  63 {
  64     return _mm_add_pd( _mm_add_pd( _mm_mul_pd(dx, dx), _mm_mul_pd(dy, dy) ), _mm_mul_pd(dz, dz) );
  65 }
  66
  67
  68 /* Load a double value from 1-2 places, merge into xmm register */
  69 static gmx_inline __m128d gmx_simdcall
  70 gmx_mm_load_2real_swizzle_pd(const double * gmx_restrict ptrA,
  71                              const double * gmx_restrict ptrB)
  72 {
  73     return _mm_unpacklo_pd(_mm_load_sd(ptrA), _mm_load_sd(ptrB));
  74 }
  75
  76 static gmx_inline __m128d gmx_simdcall
  77 gmx_mm_load_1real_pd(const double * gmx_restrict ptrA)
  78 {
  79     return _mm_load_sd(ptrA);
  80 }
  81
  82
  83 static gmx_inline void gmx_simdcall
  84 gmx_mm_store_2real_swizzle_pd(double * gmx_restrict ptrA,
  85                               double * gmx_restrict ptrB,
  86                               __m128d               xmm1)
  87 {
  88     __m128d t2;
  89
  90     t2       = _mm_unpackhi_pd(xmm1, xmm1);
  91     _mm_store_sd(ptrA, xmm1);
  92     _mm_store_sd(ptrB, t2);
  93 }
  94
  95 static gmx_inline void gmx_simdcall
  96 gmx_mm_store_1real_pd(double * gmx_restrict ptrA, __m128d xmm1)
  97 {
  98     _mm_store_sd(ptrA, xmm1);
  99 }
 100
 101
 102 /* Similar to store, but increments value in memory */
 103 static gmx_inline void gmx_simdcall
 104 gmx_mm_increment_2real_swizzle_pd(double * gmx_restrict ptrA,
 105                                   double * gmx_restrict ptrB, __m128d xmm1)
 106 {
 107     __m128d t1;
 108
 109     t1   = _mm_unpackhi_pd(xmm1, xmm1);
 110     xmm1 = _mm_add_sd(xmm1, _mm_load_sd(ptrA));
 111     t1   = _mm_add_sd(t1, _mm_load_sd(ptrB));
 112     _mm_store_sd(ptrA, xmm1);
 113     _mm_store_sd(ptrB, t1);
 114 }
 115
 116 static gmx_inline void gmx_simdcall
 117 gmx_mm_increment_1real_pd(double * gmx_restrict ptrA, __m128d xmm1)
 118 {
 119     __m128d tmp;
 120
 121     tmp = gmx_mm_load_1real_pd(ptrA);
 122     tmp = _mm_add_sd(tmp, xmm1);
 123     gmx_mm_store_1real_pd(ptrA, tmp);
 124 }
 125
 126
 127 static gmx_inline void gmx_simdcall
 128 gmx_mm_load_2pair_swizzle_pd(const double * gmx_restrict p1,
 129                              const double * gmx_restrict p2,
 130                              __m128d * gmx_restrict      c6,
 131                              __m128d * gmx_restrict      c12)
 132 {
 133     __m128d t1, t2, t3;
 134
 135     t1   = _mm_loadu_pd(p1);
 136     t2   = _mm_loadu_pd(p2);
 137     *c6  = _mm_unpacklo_pd(t1, t2);
 138     *c12 = _mm_unpackhi_pd(t1, t2);
 139 }
 140
 141 static gmx_inline void gmx_simdcall
 142 gmx_mm_load_1pair_swizzle_pd(const double * gmx_restrict p1,
 143                              __m128d * gmx_restrict      c6,
 144                              __m128d * gmx_restrict      c12)
 145 {
 146     *c6     = _mm_load_sd(p1);
 147     *c12    = _mm_load_sd(p1+1);
 148 }
 149
 150
 151
 152 static gmx_inline void gmx_simdcall
 153 gmx_mm_load_shift_and_1rvec_broadcast_pd(const double * gmx_restrict xyz_shift,
 154                                          const double * gmx_restrict xyz,
 155                                          __m128d * gmx_restrict      x1,
 156                                          __m128d * gmx_restrict      y1,
 157                                          __m128d * gmx_restrict      z1)
 158 {
 159     __m128d mem_xy, mem_z, mem_sxy, mem_sz;
 160
 161     mem_xy  = _mm_loadu_pd(xyz);
 162     mem_z   = _mm_load_sd(xyz+2);
 163     mem_sxy = _mm_loadu_pd(xyz_shift);
 164     mem_sz  = _mm_load_sd(xyz_shift+2);
 165
 166     mem_xy  = _mm_add_pd(mem_xy, mem_sxy);
 167     mem_z   = _mm_add_pd(mem_z, mem_sz);
 168
 169     *x1  = _mm_shuffle_pd(mem_xy, mem_xy, _MM_SHUFFLE2(0, 0));
 170     *y1  = _mm_shuffle_pd(mem_xy, mem_xy, _MM_SHUFFLE2(1, 1));
 171     *z1  = _mm_shuffle_pd(mem_z, mem_z, _MM_SHUFFLE2(0, 0));
 172 }
 173
 174
 175 static gmx_inline void gmx_simdcall
 176 gmx_mm_load_shift_and_3rvec_broadcast_pd(const double * gmx_restrict xyz_shift,
 177                                          const double * gmx_restrict xyz,
 178                                          __m128d * gmx_restrict x1, __m128d * gmx_restrict y1, __m128d * gmx_restrict z1,
 179                                          __m128d * gmx_restrict x2, __m128d * gmx_restrict y2, __m128d * gmx_restrict z2,
 180                                          __m128d * gmx_restrict x3, __m128d * gmx_restrict y3, __m128d * gmx_restrict z3)
 181 {
 182     __m128d t1, t2, t3, t4, t5, sxy, sz, szx, syz;
 183
 184     t1  = _mm_loadu_pd(xyz);
 185     t2  = _mm_loadu_pd(xyz+2);
 186     t3  = _mm_loadu_pd(xyz+4);
 187     t4  = _mm_loadu_pd(xyz+6);
 188     t5  = _mm_load_sd(xyz+8);
 189
 190     sxy = _mm_loadu_pd(xyz_shift);
 191     sz  = _mm_load_sd(xyz_shift+2);
 192     szx = _mm_shuffle_pd(sz, sxy, _MM_SHUFFLE2(0, 0));
 193     syz = _mm_shuffle_pd(sxy, sz, _MM_SHUFFLE2(0, 1));
 194
 195     t1  = _mm_add_pd(t1, sxy);
 196     t2  = _mm_add_pd(t2, szx);
 197     t3  = _mm_add_pd(t3, syz);
 198     t4  = _mm_add_pd(t4, sxy);
 199     t5  = _mm_add_sd(t5, sz);
 200
 201     *x1  = _mm_shuffle_pd(t1, t1, _MM_SHUFFLE2(0, 0));
 202     *y1  = _mm_shuffle_pd(t1, t1, _MM_SHUFFLE2(1, 1));
 203     *z1  = _mm_shuffle_pd(t2, t2, _MM_SHUFFLE2(0, 0));
 204     *x2  = _mm_shuffle_pd(t2, t2, _MM_SHUFFLE2(1, 1));
 205     *y2  = _mm_shuffle_pd(t3, t3, _MM_SHUFFLE2(0, 0));
 206     *z2  = _mm_shuffle_pd(t3, t3, _MM_SHUFFLE2(1, 1));
 207     *x3  = _mm_shuffle_pd(t4, t4, _MM_SHUFFLE2(0, 0));
 208     *y3  = _mm_shuffle_pd(t4, t4, _MM_SHUFFLE2(1, 1));
 209     *z3  = _mm_shuffle_pd(t5, t5, _MM_SHUFFLE2(0, 0));
 210 }
 211
 212
 213 static gmx_inline void gmx_simdcall
 214 gmx_mm_load_shift_and_4rvec_broadcast_pd(const double * gmx_restrict xyz_shift,
 215                                          const double * gmx_restrict xyz,
 216                                          __m128d * gmx_restrict x1, __m128d * gmx_restrict y1, __m128d * gmx_restrict z1,
 217                                          __m128d * gmx_restrict x2, __m128d * gmx_restrict y2, __m128d * gmx_restrict z2,
 218                                          __m128d * gmx_restrict x3, __m128d * gmx_restrict y3, __m128d * gmx_restrict z3,
 219                                          __m128d * gmx_restrict x4, __m128d * gmx_restrict y4, __m128d * gmx_restrict z4)
 220 {
 221     __m128d t1, t2, t3, t4, t5, t6, sxy, sz, szx, syz;
 222
 223     t1  = _mm_loadu_pd(xyz);
 224     t2  = _mm_loadu_pd(xyz+2);
 225     t3  = _mm_loadu_pd(xyz+4);
 226     t4  = _mm_loadu_pd(xyz+6);
 227     t5  = _mm_loadu_pd(xyz+8);
 228     t6  = _mm_loadu_pd(xyz+10);
 229
 230     sxy = _mm_loadu_pd(xyz_shift);
 231     sz  = _mm_load_sd(xyz_shift+2);
 232     szx = _mm_shuffle_pd(sz, sxy, _MM_SHUFFLE2(0, 0));
 233     syz = _mm_shuffle_pd(sxy, sz, _MM_SHUFFLE2(0, 1));
 234
 235     t1  = _mm_add_pd(t1, sxy);
 236     t2  = _mm_add_pd(t2, szx);
 237     t3  = _mm_add_pd(t3, syz);
 238     t4  = _mm_add_pd(t4, sxy);
 239     t5  = _mm_add_pd(t5, szx);
 240     t6  = _mm_add_pd(t6, syz);
 241
 242     *x1  = _mm_shuffle_pd(t1, t1, _MM_SHUFFLE2(0, 0));
 243     *y1  = _mm_shuffle_pd(t1, t1, _MM_SHUFFLE2(1, 1));
 244     *z1  = _mm_shuffle_pd(t2, t2, _MM_SHUFFLE2(0, 0));
 245     *x2  = _mm_shuffle_pd(t2, t2, _MM_SHUFFLE2(1, 1));
 246     *y2  = _mm_shuffle_pd(t3, t3, _MM_SHUFFLE2(0, 0));
 247     *z2  = _mm_shuffle_pd(t3, t3, _MM_SHUFFLE2(1, 1));
 248     *x3  = _mm_shuffle_pd(t4, t4, _MM_SHUFFLE2(0, 0));
 249     *y3  = _mm_shuffle_pd(t4, t4, _MM_SHUFFLE2(1, 1));
 250     *z3  = _mm_shuffle_pd(t5, t5, _MM_SHUFFLE2(0, 0));
 251     *x4  = _mm_shuffle_pd(t5, t5, _MM_SHUFFLE2(1, 1));
 252     *y4  = _mm_shuffle_pd(t6, t6, _MM_SHUFFLE2(0, 0));
 253     *z4  = _mm_shuffle_pd(t6, t6, _MM_SHUFFLE2(1, 1));
 254 }
 255
 256
 257
 258
 259 static gmx_inline void gmx_simdcall
 260 gmx_mm_load_1rvec_1ptr_swizzle_pd(const double * gmx_restrict p1,
 261                                   __m128d * gmx_restrict x, __m128d * gmx_restrict y, __m128d * gmx_restrict z)
 262 {
 263     *x            = _mm_load_sd(p1);
 264     *y            = _mm_load_sd(p1+1);
 265     *z            = _mm_load_sd(p1+2);
 266 }
 267
 268 static gmx_inline void gmx_simdcall
 269 gmx_mm_load_3rvec_1ptr_swizzle_pd(const double * gmx_restrict p1,
 270                                   __m128d * gmx_restrict x1, __m128d * gmx_restrict y1, __m128d * gmx_restrict z1,
 271                                   __m128d * gmx_restrict x2, __m128d * gmx_restrict y2, __m128d * gmx_restrict z2,
 272                                   __m128d * gmx_restrict x3, __m128d * gmx_restrict y3, __m128d * gmx_restrict z3)
 273 {
 274     *x1            = _mm_load_sd(p1);
 275     *y1            = _mm_load_sd(p1+1);
 276     *z1            = _mm_load_sd(p1+2);
 277     *x2            = _mm_load_sd(p1+3);
 278     *y2            = _mm_load_sd(p1+4);
 279     *z2            = _mm_load_sd(p1+5);
 280     *x3            = _mm_load_sd(p1+6);
 281     *y3            = _mm_load_sd(p1+7);
 282     *z3            = _mm_load_sd(p1+8);
 283 }
 284
 285 static gmx_inline void gmx_simdcall
 286 gmx_mm_load_4rvec_1ptr_swizzle_pd(const double * gmx_restrict p1,
 287                                   __m128d * gmx_restrict x1, __m128d * gmx_restrict y1, __m128d * gmx_restrict z1,
 288                                   __m128d * gmx_restrict x2, __m128d * gmx_restrict y2, __m128d * gmx_restrict z2,
 289                                   __m128d * gmx_restrict x3, __m128d * gmx_restrict y3, __m128d * gmx_restrict z3,
 290                                   __m128d * gmx_restrict x4, __m128d * gmx_restrict y4, __m128d * gmx_restrict z4)
 291 {
 292     *x1            = _mm_load_sd(p1);
 293     *y1            = _mm_load_sd(p1+1);
 294     *z1            = _mm_load_sd(p1+2);
 295     *x2            = _mm_load_sd(p1+3);
 296     *y2            = _mm_load_sd(p1+4);
 297     *z2            = _mm_load_sd(p1+5);
 298     *x3            = _mm_load_sd(p1+6);
 299     *y3            = _mm_load_sd(p1+7);
 300     *z3            = _mm_load_sd(p1+8);
 301     *x4            = _mm_load_sd(p1+9);
 302     *y4            = _mm_load_sd(p1+10);
 303     *z4            = _mm_load_sd(p1+11);
 304 }
 305
 306
 307 static gmx_inline void gmx_simdcall
 308 gmx_mm_load_1rvec_2ptr_swizzle_pd(const double * gmx_restrict ptrA,
 309                                   const double * gmx_restrict ptrB,
 310                                   __m128d * gmx_restrict x1, __m128d * gmx_restrict y1, __m128d * gmx_restrict z1)
 311 {
 312     __m128d t1, t2, t3, t4;
 313     t1           = _mm_loadu_pd(ptrA);
 314     t2           = _mm_loadu_pd(ptrB);
 315     t3           = _mm_load_sd(ptrA+2);
 316     t4           = _mm_load_sd(ptrB+2);
 317     GMX_MM_TRANSPOSE2_PD(t1, t2);
 318     *x1          = t1;
 319     *y1          = t2;
 320     *z1          = _mm_unpacklo_pd(t3, t4);
 321 }
 322
 323
 324 static gmx_inline void gmx_simdcall
 325 gmx_mm_load_3rvec_2ptr_swizzle_pd(const double * gmx_restrict ptrA, const double * gmx_restrict ptrB,
 326                                   __m128d * gmx_restrict x1, __m128d * gmx_restrict y1, __m128d * gmx_restrict z1,
 327                                   __m128d * gmx_restrict x2, __m128d * gmx_restrict y2, __m128d * gmx_restrict z2,
 328                                   __m128d * gmx_restrict x3, __m128d * gmx_restrict y3, __m128d * gmx_restrict z3)
 329 {
 330     __m128d t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
 331     t1           = _mm_loadu_pd(ptrA);
 332     t2           = _mm_loadu_pd(ptrB);
 333     t3           = _mm_loadu_pd(ptrA+2);
 334     t4           = _mm_loadu_pd(ptrB+2);
 335     t5           = _mm_loadu_pd(ptrA+4);
 336     t6           = _mm_loadu_pd(ptrB+4);
 337     t7           = _mm_loadu_pd(ptrA+6);
 338     t8           = _mm_loadu_pd(ptrB+6);
 339     t9           = _mm_load_sd(ptrA+8);
 340     t10          = _mm_load_sd(ptrB+8);
 341     GMX_MM_TRANSPOSE2_PD(t1, t2);
 342     GMX_MM_TRANSPOSE2_PD(t3, t4);
 343     GMX_MM_TRANSPOSE2_PD(t5, t6);
 344     GMX_MM_TRANSPOSE2_PD(t7, t8);
 345     *x1          = t1;
 346     *y1          = t2;
 347     *z1          = t3;
 348     *x2          = t4;
 349     *y2          = t5;
 350     *z2          = t6;
 351     *x3          = t7;
 352     *y3          = t8;
 353     *z3          = _mm_unpacklo_pd(t9, t10);
 354 }
 355
 356
 357 static gmx_inline void gmx_simdcall
 358 gmx_mm_load_4rvec_2ptr_swizzle_pd(const double * gmx_restrict ptrA, const double * gmx_restrict ptrB,
 359                                   __m128d * gmx_restrict x1, __m128d * gmx_restrict y1, __m128d * gmx_restrict z1,
 360                                   __m128d * gmx_restrict x2, __m128d * gmx_restrict y2, __m128d * gmx_restrict z2,
 361                                   __m128d * gmx_restrict x3, __m128d * gmx_restrict y3, __m128d * gmx_restrict z3,
 362                                   __m128d * gmx_restrict x4, __m128d * gmx_restrict y4, __m128d * gmx_restrict z4)
 363 {
 364     __m128d t1, t2, t3, t4, t5, t6;
 365     t1           = _mm_loadu_pd(ptrA);
 366     t2           = _mm_loadu_pd(ptrB);
 367     t3           = _mm_loadu_pd(ptrA+2);
 368     t4           = _mm_loadu_pd(ptrB+2);
 369     t5           = _mm_loadu_pd(ptrA+4);
 370     t6           = _mm_loadu_pd(ptrB+4);
 371     GMX_MM_TRANSPOSE2_PD(t1, t2);
 372     GMX_MM_TRANSPOSE2_PD(t3, t4);
 373     GMX_MM_TRANSPOSE2_PD(t5, t6);
 374     *x1          = t1;
 375     *y1          = t2;
 376     *z1          = t3;
 377     *x2          = t4;
 378     *y2          = t5;
 379     *z2          = t6;
 380     t1           = _mm_loadu_pd(ptrA+6);
 381     t2           = _mm_loadu_pd(ptrB+6);
 382     t3           = _mm_loadu_pd(ptrA+8);
 383     t4           = _mm_loadu_pd(ptrB+8);
 384     t5           = _mm_loadu_pd(ptrA+10);
 385     t6           = _mm_loadu_pd(ptrB+10);
 386     GMX_MM_TRANSPOSE2_PD(t1, t2);
 387     GMX_MM_TRANSPOSE2_PD(t3, t4);
 388     GMX_MM_TRANSPOSE2_PD(t5, t6);
 389     *x3          = t1;
 390     *y3          = t2;
 391     *z3          = t3;
 392     *x4          = t4;
 393     *y4          = t5;
 394     *z4          = t6;
 395 }
 396
 397
 398 /* Routines to decrement rvec in memory, typically use for j particle force updates */
 399 static gmx_inline void gmx_simdcall
 400 gmx_mm_decrement_1rvec_1ptr_noswizzle_pd(double * gmx_restrict ptrA,
 401                                          __m128d xy, __m128d z)
 402 {
 403     __m128d t1, t2;
 404
 405     t1 = _mm_loadu_pd(ptrA);
 406     t2 = _mm_load_sd(ptrA+2);
 407
 408     t1 = _mm_sub_pd(t1, xy);
 409     t2 = _mm_sub_sd(t2, z);
 410
 411     _mm_storeu_pd(ptrA, t1);
 412     _mm_store_sd(ptrA+2, t2);
 413 }
 414
 415
 416 static gmx_inline void gmx_simdcall
 417 gmx_mm_decrement_1rvec_1ptr_swizzle_pd(double * gmx_restrict ptrA,
 418                                        __m128d x1, __m128d y1, __m128d z1)
 419 {
 420     __m128d t1, t2, t3;
 421
 422     t1           = _mm_load_sd(ptrA);
 423     t2           = _mm_load_sd(ptrA+1);
 424     t3           = _mm_load_sd(ptrA+2);
 425
 426     t1           = _mm_sub_sd(t1, x1);
 427     t2           = _mm_sub_sd(t2, y1);
 428     t3           = _mm_sub_sd(t3, z1);
 429     _mm_store_sd(ptrA, t1);
 430     _mm_store_sd(ptrA+1, t2);
 431     _mm_store_sd(ptrA+2, t3);
 432 }
 433
 434
 435 static gmx_inline void gmx_simdcall
 436 gmx_mm_decrement_3rvec_1ptr_swizzle_pd(double * gmx_restrict ptrA,
 437                                        __m128d x1, __m128d y1, __m128d z1,
 438                                        __m128d x2, __m128d y2, __m128d z2,
 439                                        __m128d x3, __m128d y3, __m128d z3)
 440 {
 441     __m128d t1, t2, t3, t4, t5;
 442
 443     t1          = _mm_loadu_pd(ptrA);
 444     t2          = _mm_loadu_pd(ptrA+2);
 445     t3          = _mm_loadu_pd(ptrA+4);
 446     t4          = _mm_loadu_pd(ptrA+6);
 447     t5          = _mm_load_sd(ptrA+8);
 448
 449     x1          = _mm_unpacklo_pd(x1, y1);
 450     z1          = _mm_unpacklo_pd(z1, x2);
 451     y2          = _mm_unpacklo_pd(y2, z2);
 452     x3          = _mm_unpacklo_pd(x3, y3);
 453     /* nothing to be done for z3 */
 454
 455     t1          = _mm_sub_pd(t1, x1);
 456     t2          = _mm_sub_pd(t2, z1);
 457     t3          = _mm_sub_pd(t3, y2);
 458     t4          = _mm_sub_pd(t4, x3);
 459     t5          = _mm_sub_sd(t5, z3);
 460     _mm_storeu_pd(ptrA, t1);
 461     _mm_storeu_pd(ptrA+2, t2);
 462     _mm_storeu_pd(ptrA+4, t3);
 463     _mm_storeu_pd(ptrA+6, t4);
 464     _mm_store_sd(ptrA+8, t5);
 465 }
 466
 467
 468 static gmx_inline void gmx_simdcall
 469 gmx_mm_decrement_4rvec_1ptr_swizzle_pd(double * gmx_restrict ptrA,
 470                                        __m128d x1, __m128d y1, __m128d z1,
 471                                        __m128d x2, __m128d y2, __m128d z2,
 472                                        __m128d x3, __m128d y3, __m128d z3,
 473                                        __m128d x4, __m128d y4, __m128d z4)
 474 {
 475     __m128d t1, t2, t3, t4, t5, t6;
 476
 477     t1          = _mm_loadu_pd(ptrA);
 478     t2          = _mm_loadu_pd(ptrA+2);
 479     t3          = _mm_loadu_pd(ptrA+4);
 480     t4          = _mm_loadu_pd(ptrA+6);
 481     t5          = _mm_loadu_pd(ptrA+8);
 482     t6          = _mm_loadu_pd(ptrA+10);
 483
 484     x1          = _mm_unpacklo_pd(x1, y1);
 485     z1          = _mm_unpacklo_pd(z1, x2);
 486     y2          = _mm_unpacklo_pd(y2, z2);
 487     x3          = _mm_unpacklo_pd(x3, y3);
 488     z3          = _mm_unpacklo_pd(z3, x4);
 489     y4          = _mm_unpacklo_pd(y4, z4);
 490
 491     _mm_storeu_pd(ptrA,    _mm_sub_pd( t1, x1 ));
 492     _mm_storeu_pd(ptrA+2,  _mm_sub_pd( t2, z1 ));
 493     _mm_storeu_pd(ptrA+4,  _mm_sub_pd( t3, y2 ));
 494     _mm_storeu_pd(ptrA+6,  _mm_sub_pd( t4, x3 ));
 495     _mm_storeu_pd(ptrA+8,  _mm_sub_pd( t5, z3 ));
 496     _mm_storeu_pd(ptrA+10, _mm_sub_pd( t6, y4 ));
 497 }
 498
 499
 500 static gmx_inline void gmx_simdcall
 501 gmx_mm_decrement_1rvec_2ptr_swizzle_pd(double * gmx_restrict ptrA, double * gmx_restrict ptrB,
 502                                        __m128d x1, __m128d y1, __m128d z1)
 503 {
 504     __m128d t1, t2, t3, t4, t5, t6, t7;
 505
 506     t1          = _mm_loadu_pd(ptrA);
 507     t2          = _mm_load_sd(ptrA+2);
 508     t3          = _mm_loadu_pd(ptrB);
 509     t4          = _mm_load_sd(ptrB+2);
 510
 511     t5          = _mm_unpacklo_pd(x1, y1);
 512     t6          = _mm_unpackhi_pd(x1, y1);
 513     t7          = _mm_unpackhi_pd(z1, z1);
 514
 515     t1          = _mm_sub_pd(t1, t5);
 516     t2          = _mm_sub_sd(t2, z1);
 517
 518     t3          = _mm_sub_pd(t3, t6);
 519     t4          = _mm_sub_sd(t4, t7);
 520
 521     _mm_storeu_pd(ptrA, t1);
 522     _mm_store_sd(ptrA+2, t2);
 523     _mm_storeu_pd(ptrB, t3);
 524     _mm_store_sd(ptrB+2, t4);
 525 }
 526
 527
 528 static gmx_inline void gmx_simdcall
 529 gmx_mm_decrement_3rvec_2ptr_swizzle_pd(double * gmx_restrict ptrA, double * gmx_restrict ptrB,
 530                                        __m128d x1, __m128d y1, __m128d z1,
 531                                        __m128d x2, __m128d y2, __m128d z2,
 532                                        __m128d x3, __m128d y3, __m128d z3)
 533 {
 534     __m128d t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
 535     __m128d tA, tB, tC, tD, tE, tF, tG, tH, tI;
 536
 537     t1          = _mm_loadu_pd(ptrA);
 538     t2          = _mm_loadu_pd(ptrA+2);
 539     t3          = _mm_loadu_pd(ptrA+4);
 540     t4          = _mm_loadu_pd(ptrA+6);
 541     t5          = _mm_load_sd(ptrA+8);
 542     t6          = _mm_loadu_pd(ptrB);
 543     t7          = _mm_loadu_pd(ptrB+2);
 544     t8          = _mm_loadu_pd(ptrB+4);
 545     t9          = _mm_loadu_pd(ptrB+6);
 546     t10         = _mm_load_sd(ptrB+8);
 547
 548     tA          = _mm_unpacklo_pd(x1, y1);
 549     tB          = _mm_unpackhi_pd(x1, y1);
 550     tC          = _mm_unpacklo_pd(z1, x2);
 551     tD          = _mm_unpackhi_pd(z1, x2);
 552     tE          = _mm_unpacklo_pd(y2, z2);
 553     tF          = _mm_unpackhi_pd(y2, z2);
 554     tG          = _mm_unpacklo_pd(x3, y3);
 555     tH          = _mm_unpackhi_pd(x3, y3);
 556     tI          = _mm_unpackhi_pd(z3, z3);
 557
 558     t1          = _mm_sub_pd(t1, tA);
 559     t2          = _mm_sub_pd(t2, tC);
 560     t3          = _mm_sub_pd(t3, tE);
 561     t4          = _mm_sub_pd(t4, tG);
 562     t5          = _mm_sub_sd(t5, z3);
 563
 564     t6          = _mm_sub_pd(t6, tB);
 565     t7          = _mm_sub_pd(t7, tD);
 566     t8          = _mm_sub_pd(t8, tF);
 567     t9          = _mm_sub_pd(t9, tH);
 568     t10         = _mm_sub_sd(t10, tI);
 569
 570     _mm_storeu_pd(ptrA, t1);
 571     _mm_storeu_pd(ptrA+2, t2);
 572     _mm_storeu_pd(ptrA+4, t3);
 573     _mm_storeu_pd(ptrA+6, t4);
 574     _mm_store_sd(ptrA+8, t5);
 575     _mm_storeu_pd(ptrB, t6);
 576     _mm_storeu_pd(ptrB+2, t7);
 577     _mm_storeu_pd(ptrB+4, t8);
 578     _mm_storeu_pd(ptrB+6, t9);
 579     _mm_store_sd(ptrB+8, t10);
 580 }
 581
 582
 583 static gmx_inline void gmx_simdcall
 584 gmx_mm_decrement_4rvec_2ptr_swizzle_pd(double * gmx_restrict ptrA, double * gmx_restrict ptrB,
 585                                        __m128d x1, __m128d y1, __m128d z1,
 586                                        __m128d x2, __m128d y2, __m128d z2,
 587                                        __m128d x3, __m128d y3, __m128d z3,
 588                                        __m128d x4, __m128d y4, __m128d z4)
 589 {
 590     __m128d t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11, t12;
 591     __m128d tA, tB, tC, tD, tE, tF, tG, tH, tI, tJ, tK, tL;
 592
 593     t1          = _mm_loadu_pd(ptrA);
 594     t2          = _mm_loadu_pd(ptrA+2);
 595     t3          = _mm_loadu_pd(ptrA+4);
 596     t4          = _mm_loadu_pd(ptrA+6);
 597     t5          = _mm_loadu_pd(ptrA+8);
 598     t6          = _mm_loadu_pd(ptrA+10);
 599     t7          = _mm_loadu_pd(ptrB);
 600     t8          = _mm_loadu_pd(ptrB+2);
 601     t9          = _mm_loadu_pd(ptrB+4);
 602     t10         = _mm_loadu_pd(ptrB+6);
 603     t11         = _mm_loadu_pd(ptrB+8);
 604     t12         = _mm_loadu_pd(ptrB+10);
 605
 606     tA          = _mm_unpacklo_pd(x1, y1);
 607     tB          = _mm_unpackhi_pd(x1, y1);
 608     tC          = _mm_unpacklo_pd(z1, x2);
 609     tD          = _mm_unpackhi_pd(z1, x2);
 610     tE          = _mm_unpacklo_pd(y2, z2);
 611     tF          = _mm_unpackhi_pd(y2, z2);
 612     tG          = _mm_unpacklo_pd(x3, y3);
 613     tH          = _mm_unpackhi_pd(x3, y3);
 614     tI          = _mm_unpacklo_pd(z3, x4);
 615     tJ          = _mm_unpackhi_pd(z3, x4);
 616     tK          = _mm_unpacklo_pd(y4, z4);
 617     tL          = _mm_unpackhi_pd(y4, z4);
 618
 619     t1          = _mm_sub_pd(t1, tA);
 620     t2          = _mm_sub_pd(t2, tC);
 621     t3          = _mm_sub_pd(t3, tE);
 622     t4          = _mm_sub_pd(t4, tG);
 623     t5          = _mm_sub_pd(t5, tI);
 624     t6          = _mm_sub_pd(t6, tK);
 625
 626     t7          = _mm_sub_pd(t7, tB);
 627     t8          = _mm_sub_pd(t8, tD);
 628     t9          = _mm_sub_pd(t9, tF);
 629     t10         = _mm_sub_pd(t10, tH);
 630     t11         = _mm_sub_pd(t11, tJ);
 631     t12         = _mm_sub_pd(t12, tL);
 632
 633     _mm_storeu_pd(ptrA,  t1);
 634     _mm_storeu_pd(ptrA+2, t2);
 635     _mm_storeu_pd(ptrA+4, t3);
 636     _mm_storeu_pd(ptrA+6, t4);
 637     _mm_storeu_pd(ptrA+8, t5);
 638     _mm_storeu_pd(ptrA+10, t6);
 639     _mm_storeu_pd(ptrB,  t7);
 640     _mm_storeu_pd(ptrB+2, t8);
 641     _mm_storeu_pd(ptrB+4, t9);
 642     _mm_storeu_pd(ptrB+6, t10);
 643     _mm_storeu_pd(ptrB+8, t11);
 644     _mm_storeu_pd(ptrB+10, t12);
 645 }
 646
 647
 648 static gmx_inline void gmx_simdcall
 649 gmx_mm_update_iforce_1atom_swizzle_pd(__m128d fix1, __m128d fiy1, __m128d fiz1,
 650                                       double * gmx_restrict fptr,
 651                                       double * gmx_restrict fshiftptr)
 652 {
 653     fix1 = _mm_hadd_pd(fix1, fiy1);
 654     fiz1 = _mm_hadd_pd(fiz1, fiz1);
 655
 656     _mm_storeu_pd( fptr, _mm_add_pd( _mm_loadu_pd(fptr), fix1 ));
 657     _mm_store_sd( fptr+2, _mm_add_sd( _mm_load_sd(fptr+2), fiz1 ));
 658
 659     _mm_storeu_pd( fshiftptr, _mm_add_pd( _mm_loadu_pd(fshiftptr), fix1 ));
 660     _mm_store_sd( fshiftptr+2, _mm_add_sd( _mm_load_sd(fshiftptr+2), fiz1 ));
 661 }
 662
 663
 664
 665 static gmx_inline void gmx_simdcall
 666 gmx_mm_update_iforce_3atom_swizzle_pd(__m128d fix1, __m128d fiy1, __m128d fiz1,
 667                                       __m128d fix2, __m128d fiy2, __m128d fiz2,
 668                                       __m128d fix3, __m128d fiy3, __m128d fiz3,
 669                                       double * gmx_restrict fptr,
 670                                       double * gmx_restrict fshiftptr)
 671 {
 672     __m128d t1, t2;
 673
 674     fix1 = _mm_hadd_pd(fix1, fiy1);
 675     fiz1 = _mm_hadd_pd(fiz1, fix2);
 676     fiy2 = _mm_hadd_pd(fiy2, fiz2);
 677     fix3 = _mm_hadd_pd(fix3, fiy3);
 678     fiz3 = _mm_hadd_pd(fiz3, fiz3);
 679
 680     _mm_storeu_pd( fptr, _mm_add_pd( _mm_loadu_pd(fptr), fix1 ));
 681     _mm_storeu_pd( fptr+2, _mm_add_pd( _mm_loadu_pd(fptr+2), fiz1 ));
 682     _mm_storeu_pd( fptr+4, _mm_add_pd( _mm_loadu_pd(fptr+4), fiy2 ));
 683     _mm_storeu_pd( fptr+6, _mm_add_pd( _mm_loadu_pd(fptr+6), fix3 ));
 684     _mm_store_sd( fptr+8, _mm_add_sd( _mm_load_sd(fptr+8), fiz3 ));
 685
 686     fix1 = _mm_add_pd(fix1, fix3);
 687     t1   = _mm_shuffle_pd(fiz1, fiy2, _MM_SHUFFLE2(0, 1));
 688     fix1 = _mm_add_pd(fix1, t1); /* x and y sums */
 689
 690     t2   = _mm_shuffle_pd(fiy2, fiy2, _MM_SHUFFLE2(1, 1));
 691     fiz1 = _mm_add_sd(fiz1, fiz3);
 692     fiz1 = _mm_add_sd(fiz1, t2); /* z sum */
 693
 694     _mm_storeu_pd( fshiftptr, _mm_add_pd( _mm_loadu_pd(fshiftptr), fix1 ));
 695     _mm_store_sd( fshiftptr+2, _mm_add_sd( _mm_load_sd(fshiftptr+2), fiz1 ));
 696 }
 697
 698
 699 static gmx_inline void gmx_simdcall
 700 gmx_mm_update_iforce_4atom_swizzle_pd(__m128d fix1, __m128d fiy1, __m128d fiz1,
 701                                       __m128d fix2, __m128d fiy2, __m128d fiz2,
 702                                       __m128d fix3, __m128d fiy3, __m128d fiz3,
 703                                       __m128d fix4, __m128d fiy4, __m128d fiz4,
 704                                       double * gmx_restrict fptr,
 705                                       double * gmx_restrict fshiftptr)
 706 {
 707     __m128d t1, t2;
 708
 709     fix1 = _mm_hadd_pd(fix1, fiy1);
 710     fiz1 = _mm_hadd_pd(fiz1, fix2);
 711     fiy2 = _mm_hadd_pd(fiy2, fiz2);
 712     fix3 = _mm_hadd_pd(fix3, fiy3);
 713     fiz3 = _mm_hadd_pd(fiz3, fix4);
 714     fiy4 = _mm_hadd_pd(fiy4, fiz4);
 715
 716     _mm_storeu_pd( fptr, _mm_add_pd( _mm_loadu_pd(fptr),       fix1 ));
 717     _mm_storeu_pd( fptr+2, _mm_add_pd( _mm_loadu_pd(fptr+2),   fiz1 ));
 718     _mm_storeu_pd( fptr+4, _mm_add_pd( _mm_loadu_pd(fptr+4),   fiy2 ));
 719     _mm_storeu_pd( fptr+6, _mm_add_pd( _mm_loadu_pd(fptr+6),   fix3 ));
 720     _mm_storeu_pd( fptr+8, _mm_add_pd( _mm_loadu_pd(fptr+8),   fiz3 ));
 721     _mm_storeu_pd( fptr+10, _mm_add_pd( _mm_loadu_pd(fptr+10), fiy4 ));
 722
 723     t1   = _mm_shuffle_pd(fiz1, fiy2, _MM_SHUFFLE2(0, 1));
 724     fix1 = _mm_add_pd(fix1, t1);
 725     t2   = _mm_shuffle_pd(fiz3, fiy4, _MM_SHUFFLE2(0, 1));
 726     fix3 = _mm_add_pd(fix3, t2);
 727     fix1 = _mm_add_pd(fix1, fix3); /* x and y sums */
 728
 729     fiz1 = _mm_add_sd(fiz1, _mm_unpackhi_pd(fiy2, fiy2));
 730     fiz3 = _mm_add_sd(fiz3, _mm_unpackhi_pd(fiy4, fiy4));
 731     fiz1 = _mm_add_sd(fiz1, fiz3); /* z sum */
 732
 733     _mm_storeu_pd( fshiftptr, _mm_add_pd( _mm_loadu_pd(fshiftptr), fix1 ));
 734     _mm_store_sd( fshiftptr+2, _mm_add_sd( _mm_load_sd(fshiftptr+2), fiz1 ));
 735 }
 736
 737
 738 static gmx_inline void gmx_simdcall
 739 gmx_mm_update_1pot_pd(__m128d pot1, double * gmx_restrict ptrA)
 740 {
 741     pot1 = _mm_hadd_pd(pot1, pot1);
 742     _mm_store_sd(ptrA, _mm_add_sd(pot1, _mm_load_sd(ptrA)));
 743 }
 744
 745 static gmx_inline void gmx_simdcall
 746 gmx_mm_update_2pot_pd(__m128d pot1, double * gmx_restrict ptrA,
 747                       __m128d pot2, double * gmx_restrict ptrB)
 748 {
 749     pot1 = _mm_hadd_pd(pot1, pot2);
 750     pot2 = _mm_unpackhi_pd(pot1, pot1);
 751
 752     _mm_store_sd(ptrA, _mm_add_sd(pot1, _mm_load_sd(ptrA)));
 753     _mm_store_sd(ptrB, _mm_add_sd(pot2, _mm_load_sd(ptrB)));
 754 }
 755
 756
 757 #endif /* _kernelutil_x86_sse4_1_double_h_ */