wrapping a vendor dscal

[amklinv-nnl]

@georgemalerbo is interested in automatically generating BLAS wrappers for us. Thanks George!

Before we do that, @mhoemmen could you please provide feedback on the following dscal wrapper?

// Compiled with g++-16 -std=c++23 dscal_example.cpp -lblas
#include<cassert>
#include <vector>
#include <mdspan>
#include <iostream>

using namespace std;

// Should come from CMake or a config header
using blas_index_t = int; // BLAS typically uses int for sizes and indices

extern "C" {
    // We need to account for the Fortran name mangling, which often appends an underscore to the function name
    // The mangling should be automatically detected
    void dscal_(blas_index_t* n, double* a, double* x, blas_index_t* incx);
}

// Which accessors does it make sense to support?
template<class SizeType,
         ::std::size_t ext0,
         class Layout,
         class Accessor>
// scal requires strided data, so we can enforce this at compile time
// What is the difference between is_always_strided and is_strided? Should we support both?
requires ( mdspan<double, extents<SizeType, ext0>, Layout, Accessor>::is_always_strided() )
void scale(double a, mdspan<double, extents<SizeType, ext0>, Layout, Accessor> x) {
    blas_index_t n = x.extent(0);
    // This should be checked before calling the function
    // Like, in a "does it make sense to call BLAS?" function
    assert(n <= std::numeric_limits<blas_index_t>::max() && "Size of x exceeds maximum value for blas_index_t");

    // Assuming contiguous data, this will be 1 for a standard vector
    // This code will still work for non-contiguous data, but the performance may be worse due to cache misses
    blas_index_t incx = x.stride(0);
    dscal_(&n, &a, x.data_handle(), &incx);
}

int main() {
    vector<double> x_data = {1.0, 2.0, 3.0};
    auto x = mdspan(x_data.data(), x_data.size());

    double a = 2.0;
    scale(a, x);

    for(int i = 0; i < x.extent(0); ++i) {
        cout << x[i] << " ";
    }
    cout << endl;

    return 0;
}

[mhoemmen]

Thanks @georgemalerbo for helping out! : - )

using namespace std;

Putting this declaration inside main would avoid exposing it to the rest of the program.

    // We need to account for the Fortran name mangling, which often appends an underscore to the function name
    // The mangling should be automatically detected

I might suggest just starting with a C BLAS. OpenBLAS should have a C BLAS with which you can start. If you don't have a C BLAS, then the right encapsulation would be to provide a C BLAS wrapper around a Fortran BLAS. That's a separate problem that you can solve separately, later : - ) . This would also solve the "what is the index type?" problem: You'll have a C BLAS header, so you can examine the function types to make sure that the index type is int.

template<class SizeType,

We haven't gone through the reference std::linalg implementation yet to rename SizeType to IndexType, but we can fix that here : - ) .

    // This should be checked before calling the function
    // Like, in a "does it make sense to call BLAS?" function

Would you consider the following pattern?

template<class Scalar, class ElementType, class IndexType, std::size_t Ext0, class Layout, class Accessor>
void scale(Scalar a, mdspan<Scalar, extents<IndexType, Ext0>, Layout, Accessor> x) {
  bool done = false;
  if constexpr (maybe_can_blas_scale<Scalar, decltype(x)>()) { // compile-time check
    // Run-time check.
    // If it can call the BLAS, then call the BLAS and return true.
    // Otherwise, return false.
    done = try_blas_scale(a, x);
  }
  if (! done) {
    generic_scale(a, x);
  }
}

Christian Trott and I developed this pattern for kokkos-kernels. The idea is that the user-facing "outer" function is totally generic and works for all types that std::linalg supports. (That's why there are no extra "requires" clauses or other constraints on it.) The implementation of the "outer" function then uses first compile-time information, then run-time information to decide if it can call the BLAS. We fuse the run-time check with actually calling the BLAS (if the check passes) into one function, try_blas_scale. Putting the run-time check inside the if constexpr ensures that we only compile it if the types make sense.

[mhoemmen]

generic_scale is just generic straightforward code. It's the "fall-back" implementation in case we can't call the BLAS.

template<class Scalar, class ElementType, class IndexType, std::size_t Ext0, class Layout, class Accessor>
void generic_scale(Scalar a, mdspan<Scalar, extents<IndexType, Ext0>, Layout, Accessor> x) {
  std::ranges::for_each(std::views::indices(x.extent(0)), [=] (IndexType k) { x[k] *= a; });
  
  // Or, if you don't want to use ranges, you can do this.
  // IndexType e0 = x.extent(0);
  // for (IndexType k = 0; k < e0; ++k) {
  //   x[k] *= a; 
  // } 
}

[mhoemmen]

maybe_can_blas_scale only tests the template parameters of scale -- just the types. We combine all the mdspan-related template parameters of scale into a single MdspanType template parameter, just so that calling maybe_can_blas_scale looks as much like calling scale as possible (e.g., two template parameters vs. two function parameters). maybe_can_blas_scale takes template parameters instead of function parameters because that makes sure that everything is known at compile time. Another option would be to make this a trait (exposed as a variable template). Adding constexpr to the function declaration tells the compiler that this function can be called in a "constant expression" (at compile time).

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>>;

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, std::layout_left> ||
  std::is_same_v<Layout, std::layout_left_padded> ||
  std::is_same_v<Layout, std::layout_right> ||
  std::is_same_v<Layout, std::layout_right_padded> ||
  std::is_same_v<Layout, std::layout_stride>;
}

