add autotuner

CMakeLists_files.cmake

@@ -219,6 +219,7 @@ if (HAVE_CUDA)
   ADD_CUDA_OR_HIP_FILE(MAIN_SOURCE_FILES opm/simulators/linalg set_device.cpp)
 
   # HEADERS
+  ADD_CUDA_OR_HIP_FILE(PUBLIC_HEADER_FILES opm/simulators/linalg detail/autotuner.hpp)
   ADD_CUDA_OR_HIP_FILE(PUBLIC_HEADER_FILES opm/simulators/linalg detail/coloringAndReorderingUtils.hpp)
   ADD_CUDA_OR_HIP_FILE(PUBLIC_HEADER_FILES opm/simulators/linalg detail/cuda_safe_call.hpp)
   ADD_CUDA_OR_HIP_FILE(PUBLIC_HEADER_FILES opm/simulators/linalg detail/cusparse_matrix_operations.hpp)

opm/simulators/linalg/cuistl/CuDILU.cpp

@@ -25,6 +25,7 @@
 #include <opm/common/ErrorMacros.hpp>
 #include <opm/common/TimingMacros.hpp>
 #include <opm/simulators/linalg/GraphColoring.hpp>
+#include <opm/simulators/linalg/cuistl/detail/autotuner.hpp>
 #include <opm/simulators/linalg/cuistl/CuDILU.hpp>
 #include <opm/simulators/linalg/cuistl/CuSparseMatrix.hpp>
 #include <opm/simulators/linalg/cuistl/CuVector.hpp>
@@ -33,6 +34,8 @@
 #include <opm/simulators/linalg/cuistl/detail/preconditionerKernels/DILUKernels.hpp>
 #include <opm/simulators/linalg/matrixblock.hh>
 #include <tuple>
+#include <functional>
+#include <utility>
 namespace Opm::cuistl
 {
 
@@ -76,13 +79,9 @@ CuDILU<M, X, Y, l>::CuDILU(const M& A, bool splitMatrix, bool tuneKernels)
     }
     computeDiagAndMoveReorderedData();
 
-    // HIP does currently not support automtically picking thread block sizes as well as CUDA
-    // So only when tuning and using hip should we do our own manual tuning
-#ifdef USE_HIP
     if (m_tuneThreadBlockSizes) {
         tuneThreadBlockSizes();
     }
-#endif
 }
 
 template <class M, class X, class Y, int l>
@@ -111,7 +110,7 @@ CuDILU<M, X, Y, l>::apply(X& v, const Y& d)
                     m_gpuDInv.data(),
                     d.data(),
                     v.data(),
-                    m_applyThreadBlockSize);
+                    m_lowerSolveThreadBlockSize);
             } else {
                 detail::DILU::solveLowerLevelSet<field_type, blocksize_>(
                     m_gpuMatrixReordered->getNonZeroValues().data(),
@@ -123,7 +122,7 @@ CuDILU<M, X, Y, l>::apply(X& v, const Y& d)
                     m_gpuDInv.data(),
                     d.data(),
                     v.data(),
-                    m_applyThreadBlockSize);
+                    m_lowerSolveThreadBlockSize);
             }
             levelStartIdx += numOfRowsInLevel;
         }
@@ -143,7 +142,7 @@ CuDILU<M, X, Y, l>::apply(X& v, const Y& d)
                     numOfRowsInLevel,
                     m_gpuDInv.data(),
                     v.data(),
-                    m_applyThreadBlockSize);
+                    m_upperSolveThreadBlockSize);
             } else {
                 detail::DILU::solveUpperLevelSet<field_type, blocksize_>(
                     m_gpuMatrixReordered->getNonZeroValues().data(),
@@ -154,7 +153,7 @@ CuDILU<M, X, Y, l>::apply(X& v, const Y& d)
                     numOfRowsInLevel,
                     m_gpuDInv.data(),
                     v.data(),
-                    m_applyThreadBlockSize);
+                    m_upperSolveThreadBlockSize);
             }
         }
     }
@@ -202,7 +201,7 @@ CuDILU<M, X, Y, l>::computeDiagAndMoveReorderedData()
                 m_gpuMatrixReorderedDiag->data(),
                 m_gpuNaturalToReorder.data(),
                 m_gpuMatrixReorderedLower->N(),
