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/ewald/ewald-utils.h"
55 #include "gromacs/gpu_utils/gpu_utils.h"
56 #include "gromacs/math/invertmatrix.h"
57 #include "gromacs/math/units.h"
58 #include "gromacs/mdtypes/commrec.h"
59 #include "gromacs/utility/exceptions.h"
60 #include "gromacs/utility/gmxassert.h"
61 #include "gromacs/utility/logger.h"
62 #include "gromacs/utility/stringutil.h"
65 #include "pme-internal.h"
68 * Wrapper for getting a pointer to the plain C++ part of the GPU kernel parameters structure.
70 * \param[in] pmeGpu The PME GPU structure.
71 * \returns The pointer to the kernel parameters.
73 static PmeGpuKernelParamsBase *pme_gpu_get_kernel_params_base_ptr(const PmeGpu *pmeGpu)
75 // reinterpret_cast is needed because the derived CUDA structure is not known in this file
76 auto *kernelParamsPtr = reinterpret_cast<PmeGpuKernelParamsBase *>(pmeGpu->kernelParams.get());
77 return kernelParamsPtr;
80 gmx::ArrayRef<gmx::RVec> pme_gpu_get_forces(PmeGpu *pmeGpu)
82 return pmeGpu->staging.h_forces;
85 void pme_gpu_get_energy_virial(const PmeGpu *pmeGpu, real *energy, matrix virial)
87 for (int j = 0; j < c_virialAndEnergyCount; j++)
89 GMX_ASSERT(std::isfinite(pmeGpu->staging.h_virialAndEnergy[j]), "PME GPU produces incorrect energy/virial.");
92 GMX_ASSERT(energy, "Invalid energy output pointer in PME GPU");
94 virial[XX][XX] = 0.25f * pmeGpu->staging.h_virialAndEnergy[j++];
95 virial[YY][YY] = 0.25f * pmeGpu->staging.h_virialAndEnergy[j++];
96 virial[ZZ][ZZ] = 0.25f * pmeGpu->staging.h_virialAndEnergy[j++];
97 virial[XX][YY] = virial[YY][XX] = 0.25f * pmeGpu->staging.h_virialAndEnergy[j++];
98 virial[XX][ZZ] = virial[ZZ][XX] = 0.25f * pmeGpu->staging.h_virialAndEnergy[j++];
99 virial[YY][ZZ] = virial[ZZ][YY] = 0.25f * pmeGpu->staging.h_virialAndEnergy[j++];
100 *energy = 0.5f * pmeGpu->staging.h_virialAndEnergy[j++];
103 void pme_gpu_update_input_box(PmeGpu gmx_unused *pmeGpu,
104 const matrix gmx_unused box)
107 GMX_THROW(gmx::NotImplementedError("PME is implemented for single-precision only on GPU"));
110 pmeGpu->common->boxScaler->scaleBox(box, scaledBox);
111 auto *kernelParamsPtr = pme_gpu_get_kernel_params_base_ptr(pmeGpu);
112 kernelParamsPtr->current.boxVolume = scaledBox[XX][XX] * scaledBox[YY][YY] * scaledBox[ZZ][ZZ];
113 GMX_ASSERT(kernelParamsPtr->current.boxVolume != 0.0f, "Zero volume of the unit cell");
115 gmx::invertBoxMatrix(scaledBox, recipBox);
117 /* The GPU recipBox is transposed as compared to the CPU recipBox.
118 * Spread uses matrix columns (while solve and gather use rows).
119 * There is no particular reason for this; it might be further rethought/optimized for better access patterns.
121 const real newRecipBox[DIM][DIM] =
123 {recipBox[XX][XX], recipBox[YY][XX], recipBox[ZZ][XX]},
124 { 0.0, recipBox[YY][YY], recipBox[ZZ][YY]},
125 { 0.0, 0.0, recipBox[ZZ][ZZ]}
127 memcpy(kernelParamsPtr->current.recipBox, newRecipBox, sizeof(matrix));
131 /*! \brief \libinternal
132 * The PME GPU reinitialization function that is called both at the end of any PME computation and on any load balancing.
134 * \param[in] pmeGpu The PME GPU structure.
136 void pme_gpu_reinit_computation(const PmeGpu *pmeGpu)
138 pme_gpu_clear_grids(pmeGpu);
139 pme_gpu_clear_energy_virial(pmeGpu);
142 /*! \brief \libinternal
143 * (Re-)initializes all the PME GPU data related to the grid size and cut-off.
145 * \param[in] pmeGpu The PME GPU structure.
147 static void pme_gpu_reinit_grids(PmeGpu *pmeGpu)
149 auto *kernelParamsPtr = pme_gpu_get_kernel_params_base_ptr(pmeGpu);
150 kernelParamsPtr->grid.ewaldFactor = (M_PI * M_PI) / (pmeGpu->common->ewaldcoeff_q * pmeGpu->common->ewaldcoeff_q);
152 /* The grid size variants */
153 for (int i = 0; i < DIM; i++)
155 kernelParamsPtr->grid.realGridSize[i] = pmeGpu->common->nk[i];
156 kernelParamsPtr->grid.realGridSizeFP[i] = (float)kernelParamsPtr->grid.realGridSize[i];
157 kernelParamsPtr->grid.realGridSizePadded[i] = kernelParamsPtr->grid.realGridSize[i];
159 // The complex grid currently uses no padding;
160 // if it starts to do so, then another test should be added for that
161 kernelParamsPtr->grid.complexGridSize[i] = kernelParamsPtr->grid.realGridSize[i];
162 kernelParamsPtr->grid.complexGridSizePadded[i] = kernelParamsPtr->grid.realGridSize[i];
164 /* FFT: n real elements correspond to (n / 2 + 1) complex elements in minor dimension */
165 if (!pme_gpu_performs_FFT(pmeGpu))
167 // This allows for GPU spreading grid and CPU fftgrid to have the same layout, so that we can copy the data directly
168 kernelParamsPtr->grid.realGridSizePadded[ZZ] = (kernelParamsPtr->grid.realGridSize[ZZ] / 2 + 1) * 2;
171 /* GPU FFT: n real elements correspond to (n / 2 + 1) complex elements in minor dimension */
172 kernelParamsPtr->grid.complexGridSize[ZZ] /= 2;
173 kernelParamsPtr->grid.complexGridSize[ZZ]++;
174 kernelParamsPtr->grid.complexGridSizePadded[ZZ] = kernelParamsPtr->grid.complexGridSize[ZZ];
176 pme_gpu_realloc_and_copy_fract_shifts(pmeGpu);
177 pme_gpu_realloc_and_copy_bspline_values(pmeGpu);
178 pme_gpu_realloc_grids(pmeGpu);
179 pme_gpu_reinit_3dfft(pmeGpu);
182 /* Several GPU functions that refer to the CPU PME data live here.
183 * We would like to keep these away from the GPU-framework specific code for clarity,
184 * as well as compilation issues with MPI.
187 /*! \brief \libinternal
188 * Copies everything useful from the PME CPU to the PME GPU structure.
189 * The goal is to minimize interaction with the PME CPU structure in the GPU code.
191 * \param[in] pme The PME structure.
193 static void pme_gpu_copy_common_data_from(const gmx_pme_t *pme)
195 /* TODO: Consider refactoring the CPU PME code to use the same structure,
196 * so that this function becomes 2 lines */
197 PmeGpu *pmeGpu = pme->gpu;
198 pmeGpu->common->ngrids = pme->ngrids;
199 pmeGpu->common->epsilon_r = pme->epsilon_r;
200 pmeGpu->common->ewaldcoeff_q = pme->ewaldcoeff_q;
201 pmeGpu->common->nk[XX] = pme->nkx;
202 pmeGpu->common->nk[YY] = pme->nky;
203 pmeGpu->common->nk[ZZ] = pme->nkz;
204 pmeGpu->common->pmegrid_n[XX] = pme->pmegrid_nx;
205 pmeGpu->common->pmegrid_n[YY] = pme->pmegrid_ny;
206 pmeGpu->common->pmegrid_n[ZZ] = pme->pmegrid_nz;
207 pmeGpu->common->pme_order = pme->pme_order;
208 for (int i = 0; i < DIM; i++)
210 pmeGpu->common->bsp_mod[i].assign(pme->bsp_mod[i], pme->bsp_mod[i] + pmeGpu->common->nk[i]);
212 const int cellCount = c_pmeNeighborUnitcellCount;
213 pmeGpu->common->fsh.resize(0);
214 pmeGpu->common->fsh.insert(pmeGpu->common->fsh.end(), pme->fshx, pme->fshx + cellCount * pme->nkx);
215 pmeGpu->common->fsh.insert(pmeGpu->common->fsh.end(), pme->fshy, pme->fshy + cellCount * pme->nky);
216 pmeGpu->common->fsh.insert(pmeGpu->common->fsh.end(), pme->fshz, pme->fshz + cellCount * pme->nkz);
217 pmeGpu->common->nn.resize(0);
218 pmeGpu->common->nn.insert(pmeGpu->common->nn.end(), pme->nnx, pme->nnx + cellCount * pme->nkx);
219 pmeGpu->common->nn.insert(pmeGpu->common->nn.end(), pme->nny, pme->nny + cellCount * pme->nky);
220 pmeGpu->common->nn.insert(pmeGpu->common->nn.end(), pme->nnz, pme->nnz + cellCount * pme->nkz);
221 pmeGpu->common->runMode = pme->runMode;
222 pmeGpu->common->boxScaler = pme->boxScaler;
225 /*! \brief \libinternal
226 * Finds out if PME with given inputs is possible to run on GPU.
228 * \param[in] pme The PME structure.
229 * \param[out] error The error message if the input is not supported on GPU.
230 * \returns True if this PME input is possible to run on GPU, false otherwise.
232 static bool pme_gpu_check_restrictions(const gmx_pme_t *pme, std::string *error)
234 std::list<std::string> errorReasons;
235 if (pme->nnodes != 1)
237 errorReasons.push_back("PME decomposition");
239 if (pme->pme_order != 4)
241 errorReasons.push_back("interpolation orders other than 4");
245 errorReasons.push_back("free energy calculations (multiple grids)");
249 errorReasons.push_back("Lennard-Jones PME");
253 errorReasons.push_back("double precision");
256 #if GMX_GPU != GMX_GPU_CUDA
258 errorReasons.push_back("non-CUDA build of GROMACS");
262 bool inputSupported = errorReasons.empty();
263 if (!inputSupported && error)
265 std::string regressionTestMarker = "PME GPU does not support";
266 // this prefix is tested for in the regression tests script gmxtest.pl
267 *error = regressionTestMarker + ": " + gmx::joinStrings(errorReasons, "; ") + ".";
269 return inputSupported;
272 /*! \libinternal \brief
273 * Initializes the PME GPU data at the beginning of the run.
275 * \param[in,out] pme The PME structure.
276 * \param[in,out] gpuInfo The GPU information structure.
278 static void pme_gpu_init(gmx_pme_t *pme, gmx_device_info_t *gpuInfo)
280 std::string errorString;
281 bool canRunOnGpu = pme_gpu_check_restrictions(pme, &errorString);
284 GMX_THROW(gmx::NotImplementedError(errorString));
287 pme->gpu = new PmeGpu();
288 PmeGpu *pmeGpu = pme->gpu;
289 changePinningPolicy(&pmeGpu->staging.h_forces, gmx::PinningPolicy::CanBePinned);
290 pmeGpu->common = std::shared_ptr<PmeShared>(new PmeShared());
292 /* These settings are set here for the whole run; dynamic ones are set in pme_gpu_reinit() */
293 /* A convenience variable. */
294 pmeGpu->settings.useDecomposition = (pme->nnodes == 1);
295 /* TODO: CPU gather with GPU spread is broken due to different theta/dtheta layout. */
296 pmeGpu->settings.performGPUGather = true;
298 pme_gpu_set_testing(pmeGpu, false);
300 pmeGpu->deviceInfo = gpuInfo;
302 pme_gpu_init_internal(pmeGpu);
303 pme_gpu_init_sync_events(pmeGpu);
304 pme_gpu_alloc_energy_virial(pmeGpu);
306 pme_gpu_copy_common_data_from(pme);
308 GMX_ASSERT(pmeGpu->common->epsilon_r != 0.0f, "PME GPU: bad electrostatic coefficient");
310 auto *kernelParamsPtr = pme_gpu_get_kernel_params_base_ptr(pmeGpu);
311 kernelParamsPtr->constants.elFactor = ONE_4PI_EPS0 / pmeGpu->common->epsilon_r;
314 void pme_gpu_transform_spline_atom_data(const PmeGpu *pmeGpu, const pme_atomcomm_t *atc,
315 PmeSplineDataType type, int dimIndex, PmeLayoutTransform transform)
317 // The GPU atom spline data is laid out in a different way currently than the CPU one.
318 // This function converts the data from GPU to CPU layout (in the host memory).
319 // It is only intended for testing purposes so far.
320 // Ideally we should use similar layouts on CPU and GPU if we care about mixed modes and their performance
321 // (e.g. spreading on GPU, gathering on CPU).
322 GMX_RELEASE_ASSERT(atc->nthread == 1, "Only the serial PME data layout is supported");
323 const uintmax_t threadIndex = 0;
324 const auto atomCount = pme_gpu_get_kernel_params_base_ptr(pmeGpu)->atoms.nAtoms;
325 const auto atomsPerWarp = pme_gpu_get_atoms_per_warp(pmeGpu);
326 const auto pmeOrder = pmeGpu->common->pme_order;
328 real *cpuSplineBuffer;
329 float *h_splineBuffer;
332 case PmeSplineDataType::Values:
333 cpuSplineBuffer = atc->spline[threadIndex].theta[dimIndex];
334 h_splineBuffer = pmeGpu->staging.h_theta;
337 case PmeSplineDataType::Derivatives:
338 cpuSplineBuffer = atc->spline[threadIndex].dtheta[dimIndex];
339 h_splineBuffer = pmeGpu->staging.h_dtheta;
343 GMX_THROW(gmx::InternalError("Unknown spline data type"));
346 for (auto atomIndex = 0; atomIndex < atomCount; atomIndex++)
348 auto atomWarpIndex = atomIndex % atomsPerWarp;
349 auto warpIndex = atomIndex / atomsPerWarp;
350 for (auto orderIndex = 0; orderIndex < pmeOrder; orderIndex++)
352 const auto gpuValueIndex = ((pmeOrder * warpIndex + orderIndex) * DIM + dimIndex) * atomsPerWarp + atomWarpIndex;
353 const auto cpuValueIndex = atomIndex * pmeOrder + orderIndex;
354 GMX_ASSERT(cpuValueIndex < atomCount * pmeOrder, "Atom spline data index out of bounds (while transforming GPU data layout for host)");
357 case PmeLayoutTransform::GpuToHost:
358 cpuSplineBuffer[cpuValueIndex] = h_splineBuffer[gpuValueIndex];
361 case PmeLayoutTransform::HostToGpu:
362 h_splineBuffer[gpuValueIndex] = cpuSplineBuffer[cpuValueIndex];
366 GMX_THROW(gmx::InternalError("Unknown layout transform"));
372 void pme_gpu_get_real_grid_sizes(const PmeGpu *pmeGpu, gmx::IVec *gridSize, gmx::IVec *paddedGridSize)
374 GMX_ASSERT(gridSize != nullptr, "");
375 GMX_ASSERT(paddedGridSize != nullptr, "");
376 GMX_ASSERT(pmeGpu != nullptr, "");
377 auto *kernelParamsPtr = pme_gpu_get_kernel_params_base_ptr(pmeGpu);
378 for (int i = 0; i < DIM; i++)
380 (*gridSize)[i] = kernelParamsPtr->grid.realGridSize[i];
381 (*paddedGridSize)[i] = kernelParamsPtr->grid.realGridSizePadded[i];
385 void pme_gpu_reinit(gmx_pme_t *pme, gmx_device_info_t *gpuInfo)
387 if (!pme_gpu_active(pme))
394 /* First-time initialization */
395 pme_gpu_init(pme, gpuInfo);
399 /* After this call nothing in the GPU code should refer to the gmx_pme_t *pme itself - until the next pme_gpu_reinit */
400 pme_gpu_copy_common_data_from(pme);
402 /* GPU FFT will only get used for a single rank.*/
403 pme->gpu->settings.performGPUFFT = (pme->gpu->common->runMode == PmeRunMode::GPU) && !pme_gpu_uses_dd(pme->gpu);
404 pme->gpu->settings.performGPUSolve = (pme->gpu->common->runMode == PmeRunMode::GPU);
406 /* Reinit active timers */
407 pme_gpu_reinit_timings(pme->gpu);
409 pme_gpu_reinit_grids(pme->gpu);
410 pme_gpu_reinit_computation(pme->gpu);
411 /* Clear the previous box - doesn't hurt, and forces the PME CPU recipbox
412 * update for mixed mode on grid switch. TODO: use shared recipbox field.
414 std::memset(pme->gpu->common->previousBox, 0, sizeof(pme->gpu->common->previousBox));
417 void pme_gpu_destroy(PmeGpu *pmeGpu)
419 /* Free lots of data */
420 pme_gpu_free_energy_virial(pmeGpu);
421 pme_gpu_free_bspline_values(pmeGpu);
422 pme_gpu_free_forces(pmeGpu);
423 pme_gpu_free_coordinates(pmeGpu);
424 pme_gpu_free_coefficients(pmeGpu);
425 pme_gpu_free_spline_data(pmeGpu);
426 pme_gpu_free_grid_indices(pmeGpu);
427 pme_gpu_free_fract_shifts(pmeGpu);
428 pme_gpu_free_grids(pmeGpu);
430 pme_gpu_destroy_3dfft(pmeGpu);
431 pme_gpu_destroy_sync_events(pmeGpu);
433 /* Free the GPU-framework specific data last */
434 pme_gpu_destroy_specific(pmeGpu);
439 void pme_gpu_reinit_atoms(PmeGpu *pmeGpu, const int nAtoms, const real *charges)
441 auto *kernelParamsPtr = pme_gpu_get_kernel_params_base_ptr(pmeGpu);
442 kernelParamsPtr->atoms.nAtoms = nAtoms;
443 const int alignment = pme_gpu_get_atom_data_alignment(pmeGpu);
444 pmeGpu->nAtomsPadded = ((nAtoms + alignment - 1) / alignment) * alignment;
445 const int nAtomsAlloc = c_usePadding ? pmeGpu->nAtomsPadded : nAtoms;
446 const bool haveToRealloc = (pmeGpu->nAtomsAlloc < nAtomsAlloc); /* This check might be redundant, but is logical */
447 pmeGpu->nAtomsAlloc = nAtomsAlloc;
450 GMX_RELEASE_ASSERT(false, "Only single precision supported");
451 GMX_UNUSED_VALUE(charges);
453 pme_gpu_realloc_and_copy_input_coefficients(pmeGpu, reinterpret_cast<const float *>(charges));
454 /* Could also be checked for haveToRealloc, but the copy always needs to be performed */
459 pme_gpu_realloc_coordinates(pmeGpu);
460 pme_gpu_realloc_forces(pmeGpu);
461 pme_gpu_realloc_spline_data(pmeGpu);
462 pme_gpu_realloc_grid_indices(pmeGpu);