blas.hpp

Go to the documentation of this file.

 
#ifndef NUMERICAL_BLAS_HPP
#define NUMERICAL_BLAS_HPP
 
#include <algorithm>
#include <array>
#include <atomic>
#include <cmath>
#include <complex>
#include <cstddef>
#include <cstdint>
#include <cstring>
#include <limits>
#include <optional>
#include <utility>
#include <vector>
 
#include "decomposition.hpp"
#include "eigen.hpp"
#include "gemm.hpp"
#include "matrix.hpp"
#include "parallel.hpp"
 
namespace numerical::blas {
 
  inline constexpr double lapack_precision {0x1p-53};
 
  inline constexpr std::size_t syrk_block {128};
 
  inline constexpr std::size_t syrk_kept_limit {(std::size_t {32} << 20U) / sizeof(double)};
 
  namespace detail {
 
    [[nodiscard]] inline buffer<double>& syrk_operand() noexcept
    {
      thread_local buffer<double> storage;
      return storage;
    }
 
    [[nodiscard]] constexpr bool is_operation(char op) noexcept
    {
      return op == 'N' || op == 'n' || op == 'T' || op == 't' || op == 'C' || op == 'c';
    }
 
    [[nodiscard]] constexpr bool transposes(char op) noexcept
    {
      return op != 'N' && op != 'n';
    }
 
    [[nodiscard]] constexpr bool conjugates(char op) noexcept
    {
      return op == 'C' || op == 'c';
    }
 
    template <typename T>
    [[nodiscard]] bool finite(const T& x)
    {
      return std::isfinite(x);
    }
 
    template <typename U>
    [[nodiscard]] bool finite(const std::complex<U>& x)
    {
      return std::isfinite(x.real()) && std::isfinite(x.imag());
    }
 
    inline constexpr std::size_t transpose_tile {32};
 
    inline constexpr std::size_t transpose_weight {4};
 
    template <typename T>
    void row_major_into(buffer<T>&     out,
                        const T*       a,
                        std::ptrdiff_t lda,
                        std::ptrdiff_t rows,
                        std::ptrdiff_t cols,
                        T              scale,
                        bool           conjugate = false)
    {
      const auto r  {static_cast<std::size_t>(rows)};
      const auto c  {static_cast<std::size_t>(cols)};
      const auto ld {static_cast<std::size_t>(lda)};
      out.resize(r * c);
      T* const to   {out.data()};
      numerical::detail::parallel_for((r + transpose_tile - 1) / transpose_tile, transpose_weight * transpose_tile * c, [a, ld, r, c, to, scale, conjugate](std::size_t band) {
                                        const std::size_t i0 {band * transpose_tile};
                                        for (std::size_t j0 {0}; j0 < c; j0 += transpose_tile) {
                                          for (std::size_t j {j0}; j < std::min(j0 + transpose_tile, c); ++j) {
                                            for (std::size_t i {i0}; i < std::min(i0 + transpose_tile, r); ++i) {
                                              const T x {a[i + (j * ld)]};
                                              to[(i * c) + j] = scale * (conjugate ? numerical::detail::conj_of(x) : x);
                                            }
                                          }
                                        }
                                      });
    }
 
    template <typename T>
    [[nodiscard]] buffer<T> row_major(const T*       a,
                                      std::ptrdiff_t lda,
                                      std::ptrdiff_t rows,
                                      std::ptrdiff_t cols,
                                      T              scale,
                                      bool           conjugate = false)
    {
      buffer<T> out;
      row_major_into(out, a, lda, rows, cols, scale, conjugate);
      return out;
    }
 
    inline constexpr std::size_t kept_bytes {std::size_t {32} << 20U};
 
    template <typename T>
    [[nodiscard]] buffer<T>& kept_storage(std::size_t slot) noexcept
    {
      thread_local std::array<buffer<T>, 2> kept;
      return kept[slot];
    }
 
    template <typename T>
    [[nodiscard]] buffer<T> take_kept(std::size_t slot)
    {
      return std::exchange(kept_storage<T>(slot), buffer<T> {});
    }
 
    template <typename T>
    void give_back(std::size_t slot,
                   buffer<T>&& storage)
    {
      if (storage.capacity() <= kept_bytes / sizeof(T)) {
        kept_storage<T>(slot) = std::move(storage);
      }
    }
 
    template <typename T>
    [[nodiscard]] bool is_identity(const T*    b,
                                   std::size_t ldb,
                                   std::size_t n)
    {
      buffer<T> unit((2 * n) - 1, T {});
      unit[n - 1] = T {1};
      const T* const middle         {unit.data() + (n - 1)};
      if (std::memcmp(b, middle, n * sizeof(T)) != 0) {
        return false;
      }
      std::atomic<bool> holds {true};
      numerical::detail::parallel_for(n - 1, transpose_weight * n, [b, ldb, n, middle, &holds](std::size_t j) {
                                        if (std::memcmp(b + ((j + 1) * ldb), middle - (j + 1), n * sizeof(T)) != 0) {
                                          holds.store(false, std::memory_order_relaxed);
                                        }
                                      });
      return holds.load(std::memory_order_relaxed);
    }
 
    template <typename T>
    void column_major_into(const T*    from,
                           std::size_t rows,
                           std::size_t cols,
                           T*          a,
                           std::size_t lda)
    {
      numerical::detail::parallel_for((cols + transpose_tile - 1) / transpose_tile, transpose_weight * transpose_tile * rows, [from, rows, cols, a, lda](std::size_t band) {
                                        const std::size_t j0 {band * transpose_tile};
                                        for (std::size_t i0 {0}; i0 < rows; i0 += transpose_tile) {
                                          for (std::size_t i {i0}; i < std::min(i0 + transpose_tile, rows); ++i) {
                                            for (std::size_t j {j0}; j < std::min(j0 + transpose_tile, cols); ++j) {
                                              a[i + (j * lda)] = from[(i * cols) + j];
                                            }
                                          }
                                        }
                                      });
    }
 
    template <typename T>
    void column_major_permuted_into(const T*           from,
                                    std::size_t        n,
                                    const std::size_t* to,
                                    T*                 a,
                                    std::size_t        lda)
    {
      numerical::detail::parallel_for((n + transpose_tile - 1) / transpose_tile, transpose_weight * transpose_tile * n, [from, n, to, a, lda](std::size_t band) {
                                        const std::size_t j0 {band * transpose_tile};
                                        for (std::size_t i0 {0}; i0 < n; i0 += transpose_tile) {
                                          for (std::size_t i {i0}; i < std::min(i0 + transpose_tile, n); ++i) {
                                            for (std::size_t j {j0}; j < std::min(j0 + transpose_tile, n); ++j) {
                                              a[i + (to[j] * lda)] = from[(i * n) + j];
                                            }
                                          }
                                        }
                                      });
    }
 
    template <typename T>
    void product_by_op_a(char        transa,
                         std::size_t m,
                         std::size_t n,
                         std::size_t k,
                         const T*    l,
                         std::size_t ldl,
                         const T*    a,
                         std::size_t lda,
                         T*          c,
                         std::size_t ldc)
    {
      const bool      copied {transposes(transa)};
      const buffer<T> right  {copied ? row_major(a, static_cast<std::ptrdiff_t>(lda), static_cast<std::ptrdiff_t>(k), static_cast<std::ptrdiff_t>(n), T {1}, conjugates(transa))
                                     : buffer<T> {}};
      numerical::detail::gemm(m, n, k, l, ldl, copied ? right.data() : a, copied ? n : lda, c, ldc);
    }
 
    inline void product_by_op_a(char          transa,
                                std::size_t   m,
                                std::size_t   n,
                                std::size_t   k,
                                const double* l,
                                std::size_t   ldl,
                                const double* a,
                                std::size_t   lda,
                                double*       c,
                                std::size_t   ldc)
    {
      if (transposes(transa)) {
        numerical::detail::gemm_bt(m, n, k, l, ldl, a, lda, c, ldc);
        return;
      }
      numerical::detail::gemm(m, n, k, l, ldl, a, lda, c, ldc);
    }
 
