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 auto 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 = 0, c = 0;

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

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

330 if (want_uv) {

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

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

333 g = w[i];

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

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

336 if (g != 0.0) {

337 g = 1.0 / g;

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

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

340 double sum = 0.0;

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

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

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

344 }

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

346 } else {

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

348 }

349 Ui[i] += 1.0;

350 }

351 }

352

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

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

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

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

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

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

359 bool flag = true;

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

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

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

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

364 nm = l - 1;

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

366 --l;

367 }

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

369 double c = 0.0, s = 1.0;

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

371 double f = s * rv1[i];

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

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

374 double gg = w[i];

375 double h = pythag(f, gg);

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

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

378 if (want_uv) {

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

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

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

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

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

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

385 }

386 }

387 }

388 }

389 double z = w[k];

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

391 if (z < 0.0) {

392 w[k] = -z;

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

394 }

395 break;

396 }

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

398 double x = w[l];

399 nm = k - 1;

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

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

402 gg = pythag(f, 1.0);

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

404 double c = 1.0, s = 1.0;

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

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

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

408 z = pythag(f, h);

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

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

411 if (want_uv) {

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

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

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

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

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

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

418 }

419 }

420 z = pythag(f, h);

421 w[j] = z;

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

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

424 if (want_uv) {

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

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

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

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

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

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

431 }

432 }

433 }

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

435 }

436 }

437

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

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

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

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

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

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

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

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

446 if (want_uv)

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

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

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

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

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

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

453 }

454 return out;

455}

456

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

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

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

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

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

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

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

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

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

466 values.assign(n, 0.0);

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

468 if (n == 0) return;

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

470

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

472

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

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

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

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

477 double h = 0.0;

478 if (l > 0) {

479 double scale = 0.0;

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

481 if (scale == 0.0) {

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

483 } else {

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

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

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

487 }

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

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

490 e[i] = scale * g;

491 h -= f * g;

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

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

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

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

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

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

498 f = 0.0;

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

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

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

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

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

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

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

506 std::size_t k = 0;

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

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

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

510 }

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

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

513 e[j] = g / h;

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

515 }

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

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

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

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

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

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

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

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

524 }

525 }

526 } else {

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

528 }

529 d[i] = h;

530 }

531 d[0] = 0.0;

532 e[0] = 0.0;

533 if (want_vectors) {

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

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

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

537 double g = 0.0;

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

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

540 }

541 }

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

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

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

545 }

546 } else {

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

548 }

549

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

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

552 e[n - 1] = 0.0;

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

554 int iter = 0;

555 std::size_t m = 0;

556 for (;;) { // sweep until the off-diagonal at l deflates (m == l)

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

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

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

560 }

561 if (m != l) {

562 if (iter++ == 50)

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

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

565 double r = pythag(g, 1.0);

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

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

568 bool zeroed = false;

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

570 double f = s * e[i];

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

572 r = pythag(f, g);

573 e[i + 1] = r;

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

575 s = f / r;

576 c = g / r;

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

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

579 p = s * r;

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

581 g = c * r - b;

582 if (want_vectors)

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

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

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

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

587 }

588 }

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

590 }

591 if (m == l) break;

592 }

593 }

594

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

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

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

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

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

600 if (want_vectors)

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

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

603}

604

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

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

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

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

609 return false;

610 return true;

611}

612

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

632 }

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

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

635 evals.resize(n);

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

637 if (want_vectors) {

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

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

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

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

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

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

644 }

645 }

646 }

647}

648

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

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

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

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

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

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

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

656 std::size_t p = k;

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

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

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

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

661 }

662 piv[k] = p;

663 if (p != k)

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

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

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

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

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

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

670 }

671 }

672 return piv;

673}

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

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

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

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

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

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

680 Cplx s = b[i];

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

682 b[i] = s;

683 }

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

685 Cplx s = b[i];

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

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

688 }

689}

