PME spline+spread CUDA kernel and unit tests
[alexxy/gromacs.git] / src / gromacs / ewald / tests / pmesplinespreadtest.cpp
index c751ef1ebfb59815f55737735b3a9e99d41dea18..e6352d0c41a177096b807fae7e83e830b05391b5 100644 (file)
@@ -103,96 +103,134 @@ class PmeSplineAndSpreadTest : public ::testing::TestWithParam<SplineAndSpreadIn
             inputRec.nkz         = gridSize[ZZ];
             inputRec.pme_order   = pmeOrder;
             inputRec.coulombtype = eelPME;
+            inputRec.epsilon_r   = 1.0;
 
             TestReferenceData                                      refData;
 
-            const std::map<CodePath, std::string>                  modesToTest   = {{CodePath::CPU, "CPU"}};
+            const std::map<CodePath, std::string>                  modesToTest   = {{CodePath::CPU, "CPU"},
+                                                                                    {CodePath::CUDA, "CUDA"}};
+
             const std::map<PmeSplineAndSpreadOptions, std::string> optionsToTest = {{PmeSplineAndSpreadOptions::SplineAndSpreadUnified, "spline computation and charge spreading (fused)"},
                                                                                     {PmeSplineAndSpreadOptions::SplineOnly, "spline computation"},
                                                                                     {PmeSplineAndSpreadOptions::SpreadOnly, "charge spreading"}};