    template <typename T>
    void lower_triangle_rows(T*          a,
                             std::size_t lda,
                             std::size_t n,
                             T*          rows,
                             bool        in)
    {
      numerical::detail::parallel_for((n + transpose_tile - 1) / transpose_tile, transpose_weight * transpose_tile * n / 2, [a, lda, n, rows, in](std::size_t band) {
                                        const std::size_t i0 {band * transpose_tile};
                                        for (std::size_t j0 {0}; j0 <= i0; j0 += transpose_tile) {
                                          for (std::size_t j {j0}; j < std::min(j0 + transpose_tile, n); ++j) {
                                            for (std::size_t i {std::max(i0, j)}; i < std::min(i0 + transpose_tile, n); ++i) {
                                              if (in) {
                                                rows[(i * n) + j] = a[i + (j * lda)];
                                              }
                                              else {
                                                a[i + (j * lda)] = rows[(i * n) + j];
                                              }
                                            }
                                          }
                                        }
                                      });
    }
 
    template <typename T>
    [[nodiscard]] buffer<T> columns_as_rows(const T*       a,
                                            std::ptrdiff_t lda,
                                            std::ptrdiff_t rows,
                                            std::ptrdiff_t cols,
                                            T              scale)
    {
      const auto r {static_cast<std::size_t>(rows)};
      const auto c {static_cast<std::size_t>(cols)};
      buffer<T>  out(r * c);
      for (std::size_t j {0}; j < c; ++j) {
        for (std::size_t i {0}; i < r; ++i) {
          out[(j * r) + i] = scale * a[i + (j * static_cast<std::size_t>(lda))];
        }
      }
      return out;
    }
 
    template <typename T>
    void product_of_ops(char        transa,
                        char        transb,
                        T           alpha,
                        std::size_t m,
                        std::size_t n,
                        std::size_t k,
                        const T*    a,
                        std::size_t lda,
                        const T*    b,
                        std::size_t ldb,
                        T*          c,
                        std::size_t ldc)
    {
      const bool copy_b {transposes(transb) || alpha != T {1}};
      buffer<T>  left   {};
      if (transposes(transb)) {
        left = row_major(b, static_cast<std::ptrdiff_t>(ldb), static_cast<std::ptrdiff_t>(m), static_cast<std::ptrdiff_t>(k), alpha,
                         conjugates(transb));
      }
      else if (copy_b) {
        left = columns_as_rows(b, static_cast<std::ptrdiff_t>(ldb), static_cast<std::ptrdiff_t>(k), static_cast<std::ptrdiff_t>(m), alpha);
      }
      product_by_op_a(transa, m, n, k, copy_b ? left.data() : b, copy_b ? k : ldb, a, lda, c, ldc);
    }
 
    inline void product_of_ops(char          transa,
                               char          transb,
                               double        alpha,
                               std::size_t   m,
                               std::size_t   n,
                               std::size_t   k,
                               const double* a,
                               std::size_t   lda,
                               const double* b,
                               std::size_t   ldb,
                               double*       c,
                               std::size_t   ldc)
    {
      if (!transposes(transb) || alpha != 1.0) {
        product_of_ops<double>(transa, transb, alpha, m, n, k, a, lda, b, ldb, c, ldc);
        return;
      }
      if (transposes(transa)) {
        numerical::detail::gemm_atbt(m, n, k, b, ldb, a, lda, c, ldc);
        return;
      }
      numerical::detail::gemm_at(m, n, k, b, ldb, a, lda, c, ldc);
    }
  } // namespace detail
 
  namespace detail {
 
    template <typename Index>
    void interchanges(const std::vector<std::size_t>& pivots,
                      Index*                          ipiv)
    {
      const std::size_t        size {pivots.size()};
      std::vector<std::size_t> at(size);    // the original row at each position
      std::vector<std::size_t> where(size); // the position of each original row
      for (std::size_t i {0}; i < size; ++i) {
        at[i]    = i;
        where[i] = i;
      }
      for (std::size_t i {0}; i < size; ++i) {
        const std::size_t from {where[pivots[i]]};
        ipiv[i]          = static_cast<Index>(from + 1);
        where[at[i]]     = from;
        where[pivots[i]] = i;
        std::swap(at[i], at[from]);
      }
    }
 
    template <typename T, typename Index>
    int rectangular_lu(std::size_t m,
                       std::size_t n,
                       T*          a,
                       std::size_t lda,
                       Index*      ipiv)
    {
      int info {0};
      for (std::size_t j {0}; j < std::min(m, n); ++j) {
        T* const    column {a + (j * lda)};
        std::size_t pivot  {j};
        for (std::size_t i {j + 1}; i < m; ++i) {
          if (std::abs(column[i]) > std::abs(column[pivot])) {
            pivot = i;
          }
        }
        ipiv[j] = static_cast<Index>(pivot + 1);
        if (column[pivot] == T {}) {
          info = info == 0 ? static_cast<int>(j + 1) : info;
          continue;
        }
        for (std::size_t c {0}; c < n; ++c) { // the whole rows interchanged, as dgetf2's dswap
          std::swap(a[j + (c * lda)], a[pivot + (c * lda)]);
        }
        for (std::size_t i {j + 1}; i < m; ++i) {
          column[i] = column[i] / column[j];
        }
        for (std::size_t c {j + 1}; c < n; ++c) {
          T* const target {a + (c * lda)};
          const T  factor {target[j]};
          for (std::size_t i {j + 1}; i < m; ++i) {
            target[i] -= column[i] * factor;
          }
        }
      }
      return info;
    }
  } // namespace detail
 
  template <typename T>
  int gemm(char           transa,
           char           transb,
           std::ptrdiff_t m,
           std::ptrdiff_t n,
           std::ptrdiff_t k,
           T              alpha,
           const T*       a,
           std::ptrdiff_t lda,
           const T*       b,
           std::ptrdiff_t ldb,
           T              beta,
           T*             c,
           std::ptrdiff_t ldc)
  {
    const std::ptrdiff_t a_rows  {detail::transposes(transa) ? k : m};
    const std::ptrdiff_t b_rows  {detail::transposes(transb) ? n : k};
    int                  illegal {0};
    if (!detail::is_operation(transa)) {
      illegal = 1;
    }
    else if (!detail::is_operation(transb)) {
      illegal = 2;
    }
    else if (m < 0) {
      illegal = 3;
    }
    else if (n < 0) {
      illegal = 4;
    }
    else if (k < 0) {
      illegal = 5;
    }
    else if (lda < std::max<std::ptrdiff_t>(1, a_rows)) {
      illegal = 8;
    }
    else if (ldb < std::max<std::ptrdiff_t>(1, b_rows)) {
      illegal = 10;
    }
    else if (ldc < std::max<std::ptrdiff_t>(1, m)) {
      illegal = 13;
    }
    if (illegal != 0 || m == 0 || n == 0) {
      return illegal;
    }
    const auto rows   {static_cast<std::size_t>(m)};
    const auto cols   {static_cast<std::size_t>(n)};
    const auto stride {static_cast<std::size_t>(ldc)};
    for (std::size_t j {0}; j < cols; ++j) { // C ← β C, column by column
      T* const column {c + (j * stride)};
      if (beta == T {}) {
        std::fill_n(column, rows, T {});
      }
      else if (beta != T {1}) {
        std::transform(column, column + rows, column, [beta](T x) { return beta * x; });
      }
    }
    if (alpha == T {} || k == 0) {
      return 0;
    }
    detail::product_of_ops(transa, transb, alpha, cols, rows, static_cast<std::size_t>(k), a, static_cast<std::size_t>(lda), b,
                           static_cast<std::size_t>(ldb), c, stride);
    return 0;
  }
 
  namespace detail {
 
    [[nodiscard]] constexpr std::size_t strided(std::size_t    i,
                                                std::size_t    length,
                                                std::ptrdiff_t inc) noexcept
    {
      return inc > 0 ? i * static_cast<std::size_t>(inc) : (length - 1 - i) * static_cast<std::size_t>(-inc);
    }
 
    template <typename T>
    [[nodiscard]] T real_dot(std::ptrdiff_t n,
                             const T*       x,
                             std::ptrdiff_t incx,
                             const T*       y,
                             std::ptrdiff_t incy)
    {
      if (n <= 0) {
        return T {};
      }
      const auto count {static_cast<std::size_t>(n)};
      if (incx == 1 && incy == 1) {
        return numerical::detail::chunked_dot(x, y, count);
      }
      buffer<T> xs(count, T {});
      buffer<T> ys(count, T {});
      for (std::size_t i {0}; i < count; ++i) {
        xs[i] = x[strided(i, count, incx)];
        ys[i] = y[strided(i, count, incy)];
      }
      return numerical::detail::chunked_dot(xs.data(), ys.data(), count);
    }
  } // namespace detail
 