690

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

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

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

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

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

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

697 Cplx lambda) {

698 double scale = 1.0;

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

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

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

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

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

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

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

706 double nrm = 0.0;

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

708 nrm = std::sqrt(nrm);

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

710 };

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

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

713 normalize(v);

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

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

716 normalize(v);

717 }

718 std::size_t mi = 0;

719 double mb = 0.0;

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

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

722 if (m > mb) {

723 mb = m;

724 mi = i;

725 }

726 }

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

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

729 return v;

730}

731

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

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

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

735// eigenvalues are complex.

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

737 // Householder reduction to upper Hessenberg.

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

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

740 double scale = 0.0;

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

742 if (scale == 0.0) continue;

743 double h = 0.0;

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

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

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

747 }

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

749 h -= u[k] * g;

750 u[k] -= g;

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

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

753 double f = 0.0;

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

755 f /= h;

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

757 }

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

759 double f = 0.0;

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

761 f /= h;

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

763 }

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

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

766 }

767

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

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

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

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

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

773 const double eps = 1e-14;

774 int iter = 0;

775 while (hi >= 0) {

776 long long l = hi;

777 while (l > 0) {

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

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

780 --l;

781 }

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

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

784 --hi;

785 iter = 0;

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

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

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

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

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

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

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

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

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

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

796 } else { // complex conjugate pair

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

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

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

800 }

801 hi -= 2;

802 iter = 0;

803 } else { // QR sweep with Wilkinson shift

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

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

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

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

808 cs.clear();

809 sn.clear();

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

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

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

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

814 cs.push_back(c);

815 sn.push_back(s);

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

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

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

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

820 }

821 }

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

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

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

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

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

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

828 }

829 }

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

831 }

832 }

833 return w;

834}

835

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

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

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

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

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

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

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

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

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

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

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

847template <ndarray::Field T>

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

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

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

851 !ndarray::is_contiguous(out))

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

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

854}

855template <ndarray::Field T>

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

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

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

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

860}

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

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

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

864template <ndarray::Field T>

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

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

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

868}

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

870template <ndarray::Field T>

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

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

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

874}

875

876} // namespace

877

878// ================= public routines =================

879

880// ---- products ----

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

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

883// and allocate nothing but the result.

884//

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

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

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

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

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

890namespace {

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

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

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

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

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

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

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

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

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

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

901 });

902}

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

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

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

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

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

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

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

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

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

912}

913} // namespace

914

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

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

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

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

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

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

921 requires HostArray<Array<T>>

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

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

924 requires HostArray<Array<T>>

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

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

927 requires HostArray<Array<T>>

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

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

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

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

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

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

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

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

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

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

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

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

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

941

942namespace {

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

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

945template <ndarray::Field T>

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

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

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

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

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

951 }

952}

953} // namespace

954

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

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

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

958 requires HostArray<Array<T>>

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

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

961 reject_alias(out, a);

962 reject_alias(out, b);

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

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

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

966}

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

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

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

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

971

972namespace {

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

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

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

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

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

978template <ndarray::Field T>

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

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

981 std::size_t i = 0;

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

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

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

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

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

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

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

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

990 const T bkj = bk[j];

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

992 }

993 }

994 }

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

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

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

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

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

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

1001 }

1002 }

1003}

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

1005template <ndarray::Field T>

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

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

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

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

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

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

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

1013}

1014} // namespace

1015

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

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

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

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

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

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

1022 requires HostArray<Array<T>>

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

1024 if (a.ndim() == 3 || b.ndim() == 3) {

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

1026 // Validated HERE, not left to the caller: this kernel is public (declared in backend.hpp,

1027 // documented and explicitly instantiated), so it is reachable without the allocating

1028 // front. The 2-D control flow below is untouched.

1029 check_matmul_batched<T, Array>(a, b);

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

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

1032 reject_alias(out, a);

1033 reject_alias(out, b);

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

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

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

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

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

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

1040 return;

1041 }