-                m_updateThreadBlockSize);
+                m_moveThreadBlockSize);
         } else {
             detail::copyMatDataToReordered<field_type, blocksize_>(m_gpuMatrix.getNonZeroValues().data(),
                                                                    m_gpuMatrix.getRowIndices().data(),
@@ -210,7 +209,7 @@ CuDILU<M, X, Y, l>::computeDiagAndMoveReorderedData()
                                                                    m_gpuMatrixReordered->getRowIndices().data(),
                                                                    m_gpuNaturalToReorder.data(),
                                                                    m_gpuMatrixReordered->N(),
-                                                                   m_updateThreadBlockSize);
+                                                                   m_moveThreadBlockSize);
         }
 
         int levelStartIdx = 0;
@@ -230,7 +229,7 @@ CuDILU<M, X, Y, l>::computeDiagAndMoveReorderedData()
                     levelStartIdx,
                     numOfRowsInLevel,
                     m_gpuDInv.data(),
-                    m_updateThreadBlockSize);
+                    m_DILUFactorizationThreadBlockSize);
             } else {
                 detail::DILU::computeDiluDiagonal<field_type, blocksize_>(
                     m_gpuMatrixReordered->getNonZeroValues().data(),
@@ -241,7 +240,7 @@ CuDILU<M, X, Y, l>::computeDiagAndMoveReorderedData()
                     levelStartIdx,
                     numOfRowsInLevel,
                     m_gpuDInv.data(),
-                    m_updateThreadBlockSize);
+                    m_DILUFactorizationThreadBlockSize);
             }
             levelStartIdx += numOfRowsInLevel;
         }
@@ -252,54 +251,21 @@ template <class M, class X, class Y, int l>
 void
 CuDILU<M, X, Y, l>::tuneThreadBlockSizes()
 {
-    // TODO: generalize this code and put it somewhere outside of this class
-    long long bestApplyTime = std::numeric_limits<long long>::max();
-    long long bestUpdateTime = std::numeric_limits<long long>::max();
-    int bestApplyBlockSize = -1;
-    int bestUpdateBlockSize = -1;
-    int interval = 64;
+    using CuDILUType = std::remove_reference_t<decltype(*this)>;
 
-    // temporary buffers for the apply
+    // tune the thread-block size of the update function
+    auto updateFunc = std::bind(&CuDILUType::update, this);
+    detail::tuneThreadBlockSize(updateFunc, m_moveThreadBlockSize);
+    detail::tuneThreadBlockSize(updateFunc, m_DILUFactorizationThreadBlockSize);
+
+    // tune the thread-block size of the apply
+    auto applyFunc = std::bind(&CuDILUType::apply, this, std::placeholders::_1, std::placeholders::_1);
     CuVector<field_type> tmpV(m_gpuMatrix.N() * m_gpuMatrix.blockSize());
     CuVector<field_type> tmpD(m_gpuMatrix.N() * m_gpuMatrix.blockSize());
     tmpD = 1;
 
-    for (int thrBlockSize = interval; thrBlockSize <= 1024; thrBlockSize += interval) {
-        // sometimes the first kernel launch kan be slower, so take the time twice
-        for (int i = 0; i < 2; ++i) {
-
-            auto beforeUpdate = std::chrono::high_resolution_clock::now();
-            m_updateThreadBlockSize = thrBlockSize;
-            update();
-            std::ignore = cudaDeviceSynchronize();
-            auto afterUpdate = std::chrono::high_resolution_clock::now();
-            if (cudaSuccess == cudaGetLastError()) { // kernel launch was valid
-                long long durationInMicroSec
-                    = std::chrono::duration_cast<std::chrono::microseconds>(afterUpdate - beforeUpdate).count();
-                if (durationInMicroSec < bestUpdateTime) {
-                    bestUpdateTime = durationInMicroSec;
-                    bestUpdateBlockSize = thrBlockSize;
-                }
-            }
-
-            auto beforeApply = std::chrono::high_resolution_clock::now();
-            m_applyThreadBlockSize = thrBlockSize;
-            apply(tmpV, tmpD);
-            std::ignore = cudaDeviceSynchronize();
-            auto afterApply = std::chrono::high_resolution_clock::now();
-            if (cudaSuccess == cudaGetLastError()) { // kernel launch was valid
-                long long durationInMicroSec
-                    = std::chrono::duration_cast<std::chrono::microseconds>(afterApply - beforeApply).count();
-                if (durationInMicroSec < bestApplyTime) {
-                    bestApplyTime = durationInMicroSec;
-                    bestApplyBlockSize = thrBlockSize;
-                }
-            }
-        }
-    }
-
-    m_applyThreadBlockSize = bestApplyBlockSize;
-    m_updateThreadBlockSize = bestUpdateBlockSize;
+    detail::tuneThreadBlockSize(applyFunc, m_lowerSolveThreadBlockSize, tmpV, tmpD);
+    detail::tuneThreadBlockSize(applyFunc, m_upperSolveThreadBlockSize, tmpV, tmpD);
 }
 
 } // namespace Opm::cuistl