  [[nodiscard]] inline double dot(std::ptrdiff_t n,
                                  const double*  x,
                                  std::ptrdiff_t incx,
                                  const double*  y,
                                  std::ptrdiff_t incy)
  {
    return detail::real_dot(n, x, incx, y, incy);
  }
 
  [[nodiscard]] inline float dot(std::ptrdiff_t n,
                                 const float*   x,
                                 std::ptrdiff_t incx,
                                 const float*   y,
                                 std::ptrdiff_t incy)
  {
    return detail::real_dot(n, x, incx, y, incy);
  }
 
  template <typename T>
  [[nodiscard]] std::complex<T> complex_dot(bool                   conjugate,
                                            std::ptrdiff_t         n,
                                            const std::complex<T>* x,
                                            std::ptrdiff_t         incx,
                                            const std::complex<T>* y,
                                            std::ptrdiff_t         incy) noexcept
  {
    std::complex<T> sum {};
    if (n <= 0) {
      return sum;
    }
    const auto count {static_cast<std::size_t>(n)};
    if (incx == 1 && incy == 1) {
      return conjugate ? numerical::detail::conjugate_dot_product(x, y, count) : numerical::detail::complex_dot_product(x, y, count);
    }
    for (std::size_t i {0}; i < count; ++i) {
      const std::complex<T> xi {x[detail::strided(i, count, incx)]};
      sum += (conjugate ? std::conj(xi) : xi) * y[detail::strided(i, count, incy)];
    }
    return sum;
  }
 
  inline constexpr std::size_t gemv_rows {64};
 
  inline constexpr std::size_t gemv_weight {4};
 
  inline int gemv(char           trans,
                  std::ptrdiff_t m,
                  std::ptrdiff_t n,
                  double         alpha,
                  const double*  a,
                  std::ptrdiff_t lda,
                  const double*  x,
                  std::ptrdiff_t incx,
                  double         beta,
                  double*        y,
                  std::ptrdiff_t incy)
  {
    if (!detail::is_operation(trans)) {
      return 1;
    }
    if (m < 0) {
      return 2;
    }
    if (n < 0) {
      return 3;
    }
    if (lda < std::max<std::ptrdiff_t>(1, m)) {
      return 6;
    }
    if (incx == 0) {
      return 8;
    }
    if (incy == 0) {
      return 11;
    }
    if (m == 0 || n == 0 || (alpha == 0.0 && beta == 1.0)) {
      return 0;
    }
    const bool        transposed {detail::transposes(trans)};
    const auto        rows       {static_cast<std::size_t>(m)};
    const auto        cols       {static_cast<std::size_t>(n)};
    const auto        ld         {static_cast<std::size_t>(lda)};
    const std::size_t x_length   {transposed ? rows : cols};
    const std::size_t y_length   {transposed ? cols : rows};
    // y itself where its elements are contiguous, else a copy written back; x itself for Aᵀ x where it is contiguous
    // and α 1, else α x gathered (A x takes −α x, its factors): the same values either way, so the same bits.
    buffer<double>    gathered_y(incy == 1 ? 0 : y_length);
    double* const     out {incy == 1 ? y : gathered_y.data()};
    for (std::size_t i {0}; i < y_length; ++i) {
      out[i] = beta == 0.0 ? 0.0 : beta * y[detail::strided(i, y_length, incy)];
    }
    if (alpha != 0.0) {
      const bool     in_place {transposed && incx == 1 && alpha == 1.0};
      buffer<double> scaled(in_place ? 0 : x_length); // α x, gathered
      for (std::size_t i {0}; i < scaled.size(); ++i) {
        scaled[i] = alpha * x[detail::strided(i, x_length, incx)];
      }
      if (transposed) {
        const double* const factors {in_place ? x : scaled.data()};
        buffer<double>      dots(cols, 0.0);
        numerical::detail::parallel_for((cols + gemv_rows - 1) / gemv_rows, gemv_weight * gemv_rows * rows, [&](std::size_t group) {
                                          const std::size_t first {group * gemv_rows};
                                          numerical::detail::row_dots(a + (first * ld), ld, std::min(gemv_rows, cols - first),
                                                                      factors, rows, dots.data() + first);
                                        });
        for (std::size_t j {0}; j < cols; ++j) {
          out[j] += dots[j];
        }
      }
      else { // A x: each task's rows of y held through every column (numerical::detail::subtract_products, the
             // factors −α x_j), each element fma(α x_j, a_ij, y_i) for j ascending, as column by column
        for (double& factor : scaled) {
          factor = -factor;
        }
        numerical::detail::parallel_for((rows + gemv_rows - 1) / gemv_rows, gemv_weight * gemv_rows * cols, [&](std::size_t group) {
                                          const std::size_t first {group * gemv_rows};
                                          numerical::detail::subtract_products(out + first, scaled.data(), a + first, ld, cols,
                                                                               std::min(gemv_rows, rows - first));
                                        });
      }
    }
    for (std::size_t i {0}; incy != 1 && i < y_length; ++i) {
      y[detail::strided(i, y_length, incy)] = out[i];
    }
    return 0;
  }
 
  namespace detail {
 
    struct syrk_operands
    {
      const double* p;       
      std::size_t   ldp;     
      const double* q;       
      std::size_t   ldq;     
      std::size_t   depth;   
      double*       c;       
      std::size_t   ldc;     
      bool          lower;   
      bool          q_is_pt; 
      bool          p_is_qt; 
    };
 
    inline void syrk_product(const syrk_operands& operands,
                             std::size_t          c0,
                             std::size_t          cols,
                             std::size_t          r0,
                             std::size_t          rows,
                             double*              out,
                             std::size_t          ldo)
    {
      if (operands.q_is_pt) {
        numerical::detail::gemm_bt(cols, rows, operands.depth, operands.p + (c0 * operands.ldp), operands.ldp, operands.p + (r0 * operands.ldp),
                                   operands.ldp, out, ldo);
        return;
      }
      if (operands.p_is_qt) {
        numerical::detail::gemm_at(cols, rows, operands.depth, operands.q + c0, operands.ldq, operands.q + r0, operands.ldq, out, ldo);
        return;
      }
      numerical::detail::gemm(cols, rows, operands.depth, operands.p + (c0 * operands.ldp), operands.ldp, operands.q + r0, operands.ldq, out, ldo);
    }
 
    inline void syrk_triangle(syrk_operands operands,
                              std::size_t   j0,
                              std::size_t   w)
    {
      if (w <= syrk_block) {
        buffer<double> diagonal(w * w, 0.0); // row-major: (i, j) is C(j0 + j, j0 + i)
        syrk_product(operands, j0, w, j0, w, diagonal.data(), w);
        for (std::size_t j {0}; j < w; ++j) {
          for (std::size_t i {operands.lower ? j : 0}; i < (operands.lower ? w : j + 1); ++i) {
            operands.c[(j0 + i) + ((j0 + j) * operands.ldc)] += diagonal[(j * w) + i];
          }
        }
        return;
      }
      const std::size_t h    {w / 2};
      syrk_triangle(operands, j0, h);
      const std::size_t r0   {operands.lower ? j0 + h : j0}; // below the first half (L), or right of it (U)
      const std::size_t c0   {operands.lower ? j0 : j0 + h};
      const std::size_t rows {operands.lower ? w - h : h};
      const std::size_t cols {operands.lower ? h : w - h};
      syrk_product(operands, c0, cols, r0, rows, operands.c + r0 + (c0 * operands.ldc), operands.ldc);
      syrk_triangle(operands, j0 + h, w - h);
    }
 
    inline void syrk_all(const syrk_operands& operands,
                         std::size_t          n)
    {
      const std::size_t k {operands.depth};
      if (k < numerical::detail::gemm_min_depth || n * n < ((numerical::detail::gemm_packed_threshold + k) - 1) / k) {
        syrk_triangle(operands, 0, n);
        return;
      }
      const bool m_lower {!operands.lower};
      if (operands.q_is_pt) {
        numerical::detail::gemm_triangle<numerical::detail::a_rows, numerical::detail::b_transposed>(
          n, k, operands.p, operands.ldp, operands.p, operands.ldp, operands.c, operands.ldc, m_lower);
        return;
      }
      if (operands.p_is_qt) {
        numerical::detail::gemm_triangle<numerical::detail::a_transposed, numerical::detail::b_rows>(
          n, k, operands.q, operands.ldq, operands.q, operands.ldq, operands.c, operands.ldc, m_lower);
        return;
      }
      numerical::detail::gemm_triangle<numerical::detail::a_rows, numerical::detail::b_rows>(
        n, k, operands.p, operands.ldp, operands.q, operands.ldq, operands.c, operands.ldc, m_lower);
    }
  } // namespace detail
 
