decomposition.hpp

Go to the documentation of this file.

 
#ifndef NUMERICAL_DECOMPOSITION_HPP
#define NUMERICAL_DECOMPOSITION_HPP
 
#include <algorithm>
#include <array>
#include <atomic>
#include <bit>
#include <cmath>
#include <complex> // so std::abs resolves for a complex scalar, regardless of include order
#include <cstddef>
#include <cstdint>
#include <cstring>
#include <functional>
#include <limits>
#include <memory>
#include <numbers>
#include <optional>
#include <stdexcept>
#include <string>
#include <tuple>
#include <type_traits>
#include <utility>
#include <vector>
 
#include "matrix.hpp"
#include "simd.hpp"
#include "parallel.hpp"
#include "vector.hpp"
 
namespace numerical {
 
  namespace detail {
    template <typename T>
    [[nodiscard]] T pivot_divisor(const T& pivot)
    {
      return pivot;
    }
 
    template <typename R>
    [[nodiscard]] std::complex<R> pivot_divisor(const std::complex<R>& pivot)
    {
      return std::complex<R> {1} / pivot;
    }
 
    template <typename T>
    [[nodiscard]] T over_pivot(const T& x,
                               const T& divisor)
    {
      return x / divisor;
    }
 
    template <typename R>
    [[nodiscard]] std::complex<R> over_pivot(const std::complex<R>& x,
                                             const std::complex<R>& divisor)
    {
      return operator_times(x, divisor);
    }
 
    template <typename T>
    void subtract_term(T&       sum,
                       const T& a,
                       const T& b)
    {
      sum -= a * b;
    }
 
    template <typename R>
    void subtract_term(std::complex<R>&       sum,
                       const std::complex<R>& a,
                       const std::complex<R>& b)
    {
      sum -= operator_times(a, b);
    }
 
    inline constexpr std::size_t lu_block {64};
 
#ifdef NUMERICAL_X86_COMPILED
    inline constexpr std::size_t lu_small_order {384};
#else
    inline constexpr std::size_t lu_small_order {256};
#endif
 
    inline constexpr std::size_t lu_small_block {32};
 
    [[nodiscard]] constexpr std::size_t lu_block_for(std::size_t n) noexcept
    {
      static_assert(std::has_single_bit(lu_small_block) && std::has_single_bit(lu_block), "a panel is a mask away");
      return n <= lu_small_order ? lu_small_block : lu_block;
    }
 
    inline constexpr std::size_t lu_lookahead_columns {1024};
 
#ifdef NUMERICAL_X86_COMPILED
    inline constexpr std::size_t lu_update_speedup {gemm_kernel_speedup / 4};
#else
    inline constexpr std::size_t lu_update_speedup {gemm_kernel_speedup};
#endif
 
    inline constexpr std::size_t narrow_solve {16};
 
#ifdef NUMERICAL_X86_COMPILED
    inline constexpr std::size_t slice_columns {48};
#else
    inline constexpr std::size_t slice_columns {64};
#endif
 
#ifdef NUMERICAL_X86_COMPILED
    inline constexpr std::size_t slice_solve_order {384};
#endif
 
 
    inline constexpr std::size_t solve_leaf {16};
 
    inline constexpr std::size_t inverse_band {16};
 
    inline constexpr std::size_t memory_element_weight {4};
 
    inline constexpr std::size_t svd_band {32};
 
    inline constexpr std::size_t zeroed_rows {32};
 
#ifdef NUMERICAL_X86_COMPILED
    inline constexpr std::size_t panel_copy_rows {8};
 
    inline constexpr std::size_t panel_tile_rows {256};
#endif
 
    inline constexpr std::size_t cholesky_block {64};
 
    inline constexpr std::size_t cholesky_diagonal_rows {4};
 
    inline constexpr std::size_t cholesky_together_width {16};
 
#ifdef NUMERICAL_X86_COMPILED
    inline constexpr std::size_t cholesky_one_panel {40};
#else
    inline constexpr std::size_t cholesky_one_panel {48};
#endif
    static_assert(cholesky_one_panel <= cholesky_block, "one panel of cholesky_block columns holds the whole matrix");
 
    [[nodiscard]] constexpr std::size_t cholesky_block_for(std::size_t n) noexcept
    {
      if (n < cholesky_one_panel) {
        return cholesky_block;
      }
      if (n <= 128) {
        return cholesky_block / 4;
      }
      return n <= 256 ? cholesky_block / 2 : cholesky_block;
    }
 
    inline constexpr std::size_t cholesky_rows {64};
 
    inline constexpr std::size_t cholesky_update_columns {64};
 
    inline constexpr std::size_t cholesky_update_rows {256};
 
 
    // conj_of / real_part are written as overloads (not `if constexpr`) so each
    // is a single, fully-covered expression: a complex scalar selects the
    // std::complex overload, any other scalar the primary template.
 
    template <typename T>
    [[nodiscard]] T conj_of(const T& x)
    {
      return x;
    }
 
    template <typename U>
    [[nodiscard]] std::complex<U> conj_of(const std::complex<U>& x)
    {
      return std::conj(x);
    }
 
    template <typename T>
    using magnitude_t = decltype(std::abs(std::declval<const T&>()));
 
    template <typename T>
    [[nodiscard]] T from_double(double x)
    {
      return static_cast<T>(static_cast<magnitude_t<T>>(x));
    }
 
    template <typename T>
    [[nodiscard]] double real_part(const T& x)
    {
      return static_cast<double>(x);
    }
 
    template <typename U>
    [[nodiscard]] double real_part(const std::complex<U>& x)
    {
      return static_cast<double>(x.real());
    }
 
    template <typename T>
    using value_range_t =
      std::conditional_t<std::is_same_v<T, float> || std::is_same_v<T, std::complex<float>>, float, double>;
 
    template <typename T>
    inline constexpr double rescue_below {std::is_same_v<value_range_t<T>, float> ? 0x1p-60 : 0x1p-500};
 
    template <typename T>
    [[nodiscard]] double rescue_factor(double largest)
    {
      constexpr double least {static_cast<double>(std::numeric_limits<value_range_t<T>>::min())};
      return largest < rescue_below<T> ? std::ldexp(1.0, -std::ilogb(std::max(largest, least))) : 1.0;
    }
 
    template <typename T>
    [[nodiscard]] T phase_of(const T& x)
    {
      const T      scaled    {x * detail::from_double<T>(rescue_factor<T>(static_cast<double>(std::abs(x))))};
      const double magnitude {static_cast<double>(std::abs(scaled))};
      return magnitude == 0.0 ? T {1} : scaled / detail::from_double<T>(magnitude);
    }
 
    template <typename T>
    double rescale_tiny(buffer<T>& x)
    {
      double largest {0.0};
      for (const T& element : x) {
        largest = std::max(largest, static_cast<double>(std::abs(element)));
      }
      const double factor {rescue_factor<T>(largest)};
      for (T& element : x) {
        element *= detail::from_double<T>(factor);
      }
      return factor;
    }
 
    template <typename T>
    [[nodiscard]] double norm_at_scale(const T*    x,
                                       std::size_t count)
    {
      double largest {0.0};
      for (std::size_t i {0}; i < count; ++i) {
        largest = std::max(largest, static_cast<double>(std::abs(x[i])));
      }
      if (largest == 0.0) { // the entries are finite: require_finite runs where every caller starts
        return largest;
      }
      double sum {0.0};
      for (std::size_t i {0}; i < count; ++i) {
        const double ratio {static_cast<double>(std::abs(x[i])) / largest};
        sum += ratio * ratio;
      }
      return largest * std::sqrt(sum);
    }
 
    inline constexpr double plain_squares_low {0x1p-900};
 
    inline constexpr double plain_squares_high {0x1p900};
 
    template <typename T>
    [[nodiscard]] double plain_or_at_scale(double      plain,
                                           const T*    x,
                                           std::size_t count)
    {
      if (plain >= plain_squares_low && plain <= plain_squares_high) {
        return std::sqrt(plain);
      }
      return norm_at_scale(x, count);
    }
 
    template <typename T>
    [[nodiscard]] double scaled_norm(const T*    x,
                                     std::size_t count)
    {
      return norm_at_scale(x, count);
    }
 
    [[nodiscard]] inline double scaled_norm(const double* x,
                                            std::size_t   count)
    {
      return plain_or_at_scale(dot_product(x, x, count), x, count);
    }
 
    [[nodiscard]] inline double scaled_norm(const std::complex<double>* x,
                                            std::size_t                 count)
    {
      const auto* const parts {reinterpret_cast<const double*>(x)}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the layout guarantee above
      return plain_or_at_scale(dot_product(parts, parts, 2 * count), x, count);
    }
 
    template <typename T>
    [[nodiscard]] double scaled_norm(const buffer<T>& x)
    {
      return scaled_norm(x.data(), x.size());
    }
 
    [[nodiscard]] inline double rotation_length(double a,
                                                double b) noexcept
    {
      constexpr double low    {0x1p-500};
      constexpr double high   {0x1p500};
      const double     larger {std::max(std::abs(a), std::abs(b))};
      if (larger > low && larger < high) {
        return std::sqrt((a * a) + (b * b));
      }
      return std::hypot(a, b);
    }
  } // namespace detail
 
  class singular_matrix_error : public std::runtime_error
  {
  public:
 
    explicit singular_matrix_error(const std::string& message)
      : std::runtime_error(message)
    {}
  };
 
  class non_finite_error : public std::domain_error
  {
  public:
 
    explicit non_finite_error(const std::string& message)
      : std::domain_error(message)
    {}
  };
 
  namespace detail {
 
    template <typename T>
    [[nodiscard]] bool is_finite_entry(const T& x)
    {
      return std::isfinite(x);
    }
 
    template <typename U>
    [[nodiscard]] bool is_finite_entry(const std::complex<U>& x)
    {
      return std::isfinite(x.real()) && std::isfinite(x.imag());
    }
 
    template <typename T>
    void require_finite(const matrix<T>& a,
                        const char*      what)
    {
      for (std::size_t i {0}; i < a.rows(); ++i) {
        for (std::size_t j {0}; j < a.cols(); ++j) {
          if (!is_finite_entry(a(i, j))) {
            throw non_finite_error(std::string {what} + ": the matrix holds an infinite or NaN entry");
          }
        }
      }
    }
 
    [[nodiscard]] inline bool all_finite(const double* x,
                                         std::size_t   count) noexcept
    {
      constexpr std::uint64_t exponent {0x7FF0000000000000ULL};
      constexpr std::uint64_t carry    {std::uint64_t {1} << 52U};
      std::uint64_t           flagged  {0};
      std::size_t             i        {0};
#if defined(__GNUC__)
      using pair = std::uint64_t __attribute__((vector_size(16)));
      constexpr pair exponents {exponent, exponent};
      constexpr pair carries   {carry, carry};
      pair           first     {};
      pair           second    {};
      pair           third     {};
      pair           fourth    {};
      for (; i + 8 <= count; i += 8) {
        std::array<pair, 4> bits {};
        std::memcpy(bits.data(), x + i, sizeof(bits));
        first  |= (bits[0] & exponents) + carries;
        second |= (bits[1] & exponents) + carries;
        third  |= (bits[2] & exponents) + carries;
        fourth |= (bits[3] & exponents) + carries;
      }
      const pair all {first | second | third | fourth};
      flagged = all[0] | all[1];
#endif
      for (; i < count; ++i) {
        std::uint64_t bits {};
        std::memcpy(&bits, x + i, sizeof(bits));
        flagged |= (bits & exponent) + carry;
      }
      return (flagged >> 63U) == 0;
    }
 
    template <typename T>
    [[nodiscard]] matrix<T> parallel_copy_of(const matrix<T>& a)
    {
      buffer<T> elements(a.rows() * a.cols());
      parallel_copy(std::to_address(a.begin()), elements.size(), elements.data());
      return matrix<T> {a.rows(), a.cols(), std::move(elements)};
    }
 
    inline constexpr std::size_t finite_chunk {std::size_t {1} << 15U};
 
    [[nodiscard]] inline bool all_finite_shared(const double* x,
                                                std::size_t   count)
    {
      if (count < 2 * finite_chunk) {
        return all_finite(x, count);
      }
      std::atomic<bool> finite {true};
      parallel_for((count + finite_chunk - 1) / finite_chunk, memory_element_weight * finite_chunk / 2, [x, count, &finite](std::size_t chunk) {
                     const std::size_t first {chunk * finite_chunk};
                     if (!all_finite(x + first, std::min(finite_chunk, count - first))) {
                       finite.store(false, std::memory_order_relaxed);
                     }
                   });
      return finite.load(std::memory_order_relaxed);
    }
 
    inline void require_finite(const matrix<double>& a,
                               const char*           what)
    {
      if (!all_finite_shared(std::to_address(a.begin()), a.rows() * a.cols())) {
        throw non_finite_error(std::string {what} + ": the matrix holds an infinite or NaN entry");
      }
    }
 
    inline void require_finite(const matrix<std::complex<double>>& a,
                               const char*                         what)
    {
      const auto* const parts {reinterpret_cast<const double*>(std::to_address(a.begin()))}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the layout guarantee above
      if (!all_finite(parts, 2 * a.rows() * a.cols())) {
        throw non_finite_error(std::string {what} + ": the matrix holds an infinite or NaN entry");
      }
    }
 
    template <typename T>
    NUMERICAL_OUT_OF_LINE bool lower_is_finite(const matrix<T>& a)
    {
      for (std::size_t i {0}; i < a.rows(); ++i) {
        bool finite {std::isfinite(real_part(a(i, i)))};
        for (std::size_t j {0}; j < i; ++j) {
          finite = finite && is_finite_entry(a(i, j));
        }
        if (!finite) {
          return false;
        }
      }
      return true;
    }
 
    template <typename T>
    void require_finite_lower(const matrix<T>& a,
                              const char*      what)
    {
      if (!lower_is_finite(a)) {
        throw non_finite_error(std::string {what} + ": the matrix holds an infinite or NaN entry");
      }
    }
 
    NUMERICAL_OUT_OF_LINE inline void require_finite_lower(const matrix<double>& a,
                                                           const char*           what)
    {
      if (!all_finite_shared(std::to_address(a.begin()), a.rows() * a.cols()) && !lower_is_finite(a)) {
        throw non_finite_error(std::string {what} + ": the matrix holds an infinite or NaN entry");
      }
    }
 
    NUMERICAL_OUT_OF_LINE inline void require_finite_lower(const matrix<std::complex<double>>& a,
                                                           const char*                         what)
    {
      const auto* const parts {reinterpret_cast<const double*>(std::to_address(a.begin()))}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the layout guarantee above
      if (!all_finite(parts, 2 * a.rows() * a.cols()) && !lower_is_finite(a)) {
        throw non_finite_error(std::string {what} + ": the matrix holds an infinite or NaN entry");
      }
    }
 
    template <typename T>
    void subtract_panel_columns(T*          column,
                                const T*    columns,
                                std::size_t height,
                                const T*    l,
                                std::size_t terms)
    {
      for (std::size_t t {0}; t < terms; ++t) {
        subtract_product(column, columns + (t * height), conj_of(l[t]), height);
      }
    }
 
    inline void subtract_panel_columns(double*       column,
                                       const double* columns,
                                       std::size_t   height,
                                       const double* l,
                                       std::size_t   terms)
    {
      subtract_products(column, l, columns, height, terms, height);
    }
  } // namespace detail
 
  template <typename T>
  class lu_decomposition
  {
  public:
 
    explicit lu_decomposition(const matrix<T>& a)
      : lu_decomposition(detail::parallel_copy_of(a))
    {}
 
    explicit lu_decomposition(matrix<T>&& a)
      : n_(a.rows())
    {
      if (a.cols() != n_) {
        throw dimension_error("LU decomposition requires a square matrix");
      }
      detail::require_finite(a, "LU decomposition");
      lu_ = std::move(a).elements();
      pivots_.resize(n_);
      for (std::size_t i {0}; i < n_; ++i) {
        pivots_[i] = i;
      }
      factorize();
    }
 
    [[nodiscard]] vector<T> solve(const vector<T>& b) const
    {
      if (b.size() != n_) {
        throw dimension_error("solve: right-hand side size differs from the matrix");
      }
      if (singular_) {
        throw singular_matrix_error("solve: matrix is singular");
      }
      buffer<T> y(n_, T {});
      for (std::size_t i {0}; i < n_; ++i) { // forward substitution (unit lower L)
        T sum {b[pivots_[i]]};
        for (std::size_t j {0}; j < i; ++j) {
          detail::subtract_term(sum, lu_[(i * n_) + j], y[j]);
        }
        y[i] = sum;
      }
      for (std::size_t i {n_}; i-- > 0;) { // back substitution (upper U)
        T sum {y[i]};
        for (std::size_t j {i + 1}; j < n_; ++j) {
          detail::subtract_term(sum, lu_[(i * n_) + j], y[j]);
        }
        y[i] = detail::over_pivot(sum, detail::pivot_divisor(lu_[(i * n_) + i]));
      }
      return vector<T> {std::move(y)};
    }
 
    [[nodiscard]] T determinant() const
    {
      T value {detail::from_double<T>(sign_)};
      for (std::size_t i {0}; i < n_; ++i) {
        value *= lu_[(i * n_) + i];
      }
      return value;
    }
 
    [[nodiscard]] std::pair<T, double> log_determinant() const
    {
      if (singular_) {
        return {T {0}, -std::numeric_limits<double>::infinity()};
      }
      T      sign      {static_cast<T>(sign_)};
      double logarithm {0.0};
      for (std::size_t i {0}; i < n_; ++i) {
        const T pivot {lu_[(i * n_) + i]};
        sign      *= detail::phase_of(pivot);
        logarithm += std::log(static_cast<double>(std::abs(pivot)));
      }
      return {sign, logarithm};
    }
 
    [[nodiscard]] matrix<T> inverse() const
    {
      return inverse(buffer<T> {});
    }
 
    [[nodiscard]] matrix<T> inverse(buffer<T> storage) const
    {
      buffer<T>                rows {pivoted_inverse(std::move(storage)).elements()};
      const std::size_t        n    {n_};
      T* const                 at   {rows.data()};
      const std::size_t* const to   {pivots_.data()};
      // Bounds captured by value: by reference, the stores might have changed them, and each would be read again.
      detail::parallel_for((n + detail::inverse_band - 1) / detail::inverse_band, detail::memory_element_weight * detail::inverse_band * n, [at, n, to](std::size_t band) {
                             buffer<T> held(n); // column to[i] of A⁻¹ is column i of U⁻¹ L⁻¹
                             for (std::size_t r {band * detail::inverse_band}; r < std::min(n, (band + 1) * detail::inverse_band); ++r) {
                               std::copy_n(at + (r * n), n, held.data());
                               for (std::size_t i {0}; i < n; ++i) {
                                 at[(r * n) + to[i]] = held[i];
                               }
                             }
                           });
      return matrix<T> {n, n, std::move(rows)};
    }
 
    [[nodiscard]] matrix<T> pivoted_inverse(buffer<T> storage) const
    {
      if (singular_) {
        throw singular_matrix_error("inverse: matrix is singular");
      }
      const std::size_t n    {n_};
      storage.resize(n * n);
      T* const          rows {storage.data()};
      detail::parallel_for((n + detail::inverse_band - 1) / detail::inverse_band, detail::memory_element_weight * detail::inverse_band * n, [rows, n](std::size_t band) {
                             for (std::size_t r {band * detail::inverse_band}; r < std::min(n, (band + 1) * detail::inverse_band); ++r) {
                               std::fill_n(rows + (r * n), n, T {});
                               rows[(r * n) + r] = T {1};
                             }
                           });
      solve_rows(storage, n, true);
      return matrix<T> {n, n, std::move(storage)};
    }
 
    [[nodiscard]] matrix<T> solve(const matrix<T>& b) const
    {
      if (b.rows() != n_) {
        throw dimension_error("solve: right-hand side row count differs from the matrix");
      }
      if (singular_) {
        throw singular_matrix_error("solve: matrix is singular");
      }
      const std::size_t cols              {b.cols()};
      buffer<T>         rows(n_ * cols, T {});
      for (std::size_t i {0}; i < n_; ++i) {
        for (std::size_t c {0}; c < cols; ++c) {
          rows[(i * cols) + c] = b(pivots_[i], c);
        }
      }
      solve_rows(rows, cols);
      return matrix<T> {n_, cols, std::move(rows)};
    }
 
    [[nodiscard]] matrix<T> solve(matrix<T>&& b) const
    {
      if (b.rows() != n_) {
        throw dimension_error("solve: right-hand side row count differs from the matrix");
      }
      if (singular_) {
        throw singular_matrix_error("solve: matrix is singular");
      }
      const std::size_t cols {b.cols()};
      buffer<T>         rows {std::move(b).elements()};
      // Row i of P B is row pivots_[i] of B: each cycle i → pivots_[i] → … walked once, its first row held aside.
      std::vector<bool> placed(n_, false);
      buffer<T>         held(cols);
      for (std::size_t start {0}; start < n_; ++start) {
        if (placed[start] || pivots_[start] == start) {
          continue;
        }
        std::copy_n(rows.data() + (start * cols), cols, held.data());
        std::size_t i {start};
        for (std::size_t from {pivots_[i]}; from != start; i = from, from = pivots_[i]) {
          std::copy_n(rows.data() + (from * cols), cols, rows.data() + (i * cols));
          placed[i] = true;
        }
        std::copy_n(held.data(), cols, rows.data() + (i * cols));
        placed[i] = true;
      }
      solve_rows(rows, cols);
      return matrix<T> {n_, cols, std::move(rows)};
    }
 
    [[nodiscard]] matrix<T> factors() const&
    {
      return matrix<T> {n_, n_, lu_};
    }
 
    [[nodiscard]] matrix<T> factors() &&
    {
      return matrix<T> {n_, n_, std::move(lu_)};
    }
 
    [[nodiscard]] const std::vector<std::size_t>& pivots() const noexcept
    {
      return pivots_;
    }
 
    [[nodiscard]] bool singular() const noexcept
    {
      return singular_;
    }
 
  private:
 
 
    void solve_rows(buffer<T>&      y,
                    std::size_t     cols,
                    bool            lower = false) const
    {
      if constexpr (std::is_same_v<T, double>) { // there subtract_products holds its row in registers through the terms
        // The rows beat the blocked passes on the x86-64 units at any order (the portable path's products there fuse
        // nothing, the rows would), and at any order elsewhere, by groups (arm64, four columns: order 256 0.65 times the
        // blocked passes' time under Clang, 0.45 under GCC; order 1024 on every thread 0.91 and 0.62).
#ifdef NUMERICAL_X86_COMPILED
        const bool narrow {!lower && !y.empty() && cols <= detail::narrow_solve && detail::gemm_path_current() != gemm_path::portable};
#else
        const bool narrow {!lower && !y.empty() && cols <= detail::narrow_solve};
#endif
        if (narrow) {
          solve_narrow(y, cols);
          return;
        }
      }
      const std::size_t slices {column_slices(cols)};
      if (slices > 1) { // independent right-hand sides: each thread its own columns, all the blocks, no join per block
        const std::size_t slice {(((cols + slices - 1) / slices) + 7) / 8 * 8};
        const auto        start {[this, slice, slices, cols, lower](std::size_t s) {
                                   return lower ? lower_slice_start(s, slices, cols) : std::min(cols, s * slice);
                                 }};
        detail::parallel_for(lower ? slices : (cols + slice - 1) / slice, n_ * n_ * slice / detail::gemm_kernel_speedup, [&](std::size_t s) {
                               solve_columns(y, cols, start(s), start(s + 1), lower);
                             });
        return;
      }
      solve_columns(y, cols, 0, cols, lower);
    }
 
    [[nodiscard]] std::size_t column_slices(std::size_t cols) const
    {
#ifdef NUMERICAL_X86_COMPILED
      if (n_ > detail::slice_solve_order) {
        return 1;
      }
#endif
      if (cols < 2 * detail::slice_columns) {
        return 1;
      }
      return std::min(max_threads(), cols / detail::slice_columns);
    }
 
    [[nodiscard]] std::size_t lower_slice_start(std::size_t s,
                                                std::size_t slices,
                                                std::size_t cols) const
    {
      const auto   n {static_cast<double>(n_)};
      const auto   weight {[n](double c) { return (n * n * c) + (((n * n * n) - ((n - c) * (n - c) * (n - c))) / 3.0); }};
      const double target {weight(static_cast<double>(cols)) * static_cast<double>(s) / static_cast<double>(slices)};
      std::size_t  low    {0};
      std::size_t  high   {cols};
      while (low < high) { // the fewest columns that weigh the target
        const std::size_t middle {low + ((high - low) / 2)};
        if (weight(static_cast<double>(middle)) < target) {
          low = middle + 1;
        }
        else {
          high = middle;
        }
      }
      return std::min(cols, (low + 7) / 8 * 8);
    }
 
    void solve_columns(buffer<T>&      y,
                       std::size_t     cols,
                       std::size_t     c0,
                       std::size_t     c1,
                       bool            lower) const
    {
      constexpr std::size_t width  {64};
      const std::size_t     chunks {(c1 - c0 + width - 1) / width};
      for (std::size_t k0 {0}; k0 < n_; k0 += block_) {           // forward: unit-lower L
        const std::size_t end   {std::min(k0 + block_, n_)};
        const std::size_t reach {lower ? std::min(end, c1) : c1}; // the block's rows' columns that can be nonzero
        if (reach <= c0) {
          continue; // the slice's columns are zero in these rows and every row past them
        }
        detail::parallel_for((reach - c0 + width - 1) / width, (end - k0) * (end - k0) * width / 2, [&](std::size_t chunk) {
                               const std::size_t first {c0 + (chunk * width)};
                               const std::size_t last {std::min(first + width, reach)};
                               for (std::size_t i {k0 + 1}; i < end; ++i) {
                                 detail::subtract_products(y.data() + (i * cols) + first, lu_.data() + (i * n_) + k0, y.data() + (k0 * cols) + first, cols, i - k0, last - first);
                               }
                             });
        subtract_block(y, cols, end, n_, k0, end, c0, reach);
      }
      for (std::size_t k1 {n_}; k1 > 0;) { // backward: upper U, the blocks of the forward pass, last first
        const std::size_t k0 {panel_of(k1 - 1)};
        detail::parallel_for(chunks, (k1 - k0) * (k1 - k0) * width / 2, [&](std::size_t chunk) {
                               const std::size_t first {c0 + (chunk * width)};
                               const std::size_t last {std::min(first + width, c1)};
                               for (std::size_t i {k1}; i-- > k0;) {
                                 detail::subtract_products(y.data() + (i * cols) + first, lu_.data() + (i * n_) + i + 1, y.data() + ((i + 1) * cols) + first, cols, k1 - (i + 1), last - first);
                                 const T divisor {detail::pivot_divisor(lu_[(i * n_) + i])};
                                 for (std::size_t c {first}; c < last; ++c) {
                                   y[(i * cols) + c] = detail::over_pivot(y[(i * cols) + c], divisor);
                                 }
                               }
                             });
        subtract_block(y, cols, 0, k0, k0, k1, c0, c1);
        k1 = k0;
      }
    }
 
    void solve_narrow(buffer<T>&      y,
                      std::size_t     cols) const
    {
      // Forward, by groups of rows: their terms before the group are independent of one another, and go together
      // (\ref detail::subtract_rows), the group's own triangle after; each element's terms in ascending order still.
      constexpr std::size_t group {detail::product_rows};
      static_assert(detail::lu_small_block % group == 0 && detail::lu_block % group == 0, "a panel is a whole number of groups");
      std::size_t           i0    {0};
      for (; i0 + group <= n_; i0 += group) {
        if (i0 > 0) {
          detail::subtract_rows(y.data() + (i0 * cols), cols, lu_.data() + (i0 * n_), n_, y.data(), cols, i0, cols);
        }
        for (std::size_t i {i0 + 1}; i < i0 + group; ++i) {
          detail::subtract_products(y.data() + (i * cols), lu_.data() + (i * n_) + i0, y.data() + (i0 * cols), cols, i - i0, cols);
        }
      }
      for (std::size_t i {std::max<std::size_t>(i0, 1)}; i < n_; ++i) {
        detail::subtract_products(y.data() + (i * cols), lu_.data() + (i * n_), y.data(), cols, i, cols);
      }
      // Backward, each row's terms block by block from the last down to its own, ascending within each: a group's
      // terms in the blocks past its own are independent of one another and go together (\ref detail::subtract_rows),
      // in that same order; its own block's after, row by row (\p finish). The rows past the last whole group lie in
      // the last block, a panel being a whole number of groups: they have their own block's terms alone.
      const std::size_t last {panel_of(n_ - 1)};
      const auto        finish {[this, &y, cols](std::size_t i) {
                                  const std::size_t to {std::min(panel_of(i) + block_, n_)};
                                  if (i + 1 < to) {
                                    detail::subtract_products(y.data() + (i * cols), lu_.data() + (i * n_) + i + 1, y.data() + ((i + 1) * cols), cols, to - (i + 1), cols);
                                  }
                                  const T divisor {detail::pivot_divisor(lu_[(i * n_) + i])};
                                  for (std::size_t c {0}; c < cols; ++c) {
                                    y[(i * cols) + c] = detail::over_pivot(y[(i * cols) + c], divisor);
                                  }
                                }};
      std::size_t i1 {n_};
      for (; i1 % group != 0; --i1) {
        finish(i1 - 1);
      }
      for (; i1 > 0; i1 -= group) {
        const std::size_t top {i1 - group};
        const std::size_t own {panel_of(top)};
        for (std::size_t b0 {last}; b0 > own; b0 -= block_) {
          detail::subtract_rows(y.data() + (top * cols), cols, lu_.data() + (top * n_) + b0, n_, y.data() + (b0 * cols), cols, std::min(block_, n_ - b0), cols);
        }
        for (std::size_t i {i1}; i-- > top;) {
          finish(i);
        }
      }
    }
 
