From 32e4b0d2957f62798eb341ee4813ab0440380620 Mon Sep 17 00:00:00 2001 From: Chris Date: Thu, 18 Jun 2026 14:43:23 +1200 Subject: [PATCH 1/2] chore: move matrix code to custom --- Cargo.lock | 285 ------------------ crates/cnvx-math/Cargo.toml | 1 - crates/cnvx-math/src/lib.rs | 6 +- crates/cnvx-math/src/matrix/dense.rs | 287 +++++++++++++++++-- crates/cnvx-math/src/matrix/mod.rs | 136 ++++++--- crates/cnvx-math/src/matrix/sparse.rs | 73 ++++- crates/cnvx-math/tests/dense_matrix_tests.rs | 130 +++++++++ 7 files changed, 567 insertions(+), 351 deletions(-) create mode 100644 crates/cnvx-math/tests/dense_matrix_tests.rs diff --git a/Cargo.lock b/Cargo.lock index 4e6301e..c1481bf 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -52,27 +52,12 @@ dependencies = [ "windows-sys 0.61.2", ] -[[package]] -name = "approx" -version = "0.5.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "cab112f0a86d568ea0e627cc1d6be74a1e9cd55214684db5561995f6dad897c6" -dependencies = [ - "num-traits", -] - [[package]] name = "atomic-waker" version = "1.1.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "1505bd5d3d116872e7271a6d4e16d81d0c8570876c8de68093a09ac269d8aac0" -[[package]] -name = "autocfg" -version = "1.5.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "c08606f8c3cbf4ce6ec8e28fb0014a2c086708fe954eaa885384a6165172e7e8" - [[package]] name = "aws-lc-rs" version = "1.15.4" @@ -113,12 +98,6 @@ version = "3.19.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "5dd9dc738b7a8311c7ade152424974d8115f2cdad61e8dab8dac9f2362298510" -[[package]] -name = "bytemuck" -version = "1.25.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "c8efb64bd706a16a1bdde310ae86b351e4d21550d98d056f22f8a7f7a2183fec" - [[package]] name = "bytes" version = "1.11.1" @@ -248,9 +227,6 @@ dependencies = [ [[package]] name = "cnvx-math" version = "0.0.1" -dependencies = [ - "nalgebra", -] [[package]] name = "cnvx-parse" @@ -491,114 +467,6 @@ dependencies = [ "wasm-bindgen", ] -[[package]] -name = "glam" -version = "0.14.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "333928d5eb103c5d4050533cec0384302db6be8ef7d3cebd30ec6a35350353da" - -[[package]] -name = "glam" -version = "0.15.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "3abb554f8ee44336b72d522e0a7fe86a29e09f839a36022fa869a7dfe941a54b" - -[[package]] -name = "glam" -version = "0.16.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "4126c0479ccf7e8664c36a2d719f5f2c140fbb4f9090008098d2c291fa5b3f16" - -[[package]] -name = "glam" -version = "0.17.3" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e01732b97afd8508eee3333a541b9f7610f454bb818669e66e90f5f57c93a776" - -[[package]] -name = "glam" -version = "0.18.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "525a3e490ba77b8e326fb67d4b44b4bd2f920f44d4cc73ccec50adc68e3bee34" - -[[package]] -name = "glam" -version = "0.19.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "2b8509e6791516e81c1a630d0bd7fbac36d2fa8712a9da8662e716b52d5051ca" - -[[package]] -name = "glam" -version = "0.20.5" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f43e957e744be03f5801a55472f593d43fabdebf25a4585db250f04d86b1675f" - -[[package]] -name = "glam" -version = "0.21.3" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "518faa5064866338b013ff9b2350dc318e14cc4fcd6cb8206d7e7c9886c98815" - -[[package]] -name = "glam" -version = "0.22.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "12f597d56c1bd55a811a1be189459e8fad2bbc272616375602443bdfb37fa774" - -[[package]] -name = "glam" -version = "0.23.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "8e4afd9ad95555081e109fe1d21f2a30c691b5f0919c67dfa690a2e1eb6bd51c" - -[[package]] -name = "glam" -version = "0.24.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b5418c17512bdf42730f9032c74e1ae39afc408745ebb2acf72fbc4691c17945" - -[[package]] -name = "glam" -version = "0.25.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "151665d9be52f9bb40fc7966565d39666f2d1e69233571b71b87791c7e0528b3" - -[[package]] -name = "glam" -version = "0.27.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "9e05e7e6723e3455f4818c7b26e855439f7546cf617ef669d1adedb8669e5cb9" - -[[package]] -name = "glam" -version = "0.28.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "779ae4bf7e8421cf91c0b3b64e7e8b40b862fba4d393f59150042de7c4965a94" - -[[package]] -name = "glam" -version = "0.29.3" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "8babf46d4c1c9d92deac9f7be466f76dfc4482b6452fc5024b5e8daf6ffeb3ee" - -[[package]] -name = "glam" -version = "0.30.10" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "19fc433e8437a212d1b6f1e68c7824af3aed907da60afa994e7f542d18d12aa9" - -[[package]] -name = "glam" -version = "0.31.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "556f6b2ea90b8d15a74e0e7bb41671c9bdf38cd9f78c284d750b9ce58a2b5be7" - -[[package]] -name = "glam" -version = "0.32.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f70749695b063ecbf6b62949ccccde2e733ec3ecbbd71d467dca4e5c6c97cca0" - [[package]] name = "h2" version = "0.4.13" @@ -944,16 +812,6 @@ version = "0.1.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "112b39cec0b298b6c1999fee3e31427f74f676e4cb9879ed1a121b43661a4154" -[[package]] -name = "matrixmultiply" -version = "0.3.10" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a06de3016e9fae57a36fd14dba131fccf49f74b40b7fbdb472f96e361ec71a08" -dependencies = [ - "autocfg", - "rawpointer", -] - [[package]] name = "memchr" version = "2.8.0" @@ -977,99 +835,6 @@ dependencies = [ "windows-sys 0.61.2", ] -[[package]] -name = "nalgebra" -version = "0.34.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "df76ea0ff5c7e6b88689085804d6132ded0ddb9de5ca5b8aeb9eeadc0508a70a" -dependencies = [ - "approx", - "glam 0.14.0", - "glam 0.15.2", - "glam 0.16.0", - "glam 0.17.3", - "glam 0.18.0", - "glam 0.19.0", - "glam 0.20.5", - "glam 0.21.3", - "glam 0.22.0", - "glam 0.23.0", - "glam 0.24.2", - "glam 0.25.0", - "glam 0.27.0", - "glam 0.28.0", - "glam 0.29.3", - "glam 0.30.10", - "glam 0.31.1", - "glam 0.32.1", - "matrixmultiply", - "nalgebra-macros", - "num-complex", - "num-rational", - "num-traits", - "simba", - "typenum", -] - -[[package]] -name = "nalgebra-macros" -version = "0.3.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "973e7178a678cfd059ccec50887658d482ce16b0aa9da3888ddeab5cd5eb4889" -dependencies = [ - "proc-macro2", - "quote", - "syn", -] - -[[package]] -name = "num-bigint" -version = "0.4.6" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a5e44f723f1133c9deac646763579fdb3ac745e418f2a7af9cd0c431da1f20b9" -dependencies = [ - "num-integer", - "num-traits", -] - -[[package]] -name = "num-complex" -version = "0.4.6" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "73f88a1307638156682bada9d7604135552957b7818057dcef22705b4d509495" -dependencies = [ - "num-traits", -] - -[[package]] -name = "num-integer" -version = "0.1.46" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "7969661fd2958a5cb096e56c8e1ad0444ac2bbcd0061bd28660485a44879858f" -dependencies = [ - "num-traits", -] - -[[package]] -name = "num-rational" -version = "0.4.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f83d14da390562dca69fc84082e73e548e1ad308d24accdedd2720017cb37824" -dependencies = [ - "num-bigint", - "num-integer", - "num-traits", -] - -[[package]] -name = "num-traits" -version = "0.2.19" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841" -dependencies = [ - "autocfg", -] - [[package]] name = "once_cell" version = "1.21.3" @@ -1088,12 +853,6 @@ version = "0.2.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "7c87def4c32ab89d880effc9e097653c8da5d6ef28e6b539d313baaacfbafcbe" -[[package]] -name = "paste" -version = "1.0.15" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "57c0d7b74b563b49d38dae00a0c37d4d6de9b432382b2892f0574ddcae73fd0a" - [[package]] name = "percent-encoding" version = "2.3.2" @@ -1303,12 +1062,6 @@ dependencies = [ "getrandom 0.3.4", ] -[[package]] -name = "rawpointer" -version = "0.2.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3" - [[package]] name = "reqwest" version = "0.13.2" @@ -1463,15 +1216,6 @@ version = "1.0.22" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "b39cdef0fa800fc44525c84ccb54a029961a8215f9619753635a9c0d2538d46d" -[[package]] -name = "safe_arch" -version = "0.7.4" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "96b02de82ddbe1b636e6170c21be622223aea188ef2e139be0a5b219ec215323" -dependencies = [ - "bytemuck", -] - [[package]] name = "same-file" version = "1.0.6" @@ -1571,19 +1315,6 @@ dependencies = [ "libc", ] -[[package]] -name = "simba" -version = "0.9.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "c99284beb21666094ba2b75bbceda012e610f5479dfcc2d6e2426f53197ffd95" -dependencies = [ - "approx", - "num-complex", - "num-traits", - "paste", - "wide", -] - [[package]] name = "slab" version = "0.4.12" @@ -1900,12 +1631,6 @@ version = "0.2.5" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "e421abadd41a4225275504ea4d6566923418b7f05506fbc9c0fe86ba7396114b" -[[package]] -name = "typenum" -version = "1.19.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "562d481066bde0658276a35467c4af00bdc6ee726305698a55b86e61d7ad82bb" - [[package]] name = "unicode-ident" version = "1.0.22" @@ -2085,16 +1810,6 @@ dependencies = [ "rustls-pki-types", ] -[[package]] -name = "wide" -version = "0.7.33" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "0ce5da8ecb62bcd8ec8b7ea19f69a51275e91299be594ea5cc6ef7819e16cd03" -dependencies = [ - "bytemuck", - "safe_arch", -] - [[package]] name = "winapi-util" version = "0.1.11" diff --git a/crates/cnvx-math/Cargo.toml b/crates/cnvx-math/Cargo.toml index f2918e2..fc7586d 100644 --- a/crates/cnvx-math/Cargo.toml +++ b/crates/cnvx-math/Cargo.toml @@ -10,7 +10,6 @@ keywords = { workspace = true } readme = { workspace = true } [dependencies] -nalgebra = { workspace = true } [lints] workspace = true diff --git a/crates/cnvx-math/src/lib.rs b/crates/cnvx-math/src/lib.rs index 5769665..f065436 100644 --- a/crates/cnvx-math/src/lib.rs +++ b/crates/cnvx-math/src/lib.rs @@ -1,8 +1,8 @@ //! # CNVX Math //! //! Linear algebra utilities for LP solvers and numerical algorithms. -//! Provides matrix types and traits used in simplex computations and -//! other numerical routines. +//! Provides native matrix types and traits used in simplex computations and +//! other numerical routines without relying on external heavy dependencies. //! //! # Modules //! @@ -10,4 +10,4 @@ pub mod matrix; -pub use matrix::{DenseMatrix, MatrixWrapper as Matrix, SparseMatrix}; +pub use matrix::{DenseMatrix, Matrix, SparseMatrix}; diff --git a/crates/cnvx-math/src/matrix/dense.rs b/crates/cnvx-math/src/matrix/dense.rs index e30e2dd..e592fd2 100644 --- a/crates/cnvx-math/src/matrix/dense.rs +++ b/crates/cnvx-math/src/matrix/dense.rs @@ -1,45 +1,292 @@ -use nalgebra::DMatrix; +use crate::matrix::Matrix; -use crate::matrix::MatrixWrapper; +/// A dense matrix stored in row-major order using a flat 1D vector. +#[derive(Debug, Clone, PartialEq)] +pub struct DenseMatrix { + rows: usize, + cols: usize, + data: Vec, +} -#[derive(Debug, Clone)] -pub struct ExposedDenseMatrix { - pub inner: DMatrix, +impl DenseMatrix { + /// Helper to calculate the 1D index from 2D coordinates. + #[inline(always)] + fn index(&self, row: usize, col: usize) -> usize { + row * self.cols + col + } } -impl MatrixWrapper for ExposedDenseMatrix { +impl Matrix for DenseMatrix { fn new(rows: usize, cols: usize) -> Self { - Self { inner: DMatrix::zeros(rows, cols) } + Self { rows, cols, data: vec![0.0; rows * cols] } } + #[inline(always)] fn rows(&self) -> usize { - self.inner.nrows() + self.rows } + #[inline(always)] fn cols(&self) -> usize { - self.inner.ncols() + self.cols } + #[inline(always)] fn get(&self, row: usize, col: usize) -> f64 { - self.inner[(row, col)] + assert!(row < self.rows && col < self.cols, "Index out of bounds"); + self.data[self.index(row, col)] } + #[inline(always)] fn set(&mut self, row: usize, col: usize, value: f64) { - self.inner[(row, col)] = value; + assert!(row < self.rows && col < self.cols, "Index out of bounds"); + let idx = self.index(row, col); + self.data[idx] = value; + } + + fn add(&self, other: &Self) -> Result { + if self.rows != other.rows || self.cols != other.cols { + return Err("Matrix dimensions must match for addition".to_string()); + } + let data = self.data.iter().zip(other.data.iter()).map(|(a, b)| a + b).collect(); + Ok(Self { rows: self.rows, cols: self.cols, data }) + } + + fn sub(&self, other: &Self) -> Result { + if self.rows != other.rows || self.cols != other.cols { + return Err("Matrix dimensions must match for subtraction".to_string()); + } + let data = self.data.iter().zip(other.data.iter()).map(|(a, b)| a - b).collect(); + Ok(Self { rows: self.rows, cols: self.cols, data }) + } + + fn mul(&self, other: &Self) -> Result { + if self.cols != other.rows { + return Err("Incompatible dimensions for matrix multiplication".to_string()); + } + + let mut result = Self::new(self.rows, other.cols); + + // IKJ loop ordering for cache-friendly memory access + for i in 0..self.rows { + for k in 0..self.cols { + let a_ik = self.get(i, k); + for j in 0..other.cols { + let b_kj = other.get(k, j); + let current = result.get(i, j); + result.set(i, j, current + a_ik * b_kj); + } + } + } + Ok(result) + } + + fn mul_vec(&self, rhs: &[f64]) -> Result, String> { + if self.cols != rhs.len() { + return Err("Matrix columns must match vector length".to_string()); + } + let mut result = vec![0.0; self.rows]; + for i in 0..self.rows { + let mut sum = 0.0; + for j in 0..self.cols { + sum += self.get(i, j) * rhs[j]; + } + result[i] = sum; + } + Ok(result) + } + + fn add_scalar(&self, scalar: f64) -> Self { + let data = self.data.iter().map(|&x| x + scalar).collect(); + Self { rows: self.rows, cols: self.cols, data } + } + + fn sub_scalar(&self, scalar: f64) -> Self { + let data = self.data.iter().map(|&x| x - scalar).collect(); + Self { rows: self.rows, cols: self.cols, data } + } + + fn mul_scalar(&self, scalar: f64) -> Self { + let data = self.data.iter().map(|&x| x * scalar).collect(); + Self { rows: self.rows, cols: self.cols, data } + } + + fn div_scalar(&self, scalar: f64) -> Self { + let data = self.data.iter().map(|&x| x / scalar).collect(); + Self { rows: self.rows, cols: self.cols, data } + } + + fn mul_elementwise(&self, other: &Self) -> Result { + if self.rows != other.rows || self.cols != other.cols { + return Err("Matrix dimensions must match for element-wise multiplication" + .to_string()); + } + let data = self + .data + .iter() + .zip(other.data.iter()) + .map(|(&a, &b)| a * b) + .collect(); + Ok(Self { rows: self.rows, cols: self.cols, data }) + } + + fn get_row(&self, row: usize) -> Vec { + assert!(row < self.rows, "Row index out of bounds"); + let start = row * self.cols; + self.data[start..start + self.cols].to_vec() + } + + fn get_col(&self, col: usize) -> Vec { + assert!(col < self.cols, "Column index out of bounds"); + (0..self.rows).map(|i| self.data[i * self.cols + col]).collect() + } + + fn set_row(&mut self, row: usize, vec: &[f64]) -> Result<(), String> { + if row >= self.rows { + return Err("Row index out of bounds".to_string()); + } + if vec.len() != self.cols { + return Err("Vector length must match matrix columns".to_string()); + } + let start = row * self.cols; + self.data[start..start + self.cols].copy_from_slice(vec); + Ok(()) + } + + fn set_col(&mut self, col: usize, vec: &[f64]) -> Result<(), String> { + if col >= self.cols { + return Err("Column index out of bounds".to_string()); + } + if vec.len() != self.rows { + return Err("Vector length must match matrix rows".to_string()); + } + for i in 0..self.rows { + self.data[i * self.cols + col] = vec[i]; + } + Ok(()) + } + + fn swap_rows(&mut self, row1: usize, row2: usize) { + assert!(row1 < self.rows && row2 < self.rows, "Row index out of bounds"); + if row1 == row2 { + return; + } + // Ensure row1 is the smaller index to safely split the slice + let (r1, r2) = if row1 < row2 { (row1, row2) } else { (row2, row1) }; + let (first, second) = self.data.split_at_mut(r2 * self.cols); + + let row1_slice = &mut first[r1 * self.cols..(r1 + 1) * self.cols]; + let row2_slice = &mut second[0..self.cols]; + row1_slice.swap_with_slice(row2_slice); + } + + fn norm_inf(&self) -> f64 { + let mut max_sum = 0.0_f64; + for i in 0..self.rows { + let start = i * self.cols; + let row_sum: f64 = + self.data[start..start + self.cols].iter().map(|&x| x.abs()).sum(); + if row_sum > max_sum { + max_sum = row_sum; + } + } + max_sum + } + + fn transpose(&self) -> Self { + let mut result = Self::new(self.cols, self.rows); + for i in 0..self.rows { + for j in 0..self.cols { + result.set(j, i, self.get(i, j)); + } + } + result + } + + fn identity(size: usize) -> Self { + let mut m = Self::new(size, size); + for i in 0..size { + m.set(i, i, 1.0); + } + m + } + + fn from_diagonal(vec: &[f64]) -> Self { + let size = vec.len(); + let mut m = Self::new(size, size); + for i in 0..size { + m.set(i, i, vec[i]); + } + m + } + + fn diagonal(&self) -> Vec { + let size = self.rows.min(self.cols); + (0..size).map(|i| self.get(i, i)).collect() } - // TODO: This is a very naive implementation, and will need to be improved for performance depending on the matrix shape and sparsity. fn mldivide(&self, rhs: &mut [f64]) -> Result<(), String> { - let a = &self.inner; - let b = DMatrix::from_column_slice(rhs.len(), 1, rhs); - match a.clone().lu().solve(&b) { - Some(solution) => { - for i in 0..rhs.len() { - rhs[i] = solution[(i, 0)]; + if self.rows != self.cols { + return Err("Matrix must be square to solve linear systems".to_string()); + } + if self.rows != rhs.len() { + return Err("RHS vector length must match matrix dimensions".to_string()); + } + + let n = self.rows; + let mut lu = self.data.clone(); + let mut p: Vec = (0..n).collect(); + + // LU Decomposition with Partial Pivoting + for i in 0..n { + let mut max_a = 0.0; + let mut imax = i; + for k in i..n { + let abs_a = lu[k * n + i].abs(); + if abs_a > max_a { + max_a = abs_a; + imax = k; + } + } + + if max_a < 1e-12 { + return Err("Matrix is singular or nearly singular".to_string()); + } + + if imax != i { + for k in 0..n { + lu.swap(i * n + k, imax * n + k); } - Ok(()) + p.swap(i, imax); + } + + for j in (i + 1)..n { + lu[j * n + i] /= lu[i * n + i]; + for k in (i + 1)..n { + let factor = lu[j * n + i] * lu[i * n + k]; + lu[j * n + k] -= factor; + } + } + } + + // Forward substitution + let mut x = vec![0.0; n]; + for i in 0..n { + x[i] = rhs[p[i]]; + for k in 0..i { + x[i] -= lu[i * n + k] * x[k]; } - None => Err("Matrix is singular or not square".to_string()), } + + // Backward substitution + for i in (0..n).rev() { + for k in (i + 1)..n { + x[i] -= lu[i * n + k] * x[k]; + } + x[i] /= lu[i * n + i]; + } + + rhs.copy_from_slice(&x); + Ok(()) } } diff --git a/crates/cnvx-math/src/matrix/mod.rs b/crates/cnvx-math/src/matrix/mod.rs index d878895..03874be 100644 --- a/crates/cnvx-math/src/matrix/mod.rs +++ b/crates/cnvx-math/src/matrix/mod.rs @@ -1,27 +1,18 @@ mod dense; mod sparse; -pub use dense::ExposedDenseMatrix as DenseMatrix; -pub use sparse::ExposedSparseMatrix as SparseMatrix; +pub use dense::DenseMatrix; +pub use sparse::SparseMatrix; /// A generic matrix trait for linear algebra operations. /// /// This trait defines the interface for matrices used in the LP solver, -/// including creation, element access, and solving linear systems. -/// -/// Implementors must provide row/column indexing, element access, and -/// a method for solving linear systems. -pub trait MatrixWrapper: Clone { +/// including creation, element access, arithmetic, and solving linear systems. +pub trait Matrix: Clone { + // --- 1. Core Constructors and Accessors --- + /// Create a new matrix with the given number of rows and columns, /// initialized with zeros. - /// - /// # Example - /// ``` - /// # use cnvx_math::{DenseMatrix, Matrix}; - /// let m = DenseMatrix::new(2, 3); - /// assert_eq!(m.rows(), 2); - /// assert_eq!(m.cols(), 3); - /// ``` fn new(rows: usize, cols: usize) -> Self where Self: Sized; @@ -44,29 +35,106 @@ pub trait MatrixWrapper: Clone { /// Panics if `row` or `col` are out of bounds. fn set(&mut self, row: usize, col: usize, value: f64); + // --- 2. Standard Arithmetic --- + + /// Adds another matrix to this matrix. + fn add(&self, other: &Self) -> Result + where + Self: Sized; + + /// Subtracts another matrix from this matrix. + fn sub(&self, other: &Self) -> Result + where + Self: Sized; + + /// Multiplies this matrix by another matrix. + fn mul(&self, other: &Self) -> Result + where + Self: Sized; + + /// Multiplies this matrix by a column vector. + fn mul_vec(&self, rhs: &[f64]) -> Result, String>; + + // --- 3. Scalar and Element-wise Operations --- + + /// Adds a scalar to every element in the matrix. + fn add_scalar(&self, scalar: f64) -> Self + where + Self: Sized; + + /// Subtracts a scalar from every element in the matrix. + fn sub_scalar(&self, scalar: f64) -> Self + where + Self: Sized; + + /// Multiplies every element in the matrix by a scalar. + fn mul_scalar(&self, scalar: f64) -> Self + where + Self: Sized; + + /// Divides every element in the matrix by a scalar. + fn div_scalar(&self, scalar: f64) -> Self + where + Self: Sized; + + /// Performs element-wise (Hadamard) multiplication with another matrix. + fn mul_elementwise(&self, other: &Self) -> Result + where + Self: Sized; + + // --- 4. Row and Column Manipulations --- + + /// Extracts a specific row as a standard vector. + fn get_row(&self, row: usize) -> Vec; + + /// Extracts a specific column as a standard vector. + fn get_col(&self, col: usize) -> Vec; + + /// Overwrites a specific row with a new vector. + fn set_row(&mut self, row: usize, vec: &[f64]) -> Result<(), String>; + + /// Overwrites a specific column with a new vector. + fn set_col(&mut self, col: usize, vec: &[f64]) -> Result<(), String>; + + /// Swaps two rows in place. + fn swap_rows(&mut self, row1: usize, row2: usize); + + // --- 5. Norms and Convergence Checks --- + + /// Calculates the L_infinity norm (maximum absolute row sum). + fn norm_inf(&self) -> f64; + + // --- 6. Constructors and Utilities --- + + /// Returns the transpose of the matrix. + fn transpose(&self) -> Self + where + Self: Sized; + + /// Creates a square identity matrix of the given size. + fn identity(size: usize) -> Self + where + Self: Sized; + + /// Creates a square diagonal matrix from a vector of diagonal elements. + fn from_diagonal(vec: &[f64]) -> Self + where + Self: Sized; + + /// Extracts the diagonal elements of the matrix as a vector. + fn diagonal(&self) -> Vec; + + // --- 7. Solvers --- + /// Solve a square linear system `Ax = rhs` /// /// On success, `rhs` is overwritten with the solution vector `x`. - /// - /// # Errors - /// Returns an `Err(String)` if the system cannot be solved, e.g., if the matrix - /// is singular. - /// - /// # Example - /// ``` - /// # use cnvx_math::{DenseMatrix, Matrix}; - /// let mut a = DenseMatrix::new(2, 2); - /// a.set(0, 0, 2.0); - /// a.set(0, 1, 1.0); - /// a.set(1, 0, 1.0); - /// a.set(1, 1, 3.0); - /// - /// let mut rhs = vec![3.0, 7.0]; - /// a.mldivide(&mut rhs).unwrap(); - /// assert!((rhs[0] - 0.4).abs() < 1e-6); - /// assert!((rhs[1] - 2.2).abs() < 1e-6); - /// ``` fn mldivide(&self, rhs: &mut [f64]) -> Result<(), String> where Self: Sized; } + +/// Helper function to calculate the L2 (Euclidean) norm of a standard vector. +pub fn vector_norm_l2(vec: &[f64]) -> f64 { + vec.iter().map(|&x| x * x).sum::().sqrt() +} diff --git a/crates/cnvx-math/src/matrix/sparse.rs b/crates/cnvx-math/src/matrix/sparse.rs index efe6207..b5d9996 100644 --- a/crates/cnvx-math/src/matrix/sparse.rs +++ b/crates/cnvx-math/src/matrix/sparse.rs @@ -1,31 +1,88 @@ -use crate::matrix::MatrixWrapper; +use crate::matrix::Matrix; -// TODO: This will either use sprs or nalgebra_sparse csr format #[derive(Debug, Clone)] -pub struct ExposedSparseMatrix {} +pub struct SparseMatrix {} #[allow(unused)] -impl MatrixWrapper for ExposedSparseMatrix { +impl Matrix for SparseMatrix { fn new(rows: usize, cols: usize) -> Self { unimplemented!() } - fn rows(&self) -> usize { unimplemented!() } - fn cols(&self) -> usize { unimplemented!() } - fn get(&self, row: usize, col: usize) -> f64 { unimplemented!() } - fn set(&mut self, row: usize, col: usize, value: f64) { unimplemented!() } + fn add(&self, other: &Self) -> Result { + unimplemented!() + } + fn sub(&self, other: &Self) -> Result { + unimplemented!() + } + fn mul(&self, other: &Self) -> Result { + unimplemented!() + } + fn mul_vec(&self, rhs: &[f64]) -> Result, String> { + unimplemented!() + } + + fn add_scalar(&self, scalar: f64) -> Self { + unimplemented!() + } + fn sub_scalar(&self, scalar: f64) -> Self { + unimplemented!() + } + fn mul_scalar(&self, scalar: f64) -> Self { + unimplemented!() + } + fn div_scalar(&self, scalar: f64) -> Self { + unimplemented!() + } + fn mul_elementwise(&self, other: &Self) -> Result { + unimplemented!() + } + + fn get_row(&self, row: usize) -> Vec { + unimplemented!() + } + fn get_col(&self, col: usize) -> Vec { + unimplemented!() + } + fn set_row(&mut self, row: usize, vec: &[f64]) -> Result<(), String> { + unimplemented!() + } + fn set_col(&mut self, col: usize, vec: &[f64]) -> Result<(), String> { + unimplemented!() + } + fn swap_rows(&mut self, row1: usize, row2: usize) { + unimplemented!() + } + + fn norm_inf(&self) -> f64 { + unimplemented!() + } + + fn transpose(&self) -> Self { + unimplemented!() + } + fn identity(size: usize) -> Self { + unimplemented!() + } + fn from_diagonal(vec: &[f64]) -> Self { + unimplemented!() + } + fn diagonal(&self) -> Vec { + unimplemented!() + } + fn mldivide(&self, rhs: &mut [f64]) -> Result<(), String> { unimplemented!() } diff --git a/crates/cnvx-math/tests/dense_matrix_tests.rs b/crates/cnvx-math/tests/dense_matrix_tests.rs new file mode 100644 index 0000000..78d2d05 --- /dev/null +++ b/crates/cnvx-math/tests/dense_matrix_tests.rs @@ -0,0 +1,130 @@ +use cnvx_math::{DenseMatrix, Matrix, matrix::vector_norm_l2}; + +const EPSILON: f64 = 1e-9; + +#[test] +fn test_scalar_operations() { + let mut m = DenseMatrix::new(2, 2); + m.set(0, 0, 1.0); + m.set(0, 1, 2.0); + m.set(1, 0, 3.0); + m.set(1, 1, 4.0); + + let add_m = m.add_scalar(2.0); + assert_eq!(add_m.get(0, 0), 3.0); + assert_eq!(add_m.get(1, 1), 6.0); + + let mul_m = m.mul_scalar(3.0); + assert_eq!(mul_m.get(0, 1), 6.0); + assert_eq!(mul_m.get(1, 0), 9.0); +} + +#[test] +fn test_elementwise_multiplication() { + let mut a = DenseMatrix::new(2, 2); + a.set(0, 0, 1.0); + a.set(0, 1, 2.0); + a.set(1, 0, 3.0); + a.set(1, 1, 4.0); + + let mut b = DenseMatrix::new(2, 2); + b.set(0, 0, 2.0); + b.set(0, 1, 3.0); + b.set(1, 0, 4.0); + b.set(1, 1, 5.0); + + let c = a.mul_elementwise(&b).unwrap(); + assert_eq!(c.get(0, 0), 2.0); + assert_eq!(c.get(0, 1), 6.0); + assert_eq!(c.get(1, 0), 12.0); + assert_eq!(c.get(1, 1), 20.0); +} + +#[test] +fn test_row_col_manipulations() { + let mut m = DenseMatrix::new(3, 3); + for i in 0..3 { + for j in 0..3 { + m.set(i, j, (i * 3 + j) as f64); + } + } + // Matrix is: + // [0, 1, 2] + // [3, 4, 5] + // [6, 7, 8] + + assert_eq!(m.get_row(1), vec![3.0, 4.0, 5.0]); + assert_eq!(m.get_col(2), vec![2.0, 5.0, 8.0]); + + m.set_row(0, &[9.0, 9.0, 9.0]).unwrap(); + assert_eq!(m.get(0, 1), 9.0); + + m.set_col(1, &[0.0, 0.0, 0.0]).unwrap(); + assert_eq!(m.get(2, 1), 0.0); + + m.swap_rows(1, 2); + assert_eq!(m.get_row(1), vec![6.0, 0.0, 8.0]); + assert_eq!(m.get_row(2), vec![3.0, 0.0, 5.0]); +} + +#[test] +fn test_norms() { + let mut m = DenseMatrix::new(2, 2); + m.set(0, 0, -3.0); + m.set(0, 1, 4.0); // abs sum = 7 + m.set(1, 0, 1.0); + m.set(1, 1, -8.0); // abs sum = 9 + + assert_eq!(m.norm_inf(), 9.0); + + let v = vec![3.0, -4.0, 0.0]; + assert_eq!(vector_norm_l2(&v), 5.0); +} + +#[test] +fn test_constructors() { + let id = DenseMatrix::identity(3); + assert_eq!(id.get(0, 0), 1.0); + assert_eq!(id.get(1, 1), 1.0); + assert_eq!(id.get(0, 1), 0.0); + + let diag = DenseMatrix::from_diagonal(&[1.5, 2.5]); + assert_eq!(diag.rows(), 2); + assert_eq!(diag.cols(), 2); + assert_eq!(diag.get(0, 0), 1.5); + assert_eq!(diag.get(1, 1), 2.5); + + let mut m = DenseMatrix::new(2, 3); + m.set(0, 0, 7.0); + m.set(1, 1, 8.0); + assert_eq!(m.diagonal(), vec![7.0, 8.0]); +} + +#[test] +fn test_mldivide_success() { + let mut a = DenseMatrix::new(2, 2); + a.set(0, 0, 2.0); + a.set(0, 1, 1.0); + a.set(1, 0, 1.0); + a.set(1, 1, 3.0); + + let mut rhs = vec![3.0, 7.0]; + a.mldivide(&mut rhs).unwrap(); + + assert!((rhs[0] - 0.4).abs() < EPSILON); + assert!((rhs[1] - 2.2).abs() < EPSILON); +} + +#[test] +fn test_mldivide_singular_matrix() { + let mut a = DenseMatrix::new(2, 2); + a.set(0, 0, 1.0); + a.set(0, 1, 1.0); + a.set(1, 0, 1.0); + a.set(1, 1, 1.0); // Singular + + let mut rhs = vec![2.0, 2.0]; + let result = a.mldivide(&mut rhs); + + assert!(result.is_err()); +} From d09ad81b4220cfe6ba6801515b9f42eb5e2ad4e7 Mon Sep 17 00:00:00 2001 From: Chris Date: Thu, 18 Jun 2026 17:39:27 +1200 Subject: [PATCH 2/2] fix: resolve clippy --- crates/cnvx-math/src/matrix/dense.rs | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/crates/cnvx-math/src/matrix/dense.rs b/crates/cnvx-math/src/matrix/dense.rs index e592fd2..b8bf00f 100644 --- a/crates/cnvx-math/src/matrix/dense.rs +++ b/crates/cnvx-math/src/matrix/dense.rs @@ -86,12 +86,12 @@ impl Matrix for DenseMatrix { return Err("Matrix columns must match vector length".to_string()); } let mut result = vec![0.0; self.rows]; - for i in 0..self.rows { + for (i, res_i) in result.iter_mut().enumerate() { let mut sum = 0.0; - for j in 0..self.cols { - sum += self.get(i, j) * rhs[j]; + for (j, &rhs_j) in rhs.iter().enumerate() { + sum += self.get(i, j) * rhs_j; } - result[i] = sum; + *res_i = sum; } Ok(result) } @@ -160,8 +160,8 @@ impl Matrix for DenseMatrix { if vec.len() != self.rows { return Err("Vector length must match matrix rows".to_string()); } - for i in 0..self.rows { - self.data[i * self.cols + col] = vec[i]; + for (i, &val) in vec.iter().enumerate() { + self.data[i * self.cols + col] = val; } Ok(()) } @@ -214,8 +214,8 @@ impl Matrix for DenseMatrix { fn from_diagonal(vec: &[f64]) -> Self { let size = vec.len(); let mut m = Self::new(size, size); - for i in 0..size { - m.set(i, i, vec[i]); + for (i, &val) in vec.iter().enumerate() { + m.set(i, i, val); } m }