  inline int syrk(char           uplo,
                  char           trans,
                  std::ptrdiff_t n,
                  std::ptrdiff_t k,
                  double         alpha,
                  const double*  a,
                  std::ptrdiff_t lda,
                  double         beta,
                  double*        c,
                  std::ptrdiff_t ldc)
  {
    const bool lower      {uplo == 'L' || uplo == 'l'};
    const bool transposed {detail::transposes(trans)};
    if (!lower && uplo != 'U' && uplo != 'u') {
      return 1;
    }
    if (!detail::is_operation(trans)) {
      return 2;
    }
    if (n < 0) {
      return 3;
    }
    if (k < 0) {
      return 4;
    }
    if (lda < std::max<std::ptrdiff_t>(1, transposed ? k : n)) {
      return 7;
    }
    if (ldc < std::max<std::ptrdiff_t>(1, n)) {
      return 10;
    }
    if (n == 0) {
      return 0;
    }
    const auto order {static_cast<std::size_t>(n)};
    const auto ld    {static_cast<std::size_t>(ldc)};
    if (alpha == 0.0 || k == 0) { // C ← β C on the triangle
      for (std::size_t j {0}; j < order && beta != 1.0; ++j) {
        for (std::size_t i {lower ? j : 0}; i < (lower ? order : j + 1); ++i) {
          c[i + (j * ld)] = beta == 0.0 ? 0.0 : beta * c[i + (j * ld)];
        }
      }
      return 0;
    }
    // C(r, c) = α Σ_l op(A)(r, l) op(A)(c, l). The kernel takes row-major operands: P = op(A) (n × k) and
    // Q = op(A)ᵀ (k × n). A column-major A is one of them read in place (Q for N, P for T). With α = 1 the other is
    // read from it as its transpose (detail::syrk_product; NumPy's a @ a.T is T): no copy, the same bits. Otherwise
    // the other is copied once, α folded into it.
    const auto depth {static_cast<std::size_t>(k)};
    const auto lda_s {static_cast<std::size_t>(lda)};
    for (std::size_t j {0}; j < order; ++j) { // C ← β C on the triangle; the products then add into it
      for (std::size_t i {lower ? j : 0}; i < (lower ? order : j + 1); ++i) {
        double& entry {c[i + (j * ld)]};
        entry = beta == 0.0 ? 0.0 : beta * entry;
      }
    }
    if (alpha == 1.0) {
      detail::syrk_all({transposed ? a : nullptr, lda_s, transposed ? nullptr : a, lda_s, depth, c, ld, lower, transposed, !transposed}, order);
      return 0;
    }
    // The copy in storage the thread keeps between calls (as the kernel's packing_storage): a fresh one of a megabyte
    // is mapped anew and faulted in page by page on every call.
    buffer<double>    copy      {std::exchange(detail::syrk_operand(), buffer<double> {})};
    const std::size_t copy_rows {transposed ? depth : order}; // row-major: α A as it lies, k × n (T) or n × k (N)
    const std::size_t copy_cols {transposed ? order : depth};
    if (copy.size() < copy_rows * copy_cols) {
      copy = buffer<double>(copy_rows * copy_cols);
    }
    for (std::size_t i {0}; i < copy_rows; ++i) {
      for (std::size_t j {0}; j < copy_cols; ++j) {
        copy[(i * copy_cols) + j] = alpha * a[i + (j * lda_s)];
      }
    }
    const double* const p_rows {transposed ? a : copy.data()};
    const std::size_t   ldp    {transposed ? lda_s : depth};
    const double* const q_cols {transposed ? copy.data() : a};
    const std::size_t   ldq    {transposed ? order : lda_s};
    detail::syrk_all({p_rows, ldp, q_cols, ldq, depth, c, ld, lower, false, false}, order);
    if (copy.size() <= syrk_kept_limit) {
      detail::syrk_operand() = std::move(copy);
    }
    return 0;
  }
 
  inline constexpr std::size_t trsm_block {64};
 
  namespace detail {
 
    inline void triangular_solve(bool          lower,
                                 bool          unit,
                                 std::size_t   k,
                                 std::size_t   n,
                                 const double* a,
                                 std::size_t   lda,
                                 double*       b,
                                 std::size_t   ldb)
    {
      for (std::size_t done {0}; done < k; done += trsm_block) {
        const std::size_t w     {std::min(trsm_block, k - done)};
        const std::size_t k0    {lower ? done : k - done - w};
        const std::size_t k1    {k0 + w};
        // The block's triangle, each column of B alone, the columns in parallel.
        numerical::detail::parallel_for(n, w * w, [a, b, lda, ldb, k0, k1, lower, unit](std::size_t j) {
                                          double* x {b + (j * ldb)};
                                          for (std::size_t s {0}; s < k1 - k0; ++s) {
                                            const std::size_t p {lower ? k0 + s : k1 - 1 - s};
                                            if (!unit) {
                                              x[p] /= a[p + (p * lda)];
                                            }
                                            const std::size_t from {lower ? p + 1 : k0};
                                            const std::size_t to {lower ? k1 : p};
                                            numerical::detail::add_product(x + from, -x[p], a + from + (p * lda), to - from);
                                          }
                                        });
        const std::size_t first {lower ? k1 : 0};
        const std::size_t count {lower ? k - k1 : k0};
        if (count > 0) { // B(rows still to solve) −= A(those rows, block) X(block)
          (void)gemm('N', 'N', static_cast<std::ptrdiff_t>(count), static_cast<std::ptrdiff_t>(n), static_cast<std::ptrdiff_t>(w), -1.0,
                     a + first + (k0 * lda), static_cast<std::ptrdiff_t>(lda), b + k0, static_cast<std::ptrdiff_t>(ldb), 1.0, b + first,
                     static_cast<std::ptrdiff_t>(ldb));
        }
      }
    }
  } // namespace detail
 
  inline int trsm(char           side,
                  char           uplo,
                  char           transa,
                  char           diag,
                  std::ptrdiff_t m,
                  std::ptrdiff_t n,
                  double         alpha,
                  const double*  a,
                  std::ptrdiff_t lda,
                  double*        b,
                  std::ptrdiff_t ldb)
  {
    const bool left  {side == 'L' || side == 'l'};
    const bool lower {uplo == 'L' || uplo == 'l'};
    const bool unit  {diag == 'U' || diag == 'u'};
    if (!left && side != 'R' && side != 'r') {
      return 1;
    }
    if (!lower && uplo != 'U' && uplo != 'u') {
      return 2;
    }
    if (!detail::is_operation(transa)) {
      return 3;
    }
    if (!unit && diag != 'N' && diag != 'n') {
      return 4;
    }
    if (m < 0) {
      return 5;
    }
    if (n < 0) {
      return 6;
    }
    if (lda < std::max<std::ptrdiff_t>(1, left ? m : n)) {
      return 9;
    }
    if (ldb < std::max<std::ptrdiff_t>(1, m)) {
      return 11;
    }
    if (m == 0 || n == 0) {
      return 0;
    }
    const auto rows {static_cast<std::size_t>(m)};
    const auto cols {static_cast<std::size_t>(n)};
    const auto ld   {static_cast<std::size_t>(ldb)};
    for (std::size_t j {0}; j < cols; ++j) {
      for (std::size_t i {0}; i < rows; ++i) {
        double& entry {b[i + (j * ld)]};
        entry = alpha == 0.0 ? 0.0 : alpha * entry;
      }
    }
    if (alpha == 0.0) {
      return 0;
    }
    const std::size_t k            {left ? rows : cols};
    const bool        flip_a       {left == detail::transposes(transa)}; // A' = Aᵀ: op(A) = Aᵀ on the left, A on the right
    buffer<double>    transposed_a;
    const double*     a_used       {a};
    std::size_t       lda_used     {static_cast<std::size_t>(lda)};
    if (flip_a) {
      transposed_a = buffer<double>(k * k, 0.0);
      for (std::size_t j {0}; j < k; ++j) {
        for (std::size_t i {lower ? j : 0}; i < (lower ? k : j + 1); ++i) { // A's triangle, transposed
          transposed_a[j + (i * k)] = a[i + (j * static_cast<std::size_t>(lda))];
        }
      }
      a_used   = transposed_a.data();
      lda_used = k;
    }
    if (left) {
      detail::triangular_solve(lower != flip_a, unit, k, cols, a_used, lda_used, b, ld);
      return 0;
    }
    buffer<double> bt(cols * rows, 0.0); // Bᵀ, column-major n × m
    for (std::size_t j {0}; j < cols; ++j) {
      for (std::size_t i {0}; i < rows; ++i) {
        bt[j + (i * cols)] = b[i + (j * ld)];
      }
    }
    detail::triangular_solve(lower != flip_a, unit, k, rows, a_used, lda_used, bt.data(), cols);
    for (std::size_t j {0}; j < cols; ++j) {
      for (std::size_t i {0}; i < rows; ++i) {
        b[i + (j * ld)] = bt[j + (i * cols)];
      }
    }
    return 0;
  }
 
