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]

# the wrappers that do nothing but forward the call straight into cblas_*
# this is not the full cblas backend
# the inside of each wrapper i just copy off the cblas.h signature - args, order, types are all right there.


# the struct convention & shape itself i.e. BlasDot<float>::dot I used around it is
# not dictated by cblas.h so I need to check if it's what stdblas actually wants?

# Is this a good convention for how stdBLAS should find these wrappers?
# It's C++17
# idea: each routine as its own little class template, Blas<Op><T>,
# with one static method inside that just forwards to the matching cblas_ call which this script generates

# a clean separation of concerns :)

# The cblas_* call = the what which vendor function, with which args. 
# The Blas<Op><T> struct shell = the how it's exposed the handle the dispatch layer finds and calls. 

# looking at the types and args across cblas.h header realized most of it is
# a big repeating grid. So this script is that grid with each routine listed once in shorthand

# right now it only writes two things: the structs that declare each routine,
# and the little functions that just pass the call on to cblas.
# all the setup that has to go around them - the #includes, the integer type
# used for sizes, the struct templates these plug into, and the namespace and
# include guard - i haven't worked out yet, so i left it out for now.
# to add a routine you just add one more line to the list below.

#
# A few routines are weird, so there is a small override
# for each as I hit it those are the odd cases flagged in the comments.

from dataclasses import dataclass
from math import ceil
from pathlib import Path


# BLAS scalar types

@dataclass(frozen=True) 
class Elem:
    # One of the four BLAS scalar types. Each carries its single-letter prefix
    # (s, d, c, z) plus a real_prefix tracking the prefix of its underlying real
    # type. For the real types that's just the same letter.
    # For complex it's the s or d of the real part - a few odd routines stitch
    # that onto the regular prefix to name their function (sc/dz for nrm2/asum,
    # cs/zd for rot), which is the only reason real_prefix exists.
    cpp: str; prefix: str; is_complex: bool; real: str; real_prefix: str

S = Elem("float",                "s", False, "float",  "s")
D = Elem("double",               "d", False, "double", "d")
C = Elem("std::complex<float>",  "c", True,  "float",  "s")
Z = Elem("std::complex<double>", "z", True,  "double", "d")
REAL, COMPLEX = (S, D), (C, Z)   # the real ones and the complex ones
ALL = REAL + COMPLEX  # all four


# Parameters

@dataclass(frozen=True)
class Param:
    # One argument of a CBLAS routine. decl(t) builds how it looks in the
    # signature, arg(t) how it looks at the call site.
    # The format string fills in the arg name, the scalar type, or its real type.
    # arg_fmt_complex is the one wrinkle - ALPHA/BETA go by value for the real
    # types but need an & for the complex ones, so when its set it takes over
    # for complex only.
    name: str
    decl_fmt: str
    arg_fmt: str = "{n}"
    arg_fmt_complex: str | None = None

    def decl(self, t): return self.decl_fmt.format(n=self.name, cpp=t.cpp, real=t.real)
    def arg(self, t):  return (self.arg_fmt_complex if t.is_complex and self.arg_fmt_complex else self.arg_fmt).format(n=self.name)


# Definition builders for the args. Most params in the cblas.h signatures fall
# into a few shapes and each builder stamps one out:
# Int is a plain integer index; Scalar is a value like ALPHA/BETA (by value for
# real, by pointer for complex); RealScalar stays real even when the routine is
# complex (a few do); Const/Mut/MutReal are the pointer args (read-only,
# writable, writable-real); and the rest are the little CBLAS option enums -
#layout, transpose, uplo, side, diag.
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}")
# _enum is another one of those internal helpers (the _ thing again)
def _enum(k): return lambda n: Param(n, f"CBLAS_{k} {{n}}")
Order, Trans, Uplo, Side, Diag = map(_enum, ("ORDER", "TRANSPOSE", "UPLO", "SIDE", "DIAG"))


# Routines

# the _ in front means its meant to be internal/private i think? like dont
# import this from outside, its only a helper. 
def _wrap(prefix, items, indent=2, limit=80):
    # Little formatting helper that keeps the generated calls at 80 cols or
    # under. Fits on one line if it can.
    # If not, it slices the arg list into even blocks and wraps, trying the fewest
    # lines first and adding more only until it fits. Even slices keep the blocks
    # balanced instead of leaving a lonely orphan like `      LDC);` hanging off
    # the end the way greedy wrapping does.
    head = " " * indent + prefix
    cont = " " * len(head)
    single = head + ", ".join(items) + ");"
    n = len(items)
    if len(single) <= limit or n < 2:
        return single  # easy, it all fits on one line
    for k in range(2, n + 1):
        per = ceil(n / k)
        groups = [items[i:i + per] for i in range(0, n, per)]
        lines  = [head + ", ".join(groups[0]) + ","]
        lines += [cont + ", ".join(g) + "," for g in groups[1:-1]]
        lines.append(cont + ", ".join(groups[-1]) + ");")
        if k == n or all(len(l) <= limit for l in lines):
            return "\n".join(lines)