    [[nodiscard]] std::size_t panel_of(std::size_t i) const noexcept
    {
      return i & ~(block_ - 1);
    }
 
    void subtract_block(buffer<T>&      y,
                        std::size_t     cols,
                        std::size_t     r0,
                        std::size_t     r1,
                        std::size_t     c0,
                        std::size_t     c1,
                        std::size_t     first,
                        std::size_t     last) const
    {
      if (r0 < r1) {
        const buffer<T> minus {negated_block(r0, r1, c0, c1)};
        detail::gemm(r1 - r0, last - first, c1 - c0, minus.data(), c1 - c0, y.data() + (c0 * cols) + first, cols,
                     y.data() + (r0 * cols) + first, cols);
      }
    }
 
    void negated_block_into(buffer<T>&  out,
                            std::size_t r0,
                            std::size_t r1,
                            std::size_t c0,
                            std::size_t c1) const
    {
      out.resize((r1 - r0) * (c1 - c0));
      for (std::size_t r {r0}; r < r1; ++r) {
        for (std::size_t c {c0}; c < c1; ++c) {
          out[((r - r0) * (c1 - c0)) + (c - c0)] = -lu_[(r * n_) + c];
        }
      }
    }
 
    [[nodiscard]] buffer<T> negated_block(std::size_t r0,
                                          std::size_t r1,
                                          std::size_t c0,
                                          std::size_t c1) const
    {
      buffer<T> out((r1 - r0) * (c1 - c0), T {});
      for (std::size_t r {r0}; r < r1; ++r) {
        for (std::size_t c {c0}; c < c1; ++c) {
          out[((r - r0) * (c1 - c0)) + (c - c0)] = -lu_[(r * n_) + c];
        }
      }
      return out;
    }
 
    // Right-looking, by panels of block_ columns, as LAPACK's dgetrf: each panel is factorised
    // with partial pivoting (whole rows swapped), its rows of U to the right are solved against the panel's
    // unit-lower L, and the trailing block takes the product of the two through detail::gemm, where the
    // work of a large factorisation is. Each element still receives its updates in ascending order of the
    // pivot, so without multiply-add fusion the factors are the unblocked algorithm's to the bit; an order
    // up to one panel is that algorithm.
    //
    // With a look-ahead of one panel: the trailing update reaches the next panel's columns first, then the
    // next panel is factorised on the calling thread while the workers update the columns right of it (the
    // panel, a column at a time, would otherwise hold every other thread idle). Its row swaps reach those
    // columns once both are done: an update reads its rows of L from a copy and writes row by row, so a swap
    // after it is the swap before it, and the factors stay the ones above, to the bit.
    //
    // A pivot's swap reaches the columns left of its panel only at the end (\ref swap_left_columns): nothing reads
    // them before, and on the calling thread each swap held up the panel the other threads wait for.
    void factorize()
    {
      // Every panel in the one column-major storage, of the first (and tallest) one's size: a fresh one each panel was
      // mapped anew and faulted in page by page, on the critical path (order 1000 on six AVX2 cores: 24 µs of a panel's
      // 470, and its copy in 97 µs instead of about 40).
      buffer<T>         panel(n_ * std::min(block_, n_));
      const std::size_t end {std::min(block_, n_)};
      factorize_panel(panel, 0, end);
      if (end < n_) {
        factorize_trailing(panel, end);
      }
      swap_left_columns();
    }
 
    NUMERICAL_OUT_OF_LINE void factorize_trailing(buffer<T>&  panel,
                                                  std::size_t first_end)
    {
      // −L₂₁ and its packing for the product's tile, in storage kept from panel to panel: packed once, it serves every
      // block of the trailing update, where each block's product packed it again (6.9 % of the time on six AVX2 cores
      // at order 1000).
      buffer<T>   minus_l;
      buffer<T>   packing;
      std::size_t end {first_end};
      for (std::size_t k0 {0}; end < n_;) {
        const std::size_t next       {std::min(end + block_, n_)};
        solve_panel_rows(k0, end);
        negated_block_into(minus_l, end, n_, k0, end);
        const detail::packed_left<T> left {detail::pack_left(n_ - end, block_, end - k0, minus_l.data(), end - k0, packing)};
        update_columns(left, k0, end, next);
        const std::size_t right           {n_ - next};
        if (right > detail::lu_lookahead_columns) { // one product of the whole trailing block, then the panel
          update_columns(left, k0, next, n_);
          factorize_panel(panel, end, next);
          k0  = end;
          end = next;
          continue;
        }
        // About two blocks of columns for each thread the work pays for: a product too small to share ran one block for
        // each thread the machine has (order 256 on eight arm64 threads: 329 → 311 µs); two, not one, as on x86-64 a
        // block of half the columns stays in its cache, and the balance beat one (six AVX2 cores, order 1000: 12 %).
        const std::size_t work    {(n_ - end) * right * (end - k0) / detail::lu_update_speedup};
        const std::size_t helpers {std::clamp(work / detail::work_per_thread, std::size_t {1}, max_threads())};
        const std::size_t width   {std::max(block_, (right + (2 * helpers) - 1) / (2 * helpers))};
        const std::size_t blocks  {(right + width - 1) / width};
        swap_limit_ = next;
        detail::parallel_for_beside(blocks, (n_ - end) * width * (end - k0) / detail::lu_update_speedup,
                                    [&](std::size_t block) {
                                      const std::size_t first {next + (block * width)};
                                      update_columns(left, k0, first, std::min(first + width, n_));
                                    },
                                    [&] { factorize_panel(panel, end, next); });
        swap_limit_ = n_;
        for (const auto& [a, b] : deferred_swaps_) {
          for (std::size_t c {next}; c < n_; ++c) {
            std::swap(lu_[(a * n_) + c], lu_[(b * n_) + c]);
          }
        }
        deferred_swaps_.clear();
        k0  = end;
        end = next;
      }
    }
 
    void factorize_panel(buffer<T>&      panel,
                         std::size_t     k0,
                         std::size_t     end)
    {
      const std::size_t rows                  {n_ - k0};
      const std::size_t width                 {end - k0};
      transpose_panel(panel, rows, k0, end, true); // column j of the panel at panel[j * rows], every element written
      factorize_columns(panel, rows, 0, width, k0, end);
      transpose_panel(panel, rows, k0, end, false);
    }
 
    void transpose_panel(buffer<T>&      panel,
                         std::size_t     rows,
                         std::size_t     k0,
                         std::size_t     end,
                         bool            in)
    {
      const std::size_t width {end - k0};
#ifdef NUMERICAL_X86_COMPILED
      if (rows >= detail::panel_tile_rows) {
        T* const base {lu_.data() + (k0 * n_) + k0};
        for (std::size_t r0 {0}; r0 < rows; r0 += detail::panel_copy_rows) {
          const std::size_t r1 {std::min(r0 + detail::panel_copy_rows, rows)};
          if (in) {
            for (std::size_t j {0}; j < width; ++j) {
              for (std::size_t r {r0}; r < r1; ++r) {
                panel[(j * rows) + r] = base[(r * n_) + j];
              }
            }
          }
          else {
            for (std::size_t j {0}; j < width; ++j) {
              for (std::size_t r {r0}; r < r1; ++r) {
                base[(r * n_) + j] = panel[(j * rows) + r];
              }
            }
          }
        }
        return;
      }
#endif
      for (std::size_t r {0}; r < rows; ++r) {
        for (std::size_t j {0}; j < width; ++j) {
          if (in) {
            panel[(j * rows) + r] = lu_[((k0 + r) * n_) + k0 + j];
          }
          else {
            lu_[((k0 + r) * n_) + k0 + j] = panel[(j * rows) + r];
          }
        }
      }
    }
 
    void factorize_columns(buffer<T>&      panel,
                           std::size_t     rows,
                           std::size_t     first,
                           std::size_t     last,
                           std::size_t     k0,
                           std::size_t     end)
    {
      // A leaf, or a panel too short for its products to pay for packing (order 64 on arm64: 30 µs unblocked,
      // 37 recursive).
      if (last - first <= detail::solve_leaf || rows - first < 2 * block_) {
        for (std::size_t j {first}; j < last; ++j) {
          T           *          column {panel.data() + (j * rows)};
          std::size_t            pivot  {0}; // from row j
          detail::magnitude_t<T> best   {std::abs(column[j])};
          for (std::size_t r {detail::largest_magnitude_vectors(column + j, rows - j, best, pivot)}; r < rows - j; ++r) {
            const detail::magnitude_t<T> mag {std::abs(column[j + r])};
            if (mag > best) {
              best  = mag;
              pivot = r;
            }
          }
          pivot += j;
          if (best == detail::magnitude_t<T> {}) {
            singular_ = true;
            continue; // a zero pivot column: leave it, the determinant becomes 0
          }
          if (pivot != j) {
            for (std::size_t c {0}; c < end - k0; ++c) { // the whole panel row
              std::swap(panel[(c * rows) + j], panel[(c * rows) + pivot]);
            }
            swap_rows_outside(k0 + j, k0 + pivot, k0, end);
            std::swap(pivots_[k0 + j], pivots_[k0 + pivot]);
            sign_ = -sign_;
          }
          const T divisor {detail::pivot_divisor(column[j])};
          for (std::size_t r {j + 1}; r < rows; ++r) {
            column[r] = detail::over_pivot(column[r], divisor); // the multiplier, in L's slot
          }
          for (std::size_t c {j + 1}; c < last; ++c) {
            T*      target {panel.data() + (c * rows)};
            detail::subtract_product(target + j + 1, column + j + 1, target[j], rows - (j + 1));
          }
        }
        return;
      }
      const std::size_t mid {first + ((last - first) / 2)};
      factorize_columns(panel, rows, first, mid, k0, end);
      for (std::size_t c {mid}; c < last; ++c) { // the top of the right half: rows [first, mid)
        T* target {panel.data() + (c * rows)};
        for (std::size_t j {first}; j < mid; ++j) {
          detail::subtract_product(target + j + 1, panel.data() + (j * rows) + j + 1, target[j], mid - (j + 1));
        }
      }
      // The bottom, rows [mid, rows): A₂₂ −= L₂₁ U₁₂, column-major, so as its transpose, row-major:
      // A₂₂ᵀ += (−U₁₂ᵀ) L₂₁ᵀ.
      const std::size_t depth                           {mid - first};
      buffer<T>         minus_u((last - mid) * depth, T {});
      for (std::size_t c {mid}; c < last; ++c) {
        for (std::size_t j {first}; j < mid; ++j) {
          minus_u[((c - mid) * depth) + (j - first)] = -panel[(c * rows) + j];
        }
      }
      // The plain loop, which rounds as the tiles do (the same bits): a product this narrow packed a sliver of the
      // panel's whole height for each few columns (six x86-64 AVX2 cores, LU of order 1000: 6.86 → 6.36 ms).
      detail::gemm_rows(last - mid, rows - mid, depth, minus_u.data(), depth, panel.data() + (first * rows) + mid, rows,
                        panel.data() + (mid * rows) + mid, rows);
      factorize_columns(panel, rows, mid, last, k0, end);
    }
 
    void swap_rows_outside(std::size_t a,
                           std::size_t b,
                           std::size_t k0,
                           std::size_t end)
    {
      if (k0 > 0) {
        left_swaps_.emplace_back(a, b);
      }
      for (std::size_t c {end}; c < swap_limit_; ++c) {
        std::swap(lu_[(a * n_) + c], lu_[(b * n_) + c]);
      }
      if (swap_limit_ < n_) {
        deferred_swaps_.emplace_back(a, b);
      }
    }
 
    void swap_left_columns()
    {
      if (left_swaps_.empty()) {
        return;
      }
      const std::size_t panels {(n_ - 1) / block_}; // every panel but the last, which no later swap reaches
      detail::parallel_for(panels, left_swaps_.size() * block_, [&](std::size_t p) {
                             const std::size_t c0 {p * block_};
                             const std::size_t c1 {c0 + block_};
                             const auto first {std::lower_bound(left_swaps_.begin(), left_swaps_.end(), c1,
                                                                [](const auto& swap, std::size_t row) { return swap.first < row; })};
                             for (auto swap {first}; swap != left_swaps_.end(); ++swap) {
                               std::swap_ranges(lu_.data() + (swap->first * n_) + c0, lu_.data() + (swap->first * n_) + c1,
                                                lu_.data() + (swap->second * n_) + c0);
                             }
                           });
      left_swaps_.clear();
    }
 
    void solve_panel_rows(std::size_t k0,
                          std::size_t end)
    {
      constexpr std::size_t width  {64};
      const std::size_t     blocks {((n_ - end) + width - 1) / width};
      detail::parallel_for(blocks, (end - k0) * (end - k0) * width / 2, [&](std::size_t block) {
                             const std::size_t first {end + (block * width)};
                             const std::size_t last {std::min(first + width, n_)};
                             for (std::size_t i {k0 + 1}; i < end; ++i) {
                               detail::subtract_products(lu_.data() + (i * n_) + first, lu_.data() + (i * n_) + k0, lu_.data() + (k0 * n_) + first, n_, i - k0, last - first);
                             }
                           });
    }
 
    void update_columns(const detail::packed_left<T>& left,
                        std::size_t                   k0,
                        std::size_t                   c0,
                        std::size_t                   c1)
    {
      const std::size_t end {k0 + left.k};
      detail::gemm_with(left, c1 - c0, lu_.data() + (k0 * n_) + c0, n_, lu_.data() + (end * n_) + c0, n_);
    }
 
    std::size_t                                      n_;                                
    std::size_t                                      block_ {detail::lu_block_for(n_)}; 
    buffer<T>                                        lu_;                               
    std::vector<std::size_t>                         pivots_;                           
    int                                              sign_       {1};                   
    bool                                             singular_   {false};               
    std::size_t                                      swap_limit_ {n_};                  
    std::vector<std::pair<std::size_t, std::size_t>> deferred_swaps_;                   
    std::vector<std::pair<std::size_t, std::size_t>> left_swaps_;                       
  };
 
  template <typename T>
  [[nodiscard]] vector<T> solve(matrix<T>        a,
                                const vector<T>& b)
  {
    return lu_decomposition<T> {std::move(a)}.solve(b);
  }
 
  template <typename T>
  [[nodiscard]] matrix<T> solve(matrix<T>        a,
                                const matrix<T>& b)
  {
    return lu_decomposition<T> {std::move(a)}.solve(b);
  }
 
  template <typename T>
  [[nodiscard]] T determinant(matrix<T> a)
  {
    return lu_decomposition<T> {std::move(a)}.determinant();
  }
 
  template <typename T>
  [[nodiscard]] matrix<T> inverse(matrix<T> a)
  {
    return lu_decomposition<T> {std::move(a)}.inverse();
  }
 
  class not_positive_definite_error : public std::runtime_error
  {
  public:
 
    explicit not_positive_definite_error(const std::string& message,
                                         std::size_t        order = 0)
      : std::runtime_error(message), order_ {order}
    {}
 
    [[nodiscard]] std::size_t order() const noexcept
    {
      return order_;
    }
 
  private:
 
    std::size_t order_; 
  };
 
  template <typename T>
  class cholesky_decomposition
  {
  public:
 
    explicit cholesky_decomposition(const matrix<T>& a)
      : cholesky_decomposition(detail::parallel_copy_of(a))
    {}
 
    explicit cholesky_decomposition(matrix<T>&& a)
      : n_(a.rows())
    {
      if (a.cols() != n_) {
        throw dimension_error("Cholesky decomposition requires a square matrix");
      }
      detail::require_finite_lower(a, "Cholesky decomposition");
      lower_ = std::move(a).elements(); // the upper triangle is never read, and zeroed once factorised
      factorize();
    }
 
    [[nodiscard]] matrix<T> lower() const&
    {
      return matrix<T> {n_, n_, lower_};
    }
 
    [[nodiscard]] matrix<T> lower() &&
    {
      return matrix<T> {n_, n_, std::move(lower_)};
    }
 
    [[nodiscard]] vector<T> solve(const vector<T>& b) const
    {
      if (b.size() != n_) {
        throw dimension_error("solve: right-hand side size differs from the matrix");
      }
      buffer<T> y(n_, T {});
      for (std::size_t i {0}; i < n_; ++i) { // forward: L y = b
        T sum {b[i]};
        for (std::size_t k {0}; k < i; ++k) {
          sum -= lower_[(i * n_) + k] * y[k];
        }
        y[i] = sum / lower_[(i * n_) + i];
      }
      buffer<T> x(n_, T {});
      for (std::size_t i {n_}; i-- > 0;) { // back: Lᴴ x = y, (Lᴴ)[i][k] = conj(L[k][i])
        T sum {y[i]};
        for (std::size_t k {i + 1}; k < n_; ++k) {
          sum -= detail::conj_of(lower_[(k * n_) + i]) * x[k];
        }
        x[i] = sum / detail::conj_of(lower_[(i * n_) + i]);
      }
      return vector<T> {std::move(x)};
    }
 
    [[nodiscard]] T determinant() const
    {
      T value {static_cast<T>(1)};
      for (std::size_t i {0}; i < n_; ++i) {
        value *= lower_[(i * n_) + i] * lower_[(i * n_) + i];
      }
      return value;
    }
 
  private:
 
    // By panels of block_ columns: the panel's diagonal block row by row, the rows below it
    // by leaves (each depends only on the diagonal block; factor_rows_below), then the lower part of the trailing block takes
    // A22 += (-L21) L21ᴴ from the matrix product kernel. Every element still receives its terms in ascending
    // order, so without multiply-add fusion the factor is the unblocked algorithm's to the bit.
    void factorize()
    {
      buffer<T> packing;        // L₂₁ packed for the trailing update's tiles, kept from panel to panel
      buffer<T> minus_adjoint;  // −L₂₁ᴴ, row-major, the same
      for (std::size_t k0 {0}; k0 < n_; k0 += block_) {
        const std::size_t end {std::min(k0 + block_, n_)};
        std::size_t       i   {k0};
        for (; end - k0 >= detail::cholesky_together_width && i + detail::cholesky_diagonal_rows <= end; i += detail::cholesky_diagonal_rows) {
          factor_rows_together(i, k0);
          for (std::size_t r {0}; r < detail::cholesky_diagonal_rows; ++r) { // their triangle, row after row
            factor_row(i + r, k0, i, i + r + 1);
          }
        }
        for (; i < end; ++i) {
          factor_row(i, k0, k0, i + 1);
        }
        if (end < n_) {
          minus_adjoint.resize((end - k0) * (n_ - end));
          factor_rows_below(k0, end, minus_adjoint);
          update_trailing(k0, end, minus_adjoint, packing);
        }
      }
      T* const          rows {lower_.data()};
      const std::size_t n    {n_};
      // The upper triangle zeroed by bands of rows on the threads: a pass bound by the memory, half a row an element.
      detail::parallel_for((n + detail::zeroed_rows - 1) / detail::zeroed_rows, detail::memory_element_weight * detail::zeroed_rows * n / 2, [rows, n](std::size_t band) {
                             for (std::size_t i {band * detail::zeroed_rows}; i < std::min(n, (band + 1) * detail::zeroed_rows); ++i) {
                               std::fill(rows + (i * n) + i + 1, rows + ((i + 1) * n), T {});
                             }
                           });
    }
 
    void factor_rows_together(std::size_t i,
                              std::size_t k0)
    {
      std::array<T, detail::cholesky_diagonal_rows> sums {};
      for (std::size_t j {k0}; j < i; ++j) {
        const T* const above {lower_.data() + (j * n_)};
        NUMERICAL_UNROLL_BLOCK
        for (std::size_t r {0}; r < detail::cholesky_diagonal_rows; ++r) {
          sums[r] = lower_[((i + r) * n_) + j];
        }
        for (std::size_t k {k0}; k < j; ++k) {
          const T other {detail::conj_of(above[k])};
          NUMERICAL_UNROLL_BLOCK
          for (std::size_t r {0}; r < detail::cholesky_diagonal_rows; ++r) {
            sums[r] -= lower_[((i + r) * n_) + k] * other;
          }
        }
        NUMERICAL_UNROLL_BLOCK
        for (std::size_t r {0}; r < detail::cholesky_diagonal_rows; ++r) {
          lower_[((i + r) * n_) + j] = sums[r] / above[j];
        }
      }
    }
 
    void factor_row(std::size_t i,
                    std::size_t k0,
                    std::size_t first,
                    std::size_t last)
    {
      for (std::size_t j {first}; j < last; ++j) {
        T sum {lower_[(i * n_) + j]};
        for (std::size_t k {k0}; k < j; ++k) {
          sum -= lower_[(i * n_) + k] * detail::conj_of(lower_[(j * n_) + k]);
        }
        if (i == j) {
          const double diagonal {detail::real_part(sum)};
          if (diagonal <= 0.0) {
            throw not_positive_definite_error("Cholesky: matrix is not positive-definite", i + 1);
          }
          lower_[(i * n_) + i] = detail::from_double<T>(std::sqrt(diagonal));
        }
        else {
          lower_[(i * n_) + j] = sum / lower_[(j * n_) + j];
        }
      }
    }
 
    void factor_rows_below(std::size_t k0,
                           std::size_t end,
                           buffer<T>&  minus_adjoint)
    {
      const std::size_t width  {end - k0};
      const std::size_t chunks {((n_ - end) + detail::cholesky_rows - 1) / detail::cholesky_rows};
      detail::parallel_for(chunks, detail::cholesky_rows * width * width / 2, [&](std::size_t chunk) {
                             const std::size_t first {end + (chunk * detail::cholesky_rows)};
                             const std::size_t height {std::min(detail::cholesky_rows, n_ - first)};
                             buffer<T>    panel(height * width); // column j at panel[j * height], every element written
                             for (std::size_t r {0}; r < height; ++r) {
                               for (std::size_t j {0}; j < width; ++j) {
                                 panel[(j * height) + r] = lower_[((first + r) * n_) + k0 + j];
                               }
                             }
                             constexpr std::size_t sub {16};
                             buffer<T>             minus_conj;
                             for (std::size_t s0 {0}; s0 < width; s0 += sub) {
                               const std::size_t s1 {std::min(s0 + sub, width)};
                               if (s0 > 0) { // the earlier sub-blocks' terms, through the product kernel: Pᵀ[s0, s1) += (−conj L₁₁[s0, s1) × [0, s0)) Pᵀ[0, s0)
                                 minus_conj.resize((s1 - s0) * s0);
                                 for (std::size_t j {s0}; j < s1; ++j) {
                                   for (std::size_t k {0}; k < s0; ++k) {
                                     minus_conj[((j - s0) * s0) + k] = -detail::conj_of(lower_[((k0 + j) * n_) + k0 + k]);
                                   }
                                 }
                                 // The plain loop, which rounds as the tiles do (the same bits), where the packed
                                 // product packed slivers too small to pay for it (six AVX2 cores' one, order 1000:
                                 // this pass 2.82 → 2.23 ms).
                                 detail::gemm_rows(s1 - s0, height, s0, minus_conj.data(), s0, panel.data(), height, panel.data() + (s0 * height), height);
                               }
                               for (std::size_t j {s0}; j < s1; ++j) {
                                 T* const column {panel.data() + (j * height)};
                                 detail::subtract_panel_columns(column, panel.data() + (s0 * height), height, lower_.data() + ((k0 + j) * n_) + k0 + s0, j - s0);
                                 const T pivot {lower_[((k0 + j) * n_) + k0 + j]};
                                 for (std::size_t r {0}; r < height; ++r) {
                                   column[r] = column[r] / pivot;
                                 }
                               }
                             }
                             const std::size_t below {n_ - end};
#ifdef NUMERICAL_X86_COMPILED
                             // On the x86-64 units the factor's rows, then −L₂₁ᴴ's columns, each written along: written
                             // together, one of them went across (one AVX2 core, order 1000: this pass 2.23 → 1.97 ms).
                             // Together on arm64, where apart was 3 % slower under GCC 14 (orders 48 to 128).
                             for (std::size_t r {0}; r < height; ++r) {
                               for (std::size_t j {0}; j < width; ++j) {
                                 lower_[((first + r) * n_) + k0 + j] = panel[(j * height) + r];
                               }
                             }
                             for (std::size_t j {0}; j < width; ++j) {
                               const T* const from {panel.data() + (j * height)};
                               T* const       to {minus_adjoint.data() + (j * below) + (first - end)};
                               for (std::size_t r {0}; r < height; ++r) {
                                 to[r] = -detail::conj_of(from[r]);
                               }
                             }
#else
                             for (std::size_t r {0}; r < height; ++r) {
                               for (std::size_t j {0}; j < width; ++j) {
                                 const T value {panel[(j * height) + r]};
                                 lower_[((first + r) * n_) + k0 + j]            = value;
                                 minus_adjoint[(j * below) + (first - end + r)] = -detail::conj_of(value);
                               }
                             }
#endif
                           });
    }
 
    void update_trailing(std::size_t      k0,
                         std::size_t      end,
                         const buffer<T>& minus_adjoint,
                         buffer<T>&       packing)
    {
      const std::size_t                                rows    {n_ - end};
      const std::size_t                                depth   {end - k0};
#ifdef NUMERICAL_X86_COMPILED
      // On the x86-64 units a real product the kernel packs runs on the triangle's tiles alone (detail::gemm_triangle,
      // each element the tiles' sum below), where the blocks below computed the part of each diagonal block past the
      // diagonal: six AVX2 cores, order 500 787 → 721 µs. On arm64 the blocks stay, faster there (all cores, order 500:
      // 0.85 ms, by the triangle 1.03: the AMX takes −L₂₁ packed once for every block). A smaller or complex product by
      // the blocks.
      if constexpr (std::is_same_v<T, double>) {
        if (depth >= detail::gemm_min_depth && rows * rows >= ((detail::gemm_packed_threshold + depth) - 1) / depth) {
          detail::gemm_triangle<detail::a_rows, detail::b_rows>(rows, depth, lower_.data() + (end * n_) + k0, n_, minus_adjoint.data(), rows,
                                                                lower_.data() + (end * n_) + end, n_, true);
          return;
        }
      }
#endif
      const std::size_t                                step    {std::max(block_, detail::cholesky_update_columns)};
      // Packed once where a block of rows serves several tiles: past two tiles' rows (one AVX2 core, the shared packing
      // cost order 300 5 % and 500 3 %, where a block of columns has one or two tiles), else each tile packs its own.
      const bool                                       shared  {rows > 2 * detail::cholesky_update_rows};
      const T* const                                   l21     {lower_.data() + (end * n_) + k0};
      const detail::packed_left<T>                     left    {shared ? detail::pack_left(rows, step, depth, l21, n_, packing)
                                                                       : detail::read_in_place(rows, depth, l21, n_)};
      const std::size_t                                granule {std::max<std::size_t>(1, left.granule)}; // a row at least
      const std::size_t                                height  {std::max(granule, detail::cholesky_update_rows / granule * granule)};
      std::vector<std::pair<std::size_t, std::size_t>> tiles;                                            // each one's first column and first row, below the panel
      for (std::size_t jb {0}; jb < rows; jb += step) {
        for (std::size_t r0 {jb / granule * granule}; r0 < rows; r0 += height) {
          tiles.emplace_back(jb, r0);
        }
      }
      detail::parallel_for(tiles.size(), height * step * depth / detail::gemm_kernel_speedup, [&](std::size_t t) {
                             const auto [jb, r0] {tiles[t]};
                             detail::gemm_with(detail::left_rows(left, r0, std::min(height, rows - r0)), std::min(step, rows - jb),
                                               minus_adjoint.data() + jb, rows, lower_.data() + ((end + r0) * n_) + end + jb, n_);
                           });
    }
 
    std::size_t    n_;                                      
    std::size_t    block_ {detail::cholesky_block_for(n_)}; 
    buffer<T>      lower_;                                  
  };
 
  template <typename T>
  [[nodiscard]] cholesky_decomposition<T> cholesky(matrix<T> a)
  {
    return cholesky_decomposition<T> {std::move(a)};
  }
 