// The BLAS accepts accessors that deal with pointers to memory:
// default_accessor and aligned_accessor.
// Those are both templated classes, 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<std::default_accessor<ElementType>> = true;

template<class Accessor>
constexpr bool is_aligned_accessor_v = false;

template<class ElementType, std::size_t ByteAlignment>
constexpr bool is_aligned_accessor_v<std::aligned_accessor<ElementType, ByteAlignment>> = true;

template<class Accessor>
constexpr bool is_blas_accessor_type_v =
  is_default_accessor_v<Accessor> ||
  is_aligned_accessor_v<Accessor>;

// 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.

template<class Scalar, class MdspanType>
constexpr void 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> &&
    is_blas_value_type_v<value_type>();

  constexpr bool blas_layout =
    is_blas_layout_type_v<typename MdspanType::layout_type>;

  constexpr bool blas_accessor =
    is_blas_accessor_v<typename MdspanType::accessor_type>;

  return blas_value_type && blas_layout && blas_accessor;
}  

[amklinv-nnl]

provide a C BLAS wrapper around a Fortran BLAS. That's a separate problem that you can solve separately, later

My OCD is sufficiently satisfied with that response.

We haven't gone through the reference std::linalg implementation yet to rename SizeType to IndexType, but we can fix that here.

I was just about to do that anyway!

Would you consider the following pattern?

Yup! It almost sounds like y'all already put a lot of thought into it!

Or, if you don't want to use ranges, you can do this.

What is the oldest standard version the reference implementation has to support? Just want to know if I need to guard ranges with an ifdef.

Thanks so much for the detailed explanation!

[mhoemmen]

@amklinv-nnl wrote:

What is the oldest standard version the reference implementation has to support? Just want to know if I need to guard ranges with an ifdef.

Ah, sorry, I forgot that we should still support C++17. For loops it is, then.

Thanks for working on this! : - D

[mhoemmen]

@amklinv-nnl Btw, the implementation uses a slightly different customization mechanism. It permits customization on the execution policy. I haven't looked at this for a while, but the important thing is not to get in the way of that. The technique I showed would happen in the existing serial "fall-back" implementation.

[georgemalerbo]

Thanks for your help! I'm really interested in contributing in any way I can. As soon as I started to really look closely at how those types and arguments are structured across the parts of the cblas.h header I needed, I realized a lot of it could be broken down into a giant repeating grid. I recently created a python script as an experiment.

I haven't uploaded anything to any repository yet but this is what I've been researching and working on recently. Let me know what you think.

I wanted to create a short and concise script as clean and readable as possible, and most of the repeated parameters from the cblas.h signatures map nicely to tiny two letter shorthand variables at the top of my file with definition builders for the routine arguments

Int        = lambda n: Param(n, "blas_int_t {n}")
Scalar     = lambda n: Param(n, "{cpp} {n}", arg_fmt_complex="&{n}")
RealScalar = lambda n: Param(n, "{real} {n}")
Const      = lambda n: Param(n, "const {cpp}* {n}")
Mut        = lambda n: Param(n, "{cpp}* {n}")
MutReal    = lambda n: Param(n, "{real}* {n}")

LO, UP, TR, DI, SD   = Order("LAYOUT"), Uplo("UPLO"), Trans("TRANS"), Diag("DIAG"), Side("SIDE")
TRA, TRB             = Trans("TRANSA"), Trans("TRANSB")
M, N, K              = Int("M"), Int("N"), Int("K")
LDA, LDB, LDC        = Int("LDA"), Int("LDB"), Int("LDC")
INCX, INCY           = Int("INCX"), Int("INCY")
Xc, Xm, Yc, Ym       = Const("X"), Mut("X"), Const("Y"), Mut("Y")
Ac, Am, Bc, Bm, Cm   = Const("A"), Mut("A"), Const("B"), Mut("B"), Mut("C")
AL, BE, ALr, BEr     = Scalar("ALPHA"), Scalar("BETA"), RealScalar("ALPHA"), RealScalar("BETA")

With those in place a lot of it collapses into a clean and readable list where the shorthand blocks plug in like this:

Routine("gemm", ALL, (LO, TRA, TRB, M, N, K, AL, Ac, LDA, Bc, LDB, BE, Cm, LDC), level=3)