  template <typename T, typename Index>
  int getrs(char           trans,
            std::ptrdiff_t n,
            std::ptrdiff_t nrhs,
            const T*       a,
            std::ptrdiff_t lda,
            const Index*   ipiv,
            T*             b,
            std::ptrdiff_t ldb);
 
  template <typename T, typename Index>
  int gesv(std::ptrdiff_t n,
           std::ptrdiff_t nrhs,
           T*             a,
           std::ptrdiff_t lda,
           Index*         ipiv,
           T*             b,
           std::ptrdiff_t ldb)
  {
    if (n < 0) {
      return -1;
    }
    if (nrhs < 0) {
      return -2;
    }
    if (lda < std::max<std::ptrdiff_t>(1, n)) {
      return -4;
    }
    if (ldb < std::max<std::ptrdiff_t>(1, n)) {
      return -7;
    }
    if (n == 0) {
      return 0;
    }
    const auto               size {static_cast<std::size_t>(n)};
    const auto               ld   {static_cast<std::size_t>(lda)};
    std::vector<std::size_t> pivots;
    matrix<T>                factors;
    bool                     singular {false};
    try { // B and A each transposed in and out once, in storage the thread keeps: the LU solves in B's copy and gives
          // its factors up
      buffer<T> a_rows       {detail::take_kept<T>(0)};
      detail::row_major_into(a_rows, a, lda, n, n, T {1});
      lu_decomposition<T> lu {matrix<T> {size, size, std::move(a_rows)}};
      pivots   = lu.pivots();
      singular = lu.singular();
      if (!singular && nrhs > 0) {
        // B = I (NumPy's inv) as the LU's inverse, its forward pass without the identity's zeros: the same terms
        // otherwise, so the same bits, where the right-hand sides are too many to go row by row.
        const bool inverting {nrhs == n && size > numerical::detail::narrow_solve && detail::is_identity(b, static_cast<std::size_t>(ldb), size)};
        buffer<T>  b_rows    {detail::take_kept<T>(1)};
        if (!inverting) {
          detail::row_major_into(b_rows, b, ldb, n, nrhs, T {1});
        }
        // (P A)⁻¹'s columns moved to A⁻¹'s as they are copied out.
        matrix<T> x {inverting ? lu.pivoted_inverse(std::move(b_rows)) : lu.solve(matrix<T> {size, static_cast<std::size_t>(nrhs), std::move(b_rows)})};
        if (inverting) {
          detail::column_major_permuted_into(&*x.begin(), size, lu.pivots().data(), b, static_cast<std::size_t>(ldb));
        }
        else {
          detail::column_major_into(&*x.begin(), size, static_cast<std::size_t>(nrhs), b, static_cast<std::size_t>(ldb));
        }
        detail::give_back(1, std::move(x).elements());
      }
      factors = std::move(lu).factors();
    }
    catch (const non_finite_error&) { // A and B untouched: dgetf2's factors, then getrs, the NaN carried through
      const int info {detail::rectangular_lu(size, size, a, ld, ipiv)};
      if (info == 0) {
        (void)getrs('N', n, nrhs, a, lda, ipiv, b, ldb); // its arguments already checked
      }
      return info;
    }
    detail::column_major_into(&*factors.begin(), size, size, a, ld);
    detail::give_back(0, std::move(factors).elements());
    detail::interchanges(pivots, ipiv);
    for (std::size_t i {0}; i < size; ++i) {
      if (a[i + (i * ld)] == T {}) { // U(i, i), as just written
        return static_cast<int>(i + 1);
      }
    }
    return 0;
  }
 
  namespace detail {
    template <typename T>
    int column_cholesky(bool        lower,
                        std::size_t size,
                        T*          a,
                        std::size_t ld)
    {
      // Element (i, j) of L: below the diagonal for lower, the conjugate of U's (j, i) above it otherwise.
      const auto l = [&](std::size_t i, std::size_t j) -> T& { return lower ? a[i + (j * ld)] : a[j + (i * ld)]; };
      const auto from_l = [&](T x) { return lower ? x : numerical::detail::conj_of(x); }; // L's value as stored
      for (std::size_t j {0}; j < size; ++j) {
        auto diagonal {std::real(a[j + (j * ld)])};
        for (std::size_t k {0}; k < j; ++k) {
          const T ljk {from_l(l(j, k))};
          diagonal -= std::real(ljk * numerical::detail::conj_of(ljk));
        }
        if (!(diagonal > 0)) { // not positive, or NaN
          a[j + (j * ld)] = T {diagonal};
          return static_cast<int>(j + 1);
        }
        diagonal        = std::sqrt(diagonal);
        a[j + (j * ld)] = T {diagonal};
        for (std::size_t i {j + 1}; i < size; ++i) {
          T sum {from_l(l(i, j))};
          for (std::size_t k {0}; k < j; ++k) {
            sum -= from_l(l(i, k)) * numerical::detail::conj_of(from_l(l(j, k)));
          }
          l(i, j) = from_l(sum / diagonal);
        }
      }
      return 0;
    }
  } // namespace detail
 
  template <typename T>
  int potrf(char           uplo,
            std::ptrdiff_t n,
            T*             a,
            std::ptrdiff_t lda)
  {
    const bool lower {uplo == 'L' || uplo == 'l'};
    if (!lower && uplo != 'U' && uplo != 'u') {
      return -1;
    }
    if (n < 0) {
      return -2;
    }
    if (lda < std::max<std::ptrdiff_t>(1, n)) {
      return -4;
    }
    if (n == 0) {
      return 0;
    }
    const auto size {static_cast<std::size_t>(n)};
    const auto ld   {static_cast<std::size_t>(lda)};
    // A's lower triangle, row-major (the upper one's conjugate, Hermitian), in storage the thread keeps: the decomposition
    // reads no more of it and zeroes the rest. L's triangle by tiles (a row of it read across columns lda apart missed on
    // every element), U's columns as they lie; the factor written back the same ways.
    buffer<T> lower_rows {detail::take_kept<T>(0)};
    lower_rows.resize(size * size);
    if (lower) {
      detail::lower_triangle_rows(a, ld, size, lower_rows.data(), true);
    }
    else {
      for (std::size_t i {0}; i < size; ++i) {
        for (std::size_t j {0}; j <= i; ++j) {
          lower_rows[(i * size) + j] = numerical::detail::conj_of(a[j + (i * ld)]);
        }
      }
    }
    try {
      buffer<T> l {cholesky_decomposition<T> {matrix<T> {size, size, std::move(lower_rows)}}.lower().elements()};
      if (lower) {
        detail::lower_triangle_rows(a, ld, size, l.data(), false);
      }
      else { // U = Lᴴ
        for (std::size_t i {0}; i < size; ++i) {
          for (std::size_t j {0}; j <= i; ++j) {
            a[j + (i * ld)] = numerical::detail::conj_of(l[(i * size) + j]);
          }
        }
      }
      detail::give_back(0, std::move(l));
    }
    catch (const not_positive_definite_error& failed) {
      return static_cast<int>(failed.order());
    }
    catch (const non_finite_error&) { // A untouched: factorised as the reference's dpotf2 does, the NaN carried through
      return detail::column_cholesky(lower, size, a, ld);
    }
    return 0;
  }
 