1042 std::size_t ar = 0, ac = 0, bc = 0;

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

1044 reject_alias(out, a);

1045 reject_alias(out, b);

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

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

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

1049}

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

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

1052// attributes their bodies to this TU.

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

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

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

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

1057

1058namespace {

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

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

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

1062template <ndarray::Field T>

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

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

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

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

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

1068 }

1069}

1070} // namespace

1071

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

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

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

1075 requires HostArray<Array<T>>

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

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

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

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

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

1081 std::vector<T> sa;

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

1083}

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

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

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

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

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

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

1090

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

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

1093void matrix_power(Array<T>& out, const Array<T>& a, long long n) {

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

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

1096 require_square(r, c);

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

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

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

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

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

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

1103 NDArray base = (n < 0) ? inv(a) : a;

1104 long long e = n < 0 ? -n : n;

1105 while (e > 0) {

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

1107 base = matmul(base, base);

1108 e >>= 1;

1109 }

1110 copy_into(out, acc);

1111}

1112

1113namespace {

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

1115// whole of B.

1116template <ndarray::Field T>

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

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

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

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

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

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

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

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

1125}

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

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

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

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

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

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

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

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

1134}

1135} // namespace

1136

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

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

1139 requires HostArray<Array<T>>

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

1141 std::size_t ar = 0, ac = 0, br = 0, bc = 0;

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

1143 reject_alias(out, a);

1144 reject_alias(out, b);

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

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

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

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

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

1150}

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

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

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

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

1155

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

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

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

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

1160 requires HostArray<Array<T>>

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

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

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

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

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

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

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

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

1169 });

1170}

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

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

1173

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

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

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

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

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

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

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

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

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

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

1184 return;

1185 }

1186 std::vector<T> scratch;

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

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

1189}

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

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

1192

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

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

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

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

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

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

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

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

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

1202 const LU lu = lu_prepare(a);

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

1204 std::size_t bn = 0;

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

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

1207 lu_solve(lu, x);

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

1209}

1210

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

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

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

1214 const LU lu = lu_prepare(a);

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

1216 T d = lu.sign;

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

1218 out = d;

1219}

1220

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

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

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

1224 const LU lu = lu_prepare(a);

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

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

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

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

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

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

1231 }

1232 out = {sign, logabs};

1233}

1234

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

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

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

1238 const LU lu = lu_prepare(a);

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

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

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

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

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

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

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

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

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

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

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

1250 for (std::size_t col = 0; col < n; ++col) std::swap(X[k * n + col], X[lu.piv[k] * n + col]);

1251 for (std::size_t i = 0; i < n; ++i) // forward: unit-lower L·Y = P

1252 for (std::size_t j = 0; j < i; ++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 for (std::size_t i = n; i-- > 0;) { // back: upper U·X = Y

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

1258 const double f = M[i * n + j];

1259 for (std::size_t col = 0; col < n; ++col) X[i * n + col] -= f * X[j * n + col];

1260 }

1261 const double d = M[i * n + i];

1262 for (std::size_t col = 0; col < n; ++col) X[i * n + col] /= d;

1263 }

1264 copy_into(out, make_matrix(n, n, std::move(X)));

1265}

1266// Explicit instantiations of the LU family: the host kernels AND the (now header-inline)

1267// allocating fronts — instantiating the fronts here keeps the exported symbols the library ships.

1268template void solve<double, ndarray::basic_ndarray>(NDArray&, const NDArray&, const NDArray&);

1269template void det<double, ndarray::basic_ndarray>(double&, const NDArray&);

1270template void slogdet<double, ndarray::basic_ndarray>(SLogDet&, const NDArray&);

1271template void inv<double, ndarray::basic_ndarray>(NDArray&, const NDArray&);

1272template NDArray solve<double, ndarray::basic_ndarray>(const NDArray&, const NDArray&);

