feat: Heisenberg matrix reduction and workspace refactor

.github/workflows/clippy-lint.yml

@@ -16,12 +16,9 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        crate-target:
-          - { crate: spindalis }
-          - { crate: spindalis_core }
-          - { crate: spindalis_macros }
+        package: [spindalis, spindalis_core, 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 --all-targets --all-features -p ${{ matrix.package }}

.github/workflows/rust_tests.yml

@@ -10,22 +10,16 @@ env:
   CARGO_TERM_COLOR: always
 
 jobs:
-  build:
-
+  build_and_test:
     runs-on: ubuntu-latest
     strategy:
       fail-fast: false
       matrix:
-        crate-target:
-          - { crate: spindalis }
-          - { crate: spindalis_core }
-          - { crate: spindalis_macros }
+        package: [spindalis, spindalis_core, spindalis_macros]
 
     steps:
-    - uses: actions/checkout@v4
-    - name: Build
-      working-directory: ./${{ matrix.crate-target.crate }}
-      run: cargo build --verbose
-    - name: Run tests
-      working-directory: ./${{ matrix.crate-target.crate }}
-      run: cargo test --verbose
+      - uses: actions/checkout@v5
+      - name: Build
+        run: cargo build --verbose -p ${{ matrix.package }}
+      - name: Run tests
+        run: cargo test --verbose -p ${{ matrix.package }}

.justfile

@@ -4,56 +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" \
-}
-
 
 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}}
+    cargo test --workspace
 
 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}}
+    cargo clippy --workspace --all-targets --all-features
 
 format FORMAT:
-    cargo fmt --mainfest-path {{ROOT_DIR}}/{{FORMAT}}/Cargo.toml
+    cargo fmt -p {{FORMAT}}
 
 [doc('Format all modules')]
 format-all:
-    {{fmt_cmd}}
+    cargo fmt --all
 
 latest:
     git pull origin main

Cargo.toml

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

spindalis/src/reduction/matrix/hessenberg.rs

@@ -0,0 +1,148 @@
+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 = 0.0;
+        for i in k + 1..n {
+            x += h[(i, k)] * h[(i, k)];
+        }
+        let norm_x = x.sqrt();
+
+        if norm_x == 0.0 {
+            continue;
+        }
+
+        // v = x + sign(x[0]) * ||x|| * e1
+        let h_first = h[(k + 1, k)];
+        let sign = if h_first >= 0.0 { -1.0 } else { 1.0 };
+        let u1 = h_first - sign * norm_x;
+        // v[0] is implicitly 1.0, we store the rest of v in a temporary slice or vec
+        let mut v = vec![0.0; n - (k + 1)];
+        v[0] = 1.0;
+        for i in 1..v.len() {
+            v[i] = h[(k + 1 + i, k)] / u1;
+        }
+
+        let tau = -sign * u1 / norm_x;
+
+        // Apply H_k = I - tau * v * v^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)] -= tau * 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)] -= tau * 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)] -= tau * 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;

spindalis_macros/src/lib.rs

@@ -41,7 +41,8 @@ pub fn parse_intermediate_polynomial(input: TokenStream) -> TokenStream {
         }
     };
 
-    let mut tokens = String::from("::spindalis_core::polynomials::structs::IntermediatePolynomial { ");
+    let mut tokens =
+        String::from("::spindalis_core::polynomials::structs::IntermediatePolynomial { ");
     tokens.push_str("terms: vec![");
     for term in output.terms {
         tokens.push_str(&format!(