Source

`stdlib/linalg/routines.cpp`

2// Original work; see ACKNOWLEDGMENTS.md for the open-source ideas we build upon.

3#include "routines.hpp"

5#include <algorithm>

6#include <cmath>

7#include <complex>

8#include <limits>

9#include <stdexcept>

10#include <vector>

12// Dense linear-algebra routines on ndarray::NDArray. Algorithms reimplemented from

13// standard numerical methods (LU w/ partial pivoting, Cholesky, Householder QR,

// Golub–Reinsch SVD (bidiagonalization + implicit QR), Householder-tridiagonal + QL symmetric eigen, Hessenberg+shifted-QR for

15// the general real spectrum). Hot loops are contiguous so -O3 -march=native

16// auto-vectorizes them (SIMD). The matrices are real (double) but the general

17// eigensolvers return a COMPLEX spectrum (CNDArray) — a real matrix can have

18// complex conjugate eigenvalue pairs — built from the real arithmetic below.

19namespace cheatah::linalg {

21/// @cond INTERNAL

22using ndarray::NDArray;

23/// @endcond

25namespace {

27// ---- extract / build contiguous row-major matrices & vectors ----

28//

29// IMPORTANT: read the shared buffer directly. The element accessor `a.at({i, j})`

30// constructs a `std::vector` index *per call* (one heap allocation per element), so

31// the old extractors did rows*cols allocations just to read a matrix. These pack

32// C-order with a flat `copy_n` when the array is already contiguous (the common

33// case — a freshly built matrix/vector), and a direct strided walk otherwise. No

34// per-element allocation either way.

36// Pack `a`'s elements into `out` (size a.size()) in C-order via direct buffer

37// indexing. Used only for the non-contiguous (view/broadcast/permuted) fallback.

38template <ndarray::Field T>

39void pack_corder(const ndarray::basic_ndarray<T>& a, T* out) {

40 const T* base = a.buffer()->data();

41 const auto& shp = a.shape();

42 const auto& st = a.strides();

43 const std::size_t nd = shp.size();

44 const std::ptrdiff_t off0 = static_cast<std::ptrdiff_t>(a.offset());

45 std::vector<std::size_t> idx(nd, 0);

46 const std::size_t total = a.size();

47 for (std::size_t lin = 0; lin < total; ++lin) {

48 std::ptrdiff_t off = off0;

49 for (std::size_t d = 0; d < nd; ++d)

50 off += static_cast<std::ptrdiff_t>(idx[d]) * st[d];

51 out[lin] = base[static_cast<std::size_t>(off)];

52 for (std::size_t d = nd; d-- > 0;) { // C-order increment

53 if (++idx[d] < shp[d]) break;

54 idx[d] = 0;

55 }

56 }

57}

59// A read-only contiguous C-order pointer to `a`'s data. Zero-copy when `a` is

60// already contiguous (returns straight into its buffer); otherwise packs into

61// `scratch`. Use for routines that only READ their operands (the products).

62template <ndarray::Field T>

63const T* contig(const ndarray::basic_ndarray<T>& a, std::vector<T>& scratch) {

64 if (ndarray::is_contiguous(a)) return a.buffer()->data() + a.offset();

65 scratch.resize(a.size());

66 pack_corder(a, scratch.data());

67 return scratch.data();

68}

70template <ndarray::Field T>

71std::vector<T> as_matrix(const ndarray::basic_ndarray<T>& a, std::size_t& rows, std::size_t& cols) {

72 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");

73 rows = a.shape()[0];

74 cols = a.shape()[1];

75 std::vector<T> m(rows * cols);

76 if (ndarray::is_contiguous(a))

77 std::copy_n(a.buffer()->data() + a.offset(), rows * cols, m.data());

78 else

79 pack_corder(a, m.data());

80 return m;

81}

82// (vector_len — validate a vector shape and return its flattened length — moved to backend.hpp:

83// the generic fronts there share it, and it reads only shape metadata so it serves device

84// containers too.)

85template <ndarray::Field T>

86std::vector<T> as_vector(const ndarray::basic_ndarray<T>& a, std::size_t& n) {

87 n = vector_len(a);

88 std::vector<T> v(n);

89 if (ndarray::is_contiguous(a))

90 std::copy_n(a.buffer()->data() + a.offset(), n, v.data());

91 else

92 pack_corder(a, v.data());

93 return v;

94}

95// Wrap an already-computed buffer as a contiguous NDArray WITHOUT the throwaway

96// zero-init that `NDArray(shape)` would do (it value-fills `product(shape)` elements

97// that we then immediately overwrite — a full wasted pass, ruinous for big results

98// like `outer`). Build straight from the buffer + C-order strides instead.

99// Zero-copy: the result is ALREADY in the ndarray storage type (see ndarray::buffer_t),

100// so move its buffer straight in — no element copy, no second large allocation. Use this

101// for memory-bound results (outer, transpose, kron) where the result is as big as the

102// work and an extra copy would dominate (and, for >128 KiB results, trip glibc's mmap

103// threshold so the copy's fresh pages fault in one by one).

104template <typename T>

105ndarray::basic_ndarray<T> wrap_buffer(std::vector<std::size_t> shape, ndarray::buffer_t<T> data) {

106 auto strides = ndarray::detail::contiguous_strides(shape);

107 auto buf = std::make_shared<ndarray::buffer_t<T>>(std::move(data));

108 return ndarray::basic_ndarray<T>(std::move(buf), std::move(shape), std::move(strides), 0);

109}

110// Plain-std::vector result: one bulk copy into the ndarray storage type. resize

111// (default-init: no zero pass) + std::copy keeps libstdc++'s memmove fast path, so it is a

112// single contiguous pass — negligible next to the O(n³) work of the routines that use it

113// (matmul, inv, the SVD/eig family). (vector::assign through the default-init allocator

114// would instead force an element-by-element copy, which is much slower.)

115template <typename T>

116ndarray::basic_ndarray<T> wrap_buffer(std::vector<std::size_t> shape, std::vector<T> data) {

117 ndarray::buffer_t<T> buf;

118 buf.resize(data.size());

119 std::copy(data.begin(), data.end(), buf.begin());

120 return wrap_buffer<T>(std::move(shape), std::move(buf));

121}

122template <ndarray::Field T>

123ndarray::basic_ndarray<T> make_matrix(std::size_t rows, std::size_t cols, std::vector<T> data) {

124 return wrap_buffer<T>({rows, cols}, std::move(data));

125}

126template <ndarray::Field T>

127ndarray::basic_ndarray<T> make_vector(std::vector<T> data) {

128 const std::size_t n = data.size();

129 return wrap_buffer<T>({n}, std::move(data));

130}

131// Promote a freshly-built (contiguous, offset-0) real result to complex (imag 0).

132CNDArray to_complex(const NDArray& a) {

133 const auto& src = *a.buffer();

134 return wrap_buffer<Cplx>(a.shape(), std::vector<Cplx>(src.begin(), src.end()));

135}

136// Descending order for a complex spectrum: by real part, then imaginary part.

137bool cgreater(const Cplx& x, const Cplx& y) {

138 if (x.real() != y.real()) return x.real() > y.real();

139 return x.imag() > y.imag();

140}

141// (as_matrix / as_vector / make_matrix / make_vector above are templated over Field T, so they

142// serve both real and complex — the former as_cmatrix / make_cmatrix / make_cvector are gone.)

143void require_square(std::size_t r, std::size_t c) {

144 if (r != c) throw std::runtime_error("linalg: expected a square matrix");

145}

146

147// ---- LU decomposition with partial pivoting (in place on a copy) ----

148struct LU {

149 std::vector<double> a; // L (below diag, unit) + U (diag/above), row-major n×n

150 std::vector<std::size_t> piv;

151 double sign;

152 std::size_t n;

153 bool singular;

154};

155LU lu_decompose(std::vector<double> a, std::size_t n) {

156 std::vector<std::size_t> piv(n);

157 std::vector<double> vv(n);

158 double sign = 1.0;

159 bool singular = false;

160 for (std::size_t i = 0; i < n; ++i) {

161 double big = 0.0;

162 for (std::size_t j = 0; j < n; ++j) big = std::max(big, std::fabs(a[i * n + j]));

163 if (big == 0.0) { singular = true; big = 1.0; }

164 vv[i] = 1.0 / big;

165 }

166 for (std::size_t k = 0; k < n; ++k) {

167 double big = 0.0;

168 std::size_t imax = k;

169 for (std::size_t i = k; i < n; ++i) {

170 const double t = vv[i] * std::fabs(a[i * n + k]);

171 if (t > big) { big = t; imax = i; }

172 }

173 if (k != imax) {

174 for (std::size_t j = 0; j < n; ++j) std::swap(a[imax * n + j], a[k * n + j]);

175 sign = -sign;

176 vv[imax] = vv[k];

177 }

178 piv[k] = imax;

179 if (a[k * n + k] == 0.0) { a[k * n + k] = 1e-300; singular = true; }

180 for (std::size_t i = k + 1; i < n; ++i) {

181 const double f = a[i * n + k] / a[k * n + k];

182 a[i * n + k] = f;

183 for (std::size_t j = k + 1; j < n; ++j) a[i * n + j] -= f * a[k * n + j];

184 }

185 }

186 return {std::move(a), std::move(piv), sign, n, singular};

187}

188void lu_solve(const LU& lu, std::vector<double>& b) {

189 const std::size_t n = lu.n;

190 for (std::size_t k = 0; k < n; ++k) std::swap(b[k], b[lu.piv[k]]);

191 for (std::size_t i = 0; i < n; ++i) { // forward (unit L)

192 double s = b[i];

193 for (std::size_t j = 0; j < i; ++j) s -= lu.a[i * n + j] * b[j];

194 b[i] = s;

195 }

196 for (std::size_t i = n; i-- > 0;) { // back (U)

197 double s = b[i];

198 for (std::size_t j = i + 1; j < n; ++j) s -= lu.a[i * n + j] * b[j];

199 b[i] = s / lu.a[i * n + i];

200 }

201}

202// Shared preamble for the LU-based routines (solve/det/slogdet/inv): unpack the operand to a

203// square real workspace and factor it. `LU::n` carries the dimension, so callers need only the tail.

204template <ndarray::Field T, template <typename> class Array>

205LU lu_prepare(const Array<T>& a) {

206 std::size_t n, c;

207 std::vector<double> A = as_matrix(a, n, c);

208 require_square(n, c);

209 return lu_decompose(std::move(A), n);

210}

211

212// ---- Golub–Reinsch SVD: A(m×n) = U(m×n) diag(w) V(n×n)ᵀ, requires m ≥ n ----

213struct SVDc {

214 std::vector<double> u, w, v;

215 std::size_t m, n;

216};

217// Overflow-safe √(a²+b²) for the QR sweeps. std::hypot is correctly-rounded and several

218// times slower; called once per Givens rotation (O(n²) of them) it dominated the

219// values-only SVD. This EISPACK form is plenty accurate and much faster.

220inline double pythag(double a, double b) {

221 const double aa = std::fabs(a), ab = std::fabs(b);

222 if (aa > ab) { const double r = ab / aa; return aa * std::sqrt(1.0 + r * r); }

223 if (ab == 0.0) return 0.0;

224 const double r = aa / ab;

225 return ab * std::sqrt(1.0 + r * r);

226}

227// The world-standard dense SVD (what LAPACK's dgesvd reduces to): Householder

228// bidiagonalization to an upper-bidiagonal B = Uᵦᵀ A Vᵦ, then diagonalization of B by

229// implicit-shift QR, accumulating the orthogonal factors. One reduction plus a

230// quadratically-converging QR sweep — vastly fewer flops than one-sided Jacobi's

231// repeated full passes. On input `a` is m×n row-major; on output it holds U (m×n).

232SVDc svd_golub_reinsch(std::vector<double> a_rm, std::size_t m, std::size_t n,

233 bool want_uv = true) {

234 // When @p want_uv is false only the singular values are produced: the U/V

235 // accumulation and the (dominant) U/V Givens rotations in the QR sweep are skipped

236 // — the same values-only fast path NumPy's `svd(compute_uv=False)` / `svdvals` take,

237 // and what `cond`/`matrix_rank` need.

238 // Work entirely COLUMN-MAJOR: U(r,c) = uc[c*m + r], V(r,c) = vc[c*n + r]. The bulk

239 // of Golub–Reinsch is the length-m LEFT Householder reflectors — the column

240 // reductions and their trailing-column updates. Column-major makes those unit-stride

241 // so -O3 -march=native vectorizes them (FMA over contiguous columns); in row-major

242 // they were stride-n and ran scalar, which is what left the bare SVD behind LAPACK.

243 // The QR sweep (rotating whole U/V columns) is contiguous for the same reason.

244 // Input arrives row-major; transpose it in once; uc holds U on output.

245 std::vector<double> uc(n * m), vc(n * n, 0.0), w(n, 0.0), rv1(n, 0.0), tbuf(m, 0.0);

246 for (std::size_t r = 0; r < m; ++r)

247 for (std::size_t c = 0; c < n; ++c) uc[c * m + r] = a_rm[r * n + c];

248 auto sign = [](double x, double s) { return s >= 0.0 ? std::fabs(x) : -std::fabs(x); };

249 // `g` and `scale` carry across iterations: the super-diagonal rv1[i] is the previous

250 // row-reflector's `scale * g`.

251 double g = 0.0, scale = 0.0, anorm = 0.0;

252

253 // --- Householder reduction to bidiagonal form (diagonal w, super-diagonal rv1) ---

254 for (std::size_t i = 0; i < n; ++i) {

255 const std::size_t l = i + 1;

256 rv1[i] = scale * g;

257 g = 0.0; scale = 0.0;

258 double s = 0.0;

259 double* Ui = &uc[i * m]; // column i — contiguous

260 for (std::size_t k = i; k < m; ++k) scale += std::fabs(Ui[k]);

261 if (scale != 0.0) { // left (column) reflector -> w[i]

262 for (std::size_t k = i; k < m; ++k) { Ui[k] /= scale; s += Ui[k] * Ui[k]; }

263 double f = Ui[i];

264 g = -sign(std::sqrt(s), f);

265 const double h = f * g - s;

266 Ui[i] = f - g;

267 for (std::size_t j = l; j < n; ++j) { // apply to trailing columns

268 double* Uj = &uc[j * m];

269 double sum = 0.0;

270 for (std::size_t k = i; k < m; ++k) sum += Ui[k] * Uj[k]; // contiguous → SIMD

271 const double fr = sum / h;

272 for (std::size_t k = i; k < m; ++k) Uj[k] += fr * Ui[k]; // contiguous → SIMD

273 }

274 for (std::size_t k = i; k < m; ++k) Ui[k] *= scale;

275 }

276 w[i] = scale * g;

277 g = 0.0; scale = 0.0; s = 0.0;

278 // right (row) reflector over columns l..n — length n, the minor half (strided)

279 if (l < n) {

280 for (std::size_t k = l; k < n; ++k) scale += std::fabs(uc[k * m + i]);

281 if (scale != 0.0) {

282 for (std::size_t k = l; k < n; ++k) { uc[k * m + i] /= scale; s += uc[k * m + i] * uc[k * m + i]; }

283 double f = uc[l * m + i];

284 g = -sign(std::sqrt(s), f);

285 const double h = f * g - s;

286 uc[l * m + i] = f - g;

287 for (std::size_t k = l; k < n; ++k) rv1[k] = uc[k * m + i] / h;

288 // Trailing update A(l:m, l:n) += (A·u)·rv1ᵀ, done COLUMN-by-column so the

289 // inner loops sweep contiguous rows of a column (vectorize) — the naive

290 // row-by-row form strode across columns (stride m) and ran scalar.

291 for (std::size_t j = l; j < m; ++j) tbuf[j] = 0.0;

292 for (std::size_t k = l; k < n; ++k) { // t[j] = Σ_k A(j,k)·u[k]

293 const double uk = uc[k * m + i];

294 const double* Ck = &uc[k * m];

295 for (std::size_t j = l; j < m; ++j) tbuf[j] += Ck[j] * uk;

296 }

297 for (std::size_t k = l; k < n; ++k) { // A(j,k) += t[j]·rv1[k]

298 const double r = rv1[k];

299 double* Ck = &uc[k * m];

300 for (std::size_t j = l; j < m; ++j) Ck[j] += tbuf[j] * r;

301 }

302 for (std::size_t k = l; k < n; ++k) uc[k * m + i] *= scale;

303 }

304 }

305 anorm = std::max(anorm, std::fabs(w[i]) + std::fabs(rv1[i]));

306 }

307

308 // --- accumulate the right-hand transformations into V (column-major vc) ---

309 if (want_uv)

310 for (std::size_t i = n; i-- > 0;) {

311 const std::size_t l = i + 1;

312 if (l < n) {

313 if (g != 0.0) {

314 for (std::size_t j = l; j < n; ++j) // V(j,i); double division guards overflow

315 vc[i * n + j] = (uc[j * m + i] / uc[l * m + i]) / g;

316 for (std::size_t j = l; j < n; ++j) {

317 double sum = 0.0;

318 for (std::size_t k = l; k < n; ++k) sum += uc[k * m + i] * vc[j * n + k];

319 for (std::size_t k = l; k < n; ++k) vc[j * n + k] += sum * vc[i * n + k];

320 }

321 }

322 for (std::size_t j = l; j < n; ++j) { vc[j * n + i] = 0.0; vc[i * n + j] = 0.0; }

323 }

324 vc[i * n + i] = 1.0;

325 g = rv1[i];

326 }

327

328 // --- accumulate the left-hand transformations into U (held in uc) ---

329 if (want_uv)

330 for (std::size_t i = n; i-- > 0;) { // min(m,n) == n since m >= n

331 const std::size_t l = i + 1;

332 g = w[i];

333 for (std::size_t j = l; j < n; ++j) uc[j * m + i] = 0.0;

334 double* Ui = &uc[i * m];

335 if (g != 0.0) {

336 g = 1.0 / g;

337 for (std::size_t j = l; j < n; ++j) {

338 double* Uj = &uc[j * m];

339 double sum = 0.0;

340 for (std::size_t k = l; k < m; ++k) sum += Ui[k] * Uj[k]; // contiguous → SIMD

341 const double f = (sum / Ui[i]) * g;

342 for (std::size_t k = i; k < m; ++k) Uj[k] += f * Ui[k]; // contiguous → SIMD

343 }

344 for (std::size_t k = i; k < m; ++k) Ui[k] *= g;

345 } else {

346 for (std::size_t k = i; k < m; ++k) Ui[k] = 0.0;

347 }

348 Ui[i] += 1.0;

349 }

350

351 // U (uc) and V (vc) are already column-major, so the QR sweep's whole-column

352 // rotations below are contiguous and vectorizable — no repacking needed.

353 // --- diagonalize the bidiagonal form: implicit-shift QR with deflation ---

354 const double eps = std::numeric_limits<double>::epsilon();

355 for (std::size_t k = n; k-- > 0;) {

356 for (int its = 0; its < 60; ++its) {

357 bool flag = true;

358 std::size_t l = k, nm = 0;

359 while (true) { // find a negligible super-diagonal to split at

360 if (l == 0) { flag = false; break; } // rv1[0] is structurally 0

361 if (std::fabs(rv1[l]) <= eps * anorm) { flag = false; break; }

362 nm = l - 1;

363 if (std::fabs(w[nm]) <= eps * anorm) break;

364 --l;

365 }

366 if (flag) { // cancel rv1[l] via Givens rotations in U

367 double c = 0.0, s = 1.0;

368 for (std::size_t i = l; i <= k; ++i) {

369 double f = s * rv1[i];

370 rv1[i] = c * rv1[i];

371 if (std::fabs(f) <= eps * anorm) break;

372 double gg = w[i];

373 double h = pythag(f, gg);

374 w[i] = h; h = 1.0 / h;

375 c = gg * h; s = -f * h;

376 if (want_uv) {

377 double* Unm = &uc[nm * m];

378 double* Ui = &uc[i * m];

379 for (std::size_t j = 0; j < m; ++j) {

380 const double y = Unm[j], z = Ui[j];

381 Unm[j] = y * c + z * s;

382 Ui[j] = z * c - y * s;

383 }

384 }

385 }

386 }

387 double z = w[k];

388 if (l == k) { // converged: make the singular value non-negative

389 if (z < 0.0) {

390 w[k] = -z;

391 if (want_uv) { double* Vk = &vc[k * n]; for (std::size_t j = 0; j < n; ++j) Vk[j] = -Vk[j]; }

392 }

393 break;

394 }

395 if (its == 59) throw std::runtime_error("linalg: SVD did not converge");

396 double x = w[l];

397 nm = k - 1;

398 double y = w[nm], gg = rv1[nm], h = rv1[k];

399 double f = ((y - z) * (y + z) + (gg - h) * (gg + h)) / (2.0 * h * y);

400 gg = pythag(f, 1.0);

401 f = ((x - z) * (x + z) + h * ((y / (f + sign(gg, f))) - h)) / x;

402 double c = 1.0, s = 1.0;

403 for (std::size_t j = l; j <= nm; ++j) { // QR sweep: chase the bulge

404 const std::size_t i = j + 1;

405 gg = rv1[i]; y = w[i]; h = s * gg; gg = c * gg;

406 z = pythag(f, h);

407 rv1[j] = z; c = f / z; s = h / z;

408 f = x * c + gg * s; gg = gg * c - x * s; h = y * s; y *= c;

409 if (want_uv) {

410 double* Vj = &vc[j * n];

411 double* Vi = &vc[i * n];

412 for (std::size_t jj = 0; jj < n; ++jj) { // rotate V columns j, i (contiguous)

413 const double vx = Vj[jj], vz = Vi[jj];

414 Vj[jj] = vx * c + vz * s;

415 Vi[jj] = vz * c - vx * s;

416 }

417 }

418 z = pythag(f, h);

419 w[j] = z;

420 if (z != 0.0) { z = 1.0 / z; c = f * z; s = h * z; }

421 f = c * gg + s * y; x = c * y - s * gg;

422 if (want_uv) {

423 double* Uj = &uc[j * m];

424 double* Ui = &uc[i * m];

425 for (std::size_t jj = 0; jj < m; ++jj) { // rotate U columns j, i (contiguous)

426 const double uy = Uj[jj], uz = Ui[jj];

427 Uj[jj] = uy * c + uz * s;

428 Ui[jj] = uz * c - uy * s;

429 }

430 }

431 }

432 rv1[l] = 0.0; rv1[k] = f; w[k] = x;

433 }

434 }

435

436 // singular values come out non-negative but unordered — sort descending, carrying

437 // the matching columns of U and V (read straight from the column-major buffers).

438 std::vector<std::size_t> idx(n);

439 for (std::size_t i = 0; i < n; ++i) idx[i] = i;

440 std::sort(idx.begin(), idx.end(), [&](std::size_t x, std::size_t y) { return w[x] > w[y]; });

441 SVDc out{std::vector<double>(want_uv ? m * n : 0), std::vector<double>(n),

442 std::vector<double>(want_uv ? n * n : 0), m, n};

443 for (std::size_t j = 0; j < n; ++j) out.w[j] = w[idx[j]];

444 if (want_uv)

445 for (std::size_t j = 0; j < n; ++j) {

446 const std::size_t src = idx[j];

447 const double* Uc = &uc[src * m];

448 for (std::size_t i = 0; i < m; ++i) out.u[i * n + j] = Uc[i];

449 const double* Vc = &vc[src * n];

450 for (std::size_t i = 0; i < n; ++i) out.v[i * n + j] = Vc[i];

451 }

452 return out;

453}

454

455// ---- real symmetric eigensolver: Householder tridiagonalization (tred2) + ----

456// ---- implicit-shift QL (tql2). ------------------------------------------------

457// The O(n³) method LAPACK uses (one reduction + a QL sweep that converges in O(n)

458// rotations), far cheaper than cyclic Jacobi's repeated full-matrix sweeps. `a` is a

459// row-major n×n matrix ASSUMED symmetric (only the working triangle is used). Returns

460// eigenvalues DESCENDING in `values`, with the matching orthonormal eigenvector as

461// column j of the row-major `vectors` (vectors[i*n+j] = component i of eigenvector j).

462void symmetric_eig(std::vector<double> z, std::size_t n, std::vector<double>& values,

463 std::vector<double>& vectors, bool want_vectors = true) {

464 values.assign(n, 0.0);

465 vectors.assign(want_vectors ? n * n : 0, 0.0);

466 if (n == 0) return;

467 if (n == 1) { values[0] = z[0]; if (want_vectors) vectors[0] = 1.0; return; }

468

469 std::vector<double> d(n, 0.0), e(n, 0.0);

470

471 // --- tred2: reduce symmetric z -> tridiagonal (d diagonal, e subdiagonal),

472 // leaving the accumulated orthogonal transform in z. ---

473 for (std::size_t i = n - 1; i >= 1; --i) {

474 const std::size_t l = i - 1;

475 double h = 0.0;

476 if (l > 0) {

477 double scale = 0.0;

478 for (std::size_t k = 0; k <= l; ++k) scale += std::fabs(z[i * n + k]);

479 if (scale == 0.0) {

480 e[i] = z[i * n + l];

481 } else {

482 for (std::size_t k = 0; k <= l; ++k) {

483 z[i * n + k] /= scale;

484 h += z[i * n + k] * z[i * n + k];

485 }

486 double f = z[i * n + l];

487 double g = (f >= 0.0) ? -std::sqrt(h) : std::sqrt(h);

488 e[i] = scale * g;

489 h -= f * g;

490 z[i * n + l] = f - g;

491 // The active block [0..l]×[0..l] is kept FULL-symmetric (the rank-2

492 // update below writes both triangles), so the matrix–vector product

493 // p = A·u is a single contiguous, vectorizing row·u dot — no

494 // column-stride walk. u is the Householder vector (row i). 2× the

495 // update flops vs the packed form, but both phases now hit SIMD.

496 f = 0.0;

497 const double* ui = &z[i * n]; // Householder vector u (= row i)

498 for (std::size_t j = 0; j <= l; ++j) {

499 z[j * n + i] = z[i * n + j] / h; // store u/h in column i (for Q)

500 const double* zj = &z[j * n];

501 // full row · u, four independent accumulators so it vectorizes

502 // (a single running sum is FMA-latency-bound — the dot mistake).

503 double g0 = 0, g1 = 0, g2 = 0, g3 = 0;

504 std::size_t k = 0;

505 for (; k + 4 <= l + 1; k += 4) {

506 g0 += zj[k] * ui[k]; g1 += zj[k + 1] * ui[k + 1];

507 g2 += zj[k + 2] * ui[k + 2]; g3 += zj[k + 3] * ui[k + 3];

508 }

509 g = (g0 + g1) + (g2 + g3);

510 for (; k <= l; ++k) g += zj[k] * ui[k];

511 e[j] = g / h;

512 f += e[j] * ui[j];

513 }

514 const double hh = f / (h + h);

515 for (std::size_t j = 0; j <= l; ++j) e[j] -= hh * ui[j]; // e := w = p/h − hh·u

516 // Symmetric rank-2 update A −= u·wᵀ + w·uᵀ over the full block (w fully

517 // formed above, so no in-place hazard); contiguous inner loop.

518 for (std::size_t j = 0; j <= l; ++j) {

519 const double uj = ui[j], wj = e[j];

520 double* zj = &z[j * n];

521 for (std::size_t k = 0; k <= l; ++k) zj[k] -= uj * e[k] + wj * ui[k];

522 }

523 }

524 } else {

525 e[i] = z[i * n + l];

526 }

527 d[i] = h;

528 }

529 d[0] = 0.0;

530 e[0] = 0.0;

531 if (want_vectors) {

532 for (std::size_t i = 0; i < n; ++i) { // accumulate the transform into z

533 if (d[i] != 0.0) {

534 for (std::size_t j = 0; j < i; ++j) {

535 double g = 0.0;

536 for (std::size_t k = 0; k < i; ++k) g += z[i * n + k] * z[k * n + j];

537 for (std::size_t k = 0; k < i; ++k) z[k * n + j] -= g * z[k * n + i];

538 }

539 }

540 d[i] = z[i * n + i];

541 z[i * n + i] = 1.0;

542 for (std::size_t j = 0; j < i; ++j) { z[j * n + i] = 0.0; z[i * n + j] = 0.0; }

543 }

544 } else {

545 for (std::size_t i = 0; i < n; ++i) d[i] = z[i * n + i]; // values only — skip Q

546 }

547

548 // --- tql2: implicit-shift QL on the tridiagonal (d, e), rotating z alongside. ---

549 for (std::size_t i = 1; i < n; ++i) e[i - 1] = e[i];

550 e[n - 1] = 0.0;

551 for (std::size_t l = 0; l < n; ++l) {

552 int iter = 0;

553 std::size_t m;

554 do {

555 for (m = l; m + 1 < n; ++m) {

556 const double dd = std::fabs(d[m]) + std::fabs(d[m + 1]);

557 if (std::fabs(e[m]) <= 2.2e-16 * dd) break;

558 }

559 if (m != l) {

560 if (iter++ == 50)

561 throw std::runtime_error("linalg: symmetric eigen QL did not converge");

562 double g = (d[l + 1] - d[l]) / (2.0 * e[l]);

563 double r = pythag(g, 1.0);

564 g = d[m] - d[l] + e[l] / (g + (g >= 0.0 ? std::fabs(r) : -std::fabs(r)));

565 double s = 1.0, c = 1.0, p = 0.0;

566 bool zeroed = false;

567 for (std::size_t i = m; i-- > l;) { // i = m-1 … l

568 double f = s * e[i];

569 const double b = c * e[i];

570 r = pythag(f, g);

571 e[i + 1] = r;

572 if (r == 0.0) { d[i + 1] -= p; e[m] = 0.0; zeroed = true; break; }

573 s = f / r;

574 c = g / r;

575 g = d[i + 1] - p;

576 r = (d[i] - g) * s + 2.0 * c * b;

577 p = s * r;

578 d[i + 1] = g + p;

579 g = c * r - b;

580 if (want_vectors)

581 for (std::size_t k = 0; k < n; ++k) { // rotate eigenvector columns

582 f = z[k * n + i + 1];

583 z[k * n + i + 1] = s * z[k * n + i] + c * f;

584 z[k * n + i] = c * z[k * n + i] - s * f;

585 }

586 }

587 if (!zeroed) { d[l] -= p; e[l] = g; e[m] = 0.0; }

588 }

589 } while (m != l);

590 }

591

592 // sort DESCENDING, carrying the matching eigenvector columns.

593 std::vector<std::size_t> idx(n);

594 for (std::size_t i = 0; i < n; ++i) idx[i] = i;

595 std::sort(idx.begin(), idx.end(), [&](std::size_t x, std::size_t y) { return d[x] > d[y]; });

596 for (std::size_t j = 0; j < n; ++j) values[j] = d[idx[j]];

597 if (want_vectors)

598 for (std::size_t j = 0; j < n; ++j)

599 for (std::size_t i = 0; i < n; ++i) vectors[i * n + j] = z[i * n + idx[j]];

600}

601

602bool is_symmetric(const std::vector<double>& a, std::size_t n) {

603 for (std::size_t i = 0; i < n; ++i)

604 for (std::size_t j = i + 1; j < n; ++j)

605 if (std::fabs(a[i * n + j] - a[j * n + i]) > 1e-12 * (1 + std::fabs(a[i * n + j])))

606 return false;

607 return true;

608}

609

610// Complex Hermitian eigensolver with REAL eigenvalues and COMPLEX eigenvectors,

611// reusing the real symmetric tridiagonal-QL solver via the standard 2n×2n real embedding:

612// for H = A + iB (A symmetric, B antisymmetric), the real symmetric matrix

613// M = [[A, -B], [B, A]]

614// has each eigenvalue of H twice, and a real eigenvector [x; y] of M corresponds to

615// the complex eigenvector x + iy of H (already unit-norm: |x|²+|y|² = 1). We take

616// one representative per duplicated pair. @p evecs (when requested) is row-major n×n

617// with column k the eigenvector for evals[k]; both come out sorted descending.

618void hermitian_eig(const std::vector<Cplx>& H, std::size_t n, std::vector<double>& evals,

619 std::vector<Cplx>& evecs, bool want_vectors) {

620 const std::size_t N = 2 * n;

621 std::vector<double> M(N * N, 0.0);

622 for (std::size_t i = 0; i < n; ++i)

623 for (std::size_t j = 0; j < n; ++j) {

624 const double re = H[i * n + j].real(), im = H[i * n + j].imag();

625 M[i * N + j] = re; // top-left A

626 M[(i + n) * N + (j + n)] = re; // bottom-right A

627 M[i * N + (j + n)] = -im; // top-right -B

628 M[(i + n) * N + j] = im; // bottom-left B

629 }

630 std::vector<double> w, V;

631 symmetric_eig(M, N, w, V, want_vectors); // 2n eigenvalues (desc, paired) + vectors

632 evals.resize(n);

633 for (std::size_t k = 0; k < n; ++k) evals[k] = w[2 * k]; // one of each equal pair

634 if (want_vectors) {

635 evecs.assign(n * n, Cplx{});

636 for (std::size_t k = 0; k < n; ++k) {

637 const std::size_t col = 2 * k;

638 for (std::size_t p = 0; p < n; ++p) {

639 const double x = V[p * N + col], y = V[(p + n) * N + col];

640 evecs[p * n + k] = Cplx(x, y); // column k = eigenvector for evals[k]

641 }

642 }

643 }

644}

645

646// Complex LU with partial pivoting, factored in place on M (row-major n×n): the unit

647// lower factor's multipliers are stored below the diagonal, U on/above it. Returns the

648// pivot vector. Factor ONCE, then `complex_lu_solve` for each right-hand side — inverse

649// iteration reuses the same (deliberately near-singular) M across several RHS.

650std::vector<std::size_t> complex_lu(std::vector<Cplx>& M, std::size_t n) {

651 std::vector<std::size_t> piv(n);

652 for (std::size_t k = 0; k < n; ++k) {

653 std::size_t p = k;

654 double best = std::abs(M[k * n + k]);

655 for (std::size_t i = k + 1; i < n; ++i) {

656 const double m = std::abs(M[i * n + k]);

657 if (m > best) { best = m; p = i; }

658 }

659 piv[k] = p;

660 if (p != k)

661 for (std::size_t j = 0; j < n; ++j) std::swap(M[k * n + j], M[p * n + j]);

662 const Cplx d = M[k * n + k];

663 for (std::size_t i = k + 1; i < n; ++i) {

664 const Cplx f = M[i * n + k] / d;

665 M[i * n + k] = f;

666 for (std::size_t j = k + 1; j < n; ++j) M[i * n + j] -= f * M[k * n + j];

667 }

668 }

669 return piv;

670}

671// Solve (already-factored) M·x = b in place on @p b (forward unit-L, then back-U).

672void complex_lu_solve(const std::vector<Cplx>& M, const std::vector<std::size_t>& piv,

673 std::vector<Cplx>& b, std::size_t n) {

674 for (std::size_t k = 0; k < n; ++k)

675 if (piv[k] != k) std::swap(b[k], b[piv[k]]);

676 for (std::size_t i = 0; i < n; ++i) {

677 Cplx s = b[i];

678 for (std::size_t j = 0; j < i; ++j) s -= M[i * n + j] * b[j];

679 b[i] = s;

680 }

681 for (std::size_t i = n; i-- > 0;) {

682 Cplx s = b[i];

683 for (std::size_t j = i + 1; j < n; ++j) s -= M[i * n + j] * b[j];

684 b[i] = s / M[i * n + i];

685 }

686}

687

688// Eigenvector of the real matrix @p A for (complex) eigenvalue @p lambda, by inverse

689// iteration. C = A − (λ + tiny complex shift)·I is made just non-singular by the

690// shift, then a few inverse-iteration steps converge to the eigenvector; the phase

691// is fixed so the largest-magnitude component is real-positive (a stable, if

692// arbitrary, choice — eigenvectors are only defined up to phase).

693std::vector<Cplx> eigvector_inverse_iteration(const std::vector<double>& A, std::size_t n,

694 Cplx lambda) {

695 double scale = 1.0;

696 for (double a : A) scale = std::max(scale, std::fabs(a));

697 const Cplx shifted = lambda + Cplx(scale * 1e-10, scale * 1e-10);

698 std::vector<Cplx> C(n * n);

699 for (std::size_t i = 0; i < n; ++i)

700 for (std::size_t j = 0; j < n; ++j)

701 C[i * n + j] = Cplx(A[i * n + j], 0.0) - (i == j ? shifted : Cplx{});

702 const auto normalize = [&](std::vector<Cplx>& x) {

703 double nrm = 0.0;

704 for (const Cplx& z : x) nrm += std::norm(z);

705 nrm = std::sqrt(nrm);

706 for (Cplx& z : x) z /= nrm;

707 };

708 const std::vector<std::size_t> piv = complex_lu(C, n); // factor ONCE, reuse per step

709 std::vector<Cplx> v(n, Cplx(1.0, 0.0));

710 normalize(v);

711 for (int it = 0; it < 5; ++it) {

712 complex_lu_solve(C, piv, v, n); // in place on v — no per-step copy or re-factor

713 normalize(v);

714 }

715 std::size_t mi = 0;

716 double mb = 0.0;

717 for (std::size_t i = 0; i < n; ++i) {

718 const double m = std::abs(v[i]);

719 if (m > mb) {

720 mb = m;

721 mi = i;

722 }

723 }

724 const Cplx phase = v[mi] / std::abs(v[mi]); // unit-norm v -> mb > 0

725 for (Cplx& z : v) z /= phase;

726 return v;

727}

728

729// ---- general eigenvalues: Hessenberg reduction + shifted QR ----

730// Real matrix in; COMPLEX spectrum out (a 2×2 block with negative discriminant is a

731// conjugate pair, not an error). The arithmetic stays real; only the extracted

732// eigenvalues are complex.

733std::vector<Cplx> eigvals_general(std::vector<double> a, std::size_t n) {

734 // Householder reduction to upper Hessenberg.

735 std::vector<double> u(n); // reflector, reused per column (entries < k unused)

736 for (std::size_t k = 1; k + 1 < n; ++k) {

737 double scale = 0.0;

738 for (std::size_t i = k; i < n; ++i) scale += std::fabs(a[i * n + (k - 1)]);

739 if (scale == 0.0) continue;

740 double h = 0.0;

741 for (std::size_t i = k; i < n; ++i) {

742 u[i] = a[i * n + (k - 1)] / scale;

743 h += u[i] * u[i];

744 }

745 double g = (u[k] >= 0 ? -std::sqrt(h) : std::sqrt(h));

746 h -= u[k] * g;

747 u[k] -= g;

748 // A = (I - uuᵀ/h) A (I - uuᵀ/h)

749 for (std::size_t j = 0; j < n; ++j) { // right: columns

750 double f = 0.0;

751 for (std::size_t i = k; i < n; ++i) f += u[i] * a[j * n + i];

752 f /= h;

753 for (std::size_t i = k; i < n; ++i) a[j * n + i] -= f * u[i];

754 }

755 for (std::size_t i = 0; i < n; ++i) { // left: rows

756 double f = 0.0;

757 for (std::size_t j = k; j < n; ++j) f += u[j] * a[j * n + i];

758 f /= h;

759 for (std::size_t j = k; j < n; ++j) a[j * n + i] -= f * u[j];

760 }

761 a[k * n + (k - 1)] = scale * g;

762 for (std::size_t i = k + 1; i < n; ++i) a[i * n + (k - 1)] = 0.0;

763 }

764

765 // Shifted QR on the Hessenberg matrix. A 1×1 block is a real eigenvalue; a 2×2

766 // block is two reals (disc ≥ 0) or a complex conjugate pair (disc < 0).

767 std::vector<Cplx> w(n);

768 std::vector<double> cs, sn; // Givens rotations, reused per QR sweep (clear keeps capacity)

769 long long hi = static_cast<long long>(n) - 1;

770 const double eps = 1e-14;

771 int iter = 0;

772 while (hi >= 0) {

773 long long l = hi;

774 while (l > 0) {

775 const double s = std::fabs(a[(l - 1) * n + (l - 1)]) + std::fabs(a[l * n + l]);

776 if (std::fabs(a[l * n + (l - 1)]) <= eps * (s == 0 ? 1.0 : s)) break;

777 --l;

778 }

779 if (l == hi) { // 1×1 block -> real eigenvalue

780 w[hi] = a[hi * n + hi];

781 --hi;

782 iter = 0;

783 } else if (l == hi - 1) { // 2×2 block

784 const std::size_t p = static_cast<std::size_t>(hi - 1), q = static_cast<std::size_t>(hi);

785 const double app = a[p * n + p], aqq = a[q * n + q];

786 const double apq = a[p * n + q], aqp = a[q * n + p];

787 const double tr = app + aqq, det = app * aqq - apq * aqp;

788 const double disc = tr * tr - 4.0 * det;

789 if (disc >= 0.0) { // two real eigenvalues

790 const double sq = std::sqrt(disc);

791 w[p] = (tr + sq) / 2.0;

792 w[q] = (tr - sq) / 2.0;

793 } else { // complex conjugate pair

794 const double im = std::sqrt(-disc) / 2.0;

795 w[p] = Cplx(tr / 2.0, im);

796 w[q] = Cplx(tr / 2.0, -im);

797 }

798 hi -= 2;

799 iter = 0;

800 } else { // QR sweep with Wilkinson shift

801 if (++iter > 200) throw std::runtime_error("linalg: eigenvalue iteration did not converge");

802 const double shift = a[hi * n + hi];

803 for (long long i = l; i <= hi; ++i) a[i * n + i] -= shift;

804 // one explicit QR step via Givens rotations on the Hessenberg block

805 cs.clear();

806 sn.clear();

807 for (long long i = l; i < hi; ++i) {

808 const double x = a[i * n + i], y = a[(i + 1) * n + i];

809 const double r = pythag(x, y);

810 const double c = r == 0 ? 1.0 : x / r, s = r == 0 ? 0.0 : y / r;

811 cs.push_back(c);

812 sn.push_back(s);

813 for (long long j = i; j <= hi; ++j) {

814 const double t1 = a[i * n + j], t2 = a[(i + 1) * n + j];

815 a[i * n + j] = c * t1 + s * t2;

816 a[(i + 1) * n + j] = -s * t1 + c * t2;

817 }

818 }

819 for (long long i = l; i < hi; ++i) { // RQ: post-multiply

820 const double c = cs[static_cast<std::size_t>(i - l)], s = sn[static_cast<std::size_t>(i - l)];

821 for (long long j = l; j <= i + 1; ++j) {

822 const double t1 = a[j * n + i], t2 = a[j * n + (i + 1)];

823 a[j * n + i] = c * t1 + s * t2;

824 a[j * n + (i + 1)] = -s * t1 + c * t2;

825 }

826 }

827 for (long long i = l; i <= hi; ++i) a[i * n + i] += shift;

828 }

829 }

830 return w;

831}

832

833// ---- user-provided output buffers: reuse the caller's storage, no result NDArray allocation ----

834// Every array-returning routine below also has an `out`-FIRST overload (like matmul and the ndarray

835// elementwise ops) that writes into a caller-supplied array instead of allocating a fresh one, so a

836// hot loop can hand the same scratch every call. `out_buf` validates the destination and returns a

837// writable pointer into it; the memory-bound products/transposes write their kernel STRAIGHT into

838// that pointer (genuinely zero result allocation). `copy_into` places an already-built result into

839// it — used by the O(n³) factorizations, whose internal workspace is allocated regardless and for

840// which the single O(n²) copy is negligible next to the decomposition.

841// The hot-loop overload takes the expected shape as an initializer_list: a braced `{r, c}` at a

842// call site would otherwise materialize a temporary std::vector — ONE HEAP ALLOCATION PER CALL —

843// which is exactly what the out-param forms exist to avoid.

844template <ndarray::Field T>

845T* out_buf(ndarray::basic_ndarray<T>& out, std::initializer_list<std::size_t> shape) {

846 const std::vector<std::size_t>& os = out.shape();

847 if (os.size() != shape.size() || !std::equal(os.begin(), os.end(), shape.begin()) ||

848 !ndarray::is_contiguous(out))

849 throw std::runtime_error("linalg: out must be a contiguous array of the result's shape");

850 return out.buffer()->data() + out.offset();

851}

852template <ndarray::Field T>

853T* out_buf(ndarray::basic_ndarray<T>& out, const std::vector<std::size_t>& shape) {

854 if (out.shape() != shape || !ndarray::is_contiguous(out))

855 throw std::runtime_error("linalg: out must be a contiguous array of the result's shape");

856 return out.buffer()->data() + out.offset();

857}

858// Reject an out that aliases an operand still being READ through a zero-copy `contig` pointer

859// (the products and the transpose). The factorizations copy their inputs out first, so they never

860// call this — out may safely alias the input there.

861template <ndarray::Field T>

862void reject_alias(const ndarray::basic_ndarray<T>& out, const ndarray::basic_ndarray<T>& a) {

863 if (out.buffer().get() == a.buffer().get())

864 throw std::runtime_error("linalg: out must not alias an input (it is not computed in place)");

865}

866// Copy a freshly-built contiguous result into the caller's out buffer (validated, reuses its storage).

867template <ndarray::Field T>

868void copy_into(ndarray::basic_ndarray<T>& out, const ndarray::basic_ndarray<T>& result) {

869 T* dst = out_buf(out, result.shape());

870 std::copy_n(result.buffer()->data() + result.offset(), result.size(), dst);

871}

872

873} // namespace

874

875// ================= public routines =================

876

877// ---- products ----

878// The products only READ their operands, so they take a zero-copy `contig` pointer

879// (straight into the array's own buffer when it is contiguous — the common case)

880// and allocate nothing but the result.

881//

882// Reduction kernels use several independent accumulators. A single running sum

883// serializes the loop on floating-point-add latency (the compiler may not reassociate

884// FP adds without -ffast-math), so a plain `s += x[i]*y[i]` runs at ~one element per

885// FADD latency. Independent lanes break that dependency chain, letting -O3

886// -march=native issue SIMD + FMA and hit memory bandwidth instead of add latency.

887namespace {

888// The reduction kernel over any Field T with a compile-time conjugation choice (@ref Conj) —

889// ONE kernel replacing the former real `ddot` and complex `cdot`, now routed through the shared

890// multi-accumulator reduction @ref cheatah::ndarray::detail::reduce_lanes. The per-element term

891// conjugates the first operand only for a complex element under Conj::Conjugate (Hermitian inner

892// product); for real T or Conj::None the conjugation branch is compiled OUT by `if constexpr`.

893template <ndarray::Field T, Conj C>

894T dot_kernel(const T* x, const T* y, std::size_t n) {

895 return ndarray::detail::reduce_lanes<T>(n, [x, y](std::size_t i) -> T {

896 if constexpr (ndarray::is_complex_v<T> && C == Conj::Conjugate) return std::conj(x[i]) * y[i];

897 else return x[i] * y[i];

898 });

899}

900// Read two operands as contiguous pointers (zero-copy when contiguous, else pack once) and reduce.

901template <ndarray::Field T, Conj C, template <typename> class Array>

902T dot_reduce(const Array<T>& a, const Array<T>& b) {

903 const std::size_t n = vector_len(a), m = vector_len(b);

904 if (n != m) throw std::runtime_error("linalg: dot dimension mismatch");

905 if (ndarray::is_contiguous(a) && ndarray::is_contiguous(b))

906 return dot_kernel<T, C>(a.buffer()->data() + a.offset(), b.buffer()->data() + b.offset(), n);

907 std::vector<T> sa, sb;

908 return dot_kernel<T, C>(contig(a, sa), contig(b, sb), n);

909}

910} // namespace

911

912// dot / vdot / inner over any Field T and (host) container Array — the HOST scalar-out kernels

913// of the backend.hpp reduction pattern (the length validation lives in the generic fronts there;

914// dot_reduce's own check is a harmless second line of defense). `dot`/`inner` are bilinear

915// (Σ aᵢbᵢ); `vdot` is the conjugate-linear Hermitian inner product Σ conj(aᵢ)·bᵢ (identical to dot

916// for a real element). Both operands are Array<T> (the deduction firewall).

917template <ndarray::Field T, template <typename> class Array>

918 requires HostArray<Array<T>>

919void dot(T& out, const Array<T>& a, const Array<T>& b) { out = dot_reduce<T, Conj::None>(a, b); }

920template <ndarray::Field T, template <typename> class Array>

921 requires HostArray<Array<T>>

922void vdot(T& out, const Array<T>& a, const Array<T>& b) { out = dot_reduce<T, Conj::Conjugate>(a, b); }

923template <ndarray::Field T, template <typename> class Array>

924 requires HostArray<Array<T>>

925void inner(T& out, const Array<T>& a, const Array<T>& b) { out = dot_reduce<T, Conj::None>(a, b); }

926// Explicit instantiations: the host real + complex kernels, AND the (now header-inline) allocating

927// fronts — instantiating the fronts here keeps the exported symbols the library always shipped.

928template void dot<double, ndarray::basic_ndarray>(double&, const NDArray&, const NDArray&);

929template void dot<Cplx, ndarray::basic_ndarray>(Cplx&, const CNDArray&, const CNDArray&);

930template void vdot<double, ndarray::basic_ndarray>(double&, const NDArray&, const NDArray&);

931template void vdot<Cplx, ndarray::basic_ndarray>(Cplx&, const CNDArray&, const CNDArray&);

932template void inner<double, ndarray::basic_ndarray>(double&, const NDArray&, const NDArray&);

933template double dot<double, ndarray::basic_ndarray>(const NDArray&, const NDArray&);

934template Cplx dot<Cplx, ndarray::basic_ndarray>(const CNDArray&, const CNDArray&);

935template double vdot<double, ndarray::basic_ndarray>(const NDArray&, const NDArray&);

936template Cplx vdot<Cplx, ndarray::basic_ndarray>(const CNDArray&, const CNDArray&);

937template double inner<double, ndarray::basic_ndarray>(const NDArray&, const NDArray&);

938

939namespace {

940// outer-product kernel over any Field T: writes rp[n×m] = x[i]·y[j]. Loop-invariant xi + a clean

941// row pointer keep the inner store contiguous so it vectorizes.

942template <ndarray::Field T>

943void outer_kernel(T* rp, const T* x, const T* y, std::size_t n, std::size_t m) {

944 for (std::size_t i = 0; i < n; ++i) {

945 const T xi = x[i]; // loop-invariant scalar…

946 T* ri = rp + i * m; // …and a clean row pointer, so the inner

947 for (std::size_t j = 0; j < m; ++j) ri[j] = xi * y[j]; // store vectorizes

948 }

949}

950} // namespace

951

952// Outer product a⊗b (rank-1 n×m matrix) into the caller's buffer — the HOST out-parameter form

953// (two-layer over element T and container Array). Writes the kernel straight into @p out.

954template <ndarray::Field T, template <typename> class Array>

955 requires HostArray<Array<T>>

956void outer(Array<T>& out, const Array<T>& a, const Array<T>& b) {

957 const std::size_t n = vector_len(a), m = vector_len(b);

958 reject_alias(out, a);

959 reject_alias(out, b);

960 T* rp = out_buf(out, {n, m});

961 std::vector<T> sa, sb;

962 outer_kernel<T>(rp, contig(a, sa), contig(b, sb), n, m);

963}

964// (The allocating front is inline in backend.hpp — the matmul pattern; instantiating it here

965// keeps the exported symbol the library always shipped.)

966template void outer<double, ndarray::basic_ndarray>(NDArray&, const NDArray&, const NDArray&);

967template NDArray outer<double, ndarray::basic_ndarray>(const NDArray&, const NDArray&);

968

969namespace {

970// The matmul kernel over ANY Field T (real or complex). The loop is element-generic — the

971// only element-specific step is the `T{}` zero-fill — so ONE kernel now serves what used to be

972// a `double*` and a `Cplx*` overload. Writes C[ar×bc] = A[ar×ac]·B[ac×bc] into the caller's @p C

973// (zeroed, then accumulated). ikj keeps the inner (j) loop contiguous so it vectorizes; blocking

974// FOUR rows of A reuses each B[k][j] load across four C rows (4 FMAs per B load instead of 1).

975template <ndarray::Field T>

976void matmul_kernel(T* C, const T* A, const T* B, std::size_t ar, std::size_t ac, std::size_t bc) {

977 std::fill(C, C + ar * bc, T{});

978 std::size_t i = 0;

979 for (; i + 4 <= ar; i += 4) {

980 T* c0 = &C[(i + 0) * bc]; T* c1 = &C[(i + 1) * bc];

981 T* c2 = &C[(i + 2) * bc]; T* c3 = &C[(i + 3) * bc];

982 for (std::size_t k = 0; k < ac; ++k) {

983 const T a0 = A[(i + 0) * ac + k], a1 = A[(i + 1) * ac + k];

984 const T a2 = A[(i + 2) * ac + k], a3 = A[(i + 3) * ac + k];

985 const T* bk = &B[k * bc];

986 for (std::size_t j = 0; j < bc; ++j) {

987 const T bkj = bk[j];

988 c0[j] += a0 * bkj; c1[j] += a1 * bkj; c2[j] += a2 * bkj; c3[j] += a3 * bkj;

989 }

990 }

991 }

992 for (; i < ar; ++i) { // remainder rows (ar not a multiple of 4)

993 T* ci = &C[i * bc];

994 for (std::size_t k = 0; k < ac; ++k) {

995 const T aik = A[i * ac + k];

996 const T* bk = &B[k * bc];

997 for (std::size_t j = 0; j < bc; ++j) ci[j] += aik * bk[j];

998 }

999 }

1000}

1001// Shared 2-D shape validation → (ar, ac, bc); throws on a non-2-D input or inner-dim mismatch.

1002template <ndarray::Field T>

1003void check_matmul(const ndarray::basic_ndarray<T>& a, const ndarray::basic_ndarray<T>& b,

1004 std::size_t& ar, std::size_t& ac, std::size_t& bc) {

1005 if (a.ndim() != 2 || b.ndim() != 2)

1006 throw std::runtime_error("linalg: matmul expects 2-D matrices");

1007 ar = a.shape()[0]; ac = a.shape()[1];

1008 const std::size_t br = b.shape()[0]; bc = b.shape()[1];

1009 if (ac != br) throw std::runtime_error("linalg: matmul inner dimension mismatch");

1010}

1011} // namespace

1012

1013// Matmul into the caller's buffer @p out (out FIRST) — the HOST out-parameter kernel (the two-layer

1014// `template <Field T, template<typename> class Array> requires HostArray<Array<T>>` overload declared

1015// in backend.hpp). ONE definition unifying the former real and complex out-param functions. Validates

1016// shapes, rejects aliasing (out reads all of A and B while writing, so it is not in-place), packs a

1017// strided operand once, and runs the single matmul_kernel. The allocating matmul(a,b) front calls it.

1018template <ndarray::Field T, template <typename> class Array>

1019 requires HostArray<Array<T>>

1020void matmul(Array<T>& out, const Array<T>& a, const Array<T>& b) {

1021 if (a.ndim() == 3) {

1022 // Batched [B,M,K] @ [B,K,N]: the same single-matrix kernel per contiguous batch slice

1023 // (validated by the front — equal batch counts, matching inner dims). The 2-D control

1024 // flow below is untouched.

1025 const std::size_t B = a.shape()[0], M = a.shape()[1], K = a.shape()[2];

1026 const std::size_t N = b.shape()[2];

1027 reject_alias(out, a);

1028 reject_alias(out, b);

1029 T* C = out_buf(out, {B, M, N});

1030 std::vector<T> sa, sb;

1031 const T* A = contig(a, sa);

1032 const T* Bp = contig(b, sb);

1033 for (std::size_t z = 0; z < B; ++z)

1034 matmul_kernel<T>(C + z * M * N, A + z * M * K, Bp + z * K * N, M, K, N);

1035 return;

1036 }

1037 std::size_t ar, ac, bc;

1038 check_matmul(a, b, ar, ac, bc);

1039 reject_alias(out, a);

1040 reject_alias(out, b);

1041 T* C = out_buf(out, {ar, bc});

1042 std::vector<T> sa, sb;

1043 matmul_kernel<T>(C, contig(a, sa), contig(b, sb), ar, ac, bc);

1044}

1045// Explicit instantiations for the two host element types the library ships — both the out-param

1046// kernel and the allocating front — so the header templates link from other TUs and llvm coverage

1047// attributes their bodies to this TU.

1048template void matmul<double, ndarray::basic_ndarray>(NDArray&, const NDArray&, const NDArray&);

1049template void matmul<Cplx, ndarray::basic_ndarray>(CNDArray&, const CNDArray&, const CNDArray&);

1050template NDArray matmul<double, ndarray::basic_ndarray>(const NDArray&, const NDArray&);

1051template CNDArray matmul<Cplx, ndarray::basic_ndarray>(const CNDArray&, const CNDArray&);

1052

1053namespace {

1054// (Conjugate-)transpose kernel over any Field T: D[c×r] = A[r×c]ᵀ, conjugated for a complex

1055// element (Hermitian adjoint). The conjugation is an `if constexpr` branch — a real element

1056// gets a plain transpose, a complex element the adjoint, from ONE kernel.

1057template <ndarray::Field T>

1058void transpose_kernel(T* D, const T* A, std::size_t r, std::size_t c) {

1059 for (std::size_t i = 0; i < r; ++i)

1060 for (std::size_t j = 0; j < c; ++j) {

1061 if constexpr (ndarray::is_complex_v<T>) D[j * r + i] = std::conj(A[i * c + j]);

1062 else D[j * r + i] = A[i * c + j];

1063 }

1064}

1065} // namespace

1066

1067// Conjugate transpose (Hermitian adjoint) Aᴴ into the caller's buffer — the HOST out-parameter

1068// form (two-layer). For a real element this is a plain transpose (conjugation compiled out).

1069template <ndarray::Field T, template <typename> class Array>

1070 requires HostArray<Array<T>>

1071void conj_transpose(Array<T>& out, const Array<T>& a) {

1072 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");

1073 const std::size_t r = a.shape()[0], c = a.shape()[1];

1074 reject_alias(out, a); // reads A while writing the transposed out — not in place

1075 T* D = out_buf(out, {c, r});

1076 std::vector<T> sa;

1077 transpose_kernel<T>(D, contig(a, sa), r, c);

1078}

1079// (The allocating front is inline in backend.hpp — the matmul pattern; instantiating it here

1080// keeps the exported symbol the library always shipped.)

1081template void conj_transpose<Cplx, ndarray::basic_ndarray>(CNDArray&, const CNDArray&);

1082template CNDArray conj_transpose<Cplx, ndarray::basic_ndarray>(const CNDArray&);

1083template void conj_transpose<double, ndarray::basic_ndarray>(NDArray&, const NDArray&);

1084template NDArray conj_transpose<double, ndarray::basic_ndarray>(const NDArray&);

1085

1086template <ndarray::Field T, template <typename> class Array>

1087 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>

1088void matrix_power(Array<T>& out, const Array<T>& a, long long p) {

1089 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");

1090 const std::size_t r = a.shape()[0], c = a.shape()[1]; // dims only — no copy

1091 require_square(r, c);

1092 // r*r is a product of two dims sized BEFORE any overflow-checked path (a.size()/as_matrix)

1093 // runs, so guard it here: product({r, r}) throws on a size_t wrap instead of silently

1094 // under-allocating `result` and letting the identity-fill write out of bounds.

1095 std::vector<double> result(ndarray::detail::product({r, r}), 0.0);

1096 for (std::size_t i = 0; i < r; ++i) result[i * r + i] = 1.0; // identity

1097 NDArray acc = make_matrix(r, r, std::move(result));

1098 NDArray base = (p < 0) ? inv(a) : a;

1099 long long e = p < 0 ? -p : p;

1100 while (e > 0) {

1101 if (e & 1) acc = matmul(acc, base);

1102 base = matmul(base, base);

1103 e >>= 1;

1104 }

1105 copy_into(out, acc);

1106}

1107

1108namespace {

1109// Kronecker-product kernel over any Field T: K[(ar·br)×(ac·bc)] = A⊗B, each A entry scaling the

1110// whole of B.

1111template <ndarray::Field T>

1112void kron_kernel(T* K, const T* A, const T* B, std::size_t ar, std::size_t ac,

1113 std::size_t br, std::size_t bc) {

1114 const std::size_t kc = ac * bc;

1115 for (std::size_t i = 0; i < ar; ++i)

1116 for (std::size_t j = 0; j < ac; ++j)

1117 for (std::size_t p = 0; p < br; ++p)

1118 for (std::size_t q = 0; q < bc; ++q)

1119 K[(i * br + p) * kc + (j * bc + q)] = A[i * ac + j] * B[p * bc + q];

1120}

1121// Shared 2-D validation → (ar, ac, br, bc); throws on a non-2-D operand.

1122template <ndarray::Field T, template <typename> class Array>

1123void kron_dims(const Array<T>& a, const Array<T>& b, std::size_t& ar, std::size_t& ac,

1124 std::size_t& br, std::size_t& bc) {

1125 if (a.ndim() != 2 || b.ndim() != 2)

1126 throw std::runtime_error("linalg: kron expects 2-D matrices");

1127 ar = a.shape()[0]; ac = a.shape()[1];

1128 br = b.shape()[0]; bc = b.shape()[1];

1129}

1130} // namespace

1131

1132// Kronecker product A⊗B into the caller's buffer — the HOST out-parameter form (two-layer).

1133template <ndarray::Field T, template <typename> class Array>

1134 requires HostArray<Array<T>>

1135void kron(Array<T>& out, const Array<T>& a, const Array<T>& b) {

1136 std::size_t ar, ac, br, bc;

1137 kron_dims(a, b, ar, ac, br, bc);

1138 reject_alias(out, a);

1139 reject_alias(out, b);

1140 // Overflow-check each output dim (a product of two input dims) before it collapses into the

1141 // shape — mirrors the allocating front in backend.hpp so the direct out-param path is guarded too.

1142 T* K = out_buf(out, {ndarray::detail::product({ar, br}), ndarray::detail::product({ac, bc})});

1143 std::vector<T> sa, sb;

1144 kron_kernel<T>(K, contig(a, sa), contig(b, sb), ar, ac, br, bc);

1145}

1146// (The allocating front is inline in backend.hpp — the matmul pattern; instantiating it here

1147// keeps the exported symbol the library always shipped.)

1148template void kron<double, ndarray::basic_ndarray>(NDArray&, const NDArray&, const NDArray&);

1149template NDArray kron<double, ndarray::basic_ndarray>(const NDArray&, const NDArray&);

1150

1151// Trace (sum of the diagonal) — the HOST scalar-out kernel of the backend.hpp reduction pattern

1152// (the 2-D validation lives in the generic front there). Reads the diagonal straight from the

1153// buffer, no copy, even for a strided view.

1154template <ndarray::Field T, template <typename> class Array>

1155 requires HostArray<Array<T>>

1156void trace(T& out, const Array<T>& a) {

1157 const std::size_t r = a.shape()[0], c = a.shape()[1];

1158 const T* base = a.buffer()->data();

1159 const std::ptrdiff_t off = static_cast<std::ptrdiff_t>(a.offset());

1160 const std::ptrdiff_t step = a.strides()[0] + a.strides()[1]; // (i,i) advances by s0+s1

1161 // Diagonal sum through the shared multi-accumulator reduction — the term is a strided read.

1162 out = ndarray::detail::reduce_lanes<T>(std::min(r, c), [base, off, step](std::size_t i) {

1163 return base[static_cast<std::size_t>(off + static_cast<std::ptrdiff_t>(i) * step)];

1164 });

1165}

1166template void trace<double, ndarray::basic_ndarray>(double&, const NDArray&);

1167template double trace<double, ndarray::basic_ndarray>(const NDArray&);

1168

1169// norm — the HOST scalar-out kernel (Frobenius for matrices / L2 for vectors — same flat sum).

1170template <ndarray::Field T, template <typename> class Array>

1171 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>

1172void norm(T& out, const Array<T>& a) {

1173 // Frobenius/L2 norm is sqrt(x·x); reuse the multi-accumulator dot_kernel so the

1174 // squared-sum reaches memory bandwidth instead of serializing on FP-add latency.

1175 // Contiguous fast path reads straight from the buffer (no scratch allocation).

1176 if (ndarray::is_contiguous(a)) {

1177 const T* p = a.buffer()->data() + a.offset();

1178 out = std::sqrt(dot_kernel<T, Conj::None>(p, p, a.size()));

1179 return;

1180 }

1181 std::vector<T> scratch;

1182 const T* p = contig(a, scratch);

1183 out = std::sqrt(dot_kernel<T, Conj::None>(p, p, a.size()));

1184}

1185template void norm<double, ndarray::basic_ndarray>(double&, const NDArray&);

1186template double norm<double, ndarray::basic_ndarray>(const NDArray&);

1187

1188// ---- LU-based: solve / det / slogdet / inv / lstsq ----

1189// LU-based solve / det / slogdet / inv — the HOST out-param/scalar-out kernels of the routines.hpp

1190// seam pattern (the allocating fronts are inline in routines.hpp; they validate metadata and call

1191// these unqualified, so a device extension's DeviceArray overloads are found by ADL). Constrained

1192// to a real floating element (the LU core is real double). Only `double` is shipped; the internal

1193// helpers (as_matrix/as_vector/make_vector, all templated over T) deduce the element.

1194template <ndarray::Field T, template <typename> class Array>

1195 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>

1196void solve(Array<T>& out, const Array<T>& a, const Array<T>& b) {

1197 const LU lu = lu_prepare(a);

1198 const std::size_t n = lu.n;

1199 std::size_t bn;

1200 std::vector<T> x = as_vector(b, bn);

1201 if (bn != n) throw std::runtime_error("linalg: solve dimension mismatch");

1202 lu_solve(lu, x);

1203 copy_into(out, make_vector(std::move(x)));

1204}

1205

1206template <ndarray::Field T, template <typename> class Array>

1207 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>

1208void det(T& out, const Array<T>& a) {

1209 const LU lu = lu_prepare(a);

1210 const std::size_t n = lu.n;

1211 T d = lu.sign;

1212 for (std::size_t i = 0; i < n; ++i) d *= lu.a[i * n + i];

1213 out = d;

1214}

1215

1216template <ndarray::Field T, template <typename> class Array>

1217 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>

1218void slogdet(SLogDet& out, const Array<T>& a) {

1219 const LU lu = lu_prepare(a);

1220 const std::size_t n = lu.n;

1221 double sign = lu.sign, logabs = 0.0;

1222 for (std::size_t i = 0; i < n; ++i) {

1223 const double d = lu.a[i * n + i];

1224 if (d < 0) sign = -sign;

1225 logabs += std::log(std::fabs(d));

1226 }

1227 out = {sign, logabs};

1228}

1229

1230template <ndarray::Field T, template <typename> class Array>

1231 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>

1232void inv(Array<T>& out, const Array<T>& a) {

1233 const LU lu = lu_prepare(a);

1234 const std::size_t n = lu.n;

1235 const std::vector<double>& M = lu.a; // L (unit, below diag) + U (on/above), row-major

1236 // Invert by solving L·U·X = P·I for the WHOLE identity at once. Doing the forward

1237 // and back substitution across all n columns turns each inner loop into a SAXPY

1238 // over a contiguous row (`X[i,:] -= M[i,j]·X[j,:]`), which auto-vectorizes — unlike

1239 // n separate single-RHS solves, whose substitution is a serial-reduction dot that

1240 // cannot vectorize (the reason a naive `inv` lost to LAPACK while `det` won).

1241 std::vector<double> X(n * n, 0.0);

1242 for (std::size_t i = 0; i < n; ++i) X[i * n + i] = 1.0; // identity

1243 for (std::size_t k = 0; k < n; ++k) // apply LU's row pivots: X = P·I

1244 if (lu.piv[k] != k)

1245 for (std::size_t col = 0; col < n; ++col) std::swap(X[k * n + col], X[lu.piv[k] * n + col]);

1246 for (std::size_t i = 0; i < n; ++i) // forward: unit-lower L·Y = P

1247 for (std::size_t j = 0; j < i; ++j) {

1248 const double f = M[i * n + j];

1249 for (std::size_t col = 0; col < n; ++col) X[i * n + col] -= f * X[j * n + col];

1250 }

1251 for (std::size_t i = n; i-- > 0;) { // back: upper U·X = Y

1252 for (std::size_t j = i + 1; j < n; ++j) {

1253 const double f = M[i * n + j];

1254 for (std::size_t col = 0; col < n; ++col) X[i * n + col] -= f * X[j * n + col];

1255 }

1256 const double d = M[i * n + i];

1257 for (std::size_t col = 0; col < n; ++col) X[i * n + col] /= d;

1258 }

1259 copy_into(out, make_matrix(n, n, std::move(X)));

1260}

1261// Explicit instantiations of the LU family: the host kernels AND the (now header-inline)

1262// allocating fronts — instantiating the fronts here keeps the exported symbols the library ships.

1263template void solve<double, ndarray::basic_ndarray>(NDArray&, const NDArray&, const NDArray&);

1264template void det<double, ndarray::basic_ndarray>(double&, const NDArray&);

1265template void slogdet<double, ndarray::basic_ndarray>(SLogDet&, const NDArray&);

1266template void inv<double, ndarray::basic_ndarray>(NDArray&, const NDArray&);

1267template NDArray solve<double, ndarray::basic_ndarray>(const NDArray&, const NDArray&);

1268template double det<double, ndarray::basic_ndarray>(const NDArray&);

1269template SLogDet slogdet<double, ndarray::basic_ndarray>(const NDArray&);

1270template NDArray inv<double, ndarray::basic_ndarray>(const NDArray&);

1271

1272template <ndarray::Field T, template <typename> class Array>

1273 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>

1274void lstsq(Array<T>& out, const Array<T>& a, const Array<T>& b) { // min ‖Ax−b‖ via the pseudo-inverse

1275 Array<T> p = Array<T>::uninitialized({a.shape()[1], a.shape()[0]});

1276 pinv(p, a);

1277 matmul(out, p, b);

1278}

1279

1280// ---- Cholesky ----

1281template <ndarray::Field T, template <typename> class Array>

1282 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>

1283void cholesky(Array<T>& out, const Array<T>& a) {

1284 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");

1285 const std::size_t n = a.shape()[0], c = a.shape()[1];

1286 require_square(n, c);

1287 std::vector<double> scratch;

1288 const double* A = contig(a, scratch); // read-only — zero-copy when contiguous

1289 std::vector<double> L(n * n, 0.0);

1290 for (std::size_t i = 0; i < n; ++i) {

1291 const double* Li = &L[i * n];

1292 for (std::size_t j = 0; j <= i; ++j) {

1293 // s = A[i][j] − (row i · row j over k<j): four accumulators so the O(n³)

1294 // inner dot vectorizes instead of serializing on FP-sub latency.

1295 const double* Lj = &L[j * n];

1296 double d0 = 0, d1 = 0, d2 = 0, d3 = 0;

1297 std::size_t k = 0;

1298 for (; k + 4 <= j; k += 4) {

1299 d0 += Li[k] * Lj[k]; d1 += Li[k + 1] * Lj[k + 1];

1300 d2 += Li[k + 2] * Lj[k + 2]; d3 += Li[k + 3] * Lj[k + 3];

1301 }

1302 double s = A[i * n + j] - ((d0 + d1) + (d2 + d3));

1303 for (; k < j; ++k) s -= Li[k] * Lj[k];

1304 if (i == j) {

1305 if (s <= 0.0) throw std::runtime_error("linalg: matrix is not positive-definite");

1306 L[i * n + i] = std::sqrt(s);

1307 } else {

1308 L[i * n + j] = s / L[j * n + j];

1309 }

1310 }

1311 }

1312 copy_into(out, make_matrix(n, n, std::move(L)));

1313}

1314

1315// ---- Householder QR (reduced: Q is m×n, R is n×n) ----

1316template <ndarray::Field T, template <typename> class Array>

1317 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>

1318void qr(Array<T>& q, Array<T>& r, const Array<T>& a) {

1319 std::size_t m, n;

1320 std::vector<double> A = as_matrix(a, m, n);

1321 if (m < n) throw std::runtime_error("linalg: qr requires rows >= cols");

1322 // Work on the TRANSPOSE At (n×m, row-major). A Householder QR repeatedly reads and

1323 // updates COLUMNS of A, which stride by n in row-major and don't vectorize (the

1324 // original cost ~3× Eigen); as ROWS of At those same operations are contiguous, and

1325 // the reductions are multi-accumulated like ddot.

1326 std::vector<double> At(n * m);

1327 for (std::size_t i = 0; i < m; ++i)

1328 for (std::size_t j = 0; j < n; ++j) At[j * m + i] = A[i * n + j];

1329 std::vector<double> Q(m * m, 0.0);

1330 for (std::size_t i = 0; i < m; ++i) Q[i * m + i] = 1.0;

1331 std::vector<double> u(m); // Householder vector, reused per column (entries < k unused)

1332 // Reflect: s = u · row over [k, m) via the shared multi-accumulator reduction, then

1333 // row -= (2 s / ‖u‖²) u — contiguous.

1334 auto reflect = [&u](double* row, std::size_t k, std::size_t m, double inv) {

1335 double s = ndarray::detail::reduce_lanes<double>(

1336 m - k, [&u, row, k](std::size_t i) { return u[k + i] * row[k + i]; });

1337 s *= inv;

1338 for (std::size_t i = k; i < m; ++i) row[i] -= s * u[i];

1339 };

1340 for (std::size_t k = 0; k < n; ++k) {

1341 double* Atk = &At[k * m]; // column k of A == row k of At

1342 double nrm = 0.0;

1343 for (std::size_t i = k; i < m; ++i) nrm += Atk[i] * Atk[i];

1344 nrm = std::sqrt(nrm);

1345 if (nrm == 0.0) continue;

1346 const double alpha = Atk[k] >= 0 ? -nrm : nrm;

1347 for (std::size_t i = k; i < m; ++i) u[i] = Atk[i];

1348 u[k] -= alpha;

1349 double unorm2 = 0.0;

1350 for (std::size_t i = k; i < m; ++i) unorm2 += u[i] * u[i];

1351 if (unorm2 == 0.0) continue;

1352 const double inv = 2.0 / unorm2;

1353 for (std::size_t j = k; j < n; ++j) reflect(&At[j * m], k, m, inv); // A's cols j≥k

1354 for (std::size_t j = 0; j < m; ++j) reflect(&Q[j * m], k, m, inv); // Q = Q · Hₖ

1355 }

1356 std::vector<double> Qr(m * n), Rr(n * n, 0.0); // reduced

1357 for (std::size_t i = 0; i < m; ++i)

1358 for (std::size_t j = 0; j < n; ++j) Qr[i * n + j] = Q[i * m + j];

1359 for (std::size_t i = 0; i < n; ++i)

1360 for (std::size_t j = i; j < n; ++j) Rr[i * n + j] = At[j * m + i]; // R[i][j]=A[i][j]=At[j][i]

1361 copy_into(q, make_matrix(m, n, std::move(Qr)));

1362 copy_into(r, make_matrix(n, n, std::move(Rr)));

1363}

1364

1365// ---- SVD and its derived quantities ----

1366template <ndarray::Field T, template <typename> class Array>

1367 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>

1368void svd(Array<T>& u, Array<T>& sv, Array<T>& vhh, const Array<T>& a) {

1369 std::size_t m, n;

1370 std::vector<double> A = as_matrix(a, m, n);

1371 if (m < n) throw std::runtime_error("linalg: svd requires rows >= cols (transpose otherwise)");

1372 const SVDc s = svd_golub_reinsch(std::move(A), m, n);

1373 // vh = Vᵀ

1374 std::vector<double> vh(n * n);

1375 for (std::size_t i = 0; i < n; ++i)

1376 for (std::size_t j = 0; j < n; ++j) vh[i * n + j] = s.v[j * n + i];

1377 copy_into(u, make_matrix(m, n, s.u));

1378 copy_into(sv, make_vector(s.w));

1379 copy_into(vhh, make_matrix(n, n, std::move(vh)));

1380}

1381

1382template <ndarray::Field T, template <typename> class Array>

1383 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>

1384void svdvals(Array<T>& out, const Array<T>& a) { // singular values only — skips the U/V work entirely

1385 std::size_t m, n;

1386 std::vector<double> A = as_matrix(a, m, n);

1387 SVDc s;

1388 if (m >= n) {

1389 s = svd_golub_reinsch(std::move(A), m, n, /*want_uv=*/false);

1390 } else { // A and Aᵀ share singular values; reduce the tall one

1391 std::vector<double> At(n * m);

1392 for (std::size_t i = 0; i < m; ++i)

1393 for (std::size_t j = 0; j < n; ++j) At[j * m + i] = A[i * n + j];

1394 s = svd_golub_reinsch(std::move(At), n, m, /*want_uv=*/false);

1395 }

1396 copy_into(out, make_vector(std::move(s.w)));

1397}

1398

1399template <ndarray::Field T, template <typename> class Array>

1400 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>

1401void pinv(Array<T>& out, const Array<T>& a) {

1402 std::size_t m, n;

1403 std::vector<double> A = as_matrix(a, m, n);

1404 if (m >= n) {

1405 const SVDc s = svd_golub_reinsch(std::move(A), m, n); // A = U(m×n) diag(w) V(n×n)ᵀ

1406 const double tsh = 0.5 * std::sqrt(double(m + n + 1)) * (s.w.empty() ? 0 : s.w[0]) * 1e-15;

1407 std::vector<double> p(n * m, 0.0); // pinv = V diag(1/w) Uᵀ -> n×m

1408 for (std::size_t i = 0; i < n; ++i)

1409 for (std::size_t j = 0; j < m; ++j) {

1410 double acc = 0.0;

1411 for (std::size_t k = 0; k < n; ++k)

1412 if (s.w[k] > tsh) acc += s.v[i * n + k] * (s.u[j * n + k] / s.w[k]);

1413 p[i * m + j] = acc;

1414 }

1415 copy_into(out, make_matrix(n, m, std::move(p)));

1416 return;

1417 }

1418 // m < n: compute on Aᵀ (n×m, rows>=cols) then transpose the result.

1419 std::vector<double> At(n * m);

1420 for (std::size_t i = 0; i < m; ++i)

1421 for (std::size_t j = 0; j < n; ++j) At[j * m + i] = A[i * n + j];

1422 const SVDc s = svd_golub_reinsch(std::move(At), n, m); // Aᵀ = U(n×m) diag(w) V(m×m)ᵀ

1423 const double tsh = 0.5 * std::sqrt(double(n + m + 1)) * (s.w.empty() ? 0 : s.w[0]) * 1e-15;

1424 std::vector<double> res(n * m, 0.0); // pinv(A) = (V diag(1/w) Uᵀ)ᵀ -> n×m

1425 for (std::size_t i = 0; i < m; ++i)

1426 for (std::size_t j = 0; j < n; ++j) {

1427 double acc = 0.0;

1428 for (std::size_t k = 0; k < m; ++k)

1429 if (s.w[k] > tsh) acc += s.v[i * m + k] * (s.u[j * m + k] / s.w[k]);

1430 res[j * m + i] = acc; // transpose into the n×m result

1431 }

1432 copy_into(out, make_matrix(n, m, std::move(res)));

1433}

1434

1435template <ndarray::Field T, template <typename> class Array>

1436 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>

1437void cond(T& out, const Array<T>& a) {

1438 std::size_t m, n;

1439 std::vector<double> A = as_matrix(a, m, n);

1440 SVDc s;

1441 if (m >= n) {

1442 s = svd_golub_reinsch(std::move(A), m, n, /*want_uv=*/false);

1443 } else {

1444 std::vector<double> At(n * m);

1445 for (std::size_t i = 0; i < m; ++i)

1446 for (std::size_t j = 0; j < n; ++j) At[j * m + i] = A[i * n + j];

1447 s = svd_golub_reinsch(std::move(At), n, m, /*want_uv=*/false);

1448 }

1449 const double wmin = s.w.empty() ? 0 : s.w.back();

1450 out = wmin == 0 ? std::numeric_limits<double>::infinity() : s.w.front() / wmin;

1451}

1452

1453template <ndarray::Field T, template <typename> class Array>

1454 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>

1455void matrix_rank(long long& out, const Array<T>& a) {

1456 std::size_t m, n;

1457 std::vector<double> A = as_matrix(a, m, n);

1458 const bool tr = m < n;

1459 SVDc s;

1460 if (tr) {

1461 std::vector<double> At(n * m);

1462 for (std::size_t i = 0; i < m; ++i)

1463 for (std::size_t j = 0; j < n; ++j) At[j * m + i] = A[i * n + j];

1464 s = svd_golub_reinsch(std::move(At), n, m, /*want_uv=*/false);

1465 } else {

1466 s = svd_golub_reinsch(std::move(A), m, n, /*want_uv=*/false);

1467 }

1468 const double tsh = 0.5 * std::sqrt(double(m + n + 1)) * (s.w.empty() ? 0 : s.w[0]) * 1e-15;

1469 long long r = 0;

1470 for (double w : s.w)

1471 if (w > tsh) ++r;

1472 out = r;

1473}

1474

1475// ---- eigenvalues ----

1476// eigh — eigen-decomposition of a symmetric (real) / Hermitian (complex) matrix. ONE two-layer

1477// KERNEL collapsing the former real and complex overloads: values are always the real spectrum

1478// (Array<real_base_t<T>>), vectors match the input element, and the Hermitian complex path (2n

1479// real embedding) vs the symmetric real path is an `if constexpr` branch on the element.

1480template <ndarray::Field T, template <typename> class Array>

1481 requires HostArray<Array<T>> && ndarray::FloatingPoint<ndarray::real_base_t<T>>

1482void eigh(Array<ndarray::real_base_t<T>>& values, Array<T>& vectors, const Array<T>& a) {

1483 std::size_t n, c;

1484 std::vector<T> A = as_matrix(a, n, c);

1485 require_square(n, c);

1486 std::vector<ndarray::real_base_t<T>> vals;

1487 std::vector<T> vecs;

1488 if constexpr (ndarray::is_complex_v<T>)

1489 hermitian_eig(A, n, vals, vecs, /*want_vectors=*/true);

1490 else

1491 symmetric_eig(std::move(A), n, vals, vecs); // solver owns the copy — no second one

1492 copy_into(values, make_vector(std::move(vals)));

1493 copy_into(vectors, make_matrix(n, n, std::move(vecs)));

1494}

1495// eigvalsh — eigenvalues of a symmetric (real) / Hermitian (complex) matrix; ALWAYS real. Same

1496// unified two-layer kernel shape as eigh, values only.

1497template <ndarray::Field T, template <typename> class Array>

1498 requires HostArray<Array<T>> && ndarray::FloatingPoint<ndarray::real_base_t<T>>

1499void eigvalsh(Array<ndarray::real_base_t<T>>& out, const Array<T>& a) {

1500 std::size_t n, c;

1501 std::vector<T> A = as_matrix(a, n, c);

1502 require_square(n, c);

1503 std::vector<ndarray::real_base_t<T>> vals;

1504 std::vector<T> vecs;

1505 if constexpr (ndarray::is_complex_v<T>)

1506 hermitian_eig(A, n, vals, vecs, /*want_vectors=*/false);

1507 else

1508 symmetric_eig(std::move(A), n, vals, vecs, /*want_vectors=*/false);

1509 copy_into(out, make_vector(std::move(vals)));

1510}

1511// (complex Hermitian eigh/eigvalsh are the SAME two-layer kernels above at T = std::complex<double>.)

1512template void eigvalsh<double, ndarray::basic_ndarray>(NDArray&, const NDArray&);

1513template void eigvalsh<Cplx, ndarray::basic_ndarray>(NDArray&, const CNDArray&);

1514template void eigh<double, ndarray::basic_ndarray>(NDArray&, NDArray&, const NDArray&);

1515template void eigh<Cplx, ndarray::basic_ndarray>(NDArray&, CNDArray&, const CNDArray&);

1516template NDArray eigvalsh<double, ndarray::basic_ndarray>(const NDArray&);

1517template NDArray eigvalsh<Cplx, ndarray::basic_ndarray>(const CNDArray&);

1518template Eig<NDArray> eigh<double, ndarray::basic_ndarray>(const NDArray&);

1519template EighC<NDArray, CNDArray> eigh<Cplx, ndarray::basic_ndarray>(const CNDArray&);

1520

1521// eig — general eigen-decomposition; a real matrix can have a COMPLEX conjugate spectrum, so the

1522// outputs are complex (Array<complex_of_t<T>>) for any real input.

1523template <ndarray::Field T, template <typename> class Array>

1524 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>

1525void eig(Array<ndarray::complex_of_t<T>>& values, Array<ndarray::complex_of_t<T>>& vectors,

1526 const Array<T>& a) {

1527 std::size_t n, c;

1528 std::vector<double> A = as_matrix(a, n, c);

1529 require_square(n, c);

1530 if (is_symmetric(A, n)) { // symmetric -> real spectrum + eigenvectors

1531 // Reuse the matrix we already extracted (eigh(a) would re-extract it — a

1532 // wasted O(n²) copy); the symmetric branch returns here, so moving A is safe.

1533 std::vector<double> rvals, rvecs;

1534 symmetric_eig(std::move(A), n, rvals, rvecs, /*want_vectors=*/true);

1535 copy_into(values, to_complex(make_vector(std::move(rvals))));

1536 copy_into(vectors, to_complex(make_matrix(n, n, std::move(rvecs))));

1537 return;

1538 }

1539 std::vector<Cplx> vals = eigvals_general(A, n);

1540 std::sort(vals.begin(), vals.end(), cgreater);

1541 // Complex eigenvectors via inverse iteration: column k is the eigenvector for vals[k].

1542 std::vector<Cplx> vecs(n * n, Cplx{});

1543 for (std::size_t k = 0; k < n; ++k) {

1544 const std::vector<Cplx> vk = eigvector_inverse_iteration(A, n, vals[k]);

1545 for (std::size_t i = 0; i < n; ++i) vecs[i * n + k] = vk[i];

1546 }

1547 copy_into(values, make_vector(std::move(vals)));

1548 copy_into(vectors, make_matrix(n, n, std::move(vecs)));

1549}

1550// eigvals — the general spectrum (values only); complex for any real input (conjugate pairs).

1551template <ndarray::Field T, template <typename> class Array>

1552 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>

1553void eigvals(Array<ndarray::complex_of_t<T>>& out, const Array<T>& a) {

1554 std::size_t n, c;

1555 std::vector<double> A = as_matrix(a, n, c);

1556 require_square(n, c);

1557 std::vector<Cplx> vals;

1558 if (is_symmetric(A, n)) { // symmetric -> real spectrum, values only

1559 std::vector<double> rvals, vecs;

1560 symmetric_eig(std::move(A), n, rvals, vecs, /*want_vectors=*/false);

1561 vals.assign(rvals.begin(), rvals.end());

1562 } else {

1563 vals = eigvals_general(std::move(A), n);

1564 }

1565 std::sort(vals.begin(), vals.end(), cgreater);

1566 copy_into(out, make_vector(std::move(vals)));

1567}

1568

1569// Explicit instantiations of the remaining kernels AND their (header-inline) allocating fronts —

1570// instantiating the fronts here keeps the exported symbols the library always shipped; the kernel

1571// instantiations provide the exact out-param signatures the buffer-reuse tests call.

1572template void matrix_power<double, ndarray::basic_ndarray>(NDArray&, const NDArray&, long long);

1573template void lstsq<double, ndarray::basic_ndarray>(NDArray&, const NDArray&, const NDArray&);

1574template void cholesky<double, ndarray::basic_ndarray>(NDArray&, const NDArray&);

1575template void svdvals<double, ndarray::basic_ndarray>(NDArray&, const NDArray&);

1576template void pinv<double, ndarray::basic_ndarray>(NDArray&, const NDArray&);

1577template void cond<double, ndarray::basic_ndarray>(double&, const NDArray&);

1578template void matrix_rank<double, ndarray::basic_ndarray>(long long&, const NDArray&);

1579template void qr<double, ndarray::basic_ndarray>(NDArray&, NDArray&, const NDArray&);

1580template void svd<double, ndarray::basic_ndarray>(NDArray&, NDArray&, NDArray&, const NDArray&);

1581template void eig<double, ndarray::basic_ndarray>(CNDArray&, CNDArray&, const NDArray&);

1582template void eigvals<double, ndarray::basic_ndarray>(CNDArray&, const NDArray&);

1583template NDArray matrix_power<double, ndarray::basic_ndarray>(const NDArray&, long long);

1584template NDArray lstsq<double, ndarray::basic_ndarray>(const NDArray&, const NDArray&);

1585template NDArray cholesky<double, ndarray::basic_ndarray>(const NDArray&);

1586template NDArray svdvals<double, ndarray::basic_ndarray>(const NDArray&);

1587template NDArray pinv<double, ndarray::basic_ndarray>(const NDArray&);

1588template double cond<double, ndarray::basic_ndarray>(const NDArray&);

1589template long long matrix_rank<double, ndarray::basic_ndarray>(const NDArray&);

1590template QR<NDArray> qr<double, ndarray::basic_ndarray>(const NDArray&);

1591template SVD<NDArray> svd<double, ndarray::basic_ndarray>(const NDArray&);

1592template EigC<CNDArray> eig<double, ndarray::basic_ndarray>(const NDArray&);

1593template CNDArray eigvals<double, ndarray::basic_ndarray>(const NDArray&);

1594

1595} // namespace cheatah::linalg