1273template double det<double, ndarray::basic_ndarray>(const NDArray&);

1274template SLogDet slogdet<double, ndarray::basic_ndarray>(const NDArray&);

1275template NDArray inv<double, ndarray::basic_ndarray>(const NDArray&);

1276

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

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

1279void lstsq(Array<T>& out, const Array<T>& a, const Array<T>& b) { // min ‖Ax−b‖ via the pseudo-inverse

1280 Array<T> p = Array<T>::uninitialized({a.shape()[1], a.shape()[0]});

1281 pinv(p, a);

1282 matmul(out, p, b);

1283}

1284

1285// ---- Cholesky ----

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

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

1288void cholesky(Array<T>& out, const Array<T>& a) {

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

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

1291 require_square(n, c);

1292 std::vector<double> scratch;

1293 const double* A = contig(a, scratch); // read-only — zero-copy when contiguous

1294 std::vector<double> L(n * n, 0.0);

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

1296 const double* Li = &L[i * n];

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

1298 // s = A[i][j] − (row i · row j over k<j): four accumulators so the O(n³)

1299 // inner dot vectorizes instead of serializing on FP-sub latency.

1300 const double* Lj = &L[j * n];

1301 double d0 = 0, d1 = 0, d2 = 0, d3 = 0;

1302 std::size_t k = 0;

1303 for (; k + 4 <= j; k += 4) {

1304 d0 += Li[k] * Lj[k]; d1 += Li[k + 1] * Lj[k + 1];

1305 d2 += Li[k + 2] * Lj[k + 2]; d3 += Li[k + 3] * Lj[k + 3];

1306 }

1307 double s = A[i * n + j] - ((d0 + d1) + (d2 + d3));

1308 for (; k < j; ++k) s -= Li[k] * Lj[k];

1309 if (i == j) {

1310 if (s <= 0.0) throw std::runtime_error("linalg: matrix is not positive-definite");

1311 L[i * n + i] = std::sqrt(s);

1312 } else {

1313 L[i * n + j] = s / L[j * n + j];

1314 }

1315 }

1316 }

1317 copy_into(out, make_matrix(n, n, std::move(L)));

1318}

1319

1320// ---- Householder QR (reduced: Q is m×n, R is n×n) ----

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

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

1323void qr(Array<T>& q, Array<T>& r, const Array<T>& a) {

1324 std::size_t m = 0, n = 0;

1325 std::vector<double> A = as_matrix(a, m, n);

1326 if (m < n) throw std::runtime_error("linalg: qr requires rows >= cols");

1327 // Work on the TRANSPOSE At (n×m, row-major). A Householder QR repeatedly reads and

1328 // updates COLUMNS of A, which stride by n in row-major and don't vectorize (the

1329 // original cost ~3× Eigen); as ROWS of At those same operations are contiguous, and

1330 // the reductions are multi-accumulated like ddot.

1331 std::vector<double> At(n * m);

1332 for (std::size_t i = 0; i < m; ++i)

1333 for (std::size_t j = 0; j < n; ++j) At[j * m + i] = A[i * n + j];

1334 std::vector<double> Q(m * m, 0.0);

1335 for (std::size_t i = 0; i < m; ++i) Q[i * m + i] = 1.0;

1336 std::vector<double> u(m); // Householder vector, reused per column (entries < k unused)

1337 // Reflect: s = u · row over [k, m) via the shared multi-accumulator reduction, then

1338 // row -= (2 s / ‖u‖²) u — contiguous.

1339 auto reflect = [&u](double* row, std::size_t k, std::size_t m, double inv) {

1340 auto s = ndarray::detail::reduce_lanes<double>(

1341 m - k, [&u, row, k](std::size_t i) { return u[k + i] * row[k + i]; });

1342 s *= inv;

1343 for (std::size_t i = k; i < m; ++i) row[i] -= s * u[i];

1344 };

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

1346 double* Atk = &At[k * m]; // column k of A == row k of At

1347 double nrm = 0.0;

1348 for (std::size_t i = k; i < m; ++i) nrm += Atk[i] * Atk[i];

1349 nrm = std::sqrt(nrm);

1350 if (nrm == 0.0) continue;

1351 const double alpha = Atk[k] >= 0 ? -nrm : nrm;

1352 for (std::size_t i = k; i < m; ++i) u[i] = Atk[i];

1353 u[k] -= alpha;

1354 double unorm2 = 0.0;

1355 for (std::size_t i = k; i < m; ++i) unorm2 += u[i] * u[i];

1356 if (unorm2 == 0.0) continue;

1357 const double inv = 2.0 / unorm2;

1358 for (std::size_t j = k; j < n; ++j) reflect(&At[j * m], k, m, inv); // A's cols j≥k

1359 for (std::size_t j = 0; j < m; ++j) reflect(&Q[j * m], k, m, inv); // Q = Q · Hₖ

1360 }