  template <typename T, typename Index>
  int getrf(std::ptrdiff_t m,
            std::ptrdiff_t n,
            T*             a,
            std::ptrdiff_t lda,
            Index*         ipiv)
  {
    if (m < 0) {
      return -1;
    }
    if (n < 0) {
      return -2;
    }
    if (lda < std::max<std::ptrdiff_t>(1, m)) {
      return -4;
    }
    if (m == 0 || n == 0) {
      return 0;
    }
    const auto rows {static_cast<std::size_t>(m)};
    const auto cols {static_cast<std::size_t>(n)};
    const auto ld   {static_cast<std::size_t>(lda)};
    if (rows != cols) {
      return detail::rectangular_lu(rows, cols, a, ld, ipiv);
    }
    buffer<T> a_rows {detail::take_kept<T>(0)};  // transposed in and out once, in storage the thread keeps
    detail::row_major_into(a_rows, a, lda, m, n, T {1});
    std::optional<lu_decomposition<T>> factored; // a square A checked once, by the decomposition, where it was also
    try {                             // checked here, entry by entry across its columns
      factored.emplace(matrix<T> {rows, cols, std::move(a_rows)});
    }
    catch (const non_finite_error&) { // A untouched: factorised as the reference's dgetf2 does, the NaN carried through
      return detail::rectangular_lu(rows, cols, a, ld, ipiv);
    }
    lu_decomposition<T>& lu {*factored};
    detail::interchanges(lu.pivots(), ipiv);
    matrix<T> factors       {std::move(lu).factors()};
    detail::column_major_into(&*factors.begin(), rows, cols, a, ld);
    detail::give_back(0, std::move(factors).elements());
    for (std::size_t i {0}; i < rows; ++i) {
      if (a[i + (i * ld)] == T {}) { // U(i, i), as just written
        return static_cast<int>(i + 1);
      }
    }
    return 0;
  }
 
  template <typename T, typename Index>
  int getrs(char           trans,
            std::ptrdiff_t n,
            std::ptrdiff_t nrhs,
            const T*       a,
            std::ptrdiff_t lda,
            const Index*   ipiv,
            T*             b,
            std::ptrdiff_t ldb)
  {
    if (!detail::is_operation(trans)) {
      return -1;
    }
    if (n < 0) {
      return -2;
    }
    if (nrhs < 0) {
      return -3;
    }
    if (lda < std::max<std::ptrdiff_t>(1, n)) {
      return -5;
    }
    if (ldb < std::max<std::ptrdiff_t>(1, n)) {
      return -8;
    }
    const auto size {static_cast<std::size_t>(n)};
    const auto ld   {static_cast<std::size_t>(lda)};
    const bool conjugate {detail::conjugates(trans)};
    const auto at        = [&](std::size_t i, std::size_t j) { return a[i + (j * ld)]; };
    // An element of op(A)'s factor read transposed: conjugated for Aᴴ.
    const auto across    = [&](std::size_t i, std::size_t j) { return conjugate ? numerical::detail::conj_of(at(i, j)) : at(i, j); };
    for (std::size_t c {0}; c < static_cast<std::size_t>(nrhs); ++c) {
      T* const x {b + (c * static_cast<std::size_t>(ldb))};
      if (!detail::transposes(trans)) { // L U x = P b
        for (std::size_t i {0}; i < size; ++i) {
          std::swap(x[i], x[static_cast<std::size_t>(ipiv[i]) - 1]);
        }
        for (std::size_t j {0}; j < size; ++j) { // L, unit lower, by columns
          for (std::size_t i {j + 1}; i < size; ++i) {
            x[i] -= at(i, j) * x[j];
          }
        }
        for (std::size_t j {size}; j-- > 0;) { // U, by columns
          x[j] /= at(j, j);
          for (std::size_t i {0}; i < j; ++i) {
            x[i] -= at(i, j) * x[j];
          }
        }
        continue;
      }
      for (std::size_t i {0}; i < size; ++i) { // Uᵀ (Uᴴ), by rows of U's columns
        T sum {x[i]};
        for (std::size_t k {0}; k < i; ++k) {
          sum -= across(k, i) * x[k];
        }
        x[i] = sum / across(i, i);
      }
      for (std::size_t i {size}; i-- > 0;) { // Lᵀ (Lᴴ)
        T sum {x[i]};
        for (std::size_t k {i + 1}; k < size; ++k) {
          sum -= across(k, i) * x[k];
        }
        x[i] = sum;
      }
      for (std::size_t i {size}; i-- > 0;) { // Pᵀ: the interchanges undone in reverse
        std::swap(x[i], x[static_cast<std::size_t>(ipiv[i]) - 1]);
      }
    }
    return 0;
  }
 
  template <typename T>
  int potrs(char           uplo,
            std::ptrdiff_t n,
            std::ptrdiff_t nrhs,
            const T*       a,
            std::ptrdiff_t lda,
            T*             b,
            std::ptrdiff_t ldb)
  {
    const bool lower {uplo == 'L' || uplo == 'l'};
    if (!lower && uplo != 'U' && uplo != 'u') {
      return -1;
    }
    if (n < 0) {
      return -2;
    }
    if (nrhs < 0) {
      return -3;
    }
    if (lda < std::max<std::ptrdiff_t>(1, n)) {
      return -5;
    }
    if (ldb < std::max<std::ptrdiff_t>(1, n)) {
      return -7;
    }
    const auto size {static_cast<std::size_t>(n)};
    const auto ld   {static_cast<std::size_t>(lda)};
    // L(i, j), i ≥ j: U(j, i) conjugated where U = Lᴴ is stored.
    const auto l = [&](std::size_t i, std::size_t j) { return lower ? a[i + (j * ld)] : numerical::detail::conj_of(a[j + (i * ld)]); };
    for (std::size_t c {0}; c < static_cast<std::size_t>(nrhs); ++c) {
      T* const x {b + (c * static_cast<std::size_t>(ldb))};
      for (std::size_t i {0}; i < size; ++i) { // L y = b
        T sum {x[i]};
        for (std::size_t k {0}; k < i; ++k) {
          sum -= l(i, k) * x[k];
        }
        x[i] = sum / l(i, i);
      }
      for (std::size_t i {size}; i-- > 0;) { // Lᴴ x = y
        T sum {x[i]};
        for (std::size_t k {i + 1}; k < size; ++k) {
          sum -= numerical::detail::conj_of(l(k, i)) * x[k];
        }
        x[i] = sum / numerical::detail::conj_of(l(i, i));
      }
    }
    return 0;
  }
 
  template <typename Index>
  int syevd(char           jobz,
            char           uplo,
            std::ptrdiff_t n,
            double*        a,
            std::ptrdiff_t lda,
            double*        w,
            double*        work,
            std::ptrdiff_t lwork,
            Index*         iwork,
            std::ptrdiff_t liwork)
  {
    const bool vectors {jobz == 'V' || jobz == 'v'};
    const bool lower   {uplo == 'L' || uplo == 'l'};
    if (!vectors && jobz != 'N' && jobz != 'n') {
      return -1;
    }
    if (!lower && uplo != 'U' && uplo != 'u') {
      return -2;
    }
    if (n < 0) {
      return -3;
    }
    if (lda < std::max<std::ptrdiff_t>(1, n)) {
      return -5;
    }
    std::ptrdiff_t lwork_least  {1};
    std::ptrdiff_t liwork_least {1};
    if (n > 1) {
      lwork_least  = vectors ? 1 + (6 * n) + (2 * n * n) : (2 * n) + 1;
      liwork_least = vectors ? 3 + (5 * n) : 1;
    }
    if (lwork == -1 || liwork == -1) {
      work[0]  = static_cast<double>(lwork_least);
      iwork[0] = static_cast<Index>(liwork_least);
      return 0;
    }
    if (lwork < lwork_least) {
      return -8;
    }
    if (liwork < liwork_least) {
      return -10;
    }
    if (n == 0) {
      return 0;
    }
    const auto     size {static_cast<std::size_t>(n)};
    const auto     ld   {static_cast<std::size_t>(lda)};
    buffer<double> lower_rows(size * size, 0.0); // the triangle read, as the lower one, row-major
    for (std::size_t i {0}; i < size; ++i) {
      for (std::size_t j {0}; j <= i; ++j) {
        lower_rows[(i * size) + j] = lower ? a[i + (j * ld)] : a[j + (i * ld)];
      }
    }
    try {
      const matrix<double> symmetric {size, size, std::move(lower_rows)};
      if (!vectors) {
        const vector<double> values {symmetric_eigenvalues(symmetric)};
        std::copy_n(&values[0], size, w);
        return 0;
      }
      const symmetric_eigen eigen {symmetric};
      std::copy_n(&eigen.values()[0], size, w);
      for (std::size_t j {0}; j < size; ++j) {
        for (std::size_t i {0}; i < size; ++i) {
          a[i + (j * ld)] = eigen.vectors()(i, j);
        }
      }
    }
    catch (const non_finite_error&) { // as the reference's divide and conquer on a NaN: no eigenvalue converges
      std::fill_n(w, size, std::numeric_limits<double>::quiet_NaN());
      if (vectors) {
        for (std::size_t j {0}; j < size; ++j) {
          std::fill_n(a + (j * ld), size, std::numeric_limits<double>::quiet_NaN());
        }
      }
      return static_cast<int>(size);
    }
    return 0;
  }
 
