-/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+/*
+ * 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.
*
- * This source code is part of
- *
- * G R O M A C S
+ * 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.
*
- * GROningen MAchine for Chemical Simulations
+ * 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.
*
- * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2012, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
-
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
- * of the License, or (at your option) any later version.
+ * 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, 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 www.gromacs.org.
+ * 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 papers on the package - you can find them in the top README file.
- *
- * For more info, check our website at http://www.gromacs.org
+ * the research papers on the package. Check out http://www.gromacs.org.
*/
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "gmxpre.h"
-#include <math.h>
-#include <string.h>
-#include "sysstuff.h"
-#include "smalloc.h"
-#include "macros.h"
-#include "maths.h"
-#include "vec.h"
-#include "pbc.h"
-#include "nbnxn_consts.h"
-#include "nbnxn_internal.h"
-#include "nbnxn_atomdata.h"
#include "nbnxn_search.h"
-#include "gmx_cyclecounter.h"
-#include "gmxfio.h"
-#include "gmx_omp_nthreads.h"
-#include "nrnb.h"
+#include "config.h"
-/* Pair search box lower and upper corner in x,y,z.
- * Store this in 4 iso 3 reals, which is useful with SSE.
- * To avoid complicating the code we also use 4 without SSE.
- */
-#define NNBSBB_C 4
-#define NNBSBB_B (2*NNBSBB_C)
-/* Pair search box lower and upper bound in z only. */
-#define NNBSBB_D 2
-/* Pair search box lower and upper corner x,y,z indices */
-#define BBL_X 0
-#define BBL_Y 1
-#define BBL_Z 2
-#define BBU_X 4
-#define BBU_Y 5
-#define BBU_Z 6
-
-
-#ifdef NBNXN_SEARCH_BB_SSE
-/* We use SSE or AVX-128bit for bounding box calculations */
+#include <assert.h>
+#include <math.h>
+#include <string.h>
-#ifndef GMX_DOUBLE
-/* Single precision BBs + coordinates, we can also load coordinates using SSE */
-#define NBNXN_SEARCH_SSE_SINGLE
+#include "gromacs/legacyheaders/gmx_omp_nthreads.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/legacyheaders/nrnb.h"
+#include "gromacs/legacyheaders/ns.h"
+#include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/math/utilities.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/mdlib/nb_verlet.h"
+#include "gromacs/mdlib/nbnxn_atomdata.h"
+#include "gromacs/mdlib/nbnxn_consts.h"
+#include "gromacs/pbcutil/ishift.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/utility/smalloc.h"
+
+/* nbnxn_internal.h included gromacs/simd/macros.h */
+#include "gromacs/mdlib/nbnxn_internal.h"
+#ifdef GMX_SIMD
+#include "gromacs/simd/vector_operations.h"
#endif
-/* Include basic SSE2 stuff */
-#include <emmintrin.h>
+#ifdef NBNXN_SEARCH_BB_SIMD4
+/* Always use 4-wide SIMD for bounding box calculations */
+
+# ifndef GMX_DOUBLE
+/* Single precision BBs + coordinates, we can also load coordinates with SIMD */
+# define NBNXN_SEARCH_SIMD4_FLOAT_X_BB
+# endif
-#if defined NBNXN_SEARCH_SSE_SINGLE && (GPU_NSUBCELL == 4 || GPU_NSUBCELL == 8)
+# if defined NBNXN_SEARCH_SIMD4_FLOAT_X_BB && (GPU_NSUBCELL == 4 || GPU_NSUBCELL == 8)
/* Store bounding boxes with x, y and z coordinates in packs of 4 */
-#define NBNXN_PBB_SSE
-#endif
+# define NBNXN_PBB_SIMD4
+# endif
-/* The width of SSE/AVX128 with single precision for bounding boxes with GPU.
- * Here AVX-256 turns out to be slightly slower than AVX-128.
+/* The packed bounding box coordinate stride is always set to 4.
+ * With AVX we could use 8, but that turns out not to be faster.
*/
-#define STRIDE_PBB 4
-#define STRIDE_PBB_2LOG 2
+# define STRIDE_PBB 4
+# define STRIDE_PBB_2LOG 2
-#endif /* NBNXN_SEARCH_BB_SSE */
+#endif /* NBNXN_SEARCH_BB_SIMD4 */
#ifdef GMX_NBNXN_SIMD
#define X_IND_CJ_J8(cj) ((cj)*STRIDE_P8)
/* The j-cluster size is matched to the SIMD width */
-#if GMX_NBNXN_SIMD_BITWIDTH == 128
-#ifdef GMX_DOUBLE
+#if GMX_SIMD_REAL_WIDTH == 2
#define CI_TO_CJ_SIMD_4XN(ci) CI_TO_CJ_J2(ci)
#define X_IND_CI_SIMD_4XN(ci) X_IND_CI_J2(ci)
#define X_IND_CJ_SIMD_4XN(cj) X_IND_CJ_J2(cj)
#else
-#define CI_TO_CJ_SIMD_4XN(ci) CI_TO_CJ_J4(ci)
-#define X_IND_CI_SIMD_4XN(ci) X_IND_CI_J4(ci)
-#define X_IND_CJ_SIMD_4XN(cj) X_IND_CJ_J4(cj)
-#endif
-#else
-#if GMX_NBNXN_SIMD_BITWIDTH == 256
-#ifdef GMX_DOUBLE
+#if GMX_SIMD_REAL_WIDTH == 4
#define CI_TO_CJ_SIMD_4XN(ci) CI_TO_CJ_J4(ci)
#define X_IND_CI_SIMD_4XN(ci) X_IND_CI_J4(ci)
#define X_IND_CJ_SIMD_4XN(cj) X_IND_CJ_J4(cj)
#else
+#if GMX_SIMD_REAL_WIDTH == 8
#define CI_TO_CJ_SIMD_4XN(ci) CI_TO_CJ_J8(ci)
#define X_IND_CI_SIMD_4XN(ci) X_IND_CI_J8(ci)
#define X_IND_CJ_SIMD_4XN(cj) X_IND_CJ_J8(cj)
#define CI_TO_CJ_SIMD_2XNN(ci) CI_TO_CJ_J4(ci)
#define X_IND_CI_SIMD_2XNN(ci) X_IND_CI_J4(ci)
#define X_IND_CJ_SIMD_2XNN(cj) X_IND_CJ_J4(cj)
-#endif
#else
-#error "unsupported GMX_NBNXN_SIMD_WIDTH"
+#if GMX_SIMD_REAL_WIDTH == 16
+#define CI_TO_CJ_SIMD_2XNN(ci) CI_TO_CJ_J8(ci)
+#define X_IND_CI_SIMD_2XNN(ci) X_IND_CI_J8(ci)
+#define X_IND_CJ_SIMD_2XNN(cj) X_IND_CJ_J8(cj)
+#else
+#error "unsupported GMX_SIMD_REAL_WIDTH"
+#endif
+#endif
#endif
#endif
#endif /* GMX_NBNXN_SIMD */
-/* Interaction masks for 4xN atom interactions.
- * Bit i*CJ_SIZE + j tells if atom i and j interact.
- */
-/* All interaction mask is the same for all kernels */
-#define NBNXN_INT_MASK_ALL 0xffffffff
-/* 4x4 kernel diagonal mask */
-#define NBNXN_INT_MASK_DIAG 0x08ce
-/* 4x2 kernel diagonal masks */
-#define NBNXN_INT_MASK_DIAG_J2_0 0x0002
-#define NBNXN_INT_MASK_DIAG_J2_1 0x002F
-/* 4x8 kernel diagonal masks */
-#define NBNXN_INT_MASK_DIAG_J8_0 0xf0f8fcfe
-#define NBNXN_INT_MASK_DIAG_J8_1 0x0080c0e0
-
-
-#ifdef NBNXN_SEARCH_BB_SSE
+#ifdef NBNXN_SEARCH_BB_SIMD4
/* Store bounding boxes corners as quadruplets: xxxxyyyyzzzz */
#define NBNXN_BBXXXX
/* Size of bounding box corners quadruplet */
{
int i;
- for(i=0; i<enbsCCnr; i++)
+ for (i = 0; i < enbsCCnr; i++)
{
cc[i].count = 0;
cc[i].c = 0;
return (double)cc->c*1e-6/cc->count;
}
-static void nbs_cycle_print(FILE *fp,const nbnxn_search_t nbs)
+static void nbs_cycle_print(FILE *fp, const nbnxn_search_t nbs)
{
int n;
int t;
- fprintf(fp,"\n");
- fprintf(fp,"ns %4d grid %4.1f search %4.1f red.f %5.3f",
+ fprintf(fp, "\n");
+ fprintf(fp, "ns %4d grid %4.1f search %4.1f red.f %5.3f",
nbs->cc[enbsCCgrid].count,
Mcyc_av(&nbs->cc[enbsCCgrid]),
Mcyc_av(&nbs->cc[enbsCCsearch]),
{
if (nbs->cc[enbsCCcombine].count > 0)
{
- fprintf(fp," comb %5.2f",
+ fprintf(fp, " comb %5.2f",
Mcyc_av(&nbs->cc[enbsCCcombine]));
}
- fprintf(fp," s. th");
- for(t=0; t<nbs->nthread_max; t++)
+ fprintf(fp, " s. th");
+ for (t = 0; t < nbs->nthread_max; t++)
{
- fprintf(fp," %4.1f",
+ fprintf(fp, " %4.1f",
Mcyc_av(&nbs->work[t].cc[enbsCCsearch]));
}
}
- fprintf(fp,"\n");
+ fprintf(fp, "\n");
}
static void nbnxn_grid_init(nbnxn_grid_t * grid)
}
if ((1<<log2) != n)
{
- gmx_fatal(FARGS,"nbnxn na_c (%d) is not a power of 2",n);
+ gmx_fatal(FARGS, "nbnxn na_c (%d) is not a power of 2", n);
}
return log2;
{
switch (nb_kernel_type)
{
- case nbnxnk4x4_PlainC:
- case nbnxnk4xN_SIMD_4xN:
- case nbnxnk4xN_SIMD_2xNN:
- return NBNXN_CPU_CLUSTER_I_SIZE;
- case nbnxnk8x8x8_CUDA:
- case nbnxnk8x8x8_PlainC:
- /* The cluster size for super/sub lists is only set here.
- * Any value should work for the pair-search and atomdata code.
- * The kernels, of course, might require a particular value.
- */
- return NBNXN_GPU_CLUSTER_SIZE;
- default:
- gmx_incons("unknown kernel type");
+ case nbnxnk4x4_PlainC:
+ case nbnxnk4xN_SIMD_4xN:
+ case nbnxnk4xN_SIMD_2xNN:
+ return NBNXN_CPU_CLUSTER_I_SIZE;
+ case nbnxnk8x8x8_CUDA:
+ case nbnxnk8x8x8_PlainC:
+ /* The cluster size for super/sub lists is only set here.
+ * Any value should work for the pair-search and atomdata code.
+ * The kernels, of course, might require a particular value.
+ */
+ return NBNXN_GPU_CLUSTER_SIZE;
+ default:
+ gmx_incons("unknown kernel type");
}
return 0;
int nbnxn_kernel_to_cj_size(int nb_kernel_type)
{
- int nbnxn_simd_width=0;
- int cj_size=0;
+ int nbnxn_simd_width = 0;
+ int cj_size = 0;
#ifdef GMX_NBNXN_SIMD
- nbnxn_simd_width = GMX_NBNXN_SIMD_BITWIDTH/(sizeof(real)*8);
+ nbnxn_simd_width = GMX_SIMD_REAL_WIDTH;
#endif
switch (nb_kernel_type)
{
- case nbnxnk4x4_PlainC:
- cj_size = NBNXN_CPU_CLUSTER_I_SIZE;
- break;
- case nbnxnk4xN_SIMD_4xN:
- cj_size = nbnxn_simd_width;
- break;
- case nbnxnk4xN_SIMD_2xNN:
- cj_size = nbnxn_simd_width/2;
- break;
- case nbnxnk8x8x8_CUDA:
- case nbnxnk8x8x8_PlainC:
- cj_size = nbnxn_kernel_to_ci_size(nb_kernel_type);
- break;
- default:
- gmx_incons("unknown kernel type");
+ case nbnxnk4x4_PlainC:
+ cj_size = NBNXN_CPU_CLUSTER_I_SIZE;
+ break;
+ case nbnxnk4xN_SIMD_4xN:
+ cj_size = nbnxn_simd_width;
+ break;
+ case nbnxnk4xN_SIMD_2xNN:
+ cj_size = nbnxn_simd_width/2;
+ break;
+ case nbnxnk8x8x8_CUDA:
+ case nbnxnk8x8x8_PlainC:
+ cj_size = nbnxn_kernel_to_ci_size(nb_kernel_type);
+ break;
+ default:
+ gmx_incons("unknown kernel type");
}
return cj_size;
}
-static int ci_to_cj(int na_cj_2log,int ci)
+static int ci_to_cj(int na_cj_2log, int ci)
{
switch (na_cj_2log)
{
- case 2: return ci; break;
- case 1: return (ci<<1); break;
- case 3: return (ci>>1); break;
+ case 2: return ci; break;
+ case 1: return (ci<<1); break;
+ case 3: return (ci>>1); break;
}
return 0;
switch (nb_kernel_type)
{
- case nbnxnk8x8x8_CUDA:
- case nbnxnk8x8x8_PlainC:
- return FALSE;
+ case nbnxnk8x8x8_CUDA:
+ case nbnxnk8x8x8_PlainC:
+ return FALSE;
- case nbnxnk4x4_PlainC:
- case nbnxnk4xN_SIMD_4xN:
- case nbnxnk4xN_SIMD_2xNN:
- return TRUE;
+ case nbnxnk4x4_PlainC:
+ case nbnxnk4xN_SIMD_4xN:
+ case nbnxnk4xN_SIMD_2xNN:
+ return TRUE;
- default:
- gmx_incons("Invalid nonbonded kernel type passed!");
- return FALSE;
+ default:
+ gmx_incons("Invalid nonbonded kernel type passed!");
+ return FALSE;
}
}
-void nbnxn_init_search(nbnxn_search_t * nbs_ptr,
- ivec *n_dd_cells,
+/* Initializes a single nbnxn_pairlist_t data structure */
+static void nbnxn_init_pairlist_fep(t_nblist *nl)
+{
+ nl->type = GMX_NBLIST_INTERACTION_FREE_ENERGY;
+ nl->igeometry = GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE;
+ /* The interaction functions are set in the free energy kernel fuction */
+ nl->ivdw = -1;
+ nl->ivdwmod = -1;
+ nl->ielec = -1;
+ nl->ielecmod = -1;
+
+ nl->maxnri = 0;
+ nl->maxnrj = 0;
+ nl->nri = 0;
+ nl->nrj = 0;
+ nl->iinr = NULL;
+ nl->gid = NULL;
+ nl->shift = NULL;
+ nl->jindex = NULL;
+ nl->jjnr = NULL;
+ nl->excl_fep = NULL;
+
+}
+
+void nbnxn_init_search(nbnxn_search_t * nbs_ptr,
+ ivec *n_dd_cells,
gmx_domdec_zones_t *zones,
- int nthread_max)
+ gmx_bool bFEP,
+ int nthread_max)
{
nbnxn_search_t nbs;
- int d,g,t;
+ int d, g, t;
- snew(nbs,1);
+ snew(nbs, 1);
*nbs_ptr = nbs;
+ nbs->bFEP = bFEP;
+
nbs->DomDec = (n_dd_cells != NULL);
clear_ivec(nbs->dd_dim);
{
nbs->zones = zones;
- for(d=0; d<DIM; d++)
+ for (d = 0; d < DIM; d++)
{
if ((*n_dd_cells)[d] > 1)
{
}
}
- snew(nbs->grid,nbs->ngrid);
- for(g=0; g<nbs->ngrid; g++)
+ snew(nbs->grid, nbs->ngrid);
+ for (g = 0; g < nbs->ngrid; g++)
{
nbnxn_grid_init(&nbs->grid[g]);
}
nbs->nthread_max = nthread_max;
/* Initialize the work data structures for each thread */
- snew(nbs->work,nbs->nthread_max);
- for(t=0; t<nbs->nthread_max; t++)
+ snew(nbs->work, nbs->nthread_max);
+ for (t = 0; t < nbs->nthread_max; t++)
{
nbs->work[t].cxy_na = NULL;
nbs->work[t].cxy_na_nalloc = 0;
nbs->work[t].sort_work = NULL;
nbs->work[t].sort_work_nalloc = 0;
+
+ snew(nbs->work[t].nbl_fep, 1);
+ nbnxn_init_pairlist_fep(nbs->work[t].nbl_fep);
}
/* Initialize detailed nbsearch cycle counting */
nbs->print_cycles = (getenv("GMX_NBNXN_CYCLE") != 0);
nbs->search_count = 0;
nbs_cycle_clear(nbs->cc);
- for(t=0; t<nbs->nthread_max; t++)
+ for (t = 0; t < nbs->nthread_max; t++)
{
nbs_cycle_clear(nbs->work[t].cc);
}
}
-static real grid_atom_density(int n,rvec corner0,rvec corner1)
+static real grid_atom_density(int n, rvec corner0, rvec corner1)
{
rvec size;
- rvec_sub(corner1,corner0,size);
+ if (n == 0)
+ {
+ /* To avoid zero density we use a minimum of 1 atom */
+ n = 1;
+ }
+
+ rvec_sub(corner1, corner0, size);
return n/(size[XX]*size[YY]*size[ZZ]);
}
static int set_grid_size_xy(const nbnxn_search_t nbs,
nbnxn_grid_t *grid,
int dd_zone,
- int n,rvec corner0,rvec corner1,
- real atom_density,
- int XFormat)
+ int n, rvec corner0, rvec corner1,
+ real atom_density)
{
rvec size;
int na_c;
- real adens,tlen,tlen_x,tlen_y,nc_max;
+ real adens, tlen, tlen_x, tlen_y, nc_max;
int t;
- rvec_sub(corner1,corner0,size);
+ rvec_sub(corner1, corner0, size);
if (n > grid->na_sc)
{
+ assert(atom_density > 0);
+
/* target cell length */
if (grid->bSimple)
{
/* To minimize the zero interactions, we should make
* the largest of the i/j cell cubic.
*/
- na_c = max(grid->na_c,grid->na_cj);
+ na_c = max(grid->na_c, grid->na_cj);
/* Approximately cubic cells */
- tlen = pow(na_c/atom_density,1.0/3.0);
+ tlen = pow(na_c/atom_density, 1.0/3.0);
tlen_x = tlen;
tlen_y = tlen;
}
else
{
/* Approximately cubic sub cells */
- tlen = pow(grid->na_c/atom_density,1.0/3.0);
+ tlen = pow(grid->na_c/atom_density, 1.0/3.0);
tlen_x = tlen*GPU_NSUBCELL_X;
tlen_y = tlen*GPU_NSUBCELL_Y;
}
* in the nbsist when the fixed cell dimensions (x,y) are
* larger than the variable one (z) than the other way around.
*/
- grid->ncx = max(1,(int)(size[XX]/tlen_x));
- grid->ncy = max(1,(int)(size[YY]/tlen_y));
+ grid->ncx = max(1, (int)(size[XX]/tlen_x));
+ grid->ncy = max(1, (int)(size[YY]/tlen_y));
}
else
{
grid->ncy = 1;
}
- grid->sx = size[XX]/grid->ncx;
- grid->sy = size[YY]/grid->ncy;
+ grid->sx = size[XX]/grid->ncx;
+ grid->sy = size[YY]/grid->ncy;
grid->inv_sx = 1/grid->sx;
grid->inv_sy = 1/grid->sy;
if (grid->ncx*grid->ncy+1 > grid->cxy_nalloc)
{
grid->cxy_nalloc = over_alloc_large(grid->ncx*grid->ncy+1);
- srenew(grid->cxy_na,grid->cxy_nalloc);
- srenew(grid->cxy_ind,grid->cxy_nalloc+1);
+ srenew(grid->cxy_na, grid->cxy_nalloc);
+ srenew(grid->cxy_ind, grid->cxy_nalloc+1);
}
- for(t=0; t<nbs->nthread_max; t++)
+ for (t = 0; t < nbs->nthread_max; t++)
{
if (grid->ncx*grid->ncy+1 > nbs->work[t].cxy_na_nalloc)
{
nbs->work[t].cxy_na_nalloc = over_alloc_large(grid->ncx*grid->ncy+1);
- srenew(nbs->work[t].cxy_na,nbs->work[t].cxy_na_nalloc);
+ srenew(nbs->work[t].cxy_na, nbs->work[t].cxy_na_nalloc);
}
}
if (nc_max > grid->nc_nalloc)
{
- int bb_nalloc;
-
grid->nc_nalloc = over_alloc_large(nc_max);
- srenew(grid->nsubc,grid->nc_nalloc);
- srenew(grid->bbcz,grid->nc_nalloc*NNBSBB_D);
-#ifdef NBNXN_PBB_SSE
- bb_nalloc = grid->nc_nalloc*GPU_NSUBCELL/STRIDE_PBB*NNBSBB_XXXX;
-#else
- bb_nalloc = grid->nc_nalloc*GPU_NSUBCELL*NNBSBB_B;
-#endif
+ srenew(grid->nsubc, grid->nc_nalloc);
+ srenew(grid->bbcz, grid->nc_nalloc*NNBSBB_D);
+
sfree_aligned(grid->bb);
/* This snew also zeros the contents, this avoid possible
- * floating exceptions in SSE with the unused bb elements.
+ * floating exceptions in SIMD with the unused bb elements.
*/
- snew_aligned(grid->bb,bb_nalloc,16);
+ if (grid->bSimple)
+ {
+ snew_aligned(grid->bb, grid->nc_nalloc, 16);
+ }
+ else
+ {
+#ifdef NBNXN_BBXXXX
+ int pbb_nalloc;
+
+ pbb_nalloc = grid->nc_nalloc*GPU_NSUBCELL/STRIDE_PBB*NNBSBB_XXXX;
+ snew_aligned(grid->pbb, pbb_nalloc, 16);
+#else
+ snew_aligned(grid->bb, grid->nc_nalloc*GPU_NSUBCELL, 16);
+#endif
+ }
if (grid->bSimple)
{
else
{
sfree_aligned(grid->bbj);
- snew_aligned(grid->bbj,bb_nalloc*grid->na_c/grid->na_cj,16);
+ snew_aligned(grid->bbj, grid->nc_nalloc*grid->na_c/grid->na_cj, 16);
}
}
- srenew(grid->flags,grid->nc_nalloc);
+ srenew(grid->flags, grid->nc_nalloc);
+ if (nbs->bFEP)
+ {
+ srenew(grid->fep, grid->nc_nalloc*grid->na_sc/grid->na_c);
+ }
}
- copy_rvec(corner0,grid->c0);
- copy_rvec(corner1,grid->c1);
+ copy_rvec(corner0, grid->c0);
+ copy_rvec(corner1, grid->c1);
return nc_max;
}
/* Sort particle index a on coordinates x along dim.
* Backwards tells if we want decreasing iso increasing coordinates.
* h0 is the minimum of the coordinate range.
- * invh is the inverse hole spacing.
- * nsort, the theortical hole limit, is only used for debugging.
+ * invh is the 1/length of the sorting range.
+ * n_per_h (>=n) is the expected average number of particles per 1/invh
* sort is the sorting work array.
+ * sort should have a size of at least n_per_h*SORT_GRID_OVERSIZE + n,
+ * or easier, allocate at least n*SGSF elements.
*/
-static void sort_atoms(int dim,gmx_bool Backwards,
- int *a,int n,rvec *x,
- real h0,real invh,int nsort,int *sort)
+static void sort_atoms(int dim, gmx_bool Backwards,
+ int gmx_unused dd_zone,
+ int *a, int n, rvec *x,
+ real h0, real invh, int n_per_h,
+ int *sort)
{
- int i,c;
- int zi,zim,zi_min,zi_max;
- int cp,tmp;
+ int nsort, i, c;
+ int zi, zim, zi_min, zi_max;
+ int cp, tmp;
if (n <= 1)
{
return;
}
+#ifndef NDEBUG
+ if (n > n_per_h)
+ {
+ gmx_incons("n > n_per_h");
+ }
+#endif
+
+ /* Transform the inverse range height into the inverse hole height */
+ invh *= n_per_h*SORT_GRID_OVERSIZE;
+
+ /* Set nsort to the maximum possible number of holes used.
+ * In worst case all n elements end up in the last bin.
+ */
+ nsort = n_per_h*SORT_GRID_OVERSIZE + n;
+
/* Determine the index range used, so we can limit it for the second pass */
zi_min = INT_MAX;
zi_max = -1;
/* Sort the particles using a simple index sort */
- for(i=0; i<n; i++)
+ for (i = 0; i < n; i++)
{
/* The cast takes care of float-point rounding effects below zero.
* This code assumes particles are less than 1/SORT_GRID_OVERSIZE
*/
zi = (int)((x[a[i]][dim] - h0)*invh);
-#ifdef DEBUG_NBNXN_GRIDDING
- if (zi < 0 || zi >= nsort)
+#ifndef NDEBUG
+ /* As we can have rounding effect, we use > iso >= here */
+ if (zi < 0 || (dd_zone == 0 && zi > n_per_h*SORT_GRID_OVERSIZE))
{
- gmx_fatal(FARGS,"(int)((x[%d][%c]=%f - %f)*%f) = %d, not in 0 - %d\n",
- a[i],'x'+dim,x[a[i]][dim],h0,invh,zi,nsort);
+ gmx_fatal(FARGS, "(int)((x[%d][%c]=%f - %f)*%f) = %d, not in 0 - %d*%d\n",
+ a[i], 'x'+dim, x[a[i]][dim], h0, invh, zi,
+ n_per_h, SORT_GRID_OVERSIZE);
}
#endif
+ /* In a non-local domain, particles communcated for bonded interactions
+ * can be far beyond the grid size, which is set by the non-bonded
+ * cut-off distance. We sort such particles into the last cell.
+ */
+ if (zi > n_per_h*SORT_GRID_OVERSIZE)
+ {
+ zi = n_per_h*SORT_GRID_OVERSIZE;
+ }
+
/* Ideally this particle should go in sort cell zi,
* but that might already be in use,
* in that case find the first empty cell higher up
if (sort[zi] < 0)
{
sort[zi] = a[i];
- zi_min = min(zi_min,zi);
- zi_max = max(zi_max,zi);
+ zi_min = min(zi_min, zi);
+ zi_max = max(zi_max, zi);
}
else
{
* well-defined output order, independent of input order
* to ensure binary reproducibility after restarts.
*/
- while(sort[zi] >= 0 && ( x[a[i]][dim] > x[sort[zi]][dim] ||
- (x[a[i]][dim] == x[sort[zi]][dim] &&
- a[i] > sort[zi])))
+ while (sort[zi] >= 0 && ( x[a[i]][dim] > x[sort[zi]][dim] ||
+ (x[a[i]][dim] == x[sort[zi]][dim] &&
+ a[i] > sort[zi])))
{
zi++;
}
if (sort[zi] >= 0)
{
/* Shift all elements by one slot until we find an empty slot */
- cp = sort[zi];
+ cp = sort[zi];
zim = zi + 1;
while (sort[zim] >= 0)
{
- tmp = sort[zim];
+ tmp = sort[zim];
sort[zim] = cp;
- cp = tmp;
+ cp = tmp;
zim++;
}
sort[zim] = cp;
- zi_max = max(zi_max,zim);
+ zi_max = max(zi_max, zim);
}
sort[zi] = a[i];
- zi_max = max(zi_max,zi);
+ zi_max = max(zi_max, zi);
}
}
c = 0;
if (!Backwards)
{
- for(zi=0; zi<nsort; zi++)
+ for (zi = 0; zi < nsort; zi++)
{
if (sort[zi] >= 0)
{
- a[c++] = sort[zi];
+ a[c++] = sort[zi];
sort[zi] = -1;
}
}
}
else
{
- for(zi=zi_max; zi>=zi_min; zi--)
+ for (zi = zi_max; zi >= zi_min; zi--)
{
if (sort[zi] >= 0)
{
- a[c++] = sort[zi];
+ a[c++] = sort[zi];
sort[zi] = -1;
}
}
#endif
/* Coordinate order x,y,z, bb order xyz0 */
-static void calc_bounding_box(int na,int stride,const real *x,float *bb)
+static void calc_bounding_box(int na, int stride, const real *x, nbnxn_bb_t *bb)
{
- int i,j;
- real xl,xh,yl,yh,zl,zh;
+ int i, j;
+ real xl, xh, yl, yh, zl, zh;
- i = 0;
+ i = 0;
xl = x[i+XX];
xh = x[i+XX];
yl = x[i+YY];
zl = x[i+ZZ];
zh = x[i+ZZ];
i += stride;
- for(j=1; j<na; j++)
- {
- xl = min(xl,x[i+XX]);
- xh = max(xh,x[i+XX]);
- yl = min(yl,x[i+YY]);
- yh = max(yh,x[i+YY]);
- zl = min(zl,x[i+ZZ]);
- zh = max(zh,x[i+ZZ]);
+ for (j = 1; j < na; j++)
+ {
+ xl = min(xl, x[i+XX]);
+ xh = max(xh, x[i+XX]);
+ yl = min(yl, x[i+YY]);
+ yh = max(yh, x[i+YY]);
+ zl = min(zl, x[i+ZZ]);
+ zh = max(zh, x[i+ZZ]);
i += stride;
}
/* Note: possible double to float conversion here */
- bb[BBL_X] = R2F_D(xl);
- bb[BBL_Y] = R2F_D(yl);
- bb[BBL_Z] = R2F_D(zl);
- bb[BBU_X] = R2F_U(xh);
- bb[BBU_Y] = R2F_U(yh);
- bb[BBU_Z] = R2F_U(zh);
+ bb->lower[BB_X] = R2F_D(xl);
+ bb->lower[BB_Y] = R2F_D(yl);
+ bb->lower[BB_Z] = R2F_D(zl);
+ bb->upper[BB_X] = R2F_U(xh);
+ bb->upper[BB_Y] = R2F_U(yh);
+ bb->upper[BB_Z] = R2F_U(zh);
}
/* Packed coordinates, bb order xyz0 */
-static void calc_bounding_box_x_x4(int na,const real *x,float *bb)
+static void calc_bounding_box_x_x4(int na, const real *x, nbnxn_bb_t *bb)
{
int j;
- real xl,xh,yl,yh,zl,zh;
+ real xl, xh, yl, yh, zl, zh;
xl = x[XX*PACK_X4];
xh = x[XX*PACK_X4];
yh = x[YY*PACK_X4];
zl = x[ZZ*PACK_X4];
zh = x[ZZ*PACK_X4];
- for(j=1; j<na; j++)
+ for (j = 1; j < na; j++)
{
- xl = min(xl,x[j+XX*PACK_X4]);
- xh = max(xh,x[j+XX*PACK_X4]);
- yl = min(yl,x[j+YY*PACK_X4]);
- yh = max(yh,x[j+YY*PACK_X4]);
- zl = min(zl,x[j+ZZ*PACK_X4]);
- zh = max(zh,x[j+ZZ*PACK_X4]);
+ xl = min(xl, x[j+XX*PACK_X4]);
+ xh = max(xh, x[j+XX*PACK_X4]);
+ yl = min(yl, x[j+YY*PACK_X4]);
+ yh = max(yh, x[j+YY*PACK_X4]);
+ zl = min(zl, x[j+ZZ*PACK_X4]);
+ zh = max(zh, x[j+ZZ*PACK_X4]);
}
/* Note: possible double to float conversion here */
- bb[BBL_X] = R2F_D(xl);
- bb[BBL_Y] = R2F_D(yl);
- bb[BBL_Z] = R2F_D(zl);
- bb[BBU_X] = R2F_U(xh);
- bb[BBU_Y] = R2F_U(yh);
- bb[BBU_Z] = R2F_U(zh);
+ bb->lower[BB_X] = R2F_D(xl);
+ bb->lower[BB_Y] = R2F_D(yl);
+ bb->lower[BB_Z] = R2F_D(zl);
+ bb->upper[BB_X] = R2F_U(xh);
+ bb->upper[BB_Y] = R2F_U(yh);
+ bb->upper[BB_Z] = R2F_U(zh);
}
/* Packed coordinates, bb order xyz0 */
-static void calc_bounding_box_x_x8(int na,const real *x,float *bb)
+static void calc_bounding_box_x_x8(int na, const real *x, nbnxn_bb_t *bb)
{
int j;
- real xl,xh,yl,yh,zl,zh;
+ real xl, xh, yl, yh, zl, zh;
xl = x[XX*PACK_X8];
xh = x[XX*PACK_X8];
yh = x[YY*PACK_X8];
zl = x[ZZ*PACK_X8];
zh = x[ZZ*PACK_X8];
- for(j=1; j<na; j++)
+ for (j = 1; j < na; j++)
{
- xl = min(xl,x[j+XX*PACK_X8]);
- xh = max(xh,x[j+XX*PACK_X8]);
- yl = min(yl,x[j+YY*PACK_X8]);
- yh = max(yh,x[j+YY*PACK_X8]);
- zl = min(zl,x[j+ZZ*PACK_X8]);
- zh = max(zh,x[j+ZZ*PACK_X8]);
+ xl = min(xl, x[j+XX*PACK_X8]);
+ xh = max(xh, x[j+XX*PACK_X8]);
+ yl = min(yl, x[j+YY*PACK_X8]);
+ yh = max(yh, x[j+YY*PACK_X8]);
+ zl = min(zl, x[j+ZZ*PACK_X8]);
+ zh = max(zh, x[j+ZZ*PACK_X8]);
}
/* Note: possible double to float conversion here */
- bb[BBL_X] = R2F_D(xl);
- bb[BBL_Y] = R2F_D(yl);
- bb[BBL_Z] = R2F_D(zl);
- bb[BBU_X] = R2F_U(xh);
- bb[BBU_Y] = R2F_U(yh);
- bb[BBU_Z] = R2F_U(zh);
+ bb->lower[BB_X] = R2F_D(xl);
+ bb->lower[BB_Y] = R2F_D(yl);
+ bb->lower[BB_Z] = R2F_D(zl);
+ bb->upper[BB_X] = R2F_U(xh);
+ bb->upper[BB_Y] = R2F_U(yh);
+ bb->upper[BB_Z] = R2F_U(zh);
}
-#ifdef NBNXN_SEARCH_BB_SSE
-
/* Packed coordinates, bb order xyz0 */
-static void calc_bounding_box_x_x4_halves(int na,const real *x,
- float *bb,float *bbj)
+static void calc_bounding_box_x_x4_halves(int na, const real *x,
+ nbnxn_bb_t *bb, nbnxn_bb_t *bbj)
{
- calc_bounding_box_x_x4(min(na,2),x,bbj);
+ calc_bounding_box_x_x4(min(na, 2), x, bbj);
if (na > 2)
{
- calc_bounding_box_x_x4(min(na-2,2),x+(PACK_X4>>1),bbj+NNBSBB_B);
+ calc_bounding_box_x_x4(min(na-2, 2), x+(PACK_X4>>1), bbj+1);
}
else
{
/* Set the "empty" bounding box to the same as the first one,
* so we don't need to treat special cases in the rest of the code.
*/
- _mm_store_ps(bbj+NNBSBB_B ,_mm_load_ps(bbj));
- _mm_store_ps(bbj+NNBSBB_B+NNBSBB_C,_mm_load_ps(bbj+NNBSBB_C));
+#ifdef NBNXN_SEARCH_BB_SIMD4
+ gmx_simd4_store_f(&bbj[1].lower[0], gmx_simd4_load_f(&bbj[0].lower[0]));
+ gmx_simd4_store_f(&bbj[1].upper[0], gmx_simd4_load_f(&bbj[0].upper[0]));
+#else
+ bbj[1] = bbj[0];
+#endif
}
- _mm_store_ps(bb ,_mm_min_ps(_mm_load_ps(bbj),
- _mm_load_ps(bbj+NNBSBB_B)));
- _mm_store_ps(bb+NNBSBB_C,_mm_max_ps(_mm_load_ps(bbj+NNBSBB_C),
- _mm_load_ps(bbj+NNBSBB_B+NNBSBB_C)));
+#ifdef NBNXN_SEARCH_BB_SIMD4
+ gmx_simd4_store_f(&bb->lower[0],
+ gmx_simd4_min_f(gmx_simd4_load_f(&bbj[0].lower[0]),
+ gmx_simd4_load_f(&bbj[1].lower[0])));
+ gmx_simd4_store_f(&bb->upper[0],
+ gmx_simd4_max_f(gmx_simd4_load_f(&bbj[0].upper[0]),
+ gmx_simd4_load_f(&bbj[1].upper[0])));
+#else
+ {
+ int i;
+
+ for (i = 0; i < NNBSBB_C; i++)
+ {
+ bb->lower[i] = min(bbj[0].lower[i], bbj[1].lower[i]);
+ bb->upper[i] = max(bbj[0].upper[i], bbj[1].upper[i]);
+ }
+ }
+#endif
}
+#ifdef NBNXN_SEARCH_BB_SIMD4
+
/* Coordinate order xyz, bb order xxxxyyyyzzzz */
-static void calc_bounding_box_xxxx(int na,int stride,const real *x,float *bb)
+static void calc_bounding_box_xxxx(int na, int stride, const real *x, float *bb)
{
- int i,j;
- real xl,xh,yl,yh,zl,zh;
+ int i, j;
+ real xl, xh, yl, yh, zl, zh;
- i = 0;
+ i = 0;
xl = x[i+XX];
xh = x[i+XX];
yl = x[i+YY];
zl = x[i+ZZ];
zh = x[i+ZZ];
i += stride;
- for(j=1; j<na; j++)
- {
- xl = min(xl,x[i+XX]);
- xh = max(xh,x[i+XX]);
- yl = min(yl,x[i+YY]);
- yh = max(yh,x[i+YY]);
- zl = min(zl,x[i+ZZ]);
- zh = max(zh,x[i+ZZ]);
+ for (j = 1; j < na; j++)
+ {
+ xl = min(xl, x[i+XX]);
+ xh = max(xh, x[i+XX]);
+ yl = min(yl, x[i+YY]);
+ yh = max(yh, x[i+YY]);
+ zl = min(zl, x[i+ZZ]);
+ zh = max(zh, x[i+ZZ]);
i += stride;
}
/* Note: possible double to float conversion here */
bb[5*STRIDE_PBB] = R2F_U(zh);
}
-#endif /* NBNXN_SEARCH_BB_SSE */
+#endif /* NBNXN_SEARCH_BB_SIMD4 */
-#ifdef NBNXN_SEARCH_SSE_SINGLE
+#ifdef NBNXN_SEARCH_SIMD4_FLOAT_X_BB
/* Coordinate order xyz?, bb order xyz0 */
-static void calc_bounding_box_sse(int na,const float *x,float *bb)
+static void calc_bounding_box_simd4(int na, const float *x, nbnxn_bb_t *bb)
{
- __m128 bb_0_SSE,bb_1_SSE;
- __m128 x_SSE;
+ gmx_simd4_float_t bb_0_S, bb_1_S;
+ gmx_simd4_float_t x_S;
- int i;
+ int i;
- bb_0_SSE = _mm_load_ps(x);
- bb_1_SSE = bb_0_SSE;
+ bb_0_S = gmx_simd4_load_f(x);
+ bb_1_S = bb_0_S;
- for(i=1; i<na; i++)
+ for (i = 1; i < na; i++)
{
- x_SSE = _mm_load_ps(x+i*NNBSBB_C);
- bb_0_SSE = _mm_min_ps(bb_0_SSE,x_SSE);
- bb_1_SSE = _mm_max_ps(bb_1_SSE,x_SSE);
+ x_S = gmx_simd4_load_f(x+i*NNBSBB_C);
+ bb_0_S = gmx_simd4_min_f(bb_0_S, x_S);
+ bb_1_S = gmx_simd4_max_f(bb_1_S, x_S);
}
- _mm_store_ps(bb ,bb_0_SSE);
- _mm_store_ps(bb+4,bb_1_SSE);
+ gmx_simd4_store_f(&bb->lower[0], bb_0_S);
+ gmx_simd4_store_f(&bb->upper[0], bb_1_S);
}
/* Coordinate order xyz?, bb order xxxxyyyyzzzz */
-static void calc_bounding_box_xxxx_sse(int na,const float *x,
- float *bb_work,
- real *bb)
+static void calc_bounding_box_xxxx_simd4(int na, const float *x,
+ nbnxn_bb_t *bb_work_aligned,
+ real *bb)
{
- calc_bounding_box_sse(na,x,bb_work);
-
- bb[0*STRIDE_PBB] = bb_work[BBL_X];
- bb[1*STRIDE_PBB] = bb_work[BBL_Y];
- bb[2*STRIDE_PBB] = bb_work[BBL_Z];
- bb[3*STRIDE_PBB] = bb_work[BBU_X];
- bb[4*STRIDE_PBB] = bb_work[BBU_Y];
- bb[5*STRIDE_PBB] = bb_work[BBU_Z];
+ calc_bounding_box_simd4(na, x, bb_work_aligned);
+
+ bb[0*STRIDE_PBB] = bb_work_aligned->lower[BB_X];
+ bb[1*STRIDE_PBB] = bb_work_aligned->lower[BB_Y];
+ bb[2*STRIDE_PBB] = bb_work_aligned->lower[BB_Z];
+ bb[3*STRIDE_PBB] = bb_work_aligned->upper[BB_X];
+ bb[4*STRIDE_PBB] = bb_work_aligned->upper[BB_Y];
+ bb[5*STRIDE_PBB] = bb_work_aligned->upper[BB_Z];
}
-#endif /* NBNXN_SEARCH_SSE_SINGLE */
+#endif /* NBNXN_SEARCH_SIMD4_FLOAT_X_BB */
-#ifdef NBNXN_SEARCH_BB_SSE
/* Combines pairs of consecutive bounding boxes */
-static void combine_bounding_box_pairs(nbnxn_grid_t *grid,const float *bb)
+static void combine_bounding_box_pairs(nbnxn_grid_t *grid, const nbnxn_bb_t *bb)
{
- int i,j,sc2,nc2,c2;
- __m128 min_SSE,max_SSE;
+ int i, j, sc2, nc2, c2;
- for(i=0; i<grid->ncx*grid->ncy; i++)
+ for (i = 0; i < grid->ncx*grid->ncy; i++)
{
/* Starting bb in a column is expected to be 2-aligned */
sc2 = grid->cxy_ind[i]>>1;
/* For odd numbers skip the last bb here */
nc2 = (grid->cxy_na[i]+3)>>(2+1);
- for(c2=sc2; c2<sc2+nc2; c2++)
+ for (c2 = sc2; c2 < sc2+nc2; c2++)
{
- min_SSE = _mm_min_ps(_mm_load_ps(bb+(c2*4+0)*NNBSBB_C),
- _mm_load_ps(bb+(c2*4+2)*NNBSBB_C));
- max_SSE = _mm_max_ps(_mm_load_ps(bb+(c2*4+1)*NNBSBB_C),
- _mm_load_ps(bb+(c2*4+3)*NNBSBB_C));
- _mm_store_ps(grid->bbj+(c2*2+0)*NNBSBB_C,min_SSE);
- _mm_store_ps(grid->bbj+(c2*2+1)*NNBSBB_C,max_SSE);
+#ifdef NBNXN_SEARCH_BB_SIMD4
+ gmx_simd4_float_t min_S, max_S;
+
+ min_S = gmx_simd4_min_f(gmx_simd4_load_f(&bb[c2*2+0].lower[0]),
+ gmx_simd4_load_f(&bb[c2*2+1].lower[0]));
+ max_S = gmx_simd4_max_f(gmx_simd4_load_f(&bb[c2*2+0].upper[0]),
+ gmx_simd4_load_f(&bb[c2*2+1].upper[0]));
+ gmx_simd4_store_f(&grid->bbj[c2].lower[0], min_S);
+ gmx_simd4_store_f(&grid->bbj[c2].upper[0], max_S);
+#else
+ for (j = 0; j < NNBSBB_C; j++)
+ {
+ grid->bbj[c2].lower[j] = min(bb[c2*2+0].lower[j],
+ bb[c2*2+1].lower[j]);
+ grid->bbj[c2].upper[j] = max(bb[c2*2+0].upper[j],
+ bb[c2*2+1].upper[j]);
+ }
+#endif
}
if (((grid->cxy_na[i]+3)>>2) & 1)
{
- /* Copy the last bb for odd bb count in this column */
- for(j=0; j<NNBSBB_C; j++)
+ /* The bb count in this column is odd: duplicate the last bb */
+ for (j = 0; j < NNBSBB_C; j++)
{
- grid->bbj[(c2*2+0)*NNBSBB_C+j] = bb[(c2*4+0)*NNBSBB_C+j];
- grid->bbj[(c2*2+1)*NNBSBB_C+j] = bb[(c2*4+1)*NNBSBB_C+j];
+ grid->bbj[c2].lower[j] = bb[c2*2].lower[j];
+ grid->bbj[c2].upper[j] = bb[c2*2].upper[j];
}
}
}
}
-#endif
-
/* Prints the average bb size, used for debug output */
-static void print_bbsizes_simple(FILE *fp,
+static void print_bbsizes_simple(FILE *fp,
const nbnxn_search_t nbs,
- const nbnxn_grid_t *grid)
+ const nbnxn_grid_t *grid)
{
- int c,d;
+ int c, d;
dvec ba;
clear_dvec(ba);
- for(c=0; c<grid->nc; c++)
+ for (c = 0; c < grid->nc; c++)
{
- for(d=0; d<DIM; d++)
+ for (d = 0; d < DIM; d++)
{
- ba[d] += grid->bb[c*NNBSBB_B+NNBSBB_C+d] - grid->bb[c*NNBSBB_B+d];
+ ba[d] += grid->bb[c].upper[d] - grid->bb[c].lower[d];
}
}
- dsvmul(1.0/grid->nc,ba,ba);
+ dsvmul(1.0/grid->nc, ba, ba);
- fprintf(fp,"ns bb: %4.2f %4.2f %4.2f %4.2f %4.2f %4.2f rel %4.2f %4.2f %4.2f\n",
+ fprintf(fp, "ns bb: %4.2f %4.2f %4.2f %4.2f %4.2f %4.2f rel %4.2f %4.2f %4.2f\n",
nbs->box[XX][XX]/grid->ncx,
nbs->box[YY][YY]/grid->ncy,
nbs->box[ZZ][ZZ]*grid->ncx*grid->ncy/grid->nc,
- ba[XX],ba[YY],ba[ZZ],
+ ba[XX], ba[YY], ba[ZZ],
ba[XX]*grid->ncx/nbs->box[XX][XX],
ba[YY]*grid->ncy/nbs->box[YY][YY],
ba[ZZ]*grid->nc/(grid->ncx*grid->ncy*nbs->box[ZZ][ZZ]));
}
/* Prints the average bb size, used for debug output */
-static void print_bbsizes_supersub(FILE *fp,
+static void print_bbsizes_supersub(FILE *fp,
const nbnxn_search_t nbs,
- const nbnxn_grid_t *grid)
+ const nbnxn_grid_t *grid)
{
- int ns,c,s;
+ int ns, c, s;
dvec ba;
clear_dvec(ba);
ns = 0;
- for(c=0; c<grid->nc; c++)
+ for (c = 0; c < grid->nc; c++)
{
#ifdef NBNXN_BBXXXX
- for(s=0; s<grid->nsubc[c]; s+=STRIDE_PBB)
+ for (s = 0; s < grid->nsubc[c]; s += STRIDE_PBB)
{
- int cs_w,i,d;
+ int cs_w, i, d;
cs_w = (c*GPU_NSUBCELL + s)/STRIDE_PBB;
- for(i=0; i<STRIDE_PBB; i++)
+ for (i = 0; i < STRIDE_PBB; i++)
{
- for(d=0; d<DIM; d++)
+ for (d = 0; d < DIM; d++)
{
ba[d] +=
- grid->bb[cs_w*NNBSBB_XXXX+(DIM+d)*STRIDE_PBB+i] -
- grid->bb[cs_w*NNBSBB_XXXX+ d *STRIDE_PBB+i];
+ grid->pbb[cs_w*NNBSBB_XXXX+(DIM+d)*STRIDE_PBB+i] -
+ grid->pbb[cs_w*NNBSBB_XXXX+ d *STRIDE_PBB+i];
}
}
}
#else
- for(s=0; s<grid->nsubc[c]; s++)
+ for (s = 0; s < grid->nsubc[c]; s++)
{
- int cs,d;
+ int cs, d;
cs = c*GPU_NSUBCELL + s;
- for(d=0; d<DIM; d++)
+ for (d = 0; d < DIM; d++)
{
- ba[d] +=
- grid->bb[cs*NNBSBB_B+NNBSBB_C+d] -
- grid->bb[cs*NNBSBB_B +d];
+ ba[d] += grid->bb[cs].upper[d] - grid->bb[cs].lower[d];
}
}
#endif
ns += grid->nsubc[c];
}
- dsvmul(1.0/ns,ba,ba);
+ dsvmul(1.0/ns, ba, ba);
- fprintf(fp,"ns bb: %4.2f %4.2f %4.2f %4.2f %4.2f %4.2f rel %4.2f %4.2f %4.2f\n",
+ fprintf(fp, "ns bb: %4.2f %4.2f %4.2f %4.2f %4.2f %4.2f rel %4.2f %4.2f %4.2f\n",
nbs->box[XX][XX]/(grid->ncx*GPU_NSUBCELL_X),
nbs->box[YY][YY]/(grid->ncy*GPU_NSUBCELL_Y),
nbs->box[ZZ][ZZ]*grid->ncx*grid->ncy/(grid->nc*GPU_NSUBCELL_Z),
- ba[XX],ba[YY],ba[ZZ],
+ ba[XX], ba[YY], ba[ZZ],
ba[XX]*grid->ncx*GPU_NSUBCELL_X/nbs->box[XX][XX],
ba[YY]*grid->ncy*GPU_NSUBCELL_Y/nbs->box[YY][YY],
ba[ZZ]*grid->nc*GPU_NSUBCELL_Z/(grid->ncx*grid->ncy*nbs->box[ZZ][ZZ]));
/* Potentially sorts atoms on LJ coefficients !=0 and ==0.
* Also sets interaction flags.
*/
-void sort_on_lj(nbnxn_atomdata_t *nbat,int na_c,
- int a0,int a1,const int *atinfo,
+void sort_on_lj(int na_c,
+ int a0, int a1, const int *atinfo,
int *order,
int *flags)
{
- int subc,s,a,n1,n2,a_lj_max,i,j;
- int sort1[NBNXN_NA_SC_MAX/GPU_NSUBCELL];
- int sort2[NBNXN_NA_SC_MAX/GPU_NSUBCELL];
- gmx_bool haveQ;
+ int subc, s, a, n1, n2, a_lj_max, i, j;
+ int sort1[NBNXN_NA_SC_MAX/GPU_NSUBCELL];
+ int sort2[NBNXN_NA_SC_MAX/GPU_NSUBCELL];
+ gmx_bool haveQ, bFEP;
*flags = 0;
subc = 0;
- for(s=a0; s<a1; s+=na_c)
+ for (s = a0; s < a1; s += na_c)
{
/* Make lists for this (sub-)cell on atoms with and without LJ */
- n1 = 0;
- n2 = 0;
- haveQ = FALSE;
+ n1 = 0;
+ n2 = 0;
+ haveQ = FALSE;
a_lj_max = -1;
- for(a=s; a<min(s+na_c,a1); a++)
+ for (a = s; a < min(s+na_c, a1); a++)
{
haveQ = haveQ || GET_CGINFO_HAS_Q(atinfo[order[a]]);
if (GET_CGINFO_HAS_VDW(atinfo[order[a]]))
{
sort1[n1++] = order[a];
- a_lj_max = a;
+ a_lj_max = a;
}
else
{
}
}
- /* If we don't have atom with LJ, there's nothing to sort */
+ /* If we don't have atoms with LJ, there's nothing to sort */
if (n1 > 0)
{
*flags |= NBNXN_CI_DO_LJ(subc);
*/
if (2*(a_lj_max - s) >= na_c)
{
- for(i=0; i<n1; i++)
+ for (i = 0; i < n1; i++)
{
order[a0+i] = sort1[i];
}
- for(j=0; j<n2; j++)
+ for (j = 0; j < n2; j++)
{
order[a0+n1+j] = sort2[j];
}
void fill_cell(const nbnxn_search_t nbs,
nbnxn_grid_t *grid,
nbnxn_atomdata_t *nbat,
- int a0,int a1,
+ int a0, int a1,
const int *atinfo,
rvec *x,
- int sx,int sy, int sz,
- float *bb_work)
+ int sx, int sy, int sz,
+ nbnxn_bb_t gmx_unused *bb_work_aligned)
{
- int na,a;
- size_t offset;
- float *bb_ptr;
+ int na, a;
+ size_t offset;
+ nbnxn_bb_t *bb_ptr;
+#ifdef NBNXN_BBXXXX
+ float *pbb_ptr;
+#endif
na = a1 - a0;
if (grid->bSimple)
{
- sort_on_lj(nbat,grid->na_c,a0,a1,atinfo,nbs->a,
+ sort_on_lj(grid->na_c, a0, a1, atinfo, nbs->a,
grid->flags+(a0>>grid->na_c_2log)-grid->cell0);
}
+ if (nbs->bFEP)
+ {
+ /* Set the fep flag for perturbed atoms in this (sub-)cell */
+ int c, at;
+
+ /* The grid-local cluster/(sub-)cell index */
+ c = (a0 >> grid->na_c_2log) - grid->cell0*(grid->bSimple ? 1 : GPU_NSUBCELL);
+ grid->fep[c] = 0;
+ for (at = a0; at < a1; at++)
+ {
+ if (nbs->a[at] >= 0 && GET_CGINFO_FEP(atinfo[nbs->a[at]]))
+ {
+ grid->fep[c] |= (1 << (at - a0));
+ }
+ }
+ }
+
/* Now we have sorted the atoms, set the cell indices */
- for(a=a0; a<a1; a++)
+ for (a = a0; a < a1; a++)
{
nbs->cell[nbs->a[a]] = a;
}
- copy_rvec_to_nbat_real(nbs->a+a0,a1-a0,grid->na_c,x,
- nbat->XFormat,nbat->x,a0,
- sx,sy,sz);
+ copy_rvec_to_nbat_real(nbs->a+a0, a1-a0, grid->na_c, x,
+ nbat->XFormat, nbat->x, a0,
+ sx, sy, sz);
if (nbat->XFormat == nbatX4)
{
/* Store the bounding boxes as xyz.xyz. */
- offset = ((a0 - grid->cell0*grid->na_sc)>>grid->na_c_2log)*NNBSBB_B;
+ offset = (a0 - grid->cell0*grid->na_sc) >> grid->na_c_2log;
bb_ptr = grid->bb + offset;
-#if defined GMX_DOUBLE && defined NBNXN_SEARCH_BB_SSE
+#if defined GMX_NBNXN_SIMD && GMX_SIMD_REAL_WIDTH == 2
if (2*grid->na_cj == grid->na_c)
{
- calc_bounding_box_x_x4_halves(na,nbat->x+X4_IND_A(a0),bb_ptr,
+ calc_bounding_box_x_x4_halves(na, nbat->x+X4_IND_A(a0), bb_ptr,
grid->bbj+offset*2);
}
else
#endif
{
- calc_bounding_box_x_x4(na,nbat->x+X4_IND_A(a0),bb_ptr);
+ calc_bounding_box_x_x4(na, nbat->x+X4_IND_A(a0), bb_ptr);
}
}
else if (nbat->XFormat == nbatX8)
{
/* Store the bounding boxes as xyz.xyz. */
- offset = ((a0 - grid->cell0*grid->na_sc)>>grid->na_c_2log)*NNBSBB_B;
+ offset = (a0 - grid->cell0*grid->na_sc) >> grid->na_c_2log;
bb_ptr = grid->bb + offset;
- calc_bounding_box_x_x8(na,nbat->x+X8_IND_A(a0),bb_ptr);
+ calc_bounding_box_x_x8(na, nbat->x+X8_IND_A(a0), bb_ptr);
}
#ifdef NBNXN_BBXXXX
else if (!grid->bSimple)
{
/* Store the bounding boxes in a format convenient
- * for SSE calculations: xxxxyyyyzzzz...
- */
- bb_ptr =
- grid->bb +
+ * for SIMD4 calculations: xxxxyyyyzzzz...
+ */
+ pbb_ptr =
+ grid->pbb +
((a0-grid->cell0*grid->na_sc)>>(grid->na_c_2log+STRIDE_PBB_2LOG))*NNBSBB_XXXX +
(((a0-grid->cell0*grid->na_sc)>>grid->na_c_2log) & (STRIDE_PBB-1));
-#ifdef NBNXN_SEARCH_SSE_SINGLE
+#ifdef NBNXN_SEARCH_SIMD4_FLOAT_X_BB
if (nbat->XFormat == nbatXYZQ)
{
- calc_bounding_box_xxxx_sse(na,nbat->x+a0*nbat->xstride,
- bb_work,bb_ptr);
+ calc_bounding_box_xxxx_simd4(na, nbat->x+a0*nbat->xstride,
+ bb_work_aligned, pbb_ptr);
}
else
#endif
{
- calc_bounding_box_xxxx(na,nbat->xstride,nbat->x+a0*nbat->xstride,
- bb_ptr);
+ calc_bounding_box_xxxx(na, nbat->xstride, nbat->x+a0*nbat->xstride,
+ pbb_ptr);
}
if (gmx_debug_at)
{
- fprintf(debug,"%2d %2d %2d bb %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f\n",
- sx,sy,sz,
- bb_ptr[0*STRIDE_PBB],bb_ptr[3*STRIDE_PBB],
- bb_ptr[1*STRIDE_PBB],bb_ptr[4*STRIDE_PBB],
- bb_ptr[2*STRIDE_PBB],bb_ptr[5*STRIDE_PBB]);
+ fprintf(debug, "%2d %2d %2d bb %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f\n",
+ sx, sy, sz,
+ pbb_ptr[0*STRIDE_PBB], pbb_ptr[3*STRIDE_PBB],
+ pbb_ptr[1*STRIDE_PBB], pbb_ptr[4*STRIDE_PBB],
+ pbb_ptr[2*STRIDE_PBB], pbb_ptr[5*STRIDE_PBB]);
}
}
#endif
else
{
/* Store the bounding boxes as xyz.xyz. */
- bb_ptr = grid->bb+((a0-grid->cell0*grid->na_sc)>>grid->na_c_2log)*NNBSBB_B;
+ bb_ptr = grid->bb+((a0-grid->cell0*grid->na_sc)>>grid->na_c_2log);
- calc_bounding_box(na,nbat->xstride,nbat->x+a0*nbat->xstride,
+ calc_bounding_box(na, nbat->xstride, nbat->x+a0*nbat->xstride,
bb_ptr);
if (gmx_debug_at)
{
int bbo;
bbo = (a0 - grid->cell0*grid->na_sc)/grid->na_c;
- fprintf(debug,"%2d %2d %2d bb %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f\n",
- sx,sy,sz,
- (grid->bb+bbo*NNBSBB_B)[BBL_X],
- (grid->bb+bbo*NNBSBB_B)[BBU_X],
- (grid->bb+bbo*NNBSBB_B)[BBL_Y],
- (grid->bb+bbo*NNBSBB_B)[BBU_Y],
- (grid->bb+bbo*NNBSBB_B)[BBL_Z],
- (grid->bb+bbo*NNBSBB_B)[BBU_Z]);
+ fprintf(debug, "%2d %2d %2d bb %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f\n",
+ sx, sy, sz,
+ grid->bb[bbo].lower[BB_X],
+ grid->bb[bbo].lower[BB_Y],
+ grid->bb[bbo].lower[BB_Z],
+ grid->bb[bbo].upper[BB_X],
+ grid->bb[bbo].upper[BB_Y],
+ grid->bb[bbo].upper[BB_Z]);
}
}
}
static void sort_columns_simple(const nbnxn_search_t nbs,
int dd_zone,
nbnxn_grid_t *grid,
- int a0,int a1,
+ int a0, int a1,
const int *atinfo,
rvec *x,
nbnxn_atomdata_t *nbat,
- int cxy_start,int cxy_end,
+ int cxy_start, int cxy_end,
int *sort_work)
{
int cxy;
- int cx,cy,cz,ncz,cfilled,c;
- int na,ash,ind,a;
- int na_c,ash_c;
+ int cx, cy, cz, ncz, cfilled, c;
+ int na, ash, ind, a;
+ int na_c, ash_c;
if (debug)
{
- fprintf(debug,"cell0 %d sorting columns %d - %d, atoms %d - %d\n",
- grid->cell0,cxy_start,cxy_end,a0,a1);
+ fprintf(debug, "cell0 %d sorting columns %d - %d, atoms %d - %d\n",
+ grid->cell0, cxy_start, cxy_end, a0, a1);
}
/* Sort the atoms within each x,y column in 3 dimensions */
- for(cxy=cxy_start; cxy<cxy_end; cxy++)
+ for (cxy = cxy_start; cxy < cxy_end; cxy++)
{
cx = cxy/grid->ncy;
cy = cxy - cx*grid->ncy;
ash = (grid->cell0 + grid->cxy_ind[cxy])*grid->na_sc;
/* Sort the atoms within each x,y column on z coordinate */
- sort_atoms(ZZ,FALSE,
- nbs->a+ash,na,x,
+ sort_atoms(ZZ, FALSE, dd_zone,
+ nbs->a+ash, na, x,
grid->c0[ZZ],
- ncz*grid->na_sc*SORT_GRID_OVERSIZE/nbs->box[ZZ][ZZ],
- ncz*grid->na_sc*SGSF,sort_work);
+ 1.0/nbs->box[ZZ][ZZ], ncz*grid->na_sc,
+ sort_work);
/* Fill the ncz cells in this column */
cfilled = grid->cxy_ind[cxy];
- for(cz=0; cz<ncz; cz++)
+ for (cz = 0; cz < ncz; cz++)
{
- c = grid->cxy_ind[cxy] + cz ;
+ c = grid->cxy_ind[cxy] + cz;
ash_c = ash + cz*grid->na_sc;
- na_c = min(grid->na_sc,na-(ash_c-ash));
+ na_c = min(grid->na_sc, na-(ash_c-ash));
- fill_cell(nbs,grid,nbat,
- ash_c,ash_c+na_c,atinfo,x,
+ fill_cell(nbs, grid, nbat,
+ ash_c, ash_c+na_c, atinfo, x,
grid->na_sc*cx + (dd_zone >> 2),
grid->na_sc*cy + (dd_zone & 3),
grid->na_sc*cz,
{
cfilled = c;
}
- grid->bbcz[c*NNBSBB_D ] = grid->bb[cfilled*NNBSBB_B+2];
- grid->bbcz[c*NNBSBB_D+1] = grid->bb[cfilled*NNBSBB_B+6];
+ grid->bbcz[c*NNBSBB_D ] = grid->bb[cfilled].lower[BB_Z];
+ grid->bbcz[c*NNBSBB_D+1] = grid->bb[cfilled].upper[BB_Z];
}
/* Set the unused atom indices to -1 */
- for(ind=na; ind<ncz*grid->na_sc; ind++)
+ for (ind = na; ind < ncz*grid->na_sc; ind++)
{
nbs->a[ash+ind] = -1;
}
static void sort_columns_supersub(const nbnxn_search_t nbs,
int dd_zone,
nbnxn_grid_t *grid,
- int a0,int a1,
+ int a0, int a1,
const int *atinfo,
rvec *x,
nbnxn_atomdata_t *nbat,
- int cxy_start,int cxy_end,
+ int cxy_start, int cxy_end,
int *sort_work)
{
- int cxy;
- int cx,cy,cz=-1,c=-1,ncz;
- int na,ash,na_c,ind,a;
- int subdiv_z,sub_z,na_z,ash_z;
- int subdiv_y,sub_y,na_y,ash_y;
- int subdiv_x,sub_x,na_x,ash_x;
+ int cxy;
+ int cx, cy, cz = -1, c = -1, ncz;
+ int na, ash, na_c, ind, a;
+ int subdiv_z, sub_z, na_z, ash_z;
+ int subdiv_y, sub_y, na_y, ash_y;
+ int subdiv_x, sub_x, na_x, ash_x;
- /* cppcheck-suppress unassignedVariable */
- float bb_work_array[NNBSBB_B+3],*bb_work_align;
+ nbnxn_bb_t bb_work_array[2], *bb_work_aligned;
- bb_work_align = (float *)(((size_t)(bb_work_array+3)) & (~((size_t)15)));
+ bb_work_aligned = (nbnxn_bb_t *)(((size_t)(bb_work_array+1)) & (~((size_t)15)));
if (debug)
{
- fprintf(debug,"cell0 %d sorting columns %d - %d, atoms %d - %d\n",
- grid->cell0,cxy_start,cxy_end,a0,a1);
+ fprintf(debug, "cell0 %d sorting columns %d - %d, atoms %d - %d\n",
+ grid->cell0, cxy_start, cxy_end, a0, a1);
}
subdiv_x = grid->na_c;
subdiv_z = GPU_NSUBCELL_Y*subdiv_y;
/* Sort the atoms within each x,y column in 3 dimensions */
- for(cxy=cxy_start; cxy<cxy_end; cxy++)
+ for (cxy = cxy_start; cxy < cxy_end; cxy++)
{
cx = cxy/grid->ncy;
cy = cxy - cx*grid->ncy;
ash = (grid->cell0 + grid->cxy_ind[cxy])*grid->na_sc;
/* Sort the atoms within each x,y column on z coordinate */
- sort_atoms(ZZ,FALSE,
- nbs->a+ash,na,x,
+ sort_atoms(ZZ, FALSE, dd_zone,
+ nbs->a+ash, na, x,
grid->c0[ZZ],
- ncz*grid->na_sc*SORT_GRID_OVERSIZE/nbs->box[ZZ][ZZ],
- ncz*grid->na_sc*SGSF,sort_work);
+ 1.0/nbs->box[ZZ][ZZ], ncz*grid->na_sc,
+ sort_work);
/* This loop goes over the supercells and subcells along z at once */
- for(sub_z=0; sub_z<ncz*GPU_NSUBCELL_Z; sub_z++)
+ for (sub_z = 0; sub_z < ncz*GPU_NSUBCELL_Z; sub_z++)
{
ash_z = ash + sub_z*subdiv_z;
- na_z = min(subdiv_z,na-(ash_z-ash));
+ na_z = min(subdiv_z, na-(ash_z-ash));
/* We have already sorted on z */
if (sub_z % GPU_NSUBCELL_Z == 0)
{
cz = sub_z/GPU_NSUBCELL_Z;
- c = grid->cxy_ind[cxy] + cz ;
+ c = grid->cxy_ind[cxy] + cz;
/* The number of atoms in this supercell */
- na_c = min(grid->na_sc,na-(ash_z-ash));
+ na_c = min(grid->na_sc, na-(ash_z-ash));
- grid->nsubc[c] = min(GPU_NSUBCELL,(na_c+grid->na_c-1)/grid->na_c);
+ grid->nsubc[c] = min(GPU_NSUBCELL, (na_c+grid->na_c-1)/grid->na_c);
/* Store the z-boundaries of the super cell */
grid->bbcz[c*NNBSBB_D ] = x[nbs->a[ash_z]][ZZ];
#if GPU_NSUBCELL_Y > 1
/* Sort the atoms along y */
- sort_atoms(YY,(sub_z & 1),
- nbs->a+ash_z,na_z,x,
- grid->c0[YY]+cy*grid->sy,grid->inv_sy,
- subdiv_y*SGSF,sort_work);
+ sort_atoms(YY, (sub_z & 1), dd_zone,
+ nbs->a+ash_z, na_z, x,
+ grid->c0[YY]+cy*grid->sy,
+ grid->inv_sy, subdiv_z,
+ sort_work);
#endif
- for(sub_y=0; sub_y<GPU_NSUBCELL_Y; sub_y++)
+ for (sub_y = 0; sub_y < GPU_NSUBCELL_Y; sub_y++)
{
ash_y = ash_z + sub_y*subdiv_y;
- na_y = min(subdiv_y,na-(ash_y-ash));
+ na_y = min(subdiv_y, na-(ash_y-ash));
#if GPU_NSUBCELL_X > 1
/* Sort the atoms along x */
- sort_atoms(XX,((cz*GPU_NSUBCELL_Y + sub_y) & 1),
- nbs->a+ash_y,na_y,x,
- grid->c0[XX]+cx*grid->sx,grid->inv_sx,
- subdiv_x*SGSF,sort_work);
+ sort_atoms(XX, ((cz*GPU_NSUBCELL_Y + sub_y) & 1), dd_zone,
+ nbs->a+ash_y, na_y, x,
+ grid->c0[XX]+cx*grid->sx,
+ grid->inv_sx, subdiv_y,
+ sort_work);
#endif
- for(sub_x=0; sub_x<GPU_NSUBCELL_X; sub_x++)
+ for (sub_x = 0; sub_x < GPU_NSUBCELL_X; sub_x++)
{
ash_x = ash_y + sub_x*subdiv_x;
- na_x = min(subdiv_x,na-(ash_x-ash));
+ na_x = min(subdiv_x, na-(ash_x-ash));
- fill_cell(nbs,grid,nbat,
- ash_x,ash_x+na_x,atinfo,x,
+ fill_cell(nbs, grid, nbat,
+ ash_x, ash_x+na_x, atinfo, x,
grid->na_c*(cx*GPU_NSUBCELL_X+sub_x) + (dd_zone >> 2),
grid->na_c*(cy*GPU_NSUBCELL_Y+sub_y) + (dd_zone & 3),
grid->na_c*sub_z,
- bb_work_align);
+ bb_work_aligned);
}
}
}
/* Set the unused atom indices to -1 */
- for(ind=na; ind<ncz*grid->na_sc; ind++)
+ for (ind = na; ind < ncz*grid->na_sc; ind++)
{
nbs->a[ash+ind] = -1;
}
/* Determine in which grid column atoms should go */
static void calc_column_indices(nbnxn_grid_t *grid,
- int a0,int a1,
+ int a0, int a1,
rvec *x,
- int dd_zone,const int *move,
- int thread,int nthread,
+ int dd_zone, const int *move,
+ int thread, int nthread,
int *cell,
int *cxy_na)
{
- int n0,n1,i;
- int cx,cy;
+ int n0, n1, i;
+ int cx, cy;
/* We add one extra cell for particles which moved during DD */
- for(i=0; i<grid->ncx*grid->ncy+1; i++)
+ for (i = 0; i < grid->ncx*grid->ncy+1; i++)
{
cxy_na[i] = 0;
}
if (dd_zone == 0)
{
/* Home zone */
- for(i=n0; i<n1; i++)
+ for (i = n0; i < n1; i++)
{
if (move == NULL || move[i] >= 0)
{
cx = (int)((x[i][XX] - grid->c0[XX])*grid->inv_sx);
cy = (int)((x[i][YY] - grid->c0[YY])*grid->inv_sy);
-#ifdef DEBUG_NBNXN_GRIDDING
- if (cx < 0 || cx >= grid->ncx ||
- cy < 0 || cy >= grid->ncy)
+#ifndef NDEBUG
+ if (cx < 0 || cx > grid->ncx ||
+ cy < 0 || cy > grid->ncy)
{
gmx_fatal(FARGS,
"grid cell cx %d cy %d out of range (max %d %d)\n"
"atom %f %f %f, grid->c0 %f %f",
- cx,cy,grid->ncx,grid->ncy,
- x[i][XX],x[i][YY],x[i][ZZ],grid->c0[XX],grid->c0[YY]);
+ cx, cy, grid->ncx, grid->ncy,
+ x[i][XX], x[i][YY], x[i][ZZ], grid->c0[XX], grid->c0[YY]);
}
#endif
/* Take care of potential rouding issues */
- cx = min(cx,grid->ncx - 1);
- cy = min(cy,grid->ncy - 1);
+ cx = min(cx, grid->ncx - 1);
+ cy = min(cy, grid->ncy - 1);
/* For the moment cell will contain only the, grid local,
* x and y indices, not z.
else
{
/* Non-home zone */
- for(i=n0; i<n1; i++)
+ for (i = n0; i < n1; i++)
{
cx = (int)((x[i][XX] - grid->c0[XX])*grid->inv_sx);
cy = (int)((x[i][YY] - grid->c0[YY])*grid->inv_sy);
* matter where these end up on the grid. For performance
* we put them in an extra row at the border.
*/
- cx = max(cx,0);
- cx = min(cx,grid->ncx - 1);
- cy = max(cy,0);
- cy = min(cy,grid->ncy - 1);
+ cx = max(cx, 0);
+ cx = min(cx, grid->ncx - 1);
+ cy = max(cy, 0);
+ cy = min(cy, grid->ncy - 1);
/* For the moment cell will contain only the, grid local,
* x and y indices, not z.
static void calc_cell_indices(const nbnxn_search_t nbs,
int dd_zone,
nbnxn_grid_t *grid,
- int a0,int a1,
+ int a0, int a1,
const int *atinfo,
rvec *x,
const int *move,
nbnxn_atomdata_t *nbat)
{
- int n0,n1,i;
- int cx,cy,cxy,ncz_max,ncz;
- int nthread,thread;
- int *cxy_na,cxy_na_i;
+ int n0, n1, i;
+ int cx, cy, cxy, ncz_max, ncz;
+ int nthread, thread;
+ int *cxy_na, cxy_na_i;
nthread = gmx_omp_nthreads_get(emntPairsearch);
#pragma omp parallel for num_threads(nthread) schedule(static)
- for(thread=0; thread<nthread; thread++)
+ for (thread = 0; thread < nthread; thread++)
{
- calc_column_indices(grid,a0,a1,x,dd_zone,move,thread,nthread,
- nbs->cell,nbs->work[thread].cxy_na);
+ calc_column_indices(grid, a0, a1, x, dd_zone, move, thread, nthread,
+ nbs->cell, nbs->work[thread].cxy_na);
}
/* Make the cell index as a function of x and y */
- ncz_max = 0;
- ncz = 0;
+ ncz_max = 0;
+ ncz = 0;
grid->cxy_ind[0] = 0;
- for(i=0; i<grid->ncx*grid->ncy+1; i++)
+ for (i = 0; i < grid->ncx*grid->ncy+1; i++)
{
/* We set ncz_max at the beginning of the loop iso at the end
* to skip i=grid->ncx*grid->ncy which are moved particles
ncz_max = ncz;
}
cxy_na_i = nbs->work[0].cxy_na[i];
- for(thread=1; thread<nthread; thread++)
+ for (thread = 1; thread < nthread; thread++)
{
cxy_na_i += nbs->work[thread].cxy_na[i];
}
if (debug)
{
- fprintf(debug,"ns na_sc %d na_c %d super-cells: %d x %d y %d z %.1f maxz %d\n",
- grid->na_sc,grid->na_c,grid->nc,
- grid->ncx,grid->ncy,grid->nc/((double)(grid->ncx*grid->ncy)),
+ fprintf(debug, "ns na_sc %d na_c %d super-cells: %d x %d y %d z %.1f maxz %d\n",
+ grid->na_sc, grid->na_c, grid->nc,
+ grid->ncx, grid->ncy, grid->nc/((double)(grid->ncx*grid->ncy)),
ncz_max);
if (gmx_debug_at)
{
i = 0;
- for(cy=0; cy<grid->ncy; cy++)
+ for (cy = 0; cy < grid->ncy; cy++)
{
- for(cx=0; cx<grid->ncx; cx++)
+ for (cx = 0; cx < grid->ncx; cx++)
{
- fprintf(debug," %2d",grid->cxy_ind[i+1]-grid->cxy_ind[i]);
+ fprintf(debug, " %2d", grid->cxy_ind[i+1]-grid->cxy_ind[i]);
i++;
}
- fprintf(debug,"\n");
+ fprintf(debug, "\n");
}
}
}
/* Make sure the work array for sorting is large enough */
if (ncz_max*grid->na_sc*SGSF > nbs->work[0].sort_work_nalloc)
{
- for(thread=0; thread<nbs->nthread_max; thread++)
+ for (thread = 0; thread < nbs->nthread_max; thread++)
{
nbs->work[thread].sort_work_nalloc =
over_alloc_large(ncz_max*grid->na_sc*SGSF);
srenew(nbs->work[thread].sort_work,
nbs->work[thread].sort_work_nalloc);
/* When not in use, all elements should be -1 */
- for(i=0; i<nbs->work[thread].sort_work_nalloc; i++)
+ for (i = 0; i < nbs->work[thread].sort_work_nalloc; i++)
{
nbs->work[thread].sort_work[i] = -1;
}
/* Now we know the dimensions we can fill the grid.
* This is the first, unsorted fill. We sort the columns after this.
*/
- for(i=a0; i<a1; i++)
+ for (i = a0; i < a1; i++)
{
/* At this point nbs->cell contains the local grid x,y indices */
cxy = nbs->cell[i];
n1 = grid->nc*grid->na_sc+grid->cxy_na[grid->ncx*grid->ncy];
if (dd_zone == 0)
{
- for(i=n0; i<n1; i++)
+ for (i = n0; i < n1; i++)
{
nbs->cell[nbs->a[i]] = i;
}
}
/* Sort the super-cell columns along z into the sub-cells. */
-#pragma omp parallel for num_threads(nbs->nthread_max) schedule(static)
- for(thread=0; thread<nbs->nthread_max; thread++)
+#pragma omp parallel for num_threads(nthread) schedule(static)
+ for (thread = 0; thread < nthread; thread++)
{
if (grid->bSimple)
{
- sort_columns_simple(nbs,dd_zone,grid,a0,a1,atinfo,x,nbat,
+ sort_columns_simple(nbs, dd_zone, grid, a0, a1, atinfo, x, nbat,
((thread+0)*grid->ncx*grid->ncy)/nthread,
((thread+1)*grid->ncx*grid->ncy)/nthread,
nbs->work[thread].sort_work);
}
else
{
- sort_columns_supersub(nbs,dd_zone,grid,a0,a1,atinfo,x,nbat,
+ sort_columns_supersub(nbs, dd_zone, grid, a0, a1, atinfo, x, nbat,
((thread+0)*grid->ncx*grid->ncy)/nthread,
((thread+1)*grid->ncx*grid->ncy)/nthread,
nbs->work[thread].sort_work);
}
}
-#ifdef NBNXN_SEARCH_BB_SSE
if (grid->bSimple && nbat->XFormat == nbatX8)
{
- combine_bounding_box_pairs(grid,grid->bb);
+ combine_bounding_box_pairs(grid, grid->bb);
}
-#endif
if (!grid->bSimple)
{
grid->nsubc_tot = 0;
- for(i=0; i<grid->nc; i++)
+ for (i = 0; i < grid->nc; i++)
{
grid->nsubc_tot += grid->nsubc[i];
}
{
if (grid->bSimple)
{
- print_bbsizes_simple(debug,nbs,grid);
+ print_bbsizes_simple(debug, nbs, grid);
}
else
{
- fprintf(debug,"ns non-zero sub-cells: %d average atoms %.2f\n",
- grid->nsubc_tot,(a1-a0)/(double)grid->nsubc_tot);
+ fprintf(debug, "ns non-zero sub-cells: %d average atoms %.2f\n",
+ grid->nsubc_tot, (a1-a0)/(double)grid->nsubc_tot);
- print_bbsizes_supersub(debug,nbs,grid);
+ print_bbsizes_supersub(debug, nbs, grid);
}
}
}
static void init_buffer_flags(nbnxn_buffer_flags_t *flags,
- int natoms)
+ int natoms)
{
int b;
if (flags->nflag > flags->flag_nalloc)
{
flags->flag_nalloc = over_alloc_large(flags->nflag);
- srenew(flags->flag,flags->flag_nalloc);
+ srenew(flags->flag, flags->flag_nalloc);
}
- for(b=0; b<flags->nflag; b++)
+ for (b = 0; b < flags->nflag; b++)
{
flags->flag[b] = 0;
}
* Note that without domain decomposition there is only one domain.
*/
void nbnxn_put_on_grid(nbnxn_search_t nbs,
- int ePBC,matrix box,
+ int ePBC, matrix box,
int dd_zone,
- rvec corner0,rvec corner1,
- int a0,int a1,
+ rvec corner0, rvec corner1,
+ int a0, int a1,
real atom_density,
const int *atinfo,
rvec *x,
- int nmoved,int *move,
+ int nmoved, int *move,
int nb_kernel_type,
nbnxn_atomdata_t *nbat)
{
nbnxn_grid_t *grid;
- int n;
- int nc_max_grid,nc_max;
+ int n;
+ int nc_max_grid, nc_max;
grid = &nbs->grid[dd_zone];
if (dd_zone == 0)
{
nbs->ePBC = ePBC;
- copy_mat(box,nbs->box);
+ copy_mat(box, nbs->box);
- if (atom_density >= 0)
+ /* Avoid zero density */
+ if (atom_density > 0)
{
grid->atom_density = atom_density;
}
else
{
- grid->atom_density = grid_atom_density(n-nmoved,corner0,corner1);
+ grid->atom_density = grid_atom_density(n-nmoved, corner0, corner1);
}
grid->cell0 = 0;
* for the local atoms (dd_zone=0).
*/
nbs->natoms_nonlocal = a1 - nmoved;
+
+ if (debug)
+ {
+ fprintf(debug, "natoms_local = %5d atom_density = %5.1f\n",
+ nbs->natoms_local, grid->atom_density);
+ }
}
else
{
- nbs->natoms_nonlocal = max(nbs->natoms_nonlocal,a1);
+ nbs->natoms_nonlocal = max(nbs->natoms_nonlocal, a1);
}
- nc_max_grid = set_grid_size_xy(nbs,grid,
- dd_zone,n-nmoved,corner0,corner1,
- nbs->grid[0].atom_density,
- nbat->XFormat);
+ /* We always use the home zone (grid[0]) for setting the cell size,
+ * since determining densities for non-local zones is difficult.
+ */
+ nc_max_grid = set_grid_size_xy(nbs, grid,
+ dd_zone, n-nmoved, corner0, corner1,
+ nbs->grid[0].atom_density);
nc_max = grid->cell0 + nc_max_grid;
if (a1 > nbs->cell_nalloc)
{
nbs->cell_nalloc = over_alloc_large(a1);
- srenew(nbs->cell,nbs->cell_nalloc);
+ srenew(nbs->cell, nbs->cell_nalloc);
}
/* To avoid conditionals we store the moved particles at the end of a,
if (nc_max*grid->na_sc + nmoved > nbs->a_nalloc)
{
nbs->a_nalloc = over_alloc_large(nc_max*grid->na_sc + nmoved);
- srenew(nbs->a,nbs->a_nalloc);
+ srenew(nbs->a, nbs->a_nalloc);
}
/* We need padding up to a multiple of the buffer flag size: simply add */
if (nc_max*grid->na_sc + NBNXN_BUFFERFLAG_SIZE > nbat->nalloc)
{
- nbnxn_atomdata_realloc(nbat,nc_max*grid->na_sc+NBNXN_BUFFERFLAG_SIZE);
+ nbnxn_atomdata_realloc(nbat, nc_max*grid->na_sc+NBNXN_BUFFERFLAG_SIZE);
}
- calc_cell_indices(nbs,dd_zone,grid,a0,a1,atinfo,x,move,nbat);
+ calc_cell_indices(nbs, dd_zone, grid, a0, a1, atinfo, x, move, nbat);
if (dd_zone == 0)
{
}
/* Calls nbnxn_put_on_grid for all non-local domains */
-void nbnxn_put_on_grid_nonlocal(nbnxn_search_t nbs,
+void nbnxn_put_on_grid_nonlocal(nbnxn_search_t nbs,
const gmx_domdec_zones_t *zones,
- const int *atinfo,
- rvec *x,
- int nb_kernel_type,
- nbnxn_atomdata_t *nbat)
+ const int *atinfo,
+ rvec *x,
+ int nb_kernel_type,
+ nbnxn_atomdata_t *nbat)
{
- int zone,d;
- rvec c0,c1;
+ int zone, d;
+ rvec c0, c1;
- for(zone=1; zone<zones->n; zone++)
+ for (zone = 1; zone < zones->n; zone++)
{
- for(d=0; d<DIM; d++)
+ for (d = 0; d < DIM; d++)
{
c0[d] = zones->size[zone].bb_x0[d];
c1[d] = zones->size[zone].bb_x1[d];
}
- nbnxn_put_on_grid(nbs,nbs->ePBC,NULL,
- zone,c0,c1,
+ nbnxn_put_on_grid(nbs, nbs->ePBC, NULL,
+ zone, c0, c1,
zones->cg_range[zone],
zones->cg_range[zone+1],
-1,
atinfo,
x,
- 0,NULL,
+ 0, NULL,
nb_kernel_type,
nbat);
}
}
/* Add simple grid type information to the local super/sub grid */
-void nbnxn_grid_add_simple(nbnxn_search_t nbs,
+void nbnxn_grid_add_simple(nbnxn_search_t nbs,
nbnxn_atomdata_t *nbat)
{
nbnxn_grid_t *grid;
- float *bbcz,*bb;
- int ncd,sc;
+ float *bbcz;
+ nbnxn_bb_t *bb;
+ int ncd, sc;
grid = &nbs->grid[0];
if (grid->nc*ncd > grid->nc_nalloc_simple)
{
grid->nc_nalloc_simple = over_alloc_large(grid->nc*ncd);
- srenew(grid->bbcz_simple,grid->nc_nalloc_simple*NNBSBB_D);
- srenew(grid->bb_simple,grid->nc_nalloc_simple*NNBSBB_B);
- srenew(grid->flags_simple,grid->nc_nalloc_simple);
+ srenew(grid->bbcz_simple, grid->nc_nalloc_simple*NNBSBB_D);
+ srenew(grid->bb_simple, grid->nc_nalloc_simple);
+ srenew(grid->flags_simple, grid->nc_nalloc_simple);
if (nbat->XFormat)
{
sfree_aligned(grid->bbj);
- snew_aligned(grid->bbj,grid->nc_nalloc_simple/2,16);
+ snew_aligned(grid->bbj, grid->nc_nalloc_simple/2, 16);
}
}
bb = grid->bb_simple;
#pragma omp parallel for num_threads(gmx_omp_nthreads_get(emntPairsearch)) schedule(static)
- for(sc=0; sc<grid->nc; sc++)
+ for (sc = 0; sc < grid->nc; sc++)
{
- int c,tx,na;
+ int c, tx, na;
- for(c=0; c<ncd; c++)
+ for (c = 0; c < ncd; c++)
{
tx = sc*ncd + c;
{
switch (nbat->XFormat)
{
- case nbatX4:
- /* PACK_X4==NBNXN_CPU_CLUSTER_I_SIZE, so this is simple */
- calc_bounding_box_x_x4(na,nbat->x+tx*STRIDE_P4,
- bb+tx*NNBSBB_B);
- break;
- case nbatX8:
- /* PACK_X8>NBNXN_CPU_CLUSTER_I_SIZE, more complicated */
- calc_bounding_box_x_x8(na,nbat->x+X8_IND_A(tx*NBNXN_CPU_CLUSTER_I_SIZE),
- bb+tx*NNBSBB_B);
- break;
- default:
- calc_bounding_box(na,nbat->xstride,
- nbat->x+tx*NBNXN_CPU_CLUSTER_I_SIZE*nbat->xstride,
- bb+tx*NNBSBB_B);
- break;
+ case nbatX4:
+ /* PACK_X4==NBNXN_CPU_CLUSTER_I_SIZE, so this is simple */
+ calc_bounding_box_x_x4(na, nbat->x+tx*STRIDE_P4,
+ bb+tx);
+ break;
+ case nbatX8:
+ /* PACK_X8>NBNXN_CPU_CLUSTER_I_SIZE, more complicated */
+ calc_bounding_box_x_x8(na, nbat->x+X8_IND_A(tx*NBNXN_CPU_CLUSTER_I_SIZE),
+ bb+tx);
+ break;
+ default:
+ calc_bounding_box(na, nbat->xstride,
+ nbat->x+tx*NBNXN_CPU_CLUSTER_I_SIZE*nbat->xstride,
+ bb+tx);
+ break;
}
- bbcz[tx*NNBSBB_D+0] = bb[tx*NNBSBB_B +ZZ];
- bbcz[tx*NNBSBB_D+1] = bb[tx*NNBSBB_B+NNBSBB_C+ZZ];
+ bbcz[tx*NNBSBB_D+0] = bb[tx].lower[BB_Z];
+ bbcz[tx*NNBSBB_D+1] = bb[tx].upper[BB_Z];
/* No interaction optimization yet here */
grid->flags_simple[tx] = NBNXN_CI_DO_LJ(0) | NBNXN_CI_DO_COUL(0);
}
}
-#ifdef NBNXN_SEARCH_BB_SSE
if (grid->bSimple && nbat->XFormat == nbatX8)
{
- combine_bounding_box_pairs(grid,grid->bb_simple);
+ combine_bounding_box_pairs(grid, grid->bb_simple);
}
-#endif
}
-void nbnxn_get_ncells(nbnxn_search_t nbs,int *ncx,int *ncy)
+void nbnxn_get_ncells(nbnxn_search_t nbs, int *ncx, int *ncy)
{
*ncx = nbs->grid[0].ncx;
*ncy = nbs->grid[0].ncy;
}
-void nbnxn_get_atomorder(nbnxn_search_t nbs,int **a,int *n)
+void nbnxn_get_atomorder(nbnxn_search_t nbs, int **a, int *n)
{
const nbnxn_grid_t *grid;
void nbnxn_set_atomorder(nbnxn_search_t nbs)
{
nbnxn_grid_t *grid;
- int ao,cx,cy,cxy,cz,j;
+ int ao, cx, cy, cxy, cz, j;
/* Set the atom order for the home cell (index 0) */
grid = &nbs->grid[0];
ao = 0;
- for(cx=0; cx<grid->ncx; cx++)
+ for (cx = 0; cx < grid->ncx; cx++)
{
- for(cy=0; cy<grid->ncy; cy++)
+ for (cy = 0; cy < grid->ncy; cy++)
{
cxy = cx*grid->ncy + cy;
j = grid->cxy_ind[cxy]*grid->na_sc;
- for(cz=0; cz<grid->cxy_na[cxy]; cz++)
+ for (cz = 0; cz < grid->cxy_na[cxy]; cz++)
{
nbs->a[j] = ao;
nbs->cell[ao] = j;
/* Determines the cell range along one dimension that
* the bounding box b0 - b1 sees.
*/
-static void get_cell_range(real b0,real b1,
- int nc,real c0,real s,real invs,
- real d2,real r2,int *cf,int *cl)
+static void get_cell_range(real b0, real b1,
+ int nc, real c0, real s, real invs,
+ real d2, real r2, int *cf, int *cl)
{
- *cf = max((int)((b0 - c0)*invs),0);
+ *cf = max((int)((b0 - c0)*invs), 0);
while (*cf > 0 && d2 + sqr((b0 - c0) - (*cf-1+1)*s) < r2)
{
(*cf)--;
}
- *cl = min((int)((b1 - c0)*invs),nc-1);
+ *cl = min((int)((b1 - c0)*invs), nc-1);
while (*cl < nc-1 && d2 + sqr((*cl+1)*s - (b1 - c0)) < r2)
{
(*cl)++;
}
/* Reference code calculating the distance^2 between two bounding boxes */
-static float box_dist2(float bx0,float bx1,float by0,
- float by1,float bz0,float bz1,
- const float *bb)
+static float box_dist2(float bx0, float bx1, float by0,
+ float by1, float bz0, float bz1,
+ const nbnxn_bb_t *bb)
{
float d2;
- float dl,dh,dm,dm0;
+ float dl, dh, dm, dm0;
d2 = 0;
- dl = bx0 - bb[BBU_X];
- dh = bb[BBL_X] - bx1;
- dm = max(dl,dh);
- dm0 = max(dm,0);
+ dl = bx0 - bb->upper[BB_X];
+ dh = bb->lower[BB_X] - bx1;
+ dm = max(dl, dh);
+ dm0 = max(dm, 0);
d2 += dm0*dm0;
- dl = by0 - bb[BBU_Y];
- dh = bb[BBL_Y] - by1;
- dm = max(dl,dh);
- dm0 = max(dm,0);
+ dl = by0 - bb->upper[BB_Y];
+ dh = bb->lower[BB_Y] - by1;
+ dm = max(dl, dh);
+ dm0 = max(dm, 0);
d2 += dm0*dm0;
- dl = bz0 - bb[BBU_Z];
- dh = bb[BBL_Z] - bz1;
- dm = max(dl,dh);
- dm0 = max(dm,0);
+ dl = bz0 - bb->upper[BB_Z];
+ dh = bb->lower[BB_Z] - bz1;
+ dm = max(dl, dh);
+ dm0 = max(dm, 0);
d2 += dm0*dm0;
return d2;
}
/* Plain C code calculating the distance^2 between two bounding boxes */
-static float subc_bb_dist2(int si,const float *bb_i_ci,
- int csj,const float *bb_j_all)
+static float subc_bb_dist2(int si, const nbnxn_bb_t *bb_i_ci,
+ int csj, const nbnxn_bb_t *bb_j_all)
{
- const float *bb_i,*bb_j;
- float d2;
- float dl,dh,dm,dm0;
+ const nbnxn_bb_t *bb_i, *bb_j;
+ float d2;
+ float dl, dh, dm, dm0;
- bb_i = bb_i_ci + si*NNBSBB_B;
- bb_j = bb_j_all + csj*NNBSBB_B;
+ bb_i = bb_i_ci + si;
+ bb_j = bb_j_all + csj;
d2 = 0;
- dl = bb_i[BBL_X] - bb_j[BBU_X];
- dh = bb_j[BBL_X] - bb_i[BBU_X];
- dm = max(dl,dh);
- dm0 = max(dm,0);
+ dl = bb_i->lower[BB_X] - bb_j->upper[BB_X];
+ dh = bb_j->lower[BB_X] - bb_i->upper[BB_X];
+ dm = max(dl, dh);
+ dm0 = max(dm, 0);
d2 += dm0*dm0;
- dl = bb_i[BBL_Y] - bb_j[BBU_Y];
- dh = bb_j[BBL_Y] - bb_i[BBU_Y];
- dm = max(dl,dh);
- dm0 = max(dm,0);
+ dl = bb_i->lower[BB_Y] - bb_j->upper[BB_Y];
+ dh = bb_j->lower[BB_Y] - bb_i->upper[BB_Y];
+ dm = max(dl, dh);
+ dm0 = max(dm, 0);
d2 += dm0*dm0;
- dl = bb_i[BBL_Z] - bb_j[BBU_Z];
- dh = bb_j[BBL_Z] - bb_i[BBU_Z];
- dm = max(dl,dh);
- dm0 = max(dm,0);
+ dl = bb_i->lower[BB_Z] - bb_j->upper[BB_Z];
+ dh = bb_j->lower[BB_Z] - bb_i->upper[BB_Z];
+ dm = max(dl, dh);
+ dm0 = max(dm, 0);
d2 += dm0*dm0;
return d2;
}
-#ifdef NBNXN_SEARCH_BB_SSE
+#ifdef NBNXN_SEARCH_BB_SIMD4
-/* SSE code for bb distance for bb format xyz0 */
-static float subc_bb_dist2_sse(int na_c,
- int si,const float *bb_i_ci,
- int csj,const float *bb_j_all)
+/* 4-wide SIMD code for bb distance for bb format xyz0 */
+static float subc_bb_dist2_simd4(int si, const nbnxn_bb_t *bb_i_ci,
+ int csj, const nbnxn_bb_t *bb_j_all)
{
- const float *bb_i,*bb_j;
-
- __m128 bb_i_SSE0,bb_i_SSE1;
- __m128 bb_j_SSE0,bb_j_SSE1;
- __m128 dl_SSE;
- __m128 dh_SSE;
- __m128 dm_SSE;
- __m128 dm0_SSE;
- __m128 d2_SSE;
-#ifndef GMX_X86_SSE4_1
- float d2_array[7],*d2_align;
-
- d2_align = (float *)(((size_t)(d2_array+3)) & (~((size_t)15)));
-#else
- float d2;
-#endif
-
- bb_i = bb_i_ci + si*NNBSBB_B;
- bb_j = bb_j_all + csj*NNBSBB_B;
-
- bb_i_SSE0 = _mm_load_ps(bb_i);
- bb_i_SSE1 = _mm_load_ps(bb_i+NNBSBB_C);
- bb_j_SSE0 = _mm_load_ps(bb_j);
- bb_j_SSE1 = _mm_load_ps(bb_j+NNBSBB_C);
+ gmx_simd4_float_t bb_i_S0, bb_i_S1;
+ gmx_simd4_float_t bb_j_S0, bb_j_S1;
+ gmx_simd4_float_t dl_S;
+ gmx_simd4_float_t dh_S;
+ gmx_simd4_float_t dm_S;
+ gmx_simd4_float_t dm0_S;
- dl_SSE = _mm_sub_ps(bb_i_SSE0,bb_j_SSE1);
- dh_SSE = _mm_sub_ps(bb_j_SSE0,bb_i_SSE1);
+ bb_i_S0 = gmx_simd4_load_f(&bb_i_ci[si].lower[0]);
+ bb_i_S1 = gmx_simd4_load_f(&bb_i_ci[si].upper[0]);
+ bb_j_S0 = gmx_simd4_load_f(&bb_j_all[csj].lower[0]);
+ bb_j_S1 = gmx_simd4_load_f(&bb_j_all[csj].upper[0]);
- dm_SSE = _mm_max_ps(dl_SSE,dh_SSE);
- dm0_SSE = _mm_max_ps(dm_SSE,_mm_setzero_ps());
-#ifndef GMX_X86_SSE4_1
- d2_SSE = _mm_mul_ps(dm0_SSE,dm0_SSE);
+ dl_S = gmx_simd4_sub_f(bb_i_S0, bb_j_S1);
+ dh_S = gmx_simd4_sub_f(bb_j_S0, bb_i_S1);
- _mm_store_ps(d2_align,d2_SSE);
+ dm_S = gmx_simd4_max_f(dl_S, dh_S);
+ dm0_S = gmx_simd4_max_f(dm_S, gmx_simd4_setzero_f());
- return d2_align[0] + d2_align[1] + d2_align[2];
-#else
- /* SSE4.1 dot product of components 0,1,2 */
- d2_SSE = _mm_dp_ps(dm0_SSE,dm0_SSE,0x71);
-
- _mm_store_ss(&d2,d2_SSE);
-
- return d2;
-#endif
+ return gmx_simd4_dotproduct3_f(dm0_S, dm0_S);
}
/* Calculate bb bounding distances of bb_i[si,...,si+3] and store them in d2 */
-#define SUBC_BB_DIST2_SSE_XXXX_INNER(si,bb_i,d2) \
-{ \
- int shi; \
+#define SUBC_BB_DIST2_SIMD4_XXXX_INNER(si, bb_i, d2) \
+ { \
+ int shi; \
\
- __m128 dx_0,dy_0,dz_0; \
- __m128 dx_1,dy_1,dz_1; \
+ gmx_simd4_float_t dx_0, dy_0, dz_0; \
+ gmx_simd4_float_t dx_1, dy_1, dz_1; \
\
- __m128 mx,my,mz; \
- __m128 m0x,m0y,m0z; \
+ gmx_simd4_float_t mx, my, mz; \
+ gmx_simd4_float_t m0x, m0y, m0z; \
\
- __m128 d2x,d2y,d2z; \
- __m128 d2s,d2t; \
+ gmx_simd4_float_t d2x, d2y, d2z; \
+ gmx_simd4_float_t d2s, d2t; \
\
- shi = si*NNBSBB_D*DIM; \
+ shi = si*NNBSBB_D*DIM; \
\
- xi_l = _mm_load_ps(bb_i+shi+0*STRIDE_PBB); \
- yi_l = _mm_load_ps(bb_i+shi+1*STRIDE_PBB); \
- zi_l = _mm_load_ps(bb_i+shi+2*STRIDE_PBB); \
- xi_h = _mm_load_ps(bb_i+shi+3*STRIDE_PBB); \
- yi_h = _mm_load_ps(bb_i+shi+4*STRIDE_PBB); \
- zi_h = _mm_load_ps(bb_i+shi+5*STRIDE_PBB); \
+ xi_l = gmx_simd4_load_f(bb_i+shi+0*STRIDE_PBB); \
+ yi_l = gmx_simd4_load_f(bb_i+shi+1*STRIDE_PBB); \
+ zi_l = gmx_simd4_load_f(bb_i+shi+2*STRIDE_PBB); \
+ xi_h = gmx_simd4_load_f(bb_i+shi+3*STRIDE_PBB); \
+ yi_h = gmx_simd4_load_f(bb_i+shi+4*STRIDE_PBB); \
+ zi_h = gmx_simd4_load_f(bb_i+shi+5*STRIDE_PBB); \
\
- dx_0 = _mm_sub_ps(xi_l,xj_h); \
- dy_0 = _mm_sub_ps(yi_l,yj_h); \
- dz_0 = _mm_sub_ps(zi_l,zj_h); \
+ dx_0 = gmx_simd4_sub_f(xi_l, xj_h); \
+ dy_0 = gmx_simd4_sub_f(yi_l, yj_h); \
+ dz_0 = gmx_simd4_sub_f(zi_l, zj_h); \
\
- dx_1 = _mm_sub_ps(xj_l,xi_h); \
- dy_1 = _mm_sub_ps(yj_l,yi_h); \
- dz_1 = _mm_sub_ps(zj_l,zi_h); \
+ dx_1 = gmx_simd4_sub_f(xj_l, xi_h); \
+ dy_1 = gmx_simd4_sub_f(yj_l, yi_h); \
+ dz_1 = gmx_simd4_sub_f(zj_l, zi_h); \
\
- mx = _mm_max_ps(dx_0,dx_1); \
- my = _mm_max_ps(dy_0,dy_1); \
- mz = _mm_max_ps(dz_0,dz_1); \
+ mx = gmx_simd4_max_f(dx_0, dx_1); \
+ my = gmx_simd4_max_f(dy_0, dy_1); \
+ mz = gmx_simd4_max_f(dz_0, dz_1); \
\
- m0x = _mm_max_ps(mx,zero); \
- m0y = _mm_max_ps(my,zero); \
- m0z = _mm_max_ps(mz,zero); \
+ m0x = gmx_simd4_max_f(mx, zero); \
+ m0y = gmx_simd4_max_f(my, zero); \
+ m0z = gmx_simd4_max_f(mz, zero); \
\
- d2x = _mm_mul_ps(m0x,m0x); \
- d2y = _mm_mul_ps(m0y,m0y); \
- d2z = _mm_mul_ps(m0z,m0z); \
+ d2x = gmx_simd4_mul_f(m0x, m0x); \
+ d2y = gmx_simd4_mul_f(m0y, m0y); \
+ d2z = gmx_simd4_mul_f(m0z, m0z); \
\
- d2s = _mm_add_ps(d2x,d2y); \
- d2t = _mm_add_ps(d2s,d2z); \
+ d2s = gmx_simd4_add_f(d2x, d2y); \
+ d2t = gmx_simd4_add_f(d2s, d2z); \
\
- _mm_store_ps(d2+si,d2t); \
-}
+ gmx_simd4_store_f(d2+si, d2t); \
+ }
-/* SSE code for nsi bb distances for bb format xxxxyyyyzzzz */
-static void subc_bb_dist2_sse_xxxx(const float *bb_j,
- int nsi,const float *bb_i,
- float *d2)
+/* 4-wide SIMD code for nsi bb distances for bb format xxxxyyyyzzzz */
+static void subc_bb_dist2_simd4_xxxx(const float *bb_j,
+ int nsi, const float *bb_i,
+ float *d2)
{
- __m128 xj_l,yj_l,zj_l;
- __m128 xj_h,yj_h,zj_h;
- __m128 xi_l,yi_l,zi_l;
- __m128 xi_h,yi_h,zi_h;
+ gmx_simd4_float_t xj_l, yj_l, zj_l;
+ gmx_simd4_float_t xj_h, yj_h, zj_h;
+ gmx_simd4_float_t xi_l, yi_l, zi_l;
+ gmx_simd4_float_t xi_h, yi_h, zi_h;
- __m128 zero;
+ gmx_simd4_float_t zero;
- zero = _mm_setzero_ps();
+ zero = gmx_simd4_setzero_f();
- xj_l = _mm_set1_ps(bb_j[0*STRIDE_PBB]);
- yj_l = _mm_set1_ps(bb_j[1*STRIDE_PBB]);
- zj_l = _mm_set1_ps(bb_j[2*STRIDE_PBB]);
- xj_h = _mm_set1_ps(bb_j[3*STRIDE_PBB]);
- yj_h = _mm_set1_ps(bb_j[4*STRIDE_PBB]);
- zj_h = _mm_set1_ps(bb_j[5*STRIDE_PBB]);
+ xj_l = gmx_simd4_set1_f(bb_j[0*STRIDE_PBB]);
+ yj_l = gmx_simd4_set1_f(bb_j[1*STRIDE_PBB]);
+ zj_l = gmx_simd4_set1_f(bb_j[2*STRIDE_PBB]);
+ xj_h = gmx_simd4_set1_f(bb_j[3*STRIDE_PBB]);
+ yj_h = gmx_simd4_set1_f(bb_j[4*STRIDE_PBB]);
+ zj_h = gmx_simd4_set1_f(bb_j[5*STRIDE_PBB]);
/* Here we "loop" over si (0,STRIDE_PBB) from 0 to nsi with step STRIDE_PBB.
* But as we know the number of iterations is 1 or 2, we unroll manually.
*/
- SUBC_BB_DIST2_SSE_XXXX_INNER(0,bb_i,d2);
+ SUBC_BB_DIST2_SIMD4_XXXX_INNER(0, bb_i, d2);
if (STRIDE_PBB < nsi)
{
- SUBC_BB_DIST2_SSE_XXXX_INNER(STRIDE_PBB,bb_i,d2);
+ SUBC_BB_DIST2_SIMD4_XXXX_INNER(STRIDE_PBB, bb_i, d2);
}
}
-#endif /* NBNXN_SEARCH_BB_SSE */
+#endif /* NBNXN_SEARCH_BB_SIMD4 */
/* Plain C function which determines if any atom pair between two cells
* is within distance sqrt(rl2).
*/
static gmx_bool subc_in_range_x(int na_c,
- int si,const real *x_i,
- int csj,int stride,const real *x_j,
+ int si, const real *x_i,
+ int csj, int stride, const real *x_j,
real rl2)
{
- int i,j,i0,j0;
+ int i, j, i0, j0;
real d2;
- for(i=0; i<na_c; i++)
+ for (i = 0; i < na_c; i++)
{
i0 = (si*na_c + i)*DIM;
- for(j=0; j<na_c; j++)
+ for (j = 0; j < na_c; j++)
{
j0 = (csj*na_c + j)*stride;
d2 = sqr(x_i[i0 ] - x_j[j0 ]) +
- sqr(x_i[i0+1] - x_j[j0+1]) +
- sqr(x_i[i0+2] - x_j[j0+2]);
+ sqr(x_i[i0+1] - x_j[j0+1]) +
+ sqr(x_i[i0+2] - x_j[j0+2]);
if (d2 < rl2)
{
return FALSE;
}
-/* SSE function which determines if any atom pair between two cells,
+#ifdef NBNXN_SEARCH_SIMD4_FLOAT_X_BB
+
+/* 4-wide SIMD function which determines if any atom pair between two cells,
* both with 8 atoms, is within distance sqrt(rl2).
+ * Using 8-wide AVX is not faster on Intel Sandy Bridge.
*/
-static gmx_bool subc_in_range_sse8(int na_c,
- int si,const real *x_i,
- int csj,int stride,const real *x_j,
- real rl2)
+static gmx_bool subc_in_range_simd4(int na_c,
+ int si, const real *x_i,
+ int csj, int stride, const real *x_j,
+ real rl2)
{
-#ifdef NBNXN_SEARCH_SSE_SINGLE
- __m128 ix_SSE0,iy_SSE0,iz_SSE0;
- __m128 ix_SSE1,iy_SSE1,iz_SSE1;
+ gmx_simd4_real_t ix_S0, iy_S0, iz_S0;
+ gmx_simd4_real_t ix_S1, iy_S1, iz_S1;
- __m128 rc2_SSE;
+ gmx_simd4_real_t rc2_S;
- int na_c_sse;
- int j0,j1;
+ int dim_stride;
+ int j0, j1;
- rc2_SSE = _mm_set1_ps(rl2);
+ rc2_S = gmx_simd4_set1_r(rl2);
- na_c_sse = NBNXN_GPU_CLUSTER_SIZE/STRIDE_PBB;
- ix_SSE0 = _mm_load_ps(x_i+(si*na_c_sse*DIM+0)*STRIDE_PBB);
- iy_SSE0 = _mm_load_ps(x_i+(si*na_c_sse*DIM+1)*STRIDE_PBB);
- iz_SSE0 = _mm_load_ps(x_i+(si*na_c_sse*DIM+2)*STRIDE_PBB);
- ix_SSE1 = _mm_load_ps(x_i+(si*na_c_sse*DIM+3)*STRIDE_PBB);
- iy_SSE1 = _mm_load_ps(x_i+(si*na_c_sse*DIM+4)*STRIDE_PBB);
- iz_SSE1 = _mm_load_ps(x_i+(si*na_c_sse*DIM+5)*STRIDE_PBB);
+ dim_stride = NBNXN_GPU_CLUSTER_SIZE/STRIDE_PBB*DIM;
+ ix_S0 = gmx_simd4_load_r(x_i+(si*dim_stride+0)*STRIDE_PBB);
+ iy_S0 = gmx_simd4_load_r(x_i+(si*dim_stride+1)*STRIDE_PBB);
+ iz_S0 = gmx_simd4_load_r(x_i+(si*dim_stride+2)*STRIDE_PBB);
+ ix_S1 = gmx_simd4_load_r(x_i+(si*dim_stride+3)*STRIDE_PBB);
+ iy_S1 = gmx_simd4_load_r(x_i+(si*dim_stride+4)*STRIDE_PBB);
+ iz_S1 = gmx_simd4_load_r(x_i+(si*dim_stride+5)*STRIDE_PBB);
/* We loop from the outer to the inner particles to maximize
* the chance that we find a pair in range quickly and return.
j1 = j0 + na_c - 1;
while (j0 < j1)
{
- __m128 jx0_SSE,jy0_SSE,jz0_SSE;
- __m128 jx1_SSE,jy1_SSE,jz1_SSE;
+ gmx_simd4_real_t jx0_S, jy0_S, jz0_S;
+ gmx_simd4_real_t jx1_S, jy1_S, jz1_S;
- __m128 dx_SSE0,dy_SSE0,dz_SSE0;
- __m128 dx_SSE1,dy_SSE1,dz_SSE1;
- __m128 dx_SSE2,dy_SSE2,dz_SSE2;
- __m128 dx_SSE3,dy_SSE3,dz_SSE3;
+ gmx_simd4_real_t dx_S0, dy_S0, dz_S0;
+ gmx_simd4_real_t dx_S1, dy_S1, dz_S1;
+ gmx_simd4_real_t dx_S2, dy_S2, dz_S2;
+ gmx_simd4_real_t dx_S3, dy_S3, dz_S3;
- __m128 rsq_SSE0;
- __m128 rsq_SSE1;
- __m128 rsq_SSE2;
- __m128 rsq_SSE3;
+ gmx_simd4_real_t rsq_S0;
+ gmx_simd4_real_t rsq_S1;
+ gmx_simd4_real_t rsq_S2;
+ gmx_simd4_real_t rsq_S3;
- __m128 wco_SSE0;
- __m128 wco_SSE1;
- __m128 wco_SSE2;
- __m128 wco_SSE3;
- __m128 wco_any_SSE01,wco_any_SSE23,wco_any_SSE;
+ gmx_simd4_bool_t wco_S0;
+ gmx_simd4_bool_t wco_S1;
+ gmx_simd4_bool_t wco_S2;
+ gmx_simd4_bool_t wco_S3;
+ gmx_simd4_bool_t wco_any_S01, wco_any_S23, wco_any_S;
- jx0_SSE = _mm_load1_ps(x_j+j0*stride+0);
- jy0_SSE = _mm_load1_ps(x_j+j0*stride+1);
- jz0_SSE = _mm_load1_ps(x_j+j0*stride+2);
+ jx0_S = gmx_simd4_set1_r(x_j[j0*stride+0]);
+ jy0_S = gmx_simd4_set1_r(x_j[j0*stride+1]);
+ jz0_S = gmx_simd4_set1_r(x_j[j0*stride+2]);
- jx1_SSE = _mm_load1_ps(x_j+j1*stride+0);
- jy1_SSE = _mm_load1_ps(x_j+j1*stride+1);
- jz1_SSE = _mm_load1_ps(x_j+j1*stride+2);
+ jx1_S = gmx_simd4_set1_r(x_j[j1*stride+0]);
+ jy1_S = gmx_simd4_set1_r(x_j[j1*stride+1]);
+ jz1_S = gmx_simd4_set1_r(x_j[j1*stride+2]);
/* Calculate distance */
- dx_SSE0 = _mm_sub_ps(ix_SSE0,jx0_SSE);
- dy_SSE0 = _mm_sub_ps(iy_SSE0,jy0_SSE);
- dz_SSE0 = _mm_sub_ps(iz_SSE0,jz0_SSE);
- dx_SSE1 = _mm_sub_ps(ix_SSE1,jx0_SSE);
- dy_SSE1 = _mm_sub_ps(iy_SSE1,jy0_SSE);
- dz_SSE1 = _mm_sub_ps(iz_SSE1,jz0_SSE);
- dx_SSE2 = _mm_sub_ps(ix_SSE0,jx1_SSE);
- dy_SSE2 = _mm_sub_ps(iy_SSE0,jy1_SSE);
- dz_SSE2 = _mm_sub_ps(iz_SSE0,jz1_SSE);
- dx_SSE3 = _mm_sub_ps(ix_SSE1,jx1_SSE);
- dy_SSE3 = _mm_sub_ps(iy_SSE1,jy1_SSE);
- dz_SSE3 = _mm_sub_ps(iz_SSE1,jz1_SSE);
+ dx_S0 = gmx_simd4_sub_r(ix_S0, jx0_S);
+ dy_S0 = gmx_simd4_sub_r(iy_S0, jy0_S);
+ dz_S0 = gmx_simd4_sub_r(iz_S0, jz0_S);
+ dx_S1 = gmx_simd4_sub_r(ix_S1, jx0_S);
+ dy_S1 = gmx_simd4_sub_r(iy_S1, jy0_S);
+ dz_S1 = gmx_simd4_sub_r(iz_S1, jz0_S);
+ dx_S2 = gmx_simd4_sub_r(ix_S0, jx1_S);
+ dy_S2 = gmx_simd4_sub_r(iy_S0, jy1_S);
+ dz_S2 = gmx_simd4_sub_r(iz_S0, jz1_S);
+ dx_S3 = gmx_simd4_sub_r(ix_S1, jx1_S);
+ dy_S3 = gmx_simd4_sub_r(iy_S1, jy1_S);
+ dz_S3 = gmx_simd4_sub_r(iz_S1, jz1_S);
/* rsq = dx*dx+dy*dy+dz*dz */
- rsq_SSE0 = gmx_mm_calc_rsq_ps(dx_SSE0,dy_SSE0,dz_SSE0);
- rsq_SSE1 = gmx_mm_calc_rsq_ps(dx_SSE1,dy_SSE1,dz_SSE1);
- rsq_SSE2 = gmx_mm_calc_rsq_ps(dx_SSE2,dy_SSE2,dz_SSE2);
- rsq_SSE3 = gmx_mm_calc_rsq_ps(dx_SSE3,dy_SSE3,dz_SSE3);
+ rsq_S0 = gmx_simd4_calc_rsq_r(dx_S0, dy_S0, dz_S0);
+ rsq_S1 = gmx_simd4_calc_rsq_r(dx_S1, dy_S1, dz_S1);
+ rsq_S2 = gmx_simd4_calc_rsq_r(dx_S2, dy_S2, dz_S2);
+ rsq_S3 = gmx_simd4_calc_rsq_r(dx_S3, dy_S3, dz_S3);
- wco_SSE0 = _mm_cmplt_ps(rsq_SSE0,rc2_SSE);
- wco_SSE1 = _mm_cmplt_ps(rsq_SSE1,rc2_SSE);
- wco_SSE2 = _mm_cmplt_ps(rsq_SSE2,rc2_SSE);
- wco_SSE3 = _mm_cmplt_ps(rsq_SSE3,rc2_SSE);
+ wco_S0 = gmx_simd4_cmplt_r(rsq_S0, rc2_S);
+ wco_S1 = gmx_simd4_cmplt_r(rsq_S1, rc2_S);
+ wco_S2 = gmx_simd4_cmplt_r(rsq_S2, rc2_S);
+ wco_S3 = gmx_simd4_cmplt_r(rsq_S3, rc2_S);
- wco_any_SSE01 = _mm_or_ps(wco_SSE0,wco_SSE1);
- wco_any_SSE23 = _mm_or_ps(wco_SSE2,wco_SSE3);
- wco_any_SSE = _mm_or_ps(wco_any_SSE01,wco_any_SSE23);
+ wco_any_S01 = gmx_simd4_or_b(wco_S0, wco_S1);
+ wco_any_S23 = gmx_simd4_or_b(wco_S2, wco_S3);
+ wco_any_S = gmx_simd4_or_b(wco_any_S01, wco_any_S23);
- if (_mm_movemask_ps(wco_any_SSE))
+ if (gmx_simd4_anytrue_b(wco_any_S))
{
return TRUE;
}
}
return FALSE;
-#else
- /* No SSE */
- gmx_incons("SSE function called without SSE support");
-
- return TRUE;
-#endif
}
+#endif
+
/* Returns the j sub-cell for index cj_ind */
-static int nbl_cj(const nbnxn_pairlist_t *nbl,int cj_ind)
+static int nbl_cj(const nbnxn_pairlist_t *nbl, int cj_ind)
{
return nbl->cj4[cj_ind >> NBNXN_GPU_JGROUP_SIZE_2LOG].cj[cj_ind & (NBNXN_GPU_JGROUP_SIZE - 1)];
}
/* Returns the i-interaction mask of the j sub-cell for index cj_ind */
-static unsigned nbl_imask0(const nbnxn_pairlist_t *nbl,int cj_ind)
+static unsigned int nbl_imask0(const nbnxn_pairlist_t *nbl, int cj_ind)
{
return nbl->cj4[cj_ind >> NBNXN_GPU_JGROUP_SIZE_2LOG].imei[0].imask;
}
/* Ensures there is enough space for extra extra exclusion masks */
-static void check_excl_space(nbnxn_pairlist_t *nbl,int extra)
+static void check_excl_space(nbnxn_pairlist_t *nbl, int extra)
{
if (nbl->nexcl+extra > nbl->excl_nalloc)
{
nbnxn_realloc_void((void **)&nbl->excl,
nbl->nexcl*sizeof(*nbl->excl),
nbl->excl_nalloc*sizeof(*nbl->excl),
- nbl->alloc,nbl->free);
+ nbl->alloc, nbl->free);
}
}
/* Ensures there is enough space for ncell extra j-cells in the list */
static void check_subcell_list_space_simple(nbnxn_pairlist_t *nbl,
- int ncell)
+ int ncell)
{
int cj_max;
nbnxn_realloc_void((void **)&nbl->cj,
nbl->ncj*sizeof(*nbl->cj),
nbl->cj_nalloc*sizeof(*nbl->cj),
- nbl->alloc,nbl->free);
+ nbl->alloc, nbl->free);
}
}
/* Ensures there is enough space for ncell extra j-subcells in the list */
static void check_subcell_list_space_supersub(nbnxn_pairlist_t *nbl,
- int nsupercell)
+ int nsupercell)
{
- int ncj4_max,j4,j,w,t;
+ int ncj4_max, j4, j, w, t;
#define NWARP 2
#define WARP_SIZE 32
nbnxn_realloc_void((void **)&nbl->cj4,
nbl->work->cj4_init*sizeof(*nbl->cj4),
nbl->cj4_nalloc*sizeof(*nbl->cj4),
- nbl->alloc,nbl->free);
+ nbl->alloc, nbl->free);
}
if (ncj4_max > nbl->work->cj4_init)
{
- for(j4=nbl->work->cj4_init; j4<ncj4_max; j4++)
+ for (j4 = nbl->work->cj4_init; j4 < ncj4_max; j4++)
{
/* No i-subcells and no excl's in the list initially */
- for(w=0; w<NWARP; w++)
+ for (w = 0; w < NWARP; w++)
{
nbl->cj4[j4].imei[w].imask = 0U;
nbl->cj4[j4].imei[w].excl_ind = 0;
{
int t;
- for(t=0; t<WARP_SIZE; t++)
+ for (t = 0; t < WARP_SIZE; t++)
{
/* Turn all interaction bits on */
- excl->pair[t] = NBNXN_INT_MASK_ALL;
+ excl->pair[t] = NBNXN_INTERACTION_MASK_ALL;
}
}
/* Initializes a single nbnxn_pairlist_t data structure */
static void nbnxn_init_pairlist(nbnxn_pairlist_t *nbl,
- gmx_bool bSimple,
- nbnxn_alloc_t *alloc,
- nbnxn_free_t *free)
+ gmx_bool bSimple,
+ nbnxn_alloc_t *alloc,
+ nbnxn_free_t *free)
{
if (alloc == NULL)
{
nbl->excl = NULL;
nbl->excl_nalloc = 0;
nbl->nexcl = 0;
- check_excl_space(nbl,1);
+ check_excl_space(nbl, 1);
nbl->nexcl = 1;
set_no_excls(&nbl->excl[0]);
}
- snew(nbl->work,1);
+ snew(nbl->work, 1);
+ if (nbl->bSimple)
+ {
+ snew_aligned(nbl->work->bb_ci, 1, NBNXN_SEARCH_BB_MEM_ALIGN);
+ }
+ else
+ {
#ifdef NBNXN_BBXXXX
- snew_aligned(nbl->work->bb_ci,GPU_NSUBCELL/STRIDE_PBB*NNBSBB_XXXX,NBNXN_MEM_ALIGN);
+ snew_aligned(nbl->work->pbb_ci, GPU_NSUBCELL/STRIDE_PBB*NNBSBB_XXXX, NBNXN_SEARCH_BB_MEM_ALIGN);
#else
- snew_aligned(nbl->work->bb_ci,GPU_NSUBCELL*NNBSBB_B,NBNXN_MEM_ALIGN);
+ snew_aligned(nbl->work->bb_ci, GPU_NSUBCELL, NBNXN_SEARCH_BB_MEM_ALIGN);
#endif
- snew_aligned(nbl->work->x_ci,NBNXN_NA_SC_MAX*DIM,NBNXN_MEM_ALIGN);
+ }
+ snew_aligned(nbl->work->x_ci, NBNXN_NA_SC_MAX*DIM, NBNXN_SEARCH_BB_MEM_ALIGN);
#ifdef GMX_NBNXN_SIMD
- snew_aligned(nbl->work->x_ci_simd_4xn,1,NBNXN_MEM_ALIGN);
- snew_aligned(nbl->work->x_ci_simd_2xnn,1,NBNXN_MEM_ALIGN);
+ snew_aligned(nbl->work->x_ci_simd_4xn, 1, NBNXN_MEM_ALIGN);
+ snew_aligned(nbl->work->x_ci_simd_2xnn, 1, NBNXN_MEM_ALIGN);
#endif
- snew_aligned(nbl->work->d2,GPU_NSUBCELL,NBNXN_MEM_ALIGN);
+ snew_aligned(nbl->work->d2, GPU_NSUBCELL, NBNXN_SEARCH_BB_MEM_ALIGN);
+
+ nbl->work->sort = NULL;
+ nbl->work->sort_nalloc = 0;
+ nbl->work->sci_sort = NULL;
+ nbl->work->sci_sort_nalloc = 0;
}
void nbnxn_init_pairlist_set(nbnxn_pairlist_set_t *nbl_list,
if (!nbl_list->bCombined &&
nbl_list->nnbl > NBNXN_BUFFERFLAG_MAX_THREADS)
{
- gmx_fatal(FARGS,"%d OpenMP threads were requested. Since the non-bonded force buffer reduction is prohibitively slow with more than %d threads, we do not allow this. Use %d or less OpenMP threads.",
- nbl_list->nnbl,NBNXN_BUFFERFLAG_MAX_THREADS,NBNXN_BUFFERFLAG_MAX_THREADS);
+ gmx_fatal(FARGS, "%d OpenMP threads were requested. Since the non-bonded force buffer reduction is prohibitively slow with more than %d threads, we do not allow this. Use %d or less OpenMP threads.",
+ nbl_list->nnbl, NBNXN_BUFFERFLAG_MAX_THREADS, NBNXN_BUFFERFLAG_MAX_THREADS);
}
- snew(nbl_list->nbl,nbl_list->nnbl);
+ snew(nbl_list->nbl, nbl_list->nnbl);
+ snew(nbl_list->nbl_fep, nbl_list->nnbl);
/* Execute in order to avoid memory interleaving between threads */
#pragma omp parallel for num_threads(nbl_list->nnbl) schedule(static)
- for(i=0; i<nbl_list->nnbl; i++)
+ for (i = 0; i < nbl_list->nnbl; i++)
{
/* Allocate the nblist data structure locally on each thread
* to optimize memory access for NUMA architectures.
*/
- snew(nbl_list->nbl[i],1);
+ snew(nbl_list->nbl[i], 1);
/* Only list 0 is used on the GPU, use normal allocation for i>0 */
if (i == 0)
{
- nbnxn_init_pairlist(nbl_list->nbl[i],nbl_list->bSimple,alloc,free);
+ nbnxn_init_pairlist(nbl_list->nbl[i], nbl_list->bSimple, alloc, free);
}
else
{
- nbnxn_init_pairlist(nbl_list->nbl[i],nbl_list->bSimple,NULL,NULL);
+ nbnxn_init_pairlist(nbl_list->nbl[i], nbl_list->bSimple, NULL, NULL);
}
+
+ snew(nbl_list->nbl_fep[i], 1);
+ nbnxn_init_pairlist_fep(nbl_list->nbl_fep[i]);
}
}
/* Print statistics of a pair list, used for debug output */
-static void print_nblist_statistics_simple(FILE *fp,const nbnxn_pairlist_t *nbl,
- const nbnxn_search_t nbs,real rl)
+static void print_nblist_statistics_simple(FILE *fp, const nbnxn_pairlist_t *nbl,
+ const nbnxn_search_t nbs, real rl)
{
const nbnxn_grid_t *grid;
- int cs[SHIFTS];
- int s,i,j;
- int npexcl;
+ int cs[SHIFTS];
+ int s, i, j;
+ int npexcl;
/* This code only produces correct statistics with domain decomposition */
grid = &nbs->grid[0];
- fprintf(fp,"nbl nci %d ncj %d\n",
- nbl->nci,nbl->ncj);
- fprintf(fp,"nbl na_sc %d rl %g ncp %d per cell %.1f atoms %.1f ratio %.2f\n",
- nbl->na_sc,rl,nbl->ncj,nbl->ncj/(double)grid->nc,
+ fprintf(fp, "nbl nci %d ncj %d\n",
+ nbl->nci, nbl->ncj);
+ fprintf(fp, "nbl na_sc %d rl %g ncp %d per cell %.1f atoms %.1f ratio %.2f\n",
+ nbl->na_sc, rl, nbl->ncj, nbl->ncj/(double)grid->nc,
nbl->ncj/(double)grid->nc*grid->na_sc,
nbl->ncj/(double)grid->nc*grid->na_sc/(0.5*4.0/3.0*M_PI*rl*rl*rl*grid->nc*grid->na_sc/det(nbs->box)));
- fprintf(fp,"nbl average j cell list length %.1f\n",
+ fprintf(fp, "nbl average j cell list length %.1f\n",
0.25*nbl->ncj/(double)nbl->nci);
- for(s=0; s<SHIFTS; s++)
+ for (s = 0; s < SHIFTS; s++)
{
cs[s] = 0;
}
npexcl = 0;
- for(i=0; i<nbl->nci; i++)
+ for (i = 0; i < nbl->nci; i++)
{
cs[nbl->ci[i].shift & NBNXN_CI_SHIFT] +=
nbl->ci[i].cj_ind_end - nbl->ci[i].cj_ind_start;
j = nbl->ci[i].cj_ind_start;
while (j < nbl->ci[i].cj_ind_end &&
- nbl->cj[j].excl != NBNXN_INT_MASK_ALL)
+ nbl->cj[j].excl != NBNXN_INTERACTION_MASK_ALL)
{
npexcl++;
j++;
}
}
- fprintf(fp,"nbl cell pairs, total: %d excl: %d %.1f%%\n",
- nbl->ncj,npexcl,100*npexcl/(double)nbl->ncj);
- for(s=0; s<SHIFTS; s++)
+ fprintf(fp, "nbl cell pairs, total: %d excl: %d %.1f%%\n",
+ nbl->ncj, npexcl, 100*npexcl/(double)nbl->ncj);
+ for (s = 0; s < SHIFTS; s++)
{
if (cs[s] > 0)
{
- fprintf(fp,"nbl shift %2d ncj %3d\n",s,cs[s]);
+ fprintf(fp, "nbl shift %2d ncj %3d\n", s, cs[s]);
}
}
}
/* Print statistics of a pair lists, used for debug output */
-static void print_nblist_statistics_supersub(FILE *fp,const nbnxn_pairlist_t *nbl,
- const nbnxn_search_t nbs,real rl)
+static void print_nblist_statistics_supersub(FILE *fp, const nbnxn_pairlist_t *nbl,
+ const nbnxn_search_t nbs, real rl)
{
const nbnxn_grid_t *grid;
- int i,j4,j,si,b;
- int c[GPU_NSUBCELL+1];
+ int i, j4, j, si, b;
+ int c[GPU_NSUBCELL+1];
/* This code only produces correct statistics with domain decomposition */
grid = &nbs->grid[0];
- fprintf(fp,"nbl nsci %d ncj4 %d nsi %d excl4 %d\n",
- nbl->nsci,nbl->ncj4,nbl->nci_tot,nbl->nexcl);
- fprintf(fp,"nbl na_c %d rl %g ncp %d per cell %.1f atoms %.1f ratio %.2f\n",
- nbl->na_ci,rl,nbl->nci_tot,nbl->nci_tot/(double)grid->nsubc_tot,
+ fprintf(fp, "nbl nsci %d ncj4 %d nsi %d excl4 %d\n",
+ nbl->nsci, nbl->ncj4, nbl->nci_tot, nbl->nexcl);
+ fprintf(fp, "nbl na_c %d rl %g ncp %d per cell %.1f atoms %.1f ratio %.2f\n",
+ nbl->na_ci, rl, nbl->nci_tot, nbl->nci_tot/(double)grid->nsubc_tot,
nbl->nci_tot/(double)grid->nsubc_tot*grid->na_c,
nbl->nci_tot/(double)grid->nsubc_tot*grid->na_c/(0.5*4.0/3.0*M_PI*rl*rl*rl*grid->nsubc_tot*grid->na_c/det(nbs->box)));
- fprintf(fp,"nbl average j super cell list length %.1f\n",
+ fprintf(fp, "nbl average j super cell list length %.1f\n",
0.25*nbl->ncj4/(double)nbl->nsci);
- fprintf(fp,"nbl average i sub cell list length %.1f\n",
+ fprintf(fp, "nbl average i sub cell list length %.1f\n",
nbl->nci_tot/((double)nbl->ncj4));
- for(si=0; si<=GPU_NSUBCELL; si++)
+ for (si = 0; si <= GPU_NSUBCELL; si++)
{
c[si] = 0;
}
- for(i=0; i<nbl->nsci; i++)
+ for (i = 0; i < nbl->nsci; i++)
{
- for(j4=nbl->sci[i].cj4_ind_start; j4<nbl->sci[i].cj4_ind_end; j4++)
+ for (j4 = nbl->sci[i].cj4_ind_start; j4 < nbl->sci[i].cj4_ind_end; j4++)
{
- for(j=0; j<NBNXN_GPU_JGROUP_SIZE; j++)
+ for (j = 0; j < NBNXN_GPU_JGROUP_SIZE; j++)
{
b = 0;
- for(si=0; si<GPU_NSUBCELL; si++)
+ for (si = 0; si < GPU_NSUBCELL; si++)
{
if (nbl->cj4[j4].imei[0].imask & (1U << (j*GPU_NSUBCELL + si)))
{
}
}
}
- for(b=0; b<=GPU_NSUBCELL; b++)
- {
- fprintf(fp,"nbl j-list #i-subcell %d %7d %4.1f\n",
- b,c[b],100.0*c[b]/(double)(nbl->ncj4*NBNXN_GPU_JGROUP_SIZE));
- }
-}
-
-/* Print the full pair list, used for debug output */
-static void print_supersub_nsp(const char *fn,
- const nbnxn_pairlist_t *nbl,
- int iloc)
-{
- char buf[STRLEN];
- FILE *fp;
- int i,nsp,j4,p;
-
- sprintf(buf,"%s_%s.xvg",fn,NONLOCAL_I(iloc) ? "nl" : "l");
- fp = ffopen(buf,"w");
-
- for(i=0; i<nbl->nci; i++)
+ for (b = 0; b <= GPU_NSUBCELL; b++)
{
- nsp = 0;
- for(j4=nbl->sci[i].cj4_ind_start; j4<nbl->sci[i].cj4_ind_end; j4++)
- {
- for(p=0; p<NBNXN_GPU_JGROUP_SIZE*GPU_NSUBCELL; p++)
- {
- nsp += (nbl->cj4[j4].imei[0].imask >> p) & 1;
- }
- }
- fprintf(fp,"%4d %3d %3d\n",
- i,
- nsp,
- nbl->sci[i].cj4_ind_end-nbl->sci[i].cj4_ind_start);
+ fprintf(fp, "nbl j-list #i-subcell %d %7d %4.1f\n",
+ b, c[b], 100.0*c[b]/(double)(nbl->ncj4*NBNXN_GPU_JGROUP_SIZE));
}
-
- fclose(fp);
}
/* Returns a pointer to the exclusion mask for cj4-unit cj4, warp warp */
-static void low_get_nbl_exclusions(nbnxn_pairlist_t *nbl,int cj4,
- int warp,nbnxn_excl_t **excl)
+static void low_get_nbl_exclusions(nbnxn_pairlist_t *nbl, int cj4,
+ int warp, nbnxn_excl_t **excl)
{
if (nbl->cj4[cj4].imei[warp].excl_ind == 0)
{
}
/* Returns a pointer to the exclusion mask for cj4-unit cj4, warp warp,
- * allocates extra memory, if necessary.
+ * generates a new element and allocates extra memory, if necessary.
*/
-static void get_nbl_exclusions_1(nbnxn_pairlist_t *nbl,int cj4,
- int warp,nbnxn_excl_t **excl)
+static void get_nbl_exclusions_1(nbnxn_pairlist_t *nbl, int cj4,
+ int warp, nbnxn_excl_t **excl)
{
if (nbl->cj4[cj4].imei[warp].excl_ind == 0)
{
/* We need to make a new list entry, check if we have space */
- check_excl_space(nbl,1);
+ check_excl_space(nbl, 1);
}
- low_get_nbl_exclusions(nbl,cj4,warp,excl);
+ low_get_nbl_exclusions(nbl, cj4, warp, excl);
}
/* Returns pointers to the exclusion mask for cj4-unit cj4 for both warps,
- * allocates extra memory, if necessary.
+ * generates a new element and allocates extra memory, if necessary.
*/
-static void get_nbl_exclusions_2(nbnxn_pairlist_t *nbl,int cj4,
+static void get_nbl_exclusions_2(nbnxn_pairlist_t *nbl, int cj4,
nbnxn_excl_t **excl_w0,
nbnxn_excl_t **excl_w1)
{
/* Check for space we might need */
- check_excl_space(nbl,2);
+ check_excl_space(nbl, 2);
- low_get_nbl_exclusions(nbl,cj4,0,excl_w0);
- low_get_nbl_exclusions(nbl,cj4,1,excl_w1);
+ low_get_nbl_exclusions(nbl, cj4, 0, excl_w0);
+ low_get_nbl_exclusions(nbl, cj4, 1, excl_w1);
}
/* Sets the self exclusions i=j and pair exclusions i>j */
static void set_self_and_newton_excls_supersub(nbnxn_pairlist_t *nbl,
- int cj4_ind,int sj_offset,
+ int cj4_ind, int sj_offset,
int si)
{
nbnxn_excl_t *excl[2];
- int ei,ej,w;
+ int ei, ej, w;
/* Here we only set the set self and double pair exclusions */
- get_nbl_exclusions_2(nbl,cj4_ind,&excl[0],&excl[1]);
+ get_nbl_exclusions_2(nbl, cj4_ind, &excl[0], &excl[1]);
/* Only minor < major bits set */
- for(ej=0; ej<nbl->na_ci; ej++)
+ for (ej = 0; ej < nbl->na_ci; ej++)
{
w = (ej>>2);
- for(ei=ej; ei<nbl->na_ci; ei++)
+ for (ei = ej; ei < nbl->na_ci; ei++)
{
excl[w]->pair[(ej & (NBNXN_GPU_JGROUP_SIZE-1))*nbl->na_ci + ei] &=
~(1U << (sj_offset*GPU_NSUBCELL + si));
}
/* Returns a diagonal or off-diagonal interaction mask for plain C lists */
-static unsigned int get_imask(gmx_bool rdiag,int ci,int cj)
+static unsigned int get_imask(gmx_bool rdiag, int ci, int cj)
{
- return (rdiag && ci == cj ? NBNXN_INT_MASK_DIAG : NBNXN_INT_MASK_ALL);
+ return (rdiag && ci == cj ? NBNXN_INTERACTION_MASK_DIAG : NBNXN_INTERACTION_MASK_ALL);
}
-/* Returns a diagonal or off-diagonal interaction mask for SIMD128 lists */
-static unsigned int get_imask_x86_simd128(gmx_bool rdiag,int ci,int cj)
+/* Returns a diagonal or off-diagonal interaction mask for cj-size=2 */
+static unsigned int get_imask_simd_j2(gmx_bool rdiag, int ci, int cj)
{
-#ifndef GMX_DOUBLE /* cj-size = 4 */
- return (rdiag && ci == cj ? NBNXN_INT_MASK_DIAG : NBNXN_INT_MASK_ALL);
-#else /* cj-size = 2 */
- return (rdiag && ci*2 == cj ? NBNXN_INT_MASK_DIAG_J2_0 :
- (rdiag && ci*2+1 == cj ? NBNXN_INT_MASK_DIAG_J2_1 :
- NBNXN_INT_MASK_ALL));
-#endif
+ return (rdiag && ci*2 == cj ? NBNXN_INTERACTION_MASK_DIAG_J2_0 :
+ (rdiag && ci*2+1 == cj ? NBNXN_INTERACTION_MASK_DIAG_J2_1 :
+ NBNXN_INTERACTION_MASK_ALL));
}
-/* Returns a diagonal or off-diagonal interaction mask for SIMD256 lists */
-static unsigned int get_imask_x86_simd256(gmx_bool rdiag,int ci,int cj)
+/* Returns a diagonal or off-diagonal interaction mask for cj-size=4 */
+static unsigned int get_imask_simd_j4(gmx_bool rdiag, int ci, int cj)
{
-#ifndef GMX_DOUBLE /* cj-size = 8 */
- return (rdiag && ci == cj*2 ? NBNXN_INT_MASK_DIAG_J8_0 :
- (rdiag && ci == cj*2+1 ? NBNXN_INT_MASK_DIAG_J8_1 :
- NBNXN_INT_MASK_ALL));
-#else /* cj-size = 4 */
- return (rdiag && ci == cj ? NBNXN_INT_MASK_DIAG : NBNXN_INT_MASK_ALL);
-#endif
+ return (rdiag && ci == cj ? NBNXN_INTERACTION_MASK_DIAG : NBNXN_INTERACTION_MASK_ALL);
+}
+
+/* Returns a diagonal or off-diagonal interaction mask for cj-size=8 */
+static unsigned int get_imask_simd_j8(gmx_bool rdiag, int ci, int cj)
+{
+ return (rdiag && ci == cj*2 ? NBNXN_INTERACTION_MASK_DIAG_J8_0 :
+ (rdiag && ci == cj*2+1 ? NBNXN_INTERACTION_MASK_DIAG_J8_1 :
+ NBNXN_INTERACTION_MASK_ALL));
}
#ifdef GMX_NBNXN_SIMD
-#if GMX_NBNXN_SIMD_BITWIDTH == 128
-#define get_imask_x86_simd_4xn get_imask_x86_simd128
-#else
-#if GMX_NBNXN_SIMD_BITWIDTH == 256
-#define get_imask_x86_simd_4xn get_imask_x86_simd256
-#define get_imask_x86_simd_2xnn get_imask_x86_simd128
-#else
-#error "unsupported GMX_NBNXN_SIMD_BITWIDTH"
+#if GMX_SIMD_REAL_WIDTH == 2
+#define get_imask_simd_4xn get_imask_simd_j2
#endif
+#if GMX_SIMD_REAL_WIDTH == 4
+#define get_imask_simd_4xn get_imask_simd_j4
+#endif
+#if GMX_SIMD_REAL_WIDTH == 8
+#define get_imask_simd_4xn get_imask_simd_j8
+#define get_imask_simd_2xnn get_imask_simd_j4
+#endif
+#if GMX_SIMD_REAL_WIDTH == 16
+#define get_imask_simd_2xnn get_imask_simd_j8
#endif
#endif
*/
static void make_cluster_list_simple(const nbnxn_grid_t *gridj,
nbnxn_pairlist_t *nbl,
- int ci,int cjf,int cjl,
+ int ci, int cjf, int cjl,
gmx_bool remove_sub_diag,
const real *x_j,
- real rl2,float rbb2,
+ real rl2, float rbb2,
int *ndistc)
{
const nbnxn_list_work_t *work;
- const float *bb_ci;
- const real *x_ci;
+ const nbnxn_bb_t *bb_ci;
+ const real *x_ci;
- gmx_bool InRange;
- real d2;
- int cjf_gl,cjl_gl,cj;
+ gmx_bool InRange;
+ real d2;
+ int cjf_gl, cjl_gl, cj;
work = nbl->work;
InRange = FALSE;
while (!InRange && cjf <= cjl)
{
- d2 = subc_bb_dist2(0,bb_ci,cjf,gridj->bb);
+ d2 = subc_bb_dist2(0, bb_ci, cjf, gridj->bb);
*ndistc += 2;
/* Check if the distance is within the distance where
}
else if (d2 < rl2)
{
- int i,j;
+ int i, j;
cjf_gl = gridj->cell0 + cjf;
- for(i=0; i<NBNXN_CPU_CLUSTER_I_SIZE && !InRange; i++)
+ for (i = 0; i < NBNXN_CPU_CLUSTER_I_SIZE && !InRange; i++)
{
- for(j=0; j<NBNXN_CPU_CLUSTER_I_SIZE; j++)
+ for (j = 0; j < NBNXN_CPU_CLUSTER_I_SIZE; j++)
{
InRange = InRange ||
(sqr(x_ci[i*STRIDE_XYZ+XX] - x_j[(cjf_gl*NBNXN_CPU_CLUSTER_I_SIZE+j)*STRIDE_XYZ+XX]) +
InRange = FALSE;
while (!InRange && cjl > cjf)
{
- d2 = subc_bb_dist2(0,bb_ci,cjl,gridj->bb);
+ d2 = subc_bb_dist2(0, bb_ci, cjl, gridj->bb);
*ndistc += 2;
/* Check if the distance is within the distance where
}
else if (d2 < rl2)
{
- int i,j;
+ int i, j;
cjl_gl = gridj->cell0 + cjl;
- for(i=0; i<NBNXN_CPU_CLUSTER_I_SIZE && !InRange; i++)
+ for (i = 0; i < NBNXN_CPU_CLUSTER_I_SIZE && !InRange; i++)
{
- for(j=0; j<NBNXN_CPU_CLUSTER_I_SIZE; j++)
+ for (j = 0; j < NBNXN_CPU_CLUSTER_I_SIZE; j++)
{
InRange = InRange ||
(sqr(x_ci[i*STRIDE_XYZ+XX] - x_j[(cjl_gl*NBNXN_CPU_CLUSTER_I_SIZE+j)*STRIDE_XYZ+XX]) +
if (cjf <= cjl)
{
- for(cj=cjf; cj<=cjl; cj++)
+ for (cj = cjf; cj <= cjl; cj++)
{
/* Store cj and the interaction mask */
nbl->cj[nbl->ncj].cj = gridj->cell0 + cj;
- nbl->cj[nbl->ncj].excl = get_imask(remove_sub_diag,ci,cj);
+ nbl->cj[nbl->ncj].excl = get_imask(remove_sub_diag, ci, cj);
nbl->ncj++;
}
/* Increase the closing index in i super-cell list */
}
#ifdef GMX_NBNXN_SIMD_4XN
-#include "nbnxn_search_simd_4xn.h"
+#include "gromacs/mdlib/nbnxn_search_simd_4xn.h"
#endif
#ifdef GMX_NBNXN_SIMD_2XNN
-#include "nbnxn_search_simd_2xnn.h"
+#include "gromacs/mdlib/nbnxn_search_simd_2xnn.h"
#endif
-/* Plain C or SSE code for making a pair list of super-cell sci vs scj.
+/* Plain C or SIMD4 code for making a pair list of super-cell sci vs scj.
* Checks bounding box distances and possibly atom pair distances.
*/
-static void make_cluster_list_supersub(const nbnxn_search_t nbs,
- const nbnxn_grid_t *gridi,
+static void make_cluster_list_supersub(const nbnxn_grid_t *gridi,
const nbnxn_grid_t *gridj,
nbnxn_pairlist_t *nbl,
- int sci,int scj,
+ int sci, int scj,
gmx_bool sci_equals_scj,
- int stride,const real *x,
- real rl2,float rbb2,
+ int stride, const real *x,
+ real rl2, float rbb2,
int *ndistc)
{
- int na_c;
- int npair;
- int cjo,ci1,ci,cj,cj_gl;
- int cj4_ind,cj_offset;
- unsigned imask;
- nbnxn_cj4_t *cj4;
- const float *bb_ci;
- const real *x_ci;
- float *d2l,d2;
- int w;
+ int na_c;
+ int npair;
+ int cjo, ci1, ci, cj, cj_gl;
+ int cj4_ind, cj_offset;
+ unsigned int imask;
+ nbnxn_cj4_t *cj4;
+#ifdef NBNXN_BBXXXX
+ const float *pbb_ci;
+#else
+ const nbnxn_bb_t *bb_ci;
+#endif
+ const real *x_ci;
+ float *d2l, d2;
+ int w;
#define PRUNE_LIST_CPU_ONE
#ifdef PRUNE_LIST_CPU_ONE
- int ci_last=-1;
+ int ci_last = -1;
#endif
d2l = nbl->work->d2;
- bb_ci = nbl->work->bb_ci;
- x_ci = nbl->work->x_ci;
+#ifdef NBNXN_BBXXXX
+ pbb_ci = nbl->work->pbb_ci;
+#else
+ bb_ci = nbl->work->bb_ci;
+#endif
+ x_ci = nbl->work->x_ci;
na_c = gridj->na_c;
- for(cjo=0; cjo<gridj->nsubc[scj]; cjo++)
+ for (cjo = 0; cjo < gridj->nsubc[scj]; cjo++)
{
cj4_ind = (nbl->work->cj_ind >> NBNXN_GPU_JGROUP_SIZE_2LOG);
cj_offset = nbl->work->cj_ind - cj4_ind*NBNXN_GPU_JGROUP_SIZE;
}
#ifdef NBNXN_BBXXXX
- /* Determine all ci1 bb distances in one call with SSE */
- subc_bb_dist2_sse_xxxx(gridj->bb+(cj>>STRIDE_PBB_2LOG)*NNBSBB_XXXX+(cj & (STRIDE_PBB-1)),
- ci1,bb_ci,d2l);
+ /* Determine all ci1 bb distances in one call with SIMD4 */
+ subc_bb_dist2_simd4_xxxx(gridj->pbb+(cj>>STRIDE_PBB_2LOG)*NNBSBB_XXXX+(cj & (STRIDE_PBB-1)),
+ ci1, pbb_ci, d2l);
*ndistc += na_c*2;
#endif
npair = 0;
/* We use a fixed upper-bound instead of ci1 to help optimization */
- for(ci=0; ci<GPU_NSUBCELL; ci++)
+ for (ci = 0; ci < GPU_NSUBCELL; ci++)
{
if (ci == ci1)
{
#ifndef NBNXN_BBXXXX
/* Determine the bb distance between ci and cj */
- d2l[ci] = subc_bb_dist2(ci,bb_ci,cj,gridj->bb);
+ d2l[ci] = subc_bb_dist2(ci, bb_ci, cj, gridj->bb);
*ndistc += 2;
#endif
d2 = d2l[ci];
*/
*ndistc += na_c*na_c;
if (d2 < rbb2 ||
- (d2 < rl2 && subc_in_range_x(na_c,ci,x_ci,cj_gl,stride,x,rl2)))
+ (d2 < rl2 &&
+#ifdef NBNXN_PBB_SIMD4
+ subc_in_range_simd4
+#else
+ subc_in_range_x
+#endif
+ (na_c, ci, x_ci, cj_gl, stride, x, rl2)))
#else
/* Check if the distance between the two bounding boxes
* in within the pair-list cut-off.
{
/* Avoid using function pointers here, as it's slower */
if (
-#ifdef NBNXN_PBB_SSE
- !subc_in_range_sse8
+#ifdef NBNXN_PBB_SIMD4
+ !subc_in_range_simd4
#else
!subc_in_range_x
#endif
- (na_c,ci_last,x_ci,cj_gl,stride,x,rl2))
+ (na_c, ci_last, x_ci, cj_gl, stride, x, rl2))
{
imask &= ~(1U << (cj_offset*GPU_NSUBCELL+ci_last));
npair--;
*/
if (sci_equals_scj)
{
- set_self_and_newton_excls_supersub(nbl,cj4_ind,cj_offset,cjo);
+ set_self_and_newton_excls_supersub(nbl, cj4_ind, cj_offset, cjo);
}
/* Copy the cluster interaction mask to the list */
- for(w=0; w<NWARP; w++)
+ for (w = 0; w < NWARP; w++)
{
cj4->imei[w].imask |= imask;
}
* as masks in the pair-list for simple list i-entry nbl_ci
*/
static void set_ci_top_excls(const nbnxn_search_t nbs,
- nbnxn_pairlist_t *nbl,
- gmx_bool diagRemoved,
- int na_ci_2log,
- int na_cj_2log,
- const nbnxn_ci_t *nbl_ci,
- const t_blocka *excl)
+ nbnxn_pairlist_t *nbl,
+ gmx_bool diagRemoved,
+ int na_ci_2log,
+ int na_cj_2log,
+ const nbnxn_ci_t *nbl_ci,
+ const t_blocka *excl)
{
- const int *cell;
- int ci;
- int cj_ind_first,cj_ind_last;
- int cj_first,cj_last;
- int ndirect;
- int i,ai,aj,si,eind,ge,se;
- int found,cj_ind_0,cj_ind_1,cj_ind_m;
- int cj_m;
- gmx_bool Found_si;
- int si_ind;
+ const int *cell;
+ int ci;
+ int cj_ind_first, cj_ind_last;
+ int cj_first, cj_last;
+ int ndirect;
+ int i, ai, aj, si, eind, ge, se;
+ int found, cj_ind_0, cj_ind_1, cj_ind_m;
+ int cj_m;
+ gmx_bool Found_si;
+ int si_ind;
nbnxn_excl_t *nbl_excl;
- int inner_i,inner_e;
+ int inner_i, inner_e;
cell = nbs->cell;
ndirect++;
}
}
-#ifdef NBNXN_SEARCH_BB_SSE
+#ifdef NBNXN_SEARCH_BB_SIMD4
else
{
while (cj_ind_first + ndirect <= cj_ind_last &&
- nbl->cj[cj_ind_first+ndirect].cj == ci_to_cj(na_cj_2log,ci) + ndirect)
+ nbl->cj[cj_ind_first+ndirect].cj == ci_to_cj(na_cj_2log, ci) + ndirect)
{
ndirect++;
}
#endif
/* Loop over the atoms in the i super-cell */
- for(i=0; i<nbl->na_sc; i++)
+ for (i = 0; i < nbl->na_sc; i++)
{
ai = nbs->a[ci*nbl->na_sc+i];
if (ai >= 0)
si = (i>>na_ci_2log);
/* Loop over the topology-based exclusions for this i-atom */
- for(eind=excl->index[ai]; eind<excl->index[ai+1]; eind++)
+ for (eind = excl->index[ai]; eind < excl->index[ai+1]; eind++)
{
aj = excl->a[eind];
else
{
/* Search for se using bisection */
- found = -1;
+ found = -1;
cj_ind_0 = cj_ind_first + ndirect;
cj_ind_1 = cj_ind_last + 1;
while (found == -1 && cj_ind_0 < cj_ind_1)
inner_e = ge - (se << na_cj_2log);
nbl->cj[found].excl &= ~(1U<<((inner_i<<na_cj_2log) + inner_e));
+/* The next code line is usually not needed. We do not want to version
+ * away the above line, because there is logic that relies on being
+ * able to detect easily whether any exclusions exist. */
+#if (defined GMX_SIMD_IBM_QPX)
+ nbl->cj[found].interaction_mask_indices[inner_i] &= ~(1U << inner_e);
+#endif
}
}
}
}
}
-/* Set all atom-pair exclusions from the topology stored in excl
- * as masks in the pair-list for i-super-cell entry nbl_sci
- */
-static void set_sci_top_excls(const nbnxn_search_t nbs,
- nbnxn_pairlist_t *nbl,
- gmx_bool diagRemoved,
- int na_c_2log,
- const nbnxn_sci_t *nbl_sci,
- const t_blocka *excl)
+/* Add a new i-entry to the FEP list and copy the i-properties */
+static gmx_inline void fep_list_new_nri_copy(t_nblist *nlist)
{
- const int *cell;
- int na_c;
- int sci;
- int cj_ind_first,cj_ind_last;
- int cj_first,cj_last;
- int ndirect;
- int i,ai,aj,si,eind,ge,se;
- int found,cj_ind_0,cj_ind_1,cj_ind_m;
- int cj_m;
- gmx_bool Found_si;
- int si_ind;
- nbnxn_excl_t *nbl_excl;
- int inner_i,inner_e,w;
+ /* Add a new i-entry */
+ nlist->nri++;
- cell = nbs->cell;
+ assert(nlist->nri < nlist->maxnri);
- na_c = nbl->na_ci;
+ /* Duplicate the last i-entry, except for jindex, which continues */
+ nlist->iinr[nlist->nri] = nlist->iinr[nlist->nri-1];
+ nlist->shift[nlist->nri] = nlist->shift[nlist->nri-1];
+ nlist->gid[nlist->nri] = nlist->gid[nlist->nri-1];
+ nlist->jindex[nlist->nri] = nlist->nrj;
+}
- if (nbl_sci->cj4_ind_end == nbl_sci->cj4_ind_start)
+/* For load balancing of the free-energy lists over threads, we set
+ * the maximum nrj size of an i-entry to 40. This leads to good
+ * load balancing in the worst case scenario of a single perturbed
+ * particle on 16 threads, while not introducing significant overhead.
+ * Note that half of the perturbed pairs will anyhow end up in very small lists,
+ * since non perturbed i-particles will see few perturbed j-particles).
+ */
+const int max_nrj_fep = 40;
+
+/* Exclude the perturbed pairs from the Verlet list. This is only done to avoid
+ * singularities for overlapping particles (0/0), since the charges and
+ * LJ parameters have been zeroed in the nbnxn data structure.
+ * Simultaneously make a group pair list for the perturbed pairs.
+ */
+static void make_fep_list(const nbnxn_search_t nbs,
+ const nbnxn_atomdata_t *nbat,
+ nbnxn_pairlist_t *nbl,
+ gmx_bool bDiagRemoved,
+ nbnxn_ci_t *nbl_ci,
+ const nbnxn_grid_t *gridi,
+ const nbnxn_grid_t *gridj,
+ t_nblist *nlist)
+{
+ int ci, cj_ind_start, cj_ind_end, cj_ind, cja, cjr;
+ int nri_max;
+ int ngid, gid_i = 0, gid_j, gid;
+ int egp_shift, egp_mask;
+ int gid_cj = 0;
+ int i, j, ind_i, ind_j, ai, aj;
+ int nri;
+ gmx_bool bFEP_i, bFEP_i_all;
+
+ if (nbl_ci->cj_ind_end == nbl_ci->cj_ind_start)
{
/* Empty list */
return;
}
- sci = nbl_sci->sci;
-
- cj_ind_first = nbl_sci->cj4_ind_start*NBNXN_GPU_JGROUP_SIZE;
- cj_ind_last = nbl->work->cj_ind - 1;
+ ci = nbl_ci->ci;
- cj_first = nbl->cj4[nbl_sci->cj4_ind_start].cj[0];
- cj_last = nbl_cj(nbl,cj_ind_last);
+ cj_ind_start = nbl_ci->cj_ind_start;
+ cj_ind_end = nbl_ci->cj_ind_end;
- /* Determine how many contiguous j-clusters we have starting
- * from the first i-cluster. This number can be used to directly
- * calculate j-cluster indices for excluded atoms.
+ /* In worst case we have alternating energy groups
+ * and create #atom-pair lists, which means we need the size
+ * of a cluster pair (na_ci*na_cj) times the number of cj's.
*/
- ndirect = 0;
- while (cj_ind_first + ndirect <= cj_ind_last &&
- nbl_cj(nbl,cj_ind_first+ndirect) == sci*GPU_NSUBCELL + ndirect)
+ nri_max = nbl->na_ci*nbl->na_cj*(cj_ind_end - cj_ind_start);
+ if (nlist->nri + nri_max > nlist->maxnri)
{
- ndirect++;
+ nlist->maxnri = over_alloc_large(nlist->nri + nri_max);
+ reallocate_nblist(nlist);
}
- /* Loop over the atoms in the i super-cell */
- for(i=0; i<nbl->na_sc; i++)
+ ngid = nbat->nenergrp;
+
+ if (ngid*gridj->na_cj > sizeof(gid_cj)*8)
{
- ai = nbs->a[sci*nbl->na_sc+i];
+ gmx_fatal(FARGS, "The Verlet scheme with %dx%d kernels and free-energy only supports up to %d energy groups",
+ gridi->na_c, gridj->na_cj, (sizeof(gid_cj)*8)/gridj->na_cj);
+ }
+
+ egp_shift = nbat->neg_2log;
+ egp_mask = (1<<nbat->neg_2log) - 1;
+
+ /* Loop over the atoms in the i sub-cell */
+ bFEP_i_all = TRUE;
+ for (i = 0; i < nbl->na_ci; i++)
+ {
+ ind_i = ci*nbl->na_ci + i;
+ ai = nbs->a[ind_i];
if (ai >= 0)
{
- si = (i>>na_c_2log);
+ nri = nlist->nri;
+ nlist->jindex[nri+1] = nlist->jindex[nri];
+ nlist->iinr[nri] = ai;
+ /* The actual energy group pair index is set later */
+ nlist->gid[nri] = 0;
+ nlist->shift[nri] = nbl_ci->shift & NBNXN_CI_SHIFT;
- /* Loop over the topology-based exclusions for this i-atom */
- for(eind=excl->index[ai]; eind<excl->index[ai+1]; eind++)
+ bFEP_i = gridi->fep[ci - gridi->cell0] & (1 << i);
+
+ bFEP_i_all = bFEP_i_all && bFEP_i;
+
+ if ((nlist->nrj + cj_ind_end - cj_ind_start)*nbl->na_cj > nlist->maxnrj)
{
- aj = excl->a[eind];
+ nlist->maxnrj = over_alloc_small((nlist->nrj + cj_ind_end - cj_ind_start)*nbl->na_cj);
+ srenew(nlist->jjnr, nlist->maxnrj);
+ srenew(nlist->excl_fep, nlist->maxnrj);
+ }
- if (aj == ai)
- {
- /* The self exclusion are already set, save some time */
- continue;
- }
+ if (ngid > 1)
+ {
+ gid_i = (nbat->energrp[ci] >> (egp_shift*i)) & egp_mask;
+ }
- ge = cell[aj];
+ for (cj_ind = cj_ind_start; cj_ind < cj_ind_end; cj_ind++)
+ {
+ unsigned int fep_cj;
- /* Without shifts we only calculate interactions j>i
- * for one-way pair-lists.
- */
- if (diagRemoved && ge <= sci*nbl->na_sc + i)
+ cja = nbl->cj[cj_ind].cj;
+
+ if (gridj->na_cj == gridj->na_c)
{
- continue;
+ cjr = cja - gridj->cell0;
+ fep_cj = gridj->fep[cjr];
+ if (ngid > 1)
+ {
+ gid_cj = nbat->energrp[cja];
+ }
}
-
- se = ge>>na_c_2log;
- /* Could the cluster se be in our list? */
- if (se >= cj_first && se <= cj_last)
+ else if (2*gridj->na_cj == gridj->na_c)
{
- if (se < cj_first + ndirect)
+ cjr = cja - gridj->cell0*2;
+ /* Extract half of the ci fep/energrp mask */
+ fep_cj = (gridj->fep[cjr>>1] >> ((cjr&1)*gridj->na_cj)) & ((1<<gridj->na_cj) - 1);
+ if (ngid > 1)
{
- /* We can calculate cj_ind directly from se */
- found = cj_ind_first + se - cj_first;
+ gid_cj = nbat->energrp[cja>>1] >> ((cja&1)*gridj->na_cj*egp_shift) & ((1<<(gridj->na_cj*egp_shift)) - 1);
}
- else
+ }
+ else
+ {
+ cjr = cja - (gridj->cell0>>1);
+ /* Combine two ci fep masks/energrp */
+ fep_cj = gridj->fep[cjr*2] + (gridj->fep[cjr*2+1] << gridj->na_c);
+ if (ngid > 1)
{
- /* Search for se using bisection */
- found = -1;
- cj_ind_0 = cj_ind_first + ndirect;
- cj_ind_1 = cj_ind_last + 1;
- while (found == -1 && cj_ind_0 < cj_ind_1)
- {
+ gid_cj = nbat->energrp[cja*2] + (nbat->energrp[cja*2+1] << (gridj->na_c*egp_shift));
+ }
+ }
+
+ if (bFEP_i || fep_cj != 0)
+ {
+ for (j = 0; j < nbl->na_cj; j++)
+ {
+ /* Is this interaction perturbed and not excluded? */
+ ind_j = cja*nbl->na_cj + j;
+ aj = nbs->a[ind_j];
+ if (aj >= 0 &&
+ (bFEP_i || (fep_cj & (1 << j))) &&
+ (!bDiagRemoved || ind_j >= ind_i))
+ {
+ if (ngid > 1)
+ {
+ gid_j = (gid_cj >> (j*egp_shift)) & egp_mask;
+ gid = GID(gid_i, gid_j, ngid);
+
+ if (nlist->nrj > nlist->jindex[nri] &&
+ nlist->gid[nri] != gid)
+ {
+ /* Energy group pair changed: new list */
+ fep_list_new_nri_copy(nlist);
+ nri = nlist->nri;
+ }
+ nlist->gid[nri] = gid;
+ }
+
+ if (nlist->nrj - nlist->jindex[nri] >= max_nrj_fep)
+ {
+ fep_list_new_nri_copy(nlist);
+ nri = nlist->nri;
+ }
+
+ /* Add it to the FEP list */
+ nlist->jjnr[nlist->nrj] = aj;
+ nlist->excl_fep[nlist->nrj] = (nbl->cj[cj_ind].excl >> (i*nbl->na_cj + j)) & 1;
+ nlist->nrj++;
+
+ /* Exclude it from the normal list.
+ * Note that the charge has been set to zero,
+ * but we need to avoid 0/0, as perturbed atoms
+ * can be on top of each other.
+ */
+ nbl->cj[cj_ind].excl &= ~(1U << (i*nbl->na_cj + j));
+ }
+ }
+ }
+ }
+
+ if (nlist->nrj > nlist->jindex[nri])
+ {
+ /* Actually add this new, non-empty, list */
+ nlist->nri++;
+ nlist->jindex[nlist->nri] = nlist->nrj;
+ }
+ }
+ }
+
+ if (bFEP_i_all)
+ {
+ /* All interactions are perturbed, we can skip this entry */
+ nbl_ci->cj_ind_end = cj_ind_start;
+ }
+}
+
+/* Return the index of atom a within a cluster */
+static gmx_inline int cj_mod_cj4(int cj)
+{
+ return cj & (NBNXN_GPU_JGROUP_SIZE - 1);
+}
+
+/* Convert a j-cluster to a cj4 group */
+static gmx_inline int cj_to_cj4(int cj)
+{
+ return cj >> NBNXN_GPU_JGROUP_SIZE_2LOG;
+}
+
+/* Return the index of an j-atom within a warp */
+static gmx_inline int a_mod_wj(int a)
+{
+ return a & (NBNXN_GPU_CLUSTER_SIZE/2 - 1);
+}
+
+/* As make_fep_list above, but for super/sub lists. */
+static void make_fep_list_supersub(const nbnxn_search_t nbs,
+ const nbnxn_atomdata_t *nbat,
+ nbnxn_pairlist_t *nbl,
+ gmx_bool bDiagRemoved,
+ const nbnxn_sci_t *nbl_sci,
+ real shx,
+ real shy,
+ real shz,
+ real rlist_fep2,
+ const nbnxn_grid_t *gridi,
+ const nbnxn_grid_t *gridj,
+ t_nblist *nlist)
+{
+ int sci, cj4_ind_start, cj4_ind_end, cj4_ind, gcj, cjr;
+ int nri_max;
+ int c, c_abs;
+ int i, j, ind_i, ind_j, ai, aj;
+ int nri;
+ gmx_bool bFEP_i;
+ real xi, yi, zi;
+ const nbnxn_cj4_t *cj4;
+
+ if (nbl_sci->cj4_ind_end == nbl_sci->cj4_ind_start)
+ {
+ /* Empty list */
+ return;
+ }
+
+ sci = nbl_sci->sci;
+
+ cj4_ind_start = nbl_sci->cj4_ind_start;
+ cj4_ind_end = nbl_sci->cj4_ind_end;
+
+ /* Here we process one super-cell, max #atoms na_sc, versus a list
+ * cj4 entries, each with max NBNXN_GPU_JGROUP_SIZE cj's, each
+ * of size na_cj atoms.
+ * On the GPU we don't support energy groups (yet).
+ * So for each of the na_sc i-atoms, we need max one FEP list
+ * for each max_nrj_fep j-atoms.
+ */
+ nri_max = nbl->na_sc*nbl->na_cj*(1 + ((cj4_ind_end - cj4_ind_start)*NBNXN_GPU_JGROUP_SIZE)/max_nrj_fep);
+ if (nlist->nri + nri_max > nlist->maxnri)
+ {
+ nlist->maxnri = over_alloc_large(nlist->nri + nri_max);
+ reallocate_nblist(nlist);
+ }
+
+ /* Loop over the atoms in the i super-cluster */
+ for (c = 0; c < GPU_NSUBCELL; c++)
+ {
+ c_abs = sci*GPU_NSUBCELL + c;
+
+ for (i = 0; i < nbl->na_ci; i++)
+ {
+ ind_i = c_abs*nbl->na_ci + i;
+ ai = nbs->a[ind_i];
+ if (ai >= 0)
+ {
+ nri = nlist->nri;
+ nlist->jindex[nri+1] = nlist->jindex[nri];
+ nlist->iinr[nri] = ai;
+ /* With GPUs, energy groups are not supported */
+ nlist->gid[nri] = 0;
+ nlist->shift[nri] = nbl_sci->shift & NBNXN_CI_SHIFT;
+
+ bFEP_i = (gridi->fep[c_abs - gridi->cell0] & (1 << i));
+
+ xi = nbat->x[ind_i*nbat->xstride+XX] + shx;
+ yi = nbat->x[ind_i*nbat->xstride+YY] + shy;
+ zi = nbat->x[ind_i*nbat->xstride+ZZ] + shz;
+
+ if ((nlist->nrj + cj4_ind_end - cj4_ind_start)*NBNXN_GPU_JGROUP_SIZE*nbl->na_cj > nlist->maxnrj)
+ {
+ nlist->maxnrj = over_alloc_small((nlist->nrj + cj4_ind_end - cj4_ind_start)*NBNXN_GPU_JGROUP_SIZE*nbl->na_cj);
+ srenew(nlist->jjnr, nlist->maxnrj);
+ srenew(nlist->excl_fep, nlist->maxnrj);
+ }
+
+ for (cj4_ind = cj4_ind_start; cj4_ind < cj4_ind_end; cj4_ind++)
+ {
+ cj4 = &nbl->cj4[cj4_ind];
+
+ for (gcj = 0; gcj < NBNXN_GPU_JGROUP_SIZE; gcj++)
+ {
+ unsigned int fep_cj;
+
+ if ((cj4->imei[0].imask & (1U << (gcj*GPU_NSUBCELL + c))) == 0)
+ {
+ /* Skip this ci for this cj */
+ continue;
+ }
+
+ cjr = cj4->cj[gcj] - gridj->cell0*GPU_NSUBCELL;
+
+ fep_cj = gridj->fep[cjr];
+
+ if (bFEP_i || fep_cj != 0)
+ {
+ for (j = 0; j < nbl->na_cj; j++)
+ {
+ /* Is this interaction perturbed and not excluded? */
+ ind_j = (gridj->cell0*GPU_NSUBCELL + cjr)*nbl->na_cj + j;
+ aj = nbs->a[ind_j];
+ if (aj >= 0 &&
+ (bFEP_i || (fep_cj & (1 << j))) &&
+ (!bDiagRemoved || ind_j >= ind_i))
+ {
+ nbnxn_excl_t *excl;
+ int excl_pair;
+ unsigned int excl_bit;
+ real dx, dy, dz;
+
+ get_nbl_exclusions_1(nbl, cj4_ind, j>>2, &excl);
+
+ excl_pair = a_mod_wj(j)*nbl->na_ci + i;
+ excl_bit = (1U << (gcj*GPU_NSUBCELL + c));
+
+ dx = nbat->x[ind_j*nbat->xstride+XX] - xi;
+ dy = nbat->x[ind_j*nbat->xstride+YY] - yi;
+ dz = nbat->x[ind_j*nbat->xstride+ZZ] - zi;
+
+ /* The unpruned GPU list has more than 2/3
+ * of the atom pairs beyond rlist. Using
+ * this list will cause a lot of overhead
+ * in the CPU FEP kernels, especially
+ * relative to the fast GPU kernels.
+ * So we prune the FEP list here.
+ */
+ if (dx*dx + dy*dy + dz*dz < rlist_fep2)
+ {
+ if (nlist->nrj - nlist->jindex[nri] >= max_nrj_fep)
+ {
+ fep_list_new_nri_copy(nlist);
+ nri = nlist->nri;
+ }
+
+ /* Add it to the FEP list */
+ nlist->jjnr[nlist->nrj] = aj;
+ nlist->excl_fep[nlist->nrj] = (excl->pair[excl_pair] & excl_bit) ? 1 : 0;
+ nlist->nrj++;
+ }
+
+ /* Exclude it from the normal list.
+ * Note that the charge and LJ parameters have
+ * been set to zero, but we need to avoid 0/0,
+ * as perturbed atoms can be on top of each other.
+ */
+ excl->pair[excl_pair] &= ~excl_bit;
+ }
+ }
+
+ /* Note that we could mask out this pair in imask
+ * if all i- and/or all j-particles are perturbed.
+ * But since the perturbed pairs on the CPU will
+ * take an order of magnitude more time, the GPU
+ * will finish before the CPU and there is no gain.
+ */
+ }
+ }
+ }
+
+ if (nlist->nrj > nlist->jindex[nri])
+ {
+ /* Actually add this new, non-empty, list */
+ nlist->nri++;
+ nlist->jindex[nlist->nri] = nlist->nrj;
+ }
+ }
+ }
+ }
+}
+
+/* Set all atom-pair exclusions from the topology stored in excl
+ * as masks in the pair-list for i-super-cell entry nbl_sci
+ */
+static void set_sci_top_excls(const nbnxn_search_t nbs,
+ nbnxn_pairlist_t *nbl,
+ gmx_bool diagRemoved,
+ int na_c_2log,
+ const nbnxn_sci_t *nbl_sci,
+ const t_blocka *excl)
+{
+ const int *cell;
+ int na_c;
+ int sci;
+ int cj_ind_first, cj_ind_last;
+ int cj_first, cj_last;
+ int ndirect;
+ int i, ai, aj, si, eind, ge, se;
+ int found, cj_ind_0, cj_ind_1, cj_ind_m;
+ int cj_m;
+ gmx_bool Found_si;
+ int si_ind;
+ nbnxn_excl_t *nbl_excl;
+ int inner_i, inner_e, w;
+
+ cell = nbs->cell;
+
+ na_c = nbl->na_ci;
+
+ if (nbl_sci->cj4_ind_end == nbl_sci->cj4_ind_start)
+ {
+ /* Empty list */
+ return;
+ }
+
+ sci = nbl_sci->sci;
+
+ cj_ind_first = nbl_sci->cj4_ind_start*NBNXN_GPU_JGROUP_SIZE;
+ cj_ind_last = nbl->work->cj_ind - 1;
+
+ cj_first = nbl->cj4[nbl_sci->cj4_ind_start].cj[0];
+ cj_last = nbl_cj(nbl, cj_ind_last);
+
+ /* Determine how many contiguous j-clusters we have starting
+ * from the first i-cluster. This number can be used to directly
+ * calculate j-cluster indices for excluded atoms.
+ */
+ ndirect = 0;
+ while (cj_ind_first + ndirect <= cj_ind_last &&
+ nbl_cj(nbl, cj_ind_first+ndirect) == sci*GPU_NSUBCELL + ndirect)
+ {
+ ndirect++;
+ }
+
+ /* Loop over the atoms in the i super-cell */
+ for (i = 0; i < nbl->na_sc; i++)
+ {
+ ai = nbs->a[sci*nbl->na_sc+i];
+ if (ai >= 0)
+ {
+ si = (i>>na_c_2log);
+
+ /* Loop over the topology-based exclusions for this i-atom */
+ for (eind = excl->index[ai]; eind < excl->index[ai+1]; eind++)
+ {
+ aj = excl->a[eind];
+
+ if (aj == ai)
+ {
+ /* The self exclusion are already set, save some time */
+ continue;
+ }
+
+ ge = cell[aj];
+
+ /* Without shifts we only calculate interactions j>i
+ * for one-way pair-lists.
+ */
+ if (diagRemoved && ge <= sci*nbl->na_sc + i)
+ {
+ continue;
+ }
+
+ se = ge>>na_c_2log;
+ /* Could the cluster se be in our list? */
+ if (se >= cj_first && se <= cj_last)
+ {
+ if (se < cj_first + ndirect)
+ {
+ /* We can calculate cj_ind directly from se */
+ found = cj_ind_first + se - cj_first;
+ }
+ else
+ {
+ /* Search for se using bisection */
+ found = -1;
+ cj_ind_0 = cj_ind_first + ndirect;
+ cj_ind_1 = cj_ind_last + 1;
+ while (found == -1 && cj_ind_0 < cj_ind_1)
+ {
cj_ind_m = (cj_ind_0 + cj_ind_1)>>1;
- cj_m = nbl_cj(nbl,cj_ind_m);
+ cj_m = nbl_cj(nbl, cj_ind_m);
if (se == cj_m)
{
inner_i = i - si*na_c;
inner_e = ge - se*na_c;
-/* Macro for getting the index of atom a within a cluster */
-#define AMODCJ4(a) ((a) & (NBNXN_GPU_JGROUP_SIZE - 1))
-/* Macro for converting an atom number to a cluster number */
-#define A2CJ4(a) ((a) >> NBNXN_GPU_JGROUP_SIZE_2LOG)
-/* Macro for getting the index of an i-atom within a warp */
-#define AMODWI(a) ((a) & (NBNXN_GPU_CLUSTER_SIZE/2 - 1))
-
- if (nbl_imask0(nbl,found) & (1U << (AMODCJ4(found)*GPU_NSUBCELL + si)))
+ if (nbl_imask0(nbl, found) & (1U << (cj_mod_cj4(found)*GPU_NSUBCELL + si)))
{
w = (inner_e >> 2);
- get_nbl_exclusions_1(nbl,A2CJ4(found),w,&nbl_excl);
+ get_nbl_exclusions_1(nbl, cj_to_cj4(found), w, &nbl_excl);
- nbl_excl->pair[AMODWI(inner_e)*nbl->na_ci+inner_i] &=
- ~(1U << (AMODCJ4(found)*GPU_NSUBCELL + si));
+ nbl_excl->pair[a_mod_wj(inner_e)*nbl->na_ci+inner_i] &=
+ ~(1U << (cj_mod_cj4(found)*GPU_NSUBCELL + si));
}
-
-#undef AMODCJ4
-#undef A2CJ4
-#undef AMODWI
}
}
}
}
/* Reallocate the simple ci list for at least n entries */
-static void nb_realloc_ci(nbnxn_pairlist_t *nbl,int n)
+static void nb_realloc_ci(nbnxn_pairlist_t *nbl, int n)
{
nbl->ci_nalloc = over_alloc_small(n);
nbnxn_realloc_void((void **)&nbl->ci,
nbl->nci*sizeof(*nbl->ci),
nbl->ci_nalloc*sizeof(*nbl->ci),
- nbl->alloc,nbl->free);
+ nbl->alloc, nbl->free);
}
/* Reallocate the super-cell sci list for at least n entries */
-static void nb_realloc_sci(nbnxn_pairlist_t *nbl,int n)
+static void nb_realloc_sci(nbnxn_pairlist_t *nbl, int n)
{
nbl->sci_nalloc = over_alloc_small(n);
nbnxn_realloc_void((void **)&nbl->sci,
nbl->nsci*sizeof(*nbl->sci),
nbl->sci_nalloc*sizeof(*nbl->sci),
- nbl->alloc,nbl->free);
+ nbl->alloc, nbl->free);
}
/* Make a new ci entry at index nbl->nci */
-static void new_ci_entry(nbnxn_pairlist_t *nbl,int ci,int shift,int flags,
- nbnxn_list_work_t *work)
+static void new_ci_entry(nbnxn_pairlist_t *nbl, int ci, int shift, int flags)
{
if (nbl->nci + 1 > nbl->ci_nalloc)
{
- nb_realloc_ci(nbl,nbl->nci+1);
+ nb_realloc_ci(nbl, nbl->nci+1);
}
nbl->ci[nbl->nci].ci = ci;
nbl->ci[nbl->nci].shift = shift;
}
/* Make a new sci entry at index nbl->nsci */
-static void new_sci_entry(nbnxn_pairlist_t *nbl,int sci,int shift,int flags,
- nbnxn_list_work_t *work)
+static void new_sci_entry(nbnxn_pairlist_t *nbl, int sci, int shift)
{
if (nbl->nsci + 1 > nbl->sci_nalloc)
{
- nb_realloc_sci(nbl,nbl->nsci+1);
+ nb_realloc_sci(nbl, nbl->nsci+1);
}
nbl->sci[nbl->nsci].sci = sci;
nbl->sci[nbl->nsci].shift = shift;
/* Sort the simple j-list cj on exclusions.
* Entries with exclusions will all be sorted to the beginning of the list.
*/
-static void sort_cj_excl(nbnxn_cj_t *cj,int ncj,
+static void sort_cj_excl(nbnxn_cj_t *cj, int ncj,
nbnxn_list_work_t *work)
{
- int jnew,j;
+ int jnew, j;
if (ncj > work->cj_nalloc)
{
work->cj_nalloc = over_alloc_large(ncj);
- srenew(work->cj,work->cj_nalloc);
+ srenew(work->cj, work->cj_nalloc);
}
/* Make a list of the j-cells involving exclusions */
jnew = 0;
- for(j=0; j<ncj; j++)
+ for (j = 0; j < ncj; j++)
{
- if (cj[j].excl != NBNXN_INT_MASK_ALL)
+ if (cj[j].excl != NBNXN_INTERACTION_MASK_ALL)
{
work->cj[jnew++] = cj[j];
}
}
/* Check if there are exclusions at all or not just the first entry */
if (!((jnew == 0) ||
- (jnew == 1 && cj[0].excl != NBNXN_INT_MASK_ALL)))
+ (jnew == 1 && cj[0].excl != NBNXN_INTERACTION_MASK_ALL)))
{
- for(j=0; j<ncj; j++)
+ for (j = 0; j < ncj; j++)
{
- if (cj[j].excl == NBNXN_INT_MASK_ALL)
+ if (cj[j].excl == NBNXN_INTERACTION_MASK_ALL)
{
work->cj[jnew++] = cj[j];
}
}
- for(j=0; j<ncj; j++)
+ for (j = 0; j < ncj; j++)
{
cj[j] = work->cj[j];
}
jlen = nbl->ci[nbl->nci].cj_ind_end - nbl->ci[nbl->nci].cj_ind_start;
if (jlen > 0)
{
- sort_cj_excl(nbl->cj+nbl->ci[nbl->nci].cj_ind_start,jlen,nbl->work);
+ sort_cj_excl(nbl->cj+nbl->ci[nbl->nci].cj_ind_start, jlen, nbl->work);
/* The counts below are used for non-bonded pair/flop counts
* and should therefore match the available kernel setups.
}
/* Split sci entry for load balancing on the GPU.
- * As we only now the current count on our own thread,
+ * Splitting ensures we have enough lists to fully utilize the whole GPU.
+ * With progBal we generate progressively smaller lists, which improves
+ * load balancing. As we only know the current count on our own thread,
* we will need to estimate the current total amount of i-entries.
* As the lists get concatenated later, this estimate depends
- * both on nthread and our own thread index thread.
+ * both on nthread and our own thread index.
*/
static void split_sci_entry(nbnxn_pairlist_t *nbl,
- int nsp_max_av,gmx_bool progBal,int nc_bal,
- int thread,int nthread)
+ int nsp_max_av, gmx_bool progBal, int nc_bal,
+ int thread, int nthread)
{
int nsci_est;
int nsp_max;
- int cj4_start,cj4_end,j4len,cj4;
+ int cj4_start, cj4_end, j4len, cj4;
int sci;
- int nsp,nsp_sci,nsp_cj4,nsp_cj4_e,nsp_cj4_p;
+ int nsp, nsp_sci, nsp_cj4, nsp_cj4_e, nsp_cj4_p;
int p;
- /* Estimate the total numbers of ci's of the nblist combined
- * over all threads using the target number of ci's.
- */
- nsci_est = nc_bal*thread/nthread + nbl->nsci;
if (progBal)
{
+ /* Estimate the total numbers of ci's of the nblist combined
+ * over all threads using the target number of ci's.
+ */
+ nsci_est = nc_bal*thread/nthread + nbl->nsci;
+
/* The first ci blocks should be larger, to avoid overhead.
* The last ci blocks should be smaller, to improve load balancing.
*/
cj4_start = nbl->sci[nbl->nsci-1].cj4_ind_start;
cj4_end = nbl->sci[nbl->nsci-1].cj4_ind_end;
- j4len = cj4_end - cj4_start;
+ j4len = cj4_end - cj4_start;
if (j4len > 1 && j4len*GPU_NSUBCELL*NBNXN_GPU_JGROUP_SIZE > nsp_max)
{
nbl->nsci -= 1;
sci = nbl->nsci;
- cj4 = cj4_start;
nsp = 0;
nsp_sci = 0;
nsp_cj4_e = 0;
nsp_cj4 = 0;
- while (cj4 < cj4_end)
+ for (cj4 = cj4_start; cj4 < cj4_end; cj4++)
{
nsp_cj4_p = nsp_cj4;
+ /* Count the number of cluster pairs in this cj4 group */
nsp_cj4 = 0;
- for(p=0; p<GPU_NSUBCELL*NBNXN_GPU_JGROUP_SIZE; p++)
+ for (p = 0; p < GPU_NSUBCELL*NBNXN_GPU_JGROUP_SIZE; p++)
{
nsp_cj4 += (nbl->cj4[cj4].imei[0].imask >> p) & 1;
}
- nsp += nsp_cj4;
- if (nsp > nsp_max && nsp > nsp_cj4)
+ if (nsp_cj4 > 0 && nsp + nsp_cj4 > nsp_max)
{
+ /* Split the list at cj4 */
nbl->sci[sci].cj4_ind_end = cj4;
+ /* Create a new sci entry */
sci++;
nbl->nsci++;
if (nbl->nsci+1 > nbl->sci_nalloc)
{
- nb_realloc_sci(nbl,nbl->nsci+1);
+ nb_realloc_sci(nbl, nbl->nsci+1);
}
nbl->sci[sci].sci = nbl->sci[nbl->nsci-1].sci;
nbl->sci[sci].shift = nbl->sci[nbl->nsci-1].shift;
nbl->sci[sci].cj4_ind_start = cj4;
- nsp_sci = nsp - nsp_cj4;
- nsp_cj4_e = nsp_cj4_p;
- nsp = nsp_cj4;
+ nsp_sci = nsp;
+ nsp_cj4_e = nsp_cj4_p;
+ nsp = 0;
}
-
- cj4++;
+ nsp += nsp_cj4;
}
- /* Put the remaining cj4's in a new ci entry */
+ /* Put the remaining cj4's in the last sci entry */
nbl->sci[sci].cj4_ind_end = cj4_end;
- /* Possibly balance out the last two ci's
- * by moving the last cj4 of the second last ci.
+ /* Possibly balance out the last two sci's
+ * by moving the last cj4 of the second last sci.
*/
if (nsp_sci - nsp_cj4_e >= nsp + nsp_cj4_e)
{
nbl->sci[sci].cj4_ind_start--;
}
- sci++;
nbl->nsci++;
}
}
/* Clost this super/sub list i entry */
static void close_ci_entry_supersub(nbnxn_pairlist_t *nbl,
int nsp_max_av,
- gmx_bool progBal,int nc_bal,
- int thread,int nthread)
+ gmx_bool progBal, int nc_bal,
+ int thread, int nthread)
{
- int j4len,tlen;
- int nb,b;
+ int j4len, tlen;
+ int nb, b;
/* All content of the new ci entry have already been filled correctly,
* we only need to increase the count here (for non empty lists).
if (nsp_max_av > 0)
{
- split_sci_entry(nbl,nsp_max_av,progBal,nc_bal,thread,nthread);
+ /* Measure the size of the new entry and potentially split it */
+ split_sci_entry(nbl, nsp_max_av, progBal, nc_bal, thread, nthread);
}
}
}
nbl->work->ncj_hlj = 0;
}
+/* Clears a group scheme pair list */
+static void clear_pairlist_fep(t_nblist *nl)
+{
+ nl->nri = 0;
+ nl->nrj = 0;
+ if (nl->jindex == NULL)
+ {
+ snew(nl->jindex, 1);
+ }
+ nl->jindex[0] = 0;
+}
+
/* Sets a simple list i-cell bounding box, including PBC shift */
-static void set_icell_bb_simple(const float *bb,int ci,
- real shx,real shy,real shz,
- float *bb_ci)
+static gmx_inline void set_icell_bb_simple(const nbnxn_bb_t *bb, int ci,
+ real shx, real shy, real shz,
+ nbnxn_bb_t *bb_ci)
{
- int ia;
-
- ia = ci*NNBSBB_B;
- bb_ci[BBL_X] = bb[ia+BBL_X] + shx;
- bb_ci[BBL_Y] = bb[ia+BBL_Y] + shy;
- bb_ci[BBL_Z] = bb[ia+BBL_Z] + shz;
- bb_ci[BBU_X] = bb[ia+BBU_X] + shx;
- bb_ci[BBU_Y] = bb[ia+BBU_Y] + shy;
- bb_ci[BBU_Z] = bb[ia+BBU_Z] + shz;
+ bb_ci->lower[BB_X] = bb[ci].lower[BB_X] + shx;
+ bb_ci->lower[BB_Y] = bb[ci].lower[BB_Y] + shy;
+ bb_ci->lower[BB_Z] = bb[ci].lower[BB_Z] + shz;
+ bb_ci->upper[BB_X] = bb[ci].upper[BB_X] + shx;
+ bb_ci->upper[BB_Y] = bb[ci].upper[BB_Y] + shy;
+ bb_ci->upper[BB_Z] = bb[ci].upper[BB_Z] + shz;
}
+#ifdef NBNXN_BBXXXX
/* Sets a super-cell and sub cell bounding boxes, including PBC shift */
-static void set_icell_bb_supersub(const float *bb,int ci,
- real shx,real shy,real shz,
- float *bb_ci)
+static void set_icell_bbxxxx_supersub(const float *bb, int ci,
+ real shx, real shy, real shz,
+ float *bb_ci)
{
- int ia,m,i;
+ int ia, m, i;
-#ifdef NBNXN_BBXXXX
ia = ci*(GPU_NSUBCELL>>STRIDE_PBB_2LOG)*NNBSBB_XXXX;
- for(m=0; m<(GPU_NSUBCELL>>STRIDE_PBB_2LOG)*NNBSBB_XXXX; m+=NNBSBB_XXXX)
+ for (m = 0; m < (GPU_NSUBCELL>>STRIDE_PBB_2LOG)*NNBSBB_XXXX; m += NNBSBB_XXXX)
{
- for(i=0; i<STRIDE_PBB; i++)
+ for (i = 0; i < STRIDE_PBB; i++)
{
bb_ci[m+0*STRIDE_PBB+i] = bb[ia+m+0*STRIDE_PBB+i] + shx;
bb_ci[m+1*STRIDE_PBB+i] = bb[ia+m+1*STRIDE_PBB+i] + shy;
bb_ci[m+5*STRIDE_PBB+i] = bb[ia+m+5*STRIDE_PBB+i] + shz;
}
}
-#else
- ia = ci*GPU_NSUBCELL*NNBSBB_B;
- for(i=0; i<GPU_NSUBCELL*NNBSBB_B; i+=NNBSBB_B)
+}
+#endif
+
+/* Sets a super-cell and sub cell bounding boxes, including PBC shift */
+static void set_icell_bb_supersub(const nbnxn_bb_t *bb, int ci,
+ real shx, real shy, real shz,
+ nbnxn_bb_t *bb_ci)
+{
+ int i;
+
+ for (i = 0; i < GPU_NSUBCELL; i++)
{
- bb_ci[i+BBL_X] = bb[ia+i+BBL_X] + shx;
- bb_ci[i+BBL_Y] = bb[ia+i+BBL_Y] + shy;
- bb_ci[i+BBL_Z] = bb[ia+i+BBL_Z] + shz;
- bb_ci[i+BBU_X] = bb[ia+i+BBU_X] + shx;
- bb_ci[i+BBU_Y] = bb[ia+i+BBU_Y] + shy;
- bb_ci[i+BBU_Z] = bb[ia+i+BBU_Z] + shz;
+ set_icell_bb_simple(bb, ci*GPU_NSUBCELL+i,
+ shx, shy, shz,
+ &bb_ci[i]);
}
-#endif
}
/* Copies PBC shifted i-cell atom coordinates x,y,z to working array */
static void icell_set_x_simple(int ci,
- real shx,real shy,real shz,
- int na_c,
- int stride,const real *x,
+ real shx, real shy, real shz,
+ int gmx_unused na_c,
+ int stride, const real *x,
nbnxn_list_work_t *work)
{
- int ia,i;
+ int ia, i;
ia = ci*NBNXN_CPU_CLUSTER_I_SIZE;
- for(i=0; i<NBNXN_CPU_CLUSTER_I_SIZE; i++)
+ for (i = 0; i < NBNXN_CPU_CLUSTER_I_SIZE; i++)
{
work->x_ci[i*STRIDE_XYZ+XX] = x[(ia+i)*stride+XX] + shx;
work->x_ci[i*STRIDE_XYZ+YY] = x[(ia+i)*stride+YY] + shy;
/* Copies PBC shifted super-cell atom coordinates x,y,z to working array */
static void icell_set_x_supersub(int ci,
- real shx,real shy,real shz,
+ real shx, real shy, real shz,
int na_c,
- int stride,const real *x,
+ int stride, const real *x,
nbnxn_list_work_t *work)
{
- int ia,i;
+ int ia, i;
real *x_ci;
x_ci = work->x_ci;
ia = ci*GPU_NSUBCELL*na_c;
- for(i=0; i<GPU_NSUBCELL*na_c; i++)
+ for (i = 0; i < GPU_NSUBCELL*na_c; i++)
{
x_ci[i*DIM + XX] = x[(ia+i)*stride + XX] + shx;
x_ci[i*DIM + YY] = x[(ia+i)*stride + YY] + shy;
}
}
-#ifdef NBNXN_SEARCH_BB_SSE
+#ifdef NBNXN_SEARCH_BB_SIMD4
/* Copies PBC shifted super-cell packed atom coordinates to working array */
-static void icell_set_x_supersub_sse8(int ci,
- real shx,real shy,real shz,
- int na_c,
- int stride,const real *x,
- nbnxn_list_work_t *work)
+static void icell_set_x_supersub_simd4(int ci,
+ real shx, real shy, real shz,
+ int na_c,
+ int stride, const real *x,
+ nbnxn_list_work_t *work)
{
- int si,io,ia,i,j;
+ int si, io, ia, i, j;
real *x_ci;
x_ci = work->x_ci;
- for(si=0; si<GPU_NSUBCELL; si++)
+ for (si = 0; si < GPU_NSUBCELL; si++)
{
- for(i=0; i<na_c; i+=STRIDE_PBB)
+ for (i = 0; i < na_c; i += STRIDE_PBB)
{
io = si*na_c + i;
ia = ci*GPU_NSUBCELL*na_c + io;
- for(j=0; j<STRIDE_PBB; j++)
+ for (j = 0; j < STRIDE_PBB; j++)
{
x_ci[io*DIM + j + XX*STRIDE_PBB] = x[(ia+j)*stride+XX] + shx;
x_ci[io*DIM + j + YY*STRIDE_PBB] = x[(ia+j)*stride+YY] + shy;
}
#endif
-static real nbnxn_rlist_inc_nonloc_fac = 0.6;
+static real minimum_subgrid_size_xy(const nbnxn_grid_t *grid)
+{
+ if (grid->bSimple)
+ {
+ return min(grid->sx, grid->sy);
+ }
+ else
+ {
+ return min(grid->sx/GPU_NSUBCELL_X, grid->sy/GPU_NSUBCELL_Y);
+ }
+}
+
+static real effective_buffer_1x1_vs_MxN(const nbnxn_grid_t *gridi,
+ const nbnxn_grid_t *gridj)
+{
+ const real eff_1x1_buffer_fac_overest = 0.1;
+
+ /* Determine an atom-pair list cut-off buffer size for atom pairs,
+ * to be added to rlist (including buffer) used for MxN.
+ * This is for converting an MxN list to a 1x1 list. This means we can't
+ * use the normal buffer estimate, as we have an MxN list in which
+ * some atom pairs beyond rlist are missing. We want to capture
+ * the beneficial effect of buffering by extra pairs just outside rlist,
+ * while removing the useless pairs that are further away from rlist.
+ * (Also the buffer could have been set manually not using the estimate.)
+ * This buffer size is an overestimate.
+ * We add 10% of the smallest grid sub-cell dimensions.
+ * Note that the z-size differs per cell and we don't use this,
+ * so we overestimate.
+ * With PME, the 10% value gives a buffer that is somewhat larger
+ * than the effective buffer with a tolerance of 0.005 kJ/mol/ps.
+ * Smaller tolerances or using RF lead to a smaller effective buffer,
+ * so 10% gives a safe overestimate.
+ */
+ return eff_1x1_buffer_fac_overest*(minimum_subgrid_size_xy(gridi) +
+ minimum_subgrid_size_xy(gridj));
+}
+
+/* Clusters at the cut-off only increase rlist by 60% of their size */
+static real nbnxn_rlist_inc_outside_fac = 0.6;
/* Due to the cluster size the effective pair-list is longer than
* that of a simple atom pair-list. This function gives the extra distance.
*/
-real nbnxn_get_rlist_effective_inc(int cluster_size,real atom_density)
+real nbnxn_get_rlist_effective_inc(int cluster_size_j, real atom_density)
{
- return ((0.5 + nbnxn_rlist_inc_nonloc_fac)*sqr(((cluster_size) - 1.0)/(cluster_size))*pow((cluster_size)/(atom_density),1.0/3.0));
+ int cluster_size_i;
+ real vol_inc_i, vol_inc_j;
+
+ /* We should get this from the setup, but currently it's the same for
+ * all setups, including GPUs.
+ */
+ cluster_size_i = NBNXN_CPU_CLUSTER_I_SIZE;
+
+ vol_inc_i = (cluster_size_i - 1)/atom_density;
+ vol_inc_j = (cluster_size_j - 1)/atom_density;
+
+ return nbnxn_rlist_inc_outside_fac*pow(vol_inc_i + vol_inc_j, 1.0/3.0);
}
/* Estimates the interaction volume^2 for non-local interactions */
-static real nonlocal_vol2(const gmx_domdec_zones_t *zones,rvec ls,real r)
+static real nonlocal_vol2(const gmx_domdec_zones_t *zones, rvec ls, real r)
{
- int z,d;
- real cl,ca,za;
+ int z, d;
+ real cl, ca, za;
real vold_est;
real vol2_est_tot;
* as small parts as possible.
*/
- for(z=0; z<zones->n; z++)
+ for (z = 0; z < zones->n; z++)
{
if (zones->shift[z][XX] + zones->shift[z][YY] + zones->shift[z][ZZ] == 1)
{
cl = 0;
ca = 1;
za = 1;
- for(d=0; d<DIM; d++)
+ for (d = 0; d < DIM; d++)
{
if (zones->shift[z][d] == 0)
{
/* Estimates the average size of a full j-list for super/sub setup */
static int get_nsubpair_max(const nbnxn_search_t nbs,
- int iloc,
- real rlist,
- int min_ci_balanced)
+ int iloc,
+ real rlist,
+ int min_ci_balanced)
{
const nbnxn_grid_t *grid;
- rvec ls;
- real xy_diag2,r_eff_sup,vol_est,nsp_est,nsp_est_nl;
- int nsubpair_max;
+ rvec ls;
+ real xy_diag2, r_eff_sup, vol_est, nsp_est, nsp_est_nl;
+ int nsubpair_max;
grid = &nbs->grid[0];
xy_diag2 = ls[XX]*ls[XX] + ls[YY]*ls[YY] + ls[ZZ]*ls[ZZ];
/* The formulas below are a heuristic estimate of the average nsj per si*/
- r_eff_sup = rlist + nbnxn_rlist_inc_nonloc_fac*sqr((grid->na_c - 1.0)/grid->na_c)*sqrt(xy_diag2/3);
+ r_eff_sup = rlist + nbnxn_rlist_inc_outside_fac*sqr((grid->na_c - 1.0)/grid->na_c)*sqrt(xy_diag2/3);
if (!nbs->DomDec || nbs->zones->n == 1)
{
{
nsp_est_nl =
sqr(grid->atom_density/grid->na_c)*
- nonlocal_vol2(nbs->zones,ls,r_eff_sup);
+ nonlocal_vol2(nbs->zones, ls, r_eff_sup);
}
if (LOCAL_I(iloc))
/* 12/2 quarter pie slices on the edges */
vol_est += 2*(ls[XX] + ls[YY] + ls[ZZ])*0.25*M_PI*sqr(r_eff_sup);
/* 4 octants of a sphere */
- vol_est += 0.5*4.0/3.0*M_PI*pow(r_eff_sup,3);
+ vol_est += 0.5*4.0/3.0*M_PI*pow(r_eff_sup, 3);
nsp_est = grid->nsubc_tot*vol_est*grid->atom_density/grid->na_c;
if (debug)
{
- fprintf(debug,"nsp_est local %5.1f non-local %5.1f\n",
- nsp_est,nsp_est_nl);
+ fprintf(debug, "nsp_est local %5.1f non-local %5.1f\n",
+ nsp_est, nsp_est_nl);
}
}
else
else
{
/* Thus the (average) maximum j-list size should be as follows */
- nsubpair_max = max(1,(int)(nsp_est/min_ci_balanced+0.5));
-
- /* Since the target value is a maximum (this avoid high outliers,
- * which lead to load imbalance), not average, we get more lists
- * than we ask for (to compensate we need to add GPU_NSUBCELL*4/4).
- * But more importantly, the optimal GPU performance moves
- * to lower number of block for very small blocks.
- * To compensate we add the maximum pair count per cj4.
+ nsubpair_max = max(1, (int)(nsp_est/min_ci_balanced+0.5));
+
+ /* Since the target value is a maximum (this avoids high outliers,
+ * which lead to load imbalance), not average, we add half the
+ * number of pairs in a cj4 block to get the average about right.
*/
- nsubpair_max += GPU_NSUBCELL*NBNXN_CPU_CLUSTER_I_SIZE;
+ nsubpair_max += GPU_NSUBCELL*NBNXN_GPU_JGROUP_SIZE/2;
}
if (debug)
{
- fprintf(debug,"nbl nsp estimate %.1f, nsubpair_max %d\n",
- nsp_est,nsubpair_max);
+ fprintf(debug, "nbl nsp estimate %.1f, nsubpair_max %d\n",
+ nsp_est, nsubpair_max);
}
return nsubpair_max;
}
/* Debug list print function */
-static void print_nblist_ci_cj(FILE *fp,const nbnxn_pairlist_t *nbl)
+static void print_nblist_ci_cj(FILE *fp, const nbnxn_pairlist_t *nbl)
{
- int i,j;
+ int i, j;
- for(i=0; i<nbl->nci; i++)
+ for (i = 0; i < nbl->nci; i++)
{
- fprintf(fp,"ci %4d shift %2d ncj %3d\n",
- nbl->ci[i].ci,nbl->ci[i].shift,
+ fprintf(fp, "ci %4d shift %2d ncj %3d\n",
+ nbl->ci[i].ci, nbl->ci[i].shift,
nbl->ci[i].cj_ind_end - nbl->ci[i].cj_ind_start);
- for(j=nbl->ci[i].cj_ind_start; j<nbl->ci[i].cj_ind_end; j++)
+ for (j = nbl->ci[i].cj_ind_start; j < nbl->ci[i].cj_ind_end; j++)
{
- fprintf(fp," cj %5d imask %x\n",
+ fprintf(fp, " cj %5d imask %x\n",
nbl->cj[j].cj,
nbl->cj[j].excl);
}
}
/* Debug list print function */
-static void print_nblist_sci_cj(FILE *fp,const nbnxn_pairlist_t *nbl)
+static void print_nblist_sci_cj(FILE *fp, const nbnxn_pairlist_t *nbl)
{
- int i,j4,j;
+ int i, j4, j, ncp, si;
- for(i=0; i<nbl->nsci; i++)
+ for (i = 0; i < nbl->nsci; i++)
{
- fprintf(fp,"ci %4d shift %2d ncj4 %2d\n",
- nbl->sci[i].sci,nbl->sci[i].shift,
+ fprintf(fp, "ci %4d shift %2d ncj4 %2d\n",
+ nbl->sci[i].sci, nbl->sci[i].shift,
nbl->sci[i].cj4_ind_end - nbl->sci[i].cj4_ind_start);
- for(j4=nbl->sci[i].cj4_ind_start; j4<nbl->sci[i].cj4_ind_end; j4++)
+ ncp = 0;
+ for (j4 = nbl->sci[i].cj4_ind_start; j4 < nbl->sci[i].cj4_ind_end; j4++)
{
- for(j=0; j<NBNXN_GPU_JGROUP_SIZE; j++)
+ for (j = 0; j < NBNXN_GPU_JGROUP_SIZE; j++)
{
- fprintf(fp," sj %5d imask %x\n",
+ fprintf(fp, " sj %5d imask %x\n",
nbl->cj4[j4].cj[j],
nbl->cj4[j4].imei[0].imask);
+ for (si = 0; si < GPU_NSUBCELL; si++)
+ {
+ if (nbl->cj4[j4].imei[0].imask & (1U << (j*GPU_NSUBCELL + si)))
+ {
+ ncp++;
+ }
+ }
}
}
+ fprintf(fp, "ci %4d shift %2d ncj4 %2d ncp %3d\n",
+ nbl->sci[i].sci, nbl->sci[i].shift,
+ nbl->sci[i].cj4_ind_end - nbl->sci[i].cj4_ind_start,
+ ncp);
}
}
/* Combine pair lists *nbl generated on multiple threads nblc */
-static void combine_nblists(int nnbl,nbnxn_pairlist_t **nbl,
+static void combine_nblists(int nnbl, nbnxn_pairlist_t **nbl,
nbnxn_pairlist_t *nblc)
{
- int nsci,ncj4,nexcl;
- int n,i;
+ int nsci, ncj4, nexcl;
+ int n, i;
if (nblc->bSimple)
{
nsci = nblc->nsci;
ncj4 = nblc->ncj4;
nexcl = nblc->nexcl;
- for(i=0; i<nnbl; i++)
+ for (i = 0; i < nnbl; i++)
{
nsci += nbl[i]->nsci;
ncj4 += nbl[i]->ncj4;
if (nsci > nblc->sci_nalloc)
{
- nb_realloc_sci(nblc,nsci);
+ nb_realloc_sci(nblc, nsci);
}
if (ncj4 > nblc->cj4_nalloc)
{
nbnxn_realloc_void((void **)&nblc->cj4,
nblc->ncj4*sizeof(*nblc->cj4),
nblc->cj4_nalloc*sizeof(*nblc->cj4),
- nblc->alloc,nblc->free);
+ nblc->alloc, nblc->free);
}
if (nexcl > nblc->excl_nalloc)
{
nbnxn_realloc_void((void **)&nblc->excl,
nblc->nexcl*sizeof(*nblc->excl),
nblc->excl_nalloc*sizeof(*nblc->excl),
- nblc->alloc,nblc->free);
+ nblc->alloc, nblc->free);
}
/* Each thread should copy its own data to the combined arrays,
* as otherwise data will go back and forth between different caches.
*/
#pragma omp parallel for num_threads(gmx_omp_nthreads_get(emntPairsearch)) schedule(static)
- for(n=0; n<nnbl; n++)
+ for (n = 0; n < nnbl; n++)
{
- int sci_offset;
- int cj4_offset;
- int ci_offset;
- int excl_offset;
- int i,j4;
+ int sci_offset;
+ int cj4_offset;
+ int ci_offset;
+ int excl_offset;
+ int i, j4;
const nbnxn_pairlist_t *nbli;
/* Determine the offset in the combined data for our thread */
ci_offset = nblc->nci_tot;
excl_offset = nblc->nexcl;
- for(i=0; i<n; i++)
+ for (i = 0; i < n; i++)
{
sci_offset += nbl[i]->nsci;
cj4_offset += nbl[i]->ncj4;
nbli = nbl[n];
- for(i=0; i<nbli->nsci; i++)
+ for (i = 0; i < nbli->nsci; i++)
{
nblc->sci[sci_offset+i] = nbli->sci[i];
nblc->sci[sci_offset+i].cj4_ind_start += cj4_offset;
nblc->sci[sci_offset+i].cj4_ind_end += cj4_offset;
}
- for(j4=0; j4<nbli->ncj4; j4++)
+ for (j4 = 0; j4 < nbli->ncj4; j4++)
{
- nblc->cj4[cj4_offset+j4] = nbli->cj4[j4];
+ nblc->cj4[cj4_offset+j4] = nbli->cj4[j4];
nblc->cj4[cj4_offset+j4].imei[0].excl_ind += excl_offset;
nblc->cj4[cj4_offset+j4].imei[1].excl_ind += excl_offset;
}
- for(j4=0; j4<nbli->nexcl; j4++)
+ for (j4 = 0; j4 < nbli->nexcl; j4++)
{
nblc->excl[excl_offset+j4] = nbli->excl[j4];
}
}
- for(n=0; n<nnbl; n++)
+ for (n = 0; n < nnbl; n++)
{
nblc->nsci += nbl[n]->nsci;
nblc->ncj4 += nbl[n]->ncj4;
}
}
+static void balance_fep_lists(const nbnxn_search_t nbs,
+ nbnxn_pairlist_set_t *nbl_lists)
+{
+ int nnbl, th;
+ int nri_tot, nrj_tot, nrj_target;
+ int th_dest;
+ t_nblist *nbld;
+
+ nnbl = nbl_lists->nnbl;
+
+ if (nnbl == 1)
+ {
+ /* Nothing to balance */
+ return;
+ }
+
+ /* Count the total i-lists and pairs */
+ nri_tot = 0;
+ nrj_tot = 0;
+ for (th = 0; th < nnbl; th++)
+ {
+ nri_tot += nbl_lists->nbl_fep[th]->nri;
+ nrj_tot += nbl_lists->nbl_fep[th]->nrj;
+ }
+
+ nrj_target = (nrj_tot + nnbl - 1)/nnbl;
+
+ assert(gmx_omp_nthreads_get(emntNonbonded) == nnbl);
+
+#pragma omp parallel for schedule(static) num_threads(nnbl)
+ for (th = 0; th < nnbl; th++)
+ {
+ t_nblist *nbl;
+
+ nbl = nbs->work[th].nbl_fep;
+
+ /* Note that here we allocate for the total size, instead of
+ * a per-thread esimate (which is hard to obtain).
+ */
+ if (nri_tot > nbl->maxnri)
+ {
+ nbl->maxnri = over_alloc_large(nri_tot);
+ reallocate_nblist(nbl);
+ }
+ if (nri_tot > nbl->maxnri || nrj_tot > nbl->maxnrj)
+ {
+ nbl->maxnrj = over_alloc_small(nrj_tot);
+ srenew(nbl->jjnr, nbl->maxnrj);
+ srenew(nbl->excl_fep, nbl->maxnrj);
+ }
+
+ clear_pairlist_fep(nbl);
+ }
+
+ /* Loop over the source lists and assign and copy i-entries */
+ th_dest = 0;
+ nbld = nbs->work[th_dest].nbl_fep;
+ for (th = 0; th < nnbl; th++)
+ {
+ t_nblist *nbls;
+ int i, j;
+
+ nbls = nbl_lists->nbl_fep[th];
+
+ for (i = 0; i < nbls->nri; i++)
+ {
+ int nrj;
+
+ /* The number of pairs in this i-entry */
+ nrj = nbls->jindex[i+1] - nbls->jindex[i];
+
+ /* Decide if list th_dest is too large and we should procede
+ * to the next destination list.
+ */
+ if (th_dest+1 < nnbl && nbld->nrj > 0 &&
+ nbld->nrj + nrj - nrj_target > nrj_target - nbld->nrj)
+ {
+ th_dest++;
+ nbld = nbs->work[th_dest].nbl_fep;
+ }
+
+ nbld->iinr[nbld->nri] = nbls->iinr[i];
+ nbld->gid[nbld->nri] = nbls->gid[i];
+ nbld->shift[nbld->nri] = nbls->shift[i];
+
+ for (j = nbls->jindex[i]; j < nbls->jindex[i+1]; j++)
+ {
+ nbld->jjnr[nbld->nrj] = nbls->jjnr[j];
+ nbld->excl_fep[nbld->nrj] = nbls->excl_fep[j];
+ nbld->nrj++;
+ }
+ nbld->nri++;
+ nbld->jindex[nbld->nri] = nbld->nrj;
+ }
+ }
+
+ /* Swap the list pointers */
+ for (th = 0; th < nnbl; th++)
+ {
+ t_nblist *nbl_tmp;
+
+ nbl_tmp = nbl_lists->nbl_fep[th];
+ nbl_lists->nbl_fep[th] = nbs->work[th].nbl_fep;
+ nbs->work[th].nbl_fep = nbl_tmp;
+
+ if (debug)
+ {
+ fprintf(debug, "nbl_fep[%d] nri %4d nrj %4d\n",
+ th,
+ nbl_lists->nbl_fep[th]->nri,
+ nbl_lists->nbl_fep[th]->nrj);
+ }
+ }
+}
+
/* Returns the next ci to be processes by our thread */
static gmx_bool next_ci(const nbnxn_grid_t *grid,
int conv,
- int nth,int ci_block,
- int *ci_x,int *ci_y,
- int *ci_b,int *ci)
+ int nth, int ci_block,
+ int *ci_x, int *ci_y,
+ int *ci_b, int *ci)
{
(*ci_b)++;
(*ci)++;
*/
static float boundingbox_only_distance2(const nbnxn_grid_t *gridi,
const nbnxn_grid_t *gridj,
- real rlist,
- gmx_bool simple)
+ real rlist,
+ gmx_bool simple)
{
/* If the distance between two sub-cell bounding boxes is less
* than this distance, do not check the distance between
* is only performed when only 1 out of 8 sub-cells in within range,
* this is because the GPU is much faster than the cpu.
*/
- real bbx,bby;
+ real bbx, bby;
real rbb2;
bbx = 0.5*(gridi->sx + gridj->sx);
bby /= GPU_NSUBCELL_Y;
}
- rbb2 = sqr(max(0,rlist - 0.5*sqrt(bbx*bbx + bby*bby)));
+ rbb2 = sqr(max(0, rlist - 0.5*sqrt(bbx*bbx + bby*bby)));
#ifndef GMX_DOUBLE
return rbb2;
static int get_ci_block_size(const nbnxn_grid_t *gridi,
gmx_bool bDomDec, int nth)
{
- const int ci_block_enum = 5;
- const int ci_block_denom = 11;
+ const int ci_block_enum = 5;
+ const int ci_block_denom = 11;
const int ci_block_min_atoms = 16;
- int ci_block;
+ int ci_block;
/* Here we decide how to distribute the blocks over the threads.
* We use prime numbers to try to avoid that the grid size becomes
{
ci_block = (ci_block_min_atoms + gridi->na_sc - 1)/gridi->na_sc;
}
-
+
/* Without domain decomposition
* or with less than 3 blocks per task, divide in nth blocks.
*/
int nsubpair_max,
gmx_bool progBal,
int min_ci_balanced,
- int th,int nth,
- nbnxn_pairlist_t *nbl)
+ int th, int nth,
+ nbnxn_pairlist_t *nbl,
+ t_nblist *nbl_fep)
{
- int na_cj_2log;
- matrix box;
- real rl2;
- float rbb2;
- int d;
- int ci_b,ci,ci_x,ci_y,ci_xy,cj;
- ivec shp;
- int tx,ty,tz;
- int shift;
- gmx_bool bMakeList;
- real shx,shy,shz;
- int conv_i,cell0_i;
- const float *bb_i,*bbcz_i,*bbcz_j;
- const int *flags_i;
- real bx0,bx1,by0,by1,bz0,bz1;
- real bz1_frac;
- real d2cx,d2z,d2z_cx,d2z_cy,d2zx,d2zxy,d2xy;
- int cxf,cxl,cyf,cyf_x,cyl;
- int cx,cy;
- int c0,c1,cs,cf,cl;
- int ndistc;
- int ncpcheck;
- int gridi_flag_shift=0,gridj_flag_shift=0;
- unsigned *gridj_flag=NULL;
- int ncj_old_i,ncj_old_j;
+ int na_cj_2log;
+ matrix box;
+ real rl2, rl_fep2 = 0;
+ float rbb2;
+ int d;
+ int ci_b, ci, ci_x, ci_y, ci_xy, cj;
+ ivec shp;
+ int tx, ty, tz;
+ int shift;
+ gmx_bool bMakeList;
+ real shx, shy, shz;
+ int conv_i, cell0_i;
+ const nbnxn_bb_t *bb_i = NULL;
+#ifdef NBNXN_BBXXXX
+ const float *pbb_i = NULL;
+#endif
+ const float *bbcz_i, *bbcz_j;
+ const int *flags_i;
+ real bx0, bx1, by0, by1, bz0, bz1;
+ real bz1_frac;
+ real d2cx, d2z, d2z_cx, d2z_cy, d2zx, d2zxy, d2xy;
+ int cxf, cxl, cyf, cyf_x, cyl;
+ int cx, cy;
+ int c0, c1, cs, cf, cl;
+ int ndistc;
+ int ncpcheck;
+ int gridi_flag_shift = 0, gridj_flag_shift = 0;
+ unsigned int *gridj_flag = NULL;
+ int ncj_old_i, ncj_old_j;
nbs_cycle_start(&work->cc[enbsCCsearch]);
gridj_flag = work->buffer_flags.flag;
}
- copy_mat(nbs->box,box);
+ copy_mat(nbs->box, box);
rl2 = nbl->rlist*nbl->rlist;
- rbb2 = boundingbox_only_distance2(gridi,gridj,nbl->rlist,nbl->bSimple);
+ if (nbs->bFEP && !nbl->bSimple)
+ {
+ /* Determine an atom-pair list cut-off distance for FEP atom pairs.
+ * We should not simply use rlist, since then we would not have
+ * the small, effective buffering of the NxN lists.
+ * The buffer is on overestimate, but the resulting cost for pairs
+ * beyond rlist is neglible compared to the FEP pairs within rlist.
+ */
+ rl_fep2 = nbl->rlist + effective_buffer_1x1_vs_MxN(gridi, gridj);
+
+ if (debug)
+ {
+ fprintf(debug, "nbl_fep atom-pair rlist %f\n", rl_fep2);
+ }
+ rl_fep2 = rl_fep2*rl_fep2;
+ }
+
+ rbb2 = boundingbox_only_distance2(gridi, gridj, nbl->rlist, nbl->bSimple);
if (debug)
{
- fprintf(debug,"nbl bounding box only distance %f\n",sqrt(rbb2));
+ fprintf(debug, "nbl bounding box only distance %f\n", sqrt(rbb2));
}
/* Set the shift range */
- for(d=0; d<DIM; d++)
+ for (d = 0; d < DIM; d++)
{
/* Check if we need periodicity shifts.
* Without PBC or with domain decomposition we don't need them.
else
{
conv_i = 1;
- bb_i = gridi->bb;
+#ifdef NBNXN_BBXXXX
+ if (gridi->bSimple)
+ {
+ bb_i = gridi->bb;
+ }
+ else
+ {
+ pbb_i = gridi->pbb;
+ }
+#else
+ /* We use the normal bounding box format for both grid types */
+ bb_i = gridi->bb;
+#endif
bbcz_i = gridi->bbcz;
flags_i = gridi->flags;
}
}
if (debug)
{
- fprintf(debug,"nbl nc_i %d col.av. %.1f ci_block %d\n",
- gridi->nc,gridi->nc/(double)(gridi->ncx*gridi->ncy),ci_block);
+ fprintf(debug, "nbl nc_i %d col.av. %.1f ci_block %d\n",
+ gridi->nc, gridi->nc/(double)(gridi->ncx*gridi->ncy), ci_block);
}
- ndistc = 0;
+ ndistc = 0;
ncpcheck = 0;
/* Initially ci_b and ci to 1 before where we want them to start,
ci = th*ci_block - 1;
ci_x = 0;
ci_y = 0;
- while (next_ci(gridi,conv_i,nth,ci_block,&ci_x,&ci_y,&ci_b,&ci))
+ while (next_ci(gridi, conv_i, nth, ci_block, &ci_x, &ci_y, &ci_b, &ci))
{
if (nbl->bSimple && flags_i[ci] == 0)
{
{
if (nbl->bSimple)
{
- bx1 = bb_i[ci*NNBSBB_B+NNBSBB_C+XX];
+ bx1 = bb_i[ci].upper[BB_X];
}
else
{
ci_xy = ci_x*gridi->ncy + ci_y;
/* Loop over shift vectors in three dimensions */
- for (tz=-shp[ZZ]; tz<=shp[ZZ]; tz++)
+ for (tz = -shp[ZZ]; tz <= shp[ZZ]; tz++)
{
shz = tz*box[ZZ][ZZ];
}
/* The check with bz1_frac close to or larger than 1 comes later */
- for (ty=-shp[YY]; ty<=shp[YY]; ty++)
+ for (ty = -shp[YY]; ty <= shp[YY]; ty++)
{
shy = ty*box[YY][YY] + tz*box[ZZ][YY];
if (nbl->bSimple)
{
- by0 = bb_i[ci*NNBSBB_B +YY] + shy;
- by1 = bb_i[ci*NNBSBB_B+NNBSBB_C+YY] + shy;
+ by0 = bb_i[ci].lower[BB_Y] + shy;
+ by1 = bb_i[ci].upper[BB_Y] + shy;
}
else
{
by1 = gridi->c0[YY] + (ci_y+1)*gridi->sy + shy;
}
- get_cell_range(by0,by1,
- gridj->ncy,gridj->c0[YY],gridj->sy,gridj->inv_sy,
- d2z_cx,rl2,
- &cyf,&cyl);
+ get_cell_range(by0, by1,
+ gridj->ncy, gridj->c0[YY], gridj->sy, gridj->inv_sy,
+ d2z_cx, rl2,
+ &cyf, &cyl);
if (cyf > cyl)
{
d2z_cy += sqr(by0 - gridj->c1[YY]);
}
- for (tx=-shp[XX]; tx<=shp[XX]; tx++)
+ for (tx = -shp[XX]; tx <= shp[XX]; tx++)
{
- shift = XYZ2IS(tx,ty,tz);
+ shift = XYZ2IS(tx, ty, tz);
#ifdef NBNXN_SHIFT_BACKWARD
if (gridi == gridj && shift > CENTRAL)
if (nbl->bSimple)
{
- bx0 = bb_i[ci*NNBSBB_B +XX] + shx;
- bx1 = bb_i[ci*NNBSBB_B+NNBSBB_C+XX] + shx;
+ bx0 = bb_i[ci].lower[BB_X] + shx;
+ bx1 = bb_i[ci].upper[BB_X] + shx;
}
else
{
bx1 = gridi->c0[XX] + (ci_x+1)*gridi->sx + shx;
}
- get_cell_range(bx0,bx1,
- gridj->ncx,gridj->c0[XX],gridj->sx,gridj->inv_sx,
- d2z_cy,rl2,
- &cxf,&cxl);
+ get_cell_range(bx0, bx1,
+ gridj->ncx, gridj->c0[XX], gridj->sx, gridj->inv_sx,
+ d2z_cy, rl2,
+ &cxf, &cxl);
if (cxf > cxl)
{
if (nbl->bSimple)
{
- new_ci_entry(nbl,cell0_i+ci,shift,flags_i[ci],
- nbl->work);
+ new_ci_entry(nbl, cell0_i+ci, shift, flags_i[ci]);
}
else
{
- new_sci_entry(nbl,cell0_i+ci,shift,flags_i[ci],
- nbl->work);
+ new_sci_entry(nbl, cell0_i+ci, shift);
}
#ifndef NBNXN_SHIFT_BACKWARD
if (nbl->bSimple)
{
- set_icell_bb_simple(bb_i,ci,shx,shy,shz,
+ set_icell_bb_simple(bb_i, ci, shx, shy, shz,
nbl->work->bb_ci);
}
else
{
- set_icell_bb_supersub(bb_i,ci,shx,shy,shz,
+#ifdef NBNXN_BBXXXX
+ set_icell_bbxxxx_supersub(pbb_i, ci, shx, shy, shz,
+ nbl->work->pbb_ci);
+#else
+ set_icell_bb_supersub(bb_i, ci, shx, shy, shz,
nbl->work->bb_ci);
+#endif
}
- nbs->icell_set_x(cell0_i+ci,shx,shy,shz,
- gridi->na_c,nbat->xstride,nbat->x,
+ nbs->icell_set_x(cell0_i+ci, shx, shy, shz,
+ gridi->na_c, nbat->xstride, nbat->x,
nbl->work);
- for(cx=cxf; cx<=cxl; cx++)
+ for (cx = cxf; cx <= cxl; cx++)
{
d2zx = d2z;
if (gridj->c0[XX] + cx*gridj->sx > bx1)
cyf_x = cyf;
}
- for(cy=cyf_x; cy<=cyl; cy++)
+ for (cy = cyf_x; cy <= cyl; cy++)
{
c0 = gridj->cxy_ind[cx*gridj->ncy+cy];
c1 = gridj->cxy_ind[cx*gridj->ncy+cy+1];
* be within range.
*/
cf = cs;
- while(cf > c0 &&
- (bbcz_j[cf*NNBSBB_D+1] >= bz0 ||
- d2xy + sqr(bbcz_j[cf*NNBSBB_D+1] - bz0) < rl2))
+ while (cf > c0 &&
+ (bbcz_j[cf*NNBSBB_D+1] >= bz0 ||
+ d2xy + sqr(bbcz_j[cf*NNBSBB_D+1] - bz0) < rl2))
{
cf--;
}
* be within range.
*/
cl = cs;
- while(cl < c1-1 &&
- (bbcz_j[cl*NNBSBB_D] <= bz1 ||
- d2xy + sqr(bbcz_j[cl*NNBSBB_D] - bz1) < rl2))
+ while (cl < c1-1 &&
+ (bbcz_j[cl*NNBSBB_D] <= bz1 ||
+ d2xy + sqr(bbcz_j[cl*NNBSBB_D] - bz1) < rl2))
{
cl++;
}
int k;
cf = c1;
cl = -1;
- for(k=c0; k<c1; k++)
+ for (k = c0; k < c1; k++)
{
- if (box_dist2(bx0,bx1,by0,by1,bz0,bz1,
- bb+k*NNBSBB_B) < rl2 &&
+ if (box_dist2(bx0, bx1, by0, by1, bz0, bz1, bb+k) < rl2 &&
k < cf)
{
cf = k;
}
- if (box_dist2(bx0,bx1,by0,by1,bz0,bz1,
- bb+k*NNBSBB_B) < rl2 &&
+ if (box_dist2(bx0, bx1, by0, by1, bz0, bz1, bb+k) < rl2 &&
k > cl)
{
cl = k;
* only use cj >= ci.
*/
#ifndef NBNXN_SHIFT_BACKWARD
- cf = max(cf,ci);
+ cf = max(cf, ci);
#else
if (shift == CENTRAL)
{
- cf = max(cf,ci);
+ cf = max(cf, ci);
}
#endif
}
switch (nb_kernel_type)
{
- case nbnxnk4x4_PlainC:
- check_subcell_list_space_simple(nbl,cl-cf+1);
-
- make_cluster_list_simple(gridj,
- nbl,ci,cf,cl,
- (gridi == gridj && shift == CENTRAL),
- nbat->x,
- rl2,rbb2,
- &ndistc);
- break;
+ case nbnxnk4x4_PlainC:
+ check_subcell_list_space_simple(nbl, cl-cf+1);
+
+ make_cluster_list_simple(gridj,
+ nbl, ci, cf, cl,
+ (gridi == gridj && shift == CENTRAL),
+ nbat->x,
+ rl2, rbb2,
+ &ndistc);
+ break;
#ifdef GMX_NBNXN_SIMD_4XN
- case nbnxnk4xN_SIMD_4xN:
- check_subcell_list_space_simple(nbl,ci_to_cj(na_cj_2log,cl-cf)+2);
- make_cluster_list_simd_4xn(gridj,
- nbl,ci,cf,cl,
- (gridi == gridj && shift == CENTRAL),
- nbat->x,
- rl2,rbb2,
- &ndistc);
- break;
+ case nbnxnk4xN_SIMD_4xN:
+ check_subcell_list_space_simple(nbl, ci_to_cj(na_cj_2log, cl-cf)+2);
+ make_cluster_list_simd_4xn(gridj,
+ nbl, ci, cf, cl,
+ (gridi == gridj && shift == CENTRAL),
+ nbat->x,
+ rl2, rbb2,
+ &ndistc);
+ break;
#endif
#ifdef GMX_NBNXN_SIMD_2XNN
- case nbnxnk4xN_SIMD_2xNN:
- check_subcell_list_space_simple(nbl,ci_to_cj(na_cj_2log,cl-cf)+2);
- make_cluster_list_simd_2xnn(gridj,
- nbl,ci,cf,cl,
- (gridi == gridj && shift == CENTRAL),
- nbat->x,
- rl2,rbb2,
- &ndistc);
- break;
+ case nbnxnk4xN_SIMD_2xNN:
+ check_subcell_list_space_simple(nbl, ci_to_cj(na_cj_2log, cl-cf)+2);
+ make_cluster_list_simd_2xnn(gridj,
+ nbl, ci, cf, cl,
+ (gridi == gridj && shift == CENTRAL),
+ nbat->x,
+ rl2, rbb2,
+ &ndistc);
+ break;
#endif
- case nbnxnk8x8x8_PlainC:
- case nbnxnk8x8x8_CUDA:
- check_subcell_list_space_supersub(nbl,cl-cf+1);
- for(cj=cf; cj<=cl; cj++)
- {
- make_cluster_list_supersub(nbs,gridi,gridj,
- nbl,ci,cj,
- (gridi == gridj && shift == CENTRAL && ci == cj),
- nbat->xstride,nbat->x,
- rl2,rbb2,
- &ndistc);
- }
- break;
+ case nbnxnk8x8x8_PlainC:
+ case nbnxnk8x8x8_CUDA:
+ check_subcell_list_space_supersub(nbl, cl-cf+1);
+ for (cj = cf; cj <= cl; cj++)
+ {
+ make_cluster_list_supersub(gridi, gridj,
+ nbl, ci, cj,
+ (gridi == gridj && shift == CENTRAL && ci == cj),
+ nbat->xstride, nbat->x,
+ rl2, rbb2,
+ &ndistc);
+ }
+ break;
}
ncpcheck += cl - cf + 1;
if (bFBufferFlag && nbl->ncj > ncj_old_j)
{
- int cbf,cbl,cb;
+ int cbf, cbl, cb;
cbf = nbl->cj[ncj_old_j].cj >> gridj_flag_shift;
cbl = nbl->cj[nbl->ncj-1].cj >> gridj_flag_shift;
- for(cb=cbf; cb<=cbl; cb++)
+ for (cb = cbf; cb <= cbl; cb++)
{
gridj_flag[cb] = 1U<<th;
}
na_cj_2log,
&(nbl->ci[nbl->nci]),
excl);
+
+ if (nbs->bFEP)
+ {
+ make_fep_list(nbs, nbat, nbl,
+ shift == CENTRAL && gridi == gridj,
+ &(nbl->ci[nbl->nci]),
+ gridi, gridj, nbl_fep);
+ }
}
else
{
gridj->na_c_2log,
&(nbl->sci[nbl->nsci]),
excl);
+
+ if (nbs->bFEP)
+ {
+ make_fep_list_supersub(nbs, nbat, nbl,
+ shift == CENTRAL && gridi == gridj,
+ &(nbl->sci[nbl->nsci]),
+ shx, shy, shz,
+ rl_fep2,
+ gridi, gridj, nbl_fep);
+ }
}
/* Close this ci list */
{
close_ci_entry_supersub(nbl,
nsubpair_max,
- progBal,min_ci_balanced,
- th,nth);
+ progBal, min_ci_balanced,
+ th, nth);
}
}
}
if (debug)
{
- fprintf(debug,"number of distance checks %d\n",ndistc);
- fprintf(debug,"ncpcheck %s %d\n",gridi==gridj ? "local" : "non-local",
+ fprintf(debug, "number of distance checks %d\n", ndistc);
+ fprintf(debug, "ncpcheck %s %d\n", gridi == gridj ? "local" : "non-local",
ncpcheck);
if (nbl->bSimple)
{
- print_nblist_statistics_simple(debug,nbl,nbs,rlist);
+ print_nblist_statistics_simple(debug, nbl, nbs, rlist);
}
else
{
- print_nblist_statistics_supersub(debug,nbl,nbs,rlist);
+ print_nblist_statistics_supersub(debug, nbl, nbs, rlist);
}
+ if (nbs->bFEP)
+ {
+ fprintf(debug, "nbl FEP list pairs: %d\n", nbl_fep->nrj);
+ }
}
}
-static void reduce_buffer_flags(const nbnxn_search_t nbs,
- int nsrc,
+static void reduce_buffer_flags(const nbnxn_search_t nbs,
+ int nsrc,
const nbnxn_buffer_flags_t *dest)
{
- int s,b;
- const unsigned *flag;
+ int s, b;
+ const unsigned int *flag;
- for(s=0; s<nsrc; s++)
+ for (s = 0; s < nsrc; s++)
{
flag = nbs->work[s].buffer_flags.flag;
- for(b=0; b<dest->nflag; b++)
+ for (b = 0; b < dest->nflag; b++)
{
dest->flag[b] |= flag[b];
}
}
}
-static void print_reduction_cost(const nbnxn_buffer_flags_t *flags,int nout)
+static void print_reduction_cost(const nbnxn_buffer_flags_t *flags, int nout)
{
- int nelem,nkeep,ncopy,nred,b,c,out;
+ int nelem, nkeep, ncopy, nred, b, c, out;
nelem = 0;
nkeep = 0;
ncopy = 0;
nred = 0;
- for(b=0; b<flags->nflag; b++)
+ for (b = 0; b < flags->nflag; b++)
{
if (flags->flag[b] == 1)
{
else if (flags->flag[b] > 0)
{
c = 0;
- for(out=0; out<nout; out++)
+ for (out = 0; out < nout; out++)
{
if (flags->flag[b] & (1U<<out))
{
}
}
- fprintf(debug,"nbnxn reduction: #flag %d #list %d elem %4.2f, keep %4.2f copy %4.2f red %4.2f\n",
- flags->nflag,nout,
+ fprintf(debug, "nbnxn reduction: #flag %d #list %d elem %4.2f, keep %4.2f copy %4.2f red %4.2f\n",
+ flags->nflag, nout,
nelem/(double)(flags->nflag),
nkeep/(double)(flags->nflag),
ncopy/(double)(flags->nflag),
nred/(double)(flags->nflag));
}
+/* Perform a count (linear) sort to sort the smaller lists to the end.
+ * This avoids load imbalance on the GPU, as large lists will be
+ * scheduled and executed first and the smaller lists later.
+ * Load balancing between multi-processors only happens at the end
+ * and there smaller lists lead to more effective load balancing.
+ * The sorting is done on the cj4 count, not on the actual pair counts.
+ * Not only does this make the sort faster, but it also results in
+ * better load balancing than using a list sorted on exact load.
+ * This function swaps the pointer in the pair list to avoid a copy operation.
+ */
+static void sort_sci(nbnxn_pairlist_t *nbl)
+{
+ nbnxn_list_work_t *work;
+ int m, i, s, s0, s1;
+ nbnxn_sci_t *sci_sort;
+
+ if (nbl->ncj4 <= nbl->nsci)
+ {
+ /* nsci = 0 or all sci have size 1, sorting won't change the order */
+ return;
+ }
+
+ work = nbl->work;
+
+ /* We will distinguish differences up to double the average */
+ m = (2*nbl->ncj4)/nbl->nsci;
+
+ if (m + 1 > work->sort_nalloc)
+ {
+ work->sort_nalloc = over_alloc_large(m + 1);
+ srenew(work->sort, work->sort_nalloc);
+ }
+
+ if (work->sci_sort_nalloc != nbl->sci_nalloc)
+ {
+ work->sci_sort_nalloc = nbl->sci_nalloc;
+ nbnxn_realloc_void((void **)&work->sci_sort,
+ 0,
+ work->sci_sort_nalloc*sizeof(*work->sci_sort),
+ nbl->alloc, nbl->free);
+ }
+
+ /* Count the entries of each size */
+ for (i = 0; i <= m; i++)
+ {
+ work->sort[i] = 0;
+ }
+ for (s = 0; s < nbl->nsci; s++)
+ {
+ i = min(m, nbl->sci[s].cj4_ind_end - nbl->sci[s].cj4_ind_start);
+ work->sort[i]++;
+ }
+ /* Calculate the offset for each count */
+ s0 = work->sort[m];
+ work->sort[m] = 0;
+ for (i = m - 1; i >= 0; i--)
+ {
+ s1 = work->sort[i];
+ work->sort[i] = work->sort[i + 1] + s0;
+ s0 = s1;
+ }
+
+ /* Sort entries directly into place */
+ sci_sort = work->sci_sort;
+ for (s = 0; s < nbl->nsci; s++)
+ {
+ i = min(m, nbl->sci[s].cj4_ind_end - nbl->sci[s].cj4_ind_start);
+ sci_sort[work->sort[i]++] = nbl->sci[s];
+ }
+
+ /* Swap the sci pointers so we use the new, sorted list */
+ work->sci_sort = nbl->sci;
+ nbl->sci = sci_sort;
+}
+
/* Make a local or non-local pair-list, depending on iloc */
-void nbnxn_make_pairlist(const nbnxn_search_t nbs,
- nbnxn_atomdata_t *nbat,
- const t_blocka *excl,
- real rlist,
- int min_ci_balanced,
+void nbnxn_make_pairlist(const nbnxn_search_t nbs,
+ nbnxn_atomdata_t *nbat,
+ const t_blocka *excl,
+ real rlist,
+ int min_ci_balanced,
nbnxn_pairlist_set_t *nbl_list,
- int iloc,
- int nb_kernel_type,
- t_nrnb *nrnb)
+ int iloc,
+ int nb_kernel_type,
+ t_nrnb *nrnb)
{
- nbnxn_grid_t *gridi,*gridj;
- gmx_bool bGPUCPU;
- int nzi,zi,zj0,zj1,zj;
- int nsubpair_max;
- int th;
- int nnbl;
+ nbnxn_grid_t *gridi, *gridj;
+ gmx_bool bGPUCPU;
+ int nzi, zi, zj0, zj1, zj;
+ int nsubpair_max;
+ int th;
+ int nnbl;
nbnxn_pairlist_t **nbl;
- int ci_block;
- gmx_bool CombineNBLists;
- int np_tot,np_noq,np_hlj,nap;
+ int ci_block;
+ gmx_bool CombineNBLists;
+ gmx_bool progBal;
+ int np_tot, np_noq, np_hlj, nap;
/* Check if we are running hybrid GPU + CPU nbnxn mode */
bGPUCPU = (!nbs->grid[0].bSimple && nbl_list->bSimple);
if (debug)
{
- fprintf(debug,"ns making %d nblists\n", nnbl);
+ fprintf(debug, "ns making %d nblists\n", nnbl);
}
nbat->bUseBufferFlags = (nbat->nout > 1);
/* We should re-init the flags before making the first list */
if (nbat->bUseBufferFlags && (LOCAL_I(iloc) || bGPUCPU))
{
- init_buffer_flags(&nbat->buffer_flags,nbat->natoms);
+ init_buffer_flags(&nbat->buffer_flags, nbat->natoms);
}
if (nbl_list->bSimple)
switch (nb_kernel_type)
{
#ifdef GMX_NBNXN_SIMD_4XN
- case nbnxnk4xN_SIMD_4xN:
- nbs->icell_set_x = icell_set_x_simd_4xn;
- break;
+ case nbnxnk4xN_SIMD_4xN:
+ nbs->icell_set_x = icell_set_x_simd_4xn;
+ break;
#endif
#ifdef GMX_NBNXN_SIMD_2XNN
- case nbnxnk4xN_SIMD_2xNN:
- nbs->icell_set_x = icell_set_x_simd_2xnn;
- break;
+ case nbnxnk4xN_SIMD_2xNN:
+ nbs->icell_set_x = icell_set_x_simd_2xnn;
+ break;
#endif
- default:
- nbs->icell_set_x = icell_set_x_simple;
- break;
+ default:
+ nbs->icell_set_x = icell_set_x_simple;
+ break;
}
}
else
{
-#ifdef NBNXN_SEARCH_BB_SSE
- nbs->icell_set_x = icell_set_x_supersub_sse8;
+#ifdef NBNXN_SEARCH_BB_SIMD4
+ nbs->icell_set_x = icell_set_x_supersub_simd4;
#else
nbs->icell_set_x = icell_set_x_supersub;
#endif
if (!nbl_list->bSimple && min_ci_balanced > 0)
{
- nsubpair_max = get_nsubpair_max(nbs,iloc,rlist,min_ci_balanced);
+ nsubpair_max = get_nsubpair_max(nbs, iloc, rlist, min_ci_balanced);
}
else
{
}
/* Clear all pair-lists */
- for(th=0; th<nnbl; th++)
+ for (th = 0; th < nnbl; th++)
{
clear_pairlist(nbl[th]);
+
+ if (nbs->bFEP)
+ {
+ clear_pairlist_fep(nbl_list->nbl_fep[th]);
+ }
}
- for(zi=0; zi<nzi; zi++)
+ for (zi = 0; zi < nzi; zi++)
{
gridi = &nbs->grid[zi];
zj0++;
}
}
- for(zj=zj0; zj<zj1; zj++)
+ for (zj = zj0; zj < zj1; zj++)
{
gridj = &nbs->grid[zj];
if (debug)
{
- fprintf(debug,"ns search grid %d vs %d\n",zi,zj);
+ fprintf(debug, "ns search grid %d vs %d\n", zi, zj);
}
nbs_cycle_start(&nbs->cc[enbsCCsearch]);
}
else
{
- ci_block = get_ci_block_size(gridi,nbs->DomDec,nnbl);
+ ci_block = get_ci_block_size(gridi, nbs->DomDec, nnbl);
}
+ /* With GPU: generate progressively smaller lists for
+ * load balancing for local only or non-local with 2 zones.
+ */
+ progBal = (LOCAL_I(iloc) || nbs->zones->n <= 2);
+
#pragma omp parallel for num_threads(nnbl) schedule(static)
- for(th=0; th<nnbl; th++)
+ for (th = 0; th < nnbl; th++)
{
/* Re-init the thread-local work flag data before making
* the first list (not an elegant conditional).
if (nbat->bUseBufferFlags && ((zi == 0 && zj == 0) ||
(bGPUCPU && zi == 0 && zj == 1)))
{
- init_buffer_flags(&nbs->work[th].buffer_flags,nbat->natoms);
+ init_buffer_flags(&nbs->work[th].buffer_flags, nbat->natoms);
}
if (CombineNBLists && th > 0)
}
/* Divide the i super cell equally over the nblists */
- nbnxn_make_pairlist_part(nbs,gridi,gridj,
- &nbs->work[th],nbat,excl,
+ nbnxn_make_pairlist_part(nbs, gridi, gridj,
+ &nbs->work[th], nbat, excl,
rlist,
nb_kernel_type,
ci_block,
nbat->bUseBufferFlags,
nsubpair_max,
- (LOCAL_I(iloc) || nbs->zones->n <= 2),
- min_ci_balanced,
- th,nnbl,
- nbl[th]);
+ progBal, min_ci_balanced,
+ th, nnbl,
+ nbl[th],
+ nbl_list->nbl_fep[th]);
}
nbs_cycle_stop(&nbs->cc[enbsCCsearch]);
np_tot = 0;
np_noq = 0;
np_hlj = 0;
- for(th=0; th<nnbl; th++)
+ for (th = 0; th < nnbl; th++)
{
- inc_nrnb(nrnb,eNR_NBNXN_DIST2,nbs->work[th].ndistc);
+ inc_nrnb(nrnb, eNR_NBNXN_DIST2, nbs->work[th].ndistc);
if (nbl_list->bSimple)
{
np_tot += nbl[th]->nci_tot;
}
}
- nap = nbl[0]->na_ci*nbl[0]->na_cj;
+ nap = nbl[0]->na_ci*nbl[0]->na_cj;
nbl_list->natpair_ljq = (np_tot - np_noq)*nap - np_hlj*nap/2;
nbl_list->natpair_lj = np_noq*nap;
nbl_list->natpair_q = np_hlj*nap/2;
{
nbs_cycle_start(&nbs->cc[enbsCCcombine]);
- combine_nblists(nnbl-1,nbl+1,nbl[0]);
+ combine_nblists(nnbl-1, nbl+1, nbl[0]);
nbs_cycle_stop(&nbs->cc[enbsCCcombine]);
}
}
}
+ if (!nbl_list->bSimple)
+ {
+ /* Sort the entries on size, large ones first */
+ if (CombineNBLists || nnbl == 1)
+ {
+ sort_sci(nbl[0]);
+ }
+ else
+ {
+#pragma omp parallel for num_threads(nnbl) schedule(static)
+ for (th = 0; th < nnbl; th++)
+ {
+ sort_sci(nbl[th]);
+ }
+ }
+ }
+
if (nbat->bUseBufferFlags)
{
- reduce_buffer_flags(nbs,nnbl,&nbat->buffer_flags);
+ reduce_buffer_flags(nbs, nnbl, &nbat->buffer_flags);
}
- /*
- print_supersub_nsp("nsubpair",nbl[0],iloc);
- */
+ if (nbs->bFEP)
+ {
+ /* Balance the free-energy lists over all the threads */
+ balance_fep_lists(nbs, nbl_list);
+ }
/* Special performance logging stuff (env.var. GMX_NBNXN_CYCLE) */
if (LOCAL_I(iloc))
(!nbs->DomDec || (nbs->DomDec && !LOCAL_I(iloc))) &&
nbs->search_count % 100 == 0)
{
- nbs_cycle_print(stderr,nbs);
+ nbs_cycle_print(stderr, nbs);
}
if (debug && (CombineNBLists && nnbl > 1))
{
if (nbl[0]->bSimple)
{
- print_nblist_statistics_simple(debug,nbl[0],nbs,rlist);
+ print_nblist_statistics_simple(debug, nbl[0], nbs, rlist);
}
else
{
- print_nblist_statistics_supersub(debug,nbl[0],nbs,rlist);
+ print_nblist_statistics_supersub(debug, nbl[0], nbs, rlist);
}
}
{
if (nbl[0]->bSimple)
{
- print_nblist_ci_cj(debug,nbl[0]);
+ print_nblist_ci_cj(debug, nbl[0]);
}
else
{
- print_nblist_sci_cj(debug,nbl[0]);
+ print_nblist_sci_cj(debug, nbl[0]);
}
}
if (nbat->bUseBufferFlags)
{
- print_reduction_cost(&nbat->buffer_flags,nnbl);
+ print_reduction_cost(&nbat->buffer_flags, nnbl);
}
}
}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff