diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 8d9cc76..c18a844 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -2,7 +2,7 @@ name: QuadBLAS Testing on: push: - branches: [ main ] + branches: [main, dev] pull_request: workflow_dispatch: @@ -10,78 +10,117 @@ jobs: test-ubuntu: name: Test Ubuntu runs-on: ubuntu-latest - + steps: - - uses: actions/checkout@v4 - - - name: Install dependencies - run: | - sudo apt-get update - sudo apt-get install -y cmake build-essential libomp-dev - - - name: Build SLEEF - run: | - git clone --branch 3.8 --depth 1 https://github.com/shibatch/sleef.git - cd sleef - cmake -S . -B build -DSLEEF_BUILD_QUAD=ON -DSLEEF_BUILD_SHARED_LIBS=ON - cmake --build build --parallel $(nproc) - sudo cmake --install build --prefix /usr/local - - - name: Build QuadBLAS - env: - SLEEF_ROOT: /usr/local - run: | - mkdir build && cd build - cmake .. - make -j$(nproc) - - - name: Run tests - run: | - cd build - ./quadblas_test + - uses: actions/checkout@v4 + with: + submodules: recursive + + - name: Install dependencies + run: | + sudo apt-get update + sudo apt-get install -y cmake build-essential libomp-dev + + - name: Build SLEEF + run: | + git clone --branch 3.9.0 --depth 1 https://github.com/shibatch/sleef.git + cd sleef + cmake -S . -B build \ + -DCMAKE_BUILD_TYPE=Release \ + -DSLEEF_BUILD_QUAD=ON \ + -DSLEEF_BUILD_SHARED_LIBS=ON \ + -DSLEEF_BUILD_TESTS=OFF \ + -DSLEEF_BUILD_INLINE_HEADERS=OFF \ + -DSLEEF_ENABLE_TLFLOAT=OFF \ + -DCMAKE_POSITION_INDEPENDENT_CODE=ON + cmake --build build --parallel $(nproc) + sudo cmake --install build --prefix /usr/local + + - name: Build QuadBLAS + env: + SLEEF_ROOT: /usr/local + run: | + mkdir build && cd build + cmake -DCMAKE_BUILD_TYPE=Release .. + make -j$(nproc) + + - name: Run tests + run: | + cd build + ./quadblas_test + + - name: Run benchmarks + run: | + cd build + if [ -f "./quadblas_benchmark" ]; then + echo "Running Google Benchmark tests..." + ./quadblas_benchmark --benchmark_min_time=0.1s --benchmark_filter="BM_DotProduct|BM_AXPY|BM_VectorNorm" + else + echo "Google Benchmark not available, running legacy benchmark..." + ./quadblas_benchmark_legacy + fi test-macos: name: Test macOS ${{ matrix.os }} runs-on: ${{ matrix.os }} strategy: matrix: - os: [macos-13, macos-14] + os: [macos-14, macos-15] steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v4 + with: + submodules: recursive + + - name: Install dependencies + run: | + brew install cmake libomp - - name: Install dependencies - run: | - brew install cmake libomp + - name: Install SLEEF + env: + MACOSX_DEPLOYMENT_TARGET: "11.0" + run: | + git clone --branch 3.9.0 --depth 1 https://github.com/shibatch/sleef.git + cd sleef + cmake -S . -B build \ + -DCMAKE_BUILD_TYPE=Release \ + -DSLEEF_BUILD_QUAD=ON \ + -DSLEEF_BUILD_SHARED_LIBS=ON \ + -DSLEEF_BUILD_TESTS=OFF \ + -DSLEEF_BUILD_INLINE_HEADERS=OFF \ + -DSLEEF_ENABLE_TLFLOAT=OFF \ + -DCMAKE_POSITION_INDEPENDENT_CODE=ON \ + -DCMAKE_OSX_DEPLOYMENT_TARGET=11.0 \ + -DCMAKE_INSTALL_RPATH="@loader_path/../lib" \ + -DCMAKE_BUILD_WITH_INSTALL_RPATH=ON + cmake --build build --clean-first -j + sudo cmake --install build --prefix /usr/local - - name: Install SLEEF - env: - MACOSX_DEPLOYMENT_TARGET: "11.0" - run: | - git clone --branch 3.8 --depth 1 https://github.com/shibatch/sleef.git - cd sleef - cmake -S . -B build \ - -DSLEEF_BUILD_QUAD=ON \ - -DSLEEF_BUILD_SHARED_LIBS=ON \ - -DCMAKE_POSITION_INDEPENDENT_CODE=ON \ - -DCMAKE_OSX_DEPLOYMENT_TARGET=11.0 \ - -DCMAKE_INSTALL_RPATH="@loader_path/../lib" \ - -DCMAKE_BUILD_WITH_INSTALL_RPATH=ON - cmake --build build --clean-first -j - sudo cmake --install build --prefix /usr/local + - name: Build QuadBLAS + env: + SLEEF_ROOT: /usr/local + DYLD_LIBRARY_PATH: "/usr/local/lib:$DYLD_LIBRARY_PATH" + run: | + mkdir build && cd build + cmake -DCMAKE_BUILD_TYPE=Release .. + make -j$(sysctl -n hw.ncpu) - - name: Build QuadBLAS - env: - SLEEF_ROOT: /usr/local - DYLD_LIBRARY_PATH: "/usr/local/lib:$DYLD_LIBRARY_PATH" - run: | - mkdir build && cd build - cmake .. - make -j$(sysctl -n hw.ncpu) + - name: Run tests + env: + DYLD_LIBRARY_PATH: "/usr/local/lib:$DYLD_LIBRARY_PATH" + run: | + cd build + ./quadblas_test - - name: Run tests - env: - DYLD_LIBRARY_PATH: "/usr/local/lib:$DYLD_LIBRARY_PATH" - run: | - cd build - ./quadblas_test + - name: Run benchmarks + env: + DYLD_LIBRARY_PATH: "/usr/local/lib:$DYLD_LIBRARY_PATH" + run: | + cd build + if [ -f "./quadblas_benchmark" ]; then + echo "Running Google Benchmark tests..." + ./quadblas_benchmark --benchmark_min_time=0.1s --benchmark_filter="BM_DotProduct|BM_AXPY|BM_VectorNorm" + else + echo "Google Benchmark not available, running legacy benchmark..." + ./quadblas_benchmark_legacy + fi diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..55a47ce --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "third_party/benchmark"] + path = third_party/benchmark + url = https://github.com/google/benchmark.git diff --git a/CMakeLists.txt b/CMakeLists.txt index a8f51bd..462d3ea 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,111 +1,149 @@ cmake_minimum_required(VERSION 3.15) project(QuadBLAS) -# Check for custom SLEEF installation via environment variable -set(SLEEF_CUSTOM_ROOT "$ENV{SLEEF_ROOT}") +# Set build type to Release if not specified +if(NOT CMAKE_BUILD_TYPE) + set(CMAKE_BUILD_TYPE Release CACHE STRING "Build type" FORCE) +endif() + +# Set C++ standard globally +set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD_REQUIRED ON) +# Find OpenMP (your existing code) +if(APPLE) + find_library(OpenMP_LIBRARY + NAMES omp gomp iomp5 + PATHS /opt/homebrew/lib /usr/local/lib /opt/homebrew/opt/libomp/lib + NO_DEFAULT_PATH + ) + find_path(OpenMP_INCLUDE_DIR + NAMES omp.h + PATHS /opt/homebrew/include /usr/local/include /opt/homebrew/opt/libomp/include + NO_DEFAULT_PATH + ) + + if(OpenMP_LIBRARY AND OpenMP_INCLUDE_DIR) + set(OpenMP_FOUND TRUE) + add_library(OpenMP::OpenMP_CXX UNKNOWN IMPORTED) + set_target_properties(OpenMP::OpenMP_CXX PROPERTIES + IMPORTED_LOCATION "${OpenMP_LIBRARY}" + INTERFACE_INCLUDE_DIRECTORIES "${OpenMP_INCLUDE_DIR}" + INTERFACE_COMPILE_OPTIONS "-Xpreprocessor;-fopenmp" + ) + message(STATUS "Found OpenMP: ${OpenMP_LIBRARY}") + else() + set(OpenMP_FOUND FALSE) + message(WARNING "OpenMP not found. Install with: brew install libomp") + endif() +else() + find_package(OpenMP) +endif() + +# Handle SLEEF library (fixed for your environment) +set(SLEEF_CUSTOM_ROOT "$ENV{SLEEF_ROOT}") if(SLEEF_CUSTOM_ROOT) - # Use custom SLEEF installation - message(STATUS "Using custom SLEEF installation: ${SLEEF_CUSTOM_ROOT}") set(SLEEF_INCLUDE_DIR "${SLEEF_CUSTOM_ROOT}/include") set(SLEEF_LIB_DIR "${SLEEF_CUSTOM_ROOT}/lib") - - # Verify paths exist - if(NOT EXISTS ${SLEEF_INCLUDE_DIR}) - message(FATAL_ERROR "SLEEF include directory not found: ${SLEEF_INCLUDE_DIR}") - endif() - if(NOT EXISTS ${SLEEF_LIB_DIR}) - message(FATAL_ERROR "SLEEF library directory not found: ${SLEEF_LIB_DIR}") + if(NOT EXISTS ${SLEEF_INCLUDE_DIR} OR NOT EXISTS ${SLEEF_LIB_DIR}) + message(FATAL_ERROR "SLEEF installation not found at: ${SLEEF_CUSTOM_ROOT}") endif() + find_library(SLEEF_LIBRARY + NAMES sleef + PATHS ${SLEEF_LIB_DIR} + NO_DEFAULT_PATH + ) + find_library(SLEEFQUAD_LIBRARY + NAMES sleefquad + PATHS ${SLEEF_LIB_DIR} + NO_DEFAULT_PATH + ) + + if(SLEEF_LIBRARY AND SLEEFQUAD_LIBRARY) + set(SLEEF_LIBRARIES ${SLEEF_LIBRARY} ${SLEEFQUAD_LIBRARY}) + message(STATUS "Using SLEEF from: ${SLEEF_CUSTOM_ROOT}") + message(STATUS "Found SLEEF libraries: ${SLEEF_LIBRARIES}") + else() + message(FATAL_ERROR "SLEEF libraries not found in ${SLEEF_LIB_DIR}") + endif() else() - # Use system SLEEF via pkg-config find_package(PkgConfig REQUIRED) pkg_check_modules(SLEEF REQUIRED sleef) endif() -# Find OpenMP -find_package(OpenMP) +# Add Google Benchmark with proper Release configuration +if(EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/third_party/benchmark/CMakeLists.txt") + # Configure Google Benchmark for Release mode and disable testing + set(BENCHMARK_ENABLE_TESTING OFF CACHE BOOL "Enable testing of the benchmark library.") + set(BENCHMARK_ENABLE_GTEST_TESTS OFF CACHE BOOL "Enable Google Test in benchmark.") + set(BENCHMARK_DOWNLOAD_DEPENDENCIES ON CACHE BOOL "Allow benchmark to download dependencies.") + + # Ensure benchmark is built in Release mode + set(CMAKE_BUILD_TYPE_BACKUP ${CMAKE_BUILD_TYPE}) + set(CMAKE_BUILD_TYPE Release) + + add_subdirectory(third_party/benchmark) + + # Restore original build type + set(CMAKE_BUILD_TYPE ${CMAKE_BUILD_TYPE_BACKUP}) + + set(GOOGLE_BENCHMARK_AVAILABLE TRUE) + message(STATUS "Google Benchmark added (Release mode)") +else() + set(GOOGLE_BENCHMARK_AVAILABLE FALSE) + message(STATUS "Google Benchmark not found (run: git submodule add https://github.com/google/benchmark.git third_party/benchmark)") +endif() # Create executables add_executable(quadblas_test test_quadblas.cpp) -add_executable(quadblas_benchmark quaddtype_tests.cpp) - -# Set C++ standard -target_compile_features(quadblas_test PRIVATE cxx_std_17) -target_compile_features(quadblas_benchmark PRIVATE cxx_std_17) +add_executable(quadblas_benchmark_legacy benchmarks/legacy_benchmark.cpp) -# Add compiler flags for optimization -target_compile_options(quadblas_test PRIVATE - -O3 -march=native -ffast-math - $<$:-g -fsanitize=address> -) -target_compile_options(quadblas_benchmark PRIVATE - -O3 -march=native -ffast-math - $<$:-g -fsanitize=address> -) +# Add Google Benchmark executable if available +if(GOOGLE_BENCHMARK_AVAILABLE) + add_executable(quadblas_benchmark benchmarks/benchmark.cpp) +endif() -# Include directories -if(SLEEF_CUSTOM_ROOT) - target_include_directories(quadblas_test PRIVATE - ${CMAKE_CURRENT_SOURCE_DIR} # For QuadBLAS headers - ${SLEEF_INCLUDE_DIR} - ) - target_include_directories(quadblas_benchmark PRIVATE - ${CMAKE_CURRENT_SOURCE_DIR} # For QuadBLAS headers - ${SLEEF_INCLUDE_DIR} - ) - - # Link directories and libraries for custom SLEEF - target_link_directories(quadblas_test PRIVATE ${SLEEF_LIB_DIR}) - target_link_directories(quadblas_benchmark PRIVATE ${SLEEF_LIB_DIR}) +# Common configuration function +function(configure_target target_name) + target_compile_features(${target_name} PRIVATE cxx_std_17) - target_link_libraries(quadblas_test - sleef sleefquad - $<$:OpenMP::OpenMP_CXX> - $<$:-fsanitize=address> - ) - target_link_libraries(quadblas_benchmark - sleef sleefquad - $<$:OpenMP::OpenMP_CXX> - $<$:-fsanitize=address> - ) -else() - # Use pkg-config provided paths and libraries - target_include_directories(quadblas_test PRIVATE - ${CMAKE_CURRENT_SOURCE_DIR} # For QuadBLAS headers - ${SLEEF_INCLUDE_DIRS} - ) - target_include_directories(quadblas_benchmark PRIVATE - ${CMAKE_CURRENT_SOURCE_DIR} # For QuadBLAS headers - ${SLEEF_INCLUDE_DIRS} + # Set optimization flags based on build type + target_compile_options(${target_name} PRIVATE + $<$:-O3 -march=native -ffast-math -DNDEBUG> + $<$:-g -O0 -fsanitize=address> + $<$:-O2 -g -march=native -ffast-math> ) - target_link_directories(quadblas_test PRIVATE ${SLEEF_LIBRARY_DIRS}) - target_link_directories(quadblas_benchmark PRIVATE ${SLEEF_LIBRARY_DIRS}) + target_include_directories(${target_name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) - target_compile_definitions(quadblas_test PRIVATE ${SLEEF_CFLAGS_OTHER}) - target_compile_definitions(quadblas_benchmark PRIVATE ${SLEEF_CFLAGS_OTHER}) + if(SLEEF_CUSTOM_ROOT) + target_include_directories(${target_name} PRIVATE ${SLEEF_INCLUDE_DIR}) + else() + target_include_directories(${target_name} PRIVATE ${SLEEF_INCLUDE_DIRS}) + target_link_directories(${target_name} PRIVATE ${SLEEF_LIBRARY_DIRS}) + target_compile_definitions(${target_name} PRIVATE ${SLEEF_CFLAGS_OTHER}) + endif() - target_link_libraries(quadblas_test - ${SLEEF_LIBRARIES} - $<$:OpenMP::OpenMP_CXX> - $<$:-fsanitize=address> - ) - target_link_libraries(quadblas_benchmark + target_link_libraries(${target_name} ${SLEEF_LIBRARIES} - $<$:OpenMP::OpenMP_CXX> - $<$:-fsanitize=address> + $<$:OpenMP::OpenMP_CXX> + $<$:-fsanitize=address> ) -endif() +endfunction() -# Print configuration info -message(STATUS "SLEEF configuration:") -if(SLEEF_CUSTOM_ROOT) - message(STATUS " Custom installation: ${SLEEF_CUSTOM_ROOT}") -else() - message(STATUS " System installation via pkg-config") - message(STATUS " Include dirs: ${SLEEF_INCLUDE_DIRS}") - message(STATUS " Libraries: ${SLEEF_LIBRARIES}") +# Configure targets +configure_target(quadblas_test) +configure_target(quadblas_benchmark_legacy) + +# Configure Google Benchmark target with proper linking +if(GOOGLE_BENCHMARK_AVAILABLE) + configure_target(quadblas_benchmark) + target_link_libraries(quadblas_benchmark benchmark::benchmark) endif() -message(STATUS "OpenMP found: ${OpenMP_FOUND}") \ No newline at end of file + +message(STATUS "=== Configuration ===") +message(STATUS "Build Type: ${CMAKE_BUILD_TYPE}") +message(STATUS "OpenMP: ${OpenMP_FOUND}") +message(STATUS "SLEEF: ${SLEEF_LIBRARIES}") +message(STATUS "Google Benchmark: ${GOOGLE_BENCHMARK_AVAILABLE}") \ No newline at end of file diff --git a/README.md b/README.md index 4e70e05..64d349a 100644 --- a/README.md +++ b/README.md @@ -2,8 +2,6 @@ QuadBLAS is a high-performance linear algebra library implementing BLAS-compliant routines for IEEE 754 quadruple precision (binary128) floating-point arithmetic. Built as a header-only templated library on top of the SLEEF vectorized mathematical library, QuadBLAS provides optimized implementations of fundamental linear algebra operations with significant performance improvements over naive implementations. - - ## Technical Features - **Quadruple Precision Arithmetic**: Full IEEE 754 binary128 (128-bit) floating-point support @@ -19,13 +17,13 @@ QuadBLAS is a high-performance linear algebra library implementing BLAS-complian QuadBLAS demonstrates substantial performance improvements over naive implementations through algorithmic optimizations and parallel execution: -| Operation Category | Problem Scale | Performance Improvement | Throughput | -|-------------------|---------------|------------------------|------------| -| Level 1 BLAS (DDOT) | 10^5 elements | 21× over serial implementation | Vectorized execution | -| Level 2 BLAS (DGEMV) | 1500×1500 matrix | 75× over serial implementation | 1.6 GFLOPS sustained | +| Operation Category | Problem Scale | Performance Improvement | Throughput | +| -------------------- | ------------------ | ------------------------------- | --------------------- | +| Level 1 BLAS (DDOT) | 10^5 elements | 21× over serial implementation | Vectorized execution | +| Level 2 BLAS (DGEMV) | 1500×1500 matrix | 75× over serial implementation | 1.6 GFLOPS sustained | | Level 3 BLAS (DGEMM) | 1000×1000 matrices | 2.8× over serial implementation | 0.06 GFLOPS sustained | -*Performance measurements conducted on multi-core x86-64 architecture with OpenMP threading enabled.* +_Performance measurements conducted on multi-core x86-64 architecture with OpenMP threading enabled._ ## Prerequisites @@ -36,6 +34,8 @@ QuadBLAS demonstrates substantial performance improvements over naive implementa ## Installation +We already provide a `build.sh` file, which can directly executes the steps detailed as follows: + For systems with SLEEF built from source or non-standard installation paths: ```bash @@ -46,12 +46,12 @@ export SLEEF_ROOT=/path/to/sleef/installation git clone https://github.com/your-org/QuadBLAS.git cd QuadBLAS mkdir build && cd build -cmake .. +cmake -DCMAKE_BUILD_TYPE=Release .. make -j$(nproc) # Execute test suite ./quadblas_test -./quadblas_benchmark +./quadblas_benchmark --benchmark_min_time=0.5s ``` ## Usage @@ -74,12 +74,12 @@ For direct integration into existing projects: int main() { using namespace QuadBLAS; - + // Instantiate quadruple precision containers const size_t n = 1000; DefaultVector<> x(n), y(n); DefaultMatrix<> A(n, n); - + // Initialize with quadruple precision values for (size_t i = 0; i < n; ++i) { x[i] = Sleef_cast_from_doubleq1(static_cast(i + 1)); @@ -88,12 +88,12 @@ int main() { A(i, j) = Sleef_cast_from_doubleq1((i == j) ? 2.0 : 0.1); } } - + // Execute optimized linear algebra operations Sleef_quad dot_result = x.dot(y); // Level 1 BLAS: DDOT Sleef_quad norm_result = x.norm(); // Level 1 BLAS: DNRM2 A.gemv(SLEEF_QUAD_C(1.0), x, SLEEF_QUAD_C(0.0), y); // Level 2 BLAS: DGEMV - + return 0; } ``` @@ -110,7 +110,7 @@ double result = quadblas_qdot(n, x_ptr, 1, y_ptr, 1); // y := alpha*A*x + beta*y quadblas_qgemv('R', 'N', m, n, 1.0, A_ptr, lda, x_ptr, 1, 0.0, y_ptr, 1); -// Level 3 BLAS: Matrix-matrix multiplication +// Level 3 BLAS: Matrix-matrix multiplication // C := alpha*A*B + beta*C quadblas_qgemm('R', 'N', 'N', m, n, k, 1.0, A_ptr, lda, B_ptr, ldb, 0.0, C_ptr, ldc); @@ -144,7 +144,7 @@ QuadBLAS::DefaultVector<> x(3), y(3); // Configure test case prone to catastrophic cancellation in double precision x[0] = Sleef_cast_from_doubleq1(1e20); // Large positive value -x[1] = Sleef_cast_from_doubleq1(1.0); // Unit value +x[1] = Sleef_cast_from_doubleq1(1.0); // Unit value x[2] = Sleef_cast_from_doubleq1(-1e20); // Large negative value y[0] = y[1] = y[2] = Sleef_cast_from_doubleq1(1.0); @@ -157,98 +157,6 @@ std::cout << "Double precision computation: " << double_precision_result << std: std::cout << "Quadruple precision computation: " << Sleef_cast_to_doubleq1(quad_precision_result) << std::endl; ``` -## Algorithm Implementation - -QuadBLAS implements a hierarchical BLAS architecture with optimizations at multiple levels: - -### Level 1 BLAS (Vector-Vector Operations) -- **DDOT**: Vectorized dot product with parallel reduction -- **DNRM2**: Euclidean norm computation with overflow protection -- **DAXPY**: Scaled vector addition with threading - -### Level 2 BLAS (Matrix-Vector Operations) -- **DGEMV**: Matrix-vector multiplication with cache-aware row/column processing -- Support for both row-major and column-major memory layouts -- Adaptive parallelization based on problem size - -### Level 3 BLAS (Matrix-Matrix Operations) -- **DGEMM**: Matrix-matrix multiplication with multi-level cache blocking -- Panel-panel multiplication algorithms -- Memory hierarchy optimization through data packing - -### Optimization Techniques - -1. **SIMD Vectorization**: Exploitation of 2-way quadruple precision vector units -2. **Cache Blocking**: Multi-level tiling for L1, L2, and L3 cache optimization -3. **Parallel Execution**: OpenMP-based threading with work distribution -4. **Memory Alignment**: 32-byte aligned allocations for optimal SIMD performance - -## Supported Operations - -| BLAS Level | Routine | Description | Implementation Status | -|------------|---------|-------------|----------------------| -| 1 | DDOT | Vector dot product | Optimized | -| 1 | DNRM2 | Vector Euclidean norm | Optimized | -| 1 | DAXPY | Scaled vector addition | Optimized | -| 2 | DGEMV | Matrix-vector multiplication | Optimized | -| 3 | DGEMM | Matrix-matrix multiplication | Optimized | - -## Integration with Numerical Computing Environments - -### Python/NumPy Integration - -QuadBLAS provides C interface compatibility for integration with Python numerical computing environments: - -```python -# Cython interface definition -cdef extern from "include/quadblas/quadblas.hpp": - double quadblas_qdot(int n, void* x, int incx, void* y, int incy) - void quadblas_qgemv(char layout, char trans, int m, int n, double alpha, - void* A, int lda, void* x, int incx, - double beta, void* y, int incy) - -def optimized_dot_product(x, y): - """High-performance quadruple precision dot product.""" - return quadblas_qdot(len(x), x.data, 1, y.data, 1) - -def optimized_matrix_vector_multiply(A, x): - """High-performance quadruple precision matrix-vector multiplication.""" - y = np.zeros(A.shape[0], dtype=quad128) - quadblas_qgemv(b'R', b'N', A.shape[0], A.shape[1], 1.0, - A.data, A.shape[1], x.data, 1, 0.0, y.data, 1) - return y -``` - -## Building from Source - -### Build Configuration Options - -```bash -# Optimized release build (recommended for production) -cmake -DCMAKE_BUILD_TYPE=Release .. - -# Debug build with address sanitization -cmake -DCMAKE_BUILD_TYPE=Debug .. - -# Custom SLEEF library location -export SLEEF_ROOT=/custom/sleef/path -cmake .. -``` - -### Test Suite Execution - -```bash -# Execute functionality verification tests -./quadblas_test - -# Execute performance benchmark suite -./quadblas_benchmark - -# Configure OpenMP thread count -export OMP_NUM_THREADS=16 -./quadblas_benchmark -``` - ## Contributing We welcome contributions to QuadBLAS development. Please follow established contribution guidelines: @@ -270,6 +178,7 @@ We welcome contributions to QuadBLAS development. Please follow established cont - Integration support for additional programming languages ## Citation + ```bibtex @software{quadblas2025, title = {QuadBLAS: High-Performance Linear Algebra for Quadruple Precision Computing}, @@ -282,6 +191,7 @@ We welcome contributions to QuadBLAS development. Please follow established cont ``` ## Acknowledgments + ``` @software{sleef, title = {SLEEF: A Portable Vectorized Library of Elementary Functions}, diff --git a/benchmarks/benchmark.cpp b/benchmarks/benchmark.cpp new file mode 100644 index 0000000..7423834 --- /dev/null +++ b/benchmarks/benchmark.cpp @@ -0,0 +1,432 @@ +#include "include/quadblas/quadblas.hpp" +#include +#include +#include +#include +#include + +class QuadRandom +{ +private: + std::mt19937 gen; + std::uniform_real_distribution dist; + +public: + QuadRandom(unsigned seed = 42) : gen(seed), dist(-1.0, 1.0) {} + + Sleef_quad next() + { + return Sleef_cast_from_doubleq1(dist(gen)); + } + + void fill_array(Sleef_quad *arr, size_t n) + { + for (size_t i = 0; i < n; ++i) + { + arr[i] = next(); + } + } +}; + +template +class AlignedArray +{ +private: + T *data_; + size_t size_; + +public: + AlignedArray(size_t size) : size_(size) + { + data_ = QuadBLAS::aligned_alloc(size); + if (!data_) + { + throw std::bad_alloc(); + } + } + + ~AlignedArray() + { + QuadBLAS::aligned_free(data_); + } + + T *data() { return data_; } + const T *data() const { return data_; } + size_t size() const { return size_; } + + T &operator[](size_t i) { return data_[i]; } + const T &operator[](size_t i) const { return data_[i]; } +}; + +// Level 1 BLAS Benchmarks + +static void BM_DotProduct(benchmark::State &state) +{ + const size_t n = state.range(0); + + AlignedArray x(n), y(n); + QuadRandom rng; + rng.fill_array(x.data(), n); + rng.fill_array(y.data(), n); + + Sleef_quad result; + + for (auto _ : state) + { + result = QuadBLAS::dot(n, x.data(), 1, y.data(), 1); + benchmark::DoNotOptimize(result); + } + + state.SetItemsProcessed(state.iterations() * n); + state.SetBytesProcessed(state.iterations() * n * 2 * sizeof(Sleef_quad)); + state.counters["GFLOPS"] = benchmark::Counter( + 2.0 * n * state.iterations(), + benchmark::Counter::kIsRate, + benchmark::Counter::kIs1000); +} + +static void BM_VectorNorm(benchmark::State &state) +{ + const size_t n = state.range(0); + + QuadBLAS::DefaultVector<> x(n); + QuadRandom rng; + for (size_t i = 0; i < n; ++i) + { + x[i] = rng.next(); + } + + Sleef_quad result; + + for (auto _ : state) + { + result = x.norm(); + benchmark::DoNotOptimize(result); + } + + state.SetItemsProcessed(state.iterations() * n); + state.counters["GFLOPS"] = benchmark::Counter( + 2.0 * n * state.iterations(), + benchmark::Counter::kIsRate, + benchmark::Counter::kIs1000); +} + +// y = alpha * x + y +static void BM_AXPY(benchmark::State &state) +{ + const size_t n = state.range(0); + + QuadBLAS::DefaultVector<> x(n), y(n); + QuadRandom rng; + Sleef_quad alpha = rng.next(); + + for (size_t i = 0; i < n; ++i) + { + x[i] = rng.next(); + y[i] = rng.next(); + } + + for (auto _ : state) + { + y.axpy(alpha, x); + benchmark::DoNotOptimize(y.data()); + } + + state.SetItemsProcessed(state.iterations() * n); + state.counters["GFLOPS"] = benchmark::Counter( + 2.0 * n * state.iterations(), + benchmark::Counter::kIsRate, + benchmark::Counter::kIs1000); +} + +// Level 2 BLAS Benchmarks + +static void BM_GEMV(benchmark::State &state) +{ + const size_t n = state.range(0); + const size_t m = n; + + AlignedArray A(m * n), x(n), y(m); + QuadRandom rng; + + rng.fill_array(A.data(), m * n); + rng.fill_array(x.data(), n); + rng.fill_array(y.data(), m); + + Sleef_quad alpha = Sleef_cast_from_doubleq1(1.0); + Sleef_quad beta = Sleef_cast_from_doubleq1(0.0); + + for (auto _ : state) + { + QuadBLAS::gemv(QuadBLAS::Layout::RowMajor, m, n, alpha, + A.data(), n, x.data(), 1, beta, y.data(), 1); + benchmark::DoNotOptimize(y.data()); + } + + state.SetItemsProcessed(state.iterations() * m * n); + state.counters["GFLOPS"] = benchmark::Counter( + 2.0 * m * n * state.iterations(), + benchmark::Counter::kIsRate, + benchmark::Counter::kIs1000); + state.counters["MatrixSize"] = n; +} + +// Level 3 BLAS Benchmarks + +static void BM_GEMM(benchmark::State &state) +{ + const size_t n = state.range(0); + const size_t m = n, k = n; + + AlignedArray A(m * k), B(k * n), C(m * n); + QuadRandom rng; + + rng.fill_array(A.data(), m * k); + rng.fill_array(B.data(), k * n); + rng.fill_array(C.data(), m * n); + + Sleef_quad alpha = Sleef_cast_from_doubleq1(1.0); + Sleef_quad beta = Sleef_cast_from_doubleq1(0.0); + + for (auto _ : state) + { + QuadBLAS::gemm(QuadBLAS::Layout::RowMajor, m, n, k, alpha, + A.data(), k, B.data(), n, beta, C.data(), n); + benchmark::DoNotOptimize(C.data()); + } + + state.SetItemsProcessed(state.iterations() * m * n * k); + state.counters["GFLOPS"] = benchmark::Counter( + 2.0 * m * n * k * state.iterations(), + benchmark::Counter::kIsRate, + benchmark::Counter::kIs1000); + state.counters["MatrixSize"] = n; +} + +// Threading Scalability Benchmarks + +static void BM_ThreadingScalability(benchmark::State &state) +{ + const int num_threads = state.range(0); + const size_t n = 10000; + + QuadBLAS::set_num_threads(num_threads); + + AlignedArray x(n), y(n); + QuadRandom rng; + rng.fill_array(x.data(), n); + rng.fill_array(y.data(), n); + + Sleef_quad result; + + for (auto _ : state) + { + result = QuadBLAS::dot(n, x.data(), 1, y.data(), 1); + benchmark::DoNotOptimize(result); + } + + state.counters["Threads"] = num_threads; + state.counters["GFLOPS"] = benchmark::Counter( + 2.0 * n * state.iterations(), + benchmark::Counter::kIsRate, + benchmark::Counter::kIs1000); + + static double baseline_time = 0.0; + if (num_threads == 1) + { + baseline_time = state.counters["GFLOPS"]; + } + else if (baseline_time > 0.0) + { + double current_gflops = state.counters["GFLOPS"]; + double speedup = current_gflops / baseline_time; + double efficiency = speedup / num_threads * 100.0; + state.counters["Efficiency%"] = efficiency; + } +} + +// Numerical Precision Test +static void BM_NumericalPrecision(benchmark::State &state) +{ + QuadBLAS::DefaultVector<> x(3), y(3); + + x[0] = Sleef_cast_from_doubleq1(1e20); + x[1] = Sleef_cast_from_doubleq1(1.0); + x[2] = Sleef_cast_from_doubleq1(-1e20); + + y[0] = y[1] = y[2] = Sleef_cast_from_doubleq1(1.0); + + Sleef_quad result; + + for (auto _ : state) + { + result = x.dot(y); + benchmark::DoNotOptimize(result); + } + + double final_result = Sleef_cast_to_doubleq1(result); + state.counters["Result"] = final_result; + state.counters["Error"] = std::abs(final_result - 1.0); +} + +// Memory Layout Benchmarks + +static void BM_RowMajorGEMV(benchmark::State &state) +{ + const size_t n = state.range(0); + QuadBLAS::MatrixRowMajor A(n, n); + QuadBLAS::VectorRowMajor x(n), y(n); + + QuadRandom rng; + for (size_t i = 0; i < n; ++i) + { + x[i] = rng.next(); + for (size_t j = 0; j < n; ++j) + { + A(i, j) = rng.next(); + } + } + + for (auto _ : state) + { + A.gemv(Sleef_cast_from_doubleq1(1.0), x, Sleef_cast_from_doubleq1(0.0), y); + benchmark::DoNotOptimize(y.data()); + } + + state.counters["Layout"] = 0; + state.counters["GFLOPS"] = benchmark::Counter( + 2.0 * n * n * state.iterations(), + benchmark::Counter::kIsRate, + benchmark::Counter::kIs1000); +} + +static void BM_ColMajorGEMV(benchmark::State &state) +{ + const size_t n = state.range(0); + QuadBLAS::MatrixColMajor A(n, n); + QuadBLAS::VectorColMajor x(n), y(n); + + QuadRandom rng; + for (size_t i = 0; i < n; ++i) + { + x[i] = rng.next(); + for (size_t j = 0; j < n; ++j) + { + A(i, j) = rng.next(); + } + } + + for (auto _ : state) + { + A.gemv(Sleef_cast_from_doubleq1(1.0), x, Sleef_cast_from_doubleq1(0.0), y); + benchmark::DoNotOptimize(y.data()); + } + + state.counters["Layout"] = 1; + state.counters["GFLOPS"] = benchmark::Counter( + 2.0 * n * n * state.iterations(), + benchmark::Counter::kIsRate, + benchmark::Counter::kIs1000); +} + +// =============================================================================== +// REGISTER BENCHMARKS WITH PARAMETER RANGES +// =============================================================================== + +// Level 1 BLAS - Vector sizes from 1K to 1M +BENCHMARK(BM_DotProduct) + ->RangeMultiplier(10) + ->Range(1000, 1000000) + ->Unit(benchmark::kMicrosecond); + +BENCHMARK(BM_VectorNorm) + ->RangeMultiplier(4) + ->Range(1000, 64000) + ->Unit(benchmark::kMicrosecond); + +BENCHMARK(BM_AXPY) + ->RangeMultiplier(4) + ->Range(1000, 64000) + ->Unit(benchmark::kMicrosecond); + +// Level 2 BLAS - Matrix sizes from 100x100 to 2000x2000 +BENCHMARK(BM_GEMV) + ->RangeMultiplier(2) + ->Range(100, 1600) + ->Unit(benchmark::kMillisecond); + +BENCHMARK(BM_RowMajorGEMV) + ->RangeMultiplier(2) + ->Range(100, 800) + ->Unit(benchmark::kMillisecond); + +BENCHMARK(BM_ColMajorGEMV) + ->RangeMultiplier(2) + ->Range(100, 800) + ->Unit(benchmark::kMillisecond); + +// Level 3 BLAS - Matrix sizes from 64x64 to 512x512 (practical for quad precision) +BENCHMARK(BM_GEMM) + ->RangeMultiplier(2) + ->Range(64, 512) + ->Unit(benchmark::kMillisecond) + ->MinTime(2.0); // Ensure sufficient measurement time + +// Threading scalability - range set in main() to avoid static init issues +static void RegisterThreadingBenchmarks() +{ + int max_threads = std::min(16, QuadBLAS::get_num_threads()); + ::benchmark::RegisterBenchmark("BM_ThreadingScalability", BM_ThreadingScalability) + ->DenseRange(1, max_threads, 1) + ->Unit(benchmark::kMicrosecond); +} + +// Numerical precision test +BENCHMARK(BM_NumericalPrecision) + ->Iterations(1000) + ->Unit(benchmark::kNanosecond); + +// Custom main function with QuadBLAS info +int main(int argc, char **argv) +{ + // Print system information + std::cout << "=== QuadBLAS Google Benchmark Suite ===" << std::endl; + std::cout << "Library Version: " << QuadBLAS::VERSION << std::endl; + std::cout << "Available Threads: " << QuadBLAS::get_num_threads() << std::endl; + +#ifdef QUADBLAS_X86_64 + std::cout << "Platform: x86-64 with SIMD vectorization" << std::endl; +#elif defined(QUADBLAS_AARCH64) + std::cout << "Platform: AArch64 with NEON vectorization" << std::endl; +#else + std::cout << "Platform: Generic (scalar fallback)" << std::endl; +#endif + +#ifdef _OPENMP + std::cout << "OpenMP: Enabled" << std::endl; +#else + std::cout << "OpenMP: Disabled" << std::endl; +#endif + + std::cout << "SIMD Vector Size: " << QuadBLAS::VECTOR_SIZE << " elements" << std::endl; + std::cout << "Memory Alignment: " << QuadBLAS::ALIGNMENT << " bytes" << std::endl; + + // Add macOS-specific system info +#ifdef __APPLE__ + std::cout << "Note: macOS warnings about CPU frequency/thread affinity are normal and don't affect timing accuracy" << std::endl; +#endif + + std::cout << "=========================================" << std::endl; + + // Register benchmarks that need runtime info + RegisterThreadingBenchmarks(); + + // Initialize and run benchmarks + ::benchmark::Initialize(&argc, argv); + if (::benchmark::ReportUnrecognizedArguments(argc, argv)) + return 1; + ::benchmark::RunSpecifiedBenchmarks(); + ::benchmark::Shutdown(); + + return 0; +} \ No newline at end of file diff --git a/benchmark.cpp b/benchmarks/legacy_benchmark.cpp similarity index 64% rename from benchmark.cpp rename to benchmarks/legacy_benchmark.cpp index 560e238..a348178 100644 --- a/benchmark.cpp +++ b/benchmarks/legacy_benchmark.cpp @@ -15,39 +15,72 @@ Sleef_quad from_double(double d) { return Sleef_cast_from_doubleq1(d); } -// High-precision timing -class PrecisionTimer { +// High-precision timing with warmup and averaging +class BenchTimer { private: std::chrono::high_resolution_clock::time_point start_time; public: void start() { start_time = std::chrono::high_resolution_clock::now(); } - + double stop_ms() { auto end_time = std::chrono::high_resolution_clock::now(); auto duration = std::chrono::duration_cast(end_time - start_time); - return duration.count() / 1e6; // Convert to milliseconds + return duration.count() / 1e6; } - + double stop_s() { return stop_ms() / 1000.0; } }; +// Compiler barrier to prevent dead code elimination of results +template +inline void do_not_optimize(const T &value) { +#if defined(__GNUC__) || defined(__clang__) + asm volatile("" : : "r,m"(value) : "memory"); +#else + volatile T sink = value; + (void)sink; +#endif +} + +// Run a callable multiple times and return median time +template +double benchmark_median(Func &&fn, int warmup = 2, int trials = 5) { + // Warmup + for (int i = 0; i < warmup; ++i) { + fn(); + } + + // Timed trials + std::vector times(trials); + BenchTimer timer; + for (int i = 0; i < trials; ++i) { + timer.start(); + fn(); + times[i] = timer.stop_ms(); + } + + // Return median + std::sort(times.begin(), times.end()); + return times[trials / 2]; +} + // Random number generator for consistent tests class QuadRandom { private: std::mt19937 gen; std::uniform_real_distribution dist; - + public: QuadRandom(unsigned seed = 12345) : gen(seed), dist(-1.0, 1.0) {} - + Sleef_quad next() { return from_double(dist(gen)); } - + void fill_array(Sleef_quad* arr, size_t n) { for (size_t i = 0; i < n; ++i) { arr[i] = next(); @@ -56,108 +89,92 @@ class QuadRandom { }; // ================================================================================ -// NAIVE IMPLEMENTATIONS (No optimizations) +// NAIVE IMPLEMENTATIONS (No optimizations - baseline for comparison) // ================================================================================ -// Naive dot product: O(n) with simple loop Sleef_quad naive_dot(size_t n, const Sleef_quad* x, const Sleef_quad* y) { Sleef_quad result = from_double(0.0); - for (size_t i = 0; i < n; ++i) { - result = Sleef_fmaq1_u05(x[i], y[i], result); // x[i]*y[i] + result + result = Sleef_fmaq1_u05(x[i], y[i], result); } - return result; } -// Naive matrix-vector: O(mn) with nested loops -void naive_gemv(size_t m, size_t n, Sleef_quad alpha, +void naive_gemv(size_t m, size_t n, Sleef_quad alpha, const Sleef_quad* A, size_t lda, const Sleef_quad* x, Sleef_quad beta, Sleef_quad* y) { - for (size_t i = 0; i < m; ++i) { Sleef_quad sum = from_double(0.0); - - // Compute dot product of row i with vector x for (size_t j = 0; j < n; ++j) { sum = Sleef_fmaq1_u05(A[i * lda + j], x[j], sum); } - - // y[i] = alpha * sum + beta * y[i] y[i] = Sleef_fmaq1_u05(alpha, sum, Sleef_mulq1_u05(beta, y[i])); } } -// Naive matrix-matrix: O(mnk) with triple nested loops void naive_gemm(size_t m, size_t n, size_t k, Sleef_quad alpha, const Sleef_quad* A, size_t lda, - const Sleef_quad* B, size_t ldb, + const Sleef_quad* B, size_t ldb, Sleef_quad beta, Sleef_quad* C, size_t ldc) { - for (size_t i = 0; i < m; ++i) { for (size_t j = 0; j < n; ++j) { Sleef_quad sum = from_double(0.0); - - // Compute dot product of row i of A with column j of B for (size_t l = 0; l < k; ++l) { sum = Sleef_fmaq1_u05(A[i * lda + l], B[l * ldb + j], sum); } - - // C[i,j] = alpha * sum + beta * C[i,j] C[i * ldc + j] = Sleef_fmaq1_u05(alpha, sum, Sleef_mulq1_u05(beta, C[i * ldc + j])); } } } // ================================================================================ -// BENCHMARK FUNCTIONS +// BENCHMARK FUNCTIONS (with warmup and median timing) // ================================================================================ void benchmark_dot_product() { std::cout << "\n" << std::string(80, '=') << "\n"; - std::cout << "DOT PRODUCT BENCHMARK (Vector • Vector)\n"; + std::cout << "DOT PRODUCT BENCHMARK (Vector . Vector)\n"; std::cout << std::string(80, '=') << "\n"; - + std::vector sizes = {1000, 5000, 10000, 50000, 100000}; - - std::cout << std::setw(10) << "Size" + + std::cout << std::setw(10) << "Size" << std::setw(15) << "Naive (ms)" - << std::setw(15) << "QuadBLAS (ms)" + << std::setw(15) << "QuadBLAS (ms)" << std::setw(12) << "Speedup" << std::setw(15) << "Error" << "\n"; std::cout << std::string(80, '-') << "\n"; - + QuadRandom rng; - + for (size_t n : sizes) { - // Allocate and initialize vectors Sleef_quad* x = QuadBLAS::aligned_alloc(n); Sleef_quad* y = QuadBLAS::aligned_alloc(n); - rng.fill_array(x, n); rng.fill_array(y, n); - - // Benchmark naive implementation - PrecisionTimer timer; - timer.start(); - Sleef_quad naive_result = naive_dot(n, x, y); - double naive_time = timer.stop_ms(); - - // Benchmark QuadBLAS implementation - timer.start(); - Sleef_quad quadblas_result = QuadBLAS::dot(n, x, 1, y, 1); - double quadblas_time = timer.stop_ms(); - + + Sleef_quad naive_result, quadblas_result; + + double naive_time = benchmark_median([&]() { + naive_result = naive_dot(n, x, y); + do_not_optimize(naive_result); + }); + + double quadblas_time = benchmark_median([&]() { + quadblas_result = QuadBLAS::dot(n, x, 1, y, 1); + do_not_optimize(quadblas_result); + }); + double speedup = naive_time / quadblas_time; double error = std::abs(to_double(naive_result) - to_double(quadblas_result)); - + std::cout << std::setw(10) << n << std::setw(15) << std::fixed << std::setprecision(3) << naive_time << std::setw(15) << std::fixed << std::setprecision(3) << quadblas_time << std::setw(12) << std::fixed << std::setprecision(2) << speedup << std::setw(15) << std::scientific << std::setprecision(2) << error << "\n"; - + QuadBLAS::aligned_free(x); QuadBLAS::aligned_free(y); } @@ -165,139 +182,116 @@ void benchmark_dot_product() { void benchmark_matrix_vector() { std::cout << "\n" << std::string(80, '=') << "\n"; - std::cout << "MATRIX-VECTOR BENCHMARK (Matrix × Vector)\n"; + std::cout << "MATRIX-VECTOR BENCHMARK (Matrix x Vector)\n"; std::cout << std::string(80, '=') << "\n"; - + std::vector sizes = {100, 300, 500, 1000, 1500}; - + std::cout << std::setw(10) << "Size" << std::setw(15) << "Naive (ms)" << std::setw(15) << "QuadBLAS (ms)" - << std::setw(12) << "Speedup" + << std::setw(12) << "Speedup" << std::setw(15) << "GFLOPS" << "\n"; std::cout << std::string(80, '-') << "\n"; - + QuadRandom rng; - + for (size_t n : sizes) { - size_t m = n; // Square matrices - - // Allocate matrices and vectors + size_t m = n; Sleef_quad* A = QuadBLAS::aligned_alloc(m * n); Sleef_quad* x = QuadBLAS::aligned_alloc(n); - Sleef_quad* y_naive = QuadBLAS::aligned_alloc(m); - Sleef_quad* y_quadblas = QuadBLAS::aligned_alloc(m); - + Sleef_quad* y = QuadBLAS::aligned_alloc(m); + rng.fill_array(A, m * n); rng.fill_array(x, n); - + Sleef_quad alpha = from_double(1.0); Sleef_quad beta = from_double(0.0); - - // Initialize y vectors - for (size_t i = 0; i < m; ++i) { - y_naive[i] = y_quadblas[i] = from_double(0.0); - } - - // Benchmark naive implementation - PrecisionTimer timer; - timer.start(); - naive_gemv(m, n, alpha, A, n, x, beta, y_naive); - double naive_time = timer.stop_ms(); - - // Reset y_quadblas - for (size_t i = 0; i < m; ++i) { - y_quadblas[i] = from_double(0.0); - } - - // Benchmark QuadBLAS implementation - timer.start(); - QuadBLAS::gemv(QuadBLAS::Layout::RowMajor, m, n, alpha, A, n, x, 1, beta, y_quadblas, 1); - double quadblas_time = timer.stop_ms(); - + + double naive_time = benchmark_median([&]() { + for (size_t i = 0; i < m; ++i) y[i] = from_double(0.0); + naive_gemv(m, n, alpha, A, n, x, beta, y); + }); + + double quadblas_time = benchmark_median([&]() { + for (size_t i = 0; i < m; ++i) y[i] = from_double(0.0); + QuadBLAS::gemv(QuadBLAS::Layout::RowMajor, m, n, alpha, A, n, x, 1, beta, y, 1); + }); + double speedup = naive_time / quadblas_time; - double gflops = (2.0 * m * n) / (quadblas_time * 1e6); // GFLOP/s - + double gflops = (2.0 * m * n) / (quadblas_time * 1e6); + std::cout << std::setw(10) << n << std::setw(15) << std::fixed << std::setprecision(3) << naive_time << std::setw(15) << std::fixed << std::setprecision(3) << quadblas_time << std::setw(12) << std::fixed << std::setprecision(2) << speedup << std::setw(15) << std::fixed << std::setprecision(3) << gflops << "\n"; - + QuadBLAS::aligned_free(A); QuadBLAS::aligned_free(x); - QuadBLAS::aligned_free(y_naive); - QuadBLAS::aligned_free(y_quadblas); + QuadBLAS::aligned_free(y); } } void benchmark_matrix_matrix() { std::cout << "\n" << std::string(80, '=') << "\n"; - std::cout << "MATRIX-MATRIX BENCHMARK (Matrix × Matrix)\n"; + std::cout << "MATRIX-MATRIX BENCHMARK (Matrix x Matrix)\n"; std::cout << std::string(80, '=') << "\n"; - - std::vector sizes = {100, 500, 1000, 2000, 4000}; - + + // Reduced to practical sizes for quad precision + // 500x500 GEMM = 250M quad FMAs which already takes seconds + std::vector sizes = {32, 64, 128, 200, 500}; + std::cout << std::setw(10) << "Size" << std::setw(15) << "Naive (s)" << std::setw(15) << "QuadBLAS (s)" << std::setw(12) << "Speedup" << std::setw(15) << "GFLOPS" << "\n"; std::cout << std::string(80, '-') << "\n"; - + QuadRandom rng; - + for (size_t n : sizes) { - size_t m = n, k = n; // Cubic matrices - - // Allocate matrices + size_t m = n, k = n; + Sleef_quad* A = QuadBLAS::aligned_alloc(m * k); Sleef_quad* B = QuadBLAS::aligned_alloc(k * n); - Sleef_quad* C_naive = QuadBLAS::aligned_alloc(m * n); - Sleef_quad* C_quadblas = QuadBLAS::aligned_alloc(m * n); - + Sleef_quad* C = QuadBLAS::aligned_alloc(m * n); + rng.fill_array(A, m * k); rng.fill_array(B, k * n); - + Sleef_quad alpha = from_double(1.0); Sleef_quad beta = from_double(0.0); - - // Initialize C matrices - for (size_t i = 0; i < m * n; ++i) { - C_naive[i] = C_quadblas[i] = from_double(0.0); - } - - // Benchmark naive implementation - PrecisionTimer timer; - timer.start(); - naive_gemm(m, n, k, alpha, A, k, B, n, beta, C_naive, n); - double naive_time = timer.stop_s(); - - // Reset C_quadblas - for (size_t i = 0; i < m * n; ++i) { - C_quadblas[i] = from_double(0.0); - } - - // Benchmark QuadBLAS implementation - timer.start(); - QuadBLAS::gemm(QuadBLAS::Layout::RowMajor, m, n, k, alpha, A, k, B, n, beta, C_quadblas, n); - double quadblas_time = timer.stop_s(); - + + // Fewer trials for large sizes + int trials = (n <= 128) ? 5 : 3; + int warmup = (n <= 128) ? 2 : 1; + + double naive_time = benchmark_median([&]() { + for (size_t i = 0; i < m * n; ++i) C[i] = from_double(0.0); + naive_gemm(m, n, k, alpha, A, k, B, n, beta, C, n); + }, warmup, trials) / 1000.0; // Convert to seconds + + double quadblas_time = benchmark_median([&]() { + for (size_t i = 0; i < m * n; ++i) C[i] = from_double(0.0); + QuadBLAS::gemm(QuadBLAS::Layout::RowMajor, m, n, k, alpha, A, k, B, n, beta, C, n); + }, warmup, trials) / 1000.0; + double speedup = naive_time / quadblas_time; - double gflops = (2.0 * m * n * k) / (quadblas_time * 1e9); // GFLOP/s - + double gflops = (2.0 * m * n * k) / (quadblas_time * 1e9); + std::cout << std::setw(10) << n << std::setw(15) << std::fixed << std::setprecision(3) << naive_time << std::setw(15) << std::fixed << std::setprecision(3) << quadblas_time << std::setw(12) << std::fixed << std::setprecision(2) << speedup << std::setw(15) << std::fixed << std::setprecision(3) << gflops << "\n"; - + QuadBLAS::aligned_free(A); QuadBLAS::aligned_free(B); - QuadBLAS::aligned_free(C_naive); - QuadBLAS::aligned_free(C_quadblas); + QuadBLAS::aligned_free(C); } } @@ -305,56 +299,52 @@ void benchmark_threading_scaling() { std::cout << "\n" << std::string(80, '=') << "\n"; std::cout << "THREADING SCALABILITY BENCHMARK\n"; std::cout << std::string(80, '=') << "\n"; - - const size_t n = 10000; // Large vector for threading test - - // Allocate vectors + + const size_t n = 50000; // Large enough for threading to matter + Sleef_quad* x = QuadBLAS::aligned_alloc(n); Sleef_quad* y = QuadBLAS::aligned_alloc(n); - + QuadRandom rng; rng.fill_array(x, n); rng.fill_array(y, n); - + std::cout << std::setw(12) << "Threads" << std::setw(15) << "Time (ms)" << std::setw(12) << "Speedup" << std::setw(15) << "Efficiency" << "\n"; std::cout << std::string(60, '-') << "\n"; - + double baseline_time = 0.0; int max_threads = QuadBLAS::get_num_threads(); - + for (int threads = 1; threads <= max_threads; threads *= 2) { QuadBLAS::set_num_threads(threads); - - PrecisionTimer timer; - const int iterations = 50; - - timer.start(); - for (int i = 0; i < iterations; ++i) { - volatile Sleef_quad result = QuadBLAS::dot(n, x, 1, y, 1); - (void)result; // Suppress unused variable warning - } - double avg_time = timer.stop_ms() / iterations; - + + Sleef_quad result; + double avg_time = benchmark_median([&]() { + result = QuadBLAS::dot(n, x, 1, y, 1); + do_not_optimize(result); + }, 3, 10); + + (void)result; + if (threads == 1) baseline_time = avg_time; - + double speedup = baseline_time / avg_time; double efficiency = speedup / threads * 100.0; - + std::cout << std::setw(12) << threads << std::setw(15) << std::fixed << std::setprecision(3) << avg_time << std::setw(12) << std::fixed << std::setprecision(2) << speedup << std::setw(14) << std::fixed << std::setprecision(1) << efficiency << "%" << "\n"; - - if (threads >= 16) break; // Don't test too many threads + + if (threads >= 16) break; } - - // Restore original thread count + QuadBLAS::set_num_threads(max_threads); - + QuadBLAS::aligned_free(x); QuadBLAS::aligned_free(y); } @@ -363,10 +353,10 @@ void print_system_info() { std::cout << "\n" << std::string(80, '=') << "\n"; std::cout << "SYSTEM INFORMATION\n"; std::cout << std::string(80, '=') << "\n"; - + std::cout << "QuadBLAS Version: " << quadblas_get_version() << "\n"; std::cout << "Available CPU threads: " << QuadBLAS::get_num_threads() << "\n"; - + #ifdef QUADBLAS_X86_64 std::cout << "Platform: x86-64 with SIMD vectorization\n"; #elif defined(QUADBLAS_AARCH64) @@ -380,36 +370,29 @@ void print_system_info() { #else std::cout << "OpenMP: Disabled\n"; #endif - + std::cout << "SIMD Vector Size: " << QuadBLAS::VECTOR_SIZE << " quad elements\n"; std::cout << "Memory Alignment: " << QuadBLAS::ALIGNMENT << " bytes\n"; - - // Test basic precision - Sleef_quad pi = SLEEF_M_PIq; - std::cout << "Quad precision π: " << std::setprecision(30) << to_double(pi) << "\n"; + std::cout << "Quad precision size: " << sizeof(Sleef_quad) << " bytes\n"; } int main() { print_system_info(); - + try { benchmark_dot_product(); benchmark_matrix_vector(); benchmark_matrix_matrix(); benchmark_threading_scaling(); - + std::cout << "\n" << std::string(80, '=') << "\n"; - std::cout << "BENCHMARK SUMMARY\n"; + std::cout << "All benchmarks completed.\n"; std::cout << std::string(80, '=') << "\n"; - std::cout << "✓ All benchmarks completed successfully!\n"; - std::cout << "✓ QuadBLAS demonstrates significant performance improvements\n"; - std::cout << "✓ Vectorization and threading optimizations are working\n"; - std::cout << "✓ Ready for production use with numpy_quaddtype\n"; - + } catch (const std::exception& e) { std::cerr << "Benchmark error: " << e.what() << std::endl; return 1; } - + return 0; -} \ No newline at end of file +} diff --git a/build.sh b/build.sh new file mode 100644 index 0000000..2a703fc --- /dev/null +++ b/build.sh @@ -0,0 +1,42 @@ +#!/bin/bash +set -e + +# Clean and setup +rm -rf build +mkdir build +cd build + +# Environment setup +export SLEEF_ROOT="/Users/swayam/Desktop/numpy_dtypes/sleef/build" + +# Detect available parallelism +if command -v nproc &> /dev/null; then + NJOBS=$(nproc) +elif command -v sysctl &> /dev/null; then + NJOBS=$(sysctl -n hw.ncpu) +else + NJOBS=4 +fi + +# Initialize submodules if needed +if [ ! -f "../third_party/benchmark/CMakeLists.txt" ]; then + echo "Initializing Google Benchmark submodule..." + (cd .. && git submodule update --init --recursive) +fi + +# Build +cmake -DCMAKE_BUILD_TYPE=Release .. +make -j"${NJOBS}" + +echo "=== Running Tests ===" +./quadblas_test + +# Google Benchmark +echo "" +echo "=== Running Google Benchmarks ===" +if [ -f "./quadblas_benchmark" ]; then + ./quadblas_benchmark --benchmark_min_time=0.5s +else + echo "Google Benchmark not available, running legacy..." + ./quadblas_benchmark_legacy +fi diff --git a/include/quadblas/algorithms/level1.hpp b/include/quadblas/algorithms/level1.hpp index 20cf902..fd0a5cd 100644 --- a/include/quadblas/algorithms/level1.hpp +++ b/include/quadblas/algorithms/level1.hpp @@ -10,25 +10,22 @@ namespace QuadBLAS { - // Vectorized dot product kernel for aligned data - inline Sleef_quad dot_kernel_vectorized(Sleef_quad *x, Sleef_quad *y, size_t n) + // Vectorized dot product kernel for contiguous data + inline Sleef_quad dot_kernel_vectorized(const Sleef_quad *x, const Sleef_quad *y, size_t n) { const size_t vec_n = n / VECTOR_SIZE; - const size_t remainder = n % VECTOR_SIZE; QuadVector sum_vec(SLEEF_QUAD_C(0.0)); - // Vectorized loop for (size_t i = 0; i < vec_n; ++i) { QuadVector x_vec = QuadVector::load(&x[i * VECTOR_SIZE]); QuadVector y_vec = QuadVector::load(&y[i * VECTOR_SIZE]); - sum_vec = x_vec.fma(y_vec, sum_vec); // x*y + sum + sum_vec = x_vec.fma(y_vec, sum_vec); } Sleef_quad result = sum_vec.horizontal_sum(); - // Handle remainder for (size_t i = vec_n * VECTOR_SIZE; i < n; ++i) { result = Sleef_fmaq1_u05(x[i], y[i], result); @@ -38,7 +35,7 @@ namespace QuadBLAS } // Parallel dot product for large vectors - inline Sleef_quad dot_parallel(Sleef_quad *x, Sleef_quad *y, size_t n) + inline Sleef_quad dot_parallel(const Sleef_quad *x, const Sleef_quad *y, size_t n) { if (n < PARALLEL_THRESHOLD) { @@ -49,7 +46,6 @@ namespace QuadBLAS const int num_threads = get_num_threads(); const size_t chunk_size = n / num_threads; - // Use array for partial results since OpenMP reduction doesn't work with struct types std::vector partial_results(num_threads); for (int i = 0; i < num_threads; ++i) { @@ -68,7 +64,6 @@ namespace QuadBLAS } } - // Combine partial results Sleef_quad result = SLEEF_QUAD_C(0.0); for (int i = 0; i < num_threads; ++i) { @@ -82,19 +77,17 @@ namespace QuadBLAS } // Main dot product function with stride support - inline Sleef_quad dot(size_t n, Sleef_quad *x, size_t incx, - Sleef_quad *y, size_t incy) + inline Sleef_quad dot(size_t n, const Sleef_quad *x, size_t incx, + const Sleef_quad *y, size_t incy) { if (n == 0) return SLEEF_QUAD_C(0.0); - // Fast path for unit strides if (incx == 1 && incy == 1) { return dot_parallel(x, y, n); } - // Strided access Sleef_quad result = SLEEF_QUAD_C(0.0); if (n >= PARALLEL_THRESHOLD) @@ -102,7 +95,6 @@ namespace QuadBLAS #ifdef _OPENMP const int num_threads = get_num_threads(); - // Use array for partial results std::vector partial_results(num_threads); for (int i = 0; i < num_threads; ++i) { @@ -122,7 +114,6 @@ namespace QuadBLAS } } - // Combine partial results for (int i = 0; i < num_threads; ++i) { result = Sleef_addq1_u05(result, partial_results[i]); @@ -145,6 +136,92 @@ namespace QuadBLAS return result; } + // Vectorized AXPY kernel (y = alpha * x + y) + inline void axpy_kernel_vectorized(Sleef_quad alpha, const Sleef_quad *x, Sleef_quad *y, size_t n) + { + const size_t vec_n = n / VECTOR_SIZE; + + QuadVector alpha_vec(alpha); + + for (size_t i = 0; i < vec_n; ++i) + { + QuadVector x_vec = QuadVector::load(&x[i * VECTOR_SIZE]); + QuadVector y_vec = QuadVector::load(&y[i * VECTOR_SIZE]); + QuadVector result_vec = x_vec.fma(alpha_vec, y_vec); + result_vec.store(&y[i * VECTOR_SIZE]); + } + + for (size_t i = vec_n * VECTOR_SIZE; i < n; ++i) + { + y[i] = Sleef_fmaq1_u05(alpha, x[i], y[i]); + } + } + + // Parallel AXPY for large vectors + inline void axpy_parallel(Sleef_quad alpha, const Sleef_quad *x, Sleef_quad *y, size_t n) + { + if (n < PARALLEL_THRESHOLD) + { + axpy_kernel_vectorized(alpha, x, y, n); + return; + } + +#ifdef _OPENMP + const int num_threads = get_num_threads(); + +#pragma omp parallel + { + int tid = omp_get_thread_num(); + size_t chunk_size = n / num_threads; + size_t start = tid * chunk_size; + size_t end = (tid == num_threads - 1) ? n : start + chunk_size; + + if (start < end) + { + axpy_kernel_vectorized(alpha, &x[start], &y[start], end - start); + } + } +#else + axpy_kernel_vectorized(alpha, x, y, n); +#endif + } + + // Main AXPY function with stride support + inline void axpy(size_t n, Sleef_quad alpha, const Sleef_quad *x, size_t incx, Sleef_quad *y, size_t incy) + { + if (n == 0) + return; + + if (incx == 1 && incy == 1) + { + axpy_parallel(alpha, x, y, n); + return; + } + + if (n >= PARALLEL_THRESHOLD) + { +#ifdef _OPENMP +#pragma omp parallel for + for (size_t i = 0; i < n; ++i) + { + y[i * incy] = Sleef_fmaq1_u05(alpha, x[i * incx], y[i * incy]); + } +#else + for (size_t i = 0; i < n; ++i) + { + y[i * incy] = Sleef_fmaq1_u05(alpha, x[i * incx], y[i * incy]); + } +#endif + } + else + { + for (size_t i = 0; i < n; ++i) + { + y[i * incy] = Sleef_fmaq1_u05(alpha, x[i * incx], y[i * incy]); + } + } + } + } // namespace QuadBLAS -#endif // QUADBLAS_ALGORITHMS_LEVEL1_HPP \ No newline at end of file +#endif // QUADBLAS_ALGORITHMS_LEVEL1_HPP diff --git a/include/quadblas/algorithms/level2.hpp b/include/quadblas/algorithms/level2.hpp index 844d04a..68c079d 100644 --- a/include/quadblas/algorithms/level2.hpp +++ b/include/quadblas/algorithms/level2.hpp @@ -13,8 +13,8 @@ namespace QuadBLAS // GEMV: y = alpha * A * x + beta * y // Row-major implementation inline void gemv_row_major(size_t m, size_t n, Sleef_quad alpha, - Sleef_quad *A, size_t lda, - Sleef_quad *x, size_t incx, + const Sleef_quad *A, size_t lda, + const Sleef_quad *x, size_t incx, Sleef_quad beta, Sleef_quad *y, size_t incy) { @@ -28,26 +28,22 @@ namespace QuadBLAS #endif for (size_t i = 0; i < m; ++i) { - // Compute dot product of row i with vector x Sleef_quad sum = SLEEF_QUAD_C(0.0); - Sleef_quad *row = &A[i * lda]; + const Sleef_quad *row = &A[i * lda]; - // Vectorized inner loop if (incx == 1) { sum = dot_kernel_vectorized(row, x, n); } else { - // Strided access for (size_t j = 0; j < n; ++j) { sum = Sleef_fmaq1_u05(row[j], x[j * incx], sum); } } - // y[i] = alpha * sum + beta * y[i] size_t y_idx = i * incy; y[y_idx] = Sleef_fmaq1_u05(alpha, sum, Sleef_mulq1_u05(beta, y[y_idx])); } @@ -87,8 +83,8 @@ namespace QuadBLAS // Main GEMV function inline void gemv(Layout layout, size_t m, size_t n, Sleef_quad alpha, - Sleef_quad *A, size_t lda, - Sleef_quad *x, size_t incx, + const Sleef_quad *A, size_t lda, + const Sleef_quad *x, size_t incx, Sleef_quad beta, Sleef_quad *y, size_t incy) { @@ -104,4 +100,4 @@ namespace QuadBLAS } // namespace QuadBLAS -#endif // QUADBLAS_ALGORITHMS_LEVEL2_HPP \ No newline at end of file +#endif // QUADBLAS_ALGORITHMS_LEVEL2_HPP diff --git a/include/quadblas/algorithms/level3.hpp b/include/quadblas/algorithms/level3.hpp index c55f915..b4b7bd1 100644 --- a/include/quadblas/algorithms/level3.hpp +++ b/include/quadblas/algorithms/level3.hpp @@ -13,9 +13,13 @@ namespace QuadBLAS { + // Micro-kernel tile sizes + constexpr size_t GEMM_MR = 4; + constexpr size_t GEMM_NR = 4; + inline void gemm_micro_kernel_scalar(size_t mr, size_t nr, size_t kc, Sleef_quad alpha, - Sleef_quad *A_packed, Sleef_quad *B_packed, + const Sleef_quad *A_packed, const Sleef_quad *B_packed, Sleef_quad beta, Sleef_quad *C, size_t ldc) { @@ -40,7 +44,7 @@ namespace QuadBLAS inline void gemm_micro_kernel_vectorized(size_t mr, size_t nr, size_t kc, Sleef_quad alpha, - Sleef_quad *A_packed, Sleef_quad *B_packed, + const Sleef_quad *A_packed, Sleef_quad *B_packed, Sleef_quad beta, Sleef_quad *C, size_t ldc) { @@ -50,10 +54,10 @@ namespace QuadBLAS return; } - const size_t mr_vec = mr / VECTOR_SIZE; const size_t nr_vec = nr / VECTOR_SIZE; - Sleef_quad c_acc[mr][nr]; + // Fixed-size accumulator (max GEMM_MR x GEMM_NR) - no VLA + Sleef_quad c_acc[GEMM_MR][GEMM_NR]; for (size_t i = 0; i < mr; ++i) { @@ -107,7 +111,7 @@ namespace QuadBLAS inline void gemm_micro_kernel(size_t mr, size_t nr, size_t kc, Sleef_quad alpha, - Sleef_quad *A_packed, Sleef_quad *B_packed, + const Sleef_quad *A_packed, Sleef_quad *B_packed, Sleef_quad beta, Sleef_quad *C, size_t ldc) { @@ -123,24 +127,24 @@ namespace QuadBLAS inline void gemm_macro_kernel(size_t mc, size_t nc, size_t kc, Sleef_quad alpha, - Sleef_quad *A_packed, Sleef_quad *B_packed, + const Sleef_quad *A_packed, Sleef_quad *B_packed, Sleef_quad beta, Sleef_quad *C, size_t ldc) { - constexpr size_t MR = 4; - constexpr size_t NR = 4; + // Allocate B_sub once outside the loop instead of per-tile + Sleef_quad *B_sub = aligned_alloc(kc * GEMM_NR); + const bool have_b_sub = (B_sub != nullptr); - for (size_t i = 0; i < mc; i += MR) + for (size_t i = 0; i < mc; i += GEMM_MR) { - size_t mr = std::min(MR, mc - i); + size_t mr = std::min(GEMM_MR, mc - i); - for (size_t j = 0; j < nc; j += NR) + for (size_t j = 0; j < nc; j += GEMM_NR) { - size_t nr = std::min(NR, nc - j); + size_t nr = std::min(GEMM_NR, nc - j); - Sleef_quad *B_sub = aligned_alloc(kc * nr); - if (B_sub) + if (have_b_sub) { - + // Repack B sub-tile into contiguous buffer for (size_t k = 0; k < kc; ++k) { for (size_t jj = 0; jj < nr; ++jj) @@ -152,24 +156,38 @@ namespace QuadBLAS gemm_micro_kernel(mr, nr, kc, alpha, &A_packed[i * kc], B_sub, beta, &C[i * ldc + j], ldc); - - aligned_free(B_sub); } else { - - gemm_micro_kernel_scalar(mr, nr, kc, alpha, - &A_packed[i * kc], &B_packed[j], - beta, &C[i * ldc + j], ldc); + // Fallback: work directly on B_packed with stride + for (size_t ii = 0; ii < mr; ++ii) + { + for (size_t jj = 0; jj < nr; ++jj) + { + Sleef_quad sum = SLEEF_QUAD_C(0.0); + for (size_t k = 0; k < kc; ++k) + { + sum = Sleef_fmaq1_u05(A_packed[(i + ii) * kc + k], + B_packed[k * nc + (j + jj)], sum); + } + Sleef_quad c_old = C[(i + ii) * ldc + (j + jj)]; + C[(i + ii) * ldc + (j + jj)] = Sleef_fmaq1_u05(alpha, sum, Sleef_mulq1_u05(beta, c_old)); + } + } } } } + + if (have_b_sub) + { + aligned_free(B_sub); + } } inline void gemm_simple(Layout layout, size_t m, size_t n, size_t k, Sleef_quad alpha, - Sleef_quad *A, size_t lda, - Sleef_quad *B, size_t ldb, + const Sleef_quad *A, size_t lda, + const Sleef_quad *B, size_t ldb, Sleef_quad beta, Sleef_quad *C, size_t ldc) { #ifdef _OPENMP @@ -196,8 +214,8 @@ namespace QuadBLAS inline void gemm(Layout layout, size_t m, size_t n, size_t k, Sleef_quad alpha, - Sleef_quad *A, size_t lda, - Sleef_quad *B, size_t ldb, + const Sleef_quad *A, size_t lda, + const Sleef_quad *B, size_t ldb, Sleef_quad beta, Sleef_quad *C, size_t ldc) { if (m == 0 || n == 0 || k == 0) @@ -227,29 +245,78 @@ namespace QuadBLAS { size_t kc = std::min(params.kc, k - kk); - for (size_t mm = 0; mm < m; mm += params.mc) + for (size_t nn = 0; nn < n; nn += params.nc) { - size_t mc = std::min(params.mc, m - mm); + size_t nc = std::min(params.nc, n - nn); - for (size_t i = 0; i < mc; ++i) + // Pack B panel once per (kk, nn) block + for (size_t i = 0; i < kc; ++i) { - for (size_t j = 0; j < kc; ++j) + for (size_t j = 0; j < nc; ++j) { - size_t src_idx = (layout == Layout::RowMajor) ? (mm + i) * lda + (kk + j) : (kk + j) * lda + (mm + i); - A_packed[i * kc + j] = A[src_idx]; + size_t src_idx = (layout == Layout::RowMajor) ? (kk + i) * ldb + (nn + j) : (nn + j) * ldb + (kk + i); + B_packed[i * nc + j] = B[src_idx]; } } - for (size_t nn = 0; nn < n; nn += params.nc) + // Iterate over row panels of A with OpenMP parallelism +#ifdef _OPENMP +#pragma omp parallel if (m >= PARALLEL_THRESHOLD) { - size_t nc = std::min(params.nc, n - nn); + // Each thread allocates its own A buffer once + Sleef_quad *A_local = nullptr; + bool multi_thread = (omp_get_num_threads() > 1); + if (multi_thread) + { + A_local = aligned_alloc(params.mc * kc); + } + else + { + A_local = A_packed; + } + + if (A_local) + { +#pragma omp for schedule(static) + for (size_t mm = 0; mm < m; mm += params.mc) + { + size_t mc = std::min(params.mc, m - mm); + + // Pack A panel + for (size_t i = 0; i < mc; ++i) + { + for (size_t j = 0; j < kc; ++j) + { + size_t src_idx = (layout == Layout::RowMajor) ? (mm + i) * lda + (kk + j) : (kk + j) * lda + (mm + i); + A_local[i * kc + j] = A[src_idx]; + } + } + + Sleef_quad *C_block = &C[(layout == Layout::RowMajor) ? mm * ldc + nn : nn * ldc + mm]; + + gemm_macro_kernel(mc, nc, kc, alpha, + A_local, B_packed, + (kk == 0) ? beta : SLEEF_QUAD_C(1.0), + C_block, ldc); + } + } - for (size_t i = 0; i < kc; ++i) + if (multi_thread && A_local) { - for (size_t j = 0; j < nc; ++j) + aligned_free(A_local); + } + } +#else + for (size_t mm = 0; mm < m; mm += params.mc) + { + size_t mc = std::min(params.mc, m - mm); + + for (size_t i = 0; i < mc; ++i) + { + for (size_t j = 0; j < kc; ++j) { - size_t src_idx = (layout == Layout::RowMajor) ? (kk + i) * ldb + (nn + j) : (nn + j) * ldb + (kk + i); - B_packed[i * nc + j] = B[src_idx]; + size_t src_idx = (layout == Layout::RowMajor) ? (mm + i) * lda + (kk + j) : (kk + j) * lda + (mm + i); + A_packed[i * kc + j] = A[src_idx]; } } @@ -260,6 +327,7 @@ namespace QuadBLAS (kk == 0) ? beta : SLEEF_QUAD_C(1.0), C_block, ldc); } +#endif } } @@ -269,4 +337,4 @@ namespace QuadBLAS } -#endif \ No newline at end of file +#endif diff --git a/include/quadblas/core/constants.hpp b/include/quadblas/core/constants.hpp index d424e28..45b21b6 100644 --- a/include/quadblas/core/constants.hpp +++ b/include/quadblas/core/constants.hpp @@ -11,7 +11,7 @@ namespace QuadBLAS constexpr size_t CACHE_LINE_SIZE = 64; constexpr size_t L1_CACHE_SIZE = 32768; constexpr size_t L2_CACHE_SIZE = 262144; - constexpr size_t PARALLEL_THRESHOLD = 1000; + constexpr size_t PARALLEL_THRESHOLD = 500; constexpr size_t GEMM_BLOCK_SIZE = 64; // Memory alignment for SIMD operations diff --git a/include/quadblas/detail/blocking.hpp b/include/quadblas/detail/blocking.hpp index c815d6f..c281e8e 100644 --- a/include/quadblas/detail/blocking.hpp +++ b/include/quadblas/detail/blocking.hpp @@ -8,44 +8,64 @@ namespace QuadBLAS { - // Cache-friendly blocking parameters + // Cache-friendly blocking parameters for GEMM + // + // Memory hierarchy usage (standard BLAS approach): + // L1: Micro-kernel slices: A(MR x kc) + B(kc x NR) + C(MR x NR) + // L2: A panel (mc x kc) stays resident + // L3/Memory: B panel (kc x nc) streams through struct BlockingParams { size_t mc, kc, nc; // Panel sizes for GEMM BlockingParams(size_t m = 0, size_t n = 0, size_t k = 0) { - // Calculate optimal blocking sizes based on cache hierarchy - const size_t quad_size = sizeof(Sleef_quad); - const size_t l1_quads = L1_CACHE_SIZE / (3 * quad_size); - const size_t l2_quads = L2_CACHE_SIZE / quad_size; - - if (m * n * k == 0) - { - // Default blocking - mc = std::min(static_cast(GEMM_BLOCK_SIZE), l1_quads / 4); - kc = std::min(static_cast(GEMM_BLOCK_SIZE), l1_quads / 4); - nc = std::min(static_cast(GEMM_BLOCK_SIZE), l2_quads / (mc + kc)); - } - else - { - // Adaptive blocking based on matrix sizes - mc = std::min({m, static_cast(GEMM_BLOCK_SIZE), l1_quads / 4}); - kc = std::min({k, static_cast(GEMM_BLOCK_SIZE), l1_quads / 4}); - nc = std::min({n, static_cast(GEMM_BLOCK_SIZE), l2_quads / (mc + kc)}); - } - - // Ensure sizes are multiples of vector size for better vectorization - mc = (mc / VECTOR_SIZE) * VECTOR_SIZE; - nc = (nc / VECTOR_SIZE) * VECTOR_SIZE; - - // Minimum sizes - mc = std::max(mc, VECTOR_SIZE); + constexpr size_t MR = 4; // Must match GEMM_MR + constexpr size_t NR = 4; // Must match GEMM_NR + const size_t elem_size = sizeof(Sleef_quad); + + // Determine cache sizes (detect at compile time for known platforms) + size_t l1_size = L1_CACHE_SIZE; + size_t l2_size = L2_CACHE_SIZE; + +#if defined(__APPLE__) && defined(QUADBLAS_AARCH64) + // Apple Silicon has larger caches + l1_size = 131072; // 128KB L1D per core + l2_size = 4194304; // 4MB L2 per cluster +#endif + + // kc: chosen so micro-kernel data fits in L1 + // Need: (MR * kc + kc * NR + MR * NR) * elem_size <= L1 / 2 + size_t l1_budget = l1_size / (2 * elem_size); + size_t max_kc = (l1_budget - MR * NR) / (MR + NR); + kc = std::min(max_kc, static_cast(256)); + + // mc: A panel (mc x kc) fits in L2 / 2 + size_t l2_budget = l2_size / (2 * elem_size); + size_t max_mc = l2_budget / kc; + mc = std::min(max_mc, static_cast(256)); + + // nc: B panel (kc x nc) can be larger - fits in remaining L2 or streams from L3 + size_t remaining_l2 = l2_budget - mc * kc; + size_t max_nc = remaining_l2 / kc; + nc = std::min(max_nc, static_cast(512)); + + // Clamp to actual matrix dimensions + if (m > 0) mc = std::min(mc, m); + if (k > 0) kc = std::min(kc, k); + if (n > 0) nc = std::min(nc, n); + + // Align to micro-kernel tile sizes for clean tiling + mc = (mc / MR) * MR; + nc = (nc / NR) * NR; + + // Ensure minimum sizes + mc = std::max(mc, MR); kc = std::max(kc, static_cast(4)); - nc = std::max(nc, VECTOR_SIZE); + nc = std::max(nc, NR); } }; } // namespace QuadBLAS -#endif // QUADBLAS_DETAIL_BLOCKING_HPP \ No newline at end of file +#endif // QUADBLAS_DETAIL_BLOCKING_HPP diff --git a/include/quadblas/interface/c_interface.hpp b/include/quadblas/interface/c_interface.hpp index 5088d40..b9dde00 100644 --- a/include/quadblas/interface/c_interface.hpp +++ b/include/quadblas/interface/c_interface.hpp @@ -7,16 +7,18 @@ #include "../algorithms/level2.hpp" #include "../algorithms/level3.hpp" #include "../threading/openmp_utils.hpp" +#include #include // C Interface for easy integration with Python/numpy +// All functions are inline to support header-only inclusion in multiple TUs extern "C" { // C API - LEVEL 1 BLAS (Vector operations) // QDOT: dot product of two vectors - double quadblas_qdot(int n, void *x, int incx, void *y, int incy) + inline double quadblas_qdot(int n, void *x, int incx, void *y, int incy) { Sleef_quad *qx = static_cast(x); Sleef_quad *qy = static_cast(y); @@ -29,7 +31,7 @@ extern "C" } // QNRM2: Euclidean norm of a vector - double quadblas_qnrm2(int n, void *x, int incx) + inline double quadblas_qnrm2(int n, void *x, int incx) { Sleef_quad *qx = static_cast(x); @@ -41,33 +43,25 @@ extern "C" return static_cast(Sleef_cast_to_doubleq1(result)); } - // QAXPY: y := alpha*x + y - void quadblas_qaxpy(int n, double alpha, void *x, int incx, void *y, int incy) + // QAXPY: y := alpha*x + y (uses optimized SIMD/threaded implementation) + inline void quadblas_qaxpy(int n, double alpha, void *x, int incx, void *y, int incy) { - Sleef_quad *qx = static_cast(x); - Sleef_quad *qy = static_cast(y); - Sleef_quad qalpha = Sleef_cast_from_doubleq1(alpha); - if (n <= 0) return; - size_t un = static_cast(n); - size_t uincx = static_cast(incx); - size_t uincy = static_cast(incy); + Sleef_quad *qx = static_cast(x); + Sleef_quad *qy = static_cast(y); + Sleef_quad qalpha = Sleef_cast_from_doubleq1(alpha); -#ifdef _OPENMP -#pragma omp parallel for if (un >= QuadBLAS::PARALLEL_THRESHOLD) -#endif - for (size_t i = 0; i < un; ++i) - { - qy[i * uincy] = Sleef_fmaq1_u05(qalpha, qx[i * uincx], qy[i * uincy]); - } + QuadBLAS::axpy(static_cast(n), qalpha, qx, + static_cast(incx), qy, + static_cast(incy)); } // C API - LEVEL 2 BLAS (Matrix-vector operations) // QGEMV: matrix-vector multiplication - void quadblas_qgemv(char layout, char trans, int m, int n, double alpha, + inline void quadblas_qgemv(char layout, char trans, int m, int n, double alpha, void *A, int lda, void *x, int incx, double beta, void *y, int incy) { @@ -81,18 +75,13 @@ extern "C" QuadBLAS::Layout quadblas_layout = (layout == 'C' || layout == 'c') ? QuadBLAS::Layout::ColMajor : QuadBLAS::Layout::RowMajor; - // Handle transpose (swap dimensions) + // Handle transpose (swap dimensions and flip layout) if (trans == 'T' || trans == 't' || trans == 'C' || trans == 'c') { std::swap(m, n); - if (quadblas_layout == QuadBLAS::Layout::RowMajor) - { - quadblas_layout = QuadBLAS::Layout::ColMajor; - } - else - { - quadblas_layout = QuadBLAS::Layout::RowMajor; - } + quadblas_layout = (quadblas_layout == QuadBLAS::Layout::RowMajor) + ? QuadBLAS::Layout::ColMajor + : QuadBLAS::Layout::RowMajor; } QuadBLAS::gemv(quadblas_layout, static_cast(m), static_cast(n), @@ -103,7 +92,7 @@ extern "C" // C API - LEVEL 3 BLAS (Matrix-matrix operations) // QGEMM: matrix-matrix multiplication - void quadblas_qgemm(char layout, char transa, char transb, int m, int n, int k, + inline void quadblas_qgemm(char layout, char transa, char transb, int m, int n, int k, double alpha, void *A, int lda, void *B, int ldb, double beta, void *C, int ldc) { @@ -117,16 +106,10 @@ extern "C" QuadBLAS::Layout quadblas_layout = (layout == 'C' || layout == 'c') ? QuadBLAS::Layout::ColMajor : QuadBLAS::Layout::RowMajor; - // For simplicity, we assume no transpose for now - // Full transpose support would require additional matrix manipulation - if (transa != 'N' && transa != 'n') - { - // TODO: Handle transpose A - } - if (transb != 'N' && transb != 'n') - { - // TODO: Handle transpose B - } + // Note: transpose support requires additional matrix manipulation + // For now, silently proceed with NoTranspose semantics + (void)transa; + (void)transb; QuadBLAS::gemm(quadblas_layout, static_cast(m), static_cast(n), static_cast(k), qalpha, qA, static_cast(lda), @@ -136,25 +119,25 @@ extern "C" // Utility functions // Set number of threads for OpenMP - void quadblas_set_num_threads(int num_threads) + inline void quadblas_set_num_threads(int num_threads) { QuadBLAS::set_num_threads(num_threads); } // Get number of threads - int quadblas_get_num_threads(void) + inline int quadblas_get_num_threads(void) { return QuadBLAS::get_num_threads(); } // Version information - const char *quadblas_get_version(void) + inline const char *quadblas_get_version(void) { return "QuadBLAS 1.0.0 - High Performance Quad Precision BLAS"; } // Memory alignment check - int quadblas_is_aligned(const void *ptr) + inline int quadblas_is_aligned(const void *ptr) { uintptr_t addr = reinterpret_cast(ptr); return (addr % QuadBLAS::ALIGNMENT) == 0; @@ -162,4 +145,4 @@ extern "C" } // extern "C" -#endif // QUADBLAS_INTERFACE_C_INTERFACE_HPP \ No newline at end of file +#endif // QUADBLAS_INTERFACE_C_INTERFACE_HPP diff --git a/include/quadblas/interface/cpp_classes.hpp b/include/quadblas/interface/cpp_classes.hpp index 91050cd..0a5b210 100644 --- a/include/quadblas/interface/cpp_classes.hpp +++ b/include/quadblas/interface/cpp_classes.hpp @@ -55,12 +55,12 @@ namespace QuadBLAS Vector &operator=(const Vector &) = delete; Sleef_quad &operator[](size_t i) { return data_[i * stride_]; } - Sleef_quad &operator[](size_t i) const { return data_[i * stride_]; } + const Sleef_quad &operator[](size_t i) const { return data_[i * stride_]; } size_t size() const { return size_; } size_t stride() const { return stride_; } Sleef_quad *data() { return data_; } - Sleef_quad *data() const { return data_; } + const Sleef_quad *data() const { return data_; } // Dot product Sleef_quad dot(const Vector &other) const @@ -71,10 +71,7 @@ namespace QuadBLAS // AXPY: this = alpha * other + this void axpy(Sleef_quad alpha, const Vector &other) { - for (size_t i = 0; i < size_; ++i) - { - (*this)[i] = Sleef_fmaq1_u05(alpha, other[i], (*this)[i]); - } + QuadBLAS::axpy(size_, alpha, other.data_, other.stride_, data_, stride_); } // Norm @@ -131,7 +128,7 @@ namespace QuadBLAS return layout == Layout::RowMajor ? data_[i * ld_ + j] : data_[j * ld_ + i]; } - Sleef_quad &operator()(size_t i, size_t j) const + const Sleef_quad &operator()(size_t i, size_t j) const { return layout == Layout::RowMajor ? data_[i * ld_ + j] : data_[j * ld_ + i]; } @@ -140,7 +137,7 @@ namespace QuadBLAS size_t cols() const { return cols_; } size_t leading_dimension() const { return ld_; } Sleef_quad *data() { return data_; } - Sleef_quad *data() const { return data_; } + const Sleef_quad *data() const { return data_; } // Matrix-vector multiplication void gemv(Sleef_quad alpha, const Vector &x, Sleef_quad beta, Vector &y) const diff --git a/include/quadblas/simd/quad_vector.hpp b/include/quadblas/simd/quad_vector.hpp index 4f9e74c..4e5fa4f 100644 --- a/include/quadblas/simd/quad_vector.hpp +++ b/include/quadblas/simd/quad_vector.hpp @@ -35,27 +35,18 @@ namespace QuadBLAS QuadVector(Sleef_quad a, Sleef_quad b) { - // Universal fallback: create via load from heap memory - Sleef_quad *temp = aligned_alloc(2); - if (temp) - { - temp[0] = a; - temp[1] = b; - *this = QuadVector::load(temp); // Use the working load function - aligned_free(temp); - } - else - { - // Emergency fallback if allocation fails + // Use stack-aligned buffer instead of heap allocation + alignas(ALIGNMENT) Sleef_quad temp[2]; + temp[0] = a; + temp[1] = b; #ifdef QUADBLAS_X86_64 - data_ = Sleef_splatq2_sse2(a); // At least get one value + data_ = Sleef_loadq2_sse2(temp); #elif defined(QUADBLAS_AARCH64) - data_ = Sleef_splatq2_advsimd(a); + data_ = Sleef_loadq2_advsimd(temp); #else - data_[0] = a; - data_[1] = a; // Not ideal but safe + data_[0] = a; + data_[1] = b; #endif - } } static QuadVector load(Sleef_quad *ptr) @@ -72,6 +63,12 @@ namespace QuadBLAS return result; } + // Const overload: SLEEF load takes non-const but doesn't modify memory + static QuadVector load(const Sleef_quad *ptr) + { + return load(const_cast(ptr)); + } + void store(Sleef_quad *ptr) const { #ifdef QUADBLAS_X86_64 diff --git a/perf_compare.cpp b/perf_compare.cpp new file mode 100644 index 0000000..a0992e2 --- /dev/null +++ b/perf_compare.cpp @@ -0,0 +1,185 @@ +#include "include/quadblas/quadblas.hpp" +#include +#include +#include +#include +#include +#include +#include + +static double to_double(Sleef_quad q) { + return static_cast(Sleef_cast_to_doubleq1(q)); +} + +static Sleef_quad from_double(double d) { + return Sleef_cast_from_doubleq1(d); +} + +// Compiler barrier to prevent dead-code elimination +template +static void do_not_optimize(T const &value) { + asm volatile("" : : "r,m"(value) : "memory"); +} + +static void fill_random(Sleef_quad *arr, size_t n, unsigned seed = 42) { + std::mt19937 gen(seed); + std::uniform_real_distribution dist(-1.0, 1.0); + for (size_t i = 0; i < n; ++i) + arr[i] = from_double(dist(gen)); +} + +// Returns median time in milliseconds over `trials` runs (with 2 warmup runs) +template +static double benchmark_median(Func fn, int trials = 7) { + // warmup + for (int i = 0; i < 2; ++i) fn(); + + std::vector times(trials); + for (int i = 0; i < trials; ++i) { + auto t0 = std::chrono::high_resolution_clock::now(); + fn(); + auto t1 = std::chrono::high_resolution_clock::now(); + times[i] = std::chrono::duration(t1 - t0).count(); + } + std::sort(times.begin(), times.end()); + return times[trials / 2]; // median +} + +// ---- Benchmarks ---- + +static void bench_dot() { + std::cout << "\n=== DOT PRODUCT ===\n"; + std::cout << std::setw(10) << "N" + << std::setw(15) << "Time (ms)" << "\n"; + std::cout << std::string(25, '-') << "\n"; + + for (size_t n : {1000, 5000, 10000, 50000}) { + Sleef_quad *x = QuadBLAS::aligned_alloc(n); + Sleef_quad *y = QuadBLAS::aligned_alloc(n); + fill_random(x, n, 1); + fill_random(y, n, 2); + + double ms = benchmark_median([&]() { + Sleef_quad r = QuadBLAS::dot(n, x, 1, y, 1); + do_not_optimize(r); + }); + + std::cout << std::setw(10) << n + << std::setw(15) << std::fixed << std::setprecision(3) << ms << "\n"; + + QuadBLAS::aligned_free(x); + QuadBLAS::aligned_free(y); + } +} + +static void bench_axpy() { + std::cout << "\n=== AXPY ===\n"; + std::cout << std::setw(10) << "N" + << std::setw(15) << "Time (ms)" << "\n"; + std::cout << std::string(25, '-') << "\n"; + + Sleef_quad alpha = from_double(2.5); + + for (size_t n : {1000, 5000, 10000, 50000}) { + Sleef_quad *x = QuadBLAS::aligned_alloc(n); + Sleef_quad *y = QuadBLAS::aligned_alloc(n); + Sleef_quad *y_backup = QuadBLAS::aligned_alloc(n); + fill_random(x, n, 3); + fill_random(y, n, 4); + std::copy(y, y + n, y_backup); + + double ms = benchmark_median([&]() { + std::copy(y_backup, y_backup + n, y); // restore y each run + QuadBLAS::axpy(n, alpha, x, 1, y, 1); + do_not_optimize(y[0]); + }); + + std::cout << std::setw(10) << n + << std::setw(15) << std::fixed << std::setprecision(3) << ms << "\n"; + + QuadBLAS::aligned_free(x); + QuadBLAS::aligned_free(y); + QuadBLAS::aligned_free(y_backup); + } +} + +static void bench_gemv() { + std::cout << "\n=== GEMV ===\n"; + std::cout << std::setw(10) << "N" + << std::setw(15) << "Time (ms)" << "\n"; + std::cout << std::string(25, '-') << "\n"; + + Sleef_quad alpha = from_double(1.0); + Sleef_quad beta = from_double(0.0); + + for (size_t n : {100, 300, 500, 1000}) { + Sleef_quad *A = QuadBLAS::aligned_alloc(n * n); + Sleef_quad *x = QuadBLAS::aligned_alloc(n); + Sleef_quad *y = QuadBLAS::aligned_alloc(n); + fill_random(A, n * n, 5); + fill_random(x, n, 6); + + double ms = benchmark_median([&]() { + for (size_t i = 0; i < n; ++i) + y[i] = from_double(0.0); + QuadBLAS::gemv(QuadBLAS::Layout::RowMajor, n, n, alpha, A, n, x, 1, beta, y, 1); + do_not_optimize(y[0]); + }); + + std::cout << std::setw(10) << n + << std::setw(15) << std::fixed << std::setprecision(3) << ms << "\n"; + + QuadBLAS::aligned_free(A); + QuadBLAS::aligned_free(x); + QuadBLAS::aligned_free(y); + } +} + +static void bench_gemm() { + std::cout << "\n=== GEMM ===\n"; + std::cout << std::setw(10) << "N" + << std::setw(15) << "Time (ms)" << "\n"; + std::cout << std::string(25, '-') << "\n"; + + Sleef_quad alpha = from_double(1.0); + Sleef_quad beta = from_double(0.0); + + for (size_t n : {32, 64, 128, 200, 300}) { + Sleef_quad *A = QuadBLAS::aligned_alloc(n * n); + Sleef_quad *B = QuadBLAS::aligned_alloc(n * n); + Sleef_quad *C = QuadBLAS::aligned_alloc(n * n); + fill_random(A, n * n, 7); + fill_random(B, n * n, 8); + + int trials = (n <= 128) ? 5 : 3; + + double ms = benchmark_median([&]() { + for (size_t i = 0; i < n * n; ++i) + C[i] = from_double(0.0); + QuadBLAS::gemm(QuadBLAS::Layout::RowMajor, n, n, n, alpha, A, n, B, n, beta, C, n); + do_not_optimize(C[0]); + }, trials); + + std::cout << std::setw(10) << n + << std::setw(15) << std::fixed << std::setprecision(3) << ms << "\n"; + + QuadBLAS::aligned_free(A); + QuadBLAS::aligned_free(B); + QuadBLAS::aligned_free(C); + } +} + +int main() { + std::cout << "QuadBLAS Performance Comparison Benchmark\n"; + std::cout << "==========================================\n"; + std::cout << "Version: " << quadblas_get_version() << "\n"; + std::cout << "Threads: " << QuadBLAS::get_num_threads() << "\n"; + + bench_dot(); + bench_axpy(); + bench_gemv(); + bench_gemm(); + + std::cout << "\nDone.\n"; + return 0; +} diff --git a/test.sh b/test.sh deleted file mode 100644 index c7f43e8..0000000 --- a/test.sh +++ /dev/null @@ -1,9 +0,0 @@ -rm -rf build -mkdir build -cd build - -export SLEEF_ROOT=/datadrive/guestuser/t-swsingh/temp/sleef/build -cmake .. -make -./quadblas_test -./quadblas_benchmark \ No newline at end of file diff --git a/third_party/benchmark b/third_party/benchmark new file mode 160000 index 0000000..7697796 --- /dev/null +++ b/third_party/benchmark @@ -0,0 +1 @@ +Subproject commit 7697796fe5a70edbc3394d0593ef82afeebb7ddf