  namespace detail {
 
    template <typename T>
    struct thin_qr
    {
      buffer<T> q; 
      buffer<T> r; 
    };
 
    template <typename T>
    void reflect_rows(buffer<T>&            a,
                      std::size_t           cols,
                      const buffer<T>&      v,
                      std::size_t           first,
                      std::size_t           from)
    {
      constexpr std::size_t width {64};
      const std::size_t     span  {cols - from};
      parallel_for((span + width - 1) / width, 2 * v.size() * width, [&](std::size_t chunk) {
                     const std::size_t start {from + (chunk * width)};
                     const std::size_t count {std::min(width, cols - start)};
                     buffer<T>    w(count, T {});
                     for (std::size_t i {0}; i < v.size(); ++i) {
                       const T  weight {conj_of(v[i])};
                       const T* row    {a.data() + ((first + i) * cols) + start};
                       add_product(w.data(), weight, row, count);
                     }
                     for (std::size_t i {0}; i < v.size(); ++i) {
                       const T scale {static_cast<T>(2) * v[i]};
                       T*      row   {a.data() + ((first + i) * cols) + start};
                       subtract_product(row, scale, w.data(), count);
                     }
                   });
    }
 
    template <typename T>
    struct reflector_block
    {
      std::size_t    first;  
      std::size_t    height; 
      std::size_t    width;  
      buffer<T>      v;      
      buffer<T>      t;      
    };
 
    template <typename T>
    [[nodiscard]] reflector_block<T> gather_reflectors(const std::vector<buffer<T>>&      reflectors,
                                                       std::size_t                        k0,
                                                       std::size_t                        k1,
                                                       std::size_t                        rows,
                                                       const std::vector<T>*              taus = nullptr)
    {
      reflector_block<T> block {.first = k0, .height = rows - k0, .width = k1 - k0, .v = {}, .t = {}};
      const std::size_t  b     {block.width};
      block.v.assign(block.height * b, T {});
      block.t.assign(b * b, T {});
      for (std::size_t j {0}; j < b; ++j) {
        const buffer<T>& v {reflectors[k0 + j]};
        for (std::size_t r {0}; r < v.size(); ++r) {
          block.v[((j + r) * b) + j] = v[r];
        }
      }
      // Vᴴ V once, through the matrix product kernel: its column j above the diagonal is V(:, 0…j-1)ᴴ v_j.
      buffer<T> v_adjoint(b * block.height, T {});
      for (std::size_t r {0}; r < block.height; ++r) {
        for (std::size_t j {0}; j < b; ++j) {
          v_adjoint[(j * block.height) + r] = conj_of(block.v[(r * b) + j]);
        }
      }
      buffer<T> gram(b * b, T {});
      gemm(b, b, block.height, v_adjoint.data(), block.height, block.v.data(), b, gram.data(), b);
      for (std::size_t j {0}; j < b; ++j) {
        if (reflectors[k0 + j].empty()) {
          continue;                           // the identity: a zero column of T
        }
        const T tau {taus == nullptr ? static_cast<T>(2) : (*taus)[k0 + j]};
        for (std::size_t i {0}; i < j; ++i) { // T(0…j-1, j) = −τ_j T(0…j-1, 0…j-1) z, z = V(:, 0…j-1)ᴴ v_j
          T sum {};
          for (std::size_t p {i}; p < j; ++p) {
            sum += block.t[(i * b) + p] * gram[(p * b) + j];
          }
          block.t[(i * b) + j] = -tau * sum;
        }
        block.t[(j * b) + j] = tau;
      }
      return block;
    }
 
    template <typename T>
    void apply_reflector_block(const reflector_block<T>& block,
                               bool                      adjoint,
                               buffer<T>&                a,
                               std::size_t               cols,
                               std::size_t               from,
                               std::size_t               to)
    {
      const std::size_t b                  {block.width};
      const std::size_t h                  {block.height};
      const std::size_t w                  {to - from};
      buffer<T>         v_adjoint(b * h, T {});
      buffer<T>         minus_v(h * b, T   {});
      for (std::size_t r {0}; r < h; ++r) {
        for (std::size_t j {0}; j < b; ++j) {
          v_adjoint[(j * h) + r]  = conj_of(block.v[(r * b) + j]);
          minus_v[(r * b) + j]    = -block.v[(r * b) + j];
        }
      }
      T*             top                 {a.data() + (block.first * cols) + from};
      buffer<T>      projection(b * w, T {});
      gemm(b, w, h, v_adjoint.data(), h, top, cols, projection.data(), w);
      buffer<T> scaled(b * w, T          {});
      for (std::size_t i {0}; i < b; ++i) {
        for (std::size_t p {0}; p < b; ++p) {
          const T factor {adjoint ? conj_of(block.t[(p * b) + i]) : block.t[(i * b) + p]};
          if (factor == T {}) {
            continue; // T is triangular: half of these are zero
          }
          add_product(scaled.data() + (i * w), factor, projection.data() + (p * w), w);
        }
      }
      gemm(h, w, b, minus_v.data(), b, scaled.data(), w, top, cols);
    }
 
    template <typename T>
    void apply_transposed_reflector_block(const reflector_block<T>& block,
                                          buffer<T>&                at,
                                          std::size_t               cols,
                                          std::size_t               from,
                                          std::size_t               to,
                                          bool                      adjoint = false)
    {
      const std::size_t b            {block.width};
      const std::size_t h            {block.height};
      const std::size_t r            {to - from};
      const auto        v_conjugate  {std::make_unique_for_overwrite<T[]>(h * b)}; // conj(V), h × b, written whole
      const auto        minus_v_t    {std::make_unique_for_overwrite<T[]>(b * h)}; // −Vᵀ, b × h
      const auto        t_transposed {std::make_unique_for_overwrite<T[]>(b * b)}; // Tᵀ, or conj(T)
      for (std::size_t row {0}; row < h; ++row) {
        for (std::size_t j {0}; j < b; ++j) {
          v_conjugate[(row * b) + j] = conj_of(block.v[(row * b) + j]);
          minus_v_t[(j * h) + row]   = -block.v[(row * b) + j];
        }
      }
      for (std::size_t i {0}; i < b; ++i) {
        for (std::size_t p {0}; p < b; ++p) {
          t_transposed[(p * b) + i] = adjoint ? conj_of(block.t[(p * b) + i]) : block.t[(i * b) + p];
        }
      }
      T*             left                {at.data() + (from * cols) + block.first};
      buffer<T>      projection(r * b, T {});
      gemm(r, b, h, left, cols, v_conjugate.get(), b, projection.data(), b);
      buffer<T> scaled(r * b, T          {});
      gemm(r, b, b, projection.data(), b, t_transposed.get(), b, scaled.data(), b);
      gemm(r, h, b, scaled.data(), b, minus_v_t.get(), h, left, cols);
    }
 
    template <typename T>
    requires std::is_floating_point_v<T>
    [[nodiscard]] T conjugate_dot(const T*    x,
                                  const T*    y,
                                  std::size_t count)
    {
      return dot_product(x, y, count);
    }
 
    template <typename T>
    [[nodiscard]] T conjugate_dot(const T*    x,
                                  const T*    y,
                                  std::size_t count)
    {
      return conjugate_dot_product(x, y, count);
    }
 
    template <typename T>
    void reflect_columns(T*                    panel,
                         std::size_t           height,
                         const buffer<T>&      v,
                         std::size_t           first,
                         std::size_t           from,
                         std::size_t           to,
                         const T&              s = T {2})
    {
      for (std::size_t j {from}; j < to; ++j) {
        T* const column {panel + (j * height) + first};
        const T  scaled {s * conjugate_dot(v.data(), column, v.size())};
        subtract_product(column, scaled, v.data(), v.size());
      }
    }
 
    inline constexpr std::size_t qr_block {32};
 
    template <typename T>
    void factor_qr_panel(T*                           panel,
                         std::size_t                  ld,
                         std::size_t                  top,
                         std::vector<buffer<T>>&      reflectors,
                         std::vector<T>&              taus,
                         std::size_t                  k0,
                         std::size_t                  k1)
    {
      for (std::size_t k {k0}; k < k1; ++k) {
        const std::size_t c      {k - k0};
        T* const          x      {panel + (c * ld) + (k - top)};
        buffer<T>         w(x, x + (ld + top - k));
        const double      factor {rescale_tiny(w)}; // w, τ and the signs the same at any scale; β scaled back
        const T           alpha  {w[0]};
        const double      below  {scaled_norm(w.data() + 1, w.size() - 1)};
        if (below == 0.0 && alpha == detail::from_double<T>(detail::real_part(alpha))) {
          continue; // the identity: nothing below the diagonal, and the diagonal real
        }
        const double beta  {-std::copysign(scaled_norm(w), detail::real_part(alpha))};
        const T      tau   {(detail::from_double<T>(beta) - alpha) / detail::from_double<T>(beta)};
        const T      scale {T {1} / (alpha - detail::from_double<T>(beta))}; // |α − β| ≥ ‖x‖: no cancellation
        w[0] = T {1};
        for (std::size_t r {1}; r < w.size(); ++r) {
          w[r] *= scale;
        }
        reflect_columns(panel, ld, w, k - top, c + 1, k1 - k0, conj_of(tau)); // Hᴴ on the panel's columns past k
        x[0] = detail::from_double<T>(beta / factor);
        std::copy(w.begin() + 1, w.end(), x + 1);
        reflectors[k] = std::move(w);
        taus[k]       = tau;
      }
    }
 
    template <typename T>
    struct qr_reflections
    {
      buffer<T>                       r;      
      std::vector<reflector_block<T>> blocks; 
    };
 
    template <typename T>
    [[nodiscard]] std::vector<reflector_block<T>> reflect_in_columns(buffer<T>&  columns,
                                                                     std::size_t rows,
                                                                     std::size_t cols)
    {
      const std::size_t               depth         {std::min(rows, cols)};
      std::vector<buffer<T>>          reflectors(depth); // w_k on rows k…m-1; empty for the identity
      std::vector<T>                  taus(depth, T {});
      std::vector<reflector_block<T>> blocks;
      for (std::size_t k0 {0}; k0 < depth; k0 += qr_block) {
        const std::size_t k1 {std::min(k0 + qr_block, depth)};
        factor_qr_panel(columns.data() + (k0 * rows), rows, 0, reflectors, taus, k0, k1);
        blocks.push_back(gather_reflectors(reflectors, k0, k1, rows, &taus));
        if (k1 < cols) { // right of the panel: H_{k1-1} … H_{k0} A = (I − V T Vᴴ)ᴴ A
          apply_transposed_reflector_block(blocks.back(), columns, rows, k1, cols, true);
        }
      }
      return blocks;
    }
 
    inline constexpr std::size_t qr_transpose_tile {32};
 
    template <typename T>
    void transpose_into(const T*    from,
                        std::size_t rows,
                        std::size_t cols,
                        T*          to)
    {
      parallel_for((rows + qr_transpose_tile - 1) / qr_transpose_tile, memory_element_weight * qr_transpose_tile * cols, [from, rows, cols, to](std::size_t band) {
                     const std::size_t i0 {band * qr_transpose_tile};
                     for (std::size_t j0 {0}; j0 < cols; j0 += qr_transpose_tile) {
                       for (std::size_t j {j0}; j < std::min(j0 + qr_transpose_tile, cols); ++j) {
                         for (std::size_t i {i0}; i < std::min(i0 + qr_transpose_tile, rows); ++i) {
                           to[(j * rows) + i] = from[(i * cols) + j];
                         }
                       }
                     }
                   });
    }
 
    template <typename T>
    void factor_row_major_panel(buffer<T>&              a,
                                std::size_t             rows,
                                std::size_t             cols,
                                std::vector<buffer<T>>& reflectors,
                                std::vector<T>&         taus,
                                std::size_t             k0,
                                std::size_t             k1)
    {
      const std::size_t height                  {rows - k0};
      const std::size_t width                   {k1 - k0};
      buffer<T>         panel(height * width, T {});
      for (std::size_t i {0}; i < height; ++i) {
        for (std::size_t j {0}; j < width; ++j) {
          panel[(j * height) + i] = a[((k0 + i) * cols) + k0 + j];
        }
      }
      factor_qr_panel(panel.data(), height, k0, reflectors, taus, k0, k1);
      for (std::size_t i {0}; i < height; ++i) {
        for (std::size_t j {0}; j < width; ++j) {
          a[((k0 + i) * cols) + k0 + j] = panel[(j * height) + i];
        }
      }
    }
 
    template <typename T>
    [[nodiscard]] std::vector<reflector_block<T>> reflect_in_rows(buffer<T>&  a,
                                                                  std::size_t rows,
                                                                  std::size_t cols,
                                                                  buffer<T>*  r = nullptr)
    {
      const std::size_t      depth         {std::min(rows, cols)};
      std::vector<buffer<T>> reflectors(depth); // w_k on rows k…m-1; empty for the identity
      std::vector<T>         taus(depth, T {});
      if (r != nullptr) {
        *r = buffer<T>(depth * cols, T {}); // constructed, not assigned: the generic fill of assign took 4 % of qr at order 128
      }
      std::vector<reflector_block<T>> blocks;
      for (std::size_t k0 {0}; k0 < depth; k0 += qr_block) {
        const std::size_t k1 {std::min(k0 + qr_block, depth)};
        factor_row_major_panel(a, rows, cols, reflectors, taus, k0, k1);
        blocks.push_back(gather_reflectors(reflectors, k0, k1, rows, &taus));
        if (k1 < cols) { // right of the panel: H_{k1-1} … H_{k0} A = (I − V T Vᴴ)ᴴ A
          apply_reflector_block(blocks.back(), true, a, cols, k1, cols);
        }
      }
      for (std::size_t i {0}; r != nullptr && i < depth; ++i) {
        for (std::size_t j {i}; j < cols; ++j) {
          (*r)[(i * cols) + j] = a[(i * cols) + j];
        }
      }
      return blocks;
    }
 
    template <typename T>
    [[nodiscard]] qr_reflections<T> householder_reflect(buffer<T>   a,
                                                        std::size_t rows,
                                                        std::size_t cols)
    {
#ifdef NUMERICAL_X86_COMPILED
      const std::size_t depth {std::min(rows, cols)};
      buffer<T>         columns(rows * cols);
      transpose_into(a.data(), rows, cols, columns.data());
      qr_reflections<T> out {.r = buffer<T>(depth * cols), .blocks = reflect_in_columns(columns, rows, cols)};
      for (std::size_t j0 {0}; j0 < cols; j0 += qr_transpose_tile) { // R, row-major: written whole, zeros below
        for (std::size_t i0 {0}; i0 < depth; i0 += qr_transpose_tile) {
          for (std::size_t i {i0}; i < std::min(i0 + qr_transpose_tile, depth); ++i) {
            for (std::size_t j {j0}; j < std::min(j0 + qr_transpose_tile, cols); ++j) {
              out.r[(i * cols) + j] = i <= j ? columns[(j * rows) + i] : T {};
            }
          }
        }
      }
#else
      qr_reflections<T> out;
      out.blocks = reflect_in_rows(a, rows, cols, &out.r);
#endif
      return out;
    }
 
    template <typename T>
    void form_q_into(buffer<T>&                             out,
                     const std::vector<reflector_block<T>>& blocks,
                     std::size_t                            rows,
                     std::size_t                            first,
                     std::size_t                            count,
                     bool                                   transposed)
    {
#ifdef NUMERICAL_X86_COMPILED
      buffer<T>  formed;
      buffer<T>& qt {transposed ? out : formed};
      qt = buffer<T>(count * rows, T {});
      for (std::size_t c {0}; c < count; ++c) {
        qt[(c * rows) + first + c] = T {1};
      }
      for (std::size_t b {blocks.size()}; b-- > 0;) {
        apply_transposed_reflector_block(blocks[b], qt, rows, std::max(first, blocks[b].first) - first, count);
      }
      if (!transposed) {
        out.resize(rows * count);
        transpose_into(qt.data(), count, rows, out.data());
      }
#else
      buffer<T>  formed;
      buffer<T>& q {transposed ? formed : out};
      q = buffer<T>(rows * count, T {});
      for (std::size_t c {0}; c < count; ++c) {
        q[((first + c) * count) + c] = T {1};
      }
      for (std::size_t b {blocks.size()}; b-- > 0;) {
        apply_reflector_block(blocks[b], false, q, count, std::max(first, blocks[b].first) - first, count);
      }
      if (transposed) {
        out.resize(count * rows);
        transpose_into(q.data(), rows, count, out.data());
      }
#endif
    }
 
    template <typename T>
    [[nodiscard]] thin_qr<T> householder_qr(buffer<T>      a,
                                            std::size_t    rows,
                                            std::size_t    cols,
                                            bool           complete = false)
    {
      const std::size_t depth   {std::min(rows, cols)};
      const std::size_t width   {complete ? rows : depth};
      qr_reflections<T> factors {householder_reflect(std::move(a), rows, cols)};
      // Q = H₀ H₁ … applied to the columns [first, first + count) of I, each column's own product: the first `depth`
      // formed alone whether Q is complete or not, so that the complete Q's are the reduced one's to the bit (formed
      // together, the products' tiles differ, and a fused multiply-add may round a column otherwise).
      const auto formed {[&](std::size_t first, std::size_t count) {
                           buffer<T> columns;
                           form_q_into(columns, factors.blocks, rows, first, count, false);
                           return columns;
                         }};
      thin_qr<T> out {.q = formed(0, depth), .r = std::move(factors.r)};
      if (width > depth) {
        const buffer<T> rest                  {formed(depth, width - depth)};
        buffer<T>       whole(rows * width, T {});
        for (std::size_t i {0}; i < rows; ++i) {
          std::copy_n(out.q.data() + (i * depth), depth, whole.data() + (i * width));
          std::copy_n(rest.data() + (i * (width - depth)), width - depth, whole.data() + (i * width) + depth);
        }
        out.q = std::move(whole);
      }
      for (std::size_t k {0}; k < depth; ++k) {               // R's diagonal real and non-negative: rotate row k, column k
        const T phase {phase_of(out.r[(k * cols) + k])};
        for (std::size_t j {k}; j < cols; ++j) {
          out.r[(k * cols) + j] *= conj_of(phase);
        }
        for (std::size_t i {0}; i < rows; ++i) {
          out.q[(i * width) + k] *= phase;
        }
      }
      return out;
    }
 
    [[nodiscard]] inline double rank_tolerance(double      largest,
                                               std::size_t rows,
                                               std::size_t cols)
    {
      return largest * std::numeric_limits<double>::epsilon() * static_cast<double>(std::max(rows, cols));
    }
  } // namespace detail
 
  class rank_cutoff
  {
  public:
 
    constexpr rank_cutoff() noexcept = default;
 
    [[nodiscard]] static rank_cutoff relative(double ratio)
    {
      return rank_cutoff {checked(ratio), false};
    }
 
    [[nodiscard]] static rank_cutoff absolute(double magnitude)
    {
      return rank_cutoff {checked(magnitude), true};
    }
 
    [[nodiscard]] double magnitude(double      largest,
                                   std::size_t rows,
                                   std::size_t cols) const noexcept
    {
      if (absolute_) {
        return value_;
      }
      return value_ < 0.0 ? detail::rank_tolerance(largest, rows, cols) : largest * value_;
    }
 
    [[nodiscard]] rank_cutoff scaled(int exponent) const noexcept
    {
      return absolute_ ? rank_cutoff {std::ldexp(value_, -exponent), true} : *this;
    }
 
  private:
 
    constexpr rank_cutoff(double value,
                          bool   absolute) noexcept
      : value_(value), absolute_(absolute)
    {}
 
    [[nodiscard]] static double checked(double value)
    {
      if (!(value >= 0.0)) {
        throw std::invalid_argument("rank cutoff: a tolerance is non-negative");
      }
      return value;
    }
 
    double value_    {-1.0};  
    bool   absolute_ {false}; 
  };
 
  enum class qr_mode
  {
    reduced, 
    complete 
  };
 
  template <typename T>
  class qr_decomposition
  {
  public:
 
    explicit qr_decomposition(const matrix<T>& a,
                              qr_mode          mode = qr_mode::reduced)
      : qr_decomposition(detail::parallel_copy_of(a),
                         mode)
    {}
 
    explicit qr_decomposition(matrix<T>&& a,
                              qr_mode     mode = qr_mode::reduced)
      : rows_(a.rows()), cols_(a.cols()), depth_(std::min(rows_, cols_)),
        width_(mode == qr_mode::complete ? rows_ : depth_)
    {
      detail::require_finite(a, "QR decomposition");
      detail::thin_qr<T> factors {detail::householder_qr(std::move(a).elements(), rows_, cols_, mode == qr_mode::complete)};
      q_ = std::move(factors.q);
      r_ = std::move(factors.r);
      r_.resize(width_ * cols_, T {}); // complete: R's rows past k are zero
    }
 
    [[nodiscard]] matrix<T> q() const&
    {
      return matrix<T> {rows_, width_, q_};
    }
 
    [[nodiscard]] matrix<T> q() &&
    {
      return matrix<T> {rows_, width_, std::move(q_)};
    }
 
    [[nodiscard]] matrix<T> r() const&
    {
      return matrix<T> {width_, cols_, r_};
    }
 
    [[nodiscard]] matrix<T> r() &&
    {
      return matrix<T> {width_, cols_, std::move(r_)};
    }
 
    [[nodiscard]] vector<T> solve(const vector<T>& b) const
    {
      if (b.size() != rows_) {
        throw dimension_error("solve: right-hand side size differs from the matrix");
      }
      if (rows_ < cols_) {
        throw dimension_error("QR solve: the system is underdetermined (lstsq gives its minimum-norm solution)");
      }
      double largest {0.0};
      for (std::size_t i {0}; i < cols_; ++i) {
        largest = std::max(largest, detail::real_part(r_[(i * cols_) + i]));
      }
      const double   tolerance  {detail::rank_tolerance(largest, rows_, cols_)};
      buffer<T>      y(cols_, T {}); // y = Qᴴ b
      for (std::size_t i {0}; i < cols_; ++i) {
        if (detail::real_part(r_[(i * cols_) + i]) <= tolerance) {
          throw singular_matrix_error("QR solve: the matrix is rank-deficient (lstsq gives its minimum-norm solution)");
        }
        T sum {};
        for (std::size_t k {0}; k < rows_; ++k) {
          sum += detail::conj_of(q_[(k * width_) + i]) * b[k];
        }
        y[i] = sum;
      }
      buffer<T> x(cols_, T {}); // back substitution R x = y
      for (std::size_t i {cols_}; i-- > 0;) {
        T sum {y[i]};
        for (std::size_t j {i + 1}; j < cols_; ++j) {
          sum -= r_[(i * cols_) + j] * x[j];
        }
        x[i] = sum / r_[(i * cols_) + i];
      }
      return vector<T> {std::move(x)};
    }
 
  private:
 
    std::size_t    rows_;  
    std::size_t    cols_;  
    std::size_t    depth_; 
    std::size_t    width_; 
    buffer<T>      q_;     
    buffer<T>      r_;     
  };
 
  template <typename T>
  [[nodiscard]] qr_decomposition<T> qr(const matrix<T>& a)
  {
    return qr_decomposition<T> {a};
  }
 
  class convergence_error : public std::runtime_error
  {
  public:
 
    explicit convergence_error(const std::string& message)
      : std::runtime_error(message)
    {}
  };
 
  namespace detail {
 
    inline constexpr std::size_t jacobi_sweeps {60};
 
    template <typename T>
    requires std::is_floating_point_v<T>
    void add_conjugate_product(T*          y,
                               const T*    x,
                               const T&    f,
                               std::size_t count)
    {
      add_product(y, f, x, count);
    }
 
    template <typename T>
    void add_conjugate_product(T*          y,
                               const T*    x,
                               const T&    f,
                               std::size_t count)
    {
      using R = typename T::value_type;
      const R c {f.real()};
      const R d {f.imag()};
      for (std::size_t i {0}; i < count; ++i) {
        const R a  {x[i].real()};
        const R b  {-x[i].imag()};
        const R ac {a * c};
        const R bd {b * d};
        const R ad {a * d};
        const R bc {b * c};
        const R re {ac - bd};
        const R im {ad + bc};
        y[i] = T {y[i].real() + re, y[i].imag() + im};
      }
    }
 
    using symmetric_row_kernel = double (*)(const double*, const double*, double, double*, std::size_t);
 
    using symmetric_rows_kernel = std::array<double, symmetric_row_group> (*)(const double*, std::size_t, const double*,
                                                                              double*, std::size_t);
 
    template <typename T>
    [[nodiscard]] T hermitian_row(symmetric_row_kernel /*kernel*/,
                                  const T*             row,
                                  const T*             v,
                                  T*                   out,
                                  std::size_t          i)
    {
      add_conjugate_product(out, row, v[i], i);
      return dot_product(row, v, i + 1);
    }
 
#ifdef NUMERICAL_NEON_COMPILED
    inline std::complex<double> neon_hermitian_row(const double* row,
                                                   const double* v,
                                                   double        vr,
                                                   double        vi,
                                                   double*       out,
                                                   std::size_t   count) noexcept
    {
      const float64x2_t real      {vr, -vr};
      const float64x2_t imag      {vi, vi};
      float64x2_t       same[2]   {vdupq_n_f64(0.0), vdupq_n_f64(0.0)};
      float64x2_t       across[2] {vdupq_n_f64(0.0), vdupq_n_f64(0.0)};
      std::size_t       j         {0};
      for (; j + 2 <= count; j += 2) {
        for (std::size_t l {0}; l < 2; ++l) {
          const float64x2_t r {vld1q_f64(row + (2 * (j + l)))};
          const float64x2_t x {vld1q_f64(v + (2 * (j + l)))};
          const float64x2_t s {vextq_f64(r, r, 1)};
          float64x2_t       o {vld1q_f64(out + (2 * (j + l)))};
          o         = vfmaq_f64(o, real, r);
          o         = vfmaq_f64(o, imag, s);
          same[l]   = vfmaq_f64(same[l], r, x);
          across[l] = vfmaq_f64(across[l], s, x);
          vst1q_f64(out + (2 * (j + l)), o);
        }
      }
      const float64x2_t sum_same   {vaddq_f64(same[0], same[1])};
      const float64x2_t sum_across {vaddq_f64(across[0], across[1])};
      double            re         {vgetq_lane_f64(sum_same, 0) - vgetq_lane_f64(sum_same, 1)};
      double            im         {vgetq_lane_f64(sum_across, 0) + vgetq_lane_f64(sum_across, 1)};
      for (; j < count; ++j) {
        const double rr {row[2 * j]};
        const double ri {row[(2 * j) + 1]};
        out[2 * j]       = std::fma(vi, ri, std::fma(vr, rr, out[2 * j]));
        out[(2 * j) + 1] = std::fma(vi, rr, std::fma(-vr, ri, out[(2 * j) + 1]));
        re               = std::fma(-ri, v[(2 * j) + 1], std::fma(rr, v[2 * j], re));
        im               = std::fma(rr, v[(2 * j) + 1], std::fma(ri, v[2 * j], im));
      }
      return std::complex<double> {re, im};
    }
 
#endif
 
    [[nodiscard]] inline std::complex<double> hermitian_row(symmetric_row_kernel        kernel,
                                                            const std::complex<double>* row,
                                                            const std::complex<double>* v,
                                                            std::complex<double>*       out,
                                                            std::size_t                 i)
    {
#ifdef NUMERICAL_X86_COMPILED
      // A complex number's parts are an array of two (the standard's layout guarantee): addressed as doubles.
      const gemm_path path {gemm_path_current()};
      if (path == gemm_path::avx512 || path == gemm_path::avx2) {
        const auto* const          values {reinterpret_cast<const double*>(v)};   // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the layout guarantee above
        const auto* const          parts  {reinterpret_cast<const double*>(row)}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the layout guarantee above
        auto* const                into   {reinterpret_cast<double*>(out)};       // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the layout guarantee above
        const std::complex<double> lower  {path == gemm_path::avx512
                                          ? avx512_hermitian_row(parts, values, v[i].real(), v[i].imag(), into, i)
                                          : avx2_hermitian_row(parts, values, v[i].real(), v[i].imag(), into, i)};
        return lower + (row[i] * v[i]);
      }
#endif
#ifdef NUMERICAL_NEON_COMPILED
      if (const gemm_path neon_path {gemm_path_current()}; neon_path == gemm_path::neon || neon_path == gemm_path::amx) {
        const auto* const          values {reinterpret_cast<const double*>(v)};   // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the layout guarantee above
        const auto* const          parts  {reinterpret_cast<const double*>(row)}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the layout guarantee above
        auto* const                into   {reinterpret_cast<double*>(out)};       // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the layout guarantee above
        const std::complex<double> lower  {neon_hermitian_row(parts, values, v[i].real(), v[i].imag(), into, i)};
        return lower + (row[i] * v[i]);
      }
#endif
      return hermitian_row<std::complex<double>>(kernel, row, v, out, i);
    }
 
