feat: Heisenberg matrix reduction and workspace refactor

.github/workflows/clippy-lint.yml

@@ -13,15 +13,7 @@ env:
 jobs:
   clippy_check:
     runs-on: ubuntu-latest
-    strategy:
-      fail-fast: false
-      matrix:
-        crate-target:
-          - { crate: spindalis }
-          - { crate: spindalis_core }
-          - { crate: spindalis_macros }
     steps:
       - uses: actions/checkout@v5
       - name: Run Clippy
-        working-directory: ./${{ matrix.crate-target.crate }}
-        run: cargo clippy --all-targets --all-features
+        run: cargo clippy --workspace --all-targets --all-features

.github/workflows/rust_tests.yml

@@ -13,19 +13,9 @@ jobs:
   build:
 
     runs-on: ubuntu-latest
-    strategy:
-      fail-fast: false
-      matrix:
-        crate-target:
-          - { crate: spindalis }
-          - { crate: spindalis_core }
-          - { crate: spindalis_macros }
-
     steps:
     - uses: actions/checkout@v4
     - name: Build
-      working-directory: ./${{ matrix.crate-target.crate }}
-      run: cargo build --verbose
+      run: cargo build --verbose --workspace
     - name: Run tests
-      working-directory: ./${{ matrix.crate-target.crate }}
-      run: cargo test --verbose
+      run: cargo test --verbose --workspace

.justfile

@@ -4,52 +4,30 @@ set windows-shell := ["cmd.exe", "/c"]
 ROOT_DIR := justfile_directory()
 submodules := "spindalis spindalis_core spindalis_macros"
 
-test_cmd := if os_family() == "windows" { \
-    f"cmd.exe /c FOR %f IN ({{submodules}}) DO " + \
-    f"cargo test --manifest-path {{ROOT_DIR}}/%f/Cargo.toml" \
-} else { \
-    f"for mod in {{submodules}}; do " + \
-    f"cargo test --manifest-path {{ROOT_DIR}}/$mod/Cargo.toml; done" \
-}
-
-clippy_cmd := if os_family() == "windows" { \
-    f"cmd.exe /c FOR %f IN ({{submodules}}) DO " + \
-    "cargo clippy --all-targets --all-features " + \
-    f"--manifest-path {{ROOT_DIR}}/%f/Cargo.toml" \
-} else { \
-    f"for mod in {{submodules}}; do " + \
-    "cargo clippy --all-targets --all-features " + \
-    f"--manifest-path {{ROOT_DIR}}/$mod/Cargo.toml; done" \
-}
-
-fmt_cmd := if os_family() == "windows" { \
-    f"cmd.exe /c FOR %f IN ({{submodules}}) DO " + \
-    f"cargo fmt --manifest-path {{ROOT_DIR}}/%f/Cargo.toml" \
-} else { \
-    f"for mod in {{submodules}}; do " + \
-    f"cargo fmt --manifest-path {{ROOT_DIR}}/$mod/Cargo.toml; done" \
-}
+test_cmd := "cargo test --workspace"
+clippy_cmd := "cargo clippy --workspace --all-targets --all-features"
+fmt_cmd := "cargo fmt --all"
 
 
 recipes:
     just -l
 
 test TEST:
-    cargo test --mainfest-path {{ROOT_DIR}}/{{TEST}}/Cargo.toml
+    cargo test -p {{TEST}}
 
 [doc('Test all modules')]
 test-all:
     {{test_cmd}}
 
 lint LINT:
-    cargo clippy --all-targets --all-features --mainfest-path {{ROOT_DIR}}/{{LINT}}/Cargo.toml
+    cargo clippy --all-targets --all-features -p {{LINT}}
 
 [doc('Lint all modules with clippy')]
 lint-all:
     {{clippy_cmd}}
 
 format FORMAT:
-    cargo fmt --mainfest-path {{ROOT_DIR}}/{{FORMAT}}/Cargo.toml
+    cargo fmt -p {{FORMAT}}
 
 [doc('Format all modules')]
 format-all:

Cargo.toml

@@ -0,0 +1,7 @@
+[workspace]
+members = [
+    "spindalis",
+    "spindalis_core",
+    "spindalis_macros",
+]
+resolver = "2"

spindalis/src/reduction/matrix/hessenberg.rs