+
+            // There is a subtle problem with multiple comparisons against same reference data:
+            // The subsequent (GPU) spreading runs at one point didn't actually copy the output grid into the proper buffer,
+            // but the reference data was already marked as checked (hasBeenChecked_) by the CPU run, so nothing failed.
+            // For now we will manually track that the count of the grid entries is the same on each run.
+            // This is just a hack for a single specific output though.
+            // What would be much better TODO is to split different codepaths into separate tests,
+            // while making them use the same reference files.
+            bool   gridValuesSizeAssigned = false;
+            size_t previousGridValuesSize;
+
             for (const auto &mode : modesToTest)
             {
-                for (const auto &option : optionsToTest)
+                const bool supportedInput = pmeSupportsInputForMode(&inputRec, mode.first);
+                if (!supportedInput)
                 {
-                    /* Describing the test uniquely in case it fails */
-
-                    SCOPED_TRACE(formatString("Testing %s with %s for PME grid size %d %d %d"
-                                              ", order %d, %zu atoms",
-                                              option.second.c_str(), mode.second.c_str(),
-                                              gridSize[XX], gridSize[YY], gridSize[ZZ],
-                                              pmeOrder,
-                                              atomCount));
+                    /* Testing the failure for the unsupported input */
+                    EXPECT_THROW(pmeInitAtoms(&inputRec, mode.first, nullptr, coordinates, charges, box), NotImplementedError);
+                    continue;
+                }
 
-                    /* Running the test */
+                const auto contextsToTest = pmeEnv->getHardwareContexts(mode.first);
+                for (const auto &context : contextsToTest)
+                {
+                    for (const auto &option : optionsToTest)
+                    {
+                        /* Describing the test uniquely in case it fails */
 
-                    PmeSafePointer pmeSafe = pmeInitAtoms(&inputRec, coordinates, charges, box);
+                        SCOPED_TRACE(formatString("Testing %s with %s %sfor PME grid size %d %d %d"
+                                                  ", order %d, %zu atoms",
+                                                  option.second.c_str(), mode.second.c_str(),
+                                                  context.getDescription().c_str(),
+                                                  gridSize[XX], gridSize[YY], gridSize[ZZ],
+                                                  pmeOrder,
+                                                  atomCount));
 
-                    const bool     computeSplines = (option.first == PmeSplineAndSpreadOptions::SplineOnly) || (option.first == PmeSplineAndSpreadOptions::SplineAndSpreadUnified);
-                    const bool     spreadCharges  = (option.first == PmeSplineAndSpreadOptions::SpreadOnly) || (option.first == PmeSplineAndSpreadOptions::SplineAndSpreadUnified);
+                        /* Running the test */
 
-                    if (!computeSplines)
-                    {
-                        // Here we should set up the results of the spline computation so that the spread can run.
-                        // What is lazy and works is running the separate spline so that it will set it up for us:
-                        pmePerformSplineAndSpread(pmeSafe.get(), mode.first, true, false);
-                        // We know that it is tested in another iteration.
-                        // TODO: Clean alternative: read and set the reference gridline indices, spline params
-                    }
+                        PmeSafePointer pmeSafe = pmeInitAtoms(&inputRec, mode.first, context.getDeviceInfo(), coordinates, charges, box);
 
-                    pmePerformSplineAndSpread(pmeSafe.get(), mode.first, computeSplines, spreadCharges);
+                        const bool     computeSplines = (option.first == PmeSplineAndSpreadOptions::SplineOnly) || (option.first == PmeSplineAndSpreadOptions::SplineAndSpreadUnified);
+                        const bool     spreadCharges  = (option.first == PmeSplineAndSpreadOptions::SpreadOnly) || (option.first == PmeSplineAndSpreadOptions::SplineAndSpreadUnified);
 
-                    /* Outputs correctness check */
-                    /* All tolerances were picked empirically for single precision on CPU */
+                        if (!computeSplines)
+                        {
+                            // Here we should set up the results of the spline computation so that the spread can run.
+                            // What is lazy and works is running the separate spline so that it will set it up for us:
+                            pmePerformSplineAndSpread(pmeSafe.get(), mode.first, true, false);
+                            // We know that it is tested in another iteration.
+                            // TODO: Clean alternative: read and set the reference gridline indices, spline params
+                        }
 
-                    TestReferenceChecker rootChecker(refData.rootChecker());
+                        pmePerformSplineAndSpread(pmeSafe.get(), mode.first, computeSplines, spreadCharges);
+                        pmeFinalizeTest(pmeSafe.get(), mode.first);
 
-                    const auto           maxGridSize              = std::max(std::max(gridSize[XX], gridSize[YY]), gridSize[ZZ]);
-                    const auto           ulpToleranceSplineValues = 2 * (pmeOrder - 2) * maxGridSize;
-                    /* 2 is empiric, the rest follows from the amount of operations */
+                        /* Outputs correctness check */
+                        /* All tolerances were picked empirically for single precision on CPU */
 
-                    if (computeSplines)
-                    {
-                        const char *dimString[] = { "X", "Y", "Z" };
+                        TestReferenceChecker rootChecker(refData.rootChecker());
 
-                        /* Spline values */
-                        SCOPED_TRACE(formatString("Testing spline values with tolerance of %ld", ulpToleranceSplineValues));
-                        TestReferenceChecker splineValuesChecker(rootChecker.checkCompound("Splines", "Values"));
-                        splineValuesChecker.setDefaultTolerance(relativeToleranceAsUlp(1.0, ulpToleranceSplineValues));
-                        for (int i = 0; i < DIM; i++)
-                        {
-                            auto splineValuesDim = pmeGetSplineData(pmeSafe.get(), mode.first, PmeSplineDataType::Values, i);
-                            splineValuesChecker.checkSequence(splineValuesDim.begin(), splineValuesDim.end(), dimString[i]);
-                        }
+                        const auto           maxGridSize              = std::max(std::max(gridSize[XX], gridSize[YY]), gridSize[ZZ]);
+                        const auto           ulpToleranceSplineValues = 2 * (pmeOrder - 2) * maxGridSize;
+                        /* 2 is empiric, the rest follows from the amount of operations */
 
-                        /* Spline derivatives */
-                        const auto ulpToleranceSplineDerivatives = 4 * ulpToleranceSplineValues;
-                        /* 4 is just a wild guess since the derivatives are deltas of neighbor spline values which could differ greatly */
-                        SCOPED_TRACE(formatString("Testing spline derivatives with tolerance of %ld", ulpToleranceSplineDerivatives));
-                        TestReferenceChecker splineDerivativesChecker(rootChecker.checkCompound("Splines", "Derivatives"));
-                        splineDerivativesChecker.setDefaultTolerance(relativeToleranceAsUlp(1.0, ulpToleranceSplineDerivatives));
-                        for (int i = 0; i < DIM; i++)
+                        if (computeSplines)
                         {
-                            auto splineDerivativesDim = pmeGetSplineData(pmeSafe.get(), mode.first, PmeSplineDataType::Derivatives, i);
-                            splineDerivativesChecker.checkSequence(splineDerivativesDim.begin(), splineDerivativesDim.end(), dimString[i]);
+                            const char *dimString[] = { "X", "Y", "Z" };
+
+                            /* Spline values */
+                            SCOPED_TRACE(formatString("Testing spline values with tolerance of %ld", ulpToleranceSplineValues));
+                            TestReferenceChecker splineValuesChecker(rootChecker.checkCompound("Splines", "Values"));
+                            splineValuesChecker.setDefaultTolerance(relativeToleranceAsUlp(1.0, ulpToleranceSplineValues));
+                            for (int i = 0; i < DIM; i++)
+                            {
+                                auto splineValuesDim = pmeGetSplineData(pmeSafe.get(), mode.first, PmeSplineDataType::Values, i);
+                                splineValuesChecker.checkSequence(splineValuesDim.begin(), splineValuesDim.end(), dimString[i]);
+                            }
+
+                            /* Spline derivatives */
+                            const auto ulpToleranceSplineDerivatives = 4 * ulpToleranceSplineValues;
+                            /* 4 is just a wild guess since the derivatives are deltas of neighbor spline values which could differ greatly */
+                            SCOPED_TRACE(formatString("Testing spline derivatives with tolerance of %ld", ulpToleranceSplineDerivatives));
+                            TestReferenceChecker splineDerivativesChecker(rootChecker.checkCompound("Splines", "Derivatives"));
+                            splineDerivativesChecker.setDefaultTolerance(relativeToleranceAsUlp(1.0, ulpToleranceSplineDerivatives));
+                            for (int i = 0; i < DIM; i++)
+                            {
+                                auto splineDerivativesDim = pmeGetSplineData(pmeSafe.get(), mode.first, PmeSplineDataType::Derivatives, i);
+                                splineDerivativesChecker.checkSequence(splineDerivativesDim.begin(), splineDerivativesDim.end(), dimString[i]);
+                            }
+
+                            /* Particle gridline indices */
+                            auto gridLineIndices = pmeGetGridlineIndices(pmeSafe.get(), mode.first);
+                            rootChecker.checkSequence(gridLineIndices.begin(), gridLineIndices.end(), "Gridline indices");
                         }
 
-                        /* Particle gridline indices */
-                        auto gridLineIndices = pmeGetGridlineIndices(pmeSafe.get(), mode.first);
-                        rootChecker.checkSequence(gridLineIndices.begin(), gridLineIndices.end(), "Gridline indices");
-                    }
-
-                    if (spreadCharges)
-                    {
-                        /* The wrapped grid */
-                        SparseRealGridValuesOutput nonZeroGridValues = pmeGetRealGrid(pmeSafe.get(), mode.first);
-                        TestReferenceChecker       gridValuesChecker(rootChecker.checkCompound("NonZeroGridValues", "RealSpaceGrid"));
-                        const auto                 ulpToleranceGrid = 2 * ulpToleranceSplineValues * (int)(ceil(sqrt(atomCount)));
-                        /* 2 is empiric; sqrt(atomCount) assumes all the input charges may spread onto the same cell */
-                        SCOPED_TRACE(formatString("Testing grid values with tolerance of %ld", ulpToleranceGrid));
-                        gridValuesChecker.setDefaultTolerance(relativeToleranceAsUlp(1.0, ulpToleranceGrid));
-                        for (const auto &point : nonZeroGridValues)
+                        if (spreadCharges)
                         {
-                            gridValuesChecker.checkReal(point.second, point.first.c_str());
+                            /* The wrapped grid */
+                            SparseRealGridValuesOutput nonZeroGridValues = pmeGetRealGrid(pmeSafe.get(), mode.first);
+                            TestReferenceChecker       gridValuesChecker(rootChecker.checkCompound("NonZeroGridValues", "RealSpaceGrid"));
+                            const auto                 ulpToleranceGrid = 2 * ulpToleranceSplineValues * (int)(ceil(sqrt(atomCount)));
+                            /* 2 is empiric; sqrt(atomCount) assumes all the input charges may spread onto the same cell */
+                            SCOPED_TRACE(formatString("Testing grid values with tolerance of %ld", ulpToleranceGrid));
+                            if (!gridValuesSizeAssigned)
+                            {
+                                previousGridValuesSize = nonZeroGridValues.size();
+                                gridValuesSizeAssigned = true;
+                            }
+                            else
+                            {
+                                EXPECT_EQ(previousGridValuesSize, nonZeroGridValues.size());
+                            }
+
+                            gridValuesChecker.setDefaultTolerance(relativeToleranceAsUlp(1.0, ulpToleranceGrid));
+                            for (const auto &point : nonZeroGridValues)
+                            {
+                                gridValuesChecker.checkReal(point.second, point.first.c_str());
+                            }
                         }
                     }
                 }
@@ -201,7 +239,7 @@ class PmeSplineAndSpreadTest : public ::testing::TestWithParam<SplineAndSpreadIn
 };
 
 
-/*! \brief Test for PME B-spline moduli computation */
+/*! \brief Test for spline parameter computation and charge spreading. */
 TEST_P(PmeSplineAndSpreadTest, ReproducesOutputs)
 {
     EXPECT_NO_THROW(runTest());
@@ -210,7 +248,7 @@ TEST_P(PmeSplineAndSpreadTest, ReproducesOutputs)
 /* Valid input instances */
 
 //! A couple of valid inputs for boxes.
-static std::vector<Matrix3x3> const sampleBoxes
+static std::vector<Matrix3x3> const c_sampleBoxes
 {
     // normal box
     Matrix3x3 {{
@@ -227,7 +265,7 @@ static std::vector<Matrix3x3> const sampleBoxes
 };
 
 //! A couple of valid inputs for grid sizes.
-static std::vector<IVec> const sampleGridSizes
+static std::vector<IVec> const c_sampleGridSizes
 {
     IVec {
         16, 12, 14
@@ -238,19 +276,19 @@ static std::vector<IVec> const sampleGridSizes
 };
 
 //! Random charges
-static std::vector<real> const sampleChargesFull
+static std::vector<real> const c_sampleChargesFull
 {
     4.95f, 3.11f, 3.97f, 1.08f, 2.09f, 1.1f, 4.13f, 3.31f, 2.8f, 5.83f, 5.09f, 6.1f, 2.86f, 0.24f, 5.76f, 5.19f, 0.72f
 };
 //! 1 charge
-static auto const sampleCharges1 = ChargesVector::fromVector(sampleChargesFull.begin(), sampleChargesFull.begin() + 1);
+static auto const c_sampleCharges1 = ChargesVector::fromVector(c_sampleChargesFull.begin(), c_sampleChargesFull.begin() + 1);
 //! 2 charges
-static auto const sampleCharges2 = ChargesVector::fromVector(sampleChargesFull.begin() + 1, sampleChargesFull.begin() + 3);
+static auto const c_sampleCharges2 = ChargesVector::fromVector(c_sampleChargesFull.begin() + 1, c_sampleChargesFull.begin() + 3);
 //! 13 charges
-static auto const sampleCharges13 = ChargesVector::fromVector(sampleChargesFull.begin() + 3, sampleChargesFull.begin() + 16);
+static auto const c_sampleCharges13 = ChargesVector::fromVector(c_sampleChargesFull.begin() + 3, c_sampleChargesFull.begin() + 16);
 
 //! Random coordinate vectors
-static CoordinatesVector const sampleCoordinatesFull
+static CoordinatesVector const c_sampleCoordinatesFull
 {
     {
         5.59f, 1.37f, 0.95f
@@ -289,33 +327,34 @@ static CoordinatesVector const sampleCoordinatesFull
     }
 };
 //! 1 coordinate vector
-static CoordinatesVector const sampleCoordinates1(sampleCoordinatesFull.begin(), sampleCoordinatesFull.begin() + 1);
+static CoordinatesVector const c_sampleCoordinates1(c_sampleCoordinatesFull.begin(), c_sampleCoordinatesFull.begin() + 1);
 //! 2 coordinate vectors
-static CoordinatesVector const sampleCoordinates2(sampleCoordinatesFull.begin() + 1, sampleCoordinatesFull.begin() + 3);
+static CoordinatesVector const c_sampleCoordinates2(c_sampleCoordinatesFull.begin() + 1, c_sampleCoordinatesFull.begin() + 3);
 //! 13 coordinate vectors
-static CoordinatesVector const sampleCoordinates13(sampleCoordinatesFull.begin() + 3, sampleCoordinatesFull.begin() + 16);
+static CoordinatesVector const c_sampleCoordinates13(c_sampleCoordinatesFull.begin() + 3, c_sampleCoordinatesFull.begin() + 16);
 
 //! moved out from instantiantions for readability
-auto inputBoxes     = ::testing::ValuesIn(sampleBoxes);
+auto c_inputBoxes     = ::testing::ValuesIn(c_sampleBoxes);
 //! moved out from instantiantions for readability
-auto inputPmeOrders = ::testing::Range(3, 5 + 1);
+auto c_inputPmeOrders = ::testing::Range(3, 5 + 1);
 //! moved out from instantiantions for readability
-auto inputGridSizes = ::testing::ValuesIn(sampleGridSizes);
+auto c_inputGridSizes = ::testing::ValuesIn(c_sampleGridSizes);
 
-/*! \brief Instantiation of the PME spline computation test with valid input and 1 atom */
-INSTANTIATE_TEST_CASE_P(SaneInput1, PmeSplineAndSpreadTest, ::testing::Combine(inputBoxes, inputPmeOrders, inputGridSizes,
-                                                                                   ::testing::Values(sampleCoordinates1),
-                                                                                   ::testing::Values(sampleCharges1)
+/*! \brief Instantiation of the test with valid input and 1 atom */
+INSTANTIATE_TEST_CASE_P(SaneInput1, PmeSplineAndSpreadTest, ::testing::Combine(c_inputBoxes, c_inputPmeOrders, c_inputGridSizes,
+                                                                                   ::testing::Values(c_sampleCoordinates1),
+                                                                                   ::testing::Values(c_sampleCharges1)
                                                                                ));
-/*! \brief Instantiation of the PME spline computation test with valid input and 2 atoms */
-INSTANTIATE_TEST_CASE_P(SaneInput2, PmeSplineAndSpreadTest, ::testing::Combine(inputBoxes, inputPmeOrders, inputGridSizes,
-                                                                                   ::testing::Values(sampleCoordinates2),
-                                                                                   ::testing::Values(sampleCharges2)
+
+/*! \brief Instantiation of the test with valid input and 2 atoms */
+INSTANTIATE_TEST_CASE_P(SaneInput2, PmeSplineAndSpreadTest, ::testing::Combine(c_inputBoxes, c_inputPmeOrders, c_inputGridSizes,
+                                                                                   ::testing::Values(c_sampleCoordinates2),
+                                                                                   ::testing::Values(c_sampleCharges2)
                                                                                ));
-/*! \brief Instantiation of the PME spline computation test with valid input and 13 atoms */
-INSTANTIATE_TEST_CASE_P(SaneInput13, PmeSplineAndSpreadTest, ::testing::Combine(inputBoxes, inputPmeOrders, inputGridSizes,
-                                                                                    ::testing::Values(sampleCoordinates13),
-                                                                                    ::testing::Values(sampleCharges13)
+/*! \brief Instantiation of the test with valid input and 13 atoms */
+INSTANTIATE_TEST_CASE_P(SaneInput13, PmeSplineAndSpreadTest, ::testing::Combine(c_inputBoxes, c_inputPmeOrders, c_inputGridSizes,
+                                                                                    ::testing::Values(c_sampleCoordinates13),
+                                                                                    ::testing::Values(c_sampleCharges13)
                                                                                 ));
 }
 }