1361 std::vector<double> Qr(m * n), Rr(n * n, 0.0); // reduced

1362 for (std::size_t i = 0; i < m; ++i)

1363 for (std::size_t j = 0; j < n; ++j) Qr[i * n + j] = Q[i * m + j];

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

1365 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]

1366 copy_into(q, make_matrix(m, n, std::move(Qr)));

1367 copy_into(r, make_matrix(n, n, std::move(Rr)));

1368}

1369

1370// ---- SVD and its derived quantities ----

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

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

1373void svd(Array<T>& u, Array<T>& sv, Array<T>& vhh, const Array<T>& a) {

1374 std::size_t m = 0, n = 0;

1375 std::vector<double> A = as_matrix(a, m, n);

1376 if (m < n) throw std::runtime_error("linalg: svd requires rows >= cols (transpose otherwise)");

1377 const SVDc s = svd_golub_reinsch(std::move(A), m, n);

1378 // vh = Vᵀ

1379 std::vector<double> vh(n * n);

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

1381 for (std::size_t j = 0; j < n; ++j) vh[i * n + j] = s.v[j * n + i];

1382 copy_into(u, make_matrix(m, n, s.u));

1383 copy_into(sv, make_vector(s.w));

1384 copy_into(vhh, make_matrix(n, n, std::move(vh)));

1385}

1386

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

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

1389void svdvals(Array<T>& out, const Array<T>& a) { // singular values only — skips the U/V work entirely

1390 std::size_t m = 0, n = 0;

1391 std::vector<double> A = as_matrix(a, m, n);

1392 SVDc s;

1393 if (m >= n) {

1394 s = svd_golub_reinsch(std::move(A), m, n, /*want_uv=*/false);

1395 } else { // A and Aᵀ share singular values; reduce the tall one

1396 std::vector<double> At(n * m);

1397 for (std::size_t i = 0; i < m; ++i)

1398 for (std::size_t j = 0; j < n; ++j) At[j * m + i] = A[i * n + j];

1399 s = svd_golub_reinsch(std::move(At), n, m, /*want_uv=*/false);

1400 }

1401 copy_into(out, make_vector(std::move(s.w)));

1402}

1403

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

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

1406void pinv(Array<T>& out, const Array<T>& a) {

1407 std::size_t m = 0, n = 0;

1408 std::vector<double> A = as_matrix(a, m, n);

1409 if (m >= n) {

1410 const SVDc s = svd_golub_reinsch(std::move(A), m, n); // A = U(m×n) diag(w) V(n×n)ᵀ

1411 const double tsh = 0.5 * std::sqrt(double(m + n + 1)) * (s.w.empty() ? 0 : s.w[0]) * 1e-15;

1412 std::vector<double> p(n * m, 0.0); // pinv = V diag(1/w) Uᵀ -> n×m

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

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

1415 double acc = 0.0;

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

1417 if (s.w[k] > tsh) acc += s.v[i * n + k] * (s.u[j * n + k] / s.w[k]);

1418 p[i * m + j] = acc;

1419 }

1420 copy_into(out, make_matrix(n, m, std::move(p)));

1421 return;

1422 }