@dataclass(frozen=True)
class Routine:
    # One row of the grid a single CBLAS routine for one set of scalar types.
    # declare(t) writes the struct declaration, specialize(t) the definition that
    # forwards into cblas. The name strings fill in the scalar type, its prefix,
    # real prefix, real type, and the routine name.
    #
    # The last few fields handle the routines that don't behave:
    # out_param is for the ones that hand their answer back through a pointer arg
    # instead of returning it; cast wraps the result when cblas gives back a wider
    # type than I want; note is just a comment block printed above the def when a
    # case needs it.
    name:      str
    types:     tuple
    params:    tuple
    level:     int  = 1
    returns:   str  = "void"
    cast:      bool = False
    cblas:     str  = "cblas_{prefix}{name}"
    out_param: str | None = None
    note:      str | None = None

    @property
    def struct(self): return f"Blas{self.name.capitalize()}"

    # _fmt = helper just for this class (one _ = internal). the double __ ones
    # like the __init__ that @dataclass writes for you are pythons own special
    # built-in methods, still getting the hang of when its one underscore vs two
    def _fmt(self, tmpl, t):
        return tmpl.format(T=t.cpp, prefix=t.prefix, real_prefix=t.real_prefix,
                           name=self.name, real=t.real)

    def declare(self, t):
        ret = self._fmt(self.returns, t)
        sig = ",\n      ".join(p.decl(t) for p in self.params)
        return (f"template <> struct {self.struct}<{t.cpp}> {{\n"
                f"  static {ret} {self.name}(\n      {sig});\n}};\n")

    def specialize(self, t):
        ret, fn = self._fmt(self.returns, t), self._fmt(self.cblas, t)
        args = [p.arg(t) for p in self.params]
        if self.out_param is not None:
            args.append(f"&{self.out_param}")
            body = (f"  {ret} {self.out_param}{{}};\n"
                    f"{_wrap(fn + '(', args)}\n  return {self.out_param};")
        elif self.cast:
            body = f"  return static_cast<{ret}>({fn}({', '.join(args)}));"
        elif ret != "void":
            body = f"  return {fn}({', '.join(args)});"
        else:
            body = _wrap(fn + "(", args)
        sig = ",\n    ".join(p.decl(t) for p in self.params)
        return (f"inline {ret} {self.struct}<{t.cpp}>::{self.name}(\n"
                f"    {sig}) {{\n{body}\n}}\n")


# Routine catalog
# The tiny two-letter shorthands the routine list is built from. Each is a param
# from the cblas.h signatures, named so a whole signature reads as a short block.
# LO/UP/TR/DI/SD are layout, uplo, transpose, diag, side; TRA/TRB the two
# transpose flags; M/N/K the sizes; LD and INC names are strides;
# the c/m suffix on X, Y, A, B, C, P means const vs mutable;
# AL/BE are ALPHA/BETA, ALr/BEr their always-real versions;
# RC/RS the cosine and sine for the rotations.

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, Pc   = Const("X"), Mut("X"), Const("Y"), Mut("Y"), Const("P")
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")
RC, RS               = RealScalar("C"), RealScalar("S")

