2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2016,2017, 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.
38 * \brief This file contains internal function implementations
39 * for performing the PME calculations on GPU.
41 * \author Aleksei Iupinov <a.yupinov@gmail.com>
42 * \ingroup module_ewald
47 #include "pme-gpu-internal.h"
54 #include "gromacs/gpu_utils/gpu_utils.h"
55 #include "gromacs/math/invertmatrix.h"
56 #include "gromacs/math/units.h"
57 #include "gromacs/mdtypes/commrec.h"
58 #include "gromacs/utility/exceptions.h"
59 #include "gromacs/utility/gmxassert.h"
60 #include "gromacs/utility/logger.h"
61 #include "gromacs/utility/stringutil.h"
64 #include "pme-internal.h"
67 * Wrapper for getting a pointer to the plain C++ part of the GPU kernel parameters structure.
69 * \param[in] pmeGPU The PME GPU structure.
70 * \returns The pointer to the kernel parameters.
72 static pme_gpu_kernel_params_base_t *pme_gpu_get_kernel_params_base_ptr(const pme_gpu_t *pmeGPU)
74 // reinterpret_cast is needed because the derived CUDA structure is not known in this file
75 auto *kernelParamsPtr = reinterpret_cast<pme_gpu_kernel_params_base_t *>(pmeGPU->kernelParams.get());
76 return kernelParamsPtr;
79 void pme_gpu_get_energy_virial(const pme_gpu_t *pmeGPU, real *energy, matrix virial)
81 GMX_ASSERT(energy, "Invalid energy output pointer in PME GPU");
83 virial[XX][XX] = 0.25f * pmeGPU->staging.h_virialAndEnergy[j++];
84 virial[YY][YY] = 0.25f * pmeGPU->staging.h_virialAndEnergy[j++];
85 virial[ZZ][ZZ] = 0.25f * pmeGPU->staging.h_virialAndEnergy[j++];
86 virial[XX][YY] = virial[YY][XX] = 0.25f * pmeGPU->staging.h_virialAndEnergy[j++];
87 virial[XX][ZZ] = virial[ZZ][XX] = 0.25f * pmeGPU->staging.h_virialAndEnergy[j++];
88 virial[YY][ZZ] = virial[ZZ][YY] = 0.25f * pmeGPU->staging.h_virialAndEnergy[j++];
89 *energy = 0.5f * pmeGPU->staging.h_virialAndEnergy[j++];
92 void pme_gpu_update_input_box(pme_gpu_t *pmeGPU, const matrix box)
94 auto *kernelParamsPtr = pme_gpu_get_kernel_params_base_ptr(pmeGPU);
95 kernelParamsPtr->step.boxVolume = box[XX][XX] * box[YY][YY] * box[ZZ][ZZ];
96 GMX_ASSERT(kernelParamsPtr->step.boxVolume != 0.0f, "Zero volume of the unit cell");
99 GMX_THROW(gmx::NotImplementedError("PME is implemented for single-precision only on GPU"));
102 gmx::invertBoxMatrix(box, recipBox);
103 /* The GPU recipBox is transposed as compared to the CPU recipBox.
104 * Spread uses matrix columns (while solve and gather use rows).
105 * There is no particular reason for this; it might be further rethought/optimized for better access patterns.
107 const real newRecipBox[DIM][DIM] =
109 {recipBox[XX][XX], recipBox[YY][XX], recipBox[ZZ][XX]},
110 { 0.0, recipBox[YY][YY], recipBox[ZZ][YY]},
111 { 0.0, 0.0, recipBox[ZZ][ZZ]}
113 memcpy(kernelParamsPtr->step.recipBox, newRecipBox, sizeof(matrix));
117 void pme_gpu_start_step(pme_gpu_t *pmeGPU, bool needToUpdateBox, const matrix box, const rvec *h_coordinates)
119 pme_gpu_copy_input_coordinates(pmeGPU, h_coordinates);
122 pme_gpu_update_input_box(pmeGPU, box);
126 /*! \brief \libinternal
127 * The PME GPU reinitialization function that is called both at the end of any MD step and on any load balancing step.
129 * \param[in] pmeGPU The PME GPU structure.
131 static void pme_gpu_reinit_step(const pme_gpu_t *pmeGPU)
133 pme_gpu_clear_grids(pmeGPU);
134 pme_gpu_clear_energy_virial(pmeGPU);
137 void pme_gpu_finish_step(const pme_gpu_t *pmeGPU, const bool bCalcF, const bool bCalcEnerVir)
139 /* Needed for copy back as well as timing events */
140 pme_gpu_synchronize(pmeGPU);
142 if (bCalcF && pme_gpu_performs_gather(pmeGPU))
144 pme_gpu_sync_output_forces(pmeGPU);
146 if (bCalcEnerVir && pme_gpu_performs_solve(pmeGPU))
148 pme_gpu_sync_output_energy_virial(pmeGPU);
150 pme_gpu_update_timings(pmeGPU);
151 pme_gpu_reinit_step(pmeGPU);
154 /*! \brief \libinternal
155 * (Re-)initializes all the PME GPU data related to the grid size and cut-off.
157 * \param[in] pmeGPU The PME GPU structure.
159 static void pme_gpu_reinit_grids(pme_gpu_t *pmeGPU)
161 auto *kernelParamsPtr = pme_gpu_get_kernel_params_base_ptr(pmeGPU);
162 kernelParamsPtr->grid.ewaldFactor = (M_PI * M_PI) / (pmeGPU->common->ewaldcoeff_q * pmeGPU->common->ewaldcoeff_q);
164 /* The grid size variants */
165 for (int i = 0; i < DIM; i++)
167 kernelParamsPtr->grid.realGridSize[i] = pmeGPU->common->nk[i];
168 kernelParamsPtr->grid.realGridSizeFP[i] = (float)kernelParamsPtr->grid.realGridSize[i];
169 kernelParamsPtr->grid.realGridSizePadded[i] = kernelParamsPtr->grid.realGridSize[i];
171 // The complex grid currently uses no padding;
172 // if it starts to do so, then another test should be added for that
173 kernelParamsPtr->grid.complexGridSize[i] = kernelParamsPtr->grid.realGridSize[i];
174 kernelParamsPtr->grid.complexGridSizePadded[i] = kernelParamsPtr->grid.realGridSize[i];
176 /* FFT: n real elements correspond to (n / 2 + 1) complex elements in minor dimension */
177 if (!pme_gpu_performs_FFT(pmeGPU))
179 // This allows for GPU spreading grid and CPU fftgrid to have the same layout, so that we can copy the data directly
180 kernelParamsPtr->grid.realGridSizePadded[ZZ] = (kernelParamsPtr->grid.realGridSize[ZZ] / 2 + 1) * 2;
183 /* GPU FFT: n real elements correspond to (n / 2 + 1) complex elements in minor dimension */
184 kernelParamsPtr->grid.complexGridSize[ZZ] /= 2;
185 kernelParamsPtr->grid.complexGridSize[ZZ]++;
186 kernelParamsPtr->grid.complexGridSizePadded[ZZ] = kernelParamsPtr->grid.complexGridSize[ZZ];
188 pme_gpu_realloc_and_copy_fract_shifts(pmeGPU);
189 pme_gpu_realloc_and_copy_bspline_values(pmeGPU);
190 pme_gpu_realloc_grids(pmeGPU);
191 pme_gpu_reinit_3dfft(pmeGPU);
194 /* Several GPU functions that refer to the CPU PME data live here.
195 * We would like to keep these away from the GPU-framework specific code for clarity,
196 * as well as compilation issues with MPI.
199 /*! \brief \libinternal
200 * Copies everything useful from the PME CPU to the PME GPU structure.
201 * The goal is to minimize interaction with the PME CPU structure in the GPU code.
203 * \param[in] pme The PME structure.
205 static void pme_gpu_copy_common_data_from(const gmx_pme_t *pme)
207 /* TODO: Consider refactoring the CPU PME code to use the same structure,
208 * so that this function becomes 2 lines */
209 pme_gpu_t *pmeGPU = pme->gpu;
210 pmeGPU->common->ngrids = pme->ngrids;
211 pmeGPU->common->epsilon_r = pme->epsilon_r;
212 pmeGPU->common->ewaldcoeff_q = pme->ewaldcoeff_q;
213 pmeGPU->common->nk[XX] = pme->nkx;
214 pmeGPU->common->nk[YY] = pme->nky;
215 pmeGPU->common->nk[ZZ] = pme->nkz;
216 pmeGPU->common->pmegrid_n[XX] = pme->pmegrid_nx;
217 pmeGPU->common->pmegrid_n[YY] = pme->pmegrid_ny;
218 pmeGPU->common->pmegrid_n[ZZ] = pme->pmegrid_nz;
219 pmeGPU->common->pme_order = pme->pme_order;
220 for (int i = 0; i < DIM; i++)
222 pmeGPU->common->bsp_mod[i].assign(pme->bsp_mod[i], pme->bsp_mod[i] + pmeGPU->common->nk[i]);
224 const int cellCount = c_pmeNeighborUnitcellCount;
225 pmeGPU->common->fsh.resize(0);
226 pmeGPU->common->fsh.insert(pmeGPU->common->fsh.end(), pme->fshx, pme->fshx + cellCount * pme->nkx);
227 pmeGPU->common->fsh.insert(pmeGPU->common->fsh.end(), pme->fshy, pme->fshy + cellCount * pme->nky);
228 pmeGPU->common->fsh.insert(pmeGPU->common->fsh.end(), pme->fshz, pme->fshz + cellCount * pme->nkz);
229 pmeGPU->common->nn.resize(0);
230 pmeGPU->common->nn.insert(pmeGPU->common->nn.end(), pme->nnx, pme->nnx + cellCount * pme->nkx);
231 pmeGPU->common->nn.insert(pmeGPU->common->nn.end(), pme->nny, pme->nny + cellCount * pme->nky);
232 pmeGPU->common->nn.insert(pmeGPU->common->nn.end(), pme->nnz, pme->nnz + cellCount * pme->nkz);
233 pmeGPU->common->runMode = pme->runMode;
236 /*! \brief \libinternal
237 * Finds out if PME with given inputs is possible to run on GPU.
239 * \param[in] pme The PME structure.
240 * \param[out] error The error message if the input is not supported on GPU.
241 * \returns True if this PME input is possible to run on GPU, false otherwise.
243 static bool pme_gpu_check_restrictions(const gmx_pme_t *pme, std::string *error)
245 std::list<std::string> errorReasons;
246 if (pme->nnodes != 1)
248 errorReasons.push_back("PME decomposition");
250 if (pme->pme_order != 4)
252 errorReasons.push_back("interpolation orders other than 4");
256 errorReasons.push_back("free energy calculations (multiple grids)");
260 errorReasons.push_back("Lennard-Jones PME");
264 errorReasons.push_back("double precision");
267 #if GMX_GPU != GMX_GPU_CUDA
269 errorReasons.push_back("non-CUDA build of GROMACS");
273 bool inputSupported = errorReasons.empty();
274 if (!inputSupported && error)
276 std::string regressionTestMarker = "PME GPU does not support";
277 // this prefix is tested for in the regression tests script gmxtest.pl
278 *error = regressionTestMarker + ": " + gmx::joinStrings(errorReasons, "; ") + ".";
280 return inputSupported;
283 /*! \libinternal \brief
284 * Initializes the PME GPU data at the beginning of the run.
286 * \param[in,out] pme The PME structure.
287 * \param[in,out] gpuInfo The GPU information structure.
288 * \param[in] mdlog The logger.
289 * \param[in] cr The communication structure.
291 static void pme_gpu_init(gmx_pme_t *pme, gmx_device_info_t *gpuInfo, const gmx::MDLogger &mdlog,
294 std::string errorString;
295 bool canRunOnGpu = pme_gpu_check_restrictions(pme, &errorString);
298 GMX_THROW(gmx::NotImplementedError(errorString));
301 pme->gpu = new pme_gpu_t();
302 pme_gpu_t *pmeGPU = pme->gpu;
303 pmeGPU->common = std::shared_ptr<pme_shared_t>(new pme_shared_t());
305 /* These settings are set here for the whole run; dynamic ones are set in pme_gpu_reinit() */
306 /* A convenience variable. */
307 pmeGPU->settings.useDecomposition = (pme->nnodes == 1);
308 /* TODO: CPU gather with GPU spread is broken due to different theta/dtheta layout. */
309 pmeGPU->settings.performGPUGather = true;
311 pme_gpu_set_testing(pmeGPU, false);
313 // GPU initialization
314 init_gpu(mdlog, cr->nodeid, gpuInfo);
315 pmeGPU->deviceInfo = gpuInfo;
317 pme_gpu_init_internal(pmeGPU);
318 pme_gpu_init_sync_events(pmeGPU);
319 pme_gpu_alloc_energy_virial(pmeGPU);
321 pme_gpu_copy_common_data_from(pme);
323 GMX_ASSERT(pmeGPU->common->epsilon_r != 0.0f, "PME GPU: bad electrostatic coefficient");
325 auto *kernelParamsPtr = pme_gpu_get_kernel_params_base_ptr(pmeGPU);
326 kernelParamsPtr->constants.elFactor = ONE_4PI_EPS0 / pmeGPU->common->epsilon_r;
329 void pme_gpu_transform_spline_atom_data_for_host(const pme_gpu_t *pmeGPU, const pme_atomcomm_t *atc, PmeSplineDataType type, int dimIndex)
331 // The GPU atom spline data is laid out in a different way currently than the CPU one.
332 // This function converts the data from GPU to CPU layout. Obviously, it is slow.
333 // It is only intended for testing purposes so far.
334 // Ideally we should use similar layouts on CPU and GPU if we care about mixed modes and their performance
335 // (e.g. spreading on GPU, gathering on CPU).
336 GMX_RELEASE_ASSERT(atc->nthread == 1, "Only the serial PME data layout is supported");
337 const uintmax_t threadIndex = 0;
338 const auto atomCount = pme_gpu_get_kernel_params_base_ptr(pmeGPU)->atoms.nAtoms;
339 const auto atomsPerWarp = pme_gpu_get_atom_spline_data_alignment(pmeGPU);
340 const auto pmeOrder = pmeGPU->common->pme_order;
342 real *h_splineBuffer;
343 float *d_splineBuffer;
346 case PmeSplineDataType::Values:
347 h_splineBuffer = atc->spline[threadIndex].theta[dimIndex];
348 d_splineBuffer = pmeGPU->staging.h_theta;
351 case PmeSplineDataType::Derivatives:
352 h_splineBuffer = atc->spline[threadIndex].dtheta[dimIndex];
353 d_splineBuffer = pmeGPU->staging.h_dtheta;
357 GMX_THROW(gmx::InternalError("Unknown spline data type"));
360 for (auto atomIndex = 0; atomIndex < atomCount; atomIndex++)
362 auto atomWarpIndex = atomIndex % atomsPerWarp;
363 auto warpIndex = atomIndex / atomsPerWarp;
364 for (auto orderIndex = 0; orderIndex < pmeOrder; orderIndex++)
366 const auto gpuValueIndex = ((pmeOrder * warpIndex + orderIndex) * DIM + dimIndex) * atomsPerWarp + atomWarpIndex;
367 const auto cpuValueIndex = atomIndex * pmeOrder + orderIndex;
368 GMX_ASSERT(cpuValueIndex < atomCount * pmeOrder, "Atom spline data index out of bounds (while transforming GPU data layout for host)");
369 h_splineBuffer[cpuValueIndex] = d_splineBuffer[gpuValueIndex];
374 void pme_gpu_get_real_grid_sizes(const pme_gpu_t *pmeGPU, gmx::IVec *gridSize, gmx::IVec *paddedGridSize)
376 GMX_ASSERT(gridSize != nullptr, "");
377 GMX_ASSERT(paddedGridSize != nullptr, "");
378 GMX_ASSERT(pmeGPU != nullptr, "");
379 auto *kernelParamsPtr = pme_gpu_get_kernel_params_base_ptr(pmeGPU);
380 for (int i = 0; i < DIM; i++)
382 (*gridSize)[i] = kernelParamsPtr->grid.realGridSize[i];
383 (*paddedGridSize)[i] = kernelParamsPtr->grid.realGridSizePadded[i];
387 void pme_gpu_reinit(gmx_pme_t *pme, gmx_device_info_t *gpuInfo, const gmx::MDLogger &mdlog, const t_commrec *cr)
389 if (!pme_gpu_active(pme))
396 /* First-time initialization */
397 pme_gpu_init(pme, gpuInfo, mdlog, cr);
401 /* After this call nothing in the GPU code should refer to the gmx_pme_t *pme itself - until the next pme_gpu_reinit */
402 pme_gpu_copy_common_data_from(pme);
404 /* GPU FFT will only get used for a single rank.*/
405 pme->gpu->settings.performGPUFFT = (pme->gpu->common->runMode == PmeRunMode::GPU) && !pme_gpu_uses_dd(pme->gpu);
406 pme->gpu->settings.performGPUSolve = (pme->gpu->common->runMode == PmeRunMode::GPU);
408 pme_gpu_reinit_grids(pme->gpu);
409 pme_gpu_reinit_step(pme->gpu);
412 void pme_gpu_destroy(pme_gpu_t *pmeGPU)
414 /* Free lots of data */
415 pme_gpu_free_energy_virial(pmeGPU);
416 pme_gpu_free_bspline_values(pmeGPU);
417 pme_gpu_free_forces(pmeGPU);
418 pme_gpu_free_coordinates(pmeGPU);
419 pme_gpu_free_coefficients(pmeGPU);
420 pme_gpu_free_spline_data(pmeGPU);
421 pme_gpu_free_grid_indices(pmeGPU);
422 pme_gpu_free_fract_shifts(pmeGPU);
423 pme_gpu_free_grids(pmeGPU);
425 pme_gpu_destroy_3dfft(pmeGPU);
426 pme_gpu_destroy_sync_events(pmeGPU);
428 /* Free the GPU-framework specific data last */
429 pme_gpu_destroy_specific(pmeGPU);
434 void pme_gpu_reinit_atoms(pme_gpu_t *pmeGPU, const int nAtoms, const real *charges)
436 auto *kernelParamsPtr = pme_gpu_get_kernel_params_base_ptr(pmeGPU);
437 kernelParamsPtr->atoms.nAtoms = nAtoms;
438 const int alignment = pme_gpu_get_atom_data_alignment(pmeGPU);
439 pmeGPU->nAtomsPadded = ((nAtoms + alignment - 1) / alignment) * alignment;
440 const int nAtomsAlloc = c_usePadding ? pmeGPU->nAtomsPadded : nAtoms;
441 const bool haveToRealloc = (pmeGPU->nAtomsAlloc < nAtomsAlloc); /* This check might be redundant, but is logical */
442 pmeGPU->nAtomsAlloc = nAtomsAlloc;
445 GMX_RELEASE_ASSERT(false, "Only single precision supported");
446 GMX_UNUSED_VALUE(charges);
448 pme_gpu_realloc_and_copy_input_coefficients(pmeGPU, reinterpret_cast<const float *>(charges));
449 /* Could also be checked for haveToRealloc, but the copy always needs to be performed */
454 pme_gpu_realloc_coordinates(pmeGPU);
455 pme_gpu_realloc_forces(pmeGPU);
456 pme_gpu_realloc_spline_data(pmeGPU);
457 pme_gpu_realloc_grid_indices(pmeGPU);