1 files changed, 53 insertions, 28 deletions

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