1423 // m < n: compute on Aᵀ (n×m, rows>=cols) then transpose the result.

1424 std::vector<double> At(n * m);

1425 for (std::size_t i = 0; i < m; ++i)

1426 for (std::size_t j = 0; j < n; ++j) At[j * m + i] = A[i * n + j];

1427 const SVDc s = svd_golub_reinsch(std::move(At), n, m); // Aᵀ = U(n×m) diag(w) V(m×m)ᵀ

1428 const double tsh = 0.5 * std::sqrt(double(n + m + 1)) * (s.w.empty() ? 0 : s.w[0]) * 1e-15;

1429 std::vector<double> res(n * m, 0.0); // pinv(A) = (V diag(1/w) Uᵀ)ᵀ -> n×m

1430 for (std::size_t i = 0; i < m; ++i)

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

1432 double acc = 0.0;

1433 for (std::size_t k = 0; k < m; ++k)

1434 if (s.w[k] > tsh) acc += s.v[i * m + k] * (s.u[j * m + k] / s.w[k]);

1435 res[j * m + i] = acc; // transpose into the n×m result

1436 }

1437 copy_into(out, make_matrix(n, m, std::move(res)));

1438}

1439

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

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

1442void cond(T& out, const Array<T>& a) {

1443 std::size_t m = 0, n = 0;

1444 std::vector<double> A = as_matrix(a, m, n);

1445 SVDc s;

1446 if (m >= n) {

1447 s = svd_golub_reinsch(std::move(A), m, n, /*want_uv=*/false);

1448 } else {

1449 std::vector<double> At(n * m);

1450 for (std::size_t i = 0; i < m; ++i)

1451 for (std::size_t j = 0; j < n; ++j) At[j * m + i] = A[i * n + j];

1452 s = svd_golub_reinsch(std::move(At), n, m, /*want_uv=*/false);

1453 }

1454 const double wmin = s.w.empty() ? 0 : s.w.back();

1455 out = wmin == 0 ? std::numeric_limits<double>::infinity() : s.w.front() / wmin;

1456}

1457

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

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

1460void matrix_rank(long long& out, const Array<T>& a) {

1461 std::size_t m = 0, n = 0;

1462 std::vector<double> A = as_matrix(a, m, n);

1463 const bool tr = m < n;

1464 SVDc s;

1465 if (tr) {

1466 std::vector<double> At(n * m);

1467 for (std::size_t i = 0; i < m; ++i)

1468 for (std::size_t j = 0; j < n; ++j) At[j * m + i] = A[i * n + j];

1469 s = svd_golub_reinsch(std::move(At), n, m, /*want_uv=*/false);

1470 } else {

1471 s = svd_golub_reinsch(std::move(A), m, n, /*want_uv=*/false);

1472 }

1473 const double tsh = 0.5 * std::sqrt(double(m + n + 1)) * (s.w.empty() ? 0 : s.w[0]) * 1e-15;

1474 long long r = 0;

1475 for (double w : s.w)

1476 if (w > tsh) ++r;

1477 out = r;

1478}

1479

1480// ---- eigenvalues ----

1481// eigh — eigen-decomposition of a symmetric (real) / Hermitian (complex) matrix. ONE two-layer

1482// KERNEL collapsing the former real and complex overloads: values are always the real spectrum

1483// (Array<real_base_t<T>>), vectors match the input element, and the Hermitian complex path (2n

1484// real embedding) vs the symmetric real path is an `if constexpr` branch on the element.

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

1486 requires HostArray<Array<T>> && ndarray::FloatingPoint<ndarray::real_base_t<T>>

