From 9290aba6fad3500ca9dc38c10804c55a4ac61a24 Mon Sep 17 00:00:00 2001 From: TruNcat3 <31944766+TruNcat3@users.noreply.github.com> Date: Tue, 11 Aug 2026 14:01:13 +0800 Subject: [PATCH 1/2] Add v0.6 general-shape plan API --- CHANGELOG.md | 42 + CMakeLists.txt | 28 +- README.md | 43 +- ROADMAP.md | 24 + apps/cubutterfly_plan_bench.cpp | 320 +++ config/v100_application_profiles.json | 20 + docs/README.md | 9 +- docs/api_reference.md | 17 +- docs/c_api.md | 115 ++ docs/cubutterfly_feature_coverage.md | 21 +- docs/experiments.md | 11 + docs/general_shape_results.md | 97 + docs/general_shapes.md | 76 + docs/getting_started.md | 29 + docs/programming_guide.md | 6 + docs/research_status.md | 15 +- docs/support_matrix.md | 14 +- docs/v0.6_implementation_status.md | 101 + include/cubutterfly/cubutterfly.h | 188 ++ include/cuntt/butterfly.hpp | 8 + results/README.md | 21 + ...general_shape_direct_boundary_v100_raw.csv | 26 + .../general_shape_fused_input_v100_raw.csv | 12 + results/general_shape_selection_v100_raw.csv | 51 + results/general_shapes_v100_smoke.csv | 8 + scripts/generate_application_profiles.py | 61 + scripts/generate_design_points.py | 28 + src/bluestein_fft.cu | 270 +++ src/bluestein_fft.hpp | 32 + src/bluestein_ntt.cu | 303 +++ src/bluestein_ntt.hpp | 35 + src/butterfly.cu | 49 + src/c_api.cpp | 1753 +++++++++++++++++ src/direct_fft.cu | 114 ++ src/direct_fft.hpp | 27 + src/general_kernels.cu | 183 ++ src/general_kernels.hpp | 39 + src/generated_unit_api.cuh | 12 + src/register_fwht.cuh | 104 +- tests/test_application_profiles.py | 37 + tests/test_butterfly.cpp | 55 + tests/test_c_api.c | 689 +++++++ 42 files changed, 5055 insertions(+), 38 deletions(-) create mode 100644 apps/cubutterfly_plan_bench.cpp create mode 100644 config/v100_application_profiles.json create mode 100644 docs/c_api.md create mode 100644 docs/general_shape_results.md create mode 100644 docs/general_shapes.md create mode 100644 docs/v0.6_implementation_status.md create mode 100644 include/cubutterfly/cubutterfly.h create mode 100644 results/general_shape_direct_boundary_v100_raw.csv create mode 100644 results/general_shape_fused_input_v100_raw.csv create mode 100644 results/general_shape_selection_v100_raw.csv create mode 100644 results/general_shapes_v100_smoke.csv create mode 100644 scripts/generate_application_profiles.py create mode 100644 src/bluestein_fft.cu create mode 100644 src/bluestein_fft.hpp create mode 100644 src/bluestein_ntt.cu create mode 100644 src/bluestein_ntt.hpp create mode 100644 src/c_api.cpp create mode 100644 src/direct_fft.cu create mode 100644 src/direct_fft.hpp create mode 100644 src/general_kernels.cu create mode 100644 src/general_kernels.hpp create mode 100644 tests/test_application_profiles.py create mode 100644 tests/test_c_api.c diff --git a/CHANGELOG.md b/CHANGELOG.md index 096652a..ccc32ff 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,48 @@ All notable repository and research-artifact changes are recorded here. +## 0.6.0 - Unreleased General Shapes And Plan API + +### Added + +- CUDA-library-style C handle/descriptor/plan API with streams, explicit + workspace, status codes, selection queries, and structured logging. +- Exact arbitrary-length complex FP32/FP64 FFT through Bluestein and exact + uint64 modular Bluestein NTT when the requested root domain exists. +- Explicit zero-extended power-of-two embedding for all public operators. +- Rank-two FFT, NTT, FWHT, zeta/Mobius, and Structured 2x2 composition through + per-axis plans and transpose boundaries. +- Versioned V100 application profile table generated at build time, exact-key + user cache, bounded runtime measurement policy, and profile-miss guidance. +- A `cubutterfly_plan_bench` CLI for logical/physical shape, batch, policy, + cache, workspace, and kernel-time inspection. +- Per-axis static/cache/measurement selection for embedded and rank-two + compositions, selectable Hybrid2D Bluestein NTT cores, and complete direct + cuFFT versus Bluestein measurement for exact arbitrary FFT. +- Structured 2x2 matrices and explicit composite algorithm IDs in the plan + benchmark surface. +- Direct physical-output lowering for contiguous embedded and rank-two + compositions, eliminating the final full-buffer scatter while retaining the + general strided-output path. +- Generated zero-extended input dispatch for warp-register FP32 FWHT and + Structured 2x2 units, with the runtime selecting the one-kernel fused-input + lowering for validated FWHT workloads and retaining pack fallback otherwise. + +### Current Evidence Boundary + +- All 49 CTest entries pass on V100/CUDA 11.8, including exact arbitrary + length, embedding, rank-two, strided fallback, cache, and generated boundary + candidate coverage. +- Direct packed arbitrary FFT selects cuFFT at timing parity; this is + processing-unit abstention, not a claim of improving the vendor core. +- Fused-input FWHT reaches 0.11759 ms for `32767 -> 32768,batch=256`, within + 1.01x of the archived pre-padded Dao result, and uses zero composition + workspace. +- Structured fused input is retained as a correct generated candidate but is + not selected because the measured V100 path regresses. +- Cross-GPU selection transfer and universal specialist-library superiority + remain outside the release claim. + ## 0.5.0 - 2026-08-08 Numeric-Regime Mapping ### Added diff --git a/CMakeLists.txt b/CMakeLists.txt index f700fdb..dff0ef5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -4,7 +4,7 @@ if(NOT DEFINED CMAKE_CUDA_ARCHITECTURES) set(CMAKE_CUDA_ARCHITECTURES 70 CACHE STRING "CUDA architectures to compile for") endif() -project(cuButterfly VERSION 0.5.0 LANGUAGES CXX CUDA) +project(cuButterfly VERSION 0.6.0 LANGUAGES C CXX CUDA) include(GNUInstallDirs) include(CMakePackageConfigHelpers) @@ -144,11 +144,25 @@ execute_process( if(NOT CUBUTTERFLY_RUNTIME_SELECTOR_RESULT EQUAL 0) message(FATAL_ERROR "cuButterfly runtime selector generation failed: ${CUBUTTERFLY_RUNTIME_SELECTOR_ERROR}") endif() +execute_process( + COMMAND "${Python3_EXECUTABLE}" "${CMAKE_CURRENT_SOURCE_DIR}/scripts/generate_application_profiles.py" + --spec "${CMAKE_CURRENT_SOURCE_DIR}/config/v100_application_profiles.json" + --header "${CUNTT_GENERATED_DIR}/generated_application_profiles.hpp" + RESULT_VARIABLE CUBUTTERFLY_APPLICATION_PROFILE_RESULT + ERROR_VARIABLE CUBUTTERFLY_APPLICATION_PROFILE_ERROR) +if(NOT CUBUTTERFLY_APPLICATION_PROFILE_RESULT EQUAL 0) + message(FATAL_ERROR "cuButterfly application profile generation failed: ${CUBUTTERFLY_APPLICATION_PROFILE_ERROR}") +endif() include("${CUNTT_GENERATED_DIR}/design_points.cmake") add_library(cuntt STATIC + src/bluestein_fft.cu + src/bluestein_ntt.cu + src/c_api.cpp + src/direct_fft.cu src/butterfly.cu src/butterfly_reference.cpp + src/general_kernels.cu src/mapping_selector.cpp src/ntt.cu src/reference.cpp) @@ -225,6 +239,10 @@ add_executable(cubutterfly_bench apps/cubutterfly_bench.cpp) target_link_libraries(cubutterfly_bench PRIVATE cuntt) target_compile_options(cubutterfly_bench PRIVATE -Wall -Wextra -Wpedantic) +add_executable(cubutterfly_plan_bench apps/cubutterfly_plan_bench.cpp) +target_link_libraries(cubutterfly_plan_bench PRIVATE cuntt) +target_compile_options(cubutterfly_plan_bench PRIVATE -Wall -Wextra -Wpedantic) + add_executable(cuntt_tests tests/test_ntt.cpp) target_link_libraries(cuntt_tests PRIVATE cuntt) target_compile_options(cuntt_tests PRIVATE -Wall -Wextra -Wpedantic) @@ -236,6 +254,9 @@ target_compile_options(cubutterfly_tests PRIVATE -Wall -Wextra -Wpedantic) add_executable(cubutterfly_structured_reference_tests tests/test_structured_reference.cpp) target_link_libraries(cubutterfly_structured_reference_tests PRIVATE cuntt) target_compile_options(cubutterfly_structured_reference_tests PRIVATE -Wall -Wextra -Wpedantic) +add_executable(cubutterfly_c_api_tests tests/test_c_api.c) +target_link_libraries(cubutterfly_c_api_tests PRIVATE cuntt) +target_compile_options(cubutterfly_c_api_tests PRIVATE -Wall -Wextra -Wpedantic) if(CUBUTTERFLY_ENABLE_CUFFTDX) target_compile_definitions(cubutterfly_tests PRIVATE CUBUTTERFLY_TEST_CUFFTDX=1) endif() @@ -267,6 +288,7 @@ if(BUILD_TESTING) add_test(NAME cuntt_correctness COMMAND cuntt_tests) add_test(NAME cubutterfly_correctness COMMAND cubutterfly_tests) add_test(NAME cubutterfly_structured_reference COMMAND cubutterfly_structured_reference_tests) + add_test(NAME cubutterfly_c_api COMMAND cubutterfly_c_api_tests) add_test(NAME fft_design_space COMMAND "${Python3_EXECUTABLE}" "${CMAKE_CURRENT_SOURCE_DIR}/tests/test_fft_design_space.py") add_test(NAME butterfly_design_space @@ -312,6 +334,8 @@ if(BUILD_TESTING) -P "${CMAKE_CURRENT_SOURCE_DIR}/tests/test_install_package.cmake") add_test(NAME runtime_selector_generator COMMAND "${Python3_EXECUTABLE}" "${CMAKE_CURRENT_SOURCE_DIR}/tests/test_runtime_selector_generator.py") + add_test(NAME application_profile_generator + COMMAND "${Python3_EXECUTABLE}" "${CMAKE_CURRENT_SOURCE_DIR}/tests/test_application_profiles.py") if(CUBUTTERFLY_ENABLE_CUFFTDX) add_test(NAME fft_fp64_cufftdx_block_correctness COMMAND cubutterfly_bench --operator fft --precision fp64 --backend temporal-tile @@ -450,7 +474,7 @@ install(TARGETS cuntt ARCHIVE DESTINATION "${CMAKE_INSTALL_LIBDIR}" LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}") -install(DIRECTORY include/cuntt DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}") +install(DIRECTORY include/cuntt include/cubutterfly DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}") configure_package_config_file( "${CMAKE_CURRENT_SOURCE_DIR}/cmake/cuButterflyConfig.cmake.in" diff --git a/README.md b/README.md index 0879b59..206b4f8 100644 --- a/README.md +++ b/README.md @@ -18,6 +18,27 @@ processing-unit boundary, same-machine library comparisons, Nsight profiling scripts, and the raw CSV data used in the reports. It is a research artifact, not a drop-in replacement for cuFFT or a production cryptography library. +The v0.6 development API adds a CUDA-library-style C plan surface, exact +non-power-of-two FFT/compatible NTT lowering, explicit power-of-two embedding, +rank-two axis mapping, caller-owned workspace, static application profiles, +one-time measured tuning, and exact-key cache reuse. See the +[C Plan API](docs/c_api.md) and [General Shape Mapping](docs/general_shapes.md). +The [general-shape comparison](docs/general_shape_results.md) records both the +initial unfused lowering and the new composition selector. Exact packed FFT +now selects direct cuFFT at parity, while modular NTT and embedded operators +reuse selected physical cores. Composition boundaries are reported explicitly: +contiguous output removes scatter, selected warp-register FWHT also fuses +logical input and zero fill, and all other paths retain their documented +materialized boundary. +Contiguous embedded outputs also bypass the final scatter and expose that +choice as `+direct-output` in the selected algorithm name. +Generated FP32 warp-register FWHT can also consume the logical input and +synthesize its zero tail inside the physical kernel. This one-kernel lowering +is exposed as `+fused-input`; unsupported or measured-slower combinations keep +the explicit pack path. +The complete release boundary is summarized in the +[v0.6 Implementation Status](docs/v0.6_implementation_status.md). + ## Why This Design A fast butterfly kernel needs more than a fast butterfly instruction. It must @@ -104,7 +125,10 @@ The historical `xor-zeta` API name is retained as a compatibility spelling for its original subset-zeta behavior. See the [Operator Catalog](docs/operator_catalog.md) for exact semantics and expansion priorities. -Power-of-two lengths, batches, forward/inverse execution, normalization, +The C++ kernel-level API covers power-of-two lengths. The unified C plan API +also covers exact non-power-of-two FP32/FP64 FFT, compatible uint64 NTT, +explicit power-of-two embedding for every operator, and rank-two axis +composition. Batches, forward/inverse execution, normalization, in-place/out-of-place placement, padded batches, and strided elements are covered where listed in the [Feature Matrix](docs/cubutterfly_feature_coverage.md). Unsupported generated combinations fail explicitly rather than falling back to @@ -164,6 +188,23 @@ all precisions and shapes. | 60-bit NTT, `logN=16`, natural order | 0.274104 ms | GPU-NTT 0.291133 ms | 1.062x | fused Hybrid2D | | 60-bit NTT, `logN=20`, native bit-reversed | 0.341320 ms | GPU-NTT 0.396186 ms | 1.161x | compact-stage mapping | +### General-Shape Composition + +These v0.6 rows time the complete logical-input resident composition, not only +the selected power-of-two core: + +| Workload | Initial composition | Selected composition | Result | +|:--|--:|--:|:--| +| exact FP32 FFT `N=1000,batch=64` | 0.02622 ms Bluestein | 0.00427 ms direct cuFFT | 0.998x matching cuFFT; processing-unit abstention | +| exact rank-2 FP32 FFT `30x50,batch=64` | 0.07993 ms Bluestein | 0.00981 ms direct cuFFT | 1.000x matching cuFFT; processing-unit abstention | +| embedded FP32 FWHT `32767 -> 32768,batch=256` | 0.73772 ms | **0.11759 ms** | 6.27x composition gain; within 1.01x of pre-padded Dao core | +| embedded FP32 FWHT `255 -> 256,batch=65536` | 0.60309 ms | **0.16511 ms** | 3.65x composition gain | + +The FFT rows deliberately reuse the vendor exact-length unit and therefore do +not claim to outperform cuFFT. See [General-Shape V100 Selection +Results](docs/general_shape_results.md) for NTT, zeta, Structured, boundary, +workspace, and raw-trial details. + The generated local core, fused physical grouping, and vectorized two-pass composition bring the selected long FFT shapes to 95.4%-108.5% of cuFFT. Whole-transform cuFFTDx diff --git a/ROADMAP.md b/ROADMAP.md index e128749..199c74e 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -3,6 +3,25 @@ This roadmap separates implementation goals from evidence required for the research claim. Ordering may change after profiling or access to new hardware. +## v0.6.0 Status: General Shapes And Professional Plan Surface + +Status: implementation and single-V100 correctness/performance validation +complete; release packaging and review remain. +The new C plan API separates standard exact-length semantics from explicit +zero-extended embedding, lowers rank-two workloads to per-axis mappings, and +adds exact static-profile/cache selection with actionable miss logs. Exact FFT +selects between native cuFFT and Bluestein; modular Bluestein NTT is enabled +only when its root requirements hold and now reuses the selected radix-4 +Hybrid2D core. Embedded and rank-two plans select and cache each physical axis. +The initial Bluestein-only FFT reached 0.105x-0.179x of cuFFT; the new selector +reaches parity by abstaining to direct cuFFT. Boundary fusion for embedded +operators now removes contiguous output scatter, and selected FP32 +warp-register FWHT also fuses logical input and zero fill. The complete +logical `32767 -> 32768,batch=256` FWHT path reaches 0.11759 ms, within 1.01x +of the archived pre-padded Dao core. Structured fused input remains a generated +but unselected candidate because it regresses on V100. See +[v0.6 Implementation Status](docs/v0.6_implementation_status.md). + ## v0.4.0 Status: Library And Resource-Model Integration Status: complete on V100. This release adds the public device API and package, @@ -75,6 +94,11 @@ artifacts, and decision gates. ## Near Term Backlog +- Search boundary residence jointly with processing unit, numeric contract, + logical padding ratio, batch alignment, and layout instead of promoting a + fusion family from one successful operator. +- Add a padded-layout reuse contract for adjacent plan executions so an + application can retain physical data without repeated embedding work. - Extend the fixed-V100 regime model to stride, placement, direction, normalization, coefficient policy, and output-order neighborhoods using the same screen-follow-up-counter hierarchy. diff --git a/apps/cubutterfly_plan_bench.cpp b/apps/cubutterfly_plan_bench.cpp new file mode 100644 index 0000000..bf1b77d --- /dev/null +++ b/apps/cubutterfly_plan_bench.cpp @@ -0,0 +1,320 @@ +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { + +void check(cubutterflyStatus_t status, const char* operation) { + if (status != CUBUTTERFLY_STATUS_SUCCESS) + throw std::runtime_error(std::string(operation) + ": " + cubutterflyGetStatusString(status)); +} + +void check_cuda(cudaError_t status, const char* operation) { + if (status != cudaSuccess) + throw std::runtime_error(std::string(operation) + ": " + cudaGetErrorString(status)); +} + +void check_cufft(cufftResult status, const char* operation) { + if (status != CUFFT_SUCCESS) + throw std::runtime_error(std::string(operation) + ": cuFFT status " + std::to_string(status)); +} + +std::string take(int& index, int argc, char** argv) { + if (++index >= argc) + throw std::invalid_argument("missing command-line value"); + return argv[index]; +} + +std::vector parse_shape(const std::string& text) { + std::vector extents; + std::size_t begin = 0; + while (begin <= text.size()) { + const auto end = text.find('x', begin); + const auto token = text.substr(begin, end == std::string::npos ? end : end - begin); + std::size_t consumed = 0; + const auto extent = std::stoull(token, &consumed); + if (consumed != token.size() || extent == 0) + throw std::invalid_argument("shape extents must be positive integers"); + extents.push_back(extent); + if (end == std::string::npos) + break; + begin = end + 1; + } + if (extents.empty() || extents.size() > 2) + throw std::invalid_argument("shape must be N or RowsxColumns"); + return extents; +} + +std::array parse_matrix(const std::string& text) { + std::array values{}; + std::size_t begin = 0; + for (std::size_t index = 0; index < values.size(); ++index) { + const auto end = text.find(',', begin); + const auto token = text.substr(begin, end == std::string::npos ? end : end - begin); + std::size_t consumed = 0; + values[index] = std::stod(token, &consumed); + if (consumed != token.size() || (index + 1 != values.size()) != (end != std::string::npos)) + throw std::invalid_argument("matrix must contain exactly four comma-separated values"); + begin = end == std::string::npos ? text.size() : end + 1; + } + return values; +} + +void logger(cubutterflyLogLevel_t level, const char* event, const char* message, void*) { + const char* label = level == CUBUTTERFLY_LOG_ERROR ? "error" : + (level == CUBUTTERFLY_LOG_WARNING ? "warning" : "info"); + std::cerr << "cubutterfly[" << label << "][" << event << "]: " << message << '\n'; +} + +void usage() { + std::cout << "cubutterfly_plan_bench --operator fft|ntt|fwht|subset-zeta|superset-zeta|structured-2x2 " + "--shape N|RowsxColumns [--batch N] [--precision fp32|fp64|uint32|uint64] " + "[--length-mode standard|embedding] [--policy default|measure] [--algorithm ID] " + "[--matrix m00,m01,m10,m11] [--cache PATH] " + "[--inverse] [--in-place] [--compare-cufft] [--warmup N] [--repeat N]\n"; +} + +} // namespace + +int main(int argc, char** argv) { + cubutterflyOperator_t op = CUBUTTERFLY_OPERATOR_FFT; + cubutterflyDataType_t storage = CUBUTTERFLY_DATA_COMPLEX_FP32; + cubutterflyComputeType_t compute = CUBUTTERFLY_COMPUTE_FP32; + cubutterflyLengthMode_t length_mode = CUBUTTERFLY_LENGTH_STANDARD; + cubutterflyAlgorithmPolicy_t policy = CUBUTTERFLY_ALGORITHM_DEFAULT; + cubutterflyPlacement_t placement = CUBUTTERFLY_PLACEMENT_OUT_OF_PLACE; + std::vector shape{256}; + std::size_t batch = 1; + std::uint64_t modulus = 1152921504606584833ULL; + std::string cache; + int warmup = 10; + int repeat = 100; + bool inverse = false; + bool compare_cufft = false; + std::string explicit_algorithm; + std::array matrix{1.0, 0.25, -0.5, 1.0}; + try { + for (int index = 1; index < argc; ++index) { + const std::string arg = argv[index]; + if (arg == "--help" || arg == "-h") { + usage(); + return 0; + } else if (arg == "--shape") { + shape = parse_shape(take(index, argc, argv)); + } else if (arg == "--batch") { + batch = std::stoull(take(index, argc, argv)); + } else if (arg == "--warmup") { + warmup = std::stoi(take(index, argc, argv)); + } else if (arg == "--repeat") { + repeat = std::stoi(take(index, argc, argv)); + } else if (arg == "--cache") { + cache = take(index, argc, argv); + } else if (arg == "--modulus") { + modulus = std::stoull(take(index, argc, argv)); + } else if (arg == "--algorithm") { + explicit_algorithm = take(index, argc, argv); + policy = CUBUTTERFLY_ALGORITHM_EXPLICIT; + } else if (arg == "--matrix") { + matrix = parse_matrix(take(index, argc, argv)); + } else if (arg == "--inverse") { + inverse = true; + } else if (arg == "--in-place") { + placement = CUBUTTERFLY_PLACEMENT_IN_PLACE; + } else if (arg == "--compare-cufft") { + compare_cufft = true; + } else if (arg == "--operator") { + const auto value = take(index, argc, argv); + if (value == "fft") op = CUBUTTERFLY_OPERATOR_FFT; + else if (value == "ntt") op = CUBUTTERFLY_OPERATOR_NTT; + else if (value == "fwht") op = CUBUTTERFLY_OPERATOR_FWHT; + else if (value == "subset-zeta") op = CUBUTTERFLY_OPERATOR_SUBSET_ZETA; + else if (value == "superset-zeta") op = CUBUTTERFLY_OPERATOR_SUPERSET_ZETA; + else if (value == "structured-2x2") op = CUBUTTERFLY_OPERATOR_STRUCTURED_2X2; + else throw std::invalid_argument("unknown operator"); + } else if (arg == "--precision") { + const auto value = take(index, argc, argv); + if (value == "fp32") { storage = op == CUBUTTERFLY_OPERATOR_FFT ? CUBUTTERFLY_DATA_COMPLEX_FP32 : CUBUTTERFLY_DATA_FP32; compute = CUBUTTERFLY_COMPUTE_FP32; } + else if (value == "fp64") { storage = op == CUBUTTERFLY_OPERATOR_FFT ? CUBUTTERFLY_DATA_COMPLEX_FP64 : CUBUTTERFLY_DATA_FP64; compute = CUBUTTERFLY_COMPUTE_FP64; } + else if (value == "uint32") { storage = CUBUTTERFLY_DATA_UINT32; compute = CUBUTTERFLY_COMPUTE_UINT32; } + else if (value == "uint64") { storage = CUBUTTERFLY_DATA_UINT64; compute = CUBUTTERFLY_COMPUTE_UINT64; } + else throw std::invalid_argument("unknown precision"); + } else if (arg == "--length-mode") { + const auto value = take(index, argc, argv); + if (value == "standard") length_mode = CUBUTTERFLY_LENGTH_STANDARD; + else if (value == "embedding") length_mode = CUBUTTERFLY_LENGTH_ZERO_EXTENDED_EMBEDDING; + else throw std::invalid_argument("unknown length mode"); + } else if (arg == "--policy") { + const auto value = take(index, argc, argv); + if (value == "default") policy = CUBUTTERFLY_ALGORITHM_DEFAULT; + else if (value == "measure") policy = CUBUTTERFLY_ALGORITHM_MEASURE; + else throw std::invalid_argument("unknown policy"); + } else { + throw std::invalid_argument("unknown option: " + arg); + } + } + if (op == CUBUTTERFLY_OPERATOR_NTT && storage != CUBUTTERFLY_DATA_UINT32 && storage != CUBUTTERFLY_DATA_UINT64) { + storage = CUBUTTERFLY_DATA_UINT64; + compute = CUBUTTERFLY_COMPUTE_UINT64; + } + if ((op == CUBUTTERFLY_OPERATOR_SUBSET_ZETA || op == CUBUTTERFLY_OPERATOR_SUPERSET_ZETA)) { + storage = CUBUTTERFLY_DATA_UINT32; + compute = CUBUTTERFLY_COMPUTE_UINT32; + } + if (op == CUBUTTERFLY_OPERATOR_FFT) { + if (storage == CUBUTTERFLY_DATA_FP32) storage = CUBUTTERFLY_DATA_COMPLEX_FP32; + if (storage == CUBUTTERFLY_DATA_FP64) storage = CUBUTTERFLY_DATA_COMPLEX_FP64; + } else { + if (storage == CUBUTTERFLY_DATA_COMPLEX_FP32) storage = CUBUTTERFLY_DATA_FP32; + if (storage == CUBUTTERFLY_DATA_COMPLEX_FP64) storage = CUBUTTERFLY_DATA_FP64; + } + if (warmup < 0 || repeat <= 0) + throw std::invalid_argument("warmup must be nonnegative and repeat must be positive"); + if (compare_cufft && (op != CUBUTTERFLY_OPERATOR_FFT || inverse || + length_mode != CUBUTTERFLY_LENGTH_STANDARD)) + throw std::invalid_argument("cuFFT comparison currently requires a forward standard FFT"); + + cubutterflyHandle_t handle = nullptr; + cubutterflyDescriptor_t descriptor = nullptr; + cubutterflyPlan_t plan = nullptr; + cudaStream_t stream = nullptr; + check_cuda(cudaStreamCreateWithFlags(&stream, cudaStreamNonBlocking), "create stream"); + check(cubutterflyCreate(&handle), "create handle"); + check(cubutterflySetStream(handle, stream), "set stream"); + check(cubutterflySetLogCallback(handle, logger, nullptr), "set logger"); + if (!cache.empty()) + check(cubutterflySetCachePath(handle, cache.c_str()), "set cache path"); + check(cubutterflyCreateDescriptor(&descriptor), "create descriptor"); + check(cubutterflySetOperator(descriptor, op), "set operator"); + check(cubutterflySetDataType(descriptor, storage, compute), "set data type"); + check(cubutterflySetShape(descriptor, static_cast(shape.size()), shape.data(), batch), "set shape"); + check(cubutterflySetDirection(descriptor, inverse ? CUBUTTERFLY_DIRECTION_INVERSE : CUBUTTERFLY_DIRECTION_FORWARD, 1), "set direction"); + check(cubutterflySetPlacement(descriptor, placement), "set placement"); + check(cubutterflySetLengthMode(descriptor, length_mode), "set length mode"); + check(cubutterflySetAlgorithmPolicy(descriptor, policy, + explicit_algorithm.empty() ? nullptr : explicit_algorithm.c_str()), + "set algorithm policy"); + if (op == CUBUTTERFLY_OPERATOR_NTT) + check(cubutterflySetModulus(descriptor, modulus, storage == CUBUTTERFLY_DATA_UINT32 ? 32 : 64), "set modulus"); + if (op == CUBUTTERFLY_OPERATOR_STRUCTURED_2X2) { + const cubutterflyMatrix2x2_t stage{matrix[0], matrix[1], matrix[2], matrix[3]}; + for (std::uint32_t axis = 0; axis < shape.size(); ++axis) + check(cubutterflySetStageMatrices(descriptor, axis, &stage, 1), "set stage matrix"); + } + check(cubutterflyCreatePlan(handle, descriptor, &plan), "create plan"); + + std::size_t input_bytes = 0, output_bytes = 0, workspace_bytes = 0; + check(cubutterflyPlanGetInputSize(plan, &input_bytes), "query input size"); + check(cubutterflyPlanGetOutputSize(plan, &output_bytes), "query output size"); + check(cubutterflyPlanGetWorkspaceSize(plan, &workspace_bytes), "query workspace size"); + void* input = nullptr; + void* output = nullptr; + void* workspace = nullptr; + check_cuda(cudaMalloc(&input, std::max(input_bytes, output_bytes)), "allocate input"); + if (placement == CUBUTTERFLY_PLACEMENT_OUT_OF_PLACE) + check_cuda(cudaMalloc(&output, output_bytes), "allocate output"); + else + output = input; + if (workspace_bytes != 0) { + check_cuda(cudaMalloc(&workspace, workspace_bytes), "allocate workspace"); + check(cubutterflyPlanSetWorkspace(plan, workspace, workspace_bytes), "bind workspace"); + } + check_cuda(cudaMemsetAsync(input, 0, input_bytes, stream), "initialize input"); + for (int index = 0; index < warmup; ++index) + check(cubutterflyExecute(handle, plan, input, output), "warmup execute"); + cudaEvent_t start = nullptr, stop = nullptr; + check_cuda(cudaEventCreate(&start), "create start event"); + check_cuda(cudaEventCreate(&stop), "create stop event"); + check_cuda(cudaEventRecord(start, stream), "record start"); + for (int index = 0; index < repeat; ++index) + check(cubutterflyExecute(handle, plan, input, output), "timed execute"); + check_cuda(cudaEventRecord(stop, stream), "record stop"); + check_cuda(cudaEventSynchronize(stop), "wait for benchmark"); + float total_ms = 0.0F; + check_cuda(cudaEventElapsedTime(&total_ms, start, stop), "read benchmark time"); + float cufft_ms = 0.0F; + if (compare_cufft) { + if (batch > static_cast(std::numeric_limits::max())) + throw std::invalid_argument("cuFFT comparison batch exceeds int range"); + std::vector dimensions(shape.size()); + std::size_t distance = 1; + for (std::size_t axis = 0; axis < shape.size(); ++axis) { + if (shape[axis] > static_cast(std::numeric_limits::max())) + throw std::invalid_argument("cuFFT comparison extent exceeds int range"); + dimensions[axis] = static_cast(shape[axis]); + if (distance > static_cast(std::numeric_limits::max()) / shape[axis]) + throw std::invalid_argument("cuFFT comparison transform distance exceeds int range"); + distance *= shape[axis]; + } + cufftHandle baseline = 0; + check_cufft(cufftPlanMany(&baseline, static_cast(dimensions.size()), dimensions.data(), + nullptr, 1, static_cast(distance), nullptr, 1, + static_cast(distance), + storage == CUBUTTERFLY_DATA_COMPLEX_FP64 ? CUFFT_Z2Z : CUFFT_C2C, + static_cast(batch)), "create cuFFT comparison plan"); + check_cufft(cufftSetStream(baseline, stream), "set cuFFT comparison stream"); + auto execute_cufft = [&] { + if (storage == CUBUTTERFLY_DATA_COMPLEX_FP64) + check_cufft(cufftExecZ2Z(baseline, static_cast(input), + static_cast(output), CUFFT_FORWARD), + "execute cuFFT comparison"); + else + check_cufft(cufftExecC2C(baseline, static_cast(input), + static_cast(output), CUFFT_FORWARD), + "execute cuFFT comparison"); + }; + for (int index = 0; index < warmup; ++index) + execute_cufft(); + check_cuda(cudaEventRecord(start, stream), "record cuFFT start"); + for (int index = 0; index < repeat; ++index) + execute_cufft(); + check_cuda(cudaEventRecord(stop, stream), "record cuFFT stop"); + check_cuda(cudaEventSynchronize(stop), "wait for cuFFT benchmark"); + check_cuda(cudaEventElapsedTime(&cufft_ms, start, stop), "read cuFFT benchmark time"); + cufft_ms /= repeat; + cufftDestroy(baseline); + } + std::size_t name_bytes = 0; + check(cubutterflyPlanGetAlgorithmName(plan, nullptr, &name_bytes), "query algorithm name size"); + std::vector algorithm(name_bytes); + check(cubutterflyPlanGetAlgorithmName(plan, algorithm.data(), &name_bytes), "query algorithm name"); + std::vector physical(shape.size()); + check(cubutterflyPlanGetPhysicalExtents(plan, physical.data(), static_cast(physical.size())), "query physical shape"); + std::cout << "algorithm=" << algorithm.data() << ",logical_shape="; + for (std::size_t axis = 0; axis < shape.size(); ++axis) + std::cout << (axis ? "x" : "") << shape[axis]; + std::cout << ",physical_shape="; + for (std::size_t axis = 0; axis < physical.size(); ++axis) + std::cout << (axis ? "x" : "") << physical[axis]; + std::cout << ",batch=" << batch << ",workspace_bytes=" << workspace_bytes + << ",kernel_ms=" << total_ms / repeat; + if (compare_cufft) + std::cout << ",cufft_ms=" << cufft_ms << ",speedup_vs_cufft=" + << cufft_ms / (total_ms / repeat); + std::cout << '\n'; + + cudaEventDestroy(stop); + cudaEventDestroy(start); + cudaFree(workspace); + if (placement == CUBUTTERFLY_PLACEMENT_OUT_OF_PLACE) + cudaFree(output); + cudaFree(input); + cubutterflyDestroyPlan(plan); + cubutterflyDestroyDescriptor(descriptor); + cubutterflyDestroy(handle); + cudaStreamDestroy(stream); + } catch (const std::exception& error) { + std::cerr << "error: " << error.what() << '\n'; + return 1; + } + return 0; +} diff --git a/config/v100_application_profiles.json b/config/v100_application_profiles.json new file mode 100644 index 0000000..d8f28f1 --- /dev/null +++ b/config/v100_application_profiles.json @@ -0,0 +1,20 @@ +{ + "schema_version": 1, + "target": { + "name": "v100-sm70-80sm", + "compute_major": 7, + "compute_minor": 0, + "multiprocessors": 80 + }, + "profiles": [ + {"operator": "fft", "storage": "complex-fp32", "direction": "forward", "placement": "in-place", "logN": 8, "batches": [1, 4, 16, 64, 256, 1024, 4096, 16384, 65536], "implementation": "cuButterfly-cuFFTDx"}, + {"operator": "fft", "storage": "complex-fp32", "direction": "forward", "placement": "in-place", "logN": 14, "batches": [1, 4, 16, 64, 256, 1024], "implementation": "cuButterfly-cuFFTDx-direct"}, + {"operator": "fft", "storage": "complex-fp32", "direction": "forward", "placement": "in-place", "logN": 18, "batches": [1, 2, 4, 8, 16, 32, 64], "implementation": "cuButterfly-cuFFTDx-online"}, + {"operator": "fft", "storage": "complex-fp32", "direction": "forward", "placement": "in-place", "logN": 20, "batches": [1, 2, 4, 8, 16], "implementation": "cuButterfly-cuFFTDx-online"}, + {"operator": "fwht", "storage": "fp32", "direction": "forward", "placement": "out-of-place", "logN": 8, "batches": [1, 4, 16, 64, 256, 1024, 4096, 16384, 65536], "implementation": "cuButterfly-warp-register"}, + {"operator": "fwht", "storage": "fp32", "direction": "forward", "placement": "out-of-place", "logN": 15, "batches": [1, 4, 16, 32, 128, 256, 512], "implementation": "cuButterfly-piecewise"}, + {"operator": "fwht", "storage": "fp32", "direction": "forward", "placement": "out-of-place", "logN": 20, "batches": [1, 2, 4, 8, 16], "implementation": "cuButterfly-online-radix4"}, + {"operator": "ntt", "storage": "uint32", "direction": "forward", "placement": "out-of-place", "modulus": 1073692673, "logN": 12, "batches": [1, 4, 16, 64, 256, 1024, 4096], "implementation": "cuNTT-Hybrid2D-radix4"}, + {"operator": "ntt", "storage": "uint64", "direction": "forward", "placement": "out-of-place", "modulus": 1152921504606584833, "logN": 16, "batches": [1, 4, 16, 32, 64, 128, 256], "implementation": "cuNTT-Hybrid2D-radix4"} + ] +} diff --git a/docs/README.md b/docs/README.md index c03956b..efc129a 100644 --- a/docs/README.md +++ b/docs/README.md @@ -12,8 +12,11 @@ Started](getting_started.md); architecture and paper readers should begin with | Document | Question answered | |:--|:--| | [Getting Started](getting_started.md) | How is the library configured, built, installed, tested, and called? | +| [C Plan API](c_api.md) | How do handle, descriptor, plan, workspace, arbitrary shape, tuning cache, and logs compose? | +| [General Shape Mapping](general_shapes.md) | How are exact, embedded, and rank-two workloads lowered to the same architecture mapping? | +| [v0.6 Implementation Status](v0.6_implementation_status.md) | Which general-shape, API, selection, and boundary features are complete, and what remains outside the claim? | | [Programming Guide](programming_guide.md) | What are the plan, device, stream, layout, workspace, concurrency, and lifetime contracts? | -| [C++ API Reference](api_reference.md) | What does each public type, configuration field, and plan member mean? | +| [API Reference](api_reference.md) | What does each public type, configuration field, and plan member mean? | | [Error Handling](error_handling.md) | Which exceptions are reported synchronously and where can asynchronous CUDA errors surface? | | [Capability and Compatibility Matrix](support_matrix.md) | Which operators, types, toolchains, optional components, and runtime features are supported? | | [Operator Catalog](operator_catalog.md) | What does each transform compute, which mappings does it reuse, and which operators come next? | @@ -62,6 +65,8 @@ comprehensive and scaling reports define the paper-facing comparison protocol. | [Candidate Implementations](implementation_candidates.md) | Which processing units and kernel forms are implemented or planned? | | [Getting Started](getting_started.md) | How is the project built and used? | | [Programming Guide](programming_guide.md) | How do plans, layouts, streams, workspaces, devices, concurrency, and determinism compose? | +| [C Plan API](c_api.md) | What is the professional C integration and selection contract? | +| [General Shape Mapping](general_shapes.md) | How do non-power-of-two and rank-two lowerings preserve operator semantics? | | [C++ API Reference](api_reference.md) | What is the contract of each installed public API type and member? | | [Error Handling](error_handling.md) | How are validation, backend, CUDA launch, and asynchronous failures reported? | | [Capability and Compatibility Matrix](support_matrix.md) | Which public semantics, optional builds, and runtime integrations are supported? | @@ -79,6 +84,8 @@ comprehensive and scaling reports define the paper-facing comparison protocol. | [V100 External Baselines](v100_external_baselines.md) | What do matching-protocol Dao FHT and GPU-NTT comparisons establish? | | [V100 Library/Base/Search Comparison](v100_three_way_comparison.md) | How much comes from search, and where does the selected point stand against specialized libraries? | | [V100 Research Status](research_status.md) | Which claims are closed, what are the current performance boundaries, and what remains? | +| [General-Shape V100 Selection Results](general_shape_results.md) | How much do composition and physical-core selection recover, and which boundaries remain? | +| [v0.6 Implementation Status](v0.6_implementation_status.md) | What is implemented, validated, selected by default, retained only as a candidate, or deferred? | | [v0.3.0 Mapping-Selection Milestone](next_phase_v0.3.md) | Which falsifiable goal, execution order, and definition of done shaped this release? | | [Numeric-Regime Mapping Study](next_phase_numeric_regimes.md) | How are precision, arithmetic, length, batch, and core choice related to resource cliffs, and what did the first quick screen establish? | | [Reproducibility](reproducibility.md) | How are sweeps, external baselines, and profiler data reproduced? | diff --git a/docs/api_reference.md b/docs/api_reference.md index 82fc56f..ecbd71a 100644 --- a/docs/api_reference.md +++ b/docs/api_reference.md @@ -1,4 +1,12 @@ -# C++ API Reference +# API Reference + +The library now also installs a stable-style C plan surface in +``. It follows the handle, descriptor, immutable +plan, workspace-query, and asynchronous execute sequence used by CUDA +professional libraries. See [C Plan API](c_api.md) for the complete contract, +general shapes, selection cache, and logging behavior. + +## C++ API The installed public headers are `` and ``. All symbols are in namespace `cuntt`. The library requires @@ -137,6 +145,7 @@ The associated enums expose: | `set_workspace(ptr, bytes)` | binds aligned caller-owned scratch or removes a binding with null | | `workspace()` | active external or plan-owned workspace pointer | | `execute_async(input, output)` | typed, allocation-free device submission without synchronization | +| `execute_zero_extended_async(input, output, logical_points, input_batch_stride, input_element_stride)` | advanced FP32 warp-register boundary submission; reads a logical FWHT/Structured input and writes the full physical transform | | `execute(input, output, warmup, repeat)` | synchronous host-vector convenience and timing path | The nine `execute_async` overloads accept `Fp16`, `Bf16`, `float`, `double`, @@ -148,6 +157,12 @@ explicit FP32-compute-and-narrow emulation and benchmark CSV marks it with `ButterflyStats` reports H2D, kernel and D2H milliseconds plus transforms, butterflies, and points per second. +`execute_zero_extended_async` is a processing-unit interface used by the C +plan composition layer. It requires a forward FP32 warp-register FWHT or +Structured 2x2 plan, distinct device pointers, and a logical extent no larger +than `2^log_n`. Most applications should request embedding through the C plan +API so selection can keep the explicit pack path when fusion is slower. + CPU helpers `reference_fwht`, `reference_fft`, `reference_subset_zeta`, `reference_superset_zeta`, `reference_structured_2x2`, and legacy `reference_xor_zeta` provide small-scale diff --git a/docs/c_api.md b/docs/c_api.md new file mode 100644 index 0000000..5d41933 --- /dev/null +++ b/docs/c_api.md @@ -0,0 +1,115 @@ +# C Plan API + +The installed `` header provides a CUDA-library-style +C ABI. Applications create one handle per stream context, describe an operator +and workload, create an immutable plan, query and bind device workspace, then +submit allocation-free executions. + +```c +cubutterflyHandle_t handle; +cubutterflyDescriptor_t descriptor; +cubutterflyPlan_t plan; +size_t shape[2] = {3, 5}; +size_t workspace_bytes; + +cubutterflyCreate(&handle); +cubutterflySetStream(handle, stream); +cubutterflyCreateDescriptor(&descriptor); +cubutterflySetOperator(descriptor, CUBUTTERFLY_OPERATOR_FFT); +cubutterflySetDataType(descriptor, CUBUTTERFLY_DATA_COMPLEX_FP32, + CUBUTTERFLY_COMPUTE_FP32); +cubutterflySetShape(descriptor, 2, shape, batch); +cubutterflyCreatePlan(handle, descriptor, &plan); +cubutterflyPlanGetWorkspaceSize(plan, &workspace_bytes); +cudaMalloc(&workspace, workspace_bytes); +cubutterflyPlanSetWorkspace(plan, workspace, workspace_bytes); +cubutterflyExecute(handle, plan, input, output); +``` + +Plan creation may initialize coefficient tables and vendor plans. Execution +uses the handle's current stream, performs no allocation, and does not +synchronize the host. One plan must not be submitted concurrently from +multiple host threads. Independent plans and workspaces may use independent +streams. + +## Length Semantics + +`CUBUTTERFLY_LENGTH_STANDARD` preserves the requested mathematical length. +Power-of-two axes use native butterfly plans. For a packed non-power-of-two +FP32/FP64 FFT, direct cuFFT and Bluestein with `M = next_pow2(2N-1)` are legal +processing-unit candidates; the default currently prefers direct cuFFT and +`MEASURE` times both complete compositions. Strided data uses Bluestein. A +uint64 NTT uses modular Bluestein only when the prime modulus admits both the +required `2N`-th root and the power-of-two convolution root. Plan creation +rejects an invalid root domain. + +`CUBUTTERFLY_LENGTH_ZERO_EXTENDED_EMBEDDING` is deliberately different: each +axis is padded to `next_pow2(N)` and the power-of-two transform is returned. +Forward input has the logical shape and forward output has the embedded shape; +inverse reverses those extents. The physical extents and byte requirements are +queryable, so callers never infer allocation sizes. + +Rank two is lowered as row transform, tiled transpose, column transform, and +transpose back. This applies to FFT, NTT, FWHT, subset/superset zeta, and +Structured 2x2. The transpose is a mapping boundary, not a change to operator +arithmetic. + +Contiguous physical outputs use a direct-output boundary: the final physical +core or transpose lands in the caller buffer and no scatter kernel is issued. +Strided outputs retain scatter semantics. `cubutterflyPlanGetAlgorithmName` +marks the former with `+direct-output`. +Workspace queries reflect the lowering: a rank-one direct-output butterfly may +need no composition buffer, a direct-output NTT keeps one packed input buffer, +and rank two keeps both transpose buffers. +Forward FP32 FWHT embedding can additionally fuse logical reads and zero fill +into a selected warp-register unit. Its algorithm name includes +`+fused-input`, it launches one kernel, and its workspace remains zero. +The suffix is deliberately absent when the selected core, operator, direction, +or layout uses the materialized pack fallback. + +## Selection And Cache + +The default policy uses this order: + +1. exact workload key in the caller-selected cache; +2. exact device/semantics/length/batch entry in the generated application table; +3. legal static resource-model fallback. + +`CUBUTTERFLY_ALGORITHM_MEASURE` times a bounded legal candidate set. Direct +power-of-two plans tune one plan, embedded/rank-two plans tune every physical +axis, modular Bluestein tunes its power-of-two NTT core, and exact arbitrary +FFT compares the complete direct and Bluestein paths. When a cache path is +configured, winners are appended and subsequent default plans reuse them. +Cache keys include compute +capability, SM count, operator, numeric contract, direction, placement, length +mode, rank, extents, batch, strides, and NTT modulus. Moving the file to a +different GPU therefore cannot cause an accidental match. + +Register a log callback with `cubutterflySetLogCallback`. A missing device +profile emits `hardware-profile-miss`; a known device with an unseen workload +emits `workload-profile-miss`. Both messages include the selected fallback and +recommend measured calibration. `cubutterflyPlanGetSelectionSource`, +`cubutterflyPlanGetAlgorithmName`, and `cubutterflyPlanGetSelectionReason` +provide the same decision programmatically. + +The generated V100 application table is sourced from +`config/v100_application_profiles.json`. It contains only archived exact +anchors and does not label interpolation or a static fallback as measured. +`cubutterfly_plan_bench --compare-cufft` records a matching forward standard +FFT baseline for rank-one or rank-two shapes and reports the resident-kernel +ratio; plan creation and host transfers remain outside both timings. + +## Current Boundaries + +- Standard arbitrary-length FFT: complex FP32/FP64, rank one or two. +- Standard arbitrary-length NTT: uint64, rank one or two, compatible prime/root domain. +- Embedded non-power-of-two shapes: all public operators and rank one or two. +- Runtime `MEASURE`: direct plans, per-axis composite plans, and exact FFT + direct-versus-Bluestein selection. Composite tuning currently minimizes each + axis independently because pack/transpose boundaries are common to its + candidates. +- Structured 2x2 requires matrices for every executed axis. + +All API calls return `cubutterflyStatus_t`. Invalid descriptors, unsupported +mathematical domains, allocation failures, and backend failures are returned +across the C boundary; C++ exceptions never escape. diff --git a/docs/cubutterfly_feature_coverage.md b/docs/cubutterfly_feature_coverage.md index 209df9e..0db6c5c 100644 --- a/docs/cubutterfly_feature_coverage.md +++ b/docs/cubutterfly_feature_coverage.md @@ -22,6 +22,10 @@ arithmetic contract only on hardware that implements the required BF16 path. ## Length and Backend Capabilities +The table below describes power-of-two physical plans exposed by the C++ API. +The general-shape C plan composes those plans with the exact/embedded and +rank-two lowerings summarized after the table. + | Backend | Length | Processing unit | FP32 | FP64 | uint32 | |:--|:--|:--|:--:|:--:|:--:| | temporal-tile | shared: `logN=1..10`; register FWHT/Structured: `logN=3..15` | radix-2, fused radix-4/radix-8; FP32 FWHT and Structured register-vector units | yes | yes | yes | @@ -110,6 +114,21 @@ different lengths, precisions, directions, normalization policies, structured processing units, and batches cannot be accidentally combined. +## General-Shape C Plan Coverage + +| Shape contract | Operators and restrictions | +|:--|:--| +| native power of two | all public operators; rank one or two | +| standard exact arbitrary length | packed complex FP32/FP64 FFT through selected direct cuFFT/Bluestein; uint64 NTT through modular Bluestein when roots exist | +| zero-extended embedding | all public operators; rank one or two; forward returns the physical shape and inverse consumes it | +| rank-two layout | per-axis physical plans with tiled transpose; contiguous final output can receive the final transpose directly | +| fused input | selected forward FP32 warp-register FWHT; generated Structured candidate is correct but not selected on V100 | + +General-shape plans accept positive input/output element and batch strides. +Direct packed cuFFT requires its documented contiguous layout; strided exact +FFT uses Bluestein. Algorithm names and workspace queries expose whether pack, +scatter, transpose, direct output, or fused input is present. + ## Verification The standalone test executables cover FP32/uint32 temporal radix-2, radix-4, @@ -147,7 +166,7 @@ multiplication; these are separate CSV fields. The current public runtime does not yet claim: - full-length FP16/BF16 Tensor Core FFT beyond the generated DFT8 point; -- non-power-of-two or multidimensional transforms; +- rank greater than two or mixed-radix native butterfly cores; - `logN>20` through the common runtime; - fused application epilogues such as convolution or quantization. diff --git a/docs/experiments.md b/docs/experiments.md index d1bfd85..1999b92 100644 --- a/docs/experiments.md +++ b/docs/experiments.md @@ -56,6 +56,17 @@ experiment tests a stage-parameterized dense processing unit across the same mapping families. Its follow-up generated register codelet separates mapping and transport gains from the remaining dense-pair arithmetic cost. +The v0.6 [General-Shape V100 Selection Results](general_shape_results.md) keep +the initial Bluestein-only diagnostic separate from the five-trial selector +follow-up. The first lowering reached only 0.105x-0.179x of direct cuFFT; the +new composition selector reaches FFT parity through direct-core abstention and +substantially improves NTT/FWHT/Structured physical-core choices. Embedded +boundaries remain explicit and must not be merged with bare-core library +comparisons. The follow-up removes contiguous output scatter and fuses input +for selected warp-register FWHT: `32767 -> 32768,batch=256` reaches 0.11759 ms, +within 1.01x of the archived pre-padded Dao core. Structured input fusion is a +measured negative candidate and retains the pack path. + The evidence chain is therefore: ```text diff --git a/docs/general_shape_results.md b/docs/general_shape_results.md new file mode 100644 index 0000000..d972f1b --- /dev/null +++ b/docs/general_shape_results.md @@ -0,0 +1,97 @@ +# General-Shape V100 Selection Results + +## Composition Selector Follow-up + +The second experiment adds direct cuFFT as an exact FFT processing-unit +candidate, propagates physical-axis selection into embedded/rank-two plans, +and lets modular Bluestein use the selected NTT core. Each row below is the +median of five independent warmed processes with 50 warmups and 100 timed +executions. Plan creation, tuning, cache I/O, and host transfers are excluded. + +| Operator | Logical shape | Batch | Selected composition | Initial ms | Selected-core ms | Direct-output ms | Fused-input ms | Total gain | +|:--|:--|--:|:--|--:|--:|--:|--:|--:| +| FFT FP32 | 1000 | 64 | direct cuFFT | 0.02622 | 0.00427 | 0.00427 | 0.00427 | 6.14x | +| FFT FP32 | 30x50 | 64 | direct cuFFT | 0.07993 | 0.00981 | 0.00981 | 0.00981 | 8.15x | +| NTT uint64 | 5419 | 1 | Bluestein + Hybrid2D radix-4 | 0.08608 | 0.03147 | 0.03147 | 0.03147 | 2.74x | +| NTT uint64 | 5419 | 64 | Bluestein + Hybrid2D radix-4 | 0.79822 | 0.21323 | 0.21323 | 0.21323 | 3.74x | +| FWHT FP32 embedding | 32767 -> 32768 | 256 | warp-register radix-2, fused load | 0.73772 | 0.32409 | 0.24271 | **0.11759** | **6.27x** | +| FWHT FP32 embedding | 255 -> 256 | 65536 | warp-register radix-2, fused load | 0.60309 | 0.58422 | 0.41231 | **0.16511** | **3.65x** | +| subset zeta embedding | 32767 -> 32768 | 256 | hierarchical radix-4 | 0.73606 | 0.73606 | 0.64882 | 0.64882 | 1.13x | +| Structured FP32 embedding | 32767 -> 32768 | 256 | measured warp-register radix-2 | 0.74460 | 0.34425 | 0.25772 | not selected | 2.89x | + +Direct FFT is at timing parity with the matching cuFFT invocation: 0.998x at +`N=1000,batch=64` and 1.000x at `30x50,batch=64`. This is intentional +processing-unit reuse, not a claim that a wrapper around cuFFT outperforms +cuFFT. Runtime `MEASURE` still compares the complete direct and Bluestein +compositions and records the winner, while explicit Bluestein remains +available for architecture experiments. + +There is no matching exact arbitrary-length GPU-NTT baseline in the archived +external suite. Relative to this repository's direct power-of-two core at the +same physical `M=16384`, modular Bluestein costs 2.97x at batch 1 and 2.72x at +batch 64. That is now close to its two transforms plus clear/pack/multiply/output +boundaries; the previous excess came from hard-coding baseline radix-2 instead +of reusing Hybrid2D radix-4. + +FWHT, subset/superset zeta, and Structured 2x2 use zero-extended embedding, +not an exact `N`-point transform. For a contiguous physical output, the second +implementation removes the final scatter: rank-one floating/integer plans +write their last core directly to the caller output, while rank-two plans make +the final transpose land there. Strided outputs retain the materialized +scatter path. Algorithm names expose this decision as `+direct-output`. + +For FP32 warp-register FWHT, the generated first processing unit now accepts +the logical length and independent input/output strides. It predicates the +tail to zero while loading the useful values, so no pack kernel or composition +buffer is required. The algorithm name exposes this as `+fused-input`. At +`P=32768,batch=256`, the complete logical-shape wrapper is 0.11759 ms, within +1.01x of the archived 0.11643 ms bare Dao core despite accepting `N=32767` +input batches rather than pre-padded input. At `P=256,batch=65536`, removing +the pack boundary gives a further 2.50x over the direct-output version. + +This is a selected lowering, not a universal rule. Misaligned logical batch +distances force scalar loads in some batches. A two-aligned-vector stitching +candidate was measured and rejected because redundant memory transactions +increased latency. Structured 2x2 also retains the pack path: its tested fused +candidate was slower than the direct-output composition. The generator keeps +the specialized entry as a candidate for future aligned/padded input contracts, +but the default C plan only emits `+fused-input` for the measured-beneficial +FWHT family. + +Subset and superset zeta remain symmetric at 0.64882/0.64900 ms and have no +semantics-matched external GPU baseline. + +The repeated raw rows are +[`results/general_shape_selection_v100_raw.csv`](../results/general_shape_selection_v100_raw.csv) +and +[`results/general_shape_direct_boundary_v100_raw.csv`](../results/general_shape_direct_boundary_v100_raw.csv). +The fused-boundary trials are +[`results/general_shape_fused_input_v100_raw.csv`](../results/general_shape_fused_input_v100_raw.csv). + +## Initial Bluestein-only Diagnostic + +The first v0.6 diagnostic forced exact arbitrary FFT through Bluestein. It used +one aggregated CUDA-event run per row with 20 warmups and 100 timed executions. + +| Shape | Batch | Bluestein physical shape | cuButterfly ms | cuFFT ms | cuButterfly/cuFFT | +|:--|--:|:--|--:|--:|--:| +| 5 | 1 | 16 | 0.01534 | 0.00274 | 0.179x | +| 5 | 4096 | 16 | 0.01878 | 0.00271 | 0.144x | +| 1000 | 1 | 2048 | 0.02048 | 0.00337 | 0.165x | +| 1000 | 64 | 2048 | 0.02622 | 0.00428 | 0.163x | +| 3x5 | 4096 | 8x16 | 0.06553 | 0.00690 | 0.105x | +| 30x50 | 1 | 64x128 | 0.04690 | 0.00788 | 0.168x | +| 30x50 | 64 | 64x128 | 0.07993 | 0.00979 | 0.122x | + +Those rows are in +[`results/general_shapes_v100_smoke.csv`](../results/general_shapes_v100_smoke.csv). +Rank one performed chirp/pack, forward power-of-two cuFFT, pointwise multiply, +inverse cuFFT, and chirp/output. Rank two repeated that sequence per axis and +added transpose boundaries. The result established that a fast physical core +does not by itself make an efficient composition. + +The follow-up implements processing-unit abstention, physical-core selection, +direct physical output, and a selected fused input boundary for warp-register +FWHT. Remaining work is to make boundary residence a searched property for +more processing units and to retain padded layouts across adjacent library +calls so applications can avoid repeated embedding work. diff --git a/docs/general_shapes.md b/docs/general_shapes.md new file mode 100644 index 0000000..a9c145e --- /dev/null +++ b/docs/general_shapes.md @@ -0,0 +1,76 @@ +# General Shape Mapping + +General shapes do not add another architecture principle. They add a lowering +step before the same two-dimensional space-time mapping is selected. + +```text +logical axes and batch + | + +-- power-of-two axis --------> native staged butterfly + | + +-- standard arbitrary FFT ---> direct vendor plan OR chirp convolution + | + +-- standard arbitrary NTT ---> root validation -> modular chirp/convolution + | + +-- embedding ----------------> zero pad to next power of two + v +axis plan 0 -> online transpose -> axis plan 1 -> transpose back + v +(Ud,Td) data mapping x (Us,Ts) stage mapping x selected processing unit +``` + +For an exact FFT axis of length `N`, direct cuFFT is a processing-unit candidate. +Bluestein converts the transform into a +convolution of length at least `2N-1`; the implementation chooses the next +power of two `M` so an existing high-performance FFT core remains a selectable +physical unit. The chirp pack, pointwise multiply, and chirp output are fixed +boundary work. They do not alter the FFT formula. + +An arbitrary NTT uses the same composition over the finite field. Unlike a +floating FFT, its domain is conditional: the selected prime must contain the +needed roots of unity. Rejecting a missing root is a semantic requirement, not +a performance limitation. Embedding remains available when the application +actually wants a padded power-of-two transform rather than an exact-length NTT. + +For rank two, rows provide independent data-space work for the first axis. +After the first axis completes, transpose changes ownership/layout so columns +become packed independent transforms for the second axis. Each axis retains +its own `(Us,Ts,Ud,Td)` and processing-unit choice. The current implementation +materializes this global boundary in caller workspace because ordinary CTAs +cannot perform a grid-wide synchronization inside a conventional kernel. + +The physical extent query has operator-specific meaning: + +| Mode | Returned physical axis extent | +|:--|:--| +| native power of two | `N` | +| standard arbitrary FFT | logical `N` for direct cuFFT; convolution length for Bluestein | +| standard arbitrary NTT | convolution length `next_pow2(2N-1)` | +| zero-extended embedding | embedded transform length `next_pow2(N)` | + +Input/output byte queries always describe the logical API contract, including +the asymmetric forward/inverse embedding contract. Workspace size includes +packed matrices, transpose buffers, convolution buffers, and underlying axis +plan scratch. + +When the physical output layout is contiguous, the composition does not +materialize a separate final scatter. A rank-one embedded butterfly packs and +zero-fills directly into the caller output before its in-place physical core; +an out-of-place NTT core writes there directly; and a rank-two plan makes its +final transpose target the caller output. Arbitrary output strides keep the +explicit scatter. The selected algorithm name carries `+direct-output` so the +boundary choice is visible in logs and result files. +For a rank-one in-place physical butterfly this also removes the packed +workspace buffer entirely; the workspace query then reports only underlying +core scratch and may return zero. A rank-one out-of-place NTT retains one pack +buffer, and rank two retains the two transpose buffers. + +For forward FP32 FWHT embedding with a generated warp-register physical unit, +the input boundary is lowered into that unit as well. The kernel receives the +logical extent plus the independent input and output layouts, reads valid +logical values, synthesizes zero for the tail in registers, and writes the +complete physical result. This eliminates both composition kernels and is +reported as `+direct-output+fused-input`. Other processing units retain the +pack path unless the selector has a validated fused candidate. In particular, +Structured 2x2 currently rejects this candidate on performance grounds; a +legal fusion is not automatically an efficient mapping. diff --git a/docs/getting_started.md b/docs/getting_started.md index cef64e1..4afbcca 100644 --- a/docs/getting_started.md +++ b/docs/getting_started.md @@ -117,6 +117,35 @@ end-to-end rates. ./build/cubutterfly_bench --list-capabilities ``` +For the unified C plan surface, arbitrary lengths, and rank-two shapes, use: + +```bash +# Exact 3x5 FP32 FFT; default selects direct cuFFT when it wins +./build/cubutterfly_plan_bench --operator fft --shape 3x5 --repeat 100 +./build/cubutterfly_plan_bench --operator fft --shape 3x5 --compare-cufft --repeat 100 + +# Retain Bluestein as an explicit research composition +./build/cubutterfly_plan_bench --operator fft --shape 3x5 \ + --algorithm bluestein-cufft-power2-core --repeat 100 + +# Exact length-3 uint64 NTT; the modulus/root domain is validated +./build/cubutterfly_plan_bench --operator ntt --shape 3 --precision uint64 + +# One-time search followed by exact-key cache reuse +./build/cubutterfly_plan_bench --operator fwht --shape 4096 --batch 64 \ + --policy measure --cache results/local_selection.cache +./build/cubutterfly_plan_bench --operator fwht --shape 4096 --batch 64 \ + --cache results/local_selection.cache + +# Embedded Structured 2x2 with one broadcast matrix per axis +./build/cubutterfly_plan_bench --operator structured-2x2 --shape 32767 \ + --length-mode embedding --matrix 1,0.25,-0.5,1 --policy measure \ + --cache results/local_selection.cache +``` + +Profile misses are printed with the hardware/workload key and selected +fallback. See [C Plan API](c_api.md) and [General Shape Mapping](general_shapes.md). + Examples: ```bash diff --git a/docs/programming_guide.md b/docs/programming_guide.md index 0a6da4a..a5d8281 100644 --- a/docs/programming_guide.md +++ b/docs/programming_guide.md @@ -5,6 +5,12 @@ This guide defines the application-facing execution model. It complements the functions, and the [Capability and Compatibility Matrix](support_matrix.md), which states the supported semantic ranges. +The plan types described below are the power-of-two C++ physical-plan APIs. +Applications requiring exact arbitrary lengths, zero-extended embedding, +rank-two composition, selection caches, or status-code error handling should +use the [C Plan API](c_api.md); both surfaces share the same stream and caller- +owned workspace principles. + ## Execution Model Both API families use an immutable-shape plan: diff --git a/docs/research_status.md b/docs/research_status.md index ca22e27..8bd9df0 100644 --- a/docs/research_status.md +++ b/docs/research_status.md @@ -6,6 +6,10 @@ cuButterfly is a single-GPU research artifact for hardware-mapped space-time parallelism across regular butterfly computations. V100 is the only fully measured hardware target. Cross-GPU transfer is explicitly future work and is not used to support any current portability claim. +Power-of-two and numeric-regime evidence below is the v0.5 authority. The +general-shape/API extension and its separate protocols are summarized in +[v0.6 Implementation Status](v0.6_implementation_status.md); those results do +not retroactively change older controlled tables. ## Research Claims And Evidence @@ -19,6 +23,7 @@ not used to support any current portability claim. | A small hardware calibration set can replace exhaustive timing | leave-one-batch-out FFT interpolation | top-3 geomean regret 1.0007x, worst 1.0020x using 6.2% of candidates | | The selection method extends beyond FFT | leave-one-complete-shape-out selector over NTT/FWHT/XOR-zeta | 93.88% top-1, 100% top-3, 1.0049x geomean regret | | Numeric and batch preference is piecewise rather than globally smooth | 736-case screen, 29 full follow-ups, 12 paired NCU profiles, 99 adaptive anchors | all 9 length crossovers reproduce; 6 batch events are unstable; bounded selector coverage is 27.74% with 1.01112x worst regret | +| General-shape composition must select both physical units and boundaries | exact FFT/NTT, embedding, rank-two, direct-output and fused-input follow-ups | direct exact FFT reaches vendor parity by abstention; embedded `32767 -> 32768` FWHT reaches 0.11759 ms after selected boundary fusion | The static-model failure is part of the result, not a row to hide. Resource legality, occupancy, waves, and issued work provide a useful shortlist, but do @@ -144,9 +149,11 @@ boundaries. established, capture a newer GPU profile and test whether the learned descriptors and predicted boundaries transfer before and after calibration. -The repository does not currently claim a production cuFFT API replacement, -non-power-of-two coverage, universal library superiority, or architecture -portability beyond V100. +The repository does not currently claim a production cuFFT replacement, +universal library superiority, or architecture portability beyond V100. +Non-power-of-two and rank-two execution are implemented through the v0.6 C +plan composition layer, but only the documented exact/embedding semantics and +single-V100 evidence are claimed. ## Reproduction Map @@ -164,6 +171,8 @@ portability beyond V100. - FP64 scalar/core/mapping validation: `results/fp64_*logN16*` - FP64 length/batch robustness: `results/fp64_robustness_*` - FP64 privileged counter command: `scripts/profile_fp64_fft_ncu.sh` +- general-shape composition: `results/general_shape_*`, + `docs/general_shape_results.md` Run `./scripts/reproduce_v100_analysis.sh` to regenerate derived artifacts from the checked-in raw measurements. It does not recollect timings or counters. diff --git a/docs/support_matrix.md b/docs/support_matrix.md index 1d2fc1e..a1b2315 100644 --- a/docs/support_matrix.md +++ b/docs/support_matrix.md @@ -74,19 +74,19 @@ silent substitution with another core. | concurrent use of one plan | not supported | | CUDA Graph capture | not currently claimed/tested | | one transform across multiple GPUs | not supported | -| non-power-of-two transforms | not supported | -| multidimensional FFT/NTT | not supported | +| non-power-of-two transforms | C plan API: exact FP32/FP64 FFT via selected direct cuFFT/Bluestein, exact uint64 NTT when roots exist, or explicit power-of-two embedding for all operators | +| multidimensional transforms | C plan API: rank-two FFT, NTT, FWHT, zeta/Mobius, and Structured 2x2 via axis plans and transpose boundaries | | fused user epilogues | not supported | ## Automatic Selector Coverage -The checked-in selector recognizes V100 `sm_70` only and accepts a calibrated +The checked-in static table recognizes V100 `sm_70` only and accepts a calibrated subset of forward contiguous workloads. Its exact FFT, FWHT, legacy XOR-name and NTT shape table is documented in [Runtime Mapping Selector](runtime_selector.md). -Unsupported hardware or semantics throw `std::invalid_argument`; explicit -backend configuration remains available outside selector coverage. The new -canonical subset/superset and structured-2x2 names require explicit mapping -until their sweeps are recorded. +Unseen workloads use the legal resource-model fallback and emit a profile-miss +log; they can be measured and cached. Invalid mathematical domains still fail +plan creation. The new canonical subset/superset and structured-2x2 names +require measurement or explicit mapping until their sweeps are recorded. ## Versioning And ABI diff --git a/docs/v0.6_implementation_status.md b/docs/v0.6_implementation_status.md new file mode 100644 index 0000000..ce1e724 --- /dev/null +++ b/docs/v0.6_implementation_status.md @@ -0,0 +1,101 @@ +# v0.6 Implementation Status + +## Release Scope + +v0.6 extends the power-of-two kernel research artifact into a general-shape, +CUDA-library-style plan system. The release keeps operator arithmetic, +architecture mapping, physical processing units, composition boundaries, and +runtime selection as separate objects: + +```text +operator contract + -> exact or embedded shape lowering + -> per-axis (Ud,Td,Us,Ts) mapping + -> generated/imported processing unit + -> boundary residence and output layout + -> static, cached, measured, or explicit selection +``` + +This separation is the implementation form of the research claim. Direct +cuFFT, cuFFTDx, Hybrid2D, and register codelets are processing-unit candidates; +using one does not transfer ownership of the surrounding shape, layout, +workspace, logging, or selection contract. + +## Implemented Surface + +| Area | Implemented in v0.6 | +|:--|:--| +| C plan API | handle, immutable descriptor-derived plan, stream, caller workspace, asynchronous execution, status codes, selection queries, logs | +| exact lengths | packed FP32/FP64 FFT through direct cuFFT or Bluestein; uint64 NTT through modular Bluestein when roots exist | +| embedding | next-power-of-two zero extension for every public operator, with asymmetric forward/inverse logical and physical extents | +| rank two | per-axis FFT, NTT, FWHT, subset/superset zeta, and Structured 2x2 with tiled transpose boundaries | +| selection | exact V100 application profiles, exact-key user cache, bounded `MEASURE`, explicit algorithm IDs, auditable miss logs | +| composition | selected physical axes, selected Bluestein NTT core, direct physical output, and selected fused input | +| generator | warp-register FWHT/Structured boundary-aware variants and existing generated FFT processing units | + +The C++ `ButterflyPlan` remains the power-of-two physical-plan API. The C plan +API owns general-shape composition and should be preferred by applications that +need exact arbitrary length, embedding, rank two, or automatic composition +selection. + +## Boundary Lowering + +The implementation does not assume that every boundary should be fused. +Contiguous physical output removes the final scatter for rank-one embedding +and makes the final rank-two transpose land in caller storage. Strided output +retains scatter semantics. + +Forward FP32 FWHT with a generated warp-register unit additionally consumes +the logical input directly and synthesizes the zero tail in registers. It is a +single-kernel, zero-composition-workspace path and is reported as +`+direct-output+fused-input`. Structured 2x2 has a functionally correct +generated fused-input candidate, but the V100 measurement was slower than the +pack path, so the default selector does not enable it. This negative result is +part of the selection methodology: legal fusion is a candidate, not a rule. + +## V100 Evidence + +The matching-protocol ten-shape power-of-two matrix reports 2.349x geometric- +mean speedup from searched mappings over the fixed radix-2 base and 1.109x +geometric-mean throughput over the matched cuFFT, Dao FHT, and GPU-NTT rows. +That result is representative, not universal. + +The v0.6 composition follow-up records: + +| Workload | Selected result | Interpretation | +|:--|--:|:--| +| exact FP32 FFT `N=1000,batch=64` | 0.00427 ms, 0.998x matching cuFFT | direct vendor unit selected at parity | +| exact rank-2 FP32 FFT `30x50,batch=64` | 0.00981 ms, 1.000x matching cuFFT | rank-2 direct vendor unit selected at parity | +| uint64 NTT `N=5419,batch=64` | 0.21323 ms | modular Bluestein reuses Hybrid2D radix-4; no exact external baseline is archived | +| embedded FP32 FWHT `32767 -> 32768,batch=256` | 0.11759 ms | complete logical-input path is within 1.01x of the archived 0.11643 ms pre-padded Dao core | +| embedded FP32 FWHT `255 -> 256,batch=65536` | 0.16511 ms | fused input is 2.50x faster than the direct-output pack composition | + +Raw rows and protocol details are in [General-Shape V100 Selection +Results](general_shape_results.md). Power-of-two library claims remain governed +by the matching-protocol reports linked from [Experimental +Results](experiments.md). + +## Validation + +- 49 CTest entries pass on V100/CUDA 11.8. +- C API tests cover exact, embedded, rank-two, strided fallback, cache, and + forward/inverse contracts. +- Direct processing-unit tests cover zero-extended warp-register FWHT and + Structured 2x2 with distinct logical and physical batch distances. +- Generated-manifest validation, installed-package consumption, and static + application-profile generation are part of CTest. + +## Claim Boundary + +v0.6 does not establish universal superiority over specialist libraries or +cross-generation transfer. Exact FFT parity obtained by selecting direct +cuFFT is evidence that the composition can abstain to an appropriate unit, not +that a wrapper makes cuFFT faster. Exact arbitrary NTT depends on the finite- +field root domain. Rank-two execution retains a global transpose because +ordinary CTAs do not provide a conventional in-kernel grid barrier. + +The checked-in automatic table is V100-specific. An unseen workload or GPU +uses a legal static fallback, emits a profile-miss log, and can be measured and +cached. Cross-GPU profile capture remains future work; the next research step +is broader fixed-GPU regime mapping over numeric contract, length, batch, +stride, direction, placement, coefficient policy, and boundary residence. diff --git a/include/cubutterfly/cubutterfly.h b/include/cubutterfly/cubutterfly.h new file mode 100644 index 0000000..f93103f --- /dev/null +++ b/include/cubutterfly/cubutterfly.h @@ -0,0 +1,188 @@ +#pragma once + +#include + +#include +#include + +#if defined(_WIN32) +#define CUBUTTERFLYAPI __declspec(dllimport) +#else +#define CUBUTTERFLYAPI +#endif + +#ifdef __cplusplus +extern "C" { +#endif + +#define CUBUTTERFLY_VERSION_MAJOR 0 +#define CUBUTTERFLY_VERSION_MINOR 6 +#define CUBUTTERFLY_VERSION_PATCH 0 + +typedef struct cubutterflyHandle* cubutterflyHandle_t; +typedef struct cubutterflyDescriptor* cubutterflyDescriptor_t; +typedef struct cubutterflyPlan* cubutterflyPlan_t; + +typedef enum cubutterflyStatus { + CUBUTTERFLY_STATUS_SUCCESS = 0, + CUBUTTERFLY_STATUS_INVALID_VALUE, + CUBUTTERFLY_STATUS_INVALID_TYPE, + CUBUTTERFLY_STATUS_NOT_SUPPORTED, + CUBUTTERFLY_STATUS_ALLOC_FAILED, + CUBUTTERFLY_STATUS_CUDA_ERROR, + CUBUTTERFLY_STATUS_BACKEND_ERROR, + CUBUTTERFLY_STATUS_INTERNAL_ERROR +} cubutterflyStatus_t; + +typedef enum cubutterflyOperator { + CUBUTTERFLY_OPERATOR_FFT = 0, + CUBUTTERFLY_OPERATOR_NTT, + CUBUTTERFLY_OPERATOR_FWHT, + CUBUTTERFLY_OPERATOR_SUBSET_ZETA, + CUBUTTERFLY_OPERATOR_SUPERSET_ZETA, + CUBUTTERFLY_OPERATOR_STRUCTURED_2X2 +} cubutterflyOperator_t; + +typedef enum cubutterflyDataType { + CUBUTTERFLY_DATA_FP16 = 0, + CUBUTTERFLY_DATA_BF16, + CUBUTTERFLY_DATA_FP32, + CUBUTTERFLY_DATA_FP64, + CUBUTTERFLY_DATA_COMPLEX_FP16, + CUBUTTERFLY_DATA_COMPLEX_BF16, + CUBUTTERFLY_DATA_COMPLEX_FP32, + CUBUTTERFLY_DATA_COMPLEX_FP64, + CUBUTTERFLY_DATA_UINT32, + CUBUTTERFLY_DATA_UINT64 +} cubutterflyDataType_t; + +typedef enum cubutterflyComputeType { + CUBUTTERFLY_COMPUTE_NATIVE = 0, + CUBUTTERFLY_COMPUTE_FP32, + CUBUTTERFLY_COMPUTE_FP64, + CUBUTTERFLY_COMPUTE_UINT32, + CUBUTTERFLY_COMPUTE_UINT64 +} cubutterflyComputeType_t; + +typedef enum cubutterflyDirection { + CUBUTTERFLY_DIRECTION_FORWARD = 0, + CUBUTTERFLY_DIRECTION_INVERSE +} cubutterflyDirection_t; + +typedef enum cubutterflyPlacement { + CUBUTTERFLY_PLACEMENT_OUT_OF_PLACE = 0, + CUBUTTERFLY_PLACEMENT_IN_PLACE +} cubutterflyPlacement_t; + +typedef enum cubutterflyLengthMode { + CUBUTTERFLY_LENGTH_STANDARD = 0, + CUBUTTERFLY_LENGTH_ZERO_EXTENDED_EMBEDDING +} cubutterflyLengthMode_t; + +typedef enum cubutterflyAlgorithmPolicy { + CUBUTTERFLY_ALGORITHM_DEFAULT = 0, + CUBUTTERFLY_ALGORITHM_MEASURE, + CUBUTTERFLY_ALGORITHM_EXPLICIT +} cubutterflyAlgorithmPolicy_t; + +typedef enum cubutterflySelectionSource { + CUBUTTERFLY_SELECTION_MEASURED_TABLE = 0, + CUBUTTERFLY_SELECTION_USER_CACHE, + CUBUTTERFLY_SELECTION_RUNTIME_MEASURED, + CUBUTTERFLY_SELECTION_STATIC_MODEL, + CUBUTTERFLY_SELECTION_EXPLICIT +} cubutterflySelectionSource_t; + +typedef enum cubutterflyLogLevel { + CUBUTTERFLY_LOG_INFO = 0, + CUBUTTERFLY_LOG_WARNING, + CUBUTTERFLY_LOG_ERROR +} cubutterflyLogLevel_t; + +typedef struct cubutterflyMatrix2x2 { + double m00; + double m01; + double m10; + double m11; +} cubutterflyMatrix2x2_t; + +typedef void (*cubutterflyLogCallback_t)(cubutterflyLogLevel_t level, + const char* event, + const char* message, + void* user_data); + +CUBUTTERFLYAPI const char* cubutterflyGetStatusString(cubutterflyStatus_t status); +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflyGetVersion(int* major, int* minor, int* patch); + +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflyCreate(cubutterflyHandle_t* handle); +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflyDestroy(cubutterflyHandle_t handle); +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflySetStream(cubutterflyHandle_t handle, cudaStream_t stream); +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflyGetStream(cubutterflyHandle_t handle, cudaStream_t* stream); +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflySetLogCallback(cubutterflyHandle_t handle, + cubutterflyLogCallback_t callback, + void* user_data); +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflySetCachePath(cubutterflyHandle_t handle, const char* path); + +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflyCreateDescriptor(cubutterflyDescriptor_t* descriptor); +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflyDestroyDescriptor(cubutterflyDescriptor_t descriptor); +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflySetOperator(cubutterflyDescriptor_t descriptor, + cubutterflyOperator_t op); +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflySetDataType(cubutterflyDescriptor_t descriptor, + cubutterflyDataType_t storage, + cubutterflyComputeType_t compute); +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflySetShape(cubutterflyDescriptor_t descriptor, + uint32_t rank, + const size_t* extents, + size_t batch); +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflySetStrides(cubutterflyDescriptor_t descriptor, + const size_t* input_strides, + size_t input_batch_stride, + const size_t* output_strides, + size_t output_batch_stride); +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflySetDirection(cubutterflyDescriptor_t descriptor, + cubutterflyDirection_t direction, + int normalize_inverse); +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflySetPlacement(cubutterflyDescriptor_t descriptor, + cubutterflyPlacement_t placement); +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflySetLengthMode(cubutterflyDescriptor_t descriptor, + cubutterflyLengthMode_t mode); +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflySetAlgorithmPolicy(cubutterflyDescriptor_t descriptor, + cubutterflyAlgorithmPolicy_t policy, + const char* explicit_algorithm); +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflySetModulus(cubutterflyDescriptor_t descriptor, + uint64_t modulus, + uint32_t word_bits); +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflySetStageMatrices(cubutterflyDescriptor_t descriptor, + uint32_t axis, + const cubutterflyMatrix2x2_t* matrices, + size_t count); + +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflyCreatePlan(cubutterflyHandle_t handle, + cubutterflyDescriptor_t descriptor, + cubutterflyPlan_t* plan); +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflyDestroyPlan(cubutterflyPlan_t plan); +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflyPlanGetWorkspaceSize(cubutterflyPlan_t plan, size_t* bytes); +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflyPlanGetInputSize(cubutterflyPlan_t plan, size_t* bytes); +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflyPlanGetOutputSize(cubutterflyPlan_t plan, size_t* bytes); +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflyPlanGetPhysicalExtents(cubutterflyPlan_t plan, + size_t* extents, + uint32_t capacity); +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflyPlanGetSelectionSource(cubutterflyPlan_t plan, + cubutterflySelectionSource_t* source); +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflyPlanGetAlgorithmName(cubutterflyPlan_t plan, + char* name, + size_t* bytes); +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflyPlanGetSelectionReason(cubutterflyPlan_t plan, + char* reason, + size_t* bytes); +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflyPlanSetWorkspace(cubutterflyPlan_t plan, + void* workspace, + size_t bytes); +CUBUTTERFLYAPI cubutterflyStatus_t cubutterflyExecute(cubutterflyHandle_t handle, + cubutterflyPlan_t plan, + const void* input, + void* output); + +#ifdef __cplusplus +} +#endif diff --git a/include/cuntt/butterfly.hpp b/include/cuntt/butterfly.hpp index d030475..5b8bbb9 100644 --- a/include/cuntt/butterfly.hpp +++ b/include/cuntt/butterfly.hpp @@ -306,6 +306,14 @@ class ButterflyPlan { void execute_async(const ComplexBf16* input, ComplexBf16* output); void execute_async(const std::uint32_t* input, std::uint32_t* output); + // Specialized boundary lowering for a packed power-of-two FP32 + // warp-register plan. Values beyond logical_points are read as zero while + // the complete physical transform is written to the plan's output layout. + void execute_zero_extended_async(const float* input, float* output, + std::size_t logical_points, + std::size_t input_batch_stride, + std::size_t input_element_stride = 1); + ButterflyStats execute(const std::vector& input, std::vector& output, std::uint32_t warmup = 1, std::uint32_t repeat = 1); ButterflyStats execute(const std::vector& input, std::vector& output, std::uint32_t warmup = 1, std::uint32_t repeat = 1); ButterflyStats execute(const std::vector& input, std::vector& output, std::uint32_t warmup = 1, std::uint32_t repeat = 1); diff --git a/results/README.md b/results/README.md index 85b0dc6..3a08d53 100644 --- a/results/README.md +++ b/results/README.md @@ -30,6 +30,27 @@ replace rows in the canonical full raw data. ## Focused Evidence +`general_shapes_v100_smoke.csv` is the first seven-row v0.6 exact-length FFT +diagnostic against direct cuFFT. It uses one aggregated CUDA-event run per row, +so it locates a large implementation gap but is not canonical repeated-trial +evidence. Its protocol and interpretation are in +`docs/general_shape_results.md`. + +`general_shape_selection_v100_raw.csv` is the five-trial follow-up after +composition selection. It covers direct-versus-Bluestein FFT, selected +Bluestein NTT cores, embedded FWHT/zeta, and default versus selected Structured +2x2. Runs are separate warmed processes and are interpreted in the same report. + +`general_shape_direct_boundary_v100_raw.csv` isolates the next lowering step: +contiguous physical output replaces the final scatter for embedded FWHT, +subset/superset zeta, and Structured 2x2. The noncontiguous fallback remains a +C API correctness test. + +`general_shape_fused_input_v100_raw.csv` records the selected first-boundary +fusion for FP32 warp-register FWHT and the rejected Structured control. The +FWHT rows use five independent warmed processes with 100 timed launches; the +algorithm suffix and zero workspace verify that the pack kernel is absent. + Files prefixed with `fft_`, `cubutterfly_`, `processing_units_`, `ntt_`, or `stage_pipeline_` are focused design-space experiments. Their protocol and interpretation live in the matching document under `docs/`; they must not be diff --git a/results/general_shape_direct_boundary_v100_raw.csv b/results/general_shape_direct_boundary_v100_raw.csv new file mode 100644 index 0000000..193b7df --- /dev/null +++ b/results/general_shape_direct_boundary_v100_raw.csv @@ -0,0 +1,26 @@ +operator,precision,logical_shape,batch,policy,trial,algorithm,physical_shape,kernel_ms +fwht,fp32,32767,256,default,1,zero-extended[temporal-warp-register-radix2]+direct-output,32768,0.242842 +fwht,fp32,32767,256,default,2,zero-extended[temporal-warp-register-radix2]+direct-output,32768,0.242708 +fwht,fp32,32767,256,default,3,zero-extended[temporal-warp-register-radix2]+direct-output,32768,0.242586 +fwht,fp32,32767,256,default,4,zero-extended[temporal-warp-register-radix2]+direct-output,32768,0.242657 +fwht,fp32,32767,256,default,5,zero-extended[temporal-warp-register-radix2]+direct-output,32768,0.242729 +fwht,fp32,255,65536,default,1,zero-extended[temporal-warp-register-radix2]+direct-output,256,0.412979 +fwht,fp32,255,65536,default,2,zero-extended[temporal-warp-register-radix2]+direct-output,256,0.412692 +fwht,fp32,255,65536,default,3,zero-extended[temporal-warp-register-radix2]+direct-output,256,0.412293 +fwht,fp32,255,65536,default,4,zero-extended[temporal-warp-register-radix2]+direct-output,256,0.412314 +fwht,fp32,255,65536,default,5,zero-extended[temporal-warp-register-radix2]+direct-output,256,0.412314 +subset-zeta,uint32,32767,256,default,1,zero-extended[hierarchical-radix4]+direct-output,32768,0.648817 +subset-zeta,uint32,32767,256,default,2,zero-extended[hierarchical-radix4]+direct-output,32768,0.648817 +subset-zeta,uint32,32767,256,default,3,zero-extended[hierarchical-radix4]+direct-output,32768,0.648837 +subset-zeta,uint32,32767,256,default,4,zero-extended[hierarchical-radix4]+direct-output,32768,0.648806 +subset-zeta,uint32,32767,256,default,5,zero-extended[hierarchical-radix4]+direct-output,32768,0.648929 +superset-zeta,uint32,32767,256,default,1,zero-extended[hierarchical-radix4]+direct-output,32768,0.649032 +superset-zeta,uint32,32767,256,default,2,zero-extended[hierarchical-radix4]+direct-output,32768,0.649001 +superset-zeta,uint32,32767,256,default,3,zero-extended[hierarchical-radix4]+direct-output,32768,0.648847 +superset-zeta,uint32,32767,256,default,4,zero-extended[hierarchical-radix4]+direct-output,32768,0.649011 +superset-zeta,uint32,32767,256,default,5,zero-extended[hierarchical-radix4]+direct-output,32768,0.649001 +structured-2x2,fp32,32767,256,explicit,1,zero-extended[temporal-warp-register-radix2]+direct-output,32768,0.257444 +structured-2x2,fp32,32767,256,explicit,2,zero-extended[temporal-warp-register-radix2]+direct-output,32768,0.258212 +structured-2x2,fp32,32767,256,explicit,3,zero-extended[temporal-warp-register-radix2]+direct-output,32768,0.257618 +structured-2x2,fp32,32767,256,explicit,4,zero-extended[temporal-warp-register-radix2]+direct-output,32768,0.257905 +structured-2x2,fp32,32767,256,explicit,5,zero-extended[temporal-warp-register-radix2]+direct-output,32768,0.25772 diff --git a/results/general_shape_fused_input_v100_raw.csv b/results/general_shape_fused_input_v100_raw.csv new file mode 100644 index 0000000..4193192 --- /dev/null +++ b/results/general_shape_fused_input_v100_raw.csv @@ -0,0 +1,12 @@ +gpu,operator,logical_shape,physical_shape,batch,algorithm,trial,kernel_ms +V100-SXM2-32GB,fwht-fp32,32767,32768,256,zero-extended-temporal-warp-register-radix2-direct-output-fused-input,1,0.117586 +V100-SXM2-32GB,fwht-fp32,32767,32768,256,zero-extended-temporal-warp-register-radix2-direct-output-fused-input,2,0.117248 +V100-SXM2-32GB,fwht-fp32,32767,32768,256,zero-extended-temporal-warp-register-radix2-direct-output-fused-input,3,0.117586 +V100-SXM2-32GB,fwht-fp32,32767,32768,256,zero-extended-temporal-warp-register-radix2-direct-output-fused-input,4,0.117647 +V100-SXM2-32GB,fwht-fp32,32767,32768,256,zero-extended-temporal-warp-register-radix2-direct-output-fused-input,5,0.117443 +V100-SXM2-32GB,fwht-fp32,255,256,65536,zero-extended-temporal-warp-register-radix2-direct-output-fused-input,1,0.165100 +V100-SXM2-32GB,fwht-fp32,255,256,65536,zero-extended-temporal-warp-register-radix2-direct-output-fused-input,2,0.165161 +V100-SXM2-32GB,fwht-fp32,255,256,65536,zero-extended-temporal-warp-register-radix2-direct-output-fused-input,3,0.165059 +V100-SXM2-32GB,fwht-fp32,255,256,65536,zero-extended-temporal-warp-register-radix2-direct-output-fused-input,4,0.165140 +V100-SXM2-32GB,fwht-fp32,255,256,65536,zero-extended-temporal-warp-register-radix2-direct-output-fused-input,5,0.164772 +V100-SXM2-32GB,structured-2x2-fp32,32767,32768,256,zero-extended-temporal-warp-register-radix2-direct-output,1,0.350198 diff --git a/results/general_shape_selection_v100_raw.csv b/results/general_shape_selection_v100_raw.csv new file mode 100644 index 0000000..2fa3c01 --- /dev/null +++ b/results/general_shape_selection_v100_raw.csv @@ -0,0 +1,51 @@ +operator,precision,logical_shape,batch,length_mode,policy,trial,algorithm,physical_shape,kernel_ms,cufft_ms +fft,fp32,1000,64,standard,default,1,cufft-direct,1000,0.00427008,0.00427008 +fft,fp32,1000,64,standard,default,2,cufft-direct,1000,0.00427008,0.00427008 +fft,fp32,1000,64,standard,default,3,cufft-direct,1000,0.00428032,0.00427008 +fft,fp32,1000,64,standard,default,4,cufft-direct,1000,0.00427008,0.00428032 +fft,fp32,1000,64,standard,default,5,cufft-direct,1000,0.00427008,0.00428032 +fft,fp32,30x50,64,standard,default,1,cufft-direct,30x50,0.00979968,0.0098304 +fft,fp32,30x50,64,standard,default,2,cufft-direct,30x50,0.00980992,0.00980992 +fft,fp32,30x50,64,standard,default,3,cufft-direct,30x50,0.00980992,0.00980992 +fft,fp32,30x50,64,standard,default,4,cufft-direct,30x50,0.00978944,0.00980992 +fft,fp32,30x50,64,standard,default,5,cufft-direct,30x50,0.00982016,0.00982016 +ntt,uint64,5419,1,standard,default,1,bluestein-ntt[ntt-hybrid-radix4],16384,0.0314573, +ntt,uint64,5419,1,standard,default,2,bluestein-ntt[ntt-hybrid-radix4],16384,0.0314778, +ntt,uint64,5419,1,standard,default,3,bluestein-ntt[ntt-hybrid-radix4],16384,0.0314675, +ntt,uint64,5419,1,standard,default,4,bluestein-ntt[ntt-hybrid-radix4],16384,0.0314675, +ntt,uint64,5419,1,standard,default,5,bluestein-ntt[ntt-hybrid-radix4],16384,0.0314368, +ntt,uint64,5419,64,standard,default,1,bluestein-ntt[ntt-hybrid-radix4],16384,0.213228, +ntt,uint64,5419,64,standard,default,2,bluestein-ntt[ntt-hybrid-radix4],16384,0.213166, +ntt,uint64,5419,64,standard,default,3,bluestein-ntt[ntt-hybrid-radix4],16384,0.21333, +ntt,uint64,5419,64,standard,default,4,bluestein-ntt[ntt-hybrid-radix4],16384,0.213228, +ntt,uint64,5419,64,standard,default,5,bluestein-ntt[ntt-hybrid-radix4],16384,0.212992, +fwht,fp32,32767,256,embedding,default,1,zero-extended[temporal-warp-register-radix2],32768,0.326543, +fwht,fp32,32767,256,embedding,default,2,zero-extended[temporal-warp-register-radix2],32768,0.323973, +fwht,fp32,32767,256,embedding,default,3,zero-extended[temporal-warp-register-radix2],32768,0.324239, +fwht,fp32,32767,256,embedding,default,4,zero-extended[temporal-warp-register-radix2],32768,0.324065, +fwht,fp32,32767,256,embedding,default,5,zero-extended[temporal-warp-register-radix2],32768,0.324086, +fwht,fp32,255,65536,embedding,default,1,zero-extended[temporal-warp-register-radix2],256,0.584264, +fwht,fp32,255,65536,embedding,default,2,zero-extended[temporal-warp-register-radix2],256,0.584172, +fwht,fp32,255,65536,embedding,default,3,zero-extended[temporal-warp-register-radix2],256,0.584182, +fwht,fp32,255,65536,embedding,default,4,zero-extended[temporal-warp-register-radix2],256,0.584264, +fwht,fp32,255,65536,embedding,default,5,zero-extended[temporal-warp-register-radix2],256,0.584223, +subset-zeta,uint32,32767,256,embedding,default,1,zero-extended[hierarchical-radix4],32768,0.735898, +subset-zeta,uint32,32767,256,embedding,default,2,zero-extended[hierarchical-radix4],32768,0.736133, +subset-zeta,uint32,32767,256,embedding,default,3,zero-extended[hierarchical-radix4],32768,0.736061, +subset-zeta,uint32,32767,256,embedding,default,4,zero-extended[hierarchical-radix4],32768,0.736031, +subset-zeta,uint32,32767,256,embedding,default,5,zero-extended[hierarchical-radix4],32768,0.736051, +superset-zeta,uint32,32767,256,embedding,default,1,zero-extended[hierarchical-radix4],32768,0.736072, +superset-zeta,uint32,32767,256,embedding,default,2,zero-extended[hierarchical-radix4],32768,0.736092, +superset-zeta,uint32,32767,256,embedding,default,3,zero-extended[hierarchical-radix4],32768,0.736133, +superset-zeta,uint32,32767,256,embedding,default,4,zero-extended[hierarchical-radix4],32768,0.73601, +superset-zeta,uint32,32767,256,embedding,default,5,zero-extended[hierarchical-radix4],32768,0.736215, +structured-2x2,fp32,32767,256,embedding,default,1,zero-extended[hierarchical-radix4],32768,0.744632, +structured-2x2,fp32,32767,256,embedding,default,2,zero-extended[hierarchical-radix4],32768,0.744643, +structured-2x2,fp32,32767,256,embedding,default,3,zero-extended[hierarchical-radix4],32768,0.744376, +structured-2x2,fp32,32767,256,embedding,default,4,zero-extended[hierarchical-radix4],32768,0.744602, +structured-2x2,fp32,32767,256,embedding,default,5,zero-extended[hierarchical-radix4],32768,0.744561, +structured-2x2,fp32,32767,256,embedding,explicit,1,zero-extended[temporal-warp-register-radix2],32768,0.345221, +structured-2x2,fp32,32767,256,embedding,explicit,2,zero-extended[temporal-warp-register-radix2],32768,0.344248, +structured-2x2,fp32,32767,256,embedding,explicit,3,zero-extended[temporal-warp-register-radix2],32768,0.344238, +structured-2x2,fp32,32767,256,embedding,explicit,4,zero-extended[temporal-warp-register-radix2],32768,0.344699, +structured-2x2,fp32,32767,256,embedding,explicit,5,zero-extended[temporal-warp-register-radix2],32768,0.343921, diff --git a/results/general_shapes_v100_smoke.csv b/results/general_shapes_v100_smoke.csv new file mode 100644 index 0000000..7a9b5ca --- /dev/null +++ b/results/general_shapes_v100_smoke.csv @@ -0,0 +1,8 @@ +device,cuda,operator,precision,logical_shape,physical_shape,batch,warmup,repeat,algorithm,workspace_bytes,cubutterfly_ms,cufft_ms,speedup_vs_cufft +Tesla V100-SXM2-16GB,11.8,fft,complex-fp32,5,16,1,20,100,bluestein-cufft-power2-core,128,0.0153395,0.00274432,0.178905 +Tesla V100-SXM2-16GB,11.8,fft,complex-fp32,5,16,4096,20,100,bluestein-cufft-power2-core,524288,0.0187802,0.00271360,0.144493 +Tesla V100-SXM2-16GB,11.8,fft,complex-fp32,1000,2048,1,20,100,bluestein-cufft-power2-core,16384,0.0204800,0.00336896,0.164500 +Tesla V100-SXM2-16GB,11.8,fft,complex-fp32,1000,2048,64,20,100,bluestein-cufft-power2-core,1048576,0.0262246,0.00428032,0.163217 +Tesla V100-SXM2-16GB,11.8,fft,complex-fp32,3x5,8x16,4096,20,100,rank2-transpose-bluestein-cufft-power2-core,2555904,0.0655258,0.00690176,0.105329 +Tesla V100-SXM2-16GB,11.8,fft,complex-fp32,30x50,64x128,1,20,100,rank2-transpose-bluestein-cufft-power2-core,54720,0.0468992,0.00788480,0.168122 +Tesla V100-SXM2-16GB,11.8,fft,complex-fp32,30x50,64x128,64,20,100,rank2-transpose-bluestein-cufft-power2-core,3502080,0.0799334,0.00978944,0.122470 diff --git a/scripts/generate_application_profiles.py b/scripts/generate_application_profiles.py new file mode 100644 index 0000000..d01fdc6 --- /dev/null +++ b/scripts/generate_application_profiles.py @@ -0,0 +1,61 @@ +#!/usr/bin/env python3 +import argparse +import json +import pathlib + + +OPERATORS = {"fft": 0, "ntt": 1, "fwht": 2} +STORAGE = {"fp32": 2, "complex-fp32": 6, "uint32": 8, "uint64": 9} +DIRECTION = {"forward": 0, "inverse": 1} +PLACEMENT = {"out-of-place": 0, "in-place": 1} + + +def generate(spec_path: pathlib.Path, header_path: pathlib.Path) -> None: + spec = json.loads(spec_path.read_text()) + if spec.get("schema_version") != 1: + raise ValueError("application profile schema_version must be 1") + target = spec["target"] + profiles = spec["profiles"] + lines = [ + "// Generated from config/v100_application_profiles.json. Do not edit.", + "#pragma once", + "#include ", + "#include ", + "#include ", + "namespace cuntt::detail {", + "struct GeneratedApplicationProfile { int op; int storage; int direction; int placement; std::uint64_t modulus; std::uint32_t log_n; const std::size_t* batches; std::size_t batch_count; const char* implementation; };", + f"inline constexpr int kApplicationProfileComputeMajor = {int(target['compute_major'])};", + f"inline constexpr int kApplicationProfileComputeMinor = {int(target['compute_minor'])};", + f"inline constexpr int kApplicationProfileMultiprocessors = {int(target['multiprocessors'])};", + f"inline constexpr const char* kApplicationProfileTarget = \"{target['name']}\";", + ] + for index, profile in enumerate(profiles): + batches = profile["batches"] + if not batches or any(int(value) <= 0 for value in batches): + raise ValueError(f"profile {index} has invalid batches") + values = ", ".join(f"{int(value)}ULL" for value in batches) + lines.append(f"inline constexpr std::array kApplicationBatches{index} = {{{{{values}}}}};") + lines.append(f"inline constexpr std::array kApplicationProfiles = {{{{") + for index, profile in enumerate(profiles): + lines.append( + " GeneratedApplicationProfile{" + + f"{OPERATORS[profile['operator']]}, {STORAGE[profile['storage']]}, " + + f"{DIRECTION[profile['direction']]}, {PLACEMENT[profile['placement']]}, " + + f"{int(profile.get('modulus', 0))}ULL, {int(profile['logN'])}, kApplicationBatches{index}.data(), kApplicationBatches{index}.size(), " + + f"\"{profile['implementation']}\"" + "}," + ) + lines.extend(["}};", "} // namespace cuntt::detail", ""]) + header_path.parent.mkdir(parents=True, exist_ok=True) + header_path.write_text("\n".join(lines)) + + +def main() -> None: + parser = argparse.ArgumentParser(description="Generate built-in application workload profiles.") + parser.add_argument("--spec", type=pathlib.Path, required=True) + parser.add_argument("--header", type=pathlib.Path, required=True) + args = parser.parse_args() + generate(args.spec, args.header) + + +if __name__ == "__main__": + main() diff --git a/scripts/generate_design_points.py b/scripts/generate_design_points.py index eabf41a..532890b 100755 --- a/scripts/generate_design_points.py +++ b/scripts/generate_design_points.py @@ -55,6 +55,9 @@ def availability(log_ns): def write_fwht(path, log_ns): cases = switch_cases(log_ns, "launch_register_fwht<{log_n}>(input, output, transforms, batch_distance, element_stride, normalize, stream)") + zero_extended_cases = switch_cases( + log_ns, + "launch_register_fwht_zero_extended<{log_n}>(input, output, transforms, input_batch_distance, output_batch_distance, input_element_stride, output_element_stride, logical_points, stream)") thread_cases = "\n".join( f" case {log_n}: return register_fwht_threads({log_n});" for log_n in sorted(log_ns)) path.write_text(f'''// Generated by scripts/generate_design_points.py. Do not edit. @@ -78,6 +81,17 @@ def write_fwht(path, log_ns): default: throw std::invalid_argument("requested FWHT warp-register design point was not generated"); }} }} +void launch_generated_register_fwht_zero_extended( + std::uint32_t log_n, const float* input, float* output, + std::uint64_t transforms, std::uint64_t input_batch_distance, + std::uint64_t output_batch_distance, std::uint64_t input_element_stride, + std::uint64_t output_element_stride, std::uint32_t logical_points, + cudaStream_t stream) {{ + switch (log_n) {{ +{zero_extended_cases} + default: throw std::invalid_argument("requested zero-extended FWHT design point was not generated"); + }} +}} }} // namespace cuntt::detail ''') @@ -86,6 +100,9 @@ def write_structured(path, log_ns): cases = switch_cases( log_ns, "launch_register_structured<{log_n}>(input, output, matrices, transforms, batch_distance, element_stride, broadcast, stream)") + zero_extended_cases = switch_cases( + log_ns, + "launch_register_structured_zero_extended<{log_n}>(input, output, matrices, transforms, input_batch_distance, output_batch_distance, input_element_stride, output_element_stride, logical_points, broadcast, stream)") thread_cases = "\n".join( f" case {log_n}: return register_fwht_threads({log_n});" for log_n in sorted(log_ns)) path.write_text(f'''// Generated by scripts/generate_design_points.py. Do not edit. @@ -110,6 +127,17 @@ def write_structured(path, log_ns): default: throw std::invalid_argument("requested structured-2x2 warp-register design point was not generated"); }} }} +void launch_generated_register_structured_zero_extended( + std::uint32_t log_n, const float* input, float* output, + const float* matrices, std::uint64_t transforms, + std::uint64_t input_batch_distance, std::uint64_t output_batch_distance, + std::uint64_t input_element_stride, std::uint64_t output_element_stride, + std::uint32_t logical_points, bool broadcast, cudaStream_t stream) {{ + switch (log_n) {{ +{zero_extended_cases} + default: throw std::invalid_argument("requested zero-extended structured design point was not generated"); + }} +}} }} // namespace cuntt::detail ''') diff --git a/src/bluestein_fft.cu b/src/bluestein_fft.cu new file mode 100644 index 0000000..105a4eb --- /dev/null +++ b/src/bluestein_fft.cu @@ -0,0 +1,270 @@ +#include "bluestein_fft.hpp" + +#include "cuntt/butterfly.hpp" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace cuntt::detail { +namespace { + +void check_cuda(cudaError_t status, const char* operation) { + if (status != cudaSuccess) + throw std::runtime_error(std::string(operation) + ": " + cudaGetErrorString(status)); +} + +void check_cufft(cufftResult status, const char* operation) { + if (status != CUFFT_SUCCESS) + throw std::runtime_error(std::string(operation) + ": cuFFT status " + std::to_string(status)); +} + +template +struct ComplexOps; + +template <> +struct ComplexOps { + __host__ __device__ static Complex32 make(double real, double imag) { + return {static_cast(real), static_cast(imag)}; + } +}; + +template <> +struct ComplexOps { + __host__ __device__ static Complex64 make(double real, double imag) { return {real, imag}; } +}; + +template +__global__ void bluestein_pack_kernel(const Complex* input, Complex* work, + std::size_t n, std::size_t m, std::size_t batch, + std::size_t input_stride, std::size_t input_batch_stride, + double sign) { + const std::size_t linear = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + const std::size_t count = n * batch; + if (linear >= count) + return; + const std::size_t transform = linear / n; + const std::size_t index = linear - transform * n; + const auto value = input[transform * input_batch_stride + index * input_stride]; + const double angle = sign * 3.141592653589793238462643383279502884 * + static_cast(index) * static_cast(index) / static_cast(n); + const double c = cos(angle); + const double s = sin(angle); + work[transform * m + index] = ComplexOps::make(value.real * c - value.imag * s, + value.real * s + value.imag * c); +} + +template +__global__ void bluestein_multiply_kernel(Complex* work, const Complex* kernel, + std::size_t m, std::size_t batch) { + const std::size_t linear = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + if (linear >= m * batch) + return; + const auto left = work[linear]; + const auto right = kernel[linear % m]; + work[linear] = ComplexOps::make(left.real * right.real - left.imag * right.imag, + left.real * right.imag + left.imag * right.real); +} + +template +__global__ void bluestein_output_kernel(const Complex* work, Complex* output, + std::size_t n, std::size_t m, std::size_t batch, + std::size_t output_stride, std::size_t output_batch_stride, + double sign, double scale) { + const std::size_t linear = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + const std::size_t count = n * batch; + if (linear >= count) + return; + const std::size_t transform = linear / n; + const std::size_t index = linear - transform * n; + const auto value = work[transform * m + index]; + const double angle = sign * 3.141592653589793238462643383279502884 * + static_cast(index) * static_cast(index) / static_cast(n); + const double c = cos(angle) * scale; + const double s = sin(angle) * scale; + output[transform * output_batch_stride + index * output_stride] = + ComplexOps::make(value.real * c - value.imag * s, + value.real * s + value.imag * c); +} + +unsigned int blocks_for(std::size_t count) { + const std::size_t blocks = (count + 255) / 256; + if (blocks > UINT32_MAX) + throw std::overflow_error("Bluestein CUDA grid exceeds grid.x"); + return static_cast(blocks); +} + +} // namespace + +class BluesteinFftPlan::Impl { + public: + Impl(std::size_t n, std::size_t m, std::size_t batch, bool inverse, + bool normalize_inverse, bool fp64) + : n_(n), m_(m), batch_(batch), inverse_(inverse), + normalize_inverse_(normalize_inverse), fp64_(fp64) { + if (n_ == 0 || m_ < 2 * n_ - 1 || batch_ == 0 || m_ > static_cast(INT32_MAX)) + throw std::invalid_argument("invalid Bluestein shape"); + int length = static_cast(m_); + check_cufft(cufftPlanMany(&plan_, 1, &length, nullptr, 1, length, nullptr, 1, length, + fp64_ ? CUFFT_Z2Z : CUFFT_C2C, static_cast(batch_)), + "create Bluestein batch plan"); + try { + build_kernel(); + } catch (...) { + cufftDestroy(plan_); + plan_ = 0; + cudaFree(kernel_fft_); + kernel_fft_ = nullptr; + throw; + } + } + + ~Impl() { + if (kernel_fft_ != nullptr) + cudaFree(kernel_fft_); + if (plan_ != 0) + cufftDestroy(plan_); + } + + std::size_t workspace_size() const noexcept { + return m_ * batch_ * (fp64_ ? sizeof(Complex64) : sizeof(Complex32)); + } + + void execute(const void* input, void* output, void* workspace, + std::size_t input_stride, std::size_t input_batch_stride, + std::size_t output_stride, std::size_t output_batch_stride, + cudaStream_t stream) { + if (input == nullptr || output == nullptr || workspace == nullptr) + throw std::invalid_argument("Bluestein input, output, and workspace must be non-null"); + check_cufft(cufftSetStream(plan_, stream), "set Bluestein stream"); + check_cuda(cudaMemsetAsync(workspace, 0, workspace_size(), stream), "clear Bluestein workspace"); + if (fp64_) + execute_typed(static_cast(input), static_cast(output), + static_cast(workspace), input_stride, input_batch_stride, + output_stride, output_batch_stride, stream); + else + execute_typed(static_cast(input), static_cast(output), + static_cast(workspace), input_stride, input_batch_stride, + output_stride, output_batch_stride, stream); + } + + private: + template + void execute_typed(const Complex* input, Complex* output, Complex* workspace, + std::size_t input_stride, std::size_t input_batch_stride, + std::size_t output_stride, std::size_t output_batch_stride, + cudaStream_t stream) { + const double sign = inverse_ ? 1.0 : -1.0; + bluestein_pack_kernel<<>>( + input, workspace, n_, m_, batch_, input_stride, input_batch_stride, sign); + check_cuda(cudaGetLastError(), "launch Bluestein pack"); + if constexpr (std::is_same_v) + check_cufft(cufftExecZ2Z(plan_, reinterpret_cast(workspace), + reinterpret_cast(workspace), CUFFT_FORWARD), + "execute Bluestein forward FFT"); + else + check_cufft(cufftExecC2C(plan_, reinterpret_cast(workspace), + reinterpret_cast(workspace), CUFFT_FORWARD), + "execute Bluestein forward FFT"); + bluestein_multiply_kernel<<>>( + workspace, static_cast(kernel_fft_), m_, batch_); + check_cuda(cudaGetLastError(), "launch Bluestein pointwise multiply"); + if constexpr (std::is_same_v) + check_cufft(cufftExecZ2Z(plan_, reinterpret_cast(workspace), + reinterpret_cast(workspace), CUFFT_INVERSE), + "execute Bluestein inverse FFT"); + else + check_cufft(cufftExecC2C(plan_, reinterpret_cast(workspace), + reinterpret_cast(workspace), CUFFT_INVERSE), + "execute Bluestein inverse FFT"); + const double scale = 1.0 / static_cast(m_) * + (inverse_ && normalize_inverse_ ? 1.0 / static_cast(n_) : 1.0); + bluestein_output_kernel<<>>( + workspace, output, n_, m_, batch_, output_stride, output_batch_stride, sign, scale); + check_cuda(cudaGetLastError(), "launch Bluestein output"); + } + + template + void build_kernel_typed() { + std::vector values(m_, ComplexOps::make(0.0, 0.0)); + const double sign = inverse_ ? 1.0 : -1.0; + for (std::size_t index = 0; index < n_; ++index) { + const double angle = -sign * 3.141592653589793238462643383279502884 * + static_cast(index) * static_cast(index) / static_cast(n_); + const auto value = ComplexOps::make(std::cos(angle), std::sin(angle)); + values[index] = value; + if (index != 0) + values[m_ - index] = value; + } + const std::size_t bytes = m_ * sizeof(Complex); + void* temporary = nullptr; + cufftHandle kernel_plan = 0; + check_cuda(cudaMalloc(&temporary, bytes), "allocate Bluestein kernel temporary"); + try { + check_cuda(cudaMalloc(&kernel_fft_, bytes), "allocate Bluestein kernel FFT"); + check_cuda(cudaMemcpy(temporary, values.data(), bytes, cudaMemcpyHostToDevice), "upload Bluestein kernel"); + int length = static_cast(m_); + check_cufft(cufftPlan1d(&kernel_plan, length, + fp64_ ? CUFFT_Z2Z : CUFFT_C2C, 1), "create Bluestein kernel plan"); + if constexpr (std::is_same_v) + check_cufft(cufftExecZ2Z(kernel_plan, static_cast(temporary), + static_cast(kernel_fft_), CUFFT_FORWARD), + "transform Bluestein kernel"); + else + check_cufft(cufftExecC2C(kernel_plan, static_cast(temporary), + static_cast(kernel_fft_), CUFFT_FORWARD), + "transform Bluestein kernel"); + check_cuda(cudaDeviceSynchronize(), "wait for Bluestein kernel initialization"); + cufftDestroy(kernel_plan); + cudaFree(temporary); + } catch (...) { + if (kernel_plan != 0) + cufftDestroy(kernel_plan); + cudaFree(temporary); + throw; + } + } + + void build_kernel() { + if (fp64_) + build_kernel_typed(); + else + build_kernel_typed(); + } + + std::size_t n_; + std::size_t m_; + std::size_t batch_; + bool inverse_; + bool normalize_inverse_; + bool fp64_; + cufftHandle plan_ = 0; + void* kernel_fft_ = nullptr; +}; + +BluesteinFftPlan::BluesteinFftPlan(std::size_t length, std::size_t convolution_length, + std::size_t batch, bool inverse, bool normalize_inverse, + bool fp64) + : impl_(std::make_unique(length, convolution_length, batch, inverse, + normalize_inverse, fp64)) {} + +BluesteinFftPlan::~BluesteinFftPlan() = default; + +std::size_t BluesteinFftPlan::workspace_size() const noexcept { return impl_->workspace_size(); } + +void BluesteinFftPlan::execute(const void* input, void* output, void* workspace, + std::size_t input_stride, std::size_t input_batch_stride, + std::size_t output_stride, std::size_t output_batch_stride, + cudaStream_t stream) { + impl_->execute(input, output, workspace, input_stride, input_batch_stride, + output_stride, output_batch_stride, stream); +} + +} // namespace cuntt::detail diff --git a/src/bluestein_fft.hpp b/src/bluestein_fft.hpp new file mode 100644 index 0000000..ecdeff6 --- /dev/null +++ b/src/bluestein_fft.hpp @@ -0,0 +1,32 @@ +#pragma once + +#include + +#include +#include +#include + +namespace cuntt::detail { + +class BluesteinFftPlan { + public: + BluesteinFftPlan(std::size_t length, std::size_t convolution_length, + std::size_t batch, bool inverse, bool normalize_inverse, + bool fp64); + ~BluesteinFftPlan(); + + BluesteinFftPlan(const BluesteinFftPlan&) = delete; + BluesteinFftPlan& operator=(const BluesteinFftPlan&) = delete; + + std::size_t workspace_size() const noexcept; + void execute(const void* input, void* output, void* workspace, + std::size_t input_stride, std::size_t input_batch_stride, + std::size_t output_stride, std::size_t output_batch_stride, + cudaStream_t stream); + + private: + class Impl; + std::unique_ptr impl_; +}; + +} // namespace cuntt::detail diff --git a/src/bluestein_ntt.cu b/src/bluestein_ntt.cu new file mode 100644 index 0000000..8a15887 --- /dev/null +++ b/src/bluestein_ntt.cu @@ -0,0 +1,303 @@ +#include "bluestein_ntt.hpp" + +#include "cuntt/ntt.hpp" + +#include + +#include +#include +#include +#include +#include + +namespace cuntt::detail { +namespace { + +void check_cuda(cudaError_t status, const char* operation) { + if (status != cudaSuccess) + throw std::runtime_error(std::string(operation) + ": " + cudaGetErrorString(status)); +} + +std::uint64_t shoup_host(std::uint64_t value, std::uint64_t modulus) { + return static_cast((static_cast(value) << 64) / modulus); +} + +std::vector prime_factors(std::uint64_t value) { + std::vector factors; + for (std::uint64_t factor = 2; factor <= value / factor; ++factor) { + if (value % factor != 0) + continue; + factors.push_back(factor); + while (value % factor == 0) + value /= factor; + } + if (value > 1) + factors.push_back(value); + return factors; +} + +std::uint64_t find_root(std::uint64_t order, std::uint64_t modulus) { + if (!cuntt::is_prime(modulus)) + throw std::invalid_argument("Bluestein NTT modulus must be prime"); + if (order < 2 || (modulus - 1) % order != 0) + throw std::invalid_argument("modulus does not admit the required 2N-th Bluestein root"); + const auto factors = prime_factors(order); + for (std::uint64_t candidate = 2; candidate < (1ULL << 20); ++candidate) { + const auto root = cuntt::mod_pow(candidate, (modulus - 1) / order, modulus); + bool primitive = cuntt::mod_pow(root, order, modulus) == 1; + for (const auto factor : factors) + primitive = primitive && cuntt::mod_pow(root, order / factor, modulus) != 1; + if (primitive) + return root; + } + throw std::runtime_error("could not find the required Bluestein NTT root"); +} + +__device__ __forceinline__ std::uint64_t mul_shoup(std::uint64_t value, + std::uint64_t multiplier, + std::uint64_t multiplier_shoup, + std::uint64_t modulus) { + const std::uint64_t quotient = __umul64hi(value, multiplier_shoup); + std::uint64_t reduced = value * multiplier - quotient * modulus; + if (reduced >= modulus) + reduced -= modulus; + return reduced; +} + +__global__ void pack_kernel(const std::uint64_t* input, std::uint64_t* work, + const std::uint64_t* chirp, const std::uint64_t* chirp_shoup, + std::size_t n, std::size_t m, std::size_t batch, + std::size_t input_stride, std::size_t input_batch_stride, + std::uint64_t modulus) { + const std::size_t linear = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + if (linear >= n * batch) + return; + const std::size_t transform = linear / n; + const std::size_t index = linear % n; + const auto value = input[transform * input_batch_stride + index * input_stride] % modulus; + work[transform * m + index] = mul_shoup(value, chirp[index], chirp_shoup[index], modulus); +} + +__global__ void multiply_kernel(std::uint64_t* values, const std::uint64_t* kernel, + const std::uint64_t* kernel_shoup, std::size_t m, + std::size_t batch, std::uint64_t modulus) { + const std::size_t linear = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + if (linear >= m * batch) + return; + const std::size_t index = linear % m; + values[linear] = mul_shoup(values[linear], kernel[index], kernel_shoup[index], modulus); +} + +__global__ void output_kernel(const std::uint64_t* work, std::uint64_t* output, + const std::uint64_t* chirp, const std::uint64_t* chirp_shoup, + std::size_t n, std::size_t m, std::size_t batch, + std::size_t output_stride, std::size_t output_batch_stride, + std::uint64_t scale, std::uint64_t scale_shoup, + std::uint64_t modulus) { + const std::size_t linear = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + if (linear >= n * batch) + return; + const std::size_t transform = linear / n; + const std::size_t index = linear % n; + auto value = mul_shoup(work[transform * m + index], chirp[index], chirp_shoup[index], modulus); + value = mul_shoup(value, scale, scale_shoup, modulus); + output[transform * output_batch_stride + index * output_stride] = value; +} + +unsigned int blocks_for(std::size_t count) { + const auto blocks = (count + 255) / 256; + if (blocks > UINT_MAX) + throw std::overflow_error("Bluestein NTT CUDA grid exceeds grid.x"); + return static_cast(blocks); +} + +std::size_t align_16(std::size_t value) { return (value + 15) & ~std::size_t{15}; } + +} // namespace + +class BluesteinNttPlan::Impl { + public: + Impl(std::size_t n, std::size_t m, std::size_t batch, std::uint64_t modulus, bool inverse) + : Impl(n, m, batch, modulus, inverse, cuntt::Backend::Baseline, + cuntt::ComputeUnit::Radix2) {} + + Impl(std::size_t n, std::size_t m, std::size_t batch, std::uint64_t modulus, + bool inverse, cuntt::Backend core_backend, cuntt::ComputeUnit core_compute_unit) + : n_(n), m_(m), batch_(batch), modulus_(modulus), inverse_(inverse), + core_backend_(core_backend), core_compute_unit_(core_compute_unit), + forward_(core_config(m, batch, modulus, false, core_backend, core_compute_unit)), + inverse_core_(core_config(m, batch, modulus, true, core_backend, core_compute_unit)) { + if (n_ < 2 || m_ < 2 * n_ - 1 || batch_ == 0) + throw std::invalid_argument("invalid Bluestein NTT shape"); + if ((modulus_ - 1) % m_ != 0) + throw std::invalid_argument("modulus does not admit the power-of-two Bluestein convolution root"); + try { + build_tables(); + } catch (...) { + cudaFree(kernel_shoup_); + cudaFree(kernel_fft_); + cudaFree(chirp_shoup_); + cudaFree(chirp_); + kernel_shoup_ = nullptr; + kernel_fft_ = nullptr; + chirp_shoup_ = nullptr; + chirp_ = nullptr; + throw; + } + } + + ~Impl() { + cudaFree(kernel_shoup_); + cudaFree(kernel_fft_); + cudaFree(chirp_shoup_); + cudaFree(chirp_); + } + + std::size_t workspace_size() const noexcept { + return 2 * data_bytes() + align_16(std::max(forward_.workspace_size(), inverse_core_.workspace_size())); + } + + void execute(const std::uint64_t* input, std::uint64_t* output, void* workspace, + std::size_t input_stride, std::size_t input_batch_stride, + std::size_t output_stride, std::size_t output_batch_stride, + cudaStream_t stream) { + if (input == nullptr || output == nullptr || workspace == nullptr) + throw std::invalid_argument("Bluestein NTT input, output, and workspace must be non-null"); + auto* first = static_cast(workspace); + auto* second = reinterpret_cast(static_cast(workspace) + data_bytes()); + auto* core_workspace = static_cast(workspace) + 2 * data_bytes(); + const auto core_bytes = std::max(forward_.workspace_size(), inverse_core_.workspace_size()); + forward_.set_stream(stream); + inverse_core_.set_stream(stream); + if (core_bytes != 0) { + forward_.set_workspace(core_workspace, core_bytes); + inverse_core_.set_workspace(core_workspace, core_bytes); + } + check_cuda(cudaMemsetAsync(first, 0, data_bytes(), stream), "clear Bluestein NTT input"); + pack_kernel<<>>( + input, first, chirp_, chirp_shoup_, n_, m_, batch_, input_stride, + input_batch_stride, modulus_); + check_cuda(cudaGetLastError(), "launch Bluestein NTT pack"); + forward_.execute_async(first, second); + multiply_kernel<<>>( + second, kernel_fft_, kernel_shoup_, m_, batch_, modulus_); + check_cuda(cudaGetLastError(), "launch Bluestein NTT multiply"); + inverse_core_.execute_async(second, first); + output_kernel<<>>( + first, output, chirp_, chirp_shoup_, n_, m_, batch_, output_stride, + output_batch_stride, scale_, scale_shoup_, modulus_); + check_cuda(cudaGetLastError(), "launch Bluestein NTT output"); + } + + private: + static cuntt::PlanConfig core_config(std::size_t m, std::size_t batch, + std::uint64_t modulus, bool inverse, + cuntt::Backend backend, + cuntt::ComputeUnit compute_unit) { + if (m == 0 || (m & (m - 1)) != 0) + throw std::invalid_argument("Bluestein NTT convolution length must be a power of two"); + std::uint32_t log_m = 0; + for (auto value = m; value > 1; value >>= 1) + ++log_m; + cuntt::PlanConfig config; + config.log_n = log_m; + config.batch = batch; + config.modulus = modulus; + config.inverse = inverse; + config.backend = backend; + config.compute_unit = compute_unit; + config.word_bits = 64; + config.auto_allocate_workspace = false; + return config; + } + + std::size_t data_bytes() const noexcept { return m_ * batch_ * sizeof(std::uint64_t); } + + void build_tables() { + auto root = find_root(2 * n_, modulus_); + if (inverse_) + root = cuntt::mod_pow(root, modulus_ - 2, modulus_); + std::vector chirp(n_); + std::vector chirp_shoup(n_); + std::vector kernel(m_, 0); + for (std::size_t index = 0; index < n_; ++index) { + const auto exponent = static_cast( + (static_cast(index) * index) % (2 * n_)); + chirp[index] = cuntt::mod_pow(root, exponent, modulus_); + chirp_shoup[index] = shoup_host(chirp[index], modulus_); + const auto inverse_chirp = cuntt::mod_pow(chirp[index], modulus_ - 2, modulus_); + kernel[index] = inverse_chirp; + if (index != 0) + kernel[m_ - index] = inverse_chirp; + } + const auto chirp_bytes = n_ * sizeof(std::uint64_t); + const auto kernel_bytes = m_ * sizeof(std::uint64_t); + check_cuda(cudaMalloc(&chirp_, chirp_bytes), "allocate Bluestein NTT chirp"); + check_cuda(cudaMalloc(&chirp_shoup_, chirp_bytes), "allocate Bluestein NTT chirp Shoup table"); + check_cuda(cudaMalloc(&kernel_fft_, kernel_bytes), "allocate Bluestein NTT kernel"); + check_cuda(cudaMalloc(&kernel_shoup_, kernel_bytes), "allocate Bluestein NTT kernel Shoup table"); + std::uint64_t* temporary = nullptr; + check_cuda(cudaMalloc(&temporary, kernel_bytes), "allocate Bluestein NTT initialization buffer"); + try { + check_cuda(cudaMemcpy(chirp_, chirp.data(), chirp_bytes, cudaMemcpyHostToDevice), "upload Bluestein NTT chirp"); + check_cuda(cudaMemcpy(chirp_shoup_, chirp_shoup.data(), chirp_bytes, cudaMemcpyHostToDevice), "upload Bluestein NTT chirp Shoup table"); + check_cuda(cudaMemcpy(temporary, kernel.data(), kernel_bytes, cudaMemcpyHostToDevice), "upload Bluestein NTT kernel"); + auto init_config = core_config(m_, 1, modulus_, false, core_backend_, core_compute_unit_); + init_config.auto_allocate_workspace = true; + cuntt::Plan init_plan(init_config); + init_plan.execute_async(temporary, kernel_fft_); + check_cuda(cudaDeviceSynchronize(), "initialize Bluestein NTT kernel transform"); + std::vector transformed(m_); + check_cuda(cudaMemcpy(transformed.data(), kernel_fft_, kernel_bytes, cudaMemcpyDeviceToHost), + "download Bluestein NTT kernel transform"); + std::vector transformed_shoup(m_); + std::transform(transformed.begin(), transformed.end(), transformed_shoup.begin(), + [&](std::uint64_t value) { return shoup_host(value, modulus_); }); + check_cuda(cudaMemcpy(kernel_shoup_, transformed_shoup.data(), kernel_bytes, cudaMemcpyHostToDevice), + "upload Bluestein NTT kernel Shoup table"); + cudaFree(temporary); + } catch (...) { + cudaFree(temporary); + throw; + } + scale_ = inverse_ ? cuntt::mod_pow(static_cast(n_), modulus_ - 2, modulus_) : 1; + scale_shoup_ = shoup_host(scale_, modulus_); + } + + std::size_t n_; + std::size_t m_; + std::size_t batch_; + std::uint64_t modulus_; + bool inverse_; + cuntt::Backend core_backend_; + cuntt::ComputeUnit core_compute_unit_; + cuntt::Plan forward_; + cuntt::Plan inverse_core_; + std::uint64_t* chirp_ = nullptr; + std::uint64_t* chirp_shoup_ = nullptr; + std::uint64_t* kernel_fft_ = nullptr; + std::uint64_t* kernel_shoup_ = nullptr; + std::uint64_t scale_ = 1; + std::uint64_t scale_shoup_ = 0; +}; + +BluesteinNttPlan::BluesteinNttPlan(std::size_t length, std::size_t convolution_length, + std::size_t batch, std::uint64_t modulus, bool inverse, + cuntt::Backend core_backend, + cuntt::ComputeUnit core_compute_unit) + : impl_(std::make_unique(length, convolution_length, batch, modulus, inverse, + core_backend, core_compute_unit)) {} + +BluesteinNttPlan::~BluesteinNttPlan() = default; + +std::size_t BluesteinNttPlan::workspace_size() const noexcept { return impl_->workspace_size(); } + +void BluesteinNttPlan::execute(const std::uint64_t* input, std::uint64_t* output, + void* workspace, std::size_t input_stride, + std::size_t input_batch_stride, std::size_t output_stride, + std::size_t output_batch_stride, cudaStream_t stream) { + impl_->execute(input, output, workspace, input_stride, input_batch_stride, + output_stride, output_batch_stride, stream); +} + +} // namespace cuntt::detail diff --git a/src/bluestein_ntt.hpp b/src/bluestein_ntt.hpp new file mode 100644 index 0000000..b9d3c95 --- /dev/null +++ b/src/bluestein_ntt.hpp @@ -0,0 +1,35 @@ +#pragma once + +#include "cuntt/ntt.hpp" + +#include + +#include +#include +#include + +namespace cuntt::detail { + +class BluesteinNttPlan { + public: + BluesteinNttPlan(std::size_t length, std::size_t convolution_length, + std::size_t batch, std::uint64_t modulus, + bool inverse, cuntt::Backend core_backend, + cuntt::ComputeUnit core_compute_unit); + ~BluesteinNttPlan(); + + BluesteinNttPlan(const BluesteinNttPlan&) = delete; + BluesteinNttPlan& operator=(const BluesteinNttPlan&) = delete; + + std::size_t workspace_size() const noexcept; + void execute(const std::uint64_t* input, std::uint64_t* output, + void* workspace, std::size_t input_stride, + std::size_t input_batch_stride, std::size_t output_stride, + std::size_t output_batch_stride, cudaStream_t stream); + + private: + class Impl; + std::unique_ptr impl_; +}; + +} // namespace cuntt::detail diff --git a/src/butterfly.cu b/src/butterfly.cu index 1222adc..f75a70b 100644 --- a/src/butterfly.cu +++ b/src/butterfly.cu @@ -1324,6 +1324,47 @@ class ButterflyPlan::Impl { launch_once(input, output, scratch_buffer(), auxiliary_buffer()); } + void execute_zero_extended_async(const float* input, float* output, + std::size_t logical_points, + std::size_t input_batch_stride, + std::size_t input_element_stride) { + if (input == nullptr || output == nullptr) + throw std::invalid_argument("zero-extended butterfly input and output must be non-null"); + if (input == output) + throw std::invalid_argument("zero-extended butterfly execution requires distinct input and output pointers"); + if (config_.precision != ButterflyPrecision::Fp32 || + (config_.op != ButterflyOperator::Fwht && config_.op != ButterflyOperator::Structured2x2) || + config_.inverse || config_.local_exchange != LocalExchange::WarpRegister) { + throw std::invalid_argument( + "zero-extended fused input requires a forward FP32 warp-register FWHT or structured-2x2 plan"); + } + if (logical_points == 0 || logical_points > points_) + throw std::invalid_argument("zero-extended logical length must be in the physical transform range"); + if (input_element_stride == 0 || input_batch_stride == 0) + throw std::invalid_argument("zero-extended input strides must be positive"); + if (logical_points > 1 && + (input_element_stride > std::numeric_limits::max() / (logical_points - 1) || + input_batch_stride <= (logical_points - 1) * input_element_stride)) + throw std::invalid_argument("zero-extended input batches overlap"); + + if (config_.op == ButterflyOperator::Structured2x2) { + detail::launch_generated_register_structured_zero_extended( + config_.log_n, input, output, + reinterpret_cast(device_operator_coefficients_.data()), + config_.batch, input_batch_stride, config_.batch_stride, + input_element_stride, config_.element_stride, + static_cast(logical_points), + config_.stage_matrices.size() == 1, stream_); + } else { + detail::launch_generated_register_fwht_zero_extended( + config_.log_n, input, output, config_.batch, + input_batch_stride, config_.batch_stride, + input_element_stride, config_.element_stride, + static_cast(logical_points), stream_); + } + CUB_CUDA_CHECK(cudaGetLastError()); + } + template ButterflyStats execute(const std::vector& input, std::vector& output, std::uint32_t warmup, std::uint32_t repeat, ButterflyOperator required_op, ButterflyPrecision required_precision) { @@ -1881,6 +1922,14 @@ void ButterflyPlan::execute_async(const std::uint32_t* input, std::uint32_t* out impl_->execute_async(input, output, ButterflyOperator::XorZeta, ButterflyPrecision::Uint32); } +void ButterflyPlan::execute_zero_extended_async(const float* input, float* output, + std::size_t logical_points, + std::size_t input_batch_stride, + std::size_t input_element_stride) { + impl_->execute_zero_extended_async(input, output, logical_points, + input_batch_stride, input_element_stride); +} + ButterflyStats ButterflyPlan::execute(const std::vector& input, std::vector& output, std::uint32_t warmup, std::uint32_t repeat) { return impl_->execute(input, output, warmup, repeat, ButterflyOperator::Fwht, ButterflyPrecision::Fp32); } diff --git a/src/c_api.cpp b/src/c_api.cpp new file mode 100644 index 0000000..a589dc2 --- /dev/null +++ b/src/c_api.cpp @@ -0,0 +1,1753 @@ +#include "cubutterfly/cubutterfly.h" + +#include "cuntt/butterfly.hpp" +#include "cuntt/ntt.hpp" +#include "bluestein_fft.hpp" +#include "bluestein_ntt.hpp" +#include "direct_fft.hpp" +#include "general_kernels.hpp" +#include "generated_application_profiles.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifndef CUBUTTERFLY_HAS_CUFFTDX +#define CUBUTTERFLY_HAS_CUFFTDX 0 +#endif + +struct cubutterflyHandle { + cudaStream_t stream = nullptr; + cubutterflyLogCallback_t logger = nullptr; + void* logger_data = nullptr; + std::string cache_path; +}; + +struct cubutterflyDescriptor { + cubutterflyOperator_t op = CUBUTTERFLY_OPERATOR_FFT; + cubutterflyDataType_t storage = CUBUTTERFLY_DATA_COMPLEX_FP32; + cubutterflyComputeType_t compute = CUBUTTERFLY_COMPUTE_FP32; + cubutterflyDirection_t direction = CUBUTTERFLY_DIRECTION_FORWARD; + cubutterflyPlacement_t placement = CUBUTTERFLY_PLACEMENT_OUT_OF_PLACE; + cubutterflyLengthMode_t length_mode = CUBUTTERFLY_LENGTH_STANDARD; + cubutterflyAlgorithmPolicy_t algorithm_policy = CUBUTTERFLY_ALGORITHM_DEFAULT; + std::string explicit_algorithm; + std::uint32_t rank = 1; + std::vector extents{256}; + std::vector input_strides; + std::vector output_strides; + std::size_t input_batch_stride = 0; + std::size_t output_batch_stride = 0; + std::size_t batch = 1; + bool normalize_inverse = true; + std::uint64_t modulus = cuntt::kDefaultModulus; + std::uint32_t word_bits = 64; + std::vector> stage_matrices{2}; +}; + +struct cubutterflyPlan { + enum class Kind { Butterfly, Ntt, DirectFft, BluesteinFft, BluesteinNtt } kind = Kind::Butterfly; + std::unique_ptr butterfly; + std::unique_ptr ntt; + std::vector> bluestein_axes; + std::vector> bluestein_ntt_axes; + std::unique_ptr direct_fft; + std::vector> butterfly_axes; + std::vector> ntt_axes; + cubutterflyDescriptor descriptor; + bool composite = false; + bool direct_output_boundary = false; + bool fused_input_boundary = false; + cubutterflySelectionSource_t selection_source = CUBUTTERFLY_SELECTION_STATIC_MODEL; + std::string algorithm; + std::string reason; + std::vector physical_extents; + std::size_t input_bytes = 0; + std::size_t output_bytes = 0; + std::size_t element_bytes = 0; + std::size_t packed_bytes = 0; + std::size_t composite_buffers = 0; + std::size_t workspace_bytes = 0; + void* workspace = nullptr; +}; + +namespace { + +template +cubutterflyStatus_t protect(Function&& function) noexcept { + try { + function(); + return CUBUTTERFLY_STATUS_SUCCESS; + } catch (const std::invalid_argument&) { + return CUBUTTERFLY_STATUS_INVALID_VALUE; + } catch (const std::overflow_error&) { + return CUBUTTERFLY_STATUS_INVALID_VALUE; + } catch (const std::bad_alloc&) { + return CUBUTTERFLY_STATUS_ALLOC_FAILED; + } catch (const std::logic_error&) { + return CUBUTTERFLY_STATUS_INTERNAL_ERROR; + } catch (const std::runtime_error&) { + return CUBUTTERFLY_STATUS_BACKEND_ERROR; + } catch (...) { + return CUBUTTERFLY_STATUS_INTERNAL_ERROR; + } +} + +void require(bool condition, const char* message) { + if (!condition) + throw std::invalid_argument(message); +} + +std::string workload_key(const cubutterflyDescriptor& descriptor); + +bool is_power_of_two(std::size_t value) { + return value != 0 && (value & (value - 1)) == 0; +} + +std::uint32_t integer_log2(std::size_t value) { + require(is_power_of_two(value), "extent must be a power of two"); + std::uint32_t result = 0; + while (value > 1) { + value >>= 1; + ++result; + } + return result; +} + +std::size_t next_power_of_two(std::size_t value) { + require(value != 0, "extent must be positive"); + if (is_power_of_two(value)) + return value; + require(value <= (static_cast(-1) >> 1) + 1, "power-of-two extent overflows size_t"); + --value; + for (std::size_t shift = 1; shift < sizeof(value) * 8; shift <<= 1) + value |= value >> shift; + return value + 1; +} + +void log_event(cubutterflyHandle_t handle, cubutterflyLogLevel_t level, + const char* event, const std::string& message) { + if (handle != nullptr && handle->logger != nullptr) + handle->logger(level, event, message.c_str(), handle->logger_data); +} + +std::size_t element_bytes(cubutterflyDataType_t type) { + switch (type) { + case CUBUTTERFLY_DATA_FP16: + case CUBUTTERFLY_DATA_BF16: + return 2; + case CUBUTTERFLY_DATA_FP32: + case CUBUTTERFLY_DATA_COMPLEX_FP16: + case CUBUTTERFLY_DATA_COMPLEX_BF16: + case CUBUTTERFLY_DATA_UINT32: + return 4; + case CUBUTTERFLY_DATA_FP64: + case CUBUTTERFLY_DATA_COMPLEX_FP32: + case CUBUTTERFLY_DATA_UINT64: + return 8; + case CUBUTTERFLY_DATA_COMPLEX_FP64: + return 16; + } + throw std::invalid_argument("unknown data type"); +} + +std::size_t strided_elements(const cubutterflyDescriptor& descriptor, bool input, + const std::vector& extents) { + const auto& strides = input ? descriptor.input_strides : descriptor.output_strides; + const auto batch_stride = input ? descriptor.input_batch_stride : descriptor.output_batch_stride; + require(strides.size() == descriptor.rank, "stride rank does not match shape rank"); + std::size_t extent = 1; + for (std::uint32_t axis = 0; axis < descriptor.rank; ++axis) { + require(strides[axis] != 0, "strides must be positive"); + const std::size_t add = (extents[axis] - 1) * strides[axis]; + require(extent <= static_cast(-1) - add, "strided extent overflows size_t"); + extent += add; + } + if (descriptor.batch > 1) { + require(batch_stride >= extent, "batch stride is smaller than one transform extent"); + const std::size_t add = (descriptor.batch - 1) * batch_stride; + require(extent <= static_cast(-1) - add, "batch extent overflows size_t"); + extent += add; + } + return extent; +} + +void resolve_default_strides(cubutterflyDescriptor& descriptor, + const std::vector& input_extents, + const std::vector& output_extents) { + auto resolve = [&](std::vector& strides, std::size_t& batch_stride, + const std::vector& extents) { + if (strides.empty()) { + strides.assign(descriptor.rank, 1); + for (std::uint32_t axis = descriptor.rank; axis-- > 1;) + strides[axis - 1] = strides[axis] * extents[axis]; + } + require(strides.size() == descriptor.rank, "stride rank does not match shape rank"); + if (batch_stride == 0) { + batch_stride = 1; + for (std::uint32_t axis = 0; axis < descriptor.rank; ++axis) + batch_stride += (extents[axis] - 1) * strides[axis]; + } + }; + resolve(descriptor.input_strides, descriptor.input_batch_stride, input_extents); + resolve(descriptor.output_strides, descriptor.output_batch_stride, output_extents); +} + +bool packed_layout(const cubutterflyDescriptor& descriptor) { + std::size_t stride = 1; + for (std::uint32_t axis = descriptor.rank; axis-- > 0;) { + if (descriptor.input_strides[axis] != stride || descriptor.output_strides[axis] != stride) + return false; + stride *= descriptor.extents[axis]; + } + return descriptor.input_batch_stride == stride && descriptor.output_batch_stride == stride; +} + +bool packed_side(const cubutterflyDescriptor& descriptor, bool input, + const std::vector& extents) { + const auto& strides = input ? descriptor.input_strides : descriptor.output_strides; + const auto batch_stride = input ? descriptor.input_batch_stride + : descriptor.output_batch_stride; + std::size_t stride = 1; + for (std::uint32_t axis = descriptor.rank; axis-- > 0;) { + if (strides[axis] != stride) + return false; + stride *= extents[axis]; + } + return batch_stride == stride; +} + +cuntt::ButterflyOperator butterfly_operator(cubutterflyOperator_t op) { + switch (op) { + case CUBUTTERFLY_OPERATOR_FFT: return cuntt::ButterflyOperator::Fft; + case CUBUTTERFLY_OPERATOR_FWHT: return cuntt::ButterflyOperator::Fwht; + case CUBUTTERFLY_OPERATOR_SUBSET_ZETA: return cuntt::ButterflyOperator::SubsetZeta; + case CUBUTTERFLY_OPERATOR_SUPERSET_ZETA: return cuntt::ButterflyOperator::SupersetZeta; + case CUBUTTERFLY_OPERATOR_STRUCTURED_2X2: return cuntt::ButterflyOperator::Structured2x2; + case CUBUTTERFLY_OPERATOR_NTT: break; + } + throw std::invalid_argument("NTT is not a common butterfly plan operator"); +} + +cuntt::ButterflyPrecision butterfly_precision(const cubutterflyDescriptor& descriptor) { + switch (descriptor.storage) { + case CUBUTTERFLY_DATA_FP16: + case CUBUTTERFLY_DATA_COMPLEX_FP16: return cuntt::ButterflyPrecision::Fp16; + case CUBUTTERFLY_DATA_BF16: + case CUBUTTERFLY_DATA_COMPLEX_BF16: return cuntt::ButterflyPrecision::Bf16; + case CUBUTTERFLY_DATA_FP32: + case CUBUTTERFLY_DATA_COMPLEX_FP32: return cuntt::ButterflyPrecision::Fp32; + case CUBUTTERFLY_DATA_FP64: + case CUBUTTERFLY_DATA_COMPLEX_FP64: return cuntt::ButterflyPrecision::Fp64; + case CUBUTTERFLY_DATA_UINT32: return cuntt::ButterflyPrecision::Uint32; + case CUBUTTERFLY_DATA_UINT64: break; + } + throw std::invalid_argument("data type is not valid for a common butterfly plan"); +} + +void validate_storage(const cubutterflyDescriptor& descriptor) { + const bool complex = descriptor.storage == CUBUTTERFLY_DATA_COMPLEX_FP16 || + descriptor.storage == CUBUTTERFLY_DATA_COMPLEX_BF16 || + descriptor.storage == CUBUTTERFLY_DATA_COMPLEX_FP32 || + descriptor.storage == CUBUTTERFLY_DATA_COMPLEX_FP64; + const bool real = descriptor.storage == CUBUTTERFLY_DATA_FP16 || + descriptor.storage == CUBUTTERFLY_DATA_BF16 || + descriptor.storage == CUBUTTERFLY_DATA_FP32 || + descriptor.storage == CUBUTTERFLY_DATA_FP64; + if (descriptor.op == CUBUTTERFLY_OPERATOR_FFT) + require(complex, "FFT requires a complex storage type"); + else if (descriptor.op == CUBUTTERFLY_OPERATOR_NTT) + require(descriptor.storage == CUBUTTERFLY_DATA_UINT32 || descriptor.storage == CUBUTTERFLY_DATA_UINT64, + "NTT requires uint32 or uint64 storage"); + else if (descriptor.op == CUBUTTERFLY_OPERATOR_SUBSET_ZETA || descriptor.op == CUBUTTERFLY_OPERATOR_SUPERSET_ZETA) + require(descriptor.storage == CUBUTTERFLY_DATA_UINT32, "zeta transforms require uint32 storage"); + else + require(real, "real butterfly operator requires a real floating-point storage type"); +} + +bool v100_profile_device() { + const auto device = cuntt::current_device_info(); + return device.compute_major == cuntt::detail::kApplicationProfileComputeMajor && + device.compute_minor == cuntt::detail::kApplicationProfileComputeMinor && + device.multiprocessors == cuntt::detail::kApplicationProfileMultiprocessors; +} + +bool measured_v100_shape(const cubutterflyDescriptor& descriptor, std::uint32_t log_n) { + if (!v100_profile_device() || descriptor.rank != 1 || + descriptor.length_mode != CUBUTTERFLY_LENGTH_STANDARD || + descriptor.direction != CUBUTTERFLY_DIRECTION_FORWARD) + return false; + if (descriptor.input_strides.size() != 1 || descriptor.output_strides.size() != 1 || + descriptor.input_strides[0] != 1 || descriptor.output_strides[0] != 1 || + descriptor.input_batch_stride != descriptor.extents[0] || + descriptor.output_batch_stride != descriptor.extents[0]) + return false; + if (descriptor.op == CUBUTTERFLY_OPERATOR_FFT && !CUBUTTERFLY_HAS_CUFFTDX) + return false; + for (const auto& profile : cuntt::detail::kApplicationProfiles) { + if (profile.op != static_cast(descriptor.op) || + profile.storage != static_cast(descriptor.storage) || + profile.direction != static_cast(descriptor.direction) || + profile.placement != static_cast(descriptor.placement) || profile.log_n != log_n) + continue; + if (profile.modulus != 0 && profile.modulus != descriptor.modulus) + continue; + if (std::find(profile.batches, profile.batches + profile.batch_count, descriptor.batch) != + profile.batches + profile.batch_count) + return true; + } + return false; +} + +bool measured_v100_axis(const cubutterflyDescriptor& descriptor, std::uint32_t log_n, + std::size_t batch) { + if (!v100_profile_device() || descriptor.direction != CUBUTTERFLY_DIRECTION_FORWARD) + return false; + if (descriptor.op == CUBUTTERFLY_OPERATOR_FFT && !CUBUTTERFLY_HAS_CUFFTDX) + return false; + for (const auto& profile : cuntt::detail::kApplicationProfiles) { + if (profile.op != static_cast(descriptor.op) || + profile.storage != static_cast(descriptor.storage) || + profile.direction != static_cast(descriptor.direction) || profile.log_n != log_n) + continue; + if (profile.modulus != 0 && profile.modulus != descriptor.modulus) + continue; + if (std::find(profile.batches, profile.batches + profile.batch_count, batch) != + profile.batches + profile.batch_count) + return true; + } + return false; +} + +void log_profile_miss(cubutterflyHandle_t handle, const cubutterflyDescriptor& descriptor, + const std::string& reason, const std::string& algorithm) { + const auto device = cuntt::current_device_info(); + std::ostringstream message; + if (!v100_profile_device()) { + message << "no static performance profile for device='" << device.name << "' sm_" + << device.compute_major << device.compute_minor << " multiprocessors=" + << device.multiprocessors; + log_event(handle, CUBUTTERFLY_LOG_WARNING, "hardware-profile-miss", + message.str() + "; using " + algorithm + "; set MEASURE with a cache path to calibrate this workload"); + return; + } + message << "V100 profile has no exact entry for " << workload_key(descriptor) + << "; " << reason << "; using " << algorithm + << "; set MEASURE with a cache path to record a reusable choice"; + log_event(handle, CUBUTTERFLY_LOG_WARNING, "workload-profile-miss", message.str()); +} + +void choose_butterfly_mapping(cuntt::ButterflyConfig& config, bool measured) { + config.compute_unit = cuntt::ComputeUnit::Radix4; + config.tile_threads = 128; + if (measured && config.op == cuntt::ButterflyOperator::Fwht && + config.precision == cuntt::ButterflyPrecision::Fp32) { + if (config.log_n == 8 || (config.log_n == 15 && config.batch >= 16)) { + config.backend = cuntt::ButterflyBackend::TemporalTile; + config.compute_unit = cuntt::ComputeUnit::Radix2; + config.local_exchange = cuntt::LocalExchange::WarpRegister; + config.tile_threads = 256; + return; + } + if (config.log_n == 15) { + config.backend = cuntt::ButterflyBackend::OnlineReorder; + config.local_stages = 8; + config.reorder_columns = 1; + config.tile_threads = 256; + return; + } + } + if (config.log_n <= 10) { + config.backend = cuntt::ButterflyBackend::TemporalTile; + return; + } + if (measured && config.log_n >= 18) { + config.backend = cuntt::ButterflyBackend::OnlineReorder; + config.local_stages = config.log_n / 2; + config.reorder_columns = 1; + config.tile_threads = 256; + return; + } + config.backend = cuntt::ButterflyBackend::Hierarchical; + config.local_stages = std::min(10, config.log_n - 1); + config.tile_threads = 256; +} + +void copy_text(const std::string& value, char* destination, std::size_t* bytes) { + require(bytes != nullptr, "size query pointer must be non-null"); + const std::size_t required = value.size() + 1; + if (destination == nullptr) { + *bytes = required; + return; + } + require(*bytes >= required, "destination string buffer is too small"); + std::memcpy(destination, value.c_str(), required); + *bytes = required; +} + +void execute_butterfly(cuntt::ButterflyPlan& plan, const void* input, void* output, + cubutterflyDataType_t storage) { + switch (storage) { + case CUBUTTERFLY_DATA_FP16: + plan.execute_async(static_cast(input), static_cast(output)); return; + case CUBUTTERFLY_DATA_BF16: + plan.execute_async(static_cast(input), static_cast(output)); return; + case CUBUTTERFLY_DATA_FP32: + plan.execute_async(static_cast(input), static_cast(output)); return; + case CUBUTTERFLY_DATA_FP64: + plan.execute_async(static_cast(input), static_cast(output)); return; + case CUBUTTERFLY_DATA_COMPLEX_FP16: + plan.execute_async(static_cast(input), static_cast(output)); return; + case CUBUTTERFLY_DATA_COMPLEX_BF16: + plan.execute_async(static_cast(input), static_cast(output)); return; + case CUBUTTERFLY_DATA_COMPLEX_FP32: + plan.execute_async(static_cast(input), static_cast(output)); return; + case CUBUTTERFLY_DATA_COMPLEX_FP64: + plan.execute_async(static_cast(input), static_cast(output)); return; + case CUBUTTERFLY_DATA_UINT32: + plan.execute_async(static_cast(input), static_cast(output)); return; + case CUBUTTERFLY_DATA_UINT64: break; + } + throw std::invalid_argument("invalid butterfly execution type"); +} + +cuntt::ButterflyConfig make_axis_butterfly_config(const cubutterflyDescriptor& descriptor, + std::uint32_t axis, std::size_t length, + std::size_t batch) { + cuntt::ButterflyConfig config; + config.op = butterfly_operator(descriptor.op); + config.precision = butterfly_precision(descriptor); + const bool low_precision = config.precision == cuntt::ButterflyPrecision::Fp16 || + config.precision == cuntt::ButterflyPrecision::Bf16; + config.accumulation = low_precision && descriptor.compute == CUBUTTERFLY_COMPUTE_FP32 + ? cuntt::ButterflyAccumulation::Fp32 + : cuntt::ButterflyAccumulation::Native; + config.placement = cuntt::ButterflyPlacement::InPlace; + config.log_n = integer_log2(length); + config.batch = batch; + config.batch_stride = length; + config.element_stride = 1; + config.inverse = descriptor.direction == CUBUTTERFLY_DIRECTION_INVERSE; + config.normalize_inverse = descriptor.normalize_inverse; + config.auto_allocate_workspace = false; + if (config.op == cuntt::ButterflyOperator::Structured2x2) { + require(axis < descriptor.stage_matrices.size() && !descriptor.stage_matrices[axis].empty(), + "structured-2x2 requires matrices for every executed axis"); + for (const auto& matrix : descriptor.stage_matrices[axis]) + config.stage_matrices.push_back({matrix.m00, matrix.m01, matrix.m10, matrix.m11}); + } + choose_butterfly_mapping(config, false); + return config; +} + +cuntt::PlanConfig make_axis_ntt_config(const cubutterflyDescriptor& descriptor, + std::size_t length, std::size_t batch) { + cuntt::PlanConfig config; + config.log_n = integer_log2(length); + config.batch = batch; + config.modulus = descriptor.modulus; + config.inverse = descriptor.direction == CUBUTTERFLY_DIRECTION_INVERSE; + config.word_bits = descriptor.word_bits; + config.auto_allocate_workspace = false; + if (config.log_n >= 12 && config.log_n <= 20) { + config.backend = cuntt::Backend::Hybrid2D; + config.compute_unit = cuntt::ComputeUnit::Radix4; + } else { + config.backend = cuntt::Backend::Baseline; + config.compute_unit = cuntt::ComputeUnit::Radix2; + } + return config; +} + +std::size_t align_16(std::size_t value) { + require(value <= static_cast(-1) - 15, "workspace size overflows size_t"); + return (value + 15) & ~std::size_t{15}; +} + +const char* operator_name(cubutterflyOperator_t op) { + switch (op) { + case CUBUTTERFLY_OPERATOR_FFT: return "fft"; + case CUBUTTERFLY_OPERATOR_NTT: return "ntt"; + case CUBUTTERFLY_OPERATOR_FWHT: return "fwht"; + case CUBUTTERFLY_OPERATOR_SUBSET_ZETA: return "subset-zeta"; + case CUBUTTERFLY_OPERATOR_SUPERSET_ZETA: return "superset-zeta"; + case CUBUTTERFLY_OPERATOR_STRUCTURED_2X2: return "structured-2x2"; + } + return "unknown"; +} + +std::string workload_key(const cubutterflyDescriptor& descriptor) { + const auto device = cuntt::current_device_info(); + std::ostringstream key; + key << "v1|sm" << device.compute_major << device.compute_minor << "-mp" << device.multiprocessors + << '|' << operator_name(descriptor.op) << "|storage" << static_cast(descriptor.storage) + << "|compute" << static_cast(descriptor.compute) << "|direction" << static_cast(descriptor.direction) + << "|placement" << static_cast(descriptor.placement) << "|length-mode" << static_cast(descriptor.length_mode) + << "|rank" << descriptor.rank << "|shape"; + for (std::uint32_t axis = 0; axis < descriptor.rank; ++axis) { + if (axis != 0) + key << 'x'; + key << descriptor.extents[axis]; + } + key << "|batch" << descriptor.batch << "|in"; + for (const auto stride : descriptor.input_strides) + key << stride << ','; + key << descriptor.input_batch_stride << "|out"; + for (const auto stride : descriptor.output_strides) + key << stride << ','; + key << descriptor.output_batch_stride; + if (descriptor.op == CUBUTTERFLY_OPERATOR_NTT) + key << "|modulus" << descriptor.modulus << "|word" << descriptor.word_bits; + if (descriptor.op == CUBUTTERFLY_OPERATOR_STRUCTURED_2X2) { + key << "|matrices" << std::setprecision(std::numeric_limits::max_digits10); + for (std::uint32_t axis = 0; axis < descriptor.rank; ++axis) { + if (axis != 0) + key << '/'; + for (const auto& matrix : descriptor.stage_matrices[axis]) + key << matrix.m00 << ',' << matrix.m01 << ',' << matrix.m10 << ',' << matrix.m11 << ';'; + } + } + return key.str(); +} + +cubutterflyDescriptor make_axis_descriptor(const cubutterflyDescriptor& descriptor, + std::uint32_t source_axis, std::size_t length, + std::size_t batch) { + auto axis = descriptor; + axis.rank = 1; + axis.extents = {length}; + axis.batch = batch; + axis.length_mode = CUBUTTERFLY_LENGTH_STANDARD; + axis.placement = CUBUTTERFLY_PLACEMENT_IN_PLACE; + axis.input_strides = {1}; + axis.output_strides = {1}; + axis.input_batch_stride = length; + axis.output_batch_stride = length; + if (descriptor.op == CUBUTTERFLY_OPERATOR_STRUCTURED_2X2) { + require(source_axis < descriptor.stage_matrices.size(), "structured axis is out of range"); + axis.stage_matrices[0] = descriptor.stage_matrices[source_axis]; + axis.stage_matrices[1].clear(); + } + return axis; +} + +std::string lookup_cache(const std::string& path, const std::string& key) { + if (path.empty()) + return {}; + std::ifstream input(path); + std::string line; + std::string match; + while (std::getline(input, line)) { + const auto separator = line.find('\t'); + if (separator != std::string::npos && line.compare(0, separator, key) == 0 && separator == key.size()) + match = line.substr(separator + 1); + } + return match; +} + +void store_cache(const std::string& path, const std::string& key, const std::string& algorithm) { + if (path.empty()) + return; + std::ofstream output(path, std::ios::app); + if (!output) + throw std::runtime_error("could not open the cuButterfly selection cache for writing"); + output << key << '\t' << algorithm << '\n'; +} + +void apply_butterfly_algorithm(cuntt::ButterflyConfig& config, const std::string& algorithm) { + if (algorithm == "temporal-radix2") { + config.backend = cuntt::ButterflyBackend::TemporalTile; + config.compute_unit = cuntt::ComputeUnit::Radix2; + config.tile_threads = 128; + } else if (algorithm == "temporal-radix4") { + config.backend = cuntt::ButterflyBackend::TemporalTile; + config.compute_unit = cuntt::ComputeUnit::Radix4; + config.tile_threads = 128; + } else if (algorithm == "temporal-warp-register-radix2") { + config.backend = cuntt::ButterflyBackend::TemporalTile; + config.compute_unit = cuntt::ComputeUnit::Radix2; + config.local_exchange = cuntt::LocalExchange::WarpRegister; + config.tile_threads = 256; + } else if (algorithm == "hierarchical-radix4") { + config.backend = cuntt::ButterflyBackend::Hierarchical; + config.compute_unit = cuntt::ComputeUnit::Radix4; + config.local_stages = std::min(10, config.log_n - 1); + config.tile_threads = 256; + } else if (algorithm == "online-radix4") { + config.backend = cuntt::ButterflyBackend::OnlineReorder; + config.compute_unit = cuntt::ComputeUnit::Radix4; + config.local_stages = std::min(10, config.log_n - 1); + config.reorder_columns = 1; + config.tile_threads = 256; + } else if (algorithm == "cufft") { + require(config.op == cuntt::ButterflyOperator::Fft, "cuFFT algorithm is valid only for FFT"); + config.backend = cuntt::ButterflyBackend::CuFft; + } else { + throw std::invalid_argument("unknown explicit butterfly algorithm"); + } +} + +std::string butterfly_algorithm_id(const cuntt::ButterflyConfig& config) { + if (config.backend == cuntt::ButterflyBackend::TemporalTile && + config.local_exchange == cuntt::LocalExchange::WarpRegister) + return "temporal-warp-register-radix2"; + if (config.backend == cuntt::ButterflyBackend::TemporalTile) + return config.compute_unit == cuntt::ComputeUnit::Radix2 ? "temporal-radix2" : "temporal-radix4"; + if (config.backend == cuntt::ButterflyBackend::Hierarchical) + return "hierarchical-radix4"; + if (config.backend == cuntt::ButterflyBackend::OnlineReorder) + return "online-radix4"; + if (config.backend == cuntt::ButterflyBackend::CuFft) + return "cufft"; + throw std::invalid_argument("butterfly configuration has no public algorithm ID"); +} + +std::string ntt_algorithm_id(const cuntt::PlanConfig& config) { + if (config.backend == cuntt::Backend::Baseline) + return "ntt-baseline-radix2"; + if (config.backend == cuntt::Backend::Tile256) + return "ntt-tile256-radix2"; + if (config.backend == cuntt::Backend::Hybrid2D) + return config.compute_unit == cuntt::ComputeUnit::Radix4 ? "ntt-hybrid-radix4" + : "ntt-hybrid-radix2"; + throw std::invalid_argument("NTT configuration has no public algorithm ID"); +} + +void apply_ntt_algorithm(cuntt::PlanConfig& config, const std::string& algorithm) { + if (algorithm == "ntt-baseline-radix2") { + config.backend = cuntt::Backend::Baseline; + config.compute_unit = cuntt::ComputeUnit::Radix2; + } else if (algorithm == "ntt-tile256-radix2") { + config.backend = cuntt::Backend::Tile256; + config.compute_unit = cuntt::ComputeUnit::Radix2; + } else if (algorithm == "ntt-hybrid-radix2") { + config.backend = cuntt::Backend::Hybrid2D; + config.compute_unit = cuntt::ComputeUnit::Radix2; + } else if (algorithm == "ntt-hybrid-radix4") { + config.backend = cuntt::Backend::Hybrid2D; + config.compute_unit = cuntt::ComputeUnit::Radix4; + } else { + throw std::invalid_argument("unknown explicit NTT algorithm"); + } +} + +struct DeviceAllocation { + void* pointer = nullptr; + explicit DeviceAllocation(std::size_t bytes) { + if (bytes != 0 && cudaMalloc(&pointer, bytes) != cudaSuccess) + throw std::bad_alloc(); + } + ~DeviceAllocation() { cudaFree(pointer); } +}; + +template +double measure_cuda_candidate(cudaStream_t stream, Execute&& execute) { + execute(); + execute(); + cudaEvent_t start = nullptr; + cudaEvent_t stop = nullptr; + if (cudaEventCreate(&start) != cudaSuccess || cudaEventCreate(&stop) != cudaSuccess) + throw std::runtime_error("could not create candidate measurement events"); + cudaEventRecord(start, stream); + for (int repeat = 0; repeat < 5; ++repeat) + execute(); + cudaEventRecord(stop, stream); + cudaEventSynchronize(stop); + float milliseconds = 0.0F; + cudaEventElapsedTime(&milliseconds, start, stop); + cudaEventDestroy(stop); + cudaEventDestroy(start); + return milliseconds / 5.0; +} + +double measure_direct_fft_candidate(const cubutterflyDescriptor& descriptor, + std::size_t bytes, cudaStream_t stream) { + cuntt::detail::DirectFftPlan candidate( + descriptor.extents, descriptor.batch, + descriptor.direction == CUBUTTERFLY_DIRECTION_INVERSE, + descriptor.normalize_inverse, + descriptor.storage == CUBUTTERFLY_DATA_COMPLEX_FP64); + DeviceAllocation input(bytes); + DeviceAllocation output(bytes); + cudaMemsetAsync(input.pointer, 0, bytes, stream); + return measure_cuda_candidate(stream, [&] { candidate.execute(input.pointer, output.pointer, stream); }); +} + +double measure_bluestein_fft_candidate(const cubutterflyDescriptor& descriptor, + const std::vector& physical_extents, + std::size_t element_size, cudaStream_t stream) { + const bool inverse = descriptor.direction == CUBUTTERFLY_DIRECTION_INVERSE; + const bool fp64 = descriptor.storage == CUBUTTERFLY_DATA_COMPLEX_FP64; + std::size_t logical_points = descriptor.batch; + for (const auto extent : descriptor.extents) + logical_points *= extent; + const std::size_t logical_bytes = logical_points * element_size; + DeviceAllocation input(logical_bytes); + DeviceAllocation output(logical_bytes); + cudaMemsetAsync(input.pointer, 0, logical_bytes, stream); + if (descriptor.rank == 1) { + cuntt::detail::BluesteinFftPlan candidate( + descriptor.extents[0], physical_extents[0], descriptor.batch, + inverse, descriptor.normalize_inverse, fp64); + DeviceAllocation workspace(candidate.workspace_size()); + return measure_cuda_candidate(stream, [&] { + candidate.execute(input.pointer, output.pointer, workspace.pointer, 1, + descriptor.extents[0], 1, descriptor.extents[0], stream); + }); + } + + cuntt::detail::BluesteinFftPlan rows( + descriptor.extents[1], physical_extents[1], + descriptor.batch * descriptor.extents[0], inverse, + descriptor.normalize_inverse, fp64); + cuntt::detail::BluesteinFftPlan columns( + descriptor.extents[0], physical_extents[0], + descriptor.batch * descriptor.extents[1], inverse, + descriptor.normalize_inverse, fp64); + const std::size_t packed_bytes = logical_bytes; + const std::size_t axis_workspace = std::max(rows.workspace_size(), columns.workspace_size()); + DeviceAllocation workspace(2 * align_16(packed_bytes) + align_16(axis_workspace)); + auto* first = static_cast(workspace.pointer); + auto* second = first + align_16(packed_bytes); + auto* core_workspace = second + align_16(packed_bytes); + return measure_cuda_candidate(stream, [&] { + cuntt::detail::launch_pack_2d( + input.pointer, first, descriptor.extents[0], descriptor.extents[1], + descriptor.extents[0], descriptor.extents[1], descriptor.batch, + descriptor.input_strides[0], descriptor.input_strides[1], + descriptor.input_batch_stride, element_size, stream); + rows.execute(first, second, core_workspace, 1, descriptor.extents[1], 1, + descriptor.extents[1], stream); + cuntt::detail::launch_transpose_2d( + second, first, descriptor.extents[0], descriptor.extents[1], + descriptor.batch, element_size, stream); + columns.execute(first, second, core_workspace, 1, descriptor.extents[0], 1, + descriptor.extents[0], stream); + cuntt::detail::launch_transpose_2d( + second, first, descriptor.extents[1], descriptor.extents[0], + descriptor.batch, element_size, stream); + cuntt::detail::launch_scatter_2d( + first, output.pointer, descriptor.extents[0], descriptor.extents[1], + descriptor.extents[0], descriptor.extents[1], descriptor.batch, + descriptor.output_strides[0], descriptor.output_strides[1], + descriptor.output_batch_stride, element_size, stream); + }); +} + +double measure_butterfly(cuntt::ButterflyConfig config, cubutterflyDataType_t storage, + cudaStream_t stream) { + config.auto_allocate_workspace = false; + cuntt::ButterflyPlan plan(config); + plan.set_stream(stream); + DeviceAllocation input(plan.data_size()); + DeviceAllocation output(config.placement == cuntt::ButterflyPlacement::InPlace ? 0 : plan.data_size()); + DeviceAllocation workspace(plan.workspace_size()); + if (plan.workspace_size() != 0) + plan.set_workspace(workspace.pointer, plan.workspace_size()); + cudaMemsetAsync(input.pointer, 0, plan.data_size(), stream); + void* destination = output.pointer == nullptr ? input.pointer : output.pointer; + for (int warmup = 0; warmup < 2; ++warmup) + execute_butterfly(plan, input.pointer, destination, storage); + cudaEvent_t start = nullptr; + cudaEvent_t stop = nullptr; + if (cudaEventCreate(&start) != cudaSuccess || cudaEventCreate(&stop) != cudaSuccess) + throw std::runtime_error("could not create tuning events"); + cudaEventRecord(start, stream); + for (int repeat = 0; repeat < 5; ++repeat) + execute_butterfly(plan, input.pointer, destination, storage); + cudaEventRecord(stop, stream); + cudaEventSynchronize(stop); + float milliseconds = 0.0F; + cudaEventElapsedTime(&milliseconds, start, stop); + cudaEventDestroy(stop); + cudaEventDestroy(start); + return milliseconds / 5.0; +} + +std::pair tune_butterfly(const cuntt::ButterflyConfig& base, + cubutterflyDataType_t storage, + cudaStream_t stream) { + std::vector candidates{"temporal-radix2", "temporal-radix4", + "temporal-warp-register-radix2", + "hierarchical-radix4", "online-radix4"}; + if (base.op == cuntt::ButterflyOperator::Fft) + candidates.push_back("cufft"); + std::pair best{"", std::numeric_limits::infinity()}; + for (const auto& name : candidates) { + try { + auto candidate = base; + apply_butterfly_algorithm(candidate, name); + const double milliseconds = measure_butterfly(candidate, storage, stream); + if (milliseconds < best.second) + best = {name, milliseconds}; + } catch (const std::exception&) { + } + } + require(!best.first.empty(), "no legal butterfly candidate survived measurement"); + return best; +} + +double measure_ntt(cuntt::PlanConfig config, cudaStream_t stream) { + config.auto_allocate_workspace = false; + cuntt::Plan plan(config); + plan.set_stream(stream); + DeviceAllocation input(plan.data_size()); + DeviceAllocation output(plan.data_size()); + DeviceAllocation workspace(plan.workspace_size()); + if (plan.workspace_size() != 0) + plan.set_workspace(workspace.pointer, plan.workspace_size()); + cudaMemsetAsync(input.pointer, 0, plan.data_size(), stream); + auto execute = [&] { + if (config.word_bits == 32) + plan.execute_async(static_cast(input.pointer), static_cast(output.pointer)); + else + plan.execute_async(static_cast(input.pointer), static_cast(output.pointer)); + }; + execute(); + execute(); + cudaEvent_t start = nullptr; + cudaEvent_t stop = nullptr; + if (cudaEventCreate(&start) != cudaSuccess || cudaEventCreate(&stop) != cudaSuccess) + throw std::runtime_error("could not create tuning events"); + cudaEventRecord(start, stream); + for (int repeat = 0; repeat < 5; ++repeat) + execute(); + cudaEventRecord(stop, stream); + cudaEventSynchronize(stop); + float milliseconds = 0.0F; + cudaEventElapsedTime(&milliseconds, start, stop); + cudaEventDestroy(stop); + cudaEventDestroy(start); + return milliseconds / 5.0; +} + +std::pair tune_ntt(const cuntt::PlanConfig& base, cudaStream_t stream) { + const std::vector candidates{"ntt-baseline-radix2", "ntt-tile256-radix2", + "ntt-hybrid-radix2", "ntt-hybrid-radix4"}; + std::pair best{"", std::numeric_limits::infinity()}; + for (const auto& name : candidates) { + try { + auto candidate = base; + apply_ntt_algorithm(candidate, name); + const double milliseconds = measure_ntt(candidate, stream); + if (milliseconds < best.second) + best = {name, milliseconds}; + } catch (const std::exception&) { + } + } + require(!best.first.empty(), "no legal NTT candidate survived measurement"); + return best; +} + +struct ButterflyAxisSelection { + cuntt::ButterflyConfig config; + std::string algorithm; + cubutterflySelectionSource_t source = CUBUTTERFLY_SELECTION_STATIC_MODEL; + std::string reason; +}; + +ButterflyAxisSelection select_butterfly_axis(cubutterflyHandle_t handle, + const cubutterflyDescriptor& descriptor, + std::uint32_t source_axis, std::size_t length, + std::size_t batch) { + auto config = make_axis_butterfly_config(descriptor, source_axis, length, batch); + const auto axis_descriptor = make_axis_descriptor(descriptor, source_axis, length, batch); + const auto key = workload_key(axis_descriptor); + const auto cached = descriptor.algorithm_policy == CUBUTTERFLY_ALGORITHM_DEFAULT + ? lookup_cache(handle->cache_path, key) + : std::string{}; + ButterflyAxisSelection selected; + if (descriptor.algorithm_policy == CUBUTTERFLY_ALGORITHM_EXPLICIT) { + apply_butterfly_algorithm(config, descriptor.explicit_algorithm); + selected.source = CUBUTTERFLY_SELECTION_EXPLICIT; + selected.reason = "applied the explicit algorithm to the physical axis"; + } else if (!cached.empty()) { + apply_butterfly_algorithm(config, cached); + selected.source = CUBUTTERFLY_SELECTION_USER_CACHE; + selected.reason = "matched the physical axis in the user selection cache"; + } else if (descriptor.algorithm_policy == CUBUTTERFLY_ALGORITHM_MEASURE) { + const auto tuned = tune_butterfly(config, descriptor.storage, handle->stream); + apply_butterfly_algorithm(config, tuned.first); + store_cache(handle->cache_path, key, tuned.first); + selected.source = CUBUTTERFLY_SELECTION_RUNTIME_MEASURED; + selected.reason = "measured physical-axis candidates; kernel_ms=" + std::to_string(tuned.second); + } else { + const bool profiled = measured_v100_axis(descriptor, integer_log2(length), batch); + choose_butterfly_mapping(config, profiled); + selected.source = profiled ? CUBUTTERFLY_SELECTION_MEASURED_TABLE + : CUBUTTERFLY_SELECTION_STATIC_MODEL; + selected.reason = profiled ? "matched the built-in V100 physical-axis profile" + : "selected the physical axis with the static resource model"; + } + selected.algorithm = butterfly_algorithm_id(config); + selected.config = std::move(config); + return selected; +} + +struct NttAxisSelection { + cuntt::PlanConfig config; + std::string algorithm; + cubutterflySelectionSource_t source = CUBUTTERFLY_SELECTION_STATIC_MODEL; + std::string reason; +}; + +NttAxisSelection select_ntt_axis(cubutterflyHandle_t handle, + const cubutterflyDescriptor& descriptor, + std::uint32_t source_axis, std::size_t length, + std::size_t batch) { + auto config = make_axis_ntt_config(descriptor, length, batch); + const auto axis_descriptor = make_axis_descriptor(descriptor, source_axis, length, batch); + const auto key = workload_key(axis_descriptor); + const auto cached = descriptor.algorithm_policy == CUBUTTERFLY_ALGORITHM_DEFAULT + ? lookup_cache(handle->cache_path, key) + : std::string{}; + NttAxisSelection selected; + if (descriptor.algorithm_policy == CUBUTTERFLY_ALGORITHM_EXPLICIT) { + apply_ntt_algorithm(config, descriptor.explicit_algorithm); + selected.source = CUBUTTERFLY_SELECTION_EXPLICIT; + selected.reason = "applied the explicit algorithm to the physical NTT axis"; + } else if (!cached.empty()) { + apply_ntt_algorithm(config, cached); + selected.source = CUBUTTERFLY_SELECTION_USER_CACHE; + selected.reason = "matched the physical NTT axis in the user selection cache"; + } else if (descriptor.algorithm_policy == CUBUTTERFLY_ALGORITHM_MEASURE) { + const auto tuned = tune_ntt(config, handle->stream); + apply_ntt_algorithm(config, tuned.first); + store_cache(handle->cache_path, key, tuned.first); + selected.source = CUBUTTERFLY_SELECTION_RUNTIME_MEASURED; + selected.reason = "measured physical NTT-axis candidates; kernel_ms=" + std::to_string(tuned.second); + } else { + const bool profiled = measured_v100_axis(descriptor, integer_log2(length), batch); + selected.source = profiled ? CUBUTTERFLY_SELECTION_MEASURED_TABLE + : CUBUTTERFLY_SELECTION_STATIC_MODEL; + selected.reason = profiled ? "matched the built-in V100 physical NTT-axis profile" + : "selected the physical NTT axis with the static resource model"; + } + selected.algorithm = ntt_algorithm_id(config); + selected.config = std::move(config); + return selected; +} + +cubutterflySelectionSource_t combined_source( + cubutterflySelectionSource_t current, cubutterflySelectionSource_t next) { + if (current == CUBUTTERFLY_SELECTION_RUNTIME_MEASURED || next == CUBUTTERFLY_SELECTION_RUNTIME_MEASURED) + return CUBUTTERFLY_SELECTION_RUNTIME_MEASURED; + if (current == CUBUTTERFLY_SELECTION_USER_CACHE || next == CUBUTTERFLY_SELECTION_USER_CACHE) + return CUBUTTERFLY_SELECTION_USER_CACHE; + if (current == CUBUTTERFLY_SELECTION_EXPLICIT || next == CUBUTTERFLY_SELECTION_EXPLICIT) + return CUBUTTERFLY_SELECTION_EXPLICIT; + if (current == CUBUTTERFLY_SELECTION_MEASURED_TABLE || next == CUBUTTERFLY_SELECTION_MEASURED_TABLE) + return CUBUTTERFLY_SELECTION_MEASURED_TABLE; + return CUBUTTERFLY_SELECTION_STATIC_MODEL; +} + +} // namespace + +const char* cubutterflyGetStatusString(cubutterflyStatus_t status) { + switch (status) { + case CUBUTTERFLY_STATUS_SUCCESS: return "success"; + case CUBUTTERFLY_STATUS_INVALID_VALUE: return "invalid value"; + case CUBUTTERFLY_STATUS_INVALID_TYPE: return "invalid type"; + case CUBUTTERFLY_STATUS_NOT_SUPPORTED: return "not supported"; + case CUBUTTERFLY_STATUS_ALLOC_FAILED: return "allocation failed"; + case CUBUTTERFLY_STATUS_CUDA_ERROR: return "CUDA error"; + case CUBUTTERFLY_STATUS_BACKEND_ERROR: return "backend error"; + case CUBUTTERFLY_STATUS_INTERNAL_ERROR: return "internal error"; + } + return "unknown status"; +} + +cubutterflyStatus_t cubutterflyGetVersion(int* major, int* minor, int* patch) { + return protect([&] { + require(major != nullptr && minor != nullptr && patch != nullptr, "version pointers must be non-null"); + *major = CUBUTTERFLY_VERSION_MAJOR; + *minor = CUBUTTERFLY_VERSION_MINOR; + *patch = CUBUTTERFLY_VERSION_PATCH; + }); +} + +cubutterflyStatus_t cubutterflyCreate(cubutterflyHandle_t* handle) { + return protect([&] { require(handle != nullptr, "handle pointer must be non-null"); *handle = new cubutterflyHandle; }); +} + +cubutterflyStatus_t cubutterflyDestroy(cubutterflyHandle_t handle) { delete handle; return CUBUTTERFLY_STATUS_SUCCESS; } + +cubutterflyStatus_t cubutterflySetStream(cubutterflyHandle_t handle, cudaStream_t stream) { + return protect([&] { require(handle != nullptr, "handle must be non-null"); handle->stream = stream; }); +} + +cubutterflyStatus_t cubutterflyGetStream(cubutterflyHandle_t handle, cudaStream_t* stream) { + return protect([&] { require(handle != nullptr && stream != nullptr, "handle and stream pointer must be non-null"); *stream = handle->stream; }); +} + +cubutterflyStatus_t cubutterflySetLogCallback(cubutterflyHandle_t handle, cubutterflyLogCallback_t callback, void* user_data) { + return protect([&] { require(handle != nullptr, "handle must be non-null"); handle->logger = callback; handle->logger_data = user_data; }); +} + +cubutterflyStatus_t cubutterflySetCachePath(cubutterflyHandle_t handle, const char* path) { + return protect([&] { require(handle != nullptr, "handle must be non-null"); handle->cache_path = path == nullptr ? "" : path; }); +} + +cubutterflyStatus_t cubutterflyCreateDescriptor(cubutterflyDescriptor_t* descriptor) { + return protect([&] { require(descriptor != nullptr, "descriptor pointer must be non-null"); *descriptor = new cubutterflyDescriptor; }); +} + +cubutterflyStatus_t cubutterflyDestroyDescriptor(cubutterflyDescriptor_t descriptor) { delete descriptor; return CUBUTTERFLY_STATUS_SUCCESS; } + +cubutterflyStatus_t cubutterflySetOperator(cubutterflyDescriptor_t descriptor, cubutterflyOperator_t op) { + return protect([&] { require(descriptor != nullptr, "descriptor must be non-null"); descriptor->op = op; }); +} + +cubutterflyStatus_t cubutterflySetDataType(cubutterflyDescriptor_t descriptor, cubutterflyDataType_t storage, cubutterflyComputeType_t compute) { + return protect([&] { require(descriptor != nullptr, "descriptor must be non-null"); descriptor->storage = storage; descriptor->compute = compute; }); +} + +cubutterflyStatus_t cubutterflySetShape(cubutterflyDescriptor_t descriptor, uint32_t rank, const size_t* extents, size_t batch) { + return protect([&] { + require(descriptor != nullptr && extents != nullptr, "descriptor and extents must be non-null"); + require(rank == 1 || rank == 2, "rank must be one or two"); + require(batch != 0, "batch must be positive"); + descriptor->rank = rank; + descriptor->extents.assign(extents, extents + rank); + require(std::all_of(descriptor->extents.begin(), descriptor->extents.end(), [](std::size_t n) { return n != 0; }), + "extents must be positive"); + descriptor->batch = batch; + descriptor->input_strides.clear(); + descriptor->output_strides.clear(); + descriptor->input_batch_stride = 0; + descriptor->output_batch_stride = 0; + }); +} + +cubutterflyStatus_t cubutterflySetStrides(cubutterflyDescriptor_t descriptor, const size_t* input_strides, size_t input_batch_stride, + const size_t* output_strides, size_t output_batch_stride) { + return protect([&] { + require(descriptor != nullptr && input_strides != nullptr && output_strides != nullptr, "descriptor and strides must be non-null"); + descriptor->input_strides.assign(input_strides, input_strides + descriptor->rank); + descriptor->output_strides.assign(output_strides, output_strides + descriptor->rank); + descriptor->input_batch_stride = input_batch_stride; + descriptor->output_batch_stride = output_batch_stride; + }); +} + +cubutterflyStatus_t cubutterflySetDirection(cubutterflyDescriptor_t descriptor, cubutterflyDirection_t direction, int normalize_inverse) { + return protect([&] { require(descriptor != nullptr, "descriptor must be non-null"); descriptor->direction = direction; descriptor->normalize_inverse = normalize_inverse != 0; }); +} + +cubutterflyStatus_t cubutterflySetPlacement(cubutterflyDescriptor_t descriptor, cubutterflyPlacement_t placement) { + return protect([&] { require(descriptor != nullptr, "descriptor must be non-null"); descriptor->placement = placement; }); +} + +cubutterflyStatus_t cubutterflySetLengthMode(cubutterflyDescriptor_t descriptor, cubutterflyLengthMode_t mode) { + return protect([&] { require(descriptor != nullptr, "descriptor must be non-null"); descriptor->length_mode = mode; }); +} + +cubutterflyStatus_t cubutterflySetAlgorithmPolicy(cubutterflyDescriptor_t descriptor, cubutterflyAlgorithmPolicy_t policy, const char* explicit_algorithm) { + return protect([&] { + require(descriptor != nullptr, "descriptor must be non-null"); + require(policy != CUBUTTERFLY_ALGORITHM_EXPLICIT || (explicit_algorithm != nullptr && explicit_algorithm[0] != '\0'), + "explicit policy requires an algorithm name"); + descriptor->algorithm_policy = policy; + descriptor->explicit_algorithm = explicit_algorithm == nullptr ? "" : explicit_algorithm; + }); +} + +cubutterflyStatus_t cubutterflySetModulus(cubutterflyDescriptor_t descriptor, uint64_t modulus, uint32_t word_bits) { + return protect([&] { require(descriptor != nullptr, "descriptor must be non-null"); require(word_bits == 32 || word_bits == 64, "word width must be 32 or 64"); descriptor->modulus = modulus; descriptor->word_bits = word_bits; }); +} + +cubutterflyStatus_t cubutterflySetStageMatrices(cubutterflyDescriptor_t descriptor, uint32_t axis, + const cubutterflyMatrix2x2_t* matrices, size_t count) { + return protect([&] { + require(descriptor != nullptr && matrices != nullptr && count != 0, "matrix arguments must be non-null and non-empty"); + require(axis < 2, "matrix axis must be zero or one"); + descriptor->stage_matrices[axis].assign(matrices, matrices + count); + }); +} + +cubutterflyStatus_t cubutterflyCreatePlan(cubutterflyHandle_t handle, cubutterflyDescriptor_t descriptor, cubutterflyPlan_t* plan) { + return protect([&] { + require(handle != nullptr && descriptor != nullptr && plan != nullptr, "handle, descriptor, and plan pointer must be non-null"); + auto resolved = *descriptor; + validate_storage(resolved); + std::vector physical_extents = resolved.extents; + bool standard_arbitrary = false; + if (resolved.length_mode == CUBUTTERFLY_LENGTH_ZERO_EXTENDED_EMBEDDING) { + require(resolved.placement == CUBUTTERFLY_PLACEMENT_OUT_OF_PLACE, + "zero-extended embedding currently requires out-of-place execution"); + for (auto& extent : physical_extents) + extent = next_power_of_two(extent); + } else { + standard_arbitrary = std::any_of(physical_extents.begin(), physical_extents.end(), + [](std::size_t extent) { return !is_power_of_two(extent); }); + if (standard_arbitrary) { + require(resolved.op == CUBUTTERFLY_OPERATOR_FFT || resolved.op == CUBUTTERFLY_OPERATOR_NTT, + "standard arbitrary length supports FFT and NTT; use embedding mode for other operators"); + if (resolved.op == CUBUTTERFLY_OPERATOR_FFT) + require(resolved.storage == CUBUTTERFLY_DATA_COMPLEX_FP32 || + resolved.storage == CUBUTTERFLY_DATA_COMPLEX_FP64, + "Bluestein FFT currently supports complex FP32 and FP64 storage"); + else + require(resolved.storage == CUBUTTERFLY_DATA_UINT64 && resolved.word_bits == 64, + "Bluestein NTT currently supports the uint64 storage contract"); + for (std::uint32_t axis = 0; axis < resolved.rank; ++axis) { + require(resolved.extents[axis] <= static_cast(-1) / 2 + 1, + "Bluestein convolution length overflows size_t"); + physical_extents[axis] = next_power_of_two(2 * resolved.extents[axis] - 1); + } + } + } + const bool embedding_inverse = resolved.length_mode == CUBUTTERFLY_LENGTH_ZERO_EXTENDED_EMBEDDING && + resolved.direction == CUBUTTERFLY_DIRECTION_INVERSE; + const auto& input_extents = embedding_inverse ? physical_extents : resolved.extents; + const auto& output_extents = resolved.length_mode == CUBUTTERFLY_LENGTH_ZERO_EXTENDED_EMBEDDING && !embedding_inverse + ? physical_extents : resolved.extents; + resolve_default_strides(resolved, input_extents, output_extents); + const std::uint32_t log_n = resolved.rank == 1 ? integer_log2(physical_extents[0]) : 0; + auto created = std::make_unique(); + created->descriptor = resolved; + created->physical_extents = physical_extents; + created->element_bytes = element_bytes(resolved.storage); + created->input_bytes = strided_elements(resolved, true, input_extents) * created->element_bytes; + created->output_bytes = strided_elements(resolved, false, output_extents) * created->element_bytes; + const bool measured = !standard_arbitrary && resolved.rank == 1 && resolved.length_mode == CUBUTTERFLY_LENGTH_STANDARD && + resolved.algorithm_policy == CUBUTTERFLY_ALGORITHM_DEFAULT && measured_v100_shape(resolved, log_n); + created->selection_source = resolved.algorithm_policy == CUBUTTERFLY_ALGORITHM_EXPLICIT + ? CUBUTTERFLY_SELECTION_EXPLICIT + : (measured ? CUBUTTERFLY_SELECTION_MEASURED_TABLE : CUBUTTERFLY_SELECTION_STATIC_MODEL); + created->composite = resolved.rank == 2 || resolved.length_mode == CUBUTTERFLY_LENGTH_ZERO_EXTENDED_EMBEDDING; + created->direct_output_boundary = created->composite && output_extents == physical_extents && + packed_side(resolved, false, output_extents); + const std::string cache_key = workload_key(resolved); + const std::string cached_algorithm = resolved.algorithm_policy == CUBUTTERFLY_ALGORITHM_DEFAULT + ? lookup_cache(handle->cache_path, cache_key) + : std::string{}; + if (standard_arbitrary) { + require(resolved.placement == CUBUTTERFLY_PLACEMENT_OUT_OF_PLACE, + "Bluestein FFT currently requires out-of-place execution"); + const bool ntt_bluestein = resolved.op == CUBUTTERFLY_OPERATOR_NTT; + const char* base_algorithm = ntt_bluestein ? "bluestein-ntt-power2-core" : "bluestein-cufft-power2-core"; + bool direct_fft = false; + std::string arbitrary_selection_reason; + cubutterflySelectionSource_t arbitrary_selection_source = CUBUTTERFLY_SELECTION_STATIC_MODEL; + if (!ntt_bluestein) { + if (resolved.algorithm_policy == CUBUTTERFLY_ALGORITHM_EXPLICIT) { + direct_fft = resolved.explicit_algorithm == "cufft-direct"; + require(!direct_fft || packed_layout(resolved), + "cufft-direct requires packed contiguous batches"); + arbitrary_selection_source = CUBUTTERFLY_SELECTION_EXPLICIT; + } else if (!cached_algorithm.empty()) { + require(cached_algorithm == "cufft-direct" || cached_algorithm == base_algorithm, + "cached arbitrary FFT algorithm is not recognized"); + direct_fft = cached_algorithm == "cufft-direct"; + require(!direct_fft || packed_layout(resolved), + "cached cufft-direct selection requires packed contiguous batches"); + arbitrary_selection_source = CUBUTTERFLY_SELECTION_USER_CACHE; + arbitrary_selection_reason = "matched the exact arbitrary FFT workload in the user cache"; + } else if (resolved.algorithm_policy == CUBUTTERFLY_ALGORITHM_MEASURE) { + const double bluestein_ms = measure_bluestein_fft_candidate( + resolved, physical_extents, created->element_bytes, handle->stream); + const double direct_ms = packed_layout(resolved) + ? measure_direct_fft_candidate( + resolved, created->input_bytes, handle->stream) + : std::numeric_limits::infinity(); + direct_fft = direct_ms <= bluestein_ms; + const std::string selected_id = direct_fft ? "cufft-direct" : base_algorithm; + store_cache(handle->cache_path, cache_key, selected_id); + arbitrary_selection_source = CUBUTTERFLY_SELECTION_RUNTIME_MEASURED; + arbitrary_selection_reason = "measured exact-length candidates; direct_ms=" + + std::to_string(direct_ms) + "; bluestein_ms=" + std::to_string(bluestein_ms); + log_event(handle, CUBUTTERFLY_LOG_INFO, "measure-complete", + arbitrary_selection_reason + "; algorithm=" + selected_id); + } else { + direct_fft = packed_layout(resolved); + arbitrary_selection_reason = "static policy prefers the vendor native exact-length candidate"; + } + if (direct_fft) { + created->kind = cubutterflyPlan::Kind::DirectFft; + created->direct_fft = std::make_unique( + resolved.extents, resolved.batch, + resolved.direction == CUBUTTERFLY_DIRECTION_INVERSE, + resolved.normalize_inverse, + resolved.storage == CUBUTTERFLY_DATA_COMPLEX_FP64); + created->physical_extents = resolved.extents; + created->workspace_bytes = 0; + created->algorithm = "cufft-direct"; + created->reason = arbitrary_selection_reason.empty() + ? "selected the vendor native exact-length FFT candidate; Bluestein remains available explicitly" + : arbitrary_selection_reason; + created->selection_source = arbitrary_selection_source; + log_event(handle, CUBUTTERFLY_LOG_INFO, "algorithm-selection", + created->reason + "; algorithm=" + created->algorithm); + *plan = created.release(); + return; + } + } + auto core_descriptor = resolved; + if (resolved.algorithm_policy == CUBUTTERFLY_ALGORITHM_EXPLICIT) { + if (resolved.explicit_algorithm == base_algorithm) { + core_descriptor.algorithm_policy = CUBUTTERFLY_ALGORITHM_DEFAULT; + core_descriptor.explicit_algorithm.clear(); + } else if (ntt_bluestein && resolved.explicit_algorithm.rfind("bluestein-", 0) == 0) { + core_descriptor.explicit_algorithm = resolved.explicit_algorithm.substr(10); + } else { + throw std::invalid_argument("unknown explicit arbitrary-length algorithm"); + } + } + std::vector ntt_core_selections; + if (ntt_bluestein) { + created->selection_source = CUBUTTERFLY_SELECTION_STATIC_MODEL; + if (resolved.rank == 1) { + ntt_core_selections.push_back(select_ntt_axis( + handle, core_descriptor, 0, physical_extents[0], resolved.batch)); + } else { + ntt_core_selections.push_back(select_ntt_axis( + handle, core_descriptor, 1, physical_extents[1], + resolved.batch * resolved.extents[0])); + ntt_core_selections.push_back(select_ntt_axis( + handle, core_descriptor, 0, physical_extents[0], + resolved.batch * resolved.extents[1])); + } + for (const auto& selection : ntt_core_selections) + created->selection_source = combined_source(created->selection_source, selection.source); + } + created->kind = ntt_bluestein ? cubutterflyPlan::Kind::BluesteinNtt + : cubutterflyPlan::Kind::BluesteinFft; + created->algorithm = base_algorithm; + if (ntt_bluestein) { + std::ostringstream algorithm; + algorithm << "bluestein-ntt["; + for (std::size_t axis = 0; axis < ntt_core_selections.size(); ++axis) + algorithm << (axis == 0 ? "" : ",") << ntt_core_selections[axis].algorithm; + algorithm << ']'; + created->algorithm = algorithm.str(); + } + created->reason = ntt_bluestein + ? "non-power-of-two standard NTT selected Bluestein after validating target and convolution roots" + : "non-power-of-two standard FFT selected Bluestein with a power-of-two convolution core"; + if (!ntt_bluestein && !arbitrary_selection_reason.empty()) + created->reason += "; " + arbitrary_selection_reason; + if (!ntt_bluestein) + created->selection_source = arbitrary_selection_source; + const bool inverse = resolved.direction == CUBUTTERFLY_DIRECTION_INVERSE; + const bool fp64 = resolved.storage == CUBUTTERFLY_DATA_COMPLEX_FP64; + if (resolved.rank == 1) { + if (ntt_bluestein) { + created->bluestein_ntt_axes.push_back(std::make_unique( + resolved.extents[0], physical_extents[0], resolved.batch, resolved.modulus, inverse, + ntt_core_selections[0].config.backend, + ntt_core_selections[0].config.compute_unit)); + created->workspace_bytes = created->bluestein_ntt_axes[0]->workspace_size(); + } else { + created->bluestein_axes.push_back(std::make_unique( + resolved.extents[0], physical_extents[0], resolved.batch, + inverse, resolved.normalize_inverse, fp64)); + created->workspace_bytes = created->bluestein_axes[0]->workspace_size(); + } + } else { + require(resolved.extents[0] <= static_cast(-1) / resolved.extents[1], + "rank-2 point count overflows size_t"); + const std::size_t matrix_points = resolved.extents[0] * resolved.extents[1]; + require(resolved.batch <= static_cast(-1) / matrix_points, + "rank-2 batch point count overflows size_t"); + const std::size_t points = matrix_points * resolved.batch; + require(points <= static_cast(-1) / created->element_bytes, + "rank-2 packed byte count overflows size_t"); + created->packed_bytes = points * created->element_bytes; + if (ntt_bluestein) { + created->bluestein_ntt_axes.push_back(std::make_unique( + resolved.extents[1], physical_extents[1], resolved.batch * resolved.extents[0], + resolved.modulus, inverse, ntt_core_selections[0].config.backend, + ntt_core_selections[0].config.compute_unit)); + created->bluestein_ntt_axes.push_back(std::make_unique( + resolved.extents[0], physical_extents[0], resolved.batch * resolved.extents[1], + resolved.modulus, inverse, ntt_core_selections[1].config.backend, + ntt_core_selections[1].config.compute_unit)); + } else { + created->bluestein_axes.push_back(std::make_unique( + resolved.extents[1], physical_extents[1], resolved.batch * resolved.extents[0], + inverse, resolved.normalize_inverse, fp64)); + created->bluestein_axes.push_back(std::make_unique( + resolved.extents[0], physical_extents[0], resolved.batch * resolved.extents[1], + inverse, resolved.normalize_inverse, fp64)); + } + const auto first_work = ntt_bluestein ? created->bluestein_ntt_axes[0]->workspace_size() + : created->bluestein_axes[0]->workspace_size(); + const auto second_work = ntt_bluestein ? created->bluestein_ntt_axes[1]->workspace_size() + : created->bluestein_axes[1]->workspace_size(); + created->workspace_bytes = 2 * align_16(created->packed_bytes) + + std::max(align_16(first_work), align_16(second_work)); + if (ntt_bluestein) + created->algorithm = "rank2-transpose-" + created->algorithm; + else + created->algorithm = "rank2-transpose-bluestein-cufft-power2-core"; + } + log_event(handle, CUBUTTERFLY_LOG_INFO, "algorithm-selection", + created->reason + "; algorithm=" + created->algorithm); + if (resolved.algorithm_policy == CUBUTTERFLY_ALGORITHM_DEFAULT) + log_profile_miss(handle, resolved, created->reason, created->algorithm); + *plan = created.release(); + return; + } + if (created->composite) { + std::size_t packed_points = resolved.batch; + for (const auto extent : physical_extents) { + require(packed_points <= static_cast(-1) / extent, "packed point count overflows size_t"); + packed_points *= extent; + } + require(packed_points <= static_cast(-1) / created->element_bytes, "packed byte count overflows size_t"); + created->packed_bytes = packed_points * created->element_bytes; + created->kind = resolved.op == CUBUTTERFLY_OPERATOR_NTT ? cubutterflyPlan::Kind::Ntt : cubutterflyPlan::Kind::Butterfly; + std::vector axis_algorithms; + std::vector axis_reasons; + created->selection_source = CUBUTTERFLY_SELECTION_STATIC_MODEL; + auto add_ntt_axis = [&](std::uint32_t source_axis, std::size_t length, std::size_t batch) { + auto selected = select_ntt_axis(handle, resolved, source_axis, length, batch); + created->selection_source = combined_source(created->selection_source, selected.source); + axis_algorithms.push_back(selected.algorithm); + axis_reasons.push_back(selected.reason); + created->ntt_axes.push_back(std::make_unique(selected.config)); + created->ntt_axes.back()->set_stream(handle->stream); + }; + auto add_butterfly_axis = [&](std::uint32_t source_axis, std::size_t length, std::size_t batch) { + auto selected = select_butterfly_axis(handle, resolved, source_axis, length, batch); + created->selection_source = combined_source(created->selection_source, selected.source); + axis_algorithms.push_back(selected.algorithm); + axis_reasons.push_back(selected.reason); + created->butterfly_axes.push_back(std::make_unique(selected.config)); + created->butterfly_axes.back()->set_stream(handle->stream); + }; + if (resolved.rank == 1) { + if (created->kind == cubutterflyPlan::Kind::Ntt) + add_ntt_axis(0, physical_extents[0], resolved.batch); + else + add_butterfly_axis(0, physical_extents[0], resolved.batch); + } else { + const std::size_t row_transforms = resolved.batch * physical_extents[0]; + const std::size_t column_transforms = resolved.batch * physical_extents[1]; + if (created->kind == cubutterflyPlan::Kind::Ntt) { + add_ntt_axis(1, physical_extents[1], row_transforms); + add_ntt_axis(0, physical_extents[0], column_transforms); + } else { + add_butterfly_axis(1, physical_extents[1], row_transforms); + add_butterfly_axis(0, physical_extents[0], column_transforms); + } + } + created->fused_input_boundary = resolved.rank == 1 && + created->direct_output_boundary && + resolved.length_mode == CUBUTTERFLY_LENGTH_ZERO_EXTENDED_EMBEDDING && + resolved.direction == CUBUTTERFLY_DIRECTION_FORWARD && + created->kind == cubutterflyPlan::Kind::Butterfly && + resolved.storage == CUBUTTERFLY_DATA_FP32 && + resolved.op == CUBUTTERFLY_OPERATOR_FWHT && + created->butterfly_axes[0]->config().local_exchange == + cuntt::LocalExchange::WarpRegister; + created->composite_buffers = resolved.rank == 2 + ? 2 + : (created->kind == cubutterflyPlan::Kind::Ntt + ? (created->direct_output_boundary ? 1 : 2) + : (created->direct_output_boundary ? 0 : 1)); + created->workspace_bytes = align_16(created->packed_bytes) * created->composite_buffers; + if (created->kind == cubutterflyPlan::Kind::Ntt) { + for (const auto& axis_plan : created->ntt_axes) + created->workspace_bytes += align_16(axis_plan->workspace_size()); + } else { + for (const auto& axis_plan : created->butterfly_axes) + created->workspace_bytes += align_16(axis_plan->workspace_size()); + } + std::ostringstream algorithm; + algorithm << (resolved.rank == 2 ? "rank2-transpose" : "zero-extended") << '['; + for (std::size_t axis = 0; axis < axis_algorithms.size(); ++axis) + algorithm << (axis == 0 ? "" : ",") << axis_algorithms[axis]; + algorithm << ']'; + if (created->direct_output_boundary) + algorithm << "+direct-output"; + if (created->fused_input_boundary) + algorithm << "+fused-input"; + created->algorithm = algorithm.str(); + std::ostringstream reason; + for (std::size_t axis = 0; axis < axis_reasons.size(); ++axis) + reason << (axis == 0 ? "" : "; ") << "axis" << axis << ": " << axis_reasons[axis]; + created->reason = reason.str(); + log_event(handle, CUBUTTERFLY_LOG_INFO, "composite-selection", + created->reason + "; algorithm=" + created->algorithm); + if (created->selection_source == CUBUTTERFLY_SELECTION_STATIC_MODEL) + log_profile_miss(handle, resolved, created->reason, created->algorithm); + *plan = created.release(); + return; + } + if (resolved.op == CUBUTTERFLY_OPERATOR_NTT) { + require(resolved.input_strides[0] == 1 && resolved.output_strides[0] == 1 && + resolved.input_batch_stride == resolved.extents[0] && + resolved.output_batch_stride == resolved.extents[0], + "rank-1 NTT direct execution requires packed contiguous batches"); + cuntt::PlanConfig config; + config.log_n = log_n; + config.batch = resolved.batch; + config.modulus = resolved.modulus; + config.inverse = resolved.direction == CUBUTTERFLY_DIRECTION_INVERSE; + config.word_bits = resolved.word_bits; + config.auto_allocate_workspace = false; + config.auto_select = measured; + if (resolved.algorithm_policy == CUBUTTERFLY_ALGORITHM_EXPLICIT) { + apply_ntt_algorithm(config, resolved.explicit_algorithm); + } else if (!cached_algorithm.empty()) { + apply_ntt_algorithm(config, cached_algorithm); + created->selection_source = CUBUTTERFLY_SELECTION_USER_CACHE; + } else if (resolved.algorithm_policy == CUBUTTERFLY_ALGORITHM_MEASURE) { + const auto tuned = tune_ntt(config, handle->stream); + apply_ntt_algorithm(config, tuned.first); + created->selection_source = CUBUTTERFLY_SELECTION_RUNTIME_MEASURED; + created->reason = "measured legal NTT candidates on the current hardware; kernel_ms=" + + std::to_string(tuned.second); + store_cache(handle->cache_path, cache_key, tuned.first); + log_event(handle, CUBUTTERFLY_LOG_INFO, "measure-complete", + created->reason + "; algorithm=" + tuned.first); + } else if (log_n >= 12 && log_n <= 20) { + config.backend = cuntt::Backend::Hybrid2D; + config.compute_unit = cuntt::ComputeUnit::Radix4; + } else { + config.backend = cuntt::Backend::Baseline; + config.compute_unit = cuntt::ComputeUnit::Radix2; + } + created->kind = cubutterflyPlan::Kind::Ntt; + if (created->reason.empty()) + created->reason = created->selection_source == CUBUTTERFLY_SELECTION_USER_CACHE + ? "matched an exact workload key in the user selection cache" + : (measured ? "matched the built-in V100 measured profile" + : "no measured profile entry; selected by the static resource model"); + created->ntt = std::make_unique(config); + created->ntt->set_stream(handle->stream); + created->algorithm = std::string("ntt-") + cuntt::backend_name(created->ntt->config().backend) + + "-" + cuntt::compute_unit_name(created->ntt->config().compute_unit); + created->workspace_bytes = created->ntt->workspace_size(); + } else { + cuntt::ButterflyConfig config; + config.op = butterfly_operator(resolved.op); + config.precision = butterfly_precision(resolved); + const bool low_precision = config.precision == cuntt::ButterflyPrecision::Fp16 || + config.precision == cuntt::ButterflyPrecision::Bf16; + config.accumulation = low_precision && resolved.compute == CUBUTTERFLY_COMPUTE_FP32 + ? cuntt::ButterflyAccumulation::Fp32 + : cuntt::ButterflyAccumulation::Native; + config.placement = resolved.placement == CUBUTTERFLY_PLACEMENT_IN_PLACE ? cuntt::ButterflyPlacement::InPlace : cuntt::ButterflyPlacement::OutOfPlace; + config.log_n = log_n; + config.batch = resolved.batch; + config.element_stride = resolved.input_strides[0]; + config.batch_stride = resolved.input_batch_stride; + require(resolved.input_strides == resolved.output_strides && resolved.input_batch_stride == resolved.output_batch_stride, + "rank-1 common plans currently require matching input and output layouts"); + config.inverse = resolved.direction == CUBUTTERFLY_DIRECTION_INVERSE; + config.normalize_inverse = resolved.normalize_inverse; + config.auto_allocate_workspace = false; + config.auto_select = measured; + if (config.op == cuntt::ButterflyOperator::Structured2x2) { + require(!resolved.stage_matrices[0].empty(), "structured-2x2 requires stage matrices"); + for (const auto& matrix : resolved.stage_matrices[0]) + config.stage_matrices.push_back({matrix.m00, matrix.m01, matrix.m10, matrix.m11}); + } + if (resolved.algorithm_policy == CUBUTTERFLY_ALGORITHM_EXPLICIT) { + apply_butterfly_algorithm(config, resolved.explicit_algorithm); + } else if (!cached_algorithm.empty()) { + apply_butterfly_algorithm(config, cached_algorithm); + created->selection_source = CUBUTTERFLY_SELECTION_USER_CACHE; + } else if (resolved.algorithm_policy == CUBUTTERFLY_ALGORITHM_MEASURE) { + const auto tuned = tune_butterfly(config, resolved.storage, handle->stream); + apply_butterfly_algorithm(config, tuned.first); + created->selection_source = CUBUTTERFLY_SELECTION_RUNTIME_MEASURED; + created->reason = "measured legal butterfly candidates on the current hardware; kernel_ms=" + + std::to_string(tuned.second); + store_cache(handle->cache_path, cache_key, tuned.first); + log_event(handle, CUBUTTERFLY_LOG_INFO, "measure-complete", + created->reason + "; algorithm=" + tuned.first); + } else { + choose_butterfly_mapping(config, measured); + } + created->kind = cubutterflyPlan::Kind::Butterfly; + if (created->reason.empty()) + created->reason = created->selection_source == CUBUTTERFLY_SELECTION_USER_CACHE + ? "matched an exact workload key in the user selection cache" + : (measured ? "matched the built-in V100 measured profile" + : "no measured profile entry; selected by the static resource model"); + created->butterfly = std::make_unique(config); + created->butterfly->set_stream(handle->stream); + created->algorithm = std::string(cuntt::butterfly_backend_name(created->butterfly->config().backend)) + + "-" + cuntt::compute_unit_name(created->butterfly->config().compute_unit); + created->workspace_bytes = created->butterfly->workspace_size(); + } + if (!measured && resolved.algorithm_policy == CUBUTTERFLY_ALGORITHM_DEFAULT && cached_algorithm.empty()) + log_profile_miss(handle, resolved, created->reason, created->algorithm); + *plan = created.release(); + }); +} + +cubutterflyStatus_t cubutterflyDestroyPlan(cubutterflyPlan_t plan) { delete plan; return CUBUTTERFLY_STATUS_SUCCESS; } + +cubutterflyStatus_t cubutterflyPlanGetWorkspaceSize(cubutterflyPlan_t plan, size_t* bytes) { + return protect([&] { require(plan != nullptr && bytes != nullptr, "plan and output must be non-null"); *bytes = plan->workspace_bytes; }); +} + +cubutterflyStatus_t cubutterflyPlanGetInputSize(cubutterflyPlan_t plan, size_t* bytes) { return protect([&] { require(plan != nullptr && bytes != nullptr, "plan and output must be non-null"); *bytes = plan->input_bytes; }); } +cubutterflyStatus_t cubutterflyPlanGetOutputSize(cubutterflyPlan_t plan, size_t* bytes) { return protect([&] { require(plan != nullptr && bytes != nullptr, "plan and output must be non-null"); *bytes = plan->output_bytes; }); } + +cubutterflyStatus_t cubutterflyPlanGetPhysicalExtents(cubutterflyPlan_t plan, size_t* extents, uint32_t capacity) { + return protect([&] { require(plan != nullptr && extents != nullptr, "plan and extents must be non-null"); require(capacity >= plan->physical_extents.size(), "extent capacity is too small"); std::copy(plan->physical_extents.begin(), plan->physical_extents.end(), extents); }); +} + +cubutterflyStatus_t cubutterflyPlanGetSelectionSource(cubutterflyPlan_t plan, cubutterflySelectionSource_t* source) { return protect([&] { require(plan != nullptr && source != nullptr, "plan and output must be non-null"); *source = plan->selection_source; }); } +cubutterflyStatus_t cubutterflyPlanGetAlgorithmName(cubutterflyPlan_t plan, char* name, size_t* bytes) { return protect([&] { require(plan != nullptr, "plan must be non-null"); copy_text(plan->algorithm, name, bytes); }); } +cubutterflyStatus_t cubutterflyPlanGetSelectionReason(cubutterflyPlan_t plan, char* reason, size_t* bytes) { return protect([&] { require(plan != nullptr, "plan must be non-null"); copy_text(plan->reason, reason, bytes); }); } + +cubutterflyStatus_t cubutterflyPlanSetWorkspace(cubutterflyPlan_t plan, void* workspace, size_t bytes) { + return protect([&] { + require(plan != nullptr, "plan must be non-null"); + if (plan->kind == cubutterflyPlan::Kind::DirectFft) { + require(bytes == 0, "direct FFT does not require caller workspace"); + plan->workspace = workspace; + return; + } + if (plan->kind == cubutterflyPlan::Kind::BluesteinFft || + plan->kind == cubutterflyPlan::Kind::BluesteinNtt) { + require(workspace != nullptr, "Bluestein plans require a non-null workspace"); + require(bytes >= plan->workspace_bytes, "workspace is smaller than the Bluestein plan requirement"); + require(reinterpret_cast(workspace) % 16 == 0, + "workspace must be at least 16-byte aligned"); + plan->workspace = workspace; + return; + } + if (!plan->composite) { + if (plan->kind == cubutterflyPlan::Kind::Ntt) + plan->ntt->set_workspace(workspace, bytes); + else + plan->butterfly->set_workspace(workspace, bytes); + plan->workspace = workspace; + return; + } + if (plan->workspace_bytes == 0) { + require(bytes == 0, "zero-workspace composite plan received nonzero workspace bytes"); + plan->workspace = workspace; + return; + } + require(workspace != nullptr, "composite plans require a non-null workspace"); + require(bytes >= plan->workspace_bytes, "workspace is smaller than the composite plan requirement"); + require(reinterpret_cast(workspace) % 16 == 0, "workspace must be at least 16-byte aligned"); + plan->workspace = workspace; + std::size_t offset = align_16(plan->packed_bytes) * plan->composite_buffers; + auto* base = static_cast(workspace); + if (plan->kind == cubutterflyPlan::Kind::Ntt) { + for (auto& axis_plan : plan->ntt_axes) { + const std::size_t required = axis_plan->workspace_size(); + if (required != 0) + axis_plan->set_workspace(base + offset, required); + offset += align_16(required); + } + } else { + for (auto& axis_plan : plan->butterfly_axes) { + const std::size_t required = axis_plan->workspace_size(); + if (required != 0) + axis_plan->set_workspace(base + offset, required); + offset += align_16(required); + } + } + }); +} + +cubutterflyStatus_t cubutterflyExecute(cubutterflyHandle_t handle, cubutterflyPlan_t plan, const void* input, void* output) { + return protect([&] { + require(handle != nullptr && plan != nullptr && input != nullptr && output != nullptr, "execution arguments must be non-null"); + const bool in_place = plan->descriptor.placement == CUBUTTERFLY_PLACEMENT_IN_PLACE; + require(in_place ? input == output : input != output, + in_place ? "in-place execution requires identical pointers" : "out-of-place execution requires distinct pointers"); + if (plan->kind == cubutterflyPlan::Kind::DirectFft) { + plan->direct_fft->execute(input, output, handle->stream); + return; + } + if (plan->kind == cubutterflyPlan::Kind::BluesteinFft || + plan->kind == cubutterflyPlan::Kind::BluesteinNtt) { + require(plan->workspace != nullptr, "Bluestein plan requires a bound workspace"); + if (plan->descriptor.rank == 1) { + if (plan->kind == cubutterflyPlan::Kind::BluesteinNtt) + plan->bluestein_ntt_axes[0]->execute(static_cast(input), + static_cast(output), plan->workspace, + plan->descriptor.input_strides[0], + plan->descriptor.input_batch_stride, + plan->descriptor.output_strides[0], + plan->descriptor.output_batch_stride, + handle->stream); + else + plan->bluestein_axes[0]->execute(input, output, plan->workspace, + plan->descriptor.input_strides[0], + plan->descriptor.input_batch_stride, + plan->descriptor.output_strides[0], + plan->descriptor.output_batch_stride, + handle->stream); + } else { + auto* first = static_cast(plan->workspace); + auto* second = first + align_16(plan->packed_bytes); + auto* work = second + align_16(plan->packed_bytes); + const std::size_t rows = plan->descriptor.extents[0]; + const std::size_t cols = plan->descriptor.extents[1]; + cuntt::detail::launch_pack_2d(input, first, rows, cols, rows, cols, + plan->descriptor.batch, + plan->descriptor.input_strides[0], + plan->descriptor.input_strides[1], + plan->descriptor.input_batch_stride, + plan->element_bytes, handle->stream); + if (plan->kind == cubutterflyPlan::Kind::BluesteinNtt) + plan->bluestein_ntt_axes[0]->execute(reinterpret_cast(first), + reinterpret_cast(second), work, + 1, cols, 1, cols, handle->stream); + else + plan->bluestein_axes[0]->execute(first, second, work, 1, cols, 1, cols, handle->stream); + cuntt::detail::launch_transpose_2d(second, first, rows, cols, + plan->descriptor.batch, plan->element_bytes, handle->stream); + if (plan->kind == cubutterflyPlan::Kind::BluesteinNtt) + plan->bluestein_ntt_axes[1]->execute(reinterpret_cast(first), + reinterpret_cast(second), work, + 1, rows, 1, rows, handle->stream); + else + plan->bluestein_axes[1]->execute(first, second, work, 1, rows, 1, rows, handle->stream); + cuntt::detail::launch_transpose_2d(second, first, cols, rows, + plan->descriptor.batch, plan->element_bytes, handle->stream); + cuntt::detail::launch_scatter_2d(first, output, rows, cols, rows, cols, + plan->descriptor.batch, + plan->descriptor.output_strides[0], + plan->descriptor.output_strides[1], + plan->descriptor.output_batch_stride, + plan->element_bytes, handle->stream); + } + return; + } + if (plan->composite) { + require(plan->workspace_bytes == 0 || plan->workspace != nullptr, + "composite plan requires a bound workspace"); + auto* first = static_cast(plan->workspace); + auto* second = plan->composite_buffers >= 2 + ? first + align_16(plan->packed_bytes) + : first; + const bool inverse_embedding = plan->descriptor.length_mode == CUBUTTERFLY_LENGTH_ZERO_EXTENDED_EMBEDDING && + plan->descriptor.direction == CUBUTTERFLY_DIRECTION_INVERSE; + const auto& logical = plan->descriptor.extents; + const auto& physical = plan->physical_extents; + const bool direct_butterfly_1d = plan->direct_output_boundary && + plan->descriptor.rank == 1 && + plan->kind == cubutterflyPlan::Kind::Butterfly; + if (plan->descriptor.rank == 1) { + const std::size_t input_n = inverse_embedding ? physical[0] : logical[0]; + if (!plan->fused_input_boundary) { + void* pack_destination = direct_butterfly_1d ? output : first; + cuntt::detail::launch_pack_1d(input, pack_destination, input_n, physical[0], plan->descriptor.batch, + plan->descriptor.input_strides[0], plan->descriptor.input_batch_stride, + plan->element_bytes, handle->stream); + } + } else { + const std::size_t input_rows = inverse_embedding ? physical[0] : logical[0]; + const std::size_t input_cols = inverse_embedding ? physical[1] : logical[1]; + cuntt::detail::launch_pack_2d(input, first, input_rows, input_cols, physical[0], physical[1], + plan->descriptor.batch, plan->descriptor.input_strides[0], + plan->descriptor.input_strides[1], plan->descriptor.input_batch_stride, + plan->element_bytes, handle->stream); + } + void* packed_result = first; + if (plan->kind == cubutterflyPlan::Kind::Ntt) { + for (auto& axis_plan : plan->ntt_axes) + axis_plan->set_stream(handle->stream); + void* ntt_destination = plan->direct_output_boundary && plan->descriptor.rank == 1 + ? output + : static_cast(second); + if (plan->descriptor.storage == CUBUTTERFLY_DATA_UINT32) + plan->ntt_axes[0]->execute_async(reinterpret_cast(first), reinterpret_cast(ntt_destination)); + else + plan->ntt_axes[0]->execute_async(reinterpret_cast(first), reinterpret_cast(ntt_destination)); + packed_result = ntt_destination; + if (plan->descriptor.rank == 2) { + cuntt::detail::launch_transpose_2d(second, first, physical[0], physical[1], + plan->descriptor.batch, plan->element_bytes, handle->stream); + if (plan->descriptor.storage == CUBUTTERFLY_DATA_UINT32) + plan->ntt_axes[1]->execute_async(reinterpret_cast(first), reinterpret_cast(second)); + else + plan->ntt_axes[1]->execute_async(reinterpret_cast(first), reinterpret_cast(second)); + void* final_destination = plan->direct_output_boundary + ? output + : static_cast(first); + cuntt::detail::launch_transpose_2d(second, final_destination, physical[1], physical[0], + plan->descriptor.batch, plan->element_bytes, handle->stream); + packed_result = final_destination; + } + } else { + for (auto& axis_plan : plan->butterfly_axes) + axis_plan->set_stream(handle->stream); + void* butterfly_first = direct_butterfly_1d ? output : static_cast(first); + if (plan->fused_input_boundary) { + plan->butterfly_axes[0]->execute_zero_extended_async( + static_cast(input), static_cast(output), + logical[0], plan->descriptor.input_batch_stride, + plan->descriptor.input_strides[0]); + } else { + execute_butterfly(*plan->butterfly_axes[0], butterfly_first, butterfly_first, + plan->descriptor.storage); + } + packed_result = butterfly_first; + if (plan->descriptor.rank == 2) { + cuntt::detail::launch_transpose_2d(first, second, physical[0], physical[1], + plan->descriptor.batch, plan->element_bytes, handle->stream); + execute_butterfly(*plan->butterfly_axes[1], second, second, plan->descriptor.storage); + void* final_destination = plan->direct_output_boundary + ? output + : static_cast(first); + cuntt::detail::launch_transpose_2d(second, final_destination, physical[1], physical[0], + plan->descriptor.batch, plan->element_bytes, handle->stream); + packed_result = final_destination; + } + } + const bool forward_embedding = plan->descriptor.length_mode == CUBUTTERFLY_LENGTH_ZERO_EXTENDED_EMBEDDING && + plan->descriptor.direction == CUBUTTERFLY_DIRECTION_FORWARD; + if (plan->direct_output_boundary) { + return; + } else if (plan->descriptor.rank == 1) { + const std::size_t output_n = forward_embedding ? physical[0] : logical[0]; + cuntt::detail::launch_scatter_1d(packed_result, output, output_n, physical[0], plan->descriptor.batch, + plan->descriptor.output_strides[0], plan->descriptor.output_batch_stride, + plan->element_bytes, handle->stream); + } else { + const std::size_t output_rows = forward_embedding ? physical[0] : logical[0]; + const std::size_t output_cols = forward_embedding ? physical[1] : logical[1]; + cuntt::detail::launch_scatter_2d(packed_result, output, output_rows, output_cols, physical[0], physical[1], + plan->descriptor.batch, plan->descriptor.output_strides[0], + plan->descriptor.output_strides[1], plan->descriptor.output_batch_stride, + plan->element_bytes, handle->stream); + } + return; + } + if (plan->kind == cubutterflyPlan::Kind::Ntt) { + plan->ntt->set_stream(handle->stream); + if (plan->ntt->config().word_bits == 32) + plan->ntt->execute_async(static_cast(input), static_cast(output)); + else + plan->ntt->execute_async(static_cast(input), static_cast(output)); + } else { + plan->butterfly->set_stream(handle->stream); + const auto precision = plan->butterfly->config().precision; + cubutterflyDataType_t storage = CUBUTTERFLY_DATA_FP32; + if (plan->butterfly->config().op == cuntt::ButterflyOperator::Fft) { + if (precision == cuntt::ButterflyPrecision::Fp16) storage = CUBUTTERFLY_DATA_COMPLEX_FP16; + else if (precision == cuntt::ButterflyPrecision::Bf16) storage = CUBUTTERFLY_DATA_COMPLEX_BF16; + else if (precision == cuntt::ButterflyPrecision::Fp64) storage = CUBUTTERFLY_DATA_COMPLEX_FP64; + else storage = CUBUTTERFLY_DATA_COMPLEX_FP32; + } else if (precision == cuntt::ButterflyPrecision::Fp16) storage = CUBUTTERFLY_DATA_FP16; + else if (precision == cuntt::ButterflyPrecision::Bf16) storage = CUBUTTERFLY_DATA_BF16; + else if (precision == cuntt::ButterflyPrecision::Fp64) storage = CUBUTTERFLY_DATA_FP64; + else if (precision == cuntt::ButterflyPrecision::Uint32) storage = CUBUTTERFLY_DATA_UINT32; + execute_butterfly(*plan->butterfly, input, output, storage); + } + }); +} diff --git a/src/direct_fft.cu b/src/direct_fft.cu new file mode 100644 index 0000000..a427301 --- /dev/null +++ b/src/direct_fft.cu @@ -0,0 +1,114 @@ +#include "direct_fft.hpp" + +#include "cuntt/butterfly.hpp" + +#include + +#include +#include +#include +#include +#include + +namespace cuntt::detail { +namespace { + +void check_cufft(cufftResult status, const char* operation) { + if (status != CUFFT_SUCCESS) + throw std::runtime_error(std::string(operation) + ": cuFFT status " + std::to_string(status)); +} + +template +__global__ void scale_kernel(Complex* values, std::size_t count, double scale) { + const std::size_t index = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + if (index >= count) + return; + values[index].real = static_cast(values[index].real * scale); + values[index].imag = static_cast(values[index].imag * scale); +} + +unsigned int blocks_for(std::size_t count) { + const auto blocks = (count + 255) / 256; + if (blocks > UINT_MAX) + throw std::overflow_error("direct FFT normalization grid exceeds grid.x"); + return static_cast(blocks); +} + +} // namespace + +class DirectFftPlan::Impl { + public: + Impl(const std::vector& extents, std::size_t batch, + bool inverse, bool normalize_inverse, bool fp64) + : inverse_(inverse), normalize_inverse_(normalize_inverse), fp64_(fp64), batch_(batch) { + if (extents.empty() || extents.size() > 2 || batch == 0 || batch > static_cast(INT_MAX)) + throw std::invalid_argument("invalid direct FFT shape"); + std::vector dimensions(extents.size()); + points_ = 1; + for (std::size_t axis = 0; axis < extents.size(); ++axis) { + if (extents[axis] == 0 || extents[axis] > static_cast(INT_MAX) || + points_ > static_cast(INT_MAX) / extents[axis]) + throw std::invalid_argument("direct FFT extent exceeds cuFFT limits"); + dimensions[axis] = static_cast(extents[axis]); + points_ *= extents[axis]; + } + check_cufft(cufftPlanMany(&plan_, static_cast(dimensions.size()), dimensions.data(), + nullptr, 1, static_cast(points_), nullptr, 1, + static_cast(points_), fp64_ ? CUFFT_Z2Z : CUFFT_C2C, + static_cast(batch_)), "create direct FFT plan"); + } + + ~Impl() { + if (plan_ != 0) + cufftDestroy(plan_); + } + + void execute(const void* input, void* output, cudaStream_t stream) { + if (input == nullptr || output == nullptr) + throw std::invalid_argument("direct FFT input and output must be non-null"); + check_cufft(cufftSetStream(plan_, stream), "set direct FFT stream"); + const int direction = inverse_ ? CUFFT_INVERSE : CUFFT_FORWARD; + if (fp64_) + check_cufft(cufftExecZ2Z(plan_, reinterpret_cast(const_cast(input)), + static_cast(output), direction), + "execute direct FP64 FFT"); + else + check_cufft(cufftExecC2C(plan_, reinterpret_cast(const_cast(input)), + static_cast(output), direction), + "execute direct FP32 FFT"); + if (inverse_ && normalize_inverse_) { + const std::size_t count = points_ * batch_; + const double scale = 1.0 / static_cast(points_); + if (fp64_) + scale_kernel<<>>( + static_cast(output), count, scale); + else + scale_kernel<<>>( + static_cast(output), count, scale); + const auto status = cudaGetLastError(); + if (status != cudaSuccess) + throw std::runtime_error(std::string("launch direct FFT normalization: ") + + cudaGetErrorString(status)); + } + } + + private: + cufftHandle plan_ = 0; + bool inverse_ = false; + bool normalize_inverse_ = false; + bool fp64_ = false; + std::size_t batch_ = 0; + std::size_t points_ = 0; +}; + +DirectFftPlan::DirectFftPlan(const std::vector& extents, std::size_t batch, + bool inverse, bool normalize_inverse, bool fp64) + : impl_(std::make_unique(extents, batch, inverse, normalize_inverse, fp64)) {} + +DirectFftPlan::~DirectFftPlan() = default; + +void DirectFftPlan::execute(const void* input, void* output, cudaStream_t stream) { + impl_->execute(input, output, stream); +} + +} // namespace cuntt::detail diff --git a/src/direct_fft.hpp b/src/direct_fft.hpp new file mode 100644 index 0000000..95963f5 --- /dev/null +++ b/src/direct_fft.hpp @@ -0,0 +1,27 @@ +#pragma once + +#include + +#include +#include +#include + +namespace cuntt::detail { + +class DirectFftPlan { + public: + DirectFftPlan(const std::vector& extents, std::size_t batch, + bool inverse, bool normalize_inverse, bool fp64); + ~DirectFftPlan(); + + DirectFftPlan(const DirectFftPlan&) = delete; + DirectFftPlan& operator=(const DirectFftPlan&) = delete; + + void execute(const void* input, void* output, cudaStream_t stream); + + private: + class Impl; + std::unique_ptr impl_; +}; + +} // namespace cuntt::detail diff --git a/src/general_kernels.cu b/src/general_kernels.cu new file mode 100644 index 0000000..58a4957 --- /dev/null +++ b/src/general_kernels.cu @@ -0,0 +1,183 @@ +#include "general_kernels.hpp" + +#include + +#include +#include +#include +#include + +namespace cuntt::detail { +namespace { + +void check_cuda(cudaError_t status, const char* operation) { + if (status != cudaSuccess) + throw std::runtime_error(std::string(operation) + ": " + cudaGetErrorString(status)); +} + +struct alignas(16) Bytes16 { std::uint64_t lo; std::uint64_t hi; }; + +template +__global__ void pack_1d_kernel(const T* input, T* packed, std::size_t logical_n, + std::size_t physical_n, std::size_t batch, + std::size_t input_stride, std::size_t input_batch_stride) { + const std::size_t linear = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + const std::size_t count = logical_n * batch; + if (linear >= count) + return; + const std::size_t transform = linear / logical_n; + const std::size_t index = linear - transform * logical_n; + packed[transform * physical_n + index] = input[transform * input_batch_stride + index * input_stride]; +} + +template +__global__ void scatter_1d_kernel(const T* packed, T* output, std::size_t output_n, + std::size_t physical_n, std::size_t batch, + std::size_t output_stride, std::size_t output_batch_stride) { + const std::size_t linear = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + const std::size_t count = output_n * batch; + if (linear >= count) + return; + const std::size_t transform = linear / output_n; + const std::size_t index = linear - transform * output_n; + output[transform * output_batch_stride + index * output_stride] = packed[transform * physical_n + index]; +} + +template +__global__ void pack_2d_kernel(const T* input, T* packed, + std::size_t logical_rows, std::size_t logical_cols, + std::size_t physical_rows, std::size_t physical_cols, + std::size_t batch, std::size_t input_row_stride, + std::size_t input_col_stride, std::size_t input_batch_stride) { + const std::size_t linear = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + const std::size_t matrix_points = logical_rows * logical_cols; + const std::size_t count = matrix_points * batch; + if (linear >= count) + return; + const std::size_t transform = linear / matrix_points; + const std::size_t within = linear - transform * matrix_points; + const std::size_t row = within / logical_cols; + const std::size_t col = within - row * logical_cols; + packed[(transform * physical_rows + row) * physical_cols + col] = + input[transform * input_batch_stride + row * input_row_stride + col * input_col_stride]; +} + +template +__global__ void scatter_2d_kernel(const T* packed, T* output, + std::size_t output_rows, std::size_t output_cols, + std::size_t physical_rows, std::size_t physical_cols, + std::size_t batch, std::size_t output_row_stride, + std::size_t output_col_stride, std::size_t output_batch_stride) { + const std::size_t linear = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + const std::size_t matrix_points = output_rows * output_cols; + const std::size_t count = matrix_points * batch; + if (linear >= count) + return; + const std::size_t transform = linear / matrix_points; + const std::size_t within = linear - transform * matrix_points; + const std::size_t row = within / output_cols; + const std::size_t col = within - row * output_cols; + output[transform * output_batch_stride + row * output_row_stride + col * output_col_stride] = + packed[(transform * physical_rows + row) * physical_cols + col]; +} + +template +__global__ void transpose_2d_kernel(const T* input, T* output, + std::size_t rows, std::size_t cols, std::size_t batch) { + const std::size_t linear = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + const std::size_t matrix_points = rows * cols; + const std::size_t count = matrix_points * batch; + if (linear >= count) + return; + const std::size_t transform = linear / matrix_points; + const std::size_t within = linear - transform * matrix_points; + const std::size_t row = within / cols; + const std::size_t col = within - row * cols; + output[(transform * cols + col) * rows + row] = input[linear]; +} + +template +void dispatch_size(std::size_t bytes, Function&& function) { + switch (bytes) { + case 2: function(static_cast(nullptr)); return; + case 4: function(static_cast(nullptr)); return; + case 8: function(static_cast(nullptr)); return; + case 16: function(static_cast(nullptr)); return; + default: throw std::invalid_argument("general transform element size must be 2, 4, 8, or 16 bytes"); + } +} + +unsigned int blocks_for(std::size_t count) { + constexpr std::size_t threads = 256; + const std::size_t blocks = (count + threads - 1) / threads; + if (blocks > static_cast(UINT32_MAX)) + throw std::overflow_error("general transform copy grid exceeds CUDA grid.x"); + return static_cast(blocks); +} + +} // namespace + +void launch_pack_1d(const void* input, void* packed, std::size_t logical_n, std::size_t physical_n, + std::size_t batch, std::size_t input_stride, std::size_t input_batch_stride, + std::size_t bytes, cudaStream_t stream) { + check_cuda(cudaMemsetAsync(packed, 0, physical_n * batch * bytes, stream), "clear packed rank-1 buffer"); + dispatch_size(bytes, [&](auto type) { + using T = std::remove_pointer_t; + pack_1d_kernel<<>>( + static_cast(input), static_cast(packed), logical_n, physical_n, + batch, input_stride, input_batch_stride); + }); + check_cuda(cudaGetLastError(), "launch rank-1 pack"); +} + +void launch_scatter_1d(const void* packed, void* output, std::size_t output_n, std::size_t physical_n, + std::size_t batch, std::size_t output_stride, std::size_t output_batch_stride, + std::size_t bytes, cudaStream_t stream) { + dispatch_size(bytes, [&](auto type) { + using T = std::remove_pointer_t; + scatter_1d_kernel<<>>( + static_cast(packed), static_cast(output), output_n, physical_n, + batch, output_stride, output_batch_stride); + }); + check_cuda(cudaGetLastError(), "launch rank-1 scatter"); +} + +void launch_pack_2d(const void* input, void* packed, std::size_t logical_rows, std::size_t logical_cols, + std::size_t physical_rows, std::size_t physical_cols, std::size_t batch, + std::size_t input_row_stride, std::size_t input_col_stride, + std::size_t input_batch_stride, std::size_t bytes, cudaStream_t stream) { + check_cuda(cudaMemsetAsync(packed, 0, physical_rows * physical_cols * batch * bytes, stream), + "clear packed rank-2 buffer"); + dispatch_size(bytes, [&](auto type) { + using T = std::remove_pointer_t; + pack_2d_kernel<<>>( + static_cast(input), static_cast(packed), logical_rows, logical_cols, + physical_rows, physical_cols, batch, input_row_stride, input_col_stride, input_batch_stride); + }); + check_cuda(cudaGetLastError(), "launch rank-2 pack"); +} + +void launch_scatter_2d(const void* packed, void* output, std::size_t output_rows, std::size_t output_cols, + std::size_t physical_rows, std::size_t physical_cols, std::size_t batch, + std::size_t output_row_stride, std::size_t output_col_stride, + std::size_t output_batch_stride, std::size_t bytes, cudaStream_t stream) { + dispatch_size(bytes, [&](auto type) { + using T = std::remove_pointer_t; + scatter_2d_kernel<<>>( + static_cast(packed), static_cast(output), output_rows, output_cols, + physical_rows, physical_cols, batch, output_row_stride, output_col_stride, output_batch_stride); + }); + check_cuda(cudaGetLastError(), "launch rank-2 scatter"); +} + +void launch_transpose_2d(const void* input, void* output, std::size_t rows, std::size_t cols, + std::size_t batch, std::size_t bytes, cudaStream_t stream) { + dispatch_size(bytes, [&](auto type) { + using T = std::remove_pointer_t; + transpose_2d_kernel<<>>( + static_cast(input), static_cast(output), rows, cols, batch); + }); + check_cuda(cudaGetLastError(), "launch rank-2 transpose"); +} + +} // namespace cuntt::detail diff --git a/src/general_kernels.hpp b/src/general_kernels.hpp new file mode 100644 index 0000000..34bab70 --- /dev/null +++ b/src/general_kernels.hpp @@ -0,0 +1,39 @@ +#pragma once + +#include + +#include + +namespace cuntt::detail { + +void launch_pack_1d(const void* input, void* packed, + std::size_t logical_n, std::size_t physical_n, + std::size_t batch, std::size_t input_stride, + std::size_t input_batch_stride, std::size_t element_bytes, + cudaStream_t stream); + +void launch_scatter_1d(const void* packed, void* output, + std::size_t output_n, std::size_t physical_n, + std::size_t batch, std::size_t output_stride, + std::size_t output_batch_stride, std::size_t element_bytes, + cudaStream_t stream); + +void launch_pack_2d(const void* input, void* packed, + std::size_t logical_rows, std::size_t logical_cols, + std::size_t physical_rows, std::size_t physical_cols, + std::size_t batch, std::size_t input_row_stride, + std::size_t input_col_stride, std::size_t input_batch_stride, + std::size_t element_bytes, cudaStream_t stream); + +void launch_scatter_2d(const void* packed, void* output, + std::size_t output_rows, std::size_t output_cols, + std::size_t physical_rows, std::size_t physical_cols, + std::size_t batch, std::size_t output_row_stride, + std::size_t output_col_stride, std::size_t output_batch_stride, + std::size_t element_bytes, cudaStream_t stream); + +void launch_transpose_2d(const void* input, void* output, + std::size_t rows, std::size_t cols, std::size_t batch, + std::size_t element_bytes, cudaStream_t stream); + +} // namespace cuntt::detail diff --git a/src/generated_unit_api.cuh b/src/generated_unit_api.cuh index 335fd71..455fe79 100644 --- a/src/generated_unit_api.cuh +++ b/src/generated_unit_api.cuh @@ -12,6 +12,12 @@ bool generated_register_fwht_available(std::uint32_t log_n) noexcept; std::uint32_t generated_register_fwht_threads(std::uint32_t log_n) noexcept; void launch_generated_register_fwht(std::uint32_t log_n, const float* input, float* output, std::uint64_t transforms, std::uint64_t batch_distance, std::uint64_t element_stride, bool normalize, cudaStream_t stream); +void launch_generated_register_fwht_zero_extended( + std::uint32_t log_n, const float* input, float* output, + std::uint64_t transforms, std::uint64_t input_batch_distance, + std::uint64_t output_batch_distance, std::uint64_t input_element_stride, + std::uint64_t output_element_stride, std::uint32_t logical_points, + cudaStream_t stream); bool generated_register_structured_available(std::uint32_t log_n) noexcept; std::uint32_t generated_register_structured_threads(std::uint32_t log_n) noexcept; @@ -19,6 +25,12 @@ void launch_generated_register_structured(std::uint32_t log_n, const float* inpu const float* matrices, std::uint64_t transforms, std::uint64_t batch_distance, std::uint64_t element_stride, bool broadcast, cudaStream_t stream); +void launch_generated_register_structured_zero_extended( + std::uint32_t log_n, const float* input, float* output, + const float* matrices, std::uint64_t transforms, + std::uint64_t input_batch_distance, std::uint64_t output_batch_distance, + std::uint64_t input_element_stride, std::uint64_t output_element_stride, + std::uint32_t logical_points, bool broadcast, cudaStream_t stream); bool generated_thread_dft8_available(std::uint32_t log_n) noexcept; void launch_generated_thread_dft8(std::uint32_t log_n, const Complex32* input, Complex32* output, const Complex32* twiddles, std::uint64_t transforms, diff --git a/src/register_fwht.cuh b/src/register_fwht.cuh index a301dd7..bbe9fbd 100644 --- a/src/register_fwht.cuh +++ b/src/register_fwht.cuh @@ -161,8 +161,12 @@ __device__ __forceinline__ void exchange_warps(float (&values)[Traits::kChunks][ template __global__ __launch_bounds__(RegisterFwhtTraits::kThreads) void register_fwht_kernel(const float* input, float* output, - std::uint64_t transforms, std::uint64_t batch_distance, - std::uint64_t element_stride, float scale) { + std::uint64_t transforms, + std::uint64_t input_batch_distance, + std::uint64_t output_batch_distance, + std::uint64_t input_element_stride, + std::uint64_t output_element_stride, + std::uint32_t logical_points, float scale) { using Traits = RegisterFwhtTraits; extern __shared__ __align__(16) unsigned char storage[]; float4* shared = reinterpret_cast(storage); @@ -170,18 +174,23 @@ __global__ __launch_bounds__(RegisterFwhtTraits::kThreads) void register_f if (transform >= transforms) return; - const std::uint64_t base = transform * batch_distance; + const std::uint64_t input_base = transform * input_batch_distance; + const std::uint64_t output_base = transform * output_batch_distance; float values[Traits::kChunks][4]; #pragma unroll for (std::uint32_t chunk = 0; chunk < Traits::kChunks; ++chunk) { const std::uint32_t vector_index = chunk * Traits::kThreads + threadIdx.x; - if (element_stride == 1 && (base & 3U) == 0) { - unpack_float4(reinterpret_cast(input + base)[vector_index], values[chunk]); + const std::uint32_t vector_base = vector_index * 4; + if (input_element_stride == 1 && (input_base & 3U) == 0 && + vector_base + 3 < logical_points) { + unpack_float4(reinterpret_cast(input + input_base)[vector_index], values[chunk]); } else { #pragma unroll for (std::uint32_t item = 0; item < 4; ++item) { - const std::uint32_t index = vector_index * 4 + item; - values[chunk][item] = input[base + static_cast(index) * element_stride]; + const std::uint32_t index = vector_base + item; + values[chunk][item] = index < logical_points + ? input[input_base + static_cast(index) * input_element_stride] + : 0.0F; } } register_wht(values[chunk]); @@ -216,18 +225,18 @@ __global__ __launch_bounds__(RegisterFwhtTraits::kThreads) void register_f #pragma unroll for (std::uint32_t chunk = 0; chunk < Traits::kChunks; ++chunk) { const std::uint32_t vector_index = chunk * Traits::kThreads + threadIdx.x; - if (element_stride == 1 && (base & 3U) == 0) { + if (output_element_stride == 1 && (output_base & 3U) == 0) { float4 result = pack_float4(values[chunk]); result.x *= scale; result.y *= scale; result.z *= scale; result.w *= scale; - reinterpret_cast(output + base)[vector_index] = result; + reinterpret_cast(output + output_base)[vector_index] = result; } else { #pragma unroll for (std::uint32_t item = 0; item < 4; ++item) { const std::uint32_t index = vector_index * 4 + item; - output[base + static_cast(index) * element_stride] = values[chunk][item] * scale; + output[output_base + static_cast(index) * output_element_stride] = values[chunk][item] * scale; } } } @@ -239,13 +248,32 @@ void launch_register_fwht(const float* input, float* output, std::uint64_t trans using Traits = RegisterFwhtTraits; const float scale = normalize ? 1.0F / static_cast(Traits::kN) : 1.0F; register_fwht_kernel<<(transforms), Traits::kThreads, Traits::kSharedBytes, stream>>>( - input, output, transforms, batch_distance, element_stride, scale); + input, output, transforms, batch_distance, batch_distance, element_stride, + element_stride, Traits::kN, scale); +} + +template +void launch_register_fwht_zero_extended(const float* input, float* output, + std::uint64_t transforms, + std::uint64_t input_batch_distance, + std::uint64_t output_batch_distance, + std::uint64_t input_element_stride, + std::uint64_t output_element_stride, + std::uint32_t logical_points, + cudaStream_t stream) { + using Traits = RegisterFwhtTraits; + register_fwht_kernel<<(transforms), Traits::kThreads, + Traits::kSharedBytes, stream>>>( + input, output, transforms, input_batch_distance, output_batch_distance, + input_element_stride, output_element_stride, logical_points, 1.0F); } template __global__ __launch_bounds__(RegisterFwhtTraits::kThreads) void register_structured_kernel( const float* input, float* output, const float* matrix_entries, std::uint64_t transforms, - std::uint64_t batch_distance, std::uint64_t element_stride) { + std::uint64_t input_batch_distance, std::uint64_t output_batch_distance, + std::uint64_t input_element_stride, std::uint64_t output_element_stride, + std::uint32_t logical_points) { using Traits = RegisterFwhtTraits; extern __shared__ __align__(16) unsigned char storage[]; float4* shared = reinterpret_cast(storage); @@ -255,18 +283,23 @@ __global__ __launch_bounds__(RegisterFwhtTraits::kThreads) void register_s if (transform >= transforms) return; - const std::uint64_t base = transform * batch_distance; + const std::uint64_t input_base = transform * input_batch_distance; + const std::uint64_t output_base = transform * output_batch_distance; float values[Traits::kChunks][4]; #pragma unroll for (std::uint32_t chunk = 0; chunk < Traits::kChunks; ++chunk) { const std::uint32_t vector_index = chunk * Traits::kThreads + threadIdx.x; - if (element_stride == 1 && (base & 3U) == 0) { - unpack_float4(reinterpret_cast(input + base)[vector_index], values[chunk]); + const std::uint32_t vector_base = vector_index * 4; + if (input_element_stride == 1 && (input_base & 3U) == 0 && + vector_base + 3 < logical_points) { + unpack_float4(reinterpret_cast(input + input_base)[vector_index], values[chunk]); } else { #pragma unroll for (std::uint32_t item = 0; item < 4; ++item) { - const std::uint32_t index = vector_index * 4 + item; - values[chunk][item] = input[base + static_cast(index) * element_stride]; + const std::uint32_t index = vector_base + item; + values[chunk][item] = index < logical_points + ? input[input_base + static_cast(index) * input_element_stride] + : 0.0F; } } register_structured(values[chunk], matrices, broadcast_matrix); @@ -307,13 +340,13 @@ __global__ __launch_bounds__(RegisterFwhtTraits::kThreads) void register_s #pragma unroll for (std::uint32_t chunk = 0; chunk < Traits::kChunks; ++chunk) { const std::uint32_t vector_index = chunk * Traits::kThreads + threadIdx.x; - if (element_stride == 1 && (base & 3U) == 0) { - reinterpret_cast(output + base)[vector_index] = pack_float4(values[chunk]); + if (output_element_stride == 1 && (output_base & 3U) == 0) { + reinterpret_cast(output + output_base)[vector_index] = pack_float4(values[chunk]); } else { #pragma unroll for (std::uint32_t item = 0; item < 4; ++item) { const std::uint32_t index = vector_index * 4 + item; - output[base + static_cast(index) * element_stride] = values[chunk][item]; + output[output_base + static_cast(index) * output_element_stride] = values[chunk][item]; } } } @@ -327,11 +360,38 @@ void launch_register_structured(const float* input, float* output, const float* if (broadcast) { register_structured_kernel<<(transforms), Traits::kThreads, Traits::kSharedBytes, stream>>>( - input, output, matrices, transforms, batch_distance, element_stride); + input, output, matrices, transforms, batch_distance, batch_distance, + element_stride, element_stride, Traits::kN); } else { register_structured_kernel<<(transforms), Traits::kThreads, Traits::kSharedBytes, stream>>>( - input, output, matrices, transforms, batch_distance, element_stride); + input, output, matrices, transforms, batch_distance, batch_distance, + element_stride, element_stride, Traits::kN); + } +} + +template +void launch_register_structured_zero_extended( + const float* input, float* output, const float* matrices, + std::uint64_t transforms, std::uint64_t input_batch_distance, + std::uint64_t output_batch_distance, std::uint64_t input_element_stride, + std::uint64_t output_element_stride, std::uint32_t logical_points, + bool broadcast, cudaStream_t stream) { + using Traits = RegisterFwhtTraits; + if (broadcast) { + register_structured_kernel + <<(transforms), Traits::kThreads, + Traits::kSharedBytes, stream>>>( + input, output, matrices, transforms, input_batch_distance, + output_batch_distance, input_element_stride, + output_element_stride, logical_points); + } else { + register_structured_kernel + <<(transforms), Traits::kThreads, + Traits::kSharedBytes, stream>>>( + input, output, matrices, transforms, input_batch_distance, + output_batch_distance, input_element_stride, + output_element_stride, logical_points); } } diff --git a/tests/test_application_profiles.py b/tests/test_application_profiles.py new file mode 100644 index 0000000..8c64afd --- /dev/null +++ b/tests/test_application_profiles.py @@ -0,0 +1,37 @@ +#!/usr/bin/env python3 +import importlib.util +import json +import pathlib +import tempfile +import unittest + + +ROOT = pathlib.Path(__file__).resolve().parents[1] +SPEC = importlib.util.spec_from_file_location( + "generate_application_profiles", ROOT / "scripts/generate_application_profiles.py") +MODULE = importlib.util.module_from_spec(SPEC) +SPEC.loader.exec_module(MODULE) + + +class ApplicationProfileGeneratorTest(unittest.TestCase): + def test_v100_profiles_generate_exact_shape_batch_and_modulus_metadata(self): + source = ROOT / "config/v100_application_profiles.json" + document = json.loads(source.read_text()) + self.assertEqual(document["schema_version"], 1) + self.assertEqual(document["target"]["multiprocessors"], 80) + self.assertTrue(any(profile["operator"] == "fft" and profile["logN"] == 20 + for profile in document["profiles"])) + ntt32 = next(profile for profile in document["profiles"] + if profile["operator"] == "ntt" and profile["storage"] == "uint32") + self.assertEqual(ntt32["modulus"], 1073692673) + with tempfile.TemporaryDirectory() as directory: + output = pathlib.Path(directory) / "profiles.hpp" + MODULE.generate(source, output) + text = output.read_text() + self.assertIn("kApplicationProfileMultiprocessors = 80", text) + self.assertIn("1073692673ULL", text) + self.assertIn("cuButterfly-cuFFTDx-online", text) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/test_butterfly.cpp b/tests/test_butterfly.cpp index 8689613..5f93599 100644 --- a/tests/test_butterfly.cpp +++ b/tests/test_butterfly.cpp @@ -444,6 +444,60 @@ void test_structured_operators() { } } +void test_zero_extended_register_units() { + constexpr std::size_t logical_n = 6; + constexpr std::size_t physical_n = 8; + constexpr std::size_t batch = 3; + const std::vector input{ + 1, 2, 3, 4, 5, 6, + -1, 0.5F, 2, -3, 4, -0.25F, + 0.125F, -0.5F, 1.5F, 2.5F, -4, 8}; + const std::vector matrices{{1.0, 0.25, -0.5, 1.0}}; + TestStream stream; + TestDeviceBuffer device_input(input.size() * sizeof(float)); + TestDeviceBuffer device_output(batch * physical_n * sizeof(float)); + check_cuda(cudaMemcpyAsync(device_input.get(), input.data(), input.size() * sizeof(float), + cudaMemcpyHostToDevice, stream.get()), + "cudaMemcpyAsync zero-extended input"); + + for (const auto op : {cuntt::ButterflyOperator::Fwht, + cuntt::ButterflyOperator::Structured2x2}) { + cuntt::ButterflyConfig config; + config.op = op; + config.backend = cuntt::ButterflyBackend::TemporalTile; + config.precision = cuntt::ButterflyPrecision::Fp32; + config.log_n = 3; + config.batch = batch; + config.local_exchange = cuntt::LocalExchange::WarpRegister; + if (op == cuntt::ButterflyOperator::Structured2x2) + config.stage_matrices = matrices; + cuntt::ButterflyPlan plan(config); + plan.set_stream(stream.get()); + plan.execute_zero_extended_async( + static_cast(device_input.get()), + static_cast(device_output.get()), logical_n, logical_n); + + std::vector output(batch * physical_n); + check_cuda(cudaMemcpyAsync(output.data(), device_output.get(), output.size() * sizeof(float), + cudaMemcpyDeviceToHost, stream.get()), + "cudaMemcpyAsync zero-extended output"); + check_cuda(cudaStreamSynchronize(stream.get()), "cudaStreamSynchronize zero-extended"); + for (std::size_t transform = 0; transform < batch; ++transform) { + std::vector expected(physical_n, 0.0F); + std::copy_n(input.begin() + transform * logical_n, logical_n, expected.begin()); + if (op == cuntt::ButterflyOperator::Fwht) + cuntt::reference_fwht(expected); + else + cuntt::reference_structured_2x2(expected, matrices); + for (std::size_t index = 0; index < physical_n; ++index) { + if (std::abs(output[transform * physical_n + index] - expected[index]) > 1.0e-4F) + throw std::runtime_error("zero-extended warp-register processing-unit mismatch"); + } + } + } + std::cout << "PASS zero-extended warp-register processing units\n"; +} + void test_fp64(std::uint32_t log_n, bool inverse) { const std::size_t n = std::size_t{1} << log_n; std::mt19937 random(0x6400U + log_n + static_cast(inverse)); @@ -728,6 +782,7 @@ int main() { } test_extended_zeta_operators(); test_structured_operators(); + test_zero_extended_register_units(); for (const std::uint32_t tile_threads : {32U, 64U, 128U, 256U}) { test_fwht(cuntt::ButterflyBackend::TemporalTile, 1, 5, 8, tile_threads); test_fft(cuntt::ButterflyBackend::TemporalTile, 1, 5, 8, tile_threads); diff --git a/tests/test_c_api.c b/tests/test_c_api.c new file mode 100644 index 0000000..35cb40d --- /dev/null +++ b/tests/test_c_api.c @@ -0,0 +1,689 @@ +#include + +#include + +#include +#include +#include +#include +#include + +#define CHECK_CUB(call) do { \ + cubutterflyStatus_t status_ = (call); \ + if (status_ != CUBUTTERFLY_STATUS_SUCCESS) { \ + fprintf(stderr, "%s failed: %s\n", #call, cubutterflyGetStatusString(status_)); \ + return 1; \ + } \ +} while (0) + +#define CHECK_CUDA(call) do { \ + cudaError_t status_ = (call); \ + if (status_ != cudaSuccess) { \ + fprintf(stderr, "%s failed: %s\n", #call, cudaGetErrorString(status_)); \ + return 1; \ + } \ +} while (0) + +static void log_message(cubutterflyLogLevel_t level, const char* event, + const char* message, void* user_data) { + int* warnings = (int*)user_data; + if (level == CUBUTTERFLY_LOG_WARNING) + ++*warnings; + (void)event; + (void)message; +} + +typedef struct test_complex32 { + float real; + float imag; +} test_complex32_t; + +static int run_arbitrary_fft(cubutterflyHandle_t handle, cudaStream_t stream) { + const size_t extent = 5; + test_complex32_t original[5] = {{0, 0}, {1, 0}, {0, 0}, {0, 0}, {0, 0}}; + test_complex32_t spectrum[5]; + test_complex32_t recovered[5]; + test_complex32_t* input = NULL; + test_complex32_t* output = NULL; + test_complex32_t* result = NULL; + void* workspace = NULL; + cubutterflyDescriptor_t descriptor = NULL; + cubutterflyPlan_t forward = NULL; + cubutterflyPlan_t inverse = NULL; + size_t physical_extent = 0; + size_t workspace_bytes = 0; + CHECK_CUB(cubutterflyCreateDescriptor(&descriptor)); + CHECK_CUB(cubutterflySetOperator(descriptor, CUBUTTERFLY_OPERATOR_FFT)); + CHECK_CUB(cubutterflySetDataType(descriptor, CUBUTTERFLY_DATA_COMPLEX_FP32, CUBUTTERFLY_COMPUTE_FP32)); + CHECK_CUB(cubutterflySetShape(descriptor, 1, &extent, 1)); + CHECK_CUB(cubutterflySetAlgorithmPolicy(descriptor, CUBUTTERFLY_ALGORITHM_EXPLICIT, + "bluestein-cufft-power2-core")); + CHECK_CUB(cubutterflyCreatePlan(handle, descriptor, &forward)); + CHECK_CUB(cubutterflyPlanGetPhysicalExtents(forward, &physical_extent, 1)); + CHECK_CUB(cubutterflyPlanGetWorkspaceSize(forward, &workspace_bytes)); + if (physical_extent != 16 || workspace_bytes < physical_extent * sizeof(test_complex32_t)) + return 1; + CHECK_CUDA(cudaMalloc((void**)&input, sizeof(original))); + CHECK_CUDA(cudaMalloc((void**)&output, sizeof(spectrum))); + CHECK_CUDA(cudaMalloc((void**)&result, sizeof(recovered))); + CHECK_CUDA(cudaMalloc(&workspace, workspace_bytes)); + CHECK_CUB(cubutterflyPlanSetWorkspace(forward, workspace, workspace_bytes)); + CHECK_CUDA(cudaMemcpyAsync(input, original, sizeof(original), cudaMemcpyHostToDevice, stream)); + CHECK_CUB(cubutterflyExecute(handle, forward, input, output)); + CHECK_CUDA(cudaMemcpyAsync(spectrum, output, sizeof(spectrum), cudaMemcpyDeviceToHost, stream)); + CHECK_CUDA(cudaStreamSynchronize(stream)); + for (size_t k = 0; k < extent; ++k) { + const double angle = -2.0 * 3.14159265358979323846 * (double)k / (double)extent; + if (fabsf(spectrum[k].real - (float)cos(angle)) > 2.0e-5f || + fabsf(spectrum[k].imag - (float)sin(angle)) > 2.0e-5f) + return 1; + } + CHECK_CUB(cubutterflySetDirection(descriptor, CUBUTTERFLY_DIRECTION_INVERSE, 1)); + CHECK_CUB(cubutterflyCreatePlan(handle, descriptor, &inverse)); + CHECK_CUB(cubutterflyPlanSetWorkspace(inverse, workspace, workspace_bytes)); + CHECK_CUB(cubutterflyExecute(handle, inverse, output, result)); + CHECK_CUDA(cudaMemcpyAsync(recovered, result, sizeof(recovered), cudaMemcpyDeviceToHost, stream)); + CHECK_CUDA(cudaStreamSynchronize(stream)); + for (size_t i = 0; i < extent; ++i) { + if (fabsf(recovered[i].real - original[i].real) > 2.0e-5f || + fabsf(recovered[i].imag - original[i].imag) > 2.0e-5f) + return 1; + } + cudaFree(workspace); + cudaFree(result); + cudaFree(output); + cudaFree(input); + cubutterflyDestroyPlan(inverse); + cubutterflyDestroyPlan(forward); + cubutterflyDestroyDescriptor(descriptor); + return 0; +} + +static int run_2d_arbitrary_fft(cubutterflyHandle_t handle, cudaStream_t stream) { + const size_t extents[2] = {3, 5}; + test_complex32_t original[15] = {{0, 0}}; + test_complex32_t spectrum[15]; + test_complex32_t recovered[15]; + test_complex32_t* input = NULL; + test_complex32_t* output = NULL; + test_complex32_t* result = NULL; + void* workspace = NULL; + cubutterflyDescriptor_t descriptor = NULL; + cubutterflyPlan_t forward = NULL; + cubutterflyPlan_t inverse = NULL; + size_t physical[2] = {0, 0}; + size_t workspace_bytes = 0; + original[1 * extents[1] + 2].real = 1.0f; + CHECK_CUB(cubutterflyCreateDescriptor(&descriptor)); + CHECK_CUB(cubutterflySetOperator(descriptor, CUBUTTERFLY_OPERATOR_FFT)); + CHECK_CUB(cubutterflySetDataType(descriptor, CUBUTTERFLY_DATA_COMPLEX_FP32, CUBUTTERFLY_COMPUTE_FP32)); + CHECK_CUB(cubutterflySetShape(descriptor, 2, extents, 1)); + CHECK_CUB(cubutterflyCreatePlan(handle, descriptor, &forward)); + CHECK_CUB(cubutterflyPlanGetPhysicalExtents(forward, physical, 2)); + CHECK_CUB(cubutterflyPlanGetWorkspaceSize(forward, &workspace_bytes)); + if (physical[0] != extents[0] || physical[1] != extents[1]) + return 1; + CHECK_CUDA(cudaMalloc((void**)&input, sizeof(original))); + CHECK_CUDA(cudaMalloc((void**)&output, sizeof(spectrum))); + CHECK_CUDA(cudaMalloc((void**)&result, sizeof(recovered))); + CHECK_CUDA(cudaMalloc(&workspace, workspace_bytes)); + CHECK_CUB(cubutterflyPlanSetWorkspace(forward, workspace, workspace_bytes)); + CHECK_CUDA(cudaMemcpyAsync(input, original, sizeof(original), cudaMemcpyHostToDevice, stream)); + CHECK_CUB(cubutterflyExecute(handle, forward, input, output)); + CHECK_CUDA(cudaMemcpyAsync(spectrum, output, sizeof(spectrum), cudaMemcpyDeviceToHost, stream)); + CHECK_CUDA(cudaStreamSynchronize(stream)); + for (size_t row = 0; row < extents[0]; ++row) { + for (size_t col = 0; col < extents[1]; ++col) { + const size_t index = row * extents[1] + col; + const double angle = -2.0 * 3.14159265358979323846 * + ((double)row / (double)extents[0] + 2.0 * (double)col / (double)extents[1]); + if (fabsf(spectrum[index].real - (float)cos(angle)) > 4.0e-5f || + fabsf(spectrum[index].imag - (float)sin(angle)) > 4.0e-5f) + return 1; + } + } + CHECK_CUB(cubutterflySetDirection(descriptor, CUBUTTERFLY_DIRECTION_INVERSE, 1)); + CHECK_CUB(cubutterflyCreatePlan(handle, descriptor, &inverse)); + CHECK_CUB(cubutterflyPlanSetWorkspace(inverse, workspace, workspace_bytes)); + CHECK_CUB(cubutterflyExecute(handle, inverse, output, result)); + CHECK_CUDA(cudaMemcpyAsync(recovered, result, sizeof(recovered), cudaMemcpyDeviceToHost, stream)); + CHECK_CUDA(cudaStreamSynchronize(stream)); + for (size_t i = 0; i < 15; ++i) { + if (fabsf(recovered[i].real - original[i].real) > 4.0e-5f || + fabsf(recovered[i].imag - original[i].imag) > 4.0e-5f) + return 1; + } + cudaFree(workspace); + cudaFree(result); + cudaFree(output); + cudaFree(input); + cubutterflyDestroyPlan(inverse); + cubutterflyDestroyPlan(forward); + cubutterflyDestroyDescriptor(descriptor); + return 0; +} + +static int run_2d_power_fft(cubutterflyHandle_t handle, cudaStream_t stream) { + const size_t extents[2] = {4, 8}; + test_complex32_t original[32]; + test_complex32_t recovered[32]; + test_complex32_t* input = NULL; + test_complex32_t* spectrum = NULL; + test_complex32_t* output = NULL; + void* workspace = NULL; + cubutterflyDescriptor_t descriptor = NULL; + cubutterflyPlan_t forward = NULL; + cubutterflyPlan_t inverse = NULL; + size_t workspace_bytes = 0; + size_t inverse_workspace_bytes = 0; + for (size_t i = 0; i < 32; ++i) { + original[i].real = (float)((i * 7) % 13) / 13.0f; + original[i].imag = (float)((i * 5) % 11) / 17.0f; + } + CHECK_CUB(cubutterflyCreateDescriptor(&descriptor)); + CHECK_CUB(cubutterflySetOperator(descriptor, CUBUTTERFLY_OPERATOR_FFT)); + CHECK_CUB(cubutterflySetDataType(descriptor, CUBUTTERFLY_DATA_COMPLEX_FP32, CUBUTTERFLY_COMPUTE_FP32)); + CHECK_CUB(cubutterflySetShape(descriptor, 2, extents, 1)); + CHECK_CUB(cubutterflyCreatePlan(handle, descriptor, &forward)); + CHECK_CUB(cubutterflyPlanGetWorkspaceSize(forward, &workspace_bytes)); + CHECK_CUDA(cudaMalloc((void**)&input, sizeof(original))); + CHECK_CUDA(cudaMalloc((void**)&spectrum, sizeof(original))); + CHECK_CUDA(cudaMalloc((void**)&output, sizeof(original))); + CHECK_CUDA(cudaMalloc(&workspace, workspace_bytes)); + CHECK_CUB(cubutterflyPlanSetWorkspace(forward, workspace, workspace_bytes)); + CHECK_CUDA(cudaMemcpyAsync(input, original, sizeof(original), cudaMemcpyHostToDevice, stream)); + CHECK_CUB(cubutterflyExecute(handle, forward, input, spectrum)); + CHECK_CUB(cubutterflySetDirection(descriptor, CUBUTTERFLY_DIRECTION_INVERSE, 1)); + CHECK_CUB(cubutterflyCreatePlan(handle, descriptor, &inverse)); + CHECK_CUB(cubutterflyPlanGetWorkspaceSize(inverse, &inverse_workspace_bytes)); + if (inverse_workspace_bytes > workspace_bytes) + return 1; + CHECK_CUB(cubutterflyPlanSetWorkspace(inverse, workspace, workspace_bytes)); + CHECK_CUB(cubutterflyExecute(handle, inverse, spectrum, output)); + CHECK_CUDA(cudaMemcpyAsync(recovered, output, sizeof(recovered), cudaMemcpyDeviceToHost, stream)); + CHECK_CUDA(cudaStreamSynchronize(stream)); + for (size_t i = 0; i < 32; ++i) { + if (fabsf(recovered[i].real - original[i].real) > 2.0e-4f || + fabsf(recovered[i].imag - original[i].imag) > 2.0e-4f) + return 1; + } + cudaFree(workspace); + cudaFree(output); + cudaFree(spectrum); + cudaFree(input); + cubutterflyDestroyPlan(inverse); + cubutterflyDestroyPlan(forward); + cubutterflyDestroyDescriptor(descriptor); + return 0; +} + +static int run_2d_ntt(cubutterflyHandle_t handle, cudaStream_t stream) { + const size_t extents[2] = {4, 8}; + uint64_t original[32]; + uint64_t recovered[32]; + uint64_t* input = NULL; + uint64_t* spectrum = NULL; + uint64_t* output = NULL; + void* workspace = NULL; + cubutterflyDescriptor_t descriptor = NULL; + cubutterflyPlan_t forward = NULL; + cubutterflyPlan_t inverse = NULL; + size_t workspace_bytes = 0; + size_t inverse_workspace_bytes = 0; + for (size_t i = 0; i < 32; ++i) + original[i] = (uint64_t)(i * i + 3 * i + 1); + CHECK_CUB(cubutterflyCreateDescriptor(&descriptor)); + CHECK_CUB(cubutterflySetOperator(descriptor, CUBUTTERFLY_OPERATOR_NTT)); + CHECK_CUB(cubutterflySetDataType(descriptor, CUBUTTERFLY_DATA_UINT64, CUBUTTERFLY_COMPUTE_UINT64)); + CHECK_CUB(cubutterflySetShape(descriptor, 2, extents, 1)); + CHECK_CUB(cubutterflyCreatePlan(handle, descriptor, &forward)); + CHECK_CUB(cubutterflyPlanGetWorkspaceSize(forward, &workspace_bytes)); + CHECK_CUDA(cudaMalloc((void**)&input, sizeof(original))); + CHECK_CUDA(cudaMalloc((void**)&spectrum, sizeof(original))); + CHECK_CUDA(cudaMalloc((void**)&output, sizeof(original))); + CHECK_CUDA(cudaMalloc(&workspace, workspace_bytes)); + CHECK_CUB(cubutterflyPlanSetWorkspace(forward, workspace, workspace_bytes)); + CHECK_CUDA(cudaMemcpyAsync(input, original, sizeof(original), cudaMemcpyHostToDevice, stream)); + CHECK_CUB(cubutterflyExecute(handle, forward, input, spectrum)); + CHECK_CUB(cubutterflySetDirection(descriptor, CUBUTTERFLY_DIRECTION_INVERSE, 1)); + CHECK_CUB(cubutterflyCreatePlan(handle, descriptor, &inverse)); + CHECK_CUB(cubutterflyPlanGetWorkspaceSize(inverse, &inverse_workspace_bytes)); + if (inverse_workspace_bytes > workspace_bytes) + return 1; + CHECK_CUB(cubutterflyPlanSetWorkspace(inverse, workspace, workspace_bytes)); + CHECK_CUB(cubutterflyExecute(handle, inverse, spectrum, output)); + CHECK_CUDA(cudaMemcpyAsync(recovered, output, sizeof(recovered), cudaMemcpyDeviceToHost, stream)); + CHECK_CUDA(cudaStreamSynchronize(stream)); + for (size_t i = 0; i < 32; ++i) + if (recovered[i] != original[i]) + return 1; + cudaFree(workspace); + cudaFree(output); + cudaFree(spectrum); + cudaFree(input); + cubutterflyDestroyPlan(inverse); + cubutterflyDestroyPlan(forward); + cubutterflyDestroyDescriptor(descriptor); + return 0; +} + +static int run_arbitrary_ntt(cubutterflyHandle_t handle, cudaStream_t stream) { + const size_t extent = 3; + uint64_t original[3] = {7, 11, 19}; + uint64_t recovered[3]; + uint64_t* input = NULL; + uint64_t* spectrum = NULL; + uint64_t* output = NULL; + void* workspace = NULL; + cubutterflyDescriptor_t descriptor = NULL; + cubutterflyPlan_t forward = NULL; + cubutterflyPlan_t inverse = NULL; + size_t workspace_bytes = 0; + size_t inverse_workspace_bytes = 0; + CHECK_CUB(cubutterflyCreateDescriptor(&descriptor)); + CHECK_CUB(cubutterflySetOperator(descriptor, CUBUTTERFLY_OPERATOR_NTT)); + CHECK_CUB(cubutterflySetDataType(descriptor, CUBUTTERFLY_DATA_UINT64, CUBUTTERFLY_COMPUTE_UINT64)); + CHECK_CUB(cubutterflySetShape(descriptor, 1, &extent, 1)); + CHECK_CUB(cubutterflyCreatePlan(handle, descriptor, &forward)); + CHECK_CUB(cubutterflyPlanGetWorkspaceSize(forward, &workspace_bytes)); + CHECK_CUDA(cudaMalloc((void**)&input, sizeof(original))); + CHECK_CUDA(cudaMalloc((void**)&spectrum, sizeof(original))); + CHECK_CUDA(cudaMalloc((void**)&output, sizeof(original))); + CHECK_CUDA(cudaMalloc(&workspace, workspace_bytes)); + CHECK_CUB(cubutterflyPlanSetWorkspace(forward, workspace, workspace_bytes)); + CHECK_CUDA(cudaMemcpyAsync(input, original, sizeof(original), cudaMemcpyHostToDevice, stream)); + CHECK_CUB(cubutterflyExecute(handle, forward, input, spectrum)); + CHECK_CUB(cubutterflySetDirection(descriptor, CUBUTTERFLY_DIRECTION_INVERSE, 1)); + CHECK_CUB(cubutterflyCreatePlan(handle, descriptor, &inverse)); + CHECK_CUB(cubutterflyPlanGetWorkspaceSize(inverse, &inverse_workspace_bytes)); + if (inverse_workspace_bytes > workspace_bytes) + return 1; + CHECK_CUB(cubutterflyPlanSetWorkspace(inverse, workspace, workspace_bytes)); + CHECK_CUB(cubutterflyExecute(handle, inverse, spectrum, output)); + CHECK_CUDA(cudaMemcpyAsync(recovered, output, sizeof(recovered), cudaMemcpyDeviceToHost, stream)); + CHECK_CUDA(cudaStreamSynchronize(stream)); + for (size_t i = 0; i < extent; ++i) + if (recovered[i] != original[i]) + return 1; + cudaFree(workspace); + cudaFree(output); + cudaFree(spectrum); + cudaFree(input); + cubutterflyDestroyPlan(inverse); + cubutterflyDestroyPlan(forward); + cubutterflyDestroyDescriptor(descriptor); + return 0; +} + +static int run_measure_cache(cubutterflyHandle_t handle) { + const char* cache_path = "/tmp/cubutterfly_c_api_test.cache"; + size_t extent = 16; + cubutterflyDescriptor_t descriptor = NULL; + cubutterflyPlan_t measured = NULL; + cubutterflyPlan_t cached = NULL; + cubutterflySelectionSource_t source = CUBUTTERFLY_SELECTION_STATIC_MODEL; + remove(cache_path); + CHECK_CUB(cubutterflySetCachePath(handle, cache_path)); + CHECK_CUB(cubutterflyCreateDescriptor(&descriptor)); + CHECK_CUB(cubutterflySetOperator(descriptor, CUBUTTERFLY_OPERATOR_FWHT)); + CHECK_CUB(cubutterflySetDataType(descriptor, CUBUTTERFLY_DATA_FP32, CUBUTTERFLY_COMPUTE_FP32)); + CHECK_CUB(cubutterflySetShape(descriptor, 1, &extent, 8)); + CHECK_CUB(cubutterflySetAlgorithmPolicy(descriptor, CUBUTTERFLY_ALGORITHM_MEASURE, NULL)); + CHECK_CUB(cubutterflyCreatePlan(handle, descriptor, &measured)); + CHECK_CUB(cubutterflyPlanGetSelectionSource(measured, &source)); + if (source != CUBUTTERFLY_SELECTION_RUNTIME_MEASURED) + return 1; + cubutterflyDestroyPlan(measured); + CHECK_CUB(cubutterflySetAlgorithmPolicy(descriptor, CUBUTTERFLY_ALGORITHM_DEFAULT, NULL)); + CHECK_CUB(cubutterflyCreatePlan(handle, descriptor, &cached)); + CHECK_CUB(cubutterflyPlanGetSelectionSource(cached, &source)); + if (source != CUBUTTERFLY_SELECTION_USER_CACHE) + return 1; + cubutterflyDestroyPlan(cached); + + extent = 15; + measured = NULL; + cached = NULL; + CHECK_CUB(cubutterflySetShape(descriptor, 1, &extent, 8)); + CHECK_CUB(cubutterflySetLengthMode(descriptor, CUBUTTERFLY_LENGTH_ZERO_EXTENDED_EMBEDDING)); + CHECK_CUB(cubutterflySetAlgorithmPolicy(descriptor, CUBUTTERFLY_ALGORITHM_MEASURE, NULL)); + CHECK_CUB(cubutterflyCreatePlan(handle, descriptor, &measured)); + CHECK_CUB(cubutterflyPlanGetSelectionSource(measured, &source)); + if (source != CUBUTTERFLY_SELECTION_RUNTIME_MEASURED) + return 1; + cubutterflyDestroyPlan(measured); + CHECK_CUB(cubutterflySetAlgorithmPolicy(descriptor, CUBUTTERFLY_ALGORITHM_DEFAULT, NULL)); + CHECK_CUB(cubutterflyCreatePlan(handle, descriptor, &cached)); + CHECK_CUB(cubutterflyPlanGetSelectionSource(cached, &source)); + if (source != CUBUTTERFLY_SELECTION_USER_CACHE) + return 1; + cubutterflyDestroyPlan(cached); + + cubutterflyDestroyDescriptor(descriptor); + CHECK_CUB(cubutterflySetCachePath(handle, NULL)); + remove(cache_path); + return 0; +} + +static int run_2d_zeta(cubutterflyHandle_t handle, cudaStream_t stream) { + const size_t extents[2] = {4, 8}; + uint32_t original[32]; + uint32_t recovered[32]; + uint32_t* input = NULL; + uint32_t* transformed = NULL; + uint32_t* output = NULL; + void* workspace = NULL; + cubutterflyDescriptor_t descriptor = NULL; + cubutterflyPlan_t forward = NULL; + cubutterflyPlan_t inverse = NULL; + size_t workspace_bytes = 0; + size_t inverse_workspace_bytes = 0; + for (size_t i = 0; i < 32; ++i) + original[i] = (uint32_t)(i * 3 + 1); + CHECK_CUB(cubutterflyCreateDescriptor(&descriptor)); + CHECK_CUB(cubutterflySetOperator(descriptor, CUBUTTERFLY_OPERATOR_SUBSET_ZETA)); + CHECK_CUB(cubutterflySetDataType(descriptor, CUBUTTERFLY_DATA_UINT32, CUBUTTERFLY_COMPUTE_UINT32)); + CHECK_CUB(cubutterflySetShape(descriptor, 2, extents, 1)); + CHECK_CUB(cubutterflyCreatePlan(handle, descriptor, &forward)); + CHECK_CUB(cubutterflyPlanGetWorkspaceSize(forward, &workspace_bytes)); + CHECK_CUDA(cudaMalloc((void**)&input, sizeof(original))); + CHECK_CUDA(cudaMalloc((void**)&transformed, sizeof(original))); + CHECK_CUDA(cudaMalloc((void**)&output, sizeof(original))); + CHECK_CUDA(cudaMalloc(&workspace, workspace_bytes)); + CHECK_CUB(cubutterflyPlanSetWorkspace(forward, workspace, workspace_bytes)); + CHECK_CUDA(cudaMemcpyAsync(input, original, sizeof(original), cudaMemcpyHostToDevice, stream)); + CHECK_CUB(cubutterflyExecute(handle, forward, input, transformed)); + CHECK_CUB(cubutterflySetDirection(descriptor, CUBUTTERFLY_DIRECTION_INVERSE, 0)); + CHECK_CUB(cubutterflyCreatePlan(handle, descriptor, &inverse)); + CHECK_CUB(cubutterflyPlanGetWorkspaceSize(inverse, &inverse_workspace_bytes)); + if (inverse_workspace_bytes > workspace_bytes) + return 1; + CHECK_CUB(cubutterflyPlanSetWorkspace(inverse, workspace, workspace_bytes)); + CHECK_CUB(cubutterflyExecute(handle, inverse, transformed, output)); + CHECK_CUDA(cudaMemcpyAsync(recovered, output, sizeof(recovered), cudaMemcpyDeviceToHost, stream)); + CHECK_CUDA(cudaStreamSynchronize(stream)); + for (size_t i = 0; i < 32; ++i) + if (recovered[i] != original[i]) + return 1; + cudaFree(workspace); + cudaFree(output); + cudaFree(transformed); + cudaFree(input); + cubutterflyDestroyPlan(inverse); + cubutterflyDestroyPlan(forward); + cubutterflyDestroyDescriptor(descriptor); + return 0; +} + +static int run_2d_structured(cubutterflyHandle_t handle, cudaStream_t stream) { + const size_t extents[2] = {4, 8}; + const cubutterflyMatrix2x2_t matrix = {1.0, 0.25, -0.5, 1.0}; + float original[32]; + float recovered[32]; + float* input = NULL; + float* transformed = NULL; + float* output = NULL; + void* workspace = NULL; + cubutterflyDescriptor_t descriptor = NULL; + cubutterflyPlan_t forward = NULL; + cubutterflyPlan_t inverse = NULL; + size_t workspace_bytes = 0; + size_t inverse_workspace_bytes = 0; + for (size_t i = 0; i < 32; ++i) + original[i] = (float)((i * 5) % 17) / 17.0f; + CHECK_CUB(cubutterflyCreateDescriptor(&descriptor)); + CHECK_CUB(cubutterflySetOperator(descriptor, CUBUTTERFLY_OPERATOR_STRUCTURED_2X2)); + CHECK_CUB(cubutterflySetDataType(descriptor, CUBUTTERFLY_DATA_FP32, CUBUTTERFLY_COMPUTE_FP32)); + CHECK_CUB(cubutterflySetShape(descriptor, 2, extents, 1)); + CHECK_CUB(cubutterflySetStageMatrices(descriptor, 0, &matrix, 1)); + CHECK_CUB(cubutterflySetStageMatrices(descriptor, 1, &matrix, 1)); + CHECK_CUB(cubutterflyCreatePlan(handle, descriptor, &forward)); + CHECK_CUB(cubutterflyPlanGetWorkspaceSize(forward, &workspace_bytes)); + CHECK_CUDA(cudaMalloc((void**)&input, sizeof(original))); + CHECK_CUDA(cudaMalloc((void**)&transformed, sizeof(original))); + CHECK_CUDA(cudaMalloc((void**)&output, sizeof(original))); + CHECK_CUDA(cudaMalloc(&workspace, workspace_bytes)); + CHECK_CUB(cubutterflyPlanSetWorkspace(forward, workspace, workspace_bytes)); + CHECK_CUDA(cudaMemcpyAsync(input, original, sizeof(original), cudaMemcpyHostToDevice, stream)); + CHECK_CUB(cubutterflyExecute(handle, forward, input, transformed)); + CHECK_CUB(cubutterflySetDirection(descriptor, CUBUTTERFLY_DIRECTION_INVERSE, 0)); + CHECK_CUB(cubutterflyCreatePlan(handle, descriptor, &inverse)); + CHECK_CUB(cubutterflyPlanGetWorkspaceSize(inverse, &inverse_workspace_bytes)); + if (inverse_workspace_bytes > workspace_bytes) + return 1; + CHECK_CUB(cubutterflyPlanSetWorkspace(inverse, workspace, workspace_bytes)); + CHECK_CUB(cubutterflyExecute(handle, inverse, transformed, output)); + CHECK_CUDA(cudaMemcpyAsync(recovered, output, sizeof(recovered), cudaMemcpyDeviceToHost, stream)); + CHECK_CUDA(cudaStreamSynchronize(stream)); + for (size_t i = 0; i < 32; ++i) + if (fabsf(recovered[i] - original[i]) > 3.0e-4f) + return 1; + cudaFree(workspace); + cudaFree(output); + cudaFree(transformed); + cudaFree(input); + cubutterflyDestroyPlan(inverse); + cubutterflyDestroyPlan(forward); + cubutterflyDestroyDescriptor(descriptor); + return 0; +} + +static int run_2d_fwht(cubutterflyHandle_t handle, cudaStream_t stream) { + const size_t extents[2] = {4, 8}; + float host_input[32]; + float host_output[32]; + cubutterflyDescriptor_t descriptor = NULL; + cubutterflyPlan_t plan = NULL; + float* input = NULL; + float* output = NULL; + void* workspace = NULL; + size_t workspace_bytes = 0; + for (size_t i = 0; i < 32; ++i) + host_input[i] = 1.0f; + CHECK_CUB(cubutterflyCreateDescriptor(&descriptor)); + CHECK_CUB(cubutterflySetOperator(descriptor, CUBUTTERFLY_OPERATOR_FWHT)); + CHECK_CUB(cubutterflySetDataType(descriptor, CUBUTTERFLY_DATA_FP32, CUBUTTERFLY_COMPUTE_FP32)); + CHECK_CUB(cubutterflySetShape(descriptor, 2, extents, 1)); + CHECK_CUB(cubutterflyCreatePlan(handle, descriptor, &plan)); + CHECK_CUB(cubutterflyPlanGetWorkspaceSize(plan, &workspace_bytes)); + CHECK_CUDA(cudaMalloc((void**)&input, sizeof(host_input))); + CHECK_CUDA(cudaMalloc((void**)&output, sizeof(host_output))); + CHECK_CUDA(cudaMalloc(&workspace, workspace_bytes)); + CHECK_CUB(cubutterflyPlanSetWorkspace(plan, workspace, workspace_bytes)); + CHECK_CUDA(cudaMemcpyAsync(input, host_input, sizeof(host_input), cudaMemcpyHostToDevice, stream)); + CHECK_CUB(cubutterflyExecute(handle, plan, input, output)); + CHECK_CUDA(cudaMemcpyAsync(host_output, output, sizeof(host_output), cudaMemcpyDeviceToHost, stream)); + CHECK_CUDA(cudaStreamSynchronize(stream)); + if (fabsf(host_output[0] - 32.0f) > 1.0e-5f) + return 1; + for (size_t i = 1; i < 32; ++i) + if (fabsf(host_output[i]) > 1.0e-5f) + return 1; + cudaFree(workspace); + cudaFree(output); + cudaFree(input); + cubutterflyDestroyPlan(plan); + cubutterflyDestroyDescriptor(descriptor); + return 0; +} + +static int run_embedding_fwht(cubutterflyHandle_t handle, cudaStream_t stream) { + const size_t extent = 6; + float original[6] = {1, 2, 3, 4, 5, 6}; + float spectrum[8]; + float recovered[6]; + float* logical = NULL; + float* physical = NULL; + float* result = NULL; + void* workspace = NULL; + cubutterflyDescriptor_t descriptor = NULL; + cubutterflyPlan_t forward = NULL; + cubutterflyPlan_t inverse = NULL; + size_t workspace_bytes = 0; + size_t inverse_workspace_bytes = 0; + char algorithm[256]; + size_t algorithm_bytes = sizeof(algorithm); + CHECK_CUB(cubutterflyCreateDescriptor(&descriptor)); + CHECK_CUB(cubutterflySetOperator(descriptor, CUBUTTERFLY_OPERATOR_FWHT)); + CHECK_CUB(cubutterflySetDataType(descriptor, CUBUTTERFLY_DATA_FP32, CUBUTTERFLY_COMPUTE_FP32)); + CHECK_CUB(cubutterflySetShape(descriptor, 1, &extent, 1)); + CHECK_CUB(cubutterflySetLengthMode(descriptor, CUBUTTERFLY_LENGTH_ZERO_EXTENDED_EMBEDDING)); + CHECK_CUB(cubutterflySetAlgorithmPolicy(descriptor, CUBUTTERFLY_ALGORITHM_EXPLICIT, + "temporal-warp-register-radix2")); + CHECK_CUB(cubutterflyCreatePlan(handle, descriptor, &forward)); + CHECK_CUB(cubutterflyPlanGetAlgorithmName(forward, algorithm, &algorithm_bytes)); + if (strstr(algorithm, "+fused-input") == NULL) + return 1; + CHECK_CUB(cubutterflyPlanGetWorkspaceSize(forward, &workspace_bytes)); + CHECK_CUDA(cudaMalloc((void**)&logical, sizeof(original))); + CHECK_CUDA(cudaMalloc((void**)&physical, sizeof(spectrum))); + CHECK_CUDA(cudaMalloc((void**)&result, sizeof(recovered))); + CHECK_CUDA(cudaMalloc(&workspace, workspace_bytes)); + CHECK_CUB(cubutterflyPlanSetWorkspace(forward, workspace, workspace_bytes)); + CHECK_CUDA(cudaMemcpyAsync(logical, original, sizeof(original), cudaMemcpyHostToDevice, stream)); + CHECK_CUB(cubutterflyExecute(handle, forward, logical, physical)); + CHECK_CUDA(cudaMemcpyAsync(spectrum, physical, sizeof(spectrum), cudaMemcpyDeviceToHost, stream)); + CHECK_CUDA(cudaStreamSynchronize(stream)); + + CHECK_CUB(cubutterflySetDirection(descriptor, CUBUTTERFLY_DIRECTION_INVERSE, 1)); + CHECK_CUB(cubutterflyCreatePlan(handle, descriptor, &inverse)); + CHECK_CUB(cubutterflyPlanGetWorkspaceSize(inverse, &inverse_workspace_bytes)); + if (inverse_workspace_bytes > workspace_bytes) { + cudaFree(workspace); + CHECK_CUDA(cudaMalloc(&workspace, inverse_workspace_bytes)); + workspace_bytes = inverse_workspace_bytes; + } + CHECK_CUB(cubutterflyPlanSetWorkspace(inverse, workspace, workspace_bytes)); + CHECK_CUB(cubutterflyExecute(handle, inverse, physical, result)); + CHECK_CUDA(cudaMemcpyAsync(recovered, result, sizeof(recovered), cudaMemcpyDeviceToHost, stream)); + CHECK_CUDA(cudaStreamSynchronize(stream)); + for (size_t i = 0; i < 6; ++i) + if (fabsf(recovered[i] - original[i]) > 1.0e-5f) + return 1; + cudaFree(workspace); + cudaFree(result); + cudaFree(physical); + cudaFree(logical); + cubutterflyDestroyPlan(inverse); + cubutterflyDestroyPlan(forward); + cubutterflyDestroyDescriptor(descriptor); + return 0; +} + +static int run_strided_embedding_fwht(cubutterflyHandle_t handle, cudaStream_t stream) { + const size_t extent = 6; + const size_t input_stride[1] = {1}; + const size_t output_stride[1] = {2}; + float original[6] = {1, 2, 3, 4, 5, 6}; + float expected[8] = {1, 2, 3, 4, 5, 6, 0, 0}; + float strided[15]; + float* input = NULL; + float* output = NULL; + void* workspace = NULL; + cubutterflyDescriptor_t descriptor = NULL; + cubutterflyPlan_t plan = NULL; + size_t workspace_bytes = 0; + for (size_t width = 1; width < 8; width *= 2) { + for (size_t base = 0; base < 8; base += 2 * width) { + for (size_t offset = 0; offset < width; ++offset) { + const float left = expected[base + offset]; + const float right = expected[base + offset + width]; + expected[base + offset] = left + right; + expected[base + offset + width] = left - right; + } + } + } + CHECK_CUB(cubutterflyCreateDescriptor(&descriptor)); + CHECK_CUB(cubutterflySetOperator(descriptor, CUBUTTERFLY_OPERATOR_FWHT)); + CHECK_CUB(cubutterflySetDataType(descriptor, CUBUTTERFLY_DATA_FP32, CUBUTTERFLY_COMPUTE_FP32)); + CHECK_CUB(cubutterflySetShape(descriptor, 1, &extent, 1)); + CHECK_CUB(cubutterflySetLengthMode(descriptor, CUBUTTERFLY_LENGTH_ZERO_EXTENDED_EMBEDDING)); + CHECK_CUB(cubutterflySetStrides(descriptor, input_stride, extent, + output_stride, 15)); + CHECK_CUB(cubutterflyCreatePlan(handle, descriptor, &plan)); + CHECK_CUB(cubutterflyPlanGetWorkspaceSize(plan, &workspace_bytes)); + CHECK_CUDA(cudaMalloc((void**)&input, sizeof(original))); + CHECK_CUDA(cudaMalloc((void**)&output, sizeof(strided))); + CHECK_CUDA(cudaMalloc(&workspace, workspace_bytes)); + CHECK_CUB(cubutterflyPlanSetWorkspace(plan, workspace, workspace_bytes)); + CHECK_CUDA(cudaMemcpyAsync(input, original, sizeof(original), cudaMemcpyHostToDevice, stream)); + CHECK_CUDA(cudaMemsetAsync(output, 0, sizeof(strided), stream)); + CHECK_CUB(cubutterflyExecute(handle, plan, input, output)); + CHECK_CUDA(cudaMemcpyAsync(strided, output, sizeof(strided), cudaMemcpyDeviceToHost, stream)); + CHECK_CUDA(cudaStreamSynchronize(stream)); + for (size_t index = 0; index < 8; ++index) + if (fabsf(strided[2 * index] - expected[index]) > 1.0e-5f) + return 1; + cudaFree(workspace); + cudaFree(output); + cudaFree(input); + cubutterflyDestroyPlan(plan); + cubutterflyDestroyDescriptor(descriptor); + return 0; +} + +int main(void) { + cubutterflyHandle_t handle = NULL; + cubutterflyDescriptor_t descriptor = NULL; + cubutterflyPlan_t plan = NULL; + cudaStream_t stream = NULL; + float* device_input = NULL; + float* device_output = NULL; + float host_input[8]; + float host_output[8]; + size_t extent = 8; + size_t input_bytes = 0; + size_t output_bytes = 0; + int warnings = 0; + + for (size_t i = 0; i < extent; ++i) + host_input[i] = 1.0f; + + CHECK_CUDA(cudaStreamCreateWithFlags(&stream, cudaStreamNonBlocking)); + CHECK_CUB(cubutterflyCreate(&handle)); + CHECK_CUB(cubutterflySetStream(handle, stream)); + CHECK_CUB(cubutterflySetLogCallback(handle, log_message, &warnings)); + CHECK_CUB(cubutterflyCreateDescriptor(&descriptor)); + CHECK_CUB(cubutterflySetOperator(descriptor, CUBUTTERFLY_OPERATOR_FWHT)); + CHECK_CUB(cubutterflySetDataType(descriptor, CUBUTTERFLY_DATA_FP32, CUBUTTERFLY_COMPUTE_FP32)); + CHECK_CUB(cubutterflySetShape(descriptor, 1, &extent, 1)); + CHECK_CUB(cubutterflyCreatePlan(handle, descriptor, &plan)); + CHECK_CUB(cubutterflyPlanGetInputSize(plan, &input_bytes)); + CHECK_CUB(cubutterflyPlanGetOutputSize(plan, &output_bytes)); + if (input_bytes != sizeof(host_input) || output_bytes != sizeof(host_output)) + return 1; + + CHECK_CUDA(cudaMalloc((void**)&device_input, input_bytes)); + CHECK_CUDA(cudaMalloc((void**)&device_output, output_bytes)); + CHECK_CUDA(cudaMemcpyAsync(device_input, host_input, input_bytes, cudaMemcpyHostToDevice, stream)); + CHECK_CUB(cubutterflyExecute(handle, plan, device_input, device_output)); + CHECK_CUDA(cudaMemcpyAsync(host_output, device_output, output_bytes, cudaMemcpyDeviceToHost, stream)); + CHECK_CUDA(cudaStreamSynchronize(stream)); + if (fabsf(host_output[0] - 8.0f) > 1.0e-5f) + return 1; + + if (run_2d_fwht(handle, stream) != 0 || run_embedding_fwht(handle, stream) != 0 || + run_strided_embedding_fwht(handle, stream) != 0 || + run_arbitrary_fft(handle, stream) != 0) + return 1; + if (run_2d_arbitrary_fft(handle, stream) != 0) + return 1; + if (run_2d_power_fft(handle, stream) != 0 || run_2d_ntt(handle, stream) != 0) + return 1; + if (run_arbitrary_ntt(handle, stream) != 0) + return 1; + if (run_measure_cache(handle) != 0) + return 1; + if (run_2d_zeta(handle, stream) != 0 || run_2d_structured(handle, stream) != 0) + return 1; + if (warnings == 0) + return 1; + + cubutterflyDestroyPlan(plan); + cubutterflyDestroyDescriptor(descriptor); + cubutterflyDestroy(handle); + cudaFree(device_output); + cudaFree(device_input); + cudaStreamDestroy(stream); + return 0; +} From c0023ba5b18ae51dea1226d40b36d2b5a084592a Mon Sep 17 00:00:00 2001 From: TruNcat3 <31944766+TruNcat3@users.noreply.github.com> Date: Tue, 11 Aug 2026 14:08:50 +0800 Subject: [PATCH 2/2] Prepare v0.6.0 release --- CHANGELOG.md | 2 +- CITATION.cff | 4 ++-- README.md | 4 ++-- ROADMAP.md | 5 +++-- 4 files changed, 8 insertions(+), 7 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ccc32ff..45b00ca 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,7 +2,7 @@ All notable repository and research-artifact changes are recorded here. -## 0.6.0 - Unreleased General Shapes And Plan API +## 0.6.0 - 2026-08-11 General Shapes And Plan API ### Added diff --git a/CITATION.cff b/CITATION.cff index d5bd119..564b082 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -6,8 +6,8 @@ authors: - name: "TruNcat3" repository-code: "https://github.com/TruNcat3/cuButterfly" url: "https://github.com/TruNcat3/cuButterfly" -version: 0.5.0 -date-released: 2026-08-08 +version: 0.6.0 +date-released: 2026-08-11 keywords: - butterfly computation - GPU diff --git a/README.md b/README.md index 206b4f8..bea7c96 100644 --- a/README.md +++ b/README.md @@ -18,7 +18,7 @@ processing-unit boundary, same-machine library comparisons, Nsight profiling scripts, and the raw CSV data used in the reports. It is a research artifact, not a drop-in replacement for cuFFT or a production cryptography library. -The v0.6 development API adds a CUDA-library-style C plan surface, exact +The v0.6 API adds a CUDA-library-style C plan surface, exact non-power-of-two FFT/compatible NTT lowering, explicit power-of-two embedding, rank-two axis mapping, caller-owned workspace, static application profiles, one-time measured tuning, and exact-key cache reuse. See the @@ -397,7 +397,7 @@ software release as: author = {TruNcat3}, title = {cuButterfly: Hardware-Mapped Space-Time Parallelism for Butterfly Computations on GPUs}, year = {2026}, - version = {0.5.0}, + version = {0.6.0}, url = {https://github.com/TruNcat3/cuButterfly} } ``` diff --git a/ROADMAP.md b/ROADMAP.md index 199c74e..b544f6c 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -5,8 +5,9 @@ research claim. Ordering may change after profiling or access to new hardware. ## v0.6.0 Status: General Shapes And Professional Plan Surface -Status: implementation and single-V100 correctness/performance validation -complete; release packaging and review remain. +Status: complete on the single-V100 release target. General-shape semantics, +the professional plan surface, boundary selection, tests, raw evidence, and +release documentation are checked in. The new C plan API separates standard exact-length semantics from explicit zero-extended embedding, lowers rank-two workloads to per-axis mappings, and adds exact static-profile/cache selection with actionable miss logs. Exact FFT