diff --git a/src/nbnxm-cuda/Makefile b/src/nbnxm-cuda/Makefile index ecf5ff9fd..de735862d 100644 --- a/src/nbnxm-cuda/Makefile +++ b/src/nbnxm-cuda/Makefile @@ -46,7 +46,7 @@ endif $(program): $(obj) Makefile $(CC) $(CFLAGS) $(obj) -o $@ $(LDFLAGS) -%.o: %.cu Makefile +%.o: %.cu reference.h Makefile $(CC) $(CFLAGS) -c $< -o $@ clean: diff --git a/src/nbnxm-cuda/main.cu b/src/nbnxm-cuda/main.cu index 062ecc2e6..109cba75a 100644 --- a/src/nbnxm-cuda/main.cu +++ b/src/nbnxm-cuda/main.cu @@ -1,5 +1,6 @@ #include #include +#include #include #include "vectypes.h" @@ -8,6 +9,7 @@ typedef gmx::BasicVector Float3; typedef float4 Float4; #include "constants.h" +#include "reference.h" #if (CUDART_VERSION >= 9000) #define __shfl_up(v, d) __shfl_up_sync(0xffffffff, v, d) @@ -376,21 +378,16 @@ nbnxn_cj4_t get_cj4(int id) { value.cj[i] = i + id; } for (int i = 0; i < c_nbnxnGpuClusterpairSplit; ++i) { - value.imei[i].imask = 0U; - value.imei[i].excl_ind = 0; + value.imei[i].imask = 0xFFFFFFFFu; + value.imei[i].excl_ind = id % 19205; } return value; } -nbnxn_sci_t get_sci(int id) { - return {id, 0, 8 * id, 8 * id + 7}; -} - nbnxn_excl_t get_excl(int id) { nbnxn_excl_t value; - for (int i = 0; i < c_nbnxnGpuExclSize; ++i) { - value.pair[i] = 7; - } + std::mt19937 rng(id); + for (int i = 0; i < c_nbnxnGpuExclSize; ++i) value.pair[i] = rng(); return value; } @@ -404,92 +401,144 @@ int main(int argc, char* argv[]) { const dim3 blocks ( block_x, block_y, 1 ); const dim3 grids ( 1, 1, grid_z ); - Float4* a_xq; - cudaMallocManaged(&a_xq, sizeof(Float4) * NUM_ATOMS); - - Float3* a_f; - cudaMallocManaged(&a_f, sizeof(Float3) * NUM_ATOMS); - - Float3* shiftVec; - cudaMallocManaged(&shiftVec, sizeof(Float3) * 45); - - Float3* fShift; - cudaMallocManaged(&fShift, sizeof(Float3) * 45); - - nbnxn_cj4_t* cj4; - cudaMallocManaged(&cj4, sizeof(nbnxn_cj4_t) * 56881); - - nbnxn_sci_t* sci; - cudaMallocManaged(&sci, sizeof(nbnxn_sci_t) * 4806); - - nbnxn_excl_t* excl; - cudaMallocManaged(&excl, sizeof(nbnxn_excl_t) * 19205); - - int* atomTypes; - cudaMallocManaged(&atomTypes, sizeof(int) * NUM_ATOMS); - - Float2* nbfp; - cudaMallocManaged(&nbfp, sizeof(Float2) * 1024); + // ----------------------------------------------------------------- + // Host buffers (h_ prefix) + // ----------------------------------------------------------------- + // These hold the "source of truth" data on the CPU. They are used to: + // (1) stage initial values into the device buffers via cudaMemcpy, and + // (2) feed NbnxmReference (reference.h), which is host-only code and + // must be able to dereference xq/shiftVec/cj4/sci/excl/atomTypes/ + // nbfp directly -- it can no longer read the d_* buffers now that + // they are plain (non-managed) device memory. + // a_f / fShift do not need a persistent host mirror for staging (their + // initial values are trivial to re-write directly into h_f / h_fShift + // each time we reset them), but h_f / h_fShift ARE needed as the landing + // buffer for the periodic cudaMemcpyDeviceToHost reads used to validate + // against the CPU reference. + Float4* h_xq = new Float4[NUM_ATOMS]; + Float3* h_f = new Float3[NUM_ATOMS]; + Float3* h_shiftVec = new Float3[45]; + Float3* h_fShift = new Float3[45]; + nbnxn_cj4_t* h_cj4 = new nbnxn_cj4_t[56881]; + nbnxn_sci_t* h_sci = new nbnxn_sci_t[4806]; + nbnxn_excl_t* h_excl = new nbnxn_excl_t[19205]; + int* h_atomTypes = new int[NUM_ATOMS]; + Float2* h_nbfp = new Float2[1024]; + + // ----------------------------------------------------------------- + // Device buffers (d_ prefix) + // ----------------------------------------------------------------- + // Plain cudaMalloc (NOT cudaMallocManaged): per PR #326's rationale, + // atomics on cudaMallocManaged memory can be pathologically slow on + // unified-memory architectures because they require system-scope + // (CPU-GPU) coherency, even though the host side is never touched + // during the timed kernel loops here. nbnxmKernelTest() only + // atomicAdd()s into d_f / d_fShift and only reads the rest, so plain + // device memory + explicit H2D/D2H staging is strictly faster with + // identical results. + Float4* d_xq; + Float3* d_f; + Float3* d_shiftVec; + Float3* d_fShift; + nbnxn_cj4_t* d_cj4; + nbnxn_sci_t* d_sci; + nbnxn_excl_t* d_excl; + int* d_atomTypes; + Float2* d_nbfp; + + cudaMalloc(&d_xq, sizeof(Float4) * NUM_ATOMS); + cudaMalloc(&d_f, sizeof(Float3) * NUM_ATOMS); + cudaMalloc(&d_shiftVec, sizeof(Float3) * 45); + cudaMalloc(&d_fShift, sizeof(Float3) * 45); + cudaMalloc(&d_cj4, sizeof(nbnxn_cj4_t) * 56881); + cudaMalloc(&d_sci, sizeof(nbnxn_sci_t) * 4806); + cudaMalloc(&d_excl, sizeof(nbnxn_excl_t) * 19205); + cudaMalloc(&d_atomTypes, sizeof(int) * NUM_ATOMS); + cudaMalloc(&d_nbfp, sizeof(Float2) * 1024); + + // ----------------------------------------------------------------- + // Populate host buffers with the same deterministic/synthetic data as + // before, then copy H2D once (these inputs are re-used, unchanged, for + // every launch in both the "w/o shift" and "w/ shift" blocks below). + // ----------------------------------------------------------------- + std::mt19937 rng(1337); + std::uniform_real_distribution posDist(-20.f, 20.f); + std::uniform_real_distribution qDist(-1.f, 1.f); for (int i = 0; i < NUM_ATOMS; ++i) { - a_xq[i] = make_float4(1.0f, 0.5f, 0.25f, 0.125f); + h_xq[i] = make_float4(posDist(rng), posDist(rng), posDist(rng), qDist(rng)); } for (int i = 0; i < NUM_ATOMS; ++i) { - a_f[i] = Float3(1.0f, 0.5f, 0.25f); + h_f[i] = Float3(1.0f, 0.5f, 0.25f); } for (int i = 0; i < 45; ++i) { - shiftVec[i] = Float3(1.0f, 0.5f, 0.25f); + h_shiftVec[i] = Float3(posDist(rng)*0.1f, posDist(rng)*0.1f, posDist(rng)*0.1f); } for (int i = 0; i < 45; ++i) { - fShift[i] = Float3(1.0f, 0.5f, 0.25f); + h_fShift[i] = Float3(1.0f, 0.5f, 0.25f); } for (int i = 0; i < 56881; ++i) { - cj4[i] = get_cj4(i); + h_cj4[i] = get_cj4(i % 200); } for (int i = 0; i < 4806; ++i) { - sci[i] = get_sci(i); + h_sci[i] = {i % 400, i % c_numIvecs, (2*i) % 200, (2*i) % 200 + 1}; } for (int i = 0; i < 19205; ++i) { - excl[i] = get_excl(i); + h_excl[i] = get_excl(i); } for (int i = 0; i < NUM_ATOMS; ++i) { - atomTypes[i] = (i % 2); + h_atomTypes[i] = (i % 32); } for (int i = 0; i < 1024; ++i) { - nbfp[i] = make_float2(0.5f, 0.25f); + h_nbfp[i] = make_float2(0.5f, 0.25f); } + cudaMemcpy(d_xq, h_xq, sizeof(Float4) * NUM_ATOMS, cudaMemcpyHostToDevice); + cudaMemcpy(d_f, h_f, sizeof(Float3) * NUM_ATOMS, cudaMemcpyHostToDevice); + cudaMemcpy(d_shiftVec, h_shiftVec, sizeof(Float3) * 45, cudaMemcpyHostToDevice); + cudaMemcpy(d_fShift, h_fShift, sizeof(Float3) * 45, cudaMemcpyHostToDevice); + cudaMemcpy(d_cj4, h_cj4, sizeof(nbnxn_cj4_t) * 56881, cudaMemcpyHostToDevice); + cudaMemcpy(d_sci, h_sci, sizeof(nbnxn_sci_t) * 4806, cudaMemcpyHostToDevice); + cudaMemcpy(d_excl, h_excl, sizeof(nbnxn_excl_t) * 19205, cudaMemcpyHostToDevice); + cudaMemcpy(d_atomTypes, h_atomTypes, sizeof(int) * NUM_ATOMS, cudaMemcpyHostToDevice); + cudaMemcpy(d_nbfp, h_nbfp, sizeof(Float2) * 1024, cudaMemcpyHostToDevice); + + // NbnxmReference reads its inputs directly from host memory -- pass the + // h_* buffers, not the d_* ones (it cannot dereference device pointers). + NbnxmReference ref(h_xq, h_shiftVec, h_cj4, h_sci, h_excl, h_atomTypes, h_nbfp, + 32, 1, 3.12341f, 138.935f); + // Warming-up nbnxmKernelTest<<>>( - a_xq, - a_f, - shiftVec, - fShift, - cj4, - sci, - excl, - atomTypes, - nbfp, + d_xq, + d_f, + d_shiftVec, + d_fShift, + d_cj4, + d_sci, + d_excl, + d_atomTypes, + d_nbfp, 32, 1, 3.12341, 138.935, 0); - cudaDeviceSynchronize(); + cudaDeviceSynchronize(); auto start = std::chrono::steady_clock::now(); for (int i = 0; i < repeat; ++i) { nbnxmKernelTest<<>>( - a_xq, - a_f, - shiftVec, - fShift, - cj4, - sci, - excl, - atomTypes, - nbfp, + d_xq, + d_f, + d_shiftVec, + d_fShift, + d_cj4, + d_sci, + d_excl, + d_atomTypes, + d_nbfp, 32, 1, 3.12341, @@ -502,65 +551,42 @@ int main(int argc, char* argv[]) { auto time = std::chrono::duration_cast(end - start).count(); printf("Average kernel execution time (w/o shift): %f (us)\n", (time * 1e-3f) / repeat); -#ifdef DEBUG - float f0 = 0, f1 = 0, f2 = 0; - for (int i = 0; i < NUM_ATOMS; ++i) { - f0 += a_f[i][0]; - f1 += a_f[i][1]; - f2 += a_f[i][2]; - } - printf("Checksum (a_f): %f %f %f\n", f0, f1, f2); + // Pull the accumulated forces back to host for validation against the + // CPU reference (d_f / d_fShift are not host-readable any more). + cudaMemcpy(h_f, d_f, sizeof(Float3) * NUM_ATOMS, cudaMemcpyDeviceToHost); + cudaMemcpy(h_fShift, d_fShift, sizeof(Float3) * 45, cudaMemcpyDeviceToHost); - f0 = 0, f1 = 0, f2 = 0; - for (int i = 0; i < 45; ++i) { - f0 += fShift[i][0]; - f1 += fShift[i][1]; - f2 += fShift[i][2]; - } - printf("Checksum (fShift): %f %f %f\n", f0, f1, f2); -#endif + ref.computeDelta(/*calcShift=*/false); + ref.validate(h_f, h_fShift, /*launchCount=*/repeat + 1, + 1.0f, 0.5f, 0.25f, 1.0f, 0.5f, 0.25f, + /*absTol=*/1e-3f, "w/o shift"); + // Reset the force accumulators (host copies), then push the reset back + // to device before the second (w/ shift) timed block, exactly mirroring + // what the original cudaMallocManaged version did in-place. for (int i = 0; i < NUM_ATOMS; ++i) { - a_xq[i] = make_float4(1.0f, 0.5f, 0.25f, 0.125f); - } - for (int i = 0; i < NUM_ATOMS; ++i) { - a_f[i] = Float3(1.0f, 0.5f, 0.25f); + h_f[i] = Float3(1.0f, 0.5f, 0.25f); } for (int i = 0; i < 45; ++i) { - shiftVec[i] = Float3(1.0f, 0.5f, 0.25f); - } - for (int i = 0; i < 45; ++i) { - fShift[i] = Float3(1.0f, 0.5f, 0.25f); - } - for (int i = 0; i < 56881; ++i) { - cj4[i] = get_cj4(i); - } - for (int i = 0; i < 4806; ++i) { - sci[i] = get_sci(i); - } - for (int i = 0; i < 19205; ++i) { - excl[i] = get_excl(i); - } - for (int i = 0; i < NUM_ATOMS; ++i) { - atomTypes[i] = (i % 2); - } - for (int i = 0; i < 1024; ++i) { - nbfp[i] = make_float2(0.5f, 0.25f); + h_fShift[i] = Float3(1.0f, 0.5f, 0.25f); } + cudaMemcpy(d_f, h_f, sizeof(Float3) * NUM_ATOMS, cudaMemcpyHostToDevice); + cudaMemcpy(d_fShift, h_fShift, sizeof(Float3) * 45, cudaMemcpyHostToDevice); + cudaDeviceSynchronize(); start = std::chrono::steady_clock::now(); for (int i = 0; i < repeat; ++i) { nbnxmKernelTest<<>>( - a_xq, - a_f, - shiftVec, - fShift, - cj4, - sci, - excl, - atomTypes, - nbfp, + d_xq, + d_f, + d_shiftVec, + d_fShift, + d_cj4, + d_sci, + d_excl, + d_atomTypes, + d_nbfp, 32, 1, 3.12341, @@ -573,34 +599,35 @@ int main(int argc, char* argv[]) { time = std::chrono::duration_cast(end - start).count(); printf("Average kernel execution time (w/ shift): %f (us)\n", (time * 1e-3f) / repeat); -#ifdef DEBUG - f0 = 0, f1 = 0, f2 = 0; - for (int i = 0; i < NUM_ATOMS; ++i) { - f0 += a_f[i][0]; - f1 += a_f[i][1]; - f2 += a_f[i][2]; - } - printf("Checksum (a_f): %f %f %f\n", f0, f1, f2); - - f0 = 0, f1 = 0, f2 = 0; - for (int i = 0; i < 45; ++i) { - f0 += fShift[i][0]; - f1 += fShift[i][1]; - f2 += fShift[i][2]; - } - printf("Checksum (fShift): %f %f %f\n", f0, f1, f2); -#endif - - cudaFree(nbfp); - cudaFree(atomTypes); - cudaFree(excl); - cudaFree(sci); - cudaFree(cj4); - cudaFree(fShift); - cudaFree(shiftVec); - cudaFree(a_f); - cudaFree(a_xq); + cudaMemcpy(h_f, d_f, sizeof(Float3) * NUM_ATOMS, cudaMemcpyDeviceToHost); + cudaMemcpy(h_fShift, d_fShift, sizeof(Float3) * 45, cudaMemcpyDeviceToHost); + + ref.computeDelta(/*calcShift=*/true); + // Second block's `repeat` launches all pass calcShift=1; there is + // no extra warm-up launch before this block. + ref.validate(h_f, h_fShift, /*launchCount=*/repeat, + 1.0f, 0.5f, 0.25f, 1.0f, 0.5f, 0.25f, + /*absTol=*/1e-3f, "w/ shift"); + + cudaFree(d_nbfp); + cudaFree(d_atomTypes); + cudaFree(d_excl); + cudaFree(d_sci); + cudaFree(d_cj4); + cudaFree(d_fShift); + cudaFree(d_shiftVec); + cudaFree(d_f); + cudaFree(d_xq); + + delete[] h_nbfp; + delete[] h_atomTypes; + delete[] h_excl; + delete[] h_sci; + delete[] h_cj4; + delete[] h_fShift; + delete[] h_shiftVec; + delete[] h_f; + delete[] h_xq; return 0; } - diff --git a/src/nbnxm-cuda/reference.h b/src/nbnxm-cuda/reference.h new file mode 100644 index 000000000..6dc984a63 --- /dev/null +++ b/src/nbnxm-cuda/reference.h @@ -0,0 +1,234 @@ +#pragma once + +#include +#include +#include +#include + +static inline float pmeCorrF_host(const float z2) +{ + constexpr float FN6 = -1.7357322914161492954e-8F; + constexpr float FN5 = 1.4703624142580877519e-6F; + constexpr float FN4 = -0.000053401640219807709149F; + constexpr float FN3 = 0.0010054721316683106153F; + constexpr float FN2 = -0.019278317264888380590F; + constexpr float FN1 = 0.069670166153766424023F; + constexpr float FN0 = -0.75225204789749321333F; + + constexpr float FD4 = 0.0011193462567257629232F; + constexpr float FD3 = 0.014866955030185295499F; + constexpr float FD2 = 0.11583842382862377919F; + constexpr float FD1 = 0.50736591960530292870F; + constexpr float FD0 = 1.0F; + + const float z4 = z2 * z2; + float polyFD0 = FD4 * z4 + FD2; + const float polyFD1 = FD3 * z4 + FD1; + polyFD0 = polyFD0 * z4 + FD0; + polyFD0 = polyFD1 * z2 + polyFD0; + polyFD0 = 1.0F / polyFD0; + + float polyFN0 = FN6 * z4 + FN4; + float polyFN1 = FN5 * z4 + FN3; + polyFN0 = polyFN0 * z4 + FN2; + polyFN1 = polyFN1 * z4 + FN1; + polyFN0 = polyFN0 * z4 + FN0; + polyFN0 = polyFN1 * z2 + polyFN0; + + return polyFN0 * polyFD0; +} + +struct HostVec3 { float x = 0.f, y = 0.f, z = 0.f; }; + +class NbnxmReference +{ +public: + NbnxmReference(const Float4* xq, + const Float3* shiftVec, + const nbnxn_cj4_t* cj4, + const nbnxn_sci_t* sci, + const nbnxn_excl_t* excl, + const int* atomTypes, + const Float2* nbfp, + int numTypes, + float rCoulombSq, + float ewaldBeta, + float epsFac) + : xq_(xq), shiftVec_(shiftVec), cj4_(cj4), sci_(sci), excl_(excl), + atomTypes_(atomTypes), nbfp_(nbfp), numTypes_(numTypes), + rCoulombSq_(rCoulombSq), ewaldBeta_(ewaldBeta), epsFac_(epsFac), + deltaF_(NUM_ATOMS), deltaFShift_(c_numIvecs) + {} + + // Computes the force / shift-force contribution of a SINGLE kernel + // launch. + void computeDelta(bool calcShift) + { + std::fill(deltaF_.begin(), deltaF_.end(), HostVec3{}); + std::fill(deltaFShift_.begin(), deltaFShift_.end(), HostVec3{}); + + const float beta2 = ewaldBeta_ * ewaldBeta_; + const float beta3 = ewaldBeta_ * ewaldBeta_ * ewaldBeta_; + constexpr unsigned superClMask = + ((1U << c_nbnxnGpuNumClusterPerSupercluster) - 1U); + constexpr int prunedClusterPairSize = c_clSize * c_splitClSize; + + for (int bidx = 0; bidx < grid_z; ++bidx) { + const nbnxn_sci_t nbSci = sci_[bidx]; + const int sciIdx = nbSci.sci; + const int cij4Start = nbSci.cj4_ind_start; + const int cij4End = nbSci.cj4_ind_end; + const int shift = nbSci.shift; + const bool doCalcShift = (calcShift && shift != c_centralShiftIndex); + + for (int j4 = cij4Start; j4 < cij4End; ++j4) { + for (int jm = 0; jm < c_nbnxnGpuJgroupSize; ++jm) { + const int cj = cj4_[j4].cj[jm]; + + for (int tidxj = 0; tidxj < c_clSize; ++tidxj) { + const int aj = cj * c_clSize + tidxj; + const int imeiIdx = tidxj / c_splitClSize; + const unsigned imask = cj4_[j4].imei[imeiIdx].imask; + + // Fast skip: none of this jm's bits are set for this imeiIdx. + if (!(imask & (superClMask << (jm * c_nbnxnGpuNumClusterPerSupercluster)))) + continue; + + const int wexclIdx = cj4_[j4].imei[imeiIdx].excl_ind; + + const Float4 xqjRaw = xq_[aj]; + const float xj = xqjRaw.x, yj = xqjRaw.y, zj = xqjRaw.z, qj = xqjRaw.w; + + for (int tidxi = 0; tidxi < c_clSize; ++tidxi) { + const int tidx = tidxi + tidxj * c_clSize; + const unsigned wexcl = + excl_[wexclIdx].pair[tidx & (prunedClusterPairSize - 1)]; + const bool nonSelf = + !(shift == c_centralShiftIndex && tidxj <= (unsigned)tidxi); + + for (int i = 0; i < c_nbnxnGpuNumClusterPerSupercluster; ++i) { + const unsigned maskJI = + (1U << (jm * c_nbnxnGpuNumClusterPerSupercluster + i)); + if (!(imask & maskJI)) continue; + + const int ci = sciIdx * c_nbnxnGpuNumClusterPerSupercluster + i; + const int ai = ci * c_clSize + tidxi; + + const Float3 shiftV = shiftVec_[shift]; + const Float4 xqiRaw = xq_[ai]; + const float xi = xqiRaw.x + shiftV[0]; + const float yi = xqiRaw.y + shiftV[1]; + const float zi = xqiRaw.z + shiftV[2]; + const float qi = xqiRaw.w * epsFac_; + + const float rvx = xi - xj, rvy = yi - yj, rvz = zi - zj; + float r2 = rvx * rvx + rvy * rvy + rvz * rvz; + + const float pairExclMask = (wexcl & maskJI) ? 1.0f : 0.0f; + const bool notExcluded = nonSelf || (ci != cj); + + if (!(r2 < rCoulombSq_) || !notExcluded) continue; + + const int atomTypeI = atomTypes_[ai]; + const int atomTypeJ = atomTypes_[aj]; + const Float2 c6c12 = nbfp_[numTypes_ * atomTypeI + atomTypeJ]; + const float c6 = c6c12.x, c12 = c6c12.y; + + r2 = std::max(r2, c_nbnxnMinDistanceSquared); + const float rInv = 1.0f / std::sqrt(r2); + const float r2Inv = rInv * rInv; + float r6Inv = r2Inv * r2Inv * r2Inv; + r6Inv *= pairExclMask; + float fInvR = r6Inv * (c12 * r6Inv - c6) * r2Inv; + fInvR += qi * qj * + (pairExclMask * r2Inv * rInv + pmeCorrF_host(beta2 * r2) * beta3); + + const float fx = rvx * fInvR, fy = rvy * fInvR, fz = rvz * fInvR; + + deltaF_[ai].x += fx; deltaF_[ai].y += fy; deltaF_[ai].z += fz; + deltaF_[aj].x -= fx; deltaF_[aj].y -= fy; deltaF_[aj].z -= fz; + + if (doCalcShift) { + deltaFShift_[shift].x += fx; + deltaFShift_[shift].y += fy; + deltaFShift_[shift].z += fz; + } + } + } + } + } + } + } + } + + // Compares initValue + launchCount * delta against the GPU buffers using + // a combined absolute+relative tolerance: + // |got - expected| <= absTol + relTol * max(|got|, |expected|) + bool validate(const Float3* gpu_f, const Float3* gpu_fShift, + int launchCount, + float initFx, float initFy, float initFz, + float initFShiftX, float initFShiftY, float initFShiftZ, + float absTol, + const char* label, + float relTol = 1e-3f, + int maxReportedMismatches = 10) const + { + bool ok = true; + int reported = 0; + + auto withinTol = [&](float got, float expected) { + const float diff = std::fabs(got - expected); + const float bound = absTol + relTol * std::max(std::fabs(got), std::fabs(expected)); + return diff <= bound; + }; + + for (int i = 0; i < NUM_ATOMS; ++i) { + const float ex = initFx + launchCount * deltaF_[i].x; + const float ey = initFy + launchCount * deltaF_[i].y; + const float ez = initFz + launchCount * deltaF_[i].z; + const Float3 gv = gpu_f[i]; + if (!withinTol(gv[0], ex) || !withinTol(gv[1], ey) || !withinTol(gv[2], ez)) { + ok = false; + if (reported < maxReportedMismatches) { + printf("[%s] a_f[%d] mismatch: got (%f, %f, %f) expected (%f, %f, %f)\n", + label, i, gv[0], gv[1], gv[2], ex, ey, ez); + ++reported; + } + } + } + + for (int s = 0; s < c_numIvecs; ++s) { + const float ex = initFShiftX + launchCount * deltaFShift_[s].x; + const float ey = initFShiftY + launchCount * deltaFShift_[s].y; + const float ez = initFShiftZ + launchCount * deltaFShift_[s].z; + const Float3 gv = gpu_fShift[s]; + if (!withinTol(gv[0], ex) || !withinTol(gv[1], ey) || !withinTol(gv[2], ez)) { + ok = false; + if (reported < maxReportedMismatches) { + printf("[%s] fShift[%d] mismatch: got (%f, %f, %f) expected (%f, %f, %f)\n", + label, s, gv[0], gv[1], gv[2], ex, ey, ez); + ++reported; + } + } + } + + printf("[%s] validation: %s\n", label, ok ? "PASS" : "FAIL"); + return ok; + } + +private: + const Float4* xq_; + const Float3* shiftVec_; + const nbnxn_cj4_t* cj4_; + const nbnxn_sci_t* sci_; + const nbnxn_excl_t* excl_; + const int* atomTypes_; + const Float2* nbfp_; + int numTypes_; + float rCoulombSq_; + float ewaldBeta_; + float epsFac_; + + std::vector deltaF_; + std::vector deltaFShift_; +}; diff --git a/src/nbnxm-hip/Makefile b/src/nbnxm-hip/Makefile index 64398b6dd..e8f3e77ac 100644 --- a/src/nbnxm-hip/Makefile +++ b/src/nbnxm-hip/Makefile @@ -45,7 +45,8 @@ endif $(program): $(obj) Makefile $(CC) $(CFLAGS) $(obj) -o $@ $(LDFLAGS) -%.o: %.cu ../nbnxm-cuda/constants.h ../nbnxm-cuda/vectypes.h Makefile +%.o: %.cu ../nbnxm-cuda/constants.h ../nbnxm-cuda/vectypes.h \ + ../nbnxm-cuda/reference.h Makefile $(CC) $(CFLAGS) -c $< -o $@ clean: diff --git a/src/nbnxm-hip/main.cu b/src/nbnxm-hip/main.cu index 087303ba4..4b60dc311 100644 --- a/src/nbnxm-hip/main.cu +++ b/src/nbnxm-hip/main.cu @@ -1,5 +1,6 @@ #include #include +#include #include #include "vectypes.h" @@ -8,6 +9,7 @@ typedef gmx::BasicVector Float3; typedef float4 Float4; #include "constants.h" +#include "reference.h" inline __device__ void operator+=(float4 &a, float4 b) { @@ -371,21 +373,16 @@ nbnxn_cj4_t get_cj4(int id) { value.cj[i] = i + id; } for (int i = 0; i < c_nbnxnGpuClusterpairSplit; ++i) { - value.imei[i].imask = 0U; - value.imei[i].excl_ind = 0; + value.imei[i].imask = 0xFFFFFFFFu; + value.imei[i].excl_ind = id % 19205; } return value; } -nbnxn_sci_t get_sci(int id) { - return {id, 0, 8 * id, 8 * id + 7}; -} - nbnxn_excl_t get_excl(int id) { nbnxn_excl_t value; - for (int i = 0; i < c_nbnxnGpuExclSize; ++i) { - value.pair[i] = 7; - } + std::mt19937 rng(id); + for (int i = 0; i < c_nbnxnGpuExclSize; ++i) value.pair[i] = rng(); return value; } @@ -399,92 +396,144 @@ int main(int argc, char* argv[]) { const dim3 blocks ( block_x, block_y, 1 ); const dim3 grids ( 1, 1, grid_z ); - Float4* a_xq; - hipMallocManaged(&a_xq, sizeof(Float4) * NUM_ATOMS); - - Float3* a_f; - hipMallocManaged(&a_f, sizeof(Float3) * NUM_ATOMS); - - Float3* shiftVec; - hipMallocManaged(&shiftVec, sizeof(Float3) * 45); - - Float3* fShift; - hipMallocManaged(&fShift, sizeof(Float3) * 45); - - nbnxn_cj4_t* cj4; - hipMallocManaged(&cj4, sizeof(nbnxn_cj4_t) * 56881); - - nbnxn_sci_t* sci; - hipMallocManaged(&sci, sizeof(nbnxn_sci_t) * 4806); - - nbnxn_excl_t* excl; - hipMallocManaged(&excl, sizeof(nbnxn_excl_t) * 19205); - - int* atomTypes; - hipMallocManaged(&atomTypes, sizeof(int) * NUM_ATOMS); - - Float2* nbfp; - hipMallocManaged(&nbfp, sizeof(Float2) * 1024); + // ----------------------------------------------------------------- + // Host buffers (h_ prefix) + // ----------------------------------------------------------------- + // These hold the "source of truth" data on the CPU. They are used to: + // (1) stage initial values into the device buffers via hipMemcpy, and + // (2) feed NbnxmReference (reference.h), which is host-only code and + // must be able to dereference xq/shiftVec/cj4/sci/excl/atomTypes/ + // nbfp directly -- it can no longer read the d_* buffers now that + // they are plain (non-managed) device memory. + // a_f / fShift do not need a persistent host mirror for staging (their + // initial values are trivial to re-write directly into h_f / h_fShift + // each time we reset them), but h_f / h_fShift ARE needed as the landing + // buffer for the periodic hipMemcpyDeviceToHost reads used to validate + // against the CPU reference. + Float4* h_xq = new Float4[NUM_ATOMS]; + Float3* h_f = new Float3[NUM_ATOMS]; + Float3* h_shiftVec = new Float3[45]; + Float3* h_fShift = new Float3[45]; + nbnxn_cj4_t* h_cj4 = new nbnxn_cj4_t[56881]; + nbnxn_sci_t* h_sci = new nbnxn_sci_t[4806]; + nbnxn_excl_t* h_excl = new nbnxn_excl_t[19205]; + int* h_atomTypes = new int[NUM_ATOMS]; + Float2* h_nbfp = new Float2[1024]; + + // ----------------------------------------------------------------- + // Device buffers (d_ prefix) + // ----------------------------------------------------------------- + // Plain hipMalloc (NOT hipMallocManaged): per PR #326's rationale, + // atomics on hipMallocManaged memory can be pathologically slow on + // unified-memory architectures because they require system-scope + // (CPU-GPU) coherency, even though the host side is never touched + // during the timed kernel loops here. nbnxmKernelTest() only + // atomicAdd()s into d_f / d_fShift and only reads the rest, so plain + // device memory + explicit H2D/D2H staging is strictly faster with + // identical results. + Float4* d_xq; + Float3* d_f; + Float3* d_shiftVec; + Float3* d_fShift; + nbnxn_cj4_t* d_cj4; + nbnxn_sci_t* d_sci; + nbnxn_excl_t* d_excl; + int* d_atomTypes; + Float2* d_nbfp; + + hipMalloc(&d_xq, sizeof(Float4) * NUM_ATOMS); + hipMalloc(&d_f, sizeof(Float3) * NUM_ATOMS); + hipMalloc(&d_shiftVec, sizeof(Float3) * 45); + hipMalloc(&d_fShift, sizeof(Float3) * 45); + hipMalloc(&d_cj4, sizeof(nbnxn_cj4_t) * 56881); + hipMalloc(&d_sci, sizeof(nbnxn_sci_t) * 4806); + hipMalloc(&d_excl, sizeof(nbnxn_excl_t) * 19205); + hipMalloc(&d_atomTypes, sizeof(int) * NUM_ATOMS); + hipMalloc(&d_nbfp, sizeof(Float2) * 1024); + + // ----------------------------------------------------------------- + // Populate host buffers with the same deterministic/synthetic data as + // before, then copy H2D once (these inputs are re-used, unchanged, for + // every launch in both the "w/o shift" and "w/ shift" blocks below). + // ----------------------------------------------------------------- + std::mt19937 rng(1337); + std::uniform_real_distribution posDist(-20.f, 20.f); + std::uniform_real_distribution qDist(-1.f, 1.f); for (int i = 0; i < NUM_ATOMS; ++i) { - a_xq[i] = make_float4(1.0f, 0.5f, 0.25f, 0.125f); + h_xq[i] = make_float4(posDist(rng), posDist(rng), posDist(rng), qDist(rng)); } for (int i = 0; i < NUM_ATOMS; ++i) { - a_f[i] = Float3(1.0f, 0.5f, 0.25f); + h_f[i] = Float3(1.0f, 0.5f, 0.25f); } for (int i = 0; i < 45; ++i) { - shiftVec[i] = Float3(1.0f, 0.5f, 0.25f); + h_shiftVec[i] = Float3(posDist(rng)*0.1f, posDist(rng)*0.1f, posDist(rng)*0.1f); } for (int i = 0; i < 45; ++i) { - fShift[i] = Float3(1.0f, 0.5f, 0.25f); + h_fShift[i] = Float3(1.0f, 0.5f, 0.25f); } for (int i = 0; i < 56881; ++i) { - cj4[i] = get_cj4(i); + h_cj4[i] = get_cj4(i % 200); } for (int i = 0; i < 4806; ++i) { - sci[i] = get_sci(i); + h_sci[i] = {i % 400, i % c_numIvecs, (2*i) % 200, (2*i) % 200 + 1}; } for (int i = 0; i < 19205; ++i) { - excl[i] = get_excl(i); + h_excl[i] = get_excl(i); } for (int i = 0; i < NUM_ATOMS; ++i) { - atomTypes[i] = (i % 2); + h_atomTypes[i] = (i % 32); } for (int i = 0; i < 1024; ++i) { - nbfp[i] = make_float2(0.5f, 0.25f); + h_nbfp[i] = make_float2(0.5f, 0.25f); } + hipMemcpy(d_xq, h_xq, sizeof(Float4) * NUM_ATOMS, hipMemcpyHostToDevice); + hipMemcpy(d_f, h_f, sizeof(Float3) * NUM_ATOMS, hipMemcpyHostToDevice); + hipMemcpy(d_shiftVec, h_shiftVec, sizeof(Float3) * 45, hipMemcpyHostToDevice); + hipMemcpy(d_fShift, h_fShift, sizeof(Float3) * 45, hipMemcpyHostToDevice); + hipMemcpy(d_cj4, h_cj4, sizeof(nbnxn_cj4_t) * 56881, hipMemcpyHostToDevice); + hipMemcpy(d_sci, h_sci, sizeof(nbnxn_sci_t) * 4806, hipMemcpyHostToDevice); + hipMemcpy(d_excl, h_excl, sizeof(nbnxn_excl_t) * 19205, hipMemcpyHostToDevice); + hipMemcpy(d_atomTypes, h_atomTypes, sizeof(int) * NUM_ATOMS, hipMemcpyHostToDevice); + hipMemcpy(d_nbfp, h_nbfp, sizeof(Float2) * 1024, hipMemcpyHostToDevice); + + // NbnxmReference reads its inputs directly from host memory -- pass the + // h_* buffers, not the d_* ones (it cannot dereference device pointers). + NbnxmReference ref(h_xq, h_shiftVec, h_cj4, h_sci, h_excl, h_atomTypes, h_nbfp, + 32, 1, 3.12341f, 138.935f); + // Warming-up - hipLaunchKernelGGL(nbnxmKernelTest, grids, blocks, 0, 0, - a_xq, - a_f, - shiftVec, - fShift, - cj4, - sci, - excl, - atomTypes, - nbfp, + nbnxmKernelTest<<>>( + d_xq, + d_f, + d_shiftVec, + d_fShift, + d_cj4, + d_sci, + d_excl, + d_atomTypes, + d_nbfp, 32, 1, 3.12341, 138.935, 0); - hipDeviceSynchronize(); + hipDeviceSynchronize(); auto start = std::chrono::steady_clock::now(); for (int i = 0; i < repeat; ++i) { - hipLaunchKernelGGL(nbnxmKernelTest, grids, blocks, 0, 0, - a_xq, - a_f, - shiftVec, - fShift, - cj4, - sci, - excl, - atomTypes, - nbfp, + nbnxmKernelTest<<>>( + d_xq, + d_f, + d_shiftVec, + d_fShift, + d_cj4, + d_sci, + d_excl, + d_atomTypes, + d_nbfp, 32, 1, 3.12341, @@ -497,65 +546,42 @@ int main(int argc, char* argv[]) { auto time = std::chrono::duration_cast(end - start).count(); printf("Average kernel execution time (w/o shift): %f (us)\n", (time * 1e-3f) / repeat); -#ifdef DEBUG - float f0 = 0, f1 = 0, f2 = 0; - for (int i = 0; i < NUM_ATOMS; ++i) { - f0 += a_f[i][0]; - f1 += a_f[i][1]; - f2 += a_f[i][2]; - } - printf("Checksum (a_f): %f %f %f\n", f0, f1, f2); + // Pull the accumulated forces back to host for validation against the + // CPU reference (d_f / d_fShift are not host-readable any more). + hipMemcpy(h_f, d_f, sizeof(Float3) * NUM_ATOMS, hipMemcpyDeviceToHost); + hipMemcpy(h_fShift, d_fShift, sizeof(Float3) * 45, hipMemcpyDeviceToHost); - f0 = 0, f1 = 0, f2 = 0; - for (int i = 0; i < 45; ++i) { - f0 += fShift[i][0]; - f1 += fShift[i][1]; - f2 += fShift[i][2]; - } - printf("Checksum (fShift): %f %f %f\n", f0, f1, f2); -#endif + ref.computeDelta(/*calcShift=*/false); + ref.validate(h_f, h_fShift, /*launchCount=*/repeat + 1, + 1.0f, 0.5f, 0.25f, 1.0f, 0.5f, 0.25f, + /*absTol=*/1e-3f, "w/o shift"); + // Reset the force accumulators (host copies), then push the reset back + // to device before the second (w/ shift) timed block, exactly mirroring + // what the original hipMallocManaged version did in-place. for (int i = 0; i < NUM_ATOMS; ++i) { - a_xq[i] = make_float4(1.0f, 0.5f, 0.25f, 0.125f); - } - for (int i = 0; i < NUM_ATOMS; ++i) { - a_f[i] = Float3(1.0f, 0.5f, 0.25f); - } - for (int i = 0; i < 45; ++i) { - shiftVec[i] = Float3(1.0f, 0.5f, 0.25f); + h_f[i] = Float3(1.0f, 0.5f, 0.25f); } for (int i = 0; i < 45; ++i) { - fShift[i] = Float3(1.0f, 0.5f, 0.25f); - } - for (int i = 0; i < 56881; ++i) { - cj4[i] = get_cj4(i); - } - for (int i = 0; i < 4806; ++i) { - sci[i] = get_sci(i); - } - for (int i = 0; i < 19205; ++i) { - excl[i] = get_excl(i); - } - for (int i = 0; i < NUM_ATOMS; ++i) { - atomTypes[i] = (i % 2); - } - for (int i = 0; i < 1024; ++i) { - nbfp[i] = make_float2(0.5f, 0.25f); + h_fShift[i] = Float3(1.0f, 0.5f, 0.25f); } + hipMemcpy(d_f, h_f, sizeof(Float3) * NUM_ATOMS, hipMemcpyHostToDevice); + hipMemcpy(d_fShift, h_fShift, sizeof(Float3) * 45, hipMemcpyHostToDevice); + hipDeviceSynchronize(); start = std::chrono::steady_clock::now(); for (int i = 0; i < repeat; ++i) { - hipLaunchKernelGGL(nbnxmKernelTest, grids, blocks, 0, 0, - a_xq, - a_f, - shiftVec, - fShift, - cj4, - sci, - excl, - atomTypes, - nbfp, + nbnxmKernelTest<<>>( + d_xq, + d_f, + d_shiftVec, + d_fShift, + d_cj4, + d_sci, + d_excl, + d_atomTypes, + d_nbfp, 32, 1, 3.12341, @@ -568,33 +594,35 @@ int main(int argc, char* argv[]) { time = std::chrono::duration_cast(end - start).count(); printf("Average kernel execution time (w/ shift): %f (us)\n", (time * 1e-3f) / repeat); -#ifdef DEBUG - f0 = 0, f1 = 0, f2 = 0; - for (int i = 0; i < NUM_ATOMS; ++i) { - f0 += a_f[i][0]; - f1 += a_f[i][1]; - f2 += a_f[i][2]; - } - printf("Checksum (a_f): %f %f %f\n", f0, f1, f2); - - f0 = 0, f1 = 0, f2 = 0; - for (int i = 0; i < 45; ++i) { - f0 += fShift[i][0]; - f1 += fShift[i][1]; - f2 += fShift[i][2]; - } - printf("Checksum (fShift): %f %f %f\n", f0, f1, f2); -#endif - - hipFree(nbfp); - hipFree(atomTypes); - hipFree(excl); - hipFree(sci); - hipFree(cj4); - hipFree(fShift); - hipFree(shiftVec); - hipFree(a_f); - hipFree(a_xq); + hipMemcpy(h_f, d_f, sizeof(Float3) * NUM_ATOMS, hipMemcpyDeviceToHost); + hipMemcpy(h_fShift, d_fShift, sizeof(Float3) * 45, hipMemcpyDeviceToHost); + + ref.computeDelta(/*calcShift=*/true); + // Second block's `repeat` launches all pass calcShift=1; there is + // no extra warm-up launch before this block. + ref.validate(h_f, h_fShift, /*launchCount=*/repeat, + 1.0f, 0.5f, 0.25f, 1.0f, 0.5f, 0.25f, + /*absTol=*/1e-3f, "w/ shift"); + + hipFree(d_nbfp); + hipFree(d_atomTypes); + hipFree(d_excl); + hipFree(d_sci); + hipFree(d_cj4); + hipFree(d_fShift); + hipFree(d_shiftVec); + hipFree(d_f); + hipFree(d_xq); + + delete[] h_nbfp; + delete[] h_atomTypes; + delete[] h_excl; + delete[] h_sci; + delete[] h_cj4; + delete[] h_fShift; + delete[] h_shiftVec; + delete[] h_f; + delete[] h_xq; return 0; } diff --git a/src/nbnxm-sycl/COPYRIGHT b/src/nbnxm-sycl/COPYRIGHT new file mode 100644 index 000000000..12f63792b --- /dev/null +++ b/src/nbnxm-sycl/COPYRIGHT @@ -0,0 +1,35 @@ +This file is part of the GROMACS molecular simulation package. + +Copyright (c) 1991-2000, University of Groningen, The Netherlands. +Copyright (c) 2001-2004, The GROMACS development team. +Copyright (c) 2013,2014,2016,2017,2018 by the GROMACS development team. +Copyright (c) 2019,2020,2021, 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. + +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. + +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. + +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 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 research papers on the package. Check out http://www.gromacs.org. diff --git a/src/nbnxm-sycl/Makefile b/src/nbnxm-sycl/Makefile index 029c87474..f0c69fcfd 100644 --- a/src/nbnxm-sycl/Makefile +++ b/src/nbnxm-sycl/Makefile @@ -72,7 +72,7 @@ endif $(program): $(obj) $(CC) $(CFLAGS) $(obj) -o $@ $(LDFLAGS) -%.o: %.cpp ../nbnxm-cuda/constants.h +%.o: %.cpp ../nbnxm-cuda/constants.h reference.h $(CC) $(CFLAGS) -c $< -o $@ clean: diff --git a/src/nbnxm-sycl/main.cpp b/src/nbnxm-sycl/main.cpp index d9663c905..ad8c557c6 100644 --- a/src/nbnxm-sycl/main.cpp +++ b/src/nbnxm-sycl/main.cpp @@ -1,5 +1,6 @@ #include #include +#include #include #include "vectypes.h" @@ -8,6 +9,7 @@ typedef gmx::BasicVector Float3; typedef sycl::float4 Float4; #include "constants.h" +#include "reference.h" template static inline void atomicAdd(T& val, const T delta) @@ -207,7 +209,7 @@ auto nbnxmKernelTest( const float ewaldBeta, const float epsFac, const bool calcShift) { - return [=](sycl::nd_item<3> itemIdx) [[intel::reqd_sub_group_size(32)]] { + return [=](sycl::nd_item<3> itemIdx) { constexpr int prunedClusterPairSize = c_clSize * c_splitClSize; @@ -223,9 +225,10 @@ auto nbnxmKernelTest( c_nbnxnGpuNumClusterPerSupercluster * c_clSize * sizeof(Float4) + // sm_xq c_clSize * c_clSize * DIM * sizeof(float) + // sm_reductionBuffer c_nbnxnGpuNumClusterPerSupercluster * c_clSize * sizeof(int); // sm_atomTypeI - sycl::multi_ptr localPtr = - sycl::ext::oneapi::group_local_memory_for_overwrite(itemIdx.get_group()); - uint8_t* ptr = *localPtr; + + auto group = itemIdx.get_group(); + uint8_t *ptr = + *sycl::ext::oneapi::group_local_memory_for_overwrite(group); Float4* sm_xq = reinterpret_cast(ptr); ptr += c_nbnxnGpuNumClusterPerSupercluster * c_clSize * sizeof(Float4); @@ -385,21 +388,16 @@ nbnxn_cj4_t get_cj4(int id) { value.cj[i] = i + id; } for (int i = 0; i < c_nbnxnGpuClusterpairSplit; ++i) { - value.imei[i].imask = 0U; - value.imei[i].excl_ind = 0; + value.imei[i].imask = 0xFFFFFFFFu; + value.imei[i].excl_ind = id % 19205; } return value; } -nbnxn_sci_t get_sci(int id) { - return {id, 0, 8 * id, 8 * id + 7}; -} - nbnxn_excl_t get_excl(int id) { nbnxn_excl_t value; - for (int i = 0; i < c_nbnxnGpuExclSize; ++i) { - value.pair[i] = 7; - } + std::mt19937 rng(id); + for (int i = 0; i < c_nbnxnGpuExclSize; ++i) value.pair[i] = rng(); return value; } @@ -416,57 +414,96 @@ int main(int argc, char* argv[]) { const sycl::range<3> globalSize{ grid_z, block_y, block_x }; const sycl::nd_range<3> range{ globalSize, blockSize }; - Float4* a_xq = sycl::malloc_shared(NUM_ATOMS, q); - Float3* a_f = sycl::malloc_shared(NUM_ATOMS, q); - Float3* shiftVec = sycl::malloc_shared(45, q); - Float3* fShift = sycl::malloc_shared(45, q); - nbnxn_cj4_t* cj4 = sycl::malloc_shared(56881, q); - nbnxn_sci_t* sci = sycl::malloc_shared(4806, q); - nbnxn_excl_t* excl = sycl::malloc_shared(19205, q); - int* atomTypes = sycl::malloc_shared(NUM_ATOMS, q); - Float2* nbfp = sycl::malloc_shared(1024, q); + Float4* h_xq = new Float4[NUM_ATOMS]; + Float3* h_f = new Float3[NUM_ATOMS]; + Float3* h_shiftVec = new Float3[45]; + Float3* h_fShift = new Float3[45]; + nbnxn_cj4_t* h_cj4 = new nbnxn_cj4_t[56881]; + nbnxn_sci_t* h_sci = new nbnxn_sci_t[4806]; + nbnxn_excl_t* h_excl = new nbnxn_excl_t[19205]; + int* h_atomTypes = new int[NUM_ATOMS]; + Float2* h_nbfp = new Float2[1024]; + + Float4* d_xq; + Float3* d_f; + Float3* d_shiftVec; + Float3* d_fShift; + nbnxn_cj4_t* d_cj4; + nbnxn_sci_t* d_sci; + nbnxn_excl_t* d_excl; + int* d_atomTypes; + Float2* d_nbfp; + + d_xq = sycl::malloc_device(NUM_ATOMS, q); + d_f = sycl::malloc_device(NUM_ATOMS, q); + d_shiftVec = sycl::malloc_device(45, q); + d_fShift = sycl::malloc_device(45, q); + d_cj4 = sycl::malloc_device(56881, q); + d_sci = sycl::malloc_device(4806, q); + d_excl = sycl::malloc_device(19205, q); + d_atomTypes = sycl::malloc_device(NUM_ATOMS, q); + d_nbfp = sycl::malloc_device(1024, q); + + std::mt19937 rng(1337); + std::uniform_real_distribution posDist(-20.f, 20.f); + std::uniform_real_distribution qDist(-1.f, 1.f); for (int i = 0; i < NUM_ATOMS; ++i) { - a_xq[i] = Float4(1.0f, 0.5f, 0.25f, 0.125f); + h_xq[i] = Float4(posDist(rng), posDist(rng), posDist(rng), qDist(rng)); } for (int i = 0; i < NUM_ATOMS; ++i) { - a_f[i] = Float3(1.0f, 0.5f, 0.25f); + h_f[i] = Float3(1.0f, 0.5f, 0.25f); } for (int i = 0; i < 45; ++i) { - shiftVec[i] = Float3(1.0f, 0.5f, 0.25f); + h_shiftVec[i] = Float3(posDist(rng)*0.1f, posDist(rng)*0.1f, posDist(rng)*0.1f); } for (int i = 0; i < 45; ++i) { - fShift[i] = Float3(1.0f, 0.5f, 0.25f); + h_fShift[i] = Float3(1.0f, 0.5f, 0.25f); } for (int i = 0; i < 56881; ++i) { - cj4[i] = get_cj4(i); + h_cj4[i] = get_cj4(i % 200); } for (int i = 0; i < 4806; ++i) { - sci[i] = get_sci(i); + h_sci[i] = {i % 400, i % c_numIvecs, (2*i) % 200, (2*i) % 200 + 1}; } for (int i = 0; i < 19205; ++i) { - excl[i] = get_excl(i); + h_excl[i] = get_excl(i); } for (int i = 0; i < NUM_ATOMS; ++i) { - atomTypes[i] = (i % 2); + h_atomTypes[i] = (i % 32); } for (int i = 0; i < 1024; ++i) { - nbfp[i] = Float2(0.5f, 0.25f); + h_nbfp[i] = Float2(0.5f, 0.25f); } + q.memcpy(d_xq, h_xq, sizeof(Float4) * NUM_ATOMS); + q.memcpy(d_f, h_f, sizeof(Float3) * NUM_ATOMS); + q.memcpy(d_shiftVec, h_shiftVec, sizeof(Float3) * 45); + q.memcpy(d_fShift, h_fShift, sizeof(Float3) * 45); + q.memcpy(d_cj4, h_cj4, sizeof(nbnxn_cj4_t) * 56881); + q.memcpy(d_sci, h_sci, sizeof(nbnxn_sci_t) * 4806); + q.memcpy(d_excl, h_excl, sizeof(nbnxn_excl_t) * 19205); + q.memcpy(d_atomTypes, h_atomTypes, sizeof(int) * NUM_ATOMS); + q.memcpy(d_nbfp, h_nbfp, sizeof(Float2) * 1024); + + // NbnxmReference reads its inputs directly from host memory -- pass the + // h_* buffers, not the d_* ones (it cannot dereference device pointers). + NbnxmReference ref(h_xq, h_shiftVec, h_cj4, h_sci, h_excl, h_atomTypes, h_nbfp, + 32, 1, 3.12341f, 138.935f); + // Warming-up q.submit([&](sycl::handler& cgh) { auto kernel = nbnxmKernelTest( cgh, - a_xq, - a_f, - shiftVec, - fShift, - cj4, - sci, - excl, - atomTypes, - nbfp, + d_xq, + d_f, + d_shiftVec, + d_fShift, + d_cj4, + d_sci, + d_excl, + d_atomTypes, + d_nbfp, 32, 1, 3.12341, @@ -481,15 +518,15 @@ int main(int argc, char* argv[]) { q.submit([&](sycl::handler& cgh) { auto kernel = nbnxmKernelTest( cgh, - a_xq, - a_f, - shiftVec, - fShift, - cj4, - sci, - excl, - atomTypes, - nbfp, + d_xq, + d_f, + d_shiftVec, + d_fShift, + d_cj4, + d_sci, + d_excl, + d_atomTypes, + d_nbfp, 32, 1, 3.12341, @@ -504,51 +541,29 @@ int main(int argc, char* argv[]) { auto time = std::chrono::duration_cast(end - start).count(); printf("Average kernel execution time (w/o shift): %f (us)\n", (time * 1e-3f) / repeat); -#ifdef DEBUG - float f0 = 0, f1 = 0, f2 = 0; - for (int i = 0; i < NUM_ATOMS; ++i) { - f0 += a_f[i][0]; - f1 += a_f[i][1]; - f2 += a_f[i][2]; - } - printf("Checksum (a_f): %f %f %f\n", f0, f1, f2); + // Pull the accumulated forces back to host for validation against the + // CPU reference (d_f / d_fShift are device-only USM, not host-readable). + q.memcpy(h_f, d_f, sizeof(Float3) * NUM_ATOMS); + q.memcpy(h_fShift, d_fShift, sizeof(Float3) * 45); + q.wait(); - f0 = 0, f1 = 0, f2 = 0; - for (int i = 0; i < 45; ++i) { - f0 += fShift[i][0]; - f1 += fShift[i][1]; - f2 += fShift[i][2]; - } - printf("Checksum (fShift): %f %f %f\n", f0, f1, f2); -#endif + ref.computeDelta(/*calcShift=*/false); + ref.validate(h_f, h_fShift, /*launchCount=*/repeat + 1, + 1.0f, 0.5f, 0.25f, 1.0f, 0.5f, 0.25f, + /*absTol=*/1e-3f, "w/o shift"); + // Reset the force accumulators (host copies), then push the reset back + // to device before the second (w/ shift) timed block, exactly mirroring + // what the original sycl::malloc_shared version did in-place. for (int i = 0; i < NUM_ATOMS; ++i) { - a_xq[i] = Float4(1.0f, 0.5f, 0.25f, 0.125f); - } - for (int i = 0; i < NUM_ATOMS; ++i) { - a_f[i] = Float3(1.0f, 0.5f, 0.25f); - } - for (int i = 0; i < 45; ++i) { - shiftVec[i] = Float3(1.0f, 0.5f, 0.25f); + h_f[i] = Float3(1.0f, 0.5f, 0.25f); } for (int i = 0; i < 45; ++i) { - fShift[i] = Float3(1.0f, 0.5f, 0.25f); - } - for (int i = 0; i < 56881; ++i) { - cj4[i] = get_cj4(i); - } - for (int i = 0; i < 4806; ++i) { - sci[i] = get_sci(i); - } - for (int i = 0; i < 19205; ++i) { - excl[i] = get_excl(i); - } - for (int i = 0; i < NUM_ATOMS; ++i) { - atomTypes[i] = (i % 2); - } - for (int i = 0; i < 1024; ++i) { - nbfp[i] = Float2(0.5f, 0.25f); + h_fShift[i] = Float3(1.0f, 0.5f, 0.25f); } + q.memcpy(d_f, h_f, sizeof(Float3) * NUM_ATOMS); + q.memcpy(d_fShift, h_fShift, sizeof(Float3) * 45); + q.wait(); start = std::chrono::steady_clock::now(); @@ -556,15 +571,15 @@ int main(int argc, char* argv[]) { q.submit([&](sycl::handler& cgh) { auto kernel = nbnxmKernelTest( cgh, - a_xq, - a_f, - shiftVec, - fShift, - cj4, - sci, - excl, - atomTypes, - nbfp, + d_xq, + d_f, + d_shiftVec, + d_fShift, + d_cj4, + d_sci, + d_excl, + d_atomTypes, + d_nbfp, 32, 1, 3.12341, @@ -579,33 +594,36 @@ int main(int argc, char* argv[]) { time = std::chrono::duration_cast(end - start).count(); printf("Average kernel execution time (w/ shift): %f (us)\n", (time * 1e-3f) / repeat); -#ifdef DEBUG - f0 = 0, f1 = 0, f2 = 0; - for (int i = 0; i < NUM_ATOMS; ++i) { - f0 += a_f[i][0]; - f1 += a_f[i][1]; - f2 += a_f[i][2]; - } - printf("Checksum (a_f): %f %f %f\n", f0, f1, f2); + q.memcpy(h_f, d_f, sizeof(Float3) * NUM_ATOMS); + q.memcpy(h_fShift, d_fShift, sizeof(Float3) * 45); + q.wait(); - f0 = 0, f1 = 0, f2 = 0; - for (int i = 0; i < 45; ++i) { - f0 += fShift[i][0]; - f1 += fShift[i][1]; - f2 += fShift[i][2]; - } - printf("Checksum (fShift): %f %f %f\n", f0, f1, f2); -#endif - - sycl::free(nbfp, q); - sycl::free(atomTypes, q); - sycl::free(excl, q); - sycl::free(sci, q); - sycl::free(cj4, q); - sycl::free(fShift, q); - sycl::free(shiftVec, q); - sycl::free(a_f, q); - sycl::free(a_xq, q); + ref.computeDelta(/*calcShift=*/true); + // Second block's `repeat` launches all pass calcShift=1; there is + // no extra warm-up launch before this block. + ref.validate(h_f, h_fShift, /*launchCount=*/repeat, + 1.0f, 0.5f, 0.25f, 1.0f, 0.5f, 0.25f, + /*absTol=*/1e-3f, "w/ shift"); + + sycl::free(d_nbfp, q); + sycl::free(d_atomTypes, q); + sycl::free(d_excl, q); + sycl::free(d_sci, q); + sycl::free(d_cj4, q); + sycl::free(d_fShift, q); + sycl::free(d_shiftVec, q); + sycl::free(d_f, q); + sycl::free(d_xq, q); + + delete[] h_nbfp; + delete[] h_atomTypes; + delete[] h_excl; + delete[] h_sci; + delete[] h_cj4; + delete[] h_fShift; + delete[] h_shiftVec; + delete[] h_f; + delete[] h_xq; return 0; } diff --git a/src/nbnxm-sycl/reference.h b/src/nbnxm-sycl/reference.h new file mode 100644 index 000000000..a86b07264 --- /dev/null +++ b/src/nbnxm-sycl/reference.h @@ -0,0 +1,235 @@ +#pragma once + +#include +#include +#include +#include + +static inline float pmeCorrF_host(const float z2) +{ + constexpr float FN6 = -1.7357322914161492954e-8F; + constexpr float FN5 = 1.4703624142580877519e-6F; + constexpr float FN4 = -0.000053401640219807709149F; + constexpr float FN3 = 0.0010054721316683106153F; + constexpr float FN2 = -0.019278317264888380590F; + constexpr float FN1 = 0.069670166153766424023F; + constexpr float FN0 = -0.75225204789749321333F; + + constexpr float FD4 = 0.0011193462567257629232F; + constexpr float FD3 = 0.014866955030185295499F; + constexpr float FD2 = 0.11583842382862377919F; + constexpr float FD1 = 0.50736591960530292870F; + constexpr float FD0 = 1.0F; + + const float z4 = z2 * z2; + float polyFD0 = FD4 * z4 + FD2; + const float polyFD1 = FD3 * z4 + FD1; + polyFD0 = polyFD0 * z4 + FD0; + polyFD0 = polyFD1 * z2 + polyFD0; + polyFD0 = 1.0F / polyFD0; + + float polyFN0 = FN6 * z4 + FN4; + float polyFN1 = FN5 * z4 + FN3; + polyFN0 = polyFN0 * z4 + FN2; + polyFN1 = polyFN1 * z4 + FN1; + polyFN0 = polyFN0 * z4 + FN0; + polyFN0 = polyFN1 * z2 + polyFN0; + + return polyFN0 * polyFD0; +} + +struct HostVec3 { float x = 0.f, y = 0.f, z = 0.f; }; + +class NbnxmReference +{ +public: + NbnxmReference(const Float4* xq, + const Float3* shiftVec, + const nbnxn_cj4_t* cj4, + const nbnxn_sci_t* sci, + const nbnxn_excl_t* excl, + const int* atomTypes, + const Float2* nbfp, + int numTypes, + float rCoulombSq, + float ewaldBeta, + float epsFac) + : xq_(xq), shiftVec_(shiftVec), cj4_(cj4), sci_(sci), excl_(excl), + atomTypes_(atomTypes), nbfp_(nbfp), numTypes_(numTypes), + rCoulombSq_(rCoulombSq), ewaldBeta_(ewaldBeta), epsFac_(epsFac), + deltaF_(NUM_ATOMS), deltaFShift_(c_numIvecs) + {} + + // Computes the force / shift-force contribution of a SINGLE kernel + // launch. + void computeDelta(bool calcShift) + { + std::fill(deltaF_.begin(), deltaF_.end(), HostVec3{}); + std::fill(deltaFShift_.begin(), deltaFShift_.end(), HostVec3{}); + + const float beta2 = ewaldBeta_ * ewaldBeta_; + const float beta3 = ewaldBeta_ * ewaldBeta_ * ewaldBeta_; + constexpr unsigned superClMask = + ((1U << c_nbnxnGpuNumClusterPerSupercluster) - 1U); + constexpr int prunedClusterPairSize = c_clSize * c_splitClSize; + + for (int bidx = 0; bidx < grid_z; ++bidx) { + const nbnxn_sci_t nbSci = sci_[bidx]; + const int sciIdx = nbSci.sci; + const int cij4Start = nbSci.cj4_ind_start; + const int cij4End = nbSci.cj4_ind_end; + const int shift = nbSci.shift; + const bool doCalcShift = (calcShift && shift != c_centralShiftIndex); + + for (int j4 = cij4Start; j4 < cij4End; ++j4) { + for (int jm = 0; jm < c_nbnxnGpuJgroupSize; ++jm) { + const int cj = cj4_[j4].cj[jm]; + + for (int tidxj = 0; tidxj < c_clSize; ++tidxj) { + const int aj = cj * c_clSize + tidxj; + const int imeiIdx = tidxj / c_splitClSize; + const unsigned imask = cj4_[j4].imei[imeiIdx].imask; + + // Fast skip: none of this jm's bits are set for this imeiIdx. + if (!(imask & (superClMask << (jm * c_nbnxnGpuNumClusterPerSupercluster)))) + continue; + + const int wexclIdx = cj4_[j4].imei[imeiIdx].excl_ind; + + const Float4 xqjRaw = xq_[aj]; + const float xj = xqjRaw[0], yj = xqjRaw[1], zj = xqjRaw[2], qj = xqjRaw[3]; + + for (int tidxi = 0; tidxi < c_clSize; ++tidxi) { + const int tidx = tidxi + tidxj * c_clSize; + const unsigned wexcl = + excl_[wexclIdx].pair[tidx & (prunedClusterPairSize - 1)]; + const bool nonSelf = + !(shift == c_centralShiftIndex && tidxj <= (unsigned)tidxi); + + for (int i = 0; i < c_nbnxnGpuNumClusterPerSupercluster; ++i) { + const unsigned maskJI = + (1U << (jm * c_nbnxnGpuNumClusterPerSupercluster + i)); + if (!(imask & maskJI)) continue; + + const int ci = sciIdx * c_nbnxnGpuNumClusterPerSupercluster + i; + const int ai = ci * c_clSize + tidxi; + + const Float3 shiftV = shiftVec_[shift]; + const Float4 xqiRaw = xq_[ai]; + const float xi = xqiRaw[0] + shiftV[0]; + const float yi = xqiRaw[1] + shiftV[1]; + const float zi = xqiRaw[2] + shiftV[2]; + const float qi = xqiRaw[3] * epsFac_; + + const float rvx = xi - xj, rvy = yi - yj, rvz = zi - zj; + float r2 = rvx * rvx + rvy * rvy + rvz * rvz; + + const float pairExclMask = (wexcl & maskJI) ? 1.0f : 0.0f; + const bool notExcluded = nonSelf || (ci != cj); + + if (!(r2 < rCoulombSq_) || !notExcluded) continue; + + const int atomTypeI = atomTypes_[ai]; + const int atomTypeJ = atomTypes_[aj]; + const Float2 c6c12 = nbfp_[numTypes_ * atomTypeI + atomTypeJ]; + const float c6 = c6c12[0], c12 = c6c12[1]; + + r2 = std::max(r2, c_nbnxnMinDistanceSquared); + const float rInv = 1.0f / std::sqrt(r2); + const float r2Inv = rInv * rInv; + float r6Inv = r2Inv * r2Inv * r2Inv; + r6Inv *= pairExclMask; + float fInvR = r6Inv * (c12 * r6Inv - c6) * r2Inv; + fInvR += qi * qj * + (pairExclMask * r2Inv * rInv + pmeCorrF_host(beta2 * r2) * beta3); + + const float fx = rvx * fInvR, fy = rvy * fInvR, fz = rvz * fInvR; + + deltaF_[ai].x += fx; deltaF_[ai].y += fy; deltaF_[ai].z += fz; + deltaF_[aj].x -= fx; deltaF_[aj].y -= fy; deltaF_[aj].z -= fz; + + if (doCalcShift) { + deltaFShift_[shift].x += fx; + deltaFShift_[shift].y += fy; + deltaFShift_[shift].z += fz; + } + } + } + } + } + } + } + } + + // Compares initValue + launchCount * delta against the GPU buffers using + // a combined absolute+relative tolerance: + // + // |got - expected| <= absTol + relTol * max(|got|, |expected|) + bool validate(const Float3* gpu_f, const Float3* gpu_fShift, + int launchCount, + float initFx, float initFy, float initFz, + float initFShiftX, float initFShiftY, float initFShiftZ, + float absTol, + const char* label, + float relTol = 1e-3f, + int maxReportedMismatches = 10) const + { + bool ok = true; + int reported = 0; + + auto withinTol = [&](float got, float expected) { + const float diff = std::fabs(got - expected); + const float bound = absTol + relTol * std::max(std::fabs(got), std::fabs(expected)); + return diff <= bound; + }; + + for (int i = 0; i < NUM_ATOMS; ++i) { + const float ex = initFx + launchCount * deltaF_[i].x; + const float ey = initFy + launchCount * deltaF_[i].y; + const float ez = initFz + launchCount * deltaF_[i].z; + const Float3 gv = gpu_f[i]; + if (!withinTol(gv[0], ex) || !withinTol(gv[1], ey) || !withinTol(gv[2], ez)) { + ok = false; + if (reported < maxReportedMismatches) { + printf("[%s] a_f[%d] mismatch: got (%f, %f, %f) expected (%f, %f, %f)\n", + label, i, gv[0], gv[1], gv[2], ex, ey, ez); + ++reported; + } + } + } + + for (int s = 0; s < c_numIvecs; ++s) { + const float ex = initFShiftX + launchCount * deltaFShift_[s].x; + const float ey = initFShiftY + launchCount * deltaFShift_[s].y; + const float ez = initFShiftZ + launchCount * deltaFShift_[s].z; + const Float3 gv = gpu_fShift[s]; + if (!withinTol(gv[0], ex) || !withinTol(gv[1], ey) || !withinTol(gv[2], ez)) { + ok = false; + if (reported < maxReportedMismatches) { + printf("[%s] fShift[%d] mismatch: got (%f, %f, %f) expected (%f, %f, %f)\n", + label, s, gv[0], gv[1], gv[2], ex, ey, ez); + ++reported; + } + } + } + + printf("[%s] validation: %s\n", label, ok ? "PASS" : "FAIL"); + return ok; + } + +private: + const Float4* xq_; + const Float3* shiftVec_; + const nbnxn_cj4_t* cj4_; + const nbnxn_sci_t* sci_; + const nbnxn_excl_t* excl_; + const int* atomTypes_; + const Float2* nbfp_; + int numTypes_; + float rCoulombSq_; + float ewaldBeta_; + float epsFac_; + + std::vector deltaF_; + std::vector deltaFShift_; +};