diff --git a/CMakeLists.txt b/CMakeLists.txt index 24be243c..6664348d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -252,22 +252,33 @@ install(FILES cppuddle-config.cmake DESTINATION ${CMAKE_INSTALL_PREFIX}/lib/cmak install(EXPORT CPPuddle NAMESPACE CPPuddle:: DESTINATION ${CMAKE_INSTALL_PREFIX}/lib/cmake/CPPuddle/) #Boost::boost Boost::program_options HPX::hpx Kokkos::kokkos HPXKokkos::hpx_kokkos buffer_manager stream_manager +if (CPPUDDLE_WITH_CUDA) + add_hpx_executable( + cuda_gemm_fusion_experiments + DEPENDENCIES + Boost::boost Boost::program_options atomic Kokkos::kokkos HPXKokkos::hpx_kokkos HPX::hpx buffer_manager stream_manager + COMPONENT_DEPENDENCIES iostreams + SOURCES + benchmarks/cuda_gemm_fusion_experiments.cpp + ) +endif() + if (CPPUDDLE_WITH_TESTS) if (CPPUDDLE_WITH_CUDA) add_hpx_executable( recycling-with-hpx-cuda DEPENDENCIES - Boost::boost Boost::program_options HPX::hpx buffer_manager stream_manager + Boost::boost Boost::program_options atomic Kokkos::kokkos HPXKokkos::hpx_kokkos HPX::hpx buffer_manager stream_manager COMPONENT_DEPENDENCIES iostreams SOURCES - examples/recycling-with-hpx-cuda.cu + examples/recycling-with-hpx-cuda.cpp ) if (CPPUDDLE_WITH_KOKKOS) add_hpx_executable( recycling-with-hpx-kokkos DEPENDENCIES - Boost::boost Boost::program_options Kokkos::kokkos HPXKokkos::hpx_kokkos HPX::hpx buffer_manager stream_manager + Boost::boost Boost::program_options atomic Kokkos::kokkos HPXKokkos::hpx_kokkos HPX::hpx buffer_manager stream_manager COMPONENT_DEPENDENCIES iostreams SOURCES examples/recycling-with-hpx-kokkos.cpp @@ -275,7 +286,7 @@ if (CPPUDDLE_WITH_TESTS) add_hpx_executable( kernel-aggregation-with-hpx-kokkos DEPENDENCIES - Boost::boost Boost::program_options Kokkos::kokkos HPXKokkos::hpx_kokkos HPX::hpx buffer_manager stream_manager + Boost::boost Boost::program_options atomic Kokkos::kokkos HPXKokkos::hpx_kokkos HPX::hpx buffer_manager stream_manager COMPONENT_DEPENDENCIES iostreams SOURCES examples/kernel-aggregation-with-hpx-kokkos.cpp @@ -320,7 +331,7 @@ if (CPPUDDLE_WITH_TESTS) add_executable(allocator_cuda_test tests/allocator_cuda_test.cu) target_link_libraries(allocator_cuda_test - PRIVATE Kokkos::kokkos HPXKokkos::hpx_kokkos buffer_manager) + PRIVATE atomic Kokkos::kokkos HPXKokkos::hpx_kokkos buffer_manager) target_compile_definitions(allocator_cuda_test PRIVATE CPPUDDLE_HAVE_CUDA) add_executable( @@ -338,13 +349,13 @@ if (CPPUDDLE_WITH_TESTS) tests/allocator_kokkos_test.cpp include/buffer_manager.hpp) target_link_libraries(allocator_kokkos_test - PRIVATE Boost::boost Boost::program_options HPX::hpx Kokkos::kokkos HPXKokkos::hpx_kokkos buffer_manager) + PRIVATE Boost::boost Boost::program_options HPX::hpx atomic Kokkos::kokkos HPXKokkos::hpx_kokkos buffer_manager) target_compile_definitions(allocator_kokkos_test PRIVATE HPX_WITH_CUDA CPPUDDLE_HAVE_CUDA) add_hpx_executable( allocator_kokkos_executor_for_loop_test DEPENDENCIES - HPX::hpx Kokkos::kokkos HPXKokkos::hpx_kokkos buffer_manager + HPX::hpx atomic Kokkos::kokkos HPXKokkos::hpx_kokkos buffer_manager SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/tests/allocator_kokkos_executor_for_loop_test.cpp) # target_link_libraries(allocator_kokkos_executor_for_loop_test @@ -354,7 +365,7 @@ if (CPPUDDLE_WITH_TESTS) add_hpx_executable( work_aggregation_test DEPENDENCIES - Boost::boost Boost::program_options HPX::hpx Kokkos::kokkos HPXKokkos::hpx_kokkos buffer_manager stream_manager + Boost::boost Boost::program_options HPX::hpx atomic Kokkos::kokkos HPXKokkos::hpx_kokkos buffer_manager stream_manager COMPONENT_DEPENDENCIES iostreams SOURCES tests/work_aggregation_test.cpp @@ -363,7 +374,7 @@ if (CPPUDDLE_WITH_TESTS) add_hpx_executable( work_aggregation_cpu_triad DEPENDENCIES - Boost::boost Boost::program_options HPX::hpx Kokkos::kokkos HPXKokkos::hpx_kokkos buffer_manager stream_manager + Boost::boost Boost::program_options HPX::hpx atomic Kokkos::kokkos HPXKokkos::hpx_kokkos buffer_manager stream_manager COMPONENT_DEPENDENCIES iostreams SOURCES tests/work_aggregation_cpu_triad.cpp @@ -372,7 +383,7 @@ if (CPPUDDLE_WITH_TESTS) add_hpx_executable( work_aggregation_cuda_triad DEPENDENCIES - Boost::boost Boost::program_options HPX::hpx Kokkos::kokkos HPXKokkos::hpx_kokkos buffer_manager stream_manager + Boost::boost Boost::program_options HPX::hpx atomic Kokkos::kokkos HPXKokkos::hpx_kokkos buffer_manager stream_manager COMPONENT_DEPENDENCIES iostreams SOURCES tests/work_aggregation_cuda_triad.cpp diff --git a/README.md b/README.md index 288f6332..e6aa1c66 100644 --- a/README.md +++ b/README.md @@ -48,4 +48,8 @@ find_package(CPPuddle REQUIRED) cmake -H/path/to/source -B$/path/to/build -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=/path/to/install/cppuddle -DCPPUDDLE_WITH_HPX=ON -DCPPUDDLE_WITH_HPX_AWARE_ALLOCATORS=ON -DCPPUDDLE_WITH_TESTS=OFF -DCPPUDDLE_WITH_COUNTERS=OFF ``` +#### Disclosure + +CPPuddle contains contributions developed at Los Alamos National Laboratory (LANL). +LANL software approval number: O5227 diff --git a/benchmarks/cuda_gemm_fusion_experiments.cpp b/benchmarks/cuda_gemm_fusion_experiments.cpp new file mode 100644 index 00000000..9eccb974 --- /dev/null +++ b/benchmarks/cuda_gemm_fusion_experiments.cpp @@ -0,0 +1,909 @@ +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "gemm_cpu_utils.hpp" +#include "cuda_gemm_kernels.hpp" + +// ------------------------------------------------------------ +// Computes: +// C = alpha * A * B + beta * C +// Matrix layouts: +// A: M x K, row-major +// B: K x N, row-major +// C: M x N, row-major +// ------------------------------------------------------------ + +// serial / parallel extra argument: Number tasks +// Top-down, updown extra argument: Depth and structure (number childs per node -- lets fix this to 2 for now) +enum class graph_mode {scalar = 0, serial, parallel, top_down, bottom_up}; +enum class fusion_mode {none = 0, cuda_graph, cppuddle}; +enum class futurization_mode {none = 0, async, then, dataflow, sender_receiver}; + +struct experiment_setup { + size_t utilized_implementation_id; + size_t num_repetitions; + graph_mode utilized_graph_mode; + size_t extra_graph_parameter; + bool postprocessing_active; + fusion_mode utilized_fusion_approach; + futurization_mode utilized_futurization_mode; + futurization_mode futurized; +}; + + +struct empty_t {}; + +template +inline constexpr bool is_void_or_future_void_v = + std::is_void_v || std::is_same_v>; + + +template >> +R kernel_fusion_benchmark_iteration(const size_t utilized_implementation_id, + const double alpha, const double beta, const std::vector &hA, const std::vector &hB, + std::vector &hC, double *dA, double *dB, double *dC, const std::vector &hC_comparison, + const size_t M, const size_t N, const size_t K, cublasHandle_t &handle, const graph_mode utilized_graph_mode, + const size_t extra_graph_parameter, const bool postprocessing_active) { + + // auto policy = hpx::launch::async_policy(hpx::threads::thread_priority::default_, + // hpx::threads::thread_stacksize::large); + + using result_t = std::conditional_t, empty_t, R>; + result_t ret_fut; + if constexpr (!std::is_void_v) + ret_fut = hpx::make_ready_future(); + + // if graph == parallel + // Launch one other kernel recursively for + // return when-all future + // if (utilized_graph_mode == graph_mode::parallel) { + // if (extra_graph_parameter > 1) { + // ret_fut = hpx::async(policy, &kernel_fusion_benchmark_iteration, utilized_implementation_id, alpha, beta, hA, hB, std::ref(hC), + // dA, dB, dC, hC_comparison, M, N, K, std::ref(handle), utilized_graph_mode, extra_graph_parameter - 1, postprocessing_active); + // } + // } + + // if graph == bottom-up + // Launch two more kernels -- wait until they're done (get), and only then launch this one) + // do not return future + // if (utilized_graph_mode == graph_mode::bottom_up) { + // if (extra_graph_parameter > 1) { + // auto fut1 = hpx::async(policy, &kernel_fusion_benchmark_iteration, utilized_implementation_id, alpha, beta, hA, hB, std::ref(hC), + // dA, dB, dC, hC_comparison, M, N, K, std::ref(handle), utilized_graph_mode, extra_graph_parameter - 1, postprocessing_active); + // auto fut2 = hpx::async(policy, &kernel_fusion_benchmark_iteration, utilized_implementation_id, alpha, beta, hA, hB, std::ref(hC), + // dA, dB, dC, hC_comparison, M, N, K, std::ref(handle), utilized_graph_mode, extra_graph_parameter - 1, postprocessing_active); + // ret_fut = hpx::when_all(fut1, fut2); + // } + // } + + if (utilized_implementation_id == 0) { + dim3 block(16, 16); + dim3 grid( + (N + block.x - 1) / block.x, + (M + block.y - 1) / block.y + ); + if (postprocessing_active) + reset_C_gpu(M, N, dC); + // start_timer = std::chrono::high_resolution_clock::now(); + // preprocessing? + dgemm_naive_unoptimized_v0<<>>( + M, N, K, alpha, dA, dB, beta, dC); + // postprocessing + if (postprocessing_active && !check_gpu_results(M, N, dC, hC_comparison.data())) + std::cout << "-> ERROR: Variant 0 faulty!" << std::endl; + } else if (utilized_implementation_id == 1) { + dim3 block(16, 16); + dim3 grid( + (N + block.x - 1) / block.x, + (M + block.y - 1) / block.y + ); + if (postprocessing_active) + reset_C_gpu(M, N, dC); + dgemm_naive_v1<<>>( + M, N, K, alpha, dA, dB, beta, dC); + if (postprocessing_active && !check_gpu_results(M, N, dC, hC_comparison.data())) + std::cout << "-> ERROR: Variant 1 faulty!" << std::endl; + } else if (utilized_implementation_id == 2) { + constexpr int TILE = 32; + dim3 block(TILE, TILE); + dim3 grid( + (N + TILE - 1) / TILE, + (M + TILE - 1) / TILE + ); + if (postprocessing_active) + reset_C_gpu(M, N, dC); + dgemm_shared_tiled_v2<<>>( + M, N, K, alpha, dA, dB, beta, dC); + if (postprocessing_active && !check_gpu_results(M, N, dC, hC_comparison.data())) + std::cout << "-> ERROR: Variant 2 faulty!" << std::endl; + } else if (utilized_implementation_id == 3) { + constexpr int BM = 64, BN = 64, BK = 8, TM = 4, TN = 4; + static_assert(BM % TM == 0, "BM must be divisible by TM"); + static_assert(BN % TN == 0, "BN must be divisible by TN"); + dim3 block(BN / TN, BM / TM); + dim3 grid((N + BN - 1) / BN, (M + BM - 1) / BM); + if (postprocessing_active) + reset_C_gpu(M, N, dC); + dgemm_register_tiled_v3<<>>( + M, N, K, alpha, dA, dB, beta, dC); + if (postprocessing_active && !check_gpu_results(M, N, dC, hC_comparison.data())) + std::cout << "-> ERROR: Variant 3 faulty!" << std::endl; + } else if (utilized_implementation_id == 4) { + constexpr int BM = 64, BN = 32, BK = 16, TM = 4, TN = 2; + dim3 block(BN / TN, BM / TM); + dim3 grid( + (N + BN - 1) / BN, + (M + BM - 1) / BM + ); + if (postprocessing_active) + reset_C_gpu(M, N, dC); + dgemm_v4_shared_register_tiled + <<>>( + M, N, K, alpha, dA, dB, beta, dC); + if (postprocessing_active && !check_gpu_results(M, N, dC, hC_comparison.data())) + std::cout << "-> ERROR: Variant 4 faulty!" << std::endl; + } else if (utilized_implementation_id == 5) { + constexpr int BM = 64, BN = 64, BK = 8, TM = 4, TN = 4; + dim3 block(BN / TN, BM / TM); + dim3 grid((N + BN - 1) / BN, (M + BM - 1) / BM); + if (postprocessing_active) + reset_C_gpu(M, N, dC); + dgemm_v5_vectorized_register_tiled + <<>>( + M, N, K, alpha, dA, dB, beta, dC); + if (postprocessing_active && !check_gpu_results(M, N, dC, hC_comparison.data())) + std::cout << "-> ERROR: Variant 5 faulty!" << std::endl; + } else if (utilized_implementation_id == 6) { + constexpr int BM = 64, BN = 64, BK = 16; + if ((M % BM) != 0 || (N % BN) != 0 || (K % BK) != 0) { + std::cerr << "dgemm_v5b_fp64_wmma_tiled requires " + << "M %% 64 == 0, N %% 64 == 0, K %% 16 == 0\n" + << std::endl; + exit(0); + } + dim3 block(512); + dim3 grid(N / BN, M / BM); + if (postprocessing_active) + reset_C_gpu(M, N, dC); + dgemm_v6_fp64_wmma_tiled<<>>( + M, N, K, alpha, dA, dB, beta, dC); + if (postprocessing_active && !check_gpu_results(M, N, dC, hC_comparison.data())) + std::cout << "-> ERROR: Variant 6 faulty!" << std::endl; + } else if (utilized_implementation_id == 7) { + constexpr int BM = 64, BN = 64, BK = 16; + if ((M % BM) != 0 || (N % BN) != 0 || (K % BK) != 0) { + std::cerr << "dgemm_v7_fp64_wmma_64x64_8warps requires " + << "M %% 64 == 0, N %% 64 == 0, K %% 16 == 0\n" + << std::endl; + exit(0); + } + dim3 block(256); + dim3 grid(N / BN, M / BM); + if (postprocessing_active) + reset_C_gpu(M, N, dC); + dgemm_v7_fp64_wmma_64x64_8warps + <<>>(M, N, K, alpha, dA, dB, beta, dC); + if (postprocessing_active && !check_gpu_results(M, N, dC, hC_comparison.data())) + std::cout << "-> ERROR: Variant 7 faulty!" << std::endl; + } else if (utilized_implementation_id == 8) { // CUBLAS version + if (postprocessing_active) + reset_C_gpu(M, N, dC); + cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, &alpha, dB, N, dA, K, &beta, dC, N); + if (postprocessing_active && !check_gpu_results(M, N, dC, hC_comparison.data())) + std::cout << "-> ERROR: Variant CUBLAS faulty!" << std::endl; + } else if (utilized_implementation_id == 9) { // CPU version + reset_C_cpu(M, N, hC.data()); + dgemm_naive_cpu( + M, N, K, alpha, hA.data(), hB.data(), beta, hC.data()); + if (postprocessing_active && !compare_matrices(M, N, hC.data(), hC_comparison.data())) + std::cout << "-> ERROR: Variant CPU faulty!" << std::endl; + } + + if (utilized_graph_mode == graph_mode::serial) { + if (extra_graph_parameter > 1) { + + if constexpr (!std::is_void_v) { + ret_fut = ret_fut.then( + [=, &handle, &hA, &hB, &hC, &hC_comparison](auto && fut) { + kernel_fusion_benchmark_iteration(utilized_implementation_id, alpha, beta, hA, hB, std::ref(hC), + dA, dB, dC, hC_comparison, M, N, K, std::ref(handle), utilized_graph_mode, extra_graph_parameter - 1, postprocessing_active); + }); + } else { + kernel_fusion_benchmark_iteration(utilized_implementation_id, alpha, beta, hA, hB, std::ref(hC), + dA, dB, dC, hC_comparison, M, N, K, std::ref(handle), utilized_graph_mode, extra_graph_parameter - 1, postprocessing_active); + } + } + } + + // if (utilized_graph_mode == graph_mode::top_down) { + // if (extra_graph_parameter > 1) { + // auto fut1 = hpx::async(policy, &kernel_fusion_benchmark_iteration, utilized_implementation_id, alpha, beta, hA, hB, std::ref(hC), + // dA, dB, dC, hC_comparison, M, N, K, std::ref(handle), utilized_graph_mode, extra_graph_parameter - 1, postprocessing_active); + // auto fut2 = hpx::async(policy, &kernel_fusion_benchmark_iteration, utilized_implementation_id, alpha, beta, hA, hB, std::ref(hC), + // dA, dB, dC, hC_comparison, M, N, K, std::ref(handle), utilized_graph_mode, extra_graph_parameter - 1, postprocessing_active); + // ret_fut = hpx::when_all(fut1, fut2); + // } + // } + if constexpr (!std::is_void_v) + return ret_fut; +} + +void kernel_fusion_benchmark(const experiment_setup &experiment_parameters, + const std::vector &hA, const std::vector &hB, + std::vector &hC, const std::vector &hC_comparison, + const size_t M, const size_t N, const size_t K) { + + const size_t &utilized_implementation_id = experiment_parameters.utilized_implementation_id; + const size_t &num_repetitions = experiment_parameters.num_repetitions; + const graph_mode &utilized_graph_mode = experiment_parameters.utilized_graph_mode; + const size_t &extra_graph_parameter = experiment_parameters.extra_graph_parameter; + const bool &postprocessing_active = experiment_parameters.postprocessing_active; + const fusion_mode &utilized_fusion_approach = experiment_parameters.utilized_fusion_approach; + const futurization_mode &futurized = experiment_parameters.futurized; + + + const double flops = 2.0 * static_cast(M) + * static_cast(N) + * static_cast(K); + const double alpha = 1.0; + const double beta = 0.0; + const size_t sizeA = M * K; + const size_t sizeB = K * N; + const size_t sizeC = M * N; + const size_t bytesA = sizeA * sizeof(double); + const size_t bytesB = sizeB * sizeof(double); + const size_t bytesC = sizeC * sizeof(double); + double* dA = nullptr; + double* dB = nullptr; + double* dC = nullptr; + cudaMalloc(&dA, bytesA); + cudaMalloc(&dB, bytesB); + cudaMalloc(&dC, bytesC); + cudaMemcpy(dA, hA.data(), bytesA, cudaMemcpyHostToDevice); + cudaMemcpy(dB, hB.data(), bytesB, cudaMemcpyHostToDevice); + cudaMemcpy(dC, hC.data(), bytesC, cudaMemcpyHostToDevice); + cublasHandle_t handle; + cublasCreate(&handle); + cublasSetMathMode(handle, CUBLAS_TENSOR_OP_MATH); + cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, &alpha, dB, N, dA, K, &beta, dC, N); // first pass required -- I think there is some setup happening + + std::chrono::time_point start_timer = std::chrono::high_resolution_clock::now(); + // setup cuda graph + // setup executor pool + // setup multiple cuda streams + // setup buffers and copu initial data + // repeat experiment n times (otherwise the cuda graphs will lose badly) + // ... needs separate methods (cannot extend the recursion over the loop + + if (futurized == futurization_mode::then) { + auto loop_fut = hpx::make_ready_future(); + for (size_t it = 0; it < num_repetitions; it++) { + loop_fut = loop_fut.then( + [&](auto && fut) { + // fut.get(); + return kernel_fusion_benchmark_iteration>(utilized_implementation_id, alpha, beta, hA, hB, std::ref(hC), + dA, dB, dC, hC_comparison, M, N, K, std::ref(handle), utilized_graph_mode, extra_graph_parameter, postprocessing_active); + }); + // loop_fut = hpx::dataflow([&](auto void_val) { + // kernel_fusion_benchmark_iteration(utilized_implementation_id, alpha, beta, hA, hB, std::ref(hC), + // dA, dB, dC, hC_comparison, M, N, K, utilized_graph_mode, extra_graph_parameter, postprocessing_active); + // }, loop_fut); + // + // stack size is a problem with the hpx task + // can be fixed with: + // --hpx:ini=hpx.stacks.small_size=0x100000. + } + loop_fut.get(); + } else { + for (size_t it = 0; it < num_repetitions; it++) { + kernel_fusion_benchmark_iteration(utilized_implementation_id, alpha, beta, hA, hB, std::ref(hC), + dA, dB, dC, hC_comparison, M, N, K, std::ref(handle), utilized_graph_mode, extra_graph_parameter, postprocessing_active); + } + } + cudaGetLastError(); + cudaDeviceSynchronize(); + // cudaGetLastError(); + // cudaDeviceSynchronize(); + + auto end_timer = std::chrono::high_resolution_clock::now(); + std::chrono::duration elapsed = end_timer - start_timer; + double total_seconds = elapsed.count(); + double avg_seconds = total_seconds / static_cast(num_repetitions); + if (utilized_graph_mode == graph_mode::serial) + avg_seconds = avg_seconds / static_cast(extra_graph_parameter); + double gflops = flops / avg_seconds / 1.0e9; + std::cout << "-> Average Runtime: " << avg_seconds * 1.0e3 << "ms\n"; + std::cout << "-> Performance: " << gflops << " GFLOP/s"; + std::cout << std::endl; + +} + +int hpx_main(int argc, char *argv[]) { + std::cout << "Running benchmark with " << hpx::get_num_worker_threads() << " HPX worker threads!" << std::endl; + size_t M = 1024; + size_t N = 1024; + size_t K = 1024; + + if (argc >= 4) { + M = std::stoull(argv[1]); + N = std::stoull(argv[2]); + K = std::stoull(argv[3]); + } + + std::cout << "Using Matrix Dimensions: "; + std::cout << "M = " << M << ", N = " << N << ", K = " << K << "\n"; + + const double alpha = 1.0; + const double beta = 0.0; + + size_t sizeA = static_cast(M) * K; + size_t sizeB = static_cast(K) * N; + size_t sizeC = static_cast(M) * N; + + size_t bytesA = sizeA * sizeof(double); + size_t bytesB = sizeB * sizeof(double); + size_t bytesC = sizeC * sizeof(double); + + std::vector hA(sizeA); + std::vector hB(sizeB); + std::vector hC(sizeC); + std::vector hC_comparison(sizeC); + + const unsigned long long seed = 42; + std::mt19937_64 rng(seed); + fill_random(hA, rng); + fill_random(hB, rng); + fill_random(hC, rng); + reset_C_cpu(M, N, hC_comparison.data()); + + std::string alpha_string = std::to_string(alpha); + std::replace(alpha_string.begin(), alpha_string.end(), '.', '-'); + std::string beta_string = std::to_string(beta); + std::replace(beta_string.begin(), beta_string.end(), '.', '-'); + std::string reference_data_filename = std::string("reference_") + std::to_string(seed) + + std::string("_") + std::to_string(M) + std::string("_") + std::to_string(N) + std::string("_") + std::to_string(K) + + std::string("_") + alpha_string + std::string("_") + beta_string + + std::string(".txt"); + std::cout << "Checking if reference file for this the given M, N, K, alpha, beta and seed " << seed << " exists..." << std::endl; + std::cout << "Filename format: reference______.txt (alpha and beta use - instead of the decimal point)" << std::endl; + bool generate_results = true; + if (std::filesystem::exists(reference_data_filename)) { + std::cout << "Reference file " << reference_data_filename << " found!" << std::endl; + std::cout << "Loading reference file ... " << std::endl; + try { + std::ifstream in(reference_data_filename); + // C: M x N, row-major + for (auto row = 0; row < M; row++) { + for (auto column = 0; column < N; column++) { + in >> hC_comparison[row * N + column]; + } + } + in.close(); + generate_results = false; + std::cout << "Reference data loaded! Continuing with benchmarks..." << std::endl; + } catch (...) { + std::cerr << "Unknown error when reading reference data" << reference_data_filename << " !" << std::endl; + std::cout << "Continuing with generating input data after failed reading!" << std::endl; + } + + } else { + std::cout << "Reference file " << reference_data_filename << " not found in current directory!" << std::endl; + } + + if (generate_results) { + std::cout << "Generating reference data on CPU... " << std::endl; + dgemm_naive_cpu( + M, N, K, alpha, hA.data(), hB.data(), beta, hC_comparison.data()); + + std::cout << "Reference data generated! Storing in " << reference_data_filename << "..." << std::endl; + try { + std::ofstream out(reference_data_filename); + out << std::fixed << std::setprecision(std::numeric_limits::digits10 + 1); + // C: M x N, row-major + for (auto row = 0; row < M; row++) { + for (auto column = 0; column < N; column++) { + out << hC_comparison[row * N + column] << " "; + } + out << std::endl; + } + out.close(); + std::cout << "Reference data stored! Continuing with benchmarks..." << std::endl; + } catch (...) { + std::cerr << "Unknown error when writing reference data" << reference_data_filename << " !" << std::endl; + } + + } + + double* dA = nullptr; + double* dB = nullptr; + double* dC = nullptr; + + cudaMalloc(&dA, bytesA); + cudaMalloc(&dB, bytesB); + cudaMalloc(&dC, bytesC); + + cudaMemcpy(dA, hA.data(), bytesA, cudaMemcpyHostToDevice); + cudaMemcpy(dB, hB.data(), bytesB, cudaMemcpyHostToDevice); + cudaMemcpy(dC, hC.data(), bytesC, cudaMemcpyHostToDevice); + + dim3 block_v1(16, 16); + dim3 grid_v1( + (N + block_v1.x - 1) / block_v1.x, + (M + block_v1.y - 1) / block_v1.y + ); + + // ------------------------------------------------------------ + // FLOP estimate + // For GEMM: + // C = A * B + // FLOPs = 2 * M * N * K + // ------------------------------------------------------------ + double flops = 2.0 * static_cast(M) + * static_cast(N) + * static_cast(K); + std::cout << "Estimated FLOPs: " << flops << "\n"; + + // ------------------------------------------------------------ + // Timed execution + // Results + // ------------------------------------------------------------ + std::cout << "\nResults\n"; + std::cout << "-------\n\n"; + const size_t num_repetitions = 10; + const graph_mode utilized_graph_mode = graph_mode::serial; + const size_t extra_graph_parameter = 10; + const bool postprocessing_active = false; + + experiment_setup experiment_parameters; + experiment_parameters.num_repetitions = 10; + experiment_parameters.utilized_implementation_id = 9; + experiment_parameters.utilized_fusion_approach = fusion_mode::none; + experiment_parameters.utilized_graph_mode = graph_mode::serial; + experiment_parameters.extra_graph_parameter = 10; + experiment_parameters.postprocessing_active = false; + experiment_parameters.utilized_futurization_mode = futurization_mode::then; + + + // ============================================================ + // CPU Comparison variant (datA) + // ============================================================ + std::cout << "Variant Benchmark CPU: \n"; + experiment_parameters.utilized_implementation_id = 9; + kernel_fusion_benchmark(experiment_parameters, hA, hB, hC, hC_comparison, N, M, K); + std::cout << "Variant Benchmark CPU end \n" << std::endl; + // std::cout << "Variant CPU: Naive \n"; + + // auto start_cpu = std::chrono::high_resolution_clock::now(); + // dgemm_naive_cpu( + // M, N, K, alpha, hA.data(), hB.data(), beta, hC_comparison.data()); + // auto end_cpu = std::chrono::high_resolution_clock::now(); + // std::chrono::duration elapsed_cpu = end_cpu - start_cpu; + // double total_seconds_cpu = elapsed_cpu.count(); + // double avg_seconds_cpu = total_seconds_cpu / static_cast(1); + // double gflops_cpu = flops / avg_seconds_cpu / 1.0e9; + // std::cout << "-> Runtime: " << avg_seconds_cpu * 1.0e3 << "ms\n"; + // std::cout << "-> Variant CPU performance: " << gflops_cpu << " GFLOP/s\n"; + // std::cout << std::endl; + + // ============================================================ + // Variant 0 launch configuration + // ============================================================ + // std::cout << "Variant Benchmark 0: Naive \n"; + // kernel_fusion_benchmark(0, hA, hB, hC, hC_comparison, N, M, K); + // std::cout << "Variant Benchmark 0 end \n"; + std::cout << "Variant 0 Naive: \n"; + experiment_parameters.utilized_implementation_id = 0; + kernel_fusion_benchmark(experiment_parameters, hA, hB, hC, hC_comparison, N, M, K); + std::cout << "Variant 0 Naive \n" << std::endl; + + // std::cout << "Variant 0: Naive \n"; + // auto start_v0 = std::chrono::high_resolution_clock::now(); + // for (int repeat = 0; repeat < num_repeats; ++repeat) { + // dgemm_naive_unoptimized_v0<<>>( + // M, N, K, alpha, dA, dB, beta, dC); + // } + // cudaGetLastError(); + // cudaDeviceSynchronize(); + // auto end_v0 = std::chrono::high_resolution_clock::now(); + // std::chrono::duration elapsed_v0 = end_v0 - start_v0; + // double total_seconds_v0 = elapsed_v0.count(); + // double avg_seconds_v0 = total_seconds_v0 / static_cast(num_repeats); + // double gflops_v0 = flops / avg_seconds_v0 / 1.0e9; + // if (check_gpu_results(M, N, dC, hC_comparison.data())) + // std::cout << "-> Results correct!" << std::endl; + // else + // std::cout << "-> ERROR: Variant 0 faulty!" << std::endl; + // std::cout << "-> Average Runtime: " << avg_seconds_v0 * 1.0e3 << "ms\n"; + // std::cout << "-> Performance: " << gflops_v0 << " GFLOP/s\n"; + // std::cout << std::endl; + + + // ============================================================ + // Variant 1 launch configuration + // ============================================================ + std::cout << "Variant Benchmark 1: \n"; + std::cout << "-> Description: Naive but with optimized access pattern\n"; + experiment_parameters.utilized_implementation_id = 1; + kernel_fusion_benchmark(experiment_parameters, hA, hB, hC, hC_comparison, N, M, K); + std::cout << "Variant Benchmark 1 \n" << std::endl; + + // std::cout << "Variant 1: Naive but with optimized access pattern\n"; + // auto start_v1 = std::chrono::high_resolution_clock::now(); + // for (int repeat = 0; repeat < num_repeats; ++repeat) { + // dgemm_naive_v1<<>>( + // M, N, K, alpha, dA, dB, beta, dC); + // } + // cudaGetLastError(); + // cudaDeviceSynchronize(); + // auto end_v1 = std::chrono::high_resolution_clock::now(); + + // std::chrono::duration elapsed_v1 = end_v1 - start_v1; + // double total_seconds_v1 = elapsed_v1.count(); + // double avg_seconds_v1 = total_seconds_v1 / static_cast(num_repeats); + // double gflops_v1 = flops / avg_seconds_v1 / 1.0e9; + // std::cout << "-> Block: " << block_v1.x << " x " << block_v1.y << "\n"; + // std::cout << "-> Grid: " << grid_v1.x << " x " << grid_v1.y << "\n"; + // if (check_gpu_results(M, N, dC, hC_comparison.data())) + // std::cout << "-> Results correct!" << std::endl; + // else + // std::cout << "-> ERROR: Variant 1 faulty!" << std::endl; + // std::cout << "-> average runtime: " << avg_seconds_v1 * 1.0e3 << "ms\n"; + // std::cout << "-> performance: " << gflops_v1 << " GLOP/s\n"; + // std::cout << std::endl; + + // // Optionally copy result back to host + // cudaMemcpy(hC.data(), dC, bytesC, cudaMemcpyDeviceToHost); + + + // ============================================================ + // Variant 2 launch configuration + // ============================================================ + std::cout << "Variant Benchmark 2: \n"; + std::cout << "-> Description: Shared-memory tiled\n"; + experiment_parameters.utilized_implementation_id = 2; + kernel_fusion_benchmark(experiment_parameters, hA, hB, hC, hC_comparison, N, M, K); + std::cout << "Variant Benchmark 2 \n" << std::endl; + + // std::cout << "Variant 2: shared-memory tiled\n"; + + // constexpr int TILE_V2 = 32; + // dim3 block_v2(TILE_V2, TILE_V2); + // dim3 grid_v2( + // (N + TILE_V2 - 1) / TILE_V2, + // (M + TILE_V2 - 1) / TILE_V2 + // ); + // auto start_v2 = std::chrono::high_resolution_clock::now(); + // for (int repeat = 0; repeat < num_repeats; ++repeat) { + // dgemm_shared_tiled_v2<<>>( + // M, N, K, alpha, dA, dB, beta, dC); + // } + // cudaGetLastError(); + // cudaDeviceSynchronize(); + // auto end_v2 = std::chrono::high_resolution_clock::now(); + // std::chrono::duration elapsed_v2 = end_v2 - start_v2; + // double total_seconds_v2 = elapsed_v2.count(); + // double avg_seconds_v2 = total_seconds_v2 / static_cast(num_repeats); + // double gflops_v2 = flops / avg_seconds_v2 / 1.0e9; + // std::cout << "-> Block: " << block_v2.x << " x " << block_v2.y << "\n"; + // std::cout << "-> Grid: " << grid_v2.x << " x " << grid_v2.y << "\n"; + // if (check_gpu_results(M, N, dC, hC_comparison.data())) + // std::cout << "-> Results correct!" << std::endl; + // else + // std::cout << "-> ERROR: Variant 2 faulty!" << std::endl; + // std::cout << "-> average runtime: " << avg_seconds_v2 * 1.0e3 << "ms\n"; + // std::cout << "-> performance: " << gflops_v2 << " GLOP/s\n"; + // std::cout << std::endl; + + + // ============================================================ + // Variant 3 launch configuration + // ============================================================ + std::cout << "Variant Benchmark 3: \n"; + std::cout << "-> Description: Register tiled\n"; + experiment_parameters.utilized_implementation_id = 3; + kernel_fusion_benchmark(experiment_parameters, hA, hB, hC, hC_comparison, N, M, K); + std::cout << "Variant Benchmark end 3 \n" << std::endl; + + // std::cout << "Variant 3: register tiled\n"; + // constexpr int BM_v3 = 64; + // constexpr int BN_v3 = 64; + // constexpr int BK_v3 = 8; + // constexpr int TM_v3 = 4; + // constexpr int TN_v3 = 4; + + // static_assert(BM_v3 % TM_v3 == 0, "BM must be divisible by TM"); + // static_assert(BN_v3 % TN_v3 == 0, "BN must be divisible by TN"); + + // dim3 block_v3(BN_v3 / TN_v3, BM_v3 / TM_v3); + // dim3 grid_v3( + // (N + BN_v3 - 1) / BN_v3, + // (M + BM_v3 - 1) / BM_v3 + // ); + // auto start_v3 = std::chrono::high_resolution_clock::now(); + // for (int repeat = 0; repeat < num_repeats; ++repeat) { + // dgemm_register_tiled_v3<<>>( + // M, N, K, alpha, dA, dB, beta, dC); + // } + // cudaGetLastError(); + // cudaDeviceSynchronize(); + // auto end_v3 = std::chrono::high_resolution_clock::now(); + // std::chrono::duration elapsed_v3 = end_v3 - start_v3; + // double total_seconds_v3 = elapsed_v3.count(); + // double avg_seconds_v3 = total_seconds_v3 / static_cast(num_repeats); + // double gflops_v3 = flops / avg_seconds_v3 / 1.0e9; + + // // std::cout << "-> BM = " << BM_v3 << ", BN = " << BN_v3 + // // << ", BK = " << BK_v3 << "\n"; + // // std::cout << "-> TM = " << TM_v3 << ", TN = " << TN_v3 << "\n"; + // std::cout << "-> Block: " << block_v3.x << " x " << block_v3.y << "\n"; + // std::cout << "-> Grid: " << grid_v3.x << " x " << grid_v3.y << "\n"; + // if (check_gpu_results(M, N, dC, hC_comparison.data())) + // std::cout << "-> Results correct!" << std::endl; + // else + // std::cout << "-> ERROR: Variant 3 faulty!" << std::endl; + // std::cout << "-> average runtime: " << avg_seconds_v3 * 1.0e3 << "ms\n"; + // std::cout << "-> performance: " << gflops_v3 << " GLOP/s\n"; + // std::cout << std::endl; + + // ============================================================ + // Variant 4 launch configuration + // ============================================================ + std::cout << "Variant Benchmark 4: \n"; + std::cout << "-> Description: Shared register tiled with better ccnfig\n"; + experiment_parameters.utilized_implementation_id = 4; + kernel_fusion_benchmark(experiment_parameters, hA, hB, hC, hC_comparison, N, M, K); + std::cout << "Variant Benchmark 4 end \n" << std::endl; + // std::cout << "Variant 4: shared register tiled with better ccnfig\n"; + + // constexpr int BM_v4 = 64; + // constexpr int BN_v4 = 32; + // constexpr int BK_v4 = 16; + // constexpr int TM_v4 = 4; + // constexpr int TN_v4 = 2; + + // dim3 block_v4(BN_v4 / TN_v4, BM_v4 / TM_v4); + // dim3 grid_v4( + // (N + BN_v4 - 1) / BN_v4, + // (M + BM_v4 - 1) / BM_v4 + // ); + // auto start_v4 = std::chrono::high_resolution_clock::now(); + // for (int repeat = 0; repeat < num_repeats; ++repeat) { + // dgemm_v4_shared_register_tiled + // <<>>( + // M, N, K, alpha, dA, dB, beta, dC); + // } + // cudaGetLastError(); + // cudaDeviceSynchronize(); + // auto end_v4 = std::chrono::high_resolution_clock::now(); + // std::chrono::duration elapsed_v4 = end_v4 - start_v4; + // double total_seconds_v4 = elapsed_v4.count(); + // double avg_seconds_v4 = total_seconds_v4 / static_cast(num_repeats); + // double gflops_v4 = flops / avg_seconds_v4 / 1.0e9; + + // // std::cout << "-> BM = " << BM_v4 << ", BN = " << BN_v4 + // // << ", BK = " << BK_v4 << "\n"; + // // std::cout << "-> TM = " << TM_v4 << ", TN = " << TN_v4 << "\n"; + // std::cout << "-> Block: " << block_v4.x << " x " << block_v4.y << "\n"; + // std::cout << "-> Grid: " << grid_v4.x << " x " << grid_v4.y << "\n"; + // if (check_gpu_results(M, N, dC, hC_comparison.data())) + // std::cout << "-> Results correct!" << std::endl; + // else + // std::cout << "-> ERROR: Variant 4 faulty!" << std::endl; + // std::cout << "-> average runtime: " << avg_seconds_v4 * 1.0e3 << "ms\n"; + // std::cout << "-> performance: " << gflops_v4 << " GLOP/s\n"; + // std::cout << std::endl; + + // ============================================================ + // Variant 5 launch configuration + // ============================================================ + std::cout << "Variant Benchmark 5: \n"; + std::cout << "-> Description: Vectorized shared register tiled\n"; + experiment_parameters.utilized_implementation_id = 5; + kernel_fusion_benchmark(experiment_parameters, hA, hB, hC, hC_comparison, N, M, K); + std::cout << "Variant Benchmark 5 end \n" << std::endl; + // std::cout << "Variant 5: vectorized shared register tiled\n"; + + // constexpr int BM_v5 = 64; + // constexpr int BN_v5 = 64; + // constexpr int BK_v5 = 8; + // constexpr int TM_v5 = 4; + // constexpr int TN_v5 = 4; + + // dim3 block_v5(BN_v5 / TN_v5, BM_v5 / TM_v5); + // dim3 grid_v5( + // (N + BN_v5 - 1) / BN_v5, + // (M + BM_v5 - 1) / BM_v5 + // ); + // auto start_v5 = std::chrono::high_resolution_clock::now(); + // for (int repeat = 0; repeat < num_repeats; ++repeat) { + // dgemm_v5_vectorized_register_tiled + // <<>>( + // M, N, K, alpha, dA, dB, beta, dC); + // } + // cudaGetLastError(); + // cudaDeviceSynchronize(); + // auto end_v5 = std::chrono::high_resolution_clock::now(); + // std::chrono::duration elapsed_v5 = end_v5 - start_v5; + // double total_seconds_v5 = elapsed_v5.count(); + // double avg_seconds_v5 = total_seconds_v5 / static_cast(num_repeats); + // double gflops_v5 = flops / avg_seconds_v5 / 1.0e9; + + // // std::cout << "-> BM = " << BM_v5 << ", BN = " << BN_v5 + // // << ", BK = " << BK_v5 << "\n"; + // // std::cout << "-> TM = " << TM_v5 << ", TN = " << TN_v5 << "\n"; + // std::cout << "-> Block: " << block_v5.x << " x " << block_v5.y << "\n"; + // std::cout << "-> Grid: " << grid_v5.x << " x " << grid_v5.y << "\n"; + // if (check_gpu_results(M, N, dC, hC_comparison.data())) + // std::cout << "-> Results correct!" << std::endl; + // else + // std::cout << "-> ERROR: Variant 5 faulty!" << std::endl; + // std::cout << "-> average runtime: " << avg_seconds_v5 * 1.0e3 << "ms\n"; + // std::cout << "-> performance: " << gflops_v5 << " GLOP/s\n"; + // std::cout << std::endl; + + // ============================================================ + // Variant 6 launch configuration + // ============================================================ + std::cout << "Variant Benchmark 6: \n"; + std::cout << "-> Description: WMMA\n"; + experiment_parameters.utilized_implementation_id = 6; + kernel_fusion_benchmark(experiment_parameters, hA, hB, hC, hC_comparison, N, M, K); + std::cout << "Variant Benchmark 6 end \n" << std::endl; + // std::cout << "Variant 6: WMMA\n"; + + // constexpr int BM_v6 = 64; + // constexpr int BN_v6 = 64; + // constexpr int BK_v6 = 16; + + // if ((M % BM_v6) != 0 || (N % BN_v6) != 0 || (K % BK_v6) != 0) { + // fprintf(stderr, + // "dgemm_v5b_fp64_wmma_tiled requires " + // "M %% 64 == 0, N %% 64 == 0, K %% 16 == 0\n"); + // return 0; + // } + + // dim3 block_v6(512); + // dim3 grid_v6(N / BN_v6, M / BM_v6); + + // auto start_v6 = std::chrono::high_resolution_clock::now(); + // for (int repeat = 0; repeat < num_repeats; ++repeat) { + // dgemm_v6_fp64_wmma_tiled<<>>( + // M, N, K, alpha, dA, dB, beta, dC); + // } + // cudaGetLastError(); + // cudaDeviceSynchronize(); + // auto end_v6 = std::chrono::high_resolution_clock::now(); + // std::chrono::duration elapsed_v6 = end_v6 - start_v6; + // double total_seconds_v6 = elapsed_v6.count(); + // double avg_seconds_v6 = total_seconds_v6 / static_cast(num_repeats); + // double gflops_v6 = flops / avg_seconds_v6 / 1.0e9; + // std::cout << "-> Block: " << block_v6.x << " x " << block_v6.y << "\n"; + // std::cout << "-> Grid: " << grid_v6.x << " x " << grid_v6.y << "\n"; + // if (check_gpu_results(M, N, dC, hC_comparison.data())) + // std::cout << "-> Results correct!" << std::endl; + // else + // std::cout << "-> ERROR: Variant 6 faulty!" << std::endl; + // std::cout << "-> average runtime: " << avg_seconds_v6 * 1.0e3 << "ms\n"; + // std::cout << "-> performance: " << gflops_v6 << " GLOP/s\n"; + // std::cout << std::endl; + + // ============================================================ + // Variant 7 launch configuration + // ============================================================ + std::cout << "Variant Benchmark 7: \n"; + std::cout << "-> Description: Better WMMA\n"; + experiment_parameters.utilized_implementation_id = 7; + kernel_fusion_benchmark(experiment_parameters, hA, hB, hC, hC_comparison, N, M, K); + std::cout << "Variant Benchmark 7 end \n" << std::endl; + // std::cout << "Variant 7: Better WMMA\n"; + // constexpr int BM_v7 = 64; + // constexpr int BN_v7 = 64; + // constexpr int BK_v7 = 16; + + // if ((M % BM_v7) != 0 || (N % BN_v7) != 0 || (K % BK_v7) != 0) { + // fprintf(stderr, + // "dgemm_v7_fp64_wmma_64x64_8warps requires " + // "M %% 64 == 0, N %% 64 == 0, K %% 16 == 0\n"); + // return 0; + // } + + // dim3 block_v7(256); + // dim3 grid_v7(N / BN_v7, M / BM_v7); + + // auto start_v7 = std::chrono::high_resolution_clock::now(); + // for (int repeat = 0; repeat < num_repeats; ++repeat) { + // dgemm_v7_fp64_wmma_64x64_8warps + // <<>>(M, N, K, alpha, dA, dB, beta, dC); + // } + // cudaGetLastError(); + // cudaDeviceSynchronize(); + // auto end_v7 = std::chrono::high_resolution_clock::now(); + // std::chrono::duration elapsed_v7 = end_v7 - start_v7; + // double total_seconds_v7 = elapsed_v7.count(); + // double avg_seconds_v7 = total_seconds_v7 / static_cast(num_repeats); + // double gflops_v7 = flops / avg_seconds_v7 / 1.0e9; + // std::cout << "Block: " << block_v7.x << " x " << block_v7.y << "\n"; + // std::cout << "Grid: " << grid_v7.x << " x " << grid_v7.y << "\n"; + // if (check_gpu_results(M, N, dC, hC_comparison.data())) + // std::cout << "-> Results correct!" << std::endl; + // else + // std::cout << "-> ERROR: Variant 7 faulty!" << std::endl; + // std::cout << "-> average runtime: " << avg_seconds_v7 * 1.0e3 << "ms\n"; + // std::cout << "-> performance: " << gflops_v7 << " GLOP/s\n"; + // std::cout << std::endl; + + + std::cout << "Variant Benchmark cublas: \n"; + std::cout << "-> Description: CUBLAS: \n"; + experiment_parameters.utilized_implementation_id = 8; + kernel_fusion_benchmark(experiment_parameters, hA, hB, hC, hC_comparison, N, M, K); + std::cout << "Variant Benchmark cublas end \n" << std::endl; + + //std::cout << "Variant CUBLAS: \n"; + //cublasHandle_t handle; + //cublasCreate(&handle); + //cublasSetMathMode(handle, CUBLAS_TENSOR_OP_MATH); + //// Optional if you use a specific stream: + //// cudaStream_t stream; + //// CHECK_CUDA(cudaStreamCreate(&stream)); + //// CHECK_CUBLAS(cublasSetStream(handle, stream)); + ////cublasSetPointerMode(handle, CUBLAS_POINTER_MODE_HOST); + //// warmup + //cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, &alpha, dB, N, dA, K, &beta, dC, N); + //cudaGetLastError(); + //cudaDeviceSynchronize(); + //auto start_cublas = std::chrono::high_resolution_clock::now(); + //for (int repeat = 0; repeat < num_repeats; ++repeat) { + // // A in cuBLAS sense: B_row interpreted as B^T + // // B in cuBLAS sense: A_row interpreted as A^T + // cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, &alpha, dB, N, dA, K, &beta, dC, N); + //} + //cudaGetLastError(); + //cudaDeviceSynchronize(); + //auto end_cublas = std::chrono::high_resolution_clock::now(); + //std::chrono::duration elapsed_cublas = end_cublas - start_cublas; + //double total_seconds_cublas = elapsed_cublas.count(); + //double avg_seconds_cublas = total_seconds_cublas / static_cast(num_repeats); + //double gflops_cublas = flops / avg_seconds_cublas / 1.0e9; + //std::cout << "Variant CUBLAS: \n"; + //if (check_gpu_results(M, N, dC, hC_comparison.data())) + // std::cout << "-> Results correct!" << std::endl; + //else + // std::cout << "-> ERROR: Variant CUBLAS faulty!" << std::endl; + //std::cout << "-> average runtime: " << avg_seconds_cublas * 1.0e3 << "ms\n"; + //std::cout << "-> performance: " << gflops_cublas << " GLOP/s\n"; + //std::cout << std::endl; + + // ------------------------------------------------------------ + // Cleanup + // ------------------------------------------------------------ + cudaFree(dA); + cudaFree(dB); + cudaFree(dC); + + return hpx::finalize(); +} + +int main(int argc, char *argv[]) { + hpx::init_params p; + p.cfg = {"hpx.commandline.allow_unknown=1"}; + return hpx::init(argc, argv, p); +} diff --git a/benchmarks/cuda_gemm_kernels.hpp b/benchmarks/cuda_gemm_kernels.hpp new file mode 100644 index 00000000..1d19f015 --- /dev/null +++ b/benchmarks/cuda_gemm_kernels.hpp @@ -0,0 +1,860 @@ +// LANL Copyright Notice: +// +//© 2026. Triad National Security, LLC. All rights reserved. +// +// This program was produced under U.S. Government contract 89233218CNA000001 +// for Los Alamos National Laboratory (LANL), which is operated by Triad +// National Security, LLC for the U.S. Department of Energy/National Nuclear +// Security Administration. All rights in the program are reserved by Triad +// National Security, LLC, and the U.S. Department of Energy/National Nuclear +// Security Administration. The Government is granted for itself and others +// acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license +// in this material to reproduce, prepare. derivative works, distribute copies +// to the public, perform publicly and display publicly, and to permit others +// to do so. +// +// +// Distributed under the Boost Software License, Version 1.0. (See accompanying +// file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include +#include + +__global__ void dgemm_naive_unoptimized_v0( + int M, int N, int K, double alpha, const double* A, const double* B, double beta, double* C) { + // Deliberately bad mapping: + // threadIdx.x maps to row instead of column. + int row = blockIdx.x * blockDim.x + threadIdx.x; + int col = blockIdx.y * blockDim.y + threadIdx.y; + if (row < M && col < N) { + double acc = 0.0; + for (int k = 0; k < K; ++k) { + double a = A[row * K + k]; + double b = B[k * N + col]; + acc = acc + a * b; + } + C[row * N + col] = alpha * acc + beta * C[row * N + col]; + } +} +__global__ void dgemm_naive_v1( + int M, int N, int K, double alpha, const double* __restrict__ A, const double* __restrict__ B, + double beta, double* __restrict__ C) { + int row = blockIdx.y * blockDim.y + threadIdx.y; + int col = blockIdx.x * blockDim.x + threadIdx.x; + if (row < M && col < N) { + double acc = 0.0; + for (int k = 0; k < K; ++k) { + double a = A[row * K + k]; + double b = B[k * N + col]; + acc = fma(a, b, acc); + } + int idx = row * N + col; + C[idx] = alpha * acc + beta * C[idx]; + } +} + +template +__global__ void dgemm_shared_tiled_v2( + int M, int N, int K, double alpha, const double* __restrict__ A, const double* __restrict__ B, + double beta, double* __restrict__ C) { + + __shared__ double As[TILE][TILE + 1]; + __shared__ double Bs[TILE][TILE + 1]; + int tx = threadIdx.x; + int ty = threadIdx.y; + int row = blockIdx.y * TILE + ty; + int col = blockIdx.x * TILE + tx; + double acc = 0.0; + + for (int kb = 0; kb < K; kb += TILE) { + int a_col = kb + tx; + int b_row = kb + ty; + if (row < M && a_col < K) { + As[ty][tx] = A[row * K + a_col]; + } else { + As[ty][tx] = 0.0; + } + if (b_row < K && col < N) { + Bs[ty][tx] = B[b_row * N + col]; + } else { + Bs[ty][tx] = 0.0; + } + __syncthreads(); + #pragma unroll + for (int kk = 0; kk < TILE; ++kk) { + acc = fma(As[ty][kk], Bs[kk][tx], acc); + } + + __syncthreads(); + } + if (row < M && col < N) { + int idx = row * N + col; + C[idx] = alpha * acc + beta * C[idx]; + } +} + + +template +__global__ void dgemm_register_tiled_v3( + int M, int N, int K, double alpha, const double* __restrict__ A, const double* __restrict__ B, + double beta, double* __restrict__ C) { + + __shared__ double As[BM][BK + 1]; + __shared__ double Bs[BK][BN + 1]; + int tx = threadIdx.x; + int ty = threadIdx.y; + int tid = ty * blockDim.x + tx; + int num_threads = blockDim.x * blockDim.y; + int block_row = blockIdx.y * BM; + int block_col = blockIdx.x * BN; + int thread_row = ty * TM; + int thread_col = tx * TN; + double acc[TM][TN]; + + #pragma unroll + for (int i = 0; i < TM; ++i) { + #pragma unroll + for (int j = 0; j < TN; ++j) { + acc[i][j] = 0.0; + } + } + + for (int kb = 0; kb < K; kb += BK) { + for (int idx = tid; idx < BM * BK; idx += num_threads) { + int r = idx / BK; + int c = idx % BK; + int global_row = block_row + r; + int global_col = kb + c; + + if (global_row < M && global_col < K) { + As[r][c] = A[global_row * K + global_col]; + } else { + As[r][c] = 0.0; + } + } + + for (int idx = tid; idx < BK * BN; idx += num_threads) { + int r = idx / BN; + int c = idx % BN; + int global_row = kb + r; + int global_col = block_col + c; + + if (global_row < K && global_col < N) { + Bs[r][c] = B[global_row * N + global_col]; + } else { + Bs[r][c] = 0.0; + } + } + __syncthreads(); + + #pragma unroll + for (int kk = 0; kk < BK; ++kk) { + double a_frag[TM]; + double b_frag[TN]; + + #pragma unroll + for (int i = 0; i < TM; ++i) { + a_frag[i] = As[thread_row + i][kk]; + } + #pragma unroll + for (int j = 0; j < TN; ++j) { + b_frag[j] = Bs[kk][thread_col + j]; + } + #pragma unroll + for (int i = 0; i < TM; ++i) { + #pragma unroll + for (int j = 0; j < TN; ++j) { + acc[i][j] = fma(a_frag[i], b_frag[j], acc[i][j]); + } + } + } + __syncthreads(); + } + + #pragma unroll + for (int i = 0; i < TM; ++i) { + int global_row = block_row + thread_row + i; + + #pragma unroll + for (int j = 0; j < TN; ++j) { + int global_col = block_col + thread_col + j; + + if (global_row < M && global_col < N) { + int idx = global_row * N + global_col; + C[idx] = alpha * acc[i][j] + beta * C[idx]; + } + } + } +} + +template < + int BM = 64, // block tile size in M dimension + int BN = 32, // block tile size in N dimension + int BK = 16, // tile size in K dimension + int TM = 4, // per-thread tile size in M + int TN = 2 // per-thread tile size in N +> +__global__ void dgemm_v4_shared_register_tiled( + int M, int N, int K, double alpha, const double* __restrict__ A, const double* __restrict__ B, + double beta, double* __restrict__ C) { + + static_assert(BM % TM == 0, "BM must be divisible by TM"); + static_assert(BN % TN == 0, "BN must be divisible by TN"); + constexpr int blockDimX = BN / TN; + constexpr int blockDimY = BM / TM; + __shared__ double As[BM][BK + 1]; + __shared__ double Bs[BK][BN + 1]; + int tx = threadIdx.x; + int ty = threadIdx.y; + int blockRow = blockIdx.y * BM; + int blockCol = blockIdx.x * BN; + int localRow = ty * TM; + int localCol = tx * TN; + double acc[TM][TN]; + #pragma unroll + for (int i = 0; i < TM; ++i) { + #pragma unroll + for (int j = 0; j < TN; ++j) { + acc[i][j] = 0.0; + } + } + + int tid = ty * blockDim.x + tx; + int numThreads = blockDim.x * blockDim.y; + for (int k0 = 0; k0 < K; k0 += BK) { + for (int idx = tid; idx < BM * BK; idx += numThreads) { + int r = idx / BK; + int k = idx % BK; + int globalRow = blockRow + r; + int globalK = k0 + k; + + if (globalRow < M && globalK < K) { + As[r][k] = A[globalRow * K + globalK]; + } else { + As[r][k] = 0.0; + } + } + + for (int idx = tid; idx < BK * BN; idx += numThreads) { + int k = idx / BN; + int c = idx % BN; + int globalK = k0 + k; + int globalCol = blockCol + c; + + if (globalK < K && globalCol < N) { + Bs[k][c] = B[globalK * N + globalCol]; + } else { + Bs[k][c] = 0.0; + } + } + __syncthreads(); + + #pragma unroll + for (int k = 0; k < BK; ++k) { + double aFrag[TM]; + double bFrag[TN]; + + #pragma unroll + for (int i = 0; i < TM; ++i) { + aFrag[i] = As[localRow + i][k]; + } + #pragma unroll + for (int j = 0; j < TN; ++j) { + bFrag[j] = Bs[k][localCol + j]; + } + #pragma unroll + for (int i = 0; i < TM; ++i) { + #pragma unroll + for (int j = 0; j < TN; ++j) { + acc[i][j] = fma(aFrag[i], bFrag[j], acc[i][j]); + } + } + } + __syncthreads(); + } + + #pragma unroll + for (int i = 0; i < TM; ++i) { + int globalRow = blockRow + localRow + i; + + #pragma unroll + for (int j = 0; j < TN; ++j) { + int globalCol = blockCol + localCol + j; + + if (globalRow < M && globalCol < N) { + int idx = globalRow * N + globalCol; + C[idx] = alpha * acc[i][j] + beta * C[idx]; + } + } + } +} + +template < + int BM = 64, + int BN = 64, + int BK = 8, + int TM = 4, + int TN = 4 +> +__global__ void dgemm_v5_vectorized_register_tiled( + int M, int N, int K, double alpha, const double* __restrict__ A, const double* __restrict__ B, + double beta, double* __restrict__ C) { + + static_assert(BM % TM == 0, "BM must be divisible by TM"); + static_assert(BN % TN == 0, "BN must be divisible by TN"); + static_assert(BK % 2 == 0, "BK must be even for double2 loads"); + static_assert(BN % 2 == 0, "BN must be even for double2 loads"); + constexpr int THREADS_X = BN / TN; // 16 + constexpr int THREADS_Y = BM / TM; // 16 + __shared__ double As[BM][BK]; + __shared__ double Bs[BK][BN]; + const int tx = threadIdx.x; + const int ty = threadIdx.y; + const int tid = ty * blockDim.x + tx; + const int numThreads = blockDim.x * blockDim.y; + const int blockRow = blockIdx.y * BM; + const int blockCol = blockIdx.x * BN; + const int localRow = ty * TM; + const int localCol = tx * TN; + double acc[TM][TN]; + + #pragma unroll + for (int i = 0; i < TM; ++i) { + #pragma unroll + for (int j = 0; j < TN; ++j) { + acc[i][j] = 0.0; + } + } + + const bool canVectorize = ((K & 1) == 0) && ((N & 1) == 0); + + for (int k0 = 0; k0 < K; k0 += BK) { + constexpr int A_VEC_ELEMS = BM * BK / 2; + + for (int idx = tid; idx < A_VEC_ELEMS; idx += numThreads) { + int elem = idx * 2; + int r = elem / BK; + int k = elem - r * BK; + int globalRow = blockRow + r; + int globalK = k0 + k; + + if (canVectorize && globalRow < M && globalK + 1 < K) { + double2 v = *reinterpret_cast( + &A[globalRow * K + globalK] + ); + + *reinterpret_cast(&As[r][k]) = v; + } else { + double v0 = 0.0; + double v1 = 0.0; + + if (globalRow < M && globalK < K) { + v0 = A[globalRow * K + globalK]; + } + + if (globalRow < M && globalK + 1 < K) { + v1 = A[globalRow * K + globalK + 1]; + } + + As[r][k] = v0; + As[r][k + 1] = v1; + } + } + + constexpr int B_VEC_ELEMS = BK * BN / 2; + + for (int idx = tid; idx < B_VEC_ELEMS; idx += numThreads) { + int elem = idx * 2; + int k = elem / BN; + int c = elem - k * BN; + int globalK = k0 + k; + int globalCol = blockCol + c; + + if (canVectorize && globalK < K && globalCol + 1 < N) { + double2 v = *reinterpret_cast( + &B[globalK * N + globalCol] + ); + + *reinterpret_cast(&Bs[k][c]) = v; + } else { + double v0 = 0.0; + double v1 = 0.0; + + if (globalK < K && globalCol < N) { + v0 = B[globalK * N + globalCol]; + } + if (globalK < K && globalCol + 1 < N) { + v1 = B[globalK * N + globalCol + 1]; + } + Bs[k][c] = v0; + Bs[k][c + 1] = v1; + } + } + __syncthreads(); + + #pragma unroll + for (int k = 0; k < BK; ++k) { + double aFrag[TM]; + double bFrag[TN]; + + #pragma unroll + for (int i = 0; i < TM; ++i) { + aFrag[i] = As[localRow + i][k]; + } + #pragma unroll + for (int j = 0; j < TN; ++j) { + bFrag[j] = Bs[k][localCol + j]; + } + #pragma unroll + for (int i = 0; i < TM; ++i) { + #pragma unroll + for (int j = 0; j < TN; ++j) { + acc[i][j] = fma(aFrag[i], bFrag[j], acc[i][j]); + } + } + } + __syncthreads(); + } + + #pragma unroll + for (int i = 0; i < TM; ++i) { + int globalRow = blockRow + localRow + i; + + #pragma unroll + for (int j = 0; j < TN; ++j) { + int globalCol = blockCol + localCol + j; + + if (globalRow < M && globalCol < N) { + int cIdx = globalRow * N + globalCol; + C[cIdx] = alpha * acc[i][j] + beta * C[cIdx]; + } + } + } +} + +using namespace nvcuda; +__global__ void dgemm_v6_fp64_tensorcore_wmma( + int M, int N, int K, double alpha, const double* __restrict__ A, const double* __restrict__ B, + double beta, double* __restrict__ C) { + + constexpr int WMMA_M = 8; + constexpr int WMMA_N = 8; + constexpr int WMMA_K = 4; + constexpr int WARP_SIZE = 32; + int warpIdInBlock = threadIdx.x / WARP_SIZE; + int laneId = threadIdx.x % WARP_SIZE; + int warpsPerBlock = blockDim.x / WARP_SIZE; + int tilesM = M / WMMA_M; + int tilesN = N / WMMA_N; + int globalWarpId = blockIdx.x * warpsPerBlock + warpIdInBlock; + int tileId = globalWarpId; + if (tileId >= tilesM * tilesN) { + return; + } + int tileM = tileId / tilesN; + int tileN = tileId - tileM * tilesN; + int row = tileM * WMMA_M; + int col = tileN * WMMA_N; + + wmma::fragment< + wmma::matrix_a, + WMMA_M, + WMMA_N, + WMMA_K, + double, + wmma::row_major + > aFrag; + + wmma::fragment< + wmma::matrix_b, + WMMA_M, + WMMA_N, + WMMA_K, + double, + wmma::row_major + > bFrag; + + wmma::fragment< + wmma::accumulator, + WMMA_M, + WMMA_N, + WMMA_K, + double + > accFrag; + wmma::fill_fragment(accFrag, 0.0); + + for (int k0 = 0; k0 < K; k0 += WMMA_K) { + const double* aTile = A + row * K + k0; + const double* bTile = B + k0 * N + col; + + wmma::load_matrix_sync(aFrag, aTile, K); + wmma::load_matrix_sync(bFrag, bTile, N); + wmma::mma_sync(accFrag, aFrag, bFrag, accFrag); + } + + if (beta != 0.0) { + wmma::fragment< + wmma::accumulator, + WMMA_M, + WMMA_N, + WMMA_K, + double + > cOldFrag; + const double* cTileIn = C + row * N + col; + wmma::load_matrix_sync(cOldFrag, cTileIn, N, wmma::mem_row_major); + + #pragma unroll + for (int i = 0; i < accFrag.num_elements; ++i) { + accFrag.x[i] = alpha * accFrag.x[i] + beta * cOldFrag.x[i]; + } + } else { + #pragma unroll + for (int i = 0; i < accFrag.num_elements; ++i) { + accFrag.x[i] = alpha * accFrag.x[i]; + } + } + double* cTileOut = C + row * N + col; + wmma::store_matrix_sync(cTileOut, accFrag, N, wmma::mem_row_major); +} + +template < + int BM = 64, + int BN = 64, + int BK = 16 +> +__global__ void dgemm_v6_fp64_wmma_tiled( + int M, int N, int K, double alpha, const double* __restrict__ A, const double* __restrict__ B, + double beta, double* __restrict__ C) { + + constexpr int WMMA_M = 8; + constexpr int WMMA_N = 8; + constexpr int WMMA_K = 4; + constexpr int WARPS_PER_BLOCK = 16; + constexpr int THREADS_PER_BLOCK = WARPS_PER_BLOCK * 32; + static_assert(BM == 64, "This kernel expects BM = 64"); + static_assert(BN == 64, "This kernel expects BN = 64"); + static_assert(BK == 16, "This kernel expects BK = 16"); + + constexpr int AS_LD = BK + 4; // 20 + constexpr int BS_LD = BN + 8; // 72 + __shared__ __align__(16) double As[BM][AS_LD]; + __shared__ __align__(16) double Bs[BK][BS_LD]; + const int tid = threadIdx.x; + const int warpId = tid / 32; + const int blockRow = blockIdx.y * BM; + const int blockCol = blockIdx.x * BN; + + const int warp_m = warpId / 2; // 0..7 + const int warp_n_group = warpId & 1; // 0..1 + const int warp_n_base = warp_n_group * 4; + const int cRow = warp_m * WMMA_M; + + wmma::fragment< + wmma::accumulator, + WMMA_M, + WMMA_N, + WMMA_K, + double + > acc[4]; + + #pragma unroll + for (int t = 0; t < 4; ++t) { + wmma::fill_fragment(acc[t], 0.0); + } + + for (int k0 = 0; k0 < K; k0 += BK) { + constexpr int A_VEC_LOADS = BM * BK / 2; // 512 double2 loads + for (int idx = tid; idx < A_VEC_LOADS; idx += THREADS_PER_BLOCK) { + int linear = idx * 2; + int r = linear / BK; + int k = linear - r * BK; + + const double2 v = *reinterpret_cast( + &A[(blockRow + r) * K + (k0 + k)] + ); + + As[r][k] = v.x; + As[r][k + 1] = v.y; + } + constexpr int B_VEC_LOADS = BK * BN / 2; // 512 double2 loads + for (int idx = tid; idx < B_VEC_LOADS; idx += THREADS_PER_BLOCK) { + int linear = idx * 2; + int k = linear / BN; + int c = linear - k * BN; + + const double2 v = *reinterpret_cast( + &B[(k0 + k) * N + (blockCol + c)] + ); + + Bs[k][c] = v.x; + Bs[k][c + 1] = v.y; + } + __syncthreads(); + + #pragma unroll + for (int kk = 0; kk < BK; kk += WMMA_K) { + wmma::fragment< + wmma::matrix_a, + WMMA_M, + WMMA_N, + WMMA_K, + double, + wmma::row_major + > aFrag; + + wmma::load_matrix_sync( + aFrag, + &As[cRow][kk], + AS_LD + ); + + #pragma unroll + for (int t = 0; t < 4; ++t) { + const int cCol = (warp_n_base + t) * WMMA_N; + + wmma::fragment< + wmma::matrix_b, + WMMA_M, + WMMA_N, + WMMA_K, + double, + wmma::row_major + > bFrag; + + wmma::load_matrix_sync( + bFrag, + &Bs[kk][cCol], + BS_LD + ); + + wmma::mma_sync(acc[t], aFrag, bFrag, acc[t]); + } + } + + __syncthreads(); + } + + #pragma unroll + for (int t = 0; t < 4; ++t) { + const int cCol = (warp_n_base + t) * WMMA_N; + double* cPtr = C + (blockRow + cRow) * N + (blockCol + cCol); + + if (beta != 0.0) { + wmma::fragment< + wmma::accumulator, + WMMA_M, + WMMA_N, + WMMA_K, + double + > cOld; + + wmma::load_matrix_sync( + cOld, + cPtr, + N, + wmma::mem_row_major + ); + + #pragma unroll + for (int i = 0; i < acc[t].num_elements; ++i) { + acc[t].x[i] = alpha * acc[t].x[i] + beta * cOld.x[i]; + } + } else { + #pragma unroll + for (int i = 0; i < acc[t].num_elements; ++i) { + acc[t].x[i] = alpha * acc[t].x[i]; + } + } + wmma::store_matrix_sync( + cPtr, + acc[t], + N, + wmma::mem_row_major + ); + } +} + +template < + int BM = 64, + int BN = 64, + int BK = 16 +> +__global__ __launch_bounds__(256) +void dgemm_v7_fp64_wmma_64x64_8warps( + int M, int N, int K, double alpha, const double* __restrict__ A, const double* __restrict__ B, + double beta, double* __restrict__ C) { + + constexpr int WMMA_M = 8; + constexpr int WMMA_N = 8; + constexpr int WMMA_K = 4; + constexpr int WARPS_PER_BLOCK = 8; + constexpr int THREADS_PER_BLOCK = WARPS_PER_BLOCK * 32; + static_assert(BM == 64, "BM must be 64"); +//static_assert(BN == 64, "BN must be 64"); + static_assert(BK == 16, "BK must be 16"); + constexpr int AS_LD = BK + 4; // 20 + constexpr int BS_LD = BN + 8; // 72 + __shared__ __align__(16) double As[BM][AS_LD]; + __shared__ __align__(16) double Bs[BK][BS_LD]; + const int tid = threadIdx.x; + const int warpId = tid >> 5; + const int blockRow = blockIdx.y * BM; + const int blockCol = blockIdx.x * BN; + /* + 64 x 64 C block = 8 x 8 WMMA tiles. + + 8 warps/block. + + Each warp computes: + 2 WMMA tiles in M + 4 WMMA tiles in N + + So each warp computes a 16 x 32 C region. + */ + const int warp_m_group = warpId >> 1; // 0..3 + const int warp_n_group = warpId & 1; // 0..1 + const int warpRowBase = warp_m_group * 16; // 0, 16, 32, 48 + const int warpColBase = warp_n_group * 32; // 0, 32 + + wmma::fragment< + wmma::accumulator, + WMMA_M, + WMMA_N, + WMMA_K, + double + > acc[2][4]; + + #pragma unroll + for (int i = 0; i < 2; ++i) { + #pragma unroll + for (int j = 0; j < 4; ++j) { + wmma::fill_fragment(acc[i][j], 0.0); + } + } + + for (int k0 = 0; k0 < K; k0 += BK) { + // Load A tile: 64 x 16 + constexpr int A_VEC_LOADS = BM * BK / 2; + for (int idx = tid; idx < A_VEC_LOADS; idx += THREADS_PER_BLOCK) { + int linear = idx * 2; + int r = linear / BK; + int k = linear - r * BK; + + const double2 v = *reinterpret_cast( + &A[(blockRow + r) * K + (k0 + k)] + ); + As[r][k] = v.x; + As[r][k + 1] = v.y; + } + + // Load B tile: 16 x 64 + constexpr int B_VEC_LOADS = BK * BN / 2; + + for (int idx = tid; idx < B_VEC_LOADS; idx += THREADS_PER_BLOCK) { + int linear = idx * 2; + int k = linear / BN; + int c = linear - k * BN; + const double2 v = *reinterpret_cast( + &B[(k0 + k) * N + (blockCol + c)] + ); + Bs[k][c] = v.x; + Bs[k][c + 1] = v.y; + } + __syncthreads(); + + // Compute BK slice using FP64 Tensor Cores. + #pragma unroll + for (int kk = 0; kk < BK; kk += WMMA_K) { + wmma::fragment< + wmma::matrix_a, + WMMA_M, + WMMA_N, + WMMA_K, + double, + wmma::row_major + > aFrag[2]; + + wmma::fragment< + wmma::matrix_b, + WMMA_M, + WMMA_N, + WMMA_K, + double, + wmma::row_major + > bFrag[4]; + + #pragma unroll + for (int i = 0; i < 2; ++i) { + wmma::load_matrix_sync( + aFrag[i], + &As[warpRowBase + i * WMMA_M][kk], + AS_LD + ); + } + + #pragma unroll + for (int j = 0; j < 4; ++j) { + wmma::load_matrix_sync( + bFrag[j], + &Bs[kk][warpColBase + j * WMMA_N], + BS_LD + ); + } + + #pragma unroll + for (int i = 0; i < 2; ++i) { + #pragma unroll + for (int j = 0; j < 4; ++j) { + wmma::mma_sync(acc[i][j], aFrag[i], bFrag[j], acc[i][j]); + } + } + } + __syncthreads(); + } + + // Store C + #pragma unroll + for (int i = 0; i < 2; ++i) { + #pragma unroll + for (int j = 0; j < 4; ++j) { + const int cRow = warpRowBase + i * WMMA_M; + const int cCol = warpColBase + j * WMMA_N; + + double* cPtr = C + (blockRow + cRow) * N + (blockCol + cCol); + + if (beta != 0.0) { + wmma::fragment< + wmma::accumulator, + WMMA_M, + WMMA_N, + WMMA_K, + double + > cOld; + + wmma::load_matrix_sync( + cOld, + cPtr, + N, + wmma::mem_row_major + ); + + #pragma unroll + for (int e = 0; e < acc[i][j].num_elements; ++e) { + acc[i][j].x[e] = alpha * acc[i][j].x[e] + beta * cOld.x[e]; + } + } else { + #pragma unroll + for (int e = 0; e < acc[i][j].num_elements; ++e) { + acc[i][j].x[e] *= alpha; + } + } + wmma::store_matrix_sync( + cPtr, + acc[i][j], + N, + wmma::mem_row_major + ); + } + } +} diff --git a/benchmarks/gemm_cpu_utils.hpp b/benchmarks/gemm_cpu_utils.hpp new file mode 100644 index 00000000..8b784d0b --- /dev/null +++ b/benchmarks/gemm_cpu_utils.hpp @@ -0,0 +1,77 @@ +// LANL Copyright Notice: +// +//© 2026. Triad National Security, LLC. All rights reserved. +// +// This program was produced under U.S. Government contract 89233218CNA000001 +// for Los Alamos National Laboratory (LANL), which is operated by Triad +// National Security, LLC for the U.S. Department of Energy/National Nuclear +// Security Administration. All rights in the program are reserved by Triad +// National Security, LLC, and the U.S. Department of Energy/National Nuclear +// Security Administration. The Government is granted for itself and others +// acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license +// in this material to reproduce, prepare. derivative works, distribute copies +// to the public, perform publicly and display publicly, and to permit others +// to do so. +// +// +// Distributed under the Boost Software License, Version 1.0. (See accompanying +// file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +void dgemm_naive_cpu(int M, int N, int K, double alpha, const double* A, + const double* B, double beta, double* C) { + for (auto row = 0; row < M; row++) { + for (auto column = 0; column < N; column++) { + double acc = 0.0; + for (auto entry = 0; entry < K; entry++) { + acc += A[row * K + entry] * B[column + entry * N]; + } + C[row * N + column] = alpha * acc + beta * C[row * N + column]; + } + } +} + +bool compare_matrices(int M, int N, const double *C_1, const double *C_2, double epsilon = 1e-12) { + bool correct = true; + for (auto entry = 0; entry < M * N; entry ++) { + if (std::abs(C_1[entry] - C_2[entry]) > epsilon) { + correct = false; + std::cerr << "Matrices C_1 and C_2 differ at row " << entry / N << " col " + << entry % N << " with C_1 " << C_1[entry] << " and C_2 " << C_2[entry] + << std::endl; + } + } + return correct; +} + +bool check_gpu_results(int M, int N, double *C_device_results, const double *C_reference) { + std::cout << "checking results ... " << std::endl; + std::vector c_host(M*N); + cudaMemcpy(c_host.data(), C_device_results, M * N * sizeof(double), cudaMemcpyDeviceToHost); + cudaMemset(C_device_results, 0, M * N * sizeof(double)); + return compare_matrices(M, N, c_host.data(), C_reference); +} + +void reset_C_cpu(int M, int N, double *C) { + std::memset(C, 0, M * N * sizeof(double)); +} +void reset_C_gpu(int M, int N, double *C) { + cudaMemset(C, 0, M * N * sizeof(double)); +} + +void fill_random(std::vector& x, std::mt19937_64 &rng) +{ + std::uniform_real_distribution dist(-1.0, 1.0); + + for (double& v : x) { + v = dist(rng); + } +} + +int parse_arg(char** argv, int index, int default_value) +{ + if (argv[index] == nullptr) { + return default_value; + } + + return std::stoi(argv[index]); +} diff --git a/benchmarks/kokkos_gemm_fusion_experiments.cpp b/benchmarks/kokkos_gemm_fusion_experiments.cpp new file mode 100644 index 00000000..e69de29b diff --git a/examples/kernel-aggregation-with-hpx-kokkos.cpp b/examples/kernel-aggregation-with-hpx-kokkos.cpp index a0695b7b..477abc97 100644 --- a/examples/kernel-aggregation-with-hpx-kokkos.cpp +++ b/examples/kernel-aggregation-with-hpx-kokkos.cpp @@ -280,7 +280,8 @@ void launch_gpu_kernel_task( // correctness)) for (size_t i = 0; i < entries_per_task; i++) { if (host_c_slice[i] != 1.0 + 2.0) { - std::cerr << "Task " << task_id << " contained wrong results!!" + std::cerr << "Task " << task_id << " contained wrong results!" + << "Expected: " << 1.0 + 2.0 << " Got: " << host_c_slice[i] << std::endl; break; } diff --git a/examples/recycling-with-hpx-cuda.cu b/examples/recycling-with-hpx-cuda.cpp similarity index 96% rename from examples/recycling-with-hpx-cuda.cu rename to examples/recycling-with-hpx-cuda.cpp index 81670172..05edd33b 100644 --- a/examples/recycling-with-hpx-cuda.cu +++ b/examples/recycling-with-hpx-cuda.cpp @@ -6,13 +6,13 @@ // Developer TODOs regarding CPPuddle usability: // TODO(daissgr) Simplify specifying an executor pool (at least when using the // default round_robin_pool_impl). The current way seems awfully verbose - #include #include -#include + #include #include #include +#include #include #include @@ -28,6 +28,7 @@ #include + /** \file This example shows how to use HPX + CPPuddle with GPU-accelerated * applications. Particulary we focus on how to use a) recycled pinned host * memory, b) recycled device memory, c) the executor pool, d) the HPX-CUDA @@ -52,6 +53,12 @@ using float_t = float; // ... and we will use the HPX CUDA executor inside the executor pool later on using device_executor_t = hpx::cuda::experimental::cuda_executor; +// clang + newer cuda versions cannot deduce the type of cudaMemcpyAsync otherwise... +// deducing is necessary for the post/async calls below +auto memcpy_wrapper(void *dst, const void *src, size_t count, enum cudaMemcpyKind kind, cudaStream_t stream) { + return cudaMemcpyAsync(dst, src, count, kind, stream); +} + /** Just some example CUDA kernel. For simplicity it just adds two vectors. */ __global__ void kernel_add(const float_t *input_a, const float_t *input_b, float_t *output_c) { const int index = blockIdx.x * blockDim.x + threadIdx.x; @@ -126,19 +133,19 @@ void launch_gpu_kernel_task(const size_t task_id, const size_t entries_per_task, entries_per_task); // 4c. Launch data transfers and kernel - hpx::apply(static_cast(executor), cudaMemcpyAsync, + hpx::post(static_cast(executor), memcpy_wrapper, device_a.device_side_buffer, host_a.data(), entries_per_task * sizeof(float_t), cudaMemcpyHostToDevice); - hpx::apply(static_cast(executor), cudaMemcpyAsync, + hpx::post(static_cast(executor), memcpy_wrapper, device_b.device_side_buffer, host_b.data(), entries_per_task * sizeof(float_t), cudaMemcpyHostToDevice); void *args[] = {&device_a.device_side_buffer, &device_b.device_side_buffer, &device_c.device_side_buffer}; - hpx::apply(static_cast(executor), + hpx::post(static_cast(executor), cudaLaunchKernel, kernel_add, entries_per_task / 128, 128, args, 0); auto fut = - hpx::async(static_cast(executor), cudaMemcpyAsync, + hpx::async(static_cast(executor), memcpy_wrapper, host_c.data(), device_c.device_side_buffer, entries_per_task * sizeof(float_t), cudaMemcpyDeviceToHost); fut.get(); // Allow worker thread to jump away until the kernel and diff --git a/include/cppuddle/kernel_aggregation/detail/aggregation_executors_and_allocators.hpp b/include/cppuddle/kernel_aggregation/detail/aggregation_executors_and_allocators.hpp index b3b582cc..cb0fb05e 100644 --- a/include/cppuddle/kernel_aggregation/detail/aggregation_executors_and_allocators.hpp +++ b/include/cppuddle/kernel_aggregation/detail/aggregation_executors_and_allocators.hpp @@ -1,4 +1,18 @@ -// Copyright (c) 2022-2024 Gregor Daiß +// LANL Copyright Notice: +// +//© 2026. Triad National Security, LLC. All rights reserved. +// +// This program was produced under U.S. Government contract 89233218CNA000001 +// for Los Alamos National Laboratory (LANL), which is operated by Triad +// National Security, LLC for the U.S. Department of Energy/National Nuclear +// Security Administration. All rights in the program are reserved by Triad +// National Security, LLC, and the U.S. Department of Energy/National Nuclear +// Security Administration. The Government is granted for itself and others +// acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license +// in this material to reproduce, prepare. derivative works, distribute copies +// to the public, perform publicly and display publicly, and to permit others +// to do so. +// // // Distributed under the Boost Software License, Version 1.0. (See accompanying // file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) @@ -557,7 +571,7 @@ template class aggregated_executor { //=============================================================================== - hpx::promise slices_full_promise; + // hpx::promise slices_full_promise; /// Promises with the slice executors -- to be set when the starting criteria /// is met std::vector> executor_slices; @@ -837,10 +851,12 @@ template class aggregated_executor { executor_slices_alive = true; buffers_in_use = false; dealloc_counter = 0; + launched_slices = 0; - if (mode == aggregated_executor_modes::STRICT ) { - slices_full_promise = hpx::promise{}; - } + + // if (mode == aggregated_executor_modes::STRICT ) { + // slices_full_promise = hpx::promise{}; + // } } // Create Executor Slice future -- that will be returned later @@ -883,14 +899,10 @@ template class aggregated_executor { auto exec_fut = (*executor_wrapper).get_future(); /* auto fut = hpx::when_any(exec_fut, slices_full_fut); */ fut = std::move(exec_fut); - } else { - auto slices_full_fut = slices_full_promise.get_shared_future(); - // Just use the slices launch condition - fut = std::move(slices_full_fut); - } // Launch all executor slices within this continuation current_continuation = fut.then([this](auto &&fut) { std::lock_guard guard(mut); + // if (!slices_exhausted) { slices_exhausted = true; launched_slices = current_slices; size_t id = 0; @@ -901,23 +913,19 @@ template class aggregated_executor { } executor_slices.clear(); }); + } } if (local_slice_id >= max_slices && mode != aggregated_executor_modes::ENDLESS) { - slices_exhausted = true; // prevents any more threads from entering - // before the continuation is launched - /* launched_slices = current_slices; */ - /* size_t id = 0; */ - /* for (auto &slice_promise : executor_slices) { */ - /* slice_promise.set_value( */ - /* executor_slice{*this, id, launched_slices}); */ - /* id++; */ - /* } */ - /* executor_slices.clear(); */ - if (mode == aggregated_executor_modes::STRICT ) { - slices_full_promise.set_value(); // Trigger slices launch condition continuation - } - // that continuation will set all executor slices so far handed out to ready + slices_exhausted = true; + launched_slices = current_slices; + size_t id = 0; + for (auto &slice_promise : executor_slices) { + slice_promise.set_value( + executor_slice{*this, id, launched_slices, max_slices}); + id++; + } + executor_slices.clear(); } return ret_fut; } else { @@ -937,6 +945,8 @@ template class aggregated_executor { // Last slice goes out scope? if (local_slice_id == 0) { // Mark executor fit for reusage + if (mode != aggregated_executor_modes::STRICT ) + current_continuation.get(); // clear remaining task in case it was not run yet std::lock_guard guard(mut); executor_slices_alive = false; if (!executor_slices_alive && !buffers_in_use) { @@ -952,9 +962,9 @@ template class aggregated_executor { assert(executor_slices_alive == false); assert(buffers_in_use == false); - if (mode != aggregated_executor_modes::STRICT ) { - slices_full_promise.set_value(); // Trigger slices launch condition continuation - } + // if (mode != aggregated_executor_modes::STRICT ) { + // slices_full_promise.set_value(); // Trigger slices launch condition continuation + // } // Cleanup leftovers from last run if any function_calls.clear(); diff --git a/include/cppuddle/kernel_aggregation/util/kokkos_aggregation_util.hpp b/include/cppuddle/kernel_aggregation/util/kokkos_aggregation_util.hpp index bf7aacf3..db940b8e 100644 --- a/include/cppuddle/kernel_aggregation/util/kokkos_aggregation_util.hpp +++ b/include/cppuddle/kernel_aggregation/util/kokkos_aggregation_util.hpp @@ -1,4 +1,18 @@ -// Copyright (c) 2022-2024 Gregor Daiß +// LANL Copyright Notice: +// +//© 2026. Triad National Security, LLC. All rights reserved. +// +// This program was produced under U.S. Government contract 89233218CNA000001 +// for Los Alamos National Laboratory (LANL), which is operated by Triad +// National Security, LLC for the U.S. Department of Energy/National Nuclear +// Security Administration. All rights in the program are reserved by Triad +// National Security, LLC, and the U.S. Department of Energy/National Nuclear +// Security Administration. The Government is granted for itself and others +// acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license +// in this material to reproduce, prepare. derivative works, distribute copies +// to the public, perform publicly and display publicly, and to permit others +// to do so. +// // // Distributed under the Boost Software License, Version 1.0. (See accompanying // file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) @@ -105,7 +119,9 @@ template < typename... Args> CPPUDDLE_HOST_DEVICE_METHOD auto map_views_to_slice(const Agg_executor_t &agg_exec, - const Agg_view_t ¤t_arg, const Args &...rest) { + const Agg_view_t ¤t_arg, const Args &...rest) -> + std::tuple { const size_t slice_id = agg_exec.id; const size_t max_slices = agg_exec.max_slices; static_assert(Kokkos::is_view::value, diff --git a/include/cppuddle/memory_recycling/cuda_recycling_allocators.hpp b/include/cppuddle/memory_recycling/cuda_recycling_allocators.hpp index 6cc9c565..6fe62423 100644 --- a/include/cppuddle/memory_recycling/cuda_recycling_allocators.hpp +++ b/include/cppuddle/memory_recycling/cuda_recycling_allocators.hpp @@ -1,4 +1,18 @@ -// Copyright (c) 2020-2024 Gregor Daiß +// LANL Copyright Notice: +// +//© 2026. Triad National Security, LLC. All rights reserved. +// +// This program was produced under U.S. Government contract 89233218CNA000001 +// for Los Alamos National Laboratory (LANL), which is operated by Triad +// National Security, LLC for the U.S. Department of Energy/National Nuclear +// Security Administration. All rights in the program are reserved by Triad +// National Security, LLC, and the U.S. Department of Energy/National Nuclear +// Security Administration. The Government is granted for itself and others +// acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license +// in this material to reproduce, prepare. derivative works, distribute copies +// to the public, perform publicly and display publicly, and to permit others +// to do so. +// // // Distributed under the Boost Software License, Version 1.0. (See accompanying // file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) @@ -7,7 +21,7 @@ #define CUDA_RECYCLING_ALLOCATORS_HPP #include "buffer_management_interface.hpp" -// import cuda_pinned_allocator and cuda_device_allocator +// import underlying cuda allocators used here #include "detail/cuda_underlying_allocators.hpp" /// \file @@ -41,6 +55,10 @@ template ::value, int> = 0> using recycle_allocator_cuda_device = detail::recycle_allocator>; +template ::value, int> = 0> +using recycle_allocator_cuda_managed = + detail::aggressive_recycle_allocator>; + } // namespace memory_recycling } // end namespace cppuddle #endif diff --git a/include/cppuddle/memory_recycling/detail/buffer_management.hpp b/include/cppuddle/memory_recycling/detail/buffer_management.hpp index 7c30c781..d979da49 100644 --- a/include/cppuddle/memory_recycling/detail/buffer_management.hpp +++ b/include/cppuddle/memory_recycling/detail/buffer_management.hpp @@ -1,4 +1,18 @@ -// Copyright (c) 2020-2024 Gregor Daiß +// LANL Copyright Notice: +// +//© 2026. Triad National Security, LLC. All rights reserved. +// +// This program was produced under U.S. Government contract 89233218CNA000001 +// for Los Alamos National Laboratory (LANL), which is operated by Triad +// National Security, LLC for the U.S. Department of Energy/National Nuclear +// Security Administration. All rights in the program are reserved by Triad +// National Security, LLC, and the U.S. Department of Energy/National Nuclear +// Security Administration. The Government is granted for itself and others +// acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license +// in this material to reproduce, prepare. derivative works, distribute copies +// to the public, perform publicly and display publicly, and to permit others +// to do so. +// // // Distributed under the Boost Software License, Version 1.0. (See accompanying // file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) @@ -190,7 +204,7 @@ For better performance configure CPPuddle with CPPUDDLE_WITH_BUFFER_RECYCLING=ON /// one buffer_manager std::list> finalize_callbacks; /// Callbacks for buffer_manager cleanups - each callback destroys all buffers within - /// one buffer_manager, both used and unsued + /// one buffer_manager, both used and unused std::list> total_cleanup_callbacks; /// Callbacks for partial buffer_manager cleanups - each callback deallocates /// all unused buffers of a manager @@ -268,14 +282,17 @@ For better performance configure CPPuddle with CPPUDDLE_WITH_BUFFER_RECYCLING=ON assert(instance() && !is_finalized); for (auto i = 0; i < number_instances * max_number_gpus; i++) { std::lock_guard guard(instance()[i].mut); - for (auto &buffer_tuple : instance()[i].unused_buffer_list) { - Host_Allocator alloc; - if (std::get<3>(buffer_tuple)) { - std::destroy_n(std::get<0>(buffer_tuple), std::get<1>(buffer_tuple)); + for (auto &unused_buffer_list_tuple : instance()[i].unused_buffer_map) { + auto &unused_buffer_list = std::get<1>(unused_buffer_list_tuple); + for (auto &buffer_tuple : unused_buffer_list) { + Host_Allocator alloc; + if (std::get<3>(buffer_tuple)) { + std::destroy_n(std::get<0>(buffer_tuple), std::get<1>(buffer_tuple)); + } + alloc.deallocate(std::get<0>(buffer_tuple), std::get<1>(buffer_tuple)); } - alloc.deallocate(std::get<0>(buffer_tuple), std::get<1>(buffer_tuple)); + unused_buffer_list.clear(); } - instance()[i].unused_buffer_list.clear(); } } #if defined(CPPUDDLE_HAVE_COUNTERS) && defined(CPPUDDLE_HAVE_HPX) @@ -383,17 +400,19 @@ For better performance configure CPPuddle with CPPUDDLE_WITH_BUFFER_RECYCLING=ON location_id = location_id + device_id * number_instances; std::lock_guard guard(instance()[location_id].mut); + std::list &unused_buffer_list = instance()[location_id].unused_buffer_map[number_of_elements]; + #ifdef CPPUDDLE_HAVE_COUNTERS instance()[location_id].number_allocation++; sum_number_allocation++; #endif // Check for unused buffers we can recycle: - for (auto iter = instance()[location_id].unused_buffer_list.begin(); - iter != instance()[location_id].unused_buffer_list.end(); iter++) { + for (auto iter = unused_buffer_list.begin(); + iter != unused_buffer_list.end(); iter++) { auto tuple = *iter; if (std::get<1>(tuple) == number_of_elements) { - instance()[location_id].unused_buffer_list.erase(iter); + unused_buffer_list.erase(iter); // handle the switch from aggressive to non aggressive reusage (or // vice-versa) @@ -502,7 +521,7 @@ For better performance configure CPPuddle with CPPUDDLE_WITH_BUFFER_RECYCLING=ON // sanity checks: assert(std::get<1>(tuple) == number_of_elements); // move to the unused_buffer list - instance()[location_id].unused_buffer_list.push_front(tuple); + instance()[location_id].unused_buffer_map.at(std::get<1>(tuple)).push_front(tuple); instance()[location_id].buffer_map.erase(memory_location); return; // Success } @@ -534,7 +553,7 @@ For better performance configure CPPuddle with CPPUDDLE_WITH_BUFFER_RECYCLING=ON // sanity checks: assert(std::get<1>(tuple) == number_of_elements); // move to the unused_buffer list - instance()[location_id].unused_buffer_list.push_front(tuple); + instance()[location_id].unused_buffer_map.at(std::get<1>(tuple)).push_front(tuple); instance()[location_id].buffer_map.erase(memory_location); return; // Success } @@ -568,7 +587,7 @@ For better performance configure CPPuddle with CPPUDDLE_WITH_BUFFER_RECYCLING=ON // sanity checks: assert(std::get<1>(tuple) == number_of_elements); // move to the unused_buffer list - instance()[location_id].unused_buffer_list.push_front(tuple); + instance()[location_id].unused_buffer_map.at(std::get<1>(tuple)).push_front(tuple); instance()[location_id].buffer_map.erase(memory_location); return; // Success } @@ -594,7 +613,7 @@ For better performance configure CPPuddle with CPPUDDLE_WITH_BUFFER_RECYCLING=ON // sanity checks: assert(std::get<1>(tuple) == number_of_elements); // move to the unused_buffer list - instance()[location_id].unused_buffer_list.push_front(tuple); + instance()[location_id].unused_buffer_map.at(std::get<1>(tuple)).push_front(tuple); instance()[location_id].buffer_map.erase(memory_location); return; // Success } @@ -621,8 +640,10 @@ For better performance configure CPPuddle with CPPUDDLE_WITH_BUFFER_RECYCLING=ON private: /// List with all buffers still in usage std::unordered_map buffer_map{}; + /// Map to lists with all buffers currently not in use + std::unordered_map> unused_buffer_map{}; /// List with all buffers currently not used - std::list unused_buffer_list{}; + // std::list unused_buffer_list{}; /// Access control mutex_t mut; #ifdef CPPUDDLE_HAVE_COUNTERS @@ -675,7 +696,8 @@ For better performance configure CPPuddle with CPPUDDLE_WITH_BUFFER_RECYCLING=ON if (number_allocation == 0) return; // Print performance counters - size_t number_cleaned = unused_buffer_list.size() + buffer_map.size(); + // size_t number_cleaned = unused_buffer_list.size() + buffer_map.size(); + // std::cout << "\nBuffer manager destructor for (Alloc: " << boost::core::demangle(typeid(Host_Allocator).name()) << ", Type: " << boost::core::demangle(typeid(T).name()) @@ -696,7 +718,7 @@ For better performance configure CPPuddle with CPPUDDLE_WITH_BUFFER_RECYCLING=ON << number_creation << std::endl << "--> Number cleaned up buffers: " " " - << number_cleaned << std::endl + << number_creation << std::endl // TODO Add actual deallocation counter across map << "--> Number wrong deallocation hints: " " " << number_wrong_hints << std::endl @@ -715,16 +737,20 @@ For better performance configure CPPuddle with CPPUDDLE_WITH_BUFFER_RECYCLING=ON #ifdef CPPUDDLE_HAVE_COUNTERS if (number_allocation == 0 && number_recycling == 0 && number_bad_alloc == 0 && number_creation == 0 && - unused_buffer_list.empty() && buffer_map.empty()) { + unused_buffer_map.empty() && buffer_map.empty()) { return; } #endif - for (auto &buffer_tuple : unused_buffer_list) { - Host_Allocator alloc; - if (std::get<3>(buffer_tuple)) { - std::destroy_n(std::get<0>(buffer_tuple), std::get<1>(buffer_tuple)); + for (auto &unused_buffer_list_tuple : unused_buffer_map) { + auto &unused_buffer_list = std::get<1>(unused_buffer_list_tuple); + for (auto &buffer_tuple : unused_buffer_list) { + Host_Allocator alloc; + if (std::get<3>(buffer_tuple)) { + std::destroy_n(std::get<0>(buffer_tuple), std::get<1>(buffer_tuple)); + } + alloc.deallocate(std::get<0>(buffer_tuple), std::get<1>(buffer_tuple)); } - alloc.deallocate(std::get<0>(buffer_tuple), std::get<1>(buffer_tuple)); + unused_buffer_list.clear(); } for (auto &map_tuple : buffer_map) { auto buffer_tuple = map_tuple.second; @@ -734,7 +760,6 @@ For better performance configure CPPuddle with CPPUDDLE_WITH_BUFFER_RECYCLING=ON } alloc.deallocate(std::get<0>(buffer_tuple), std::get<1>(buffer_tuple)); } - unused_buffer_list.clear(); buffer_map.clear(); #ifdef CPPUDDLE_HAVE_COUNTERS number_allocation = 0; diff --git a/include/cppuddle/memory_recycling/detail/cuda_underlying_allocators.hpp b/include/cppuddle/memory_recycling/detail/cuda_underlying_allocators.hpp index ab1f8681..c98efbe7 100644 --- a/include/cppuddle/memory_recycling/detail/cuda_underlying_allocators.hpp +++ b/include/cppuddle/memory_recycling/detail/cuda_underlying_allocators.hpp @@ -1,4 +1,18 @@ -// Copyright (c) 2020-2024 Gregor Daiß +// LANL Copyright Notice: +// +//© 2026. Triad National Security, LLC. All rights reserved. +// +// This program was produced under U.S. Government contract 89233218CNA000001 +// for Los Alamos National Laboratory (LANL), which is operated by Triad +// National Security, LLC for the U.S. Department of Energy/National Nuclear +// Security Administration. All rights in the program are reserved by Triad +// National Security, LLC, and the U.S. Department of Energy/National Nuclear +// Security Administration. The Government is granted for itself and others +// acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license +// in this material to reproduce, prepare. derivative works, distribute copies +// to the public, perform publicly and display publicly, and to permit others +// to do so. +// // // Distributed under the Boost Software License, Version 1.0. (See accompanying // file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) @@ -9,6 +23,7 @@ #include #include #include +#include namespace cppuddle { namespace memory_recycling { @@ -94,6 +109,45 @@ constexpr bool operator!=(cuda_device_allocator const &, cuda_device_allocator const &) noexcept { return false; } +/// Underlying allocator for CUDA managed memory +template struct cuda_managed_memory_allocator { + using value_type = T; + cuda_managed_memory_allocator() noexcept = default; + template + explicit cuda_managed_memory_allocator(cuda_managed_memory_allocator const &) noexcept {} + T *allocate(std::size_t n) { + T *data; + cudaError_t error = cudaMallocManaged(&data, n * sizeof(T)); + if (error != cudaSuccess) { + std::string msg = + std::string( + "cuda_managed_memory_allocator failed due to cudaMallocManaged failure : ") + + std::string(cudaGetErrorString(error)); + throw std::runtime_error(msg); + } + return data; + } + void deallocate(T *p, std::size_t n) { + cudaError_t error = cudaFree(p); + if (error != cudaSuccess) { + std::string msg = + std::string( + "cuda_managed_memory_allocator failed due to cudaFree failure : ") + + std::string(cudaGetErrorString(error)); + throw std::runtime_error(msg); + } + } +}; +template +constexpr bool operator==(cuda_managed_memory_allocator const &, + cuda_managed_memory_allocator const &) noexcept { + return true; +} +template +constexpr bool operator!=(cuda_managed_memory_allocator const &, + cuda_managed_memory_allocator const &) noexcept { + return false; +} } // end namespace detail } // namespace memory_recycling } // end namespace cppuddle