@@ -0,0 +1,149 @@
+use crate::solvers::SolverError;
+use crate::utils::Arr2D;
+
+/// Computes the Hessenberg reduction of a square matrix A.
+/// Returns (H, Q) such that H = Q^T * A * Q, where H is upper Hessenberg and Q is orthogonal.
+pub fn hessenberg_reduction(matrix: &Arr2D<f64>) -> Result<(Arr2D<f64>, Arr2D<f64>), SolverError> {
+    if matrix.height != matrix.width {
+        return Err(SolverError::NonSquareMatrix);
+    }
+    let n = matrix.height;
+    if n <= 2 {
+        return Ok((matrix.clone(), Arr2D::identity(n)));
+    }
+
+    let mut h = matrix.clone();
+    let mut q = Arr2D::identity(n);
+
+    for k in 0..n - 2 {
+        // x = h[k+1..n, k]
+        let mut x = Vec::with_capacity(n - (k + 1));
+        for i in k + 1..n {
+            x.push(h[(i, k)]);
+        }
+
+        let norm_x = x.iter().map(|&val| val * val).sum::<f64>().sqrt();
+        if norm_x == 0.0 {
+            continue;
+        }
+
+        // v = x + sign(x[0]) * ||x|| * e1
+        let sign = if x[0] >= 0.0 { 1.0 } else { -1.0 };
+        let mut v = x;
+        v[0] += sign * norm_x;
+
+        // normalize v
+        let norm_v = v.iter().map(|&val| val * val).sum::<f64>().sqrt();
+        if norm_v == 0.0 {
+            continue;
+        }
+        for val in v.iter_mut() {
+            *val /= norm_v;
+        }
+
+        // Apply H_k = I - 2vv^T to A from the left: A = H_k * A
+        // This affects rows k+1..n
+        for col in k..n {
+            let mut dot = 0.0;
+            for i in 0..v.len() {
+                dot += v[i] * h[(k + 1 + i, col)];
+            }
+            for i in 0..v.len() {
+                h[(k + 1 + i, col)] -= 2.0 * v[i] * dot;
+            }
+        }
+
+        // Apply H_k to A from the right: A = A * H_k
+        // This affects columns k+1..n
+        for row in 0..n {
+            let mut dot = 0.0;
+            for j in 0..v.len() {
+                dot += v[j] * h[(row, k + 1 + j)];
+            }
+            for j in 0..v.len() {
+                h[(row, k + 1 + j)] -= 2.0 * v[j] * dot;
+            }
+        }
+
+        // Accumulate Q = Q * H_k
+        // This affects columns k+1..n of Q
+        for row in 0..n {
+            let mut dot = 0.0;
+            for j in 0..v.len() {
+                dot += v[j] * q[(row, k + 1 + j)];
+            }
+            for j in 0..v.len() {
+                q[(row, k + 1 + j)] -= 2.0 * v[j] * dot;
+            }
+        }
+    }
+
+    Ok((h, q))
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::utils::Rounding;
+
+    #[test]
+    fn test_hessenberg_2x2() {
+        let mat = Arr2D::from(&[[1.0, 2.0], [3.0, 4.0]]);
+        let (h, q) = hessenberg_reduction(&mat).unwrap();
+        assert_eq!(h.round_to_decimal(10), mat.round_to_decimal(10));
+        assert_eq!(
+            q.round_to_decimal(10),
+            Arr2D::identity(2).round_to_decimal(10)
+        );
+    }
+
+    #[test]
+    fn test_hessenberg_3x3() {
+        // Simple 3x3 matrix
+        let mat = Arr2D::from(&[[1.0, 5.0, 7.0], [3.0, 0.0, 6.0], [4.0, 3.0, 1.0]]);
+        let (h, q) = hessenberg_reduction(&mat).unwrap();
+
+        // Check if H is upper Hessenberg
+        assert!(h[(2, 0)].abs() < 1e-10);
+
+        // Check if Q is orthogonal: Q * Q^T = I
+        let qt = q.transpose();
+        let i = q.dot(&qt).unwrap();
+        assert_eq!(
+            i.round_to_decimal(10),
+            Arr2D::identity(3).round_to_decimal(10)
+        );
+
+        // Check if A = Q * H * Q^T
+        let reconstructed = q.dot(&h).unwrap().dot(&qt).unwrap();
+        assert_eq!(reconstructed.round_to_decimal(10), mat.round_to_decimal(10));
+    }
+
+    #[test]
+    fn test_hessenberg_4x4() {
+        let mat = Arr2D::from(&[
+            [1.0, 2.0, 3.0, 4.0],
+            [2.0, 1.0, 2.0, 3.0],
+            [3.0, 2.0, 1.0, 2.0],
+            [4.0, 3.0, 2.0, 1.0],
+        ]);
+        let (h, q) = hessenberg_reduction(&mat).unwrap();
+
+        // Check upper Hessenberg form (zeros below subdiagonal)
+        assert!(h[(2, 0)].abs() < 1e-10);
+        assert!(h[(3, 0)].abs() < 1e-10);
+        assert!(h[(3, 1)].abs() < 1e-10);
+
+        // Check orthogonality
+        let qt = q.transpose();
+        let i = q.dot(&qt).unwrap();
+        assert_eq!(
+            i.round_to_decimal(10),
+            Arr2D::identity(4).round_to_decimal(10)
+        );
+
+        // Check A = Q H Q^T
+        let reconstructed = q.dot(&h).unwrap().dot(&qt).unwrap();
+        assert_eq!(reconstructed.round_to_decimal(10), mat.round_to_decimal(10));
+    }
+}

spindalis/src/reduction/matrix/mod.rs

@@ -1 +1,2 @@
-
+pub mod hessenberg;
+pub use hessenberg::hessenberg_reduction;