1487void eigh(Array<ndarray::real_base_t<T>>& values, Array<T>& vectors, const Array<T>& a) {

1488 std::size_t n = 0, c = 0;

1489 std::vector<T> A = as_matrix(a, n, c);

1490 require_square(n, c);

1491 std::vector<ndarray::real_base_t<T>> vals;

1492 std::vector<T> vecs;

1493 if constexpr (ndarray::is_complex_v<T>)

1494 hermitian_eig(A, n, vals, vecs, /*want_vectors=*/true);

1495 else

1496 symmetric_eig(std::move(A), n, vals, vecs); // solver owns the copy — no second one

1497 copy_into(values, make_vector(std::move(vals)));

1498 copy_into(vectors, make_matrix(n, n, std::move(vecs)));

1499}

1500// eigvalsh — eigenvalues of a symmetric (real) / Hermitian (complex) matrix; ALWAYS real. Same

1501// unified two-layer kernel shape as eigh, values only.

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

1503 requires HostArray<Array<T>> && ndarray::FloatingPoint<ndarray::real_base_t<T>>

1504void eigvalsh(Array<ndarray::real_base_t<T>>& out, const Array<T>& a) {

1505 std::size_t n = 0, c = 0;

1506 std::vector<T> A = as_matrix(a, n, c);

1507 require_square(n, c);

1508 std::vector<ndarray::real_base_t<T>> vals;

1509 std::vector<T> vecs;

1510 if constexpr (ndarray::is_complex_v<T>)

1511 hermitian_eig(A, n, vals, vecs, /*want_vectors=*/false);

1512 else

1513 symmetric_eig(std::move(A), n, vals, vecs, /*want_vectors=*/false);

1514 copy_into(out, make_vector(std::move(vals)));

1515}

1516// (complex Hermitian eigh/eigvalsh are the SAME two-layer kernels above at T = std::complex<double>.)

1517template void eigvalsh<double, ndarray::basic_ndarray>(NDArray&, const NDArray&);

1518template void eigvalsh<Cplx, ndarray::basic_ndarray>(NDArray&, const CNDArray&);

1519template void eigh<double, ndarray::basic_ndarray>(NDArray&, NDArray&, const NDArray&);

1520template void eigh<Cplx, ndarray::basic_ndarray>(NDArray&, CNDArray&, const CNDArray&);

1521template NDArray eigvalsh<double, ndarray::basic_ndarray>(const NDArray&);

1522template NDArray eigvalsh<Cplx, ndarray::basic_ndarray>(const CNDArray&);

1523template Eig<NDArray> eigh<double, ndarray::basic_ndarray>(const NDArray&);

1524template EighC<NDArray, CNDArray> eigh<Cplx, ndarray::basic_ndarray>(const CNDArray&);

1525

1526// eig — general eigen-decomposition; a real matrix can have a COMPLEX conjugate spectrum, so the

1527// outputs are complex (Array<complex_of_t<T>>) for any real input.

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

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

1530void eig(Array<ndarray::complex_of_t<T>>& values, Array<ndarray::complex_of_t<T>>& vectors,

1531 const Array<T>& a) {

1532 std::size_t n = 0, c = 0;

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

1534 require_square(n, c);

1535 if (is_symmetric(A, n)) { // symmetric -> real spectrum + eigenvectors

1536 // Reuse the matrix we already extracted (eigh(a) would re-extract it — a

1537 // wasted O(n²) copy); the symmetric branch returns here, so moving A is safe.

1538 std::vector<double> rvals, rvecs;

1539 symmetric_eig(std::move(A), n, rvals, rvecs, /*want_vectors=*/true);

1540 copy_into(values, to_complex(make_vector(std::move(rvals))));

1541 copy_into(vectors, to_complex(make_matrix(n, n, std::move(rvecs))));

1542 return;

1543 }

1544 std::vector<Cplx> vals = eigvals_general(A, n);

1545 std::sort(vals.begin(), vals.end(), cgreater);

1546 // Complex eigenvectors via inverse iteration: column k is the eigenvector for vals[k].

1547 std::vector<Cplx> vecs(n * n, Cplx{});

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

1549 const std::vector<Cplx> vk = eigvector_inverse_iteration(A, n, vals[k]);

1550 for (std::size_t i = 0; i < n; ++i) vecs[i * n + k] = vk[i];

1551 }