opm/simulators/linalg/cuistl/CuDILU.hpp

@@ -131,8 +131,10 @@ private:
     bool m_tuneThreadBlockSizes;
     //! \brief variables storing the threadblocksizes to use if using the tuned sizes and AMD cards
     //! The default value of -1 indicates that we have not calibrated and selected a value yet
-    int m_applyThreadBlockSize = -1;
-    int m_updateThreadBlockSize = -1;
+    int m_upperSolveThreadBlockSize = -1;
+    int m_lowerSolveThreadBlockSize = -1;
+    int m_moveThreadBlockSize = -1;
+    int m_DILUFactorizationThreadBlockSize = -1;
 };
 } // end namespace Opm::cuistl
 

opm/simulators/linalg/cuistl/OpmCuILU0.cpp

@@ -25,6 +25,7 @@
 #include <opm/common/ErrorMacros.hpp>
 #include <opm/common/TimingMacros.hpp>
 #include <opm/simulators/linalg/GraphColoring.hpp>
+#include <opm/simulators/linalg/cuistl/detail/autotuner.hpp>
 #include <opm/simulators/linalg/cuistl/CuSparseMatrix.hpp>
 #include <opm/simulators/linalg/cuistl/CuVector.hpp>
 #include <opm/simulators/linalg/cuistl/OpmCuILU0.hpp>
@@ -33,6 +34,8 @@
 #include <opm/simulators/linalg/cuistl/detail/preconditionerKernels/ILU0Kernels.hpp>
 #include <opm/simulators/linalg/matrixblock.hh>
 #include <tuple>
+#include <functional>
+#include <utility>
 namespace Opm::cuistl
 {
 
@@ -78,11 +81,9 @@ OpmCuILU0<M, X, Y, l>::OpmCuILU0(const M& A, bool splitMatrix, bool tuneKernels)
     }
     LUFactorizeAndMoveData();
 
-#ifdef USE_HIP
     if (m_tuneThreadBlockSizes) {
         tuneThreadBlockSizes();
     }
-#endif
 }
 
 template <class M, class X, class Y, int l>
@@ -110,7 +111,7 @@ OpmCuILU0<M, X, Y, l>::apply(X& v, const Y& d)
                     numOfRowsInLevel,
                     d.data(),
                     v.data(),
-                    m_applyThreadBlockSize);
+                    m_lowerSolveThreadBlockSize);
             } else {
                 detail::ILU0::solveLowerLevelSet<field_type, blocksize_>(m_gpuReorderedLU->getNonZeroValues().data(),
                                                                          m_gpuReorderedLU->getRowIndices().data(),
@@ -120,7 +121,7 @@ OpmCuILU0<M, X, Y, l>::apply(X& v, const Y& d)
                                                                          numOfRowsInLevel,
                                                                          d.data(),
                                                                          v.data(),
-                                                                         m_applyThreadBlockSize);
+                                                                         m_lowerSolveThreadBlockSize);
             }
             levelStartIdx += numOfRowsInLevel;
         }
@@ -139,7 +140,7 @@ OpmCuILU0<M, X, Y, l>::apply(X& v, const Y& d)
                     numOfRowsInLevel,
                     m_gpuMatrixReorderedDiag.value().data(),
                     v.data(),
-                    m_applyThreadBlockSize);
+                    m_upperSolveThreadBlockSize);
             } else {
                 detail::ILU0::solveUpperLevelSet<field_type, blocksize_>(m_gpuReorderedLU->getNonZeroValues().data(),
                                                                          m_gpuReorderedLU->getRowIndices().data(),
@@ -148,7 +149,7 @@ OpmCuILU0<M, X, Y, l>::apply(X& v, const Y& d)
                                                                          levelStartIdx,
                                                                          numOfRowsInLevel,
                                                                          v.data(),
-                                                                         m_applyThreadBlockSize);
+                                                                         m_upperSolveThreadBlockSize);
             }
         }
     }
@@ -196,7 +197,7 @@ OpmCuILU0<M, X, Y, l>::LUFactorizeAndMoveData()
                 m_gpuMatrixReorderedDiag.value().data(),
                 m_gpuNaturalToReorder.data(),
                 m_gpuMatrixReorderedLower->N(),
