Module refactor

9 files changed, 362 insertions, 379 deletions

diff --git a/src/linalgebra.cpp b/src/linalgebra.cpp
new file mode 100644
index 0000000..302810f
--- /dev/null
+++ b/src/linalgebra.cpp
@@ -0,0 +1,10 @@
+export module linalgebra;
+
+export import :error;
+export import :vector;
+export import :matrix;
+export import :norms;
+export import :triangular_solve;
+export import :lu;
+export import :qr;
+export import :qr_iteration;
diff --git a/src/linalgebra_error.cpp b/src/linalgebra_error.cpp
new file mode 100644
index 0000000..367c79d
--- /dev/null
+++ b/src/linalgebra_error.cpp
@@ -0,0 +1,29 @@
+export module linalgebra:error;
+import std;
+
+export namespace linalgebra {
+
+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 linalgebra
diff --git a/src/lu.cpp b/src/lu.cpp
index f6840ea..6d12eda 100644
--- a/src/lu.cpp
+++ b/src/lu.cpp
@@ -1,14 +1,26 @@
-#include "lu.hpp"
+export module linalgebra:lu;
+import std;
+import :error;
+import :vector;
+import :matrix;
+import :triangular_solve;
 
-#include <algorithm>
-#include <cmath>
-#include <numeric>
-#include <sstream>
+export namespace linalgebra {
 
-#include "linalg_error.hpp"
-#include "triangular_solve.hpp"
+struct LUResult {
+    Matrix L;
+    Matrix U;
+    std::vector<std::size_t> perm;
+    int sign;
+};
 
-namespace linalg {
+LUResult lu_factor(const Matrix& A, double singular_tolerance = 1e-12);
+
+Vector lu_solve(const LUResult& lu, const Vector& b);
+
+}  // namespace linalgebra
+
+namespace linalgebra {
 
 LUResult lu_factor(const Matrix& A, double singular_tolerance) {
     if (A.rows() != A.cols()) {
@@ -33,7 +45,6 @@ LUResult lu_factor(const Matrix& A, double singular_tolerance) {
     int sign = 1;
 
     for (std::size_t k = 0; k < n; ++k) {
-        // ---- Partial pivoting: find row with largest magnitude in column k ----
         std::size_t pivot_row = k;
         double max_val = std::abs(work(k, k));
         for (std::size_t i = k + 1; i < n; ++i) {
@@ -55,7 +66,6 @@ LUResult lu_factor(const Matrix& A, double singular_tolerance) {
             sign = -sign;
         }
 
-        // ---- Singularity check ----
         if (std::abs(work(k, k)) <= singular_tolerance) {
             std::ostringstream oss;
             oss << "lu_factor: near-zero pivot " << work(k, k) << " at step " << k
@@ -63,12 +73,10 @@ LUResult lu_factor(const Matrix& A, double singular_tolerance) {
             throw SingularMatrixError(oss.str());
         }
 
-        // ---- Record U row k ----
         for (std::size_t j = k; j < n; ++j) {
             U(k, j) = work(k, j);
         }
 
-        // ---- Compute multipliers and eliminate below pivot ----
         for (std::size_t i = k + 1; i < n; ++i) {
             L(i, k) = work(i, k) / work(k, k);
             for (std::size_t j = k + 1; j < n; ++j) {
@@ -101,4 +109,4 @@ Vector lu_solve(const LUResult& lu, const Vector& b) {
     return backward_substitution(lu.U, y);
 }
 
-}  // namespace linalg
+}  // namespace linalgebra
diff --git a/src/matrix.cpp b/src/matrix.cpp
index c27264b..75fcb07 100644
--- a/src/matrix.cpp
+++ b/src/matrix.cpp
@@ -1,120 +1,76 @@
-#include "matrix.hpp"
-#include "linalg_error.hpp"
+module;
 
-#include <algorithm>
-#include <cstddef>
-#include <sstream>
-#include <stdexcept>
-
-#if !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && \
-    (defined(__AVX__) || defined(__AVX2__) || defined(__AVX512F__))
-#include <immintrin.h>
-#endif
-
-#if !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__ARM_NEON) && defined(__aarch64__) && defined(__ARM_FEATURE_FP64_VECTOR_ARITHMETIC)
+#if !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__ARM_NEON)
 #include <arm_neon.h>
 #endif
 
-namespace linalg {
+export module linalgebra:matrix;
+import std;
+import :error;
+import :vector;
 
-namespace {
+export namespace linalgebra {
 
-void check_same_shape(const Matrix& lhs, const Matrix& rhs, const char* operation) {
-    if (lhs.rows() != rhs.rows() || lhs.cols() != rhs.cols()) {
-        std::ostringstream oss;
-        oss << operation << " requires equal matrix shapes, got " << lhs.rows() << "x" << lhs.cols() << " and " << rhs.rows() << "x" << rhs.cols();
-        throw DimensionMismatchError(oss.str());
-    }
-}
+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);
 
-double dot_product_scalar(const double* lhs, const double* rhs, std::size_t count) {
-    double sum = 0.0;
-    for (std::size_t i = 0; i < count; ++i) {
-        sum += lhs[i] * rhs[i];
-    }
-    return sum;
-}
-
-#if !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__AVX512F__)
-double horizontal_sum(__m512d values) {
-    alignas(64) double lanes[8];
-    _mm512_store_pd(lanes, values);
-    double sum = 0.0;
-    for (double lane : lanes) {
-        sum += lane;
-    }
-    return sum;
-}
+    [[nodiscard]] std::size_t rows() const noexcept;
+    [[nodiscard]] std::size_t cols() const noexcept;
+    [[nodiscard]] bool empty() const noexcept;
 
-double dot_product_avx512(const double* lhs, const double* rhs, std::size_t count) {
-    std::size_t i = 0;
-    __m512d acc0 = _mm512_setzero_pd();
-    __m512d acc1 = _mm512_setzero_pd();
+    double& operator()(std::size_t i, std::size_t j);
+    const double& operator()(std::size_t i, std::size_t j) const;
 
-    for (; i + 15 < count; i += 16) {
-        const __m512d lhs0 = _mm512_loadu_pd(lhs + i);
-        const __m512d rhs0 = _mm512_loadu_pd(rhs + i);
-        const __m512d lhs1 = _mm512_loadu_pd(lhs + i + 8);
-        const __m512d rhs1 = _mm512_loadu_pd(rhs + i + 8);
+    void fill(double value);
 
-        acc0 = _mm512_add_pd(acc0, _mm512_mul_pd(lhs0, rhs0));
-        acc1 = _mm512_add_pd(acc1, _mm512_mul_pd(lhs1, rhs1));
-    }
+    double* data() noexcept;
+    const double* data() const noexcept;
 
-    return horizontal_sum(acc0) + horizontal_sum(acc1) +
-           dot_product_scalar(lhs + i, rhs + i, count - i);
-}
-#endif
+    static Matrix identity(std::size_t n);
+    static Matrix zeros(std::size_t rows, std::size_t cols);
 
-#if !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__AVX2__)
-double horizontal_sum(__m256d values) {
-    alignas(32) double lanes[4];
-    _mm256_store_pd(lanes, values);
-    return lanes[0] + lanes[1] + lanes[2] + lanes[3];
-}
+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;
 
-double dot_product_avx2(const double* lhs, const double* rhs, std::size_t count) {
-    std::size_t i = 0;
-    __m256d acc0 = _mm256_setzero_pd();
-    __m256d acc1 = _mm256_setzero_pd();
+    std::size_t rows_ = 0;
+    std::size_t cols_ = 0;
+    std::vector<double> data_;
+};
 
-    for (; i + 7 < count; i += 8) {
-        const __m256d lhs0 = _mm256_loadu_pd(lhs + i);
-        const __m256d rhs0 = _mm256_loadu_pd(rhs + i);
-        const __m256d lhs1 = _mm256_loadu_pd(lhs + i + 4);
-        const __m256d rhs1 = _mm256_loadu_pd(rhs + i + 4);
+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);
 
-        acc0 = _mm256_add_pd(acc0, _mm256_mul_pd(lhs0, rhs0));
-        acc1 = _mm256_add_pd(acc1, _mm256_mul_pd(lhs1, rhs1));
-    }
+}  // namespace linalgebra
 
-    return horizontal_sum(acc0) + horizontal_sum(acc1) +
-           dot_product_scalar(lhs + i, rhs + i, count - i);
-}
-#endif
+namespace {
 
-#if !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__AVX__) && !defined(__AVX2__)
-double horizontal_sum(__m256d values) {
-    alignas(32) double lanes[4];
-    _mm256_store_pd(lanes, values);
-    return lanes[0] + lanes[1] + lanes[2] + lanes[3];
+void check_same_shape(const linalgebra::Matrix& lhs, const linalgebra::Matrix& rhs,
+                      const char* operation) {
+    if (lhs.rows() != rhs.rows() || lhs.cols() != rhs.cols()) {
+        std::ostringstream oss;
+        oss << operation << " requires equal matrix shapes, got " << lhs.rows() << "x"
+            << lhs.cols() << " and " << rhs.rows() << "x" << rhs.cols();
+        throw linalgebra::DimensionMismatchError(oss.str());
+    }
 }
 
-double dot_product_avx(const double* lhs, const double* rhs, std::size_t count) {
-    std::size_t i = 0;
-    __m256d acc = _mm256_setzero_pd();
-
-    for (; i + 3 < count; i += 4) {
-        const __m256d lhs_values = _mm256_loadu_pd(lhs + i);
-        const __m256d rhs_values = _mm256_loadu_pd(rhs + i);
-        acc = _mm256_add_pd(acc, _mm256_mul_pd(lhs_values, rhs_values));
+double dot_product_scalar(const double* lhs, const double* rhs, std::size_t count) {
+    double sum = 0.0;
+    for (std::size_t i = 0; i < count; ++i) {
+        sum += lhs[i] * rhs[i];
     }
-
-    return horizontal_sum(acc) + dot_product_scalar(lhs + i, rhs + i, count - i);
+    return sum;
 }
-#endif
 
-#if !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__ARM_NEON) && defined(__aarch64__) && defined(__ARM_FEATURE_FP64_VECTOR_ARITHMETIC)
+#if !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__ARM_NEON)
 double horizontal_sum(float64x2_t values) {
     return vgetq_lane_f64(values, 0) + vgetq_lane_f64(values, 1);
 }
@@ -140,13 +96,7 @@ double dot_product_neon(const double* lhs, const double* rhs, std::size_t count)
 #endif
 
 double dot_product_simd(const double* lhs, const double* rhs, std::size_t count) {
-#if !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__AVX512F__)
-    return dot_product_avx512(lhs, rhs, count);
-#elif !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__AVX2__)
-    return dot_product_avx2(lhs, rhs, count);
-#elif !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__AVX__) && !defined(__AVX2__)
-    return dot_product_avx(lhs, rhs, count);
-#elif !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__ARM_NEON) && defined(__aarch64__) && defined(__ARM_FEATURE_FP64_VECTOR_ARITHMETIC)
+#if !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__ARM_NEON)
     return dot_product_neon(lhs, rhs, count);
 #else
     return dot_product_scalar(lhs, rhs, count);
@@ -155,13 +105,16 @@ double dot_product_simd(const double* lhs, const double* rhs, std::size_t count)
 
 }  // namespace
 
+namespace linalgebra {
+
 Matrix::Matrix(std::size_t rows, std::size_t cols)
     : rows_(rows), cols_(cols), data_(rows * cols) {}
 
 Matrix::Matrix(std::size_t rows, std::size_t cols, double value)
     : rows_(rows), cols_(cols), data_(rows * cols, value) {}
 
-Matrix::Matrix(std::initializer_list<std::initializer_list<double>> values) : rows_(values.size()) {
+Matrix::Matrix(std::initializer_list<std::initializer_list<double>> values)
+    : rows_(values.size()) {
     if (rows_ == 0) {
         cols_ = 0;
         return;
@@ -178,7 +131,6 @@ Matrix::Matrix(std::initializer_list<std::initializer_list<double>> values) : ro
     }
 }
 
-
 std::size_t Matrix::rows() const noexcept { return rows_; }
 
 std::size_t Matrix::cols() const noexcept { return cols_; }
@@ -298,4 +250,4 @@ Matrix operator*(const Matrix& lhs, const Matrix& rhs) {
     return result;
 }
 
-}  // namespace linalg
+}  // namespace linalgebra
diff --git a/src/norms.cpp b/src/norms.cpp
index cafee69..7a40f1f 100644
--- a/src/norms.cpp
+++ b/src/norms.cpp
@@ -1,9 +1,15 @@
-#include "norms.hpp"
+export module linalgebra:norms;
+import std;
+import :vector;
 
-#include <cmath>
+export namespace linalgebra {
 
-namespace linalg {
+double norm2(const Vector& vector);
+
+}  // namespace linalgebra
+
+namespace linalgebra {
 
 double norm2(const Vector& vector) { return std::sqrt(dot(vector, vector)); }
 
-}  // namespace linalg
+}  // namespace linalgebra
diff --git a/src/qr.cpp b/src/qr.cpp
index 06770b2..7e32ce3 100644
--- a/src/qr.cpp
+++ b/src/qr.cpp
@@ -1,23 +1,55 @@
-#include "qr.hpp"
+export module linalgebra:qr;
+import std;
+import :error;
+import :vector;
+import :matrix;
 
-#include <cmath>
-#include <sstream>
+export namespace linalgebra {
 
-#include "linalg_error.hpp"
+struct QRResult {
+    Matrix Q;
+    Matrix R;
+};
 
-namespace linalg {
+// 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 linalgebra
 
 namespace {
 
-void require_tall(const Matrix& A, const char* name) {
+void require_tall(const linalgebra::Matrix& A, const char* name) {
     if (A.rows() < A.cols()) {
         std::ostringstream oss;
         oss << name << " requires rows >= cols, got " << A.rows() << "x" << A.cols();
-        throw DimensionMismatchError(oss.str());
+        throw linalgebra::DimensionMismatchError(oss.str());
     }
 }
 
-double col_norm(const Matrix& M, std::size_t j) {
+double col_norm(const linalgebra::Matrix& M, std::size_t j) {
     double s = 0.0;
     for (std::size_t i = 0; i < M.rows(); ++i) {
         s += M(i, j) * M(i, j);
@@ -25,7 +57,8 @@ double col_norm(const Matrix& M, std::size_t j) {
     return std::sqrt(s);
 }
 
-double col_dot(const Matrix& M, std::size_t j, const Matrix& N, std::size_t k) {
+double col_dot(const linalgebra::Matrix& M, std::size_t j,
+               const linalgebra::Matrix& N, std::size_t k) {
     double s = 0.0;
     for (std::size_t i = 0; i < M.rows(); ++i) {
         s += M(i, j) * N(i, k);
@@ -35,7 +68,7 @@ double col_dot(const Matrix& M, std::size_t j, const Matrix& N, std::size_t k) {
 
 }  // namespace
 
-// --- Gram-Schmidt ---
+namespace linalgebra {
 
 QRResult qr_classical_gs(const Matrix& A, double zero_tolerance) {
     require_tall(A, "qr_classical_gs");
@@ -49,7 +82,7 @@ QRResult qr_classical_gs(const Matrix& A, double zero_tolerance) {
         for (std::size_t i = 0; i < m; ++i) Q(i, j) = A(i, j);
 
         for (std::size_t k = 0; k < j; ++k) {
-            R(k, j) = col_dot(A, j, Q, k);       // <a_j, q_k>
+            R(k, j) = col_dot(A, j, Q, k);
             for (std::size_t i = 0; i < m; ++i) {
                 Q(i, j) -= R(k, j) * Q(i, k);
             }
@@ -69,8 +102,6 @@ QRResult qr_classical_gs(const Matrix& A, double zero_tolerance) {
     return QRResult{std::move(Q), std::move(R)};
 }
 
-// --- Modified Gram-Schmidt ---
-
 QRResult qr_modified_gs(const Matrix& A, double zero_tolerance) {
     require_tall(A, "qr_modified_gs");
     const std::size_t m = A.rows();
@@ -83,7 +114,7 @@ QRResult qr_modified_gs(const Matrix& A, double zero_tolerance) {
         for (std::size_t i = 0; i < m; ++i) Q(i, j) = A(i, j);
 
         for (std::size_t k = 0; k < j; ++k) {
-            R(k, j) = col_dot(Q, j, Q, k);       // <v_running, q_k>
+            R(k, j) = col_dot(Q, j, Q, k);
             for (std::size_t i = 0; i < m; ++i) {
                 Q(i, j) -= R(k, j) * Q(i, k);
             }
@@ -103,20 +134,16 @@ QRResult qr_modified_gs(const Matrix& A, double zero_tolerance) {
     return QRResult{std::move(Q), std::move(R)};
 }
 
-// --- Householder QR ---
-
 QRResult qr_householder(const Matrix& A) {
     require_tall(A, "qr_householder");
     const std::size_t m = A.rows();
     const std::size_t n = A.cols();
 
-    // Will become R.
     Matrix work = A;
-
     Matrix Q_full = Matrix::identity(m);
 
     for (std::size_t k = 0; k < n; ++k) {
-        const std::size_t p = m - k;  // length of the subvector
+        const std::size_t p = m - k;
 
         std::vector<double> u(p);
         for (std::size_t i = 0; i < p; ++i) u[i] = work(k + i, k);
@@ -139,19 +166,17 @@ QRResult qr_householder(const Matrix& A) {
         }();
         const double tau = 2.0 / utu;
 
-        // Apply H_k to work[k:, k:n] 
         for (std::size_t j = k; j < n; ++j) {
-            double dot = 0.0;
-            for (std::size_t i = 0; i < p; ++i) dot += u[i] * work(k + i, j);
-            const double coeff = tau * dot;
+            double d = 0.0;
+            for (std::size_t i = 0; i < p; ++i) d += u[i] * work(k + i, j);
+            const double coeff = tau * d;
             for (std::size_t i = 0; i < p; ++i) work(k + i, j) -= coeff * u[i];
         }
 
-        // Apply H_k to Q_full[k:, 0:m] 
         for (std::size_t j = 0; j < m; ++j) {
-            double dot = 0.0;
-            for (std::size_t i = 0; i < p; ++i) dot += u[i] * Q_full(k + i, j);
-            const double coeff = tau * dot;
+            double d = 0.0;
+            for (std::size_t i = 0; i < p; ++i) d += u[i] * Q_full(k + i, j);
+            const double coeff = tau * d;
             for (std::size_t i = 0; i < p; ++i) Q_full(k + i, j) -= coeff * u[i];
         }
     }
@@ -169,4 +194,4 @@ QRResult qr_householder(const Matrix& A) {
     return QRResult{std::move(Q), std::move(R)};
 }
 
-}  // namespace linalg
+}  // namespace linalgebra
diff --git a/src/qr_iteration.cpp b/src/qr_iteration.cpp
index e476781..64fd25f 100644
--- a/src/qr_iteration.cpp
+++ b/src/qr_iteration.cpp
@@ -1,73 +1,99 @@
-#include "qr_iteration.hpp"
-
+module;
 #include <cassert>
-#include <cmath>
-#include <sstream>
 
-#include "linalg_error.hpp"
-#include "matrix.hpp"
-#include "qr.hpp"
-#include "vector.hpp"
+export module linalgebra:qr_iteration;
+import std;
+import :error;
+import :vector;
+import :matrix;
+import :qr;
 
 // References used throughout this file:
 //   T&B  — Trefethen & Bau, "Numerical Linear Algebra"
 //   GVL  — Golub & Van Loan, "Matrix Computations" 4th ed.
 
-namespace linalg {
+export namespace linalgebra {
 
-namespace {
+// ---------------------------------------------------------------------------
+// Options
+// ---------------------------------------------------------------------------
+
+struct QRIterationOptions {
+    double tolerance = 1e-10;
+    int max_iterations = 1000;
+    bool track_convergence = false;
+};
+
+struct QRIterationResult {
+    Vector eigenvalues_real;
+    Vector eigenvalues_imag;
+    int iterations = 0;
+    // std::vector is used here because linalgebra::Vector has no push_back;
+    // convergence_history is a plain time-series container, not a math object.
+    std::vector<double> convergence_history;
+};
+
+[[nodiscard]] QRIterationResult eigenvalues_unshifted(const Matrix& A,
+                                                      QRIterationOptions opts = {});
+
+[[nodiscard]] QRIterationResult eigenvalues_shifted(const Matrix& A,
+                                                    QRIterationOptions opts = {});
+
+// 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)
+struct GivensRotation {
+    double      c;
+    double      s;
+    std::size_t i;
+
+    [[nodiscard]] static GivensRotation make(double x, double y, std::size_t row_index);
+    void apply_left(Matrix& M, std::size_t col_start = 0) const;
+    void apply_right(Matrix& M, std::size_t row_end) const;
+};
 
+struct HessenbergResult {
+    Matrix H;
+    Matrix Q;
+};
 
-// Frobenius norm of the strict lower triangle of an n×n matrix.
-// This is the standard convergence diagnostic for QR iteration: as A_k
-// approaches the real Schur form, all entries below the main diagonal
-// (excluding 2×2 block sub-diagonals) tend to zero.
-// Ref: T&B §28; used as the convergence criterion in Algorithm 28.1.
-double lower_triangle_norm(const Matrix& A) {
+[[nodiscard]] HessenbergResult hessenberg_reduction(const Matrix& A);
+
+void hessenberg_qr_step(Matrix& H, double sigma);
+
+[[nodiscard]] QRIterationResult eigenvalues_hessenberg(const Matrix& A,
+                                                       QRIterationOptions opts = {});
+
+}  // namespace linalgebra
+
+namespace {
+
+double lower_triangle_norm(const linalgebra::Matrix& A) {
     const std::size_t n = A.rows();
     double s = 0.0;
-    for (std::size_t i = 1; i < n; ++i)      // row 1 .. n-1
-        for (std::size_t j = 0; j < i; ++j)  // col 0 .. i-1  (strict lower)
+    for (std::size_t i = 1; i < n; ++i)
+        for (std::size_t j = 0; j < i; ++j)
             s += A(i, j) * A(i, j);
     return std::sqrt(s);
 }
 
-// Extract eigenvalues from a quasi-upper-triangular matrix (real Schur form).
-//
-// Scans the diagonal from top-left to bottom-right.  At each position i:
-//   — |A(i+1, i)| < tol  → 1×1 block: real eigenvalue A(i,i), imag = 0.
-//   — otherwise           → 2×2 block [A(i..i+1, i..i+1)]: eigenvalues via
-//                           quadratic formula.  When the discriminant is
-//                           negative the result is a complex-conjugate pair,
-//                           stored as (re, +im) and (re, -im) in the real
-//                           and imaginary part Vectors.
-//
-// Fills positions 0..n-1 of `real_out` and `imag_out` (pre-sized to n).
-//
-// Ref: T&B Lecture 28; GVL §7.4.1.
-void extract_eigenvalues(const Matrix& T, double tol,
-                         Vector& real_out, Vector& imag_out) {
+void extract_eigenvalues(const linalgebra::Matrix& T, double tol,
+                         linalgebra::Vector& real_out, linalgebra::Vector& imag_out) {
     const std::size_t n = T.rows();
-    std::size_t out = 0;  
-    std::size_t i   = 0; 
+    std::size_t out = 0;
+    std::size_t i   = 0;
 
     while (i < n) {
         const bool is_last   = (i + 1 == n);
         const bool sub_small = is_last || (std::abs(T(i + 1, i)) < tol);
 
         if (sub_small) {
-            // 1×1 block: real eigenvalue.
             real_out[out] = T(i, i);
             imag_out[out] = 0.0;
             ++out;
             ++i;
         } else {
-            // 2×2 block:
-            //   | a  b |
-            //   | c  d |
-            // Characteristic polynomial: lambda^2 - (a+d)*lambda + (ad - bc) = 0.
-            // Discriminant: (a-d)^2 + 4*b*c.
-            // Ref: GVL §7.4.1.
             const double a = T(i, i);
             const double b = T(i, i + 1);
             const double c = T(i + 1, i);
@@ -76,15 +102,12 @@ void extract_eigenvalues(const Matrix& T, double tol,
             const double disc = (a - d) * (a - d) + 4.0 * b * c;
 
             if (disc >= 0.0) {
-                // Real eigenvalues unusual in converged real Schur form, but
-                // handled robustly in case the block didn't fully split.
                 const double sq   = std::sqrt(disc);
                 real_out[out]     = 0.5 * (tr + sq);
                 imag_out[out]     = 0.0;
                 real_out[out + 1] = 0.5 * (tr - sq);
                 imag_out[out + 1] = 0.0;
             } else {
-                // Complex-conjugate pair: real part ± imaginary part.
                 const double re       = 0.5 * tr;
                 const double im       = 0.5 * std::sqrt(-disc);
                 real_out[out]         = re;
@@ -100,37 +123,32 @@ void extract_eigenvalues(const Matrix& T, double tol,
     assert(out == n);
 }
 
-void require_square(const Matrix& A, const char* fname) {
+void require_square(const linalgebra::Matrix& A, const char* fname) {
     if (A.rows() != A.cols()) {
         std::ostringstream oss;
         oss << fname << ": requires a square matrix, got "
             << A.rows() << "x" << A.cols();
-        throw DimensionMismatchError(oss.str());
+        throw linalgebra::DimensionMismatchError(oss.str());
     }
 }
 
+double wilkinson_shift(const linalgebra::Matrix& A) {
+    const std::size_t n = A.rows();
+    const double a = A(n - 2, n - 2);
+    const double b = A(n - 1, n - 2);
+    const double d = A(n - 1, n - 1);
+    const double delta = 0.5 * (a - d);
+    const double denom = std::abs(delta) + std::hypot(delta, b);
+    if (denom == 0.0) return d;
+    const double sgn = (delta >= 0.0) ? 1.0 : -1.0;
+    return d - sgn * (b * b) / denom;
+}
+
 }  // namespace
 
-// --- Unshifted QR iteration ---
-//
-// Each step performs an orthogonal similarity transformation:
-//   A_{k-1} = Q_k R_k   (Householder QR; backward-stable)
-//   A_k     = R_k Q_k   = Q_k^T A_{k-1} Q_k
-//
-// Similarity preserves eigenvalues (GVL §7.3.1, Theorem 7.3.1).
-// The iterates converge to the real Schur form: a quasi-upper-triangular
-// matrix whose 1×1 blocks give real eigenvalues and 2×2 blocks give
-// complex-conjugate pairs.
-//
-// Convergence rate: linear, with per-step reduction factor
-//   |lambda_{j+1} / lambda_j|  for the (j, j+1) coupling.
-// (T&B Lecture 28, Theorem 28.2; GVL §7.3.2)
-//
-// Each iteration costs O(n^3) due to full Householder QR; Hessenberg
-// reduction (Stage 3) reduces subsequent steps to O(n^2).
+namespace linalgebra {
 
-QRIterationResult eigenvalues_unshifted(const Matrix& A,
-                                        QRIterationOptions opts) {
+QRIterationResult eigenvalues_unshifted(const Matrix& A, QRIterationOptions opts) {
     require_square(A, "eigenvalues_unshifted");
     const std::size_t n = A.rows();
 
@@ -153,13 +171,9 @@ QRIterationResult eigenvalues_unshifted(const Matrix& A,
     Matrix Ak = A;
 
     for (int k = 0; k < opts.max_iterations; ++k) {
-        // Factor A_{k-1} = Q R using backward-stable Householder reflections.
         const QRResult qr = qr_householder(Ak);
-
-        // A_k = R Q  (orthogonal similarity: Q^T A_{k-1} Q)
         Ak = qr.R * qr.Q;
 
-        // --- Convergence check ---
         const double lower_norm = lower_triangle_norm(Ak);
 
         if (opts.track_convergence) {
@@ -184,52 +198,6 @@ QRIterationResult eigenvalues_unshifted(const Matrix& A,
     throw NonConvergenceError(oss.str());
 }
 
-// --- Wilkinson-shifted QR iteration ---
-//
-// The Wilkinson shift is the eigenvalue of the bottom-right 2×2 block
-//   | a  b |
-//   | c  d |
-// that is closest to d (the trailing diagonal entry).
-//
-// Exact eigenvalue formula: μ_{1,2} = (a+d)/2 ± sqrt(((a-d)/2)² + b·c)
-// We pick the one with |μ - d| smaller.
-//
-// When the discriminant is negative (complex eigenvalues), fall back to σ = d
-// (Rayleigh quotient shift), which still accelerates convergence.
-//
-// Ref: T&B Lecture 29; GVL §7.4.2.
-
-namespace {
-
-double wilkinson_shift(const Matrix& A) {
-    const std::size_t n = A.rows();
-    const double a = A(n - 2, n - 2);
-    const double b = A(n - 1, n - 2);  // subdiagonal entry only
-    const double d = A(n - 1, n - 1);
-    const double delta = 0.5 * (a - d);
-    const double denom = std::abs(delta) + std::hypot(delta, b);
-    if (denom == 0.0) return d;
-    const double sgn = (delta >= 0.0) ? 1.0 : -1.0;
-    return d - sgn * (b * b) / denom;
-}
-
-}  // namespace
-
-// eigenvalues_shifted — Wilkinson-shifted QR with trailing deflation.
-//
-// After each QR step we check whether the trailing subdiagonal entry of the
-// active block is negligible (relative criterion: GVL §7.4.1).  If so, the
-// bottom diagonal entry is accepted as a converged eigenvalue and the active
-// subproblem shrinks by one.  This "trailing deflation" enables the cubic
-// convergence promised by the Wilkinson shift to compound across successive
-// eigenvalues rather than stalling on the full lower-triangle norm.
-//
-// When the active size reaches 2 we extract both eigenvalues analytically
-// from the 2×2 block (handling real and complex-conjugate pairs) rather than
-// continuing to iterate.  For symmetric inputs this is always a real pair.
-//
-// Ref: GVL §7.5.1; T&B Lecture 29.
-
 QRIterationResult eigenvalues_shifted(const Matrix& A, QRIterationOptions opts) {
     require_square(A, "eigenvalues_shifted");
     const std::size_t n = A.rows();
@@ -250,7 +218,7 @@ QRIterationResult eigenvalues_shifted(const Matrix& A, QRIterationOptions opts)
     Matrix Ak = A;
 
     std::size_t n_found = n;
-    std::size_t active  = n;  // live subproblem is rows/cols 0..active-1
+    std::size_t active  = n;
 
     auto store_real = [&](double re) {
         --n_found;
@@ -281,16 +249,14 @@ QRIterationResult eigenvalues_shifted(const Matrix& A, QRIterationOptions opts)
     };
 
     for (int k = 0; k < opts.max_iterations; ++k) {
-        // --- Deflation sweep ---
         while (active >= 2) {
             const double sub   = std::abs(Ak(active - 1, active - 2));
             const double scale = std::abs(Ak(active - 2, active - 2))
                                + std::abs(Ak(active - 1, active - 1));
-            // Relative + absolute floor tolerance (GVL §7.4.1).
             const double deflation_tol =
                 opts.tolerance * (scale > 0.0 ? scale : 1.0);
             if (sub > deflation_tol) break;
-            Ak(active - 1, active - 2) = 0.0;  // enforce exact zero
+            Ak(active - 1, active - 2) = 0.0;
             store_real(Ak(active - 1, active - 1));
             --active;
         }
@@ -299,7 +265,6 @@ QRIterationResult eigenvalues_shifted(const Matrix& A, QRIterationOptions opts)
         if (active == 1) { store_real(Ak(0, 0)); active = 0; break; }
         if (active == 2) { close_2x2();                       break; }
 
-        // --- Wilkinson-shifted QR step on the active × active subblock ---
         Matrix sub_mat(active, active);
         for (std::size_t i = 0; i < active; ++i)
             for (std::size_t j = 0; j < active; ++j)
@@ -331,8 +296,6 @@ QRIterationResult eigenvalues_shifted(const Matrix& A, QRIterationOptions opts)
     return result;
 }
 
-// --- Givens rotation ---
-
 GivensRotation GivensRotation::make(double x, double y, std::size_t row_index) {
     const double r = std::hypot(x, y);
     if (r == 0.0) return {1.0, 0.0, row_index};
@@ -340,9 +303,6 @@ GivensRotation GivensRotation::make(double x, double y, std::size_t row_index) {
 }
 
 void GivensRotation::apply_left(Matrix& M, std::size_t col_start) const {
-    // Rows i and i+1, columns col_start..n-1.
-    // [ c  s] [x]   [cx + sy]
-    // [-s  c] [y] = [-sx + cy]
     for (std::size_t j = col_start; j < M.cols(); ++j) {
         const double xi  = M(i,     j);
         const double xi1 = M(i + 1, j);
@@ -352,10 +312,6 @@ void GivensRotation::apply_left(Matrix& M, std::size_t col_start) const {
 }
 
 void GivensRotation::apply_right(Matrix& M, std::size_t row_end) const {
-    // Columns i and i+1, rows 0..row_end-1.
-    // M * G^T where G^T = [c -s; s c]:
-    // new col i   = c * old_i + s * old_{i+1}
-    // new col i+1 = -s * old_i + c * old_{i+1}
     for (std::size_t j = 0; j < row_end; ++j) {
         const double xi  = M(j, i);
         const double xi1 = M(j, i + 1);
@@ -364,16 +320,6 @@ void GivensRotation::apply_right(Matrix& M, std::size_t row_end) const {
     }
 }
 
-// --- Hessenberg reduction ---
-// For k = 0, 1, ..., n-3:
-//   Build a Householder reflector H_k that zeros A[k+2:n, k].
-//   Apply from left:  A[k+1:n, k:n] ← H_k * A[k+1:n, k:n]
-//   Apply from right: A[0:n, k+1:n] ← A[0:n, k+1:n] * H_k
-//   Accumulate Q:     Q[0:n, k+1:n] ← Q[0:n, k+1:n] * H_k
-//
-// H_k is never formed explicitly; applied via rank-1 update with tau = 2/uᵀu.
-// Ref: GVL §7.4.2 (Algorithm 7.4.2).
-
 HessenbergResult hessenberg_reduction(const Matrix& A) {
     require_square(A, "hessenberg_reduction");
     const std::size_t n = A.rows();
@@ -382,10 +328,9 @@ HessenbergResult hessenberg_reduction(const Matrix& A) {
     Matrix Q = Matrix::identity(n);
 
     for (std::size_t k = 0; k + 2 <= n; ++k) {
-        const std::size_t p = n - k - 1;  // p = n - (k+1)
+        const std::size_t p = n - k - 1;
         if (p == 0) break;
 
-        // Build Householder vector u from H[k+1:n, k].
         std::vector<double> u(p);
         for (std::size_t i = 0; i < p; ++i) u[i] = H(k + 1 + i, k);
 
@@ -402,27 +347,24 @@ HessenbergResult hessenberg_reduction(const Matrix& A) {
         for (double v : u) utu += v * v;
         const double tau = 2.0 / utu;
 
-        // Apply H_k from the LEFT to H[k+1:n, k:n].
         for (std::size_t j = k; j < n; ++j) {
-            double dot = 0.0;
-            for (std::size_t i = 0; i < p; ++i) dot += u[i] * H(k + 1 + i, j);
-            const double coeff = tau * dot;
+            double d = 0.0;
+            for (std::size_t i = 0; i < p; ++i) d += u[i] * H(k + 1 + i, j);
+            const double coeff = tau * d;
             for (std::size_t i = 0; i < p; ++i) H(k + 1 + i, j) -= coeff * u[i];
         }
 
-        // Apply H_k from the RIGHT to H[0:n, k+1:n].
         for (std::size_t j = 0; j < n; ++j) {
-            double dot = 0.0;
-            for (std::size_t i = 0; i < p; ++i) dot += H(j, k + 1 + i) * u[i];
-            const double coeff = tau * dot;
+            double d = 0.0;
+            for (std::size_t i = 0; i < p; ++i) d += H(j, k + 1 + i) * u[i];
+            const double coeff = tau * d;
             for (std::size_t i = 0; i < p; ++i) H(j, k + 1 + i) -= coeff * u[i];
         }
 
-        // Accumulate Q: Q[0:n, k+1:n] ← Q[0:n, k+1:n] * H_k.
         for (std::size_t j = 0; j < n; ++j) {
-            double dot = 0.0;
-            for (std::size_t i = 0; i < p; ++i) dot += Q(j, k + 1 + i) * u[i];
-            const double coeff = tau * dot;
+            double d = 0.0;
+            for (std::size_t i = 0; i < p; ++i) d += Q(j, k + 1 + i) * u[i];
+            const double coeff = tau * d;
             for (std::size_t i = 0; i < p; ++i) Q(j, k + 1 + i) -= coeff * u[i];
         }
 
@@ -432,20 +374,6 @@ HessenbergResult hessenberg_reduction(const Matrix& A) {
     return HessenbergResult{std::move(H), std::move(Q)};
 }
 
-// --- Hessenberg QR step via Givens rotations ---
-//
-// One shifted QR step on the upper Hessenberg matrix H:
-//   1. Shift: H ← H - σI.
-//   2. For k = 0..n-2: compute G_k = Givens(H(k,k), H(k+1,k));
-//                        apply G_k from left to rows k,k+1 of H,
-//                        starting from column k (Hessenberg: H(k+1,j)=0, j<k).
-//   3. For k = 0..n-2: apply G_k^T from right to cols k,k+1 of H,
-//                        up to row k+2 (exploits upper-triangular structure).
-//   4. Unshift: H ← H + σI.
-//
-// After the step H is again upper Hessenberg (GVL §7.4.2, Theorem 7.4.1).
-// Total cost: O(n²).  Ref: GVL §7.4.2.
-
 void hessenberg_qr_step(Matrix& H, double sigma) {
     const std::size_t n = H.rows();
 
@@ -455,10 +383,7 @@ void hessenberg_qr_step(Matrix& H, double sigma) {
     gs.reserve(n - 1);
 
     for (std::size_t k = 0; k + 1 < n; ++k) {
-        // Eliminate H(k+1, k) via a rotation on rows k and k+1.
         GivensRotation g = GivensRotation::make(H(k, k), H(k + 1, k), k);
-        // Left application: rows k, k+1; columns k..n-1.
-        // (Hessenberg: H(k+1, j) = 0 for j < k, so starting from col k is exact.)
         g.apply_left(H, k);
         gs.push_back(g);
     }
@@ -467,23 +392,10 @@ void hessenberg_qr_step(Matrix& H, double sigma) {
         gs[k].apply_right(H, std::min(k + 2, n));
     }
 
-    // Unshift.
     for (std::size_t j = 0; j < n; ++j) H(j, j) += sigma;
 }
 
-// --- Full QR algorithm ---
-//
-// Same outer deflation loop as eigenvalues_shifted, but each QR step uses
-// hessenberg_qr_step (O(n²) Givens rotations) instead of full Householder QR
-// (O(n³)).  After Hessenberg reduction the matrix stays Hessenberg throughout,
-// so the O(n²) per-step cost applies for every step after the one-time O(n³)
-// reduction.  Total cost is thus O(n³) + O(iterations · n²), which beats
-// eigenvalues_shifted's O(iterations · n³) for large n.
-//
-// Ref: GVL §7.4.2; T&B Lecture 29.
-
-QRIterationResult eigenvalues_hessenberg(const Matrix& A,
-                                         QRIterationOptions opts) {
+QRIterationResult eigenvalues_hessenberg(const Matrix& A, QRIterationOptions opts) {
     require_square(A, "eigenvalues_hessenberg");
     const std::size_t n = A.rows();
 
@@ -535,7 +447,6 @@ QRIterationResult eigenvalues_hessenberg(const Matrix& A,
     };
 
     for (int k = 0; k < opts.max_iterations; ++k) {
-        // --- Deflation sweep ---
         while (active >= 2) {
             const double sub   = std::abs(H(active - 1, active - 2));
             const double scale = std::abs(H(active - 2, active - 2))
@@ -552,7 +463,6 @@ QRIterationResult eigenvalues_hessenberg(const Matrix& A,
         if (active == 1) { store_real(H(0, 0)); active = 0; break; }
         if (active == 2) { close_2x2();                       break; }
 
-        // Wilkinson shift from trailing 2×2 of the active block.
         const double a_w = H(active - 2, active - 2);
         const double b_w = H(active - 1, active - 2);
         const double d_w = H(active - 1, active - 1);
@@ -561,7 +471,6 @@ QRIterationResult eigenvalues_hessenberg(const Matrix& A,
         const double sigma = (denom == 0.0) ? d_w
             : d_w - ((delta >= 0.0) ? 1.0 : -1.0) * (b_w * b_w) / denom;
 
-        // O(n²) Givens step on the active×active Hessenberg subblock.
         Matrix sub_H(active, active);
         for (std::size_t ii = 0; ii < active; ++ii)
             for (std::size_t jj = 0; jj < active; ++jj)
@@ -593,4 +502,4 @@ QRIterationResult eigenvalues_hessenberg(const Matrix& A,
     return result;
 }
 
-}  // namespace linalg
+}  // namespace linalgebra
diff --git a/src/triangular_solve.cpp b/src/triangular_solve.cpp
index 6a5e8ad..58dda6b 100644
--- a/src/triangular_solve.cpp
+++ b/src/triangular_solve.cpp
@@ -1,28 +1,41 @@
-#include "triangular_solve.hpp"
+export module linalgebra:triangular_solve;
+import std;
+import :error;
+import :vector;
+import :matrix;
 
-#include "linalg_error.hpp"
+export namespace linalgebra {
 
-#include <cmath>
-#include <sstream>
-#include <stdexcept>
+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 {
+}  // namespace linalgebra
 
 namespace {
 
-void validate_square_system(const Matrix& matrix, const Vector& rhs, const char* operation) {
+void validate_square_system(const linalgebra::Matrix& matrix, const linalgebra::Vector& rhs,
+                             const char* operation) {
     if (matrix.rows() != matrix.cols()) {
         std::ostringstream oss;
         oss << operation << " requires a square matrix, got " << matrix.rows() << "x"
             << matrix.cols();
-        throw DimensionMismatchError(oss.str());
+        throw linalgebra::DimensionMismatchError(oss.str());
     }
 
     if (matrix.rows() != rhs.size()) {
         std::ostringstream oss;
         oss << operation << " requires matrix dimension to match rhs size, got "
             << matrix.rows() << " and " << rhs.size();
-        throw DimensionMismatchError(oss.str());
+        throw linalgebra::DimensionMismatchError(oss.str());
     }
 }
 
@@ -32,10 +45,8 @@ void validate_tolerance(double singular_tolerance) {
     }
 }
 
-void validate_lower_triangular(
-    const Matrix& lower,
-    double singular_tolerance,
-    bool unit_diagonal) {
+void validate_lower_triangular(const linalgebra::Matrix& lower, double singular_tolerance,
+                                bool unit_diagonal) {
     for (std::size_t i = 0; i < lower.rows(); ++i) {
         for (std::size_t j = i + 1; j < lower.cols(); ++j) {
             if (std::abs(lower(i, j)) > singular_tolerance) {
@@ -45,16 +56,14 @@ void validate_lower_triangular(
         }
 
         if (!unit_diagonal && std::abs(lower(i, i)) <= singular_tolerance) {
-            throw SingularMatrixError(
+            throw linalgebra::SingularMatrixError(
                 "Forward substitution encountered a zero or tiny diagonal entry");
         }
     }
 }
 
-void validate_upper_triangular(
-    const Matrix& upper,
-    double singular_tolerance,
-    bool unit_diagonal) {
+void validate_upper_triangular(const linalgebra::Matrix& upper, double singular_tolerance,
+                                bool unit_diagonal) {
     for (std::size_t i = 0; i < upper.rows(); ++i) {
         for (std::size_t j = 0; j < i; ++j) {
             if (std::abs(upper(i, j)) > singular_tolerance) {
@@ -64,7 +73,7 @@ void validate_upper_triangular(
         }
 
         if (!unit_diagonal && std::abs(upper(i, i)) <= singular_tolerance) {
-            throw SingularMatrixError(
+            throw linalgebra::SingularMatrixError(
                 "Backward substitution encountered a negligible diagonal entry");
         }
     }
@@ -72,11 +81,10 @@ void validate_upper_triangular(
 
 }  // namespace
 
-Vector forward_substitution(
-    const Matrix& lower,
-    const Vector& rhs,
-    double singular_tolerance,
-    bool unit_diagonal) {
+namespace linalgebra {
+
+Vector forward_substitution(const Matrix& lower, const Vector& rhs,
+                             double singular_tolerance, bool unit_diagonal) {
     validate_tolerance(singular_tolerance);
     validate_square_system(lower, rhs, "Forward substitution");
     validate_lower_triangular(lower, singular_tolerance, unit_diagonal);
@@ -98,11 +106,8 @@ Vector forward_substitution(
     return solution;
 }
 
-Vector backward_substitution(
-    const Matrix& upper,
-    const Vector& rhs,
-    double singular_tolerance,
-    bool unit_diagonal) {
+Vector backward_substitution(const Matrix& upper, const Vector& rhs,
+                              double singular_tolerance, bool unit_diagonal) {
     validate_tolerance(singular_tolerance);
     validate_square_system(upper, rhs, "Backward substitution");
     validate_upper_triangular(upper, singular_tolerance, unit_diagonal);
@@ -125,4 +130,4 @@ Vector backward_substitution(
     return solution;
 }
 
-}  // namespace linalg
+}  // namespace linalgebra
diff --git a/src/vector.cpp b/src/vector.cpp
index 8c4a8c4..ed86f88 100644
--- a/src/vector.cpp
+++ b/src/vector.cpp
@@ -1,26 +1,65 @@
-#include "vector.hpp"
-#include "linalg_error.hpp"
+export module linalgebra:vector;
+import std;
+import :error;
 
-#include <algorithm>
-#include <numeric>
-#include <sstream>
-#include <stdexcept>
+export namespace linalgebra {
 
-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 linalgebra
 
 namespace {
 
-void check_same_size(const Vector& lhs, const Vector& rhs, const char* operation) {
+void check_same_size(const linalgebra::Vector& lhs, const linalgebra::Vector& rhs,
+                     const char* operation) {
     if (lhs.size() != rhs.size()) {
         std::ostringstream oss;
         oss << operation << " requires equal vector sizes, got " << lhs.size() << " and "
             << rhs.size();
-        throw DimensionMismatchError(oss.str());
+        throw linalgebra::DimensionMismatchError(oss.str());
     }
 }
 
 }  // namespace
 
+namespace linalgebra {
+
 Vector::Vector(std::size_t n) : data_(n) {}
 
 Vector::Vector(std::size_t n, double value) : data_(n, value) {}
@@ -101,4 +140,4 @@ double dot(const Vector& lhs, const Vector& rhs) {
     return std::inner_product(lhs.begin(), lhs.end(), rhs.begin(), 0.0);
 }
 
-}  // namespace linalg
+}  // namespace linalgebra