# a tuple of any number of Routines
ROUTINES: tuple[Routine, ...] = (
    # Level 1
    # dot: real types just return the answer.
    # complex is an odd case it doesn't return, it passes the result back
    # out through an extra pointer at the end (the *_sub form), so those use
    # out_param.
    # dotc: conjugated dot product, complex-only - for reals it's the same as dot.
    Routine("dot",  REAL,    (N, Xc, INCX, Yc, INCY), returns="{T}"),
    Routine("dot",  COMPLEX, (N, Xc, INCX, Yc, INCY), returns="{T}", out_param="result", cblas="cblas_{prefix}dotu_sub"),
    Routine("dotc", COMPLEX, (N, Xc, INCX, Yc, INCY), returns="{T}", out_param="result", cblas="cblas_{prefix}dotc_sub"),
    # nrm2/asum: for complex these give back a real number, and cblas names them
    # with the stitched sc/dz prefix from real_prefix.
    Routine("nrm2", REAL,    (N, Xc, INCX), returns="{T}"),
    Routine("nrm2", COMPLEX, (N, Xc, INCX), returns="{real}", cblas="cblas_{real_prefix}{prefix}nrm2"),
    Routine("asum", REAL,    (N, Xc, INCX), returns="{T}"),
    Routine("asum", COMPLEX, (N, Xc, INCX), returns="{real}", cblas="cblas_{real_prefix}{prefix}asum"),
    Routine("scal", ALL, (N, AL, Xm, INCX)),
    Routine("axpy", ALL, (N, AL, Xc, INCX, Ym, INCY)),
    # iamax: position of the largest element. cblas folds all four types into
    # the one i?amax name, so this just slots in the prefix.
    Routine("iamax", ALL, (N, Xc, INCX),
            returns="blas_int_t", cast=True, cblas="cblas_i{prefix}amax",
            note=("iamax returns a CBLAS_INDEX, an unsigned size_t. Squeezing that into\n"
                  "blas_int_t is fine here - it's just a position inside the vector, so\n"
                  "it always fits, int or int64_t build alike.")),
    Routine("copy", ALL, (N, Xc, INCX, Ym, INCY)),
    Routine("swap", ALL, (N, Xm, INCX, Ym, INCY)),
    # rot: applies a plane rotation. its cosine and sine stay real even for the
    # complex types, which land on the stitched csrot/zdrot names (crot/zrot
    # aren't in cblas, they live in LAPACK).
    Routine("rot", REAL,    (N, Xm, INCX, Ym, INCY, RC, RS)),
    Routine("rot", COMPLEX, (N, Xm, INCX, Ym, INCY, RC, RS), cblas="cblas_{prefix}{real_prefix}rot"),
    Routine("rotg",  ALL,  (Mut("A"), Mut("B"), MutReal("C"), Mut("S"))),
    Routine("rotm",  REAL, (N, Xm, INCX, Ym, INCY, Pc)),
    Routine("rotmg", REAL, (Mut("D1"), Mut("D2"), Mut("B1"), Scalar("B2"), Mut("P"))),

    # Level 2
    Routine("gemv", ALL, (LO, TR, M, N, AL, Ac, LDA, Xc, INCX, BE, Ym, INCY), level=2),
    # the * here unpacks the generator straight into the tuple. confused me for a while
    *(Routine(n, ALL, (LO, UP, TR, DI, N, Ac, LDA, Xm, INCX), level=2) for n in ("trsv", "trmv")),
    # ger/geru/gerc all share the exact same arg shape - ger for the real types,
    # geru and gerc the two complex versions.
    *(Routine(n, ts, (LO, M, N, AL, Xc, INCX, Yc, INCY, Am, LDA), level=2)
      for n, ts in (("ger", REAL), ("geru", COMPLEX), ("gerc", COMPLEX))),
    # symv/hemv - another real+complex pair on one shape: symv real, hemv complex.
    *(Routine(n, ts, (LO, UP, N, AL, Ac, LDA, Xc, INCX, BE, Ym, INCY), level=2)
      for n, ts in (("symv", REAL), ("hemv", COMPLEX))),
    # syr/her share a shape too, except her's ALPHA stays real.
    Routine("syr", REAL,    (LO, UP, N, AL,  Xc, INCX, Am, LDA), level=2),
    Routine("her", COMPLEX, (LO, UP, N, ALr, Xc, INCX, Am, LDA), level=2),
    # syr2/her2, the two-vector pair; her2 passes ALPHA by pointer.
    *(Routine(n, ts, (LO, UP, N, AL, Xc, INCX, Yc, INCY, Am, LDA), level=2)
      for n, ts in (("syr2", REAL), ("her2", COMPLEX))),

    # Level 3
    Routine("gemm", ALL,     (LO, TRA, TRB, M, N, K, AL, Ac, LDA, Bc, LDB, BE, Cm, LDC), level=3),
    Routine("trsm", ALL,     (LO, SD, UP, TR, DI, M, N, AL, Ac, LDA, Bm, LDB), level=3),
    Routine("symm", ALL,     (LO, SD, UP, M, N, AL, Ac, LDA, Bc, LDB, BE, Cm, LDC), level=3),
    Routine("trmm", ALL,     (LO, SD, UP, TR, DI, M, N, AL, Ac, LDA, Bm, LDB), level=3),
    Routine("syrk", ALL,     (LO, UP, TR, N, K, AL, Ac, LDA, BE, Cm, LDC), level=3),
    # hemm is complex-only.
    Routine("hemm", COMPLEX, (LO, SD, UP, M, N, AL, Ac, LDA, Bc, LDB, BE, Cm, LDC), level=3),
    # herk/her2k are complex-only too, ALPHA/BETA stay real (her2k still passes
    # ALPHA by pointer like the rest). pretty sure thats right
    Routine("herk",  COMPLEX, (LO, UP, TR, N, K, ALr, Ac, LDA,         BEr, Cm, LDC), level=3),
    Routine("her2k", COMPLEX, (LO, UP, TR, N, K, AL,  Ac, LDA, Bc, LDB, BEr, Cm, LDC), level=3),
    Routine("syr2k", ALL,     (LO, UP, TR, N, K, AL,  Ac, LDA, Bc, LDB, BE,  Cm, LDC), level=3),
)
# Emit

# Only emits the declarations and definitions
# templates, namespace and guard still need adding by hand. Haven't sorted that
# out yet.


def generate() -> str:
    # Build just the declarations and definitions.
    # The plumbing around them (includes, namespace, guard, index type, struct
    # templates) is left out.
    # all the declarations first, then all the definitions.
    parts = []
    for method in ("declare", "specialize"):
        for r in ROUTINES:
            if method == "specialize" and r.note:
                parts.append("\n".join(f"// {l}" for l in r.note.splitlines()))
            parts += [getattr(r, method)(t) for t in r.types]
    return "\n".join(parts)


def main() -> None:
    # always just write the output to linalg_cblas.hpp, no args, no printing.
    Path("linalg_cblas.hpp").write_text(generate())
    # print(generate())  # was using this to eyeball the output while testing


main()