    [[nodiscard]] inline double hermitian_row(symmetric_row_kernel kernel,
                                              const double*        row,
                                              const double*        v,
                                              double*              out,
                                              std::size_t          i)
    {
      return kernel(row, v, v[i], out, i) + (row[i] * v[i]);
    }
 
    template <typename T>
    [[nodiscard]] std::size_t symmetric_row_speedup()
    {
#ifdef NUMERICAL_X86_COMPILED
      if constexpr (std::is_same_v<T, double>) {
        if (gemm_path_current() != gemm_path::portable) {
          return vector_speedup;
        }
      }
#endif
      return 1;
    }
 
    template <typename T>
    void hermitian_product(const T*        block,
                           std::size_t     n,
                           std::size_t     m,
                           const T*        v,
                           T*              w,
                           buffer<T>&      scratch)
    {
      constexpr std::size_t rows_per_task {32};
      const std::size_t     tasks         {(m + rows_per_task - 1) / rows_per_task};
      if (scratch.size() < tasks * m) {
        scratch.resize(tasks * m, T {});
      }
      symmetric_row_kernel kernel {[](const double* row, const double* x, double f, double* out, std::size_t count) {
                                     return symmetric_row(row, x, f, out, count);
                                   }};
      symmetric_rows_kernel group {[](const double* rows, std::size_t ld, const double* x, double* out, std::size_t i) {
                                     return symmetric_rows(rows, ld, x, out, i);
                                   }};
#ifdef NUMERICAL_X86_COMPILED
      if constexpr (std::is_same_v<T, double>) {
        const gemm_path path {gemm_path_current()};
        if (path == gemm_path::avx512) {
          kernel = &avx512_symmetric_row;
          group  = &avx512_symmetric_rows;
        }
        else if (path == gemm_path::avx2) {
          kernel = &avx2_symmetric_row;
          group  = &avx2_symmetric_rows;
        }
      }
#endif
      parallel_for(tasks, rows_per_task * m * symmetric_row_speedup<T>(), [&](std::size_t task) {
                     const std::size_t first {task * rows_per_task};
                     const std::size_t last {std::min(first + rows_per_task, m)};
                     T*                own {scratch.data() + (task * m)};
                     std::fill(own, own + last, T {}); // a task's rows reach w[0…last) alone
                     std::size_t i {first};
                     if constexpr (std::is_same_v<T, double>) { // rows by groups: v and own read once a group
                       for (; i + symmetric_row_group <= last; i += symmetric_row_group) {
                         const std::array<double, symmetric_row_group> sums {group(block + (i * n), n, v, own, i)};
                         for (std::size_t r {0}; r < symmetric_row_group; ++r) {
                           own[i + r] += sums[r];
                         }
                       }
                     }
                     for (; i < last; ++i) {
                       own[i] += hermitian_row(kernel, block + (i * n), v, own, i);
                     }
                   });
      std::fill(w, w + m, T {});
      for (std::size_t task {0}; task < tasks; ++task) { // the partials in order, each a contiguous update
        add_product(w, T {1}, scratch.data() + (task * m), std::min((task + 1) * rows_per_task, m));
      }
    }
 
#ifdef NUMERICAL_X86_COMPILED
    inline constexpr std::size_t corrected_product_least {256};
#else
    inline constexpr std::size_t corrected_product_least {512};
#endif
 
    template <typename T>
    requires std::is_floating_point_v<T>
    void subtract_corrections(T*          y,
                              const T*    f,
                              const T*    x,
                              std::size_t ldx,
                              std::size_t terms,
                              std::size_t count)
    {
      subtract_products(y, f, x, ldx, terms, count);
    }
 
    template <split_complex T>
    void subtract_corrections(T*          y,
                              const T*    f,
                              const T*    x,
                              std::size_t ldx,
                              std::size_t terms,
                              std::size_t count)
    {
#ifdef NUMERICAL_X86_COMPILED
      if constexpr (std::is_same_v<T, std::complex<double>>) { // its own fused kernel on the x86-64 units, as fast
        if (const gemm_path path {gemm_path_current()}; path == gemm_path::avx2 || path == gemm_path::avx512) {
          subtract_products(y, f, x, ldx, terms, count);
          return;
        }
      }
#endif
      using R = typename T::value_type;
      for (std::size_t t {0}; t < terms; ++t) {
        const R        c   {f[t].real()};
        const R        d   {f[t].imag()};
        const T* const row {x + (t * ldx)};
        for (std::size_t i {0}; i < count; ++i) {
          const R a  {row[i].real()};
          const R b  {row[i].imag()};
          const R ac {c * a};
          const R bd {d * b};
          const R ad {c * b};
          const R bc {d * a};
          const R re {ac - bd};
          const R im {ad + bc};
          y[i] = T {y[i].real() - re, y[i].imag() - im};
        }
      }
    }
 
    template <typename T>
    void corrected_hermitian_product(const T*        block,
                                     std::size_t     n,
                                     std::size_t     m,
                                     const T*        v,
                                     T*              w,
                                     buffer<T>&      scratch,
                                     const T*        corrections,
                                     std::size_t     ldc,
                                     std::size_t     terms,
                                     T*              factors,
                                     T*              next         = nullptr,
                                     const T*        next_factors = nullptr)
    {
      constexpr std::size_t rows_per_task {32};
      const std::size_t     tasks         {(m + rows_per_task - 1) / rows_per_task};
      if (scratch.size() < tasks * m) {
        scratch.resize(tasks * m, T {});
      }
      symmetric_row_kernel kernel {[](const double* row, const double* x, double f, double* out, std::size_t count) {
                                     return symmetric_row(row, x, f, out, count);
                                   }};
      symmetric_rows_kernel group {[](const double* rows, std::size_t ld, const double* x, double* out, std::size_t i) {
                                     return symmetric_rows(rows, ld, x, out, i);
                                   }};
#ifdef NUMERICAL_X86_COMPILED
      if constexpr (std::is_same_v<T, double>) {
        const gemm_path path {gemm_path_current()};
        if (path == gemm_path::avx512) {
          kernel = &avx512_symmetric_row;
          group  = &avx512_symmetric_rows;
        }
        else if (path == gemm_path::avx2) {
          kernel = &avx2_symmetric_row;
          group  = &avx2_symmetric_rows;
        }
      }
#endif
      // Below corrected_product_least rows the coefficients and the bands run on the caller: woken threads cost more
      // than they take (an eigh of order 128, six x86-64 cores, +13 %).
      const bool        share {m >= corrected_product_least};
      const std::size_t pairs {terms / 2};
      const auto        coefficients = [&](std::size_t p) {
                                         factors[2 * p]       = conjugate_dot(corrections + (((2 * p) + 1) * ldc), v, m);
                                         factors[(2 * p) + 1] = conjugate_dot(corrections + ((2 * p) * ldc), v, m);
                                       };
      parallel_for(tasks + (share ? pairs : 0), rows_per_task * m * symmetric_row_speedup<T>(), [&](std::size_t task) {
                     if (task >= tasks) { // a pair of coefficients: an earlier step's w and v against v
                       coefficients(task - tasks);
                       return;
                     }
                     const std::size_t first {task * rows_per_task};
                     const std::size_t last {std::min(first + rows_per_task, m)};
                     T*                own {scratch.data() + (task * m)};
                     std::fill(own, own + last, T {});
                     std::size_t i {first};
                     if constexpr (std::is_same_v<T, double>) {
                       for (; i + symmetric_row_group <= last; i += symmetric_row_group) {
                         const std::array<double, symmetric_row_group> sums {group(block + (i * n), n, v, own, i)};
                         for (std::size_t r {0}; r < symmetric_row_group; ++r) {
                           own[i + r] += sums[r];
                         }
                       }
                     }
                     for (; i < last; ++i) {
                       own[i] += hermitian_row(kernel, block + (i * n), v, own, i);
                     }
                   });
      for (std::size_t p {0}; !share && p < pairs; ++p) {
        coefficients(p);
      }
      constexpr std::size_t band {256};
      const auto            bands = [&](std::size_t b) {
                                      const std::size_t first {b * band};
                                      const std::size_t count {std::min(band, m - first)};
                                      std::fill(w + first, w + first + count, T {});
                                      for (std::size_t task {0}; task < tasks; ++task) { // the partials in order, as hermitian_product
                                        const std::size_t reach {std::min((task + 1) * rows_per_task, m)};
                                        if (reach > first) {
                                          add_product(w + first, T {1}, scratch.data() + (task * m) + first, std::min(count, reach - first));
                                        }
                                      }
                                      subtract_corrections(w + first, factors, corrections + first, ldc, terms, count);
                                      if (next != nullptr) { // the next step's column: these earlier steps' terms (next implies an earlier step)
                                        subtract_corrections(next + first, next_factors, corrections + first, ldc, terms, count);
                                      }
                                    };
      if (share) {
        parallel_for((m + band - 1) / band, band * (tasks + terms) * 8, bands);
        return;
      }
      for (std::size_t b {0}; b < (m + band - 1) / band; ++b) {
        bands(b);
      }
    }
 
    inline constexpr std::size_t tridiagonal_block {32};
 
    inline constexpr std::size_t tridiagonal_band {512};
 
    inline constexpr std::size_t ql_steps_per_value {30};
 
    inline double real_subdiagonal(double t,
                                   double /*delta*/,
                                   double& /*next*/)
    {
      return t;
    }
 
    inline double real_subdiagonal(std::complex<double>  t,
                                   std::complex<double>  delta,
                                   std::complex<double>& next)
    {
      next = delta * std::conj(phase_of(t));
      return std::abs(t);
    }
 
    template <typename T>
    struct tridiagonal_reduction
    {
      std::vector<buffer<T>> reflectors; 
      buffer<T>              delta;      
    };
 
    template <typename T>
    tridiagonal_reduction<T> tridiagonalize(buffer<T>&           a,
                                            std::size_t          n,
                                            buffer<double>&      d,
                                            buffer<double>&      e)
    {
      std::vector<buffer<T>> reflectors(n);                                                                  // reflectors[k + 1] zeroes column k, acting on rows k + 1…
      buffer<T>              symv_scratch;                                                                   // hermitian_product's partials, kept from step to step
      buffer<T>              sub(n, T      {});                                                              // T's subdiagonal before the phases
      // The trailing update's factors, for the first panel's size, the largest: one buffer each, mapped once.
      const auto left                      {std::make_unique_for_overwrite<T[]>(n * 2 * tridiagonal_block)}; // row i: −V's then −W's
      const auto right                     {std::make_unique_for_overwrite<T[]>(2 * tridiagonal_block * n)}; // Wᴴ's rows, then Vᴴ's
      // By reflection: row 2p holds step k0 + p's vector v over the rows (zero above), row 2p + 1 its w, so that a
      // step's terms of the panel's earlier ones are one update through the rows (subtract_products), in their order.
      buffer<T> by_row(2 * tridiagonal_block * n);
      buffer<T> columns(tridiagonal_block * n);      // the panel's columns of A, gathered row by row
      buffer<T> factors(2 * tridiagonal_block);      // a step's coefficients of those rows
      buffer<T> next_factors(2 * tridiagonal_block); // the next column's, applied ahead
      bool      column_ahead {false};                // whether the step's column took its earlier terms ahead
      for (std::size_t k0 {0}; k0 + 1 < n; k0 += tridiagonal_block) {
        // A panel of reflections, as LAPACK's dlatrd: the matrix the panel's later steps see is
        // A − V Wᴴ − W Vᴴ, applied to the one column and the one product each step needs; the whole trailing
        // block takes it after the panel, through the matrix product kernel.
        const std::size_t k1                 {std::min(k0 + tridiagonal_block, n - 1)};
        const std::size_t width              {k1 - k0};
        std::fill_n(by_row.data(), 2 * width * n, T {});
        for (std::size_t i {k0}; i < n; ++i) { // column k from row k down, each row's part read once
          for (std::size_t k {k0}; k < std::min(k1, i + 1); ++k) {
            columns[((k - k0) * n) + i] = a[(i * n) + k];
          }
        }
        for (std::size_t k {k0}; k < k1; ++k) {
          const std::size_t j               {k - k0};
          const std::size_t m               {n - k - 1};
          T* const          column          {columns.data() + (j * n) + k}; // column k from the diagonal down, the panel's earlier steps applied
          // The steps before the last were applied during it (column_ahead); the last's two terms, or all, here.
          const std::size_t from            {column_ahead ? j - 1 : 0};
          for (std::size_t p {from}; p < j; ++p) {
            factors[2 * (p - from)]       = conj_of(by_row[(((2 * p) + 1) * n) + k]); // v's row by w's element
            factors[(2 * (p - from)) + 1] = conj_of(by_row[((2 * p) * n) + k]);
          }
          subtract_corrections(column, factors.data(), by_row.data() + ((2 * from) * n) + k, n, 2 * (j - from), m + 1);
          column_ahead   = false;
          a[(k * n) + k] = column[0];
          buffer<T>      v(column + 1, column + m + 1); // column k below the diagonal
          const double   below {m > 1 ? scaled_norm(v.data() + 1, m - 1) : 0.0};
          if (below == 0.0) {
            sub[k] = v[0]; // already tridiagonal in this column: nothing to reflect
            continue;
          }
          const double up    {rescale_tiny(v)};                                  // the reflector is the same at any scale; α is scaled back
          const T      alpha {-phase_of(v[0]) * static_cast<T>(scaled_norm(v))}; // the sign that avoids cancellation
          sub[k]  = alpha / up;
          v[0]   -= alpha;
          const double length {scaled_norm(v)};                                  // real: each part divided, where a complex division scaled and divided both
          for (T& element : v) {
            element /= length;
          }
          const T*       block  {a.data() + ((k + 1) * n) + k + 1};         // the trailing block, row stride n
          buffer<T>      w(m, T {});
          // w = (A − V Wᴴ − W Vᴴ) v on the trailing block, (Wᴴ v)ₚ and (Vᴴ v)ₚ beside A v; the next column of the
          // panel takes these earlier steps' terms in the same pass (its rows are w's).
          T* const next {k + 1 < k1 && j > 0 ? columns.data() + ((j + 1) * n) + k + 1 : nullptr};
          for (std::size_t p {0}; next != nullptr && p < j; ++p) {
            next_factors[2 * p]       = conj_of(by_row[(((2 * p) + 1) * n) + k + 1]);
            next_factors[(2 * p) + 1] = conj_of(by_row[((2 * p) * n) + k + 1]);
          }
          corrected_hermitian_product(block, n, m, v.data(), w.data(), symv_scratch, by_row.data() + k + 1, n, 2 * j, factors.data(),
                                      next, next_factors.data());
          column_ahead = next != nullptr;
          const T c {conjugate_dot(v.data(), w.data(), m)};
          for (std::size_t i {0}; i < m; ++i) { // w = 2 B v − 2 (vᴴ B v) v, so H B H = B − v wᴴ − w vᴴ
            w[i]                                    = (static_cast<T>(2) * w[i]) - (static_cast<T>(2) * c * v[i]);
            by_row[((2 * j) * n) + k + 1 + i]       = v[i];
            by_row[(((2 * j) + 1) * n) + k + 1 + i] = w[i];
          }
          reflectors[k + 1] = std::move(v); // acting from row k + 1, as householder_qr's from row k
        }
        // The trailing block from row and column k1: A −= V Wᴴ + W Vᴴ, one product of depth 2 · width,
        // [−V −W] [Wᴴ; Vᴴ], on the lower triangle alone, as LAPACK's dsyr2k: the reduction reads nothing above the
        // diagonal. A real product the kernel packs runs on the triangle's tiles alone (\ref gemm_triangle, each
        // element the bands' sum); a smaller or complex one by bands of rows, each to the diagonal, past which a band
        // computes what is never read.
        const std::size_t rest  {n - k1};
        const std::size_t depth {2 * width};
        for (std::size_t i {0}; i < rest; ++i) { // each row of left written whole, in order
          for (std::size_t p {0}; p < width; ++p) {
            left[(i * depth) + p]         = -by_row[((2 * p) * n) + k1 + i];
            left[(i * depth) + width + p] = -by_row[(((2 * p) + 1) * n) + k1 + i];
          }
        }
        for (std::size_t p {0}; p < width; ++p) {
          for (std::size_t i {0}; i < rest; ++i) {
            right[(p * rest) + i]           = conj_of(by_row[(((2 * p) + 1) * n) + k1 + i]);
            right[((width + p) * rest) + i] = conj_of(by_row[((2 * p) * n) + k1 + i]);
          }
        }
        T* trailing {a.data() + (k1 * n) + k1};
        if constexpr (std::is_same_v<T, double>) {
          if (depth >= gemm_min_depth && rest * rest >= ((gemm_packed_threshold + depth) - 1) / depth) {
            gemm_triangle<a_rows, b_rows>(rest, depth, left.get(), depth, right.get(), rest, trailing, n, true);
            continue;
          }
        }
        for (std::size_t r0 {0}; r0 < rest; r0 += tridiagonal_band) {
          const std::size_t rows {std::min(tridiagonal_band, rest - r0)};
          gemm(rows, r0 + rows, depth, left.get() + (r0 * depth), depth, right.get(), rest, trailing + (r0 * n), n);
        }
      }
      for (std::size_t i {0}; i < n; ++i) {
        d[i] = real_part(a[(i * n) + i]);
      }
      e[n - 1] = 0.0;
      // D T Dᴴ is real with D = diag(δ): δ₀ = 1, δ_{k+1} = δ_k · conj(phase(t_k)); then M = Q Dᴴ.
      buffer<T> delta(n, T {1});
      for (std::size_t k {0}; k + 1 < n; ++k) {
        e[k] = real_subdiagonal(sub[k], delta[k], delta[k + 1]);
      }
      return {.reflectors = std::move(reflectors), .delta = std::move(delta)};
    }
 
    template <typename T>
    [[nodiscard]] buffer<T> tridiagonal_basis(const tridiagonal_reduction<T>& reduction,
                                              std::size_t                     n)
    {
      buffer<T> mt(n * n, T {});
      for (std::size_t i {0}; i < n; ++i) {
        mt[(i * n) + i] = T {1};
      }
      for (std::size_t k1 {n}; k1 > 1;) {
        const std::size_t k0 {k1 - std::min(k1 - 1, qr_block)};
        apply_transposed_reflector_block(gather_reflectors(reduction.reflectors, k0, k1, n), mt, n, k0, n);
        k1 = k0;
      }
      for (std::size_t j {0}; j < n; ++j) { // row j of Mᵀ = Dᴴ Qᵀ: row j of Qᵀ times conj(δⱼ)
        const T factor {conj_of(reduction.delta[j])};
        for (std::size_t i {0}; i < n; ++i) {
          mt[(j * n) + i] *= factor;
        }
      }
      return mt;
    }
 
    template <typename T>
    void to_matrix_basis(const tridiagonal_reduction<T>& reduction,
                         buffer<T>&                      at,
                         std::size_t                     n)
    {
      if constexpr (!std::is_floating_point_v<T>) { // a real reduction's phases are all 1
        parallel_for(n, n, [&](std::size_t i) {
                       for (std::size_t j {0}; j < n; ++j) {
                         at[(i * n) + j] *= conj_of(reduction.delta[j]);
                       }
                     });
      }
      for (std::size_t k1 {n}; k1 > 1;) {
        const std::size_t k0 {k1 - std::min(k1 - 1, qr_block)};
        apply_transposed_reflector_block(gather_reflectors(reduction.reflectors, k0, k1, n), at, n, 0, n);
        k1 = k0;
      }
    }
 
    struct row_operation
    {
      enum class kind : std::uint8_t
      {
        rotate, 
        negate, 
        swap,   
      };
 
      kind        what; 
      std::size_t i;    
      std::size_t j;    
      double      c;    
      double      s;    
    };
 
    struct rotation_run
    {
      std::size_t first; 
      std::size_t count; 
      bool        down;  
    };
 
    [[nodiscard]] inline std::vector<rotation_run> rotation_runs(const std::vector<row_operation>& operations)
    {
      std::vector<rotation_run> runs;
      for (std::size_t t {0}; t < operations.size();) {
        const row_operation& op    {operations[t]};
        std::size_t          count {1};
        bool                 down  {false};
        if (op.what == row_operation::kind::rotate && t + 1 < operations.size()
            && operations[t + 1].what == row_operation::kind::rotate) {
          down = operations[t + 1].j == op.i;
          for (; t + count < operations.size(); ++count) {
            const row_operation& previous {operations[t + count - 1]};
            const row_operation& next     {operations[t + count]};
            const bool           chained  {down ? next.j == previous.i : next.i == previous.j};
            if (next.what != row_operation::kind::rotate || !chained) {
              break;
            }
          }
        }
        runs.push_back({.first = t, .count = count, .down = down});
        t                     += count;
      }
      return runs;
    }
 
    template <typename T>
    void apply_row_operations(buffer<T>&                        m,
                              std::size_t                       n,
                              const std::vector<row_operation>& operations)
    {
      constexpr std::size_t           width {64};
      const std::vector<rotation_run> runs  {rotation_runs(operations)};
      parallel_for((n + width - 1) / width, operations.size() * width, [&](std::size_t chunk) {
                     const std::size_t     first {chunk * width};
                     const std::size_t     count {std::min(first + width, n) - first};
                     std::array<T, width>  carry {};
                     for (const rotation_run& run : runs) {
                       const row_operation& start {operations[run.first]};
                       if (run.count > 1) {
                         // Down: rotation u's first row is rotation u+1's second; up: its second is the next's first.
                         const T* source {m.data() + ((run.down ? start.j : start.i) * n) + first};
                         std::copy(source, source + count, carry.begin());
                         for (std::size_t u {run.first}; u < run.first + run.count; ++u) {
                           const row_operation& op {operations[u]};
                           const T              c  {detail::from_double<T>(op.c)};
                           const T              s  {detail::from_double<T>(op.s)};
                           if (run.down) { // the carried row is op.j: op.i is read, op.j written
                             const T* row   {m.data() + (op.i * n) + first};
                             T*       other {m.data() + (op.j * n) + first};
                             for (std::size_t k {0}; k < count; ++k) {
                               const T a {row[k]};
                               const T b {carry[k]};
                               other[k] = (c * b) - (s * a);
                               carry[k] = (c * a) + (s * b);
                             }
                           }
                           else { // the carried row is op.i: op.j is read, op.i written
                             T*       row   {m.data() + (op.i * n) + first};
                             const T* other {m.data() + (op.j * n) + first};
                             for (std::size_t k {0}; k < count; ++k) {
                               const T a {carry[k]};
                               const T b {other[k]};
                               row[k]   = (c * a) + (s * b);
                               carry[k] = (c * b) - (s * a);
                             }
                           }
                         }
                         const row_operation& end    {operations[run.first + run.count - 1]};
                         T*                   target {m.data() + ((run.down ? end.i : end.j) * n) + first};
                         std::copy(carry.begin(), carry.begin() + static_cast<std::ptrdiff_t>(count), target);
                         continue;
                       }
                       T* row   {m.data() + (start.i * n) + first};
                       T* other {m.data() + (start.j * n) + first};
                       // The kind is tested outside the column loop, which then vectorises.
                       if (start.what == row_operation::kind::rotate) {
                         const T c {detail::from_double<T>(start.c)};
                         const T s {detail::from_double<T>(start.s)};
                         for (std::size_t k {0}; k < count; ++k) {
                           const T a {row[k]};
                           const T b {other[k]};
                           row[k]   = (c * a) + (s * b);
                           other[k] = (c * b) - (s * a);
                         }
                       }
                       else if (start.what == row_operation::kind::negate) {
                         for (std::size_t k {0}; k < count; ++k) {
                           row[k] = -row[k];
                         }
                       }
                       else {
                         std::swap_ranges(row, row + count, other);
                       }
                     }
                   });
    }
 
    template <typename T>
    void tridiagonal_ql(buffer<double>&      d,
                        buffer<double>&      e,
                        buffer<T>&           zt,
                        std::size_t          n,
                        std::size_t          max_steps)
    {
      const double epsilon {std::numeric_limits<double>::epsilon()};
      double       shift   {0.0};
      double       largest {0.0};
      std::size_t  steps   {0};
      // The iteration reads only d and e: its rotations of zt's rows are recorded and applied once at the
      // end, in parallel by columns.
      std::vector<row_operation> on_zt;
      const bool                 record {!zt.empty()}; // the eigenvalues alone record nothing
      for (std::size_t l {0}; l < n; ++l) {
        largest = std::max(largest, std::abs(d[l]) + std::abs(e[l]));
        std::size_t m {l};
        while (std::abs(e[m]) > epsilon * largest) { // e[n-1] = 0 stops it
          ++m;
        }
        while (m > l && std::abs(e[l]) > epsilon * largest) {
          if (steps == max_steps) {
            throw convergence_error("symmetric eigendecomposition: QL did not converge within its steps");
          }
          ++steps;
          const double g0   {d[l]};
          double       p    {(d[l + 1] - g0) / (2.0 * e[l])};
          double       r    {rotation_length(p, 1.0)};
          r                  = p < 0.0 ? -r : r;
          d[l]               = e[l] / (p + r);
          d[l + 1]           = e[l] * (p + r);
          const double dl1  {d[l + 1]};
          const double h0   {g0 - d[l]};
          for (std::size_t i {l + 2}; i < n; ++i) {
            d[i] -= h0;
          }
          shift             += h0;
          p                  = d[m];
          double       c    {1.0};
          double       c2   {1.0};
          double       c3   {1.0};
          const double el1  {e[l + 1]};
          double       s    {0.0};
          double       s2   {0.0};
          for (std::size_t i {m}; i-- > l;) {
            c3       = c2;
            c2       = c;
            s2       = s;
            const double g {c * e[i]};
            const double h {c * p};
            r        = rotation_length(p, e[i]);
            e[i + 1] = s * r;
            s        = e[i] / r;
            c        = p / r;
            p        = (c * d[i]) - (s * g);
            d[i + 1] = h + (s * ((c * g) + (s * d[i])));
            // rows i and i+1 (columns i and i+1 of Z): upper ← c upper − s lower, lower ← c lower + s upper
            if (record) {
              on_zt.push_back({row_operation::kind::rotate, i, i + 1, c, -s});
            }
          }
          p    = -s * s2 * c3 * el1 * e[l] / dl1;
          e[l] = s * p;
          d[l] = c * p;
        }
        d[l] += shift;
        e[l]  = 0.0;
      }
      if (record) {
        apply_row_operations(zt, n, on_zt);
      }
    }
 
    inline constexpr std::size_t divide_leaf {32};
 
    inline constexpr std::size_t secular_iterations {100};
 
    struct secular_solution
    {
      std::size_t origin; 
      double      offset; 
    };
 
    [[nodiscard]] inline double quadratic_step(double a,
                                               double b,
                                               double c)
    {
      const double root {std::sqrt(std::abs((a * a) - (4.0 * b * c)))};
      return a <= 0.0 ? (a - root) / (2.0 * c) : (2.0 * b) / (a + root);
    }
 
    template <typename Offset>
    [[nodiscard]] secular_solution solve_secular(const buffer<double>&      z,
                                                 double                     rho,
                                                 std::size_t                i,
                                                 std::size_t                max_iterations,
                                                 double*                    delta,
                                                 const Offset&              offset)
    {
      const std::size_t k      {z.size()};
      const double      rhoinv {1.0 / rho};
      const bool        last   {i + 1 == k};
      double            weight {0.0};
      for (const double element : z) {
        weight += element * element;
      }
      std::size_t origin {i};
      double      lower  {0.0};
      double      upper  {2.0 * rho * weight}; // the last root is at most ρ‖z‖² past its pole: the equation is positive here
      if (!last) {
        const double gap {offset(i, i + 1)};
        double       f   {rhoinv};
        for (std::size_t j {0}; j < k; ++j) {
          f += z[j] * z[j] / (offset(i, j) - (gap / 2.0));
        }
        upper = gap / 2.0;
        if (f < 0.0) { // the root is nearer the pole past it
          origin = i + 1;
          lower  = -gap / 2.0;
          upper  = 0.0;
        }
      }
      double tau {(lower + upper) / 2.0};
      for (std::size_t iteration {0}; iteration < max_iterations; ++iteration) {
        double psi  {0.0}; // the terms of the poles up to i, and their derivative
        double dpsi {0.0};
        double phi  {0.0}; // those of the poles past i
        double dphi {0.0};
        // The two sides apart, each in its order: a term chose its sum by a branch (an SVD of order 128 on one arm64
        // core: 1.87 → 1.68 ms, eigh 1.02 → 0.82).
        for (std::size_t j {0}; j <= i; ++j) {
          delta[j] = offset(origin, j) - tau;
          const double t {z[j] / delta[j]};
          psi  += z[j] * t;
          dpsi += t * t;
        }
        for (std::size_t j {i + 1}; j < k; ++j) {
          delta[j] = offset(origin, j) - tau;
          const double t {z[j] / delta[j]};
          phi  += z[j] * t;
          dphi += t * t;
        }
        const double w {rhoinv + psi + phi};
        if (std::abs(w) <= std::numeric_limits<double>::epsilon()
            * ((8.0 * (phi - psi + rhoinv)) + (3.0 * std::abs(tau) * (dpsi + dphi)))) {
          return {.origin = origin, .offset = tau};
        }
        (w < 0.0 ? lower : upper) = tau; // the equation increases with x
        // The fit c + s₁/(δᵢ − η) + s₂/(δᵢ₊₁ − η) matching the value and the slope of each side: C η² − A η + B = 0
        // (the last root has no pole to its right: its δᵢ₊₁ term is absent, A = δᵢ w − δᵢ² ψ' over C = w − δᵢ ψ').
        const double a1   {delta[i]};
        const double a2   {last ? 0.0 : delta[i + 1]};
        const double c    {w - (a1 * dpsi) - (a2 * dphi)};
        const double a    {last ? (a1 * c) + (a1 * a1 * dpsi) : ((a1 + a2) * w) - (a1 * a2 * (dpsi + dphi))};
        const double b    {last ? 0.0 : a1 * a2 * w};
        const double next {tau + (last ? a / c : quadratic_step(a, b, c))};
        // A step out of the bracket (or not a number): bisect. The test above stops within a few units in the
        // last place of the root, before the bracket runs out of points.
        tau = next > lower && next < upper ? next : (lower + upper) / 2.0;
      }
      throw convergence_error("the secular equation did not converge");
    }
 
