aboutsummaryrefslogtreecommitdiff
path: root/include
diff options
context:
space:
mode:
Diffstat (limited to 'include')
-rw-r--r--include/linalg_error.hpp31
-rw-r--r--include/lu.hpp22
-rw-r--r--include/matrix.hpp48
-rw-r--r--include/norms.hpp9
-rw-r--r--include/qr.hpp39
-rw-r--r--include/qr_iteration.hpp165
-rw-r--r--include/triangular_solve.hpp22
-rw-r--r--include/vector.hpp47
8 files changed, 0 insertions, 383 deletions
diff --git a/include/linalg_error.hpp b/include/linalg_error.hpp
deleted file mode 100644
index ac9b456..0000000
--- a/include/linalg_error.hpp
+++ /dev/null
@@ -1,31 +0,0 @@
-#pragma once
-
-#include <stdexcept>
-#include <string>
-
-namespace linalg {
-
-class LinAlgError : public std::runtime_error {
-public:
- using std::runtime_error::runtime_error;
-};
-
-class DimensionMismatchError : public LinAlgError {
-public:
- explicit DimensionMismatchError(const std::string& message)
- : LinAlgError(message) {}
-};
-
-class SingularMatrixError : public LinAlgError {
-public:
- explicit SingularMatrixError(const std::string& message)
- : LinAlgError(message) {}
-};
-
-class NonConvergenceError : public LinAlgError {
-public:
- explicit NonConvergenceError(const std::string& message)
- : LinAlgError(message) {}
-};
-
-} // namespace linalg
diff --git a/include/lu.hpp b/include/lu.hpp
deleted file mode 100644
index e7fe939..0000000
--- a/include/lu.hpp
+++ /dev/null
@@ -1,22 +0,0 @@
-#pragma once
-
-#include <cstddef>
-#include <vector>
-
-#include "matrix.hpp"
-#include "vector.hpp"
-
-namespace linalg {
-
-struct LUResult {
- Matrix L;
- Matrix U;
- std::vector<std::size_t> perm;
- int sign;
-};
-
-LUResult lu_factor(const Matrix& A, double singular_tolerance = 1e-12);
-
-Vector lu_solve(const LUResult& lu, const Vector& b);
-
-} // namespace linalg
diff --git a/include/matrix.hpp b/include/matrix.hpp
deleted file mode 100644
index d78f814..0000000
--- a/include/matrix.hpp
+++ /dev/null
@@ -1,48 +0,0 @@
-#pragma once
-
-#include <cstddef>
-#include <initializer_list>
-#include <vector>
-#include "vector.hpp"
-
-namespace linalg {
-
-class Matrix {
-public:
- Matrix() = default;
- Matrix(std::size_t rows, std::size_t cols);
- Matrix(std::size_t rows, std::size_t cols, double value);
- Matrix(std::initializer_list<std::initializer_list<double>> values);
-
- [[nodiscard]] std::size_t rows() const noexcept;
- [[nodiscard]] std::size_t cols() const noexcept;
- [[nodiscard]] bool empty() const noexcept;
-
- double& operator()(std::size_t i, std::size_t j);
- const double& operator()(std::size_t i, std::size_t j) const;
-
- void fill(double value);
-
- double* data() noexcept;
- const double* data() const noexcept;
-
- static Matrix identity(std::size_t n);
- static Matrix zeros(std::size_t rows, std::size_t cols);
-
-private:
- [[nodiscard]] std::size_t index(std::size_t i, std::size_t j) const;
- void check_bounds(std::size_t i, std::size_t j) const;
-
- std::size_t rows_ = 0;
- std::size_t cols_ = 0;
- std::vector<double> data_;
-};
-
-Matrix transpose(const Matrix& matrix);
-Matrix operator+(const Matrix& lhs, const Matrix& rhs);
-Matrix operator-(const Matrix& lhs, const Matrix& rhs);
-Vector operator*(const Matrix& matrix, const Vector& vector);
-Matrix operator*(const Matrix& lhs, const Matrix& rhs);
-
-} // namespace linalg
-
diff --git a/include/norms.hpp b/include/norms.hpp
deleted file mode 100644
index 85fad41..0000000
--- a/include/norms.hpp
+++ /dev/null
@@ -1,9 +0,0 @@
-#pragma once
-
-#include "vector.hpp"
-
-namespace linalg {
-
-double norm2(const Vector& vector);
-
-}
diff --git a/include/qr.hpp b/include/qr.hpp
deleted file mode 100644
index c9e1638..0000000
--- a/include/qr.hpp
+++ /dev/null
@@ -1,39 +0,0 @@
-#pragma once
-
-#include "matrix.hpp"
-#include "vector.hpp"
-
-namespace linalg {
-
-struct QRResult {
- Matrix Q;
- Matrix R;
-};
-
-// Classical Gram-Schmidt.
-// Mathematically natural but numerically fragile: orthogonality of Q
-// degrades rapidly on ill-conditioned inputs.
-// Provided for comparison — prefer modified_gs or householder in practice.
-//
-// Throws DimensionMismatchError if rows < cols.
-// Throws SingularMatrixError if a column is (nearly) linearly dependent.
-QRResult qr_classical_gs(const Matrix& A, double zero_tolerance = 1e-14);
-
-// Modified Gram-Schmidt.
-// Subtracts each projection immediately on the running vector rather than
-// on the original column. Algebraically equivalent to classical GS but
-// numerically much better — round-off stays local instead of accumulating.
-//
-// Same exceptions as classical GS.
-QRResult qr_modified_gs(const Matrix& A, double zero_tolerance = 1e-14);
-
-// Householder QR.
-// Applies a sequence of orthogonal reflections to zero out below-diagonal
-// entries column by column. Backward-stable and the standard choice for
-// dense QR. Works correctly on rank-deficient matrices (zero pivots
-// produce zero diagonal entries in R without throwing).
-//
-// Throws DimensionMismatchError if rows < cols.
-QRResult qr_householder(const Matrix& A);
-
-} // namespace linalg
diff --git a/include/qr_iteration.hpp b/include/qr_iteration.hpp
deleted file mode 100644
index ccf48cd..0000000
--- a/include/qr_iteration.hpp
+++ /dev/null
@@ -1,165 +0,0 @@
-#pragma once
-
-// QR iteration for eigenvalue computation
-//
-// Refs:
-// Trefethen & Bau, "Numerical Linear Algebra" (T&B)
-// Lecture 25 — Eigenvalue algorithms
-// Lecture 26 — Schur factorisation
-// Lecture 28 — The QR algorithm (unshifted)
-// Lecture 29 — The QR algorithm with shifts
-// Golub & Van Loan, "Matrix Computations" 4th ed. (GVL)
-// §7.3 — The Unshifted QR Algorithm
-// §7.4 — The Shifted QR Algorithm
-// §7.4.2 — Wilkinson shift
-
-#include <vector>
-
-#include "linalg_error.hpp"
-#include "matrix.hpp"
-#include "vector.hpp"
-
-namespace linalg {
-
-// ---------------------------------------------------------------------------
-// Options
-// ---------------------------------------------------------------------------
-
-// All defaults are consistent with the recommendations in T&B Lecture 28.
-struct QRIterationOptions {
- // Convergence threshold. Iteration halts once the Frobenius norm of the
- // strict lower triangle of A_k falls below this value.
- // Ref: T&B §28; GVL §7.3.
- double tolerance = 1e-10;
-
- int max_iterations = 1000;
-
- // When true, the Frobenius norm of the strict lower triangle is recorded
- // after every QR step and returned in QRIterationResult::convergence_history.
- bool track_convergence = false;
-};
-
-
-struct QRIterationResult {
- // For symmetric inputs all imaginary parts are zero.
- // Complex-conjugate pairs from 2×2 Schur blocks appear as +/-imag entries.
- Vector eigenvalues_real;
- Vector eigenvalues_imag;
-
- int iterations = 0;
-
- // Populated only when QRIterationOptions::track_convergence is true.
- // Entry k is ||lower(A_k)||_F after the k-th QR step.
- // std::vector is used here because linalg::Vector has no push_back;
- // convergence_history is a plain time-series container, not a math object.
- std::vector<double> convergence_history;
-};
-
-// Algorithm (T&B Algorithm 28.1):
-//
-// A_0 = A
-// for k = 1, 2, ...:
-// factor A_{k-1} = Q_k R_k (Householder QR)
-// set A_k = R_k Q_k (orthogonal similarity: preserves eigenvalues)
-//
-// The iterates A_k converge to the real Schur form of A: a quasi-upper-
-// triangular matrix with 1×1 blocks (real eigenvalue) and 2×2 blocks
-// (complex-conjugate pair) on the diagonal.
-//
-// Convergence rate: linear. Per-step factor ~= |lambda_{j+1} / lambda_j|
-// for the off-diagonal entries linking eigenvalue clusters j and j+1.
-// (T&B Lecture 28, Theorem 28.2)
-//
-// Throws DimensionMismatchError if A is not square.
-// Throws NonConvergenceError if convergence is not achieved within
-// opts.max_iterations steps.
-[[nodiscard]] QRIterationResult eigenvalues_unshifted(const Matrix& A,
- QRIterationOptions opts = {});
-
-// ---------------------------------------------------------------------------
-// Wilkinson-shifted QR iteration
-// ---------------------------------------------------------------------------
-//
-// Same outer loop as Stage 1, but each step applies a shift σ chosen as
-// the eigenvalue of the bottom-right 2×2 block of A_{k-1} that is closest
-// to the (n,n) entry, then unshifts after the QR step:
-//
-// factor (A_{k-1} - σI) = Q_k R_k
-// set A_k = R_k Q_k + σI
-//
-// The Wilkinson shift (T&B Lecture 29; GVL §7.4.2):
-// Given the bottom-right 2×2 block | a b |
-// | b c |
-// δ = (a - c) / 2
-// σ = c - sign(δ) * b² / (|δ| + sqrt(δ² + b²))
-// equivalently: the eigenvalue of the block closer to c.
-//
-// Convergence rate: typically cubic near a simple eigenvalue.
-// (T&B Lecture 29; GVL §7.5.1)
-//
-// Same exceptions as eigenvalues_unshifted.
-[[nodiscard]] QRIterationResult eigenvalues_shifted(const Matrix& A,
- QRIterationOptions opts = {});
-
-// ---------------------------------------------------------------------------
-// Stage 3: Hessenberg reduction algorithm
-// ---------------------------------------------------------------------------
-
-// Givens rotation G acting on rows/columns i and i+1:
-//
-// G = | c s | chosen so that G * [x; y]^T = [r; 0]^T
-// | -s c | with c = x/r, s = y/r, r = hypot(x, y)
-//
-// T&B Lecture 10 (Givens rotations).
-struct GivensRotation {
- double c; // cos(theta)
- double s; // sin(theta)
- std::size_t i; // first row/column index (second is i+1)
-
- // Construct the rotation that maps [x, y]^T → [hypot(x,y), 0]^T.
- // Returns the identity (c=1, s=0) when x == y == 0.
- [[nodiscard]] static GivensRotation make(double x, double y,
- std::size_t row_index);
-
- // Apply G from the left to rows i and i+1 of M, columns [col_start, n).
- // M[i:i+2, col_start:] ← G * M[i:i+2, col_start:]
- void apply_left(Matrix& M, std::size_t col_start = 0) const;
-
- // Apply G^T from the right to columns i and i+1 of M, rows [0, row_end).
- // M[0:row_end, i:i+2] ← M[0:row_end, i:i+2] * G^T
- void apply_right(Matrix& M, std::size_t row_end) const;
-};
-
-// Result of reducing A to upper Hessenberg form.
-// H is upper Hessenberg: H(i,j) = 0 for all i > j+1.
-// Q is orthogonal and A = Q H Q^T.
-// Ref: GVL Algorithm 7.4.2; T&B Lecture 26.
-struct HessenbergResult {
- Matrix H; // upper Hessenberg similarity of A
- Matrix Q; // accumulated orthogonal transformation
-};
-
-// Reduce A to upper Hessenberg form via Householder reflectors applied
-// from both sides. Costs O(10n³/3) flops; done once before QR iteration.
-// Ref: GVL §7.4.2 (Algorithm 7.4.2).
-//
-// Throws DimensionMismatchError if A is not square.
-[[nodiscard]] HessenbergResult hessenberg_reduction(const Matrix& A);
-
-// Apply one shifted QR step to an upper Hessenberg matrix H in-place,
-// using n-1 Givens rotations. Costs O(n²) vs O(n³) for Householder QR.
-// H remains upper Hessenberg after the step.
-// Ref: GVL §7.4.2; T&B Lecture 29.
-void hessenberg_qr_step(Matrix& H, double sigma);
-
-// Full practical QR algorithm:
-// 1. Reduce A to Hessenberg H = Q^T A Q (O(n³), done once).
-// 2. Run Wilkinson-shifted QR on H using Givens steps (O(n²) each).
-// Substantially faster than eigenvalues_shifted for n ≥ 50.
-// Ref: T&B Lecture 29.
-//
-// Same exceptions as eigenvalues_unshifted.
-[[nodiscard]] QRIterationResult eigenvalues_hessenberg(const Matrix& A,
- QRIterationOptions opts = {});
-
-} // namespace linalg
diff --git a/include/triangular_solve.hpp b/include/triangular_solve.hpp
deleted file mode 100644
index 632f225..0000000
--- a/include/triangular_solve.hpp
+++ /dev/null
@@ -1,22 +0,0 @@
-#pragma once
-
-#include <cstddef>
-
-#include "matrix.hpp"
-#include "vector.hpp"
-
-namespace linalg {
-
-Vector forward_substitution(
- const Matrix& lower,
- const Vector& rhs,
- double singular_tolerance = 1e-12,
- bool unit_diagonal = false);
-
-Vector backward_substitution(
- const Matrix& upper,
- const Vector& rhs,
- double singular_tolerance = 1e-12,
- bool unit_diagonal = false);
-
-} // namespace linalg
diff --git a/include/vector.hpp b/include/vector.hpp
deleted file mode 100644
index 3e797d8..0000000
--- a/include/vector.hpp
+++ /dev/null
@@ -1,47 +0,0 @@
-#pragma once
-
-#include <cstddef>
-#include <initializer_list>
-#include <vector>
-
-namespace linalg {
-
-class Vector {
-public:
- Vector() = default;
- explicit Vector(std::size_t n);
- Vector(std::size_t n, double value);
- Vector(std::initializer_list<double> values);
-
- [[nodiscard]] std::size_t size() const noexcept;
- [[nodiscard]] bool empty() const noexcept;
-
- double& operator[](std::size_t i);
- const double& operator[](std::size_t i) const;
-
- void fill(double value);
-
- double* data() noexcept;
- const double* data() const noexcept;
-
- auto begin() noexcept { return data_.begin(); }
- auto end() noexcept { return data_.end(); }
- auto begin() const noexcept { return data_.begin(); }
- auto end() const noexcept { return data_.end(); }
- auto cbegin() const noexcept { return data_.cbegin(); }
- auto cend() const noexcept { return data_.cend(); }
-
-private:
- void check_index(std::size_t i) const;
-
- std::vector<double> data_;
-};
-
-Vector operator+(const Vector& lhs, const Vector& rhs);
-Vector operator-(const Vector& lhs, const Vector& rhs);
-Vector operator*(const Vector& v, double scalar);
-Vector operator*(double scalar, const Vector& v);
-Vector operator/(const Vector& v, double scalar);
-double dot(const Vector& lhs, const Vector& rhs);
-
-} // namespace linalg
generated by cgit v1.2.3 (git 2.51.0) at 2026-09-01 23:49:27 +0000