2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
37 * \brief Implements PME OpenCL spline parameter computation and charge spread kernels.
38 * When including this and other PME OpenCL kernel files, plenty of common
39 * constants/macros are expected to be defined (such as "order" which is PME interpolation order).
40 * For details, please see how pme_program.cl is compiled in pme_gpu_program_impl_ocl.cpp.
42 * This file's kernels specifically expect the following definitions:
44 * - atomsPerBlock which expresses how many atoms are processed by a single work group
45 * - order which is a PME interpolation order
46 * - computeSplines and spreadCharges must evaluate to either true or false to specify which
47 * kernel falvor is being compiled
48 * - wrapX and wrapY must evaluate to either true or false to specify whether the grid overlap
49 * in dimension X/Y is to be used
51 * \author Aleksei Iupinov <a.yupinov@gmail.com>
54 #include "gromacs/gpu_utils/vectype_ops.clh"
56 #include "pme_gpu_types.h"
57 #include "pme_gpu_calculate_splines.clh"
60 * This define affects the spline calculation behaviour in the kernel.
61 * 0: a single GPU thread handles a single dimension of a single particle (calculating and storing
62 * (order) spline values and derivatives). 1: (order) threads do redundant work on this same task,
63 * each one stores only a single theta and single dtheta into global arrays. The only efficiency
64 * difference is less global store operations, countered by more redundant spline computation.
66 * TODO: estimate if this should be a boolean parameter (and add it to the unit test if so).
68 #define PME_GPU_PARALLEL_SPLINE 0
70 #ifndef COMPILE_SPREAD_HELPERS_ONCE
71 # define COMPILE_SPREAD_HELPERS_ONCE
74 * General purpose function for loading atom-related data from global to shared memory.
76 * \param[out] sm_destination Local memory array for output.
77 * \param[in] gm_source Global memory array for input.
78 * \param[in] dataCountPerAtom Number of data elements per single atom (e.g. DIM for an rvec coordinates array).
81 inline void pme_gpu_stage_atom_data(__local float* __restrict__ sm_destination,
82 __global const float* __restrict__ gm_source,
83 const int dataCountPerAtom)
85 assert((get_local_id(2) * get_local_size(1) + get_local_id(1)) * get_local_size(0) + get_local_id(0)
87 const int threadLocalIndex =
88 (int)((get_local_id(2) * get_local_size(1) + get_local_id(1)) * get_local_size(0)
90 const int localIndex = threadLocalIndex;
91 const int globalIndexBase = (int)get_group_id(XX) * atomsPerBlock * dataCountPerAtom;
92 const int globalIndex = globalIndexBase + localIndex;
93 if (localIndex < atomsPerBlock * dataCountPerAtom)
95 assert(isfinite(float(gm_source[globalIndex])));
96 sm_destination[localIndex] = gm_source[globalIndex];
101 * PME GPU spline parameter and gridline indices calculation.
102 * This corresponds to the CPU functions calc_interpolation_idx() and make_bsplines().
103 * First stage of the whole kernel.
105 * \param[in] kernelParams Input PME GPU data in constant memory.
106 * \param[in] atomIndexOffset Starting atom index for the execution block w.r.t. global memory.
107 * \param[in] sm_coordinates Atom coordinates in the shared memory (rvec)
108 * \param[in] sm_coefficients Atom charges/coefficients in the shared memory.
109 * \param[out] sm_theta Atom spline values in the shared memory.
110 * \param[out] sm_gridlineIndices Atom gridline indices in the shared memory.
111 * \param[out] sm_fractCoords Atom fractional coordinates (rvec)
112 * \param[out] gm_theta Atom spline parameter values
113 * \param[out] gm_dtheta Atom spline parameter derivatives
114 * \param[out] gm_gridlineIndices Same as \p sm_gridlineIndices but in global memory.
115 * \param[in] gm_fractShiftsTable Atom fractional coordinates correction table
116 * \param[in] gm_gridlineIndicesTable Atom fractional coordinates correction table
118 gmx_opencl_inline void calculate_splines(const struct PmeOpenCLKernelParams kernelParams,
119 const int atomIndexOffset,
120 __local const float* __restrict__ sm_coordinates,
121 __local const float* __restrict__ sm_coefficients,
122 __local float* __restrict__ sm_theta,
123 __local int* __restrict__ sm_gridlineIndices,
124 __local float* __restrict__ sm_fractCoords,
125 __global float* __restrict__ gm_theta,
126 __global float* __restrict__ gm_dtheta,
127 __global int* __restrict__ gm_gridlineIndices,
128 __global const float* __restrict__ gm_fractShiftsTable,
129 __global const int* __restrict__ gm_gridlineIndicesTable)
131 /* Thread index w.r.t. block */
132 assert((get_local_id(2) * get_local_size(1) + get_local_id(1)) * get_local_size(0) + get_local_id(0)
134 const int threadLocalIndex =
135 (int)((get_local_id(2) * get_local_size(1) + get_local_id(1)) * get_local_size(0)
137 /* Warp index w.r.t. block - could probably be obtained easier? */
138 const int warpIndex = threadLocalIndex / warp_size;
139 /* Thread index w.r.t. warp */
140 const int threadWarpIndex = threadLocalIndex % warp_size;
141 /* Atom index w.r.t. warp - alternating 0 1 0 1 .. */
142 const int atomWarpIndex = threadWarpIndex % atomsPerWarp;
143 /* Atom index w.r.t. block/shared memory */
144 const int atomIndexLocal = warpIndex * atomsPerWarp + atomWarpIndex;
146 /* Spline contribution index in one dimension */
147 const int orderIndex = threadWarpIndex / (atomsPerWarp * DIM);
148 /* Dimension index */
149 const int dimIndex = (threadWarpIndex - orderIndex * (atomsPerWarp * DIM)) / atomsPerWarp;
151 /* Multi-purpose index of rvec/ivec atom data */
152 const int sharedMemoryIndex = atomIndexLocal * DIM + dimIndex;
154 /* Spline parameters need a working buffer.
155 * With PME_GPU_PARALLEL_SPLINE == 0 it is just a local array of (order) values for some of the threads, which is fine;
156 * With PME_GPU_PARALLEL_SPLINE == 1 (order) times more threads are involved, so the shared memory is used to avoid overhead.
157 * The buffer's size, striding and indexing are adjusted accordingly.
158 * The buffer is accessed with SPLINE_DATA_PTR and SPLINE_DATA macros.
160 # if PME_GPU_PARALLEL_SPLINE
161 # define splineDataStride (atomsPerBlock * DIM)
162 const int splineDataIndex = sharedMemoryIndex;
163 __local float sm_splineData[splineDataStride * order];
164 __local float* splineDataPtr = sm_splineData;
166 # define splineDataStride 1
167 const int splineDataIndex = 0;
168 float splineData[splineDataStride * order];
169 float* splineDataPtr = splineData;
172 # define SPLINE_DATA_PTR(i) (splineDataPtr + ((i)*splineDataStride + splineDataIndex))
173 # define SPLINE_DATA(i) (*SPLINE_DATA_PTR(i))
175 const int localCheck = (dimIndex < DIM) && (orderIndex < (PME_GPU_PARALLEL_SPLINE ? order : 1));
179 /* Indices interpolation */
186 const float3 x = vload3(atomIndexLocal, sm_coordinates);
188 /* Accessing fields in fshOffset/nXYZ/recipbox/... with dimIndex offset
189 * puts them into local memory(!) instead of accessing the constant memory directly.
190 * That's the reason for the switch, to unroll explicitly.
191 * The commented parts correspond to the 0 components of the recipbox.
196 tableIndex = kernelParams.grid.tablesOffsets[XX];
197 n = kernelParams.grid.realGridSizeFP[XX];
198 t = x.x * kernelParams.current.recipBox[dimIndex][XX]
199 + x.y * kernelParams.current.recipBox[dimIndex][YY]
200 + x.z * kernelParams.current.recipBox[dimIndex][ZZ];
204 tableIndex = kernelParams.grid.tablesOffsets[YY];
205 n = kernelParams.grid.realGridSizeFP[YY];
206 t = /*x.x * kernelParams.current.recipbox[dimIndex][XX] + */ x.y
207 * kernelParams.current.recipBox[dimIndex][YY]
208 + x.z * kernelParams.current.recipBox[dimIndex][ZZ];
212 tableIndex = kernelParams.grid.tablesOffsets[ZZ];
213 n = kernelParams.grid.realGridSizeFP[ZZ];
214 t = /*x.x * kernelParams.current.recipbox[dimIndex][XX] + x.y * kernelParams.current.recipbox[dimIndex][YY] + */ x
216 * kernelParams.current.recipBox[dimIndex][ZZ];
223 const float shift = c_pmeMaxUnitcellShift;
225 /* Fractional coordinates along box vectors, adding a positive shift to ensure t is positive for triclinic boxes */
228 sm_fractCoords[sharedMemoryIndex] = t - (float)tInt;
231 assert(tInt < c_pmeNeighborUnitcellCount * n);
232 sm_fractCoords[sharedMemoryIndex] += gm_fractShiftsTable[tableIndex];
233 sm_gridlineIndices[sharedMemoryIndex] = gm_gridlineIndicesTable[tableIndex];
234 gm_gridlineIndices[atomIndexOffset * DIM + sharedMemoryIndex] =
235 sm_gridlineIndices[sharedMemoryIndex];
238 /* B-spline calculation */
240 const int chargeCheck = pme_gpu_check_atom_charge(sm_coefficients[atomIndexLocal]);
244 int o = orderIndex; // This is an index that is set once for PME_GPU_PARALLEL_SPLINE == 1
246 const float dr = sm_fractCoords[sharedMemoryIndex];
247 assert(isfinite(dr));
249 /* dr is relative offset from lower cell limit */
250 *SPLINE_DATA_PTR(order - 1) = 0.0F;
251 *SPLINE_DATA_PTR(1) = dr;
252 *SPLINE_DATA_PTR(0) = 1.0F - dr;
254 # pragma unroll order
255 for (int k = 3; k < order; k++)
257 div = 1.0F / ((float)k - 1.0F);
258 *SPLINE_DATA_PTR(k - 1) = div * dr * SPLINE_DATA(k - 2);
260 for (int l = 1; l < (k - 1); l++)
262 *SPLINE_DATA_PTR(k - l - 1) =
264 * ((dr + (float)l) * SPLINE_DATA(k - l - 2)
265 + ((float)k - (float)l - dr) * SPLINE_DATA(k - l - 1));
267 *SPLINE_DATA_PTR(0) = div * (1.0F - dr) * SPLINE_DATA(0);
270 const int thetaIndexBase = getSplineParamIndexBase(warpIndex, atomWarpIndex);
271 const int thetaGlobalOffsetBase = atomIndexOffset * DIM * order;
273 /* Differentiation and storing the spline derivatives (dtheta) */
274 # if !PME_GPU_PARALLEL_SPLINE
275 // With PME_GPU_PARALLEL_SPLINE == 1, o is already set to orderIndex;
276 // With PME_GPU_PARALLEL_SPLINE == 0, we loop o over range(order).
277 # pragma unroll order
278 for (o = 0; o < order; o++)
281 const int thetaIndex = getSplineParamIndex(thetaIndexBase, dimIndex, o);
282 const int thetaGlobalIndex = thetaGlobalOffsetBase + thetaIndex;
284 const float dtheta = ((o > 0) ? SPLINE_DATA(o - 1) : 0.0F) - SPLINE_DATA(o);
285 assert(isfinite(dtheta));
286 gm_dtheta[thetaGlobalIndex] = dtheta;
289 div = 1.0F / (order - 1.0F);
290 *SPLINE_DATA_PTR(order - 1) = div * dr * SPLINE_DATA(order - 2);
292 for (int k = 1; k < (order - 1); k++)
294 *SPLINE_DATA_PTR(order - k - 1) = div
295 * ((dr + (float)k) * SPLINE_DATA(order - k - 2)
296 + (order - k - dr) * SPLINE_DATA(order - k - 1));
298 *SPLINE_DATA_PTR(0) = div * (1.0F - dr) * SPLINE_DATA(0);
300 /* Storing the spline values (theta) */
301 # if !PME_GPU_PARALLEL_SPLINE
302 // See comment for the loop above
303 # pragma unroll order
304 for (o = 0; o < order; o++)
307 const int thetaIndex = getSplineParamIndex(thetaIndexBase, dimIndex, o);
308 const int thetaGlobalIndex = thetaGlobalOffsetBase + thetaIndex;
310 sm_theta[thetaIndex] = SPLINE_DATA(o);
311 assert(isfinite(sm_theta[thetaIndex]));
312 gm_theta[thetaGlobalIndex] = SPLINE_DATA(o);
319 * Charge spreading onto the grid.
320 * This corresponds to the CPU function spread_coefficients_bsplines_thread().
321 * Second stage of the whole kernel.
323 * \param[in] kernelParams Input PME GPU data in constant memory.
324 * \param[in] atomIndexOffset Starting atom index for the execution block w.r.t. global
325 * memory. \param[in] sm_coefficients Atom coefficents/charges in the shared memory. \param[in]
326 * sm_gridlineIndices Atom gridline indices in the shared memory. \param[in] sm_theta Atom spline
327 * values in the shared memory. \param[out] gm_grid Global 3D grid for spreading.
329 gmx_opencl_inline void spread_charges(const struct PmeOpenCLKernelParams kernelParams,
331 __local const float* __restrict__ sm_coefficients,
332 __local const int* __restrict__ sm_gridlineIndices,
333 __local const float* __restrict__ sm_theta,
334 __global float* __restrict__ gm_grid)
336 const int nx = kernelParams.grid.realGridSize[XX];
337 const int ny = kernelParams.grid.realGridSize[YY];
338 const int nz = kernelParams.grid.realGridSize[ZZ];
339 const int pny = kernelParams.grid.realGridSizePadded[YY];
340 const int pnz = kernelParams.grid.realGridSizePadded[ZZ];
346 const int atomIndexLocal = get_local_id(ZZ);
348 const int chargeCheck = pme_gpu_check_atom_charge(sm_coefficients[atomIndexLocal]);
351 // Spline Y/Z coordinates
352 const int ithy = get_local_id(YY);
353 const int ithz = get_local_id(XX);
354 const int ixBase = sm_gridlineIndices[atomIndexLocal * DIM + XX] - offx;
355 int iy = sm_gridlineIndices[atomIndexLocal * DIM + YY] - offy + ithy;
356 if (wrapY & (iy >= ny))
360 int iz = sm_gridlineIndices[atomIndexLocal * DIM + ZZ] - offz + ithz;
366 /* Atom index w.r.t. warp - alternating 0 1 0 1 .. */
367 const int atomWarpIndex = atomIndexLocal % atomsPerWarp;
368 /* Warp index w.r.t. block - could probably be obtained easier? */
369 const int warpIndex = atomIndexLocal / atomsPerWarp;
371 const int splineIndexBase = getSplineParamIndexBase(warpIndex, atomWarpIndex);
372 const int splineIndexZ = getSplineParamIndex(splineIndexBase, ZZ, ithz);
373 const float thetaZ = sm_theta[splineIndexZ];
374 const int splineIndexY = getSplineParamIndex(splineIndexBase, YY, ithy);
375 const float thetaY = sm_theta[splineIndexY];
376 const float constVal = thetaZ * thetaY * sm_coefficients[atomIndexLocal];
377 assert(isfinite(constVal));
378 const int constOffset = iy * pnz + iz;
380 # pragma unroll order
381 for (int ithx = 0; (ithx < order); ithx++)
383 int ix = ixBase + ithx;
384 if (wrapX & (ix >= nx))
388 const int gridIndexGlobal = ix * pny * pnz + constOffset;
389 const int splineIndexX = getSplineParamIndex(splineIndexBase, XX, ithx);
390 const float thetaX = sm_theta[splineIndexX];
391 assert(isfinite(thetaX));
392 assert(isfinite(gm_grid[gridIndexGlobal]));
393 atomicAdd_g_f(gm_grid + gridIndexGlobal, thetaX * constVal);
398 #endif // COMPILE_SPREAD_HELPERS_ONCE
401 * A spline computation and/or charge spreading kernel function.
402 * Please see the file description for additional defines which this kernel expects.
404 * \param[in] kernelParams Input PME GPU data in constant memory.
405 * \param[in,out] gm_theta Atom spline parameter values
406 * \param[out] gm_dtheta Atom spline parameter derivatives
407 * \param[in,out] gm_gridlineIndices Atom gridline indices (ivec)
408 * \param[out] gm_grid Global 3D grid for charge spreading.
409 * \param[in] gm_fractShiftsTable Atom fractional coordinates correction table
410 * \param[in] gm_gridlineIndicesTable Atom fractional coordinates correction table
411 * \param[in] gm_coefficients Atom charges/coefficients.
412 * \param[in] gm_coordinates Atom coordinates (rvec)
414 __attribute__((reqd_work_group_size(order, order, atomsPerBlock))) __kernel void CUSTOMIZED_KERNEL_NAME(
415 pme_spline_and_spread_kernel)(const struct PmeOpenCLKernelParams kernelParams,
416 __global float* __restrict__ gm_theta,
417 __global float* __restrict__ gm_dtheta,
418 __global int* __restrict__ gm_gridlineIndices,
419 __global float* __restrict__ gm_grid,
420 __global const float* __restrict__ gm_fractShiftsTable,
421 __global const int* __restrict__ gm_gridlineIndicesTable,
422 __global const float* __restrict__ gm_coefficients,
423 __global const float* __restrict__ gm_coordinates)
425 // Gridline indices, ivec
426 __local int sm_gridlineIndices[atomsPerBlock * DIM];
428 __local float sm_coefficients[atomsPerBlock];
430 __local float sm_theta[atomsPerBlock * DIM * order];
431 // Fractional coordinates - only for spline computation
432 __local float sm_fractCoords[atomsPerBlock * DIM];
433 // Staging coordinates - only for spline computation
434 __local float sm_coordinates[DIM * atomsPerBlock];
436 const int atomIndexOffset = (int)get_group_id(XX) * atomsPerBlock;
438 /* Staging coefficients/charges for both spline and spread */
439 pme_gpu_stage_atom_data(sm_coefficients, gm_coefficients, 1);
443 /* Staging coordinates */
444 pme_gpu_stage_atom_data(sm_coordinates, gm_coordinates, DIM);
446 barrier(CLK_LOCAL_MEM_FENCE);
447 calculate_splines(kernelParams, atomIndexOffset, sm_coordinates, sm_coefficients, sm_theta,
448 sm_gridlineIndices, sm_fractCoords, gm_theta, gm_dtheta,
449 gm_gridlineIndices, gm_fractShiftsTable, gm_gridlineIndicesTable);
450 #if !defined(_AMD_SOURCE_) && !defined(_NVIDIA_SOURCE_)
451 /* This is only here for execution of e.g. 32-sized warps on 16-wide hardware; this was
452 * __syncwarp() in CUDA. #2519
454 barrier(CLK_LOCAL_MEM_FENCE);
459 /* Staging the data for spread
460 * (the data is assumed to be in GPU global memory with proper layout already,
461 * as in after running the spline kernel)
463 /* Spline data - only thetas (dthetas will only be needed in gather) */
464 pme_gpu_stage_atom_data(sm_theta, gm_theta, DIM * order);
465 /* Gridline indices - they're actually int and not float, but C99 is angry about overloads */
466 pme_gpu_stage_atom_data((__local float*)sm_gridlineIndices,
467 (__global const float*)gm_gridlineIndices, DIM);
469 barrier(CLK_LOCAL_MEM_FENCE);
474 spread_charges(kernelParams, atomIndexOffset, sm_coefficients, sm_gridlineIndices, sm_theta, gm_grid);