From 32e4b0d2957f62798eb341ee4813ab0440380620 Mon Sep 17 00:00:00 2001 From: Chris Date: Thu, 18 Jun 2026 14:43:23 +1200 Subject: [PATCH 01/10] 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 958dd9a8c682f2fc7ce917d45c79421c6e288f7a Mon Sep 17 00:00:00 2001 From: Chris Date: Thu, 18 Jun 2026 16:27:53 +1200 Subject: [PATCH 02/10] feat: add openblas and lapack support for dense matrix solvers --- .github/workflows/ci.yml | 13 +- Cargo.lock | 277 ++++++++++++++++++ Cargo.toml | 4 +- crates/cnvx-lp/src/primal_simplex.rs | 64 ++-- crates/cnvx-math/Cargo.toml | 3 + crates/cnvx-math/src/lib.rs | 2 + crates/cnvx-math/src/matrix/dense.rs | 129 +++----- crates/cnvx-math/src/matrix/mod.rs | 25 +- .../cnvx-math/src/matrix/solvers/cholesky.rs | 30 ++ crates/cnvx-math/src/matrix/solvers/lu.rs | 36 +++ crates/cnvx-math/src/matrix/solvers/mod.rs | 52 ++++ crates/cnvx-math/src/matrix/solvers/qr.rs | 47 +++ .../src/matrix/solvers/triangular.rs | 75 +++++ crates/cnvx-math/src/matrix/sparse.rs | 2 +- crates/cnvx-math/tests/dense_matrix_tests.rs | 30 +- crates/cnvx-math/tests/mldivide_tests.rs | 60 ++++ tests/src/netlib.rs | 2 +- 17 files changed, 689 insertions(+), 162 deletions(-) create mode 100644 crates/cnvx-math/src/matrix/solvers/cholesky.rs create mode 100644 crates/cnvx-math/src/matrix/solvers/lu.rs create mode 100644 crates/cnvx-math/src/matrix/solvers/mod.rs create mode 100644 crates/cnvx-math/src/matrix/solvers/qr.rs create mode 100644 crates/cnvx-math/src/matrix/solvers/triangular.rs create mode 100644 crates/cnvx-math/tests/mldivide_tests.rs diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index ce6d0e4..1244f30 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -50,6 +50,7 @@ jobs: echo "$filters" >> $GITHUB_OUTPUT echo "EOF" >> $GITHUB_OUTPUT + changes: name: Detect Changes runs-on: ubuntu-latest @@ -70,6 +71,7 @@ jobs: select(. != "integration")]') echo "crates=$crates" >> $GITHUB_OUTPUT + required: name: Required Checks runs-on: ubuntu-latest @@ -108,6 +110,7 @@ jobs: uses: actions/checkout@v6 - name: Build run: cargo build --workspace --all-targets + lint: name: 'Lint, Format and Docs - Workspace' runs-on: ubuntu-latest @@ -122,6 +125,7 @@ jobs: run: >- RUSTDOCFLAGS="-D warnings" cargo doc --workspace --no-deps --document-private-items + test: name: Test - Crate runs-on: ubuntu-latest @@ -136,6 +140,10 @@ jobs: - uses: actions/checkout@v6 - name: 'Run ${{ matrix.crate }} tests' run: 'cargo test -p ${{ matrix.crate }}' + + # TODO: Not needed for now, as the integration tests are minor. + # This should be added back once there are more substantial integration tests, and the project + # is more mature. integration-test: name: Integration Test - Crate runs-on: ubuntu-latest @@ -149,5 +157,6 @@ jobs: steps: - name: Checkout uses: actions/checkout@v6 - - name: Run integration tests - run: timeout 240s cargo netlib -- --nocapture + # TODO: No-op + # - name: Run integration tests + # run: timeout 240s cargo netlib -- --nocapture diff --git a/Cargo.lock b/Cargo.lock index c1481bf..ff7c9cc 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -2,6 +2,12 @@ # It is not intended for manual editing. version = 4 +[[package]] +name = "adler2" +version = "2.0.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "320119579fcad9c21884f5c4861d16174d0e06250625266f50fe6898340abefa" + [[package]] name = "anstream" version = "0.6.21" @@ -52,12 +58,24 @@ dependencies = [ "windows-sys 0.61.2", ] +[[package]] +name = "anyhow" +version = "1.0.102" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7f202df86484c868dbad7eaa557ef785d5c66295e41b460ef922eca0723b842c" + [[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.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f2032f911046de80f0a198e0901378627c33f59ea0ac00e363d481118bd70a53" + [[package]] name = "aws-lc-rs" version = "1.15.4" @@ -104,6 +122,26 @@ version = "1.11.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "1e748733b7cbc798e1434b6ac524f0c1ff2ab456fe201501e6497c8417a4fc33" +[[package]] +name = "cblas" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c7c834e5b5c32a6dbcdb917788d0c38473e1ae7076ff2f4011eb96c5fec20eef" +dependencies = [ + "cblas-sys", + "libc", + "num-complex", +] + +[[package]] +name = "cblas-sys" +version = "0.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8354a31ad0d37f735e74c1973343c99dde4461ea8626a2b8d4fea9e06b987255" +dependencies = [ + "libc", +] + [[package]] name = "cc" version = "1.2.55" @@ -227,6 +265,11 @@ dependencies = [ [[package]] name = "cnvx-math" version = "0.0.1" +dependencies = [ + "cblas", + "lapacke", + "openblas-src", +] [[package]] name = "cnvx-parse" @@ -315,6 +358,36 @@ version = "0.8.7" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "773648b94d0e5d620f64f280777445740e61fe701025087ec8b57f45c791888b" +[[package]] +name = "crc32fast" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9481c1c90cbf2ac953f07c8d4a58aa3945c425b7185c9154d67a65e4230da511" +dependencies = [ + "cfg-if", +] + +[[package]] +name = "dirs" +version = "6.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c3e8aa94d75141228480295a7d0e7feb620b1a5ad9f12bc40be62411e38cce4e" +dependencies = [ + "dirs-sys", +] + +[[package]] +name = "dirs-sys" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e01a3366d27ee9890022452ee61b2b63a67e6f13f58900b651ff5665f0bb1fab" +dependencies = [ + "libc", + "option-ext", + "redox_users", + "windows-sys 0.61.2", +] + [[package]] name = "displaydoc" version = "0.2.5" @@ -363,12 +436,32 @@ version = "2.3.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "37909eebbb50d72f9059c3b6d82c0463f2ff062c9e95845c43a6c9c0355411be" +[[package]] +name = "filetime" +version = "0.2.29" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5c287a33c7f0a620c38e641e7f60827713987b3c0f26e8ddc9462cc69cf75759" +dependencies = [ + "cfg-if", + "libc", +] + [[package]] name = "find-msvc-tools" version = "0.1.9" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "5baebc0774151f905a1a2cc41989300b1e6fbb29aff0ceffa1064fdd3088d582" +[[package]] +name = "flate2" +version = "1.1.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "843fba2746e448b37e26a819579957415c8cef339bf08564fe8b7ddbd959573c" +dependencies = [ + "crc32fast", + "miniz_oxide", +] + [[package]] name = "fnv" version = "1.0.7" @@ -782,12 +875,41 @@ dependencies = [ "wasm-bindgen", ] +[[package]] +name = "lapacke" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8425aee3cfc69f1e907f8487a291d1ddebfd0db2d74493b4e443be0618648744" +dependencies = [ + "lapacke-sys", + "libc", + "num-complex", +] + +[[package]] +name = "lapacke-sys" +version = "0.1.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2f7d0817c6f4a6029f3b153de01d6498dcf9df659a7536c58bd8df5cd3ccaa6e" +dependencies = [ + "libc", +] + [[package]] name = "libc" version = "0.2.181" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "459427e2af2b9c839b132acb702a1c654d95e10f8c326bfc2ad11310e458b1c5" +[[package]] +name = "libredox" +version = "0.1.17" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f02ab6bace2054fb888a3c16f990117b579d14a3088e472d63c6011fa185c9d3" +dependencies = [ + "libc", +] + [[package]] name = "linux-raw-sys" version = "0.11.0" @@ -824,6 +946,16 @@ version = "0.3.17" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "6877bb514081ee2a7ff5ef9de3281f14a4dd4bceac4c09388074a6b5df8a139a" +[[package]] +name = "miniz_oxide" +version = "0.8.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1fa76a2c86f704bdb222d66965fb3d63269ce38518b83cb0575fca855ebb6316" +dependencies = [ + "adler2", + "simd-adler32", +] + [[package]] name = "mio" version = "1.1.1" @@ -835,6 +967,24 @@ dependencies = [ "windows-sys 0.61.2", ] +[[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-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" @@ -847,12 +997,44 @@ version = "1.70.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "384b8ab6d37215f3c5301a95a4accb5d64aa607f1fcb26a11b5303878451b4fe" +[[package]] +name = "openblas-build" +version = "0.10.16" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bb9c85e9e7dd5acdc67b9f3f0c99656b550df716bc63540c6a224a920754a5c2" +dependencies = [ + "anyhow", + "cc", + "flate2", + "tar", + "thiserror 2.0.18", + "ureq", +] + +[[package]] +name = "openblas-src" +version = "0.10.16" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f1a81a5e467f1861ad6ac32d5ec1690ad097d19854753b7424250fe27da46b98" +dependencies = [ + "dirs", + "openblas-build", + "pkg-config", + "vcpkg", +] + [[package]] name = "openssl-probe" version = "0.2.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "7c87def4c32ab89d880effc9e097653c8da5d6ef28e6b539d313baaacfbafcbe" +[[package]] +name = "option-ext" +version = "0.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "04744f49eae99ab78e0d5c0b603ab218f515ea8cfe5a456d7629ad883a3b6e7d" + [[package]] name = "percent-encoding" version = "2.3.2" @@ -871,6 +1053,12 @@ version = "0.1.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "8b870d8c151b6f2fb93e84a13146138f05d02ed11c7e7c54f8826aaaf7c9f184" +[[package]] +name = "pkg-config" +version = "0.3.33" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "19f132c84eca552bf34cab8ec81f1c1dcc229b811638f9d283dceabe58c5569e" + [[package]] name = "portable-atomic" version = "1.13.1" @@ -1062,6 +1250,17 @@ dependencies = [ "getrandom 0.3.4", ] +[[package]] +name = "redox_users" +version = "0.5.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a4e608c6638b9c18977b00b475ac1f28d14e84b27d8d42f70e0bf1e3dec127ac" +dependencies = [ + "getrandom 0.2.17", + "libredox", + "thiserror 2.0.18", +] + [[package]] name = "reqwest" version = "0.13.2" @@ -1142,7 +1341,9 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "c665f33d38cea657d9614f766881e4d510e0eda4239891eea56b4cadcf01801b" dependencies = [ "aws-lc-rs", + "log", "once_cell", + "ring", "rustls-pki-types", "rustls-webpki", "subtle", @@ -1315,6 +1516,12 @@ dependencies = [ "libc", ] +[[package]] +name = "simd-adler32" +version = "0.3.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "703d5c7ef118737c72f1af64ad2f6f8c5e1921f818cdcb97b8fe6fc69bf66214" + [[package]] name = "slab" version = "0.4.12" @@ -1407,6 +1614,17 @@ dependencies = [ "libc", ] +[[package]] +name = "tar" +version = "0.4.46" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3f6221d9a6003c78398e3b239969f352578258df48c8eb051caadae0015bc840" +dependencies = [ + "filetime", + "libc", + "xattr", +] + [[package]] name = "target-lexicon" version = "0.13.5" @@ -1643,6 +1861,34 @@ version = "0.9.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "8ecb6da28b8a351d773b68d5825ac39017e680750f980f3a1a85cd8dd28a47c1" +[[package]] +name = "ureq" +version = "3.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dea7109cdcd5864d4eeb1b58a1648dc9bf520360d7af16ec26d0a9354bafcfc0" +dependencies = [ + "base64", + "log", + "percent-encoding", + "rustls", + "rustls-pki-types", + "ureq-proto", + "utf8-zero", + "webpki-roots", +] + +[[package]] +name = "ureq-proto" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e994ba84b0bd1b1b0cf92878b7ef898a5c1760108fe7b6010327e274917a808c" +dependencies = [ + "base64", + "http", + "httparse", + "log", +] + [[package]] name = "url" version = "2.5.8" @@ -1655,6 +1901,12 @@ dependencies = [ "serde", ] +[[package]] +name = "utf8-zero" +version = "0.8.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b8c0a043c9540bae7c578c88f91dda8bd82e59ae27c21baca69c8b191aaf5a6e" + [[package]] name = "utf8_iter" version = "1.0.4" @@ -1667,6 +1919,12 @@ version = "0.2.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "06abde3611657adf66d383f00b093d7faecc7fa57071cce2578660c9f1010821" +[[package]] +name = "vcpkg" +version = "0.2.15" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "accd4ea62f7bb7a82fe23066fb0957d48ef677f6eeb8215f372f52e48bb32426" + [[package]] name = "venial" version = "0.5.0" @@ -1810,6 +2068,15 @@ dependencies = [ "rustls-pki-types", ] +[[package]] +name = "webpki-roots" +version = "1.0.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bf85cb06032201fa7c6f829d7db5a7e5aa45bcc0655327713065f6f0576731bf" +dependencies = [ + "rustls-pki-types", +] + [[package]] name = "winapi-util" version = "0.1.11" @@ -2088,6 +2355,16 @@ version = "0.6.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "9edde0db4769d2dc68579893f2306b26c6ecfbe0ef499b013d731b7b9247e0b9" +[[package]] +name = "xattr" +version = "1.6.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "32e45ad4206f6d2479085147f02bc2ef834ac85886624a23575ae137c8aa8156" +dependencies = [ + "libc", + "rustix", +] + [[package]] name = "yoke" version = "0.8.1" diff --git a/Cargo.toml b/Cargo.toml index de768ee..9ee740e 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -47,7 +47,9 @@ cnvx-math = { path = "crates/cnvx-math", version = "0.0.1" } cnvx-graph = { path = "crates/cnvx-graph", version = "0.0.1" } cnvx-parse = { path = "crates/cnvx-parse", version = "0.0.1" } clap = { version = "4.5.57", features = ["derive"] } -nalgebra = "0.34.2" +cblas = { version = "0.5.0" } +lapacke = { version = "0.5.0" } +openblas-src = { version = "0.10.16" } [workspace.lints] diff --git a/crates/cnvx-lp/src/primal_simplex.rs b/crates/cnvx-lp/src/primal_simplex.rs index 8901962..28f9a37 100644 --- a/crates/cnvx-lp/src/primal_simplex.rs +++ b/crates/cnvx-lp/src/primal_simplex.rs @@ -264,7 +264,7 @@ impl PrimalSimplexState { /// Attempt to directly run phase 2 if the initial basis is feasible. fn try_phase2(&mut self, max_iter: usize, tol: f64) -> Result { let mut bmat = self.build_bmat(); - match self.compute_basic_solution(&mut bmat) { + match self.compute_basic_solution(&bmat) { Ok(xb) if xb.iter().all(|&v| v >= -tol) => { self.x_b = xb; self.remove_artificial_from_basis(&mut bmat, self.a.cols()) @@ -366,19 +366,16 @@ impl PrimalSimplexState { pub fn build_bmat(&self) -> A { let m = self.a.rows(); let mut bmat = A::new(m, m); - for i in 0..m { - for j in 0..m { - bmat.set(i, j, self.a.get(i, self.basis[j])); - } + for j in 0..m { + let col = self.a.get_col(self.basis[j]); + bmat.set_col(j, &col).unwrap(); } bmat } /// Compute the values of the basic variables by solving `B x_B = b`. - pub fn compute_basic_solution(&self, bmat: &mut A) -> Result, String> { - let mut xb = self.b.clone(); - bmat.mldivide(&mut xb).map_err(|e| format!("gauss failed: {e}"))?; - Ok(xb) + pub fn compute_basic_solution(&self, bmat: &A) -> Result, String> { + bmat.mldivide(&self.b).map_err(|e| format!("gauss failed: {e}")) } /// Run the main simplex iteration loop. @@ -423,19 +420,11 @@ impl PrimalSimplexState { /// Compute dual variables for the current basis. fn compute_duals(&self, bmat: &A) -> Result, SolveError> { let m = bmat.rows(); - let mut pi = (0..m).map(|i| self.c[self.basis[i]]).collect::>(); - - let mut bt = A::new(m, m); - for i in 0..m { - for j in 0..m { - bt.set(i, j, bmat.get(j, i)); - } - } - - bt.mldivide(&mut pi) - .map_err(|e| SolveError::Other(format!("dual solve failed: {e}")))?; + let pi_input: Vec = (0..m).map(|i| self.c[self.basis[i]]).collect(); - Ok(pi) + let bt = bmat.transpose(); + bt.mldivide(&pi_input) + .map_err(|e| SolveError::Other(format!("dual solve failed: {e}"))) } /// Choose entering variable using reduced costs. @@ -455,15 +444,12 @@ impl PrimalSimplexState { /// Compute the simplex direction `d = B^{-1} A_j`. fn compute_direction( &self, - bmat: &mut A, + bmat: &A, entering: usize, ) -> Result, SolveError> { - let mut d = (0..bmat.rows()).map(|i| self.a.get(i, entering)).collect::>(); - - bmat.mldivide(&mut d) - .map_err(|e| SolveError::Other(format!("direction solve failed: {e}")))?; - - Ok(d) + let d_input = self.a.get_col(entering); + bmat.mldivide(&d_input) + .map_err(|e| SolveError::Other(format!("direction solve failed: {e}"))) } /// Choose leaving variable using minimum ratio test. @@ -497,9 +483,8 @@ impl PrimalSimplexState { self.basis[leave_row] = entering; self.non_basis[enter_pos] = leaving; - for i in 0..bmat.rows() { - bmat.set(i, leave_row, self.a.get(i, entering)); - } + let col = self.a.get_col(entering); + bmat.set_col(leave_row, &col).unwrap(); } /// Update the current objective value. @@ -552,10 +537,9 @@ impl PrimalSimplexState { self.x_b = b_aug; let mut bmat = A::new(m, m); - for i in 0..m { - for j in 0..m { - bmat.set(i, j, self.a.get(i, self.basis[j])); - } + for j in 0..m { + let col = self.a.get_col(self.basis[j]); + bmat.set_col(j, &col).unwrap(); } (orig_a, orig_c, bmat) @@ -582,9 +566,8 @@ impl PrimalSimplexState { let leaving = self.basis[row]; self.basis[row] = j; self.non_basis[nb_pos] = leaving; - for i in 0..m { - bmat.set(i, row, self.a.get(i, j)); - } + let col = self.a.get_col(j); + bmat.set_col(row, &col).unwrap(); } else if self.x_b[row].abs() > 1e-12 { return Err( "artificial variable left in basis with non-zero value".into() @@ -595,9 +578,8 @@ impl PrimalSimplexState { let leaving = self.basis[row]; self.basis[row] = j; self.non_basis[nb_pos] = leaving; - for i in 0..m { - bmat.set(i, row, self.a.get(i, j)); - } + let col = self.a.get_col(j); + bmat.set_col(row, &col).unwrap(); break; } } diff --git a/crates/cnvx-math/Cargo.toml b/crates/cnvx-math/Cargo.toml index fc7586d..14a4b10 100644 --- a/crates/cnvx-math/Cargo.toml +++ b/crates/cnvx-math/Cargo.toml @@ -10,6 +10,9 @@ keywords = { workspace = true } readme = { workspace = true } [dependencies] +cblas = { workspace = true } +lapacke = { workspace = true } +openblas-src = { workspace = true } # Links with openblas, which provides BLAS and LAPACK implementations [lints] workspace = true diff --git a/crates/cnvx-math/src/lib.rs b/crates/cnvx-math/src/lib.rs index f065436..02af3bf 100644 --- a/crates/cnvx-math/src/lib.rs +++ b/crates/cnvx-math/src/lib.rs @@ -8,6 +8,8 @@ //! //! - [`matrix`]: Defines [`DenseMatrix`] and the [`Matrix`] trait for linear algebra operations. +extern crate openblas_src; + pub mod matrix; 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 e592fd2..24763ff 100644 --- a/crates/cnvx-math/src/matrix/dense.rs +++ b/crates/cnvx-math/src/matrix/dense.rs @@ -1,4 +1,5 @@ -use crate::matrix::Matrix; +use crate::matrix::{Matrix, solvers}; +use cblas::{Layout, Transpose, dgemm, dgemv}; /// A dense matrix stored in row-major order using a flat 1D vector. #[derive(Debug, Clone, PartialEq)] @@ -14,6 +15,11 @@ impl DenseMatrix { fn index(&self, row: usize, col: usize) -> usize { row * self.cols + col } + + /// Expose the underlying data slice to internal LAPACK solvers + pub(crate) fn data(&self) -> &[f64] { + &self.data + } } impl Matrix for DenseMatrix { @@ -65,18 +71,28 @@ impl Matrix for DenseMatrix { return Err("Incompatible dimensions for matrix multiplication".to_string()); } + let m = self.rows as i32; + let k = self.cols as i32; + let n = other.cols as i32; 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); - } - } + unsafe { + dgemm( + Layout::RowMajor, + Transpose::None, + Transpose::None, + m, + n, + k, + 1.0, + &self.data, + k, // lda + &other.data, + n, // ldb + 0.0, + &mut result.data, + n, // ldc + ); } Ok(result) } @@ -85,13 +101,26 @@ impl Matrix for DenseMatrix { if self.cols != rhs.len() { return Err("Matrix columns must match vector length".to_string()); } + + let m = self.rows as i32; + let n = self.cols as i32; 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; + + unsafe { + dgemv( + Layout::RowMajor, + Transpose::None, + m, + n, + 1.0, + &self.data, + n, // lda + rhs, + 1, // incx + 0.0, + &mut result, + 1, // incy + ); } Ok(result) } @@ -171,7 +200,6 @@ impl Matrix for DenseMatrix { 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); @@ -225,68 +253,7 @@ impl Matrix for DenseMatrix { (0..size).map(|i| self.get(i, i)).collect() } - fn mldivide(&self, rhs: &mut [f64]) -> Result<(), String> { - 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); - } - 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]; - } - } - - // 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(()) + fn mldivide(&self, rhs: &[f64]) -> Result, String> { + solvers::mldivide_dense(self, rhs) } } diff --git a/crates/cnvx-math/src/matrix/mod.rs b/crates/cnvx-math/src/matrix/mod.rs index 03874be..daa0e2f 100644 --- a/crates/cnvx-math/src/matrix/mod.rs +++ b/crates/cnvx-math/src/matrix/mod.rs @@ -1,4 +1,5 @@ mod dense; +mod solvers; mod sparse; pub use dense::DenseMatrix; @@ -9,8 +10,6 @@ pub use sparse::SparseMatrix; /// This trait defines the interface for matrices used in the LP solver, /// 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. fn new(rows: usize, cols: usize) -> Self @@ -35,8 +34,6 @@ pub trait Matrix: 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 @@ -55,8 +52,6 @@ pub trait Matrix: Clone { /// 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 @@ -82,8 +77,6 @@ pub trait Matrix: Clone { where Self: Sized; - // --- 4. Row and Column Manipulations --- - /// Extracts a specific row as a standard vector. fn get_row(&self, row: usize) -> Vec; @@ -99,13 +92,9 @@ pub trait Matrix: Clone { /// 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 @@ -124,14 +113,10 @@ pub trait Matrix: Clone { /// 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`. - fn mldivide(&self, rhs: &mut [f64]) -> Result<(), String> - where - Self: Sized; + /// Solves Ax = b. + /// Returns a dynamically allocated Vec to support both square + /// and rectangular (least-squares) solutions. + fn mldivide(&self, rhs: &[f64]) -> Result, String>; } /// Helper function to calculate the L2 (Euclidean) norm of a standard vector. diff --git a/crates/cnvx-math/src/matrix/solvers/cholesky.rs b/crates/cnvx-math/src/matrix/solvers/cholesky.rs new file mode 100644 index 0000000..33989fa --- /dev/null +++ b/crates/cnvx-math/src/matrix/solvers/cholesky.rs @@ -0,0 +1,30 @@ +use crate::matrix::{DenseMatrix, Matrix}; +use lapacke::{Layout, dposv}; + +pub fn solve(a: &DenseMatrix, b: &[f64]) -> Result, String> { + let n = a.rows() as i32; + + let mut a_copy = a.data().to_vec(); + let mut x = b.to_vec(); + + let info = unsafe { + dposv( + Layout::RowMajor, + b'U', // Check Upper triangle + n, + 1, // nrhs + &mut a_copy, + n, // lda + &mut x, + 1, // ldb + ) + }; + + if info > 0 { + return Err("Matrix is not strictly positive definite".into()); + } else if info < 0 { + return Err(format!("Illegal value in argument {}", -info)); + } + + Ok(x) +} diff --git a/crates/cnvx-math/src/matrix/solvers/lu.rs b/crates/cnvx-math/src/matrix/solvers/lu.rs new file mode 100644 index 0000000..40fe348 --- /dev/null +++ b/crates/cnvx-math/src/matrix/solvers/lu.rs @@ -0,0 +1,36 @@ +use crate::matrix::{DenseMatrix, Matrix}; +use lapacke::{Layout, dgesv}; + +pub fn solve(a: &DenseMatrix, b: &[f64]) -> Result, String> { + let n = a.rows() as i32; + + // LAPACK mutates A to store the LU factors, so we must clone it. + let mut a_copy = a.data().to_vec(); + + // LAPACK mutates b to store the solution x. + let mut x = b.to_vec(); + + // Pivot indices array + let mut ipiv = vec![0; n as usize]; + + let info = unsafe { + dgesv( + Layout::RowMajor, + n, + 1, // nrhs + &mut a_copy, + n, // lda + &mut ipiv, + &mut x, + 1, // ldb + ) + }; + + if info > 0 { + return Err("Matrix is singular (U(i,i) is exactly zero)".into()); + } else if info < 0 { + return Err(format!("Illegal value in argument {}", -info)); + } + + Ok(x) +} diff --git a/crates/cnvx-math/src/matrix/solvers/mod.rs b/crates/cnvx-math/src/matrix/solvers/mod.rs new file mode 100644 index 0000000..d5567e8 --- /dev/null +++ b/crates/cnvx-math/src/matrix/solvers/mod.rs @@ -0,0 +1,52 @@ +pub mod cholesky; +pub mod lu; +pub mod qr; +pub mod triangular; + +use crate::matrix::{DenseMatrix, Matrix}; + +/// MATLAB-style mldivide dispatcher for full matrices +/// See https://mathworks.com/help/matlab/ref/double.mldivide.html +pub fn mldivide_dense(a: &DenseMatrix, b: &[f64]) -> Result, String> { + let m = a.rows(); + let n = a.cols(); + + // 1. Is matrix square? If not, use QR for least squares. + if m != n { + return qr::solve(a, b); + } + + if b.len() != n { + return Err("RHS vector length must match matrix rows".to_string()); + } + + // 2. Is it perfectly triangular? + if let Some(tri_type) = triangular::check_structure(a) { + return triangular::solve(a, b, tri_type); + } + + // 3. Is it symmetric? + if is_symmetric(a) { + // Try Cholesky decomposition (Symmetric Positive Definite) + // If it succeeds, we are done. If it fails, fall through to LU. + if let Ok(x) = cholesky::solve(a, b) { + return Ok(x); + } + } + + // 4. Fallback: LU Decomposition with partial pivoting + lu::solve(a, b) +} + +/// Checks if A is symmetric: A_ij == A_ji +fn is_symmetric(a: &DenseMatrix) -> bool { + let n = a.rows(); + for i in 0..n { + for j in (i + 1)..n { + if (a.get(i, j) - a.get(j, i)).abs() > 1e-12 { + return false; + } + } + } + true +} diff --git a/crates/cnvx-math/src/matrix/solvers/qr.rs b/crates/cnvx-math/src/matrix/solvers/qr.rs new file mode 100644 index 0000000..59c30e1 --- /dev/null +++ b/crates/cnvx-math/src/matrix/solvers/qr.rs @@ -0,0 +1,47 @@ +use crate::matrix::{DenseMatrix, Matrix}; +use lapacke::{Layout, dgels}; + +/// Solves overdetermined rectangular systems (Least Squares) using QR factorization. +pub fn solve(a: &DenseMatrix, b: &[f64]) -> Result, String> { + let m = a.rows() as i32; + let n = a.cols() as i32; + + if m < n { + return Err("Underdetermined systems (m < n) are not yet supported.".into()); + } + if b.len() != m as usize { + return Err("RHS vector length must match matrix rows".into()); + } + + let mut a_copy = a.data().to_vec(); + + // LAPACK's dgels requires the RHS array to have a size of max(m, n) + // because it overwrites the RHS with the solution. + let max_mn = std::cmp::max(m, n) as usize; + let mut x = vec![0.0; max_mn]; + x[..b.len()].copy_from_slice(b); + + let info = unsafe { + dgels( + Layout::RowMajor, + b'N', // No transpose + m, + n, + 1, // nrhs + &mut a_copy, + n, // lda + &mut x, + 1, // ldb + ) + }; + + if info > 0 { + return Err("Matrix is rank deficient".into()); + } else if info < 0 { + return Err(format!("Illegal value in argument {}", -info)); + } + + // The first `n` elements contain the least squares solution + x.truncate(n as usize); + Ok(x) +} diff --git a/crates/cnvx-math/src/matrix/solvers/triangular.rs b/crates/cnvx-math/src/matrix/solvers/triangular.rs new file mode 100644 index 0000000..8dd5192 --- /dev/null +++ b/crates/cnvx-math/src/matrix/solvers/triangular.rs @@ -0,0 +1,75 @@ +use crate::matrix::{DenseMatrix, Matrix}; +use cblas::{Diagonal, Layout, Part, Transpose, dtrsv}; + +pub enum TriType { + Upper, + Lower, + Diagonal, +} + +pub fn check_structure(a: &DenseMatrix) -> Option { + let n = a.rows(); + let mut is_upper = true; + let mut is_lower = true; + + for i in 0..n { + for j in 0..n { + let val = a.get(i, j).abs(); + if i > j && val > 1e-12 { + is_upper = false; + } + if i < j && val > 1e-12 { + is_lower = false; + } + } + } + + if is_upper && is_lower { + Some(TriType::Diagonal) + } else if is_upper { + Some(TriType::Upper) + } else if is_lower { + Some(TriType::Lower) + } else { + None + } +} + +pub fn solve(a: &DenseMatrix, b: &[f64], tri_type: TriType) -> Result, String> { + let n = a.rows() as i32; + let mut x = b.to_vec(); + + // Check for exact zeros on the diagonal to avoid NaN pollution before CBLAS + for i in 0..a.rows() { + if a.get(i, i).abs() < 1e-12 { + return Err("Singular matrix".into()); + } + } + + let uplo = match tri_type { + TriType::Upper => Part::Upper, + TriType::Lower => Part::Lower, + TriType::Diagonal => { + // Pure diagonal division is faster than invoking BLAS + for i in 0..a.rows() { + x[i] /= a.get(i, i); + } + return Ok(x); + } + }; + + unsafe { + dtrsv( + Layout::RowMajor, + uplo, + Transpose::None, + Diagonal::Generic, + n, + a.data(), + n, // lda + &mut x, + 1, // incx + ); + } + Ok(x) +} diff --git a/crates/cnvx-math/src/matrix/sparse.rs b/crates/cnvx-math/src/matrix/sparse.rs index b5d9996..04764f1 100644 --- a/crates/cnvx-math/src/matrix/sparse.rs +++ b/crates/cnvx-math/src/matrix/sparse.rs @@ -83,7 +83,7 @@ impl Matrix for SparseMatrix { unimplemented!() } - fn mldivide(&self, rhs: &mut [f64]) -> Result<(), String> { + fn mldivide(&self, rhs: &[f64]) -> Result, String> { unimplemented!() } } diff --git a/crates/cnvx-math/tests/dense_matrix_tests.rs b/crates/cnvx-math/tests/dense_matrix_tests.rs index 78d2d05..d4654a9 100644 --- a/crates/cnvx-math/tests/dense_matrix_tests.rs +++ b/crates/cnvx-math/tests/dense_matrix_tests.rs @@ -1,6 +1,6 @@ use cnvx_math::{DenseMatrix, Matrix, matrix::vector_norm_l2}; -const EPSILON: f64 = 1e-9; +// const EPSILON: f64 = 1e-9; #[test] fn test_scalar_operations() { @@ -100,20 +100,20 @@ fn test_constructors() { 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_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() { diff --git a/crates/cnvx-math/tests/mldivide_tests.rs b/crates/cnvx-math/tests/mldivide_tests.rs new file mode 100644 index 0000000..ee5032a --- /dev/null +++ b/crates/cnvx-math/tests/mldivide_tests.rs @@ -0,0 +1,60 @@ +use cnvx_math::{DenseMatrix, Matrix}; + +#[test] +fn test_mldivide_triangular() { + // Upper triangular + let mut a = DenseMatrix::new(2, 2); + a.set(0, 0, 2.0); + a.set(0, 1, 1.0); + a.set(1, 0, 0.0); + a.set(1, 1, 3.0); + let b = vec![5.0, 6.0]; + let x = a.mldivide(&b).unwrap(); + assert_eq!(x, vec![1.5, 2.0]); +} + +#[test] +fn test_mldivide_cholesky() { + // Symmetric Positive Definite + let mut a = DenseMatrix::new(2, 2); + a.set(0, 0, 4.0); + a.set(0, 1, 1.0); + a.set(1, 0, 1.0); + a.set(1, 1, 3.0); + let b = vec![1.0, 2.0]; + let x = a.mldivide(&b).unwrap(); + // Validates via Cholesky decomposition automatically + assert!((x[0] - 0.09090909).abs() < 1e-6); + assert!((x[1] - 0.63636363).abs() < 1e-6); +} + +#[test] +fn test_mldivide_lu() { + // General Square (Requires LU Pivot) + let mut a = DenseMatrix::new(2, 2); + a.set(0, 0, 0.0); + a.set(0, 1, 1.0); + a.set(1, 0, 1.0); + a.set(1, 1, 0.0); + let b = vec![2.0, 3.0]; + let x = a.mldivide(&b).unwrap(); + assert_eq!(x, vec![3.0, 2.0]); +} + +#[test] +fn test_mldivide_qr_least_squares() { + // 3x2 Overdetermined system (Line fitting y = mx + c) + let mut a = DenseMatrix::new(3, 2); + a.set(0, 0, 1.0); + a.set(0, 1, 1.0); + a.set(1, 0, 2.0); + a.set(1, 1, 1.0); + a.set(2, 0, 3.0); + a.set(2, 1, 1.0); + let b = vec![1.2, 1.9, 3.2]; + let x = a.mldivide(&b).unwrap(); + + // Result x[0] == slope, x[1] == y-intercept + assert!((x[0] - 1.0).abs() < 1e-5); + assert!((x[1] - 0.1).abs() < 1e-5); +} diff --git a/tests/src/netlib.rs b/tests/src/netlib.rs index feb108e..a489f0d 100644 --- a/tests/src/netlib.rs +++ b/tests/src/netlib.rs @@ -144,7 +144,7 @@ fn run_cnvx(mps: &Path) -> Result { } #[test_case("afiro", Some(-4.6475314286E+02))] -// #[test_case("adlittle", Some(2.2549496316E+05))] +#[test_case("adlittle", Some(2.2549496316E+05))] #[test_case("sc50a", Some(-6.4575077059E+01))] #[test_case("sc50b", Some(-7.0000000000E+01))] // #[test_case("sc105", Some(-5.2202061212E+01))] From a38869c31e1a21802509106042bd55ce4db7a545 Mon Sep 17 00:00:00 2001 From: Chris Date: Thu, 18 Jun 2026 17:28:42 +1200 Subject: [PATCH 03/10] chore: use system openblas --- .github/workflows/ci.yml | 2 ++ Cargo.toml | 6 +++++- crates/cnvx-math/tests/dense_matrix_tests.rs | 2 +- 3 files changed, 8 insertions(+), 2 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 1244f30..a48bdc3 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -108,6 +108,8 @@ jobs: steps: - name: Checkout uses: actions/checkout@v6 + - name: Install OpenBLAS + run: sudo apt-get update && sudo apt-get install -y libopenblas-dev liblapacke-dev - name: Build run: cargo build --workspace --all-targets diff --git a/Cargo.toml b/Cargo.toml index 9ee740e..83ab4d0 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -49,7 +49,11 @@ cnvx-parse = { path = "crates/cnvx-parse", version = "0.0.1" } clap = { version = "4.5.57", features = ["derive"] } cblas = { version = "0.5.0" } lapacke = { version = "0.5.0" } -openblas-src = { version = "0.10.16" } +openblas-src = { version = "0.10.16", features = [ + "cblas", + "lapacke", + "system", +] } [workspace.lints] diff --git a/crates/cnvx-math/tests/dense_matrix_tests.rs b/crates/cnvx-math/tests/dense_matrix_tests.rs index d4654a9..96fa77a 100644 --- a/crates/cnvx-math/tests/dense_matrix_tests.rs +++ b/crates/cnvx-math/tests/dense_matrix_tests.rs @@ -124,7 +124,7 @@ fn test_mldivide_singular_matrix() { a.set(1, 1, 1.0); // Singular let mut rhs = vec![2.0, 2.0]; - let result = a.mldivide(&mut rhs); + let result = a.mldivide(&rhs); assert!(result.is_err()); } From f6b42bc7d5d10ac5cbed78ed46ca6f8efdad9749 Mon Sep 17 00:00:00 2001 From: Chris Date: Thu, 18 Jun 2026 17:30:35 +1200 Subject: [PATCH 04/10] ci: add openblas to ci tests --- .github/workflows/ci.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index a48bdc3..07e0bc2 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -140,6 +140,8 @@ jobs: crate: '${{ fromJSON(needs.changes.outputs.crates) }}' steps: - uses: actions/checkout@v6 + - name: Install OpenBLAS + run: sudo apt-get update && sudo apt-get install -y libopenblas-dev liblapacke-dev - name: 'Run ${{ matrix.crate }} tests' run: 'cargo test -p ${{ matrix.crate }}' From df8531469461d50eda7f5c0855f776141b2da619 Mon Sep 17 00:00:00 2001 From: Chris Date: Thu, 18 Jun 2026 17:48:13 +1200 Subject: [PATCH 05/10] ci: add explicit library path --- .github/workflows/ci.yml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 07e0bc2..1718e0a 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -111,6 +111,8 @@ jobs: - name: Install OpenBLAS run: sudo apt-get update && sudo apt-get install -y libopenblas-dev liblapacke-dev - name: Build + env: + LIBRARY_PATH: /usr/lib/x86_64-linux-gnu/openblas-pthread:/usr/lib/x86_64-linux-gnu/lapack run: cargo build --workspace --all-targets lint: @@ -143,6 +145,8 @@ jobs: - name: Install OpenBLAS run: sudo apt-get update && sudo apt-get install -y libopenblas-dev liblapacke-dev - name: 'Run ${{ matrix.crate }} tests' + env: + LIBRARY_PATH: /usr/lib/x86_64-linux-gnu/openblas-pthread:/usr/lib/x86_64-linux-gnu/lapack run: 'cargo test -p ${{ matrix.crate }}' # TODO: Not needed for now, as the integration tests are minor. From 6ff93219f2c306e1b12cf747f18be9048d111320 Mon Sep 17 00:00:00 2001 From: Chris Date: Thu, 18 Jun 2026 17:49:59 +1200 Subject: [PATCH 06/10] fix: clippy --- crates/cnvx-math/src/matrix/dense.rs | 8 ++++---- crates/cnvx-math/src/matrix/solvers/triangular.rs | 4 ++-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/crates/cnvx-math/src/matrix/dense.rs b/crates/cnvx-math/src/matrix/dense.rs index 24763ff..61a8620 100644 --- a/crates/cnvx-math/src/matrix/dense.rs +++ b/crates/cnvx-math/src/matrix/dense.rs @@ -189,8 +189,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(()) } @@ -242,8 +242,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.data[i * size + i] = val; // More efficient than set() for diagonal } m } diff --git a/crates/cnvx-math/src/matrix/solvers/triangular.rs b/crates/cnvx-math/src/matrix/solvers/triangular.rs index 8dd5192..9e58dc7 100644 --- a/crates/cnvx-math/src/matrix/solvers/triangular.rs +++ b/crates/cnvx-math/src/matrix/solvers/triangular.rs @@ -51,8 +51,8 @@ pub fn solve(a: &DenseMatrix, b: &[f64], tri_type: TriType) -> Result, TriType::Lower => Part::Lower, TriType::Diagonal => { // Pure diagonal division is faster than invoking BLAS - for i in 0..a.rows() { - x[i] /= a.get(i, i); + for (i, xi) in x.iter_mut().enumerate() { + *xi /= a.get(i, i); } return Ok(x); } From d161bf770233c2be966652cdab5dbb3ebce0c3cf Mon Sep 17 00:00:00 2001 From: Chris Date: Thu, 18 Jun 2026 18:03:24 +1200 Subject: [PATCH 07/10] chore: update build system for cnvx-math --- crates/cnvx-math/build.rs | 11 +++++++++++ 1 file changed, 11 insertions(+) create mode 100644 crates/cnvx-math/build.rs diff --git a/crates/cnvx-math/build.rs b/crates/cnvx-math/build.rs new file mode 100644 index 0000000..675b380 --- /dev/null +++ b/crates/cnvx-math/build.rs @@ -0,0 +1,11 @@ +fn main() { + // On Arch (and some other distros), LAPACKE is bundled into liblapack.so + // On Ubuntu/Debian, it's a separate liblapacke.so + if std::path::Path::new("/usr/lib/liblapacke.so").exists() + || std::path::Path::new("/usr/lib/x86_64-linux-gnu/liblapacke.so").exists() + { + println!("cargo:rustc-link-lib=lapacke"); + } else { + println!("cargo:rustc-link-lib=lapack"); + } +} From 05ea167358af8079ca50d15a0d105fa34249e487 Mon Sep 17 00:00:00 2001 From: Chris Date: Thu, 18 Jun 2026 18:12:20 +1200 Subject: [PATCH 08/10] chore: fix clippy, use pkg-config --- Cargo.lock | 1 + crates/cnvx-math/Cargo.toml | 3 +++ crates/cnvx-math/build.rs | 20 ++++++++++++-------- crates/cnvx-math/tests/dense_matrix_tests.rs | 2 +- 4 files changed, 17 insertions(+), 9 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index ff7c9cc..d0a712b 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -269,6 +269,7 @@ dependencies = [ "cblas", "lapacke", "openblas-src", + "pkg-config", ] [[package]] diff --git a/crates/cnvx-math/Cargo.toml b/crates/cnvx-math/Cargo.toml index 14a4b10..f4f758b 100644 --- a/crates/cnvx-math/Cargo.toml +++ b/crates/cnvx-math/Cargo.toml @@ -14,5 +14,8 @@ cblas = { workspace = true } lapacke = { workspace = true } openblas-src = { workspace = true } # Links with openblas, which provides BLAS and LAPACK implementations +[build-dependencies] +pkg-config = { version = "0.3.33" } + [lints] workspace = true diff --git a/crates/cnvx-math/build.rs b/crates/cnvx-math/build.rs index 675b380..e6d2002 100644 --- a/crates/cnvx-math/build.rs +++ b/crates/cnvx-math/build.rs @@ -1,11 +1,15 @@ fn main() { - // On Arch (and some other distros), LAPACKE is bundled into liblapack.so - // On Ubuntu/Debian, it's a separate liblapacke.so - if std::path::Path::new("/usr/lib/liblapacke.so").exists() - || std::path::Path::new("/usr/lib/x86_64-linux-gnu/liblapacke.so").exists() - { - println!("cargo:rustc-link-lib=lapacke"); - } else { - println!("cargo:rustc-link-lib=lapack"); + if cfg!(target_os = "windows") { + println!("cargo:error=Windows is not currently supported."); + std::process::exit(1); } + + // Try lapacke first, fall back to lapack. + if pkg_config::probe_library("lapacke").is_err() { + pkg_config::probe_library("lapack") + .expect("could not find lapack or lapacke via pkg-config"); + } + + // Add blas here as a "just in case" + pkg_config::probe_library("blas").expect("could not find blas via pkg-config"); } diff --git a/crates/cnvx-math/tests/dense_matrix_tests.rs b/crates/cnvx-math/tests/dense_matrix_tests.rs index 96fa77a..3b3bd4d 100644 --- a/crates/cnvx-math/tests/dense_matrix_tests.rs +++ b/crates/cnvx-math/tests/dense_matrix_tests.rs @@ -123,7 +123,7 @@ fn test_mldivide_singular_matrix() { a.set(1, 0, 1.0); a.set(1, 1, 1.0); // Singular - let mut rhs = vec![2.0, 2.0]; + let rhs = vec![2.0, 2.0]; let result = a.mldivide(&rhs); assert!(result.is_err()); From 8d1837448b86d8fbafa7059a0c6e3b5af050d597 Mon Sep 17 00:00:00 2001 From: Chris Date: Fri, 19 Jun 2026 15:10:45 +1200 Subject: [PATCH 09/10] chore: fix webassembly builds --- .github/workflows/ci.yml | 2 + bindings/typst/.gitignore | 1 + bindings/typst/examples/cnvx.wasm | Bin 173353 -> 174551 bytes bindings/typst/examples/power.pdf | Bin 16844 -> 16844 bytes bindings/typst/justfile | 14 +++ crates/cnvx-math/Cargo.toml | 2 + crates/cnvx-math/build.rs | 12 +++ crates/cnvx-math/src/lib.rs | 1 + crates/cnvx-math/src/matrix/dense.rs | 123 +++++++++++++++++++++++++-- crates/cnvx-math/src/matrix/mod.rs | 4 +- 10 files changed, 153 insertions(+), 6 deletions(-) create mode 100644 bindings/typst/.gitignore diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 1718e0a..3fbcf92 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -121,6 +121,8 @@ jobs: steps: - name: Checkout uses: actions/checkout@v6 + - name: Install OpenBLAS + run: sudo apt-get update && sudo apt-get install -y libopenblas-dev liblapacke-dev - name: Run Clippy run: cargo clippy --workspace --all-targets -- -D warnings - name: Check Formatting diff --git a/bindings/typst/.gitignore b/bindings/typst/.gitignore new file mode 100644 index 0000000..91d782c --- /dev/null +++ b/bindings/typst/.gitignore @@ -0,0 +1 @@ +examples/*.wasm diff --git a/bindings/typst/examples/cnvx.wasm b/bindings/typst/examples/cnvx.wasm index 6ada7b6d9e4d1deb6ccd43bd84e4e33591aa8e42..14d7d3727d4707d54781afccf3775550027b078d 100755 GIT binary patch delta 48421 zcmeIb33wIN`3F4b+%>uPCX;;+=S~8IJt!oAAQ@y86{}WJv4Vg~i0q&?`dHevQvOw%;>G|L~)V=-;KCq{pn^TuKni@#n()04#m+IUSH@4ev$ zMBgBjuhYi+C>M(I>nNLFh7hFa8}Ot`Zg6rTw?B6GJzK;trjd60od&N(Uy1avBMKBqeHFX;1lBPhYWo{UV*>t+4`>JOoo zpf?!d-YhR#^JeEzrT$zm3!oN;Y)lITxfgX(9f6=XKz;C7->{%VI9Rnuqv6yupay>9O0OY3H6Db^rfTsw(rCJWVc?%B6b|Ap)}_9UBf-StniTiJtb z8GD#L!R}#?vDIuHdzAg2ZDpI;IcwRd(?_2%^8AY~7<=Kx<1U*xY4U`d+4PyQdNz;M zPOZCY+TCmcTg2{Q+u5t^HMWB_v(&X826tM`d_;8K*Uj~s+0$H4PknxUeZ8iI%w{L{ z;MdJoh37oJ>l5n+PXXWliFLo{Y<>49Tk5sp5iUHYRzYcBc>3tOC~`!EnVQh5DCZ#W zX)&48)F|6dN*7wKfmHR*1V%cm3KLQzP>gDlRUo@5t3XJ-%7hTIBhx<0WbiX1J_^YO z3~iVS$$Aa1*=1D3W3K&pi}kJdJ{h`Mq$RPR2Qi%8=_TnT`=LP z;>Ca#pJT4wPnAbGui-T`fNh^x^ZmWT+f*fh1|nYYZ|)##HDk2~tKzB!7SY^7s^xqc zpB@V30Iq-_n*(qPu5Y#e7TA#KOcDcd0}8lSV{kQ1%q1zeMjW7sfSjs-ri|&M=EPrX z_(gDZgonK$^Z2+iH?>bh7R0Qhs@&->`$rr!B`QMS-_|S2wi6Pa z8Sx7(B0(epO~iqiD;4rt6&Y{R^n9JMf$C|@9O2dzP%sE&^`wOAtdp*g)fsG9MADrc zMho7`xi$PIle0U=2+pKz0>}gSC~9KnwQ*u%QNM!wwKI)4*!`->omoA4IY3hQh*4lr z!CHc#q7lRZVSFtJ(>=SN(?9<(B0Gn9MB^l)`WHn1{J_3BhFA43D7^9h`4jr*j|QoK z6oo-4&_;s@1GKMXpBtwag=5`vIuQ`h&3Ty~rMcG-S=qV2rqBg>^WzvEXzB3Z^32Q^ z=YhM1dB`PWz_}*>C}`?y2TgqnUbp2>&Q^1y2|k9n0*!y?Czhw6Ki~6-b!tJh-yT8* z!aqvdO92CjBG~AnDix$YsDg@x%4uHvv1J!j#`eqnz#EhoKCm^yWM37R@iZsjc5s*y z;)7LuJuys`1x%M-6$C-rPH-zd!S|C>Wh^ppSE;+5)R>c7Mk9HbO5W=vJGo^nQunCT z{Z6WrS;iuDpGrOGq&k^pEK(1s)Wc4ylQ|}%k$gxc*Z)H(E{*hzIV%UGlyQmOTyQ2Q;pWi*0E zRPsVaSZ;0^i`4lPjY)M~nDz({3buw@0rou$=t~7LVZgx_F=2%CMvtB6` z(08nGX_;`?hVT*B7zIWnBM=PMMgvgCtTCmP9FKXW{gie>+@noRe=F^D=31hsXbRY` zUTbPyplKmo%^_YO6%0+J;QeTpY-seJ)NKosPNNv;~=0W+3j$GTNqF;1~|e?tA}F(e%mT#BE75LonN zs8tzfa903PX|#Gqv`>+vWl@ugVT1SA8dB%79*lPEx{IbJ%8020{CuhrFw;kqNG2f2 zMD+?m_19=c=5}W$(L$|(iG~tiZuB4ERBfi$k|2ZtcbY)mfH&mFR3&N<#;h8wmgWP+ zA$)8@blyj`hCL)M(5j0(QaDiwG{t4!cZ}S&YPm{9{~~3bPm;rh$bwQ$E#qY>X~M^b zLa`@3B${^Qut6;EHnLlPV#Zn4 zQtF@KM>#c2uYR}vMF|_^DoJL$X3<CFf$Wnjp2xZz}Ir{_Tg^mi=BCd*;mc^C#H*-=)*`-g5Tt_grPT z6HJQbc8^Z%%lu^Ea>6Unp?BK2bZb{c+^OpIpG_r3Ps z3U+YR!H>X)BAUK}5=oQQx_I_doRK$#U01E%Q3TQT4>m*Qc!UZqRG|tL+W(PV(3N%J zhaa$qc4fWU%BHzp*-sf;WgqOu^4KrzV0YHYb~bJ94st!di^U$;la)wx%dw*Q~va@@$v)JSI#NO;0bmq^!*?H(pmbd8=twOrO4xooG+2R{V(oX!>|0J|Gup&91^Tdmvy(SHMquD! zp$mwupwP`_D#nW#Zxz>o48dBes)BEzl#KeO-c{^$M$KJc&HA%?yQ!L8IpVG5NIXTW zP_6EuP_-sh=p_}(SD_bF2-Q>nnpDW??Nj#U{aD|a3up#Fx$nmXdIBObGvz-j=d=+! zx~DLqk?+zWuc_;6XvRU1L9>A+1QxdhnAuf!+rGa)D*?hZ^=FsryA`5(6x#mbVY|}+ z*5RCFgcW&1Vhf^qiL_pH6k1r|y>+Ep@S?~Xj%w8S{W_P3gk3p;`Yu7`l)m(ow{-VoE zy|wyNwuyw~^j8)nrr+L@-bLEzMR>?OB0F4-%AWGpng?WP1I?bkUQInxpod|WO7y4e?S@lXez$cENUzp3 zMGPBN$RUR3RmdTf-`F3V$}R=jH4Poi9^`Qv`7UVqJ?8ECH4SXO8rV@Jbb`cfS~P@x z#i8Ua9L5HuzD!Ut{rWTYc62z4azNwLVXTrTgwsc|?y0W+l@O17QTx}!*=Q;>_zc#; z?mB|y^Tv96=m^%4BBzdE7sV4l9KkwMKvuV--DYKVe*3PGtOIg3tslu!85Q}_XqIcg zF^UyYw>}-Ex|w-8tEI#nPG{BBXw&Jc(RWT~Go3;;^px#M83Pv9PWMB0{u!zy#?YTy z_}Lk(93`6?&S0Iru>{b47chTcd8Oued}*7|-1pnzSPN?X_Z2gmmECb6HQmq23Oh$NF@W z9Aw=PakKEQn-ArbMr4_M0#Th#`#_y1{wOU=Q23p2t3nxzaSDfrQEWUOV1i zy%kBgaSYDilzToK#)5|q5iZ?3&z^VzJJNBV6OvFK!dQ&gnAwg>=z+EEMQBaa{ITpO zu95FC4{Ei-Q^E7;4>W)KiwoH!jt3x~SV4*1k@(a_?A1$^%@9xQM~MR+o52@O3ki>z zIzXF&2h|PtgZ?bBbD`;X&@<`_M<`F|D0%EZU(AY9;zjIlFJ{F(;z6f|RH{Kcg|*v? zJ#ri?Eo+jEdb)UOhIoDib=YGzJKfpjcE`ST9IIfq{meMlFQLVMBG|XpD1$W0uo`8; zNPN_fov0VC485t|9)AfdPid+Dj!Rf^YhdbB6S-5Gk=UnIGlCM!}lg!F4|f!o}JbvFO+*Vl)d9wzc#7y?o>`-X5n_S0?8r(C6YzJ1wvT) zH5b^oPhfrcuj=hB6WC@4TAYZLYL^4^qB{@Y!iwytE@P*)&Jk}0ISTB*zs~|$pMua2 zNDayHQ^7}t)M9(=M0P&BNl#5=t6EnR9}9L73`k#_;N=RSbmJs890k9e#LCOI(nyqt ztBThQ;n6T0G8l-{TPo}t3m#Cx3PTtfyu0xZ?3mzPGNcWj>&LJ$sXj|-%MtD zX+mvLFBZ5GD?*{}m$P00Cy#yELygKT+_ewUf_x*(UhG$qL1@0>m9>RZwwf^Na9FsdfDo;2iWUh3hY(9kz z>-vxkI_{?`K9YLrx|=)^Ut%*DV_&BO5jfWrP`m_@Urk}BrnogdV6$2_6RwKcwe0G! zbX5S;f!lZ3hiX|@{(QYHrm}HAn4{_h`FN)jlU8=Lj&*H~PlB))lx$yC$GWsmj1$?4 zI<^-xdec>GZ-xVg)F{|vL@&AawZy%fmWAILUz--o~!)CIk9B;KRoJZ8>bQ`)Iu=d8;tZTr@WH-)colx-q@Zc}h z+k6f?EuMD99B`KV?3HssllFaAv-ce;@L{Wus?)2Z%)hQ?-A-0y2(NWmxp!GKBb-hy zU6-iHp&jR&n+?r_`QJ(sPIH;@u<03BX-4I9)>3Woi+Ygq00c++qg zP_q%$gT(}fI1lj2B&T!t$Rv-$9#m3@IvXbUJ)F~8G;AJnIfUJGEh{nH>Xg5lmgA8> z1Nq&`9&Gy8wXFU`+5+CLUrE-d@XsyiOvUm z)`GTBjcW`W7O*m`5$s;TwqsboxS4fF=#87%b1qBwMX|V+?y~fn|F3ZKn-;Paaa3}& z8yA)Afw!{S^hF7rmIM|e3+&f!1z)_ACWzPAPv6FJ>>jtV4yDV;f1J257}uC;h||{I z#)g2oLSyPfO!|Y{*yw;$#O`}Ln}7v9`*v)=Y~7)_+_?S@c7>}BeRT)cF>Sm5ooung zJp7SNFb^7pWFBwd$p*0d?fkpgX>lI%pgrj>HY^%vDH6r}A~Y-#p@q0lC6Dvz58A)K zE6$T>1<|b*gYYb7RVi-92bzW~hQ9WL&GreR`2X7HVz0V~O*uYiwtHHv5E|hqi_KyW zH*L1CQP<%Ckc4n#R&jxcA8QBU;15Rii9}zoWykq$0<1~f%eF)JvE0_E=rV++e)ln( zea~xq#Y@?(lO@PxLVME%a#u?W` zkF--Z&SyV(gbigI?1}~~hdygBY+z-h*5PRVYc0Fsd-&(XQ@?Ea(&pDv`Q9Q*EIYW*vb<})a+y|2PT);3myfRxXo^ORATGKN3nvuwCSTq z*<4oS62~EO^HlJW^yduaO0i#fjPd&FD%V)bH4d}g<8ijVb>fwjScgP? z4SW`GzmHtQ#_%om_Wf&M(*Fi-`ZcUW>q0eDXebK(Z4E1Fop?4So`=Nf6IhD7%^v)O zG+TcBgrvck3^m>R1Wqn=Yc-roD%r>J3~^$CTUR$jX|FK1K44E<$9kr^;av~dkFH}E zQa@48%GO;FR1q1b&pjD8Oyl5)bA_*-Wc^x~jPsk*pJGE=C&qUPt$d1g!^n0#B}ev` zr&v|%tQFJ)qOgVe87 zo@T%Np1E$PXV_o9=PJiGvKdL*?wgJ5@KKb6@0Hyhi|39aYHfdUpjhKI!r9OqUsKi_ zXx6k)9=4>C4rrE{d!;V8)&9vw_6+9Vktdy=Whcc6x|P@hK@Vt^7RNu5L-w1`vgwz~ zJtkP(WCPT83fj`PuG+Joww?$LyBy2W6)BGS_UW9@G?&{N;qRQ4-aGOc6nvCpxe23Jw1 z9y3zx!k1Zb)925z^^C2tAAO$v%}q+Tcl?^=G_C(NyO344s#!tzc(U>9=BAovb~A0B zd1(u~kL|R7@*B2QC73{~ia*b~Y`31y0#P$F0AG2(iHd z8$cu23n~IL)@jV-UTp(F?cud(y!dC6<+G~EU;hQ$;w0Z z_t6Hk@DNRZ-9ECD6{N@r87zIBjbX3b(_UwHv)Au?gB2rC@&=~mb^FvepuoItU-Sn1 zdEx7%XJEC#J474u#vOO;dUW^xn+|AHL(`{kFfVhO@XI>Cd6N}%RZXz#i~z-{;^ty< zI}PFUMfhgwj!^l{T`UU?{dgCAylf?lXj8A%0Xj$JZG5oJG0j!$Kd!gg}R+-7e5H#Jkx zP8K)E3YjB(j2utO99hT6v88p6SmN;2Xk?bOV^s3~adRA$IgZuUBTgN$qXX@>!vv9I zb!4f`ajfw^A#=0`cn7L6^F^8E)MJeIrZ1=?*~iFX%N%`eGaqkEeia0R&Y`)`@0 zyiJyt91nhZJcv9la~umIn`Mq-x0UYg;PY;or5%j7uADJ7wfO_k|0x>i?Vk0yAS3JZ z8RygE%$Z&JmJh6-&Me8kljsL$OL<%s5FWeV9(GFQ0zx;*1U(KD#6q?ycSt(1x@3)5 zbTtKS1X#I;z0Y^l+qb-pJv%$=?bUCy!H$ALImnuE*MaS?-;OIUZd!b0_sn-#1AiU6 zd){GZwkjWAXFum%HkR62`!03@?W}Kl`CZnB@n`GpZ}(zNdcIxp9=pCm?&_p<0PIqG zl|~bVByJfegINxSr$1u9`W{Xx+xA88v&qoCn%|e!k-H$8ntr6c#n7@3U?4tXkNbd) z%W0s>X-$A6wt9cjWB;Rw7hCUi^xL0)z_xNMUu@jR>RRfr{tGK}hDV(uTR3hI&-e@b z3CWEXo$A-|tog}VU-*z+)G94*8~6H%iS#|KdY}F`tg2!uH~qY9wVJH+nq%&^HG&Dd z7C9X(E(RD6KhR7o(y+AYoJTo5xLd0O`{Cn##NM%= z8IDIWMY{gBf!ou{PLuCr_C8MU*4rQdm7NlE;rSi(f#WNCExN<;@|J#$p}Ed0GV56n zE#G9d12hntzPxG6C-9HD-szNRC?t>x{#!uyg6MJkY~S)Xwh3rZ^LL46KmI$`)v)*B z#=o=k)Atcj-Qyvs@uR=9C?%$T%AQZght78BAa;gsuzMe5quD0=`h(b((q!*B$ZkjI z%70)F@UQF@|6oS%B?nNa;5{@L0dsn>B?l42VeK9oZx!~a@B#bte}IO5WtV=&y5uVN z1|6*dvS3dN@f&O^z2q~N7px~$mbTAa_ZeH{;GpYWi<=^)|76#-N{d^<>;K6f<*(J- zV-I0dy=(L|!FaimTSv-G-|fVltUsUJjXh!o`c$+1*5E!xc)z$$=O*Ld>}{5@TI$a@ zp6beTTmOPhJb_7A@&)U`p0+Z&Eheq$yf4}7A-TB;vkhh4uKbE!lrjHvjL-ZQhw!~3 zrH3~EbNj)s*rdR|4GISSBW!RxdCot=&Wh)G;s~3WDeHj(DWA0>tnajYe$DE0?x6Z> zX=1326cLhb{Pt_Mrl{Q{m#3*wra+rkL|L zQsO$r=0u8_Qj5q#3MNh@Ii*s}B{Ji^sutLtVZ^O+;+lptUdw{bPDA!4&KC!FI-#bs zJ$wjD-LM-CnU&_=&+TPizVj$Xeu8fq-I}j@08Y$ZPOSEAIX#1*6SU=Gz0;H$A;tO; zyErMf)nGFKisaGSM|^yB9JgEf+kWNee*#X=OX1`gemcmr?dJpB$a#}`LM}aN7s_oL zkd=o5{PKL2gzX`2l3ze_?W(tD1$h^z9Ur`ZbR0cMf5Pwrq2?W}+F6srXYxPP+a;;I z_!z~q?TM*;9DfsAhf?{iB(cv3@$S84FQDkSo$%^AsLLS(5u{|V{{u}6c4x;)?(q=6 zm+#+g*M#{=aYt>Nj@}mL)ooI7UX#q_+wDJv`LgV`O`}t{M);yO<>JHdk;X^HUVDhT zl)we3!}ulidf6W-i{w156#w|Z_Q)lw_7t*smLs$Z*G9dpOKNGW>`ZD4d?*#agw8CM zEn%063&OTpR>_vy>5P4LJ_kP})AA&b{!XH8+1Z*}ejESYJZo^rK^!S9+C%38clQ1} zhWWN+m~T4_Gr_xI#k!yCld4bSzq>!;g&l48?f}!esrOL4AL$*%d+A9p6uNE_FPX7H z{UZsMhXc$Tv$ot?+W4zJr}OVb=f{<{2|VEh^QslU6U;2tzVA10w-zpT+={Ip{S3&% zv-|ZJxtaXg)i~;dot<5@4Jt&sIYMhy2--7+R;dv5Z3-<@A!z0lTC&=DvEQN!Ig?rl zv!$6~b?eAN%7^8{&=STxF~q}tbI|%}|4WZAY^^W)k78dno;Kj%54W6qTu>=Q_+L+M z*ZAbVG{u(@W7n}8*mY^KGRYsXp?a#}DVOiO#=3A&tnvQ_J?@FE!lI9e;O+%_>Q#K1 z6JAPBOwVF^rd9DpPIw_bv-^AM>6wC~H1x!Q?M+U213j@Xc0E0Di^f_fyqca^XSx;- zTA8cj*KwQ=E`_QNN6L0?q~TB-MQuhDt|j78!c1rlgl>yuHghwtt>QQ&#6C`?skR2IblCZ9mG#(&PS%6N*8FgN`m|M_ z?#^q|AGw~_V4XYST-l*<@xww4TXY}b*}7v?VF@TC!>p{U1xBVFNhRgiLwRcrq9SYSsB`*l zVo@(uxSmA|=((0fi&S%H5Js&TjN)`}l#7CB2A*m76$@`PXq|O>Rayq#&^f4)m32m% zb{@qe6lf!AqZvYn_HF19(!RYS4F~^FJ_LvAsECGQ8P=iG z`@}M!yyL7fxx+yN$eAHh(17_Z{;|=KfJievbBu`j1;T+bk${8_b{IRIEv5?~BF;$% zY}nzliR)x1JTvKffC&Ae%M~I=LI6ryv>3W@6q;%Vwb6|l@#8(Q5kE#=;!K2U3MV0> zrPj-%D`OeaQb`Qg8XZhN97R2Fa5F%F&5o83R2MS?P5Ln`?;K(W72s`v4KmD(_y}6c zn{T9wj4{}nl!p}0oLUj-=FuML??E>qoseHQXGqxL)&QLFj8jT1F+lHyaV&424%)^2 z*4+pvb)%hL$bFwFfMg7S$O3}k=xsF9BXKp3LRlhZIA9$RS#m^~-8}iI3KNovAwm5n zxM(6v;QFOZ7d)9}n36EyERj(c4Y%%D7APTH3t**lv^k0g>jS$HAfFHlplpLh%@Nvu*hoP#$(s z+wv(Om}QS9)`+1bw_)yj4C!WP^DsT##hCVEML8Og(1ncD1bKF%`4}WUvjFGHS|d{g zt+8iy?TxV#?s}tHn6gYFlY%I~?vP3rE&!vwxdk*O1t!*^Qmv=Y$_wD=IL(DMVC>bA zFcG9krJRV_Jhr=a{g}bAI^ru;d?+3`;9Y|U#P>ivLRB0mf<$o@?}G=jjKd?ZijT#k zw2GgH2UyYBctDGPrjBT*qIOm*vNx+1k zpwWerYBfo&*5Y$Zt&7gBNLz|XlZ^y~kYgF5&C_dISW}*2P5I{6iZ$g4)(v_~G?#3n zB`jJ>rG`ec97d8aCxaMC2~9@6_1fs-SU!OtaU$PgP9-#%u`?q1=vEd=fMtjCX`-^6 ziSwcaf?hL`$6+>UH4Vh2s-aWXVBI3CHd;iz=x1Hay}%onYUkXge%2As|FJ zs)xa9(@|ebiAgwgk}fccbY#Qm!r<(IYj}%pn$SK z(1`%Jd^7;Z8A9Z=L5Pb&>TBV!NP#*4FS25B=Is^v@%Qh*$lKqzaXcsS*VOO%AeduP`&z%r2LVJ z2j&agn3!+23h5@;fJmDP2}g{)^MxaFvLNXQt!GNTE)_JEhLoEwTzg4zyk^m+=2XRk z(-jMbtyYYJiqnw5al250Yh*E!5=bStUnB}sMV__jqC#ID$Su!$;-ZsdV90Z%3W;&( zhE7~Tj!17cB1jGz5hMqV2x0*va)Cq|2Syh$bjT*iHaLE9mIcZjshVgBmyD3=2e}0` zUfqZSRtDu$=+-3{8<3HtXhcRAMF9oC8JA5MIziekZ#7)pk-cdB?&2Y_;&>kgl9Qp1 zK9YDs=P<7VkSZp~c!}k@kfEL93Fw0N7At7c3PPLDXb#a+y350{p_Iq~&o;Nka0eIR zivCZ^hF2Je6ibnayLR~ z)uy^qk*ohJskVseC#_c&UM>r(T34%DVIE+BIx0h*ddTzxOVN#O+(;#~1JYner;M>u zFRh3*p_KH2;P~=-Uaife$Mj!QTS}5 zi~t=`xFO0wEx~BHGYw~8u3Sa19O6&eNaf&YlOob|1nAR=KFdXhWsi+b1c??cX0*)9e$QfEbwhV|Y zR89(V771Puyy6Zz=|wsgc!ISt)Bux=^ny&1%9*4NKpQ{}$=W&koGtA#Qi%wfs)G{9 zs)H`-4^&JFor1YWH$tb5!4W{{>%+kqutr5u=$;SEU9W}ag}>;ftXrpU!fuxxH6kEg z*_p6Nr~F}&R)?|#2P1V9E{4GXQXa912S!L2Kw)IGDkRm8G?4-zpqmk>aXR{k zc90DmP>u%u!zFqmG$um!Ie3|;uJ_64mkPk=n~fDPpOz^o%OABYgD1zqI&p}Pj8!^3cW9ZVS3 z-U&rlkP}8;{Or;w)LVwB2_@6OM3DGG*paQU9sxayDF7vg z9k&y<)f@R#fu!BBWHF$iG<94nhEqXE%otYpNilc14jS)cz~Y3uBK;z+I0y+SPViX4 z;?(K23#3Oi782#_GQ4VN0U6Fy(gN%a?Y9*62t%2pkQii+;w;g5rUuR#BlcMYN6Y9U z++Ynm?px1jP@Es)9;TmT*b#?=G?JrGCN-UXO`4%h!aZ~W77;O^bkem*ugefcQsqep zd`%CTHRk9`iBlSEF;M_a8D&T`V6&_AI%Mp`9u9!t)bDLG%R?x7wjDUoesl}qCF-9b4QFFclMob#cdX$i8lnvGk^G-jkq$MdlPKA;QY=CU0E^Xra z7p^?r%$3qGLV6gHxSkss1PkI$iKwQV-ZO~qTy{sZ38JF{umfL2WU8+aCAk}XoK!B-Ztkaf|OBG%Hw zBBfbT9orlLkuohuyAz{z`%NP}j7|Dbs?r?f9iQM@nM9(&wI@6#**4_TAzcaQ)YT~a z#BnOtuMKrW^&5&-VX6R^Xjy0=9x&YdsZbwPg~J@s*|BYE$XEOUeJ1{Z(0L9@6#=85 zxDil`6xXw^o|4_6h~}P5ts?TcxDz2c7wXTpR!hyWP8y%bt;)6Zp{L0G1(Ns*b__<7$!tes(T>H z7O=v{C3SPDE=t|K1#NK&L829cEht5Ci^MPdG4nu@3lm6U?n#P~Grkqjb8&jJ#A$O+ zLo-~P-GJ3VXmh&^&0mQ)_r!qY^M_VYjI%ALmC;1@cC5ah?4%%^U@dD+VfP&*XqZc6 zrEV21{_-SuxqCc8O&7r78leD(o=C;WAjK7V&0vt|#06-!r0^>cVQWiW?x1=M*HzUI zvue2hl*K@dsU5>DZjmUmj$G9x23w@XdLZj4ULlMkSs{_Kd6@Q<6VVmvU%LG0I%P&& zi-iY7VPz7C1;m5<5o9hr-8^t;Li^f=J9m0$OBu9VvSx5WD~sj`nROwllsmxW@xUv|KU=4NXw92GfGUzDTPl(z zTPl(zY0Ro3IU!C(byU#|RP^oip0NX1S1tw^F$ss@D};&LUv32OKtJ)-&d0qh#a(i5)}|UC#6|95 zEOGhVbu^8#3ZFC!X}L>&Im5j~7GDRUbu@pHzA5k8MpogIK4xzeP9BIQ!XZ3>OE?5g zEAa3WlfHqCytGd5MLoE7h~~9BStjUOIO?T`9u1;h2El`!O=Kc*6WfCW!rLO5IURSV zEzrOiEb12<@T5YgaHKL`%N6pEi1={f8KvbL=sUWIo|807p$pkl!PBe-tDr%UX{awQ2q*t+O8eRxI*apU%BK;@ga znv)_+*NnOI?XnuQ+;X*P5Z>WENls&25rYGfB&hJ#g_l=*D&stpo0TrhJRl1>mwvYE zuq{8z8NO|~zPuQ}4qqU`*kn=>6{KY&lh;NvbX+ckSQ@%RV|Sc=($rXKUG*a*?nP9J@JfYIKZrnSK1YOSLR zJ+n3@gF*$@s6=eLrK~Gd2&YOZ^neNhMJRNS3W1(5d2*(KQNr)o9nd_qAwfASw9{K& zq83#d15ZK?nPKAslh`^dVrK!a2=41rWAE1Qq&FeG&Z86V&J10 zfVYeSz*%PcodE5lG<{A0`{EhUNM#FCz~cnqen5!)ZBzrHH?T|v1fd<>kZdp!l9m8j z0a&G=0x;;jpciCtOUlf!;BciiS#YUJ*(4NqO*ClYu)koW8EF_8ZU&QPS`1v+L#AIk zD9%IqvFJ<-_Au4SxPj#sh2a2M5s2m@oQG#Vo(0f@%$_Y1pN5I=*=FMRwoLqOYQ{RD z`*Q*I1daj`HN_;6GRFD)b)#iD4E!HbX!Zvrnx04K%{m{ z3rSnAD5VFsb&xv>%!;~`$a<)mF_Cm2K_MZ50m9T|kKs&^5vhnElt5`bu+mHUVcQm> zd!mqVNrC~yf~l^Jb`TlpbPj+}5EX!uNX+IDf}krfD482WC;{k12mDzHe~!tbz)y*Y z*y2gPBGm0e7u|RGo#5&)Qe%SF0##vheju=7IXQ%!5#pp1L0U3b6u~5rw+xdIq7Ds@ z=yH{;V+*qRfNUMxAlpuGUx&-Qtgxg+AzNOeYuS%n<*+Mq%XdaYF<(1#;9=qp?h zo+LN`1boz^2)htQMdUe9k{NX8JIM^VZX-%&K)@)O6E7NTDH}vta<0+$qEI#{%ihGL3L&%qy2R;LUwb{8jzGyrHRHV>Oe~rVlZ9ikv1`l2mr1I zE=%o7jKO8y9pileM3?HmpT$K13Nn`^cBY*eVrJWu29<5b;Nr}N2AN5WEk~r+f|=zK zmgT}R2Q15jM-FNEa9<(~ds;{e0zfX~gaHr+=LbgsQG&f-`tmVNJv3|}mm7-`(gg8V z@JI00kmeSOQ7PhAM6QM;;V|Y99&92)Py^8dNeHyzkVKq`OP(ARU?(Nuykxy|PYrqi zzZn{WszGE-o*Jmq-~{>3auzVZ*2q%C9*T(&(I4?WNr_3s=P(=aB}s!Y_KDg63!k|* z8IS@ZDH|{c$V$*K8i2YdAq9Ck8({Phy8x#lDFWo54fX<(tBqo-3}g#wq1m##G#VHu zQsrZiRqj*4U2+t(G<^$v?LFFZc9ncLF2ZD=$v_FKf zADniWwqYKz*4%@1YKH?T?g7i=sNAE40}%5d>`P)ESPXKQ2li*Q&pl$qwTOM-aDn1G zO3y{nfZ{tp?A<fER)v|;U@nlXD(?Ew1a@;c@{AUs9P{fM1VB!Ab3?e zP%e`;snI5s%sQ!tQ>+WjH1dnk5V%DSkdyZ0sa;TrBasTTk}3408>QiFNnVJd zzM2PqktnV}&X^MN3gfK{&#v|e*N7Sn&MI6h>WD0LEGw#w5oL)5Ik7y@myzeJ$z{Ob zduUZc^A3Zj*&}W$*?A6?$dUXfla>KtLSU-f%m4~tx`C33(~Twi<6{e>&ZeKkZF-ZNOY|qVt_N-K+ky>FybyZXb5ux~`i2NlT5G zRCB)SlQm~SW$ZnAuhx111h8`Us8EFp?NTA3Lfciy(Mh+ekfW1sRw3g2sAYo+(M?%k zXMJ5xm4I$!B1m+LCE^=N=uF{*gof|#Xamp`y-H)j3qOP6ccI^O^cUk`hZi0Lj%S&4 zUg4JWsfMqGBa$ep;FD-Yns^Hw7}mu%7sDukgNL+BbP1jmvmd^NLL2maFp=pkVt@J( z%kjfpui}%WC&F4ks9P)rE-36ZCp;FoUj>#$I_U}IZx2BJx&}K|b~x-Yq8E*tKO0IN z9_L|G03Kn8XsAi_sv9}cG_c5A07`ECk<2xZmNk4Jk_?gO&>@NZG^r(~iCnU|av|)| z&nyAmvDU~YX&;D!1;Od0-w1yOxs_n;(+0$}mW?-b2rwU;LgjuNA1$I_`JocI)JB^M z6uOnf9Bp~^Ay5HBe_$4eF{e{+Y!9DlWMQck^~8*9j3foQ<*;)Bd>RQ@LAo;% zB+=Fjp@~AJd0}0LG%pA)(!Ahf)WX$L=oTAT_K>%MQ4U~$v5XI3fAd|!+U33@|=?B4`Jtk63@|%qNQaKVM zm*lG`pg(y+239OLqQW}~>uC;VIhJ)Ai3gYrxOguHM+5;L6_1Ma4Ohv5^fNHxlc^__ zQN!&JZ61x5L68Gg%L#L9qj0Cuo)B7qlvG1q%|TaloEpbOs?bMl<>;%r8it4khUd6F zPOI$}^JvkiFhCb@yUTt1paVGuY#Wd!d$f}%16ayIULB)pimGUx4cpB}VOp|4t00nW z3_IExvxnYgtKMOcA+23!iB6)k$VK{!2Go^H))q`KqA|!+j6W(mIdIB!T@N{;g!oq) zStPQ@fK-AGVxWpLGWpPPLQNQZSYbjMkrH|`DpD54nht<0i8D%X zm5KnV=aM)J6>(N5^zYJ&ILm2VqMQ_Lr_M2va-zI3kuqV#Xm~7XaM7-!Bk0m-MVB30 z=&}%WSpjnE;?iX~;1j1ygBHO~lrCef36rL&D^Vt)p99kg1k(v{eguvWR0o*lq>hv; za;E)tn8su}aArG3Kpqn@;u9Q4{Z{koOT^L>nc21zi*Aq{hAT{Wj2oXB>qpiNWMFrk z>BOTlS;F@NZ*(MlSW0*SVIG!SY3X{kBnh7{C49(yOG)_F0zTi6G8=6a!zNPX3XKEO zF;mhjoFy}j>>@}Iv}2Ncf>Z3`EqZ_<_+f zN5EI8wuy%j4rY@~R;|Y3%@T1Pt?jfHnJn1|B%PAUhLXwHniH4F&^6TU8>Tu+m@TmV znyg_+$223C_K7&sA{RIkgVIKuHX$Ks&!&+n<%YyOa3~_B3Ilufl*J7f7xW-$&x%ef zqn&^>X(Er5Q7CVaZ)h=FTF99Q5)19@;4l=!G0+eQi$RMrO>|y~j+Xh({Hko zQ6yy*4vG{-b3sN*uF;Vn;By+ezsUOr>mF^xBB$S|B4WkF0r8j#v`aH@P&@83G4yOm zN+~_#JXWNkYiNx~bVY1A^`gZ(pID^yQ-C6`r9_@n0==O^?Uu==^Ml0S@j&-MZyey1 z4m()om=!uFCK;&V9?H=)5u;w9Gp^KcuWFhIz>jK>;1`{!k8*N!0}7E;y~3cMTw-!Y z9#9OdF^`UYYeEQYeu-e!MNlGbOmp;UMTNA*1B4{0Fqd!_#6b^tu@|rXl+7|C0%L)a z971wwD3g#GJLupC!uWvKAR!l*X$1^h^fZ7&(NSg#9VzI-SzV7AckBit4(#5fbBZ8V zoTeZpgvJG>D;H2Ijav~CX{BQ0Y&C8%9+-j;_Fr5D?`okV_l_d3dr3*1GXEzC% zhH|-`mm}i zU0TJM`&QBtwtIv_7TKCqge;l0k5V4@x*<@2xj88YLj{v|2Hl1C!z!TbxB?K=ZmI3E zL_1Ea;~s$r<#jN3CAzA>r_#-Ri80>=F0Jy$=95{g#Wm?!2N%ZIoeftfLRV6|9a4s<@5`1@LJ8 z8YoBRp(Ok&hUG5J#$%L%QJGxI%ym{%L8lDHbFPx!*g^=4ddLEdXfN%h1`B}$Y7LQJ zwU(wrt5k^C3(}UYXfCHj_(**86})m3mD`*cns9vTs~GcGR__mWoPwFcc+oJprw zQ3Bp&xrxXV^;*x}D`K>2ot|jn#qT0q#90t~dmO7B4p#E;NWT-_u!+Up-$YAjU{vw_ zO+XGzt4}(VXn#{c_5qtg681QGh$T~xh6)b~N2PK>`;eeoM6qZ>idBm6!l~sak300k z6mlYAZ342AM4t~=_#z|L!?1R3v1_uh3kX9Y?*v#MG$A_ajDeZU=uElU*W5%;eJjCA zOB`;?79q4@Zj-qjwsC~x${H5{=&u6;thu8t>LX5R6D3H!od=L7bw=I~H+^4~42VY!M^uS*1eQA|`Q3g{boKzfyd{Y zmi6GE{52$8um}jM2fk_}1nb4p^vCB}`yU(~yM;E(nQ~DM>q3b)DMvlzpYw;S8{b9qO za51iPz6_cYHud75pC*fX%{v=uKyaLA++?LaR9G;dMh^x~dvm?!T{0DWE|F?JbV|=f z@;XDA2Rxpei&VP4fJ%}6?=|mEsBq0gr^C2~V{u ztRy@p423^JrXNv7xHys;j0P}TaBA8aFC##SfdQcL0KQZytB^$qptydmo{HeaPHedK z^22(n3}UQ@Ke9e}xO1ra2mmHDtRp^U^g7=vx~{JrZHHFA+c>z&KV2 z^%qn?e?&cxkB6vdKFu)(cw807YHrBOlOH5iU%P;QOSiCBpy^wzzc-Wy_N#tbS<4r> ztJd=VG6-g`(9A`Tq1)hop#>hMOJP6=6~;ykJRlCg!#-ar2NPrx;T;p9OCWS5jbd=7LjdWJVJPTB0?x}}tw(U$;6dl8jC5L>3E~t5 zdJrbr8^*pZl)()S$b!vw(mcofLn_f;U1!Hedy<)Yk?fD^Pd!RL0;)GS6P>V^H(xXV zuobgn{!=qItiX%D)p~7tPd2~t>*cu3*_vmayz-v0Ev&#HS3nLu5189`Q4!~sQmLv$ z%z1P#60jbu2Jp|2*oB5cw49U|Bw)x`NDgM))L630U?;cKE(zU8Dx`LFtM*dNJx+To zXZ?bwt)YMk&M6t6RsYx++;F(_vD{cOh=x=Ky5EqdSk1F8c}eIF&S@8Mnd}m_I!`tH z4X{p$WSKe0ON)W`F2HRJ`4M|0)3G3c^EMo?4m!-pE~e7xHK8_7znDeF4M zA&#|D$Rlgny7JPN=~Y)1O16z!)aN!D-(4f{MTdLV4&`^vvxcwjhvs5yJ7nKQbz(7q zbU%DvK~?w6TDqw6-)jY?F@IezLuH~aQ6_XNNe05Y@mue%D?pJi)^!TEAW>+f-LRe& zcA*XYm>SZ~a8pI|GcGF=ID%hR~+!>YjEjjDz1#mh$zxC}-YLa$0kqx{B zNn8F6mu7stDY@8wYwm_2MO%NFTxu6fEOnCQ+t`IKea@=f*tzun=O}9rO`8Q4$Kk+} zZJBFjBcEAUZydr8Z$ZvnYsbb6_J#G{#@v9L_|3)(T&Q46V;iVoBVQY+;7}yulA(e- zD3YM^*s~HUtm(&m@r%~jrVgcxUqmO(&>X>WcozHoOg>9J=HL>v7p#c47{+E6kobNFpp_(6gx!9>WtQduK7pP>k(UMFJ z_xb(seiGEE2#k=#C9z#|9@i~mjZ&0?P$P!Nk*Bb0r69|gZ$0z-uHcYw|Gp5cZ_Pu6 zR^}`DY_Zktm7=2U3!qNAVj>ejF%QW8=Yy3`c%?ktpcG`VcY@9PMip}4A z#yVsBApYS->-OzYdOW%P^8O#mgIiH^;j_)iD7y?zaVdq&4(6uBSRC`(OoVNoEU3ye zU%irVc%!lD)xM0Le}Anz8|>mLMu=J}T5wEZ5p=Z{yt~)YU=wmRUbyoz)@jch)CkvY z%AnlQrjR?CHFKX#^B~Q7@Aa(MK^gQVWramTnxAIQf0OD-O-|%a;#!%Q_KSWh~zLb~>ttbE7r;jY+|UGy}v@-+kk?j!T^YdJZDfIa6_(MtIFQ@}ihR)!uK#VkG7e2~*;ymc|6SHES1WHgYCV*l{gD~Mgbd-?IcQL_I}WrVYl|GPOsqHX*? zDbWIdmJ+T1p9N*P_|In%PX75%?4!o{duFnXog{V5*Xb|rmo;}jWBK1P`5Vtzquwdu z|Nh*%@}1uN;OEwZ@AN11`{O%H#nRs)?+^{+^4~nf`Fhb1W@5*bNbSWMerMhJZg;-o zch=^2pXv7pS$jwG2tDJ-xq2tHS!`~lryGl%O0R3?E*Z7XdUWra;6C}je4aJ-y_@+; zNa6RYvR4v-39fX8@sN|b>-(K4{oMB(iWLTsxb&KPap^wxN>7F7&NCCZVWoc%Em$E@ zB%T4z6r5U|sjryVIO>BlB^@l?H_*$G`;Ysh^!Q@`nRqmg{CHv@uFrC{AZ>OD5%#*NUJ4y-j7B=yrl}J=)3Fw!WlhB~^6A z{`PSIzl4TMTrWrn6t;rEe-_(#+GjuE#Tyi~;}F5kn&wugP3!qX7haqsO4|cA`k8Si znf~-*?n>-0@B_ODLD`7};2WZyYotmlnE82M-Z0Nv{rLbQhJBy+AoS1r7rsC?ukoyZ z9bqw54lY_M6DY!Aw#1qc5fCYlVQtxiO^-kW?;Hciu_WynD1rr&7V1GBvK|nty99+O zoa&>f3^xiiM*FKJjWyR>FCXqjRQBcJo&8*TnlA_Fhyc{kEKh>(PpQ=x6V`*DCC6F( z%WeUO;u`yWS;K-$WDU#bHD34C=`65HhSoK{dSp5i4sVs`>v5PM9$fFim9j)#0EOTA58dnLok zxo{Q=5pn&{=|Xl0db#qHgH#v5pMT`i}b8tmEGtPqrI-<8cOlB8RR=-S?v zjgP+{{M`N~n~{98OLACWP`G_w4t-q^#d4TY)GT!}pZ*a|A9;Vx7iA1SS!rLA&hgC# z`}%Yq<@4v+4e5M!_-T3b00+t9t2pQByk~f2V$7XwVzyC?m6(TT@cuDr<3|~;B0-al zf51>p9JFj@R zdYjs@F^HM(%Vc>RBFW@uxI`MsAIjwR`Ual>Ge48b>1`_SBkx*Z=Da`B$Z^MFo5snh zZr&-A-J#i6Wx_jolYM6OZXV1L`%8h)a_}fZp0(s9q`y>-`jzADBMIM zF-ki>F+1`CQsSQYX89)jwGO-|a(~?c3%z^oy9Mvs*ZqbV$^Ze8Pu-`Y3OF?Zc>zb^ zG2P7eLBR*L3(cD7QNa9Aln+l^Emu6qDPF}_@3;RE<=K<=d<0Y?MT(Y0c2htmt2aa) zl2NhK-K8QZ91`~#_ovI8_(k#z6tzg9Yl7V=G)b^KUm%4;D+)_!d=7rN(Ork+I4{);cqf!NaLnp)?u|mVoLHpe}5( z0~zSQy|0XqXV2Rc$}n%w+xL|53-z@G4y;r28XEhyPwcPCc(3p_x%A|nRYfsX<@~wA zEvg@I*MCO+GC$SKwaVDHLml~@IdZWp&8ul1CV{~RVExO6Lo+X?*1^yL+-gLvnc5Hn zT0OoX;#ufK;QYQ^mQr9`s$$7HUE;nvt67u*m*r4w2*Jf)&4#aCxd@q7gzUh|L|CImWN=dg|Zk1&Yq354ch znkO^M+^oip0U%cRBv$Tra@wbL;XA(P&z)b|m7m03vG3{1Kl`3*e!Uz27i!+t{d=!D z&=dEZs^vrJinsgpAm<$kG-XuXW3OJhrHl5s;H}a+n`fFrZ+)oHe-U7Dh6Q#4Im;ng$)&r4C@LSC3zh#em)a|rv}sOqeFdqx5(&cS7COlkFP@^K%Dd z+B^!jv(R8As&K=z5hh@`;Y$!cIUzh5;VS!^!MrHN&9c-sYIt?b4VIwCi8YoYoLFNz z!kGy*K86~HvxZ@F?-6iE3U4aFiT)e(xXeoRQ~>xd6d2AuwsmRg0Q@S_M1Ls(3hBe7TDz0{}2v^u=}8;MaBA;VuGOtKysWSFD!#PkrtB>&xXYMf-h z8_q_UAUkE+<75i7_DxI-9+S`$Gs z4PoLrFQ?f5oW>vKccj|ZbZm3kfbnXmZ>fu|w_UzfnP!hDC zn9VOUt^@4PKm|0pH1RJ0(=Nd?7sKy@=N?~EuQ}WZ+LbGr-t?pEcx6RQ8?0+Oeo_2R z#_u%z&c*L?{I11s0e%*KtMS`_-wXKd!tZ1J4&&#q(X|5n2IF@Pev9#Y48KPFHskj@ z{QijF-|+hizZ6r~a`Eem-$|z4MLQEO7vnb-zj^rGhu?Din(*6--wyo#55NBhb^i5V delta 46726 zcmeIb33yc1`9FTonMr1nnMrQQzMr`Xkg%^|5hOP(f@nd-wN^xk%CJLJ{A!(vB1OfD zUiD%{MMcGeiUbW+)TmggB4WjTDXyh0YHCqYss29iIrmOxlCZR2`}F^Rem{6_?mgRk z&Uw$fo%byB#;fV~zmPsno3>~U(=?4e%nHu-M59p_J)6tFEUM9$CS$d;y-1@(Mt?LX zteu@g(a4l4-}sObmuhGGR7CXpOZpk!%PvFN%OngW@k{p-dl@o%{QiK~%NawI$Csjcf{ZcFLQKo> zAU)OZ)mT~@V?O3f=PZ+HNY=Ahsvq6;P$8bmJ>Jw*6yiK}imwavYuR3J4r)VIBvX!5 zuQy$`AT>2NH;?7!0-*aZHORbvG=^H^_4s|9p=s0>zmNH;)!Gfbxbo@Nr?o&-Gk4P4 z*Uy_er+&`7g%@+pzi8g<3+t!Op0CZYM)8v9G^Uv>P&>-3?jZ&Zeu&kcH06|{C2Sd6 z&!!(c_3U-*tZnQW_87~}%P+|3RyAbku&_~4*>&`U5o5>Q!G6nbW;R>NZecgEyV-qg z4ZD}EW%sjt*bD4Mww=Aio@D9IzUkk=<})7UwY*kn=B1G;t{EP4@dCqZy~q28Uuk?g zrkQ=r*yVWZOJcE@rUlFu6hdsvoORyOQJK@!DBB7uB(#VDIO?0Bjcj2J_oB5Z zMvzp^RDdc&NPWtL0J0;~8j3df@u3t7$qEc@tO`lEhS#i~5b~INKYQBprVK2Ug;g6Q zNT&_x%jA^saZ4mPk4cU3T-adn1tC-;Jc6TTG7$;;*#h&xXQ&sg4}%h^DTp&PCpPRs z#FEdgCsHQyeVRh1PjH1-prPoS8LdIXtzV=Tb=kX@;Dk9167^tNt94as-{3Np zAC)3vYwDJA3BZg-3|7MxFc#9>LaLGi8J`&l3QiGBjF;azpVu%O50rT(#b4mK3kSqwC0o2OYrY-6D zjPKoRy%wnNBJpQlhNgKgyl4U8(>g9Vp^NM*qS5;jibjv`3z}_&$CC*U2tcD#LdUle z9(Ou)PAG}+RL*-t$52gQgo^vB0c42-$ikckP{tX|9YBr<4B)VgKJz-zER$oP&?C`I zG>ZX8^Bsx>A9O~snE_x0l-z*=6~yIo$gQk`Z0j!>19gWmJaMw_+i7A}!wvjIKT$^3o+ePx?IFqn2q$lARq9=~xR0aIC?B0Eo zh>!S0GyM_57FFU>M@k!Y^v@Z3o@z6RDDmy^4$d!W-w}R5ZqA z;?%I8+DB1Fs!IZXgn`7JIVZ)LL(x;e+-`)f({lenuZp~jh#XhtJxQT6^OwZQ+1J+K zFZ0bapKXME4f24?v#O0LXd|&hgORG*m?Bc+ZEQ!VDV-2%io)V`1@$=!BbwloNH7pq zRt2MeXJaiS&7f~Q3t)_qhLJ@3Yag7R6%Tz%nM&&T?{mfcb*pIJmwRRNd zMb}a>8kir;9Wexe{u1f9$> z7NBlasjEMaDp*c#8I9ysDtWC_*vTzpk-A2uZg5ha%rX|K>s9JzC)LR;W0AT^rM5V! z>!T_xqmhjM?00&*MWt+WaykWNEHXN&Tb)$-uHieAgWH{;eAnsUMDEu=u&cUpA)2a6g`f8p$0+5*YE>QvV3E>ozx}t z#&A19CsoFx56e~R%H+&47O6L?)YVR^lbK>+Sz)+U@>(a^$t`1%x<;jLa8jMjG8U=p zRqAFZ)yXVlk-ABxwm7LyW*LjrEnhqxr9N$QlAYX$mh`?=72fFtWpoXfu_(MKW)R+L@(`{e!jFnY8om!EdPv90{OBM) z_kU*9m!IB6$zOkh{IzzLyCy+9o6bj)ww+DqLj?3ek(lAy9UjPDse7zni}Q|B@KY1u zqs35nOmHdReM-JF;fu?6?ieCgzi^JDbHkJ!)tnF>RO}d1REHACknIsenpM{7aMvDE zZH5^MZQz^FGJIy{MAFRE>$s62t)(x*Rb8Zh0M1ZD#DJiNp!lRHR_F6CN%xm zDwMC8-`vpx9{2qa{UG>G8X38+8E*Z zb<8dZ`5P#M*&Jtcf9lU-O;`V^ z2Sd~&yJzR4busOMI1+jFztPzr0|_0pS5&i;*j3gE?+fp(sIhrZHM_&h8tr9$ zSUHQ?Yx}S<2z}Ovg%Qf^%bv2V=Z1Uw>8ldZW`E*kBUn<*)*!0}Qa=x}A*HLRCm3Ux zg` zP!CRc1-&ss%jlf}lidl&=$$jnvly?m8or5w&`lfYjqrMULrtuuH)aTHobYOTdu#Z7 zyl7%y!yAZ;L!~4JG5TxaMyjL6#n=$O8UYEsdG<&)l4T?(+RaD9_jSEJbrdU^5TmSc zQElf&hR{V2#>W9vRmHdQFcIPsKm_>ZGMjmTMbIOBo7;HAyy^B9$YzhZj7RkA?foX} zX@4<_otw9sOTU-78u7PLd4jjbo>j|EVQcMub*!7crxwn*pP&hA+Bf0){0_ATRQ0Rt zAY5hjXGR?^f`Z7>-EY4LC42|s-%RN7@5DSTDSv7_PZCU1ca>Q|VrPwICl1}f!d?Qo zmW2yt#VlN`DnylOFn0W5ZfYRg!f*!OC3s4OH{8X(bu6pN$iOGsC(=-hk!AmBEHi3l zoTpZVjC8YKBpebMXkM>KZw-diF_ADh!i1Aa4{3n_oRC?Oa2Ei8ua`{-X(*Oqo8#D^ zXod)xEyx!V@BvZQ43Re5C^3J;KQ=Ms6D6jn!3de(BJ7(K@<|kh;_WRZtGNVC4e84W zDF%m?1AADMkR2AHKitH}fdG(_7S;jZ1)YjdkC6a76@4yq20_)3kC<}1{210Rnh`FS zyfEJoCLae^mZ!J23d3WvMd30+=rU%Yx>R(>+YpvxxP$>I$gm)kLTG|vM9%OmkufRM zMdTyJ(-0A%o*wNt!#rprR1fl}oAV<^dPocU0s*s)Wul1z`Ywu-c!KjhfhId8{PxWV zCN-mzRs;x@X}}`90O0^=1MaMFwnq}G90?$G9B}Nzz+#vRA|kt|rvQLZEE^r74$DDs zu##;CDGvSQq976uDo7$AbU2U_Q8Xb0dwn12baRNga$O5?y< zx7I7elcKSdfU1HEKVhGb(d56@QQ;g+BaqW+fK1ig)Td;sYQ=QbjD%LCff-uvK?C)5 zM0l{T5}KQkkE6AI)moqj%8RC-EE4vkP4o@`QT-xDw&-GiK7sX!3L-skI0%Gi6LA-Y zv*W}agzJ%-TS#0hqe8Xn1tXFWQK>=#ufiH$f)@-87?C2ohUY-gLSgGe z9DATe8{#1w5(6AtT*FU8G-5~`v*Fqtlqmyu57%ZR;$xIMzd8LlHjJ^=&0~&dT^U<% zM}EdmERX_&I08~4CZz-*RPiNt)g)GCfB7?3m9YZRCL50tM8b`fU`B@;aZ!~=q6+R> z^Pfrz<_Xr*6Jue#k)4I}so02cfg`+9B?b9G0@8h(>WD@^1_Z6}V^xlXv}XZP5xm(% zkob@yG9wILDKMkQg`)Y$TL4H<7{Pp^x&j9dQ6GYZNpK1jNot}%!YL&7#+nb@DFsro zbRj|p4P*=7xDXmDGN=HWjp3}~Lk$R$*${jnAW74@h<1^knm|GVEEJd&DiomqLaA0H z&8I_j8acO&$@**W##`(KBWD{ec35 zC_`l+;7OsPsK^5Rwp=r?*yz9}okc23l6j`Xlc=zfTsJ}xb0L^Hq6NaNc&F!K0>(Bc-#Ivjffs}iTaRequwQIXN-M($Xo@nY#j}}tIj;; zn^=tXL5B3gnyeHh5!h-_wW*+BlH-c4)nCvjY>aUrEr0W$Phy-MD4bM zjYLKu1SBRRB2a@g86r+E5E-)|r9+~?t~{9?lUty4FSG_}lXT)uB*ji<>!Uew*acFK zC&Erw(lb*@PYBT*lvH8>5`?*pNWsE5xFQ8y;g$|U#-t1XEchtHS*|F6p#=WI$`H8L z553?I>mGaEDXd>U)~gyM!=sCSLwDpx(AH05Im4-YnCemYAlymjIVdC30$mLmI@lBn zHH4%r2mwcah(f#eR5l_7Inde0Q&}-9v~N2VV*3vJ*;83vASYf^fvhP?ssdGzszBFs zsHv_p60BW%TqjNwdU`awoq5hl{CYp zd}uo$ujX-@TEj^&x`cwXWk^e?7am$3Cy-ufC0NVJjltB(ya&y+!Jy7dtYz8rPG`f& zW>(pepeh&fx+_r(6zRf4zrZVGB4BLEL_kQO@#qqDS&a^Kto~3^PxOS=qI23>UZZ_O z-Ud&#c4w3i@wCDj9}dCDYzHe} z6ZMnyaHhzlj}b-|Jns6#)IwL zam8^g4NM0?C74TT1Tzg5CtiM$4KD5?r~N7Jw48|)`@OSR@pN)S1c<_1(t^cE_G?Cx z-6=^8S)D0((S2JS!onn2Sz6v2~!J` zFPMoK`9Vf0(I@6rVaNec0u+cpTrSD4lxVroA^374csU^}hj1-)gliR_Vw4epVo7)5 zFujI`4Oza7$|F1(4J+28%5fnGPBEnjps-lNL2r?*m(oNb3`?TKf}ABDAq|huYQBXJ z!D58`+=Nd}6+SL5@g{eWcyN_s`V=KuHyKg^R>E;WBuRiwnD2&;Wu<@|7SBZJ@Brp< z5YNsO1Y6(28Ppq8(>7yG8ZURYTA*h4nV4}%<(!2n0R)N21j0?vghF&V)4X(n^gCr? z&KkndKq-qX`_7ZuU^tSpFc~$nl(mf@DKtvP?2+b|S@z$KW!<9~(+o$-AAv;#bA?O> zblDFLj#<@@=wm8A6Z@E&)W0}A7CqGKvQzZ{y;!Z$qU>mxtyxQ?KNZy zYh)Sr@20@3bSV3(K-pIq95Ol3QnIfwm?e`AD^fZPA<(YWT^Ld&%n9!f`R(FS zBb#Wn3&t{yKyT;ZNxa|WVK}&8P?9}b#dSFcOEaWcz?=jO>0h=KLod{TrIBLlfU5GPX0WvP+rD+hV3Yy#Q zYJq%h?|MlR=?9T6Fvx^gynS12A(_@vw;HpZT0&4T-F<*%JbN%&TJue9;~`M>}Sqnb>LvP zh64B;u(>gL_LVSc$V!+Y>5vX>=yr`_Y0*iexJ}LnLzq>x%@-0~0K^ATGj?wRBYhUHW~Fp^*ECdn2m+*H)%=p%_Qb(^Fot>8HXTu3XE3P9jF6#y*o8S7m`#J z&?{A*5XEtW8*x&NmSZAz{Y;iuvjUnd=GdvrK&*?E(VR3bN;1b46Mj6CRYz$~LFTh8 zj>VHUP<)L9#Sdoic+6G}mM;ivqigEI$|XLW4`&pjVvD6!Ud0y^mkJO0-3)kfw_&?( zZ*3(`0|YsCe|V}ri^I7>l3L8ZCc=tO(whh(%{V>gQuE``BV)FZNi875N#7d7|qllwSS1PSh;L5BFIF*b|4+PQNNEskW4jO?bWkbW%`|R zn?;qm=01DJEOzvmj~UInJ!T^^g0m6yQqW;^1mzTk+NM%vvp>cJ6CTF%>CbHa7JKq+ zHk=^IqHU05(Kbl3Xd5Jwd8PwtHbI&XNIPb;K1Y3w<<=5~v%T<-RZ_;?JW`29#I#Hz zw)jFGx|U&XXd`TEa_cM}+haY^T1Kz;T2H~tKKgCegE#K6rv1IN-&O*GUGW?UQV@kF z;BL>*XnTfNr;Q?7U@FIkmU z4KBEyV=+MyUO6EmqeyN$_CSAn9Rm>#5H*naX4eUj&7{q)u>D03hMnL_v?#w6sRiB7GN!(6mv^Z^-ixz6lCi zzJ`iMpU5zkuOMC?GPKZhx$8ot!Wk;Z7h8;I2(*B^?da0kkK z<^9Ap&tvy0eZ>8o-N~IEt~ObiChDakET@!AlnQd>-(vI-)iaE;Vx-_p;G^1g;+Hs7 zsYt6-R7Q-RB1pziC2>NuP*f24F<}UxA{|S{z+Z_7Q>aV`%2O`7ik{T+UJ;`!sxzv` zXwJWBAXtn+DTV;uFr$mG-_XGFF>`XXjuO>(qcNa@NVoz`PnSfL3k(Z1WR@BA(9#?H zBMZWnkW)pogqTJ18NzEiy2*fBhAv^vM!Sn9#*9X-6fvc3F_7ZzVlV@17lSo6Cdz;n zOo!6tbf}VKOB6j#cLo>E5?MV>wK}U70MdaST44s4RuIv5C)^c_?jjGUq3`bK3tJGb ziudKR1#&+?B>-al6)70y5i}mdsfDAG3`DNz2~~-(99m&G{ia@{t2D-U-5KwFnJ9DO zQDZ!bmL7pJ(rIXDi6}KWkRW@C`4tpNcTpZeEXEOIlrDja>ZmBj(6~SXpn|*zNSJ5? z6kGvoS;{g<$w}xgv6)=J8COg=!(PvxB$Sw_girz<0HxFk3^Gv;&`J6quoi@S;ssQc zgCF3jCn6>kS5!xVhu)$>uxgK79b}?6EfQ7`Ldy_`#4^A~sKm)p;^xQ(MIi|j0szOP`q6*Hg4unfoO*M5_<-3cnPD8q?hEx&5 zRH|O1-)*EVDH|M5y-uKP+8vS}r|f8Z(k>(ZrSVY14bpBG1JZ66BT2g?k+jLL5)em1 z*-^z&Kw$!DXAswC9F)T4fF@Bmb)TeQR}6&2xqxB1%vV@yNk>p%P2YAda#GUVyMxj1Y};Lg$Q7PK;0< zXM})x-CD?H^Ht5KIGQ?P$H zI>Lww4>-ynHfRFp>S=={kJQr!GfGMtazWKcWD(S9djAV?!ftJ7{n zBv~RiB(D&v&G{soN~j03NLW=+swBQ>p+pZ*;!qi0t|SH}CrIMOaUlv)Jg6ilr9ncL zMgU@1iDPg_sZ?hmW7sCbBh3YfWsx3@Lk)<05;XbX>hcJ%0Ou-U0bQ1wU8+^GHB@2j zsg$FD<|b=vg|;v{NLvHkv2y!&w^2=E2iWEcL1dnX6^^9*295|pg%Uv-hvm)${2ZC4 zf-hL&$5m5NS5Kj%iE4@l6;?btEUDv=y*@`$5*fQh1Mqv0Q7_Ln!WaMrVRG+ zG9lzgf2Bd0j)_M1$i!$U9mA#5giH&V!m&Y-BTZ-~ye(vxWA=p|>lDJbsa0z&>%iF zo~lRH$RkbqXW(i-81<0T2J-@(TyFzZDgnWDkm*nmk|Y)S4i1N?FeMHG<1b-gz0v2u zz&N>+a!dpN6>R<+dEp{E*kuD4erKTn;bXZZD6l;|Xn2Xsyg)`+9#bj7#w+g3&nv@n zeqI%Z8(d8w@xuI^W{Ca3&3KK#JQ%N`5aoCkgzZt6ur5(fBMEVZO+@Z0Cf02LH)3=l zX3iF5YoP^lLl)-Is6G=vNwMj(;K>;TU`JTO~ zAr}ypRR=-NAkQJe5lgRJkw=0o4|`Yr2@8zMB_DNnyP#fkb(ih@5GQVI6-Bf%*UqoLVepeiGr3o6J> zKq{>Hh`>%hYBQlWmvm$X`35k~40!yEYw* zAr=)K1gFp+G#le&d~o5V6>hIWwYe>^I?N2R%*`ssJ}!!7t($0ohiJW(hFGCORVuVh zg@g*lRLIdY`w7Sp&CxThD&**yJt{hjP;88gl{%Hy2^6{^9EEBkSlY(i3`EEl9{n9nV4DfEiJ=<*QSu$h8%HO^wMb$ zyJTaW@~`1D6bI7wdFRV7E!}Kztz7Hzg>!^oD%>bIT-pV2P;trB+2FFQFg96AcwIy; zt?J|xcgyus%ok}D7VUKDf_yzFlTX8y52Klwz{r&of*?t3`Y!^q0c6|A%^B(8BDD_) z`=`c_u#nq!7k0?P?MUm1et;~%9?oJB3YUO_ zi_wOXa2|DqIv^{Dq#^-?{be2YSw!QC77Zv+`QX5JNIstq>msOn=yr`J?kqu)B$F+c z63A^)LPbkxIVA$b`6N1g;oM+q@@AMQkZSI~3*_Z82?Svmb5CLn%n83sSMezg$zS+u zC33L&m#|U?E|W)vu#zav$Rgm6HdoWxM{F3$RePz>QJf4cQN~b|jQ13ZfSX}FeO2?4 zA;4pbyN7B7MVc)d-2OxH(JNwpxfIRC=g_wOQz&0+fP@f@OPIcc z1_q)nf$cchssP^~2q|&xC5Pg(MNVS0fC{1Uh6qwtA!vcDLO5-;V8F_<`$ipjj_S{) z)4QNYATE~-(A;9_WHj;|TLe;$v{IgHixj%*AFNbjVELzrn5qC%5IA}oj9hZlA!2eE z!hu*X>=3n2dSr@3SyYl<2{B2b2-A$Ngt(zZ2+SoErx=Ft(pKFfQV~UH7gPk=Q6^nL zfZC57E~KO@=`~b3P64?c;1m~d1&Fm!&SgNKA{wtI+?6PQjHdp@;BG1ul@maZU~7zO zuD}YTkG6C8VD=K3i+udKYc(=DINx7vh~K8s8gS|I@}e< zUMPwDy69;IycQ+1OnxEy1^SS7OTTb8)z$88U9CdnN<mA!|3WW8mDSWw8Ym)=uA=KS6gQ?}KYAKI5Ee!~h$&32m_mUM zIC$Z0%`J|G!D#=TN84SqSUNaOp%_eW>11=Nx)=`lx1xH1{*5oFxlYsTlHBRgpu*2S2S+G4OIY zXP=~iAX=B>4zcJ`FhN1*iT-I6E8A$*OQYz8&LJ`-D#G%9*I+9BCM7!K_@e~kJNbF#c(P<-}me8=1 z2P-UK+$eIw>;#Iy$AlrH7#i6jWOz3SvsFgLB*lM~25D<0mV|LPZT^B_kSa2dENyc4 z6p-qK+fJJ5lJ$utas^-zFBsFn0y0f>z1!#-1*?Ed;7}wU5nb!y;UbBG0{F2a?KxTJqtcw$_fM&h&*Mb#_3kt{Q|DM0+M1R zqq?@8jF>KW`_UZ0C}}GpeVT@B(Lu`;ke8_d4ZtfGJuN05XR#E8-K4dYhbUNTxzOzd zS;;@y2y`EmN48O?&7r{n4*f*YGz7vRqLu@s8<5RvdDX~8Dow$P!Ou`5#bof%Hr3Ke zKyn#GE6t6tjY`8HfvezY4grO7c13EZjmIFUl>97`ICOce5M)#=$|r@eGjJ^a!2X{U z%76qeLc^&=B{X0N(yIudi%F_U3)RFXC_(}nn}yGFxRh$ zLH^OMAOSYZ1^N@&EC^0Q$@59sVuzQ=osTx+VyPNHPZ7p%FkYx|Hb{kpz)~dyU}?JB zA;`>B@__74#}M+HEAMKFkyP+w_d^WzYr%|wxmk*CcQ<{4)z-NfCsHt1CW?Dvu;dtX zZ6ZurVp&p&d{FYvD9u*l1|-PD+?^QHdG+zwcJ?cD&`O%q#Nkz-QEL-p^v4KIa&@v& znp~JgbYgSSFi4w=sIbuR9VqO#Dm^c(t7d<+Wy><64yJM>F%3 zEeRtM9>zKn`@3mdDo9JK)_$N(1w*t9722;teidp}Avyw!sa~FIlA{Y_+tVqmPvy;$ zXvs@5MSZYB;dN0M$pdKeh+*eSzJ2FXRv5)$Fw%J04geq7vS26JvOcn9eZ|r^Hu4=P zvQxxq4wP(#S330ZFIXOJc?TZeCk3+KsYNP*xR@5Wl`*d!u-Krfooz7T#Kh7qCw*7( z4JJZ!_hG=pOO}0S`!M2j9USh}uWBB6GxM_k*C^5hIy~+GKs@=(H^E*@BMUyHali&Mc#p8|p&cVwUzVYN z5VL!&ViV8*Y&F=(l*$bId~lf#<(7YNM<{3={Gyp(%S>{6QO8Vh5z9>fl$qd*f>*Nx zByP`N#q#?lLB~c97RGf{q>{Ym3gTH)PCL8?8D~Z+ZuaSR8w($O??y_3!z7_~n4r-H zEO|)|I+mh8P7R03r|oiH^FCRGQxW^h+t^Wk)-(Zz$^?6in`=}CeG3(%U z_Cj%c>!K%+vj!jRQUl0$5)J4l;XT1fEgpj0_MNx0?5KnSCL}hL8L7Gx2x4z-D}9Hd zUCGL%4W}{!rn7IjCWH_7#DXsYin76^r>n18=o@ZlxD7mxfr@z2rZNQ>kh#e#!Ab9o`m>HuK{krCB@4$w|=yQYs7y=xpw)H1eK!1Z! z2`AydCCv^-!>DwLf{X6*YL8AALmkacdE`Of7sYW9o$iaOLTfIktWK=DMM{{%mk+kcdxd7Mz86_xsssu&^Rs91ZT>?um)(UAT(5P@-K-C5 zY~FS^8^rjFG5d!#>;_x}b^Sf8!2j{9l;@S0{qQ~5_}gf|bPx7Z;|8gFn89+;CWoaY z=g#jv!8nq z>jK*zMltjcGa&hCE|IHN*gdsYE?PZF)oX#}Ty9dVkPtER!o_Qq@*6@~F5cm)II_}k zjFbdDHouYD^_h6(;B|mDGQy?DLDe+DfTM^eUOo}U7F4*+Z*Ghf1(`C!FUUSY1dlKSMsm|^+Jb^Lc;$FYFyR#Jr`Z6weJ?yV&q!Igl{ zA$bIoa)eOcXt(~Bh4XG8e0A98X?I`8h6Gv^sau?J`{4&!wyi&e3m5LE5@ESfmq3O+ zW?ESD)M8BF&C?=Sq6bT7F6gc8dl&7%`E8II%x&Aldv-Udie*eSo4;}0r#is5Vv z7FcKLTOIHaQk{Z7DGq-dLHKv!6F0kSYNGkAbxbf&@^|am82)_Bu6=+FMeyPWSb5I# z)F2YTBzOgGjHNr+);}PxNBjE&C>OK4JjnVQDj%Vk@}ccX>V*i|p0SQ~jT1h6=LE@K z$0~;b(J>8=8VcNy+0+7xfaF#;(FG=GoDdv~ZN(+tB(kva6z+8l;x8X$-GZ{DDt>Gm zuav&o#w+p@5|NwhN$%@U#_TU2WaIe@F?+;AtSjI5uHEtwE9bv|)&Afimd&@u>{}>) z+p8VozkY~aPMe4SX9J_X#P)ACu*vk^yMgt?W@5YSVOEOr%_ARX6)d`)It4!qMjvC9 zL9e42_B((EaxRS12##yVqbK*|4H=Sz!VTlPVnB^+9&N^o0wU=kK%EdT<*BRS-Vx-+(ZIVe<@}!+r*}3 zHEvEW7_;|pVs*vuQglMOT`Z~G0Q*;uuxh^HF?-`9tb5tRk5TU4nzk@Z>uy{XG`Gr{ z*4y7a!s^%xd-$WQs_gow51BRgC_D8;hrO^X2K`J%U1=w%qmvcwQHR5JiKw%yo7riE ztSg$?$$afr`|W0s1m52_v!dwQt!R@OXb>Fd+-S=^X5RGT4Y(!gh|`j?8`wcT4(!4w zKxP0M5)E!jV3BYNtu86Rxfh1HedA_!D&Mxw{%kY*4ZF@>_B%EnuUF!)JiFgxtO`-5 zKgKrl-H$Yf9%rX`3inBRcUJn)+bYwEX}z2+&_&#yAS6l;F>DQLHo zlA(FqaubLl9I+HMa&~ugetavYA}HGO9P7%iwqJRk<+5w+1JALV z9Je@0724|Uh5l|H_j@*@)BIuEzm4T~cMExS$9VwBL&b3t1I1O3-SzNi9AM9 zpx+y~NHP_M<#p5jIyp0t~rkqPc-WX1TdDF&0NLc9B#u2MPlF@=NS|zWtTvVLRAh zMj!KEX1&-bm$ZQd*s)b@6sh_HZP7-hFFA1uxtibE$4%jTUl z-2VKJY?gmD5gG~jX|J*!T{k!(^Z{IqX$S^QL4u@S;bxt(3-a;3nEfE$FnE2hu|;Gx zT=F_AEqagY4>PpIi*v;^FC#;Z>&U+Gb%?4)`}^0~R1nk8-(cqjw=_}bTi~jL39#c0 z_H!ikc#|0cnE)E4_*36xFCxC;Ep{gqZ0Bcp?8w`!f`ibndz&pyRLeu&QEK_Bca&Ow z{2ed|^1b&C`xzS3>s|I1yR|v>J+^>l?IAZhwWwWGn4=IVf)4Rq0FEE z&i?EBtUqu0&My5E8%7j*+Mn1;vEexstJBcU)TMQr{xaB1t7B#@dY;~wZGO(yce7sn zvFGfgce9N{Ka-VoHE-OAQgQlid6}v%G2f-P8!Lh)w$<*HQS0og57=-02V~%_nEmz# z=;u2z`T z`3^W05BfVB5>4i+E!0`FI9d9m5$PQWlziKVAX;i9xU{p92|OuXzA~HM0hlIw9${SX zPV#kxmsG7}n;uveJrm2nEL#!2F3^ky59fWTZ>h_{UdlFn<&5`ZXI) zEOPGGP01RzM82{|8nvL~)*zn}g$k{~`n_w;LqB^q9MFcOH-P z<%vJC!c>AXCxuIPT#>@N6FdGcg@=OWNgN1KDt?C57mv9cb5l|i&7WlP92Q+Ksnp#8 z20G|K^c@Hs2<9Df(Zu~o_%@iF3x@y^;TCk8(J@97edD=YmeR!FPU1St*I9E_*=Qn< zc2;f2qFH8^xhe^39m^d_jKpHp5eQD7#6VL&p3e7)?x8Tek|xG5loDp>pF_Z@P_=xregPtk=5DDo@HzR+~Ir zf-|umxV-ewo2GU@nXxaL9;qr|48JQ=J)iyE8r5Se`?Bed9-L+WgIZ;z!pc)U^Xhx8 zO+C-pq2d#_{IK5JgDytHSjjb=dq%KCM$>*ArWt!}*+Pu?NZ()6!Xe@%uhfmqOT zqrXv&#C^rMbDjYMR;SkN@Ae ziyT=0mO5hMPi%46NBi|V-y#2#dz^bvy^u$p1BOZKu4w7kSkwCtWuIHC`fohG^W6Vt zlkJyhv0Q7zfERvRu-5Yf1^dq0Gw{2DWaPqRl!U9dTHg&C!1grt9z20{ZO5zFWs8k> z`lo~*98hOw3pljZy0xZt!UkG!z}m?7lm@em)-YDfMTmSXe&7le`bG&vFMF*8OSc}O zLMv2AZSM7~-D~X{Qq5zpSw9Zx&u(w(KlEZ&ctF-iyAe@+m1HsK|6Ki_*0aN!*p8-) zhEHYT!(!k1=BSEkj@SQRA|z~NcgOqi<*Z9aT+4p2(%{lxVm&;vJGw&oo}Mg z!!?{t!|IH-okBI5#%pPL=b)%_5N8cIx}n-__~E2m2GkeJHT2qcG1x+kGW~#nhks%!q1W~^C zruD%w(V)1Fjd(-$8FOb*;{%wI$13Ze¬1xV|55)NI)&u1j;}1R?l`s&ei>)SvBlZ% zz5&SQhNbw6R<86$^c$=L#}=#9zAE*_D-%))y|~f>$Cwxjz?e7bSNXP(w>S~v$|Q&@ zCOQzE`jmsqt?z9&HLB{{*Syu#IRsX7cwO*4gNjhNr(G5tpi8-gv_)Vcyti{diKI)pAl1Pl)K1?pjj` zk(4W7T|8x?TSz}6&H7|Y0lxw;rgWv)ZYQ6Sn0V*O-Bqwx`Wv0s8nDKkLM5AKos!Nd z??tEfM?ZGv&>($d` zDe9$uMg;*Md`2ZoHvRmJ?%rtUlH;%l?I^zfy-~%!I|!%$UL$aNWQijo|M!i+Zix8N zbrbjA#D(v{(b>swejGfrX14!*kpFew=zkISg!+?WKRp3hlc&}4ci*sXn1)|^ZM2@5 zHmIlkaWuKZz)MRCyGaG%ZfsKt*lVtU5JU{-j1c$Xa}4&aHRha8qpt4!Ut{^-LvmG; zj5~z9Z+d0=SmxjQEqVcGqE%LZpz9_lAkkaL820y@IgX7-mnL`-y1Hr9xo5f4Zs)o( z9GanhCNSl#P3PU=Ov&SExB(J}fshtC|E2b6Ln&c( zF$2Rj+(0I;MT3J~Z=WYfd9ZP*(>&H=GfUFqC9FTqEIu|Ka_UH<)eYR}jg9~GYYQz} zM8G)Fm&&H#7PZw>Me-r{ z+=b(i^O1#>zU?&G)cdGVx?2f01xT#Q zN_t&1*r`2)S+bjej$Uxl(P=K6udFQ>Er1W8&*BS$(gy&B#KHPK)TcQ5YIo^%A!7co}p|5Sr0KDe2P;urt38zXk~CA+hn4yV)2zc^*@ zkkZ6UEB|9%K0>Xcx_sEoR?aUs#=E>0=e|2~Q?(Vothjxib{yr)ei?Iy&bg|Gwf(Z8 z;DnUu$L*oG$h{%fE&s*@CiBEKKyNGa@&fDRSQP+V7VF8cZM4?KCLr)tth#(F&0ZY; zZ3=!81{CUa0R|PQbYx=>pOeR$+t>}IZffktf8A(3(KtGu;Jf@3aNJp!p9yYz{PGVS zdQZWIDh17}qrk&g^n{1?y(@|e6rYeU7bkk-Sx{N#lF@*3(vnetV=dVbheJP^Lal;X z*Jz!3<;YGF(5baoMvhYKOrXirnnf3_lEZbSoTIr64NZllpi26UO!C0pd;HRp@dpwa z8`lz2MqKw(%WEl*J2Ne9(S zbmp36r>D7%`O5liS>xfgOVf2%UBcLRO@V9P_^*45Pn6Q)u&e{pbLRku&1*GX+Y{2` zPuD_vG<|jLZ~yCkw8~a=g9sbD!sORAHZ5DRpLIT^b8lI}l?rW436m*H`O*p)Pr9eptb7^|pGp7bkI>VxCn2AG&32M;Fh&u@}Ctz438ZR;OU;8Oy(I zvO3MW?55)XzsT#OZ@w$8hUVLamTsSrxivv_Qn%*YXF}G!WDhS}LtfOx6|T6+R7>Kh z>&iN09J#2#imfam*?RlRi9V-{)w*&T)~8OqlV^w&|if*D|)p5_cE=w;96Pe0Qa5g*W}V#{;iFVV2Z`*L6j-h3Zfg={Onn_jywm$A>A{`H0jCZc5jGNc@>5$R2 z{;VI?*CC_%0A)PwfsPsDv&`Ec7{fMNA3ZRYZL&r`SmD3XaaUM#AMC+*#H`;u*oXbj z`on|ez9lNfVoe2>@1Zc;U{ybKhHnRr3a_!EDqee|JUAT4LZFWDp34I6Cz zQJUPt^MK!FrMZ5yHSFPHw!$(u7Fcs0&U3z2J)G}+ZI8!&AOGsHQRVC7ayC(IWJls$vbk@63uMWWaq`Parp&XpmRhx&irOW|sdN6O{_G!3 zcWt`R71_+(hZ!y32G<%WA2H zF{okf=F;}HPo}gfNPD+q+S!yg18Ms2s@s3l__wdL`3^DuZ`GEulYyJA^{RdGyKdA^)D!X$lKSxF< zG>S4?KYMl~240GJD{t$7IRCX%9A|AEj%-RZJh*kqRJln(9kG7nz$h`}= zNj_2ufh6?=4s0EP+rOPB-^%*^wfu%gYxVC7@oIYf_kUoSyT}$ahncPP)~{=_wsq?e zX>I$%#eTU7;yxqI+VNINQ~e7Y7`xi~)oa44+5VE7mT8T8DYr?ybUGW-4lK1Av(&?m zO4EiNOBp+04cU1!+hOh4c_X{Sn*K^-^%1+stN8SJ`vQ`2bw;f&stLzO$ScBkGBA-rZs!mXd@OA5Njyg&&Tl*W(`! zBp75TYiJZ zY9r|GVTFv?Nxq0W&HRH`q~ixK8>&4r59P3id@_l>Z2jt!>+!en(-Qox`Baq3ZRaQ{ zAi?Xjx<6if?XG*@+;hcV>`X_QMW0@Rd|!N8z+P_h|Jln1s49_D?N>)ZtS5_*s?wiW zB@QwBB#fi?6z96N#E=>76)I)!o=f8?*k3AeELx>B{kUhThc@m1YVX~+1qAdc|6po| zRSKv){UOprv}KSEy)6G8++fP>M(WT@VgZ?Bjm&|2+dCXLcVRpwbI2RQ+QEprS>BRa z-U{jzKR%^maGI!NNgjNHNY4b z_#ICDw#KH4uW(@g%0}z-ud7(2b@|u7%#eH0$Qg@$H@|7LLf;g!tE}$dEQeQO(>L`_ zU7R-b!@C`KUZ}MoUzhTTxkM_<#O0tVfbi)4Geg%9M#vuv|Ls*)`L~~w`#AHv*#`r9 z{dX1Zl3MKy8xy{3E@o^bb3J}(_zg~++<Wlgekq%p5(*X+|Bl4AFCd$ zX0aOmgch9(Zy{L_ztsuCTF31IIQD_NmXNwrrs81F>pr&JHDB$;es-yY0Y;0xDCPP` z?DRBtFTb_X-k8P)xMp$lcWF4a-aIqJ-eXy-2q9)lo1nDkcW2X&0PpW);;2t!qn)m^ z>*#%}&Rz+~-4=qT<2TqBWwA4|cF`P~CY?CmoYq_0^{xG37MtNa&?LKQAK!(I>Lkx? zUD(g!dA{tz7Icwtu=X$CusZY|_66B2mir#Tr*q`8jx-UFb)@I8-xhb8gxp48qOIrP zg5SU9uv_}RAp!K)w54^pX2yIyIrv&~5Y6A1%eryA-_FI=8F+WeW9Me5UZ!E~Mmt7- z&|aFy?m4HGipEa{tBV5VRh$7c&0Muf_Fu6Iq$LjHrJ2iil6ABxRN*;`IpiitE_ z4ulOLQbyaYrK~1*H7Ow4nu*yO@eh0!w>8>Fm9cXRR1%bun}lCgi|BVY+P^Dfa1+># zg259=U~=qWIeYvN_2k-Lmb3Zgcge~CF@-!ObpIU04VWz`@SDVXBE(#KonSxbs~c^9 zm|X=eY+GU0tG{d#%n!FkUR^5CfPn}yePeIaU}0maYuMVYVRjSW{GRaCj;fvYPLVh&1#@bI6^1Qs3?`Z6vrN5ZqMo%`hViPB_|5nHc zvn6)8h?n+ia}xk{VDltH)(1UzdjV_uSgNU$+bywaR}L@cBq)cC~fXuj9t{n zZl|i^`q)F9j-Skgpvjtu3q#CzlalEK>w%;g`a#saD#k7><L;8+uIF38+OvL z3f^-eM<%J2M8}z@C#wQnU(x}^(wHRspON8 z(xodm9N3;tU?VOGHo|e(GrLOIzv;>^U|5dmhS#I^(r&!G$hGO9m$->*(Q_d2IrJBP zg!QRzybp4J-i=$h-{$h}yvJbIsEbK8ihPL$NaZLC5VJOsU;%#6osaAkoBgZe*i1(A zudDdD=pK14OYNbQ-yEkMO>MmT6#_aHw2E?Ij|u))@s$*i8!pusqBhB>=&|m7dvwv2 zU|Cm4DN27Z1w+~43XtH~BcpGjZfN+e0G|i1lNhtU_%4@!eTTiEnrBDX%B$yaA5X-T zpXL*EQ;7MrX6~jB+|MLS_tVf;qHBzkJP3B^6Q^1rgF1mL^IGUuJ{&j!L+DSdmaJgR z6RiMT)N^wW{xkNtJ*fwV^Km=YgP*SNr9MW&9$rgGSk`Lw%Iq6lMqmV;{XoQdne*)% zd-91zPbq-7Htrc>X7dZ&^g$oEYUJCV{MXUl1S-R;X&xpK!-qT93C7%;x9RFaG}PlJ zpqVv%r^K>WM@@{y*bvV)CjzJEfNEz_M+AS0v!15-xuE&3& zz-*N+hra{Eaz^3YZIeTW2|DbhKYsnz=DOZ|9%H|^H}>K4@v7*{e}&h&zINXlsGuvFAcqI%>Mk3%=C#T)NAAk9$U8As%@6QLYKiIMU{L7yh z@>2u&*MR)!z@HrQSA)3WT)wELO~zpUa_H{8jKLm^wP=Brdz*h4%oagrp8oN0^d`_7rxuyE$$J_{GwtB3LK_S8}QIC1dIGv`j7 zGi|Yb#wdP5<-v1FkW(+5S$~lo9?nPDwWD~NJ$WS0wM&NcbXiA9XQgM)q}aZ*o5gV6 z$`8|)`WZ7DX3m@2Rzovt_`NFC>k!cA*DsuNLBljEJ9Sb0f(z=KGiv!j&*{^r&7U@X zX2Zo|-i7rG&z(K*qU9l7o2+SC37!-@MR``r9*Qr;E(sBHx6S{bHtrAg-kasBFZnz#{ob8NWTH^iL4ZN(jG&u#pgc8{q*7;b&Di z8gKVgh!6=G&PM|XZEgc+BHSk-yb$4VyKujR@Q;9(*Z>MArjL;D%iAhA4iOaz08U1@ zr@di3KPrx-@$IraPuUZS?8G+>o7?7A)S1}EFA>g9NbicU(JoAF+~C(`Q%2%D2UKwN z1fDe*P1RQ8@#7JA2rw0n13EK)ZNN z+XXDlN9mzXty9jg7c=M1M`JInpWZNUp_p^Qq6RT{UV}KNUewQ--*9msZ3tRUWj??| zEkAzR;+bXf3JtUdshG!#8HMrsX@ts&c9qm3pac&!h~Wsd>q2W(r)+}d{cwm@b~yu`R%S# z`5z6ppb@->KiBGNAlsC-`s?O&YDyQPd>xNzP=d-Soq%0%jo zz~@Z7?*c(kE8Or(gbBL)4~>+Pn8Nh?Wb`lA!)_eS>vA9Hp~n>)2_vGdzV@4=`DhMd zS~`ZG!cXmKUoeJ`8+lt#McmUM4xmcsPOqN^W}7dCw5u!oqc~{|QfTf!s(z5Ye+)m0 z-=A&|9?J(r2fJx_Q!?p8_aB9E6gmEn!ldZje-s|(2Joiv@OI&&5FS5m(fOi5vUPiP z@44v3VDp`4^M#qw%QAFr7IHTM3~4VntYJimW8H8H!qh1@oQCkx3E@nHsTFSeI=~}s z;D+x-cr3zVT7!V;Q!nuSVX8h!n~jfekQi1GGWdd$B&1X2#MvLF z+8@m4MFaMN0~aAY0;~&lq$P$4jvMY5PjBwMfY))uZS9!=nArQd2p6K@4*~o7Mf`SN z7qa^|;Bto5K6_>ZM9g0^?1vlRg86Gk^WPhInTKDQYG*A*h7;5635)qTMt`9BH~=8B zB!<`r4*Uh)dFW|1-m87hFD>SVzxl@)?^hL_Z0g$CcoyKf0*{5~Zak0T*@oveJRjlN zkH=f5YngZ|@eIQ=9?vOwX5zUH&+~Y8^$z%}=ML1`(t;GZ;ZbQouLsKh511q5H4J&;VzL0gk zScaEhxR1GWS!imhQITU&c$lB5yG!8Y>DKBvq&8o%Ud3vNU7M|~rh*1`smbddWdVaF BH1z-g delta 177 zcmX@p%y_1maRZM9hk?1FrKypz#bgl+ML1`(t;GZ;ZUgfWLrW`DBP*cn4J&;VzL0gk zSXDrtcaddPWq?yDKBvq&8o%Ud3vNU7M|~rh*1`smbddWdVgu BIUxW5 diff --git a/bindings/typst/justfile b/bindings/typst/justfile index ef8c2e5..50995b5 100644 --- a/bindings/typst/justfile +++ b/bindings/typst/justfile @@ -1,3 +1,17 @@ +# List tasks available +default: + @just --list --list-prefix " - " + +# Build the project using cargo to create a WebAssembly module build: cargo build --release --target wasm32-unknown-unknown cp ../../target/wasm32-unknown-unknown/release/cnvx.wasm examples/cnvx.wasm + +# Compile an example file +compile file="examples/power.typ": build + typst compile {{ file }} + +# Clean the build artifacts +clean: + cargo clean + rm -rf examples/cnvx.wasm diff --git a/crates/cnvx-math/Cargo.toml b/crates/cnvx-math/Cargo.toml index f4f758b..7c67db5 100644 --- a/crates/cnvx-math/Cargo.toml +++ b/crates/cnvx-math/Cargo.toml @@ -10,6 +10,8 @@ keywords = { workspace = true } readme = { workspace = true } [dependencies] + +[target.'cfg(not(target_arch = "wasm32"))'.dependencies] cblas = { workspace = true } lapacke = { workspace = true } openblas-src = { workspace = true } # Links with openblas, which provides BLAS and LAPACK implementations diff --git a/crates/cnvx-math/build.rs b/crates/cnvx-math/build.rs index e6d2002..4427856 100644 --- a/crates/cnvx-math/build.rs +++ b/crates/cnvx-math/build.rs @@ -1,9 +1,21 @@ fn main() { + let target_arch = std::env::var("CARGO_CFG_TARGET_ARCH").unwrap(); + if cfg!(target_os = "windows") { println!("cargo:error=Windows is not currently supported."); std::process::exit(1); } + if target_arch == "wasm32" { + println!( + "cargo:warning=WebAssembly target detected, skipping lapack probe, using custom implementations." + ); + } else { + probe_lapack(); + } +} + +fn probe_lapack() { // Try lapacke first, fall back to lapack. if pkg_config::probe_library("lapacke").is_err() { pkg_config::probe_library("lapack") diff --git a/crates/cnvx-math/src/lib.rs b/crates/cnvx-math/src/lib.rs index 02af3bf..91cbeac 100644 --- a/crates/cnvx-math/src/lib.rs +++ b/crates/cnvx-math/src/lib.rs @@ -8,6 +8,7 @@ //! //! - [`matrix`]: Defines [`DenseMatrix`] and the [`Matrix`] trait for linear algebra operations. +#[cfg(not(target_arch = "wasm32"))] extern crate openblas_src; pub mod matrix; diff --git a/crates/cnvx-math/src/matrix/dense.rs b/crates/cnvx-math/src/matrix/dense.rs index 61a8620..3180e82 100644 --- a/crates/cnvx-math/src/matrix/dense.rs +++ b/crates/cnvx-math/src/matrix/dense.rs @@ -1,4 +1,6 @@ -use crate::matrix::{Matrix, solvers}; +use crate::matrix::Matrix; + +#[cfg(not(target_arch = "wasm32"))] use cblas::{Layout, Transpose, dgemm, dgemv}; /// A dense matrix stored in row-major order using a flat 1D vector. @@ -17,6 +19,7 @@ impl DenseMatrix { } /// Expose the underlying data slice to internal LAPACK solvers + #[cfg(not(target_arch = "wasm32"))] pub(crate) fn data(&self) -> &[f64] { &self.data } @@ -71,12 +74,14 @@ impl Matrix for DenseMatrix { return Err("Incompatible dimensions for matrix multiplication".to_string()); } - let m = self.rows as i32; - let k = self.cols as i32; - let n = other.cols as i32; let mut result = Self::new(self.rows, other.cols); + #[cfg(not(target_arch = "wasm32"))] unsafe { + let m = self.rows as i32; + let k = self.cols as i32; + let n = other.cols as i32; + dgemm( Layout::RowMajor, Transpose::None, @@ -94,9 +99,24 @@ impl Matrix for DenseMatrix { n, // ldc ); } + + #[cfg(target_arch = "wasm32")] + { + for i in 0..self.rows { + for j in 0..other.cols { + let mut sum = 0.0; + for k in 0..self.cols { + sum += self.get(i, k) * other.get(k, j); + } + result.set(i, j, sum); + } + } + } + Ok(result) } + #[cfg(not(target_arch = "wasm32"))] fn mul_vec(&self, rhs: &[f64]) -> Result, String> { if self.cols != rhs.len() { return Err("Matrix columns must match vector length".to_string()); @@ -125,6 +145,23 @@ impl Matrix for DenseMatrix { Ok(result) } + #[cfg(target_arch = "wasm32")] + 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 } @@ -254,6 +291,82 @@ impl Matrix for DenseMatrix { } fn mldivide(&self, rhs: &[f64]) -> Result, String> { - solvers::mldivide_dense(self, rhs) + // solvers::mldivide_dense(self, rhs) + #[cfg(not(target_arch = "wasm32"))] + { + use crate::matrix::solvers; + #[allow(clippy::needless_return)] + return solvers::mldivide_dense(self, rhs); + } + + #[cfg(target_arch = "wasm32")] + { + 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(); + let b = rhs.to_vec(); + + // 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); + } + p.swap(i, imax); + } + + for j in (i + 1)..n { + lu[j * n + i] /= lu[i * n + i]; + for k in (i + 1)..n { + lu[j * n + k] -= lu[j * n + i] * lu[i * n + k]; + } + } + } + + // Apply permutation to RHS + let mut x = vec![0.0; n]; + for i in 0..n { + x[i] = b[p[i]]; + } + + // Forward substitution (Ly = Pb) + for i in 0..n { + for k in 0..i { + x[i] -= lu[i * n + k] * x[k]; + } + } + + // Backward substitution (Ux = y) + 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]; + } + + Ok(x) + } } } diff --git a/crates/cnvx-math/src/matrix/mod.rs b/crates/cnvx-math/src/matrix/mod.rs index daa0e2f..78f2313 100644 --- a/crates/cnvx-math/src/matrix/mod.rs +++ b/crates/cnvx-math/src/matrix/mod.rs @@ -1,7 +1,9 @@ mod dense; -mod solvers; mod sparse; +#[cfg(not(target_arch = "wasm32"))] +mod solvers; // Only include LAPACK-based solvers on non-WASM targets + pub use dense::DenseMatrix; pub use sparse::SparseMatrix; From 0d950c2673db28aebe500fe54ca5d4191fed2d34 Mon Sep 17 00:00:00 2001 From: Chris Date: Fri, 19 Jun 2026 15:22:13 +1200 Subject: [PATCH 10/10] ci: fix linting --- crates/cnvx-math/src/matrix/mod.rs | 2 +- crates/cnvx-math/src/matrix/solvers/mod.rs | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/crates/cnvx-math/src/matrix/mod.rs b/crates/cnvx-math/src/matrix/mod.rs index 78f2313..c46597d 100644 --- a/crates/cnvx-math/src/matrix/mod.rs +++ b/crates/cnvx-math/src/matrix/mod.rs @@ -116,7 +116,7 @@ pub trait Matrix: Clone { fn diagonal(&self) -> Vec; /// Solves Ax = b. - /// Returns a dynamically allocated Vec to support both square + /// Returns a dynamically allocated `Vec` to support both square /// and rectangular (least-squares) solutions. fn mldivide(&self, rhs: &[f64]) -> Result, String>; } diff --git a/crates/cnvx-math/src/matrix/solvers/mod.rs b/crates/cnvx-math/src/matrix/solvers/mod.rs index d5567e8..93e92c4 100644 --- a/crates/cnvx-math/src/matrix/solvers/mod.rs +++ b/crates/cnvx-math/src/matrix/solvers/mod.rs @@ -6,7 +6,7 @@ pub mod triangular; use crate::matrix::{DenseMatrix, Matrix}; /// MATLAB-style mldivide dispatcher for full matrices -/// See https://mathworks.com/help/matlab/ref/double.mldivide.html +/// See pub fn mldivide_dense(a: &DenseMatrix, b: &[f64]) -> Result, String> { let m = a.rows(); let n = a.cols();