From a5ca13e44165ee8a3420a57b95bcf84437a81dce Mon Sep 17 00:00:00 2001 From: y-jan137 Date: Sat, 14 Mar 2026 22:10:21 +0300 Subject: Add LU and experiment --- include/lu.hpp | 43 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) create mode 100644 include/lu.hpp (limited to 'include/lu.hpp') diff --git a/include/lu.hpp b/include/lu.hpp new file mode 100644 index 0000000..6717733 --- /dev/null +++ b/include/lu.hpp @@ -0,0 +1,43 @@ +#pragma once + +#include +#include + +#include "matrix.hpp" +#include "vector.hpp" + +namespace linalg { + +// Result of LU factorization with partial pivoting. +// +// The factorization satisfies PA = LU, where: +// P is the permutation matrix encoded by `perm` +// L is unit lower triangular (L[i][i] == 1) +// U is upper triangular +// +// `perm[i]` = index of the original row that ended up at position i. +// Applying P to a vector b means: (Pb)[i] = b[perm[i]]. +// +// `sign` is the sign of the permutation: +1 if an even number of row +// swaps were made, -1 if odd. Useful for computing det(A) = sign * prod(diag(U)). +struct LUResult { + Matrix L; + Matrix U; + std::vector perm; + int sign; +}; + +// Compute the LU factorization of A with partial pivoting. +// +// Throws DimensionMismatchError if A is not square. +// Throws SingularMatrixError if A is (numerically) singular, i.e. any +// pivot is smaller in magnitude than `singular_tolerance`. +LUResult lu_factor(const Matrix& A, double singular_tolerance = 1e-12); + +// Solve Ax = b given a precomputed LU factorization. +// +// Applies the stored permutation, then forward / backward substitution. +// Throws DimensionMismatchError if b.size() != lu.L.rows(). +Vector lu_solve(const LUResult& lu, const Vector& b); + +} // namespace linalg -- cgit v1.2.3