-                m_updateThreadBlockSize);
+                m_moveThreadBlockSize);
         } else {
             detail::copyMatDataToReordered<field_type, blocksize_>(m_gpuMatrix.getNonZeroValues().data(),
                                                                    m_gpuMatrix.getRowIndices().data(),
@@ -204,7 +205,7 @@ OpmCuILU0<M, X, Y, l>::LUFactorizeAndMoveData()
                                                                    m_gpuReorderedLU->getRowIndices().data(),
                                                                    m_gpuNaturalToReorder.data(),
                                                                    m_gpuReorderedLU->N(),
-                                                                   m_updateThreadBlockSize);
+                                                                   m_moveThreadBlockSize);
         }
         int levelStartIdx = 0;
         for (int level = 0; level < m_levelSets.size(); ++level) {
@@ -223,7 +224,7 @@ OpmCuILU0<M, X, Y, l>::LUFactorizeAndMoveData()
                     m_gpuNaturalToReorder.data(),
                     levelStartIdx,
                     numOfRowsInLevel,
-                    m_updateThreadBlockSize);
+                    m_ILU0FactorizationThreadBlockSize);
 
             } else {
                 detail::ILU0::LUFactorization<field_type, blocksize_>(m_gpuReorderedLU->getNonZeroValues().data(),
@@ -233,7 +234,7 @@ OpmCuILU0<M, X, Y, l>::LUFactorizeAndMoveData()
                                                                       m_gpuReorderToNatural.data(),
                                                                       numOfRowsInLevel,
                                                                       levelStartIdx,
-                                                                      m_updateThreadBlockSize);
+                                                                      m_ILU0FactorizationThreadBlockSize);
             }
             levelStartIdx += numOfRowsInLevel;
         }
@@ -244,54 +245,21 @@ template <class M, class X, class Y, int l>
 void
 OpmCuILU0<M, X, Y, l>::tuneThreadBlockSizes()
 {
-    // TODO generalize this tuning process in a function separate of the class
-    long long bestApplyTime = std::numeric_limits<long long>::max();
-    long long bestUpdateTime = std::numeric_limits<long long>::max();
-    int bestApplyBlockSize = -1;
-    int bestUpdateBlockSize = -1;
-    int interval = 64;
+    using CuDILUType = std::remove_reference_t<decltype(*this)>;
 
-    // temporary buffers for the apply
+    // tune the thread-block size of the update function
+    auto updateFunc = std::bind(&CuDILUType::update, this);
+    detail::tuneThreadBlockSize(updateFunc, m_moveThreadBlockSize);
+    detail::tuneThreadBlockSize(updateFunc, m_ILU0FactorizationThreadBlockSize);
+
+    // tune the thread-block size of the apply
+    auto applyFunc = std::bind(&CuDILUType::apply, this, std::placeholders::_1, std::placeholders::_1);
     CuVector<field_type> tmpV(m_gpuMatrix.N() * m_gpuMatrix.blockSize());
     CuVector<field_type> tmpD(m_gpuMatrix.N() * m_gpuMatrix.blockSize());
     tmpD = 1;
 
-    for (int thrBlockSize = interval; thrBlockSize <= 1024; thrBlockSize += interval) {
-        // sometimes the first kernel launch kan be slower, so take the time twice
-        for (int i = 0; i < 2; ++i) {
-
-            auto beforeUpdate = std::chrono::high_resolution_clock::now();
-            m_updateThreadBlockSize = thrBlockSize;
-            update();
-            std::ignore = cudaDeviceSynchronize();
-            auto afterUpdate = std::chrono::high_resolution_clock::now();
-            if (cudaSuccess == cudaGetLastError()) { // kernel launch was valid
-                long long durationInMicroSec
-                    = std::chrono::duration_cast<std::chrono::microseconds>(afterUpdate - beforeUpdate).count();
-                if (durationInMicroSec < bestUpdateTime) {
-                    bestUpdateTime = durationInMicroSec;
-                    bestUpdateBlockSize = thrBlockSize;
-                }
-            }
-
-            auto beforeApply = std::chrono::high_resolution_clock::now();
-            m_applyThreadBlockSize = thrBlockSize;
-            apply(tmpV, tmpD);
-            std::ignore = cudaDeviceSynchronize();
-            auto afterApply = std::chrono::high_resolution_clock::now();
-            if (cudaSuccess == cudaGetLastError()) { // kernel launch was valid
-                long long durationInMicroSec
-                    = std::chrono::duration_cast<std::chrono::microseconds>(afterApply - beforeApply).count();
-                if (durationInMicroSec < bestApplyTime) {
-                    bestApplyTime = durationInMicroSec;
-                    bestApplyBlockSize = thrBlockSize;
-                }
-            }
-        }
-    }
-
-    m_applyThreadBlockSize = bestApplyBlockSize;
-    m_updateThreadBlockSize = bestUpdateBlockSize;
+    detail::tuneThreadBlockSize(applyFunc, m_lowerSolveThreadBlockSize, tmpV, tmpD);
+    detail::tuneThreadBlockSize(applyFunc, m_upperSolveThreadBlockSize, tmpV, tmpD);
 }
 
 } // namespace Opm::cuistl

