diff --git a/.github/workflows/clippy-lint.yml b/.github/workflows/clippy-lint.yml index 7d8d015..5eaf7b2 100644 --- a/.github/workflows/clippy-lint.yml +++ b/.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 }} diff --git a/.github/workflows/rust_tests.yml b/.github/workflows/rust_tests.yml index 0da91f3..e64baeb 100644 --- a/.github/workflows/rust_tests.yml +++ b/.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 }} diff --git a/.justfile b/.justfile index 9e5eb55..b4c223a 100644 --- a/.justfile +++ b/.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 diff --git a/Cargo.lock b/Cargo.lock new file mode 100644 index 0000000..ef6361d --- /dev/null +++ b/Cargo.lock @@ -0,0 +1,47 @@ +# This file is automatically @generated by Cargo. +# It is not intended for manual editing. +version = 4 + +[[package]] +name = "proc-macro2" +version = "1.0.105" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "535d180e0ecab6268a3e718bb9fd44db66bbbc256257165fc699dadf70d16fe7" +dependencies = [ + "unicode-ident", +] + +[[package]] +name = "quote" +version = "1.0.43" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dc74d9a594b72ae6656596548f56f667211f8a97b3d4c3d467150794690dc40a" +dependencies = [ + "proc-macro2", +] + +[[package]] +name = "spindalis" +version = "0.6.0" +dependencies = [ + "spindalis_core", + "spindalis_macros", +] + +[[package]] +name = "spindalis_core" +version = "0.3.0" + +[[package]] +name = "spindalis_macros" +version = "0.1.2" +dependencies = [ + "quote", + "spindalis_core", +] + +[[package]] +name = "unicode-ident" +version = "1.0.22" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9312f7c4f6ff9069b165498234ce8be658059c6728633667c526e27dc2cf1df5" diff --git a/Cargo.toml b/Cargo.toml new file mode 100644 index 0000000..3ddb9bd --- /dev/null +++ b/Cargo.toml @@ -0,0 +1,7 @@ +[workspace] +members = [ + "spindalis", + "spindalis_core", + "spindalis_macros", +] +resolver = "2" diff --git a/spindalis/src/reduction/matrix/hessenberg.rs b/spindalis/src/reduction/matrix/hessenberg.rs new file mode 100644 index 0000000..7fd12b0 --- /dev/null +++ b/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) -> Result<(Arr2D, Arr2D), 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)); + } +} diff --git a/spindalis/src/reduction/matrix/mod.rs b/spindalis/src/reduction/matrix/mod.rs index 8b13789..989e145 100644 --- a/spindalis/src/reduction/matrix/mod.rs +++ b/spindalis/src/reduction/matrix/mod.rs @@ -1 +1,2 @@ - +pub mod hessenberg; +pub use hessenberg::hessenberg_reduction; diff --git a/spindalis_macros/src/lib.rs b/spindalis_macros/src/lib.rs index 5e5b31f..63ede1c 100644 --- a/spindalis_macros/src/lib.rs +++ b/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!(