Fix stride(0) ordering, bound incx, reject incx == 0 and deduce the CBLAS index type

.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 <experimental/linalg>\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*'

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=<paths>, -DBLA_VENDOR=<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 $<BUILD_INTERFACE:${_linalg_inc}>)
+  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
   $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
+  $<BUILD_INTERFACE:${CMAKE_CURRENT_BINARY_DIR}/include/experimental>
   $<INSTALL_INTERFACE:include>
+  $<INSTALL_INTERFACE:include/experimental>
 )
 
 ################################################################################

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,

cmake/LinAlgConfig.cmake.in

@@ -1,3 +1,5 @@
 @PACKAGE_INIT@
 
+include(CMakeFindDependencyMacro)
+find_dependency(mdspan)
 include("${CMAKE_CURRENT_LIST_DIR}/linalgTargets.cmake")

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 <type_traits>
+
 namespace MDSPAN_IMPL_STANDARD_NAMESPACE {
 namespace MDSPAN_IMPL_PROPOSED_NAMESPACE {
 inline namespace __p1673_version_0 {
@@ -31,7 +34,7 @@ template<class ElementType,
          class Layout,
          class Accessor,
          class Scalar>
-void linalg_scale_rank_1(
+void generic_scale_rank_1(
   const Scalar alpha,
   mdspan<ElementType, extents<IndexType, ext0>, Layout, Accessor> x)
 {
@@ -40,6 +43,98 @@ void linalg_scale_rank_1(
   }
 }
 
+template<class Scalar, class MdspanType>
+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<Scalar, value_type> &&
+    impl::is_blas_value_type_v<value_type>;
+
+  constexpr bool blas_layout =
+    impl::is_blas_layout_type_v<typename MdspanType::layout_type>;
+
+  constexpr bool blas_accessor =
+    impl::is_blas_accessor_type_v<typename MdspanType::accessor_type>;
+
+  return blas_value_type && blas_layout && blas_accessor;
+}
+
+// Return true if a BLAS routine could be called to scale the vector, false otherwise.
+template<class ElementType,
+	 class IndexType,
+         ::std::size_t ext0,
+         class Layout,
+         class Accessor,
+         class Scalar>
+bool try_blas_scale(
+  const Scalar alpha,
+  mdspan<ElementType, extents<IndexType, ext0>, 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<impl::cblas_index>::max()) {
+    auto incx = x.stride(0);
+    if (incx > std::numeric_limits<impl::cblas_index>::max()) {
+      return false;
+    }
+    // Strided layout mappings can have stride zero; the BLAS cannot.
+    if (incx == 0) {
+      return false;
+    }
+    if constexpr (std::is_same_v<ElementType, float>) {
+      cblas_sscal(n, alpha, x.data_handle(), incx);
+      return true;
+    } else if constexpr (std::is_same_v<ElementType, double>) {
+      cblas_dscal(n, alpha, x.data_handle(), incx);
+      return true;
+    } else if constexpr (std::is_same_v<ElementType, std::complex<float>>) {
+      if constexpr(std::is_convertible_v<Scalar, float>) {
+        cblas_csscal(n, alpha, x.data_handle(), incx);
+        return true;
+      }
+      else if constexpr (std::is_convertible_v<Scalar, std::complex<float>>) {
+        auto converted_alpha = static_cast<std::complex<float>>(alpha);
+        cblas_cscal(n, &converted_alpha, x.data_handle(), incx);
+        return true;
+      }
+    } else if constexpr (std::is_same_v<ElementType, std::complex<double>>) {
+      if constexpr (std::is_convertible_v<Scalar, double>) {
+        cblas_zdscal(n, alpha, x.data_handle(), incx);
+        return true;
+      }
+      else if constexpr (std::is_convertible_v<Scalar, std::complex<double>>) {
+        auto converted_alpha = static_cast<std::complex<double>>(alpha);
+        cblas_zscal(n, &converted_alpha, x.data_handle(), incx);
+        return true;
+      }
+    }
+  }
+  #endif
+  return false;
+}
+
+template<class ElementType,
+	 class IndexType,
+         ::std::size_t ext0,
+         class Layout,
+         class Accessor,
+         class Scalar>
+void linalg_scale_rank_1(
+  const Scalar alpha,
+  mdspan<ElementType, extents<IndexType, ext0>, Layout, Accessor> x)
+{
+  bool done = false;
+  if constexpr (maybe_can_blas_scale<Scalar, decltype(x)>()) {
+    done = try_blas_scale(alpha, x);
+  }
+  if (!done) {
+    generic_scale_rank_1(alpha, x);
+  }
+}
+
 template<class ElementType,
 	 class IndexType,
          ::std::size_t numRows,

include/experimental/__p1673_bits/blas_helpers.hpp

@@ -0,0 +1,103 @@
+//@HEADER
+// ************************************************************************
+//
+//                        Kokkos v. 4.0
+//       Copyright (2022) National Technology & Engineering
+//               Solutions of Sandia, LLC (NTESS).
+//
+// Under the terms of Contract DE-NA0003525 with NTESS,
+// the U.S. Government retains certain rights in this software.
+//
+// Part of Kokkos, under the Apache License v2.0 with LLVM Exceptions.
+// See https://kokkos.org/LICENSE for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+// ************************************************************************
+//@HEADER
+
+#ifndef LINALG_INCLUDE_EXPERIMENTAL___P1673_BITS_BLAS_HELPERS_HPP_
+#define LINALG_INCLUDE_EXPERIMENTAL___P1673_BITS_BLAS_HELPERS_HPP_
+
+#include "__p1673_bits/linalg_config.h"
+#include <complex>
+#include <mdspan/mdspan.hpp>
+#include <type_traits>
+#ifdef LINALG_HAS_CBLAS
+#ifndef LINALG_CBLAS_INCLUDE
+#define LINALG_CBLAS_INCLUDE <cblas.h>
+#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<class ValueType>
+constexpr bool is_blas_value_type_v =
+  std::is_same_v<ValueType, float> ||
+  std::is_same_v<ValueType, double> ||
+  std::is_same_v<ValueType, std::complex<float>> ||
+  std::is_same_v<ValueType, std::complex<double>>;
+
+// The padded layouts are class templates taking a size_t, so detecting them
+// takes a partial specialization rather than is_same_v.
+template<class Layout>
+constexpr bool is_padded_layout_v = false;
+
+template<std::size_t PaddingValue>
+constexpr bool is_padded_layout_v<layout_left_padded<PaddingValue>> = true;
+
+template<std::size_t PaddingValue>
+constexpr bool is_padded_layout_v<layout_right_padded<PaddingValue>> = true;
+
+template<class Layout>
+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<Layout, layout_left> ||
+  std::is_same_v<Layout, layout_right> ||
+  is_padded_layout_v<Layout> ||
+  std::is_same_v<Layout, layout_stride>;
+
+// 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<class Accessor>
+constexpr bool is_default_accessor_v = false;
+
+template<class ElementType>
+constexpr bool is_default_accessor_v<default_accessor<ElementType>> = true;
+
+template<class Accessor>
+constexpr bool is_blas_accessor_type_v =
+  is_default_accessor_v<Accessor>;
+
+#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<class R, class A, class... Rest> 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_

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