diff --git a/.github/workflows/ghaction.yml b/.github/workflows/ghaction.yml index 4fcc42da..a434ed70 100644 --- a/.github/workflows/ghaction.yml +++ b/.github/workflows/ghaction.yml @@ -53,3 +53,70 @@ jobs: - name: Install stdblas run: cmake --install stdblas-build + + # The default matrix never compiles the CBLAS path; this job does. + cblas: + runs-on: ubuntu-latest + strategy: + fail-fast: false + matrix: + cxx-standard: [17,20] + + steps: + - name: Install a BLAS with a C interface + run: sudo apt-get update && sudo apt-get install -y libopenblas-dev + + - name: Check Out mdspan + uses: actions/checkout@v4 + with: + repository: kokkos/mdspan + path: mdspan-src + + - name: create directories + run: cmake -E make_directory mdspan-build stdblas-build + + - name: Configure mdspan + run: cmake -S mdspan-src -B mdspan-build -DMDSPAN_CXX_STANDARD=${{matrix.cxx-standard}} -DCMAKE_BUILD_TYPE=RelWithDebInfo -DCMAKE_INSTALL_PREFIX=mdspan-install + + - name: Build mdspan + run: cmake --build mdspan-build -j 3 + + - name: Install mdspan + run: cmake --install mdspan-build + + - name: Check Out stdblas + uses: actions/checkout@v4 + with: + path: stdblas-src + + - name: Configure stdblas with CBLAS + run: cmake -S stdblas-src -B stdblas-build -Dmdspan_DIR=mdspan-install -DLINALG_CXX_STANDARD=${{matrix.cxx-standard}} -DLINALG_ENABLE_TESTS=On -DLINALG_ENABLE_BLAS=ON -DCMAKE_BUILD_TYPE=RelWithDebInfo + + - name: Confirm the CBLAS path was actually compiled + run: grep -q '^#define LINALG_HAS_CBLAS' stdblas-build/include/experimental/__p1673_bits/linalg_config.h + + - name: Confirm experimental/linalg compiles on its own + run: | + printf '#include \nint main(){}\n' > selftest.cpp + g++ -std=c++${{matrix.cxx-standard}} -fsyntax-only \ + -I stdblas-src/include -I stdblas-build/include/experimental \ + -I mdspan-install/include selftest.cpp + + - name: Build stdblas + run: cmake --build stdblas-build -j 3 + + - name: Confirm the tests really link a CBLAS + run: nm -uC stdblas-build/tests/native/scale | grep -q cblas_dscal + + - name: Test stdblas + working-directory: stdblas-build + run: ctest --output-on-failure + + - name: Prove the header-name knob with a renamed header + run: | + mkdir renamed-include + cp /usr/include/x86_64-linux-gnu/cblas.h renamed-include/renamed_cblas.h + cmake -S stdblas-src -B knob-build -Dmdspan_DIR=mdspan-install -DLINALG_CXX_STANDARD=${{matrix.cxx-standard}} -DLINALG_ENABLE_TESTS=On -DLINALG_ENABLE_BLAS=ON -DLINALG_CBLAS_HEADER=renamed_cblas.h "-DLINALG_CBLAS_INCLUDE_DIRS=$PWD/renamed-include" + grep renamed_cblas.h knob-build/include/experimental/__p1673_bits/linalg_config.h + cmake --build knob-build -j 3 --target scale + ./knob-build/tests/native/scale --gtest_filter='*mdspan*' diff --git a/CMakeLists.txt b/CMakeLists.txt index 96977aea..15a468dc 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,5 +1,14 @@ cmake_minimum_required(VERSION 3.12) + +# Honor LINALG_* set() by a parent project before add_subdirectory. +if(POLICY CMP0077) + cmake_policy(SET CMP0077 NEW) +endif() +if(POLICY CMP0126) + cmake_policy(SET CMP0126 NEW) +endif() + project(LinAlg VERSION 0.0.1 LANGUAGES CXX @@ -117,9 +126,22 @@ if (NOT mdspan_FOUND) endif() find_package(BLAS) -option(LINALG_ENABLE_BLAS - "Assume that we are linking with a BLAS library." - ${BLAS_FOUND}) + +option(LINALG_ENABLE_BLAS "Call CBLAS from eligible algorithms." OFF) +if(LINALG_ENABLE_BLAS AND NOT BLAS_FOUND) + message(FATAL_ERROR + "LINALG_ENABLE_BLAS is ON, but find_package(BLAS) found no BLAS library. " + "Set -DBLAS_LIBRARIES=, -DBLA_VENDOR=, or -DLINALG_ENABLE_BLAS=OFF.") +endif() +set(LINALG_HAS_CBLAS ${LINALG_ENABLE_BLAS}) + +set(LINALG_CBLAS_HEADER "cblas.h" CACHE STRING + "Name of the CBLAS header, if the provider does not call it cblas.h.") +set(LINALG_CBLAS_INCLUDE_DIRS "" CACHE STRING + "Where to find the CBLAS header, as a ;-list, if not on the default include path.") +if(LINALG_HAS_CBLAS) + set(LINALG_CBLAS_INCLUDE "<${LINALG_CBLAS_HEADER}>") +endif() find_package(TBB) option(LINALG_ENABLE_TBB @@ -150,7 +172,20 @@ message(STATUS "Build include directory: ${CMAKE_CURRENT_BINARY_DIR}/include/exp add_library(linalg INTERFACE) add_library(std::linalg ALIAS linalg) -target_link_libraries(linalg INTERFACE mdspan) +if(LINALG_HAS_CBLAS) + target_link_libraries(linalg INTERFACE ${BLAS_LIBRARIES}) + foreach(_linalg_inc IN LISTS LINALG_CBLAS_INCLUDE_DIRS) + target_include_directories(linalg INTERFACE $) + endforeach() + unset(_linalg_inc) +endif() + +# FetchContent defines "mdspan"; an installed package exports "mdspan::mdspan". +if(TARGET mdspan::mdspan) + target_link_libraries(linalg INTERFACE mdspan::mdspan) +else() + target_link_libraries(linalg INTERFACE mdspan) +endif() if(LINALG_ENABLE_TBB) target_link_libraries(linalg INTERFACE TBB::tbb) @@ -166,7 +201,9 @@ endif() target_include_directories(linalg INTERFACE $ + $ $ + $ ) ################################################################################ diff --git a/README.md b/README.md index 21261893..4e61ee11 100644 --- a/README.md +++ b/README.md @@ -26,8 +26,11 @@ Other compilers, including MSVC 2019, have been tested in the past. 3. Run CMake, pointing it to your googletest and mdspan install locations - If you want to build tests, set LINALG_ENABLE_TESTS=ON - If you want to build examples, set LINALG_ENABLE_EXAMPLES=ON - - If you have a BLAS installation, set LINALG_ENABLE_BLAS=ON. - BLAS support is currently experimental. + - If you have a BLAS installation that provides CBLAS and want to call it, + set LINALG_ENABLE_BLAS=ON. It is off by default. If your provider does not + name its header cblas.h, set LINALG_CBLAS_HEADER (and, if it is not on the + default include path, LINALG_CBLAS_INCLUDE_DIRS). + CBLAS support is currently experimental. - If you have a TBB (Threading Building Blocks) installation and want to use TBB, set LINALG_ENABLE_TBB=ON (and optionally set TBB_DIR to the lib/cmake/TBB subdirectory of your TBB installation, diff --git a/cmake/LinAlgConfig.cmake.in b/cmake/LinAlgConfig.cmake.in index f36e6437..2bfdf8ca 100644 --- a/cmake/LinAlgConfig.cmake.in +++ b/cmake/LinAlgConfig.cmake.in @@ -1,3 +1,5 @@ @PACKAGE_INIT@ +include(CMakeFindDependencyMacro) +find_dependency(mdspan) include("${CMAKE_CURRENT_LIST_DIR}/linalgTargets.cmake") diff --git a/include/experimental/__p1673_bits/blas1_scale.hpp b/include/experimental/__p1673_bits/blas1_scale.hpp index dca5c2e2..64098d47 100644 --- a/include/experimental/__p1673_bits/blas1_scale.hpp +++ b/include/experimental/__p1673_bits/blas1_scale.hpp @@ -18,6 +18,9 @@ #ifndef LINALG_INCLUDE_EXPERIMENTAL___P1673_BITS_BLAS1_SCALE_HPP_ #define LINALG_INCLUDE_EXPERIMENTAL___P1673_BITS_BLAS1_SCALE_HPP_ +#include "blas_helpers.hpp" +#include + namespace MDSPAN_IMPL_STANDARD_NAMESPACE { namespace MDSPAN_IMPL_PROPOSED_NAMESPACE { inline namespace __p1673_version_0 { @@ -31,7 +34,7 @@ template -void linalg_scale_rank_1( +void generic_scale_rank_1( const Scalar alpha, mdspan, Layout, Accessor> x) { @@ -40,6 +43,98 @@ void linalg_scale_rank_1( } } +template +constexpr bool maybe_can_blas_scale() { + // It's OK if Scalar and value_type aren't the same, + // as long as we can convert Scalar to value_type. + using value_type = typename MdspanType::value_type; + constexpr bool blas_value_type = + std::is_convertible_v && + impl::is_blas_value_type_v; + + constexpr bool blas_layout = + impl::is_blas_layout_type_v; + + constexpr bool blas_accessor = + impl::is_blas_accessor_type_v; + + return blas_value_type && blas_layout && blas_accessor; +} + +// Return true if a BLAS routine could be called to scale the vector, false otherwise. +template +bool try_blas_scale( + const Scalar alpha, + mdspan, Layout, Accessor> x) +{ + #ifdef LINALG_HAS_CBLAS + auto n = x.extent(0); + // We can't call x.stride(0) until we know that x.is_strided() is true. + if (x.is_strided() && n <= std::numeric_limits::max()) { + auto incx = x.stride(0); + if (incx > std::numeric_limits::max()) { + return false; + } + // Strided layout mappings can have stride zero; the BLAS cannot. + if (incx == 0) { + return false; + } + if constexpr (std::is_same_v) { + cblas_sscal(n, alpha, x.data_handle(), incx); + return true; + } else if constexpr (std::is_same_v) { + cblas_dscal(n, alpha, x.data_handle(), incx); + return true; + } else if constexpr (std::is_same_v>) { + if constexpr(std::is_convertible_v) { + cblas_csscal(n, alpha, x.data_handle(), incx); + return true; + } + else if constexpr (std::is_convertible_v>) { + auto converted_alpha = static_cast>(alpha); + cblas_cscal(n, &converted_alpha, x.data_handle(), incx); + return true; + } + } else if constexpr (std::is_same_v>) { + if constexpr (std::is_convertible_v) { + cblas_zdscal(n, alpha, x.data_handle(), incx); + return true; + } + else if constexpr (std::is_convertible_v>) { + auto converted_alpha = static_cast>(alpha); + cblas_zscal(n, &converted_alpha, x.data_handle(), incx); + return true; + } + } + } + #endif + return false; +} + +template +void linalg_scale_rank_1( + const Scalar alpha, + mdspan, Layout, Accessor> x) +{ + bool done = false; + if constexpr (maybe_can_blas_scale()) { + done = try_blas_scale(alpha, x); + } + if (!done) { + generic_scale_rank_1(alpha, x); + } +} + template +#include +#include +#ifdef LINALG_HAS_CBLAS +#ifndef LINALG_CBLAS_INCLUDE +#define LINALG_CBLAS_INCLUDE +#endif +#include LINALG_CBLAS_INCLUDE +#endif + +namespace MDSPAN_IMPL_STANDARD_NAMESPACE { +namespace MDSPAN_IMPL_PROPOSED_NAMESPACE { +inline namespace __p1673_version_0 { +namespace linalg { +namespace impl { + +template +constexpr bool is_blas_value_type_v = + std::is_same_v || + std::is_same_v || + std::is_same_v> || + std::is_same_v>; + +// The padded layouts are class templates taking a size_t, so detecting them +// takes a partial specialization rather than is_same_v. +template +constexpr bool is_padded_layout_v = false; + +template +constexpr bool is_padded_layout_v> = true; + +template +constexpr bool is_padded_layout_v> = true; + +template +constexpr bool is_blas_layout_type_v = + // Assume that we have a C BLAS, which accepts + // both row-major and column-major layouts. + // + // This just means that the layouts COULD be valid. + // For layout_stride, we need to check the strides first. + std::is_same_v || + std::is_same_v || + is_padded_layout_v || + std::is_same_v; + +// The BLAS accepts accessors that deal with pointers to memory. +// default_accessor is a class template, so we can't just use is_same_v directly. +// +// scale doesn't accept conjugated_accessor or scaled_accessor +// because those are read-only accessors, and scale needs to +// write to the mdspan's elements. + +template +constexpr bool is_default_accessor_v = false; + +template +constexpr bool is_default_accessor_v> = true; + +template +constexpr bool is_blas_accessor_type_v = + is_default_accessor_v; + +#ifdef LINALG_HAS_CBLAS +// Deduce the BLAS integer index type from the first parameter +// (the length N) of cblas_dscal, as declared by the cblas.h in scope: +// int on an LP64 build, a 64-bit integer on an ILP64 build. +template A first_param(R (*)(A, Rest...)); +using cblas_index = decltype(first_param(cblas_dscal)); +#endif + +// We made the above queries traits, with their typical `_v` prefix. +// We make maybe_can_blas_scale() a function. +// That's a matter of taste; it could be a trait too. + +} // end namespace impl +} // end namespace linalg +} // end inline namespace __p1673_version_0 +} // end namespace MDSPAN_IMPL_PROPOSED_NAMESPACE +} // end namespace MDSPAN_IMPL_STANDARD_NAMESPACE + +#endif //LINALG_INCLUDE_EXPERIMENTAL___P1673_BITS_BLAS_HELPERS_HPP_ diff --git a/include/experimental/__p1673_bits/linalg_config.h.in b/include/experimental/__p1673_bits/linalg_config.h.in index 21ab176d..0eb63104 100644 --- a/include/experimental/__p1673_bits/linalg_config.h.in +++ b/include/experimental/__p1673_bits/linalg_config.h.in @@ -1,6 +1,8 @@ #pragma once #cmakedefine LINALG_ENABLE_ATOMIC_REF +#cmakedefine LINALG_HAS_CBLAS +#cmakedefine LINALG_CBLAS_INCLUDE @LINALG_CBLAS_INCLUDE@ #cmakedefine LINALG_ENABLE_BLAS #cmakedefine LINALG_ENABLE_CONCEPTS #cmakedefine LINALG_ENABLE_KOKKOS diff --git a/include/experimental/linalg b/include/experimental/linalg index c9b56e99..9d69d986 100644 --- a/include/experimental/linalg +++ b/include/experimental/linalg @@ -18,6 +18,8 @@ #pragma once #include "__p1673_bits/linalg_config.h" +// The headers below open MDSPAN_IMPL_STANDARD_NAMESPACE, which mdspan defines. +#include #include "__p1673_bits/macros.hpp" #include "__p1673_bits/linalg_execpolicy_mapper.hpp" #include "__p1673_bits/maybe_static_size.hpp" diff --git a/tests/native/CMakeLists.txt b/tests/native/CMakeLists.txt index efd0b4da..ac608486 100644 --- a/tests/native/CMakeLists.txt +++ b/tests/native/CMakeLists.txt @@ -1,14 +1,7 @@ macro(linalg_add_test name) add_executable(${name} ${name}.cpp) - if(BLAS_FOUND) - target_link_libraries(${name} linalg GTest::GTest GTest::Main ${BLAS_LIBRARIES}) - else() - # BLAS_LIBRARIES is literally "FALSE" if the BLAS was not found. - # Linking against that causes linker errors involving "FALSE.lib". - # Thus, we exclude BLAS_LIBRARIES completely if the BLAS was not found. - target_link_libraries(${name} linalg GTest::GTest GTest::Main) - endif() + target_link_libraries(${name} linalg GTest::GTest GTest::Main) add_test(${name} ${name}) endmacro() diff --git a/tests/native/scale.cpp b/tests/native/scale.cpp index 708e76dc..976171e4 100644 --- a/tests/native/scale.cpp +++ b/tests/native/scale.cpp @@ -77,6 +77,89 @@ namespace { } } } + + // Stride other than 1. Both halves matter: the referenced elements must be + // scaled and the ones between them must not, which is what a wrong incx gets + // wrong. + TEST(BLAS1_scale, mdspan_layout_stride) + { + using scalar_t = double; + using extents_t = extents; + using vector_t = mdspan; + + constexpr std::size_t vectorSize(5); + constexpr std::size_t stride(2); + std::vector storage(vectorSize * stride); + + layout_stride::mapping mapping{ + extents_t{vectorSize}, std::array{stride}}; + vector_t x(storage.data(), mapping); + + for (std::size_t k = 0; k < storage.size(); ++k) { + storage[k] = scalar_t (k) + 1.0; + } + const scalar_t scaleFactor = 5.0; + scale(scaleFactor, x); + for (std::size_t k = 0; k < storage.size(); ++k) { + const scalar_t storage_k = scalar_t (k) + 1.0; + const scalar_t expected = + (k % stride == 0) ? scaleFactor * storage_k : storage_k; + EXPECT_EQ( storage[k], expected ); + } + } + + TEST(BLAS1_scale, mdspan_layout_left_padded) + { + using scalar_t = double; + using extents_t = extents; + using layout_t = layout_left_padded<4>; + using vector_t = mdspan; + + constexpr std::size_t vectorSize(4); + constexpr std::size_t storageSize = vectorSize; + std::vector storage(storageSize); + + layout_t::mapping mapping{extents_t{}}; + vector_t x(storage.data(), mapping); + + for (std::size_t k = 0; k < vectorSize; ++k) { + const scalar_t x_k = scalar_t (k) + 1.0; + x(k) = x_k; + } + const scalar_t scaleFactor = 5.0; + scale(scaleFactor, x); + for (std::size_t k = 0; k < vectorSize; ++k) { + const scalar_t x_k = scalar_t (k) + 1.0; + EXPECT_EQ( x(k), scaleFactor * x_k ); + } + } + + TEST(BLAS1_scale, mdspan_layout_right_padded) + { + using scalar_t = double; + using extents_t = extents; + using layout_t = layout_right_padded; + using vector_t = mdspan; + + constexpr std::size_t vectorSize(4); + constexpr std::size_t paddingValue(4); + constexpr std::size_t storageSize = vectorSize; + std::vector storage(storageSize); + + layout_t::mapping mapping{extents_t{}, paddingValue}; + vector_t x(storage.data(), mapping); + + for (std::size_t k = 0; k < vectorSize; ++k) { + const scalar_t x_k = scalar_t (k) + 1.0; + x(k) = x_k; + } + const scalar_t scaleFactor = 5.0; + scale(scaleFactor, x); + for (std::size_t k = 0; k < vectorSize; ++k) { + const scalar_t x_k = scalar_t (k) + 1.0; + EXPECT_EQ( x(k), scaleFactor * x_k ); + } + } } // int main() {