For the edge cases it was a matter of patching the script whenever it hit an exception to the conventions of the others. To account for one of those cases, the type list includes a real_prefix tracker, and the naming rule handles it directly:

cblas="cblas_{real_prefix}{prefix}nrm2"

Complex dot product functions, those ended up being one of the odd cases because I see they don’t actually return a value directly and instead they pass the result out through an extra pointer parameter at the very end of the function call. The generator handles that by passing along a small layout override rule, out_param, so those routines can be emitted differently from the normal return value cases.

There's a small conditional formatting rule keeps everything aligned at 80 columns and under with a helper to make sure argument lists slice up and line wrap into clean balanced blocks of code.

By designing the script and mapping it out in this way, it comes out to about 280 lines of python code, which outputs about 2,600 lines of C++ declarations and definitions.

[georgemalerbo]

Thanks for putting this up I think I see it clearly now

I'd been thinking of the cblas backend as it's own self contained layer that I figured would just drop into stdBLAS once finished: a codegen layer emitting the cblas_* forwarders for all four scalar types from one catalog, and a dispatch layer on top.

flowchart LR
    G[generator] --> T["cblas_* forwarders<br/>s · d · c · z"]
    T --> D["maybe_can_blas_* / try_blas_*<br/>order · lda · aliasing · tags"]
    D -.->|drop in later| SB["stdBLAS<br/>serial fallback"]

Loading

But I get it now. The integration step is the point, so there's no need for a separate layer. The wrappers go straight into the real files, around the existing bodies.

flowchart LR
    G[generator] --> M["maybe_can_blas_scale<br/>try_blas_scale"]
    H["blas_helpers.hpp<br/>traits"] --> M
    M --> S["scale(inline_exec_t&&)<br/>serial fallback"]
    S --> P["blas1_scale.hpp"]

Loading

Is that what we're after? If so I can start reworking my generator script and experiment with this :)

[mhoemmen]

Hi @georgemalerbo ! You're making great progress! I see that you're understanding the patterns in the BLAS' abbreviated function names more and more.

The idea is that you can start with an existing C BLAS, with its cblas.h header file. You don't have to generate the C BLAS functions. They already exist for you. Just start by installing the OpenBLAS and including its cblas.h header in your library.

We're focusing on this step because OpenBLAS is an open-source library with good performance on a variety of hardware. Different BLAS libraries that provide a C BLAS might have different header file names, but we can deal with that in the build system. Some BLAS libraries only have a Fortran BLAS, but a good way to deal with calling Fortran BLAS functions from C or C++ is to build a C BLAS. At that point, your script for generating C BLAS function declarations would be helpful, but we don't have to worry about that yet.

Does this make sense? A good way to start would be with my suggested DSCAL wrapper above.

[amklinv-nnl]

Hi @mhoemmen! I tried implementing your suggested scal wrapper in PR #331. Could you review and make suggestions at your convenience?

[georgemalerbo]

just wanted to share something I discovered

I built #331 locally and tried it against a few different cblas libraries.

if (x.is_strided() && n <= std::numeric_limits<CBLAS_INDEX>::max()) {

works on OpenBLAS and MKL, but it doesn't compile against BLIS or AMD's AOCL they don't define CBLAS_INDEX:

error: 'CBLAS_INDEX' was not declared in this scope; did you mean 'CBLAS_SIDE'?
   18 |   n <= std::numeric_limits<CBLAS_INDEX>::max()
      |                            ^~~~~~~~~~~
      |                            CBLAS_SIDE
error: template argument 1 is invalid

AMD's header doesn't have it but MKL's does:

$ grep CBLAS_INDEX /usr/local/amd-blis/include/LP64/cblas.h
nothing

$ grep "#define CBLAS_INDEX" /opt/intel/oneapi/mkl/2026.0/include/mkl_cblas.h
#define CBLAS_INDEX size_t /* this may vary between platforms */

Just figured I'd share hope it helps

[amklinv-nnl]

@mhoemmen Is the "good enough" solution to temporarily enforce an (optional) dependency on OpenBLAS or MKL as a BLAS provider? There is a mechanism to do this in CMake.

[georgemalerbo]

I was looking at how kokkos kernels handles this, They define their own integer type, KK_INT:

#if defined(KOKKOSKERNELS_ENABLE_TPL_MKL) && defined(MKL_PROVIDES_BLAS_LAPACK)
using KK_INT = MKL_INT;
#else
using KK_INT = int;
#endif

We could mirror that and we call ours linalg_int or something and define it in blas_helpers.hpp:

#ifndef LINALG_INT
#define LINALG_INT int
#endif
using linalg_int = LINALG_INT;

and then that line becomes:

if (x.is_strided() && n <= std::numeric_limits<linalg_int>::max()) {

I tried this locally against OpenBLAS, MKL, BLIS, and AMD and it compiles on all four.