    [[nodiscard]] inline double secular_root(const buffer<double>&      d,
                                             const buffer<double>&      z,
                                             double                     rho,
                                             std::size_t                i,
                                             std::size_t                max_iterations,
                                             double*                    delta)
    {
      const secular_solution root {solve_secular(z, rho, i, max_iterations, delta,
                                                 [&d](std::size_t o, std::size_t j) { return d[j] - d[o]; })};
      return d[root.origin] + root.offset;
    }
 
    inline constexpr std::size_t lowner_block {32};
 
    struct workspace_buffer
    {
      std::unique_ptr<double[]> data; 
      std::size_t               size; 
    };
 
    struct merge_workspace
    {
      workspace_buffer gaps    {}; 
      workspace_buffer updated {}; 
      workspace_buffer columns {}; 
      workspace_buffer factor  {}; 
    };
 
    inline constexpr std::size_t first_touch_chunk {std::size_t {1} << 15};
 
    [[nodiscard]] inline double* workspace_of(workspace_buffer& buffer,
                                              std::size_t       size)
    {
      if (buffer.size < size) {
        buffer.data = std::make_unique_for_overwrite<double[]>(size);
        buffer.size = size;
        double* const elements {buffer.data.get()};
        parallel_for((size + first_touch_chunk - 1) / first_touch_chunk, first_touch_chunk, [&](std::size_t chunk) {
                       const std::size_t first {chunk * first_touch_chunk};
                       std::fill(elements + first, elements + std::min(first + first_touch_chunk, size), 0.0);
                     });
      }
      return buffer.data.get();
    }
 
    inline constexpr std::size_t kept_workspace {std::size_t {1} << 18};
 
    [[nodiscard]] inline merge_workspace& thread_workspace()
    {
      thread_local merge_workspace work;
      return work;
    }
 
    inline void release_large_workspace(merge_workspace& work)
    {
      for (workspace_buffer* buffer : {&work.gaps, &work.updated, &work.columns, &work.factor}) {
        if (buffer->size > kept_workspace) {
          buffer->data.reset();
          buffer->size = 0;
        }
      }
    }
 
    inline void merge_rank_one(buffer<double>&      d,
                               buffer<double>       z,
                               double               rho,
                               double*              q,
                               std::size_t          ldq,
                               std::size_t          n,
                               std::size_t          split,
                               merge_workspace&     work)
    {
      // Which halves of the rows each column reaches: 1 the upper, 2 the lower.
      std::vector<unsigned> reach(n);
      for (std::size_t j {0}; j < n; ++j) {
        reach[j] = j < split ? 1U : 2U;
      }
      for (double& element : z) { // ‖z‖ = √2: the rows of two orthogonal matrices
        element /= std::numbers::sqrt2;
      }
      rho *= 2.0;
      std::vector<std::size_t> order(n);
      for (std::size_t j {0}; j < n; ++j) {
        order[j] = j;
      }
      std::sort(order.begin(), order.end(), [&d](std::size_t x, std::size_t y) { return d[x] < d[y]; });
      double largest_d {0.0};
      double largest_z {0.0};
      for (std::size_t j {0}; j < n; ++j) {
        largest_d = std::max(largest_d, std::abs(d[j]));
        largest_z = std::max(largest_z, std::abs(z[j]));
      }
      const double             tolerance {8.0 * std::numeric_limits<double>::epsilon() * std::max(largest_d, largest_z)};
      std::vector<std::size_t> kept;        // the columns the secular equation updates, by ascending d
      std::vector<std::size_t> deflated;    // the columns that keep their eigenvalue
      std::size_t              pending {n}; // the last column kept so far, still open to a rotation with the next
      for (const std::size_t j : order) {
        if (rho * std::abs(z[j]) <= tolerance) {
          deflated.push_back(j);
          continue;
        }
        if (pending == n) {
          pending = j;
          continue;
        }
        const double length {std::hypot(z[j], z[pending])};
        const double c      {z[j] / length};
        const double s      {-z[pending] / length};
        if (std::abs((d[j] - d[pending]) * c * s) <= tolerance) {
          // Close poles: a rotation of the two columns zeroes the pending one's weight, which then deflates.
          z[j]       = length;
          z[pending] = 0.0;
          for (std::size_t r {0}; r < n; ++r) {
            const double x {q[(r * ldq) + pending]};
            const double y {q[(r * ldq) + j]};
            q[(r * ldq) + pending] = (c * x) + (s * y);
            q[(r * ldq) + j]       = (c * y) - (s * x);
          }
          reach[j] = reach[pending] = reach[j] | reach[pending];
          const double rotated {(d[pending] * c * c) + (d[j] * s * s)};
          d[j]       = (d[pending] * s * s) + (d[j] * c * c);
          d[pending] = rotated;
          deflated.push_back(pending);
        }
        else {
          kept.push_back(pending);
        }
        pending = j;
      }
      if (pending != n) {
        kept.push_back(pending);
      }
      const std::size_t   k                 {kept.size()};
      buffer<double>      poles(k, double   {});
      buffer<double>      weights(k, double {});
      double              norm              {0.0};
      for (std::size_t i {0}; i < k; ++i) {
        poles[i]    = d[kept[i]];
        weights[i]  = z[kept[i]];
        norm       += weights[i] * weights[i];
      }
      norm = std::sqrt(norm);
      for (double& weight : weights) {
        weight /= norm;
      }
      const double        scaled_rho      {rho * norm * norm};
      buffer<double>      roots(k, double {});
      // gaps[i * k + j] = poles[j] − roots[i], every element written by its root's solve.
      double* const gaps                  {workspace_of(work.gaps, k * k)};
      parallel_for(k, 8 * k, [&](std::size_t i) {
                     roots[i] = secular_root(poles, weights, scaled_rho, i, secular_iterations, gaps + (i * k));
                   });
      // Löwner: ẑⱼ² = −Πᵢ (dⱼ − λᵢ) / Πᵢ≠ⱼ (dⱼ − dᵢ), each pole paired with its root, by blocks of poles, each pole's
      // product in ascending order of the roots (a row of gaps at a time, along the block): the same operations in
      // the same order on any number of threads.
      buffer<double> signed_weights(k, double {});
      parallel_for((k + lowner_block - 1) / lowner_block, lowner_block * k, [&](std::size_t block) {
                     const std::size_t first {block * lowner_block};
                     const std::size_t last {std::min(first + lowner_block, k)};
                     std::array<double, lowner_block> product {};
                     for (std::size_t j {first}; j < last; ++j) {
                       product[j - first] = gaps[(j * k) + j];
                     }
                     for (std::size_t i {0}; i < k; ++i) {
                       for (std::size_t j {first}; j < last; ++j) {
                         if (i != j) {
                           product[j - first] *= gaps[(i * k) + j] / (poles[j] - poles[i]);
                         }
                       }
                     }
                     for (std::size_t j {first}; j < last; ++j) {
                       signed_weights[j] = std::copysign(std::sqrt(-product[j - first]), weights[j]);
                     }
                   });
      // u[j][i] = ẑⱼ / (dⱼ − λᵢ), column i roots[i]'s eigenvector in the kept basis, in gaps' place: gaps[i][j] and
      // gaps[j][i] are read, then each written as the other's quotient, tile by tile across the diagonal.
      double* const     u     {gaps};
      const std::size_t tiles {(k + lowner_block - 1) / lowner_block};
      parallel_for(tiles * tiles, lowner_block * lowner_block, [&](std::size_t tile) {
                     const std::size_t ti {tile / tiles};
                     const std::size_t tj {tile % tiles};
                     if (ti > tj) {
                       return; // the pair of tiles is (tj, ti)'s
                     }
                     for (std::size_t i {ti * lowner_block}; i < std::min((ti + 1) * lowner_block, k); ++i) {
                       for (std::size_t j {std::max(tj * lowner_block, ti == tj ? i : 0)}; j < std::min((tj + 1) * lowner_block, k); ++j) {
                         const double above {gaps[(i * k) + j]};
                         const double below {gaps[(j * k) + i]};
                         u[(j * k) + i] = signed_weights[j] / above;
                         u[(i * k) + j] = signed_weights[i] / below;
                       }
                     }
                   });
      buffer<double> length(k, 0.0);
      for (std::size_t j {0}; j < k; ++j) { // by rows of u, each a contiguous update of every column's sum
        for (std::size_t i {0}; i < k; ++i) {
          length[i] += u[(j * k) + i] * u[(j * k) + i];
        }
      }
      for (double& l : length) {
        l = std::sqrt(l);
      }
      parallel_for(k, k, [&](std::size_t j) {
                     for (std::size_t i {0}; i < k; ++i) {
                       u[(j * k) + i] /= length[i];
                     }
                   });
      // updated = q(:, kept) · u, by halves of the rows: a half's rows take the kept columns that reach it, and
      // u's rows for those, the others being zero there.
      double* const updated {workspace_of(work.updated, n * k)};
      for (const auto& [half, first, rows] : {std::tuple {1U, std::size_t {0}, split}, std::tuple {2U, split, n - split}}) {
        std::vector<std::size_t> reaching;
        for (std::size_t i {0}; i < k; ++i) {
          if ((reach[kept[i]] & half) != 0U) {
            reaching.push_back(i);
          }
        }
        const std::size_t width   {reaching.size()};
        double* const     columns {workspace_of(work.columns, rows * width)};
        double* const     factor  {workspace_of(work.factor, width * k)};
        parallel_for(rows, width + k, [&](std::size_t r) {
                       for (std::size_t c {0}; c < width; ++c) {
                         columns[(r * width) + c] = q[((first + r) * ldq) + kept[reaching[c]]];
                       }
                       std::fill_n(updated + ((first + r) * k), k, 0.0); // the product adds to it
                     });
        parallel_for(width, k, [&](std::size_t c) { std::copy_n(u + (reaching[c] * k), k, factor + (c * k)); });
        gemm(rows, k, width, columns, width, factor, k, updated + (first * k), k);
      }
      // Every eigenpair, then by ascending eigenvalue: column j of the result is column source[j] of the
      // updated ones (below k) or deflated column source[j] − k of q, written row by row.
      buffer<double> values(n, double {});
      for (std::size_t i {0}; i < k; ++i) {
        values[i] = roots[i];
      }
      for (std::size_t i {0}; i < deflated.size(); ++i) {
        values[k + i] = d[deflated[i]];
      }
      for (std::size_t j {0}; j < n; ++j) {
        order[j] = j;
      }
      std::sort(order.begin(), order.end(), [&values](std::size_t x, std::size_t y) { return values[x] < values[y]; });
      for (std::size_t j {0}; j < n; ++j) {
        d[j] = values[order[j]];
      }
      parallel_for(n, n, [&](std::size_t r) {
                     const buffer<double> before(q + (r * ldq), q + (r * ldq) + n);
                     for (std::size_t j {0}; j < n; ++j) {
                       const std::size_t source {order[j]};
                       q[(r * ldq) + j] = source < k ? updated[(r * k) + source] : before[deflated[source - k]];
                     }
                   });
    }
 
    struct divide_node
    {
      double*        d        {};  
      double*        work_d   {};  
      const double*  work_e   {};  
      std::size_t    n        {};  
      double*        q        {};  
      std::size_t    ldq      {};  
      int            exponent {0}; 
      buffer<double> scaled_d {};  
      buffer<double> scaled_e {};  
      std::size_t    n1       {0}; 
      double         beta     {};  
      std::size_t    depth    {};  
      std::size_t    first    {};  
      std::size_t    second   {};  
    };
 
    NUMERICAL_OUT_OF_LINE inline void plan_divide(std::vector<divide_node>& nodes,
                                                  double*                   d,
                                                  const double*             e,
                                                  std::size_t               n,
                                                  double*                   q,
                                                  std::size_t               ldq,
                                                  std::size_t               depth)
    {
      double largest {0.0};
      for (std::size_t i {0}; i < n; ++i) {
        largest = std::max({largest, std::abs(d[i]), i + 1 < n ? std::abs(e[i]) : 0.0});
      }
      divide_node node {.d = d, .work_d = d, .work_e = e, .n = n, .q = q, .ldq = ldq, .depth = depth};
      (void)std::frexp(largest, &node.exponent);
      if (node.exponent != 0) {
        node.scaled_d.assign(n, double {});
        node.scaled_e.assign(n, 0.0);
        for (std::size_t i {0}; i < n; ++i) {
          node.scaled_d[i] = std::ldexp(d[i], -node.exponent);
          node.scaled_e[i] = i + 1 < n ? std::ldexp(e[i], -node.exponent) : 0.0;
        }
        node.work_d = node.scaled_d.data(); // the buffers' storage stays put as the node moves into the vector
        node.work_e = node.scaled_e.data();
      }
      if (n <= divide_leaf) {
        nodes.push_back(std::move(node));
        return;
      }
      node.n1                       = n / 2;
      node.beta                     = node.work_e[node.n1 - 1];
      node.work_d[node.n1 - 1]     -= std::abs(node.beta);
      node.work_d[node.n1]         -= std::abs(node.beta);
      double* const        work_d   {node.work_d};
      const double* const  work_e   {node.work_e};
      const std::size_t    n1       {node.n1};
      const std::size_t    index    {nodes.size()};
      nodes.push_back(std::move(node));
      nodes[index].first  = nodes.size();
      plan_divide(nodes, work_d, work_e, n1, q, ldq, depth + 1);
      nodes[index].second = nodes.size();
      plan_divide(nodes, work_d + n1, work_e + n1, n - n1, q + (n1 * ldq) + n1, ldq, depth + 1);
    }
 
    NUMERICAL_OUT_OF_LINE inline void solve_divide_leaf(const divide_node& node)
    {
      const std::size_t n {node.n};
      buffer<double>    diagonal(node.work_d, node.work_d + n);
      buffer<double>    sub(n, 0.0);
      std::copy(node.work_e, node.work_e + (n - 1), sub.begin());
      buffer<double> zt(n * n, 0.0);
      for (std::size_t i {0}; i < n; ++i) {
        zt[(i * n) + i] = 1.0;
      }
      tridiagonal_ql(diagonal, sub, zt, n, ql_steps_per_value * n);
      std::vector<std::size_t> order(n);
      for (std::size_t i {0}; i < n; ++i) {
        order[i] = i;
      }
      std::sort(order.begin(), order.end(), [&diagonal](std::size_t x, std::size_t y) { return diagonal[x] < diagonal[y]; });
      for (std::size_t j {0}; j < n; ++j) {
        node.work_d[j] = diagonal[order[j]];
        for (std::size_t r {0}; r < n; ++r) {
          node.q[(r * node.ldq) + j] = zt[(order[j] * n) + r];
        }
      }
    }
 
    NUMERICAL_OUT_OF_LINE inline void join_divide_node(const divide_node& node,
                                                       merge_workspace&   work)
    {
      const std::size_t n           {node.n};
      const std::size_t n1          {node.n1};
      const std::size_t n2          {n - n1};
      double* const     q           {node.q};
      buffer<double>    z(n, double {});
      for (std::size_t i {0}; i < n1; ++i) {
        z[i] = q[((n1 - 1) * node.ldq) + i];
      }
      for (std::size_t i {0}; i < n2; ++i) {
        z[n1 + i] = std::copysign(1.0, node.beta) * q[(n1 * node.ldq) + n1 + i];
      }
      buffer<double> values(node.work_d, node.work_d + n);
      merge_rank_one(values, std::move(z), std::abs(node.beta), q, node.ldq, n, n1, work);
      std::copy(values.begin(), values.end(), node.work_d);
    }
 
    NUMERICAL_OUT_OF_LINE inline void unscale_divide_node(const divide_node& node)
    {
      if (node.exponent != 0) {
        for (std::size_t i {0}; i < node.n; ++i) {
          node.d[i] = std::ldexp(node.work_d[i], node.exponent);
        }
      }
    }
 
    NUMERICAL_OUT_OF_LINE inline void solve_divide_subtree(std::vector<divide_node>& nodes,
                                                           std::size_t               index,
                                                           merge_workspace&          work)
    {
      const divide_node& node {nodes[index]};
      if (node.n1 == 0) {
        solve_divide_leaf(node);
      }
      else {
        solve_divide_subtree(nodes, node.first, work);
        solve_divide_subtree(nodes, node.second, work);
        join_divide_node(node, work);
      }
      unscale_divide_node(node);
    }
 
    inline void divide_and_conquer(double*              d,
                                   const double*        e,
                                   std::size_t          n,
                                   buffer<double>&      q)
    {
      q.assign(n * n, 0.0);
      std::vector<divide_node> nodes;
      plan_divide(nodes, d, e, n, q.data(), n, 0);
      const auto               cut {static_cast<std::size_t>(std::bit_width(max_threads() - 1))};
      std::vector<std::size_t> tasks;
      for (std::size_t i {0}; i < nodes.size(); ++i) {
        if (nodes[i].depth == cut || (nodes[i].depth < cut && nodes[i].n1 == 0)) {
          tasks.push_back(i);
        }
      }
      const std::size_t task_order {nodes[tasks[0]].n};
      parallel_for(tasks.size(), tasks.size() > 1 ? task_order * task_order * task_order : 0, [&](std::size_t t) {
                     solve_divide_subtree(nodes, tasks[t], thread_workspace());
                     release_large_workspace(thread_workspace());
                   });
      for (std::size_t level {cut}; level-- > 0;) {
        std::vector<std::size_t> joins;
        for (std::size_t i {0}; i < nodes.size(); ++i) {
          if (nodes[i].n1 != 0 && nodes[i].depth == level) {
            joins.push_back(i);
          }
        }
        if (joins.empty()) { // a tree shallower than the cut
          continue;
        }
        if (joins.size() == 1) {
          join_divide_node(nodes[joins[0]], thread_workspace());
          release_large_workspace(thread_workspace());
          unscale_divide_node(nodes[joins[0]]);
          continue;
        }
        const std::size_t order {nodes[joins[0]].n};
        parallel_for(joins.size(), order * order * order / 4, [&](std::size_t j) {
                       join_divide_node(nodes[joins[j]], thread_workspace());
                       release_large_workspace(thread_workspace());
                       unscale_divide_node(nodes[joins[j]]);
                     });
      }
    }
 
    template <typename T>
    void multiply_real_left(const buffer<double>&      left,
                            buffer<T>&                 zt,
                            std::size_t                n)
    {
      buffer<T> converted(n * n, T {});
      for (std::size_t i {0}; i < n * n; ++i) {
        converted[i] = from_double<T>(left[i]);
      }
      buffer<T> product(n * n, T {});
      gemm(n, n, n, converted.data(), n, zt.data(), n, product.data(), n);
      zt = std::move(product);
    }
 
    inline void multiply_real_left(const buffer<double>&      left,
                                   buffer<double>&            zt,
                                   std::size_t                n)
    {
      buffer<double> product(n * n, 0.0);
      gemm(n, n, n, left.data(), n, zt.data(), n, product.data(), n);
      zt = std::move(product);
    }
 
    inline void multiply_real_left(const buffer<double>&              left,
                                   buffer<std::complex<double>>&      zt,
                                   std::size_t                        n)
    {
      buffer<double> re(n * n, double {});
      buffer<double> im(n * n, double {});
      for (std::size_t i {0}; i < n * n; ++i) {
        re[i] = zt[i].real();
        im[i] = zt[i].imag();
      }
      multiply_real_left(left, re, n);
      multiply_real_left(left, im, n);
      for (std::size_t i {0}; i < n * n; ++i) {
        zt[i] = {re[i], im[i]};
      }
    }
 
    template <typename T>
    [[nodiscard]] buffer<T> tridiagonal_divide_and_conquer(buffer<double>&                 d,
                                                           buffer<double>&                 e,
                                                           const tridiagonal_reduction<T>& reduction,
                                                           std::size_t                     n)
    {
      if (n <= divide_leaf) {
        buffer<T> zt {tridiagonal_basis(reduction, n)};
        tridiagonal_ql(d, e, zt, n, ql_steps_per_value * n);
        return zt;
      }
      buffer<double> q;
      divide_and_conquer(d.data(), e.data(), n, q);
      buffer<T> at(n * n, T {}); // row j: the eigenvector of d[j], in the tridiagonal's basis, then the matrix's
      parallel_for(n, n, [&](std::size_t j) {
                     for (std::size_t r {0}; r < n; ++r) {
                       at[(j * n) + r] = from_double<T>(q[(r * n) + j]);
                     }
                   });
      to_matrix_basis(reduction, at, n);
      return at;
    }
 
    [[nodiscard]] inline double jacobi_off_diagonal(const buffer<double>&      a,
                                                    std::size_t                n)
    {
      double sum {0.0};
      for (std::size_t p {0}; p < n; ++p) {
        for (std::size_t q {p + 1}; q < n; ++q) {
          sum += a[(p * n) + q] * a[(p * n) + q];
        }
      }
      return std::sqrt(2.0 * sum);
    }
 
    inline void jacobi_rotate(buffer<double>&      a,
                              buffer<double>&      v,
                              std::size_t          n,
                              std::size_t          p,
                              std::size_t          q)
    {
      const double apq {a[(p * n) + q]};
      if (apq == 0.0) {
        return; // already zero — nothing to rotate
      }
      const double app   {a[(p * n) + p]};
      const double aqq   {a[(q * n) + q]};
      const double theta {(aqq - app) / (2.0 * apq)};
      const double t     {(theta >= 0.0 ? 1.0 : -1.0) / (std::abs(theta) + std::sqrt((theta * theta) + 1.0))};
      const double c     {1.0 / std::sqrt((t * t) + 1.0)};
      const double s     {t * c};
      a[(p * n) + p] = (c * c * app) - (2.0 * s * c * apq) + (s * s * aqq);
      a[(q * n) + q] = (s * s * app) + (2.0 * s * c * apq) + (c * c * aqq);
      a[(p * n) + q] = 0.0;
      a[(q * n) + p] = 0.0;
      for (std::size_t i {0}; i < n; ++i) {
        if (i != p && i != q) {
          const double aip {a[(i * n) + p]};
          const double aiq {a[(i * n) + q]};
          a[(i * n) + p] = (c * aip) - (s * aiq);
          a[(p * n) + i] = a[(i * n) + p];
          a[(i * n) + q] = (s * aip) + (c * aiq);
          a[(q * n) + i] = a[(i * n) + q];
        }
      }
      for (std::size_t i {0}; i < n; ++i) {
        const double vip {v[(i * n) + p]};
        const double viq {v[(i * n) + q]};
        v[(i * n) + p] = (c * vip) - (s * viq);
        v[(i * n) + q] = (s * vip) + (c * viq);
      }
    }
 
    inline void cyclic_jacobi(buffer<double>&      a,
                              buffer<double>&      v,
                              std::size_t          n,
                              double               frobenius,
                              std::size_t          max_sweeps)
    {
      // The off-diagonal norm strictly decreases each rotation and converges to 0,
      // so this terminates; the threshold is relative to the matrix scale.
      const double threshold {1e-12 * frobenius};
      for (std::size_t sweeps {0}; jacobi_off_diagonal(a, n) > threshold; ++sweeps) {
        if (sweeps == max_sweeps) {
          throw convergence_error("symmetric eigendecomposition: Jacobi did not converge within its sweeps");
        }
        for (std::size_t p {0}; p < n; ++p) {
          for (std::size_t q {p + 1}; q < n; ++q) {
            jacobi_rotate(a, v, n, p, q);
          }
        }
      }
    }
 
    [[nodiscard]] inline std::vector<std::size_t> ascending_order(const buffer<double>& diagonal)
    {
      std::vector<std::size_t> order(diagonal.size());
      for (std::size_t i {0}; i < order.size(); ++i) {
        order[i] = i;
      }
      std::sort(order.begin(), order.end(), [&diagonal](std::size_t x, std::size_t y) { return diagonal[x] < diagonal[y]; });
      return order;
    }
 
    template <typename T>
    [[nodiscard]] buffer<T> columns_by_value(const buffer<T>&                zt,
                                             const std::vector<std::size_t>& order,
                                             std::size_t                     n)
    {
      buffer<T> out(n * n, T {});
      parallel_for((n + transpose_tile - 1) / transpose_tile, transpose_tile * n, [&](std::size_t block) {
                     const std::size_t r0 {block * transpose_tile};
                     const std::size_t r1 {std::min(r0 + transpose_tile, n)};
                     for (std::size_t c0 {0}; c0 < n; c0 += transpose_tile) {
                       for (std::size_t c {c0}; c < std::min(c0 + transpose_tile, n); ++c) {
                         const T* const source {zt.data() + (order[c] * n)};
                         for (std::size_t r {r0}; r < r1; ++r) {
                           out[(r * n) + c] = source[r];
                         }
                       }
                     }
                   });
      return out;
    }
  } // namespace detail
 
  enum class eigen_method
  {
    tridiagonal_ql,    
    jacobi,            
    divide_and_conquer 
  };
 
  class symmetric_eigen
  {
  public:
 
    explicit symmetric_eigen(const matrix<double>& input,
                             eigen_method          method = eigen_method::divide_and_conquer)
      : n_(input.rows())
    {
      if (input.cols() != n_) {
        throw dimension_error("symmetric eigendecomposition requires a square matrix");
      }
      detail::require_finite_lower(input, "symmetric eigendecomposition");
      // Work on the matrix divided by its largest entry: a sum of squares of entries near 1e170
      // overflows and of entries near 1e-170 underflows to 0, and either stops the sweeps before a
      // single rotation. The eigenvalues are scaled back; the eigenvectors do not depend on it.
      double scale {0.0};
      for (std::size_t i {0}; i < n_; ++i) {
        for (std::size_t j {0}; j <= i; ++j) {
          scale = std::max(scale, std::abs(input(i, j)));
        }
      }
      if (scale == 0.0) {
        scale = 1.0;
      }
      buffer<double> a(n_ * n_, double   {});
      buffer<double> diagonal(n_, double {});
      if (method == eigen_method::jacobi) {
        buffer<double>      v(n_ * n_, 0.0);
        double              frobenius_squared {0.0};
        for (std::size_t i {0}; i < n_; ++i) {
          v[(i * n_) + i] = 1.0;
          for (std::size_t j {0}; j < n_; ++j) {
            a[(i * n_) + j]    = input(std::max(i, j), std::min(i, j)) / scale; // the lower triangle, mirrored
            frobenius_squared += a[(i * n_) + j] * a[(i * n_) + j];
          }
        }
        detail::cyclic_jacobi(a, v, n_, std::sqrt(frobenius_squared), detail::jacobi_sweeps);
        for (std::size_t i {0}; i < n_; ++i) {
          diagonal[i] = a[(i * n_) + i] * scale;
        }
        extract_sorted(diagonal, v);
        return;
      }
      for (std::size_t i {0}; i < n_; ++i) { // the lower triangle alone: the reduction reads nothing above it
        for (std::size_t j {0}; j <= i; ++j) {
          a[(i * n_) + j] = input(i, j) / scale;
        }
      }
      buffer<double> zt;
      if (n_ > 0) {
        buffer<double>                               off(n_, double {});
        const detail::tridiagonal_reduction<double>  reduction      {detail::tridiagonalize(a, n_, diagonal, off)};
        if (method == eigen_method::tridiagonal_ql) {
          zt = detail::tridiagonal_basis(reduction, n_);
          detail::tridiagonal_ql(diagonal, off, zt, n_, detail::ql_steps_per_value * n_);
        }
        else {
          zt = detail::tridiagonal_divide_and_conquer(diagonal, off, reduction, n_);
        }
      }
      const std::vector<std::size_t> order             {detail::ascending_order(diagonal)};
      buffer<double>                 values(n_, double {});
      for (std::size_t c {0}; c < n_; ++c) {
        values[c] = diagonal[order[c]] * scale;
      }
      values_  = vector<double> {std::move(values)};
      vectors_ = matrix<double> {n_, n_, detail::columns_by_value(zt, order, n_)};
    }
 