  template <typename Index>
  int gesdd(char           jobz,
            std::ptrdiff_t m,
            std::ptrdiff_t n,
            double*        a,
            std::ptrdiff_t lda,
            double*        s,
            double*        u,
            std::ptrdiff_t ldu,
            double*        vt,
            std::ptrdiff_t ldvt,
            double*        work,
            std::ptrdiff_t lwork,
            Index*         iwork)
  {
    (void)iwork; // dgesdd's integer workspace, of 8 min(m, n): not used
    const bool all       {jobz == 'A' || jobz == 'a'};
    const bool some      {jobz == 'S' || jobz == 's'};
    const bool overwrite {jobz == 'O' || jobz == 'o'};
    const bool none      {jobz == 'N' || jobz == 'n'};
    if (!all && !some && !overwrite && !none) {
      return -1;
    }
    if (m < 0) {
      return -2;
    }
    if (n < 0) {
      return -3;
    }
    if (lda < std::max<std::ptrdiff_t>(1, m)) {
      return -5;
    }
    const std::ptrdiff_t least  {std::min(m, n)};
    const std::ptrdiff_t most   {std::max(m, n)};
    const bool           u_out  {all || some || (overwrite && m < n)};
    const bool           vt_out {all || some || (overwrite && m >= n)};
    if (ldu < std::max<std::ptrdiff_t>(1, u_out ? m : 1)) {
      return -8;
    }
    std::ptrdiff_t vt_rows {1};
    if (some) {
      vt_rows = least;
    }
    else if (vt_out) { // A, or O with m ≥ n: Vᵀ whole
      vt_rows = n;
    }
    if (ldvt < std::max<std::ptrdiff_t>(1, vt_rows)) {
      return -10;
    }
    std::ptrdiff_t lwork_least {(3 * least) + std::max(most, 7 * least)};
    if (overwrite) {
      lwork_least = (3 * least) + std::max(most, (5 * least * least) + (4 * least));
    }
    else if (all || some) {
      lwork_least = (3 * least) + std::max(most, (4 * least * least) + (4 * least));
    }
    if (lwork == -1) {
      work[0] = static_cast<double>(std::max<std::ptrdiff_t>(1, lwork_least));
      return 0;
    }
    if (lwork < std::max<std::ptrdiff_t>(1, lwork_least)) {
      return -12;
    }
    if (least == 0) {
      return 0;
    }
    const auto rows {static_cast<std::size_t>(m)};
    const auto cols {static_cast<std::size_t>(n)};
    const auto k    {static_cast<std::size_t>(least)};
    try {
      const matrix<double> source {rows, cols, detail::row_major(a, lda, m, n, 1.0)};
      if (none) {
        const vector<double> values {singular_values(source)};
        std::copy_n(&values[0], k, s);
        return 0;
      }
      const svd_decomposition<double> svd {source, svd_method::divide_and_conquer, all};
      std::copy_n(&svd.singular_values()[0], k, s);
      const matrix<double> left           {svd.u()};
      const matrix<double> right          {svd.vh()};
      // With O, the factor that fits A's storage goes there, the other to its own argument.
      double* const        u_into         {overwrite && m >= n ? a : u};
      const std::size_t    ld_u           {overwrite && m >= n ? static_cast<std::size_t>(lda) : static_cast<std::size_t>(ldu)};
      double* const        vt_into        {overwrite && m < n ? a : vt};
      const std::size_t    ld_vt          {overwrite && m < n ? static_cast<std::size_t>(lda) : static_cast<std::size_t>(ldvt)};
      for (std::size_t j {0}; j < left.cols(); ++j) {
        for (std::size_t i {0}; i < left.rows(); ++i) {
          u_into[i + (j * ld_u)] = left(i, j);
        }
      }
      for (std::size_t j {0}; j < right.cols(); ++j) {
        for (std::size_t i {0}; i < right.rows(); ++i) {
          vt_into[i + (j * ld_vt)] = right(i, j);
        }
      }
    }
    catch (const non_finite_error&) {
      return -4;
    }
    return 0;
  }
 
  inline int geqrf(std::ptrdiff_t m,
                   std::ptrdiff_t n,
                   double*        a,
                   std::ptrdiff_t lda,
                   double*        tau,
                   double*        work,
                   std::ptrdiff_t lwork)
  {
    if (m < 0) {
      return -1;
    }
    if (n < 0) {
      return -2;
    }
    if (lda < std::max<std::ptrdiff_t>(1, m)) {
      return -4;
    }
    const std::ptrdiff_t lwork_least {std::min(m, n) == 0 ? 1 : n};
    if (lwork == -1) {
      work[0] = static_cast<double>(lwork_least);
      return 0;
    }
    if (lwork < lwork_least) {
      return -7;
    }
    if (std::min(m, n) == 0) {
      return 0;
    }
    const auto     rows    {static_cast<std::size_t>(m)};
    const auto     cols    {static_cast<std::size_t>(n)};
    const auto     ld      {static_cast<std::size_t>(lda)};
    // In the layout the factorisation runs in (numerical::detail::householder_reflect), in storage the thread keeps: R and
    // the reflectors' w₁… where dgeqrf leaves them (numerical::detail::factor_qr_panel, as dlarfg); τ on each block's T
    // diagonal, 0 for the identity.
    buffer<double> entries {detail::take_kept<double>(0)};
#ifdef NUMERICAL_X86_COMPILED
    entries.resize(rows * cols); // column-major: copied in and out by columns
    for (std::size_t j {0}; j < cols; ++j) {
      std::copy_n(a + (j * ld), rows, entries.data() + (j * rows));
    }
#else
    detail::row_major_into(entries, a, lda, m, n, 1.0);
#endif
#ifdef NUMERICAL_X86_COMPILED
    const std::vector<numerical::detail::reflector_block<double>> blocks {numerical::detail::reflect_in_columns(entries, rows, cols)};
    for (std::size_t j {0}; j < cols; ++j) {
      std::copy_n(entries.data() + (j * rows), rows, a + (j * ld));
    }
#else
    const std::vector<numerical::detail::reflector_block<double>> blocks {numerical::detail::reflect_in_rows(entries, rows, cols)};
    detail::column_major_into(entries.data(), rows, cols, a, ld);
#endif
    detail::give_back(0, std::move(entries));
    for (const numerical::detail::reflector_block<double>& block : blocks) {
      for (std::size_t j {0}; j < block.width; ++j) {
        tau[block.first + j] = block.t[(j * block.width) + j];
      }
    }
    return 0;
  }
 
  inline int orgqr(std::ptrdiff_t m,
                   std::ptrdiff_t n,
                   std::ptrdiff_t k,
                   double*        a,
                   std::ptrdiff_t lda,
                   const double*  tau,
                   double*        work,
                   std::ptrdiff_t lwork)
  {
    if (m < 0) {
      return -1;
    }
    if (n < 0 || n > m) {
      return -2;
    }
    if (k < 0 || k > n) {
      return -3;
    }
    if (lda < std::max<std::ptrdiff_t>(1, m)) {
      return -5;
    }
    const std::ptrdiff_t lwork_least {std::max<std::ptrdiff_t>(1, n)};
    if (lwork == -1) {
      work[0] = static_cast<double>(lwork_least);
      return 0;
    }
    if (lwork < lwork_least) {
      return -8;
    }
    if (n == 0) {
      return 0;
    }
    const auto rows  {static_cast<std::size_t>(m)};
    const auto cols  {static_cast<std::size_t>(n)};
    const auto depth {static_cast<std::size_t>(k)};
    const auto ld    {static_cast<std::size_t>(lda)};
    // Each reflector w_j = (1, A(j+1…m-1, j)) and τ_j, gathered by blocks into I − V T Vᵀ (T as dlarft forms it).
    std::vector<buffer<double>> reflectors(depth);
    const std::vector<double>   taus(tau, tau + depth);
    for (std::size_t j {0}; j < depth; ++j) {
      reflectors[j].resize(rows - j);
      reflectors[j][0] = 1.0;
      std::copy_n(a + (j + 1) + (j * ld), rows - j - 1, reflectors[j].data() + 1);
    }
    std::vector<numerical::detail::reflector_block<double>> blocks;
    for (std::size_t k0 {0}; k0 < depth; k0 += numerical::detail::qr_block) {
      blocks.push_back(numerical::detail::gather_reflectors(reflectors, k0, std::min(k0 + numerical::detail::qr_block, depth), rows, &taus));
    }
    buffer<double> columns {detail::take_kept<double>(0)}; // Q column-major, as numerical::qr() forms it
    numerical::detail::form_q_into(columns, blocks, rows, 0, cols, true);
    for (std::size_t j {0}; j < cols; ++j) {
      std::copy_n(columns.data() + (j * rows), rows, a + (j * ld));
    }
    detail::give_back(0, std::move(columns));
    return 0;
  }
 