opm/simulators/linalg/cuistl/OpmCuILU0.hpp

@@ -131,8 +131,10 @@ private:
     bool m_tuneThreadBlockSizes;
     //! \brief variables storing the threadblocksizes to use if using the tuned sizes and AMD cards
     //! The default value of -1 indicates that we have not calibrated and selected a value yet
-    int m_applyThreadBlockSize = -1;
-    int m_updateThreadBlockSize = -1;
+    int m_upperSolveThreadBlockSize = -1;
+    int m_lowerSolveThreadBlockSize = -1;
+    int m_moveThreadBlockSize = -1;
+    int m_ILU0FactorizationThreadBlockSize = -1;
 };
 } // end namespace Opm::cuistl
 

opm/simulators/linalg/cuistl/detail/autotuner.hpp

@@ -0,0 +1,88 @@
+/*
+  Copyright 2024 SINTEF AS
+  This file is part of the Open Porous Media project (OPM).
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+  OPM 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 General Public License for more details.
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+#ifndef OPM_AUTOTUNER_HPP
+#define OPM_AUTOTUNER_HPP
+
+#include <cuda.h>
+#include <cuda_runtime.h>
+#include <opm/common/ErrorMacros.hpp>
+#include <opm/simulators/linalg/cuistl/detail/cuda_safe_call.hpp>
+#include <functional>
+#include <utility>
+#include <limits>
+
+namespace Opm::cuistl::detail
+{
+
+    /// @brief Function that tests the best thread block size, assumes updating the reference will affect runtimes
+    /// @tparam func function to tune
+    /// @tparam ...Args types of the arguments needed to call the function
+    /// @param f the function to tune
+    /// @param threadBlockSize reference to the thread block size that will affect kernel executions
+    /// @param ...args arguments needed by the function
+    template <typename func, typename... Args>
+    void tuneThreadBlockSize(func f, int& threadBlockSize, Args&&... args) {
+
+        // TODO: figure out a more rigorous way of deciding how many runs will suffice?
+        constexpr const int runs = 2;
+        cudaEvent_t events[runs+1];
+
+        // create the events
+        for (int i = 0; i < runs + 1; ++i){
+            OPM_CUDA_SAFE_CALL(cudaEventCreate(&events[i]));
+        }
+
+        // Initialize helper variables
+        float bestTime = std::numeric_limits<float>::max();
+        int bestBlockSize = -1;
+        int interval = 64;
+
+        // try each possible blocksize
+        for (int thrBlockSize = interval; thrBlockSize <= 1024; thrBlockSize += interval){
+            // update the blocksize
+            threadBlockSize = thrBlockSize;
+
+            // record a first event, and then an event after each kernel
+            OPM_CUDA_SAFE_CALL(cudaEventRecord(events[0]));
+            for (int i = 0; i < runs; ++i){
+                f(std::forward<Args>(args)...); // runs an arbitrary function with the provided arguments
+                OPM_CUDA_SAFE_CALL(cudaEventRecord(events[i+1]));
+            }
+
+            // make suret he runs are over
+            OPM_CUDA_SAFE_CALL(cudaEventSynchronize(events[runs]));
+
+            // kernel launch was valid
+            if (cudaSuccess == cudaGetLastError()){
+                // check if we beat the record for the fastest kernel
+                for (int i = 0; i < runs; ++i){
+                    float candidateBlockSizeTime;
+                    OPM_CUDA_SAFE_CALL(cudaEventElapsedTime(&candidateBlockSizeTime, events[i], events[i+1]));
+                    if (candidateBlockSizeTime < bestTime){
+                        bestTime = candidateBlockSizeTime;
+                        bestBlockSize = thrBlockSize;
+                    }
+                }
+            }
+        }
+        // TODO: have this only be printed if linear solver verbosity >0?
+        printf("best size: %d, best time %f\n", bestBlockSize, bestTime);
+
+        threadBlockSize = bestBlockSize;
+    }
+
+} // end namespace Opm::cuistl::detail
+
+#endif