    [[nodiscard]] const vector<double>& values() const noexcept
    {
      return values_;
    }
 
    [[nodiscard]] const matrix<double>& vectors() const& noexcept
    {
      return vectors_;
    }
 
    [[nodiscard]] matrix<double> vectors() &&
    {
      return std::move(vectors_);
    }
 
  private:
 
    void extract_sorted(const buffer<double>& diagonal,
                        const buffer<double>& v)
    {
      std::vector<std::size_t> order(n_);
      for (std::size_t i {0}; i < n_; ++i) {
        order[i] = i;
      }
      std::sort(order.begin(), order.end(),
                [&diagonal](std::size_t x, std::size_t y) {
                  return diagonal[x] < diagonal[y];
                });
      buffer<double> values(n_, double       {});
      buffer<double> vectors(n_ * n_, double {});
      for (std::size_t col {0}; col < n_; ++col) {
        values[col] = diagonal[order[col]];
        for (std::size_t row {0}; row < n_; ++row) {
          vectors[(row * n_) + col] = v[(row * n_) + order[col]];
        }
      }
      values_  = vector<double> {std::move(values)};
      vectors_ = matrix<double> {n_, n_, std::move(vectors)};
    }
 
    std::size_t    n_;       
    vector<double> values_;  
    matrix<double> vectors_; 
  };
 
  class hermitian_eigen
  {
  public:
 
    explicit hermitian_eigen(const matrix<std::complex<double>>& input,
                             eigen_method                        method = eigen_method::divide_and_conquer)
      : n_(input.rows())
    {
      if (method == eigen_method::jacobi) {
        throw std::invalid_argument("Hermitian eigendecomposition: the Jacobi method is for real symmetric matrices");
      }
      if (input.cols() != n_) {
        throw dimension_error("Hermitian eigendecomposition requires a square matrix");
      }
      detail::require_finite_lower(input, "Hermitian eigendecomposition");
      // Work on the matrix divided by its largest entry, as symmetric_eigen does; the vectors do not depend
      // on it. The entries read: the lower triangle, and the diagonal's real part (its imaginary part is ignored,
      // so it scales nothing: an infinite one made every value NaN, a huge one wrong).
      double scale {0.0};
      for (std::size_t i {0}; i < n_; ++i) {
        for (std::size_t j {0}; j < i; ++j) {
          scale = std::max(scale, std::abs(input(i, j)));
        }
        scale = std::max(scale, std::abs(input(i, i).real()));
      }
      if (scale == 0.0) {
        scale = 1.0;
      }
      buffer<std::complex<double>> a(n_ * n_, std::complex<double> {});
      for (std::size_t i {0}; i < n_; ++i) { // the lower triangle alone, a real diagonal: the reduction reads no more
        for (std::size_t j {0}; j <= i; ++j) {
          a[(i * n_) + j] = (i == j ? std::complex<double> {input(i, j).real(), 0.0} : input(i, j)) / scale;
        }
      }
      buffer<double>               diagonal(n_, double {});
      buffer<std::complex<double>> zt;
      if (n_ > 0) {
        buffer<double>                                             off(n_, double {});
        const detail::tridiagonal_reduction<std::complex<double>>  reduction      {detail::tridiagonalize(a, n_, diagonal, off)};
        if (method == eigen_method::tridiagonal_ql) {
          zt = detail::tridiagonal_basis(reduction, n_);
          detail::tridiagonal_ql(diagonal, off, zt, n_, detail::ql_steps_per_value * n_);
        }
        else {
          zt = detail::tridiagonal_divide_and_conquer(diagonal, off, reduction, n_);
        }
      }
      const std::vector<std::size_t> order             {detail::ascending_order(diagonal)};
      buffer<double>                 values(n_, double {});
      for (std::size_t c {0}; c < n_; ++c) {
        values[c] = diagonal[order[c]] * scale;
      }
      values_  = vector<double> {std::move(values)};
      vectors_ = matrix<std::complex<double>> {n_, n_, detail::columns_by_value(zt, order, n_)};
    }
 
    [[nodiscard]] const vector<double>& values() const noexcept
    {
      return values_;
    }
 
    [[nodiscard]] const matrix<std::complex<double>>& vectors() const& noexcept
    {
      return vectors_;
    }
 
    [[nodiscard]] matrix<std::complex<double>> vectors() &&
    {
      return std::move(vectors_);
    }
 
  private:
 
    std::size_t                  n_;       
    vector<double>               values_;  
    matrix<std::complex<double>> vectors_; 
  };
 
  namespace detail {
 
    inline constexpr double jacobi_tiny {std::numeric_limits<double>::min() / std::numeric_limits<double>::epsilon()};
 
    inline constexpr double jacobi_large {1e150};
 
    template <typename T>
    void one_sided_jacobi(buffer<T>&      w,
                          buffer<T>&      v,
                          std::size_t     n,
                          std::size_t     max_sweeps)
    {
      // The sweeps rotate columns: work on the transposes, where each column is a contiguous row. A walk
      // down the columns of the row-major arrays strides by n on every access and runs several times slower.
      const auto transpose_square {[n](buffer<T>& m) {
                                     for (std::size_t i {0}; i < n; ++i) {
                                       for (std::size_t j {i + 1}; j < n; ++j) {
                                         std::swap(m[(i * n) + j], m[(j * n) + i]);
                                       }
                                     }
                                   }};
      transpose_square(w);
      transpose_square(v);
      // Orthogonal to the rounding of the columns' storage: a float column holds no inner product to double's.
      const double epsilon {static_cast<double>(std::numeric_limits<value_range_t<T>>::epsilon())};
      for (std::size_t sweep {1};; ++sweep) {
        bool rotated {false};
        for (std::size_t p {0}; p < n; ++p) {
          for (std::size_t q {p + 1}; q < n; ++q) {
            const T* wp    {w.data() + (p * n)};
            const T* wq    {w.data() + (q * n)};
            double   alpha {0.0};
            double   beta  {0.0};
            T        gamma {};
            for (std::size_t i {0}; i < n; ++i) {
              alpha += static_cast<double>(std::norm(wp[i]));
              beta  += static_cast<double>(std::norm(wq[i]));
              gamma += conj_of(wp[i]) * wq[i];
            }
            double magnitude {static_cast<double>(std::abs(gamma))};
            double zeta      {};
            if (std::min(alpha, beta) < jacobi_tiny) {
              // A column whose squares underflow (rows of a far smaller scale than the largest): its sum of squares
              // is 0 or has lost its digits, and the test below never holds. The pair is judged on the columns over
              // their norms, which need no squares: the cosine, and ζ from the ratio of the norms.
              const double norm_p {scaled_norm(wp, n)};
              const double norm_q {scaled_norm(wq, n)};
              if (norm_p == 0.0 || norm_q == 0.0) {
                continue; // a zero column is orthogonal to any other
              }
              T cosine {};
              for (std::size_t i {0}; i < n; ++i) {
                cosine += conj_of(wp[i] / detail::from_double<T>(norm_p)) * (wq[i] / detail::from_double<T>(norm_q));
              }
              const double size {static_cast<double>(std::abs(cosine))};
              zeta = ((norm_q / norm_p) - (norm_p / norm_q)) / (2.0 * size);
              if (size <= epsilon) {
                continue; // orthogonal to rounding
              }
              gamma     = cosine; // the phase is γ's
              magnitude = size;
            }
            else {
              if (magnitude <= epsilon * std::sqrt(alpha) * std::sqrt(beta)) {
                continue; // already orthogonal to rounding (or a zero column)
              }
              zeta = (beta - alpha) / (2.0 * magnitude);
            }
            // Past jacobi_large, ζ² overflows and t would be 0 (no rotation, ever): there |ζ| + √(1 + ζ²) is 2|ζ| exactly,
            // and t is 1 / (2|ζ|). A tangent that still underflows to 0 is a rotation below the smallest double.
            const double t {(zeta >= 0.0 ? 1.0 : -1.0) * (std::abs(zeta) > jacobi_large ? 0.5 / std::abs(zeta) : 1.0 / (std::abs(zeta) + std::sqrt(1.0 + (zeta * zeta))))};
            if (t == 0.0) {
              continue; // a rotation below the smallest double: the pair is as orthogonal as it can be
            }
            rotated = true;
            // Turn column q by the phase of γ so that the pair's inner product is real, then rotate.
            const T      turn {conj_of(gamma) / detail::from_double<T>(magnitude)};
            const double c    {1.0 / std::sqrt(1.0 + (t * t))};
            const double s    {c * t};
            for (buffer<T>* side : {&w, &v}) {
              T* mp {side->data() + (p * n)};
              T* mq {side->data() + (q * n)};
              for (std::size_t i {0}; i < n; ++i) {
                const T left  {mp[i]};
                const T right {turn * mq[i]};
                mp[i] = (detail::from_double<T>(c) * left) - (detail::from_double<T>(s) * right);
                mq[i] = (detail::from_double<T>(s) * left) + (detail::from_double<T>(c) * right);
              }
            }
          }
        }
        if (!rotated) {
          transpose_square(w);
          transpose_square(v);
          return;
        }
        if (sweep == max_sweeps) {
          throw convergence_error("SVD: one-sided Jacobi did not converge within its sweeps");
        }
      }
    }
 
    inline constexpr std::size_t svd_qr_steps_per_value {30};
 
    inline constexpr std::size_t bidiagonal_crossover {128};
 
    inline constexpr std::size_t bidiagonal_basis_crossover {64};
 
    template <typename T>
    void form_bidiagonal_bases(const std::vector<buffer<T>>&      left,
                               std::vector<buffer<T>>             right,
                               const buffer<T>&                   left_phase,
                               const buffer<T>&                   right_phase,
                               std::size_t                        n,
                               buffer<T>&                         ut,
                               buffer<T>&                         vt)
    {
      // U_b = L₀ᴴ L₁ᴴ … with L_kᴴ = H_k diag(conj phase at k), and V_b = R₀ R₁ … with R_k = H'_k diag(phase
      // at k+1). A phase on row k commutes with every later reflector, which acts on rows past it, so
      // U_b = (H₀ H₁ …) diag(conj phases) and V_b = (H'₀ H'₁ …) diag(phases): the reflectors are applied to the
      // diagonal by blocks, through the matrix product kernel, as QR's Q is formed.
      ut.assign(n * n, T {}); // Uᵀ and Vᵀ, formed transposed: diagonal to start, the reflectors applied to them
      vt.assign(n * n, T {});
      for (std::size_t k {0}; k < n; ++k) {
        ut[(k * n) + k] = conj_of(left_phase[k]);
        vt[(k * n) + k] = k == 0 ? T {1} : right_phase[k - 1];
      }
      std::vector<buffer<T>> shifted(n); // right[k] acts from row k+1: as reflector k+1 of a QR
      for (std::size_t k {0}; k + 1 < n; ++k) {
        shifted[k + 1] = std::move(right[k]);
      }
      if (n < bidiagonal_basis_crossover) { // a reflector at a time, the last first, as LAPACK's dorg2r
        for (std::size_t k {n}; k-- > 0;) {
          reflect_columns(ut.data(), n, left[k], k, k, n);
          reflect_columns(vt.data(), n, shifted[k], k, k, n);
        }
        return;
      }
      for (std::size_t k1 {n}; k1 > 0;) { // blocks of qr_block reflectors, the last first
        const std::size_t k0 {(k1 - 1) / qr_block * qr_block};
        apply_transposed_reflector_block(gather_reflectors(left, k0, k1, n), ut, n, k0, n);
        apply_transposed_reflector_block(gather_reflectors(shifted, k0, k1, n), vt, n, k0, n);
        k1 = k0;
      }
    }
 
#ifdef NUMERICAL_X86_COMPILED
    inline constexpr std::size_t bidiagonal_shared_least {256};
#endif
 
    inline constexpr std::size_t bidiagonal_block {32};
 
#ifdef NUMERICAL_X86_COMPILED
    inline constexpr std::size_t bidiagonal_product_weight {4};
#else
    inline constexpr std::size_t bidiagonal_product_weight {1};
#endif
 
    template <typename T>
    void bidiagonalize_unblocked(buffer<T>&           b,
                                 std::size_t          n,
                                 buffer<double>&      d,
                                 buffer<double>&      e,
                                 buffer<T>&           ut,
                                 buffer<T>&           vt,
                                 bool                 form = true);
 
    template <typename T>
    [[nodiscard]] T householder_in_place(buffer<T>& x)
    {
      const double up     {rescale_tiny(x)}; // the reflector is the same at any scale; β is scaled back
      const auto   norm   {static_cast<T>(scaled_norm(x))};
      const T      beta   {-phase_of(x[0]) * norm};
      x[0] -= beta;
      const auto length   {static_cast<T>(scaled_norm(x))};
      for (T& element : x) {
        element /= length;
      }
      return beta / static_cast<T>(up);
    }
 
    template <typename T>
    requires std::is_floating_point_v<T>
    void bidiagonalize(buffer<T>&           b,
                       std::size_t          n,
                       buffer<double>&      d,
                       buffer<double>&      e,
                       buffer<T>&           ut,
                       buffer<T>&           vt,
                       bool                 form = true)
    {
      if (n < bidiagonal_crossover) {
        bidiagonalize_unblocked(b, n, d, e, ut, vt, form);
        return;
      }
      std::vector<buffer<T>>      left(n);  // the left reflector of step k, on rows k…, empty where none was needed
      std::vector<buffer<T>>      right(n); // the right reflector of step k, on columns k+1…, likewise
      const auto                  nonzero {[](const T& value) { return value != T {}; }};
      // By reflector, each row over the whole order: row 2p of vx holds step k0 + p's left reflector v (at its rows),
      // row 2p + 1 what its right one moves, x = 2 B u (at the rows); row 2p of yu what v moves, y = 2 Bᵀ v (at the
      // columns), row 2p + 1 the right reflector u (at its columns). A step's terms of the panel's earlier ones are then
      // one update through their rows (subtract_products, a family's rows 2n apart), in their order.
      buffer<T> vx(2 * bidiagonal_block * n);
      buffer<T> yu(2 * bidiagonal_block * n);
      buffer<T> columns(bidiagonal_block * n);      // the panel's columns of B, gathered row by row
      buffer<T> factors(2 * bidiagonal_block);      // a step's coefficients of those rows
#ifdef NUMERICAL_X86_COMPILED
      buffer<T> next_factors(2 * bidiagonal_block); // the next column's, applied ahead (a shared step)
#endif
      for (std::size_t k0 {0}; k0 < n; k0 += bidiagonal_block) {
        const std::size_t k1    {std::min(k0 + bidiagonal_block, n)};
        const std::size_t width {k1 - k0};
        std::fill_n(vx.data(), 2 * width * n, T {});
        std::fill_n(yu.data(), 2 * width * n, T {});
        for (std::size_t r {k0}; r < n; ++r) { // column k from row k down, each row's part read once
          for (std::size_t k {k0}; k < std::min(k1, r + 1); ++k) {
            columns[((k - k0) * n) + r] = b[(r * n) + k];
          }
        }
        const auto v_row {[&vx, n](std::size_t p, std::size_t from) { return vx.data() + (2 * p * n) + from; }};
        const auto x_row {[&vx, n](std::size_t p, std::size_t from) { return vx.data() + (((2 * p) + 1) * n) + from; }};
        const auto y_row {[&yu, n](std::size_t p, std::size_t from) { return yu.data() + (2 * p * n) + from; }};
        const auto u_row {[&yu, n](std::size_t p, std::size_t from) { return yu.data() + (((2 * p) + 1) * n) + from; }};
#ifdef NUMERICAL_X86_COMPILED
        bool column_ahead {false}; // whether the step's column took the steps before the last ahead (a shared step)
#endif
        for (std::size_t k {k0}; k < k1; ++k) {
          const std::size_t i    {k - k0};
#ifdef NUMERICAL_X86_COMPILED
          const std::size_t from {column_ahead ? i - 1 : 0};
          column_ahead = false;
#else
          constexpr std::size_t from {0};
#endif
          for (std::size_t p {from}; p < i; ++p) { // v's row by y's element, x's by u's
            factors[2 * (p - from)]       = *y_row(p, k);
            factors[(2 * (p - from)) + 1] = *u_row(p, k);
          }
          T* const column {columns.data() + (i * n) + k}; // column k from the diagonal down, the panel's earlier steps applied
          subtract_products(column, factors.data(), v_row(from, k), n, 2 * (i - from), n - k);
          buffer<T> x(column, column + (n - k));
          if (std::any_of(x.begin() + 1, x.end(), nonzero)) { // something below the diagonal to reflect
            d[k] = static_cast<double>(householder_in_place(x));
            // y = 2 (B − V Yᵀ − X Uᵀ)ᵀ v on the columns past k: Bᵀ v by contiguous rows, then the panel's terms.
            buffer<T> by_vx(2 * i, T {});                     // Vᵀ v and Xᵀ v over the earlier steps, interleaved as their rows
            for (std::size_t p {0}; p < i; ++p) {
              by_vx[2 * p]       = dot_product(v_row(p, k), x.data(), n - k);
              by_vx[(2 * p) + 1] = dot_product(x_row(p, k), x.data(), n - k);
            }
            const std::size_t     span        {n - k - 1};
            constexpr std::size_t chunk_width {64};
            // Bᵀ v as the terms −(−v)·B, the chunk held through them (negation is exact, so each column takes the same
            // additions): row by row, each term stored the chunk back (singular values of order 512 on one x86-64
            // core: 46 → 42 ms).
            buffer<T> minus(x.size(), T {});
            std::transform(x.begin(), x.end(), minus.begin(), std::negate<T> {});
            parallel_for((span + chunk_width - 1) / chunk_width, bidiagonal_product_weight * (n - k) * chunk_width, [&](std::size_t chunk) {
                           const std::size_t first {k + 1 + (chunk * chunk_width)};
                           const std::size_t count {std::min(chunk_width, n - first)};
                           buffer<T>    sum(count, T {});
                           subtract_products(sum.data(), minus.data(), b.data() + (k * n) + first, n, n - k, count);
                           subtract_products(sum.data(), by_vx.data(), y_row(0, first), n, 2 * i, count);
                           for (std::size_t c {0}; c < count; ++c) {
                             *y_row(i, first + c) = static_cast<T>(2) * sum[c];
                           }
                         });
            std::copy(x.begin(), x.end(), v_row(i, k));
            left[k] = std::move(x);
          }
          else {
            d[k] = static_cast<double>(x[0]);
          }
          if (k + 1 == n) {
            break;
          }
          buffer<T> y(b.begin() + static_cast<std::ptrdiff_t>((k * n) + k + 1), // row k past the diagonal, the left
                      b.begin() + static_cast<std::ptrdiff_t>((k + 1) * n));    // reflection included
          for (std::size_t p {0}; p <= i; ++p) {
            factors[p] = *v_row(p, k);
          }
          subtract_products(y.data(), factors.data(), y_row(0, k + 1), 2 * n, i + 1, n - k - 1);
          for (std::size_t p {0}; p < i; ++p) {
            factors[p] = *x_row(p, k);
          }
          subtract_products(y.data(), factors.data(), u_row(0, k + 1), 2 * n, i, n - k - 1);
          if (std::any_of(y.begin() + 1, y.end(), nonzero)) { // something past the superdiagonal to reflect
            e[k] = static_cast<double>(householder_in_place(y));
            // x = 2 (B − V Yᵀ − X Uᵀ) u on the rows past k: B u by rows, then the panel's terms.
            const std::size_t rows          {n - k - 1};
            buffer<T>         by_y(i + 1, T {}); // Yᵀ u over this and the earlier steps, Uᵀ u over the earlier
            buffer<T>         by_u(i, T     {});
            buffer<T>         moved(rows, T {});
#ifdef NUMERICAL_X86_COMPILED
            if (rows >= bidiagonal_shared_least) {
              // The coefficients run among B u's tasks, and the corrections by bands on the threads, the next column
              // taking the earlier steps' terms in the same pass (its rows are these): each element the operations,
              // in the order, of the branch above.
              parallel_for(rows + (2 * i) + 1, bidiagonal_product_weight * 2 * (n - k), [&](std::size_t task) {
                             if (task < rows) {
                               moved[task] = dot_product(b.data() + ((k + 1 + task) * n) + k + 1, y.data(), y.size());
                             }
                             else if (task - rows <= i) {
                               by_y[task - rows] = dot_product(y_row(task - rows, k + 1), y.data(), y.size());
                             }
                             else {
                               by_u[task - rows - i - 1] = dot_product(u_row(task - rows - i - 1, k + 1), y.data(), y.size());
                             }
                           });
              T* const next {k + 1 < k1 && i > 0 ? columns.data() + ((i + 1) * n) + k + 1 : nullptr};
              for (std::size_t p {0}; next != nullptr && p < i; ++p) { // the next column's terms of the steps before this
                next_factors[2 * p]       = *y_row(p, k + 1);
                next_factors[(2 * p) + 1] = *u_row(p, k + 1);
              }
              constexpr std::size_t band {256};
              parallel_for((rows + band - 1) / band, band * 4 * (i + 1) * 8, [&](std::size_t c) {
                             const std::size_t first {c * band};
                             const std::size_t count {std::min(band, rows - first)};
                             subtract_products(moved.data() + first, by_y.data(), v_row(0, k + 1 + first), 2 * n, i + 1, count);
                             subtract_products(moved.data() + first, by_u.data(), x_row(0, k + 1 + first), 2 * n, i, count);
                             for (std::size_t r {first}; r < first + count; ++r) {
                               *x_row(i, k + 1 + r) = static_cast<T>(2) * moved[r];
                             }
                             if (next != nullptr) {
                               subtract_products(next + first, next_factors.data(), v_row(0, k + 1 + first), n, 2 * i, count);
                             }
                           });
              column_ahead = next != nullptr;
            }
            else
#endif
            {
              for (std::size_t p {0}; p <= i; ++p) {
                by_y[p] = dot_product(y_row(p, k + 1), y.data(), y.size());
              }
              for (std::size_t p {0}; p < i; ++p) {
                by_u[p] = dot_product(u_row(p, k + 1), y.data(), y.size());
              }
              parallel_for(rows, bidiagonal_product_weight * 2 * (n - k), [&](std::size_t row) {
                             moved[row] = dot_product(b.data() + ((k + 1 + row) * n) + k + 1, y.data(), y.size());
                           });
              subtract_products(moved.data(), by_y.data(), v_row(0, k + 1), 2 * n, i + 1, rows);
              subtract_products(moved.data(), by_u.data(), x_row(0, k + 1), 2 * n, i, rows);
              for (std::size_t r {0}; r < rows; ++r) {
                *x_row(i, k + 1 + r) = static_cast<T>(2) * moved[r];
              }
            }
 
            std::copy(y.begin(), y.end(), u_row(i, k + 1));
            right[k] = std::move(y);
          }
          else {
            e[k] = static_cast<double>(y[0]);
          }
        }
        if (k1 < n) { // the trailing block from row and column k1: B −= V Yᵀ + X Uᵀ
          const std::size_t rest                    {n - k1};
          buffer<T>         minus_v(rest * width, T {});
          buffer<T>         minus_x(rest * width, T {});
          for (std::size_t r {0}; r < rest; ++r) {
            for (std::size_t p {0}; p < width; ++p) {
              minus_v[(r * width) + p] = -*v_row(p, k1 + r);
              minus_x[(r * width) + p] = -*x_row(p, k1 + r);
            }
          }
          T* trailing {b.data() + (k1 * n) + k1};
          gemm(rest, rest, width, minus_v.data(), width, y_row(0, k1), 2 * n, trailing, n);
          gemm(rest, rest, width, minus_x.data(), width, u_row(0, k1), 2 * n, trailing, n);
        }
      }
      e[n - 1] = 0.0;
      if (form) {
        form_bidiagonal_bases(left, std::move(right), buffer<T>(n, T {1}), buffer<T>(n, T {1}), n, ut, vt);
      }
    }
 
    template <typename T>
    void bidiagonalize_unblocked(buffer<T>&           b,
                                 std::size_t          n,
                                 buffer<double>&      d,
                                 buffer<double>&      e,
                                 buffer<T>&           ut,
                                 buffer<T>&           vt,
                                 bool                 form)
    {
      std::vector<buffer<T>>      left(n);  // the left reflector of step k, empty where none was needed
      std::vector<buffer<T>>      right(n); // the right reflector acting on columns k+1…, likewise
      buffer<T>                   left_phase(n, T {1});
      buffer<T>                   right_phase(n, T {1});
      const auto                  reflector {[](buffer<T>& v) {       // unit v with (I − 2 v vᴴ) x ∥ e₁
                                               (void)rescale_tiny(v); // the same reflector at any scale
                                               const double norm {scaled_norm(v)};
                                               v[0] += phase_of(v[0]) * detail::from_double<T>(norm);
                                               const auto length {detail::from_double<T>(scaled_norm(v))};
                                               for (T & element : v) {
                                                 element /= length;
                                               }
                                             }};
      for (std::size_t k {0}; k < n; ++k) {
        buffer<T> x(n - k, T {}); // column k from the diagonal down
        for (std::size_t i {k}; i < n; ++i) {
          x[i - k] = b[(i * n) + k];
        }
        if (scaled_norm(x.data() + 1, x.size() - 1) > 0.0) {
          reflector(x);
          reflect_rows(b, n, x, k, k);
          left[k] = std::move(x);
        }
        left_phase[k] = conj_of(phase_of(b[(k * n) + k])); // row k times it puts |b_kk| on the diagonal
        for (std::size_t j {k}; j < n; ++j) {
          b[(k * n) + j] *= left_phase[k];
        }
        d[k] = static_cast<double>(std::abs(b[(k * n) + k]));
        if (k + 1 == n) {
          break;
        }
        buffer<T> y(n - k - 1, T {}); // row k right of the diagonal, conjugated: B ← B H with H yᴴ ∥ e₁
        for (std::size_t j {k + 1}; j < n; ++j) {
          y[j - k - 1] = conj_of(b[(k * n) + j]);
        }
        if (y.size() > 1 && scaled_norm(y.data() + 1, y.size() - 1) > 0.0) {
          reflector(y);
          parallel_for(n - k, 2 * y.size(), [&](std::size_t r) { // row ← row (I − 2 y yᴴ); rows above k are zero there
                         T*      row   {b.data() + ((k + r) * n) + k + 1};
                         const T scale {static_cast<T>(2) * dot_product(row, y.data(), y.size())};
                         for (std::size_t j {0}; j < y.size(); ++j) {
                           row[j] -= scale * conj_of(y[j]);
                         }
                       });
          right[k] = std::move(y);
        }
        right_phase[k] = conj_of(phase_of(b[(k * n) + k + 1])); // column k+1 times it: |b_k,k+1| above
        for (std::size_t i {k}; i < n; ++i) {
          b[(i * n) + k + 1] *= right_phase[k];
        }
        e[k] = static_cast<double>(std::abs(b[(k * n) + k + 1]));
      }
      if (n > 0) {
        e[n - 1] = 0.0;
      }
      if (form) {
        form_bidiagonal_bases(left, std::move(right), left_phase, right_phase, n, ut, vt);
      }
    }
 
    template <typename T>
    requires(!std::is_floating_point_v<T>)
    void bidiagonalize(buffer<T>&           b,
                       std::size_t          n,
                       buffer<double>&      d,
                       buffer<double>&      e,
                       buffer<T>&           ut,
                       buffer<T>&           vt,
                       bool                 form = true)
    {
      bidiagonalize_unblocked(b, n, d, e, ut, vt, form);
    }
 