  template <typename Index>
  int gelsd(std::ptrdiff_t m,
            std::ptrdiff_t n,
            std::ptrdiff_t nrhs,
            const double*  a,
            std::ptrdiff_t lda,
            double*        b,
            std::ptrdiff_t ldb,
            double*        s,
            double         rcond,
            Index*         rank,
            double*        work,
            std::ptrdiff_t lwork,
            Index*         iwork)
  {
    if (m < 0) {
      return -1;
    }
    if (n < 0) {
      return -2;
    }
    if (nrhs < 0) {
      return -3;
    }
    if (lda < std::max<std::ptrdiff_t>(1, m)) {
      return -5;
    }
    if (ldb < std::max<std::ptrdiff_t>({1, m, n})) {
      return -7;
    }
    // dgelsd's sizes: its divide and conquer's leaves of 25 values, NLVL levels above them.
    constexpr std::ptrdiff_t leaf    {25};
    const std::ptrdiff_t     least   {std::max<std::ptrdiff_t>(1, std::min(m, n))};
    const auto               levels  {std::max<std::ptrdiff_t>(
                                        static_cast<std::ptrdiff_t>(std::log2(static_cast<double>(least) / static_cast<double>(leaf + 1))) + 1, 0)};
    const std::ptrdiff_t lwork_least {(3 * least) + (9 * least) + (2 * least * leaf) + (8 * least * levels) + (least * nrhs)
                                      + ((leaf + 1) * (leaf + 1))};
    if (lwork == -1) {
      work[0]  = static_cast<double>(lwork_least);
      iwork[0] = static_cast<Index>((3 * least * levels) + (11 * least));
      return 0;
    }
    if (lwork < lwork_least) {
      return -12;
    }
    *rank = 0;
    if (m == 0 || n == 0) {
      return 0;
    }
    const auto     rows {static_cast<std::size_t>(m)};
    const auto     cols {static_cast<std::size_t>(n)};
    const auto     rhs  {static_cast<std::size_t>(nrhs)};
    const auto     ld   {static_cast<std::size_t>(ldb)};
    buffer<double> right(rows * rhs, 0.0);
    for (std::size_t i {0}; i < rows; ++i) {
      for (std::size_t c {0}; c < rhs; ++c) {
        right[(i * rhs) + c] = b[i + (c * ld)];
      }
    }
    if (!std::ranges::all_of(right, [](double x) { return std::isfinite(x); })) {
      return -6;
    }
    const rank_cutoff cutoff {rank_cutoff::relative(rcond > 0.0 && rcond < 1.0 ? rcond : lapack_precision)};
    try {
      const matrix<double>                    left {rows, cols, detail::row_major(a, lda, m, n, 1.0)};
      const least_squares_solution<double>    fit  {least_squares(left, matrix<double> {rows, rhs, std::move(right)}, cutoff)};
      std::copy_n(&fit.singular_values[0], fit.singular_values.size(), s);
      *rank = static_cast<Index>(fit.rank);
      for (std::size_t c {0}; c < rhs; ++c) {
        for (std::size_t i {0}; i < cols; ++i) {
          b[i + (c * ld)] = fit.solution(i, c);
        }
        for (std::size_t i {cols}; i < rows; ++i) {
          b[i + (c * ld)] = i == cols && fit.residuals.size() > 0 ? std::sqrt(fit.residuals[c]) : 0.0;
        }
      }
    }
    catch (const non_finite_error&) {
      return -4;
    }
    return 0;
  }
 
  namespace detail {
 
    inline void real_form(const matrix<std::complex<double>>&       v,
                          const vector<std::complex<double>>&       values,
                          const std::vector<std::size_t>&           column_of,
                          bool                                      conjugate,
                          double*                                   out,
                          std::size_t                               ld)
    {
      const std::size_t n    {values.size()};
      const double      sign {conjugate ? -1.0 : 1.0};
      for (std::size_t j {0}; j < n; ++j) {
        const std::size_t c {column_of[j]};
        if (values[j].imag() == 0.0) {
          for (std::size_t i {0}; i < n; ++i) {
            out[i + (j * ld)] = v(i, c).real();
          }
        }
        else if (values[j].imag() > 0.0) {
          for (std::size_t i {0}; i < n; ++i) {
            out[i + (j * ld)]       = v(i, c).real();
            out[i + ((j + 1) * ld)] = sign * v(i, c).imag();
          }
        }
      }
    }
  } // namespace detail
 
  inline int geev(char           jobvl,
                  char           jobvr,
                  std::ptrdiff_t n,
                  const double*  a,
                  std::ptrdiff_t lda,
                  double*        wr,
                  double*        wi,
                  double*        vl,
                  std::ptrdiff_t ldvl,
                  double*        vr,
                  std::ptrdiff_t ldvr,
                  double*        work,
                  std::ptrdiff_t lwork)
  {
    const bool left  {jobvl == 'V' || jobvl == 'v'};
    const bool right {jobvr == 'V' || jobvr == 'v'};
    if (!left && jobvl != 'N' && jobvl != 'n') {
      return -1;
    }
    if (!right && jobvr != 'N' && jobvr != 'n') {
      return -2;
    }
    if (n < 0) {
      return -3;
    }
    if (lda < std::max<std::ptrdiff_t>(1, n)) {
      return -5;
    }
    if (ldvl < 1 || (left && ldvl < n)) {
      return -9;
    }
    if (ldvr < 1 || (right && ldvr < n)) {
      return -11;
    }
    const std::ptrdiff_t lwork_least {std::max<std::ptrdiff_t>(1, (left || right ? 4 : 3) * n)};
    if (lwork == -1) {
      work[0] = static_cast<double>(lwork_least);
      return 0;
    }
    if (lwork < lwork_least) {
      return -13;
    }
    if (n == 0) {
      return 0;
    }
    const auto size {static_cast<std::size_t>(n)};
    try {
      const matrix<double>                      source {size, size, detail::row_major(a, lda, n, n, 1.0)};
      const eigen_decomposition<double>         eigen  {source, right};
      const vector<std::complex<double>>&       values {eigen.values()};
      for (std::size_t j {0}; j < size; ++j) {
        wr[j] = values[j].real();
        wi[j] = values[j].imag();
      }
      std::vector<std::size_t> same(size);
      for (std::size_t j {0}; j < size; ++j) {
        same[j] = j;
      }
      if (right) {
        detail::real_form(eigen.vectors(), values, same, false, vr, static_cast<std::size_t>(ldvr));
      }
      if (left) {
        const eigen_decomposition<double> transposed {matrix<double> {size, size, detail::columns_as_rows(a, lda, n, n, 1.0)}, true};
        // Aᵀ w = λ w is (conj w)ᴴ A = λ (conj w)ᴴ: each eigenvalue's left vector is the conjugate of Aᵀ's right one
        // for the same eigenvalue, the nearest one not yet taken (a multiple eigenvalue takes its vectors in turn).
        std::vector<std::size_t> column_of(size);
        std::vector<bool>        taken(size, false);
        for (std::size_t j {0}; j < size; ++j) {
          std::size_t best {size};
          for (std::size_t p {0}; p < size; ++p) {
            if (!taken[p] && (best == size || std::abs(transposed.values()[p] - values[j]) < std::abs(transposed.values()[best] - values[j]))) {
              best = p;
            }
          }
          taken[best]  = true;
          column_of[j] = best;
        }
        detail::real_form(transposed.vectors(), values, column_of, true, vl, static_cast<std::size_t>(ldvl));
      }
    }
    catch (const non_finite_error&) {
      return -4;
    }
    return 0;
  }
} // namespace numerical::blas
 
#endif // NUMERICAL_BLAS_HPP