1552 copy_into(values, make_vector(std::move(vals)));

1553 copy_into(vectors, make_matrix(n, n, std::move(vecs)));

1554}

1555// eigvals — the general spectrum (values only); complex for any real input (conjugate pairs).

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

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

1558void eigvals(Array<ndarray::complex_of_t<T>>& out, const Array<T>& a) {

1559 std::size_t n = 0, c = 0;

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

1561 require_square(n, c);

1562 std::vector<Cplx> vals;

1563 if (is_symmetric(A, n)) { // symmetric -> real spectrum, values only

1564 std::vector<double> rvals, vecs;

1565 symmetric_eig(std::move(A), n, rvals, vecs, /*want_vectors=*/false);

1566 vals.assign(rvals.begin(), rvals.end());

1567 } else {

1568 vals = eigvals_general(std::move(A), n);

1569 }

1570 std::sort(vals.begin(), vals.end(), cgreater);

1571 copy_into(out, make_vector(std::move(vals)));

1572}

1573

1574// Explicit instantiations of the remaining kernels AND their (header-inline) allocating fronts —

1575// instantiating the fronts here keeps the exported symbols the library always shipped; the kernel

1576// instantiations provide the exact out-param signatures the buffer-reuse tests call.

1577template void matrix_power<double, ndarray::basic_ndarray>(NDArray&, const NDArray&, long long);

1578template void lstsq<double, ndarray::basic_ndarray>(NDArray&, const NDArray&, const NDArray&);

1579template void cholesky<double, ndarray::basic_ndarray>(NDArray&, const NDArray&);

1580template void svdvals<double, ndarray::basic_ndarray>(NDArray&, const NDArray&);

1581template void pinv<double, ndarray::basic_ndarray>(NDArray&, const NDArray&);

1582template void cond<double, ndarray::basic_ndarray>(double&, const NDArray&);

1583template void matrix_rank<double, ndarray::basic_ndarray>(long long&, const NDArray&);

1584template void qr<double, ndarray::basic_ndarray>(NDArray&, NDArray&, const NDArray&);

1585template void svd<double, ndarray::basic_ndarray>(NDArray&, NDArray&, NDArray&, const NDArray&);

1586template void eig<double, ndarray::basic_ndarray>(CNDArray&, CNDArray&, const NDArray&);

1587template void eigvals<double, ndarray::basic_ndarray>(CNDArray&, const NDArray&);

1588template NDArray matrix_power<double, ndarray::basic_ndarray>(const NDArray&, long long);

1589template NDArray lstsq<double, ndarray::basic_ndarray>(const NDArray&, const NDArray&);

1590template NDArray cholesky<double, ndarray::basic_ndarray>(const NDArray&);

1591template NDArray svdvals<double, ndarray::basic_ndarray>(const NDArray&);

1592template NDArray pinv<double, ndarray::basic_ndarray>(const NDArray&);

1593template double cond<double, ndarray::basic_ndarray>(const NDArray&);

1594template long long matrix_rank<double, ndarray::basic_ndarray>(const NDArray&);

1595template QR<NDArray> qr<double, ndarray::basic_ndarray>(const NDArray&);

1596template SVD<NDArray> svd<double, ndarray::basic_ndarray>(const NDArray&);

1597template EigC<CNDArray> eig<double, ndarray::basic_ndarray>(const NDArray&);

1598template CNDArray eigvals<double, ndarray::basic_ndarray>(const NDArray&);

1599

1600} // namespace cheatah::linalg