    template <typename T>
    void bidiagonal_qr(buffer<double>&      d,
                       buffer<double>&      e,
                       buffer<T>&           ut,
                       buffer<T>&           vt,
                       std::size_t          n,
                       std::size_t          max_steps)
    {
      // The iteration reads only d and e: what it does to the rows of vt and ut is recorded and applied once
      // at the end, in parallel by columns.
      std::vector<row_operation> on_vt;
      std::vector<row_operation> on_ut;
      const bool                 record  {!vt.empty()}; // the singular values alone record nothing
      const double               epsilon {std::numeric_limits<double>::epsilon()};
      const double               tiny    {std::numeric_limits<double>::min() / epsilon};
      constexpr std::size_t      none    {std::numeric_limits<std::size_t>::max()};
      std::size_t                p       {n};
      std::size_t                steps   {0};
      while (p > 0) {
        // The largest k < p−1 with a negligible e[k] (none: the loop ran out and wrapped).
        std::size_t k {p - 1};
        while (k-- > 0) {
          if (std::abs(e[k]) <= tiny + (epsilon * (std::abs(d[k]) + std::abs(d[k + 1])))) {
            e[k] = 0.0;
            break;
          }
        }
        const std::size_t top {k == none ? 0 : k + 1}; // the unreduced block is top…p−1
        if (top == p - 1) { // d[p−1] has converged: make it non-negative, sort it into place
          if (d[top] < 0.0) {
            d[top] = -d[top];
            if (record) {
              on_vt.push_back({row_operation::kind::negate, top, top, 0.0, 0.0});
            }
          }
          for (std::size_t j {top}; j + 1 < n && d[j] < d[j + 1]; ++j) {
            std::swap(d[j], d[j + 1]);
            if (record) {
              on_vt.push_back({row_operation::kind::swap, j, j + 1, 0.0, 0.0});
              on_ut.push_back({row_operation::kind::swap, j, j + 1, 0.0, 0.0});
            }
          }
          --p;
          continue;
        }
        // The last negligible d[zero] in the block (none: there is none).
        std::size_t zero {none};
        for (std::size_t ks {p}; ks-- > top;) {
          const double around {std::abs(e[ks]) + (ks != top ? std::abs(e[ks - 1]) : 0.0)};
          if (std::abs(d[ks]) <= tiny + (epsilon * around)) {
            d[ks] = 0.0;
            zero  = ks;
            break;
          }
        }
        if (zero == p - 1) { // chase e[p−2] up the block with rotations on the right
          double f {e[p - 2]};
          e[p - 2] = 0.0;
          for (std::size_t j {p - 1}; j-- > top;) {
            const double t  {rotation_length(d[j], f)};
            const double cs {d[j] / t};
            const double sn {f / t};
            d[j] = t;
            if (j != top) {
              f        = -sn * e[j - 1];
              e[j - 1] = cs * e[j - 1];
            }
            if (record) {
              on_vt.push_back({row_operation::kind::rotate, j, p - 1, cs, sn});
            }
          }
          continue;
        }
        if (zero != none) { // split at d[zero] = 0: chase e[zero] along the row with rotations on the left
          double f {e[zero]};
          e[zero] = 0.0;
          for (std::size_t j {zero + 1}; j < p; ++j) {
            const double t  {rotation_length(d[j], f)};
            const double cs {d[j] / t};
            const double sn {f / t};
            d[j] = t;
            f    = -sn * e[j];
            e[j] = cs * e[j];
            if (record) {
              on_ut.push_back({row_operation::kind::rotate, j, zero, cs, sn});
            }
          }
          continue;
        }
        if (steps == max_steps) {
          throw convergence_error("SVD: the bidiagonal QR iteration did not converge within its steps");
        }
        ++steps;
        // One implicit QR step on top…p−1, shifted by the eigenvalue of the trailing 2 × 2 of BᵀB nearer
        // its last diagonal entry.
        const double scale {std::max({std::abs(d[p - 1]), std::abs(d[p - 2]), std::abs(e[p - 2]), std::abs(d[top]),
                                      std::abs(e[top])})};
        const double sp    {d[p - 1] / scale};
        const double spm1  {d[p - 2] / scale};
        const double epm1  {e[p - 2] / scale};
        const double sk    {d[top] / scale};
        const double ek    {e[top] / scale};
        const double bb    {(((spm1 + sp) * (spm1 - sp)) + (epm1 * epm1)) / 2.0};
        const double cc    {(sp * epm1) * (sp * epm1)};
        double       shift {0.0};
        if (bb != 0.0 || cc != 0.0) {
          shift = std::sqrt((bb * bb) + cc);
          shift = cc / (bb < 0.0 ? bb - shift : bb + shift);
        }
        double f {((sk + sp) * (sk - sp)) + shift};
        double g {sk * ek};
        for (std::size_t j {top}; j + 1 < p; ++j) {
          double t  {rotation_length(f, g)};
          double cs {f / t};
          double sn {g / t};
          if (j != top) {
            e[j - 1] = t;
          }
          f        = (cs * d[j]) + (sn * e[j]);
          e[j]     = (cs * e[j]) - (sn * d[j]);
          g        = sn * d[j + 1];
          d[j + 1] = cs * d[j + 1];
          if (record) {
            on_vt.push_back({row_operation::kind::rotate, j, j + 1, cs, sn});
          }
          t        = rotation_length(f, g);
          cs       = f / t;
          sn       = g / t;
          d[j]     = t;
          f        = (cs * e[j]) + (sn * d[j + 1]);
          d[j + 1] = (cs * d[j + 1]) - (sn * e[j]);
          g        = sn * e[j + 1];
          e[j + 1] = cs * e[j + 1];
          if (record) {
            on_ut.push_back({row_operation::kind::rotate, j, j + 1, cs, sn});
          }
        }
        e[p - 2] = f;
      }
      if (!vt.empty()) {
        apply_row_operations(vt, n, on_vt);
        apply_row_operations(ut, n, on_ut);
      }
    }
 
    template <typename T>
    [[nodiscard]] matrix<T> adjoint(const matrix<T>& a)
    {
      buffer<T> entries(a.rows() * a.cols(), T {});
      for (std::size_t i {0}; i < a.rows(); ++i) {
        for (std::size_t j {0}; j < a.cols(); ++j) {
          entries[(j * a.rows()) + i] = conj_of(a(i, j));
        }
      }
      return matrix<T> {a.cols(), a.rows(), std::move(entries)};
    }
 
    [[nodiscard]] inline double singular_secular_root(const buffer<double>&      d,
                                                      const buffer<double>&      z,
                                                      std::size_t                i,
                                                      std::size_t                max_iterations,
                                                      double*                    below,
                                                      double*                    above)
    {
      buffer<double>         delta(d.size(), double {});
      const secular_solution root                   {solve_secular(z, 1.0, i, max_iterations, delta.data(), [&d](std::size_t o, std::size_t j) {
                                                                     return (d[j] - d[o]) * (d[j] + d[o]);
                                                                   })};
      const double           pole                   {d[root.origin]};
      const double           shift                  {root.offset / (pole + std::sqrt((pole * pole) + root.offset))}; // σ − d_o from σ² − d_o²
      for (std::size_t j {0}; j < d.size(); ++j) {
        below[j] = (d[j] - pole) - shift;
        above[j] = (d[j] + pole) + shift;
      }
      return pole + shift;
    }
 
    inline void rotate_columns(buffer<double>&      m,
                               std::size_t          cols,
                               std::size_t          a,
                               std::size_t          b,
                               double               c,
                               double               s)
    {
      for (std::size_t r {0}; r < m.size() / cols; ++r) {
        const double x {m[(r * cols) + a]};
        const double y {m[(r * cols) + b]};
        m[(r * cols) + a] = (c * x) + (s * y);
        m[(r * cols) + b] = (c * y) - (s * x);
      }
    }
 
    struct arrow_workspace
    {
      workspace_buffer below    {}; 
      workspace_buffer above    {}; 
      workspace_buffer left     {}; 
      workspace_buffer right    {}; 
      workspace_buffer new_u    {}; 
      workspace_buffer new_v    {}; 
      workspace_buffer columns  {}; 
      workspace_buffer factor   {}; 
      buffer<double>   sorted_u {}; 
      buffer<double>   sorted_v {}; 
    };
 
    [[nodiscard]] inline arrow_workspace& thread_arrow_workspace()
    {
      thread_local arrow_workspace work;
      return work;
    }
 
    inline void release_large_arrow_workspace(arrow_workspace& work)
    {
      for (workspace_buffer* buffer : {&work.below, &work.above, &work.left, &work.right, &work.new_u, &work.new_v,
                                       &work.columns, &work.factor}) {
        if (buffer->size > kept_workspace) {
          buffer->data.reset();
          buffer->size = 0;
        }
      }
      for (buffer<double>* sorted : {&work.sorted_u, &work.sorted_v}) {
        if (sorted->capacity() > kept_workspace) {
          *sorted = buffer<double> {};
        }
      }
    }
 
    inline void merge_arrow(buffer<double>&      d,
                            buffer<double>       z,
                            double               alpha,
                            double               beta,
                            buffer<double>&      u,
                            buffer<double>&      v,
                            std::size_t          n,
                            std::size_t          m,
                            std::size_t          split,
                            arrow_workspace&     work)
    {
      // The rows each column reaches: 1 the upper half's, 2 u's row `split` (the arrow's first), 4 the lower half's.
      // A rotation of two columns makes each reach the rows of both.
      std::vector<unsigned> reach_u(n, 4U);
      std::vector<unsigned> reach_v(m, 4U);
      reach_u[0] = 2U;
      reach_v[0] = m > n ? 5U : 1U; // the upper half's null vector, turned with the lower one's when there is one
      for (std::size_t j {1}; j <= split; ++j) {
        reach_u[j] = 1U;
        reach_v[j] = 1U;
      }
      const auto turned {[](std::vector<unsigned>& reach, std::size_t a, std::size_t b) {
                           reach[a] = reach[b] = reach[a] | reach[b];
                         }};
      std::vector<std::size_t> order(n - 1);
      for (std::size_t j {1}; j < n; ++j) {
        order[j - 1] = j;
      }
      std::sort(order.begin(), order.end(), [&d](std::size_t x, std::size_t y) { return d[x] < d[y]; });
      double largest {std::max(std::abs(alpha), std::abs(beta))};
      for (const double value : d) {
        largest = std::max(largest, value);
      }
      if (largest == 0.0) {
        return; // a zero arrow: its singular values are zeros, its vectors the bases as they are
      }
      const double tolerance {8.0 * std::numeric_limits<double>::epsilon() * largest};
      if (std::abs(z[0]) <= tolerance) {
        z[0] = tolerance; // the zero pole keeps a weight: the arrow moves by no more than the tolerance
      }
      std::vector<std::size_t> kept    {0}; // by ascending d, the zero pole first
      std::vector<std::size_t> deflated;    // the columns that keep their singular value
      std::size_t              pending {n}; // the last column kept so far, still open to a rotation with the next
      for (const std::size_t j : order) {
        if (std::abs(z[j]) <= tolerance) {
          deflated.push_back(j);
          continue;
        }
        if (d[j] <= tolerance) {
          // Nearly the zero pole: its weight joins the zero pole's, and dⱼ deflates.
          const double length {std::hypot(z[0], z[j])};
          rotate_columns(v, m, 0, j, z[0] / length, z[j] / length);
          turned(reach_v, 0, j);
          z[0] = length;
          z[j] = 0.0;
          deflated.push_back(j);
          continue;
        }
        if (pending != n && d[j] - d[pending] <= tolerance) {
          // Close values: one rotation of both rows and columns zeroes the pending one's weight.
          const double length {std::hypot(z[j], z[pending])};
          const double c      {z[j] / length};
          const double s      {z[pending] / length};
          rotate_columns(v, m, j, pending, c, s);
          rotate_columns(u, n, j, pending, c, s);
          turned(reach_v, j, pending);
          turned(reach_u, j, pending);
          z[j]       = length;
          z[pending] = 0.0;
          deflated.push_back(pending);
        }
        else if (pending != n) {
          kept.push_back(pending);
        }
        pending = j;
      }
      if (pending != n) {
        kept.push_back(pending);
      }
      const std::size_t   k                 {kept.size()};
      buffer<double>      poles(k, double   {});
      buffer<double>      weights(k, double {});
      for (std::size_t i {0}; i < k; ++i) {
        poles[i]   = d[kept[i]];
        weights[i] = z[kept[i]];
      }
      buffer<double> roots(k, double     {});
      double* const  below               {workspace_of(work.below, k * k)}; // below[i * k + j] = poles[j] − roots[i]
      double* const  above               {workspace_of(work.above, k * k)}; // above[i * k + j] = poles[j] + roots[i]
      parallel_for(k, 8 * k, [&](std::size_t i) {
                     roots[i] = singular_secular_root(poles, weights, i, secular_iterations, below + (i * k), above + (i * k));
                   });
      double* const left  {workspace_of(work.left, k * k)};  // column i: the left vector of roots[i], in the kept columns' basis
      double* const right {workspace_of(work.right, k * k)}; // and the right one
      // Löwner: ẑⱼ² = Πᵢ (σᵢ² − dⱼ²) / Πᵢ≠ⱼ (dᵢ² − dⱼ²), each pole paired with its root, by blocks of poles on every
      // core, each pole's product in ascending order of the roots (a row of the gaps at a time, along the block).
      buffer<double> signed_weights(k, double {});
      parallel_for((k + lowner_block - 1) / lowner_block, lowner_block * k, [&](std::size_t block) {
                     const std::size_t                first {block * lowner_block};
                     const std::size_t                last {std::min(first + lowner_block, k)};
                     std::array<double, lowner_block> product {};
                     for (std::size_t j {first}; j < last; ++j) {
                       product[j - first] = -below[(j * k) + j] * above[(j * k) + j];
                     }
                     for (std::size_t i {0}; i < k; ++i) {
                       for (std::size_t j {first}; j < last; ++j) {
                         if (i != j) {
                           product[j - first] *= (-below[(i * k) + j] * above[(i * k) + j]) / ((poles[i] - poles[j]) * (poles[i] + poles[j]));
                         }
                       }
                     }
                     for (std::size_t j {first}; j < last; ++j) {
                       signed_weights[j] = std::copysign(std::sqrt(product[j - first]), weights[j]);
                     }
                   });
      // The vectors by rows of poles, each on its own core: right[j][i] = ẑⱼ / (dⱼ² − σᵢ²), read a column of the gaps.
      parallel_for(k, 4 * k, [&](std::size_t j) {
                     for (std::size_t i {0}; i < k; ++i) {
                       const double ratio {signed_weights[j] / (below[(i * k) + j] * above[(i * k) + j])}; // ẑⱼ / (dⱼ² − σᵢ²)
                       right[(j * k) + i] = ratio;
                       left[(j * k) + i]  = j == 0 ? -1.0 : poles[j] * ratio;
                     }
                   });
      for (double* const vectors : {left, right}) {
        buffer<double> length(k, 0.0);
        for (std::size_t j {0}; j < k; ++j) { // by rows, each a contiguous update of every column's sum
          for (std::size_t i {0}; i < k; ++i) {
            length[i] += vectors[(j * k) + i] * vectors[(j * k) + i];
          }
        }
        for (double& l : length) {
          l = std::sqrt(l);
        }
        parallel_for(k, k, [&](std::size_t j) {
                       for (std::size_t i {0}; i < k; ++i) {
                         vectors[(j * k) + i] /= length[i];
                       }
                     });
      }
      // The kept columns times the arrow's vectors; the deflated columns as they are; then by ascending value.
      // By segments of rows, each taking the kept columns that reach it and their rows of the arrow's vectors: the
      // others are zero there, and an element's sum without them is the same to the bit.
      using segment = std::tuple<std::size_t, std::size_t, unsigned>; // first row, rows, the bit of reach
      const auto carried {[&](const buffer<double>& basis, std::size_t rows, const double* small,
                              const std::vector<unsigned>& reach, const std::vector<segment>& segments, double* product) {
                            for (const auto& [first, count, bit] : segments) {
                              std::vector<std::size_t> reaching;
                              for (std::size_t i {0}; i < k; ++i) {
                                if ((reach[kept[i]] & bit) != 0U) {
                                  reaching.push_back(i);
                                }
                              }
                              const std::size_t   width {reaching.size()};
                              double* const       columns {workspace_of(work.columns, count * width)};
                              double* const       factor {workspace_of(work.factor, width * k)};
                              parallel_for(count, width + k, [&](std::size_t r) {
                                             for (std::size_t c {0}; c < width; ++c) {
                                               columns[(r * width) + c] = basis[((first + r) * rows) + kept[reaching[c]]];
                                             }
                                             std::fill_n(product + ((first + r) * k), k, 0.0); // the product adds to it
                                           });
                              for (std::size_t c {0}; c < width; ++c) {
                                std::copy_n(small + (reaching[c] * k), k, factor + (c * k));
                              }
                              gemm(count, k, width, columns, width, factor, k, product + (first * k), k);
                            }
                          }};
      double* const new_u {workspace_of(work.new_u, n * k)};
      double* const new_v {workspace_of(work.new_v, m * k)};
      carried(u, n, left, reach_u, {{0, split, 1U}, {split, 1, 2U}, {split + 1, n - split - 1, 4U}}, new_u);
      carried(v, m, right, reach_v, {{0, split + 1, 1U}, {split + 1, m - split - 1, 4U}}, new_v);
      buffer<double>            values(n, double {});
      std::vector<std::size_t>  from(n); // where each value's vectors are: a kept root (< k) or a deflated column
      for (std::size_t i {0}; i < k; ++i) {
        values[i] = roots[i];
        from[i]   = i;
      }
      for (std::size_t i {0}; i < deflated.size(); ++i) {
        values[k + i] = d[deflated[i]];
        from[k + i]   = k + i;
      }
      std::vector<std::size_t> rank(n);
      for (std::size_t j {0}; j < n; ++j) {
        rank[j] = j;
      }
      std::sort(rank.begin(), rank.end(), [&values](std::size_t x, std::size_t y) { return values[x] < values[y]; });
      std::vector<std::size_t> source(n); // by sorted column: a kept root (< k) or a deflated column (k + its place)
      for (std::size_t c {0}; c < n; ++c) {
        source[c] = from[rank[c]];
        d[c]      = values[rank[c]];
      }
      // Row by row on every core, each row written whole (by columns, each write was a row apart).
      buffer<double>& sorted_u {work.sorted_u};
      buffer<double>& sorted_v {work.sorted_v};
      sorted_u.resize(n * n);              // written whole below
      sorted_v.assign(v.begin(), v.end()); // the null column, if any, stays last
      parallel_for(m, n, [&](std::size_t r) { // m ≥ n: the right basis has the rows of both
                     for (std::size_t c {0}; c < n; ++c) {
                       if (r < n) {
                         sorted_u[(r * n) + c] = source[c] < k ? new_u[(r * k) + source[c]] : u[(r * n) + deflated[source[c] - k]];
                       }
                       sorted_v[(r * m) + c] = source[c] < k ? new_v[(r * k) + source[c]] : v[(r * m) + deflated[source[c] - k]];
                     }
                   });
      std::swap(u, sorted_u); // the old bases' storage is kept for the next join
      std::swap(v, sorted_v);
    }
 
    inline constexpr std::size_t bidiagonal_leaf {32};
 
    struct bidiagonal_node
    {
      const double*  d        {};  
      const double*  e        {};  
      std::size_t    n        {};  
      std::size_t    extra    {};  
      int            exponent {0}; 
      buffer<double> scaled_d {};  
      buffer<double> scaled_e {};  
      std::size_t    k        {0}; 
      std::size_t    first    {};  
      std::size_t    second   {};  
      std::size_t    depth    {};  
      buffer<double> sigma    {};  
      buffer<double> u        {};  
      buffer<double> v        {};  
    };
 
    NUMERICAL_OUT_OF_LINE inline std::size_t plan_bidiagonal(std::vector<bidiagonal_node>& nodes,
                                                             const double*                 d,
                                                             const double*                 e,
                                                             std::size_t                   n,
                                                             std::size_t                   extra,
                                                             std::size_t                   depth)
    {
      double largest {0.0};
      for (std::size_t i {0}; i < n; ++i) {
        largest = std::max({largest, std::abs(d[i]), i + 1 < n + extra ? std::abs(e[i]) : 0.0});
      }
      bidiagonal_node node {.d = d, .e = e, .n = n, .extra = extra, .depth = depth};
      (void)std::frexp(largest, &node.exponent);
      if (node.exponent != 0) {
        node.scaled_d.assign(n, double {});
        node.scaled_e.assign(n, 0.0);
        for (std::size_t i {0}; i < n; ++i) {
          node.scaled_d[i] = std::ldexp(d[i], -node.exponent);
          node.scaled_e[i] = i + 1 < n + extra ? std::ldexp(e[i], -node.exponent) : 0.0;
        }
        node.d = node.scaled_d.data(); // the buffers' storage stays put as the node moves into the vector
        node.e = node.scaled_e.data();
      }
      const std::size_t index {nodes.size()};
      if (n <= bidiagonal_leaf) {
        nodes.push_back(std::move(node));
        return index;
      }
      node.k = n / 2;
      const double* const work_d {node.d};
      const double* const work_e {node.e};
      const std::size_t   k      {node.k};
      nodes.push_back(std::move(node));
      const std::size_t first    {plan_bidiagonal(nodes, work_d, work_e, k, 1, depth + 1)};
      const std::size_t second   {plan_bidiagonal(nodes, work_d + k + 1, work_e + k + 1, n - k - 1, extra, depth + 1)};
      nodes[index].first  = first;
      nodes[index].second = second;
      return index;
    }
 
    NUMERICAL_OUT_OF_LINE inline void start_bidiagonal_node(bidiagonal_node& node)
    {
      const std::size_t m {node.n + node.extra};
      node.sigma.assign(node.n, 0.0);
      node.u.assign(node.n * node.n, 0.0);
      node.v.assign(m * m, 0.0);
      for (std::size_t i {0}; i < m; ++i) {
        node.v[(i * m) + i] = 1.0;
      }
    }
 
    NUMERICAL_OUT_OF_LINE inline void solve_bidiagonal_leaf(bidiagonal_node& node)
    {
      start_bidiagonal_node(node);
      const double* const d     {node.d};
      const double* const e     {node.e};
      const std::size_t   n     {node.n};
      const std::size_t   extra {node.extra};
      const std::size_t   m     {n + extra};
      buffer<double>&     sigma {node.sigma};
      buffer<double>&     u     {node.u};
      buffer<double>&     v     {node.v};
      buffer<double>      diagonal(d, d + n);
      buffer<double>      super(n, 0.0);
      std::copy_n(e, n - std::min<std::size_t>(n, 1), super.begin());
      // An extra column first: rotations of it against each column, from the last row up, chase its one entry
      // out (B G = [B' 0]), so the square B' remains.
      double fill {extra == 1 ? e[n - 1] : 0.0};
      for (std::size_t j {n}; extra == 1 && j-- > 0;) {
        const double length {std::hypot(diagonal[j], fill)};
        const double c      {length == 0.0 ? 1.0 : diagonal[j] / length};
        const double s      {length == 0.0 ? 0.0 : fill / length};
        diagonal[j] = length;
        rotate_columns(v, m, j, n, c, s);
        if (j > 0) {
          fill         = -s * super[j - 1];
          super[j - 1] = c * super[j - 1];
        }
      }
      buffer<double> ut(n * n, 0.0);
      buffer<double> vt(n * n, 0.0);
      for (std::size_t i {0}; i < n; ++i) {
        ut[(i * n) + i] = 1.0;
        vt[(i * n) + i] = 1.0;
      }
      bidiagonal_qr(diagonal, super, ut, vt, n, svd_qr_steps_per_value * n);   // descending
      buffer<double> rotated(v);
      for (std::size_t c {0}; c < n; ++c) {
        const std::size_t source {n - 1 - c};
        sigma[c] = diagonal[source];
        for (std::size_t r {0}; r < n; ++r) {
          u[(r * n) + c] = ut[(source * n) + r];
        }
        for (std::size_t r {0}; r < m; ++r) {   // G times B''s right vector
          double sum {0.0};
          for (std::size_t t {0}; t < n; ++t) {
            sum += rotated[(r * m) + t] * vt[(source * n) + t];
          }
          v[(r * m) + c] = sum;
        }
      }
    }
 
    NUMERICAL_OUT_OF_LINE inline void join_bidiagonal_node(bidiagonal_node&       node,
                                                           const bidiagonal_node& one,
                                                           const bidiagonal_node& two,
                                                           arrow_workspace&       work)
    {
      const std::size_t     n      {node.n};
      const std::size_t     extra  {node.extra};
      const std::size_t     m      {n + extra};
      const std::size_t     k      {node.k};
      const std::size_t     lower  {n - k - 1};
      const double          alpha  {node.d[k]};
      const double          beta   {node.e[k]};
      const buffer<double>& sigma1 {one.sigma};
      const buffer<double>& u1     {one.u};
      const buffer<double>& v1     {one.v};
      const buffer<double>& sigma2 {two.sigma};
      const buffer<double>& u2     {two.u};
      const buffer<double>& v2     {two.v};
      buffer<double>&       u      {node.u};
      buffer<double>&       v      {node.v};
      const std::size_t     m1     {k + 1};
      const std::size_t     m2     {lower + extra};
      buffer<double>        arrow(n, 0.0);
      buffer<double>        z(n, 0.0);
      // The bases: row k is the arrow's first row; the arrow's columns are the null vector, then the halves'.
      u[k * n] = 1.0;
      for (std::size_t c {0}; c < k; ++c) {
        arrow[1 + c] = sigma1[c];
        z[1 + c]     = alpha * v1[(k * m1) + c];
        for (std::size_t r {0}; r < k; ++r) {
          u[(r * n) + 1 + c] = u1[(r * k) + c];
        }
      }
      for (std::size_t c {0}; c < lower; ++c) {
        arrow[k + 1 + c] = sigma2[c];
        z[k + 1 + c]     = beta * v2[c];
        for (std::size_t r {0}; r < lower; ++r) {
          u[((k + 1 + r) * n) + k + 1 + c] = u2[(r * lower) + c];
        }
      }
      std::fill(v.begin(), v.end(), 0.0);
      for (std::size_t r {0}; r < m1; ++r) {
        v[(r * m) + 0] = v1[(r * m1) + k];
        for (std::size_t c {0}; c < k; ++c) {
          v[(r * m) + 1 + c] = v1[(r * m1) + c];
        }
      }
      for (std::size_t r {0}; r < m2; ++r) {
        for (std::size_t c {0}; c < lower; ++c) {
          v[((m1 + r) * m) + k + 1 + c] = v2[(r * m2) + c];
        }
      }
      z[0] = alpha * v1[(k * m1) + k];
      if (extra == 1) { // the two null vectors: one joins the arrow, the other is the whole's
        for (std::size_t r {0}; r < m2; ++r) {
          v[((m1 + r) * m) + n] = v2[(r * m2) + lower];
        }
        const double other  {beta * v2[lower]};
        const double length {std::hypot(z[0], other)};
        const double c      {length == 0.0 ? 1.0 : z[0] / length};
        const double s      {length == 0.0 ? 0.0 : other / length};
        rotate_columns(v, m, 0, n, c, s);
        z[0] = length;
      }
      merge_arrow(arrow, std::move(z), alpha, beta, u, v, n, m, k, work);
      node.sigma = std::move(arrow);
    }
 
    NUMERICAL_OUT_OF_LINE inline void unscale_bidiagonal_node(bidiagonal_node& node)
    {
      if (node.exponent != 0) {
        for (double& value : node.sigma) {
          value = std::ldexp(value, node.exponent);
        }
      }
    }
 
    NUMERICAL_OUT_OF_LINE inline void join_bidiagonal_halves(std::vector<bidiagonal_node>& nodes,
                                                             std::size_t                   index,
                                                             arrow_workspace&              work)
    {
      bidiagonal_node& node {nodes[index]};
      join_bidiagonal_node(node, nodes[node.first], nodes[node.second], work);
      unscale_bidiagonal_node(node);
      nodes[node.first]  = bidiagonal_node {}; // read by no other join
      nodes[node.second] = bidiagonal_node {};
    }
 
    NUMERICAL_OUT_OF_LINE inline void solve_bidiagonal_subtree(std::vector<bidiagonal_node>& nodes,
                                                               std::size_t                   index,
                                                               arrow_workspace&              work)
    {
      if (nodes[index].k == 0) {
        solve_bidiagonal_leaf(nodes[index]);
        unscale_bidiagonal_node(nodes[index]);
        return;
      }
      start_bidiagonal_node(nodes[index]); // before its halves, as the recursion did: their buffers then come back to it
      solve_bidiagonal_subtree(nodes, nodes[index].first, work);
      solve_bidiagonal_subtree(nodes, nodes[index].second, work);
      join_bidiagonal_halves(nodes, index, work);
    }
 
    inline void bidiagonal_divide(const double*        d,
                                  const double*        e,
                                  std::size_t          n,
                                  std::size_t          extra,
                                  buffer<double>&      sigma,
                                  buffer<double>&      u,
                                  buffer<double>&      v)
    {
      std::vector<bidiagonal_node> nodes;
      (void)plan_bidiagonal(nodes, d, e, n, extra, 0);
      const auto               cut {static_cast<std::size_t>(std::bit_width(max_threads() - 1))};
      std::vector<std::size_t> tasks;
      for (std::size_t i {0}; i < nodes.size(); ++i) {
        if (nodes[i].depth == cut || (nodes[i].depth < cut && nodes[i].k == 0)) {
          tasks.push_back(i);
        }
      }
      const std::size_t task_order {nodes[tasks[0]].n};
      parallel_for(tasks.size(), tasks.size() > 1 ? task_order * task_order * task_order : 0, [&](std::size_t t) {
                     solve_bidiagonal_subtree(nodes, tasks[t], thread_arrow_workspace());
                     release_large_arrow_workspace(thread_arrow_workspace());
                   });
      for (std::size_t level {cut}; level-- > 0;) {
        std::vector<std::size_t> joins;
        for (std::size_t i {0}; i < nodes.size(); ++i) {
          if (nodes[i].k != 0 && nodes[i].depth == level) {
            joins.push_back(i);
          }
        }
        const std::size_t order {joins.empty() ? 0 : nodes[joins[0]].n};
        parallel_for(joins.size(), joins.size() > 1 ? order * order * order / 4 : 0, [&](std::size_t j) {
                       start_bidiagonal_node(nodes[joins[j]]);
                       join_bidiagonal_halves(nodes, joins[j], thread_arrow_workspace());
                       release_large_arrow_workspace(thread_arrow_workspace());
                     });
      }
      sigma = std::move(nodes[0].sigma);
      u     = std::move(nodes[0].u);
      v     = std::move(nodes[0].v);
    }
 
    template <typename T>
    void bidiagonal_divide_and_conquer(buffer<double>&      d,
                                       buffer<double>&      e,
                                       buffer<T>&           ut,
                                       buffer<T>&           vt,
                                       std::size_t          n)
    {
      if (n == 0) {
        return;
      }
      buffer<double> sigma;
      buffer<double> u;
      buffer<double> v;
      bidiagonal_divide(d.data(), e.data(), n, 0, sigma, u, v);
      buffer<double> left(n * n, double  {}); // row c: the vectors of the c-th largest value
      buffer<double> right(n * n, double {});
      for (std::size_t c {0}; c < n; ++c) {
        d[c] = sigma[n - 1 - c];
        for (std::size_t r {0}; r < n; ++r) {
          left[(c * n) + r]  = u[(r * n) + (n - 1 - c)];
          right[(c * n) + r] = v[(r * n) + (n - 1 - c)];
        }
      }
      multiply_real_left(left, ut, n);
      multiply_real_left(right, vt, n);
    }
  } // namespace detail
 
  enum class svd_method
  {
    bidiagonal_qr,     
    jacobi,            
    divide_and_conquer 
  };
 
  template <typename T>
  class svd_decomposition
  {
  public:
 
    explicit svd_decomposition(const matrix<T>& a,
                               svd_method       method = svd_method::divide_and_conquer,
                               bool             full   = false)
    {
      detail::require_finite(a, "SVD");
      if (a.rows() >= a.cols()) {
        factor(a, method, full);
      }
      else { // A = (Aᴴ)ᴴ = (U' Σ V'ᴴ)ᴴ = V' Σ U'ᴴ
        factor(detail::adjoint(a), method, full);
        std::swap(u_, v_);
      }
    }
 
    [[nodiscard]] const vector<double>& singular_values() const noexcept
    {
      return singular_values_;
    }
 
    [[nodiscard]] const matrix<T>& u() const& noexcept
    {
      return u_;
    }
 
    [[nodiscard]] matrix<T> u() &&
    {
      return std::move(u_);
    }
 
    [[nodiscard]] const matrix<T>& v() const& noexcept
    {
      return v_;
    }
 
    [[nodiscard]] matrix<T> vh() const
    {
      return detail::adjoint(v_);
    }
 
  private:
 
    [[nodiscard]] static matrix<T> left_vectors(buffer<T>   q,
                                                buffer<T>   left,
                                                std::size_t rows,
                                                std::size_t cols,
                                                std::size_t width)
    {
      if (width == cols) {
        return matrix<T> {rows, cols, std::move(q)} *matrix<T> {cols, cols, std::move(left)};
      }
      buffer<T> u(rows * rows, T {});
      detail::gemm(rows, cols, cols, q.data(), rows, left.data(), cols, u.data(), rows);
      for (std::size_t i {0}; i < rows; ++i) {
        for (std::size_t j {cols}; j < rows; ++j) {
          u[(i * rows) + j] = q[(i * rows) + j];
        }
      }
      return matrix<T> {rows, rows, std::move(u)};
    }
 
    void factor(const matrix<T>& a,
                svd_method       method,
                bool             full)
    {
      const std::size_t rows {a.rows()};
      const std::size_t cols {a.cols()};
      // Work on A divided by its largest magnitude: the sums of squares in the sweeps overflow near
      // 1e155 and underflow near 1e-155. The singular values are scaled back; the vectors do not change.
      // Both passes by bands of rows on the threads, bound by the memory: a band's largest magnitude, then the
      // largest of those (exact in any order), then the quotients.
      const T* const      source {std::to_address(a.begin())};
      const std::size_t   bands  {(rows + detail::svd_band - 1) / detail::svd_band};
      std::vector<double> largest(bands, 0.0);
      detail::parallel_for(bands, detail::memory_element_weight * detail::svd_band * cols, [source, rows, cols, &largest](std::size_t band) {
                             double most {0.0};
                             for (std::size_t e {band * detail::svd_band * cols}; e < std::min(rows, (band + 1) * detail::svd_band) * cols; ++e) {
                               most = std::max(most, static_cast<double>(std::abs(source[e])));
                             }
                             largest[band] = most;
                           });
      double scale {*std::max_element(largest.begin(), largest.end())};
      if (scale == 0.0) {
        scale = 1.0;
      }
      buffer<T>         entries(rows * cols);
      T* const          into    {entries.data()};
      const T           divisor {detail::from_double<T>(scale)};
      detail::parallel_for(bands, detail::memory_element_weight * detail::svd_band * cols, [source, into, rows, cols, divisor](std::size_t band) {
                             for (std::size_t e {band * detail::svd_band * cols}; e < std::min(rows, (band + 1) * detail::svd_band) * cols; ++e) {
                               into[e] = source[e] / divisor;
                             }
                           });
      // A square matrix is reduced to bidiagonal form as it is (as LAPACK's dgesdd does short of a tall one): a QR
      // first, its Q formed and multiplied in at the end, took a third of the decomposition's time at order 1024.
      const bool         direct  {rows == cols && method != svd_method::jacobi};
      const std::size_t  width   {full ? rows : cols};
      detail::thin_qr<T> factors {direct ? detail::thin_qr<T> {.q = {}, .r = std::move(entries)}
                                         : detail::householder_qr(std::move(entries), rows, cols, full)};
      if (method != svd_method::jacobi) {
        buffer<double> d(cols, double {});
        buffer<double> e(cols, double {});
        buffer<T>      ut;
        buffer<T>      vt;
        detail::bidiagonalize(factors.r, cols, d, e, ut, vt);
        if (method == svd_method::bidiagonal_qr) {
          detail::bidiagonal_qr(d, e, ut, vt, cols, detail::svd_qr_steps_per_value * cols);
        }
        else {
          detail::bidiagonal_divide_and_conquer(d, e, ut, vt, cols);
        }
        buffer<T> left(cols * cols);  // rows of ut and vt are the vectors: made columns
        buffer<T> right(cols * cols);
        detail::transpose_into(ut.data(), cols, cols, left.data());
        detail::transpose_into(vt.data(), cols, cols, right.data());
        for (double& value : d) {
          value *= scale;
        }
        singular_values_ = vector<double> {std::move(d)};
        u_               = direct ? matrix<T> {cols, cols, std::move(left)}
                                  : left_vectors(std::move(factors.q), std::move(left), rows, cols, width);
        v_               = matrix<T> {cols, cols, std::move(right)};
        return;
      }
      buffer<T>     w                    {std::move(factors.r)}; // A = Q R, and R V = W with orthogonal columns
      buffer<T>     right(cols * cols, T {});
      for (std::size_t i {0}; i < cols; ++i) {
        right[(i * cols) + i] = T {1};
      }
      detail::one_sided_jacobi(w, right, cols, detail::jacobi_sweeps);
 
      buffer<double> norms(cols, double {});
      for (std::size_t j {0}; j < cols; ++j) {
        buffer<T> column(cols, T {});
        for (std::size_t i {0}; i < cols; ++i) {
          column[i] = w[(i * cols) + j];
        }
        norms[j] = detail::scaled_norm(column);
      }
      std::vector<std::size_t> order(cols);
      for (std::size_t j {0}; j < cols; ++j) {
        order[j] = j;
      }
      std::stable_sort(order.begin(), order.end(), [&norms](std::size_t x, std::size_t y) { return norms[x] > norms[y]; });
 
      // The left vectors of R: each nonzero column of W over its norm; a zero one has none, and takes
      // the unit vector that leaves the largest residual against the columns already set, made
      // orthogonal to them.
      buffer<T>      left(cols * cols, T         {});
      buffer<T>      sorted_right(cols * cols, T {});
      buffer<double> values(cols, double         {});
      for (std::size_t d {0}; d < cols; ++d) {
        const std::size_t j {order[d]};
        values[d] = norms[j] * scale;
        for (std::size_t i {0}; i < cols; ++i) {
          sorted_right[(i * cols) + d] = right[(i * cols) + j];
        }
        if (norms[j] > 0.0) {
          for (std::size_t i {0}; i < cols; ++i) {
            left[(i * cols) + d] = w[(i * cols) + j] / detail::from_double<T>(norms[j]);
          }
          continue;
        }
        buffer<T>      best;
        double         best_norm {-1.0};
        for (std::size_t k {0}; k < cols; ++k) {
          buffer<T> candidate(cols, T {});
          candidate[k] = T {1};
          orthogonalise(candidate, left, cols, d);
          const double norm {detail::scaled_norm(candidate)};
          if (norm > best_norm) {
            best      = std::move(candidate);
            best_norm = norm;
          }
        }
        orthogonalise(best, left, cols, d); // a second pass restores orthogonality lost to rounding
        const double norm {detail::scaled_norm(best)};
        for (std::size_t i {0}; i < cols; ++i) {
          left[(i * cols) + d] = best[i] / detail::from_double<T>(norm);
        }
      }
      singular_values_ = vector<double> {std::move(values)};
      u_               = left_vectors(std::move(factors.q), std::move(left), rows, cols, width);
      v_               = matrix<T> {cols, cols, std::move(sorted_right)};
    }
 
    static void orthogonalise(buffer<T>&            x,
                              const buffer<T>&      basis,
                              std::size_t           n,
                              std::size_t           count)
    {
      for (std::size_t c {0}; c < count; ++c) {
        T projection {};
        for (std::size_t i {0}; i < n; ++i) {
          projection += detail::conj_of(basis[(i * n) + c]) * x[i];
        }
        for (std::size_t i {0}; i < n; ++i) {
          x[i] -= projection * basis[(i * n) + c];
        }
      }
    }
 
    vector<double> singular_values_; 
    matrix<T>      u_;               
    matrix<T>      v_;               
  };
 
  namespace detail {
 
    [[nodiscard]] inline std::size_t numerical_rank(const vector<double>& sigma,
                                                    std::size_t           rows,
                                                    std::size_t           cols,
                                                    rank_cutoff           cutoff = {})
    {
      double largest {0.0};
      for (std::size_t d {0}; d < sigma.size(); ++d) {
        largest = std::max(largest, sigma[d]);
      }
      const double tolerance {cutoff.magnitude(largest, rows, cols)};
      std::size_t  rank      {0};
      for (std::size_t d {0}; d < sigma.size(); ++d) {
        rank += static_cast<std::size_t>(sigma[d] > tolerance);
      }
      return rank;
    }
 
    inline constexpr double qr_rank_margin {100.0};
 
    template <typename T>
    [[nodiscard]] T scaled_by_two(const T& x,
                                  int      exponent)
    {
      return std::scalbn(x, exponent);
    }
 
    template <typename U>
    [[nodiscard]] std::complex<U> scaled_by_two(const std::complex<U>& x,
                                                int                    exponent)
    {
      return std::complex<U> {std::scalbn(x.real(), exponent), std::scalbn(x.imag(), exponent)};
    }
 
    template <typename T>
    [[nodiscard]] buffer<T> upper_inverse(const buffer<T>&      r,
                                          std::size_t           n)
    {
      constexpr std::size_t width      {64};
      buffer<T>             x(n * n, T {});
      parallel_for((n + width - 1) / width, n * n * width / 2, [&](std::size_t chunk) {
                     const std::size_t first {chunk * width};
                     const std::size_t last {std::min(first + width, n)};
                     for (std::size_t i {last}; i-- > 0;) {
                       T* const row {x.data() + (i * n) + first};
                       if (i >= first) {
                         row[i - first] = T {1};
                       }
                       subtract_products(row, r.data() + (i * n) + i + 1, x.data() + ((i + 1) * n) + first, n, last - (i + 1), last - first);
                       for (std::size_t c {0}; c < last - first; ++c) {
                         row[c] = row[c] / r[(i * n) + i];
                       }
                     }
                   });
      return x;
    }
 
    template <typename T>
    [[nodiscard]] std::optional<matrix<T>> lstsq_by_qr(const matrix<T>&                       a,
                                                       const matrix<T>&                       b,
                                                       rank_cutoff                            cutoff,
                                                       std::pair<buffer<T>, int>* const       r = nullptr)
    {
      const std::size_t m       {a.rows()};
      const std::size_t n       {a.cols()};
      const bool        wide    {m < n};
      const std::size_t rows    {wide ? n : m}; // of the matrix factorised
      const std::size_t cols    {wide ? m : n};
      double            largest {0.0};
      for (std::size_t i {0}; i < m; ++i) {
        for (std::size_t j {0}; j < n; ++j) {
          largest = std::max(largest, static_cast<double>(std::abs(a(i, j))));
        }
      }
      if (largest == 0.0) {
        return std::nullopt;
      }
      const int      exponent               {std::ilogb(largest)};
      buffer<T>      entries(rows * cols, T {});
      for (std::size_t i {0}; i < m; ++i) {
        for (std::size_t j {0}; j < n; ++j) {
          if (wide) {
            entries[(j * cols) + i] = scaled_by_two(conj_of(a(i, j)), -exponent);
          }
          else {
            entries[(i * cols) + j] = scaled_by_two(a(i, j), -exponent);
          }
        }
      }
      qr_reflections<T>       factors         {householder_reflect(std::move(entries), rows, cols)};
      const buffer<T>         inverse         {upper_inverse(factors.r, cols)};
      double                  r_squares       {0.0};
      double                  inverse_squares {0.0};
      for (std::size_t e {0}; e < cols * cols; ++e) {
        r_squares       += static_cast<double>(std::norm(factors.r[e]));
        inverse_squares += static_cast<double>(std::norm(inverse[e]));
      }
      const double smallest_bound {1.0 / std::sqrt(inverse_squares)};
      if (!(smallest_bound > qr_rank_margin * cutoff.scaled(exponent).magnitude(std::sqrt(r_squares), m, n))) {
        return std::nullopt;
      }
      const std::size_t k          {b.cols()};
      buffer<T>         x(n * k, T {});
      if (wide) { // x = Q [R⁻ᴴ B; 0]: the blocks applied last first
        const matrix<T>      inverse_adjoint {adjoint(matrix<T> {cols, cols, inverse})};
        const buffer<T>      right(b.begin(), b.end());
        gemm(cols, k, cols, &*inverse_adjoint.begin(), cols, right.data(), k, x.data(), k);
        for (std::size_t block {factors.blocks.size()}; block-- > 0;) {
          apply_reflector_block(factors.blocks[block], false, x, k, 0, k);
        }
      }
      else { // x = R⁻¹ times the first rows of Qᴴ B: the blocks' adjoints applied first first
        buffer<T> y(b.begin(), b.end());
        for (const reflector_block<T>& block : factors.blocks) {
          apply_reflector_block(block, true, y, k, 0, k);
        }
        gemm(n, k, cols, inverse.data(), cols, y.data(), k, x.data(), k);
      }
      for (T& value : x) {
        value = scaled_by_two(value, -exponent);
      }
      if (r != nullptr) {
        *r = {std::move(factors.r), exponent};
      }
      return matrix<T> {n, k, std::move(x)};
    }
 
    template <typename T>
    [[nodiscard]] buffer<double> square_singular_values(buffer<T>      b,
                                                        std::size_t    n)
    {
      buffer<double> d(n, double {});
      buffer<double> e(n, double {});
      buffer<T>      none;
      buffer<T>      none_too;
      bidiagonalize(b, n, d, e, none, none_too, false);
      bidiagonal_qr(d, e, none, none_too, n, svd_qr_steps_per_value * n);
      return d;
    }
 
    template <typename T>
    [[nodiscard]] matrix<T> lstsq_from_svd(const svd_decomposition<T>& s,
                                           std::size_t                 rank,
                                           const matrix<T>&            b)
    {
      const vector<double>& sigma                {s.singular_values()};
      const std::size_t     rows                 {s.u().rows()};
      const std::size_t     cols                 {s.v().rows()};
      buffer<T>             x(cols * b.cols(), T {});
      for (std::size_t c {0}; c < b.cols(); ++c) {
        for (std::size_t d {0}; d < rank; ++d) {
          T coefficient {}; // (u_dᴴ b_c) / σ_d
          for (std::size_t i {0}; i < rows; ++i) {
            coefficient += conj_of(s.u()(i, d)) * b(i, c);
          }
          coefficient /= from_double<T>(sigma[d]);
          for (std::size_t j {0}; j < cols; ++j) {
            x[(j * b.cols()) + c] += coefficient * s.v()(j, d);
          }
        }
      }
      return matrix<T> {cols, b.cols(), std::move(x)};
    }
  } // namespace detail
 
  template <typename T>
  [[nodiscard]] vector<T> lstsq(const matrix<T>& a,
                                const vector<T>& b,
                                rank_cutoff      cutoff = {})
  {
    const matrix<T> column {lstsq(a, matrix<T> {b.size(), 1, buffer<T>(b.begin(), b.end())}, cutoff)};
    return vector<T> {buffer<T>(column.begin(), column.end())};
  }
 
  template <typename T>
  [[nodiscard]] matrix<T> lstsq(const matrix<T>& a,
                                const matrix<T>& b,
                                rank_cutoff      cutoff = {})
  {
    if (b.rows() != a.rows()) {
      throw dimension_error("lstsq: right-hand side size differs from the matrix");
    }
    detail::require_finite(a, "lstsq");
    if (std::optional<matrix<T>> fast {detail::lstsq_by_qr(a, b, cutoff)}) {
      return *std::move(fast);
    }
    const svd_decomposition<T> s {a};
    return detail::lstsq_from_svd(s, detail::numerical_rank(s.singular_values(), a.rows(), a.cols(), cutoff), b);
  }
 
  template <typename T>
  [[nodiscard]] vector<double> singular_values(const matrix<T>& a)
  {
    detail::require_finite(a, "singular values");
    const bool        tall {a.rows() >= a.cols()};
    const std::size_t rows {tall ? a.rows() : a.cols()};
    const std::size_t cols {tall ? a.cols() : a.rows()};
    // Divided by its largest magnitude, as the decomposition is: the sums of squares overflow near 1e155.
    double scale           {0.0};
    for (std::size_t i {0}; i < a.rows(); ++i) {
      for (std::size_t j {0}; j < a.cols(); ++j) {
        scale = std::max(scale, static_cast<double>(std::abs(a(i, j))));
      }
    }
    scale = scale == 0.0 ? 1.0 : scale;
    buffer<T> entries(rows * cols, T {});
    for (std::size_t i {0}; i < rows; ++i) {
      for (std::size_t j {0}; j < cols; ++j) {
        entries[(i * cols) + j] = (tall ? a(i, j) : detail::conj_of(a(j, i))) / detail::from_double<T>(scale);
      }
    }
    buffer<T>      r {rows == cols ? std::move(entries) : detail::householder_reflect(std::move(entries), rows, cols).r};
    buffer<double> d {detail::square_singular_values(std::move(r), cols)};
    for (double& value : d) {
      value *= scale;
    }
    return vector<double> {std::move(d)};
  }
 
  template <typename T>
  [[nodiscard]] vector<double> symmetric_eigenvalues(const matrix<T>& a)
  {
    const std::size_t n {a.rows()};
    if (a.cols() != n) {
      throw dimension_error("symmetric eigenvalues require a square matrix");
    }
    detail::require_finite_lower(a, "symmetric eigenvalues");
    // Divided by its largest magnitude read, as the decompositions are: a sum of squares overflows near 1e155.
    double scale {0.0};
    for (std::size_t i {0}; i < n; ++i) {
      for (std::size_t j {0}; j < i; ++j) {
        scale = std::max(scale, static_cast<double>(std::abs(a(i, j))));
      }
      scale = std::max(scale, std::abs(detail::real_part(a(i, i))));
    }
    scale = scale == 0.0 ? 1.0 : scale;
    buffer<T> entries(n * n, T {});
    for (std::size_t i {0}; i < n; ++i) {
      for (std::size_t j {0}; j < i; ++j) {
        entries[(i * n) + j] = a(i, j) / detail::from_double<T>(scale);
      }
      entries[(i * n) + i] = detail::from_double<T>(detail::real_part(a(i, i)) / scale);
    }
    buffer<double> d(n, double {});
    if (n > 0) {
      buffer<double> off(n, double {});
      buffer<double> none;
      static_cast<void>(detail::tridiagonalize(entries, n, d, off));
      detail::tridiagonal_ql(d, off, none, n, detail::ql_steps_per_value * n);
    }
    std::sort(d.begin(), d.end());
    for (double& value : d) {
      value *= scale;
    }
    return vector<double> {std::move(d)};
  }
 
  template <typename T>
  struct least_squares_solution
  {
    matrix<T>      solution;        
    vector<double> residuals;       
    std::size_t    rank {0};        
    vector<double> singular_values; 
  };
 
  template <typename T>
  [[nodiscard]] least_squares_solution<T> least_squares(const matrix<T>& a,
                                                        const matrix<T>& b,
                                                        rank_cutoff      cutoff = {})
  {
    if (b.rows() != a.rows()) {
      throw dimension_error("lstsq: right-hand side size differs from the matrix");
    }
    detail::require_finite(a, "lstsq");
    // The singular values from what the solution took: the QR's R (A's over a power of two), else the decomposition.
    std::pair<buffer<T>, int>           r;
    std::optional<matrix<T>>            fast {detail::lstsq_by_qr(a, b, cutoff, &r)};
    std::optional<svd_decomposition<T>> s;
    if (!fast) {
      s.emplace(a);
    }
    buffer<double> values {fast ? detail::square_singular_values(std::move(r.first), std::min(a.rows(), a.cols()))
                                     : buffer<double>(s->singular_values().begin(), s->singular_values().end())};
    if (fast) {
      for (double& value : values) {
        value = std::ldexp(value, r.second);
      }
    }
    const vector<double> sigma {std::move(values)};
    const std::size_t    rank  {detail::numerical_rank(sigma, a.rows(), a.cols(), cutoff)};
    matrix<T>            x     {fast ? *std::move(fast) : detail::lstsq_from_svd(*s, rank, b)};
    buffer<double>       residuals;
    if (rank == a.cols() && a.rows() > a.cols()) {
      const matrix<T> fitted {a * x};
      residuals.assign(b.cols(), 0.0);
      for (std::size_t i {0}; i < b.rows(); ++i) {
        for (std::size_t c {0}; c < b.cols(); ++c) {
          residuals[c] += std::norm(b(i, c) - fitted(i, c)); // |·|², real or complex
        }
      }
    }
    return least_squares_solution<T> {.solution        = std::move(x),
                                      .residuals       = vector<double> {std::move(residuals)},
                                      .rank            = rank,
                                      .singular_values = sigma};
  }
 
  namespace detail {
 
    template <typename T>
    [[nodiscard]] matrix<T> hermitian_from_lower(const matrix<T>& a)
    {
      const std::size_t n                {a.rows()};
      buffer<T>         entries(n * n, T {});
      for (std::size_t i {0}; i < n; ++i) {
        for (std::size_t j {0}; j < i; ++j) {
          entries[(i * n) + j] = a(i, j);
          entries[(j * n) + i] = conj_of(a(i, j));
        }
        entries[(i * n) + i] = from_double<T>(real_part(a(i, i)));
      }
      return matrix<T> {n, n, std::move(entries)};
    }
 
    template <typename T>
    [[nodiscard]] matrix<T> inverse_on_kept(const vector<double>& sigma,
                                            const matrix<T>&      right,
                                            const matrix<T>&      left,
                                            rank_cutoff           cutoff)
    {
      const std::size_t rows    {left.rows()};
      const std::size_t cols    {right.rows()};
      double            largest {0.0};
      for (std::size_t d {0}; d < sigma.size(); ++d) {
        largest = std::max(largest, std::abs(sigma[d]));
      }
      const double             tolerance {cutoff.magnitude(largest, rows, cols)};
      std::vector<std::size_t> kept;
      for (std::size_t d {0}; d < sigma.size(); ++d) {
        if (std::abs(sigma[d]) > tolerance) {
          kept.push_back(d);
        }
      }
      const std::size_t rank                        {kept.size()};
      buffer<T>         scaled(cols * rank, T       {});
      buffer<T>         left_adjoint(rank * rows, T {});
      for (std::size_t e {0}; e < rank; ++e) {
        const std::size_t d       {kept[e]};
        const auto        inverse {from_double<T>(1.0 / sigma[d])};
        for (std::size_t i {0}; i < cols; ++i) {
          scaled[(i * rank) + e] = right(i, d) * inverse;
        }
        for (std::size_t j {0}; j < rows; ++j) {
          left_adjoint[(e * rows) + j] = conj_of(left(j, d));
        }
      }
      buffer<T> out(cols * rows, T {});
      gemm(cols, rows, rank, scaled.data(), rank, left_adjoint.data(), rows, out.data(), rows);
      return matrix<T> {cols, rows, std::move(out)};
    }
 
    template <typename T>
    [[nodiscard]] vector<double> hermitian_singular_values(const matrix<T>& a)
    {
      return singular_values(hermitian_from_lower(a));
    }
 
    template <typename T>
    [[nodiscard]] vector<double> eigenvalue_magnitudes(const matrix<T>& a)
    {
      const vector<double> signed_values {symmetric_eigenvalues(a)};
      buffer<double>       values(signed_values.begin(), signed_values.end());
      for (double& value : values) {
        value = std::abs(value);
      }
      return vector<double> {std::move(values)};
    }
 
    [[nodiscard]] inline vector<double> hermitian_singular_values(const matrix<double>& a)
    {
      return eigenvalue_magnitudes(a);
    }
 
    [[nodiscard]] inline vector<double> hermitian_singular_values(const matrix<std::complex<double>>& a)
    {
      return eigenvalue_magnitudes(a);
    }
 
    template <typename T>
    [[nodiscard]] matrix<T> hermitian_pseudo_inverse(const matrix<T>& a,
                                                     rank_cutoff      cutoff)
    {
      const svd_decomposition<T> s {hermitian_from_lower(a)};
      return inverse_on_kept(s.singular_values(), s.v(), s.u(), cutoff);
    }
 
    [[nodiscard]] inline matrix<double> hermitian_pseudo_inverse(const matrix<double>& a,
                                                                 rank_cutoff           cutoff)
    {
      const symmetric_eigen decomposition {a};
      const matrix<double>& vectors       {decomposition.vectors()};
      return inverse_on_kept(decomposition.values(), vectors, vectors, cutoff);
    }
 
    [[nodiscard]] inline matrix<std::complex<double>> hermitian_pseudo_inverse(const matrix<std::complex<double>>& a,
                                                                               rank_cutoff                         cutoff)
    {
      const hermitian_eigen               decomposition {a};
      const matrix<std::complex<double>>& vectors       {decomposition.vectors()};
      return inverse_on_kept(decomposition.values(), vectors, vectors, cutoff);
    }
  } // namespace detail
 
  template <typename T>
  [[nodiscard]] std::size_t matrix_rank(const matrix<T>& a,
                                        rank_cutoff      cutoff    = {},
                                        bool             hermitian = false)
  {
    if (hermitian) {
      if (a.cols() != a.rows()) {
        throw dimension_error("matrix_rank: a Hermitian matrix is square");
      }
      detail::require_finite_lower(a, "matrix_rank");
      return detail::numerical_rank(detail::hermitian_singular_values(a), a.rows(), a.cols(), cutoff);
    }
    return detail::numerical_rank(singular_values(a), a.rows(), a.cols(), cutoff);
  }
 
  template <typename T>
  [[nodiscard]] matrix<T> pinv(const matrix<T>& a,
                               rank_cutoff      cutoff    = {},
                               bool             hermitian = false)
  {
    if (hermitian) {
      if (a.cols() != a.rows()) {
        throw dimension_error("pinv: a Hermitian matrix is square");
      }
      detail::require_finite_lower(a, "pinv");
      return detail::hermitian_pseudo_inverse(a, cutoff);
    }
    detail::require_finite(a, "pinv");
    const svd_decomposition<T> s {a};
    return detail::inverse_on_kept(s.singular_values(), s.v(), s.u(), cutoff);
  }
} // namespace numerical
 
#endif // NUMERICAL_DECOMPOSITION_HPP