From dc7ba40d89b680ec6c5aa5a3f3750c6bb26d5b91 Mon Sep 17 00:00:00 2001 From: Timo Betcke Date: Mon, 6 Nov 2023 05:33:05 -0600 Subject: [PATCH] Multiple interface updates for LU decomposition --- algorithms/src/lib.rs | 3 - common/src/traits/operations.rs | 6 +- dense/Cargo.toml | 17 ++- dense/src/array.rs | 7 ++ dense/src/array/mult_into.rs | 46 ++++--- dense/src/lib.rs | 3 + dense/src/linalg/lu.rs | 207 ++++++++++++++++++++++++++++---- 7 files changed, 233 insertions(+), 56 deletions(-) diff --git a/algorithms/src/lib.rs b/algorithms/src/lib.rs index f29756ec..ea6de332 100644 --- a/algorithms/src/lib.rs +++ b/algorithms/src/lib.rs @@ -1,9 +1,6 @@ //! Collection of Linear Solver Algorithms and Interfaces #![cfg_attr(feature = "strict", deny(warnings))] -extern crate rlst_blis_src; -extern crate rlst_netlib_lapack_src; - //pub mod dense; // pub mod iterative_solvers; // pub mod lapack; diff --git a/common/src/traits/operations.rs b/common/src/traits/operations.rs index 13849026..4b607574 100644 --- a/common/src/traits/operations.rs +++ b/common/src/traits/operations.rs @@ -133,17 +133,17 @@ pub trait PermuteRows { /// Multiply First * Second and sum into Self pub trait MultInto { type Item: Scalar; - fn mult_into(&mut self, alpha: Self::Item, arr_a: First, arr_b: Second, beta: Self::Item); + fn mult_into(self, alpha: Self::Item, arr_a: First, arr_b: Second, beta: Self::Item) -> Self; } /// Multiply First * Second and sum into Self. Allow to resize Self if necessary pub trait MultIntoResize { type Item: Scalar; fn mult_into_resize( - &mut self, + self, alpha: Self::Item, arr_a: First, arr_b: Second, beta: Self::Item, - ); + ) -> Self; } diff --git a/dense/Cargo.toml b/dense/Cargo.toml index 588b3f75..0f3d6e6a 100644 --- a/dense/Cargo.toml +++ b/dense/Cargo.toml @@ -23,19 +23,18 @@ cauchy = "0.4" rand = "0.8" itertools = "0.10" rand_distr = "0.4" -rlst-blis = { path = "../blis"} -approx = { version = "0.5", features=["num-complex"] } -rlst-operator = {path = "../operator"} -rlst-common = {path = "../common"} -rlst-lapack = {path = "../lapack"} +rlst-blis = { path = "../blis" } +approx = { version = "0.5", features = ["num-complex"] } +rlst-operator = { path = "../operator" } +rlst-common = { path = "../common" } +rlst-lapack = { path = "../lapack" } paste = "1" rand_chacha = "0.3" +rlst-blis-src = { path = "../blis-src" } +rlst-netlib-lapack-src = { path = "../netlib-lapack-src" } [dev-dependencies] criterion = { version = "0.3", features = ["html_reports"] } [package.metadata.docs.rs] -rustdoc-args = [ "--html-in-header", "katex-header.html" ] - - - +rustdoc-args = ["--html-in-header", "katex-header.html"] diff --git a/dense/src/array.rs b/dense/src/array.rs index 05833d24..7eae1961 100644 --- a/dense/src/array.rs +++ b/dense/src/array.rs @@ -239,3 +239,10 @@ impl< self.0.resize_in_place(shape) } } + +/// Create an empty array of given type and dimension. +pub fn empty_array() -> DynamicArray { + let shape = [0; NDIM]; + let container = VectorContainer::new(0); + Array::new(BaseArray::new(container, shape)) +} diff --git a/dense/src/array/mult_into.rs b/dense/src/array/mult_into.rs index eca599be..d4fab1d5 100644 --- a/dense/src/array/mult_into.rs +++ b/dense/src/array/mult_into.rs @@ -29,15 +29,18 @@ impl< type Item = Item; fn mult_into( - &mut self, + mut self, alpha: Item, arr_a: Array, arr_b: Array, beta: Item, - ) { + ) -> Self { let transa = TransMode::NoTrans; let transb = TransMode::NoTrans; - crate::matrix_multiply::matrix_multiply(transa, transb, alpha, &arr_a, &arr_b, beta, self) + crate::matrix_multiply::matrix_multiply( + transa, transb, alpha, &arr_a, &arr_b, beta, &mut self, + ); + self } } @@ -64,12 +67,12 @@ impl< type Item = Item; fn mult_into( - &mut self, + mut self, alpha: Item, arr_a: Array, arr_b: Array, beta: Item, - ) { + ) -> Self { let transa = TransMode::NoTrans; let transb = TransMode::NoTrans; @@ -84,7 +87,8 @@ impl< &arr_with_padded_dim, beta, &mut self_with_padded_dim, - ) + ); + self } } @@ -111,12 +115,12 @@ impl< type Item = Item; fn mult_into( - &mut self, + mut self, alpha: Item, arr_a: Array, arr_b: Array, beta: Item, - ) { + ) -> Self { let transa = TransMode::NoTrans; let transb = TransMode::NoTrans; @@ -131,7 +135,8 @@ impl< &arr_b, beta, &mut self_with_padded_dim, - ) + ); + self } } @@ -161,12 +166,12 @@ impl< type Item = Item; fn mult_into_resize( - &mut self, + mut self, alpha: Item, arr_a: Array, arr_b: Array, beta: Item, - ) { + ) -> Self { let transa = TransMode::NoTrans; let transb = TransMode::NoTrans; @@ -175,7 +180,10 @@ impl< self.resize_in_place(expected_shape); } - crate::matrix_multiply::matrix_multiply(transa, transb, alpha, &arr_a, &arr_b, beta, self) + crate::matrix_multiply::matrix_multiply( + transa, transb, alpha, &arr_a, &arr_b, beta, &mut self, + ); + self } } @@ -203,12 +211,12 @@ impl< type Item = Item; fn mult_into_resize( - &mut self, + mut self, alpha: Item, arr_a: Array, arr_b: Array, beta: Item, - ) { + ) -> Self { let transa = TransMode::NoTrans; let transb = TransMode::NoTrans; @@ -229,7 +237,8 @@ impl< &arr_with_padded_dim, beta, &mut self_with_padded_dim, - ) + ); + self } } @@ -257,12 +266,12 @@ impl< type Item = Item; fn mult_into_resize( - &mut self, + mut self, alpha: Item, arr_a: Array, arr_b: Array, beta: Item, - ) { + ) -> Self { let transa = TransMode::NoTrans; let transb = TransMode::NoTrans; @@ -283,6 +292,7 @@ impl< &arr_b, beta, &mut self_with_padded_dim, - ) + ); + self } } diff --git a/dense/src/lib.rs b/dense/src/lib.rs index 86234df4..6251a118 100644 --- a/dense/src/lib.rs +++ b/dense/src/lib.rs @@ -33,6 +33,9 @@ //! - [Examples](crate::examples) #![cfg_attr(feature = "strict", deny(warnings))] +extern crate rlst_blis_src; +extern crate rlst_netlib_lapack_src; + pub mod base_array; pub mod data_container; pub mod linalg; diff --git a/dense/src/linalg/lu.rs b/dense/src/linalg/lu.rs index 1b2df4ea..41f5e2c8 100644 --- a/dense/src/linalg/lu.rs +++ b/dense/src/linalg/lu.rs @@ -46,10 +46,6 @@ impl< } } - pub fn ipiv(&self) -> &[i32] { - self.ipiv.as_slice() - } - pub fn solve_into< ArrayImplMut: RawAccessMut + UnsafeRandomAccessByValue<2, Item = Item> @@ -92,6 +88,24 @@ impl< } } + pub fn get_l_resize< + ArrayImplMut: UnsafeRandomAccessByValue<2, Item = Item> + + Shape<2> + + UnsafeRandomAccessMut<2, Item = Item> + + UnsafeRandomAccessByRef<2, Item = Item> + + ResizeInPlace<2>, + >( + &self, + mut arr: Array, + ) { + let m = self.arr.shape()[0]; + let n = self.arr.shape()[1]; + let k = std::cmp::min(m, n); + + arr.resize_in_place([m, k]); + self.get_l(arr); + } + pub fn get_l< ArrayImplMut: UnsafeRandomAccessByValue<2, Item = Item> + Shape<2> @@ -100,7 +114,7 @@ impl< >( &self, mut arr: Array, - ) -> Array { + ) { let m = self.arr.shape()[0]; let n = self.arr.shape()[1]; let k = std::cmp::min(m, n); @@ -124,7 +138,24 @@ impl< } } } - arr + } + + pub fn get_r_resize< + ArrayImplMut: UnsafeRandomAccessByValue<2, Item = Item> + + Shape<2> + + UnsafeRandomAccessMut<2, Item = Item> + + UnsafeRandomAccessByRef<2, Item = Item> + + ResizeInPlace<2>, + >( + &self, + mut arr: Array, + ) { + let m = self.arr.shape()[0]; + let n = self.arr.shape()[1]; + let k = std::cmp::min(m, n); + + arr.resize_in_place([k, n]); + self.get_r(arr); } pub fn get_r< @@ -135,7 +166,7 @@ impl< >( &self, mut arr: Array, - ) -> Array { + ) { let m = self.arr.shape()[0]; let n = self.arr.shape()[1]; let k = std::cmp::min(m, n); @@ -151,11 +182,41 @@ impl< arr.set_zero(); for col in 0..n { - for row in 0..=col { + for row in 0..=std::cmp::min(col, k - 1) { arr[[row, col]] = self.arr.get_value([row, col]).unwrap(); } } - arr + } + + pub fn get_p_resize< + ArrayImplMut: UnsafeRandomAccessByValue<2, Item = Item> + + Shape<2> + + UnsafeRandomAccessMut<2, Item = Item> + + UnsafeRandomAccessByRef<2, Item = Item> + + ResizeInPlace<2>, + >( + &self, + mut arr: Array, + ) { + let m = self.arr.shape()[0]; + + arr.resize_in_place([m, m]); + self.get_p(arr); + } + + fn get_perm(&self) -> Vec { + let m = self.arr.shape()[0]; + // let n = self.arr.shape()[1]; + // let k = std::cmp::min(m, n); + let ipiv: Vec = self.ipiv.iter().map(|&elem| (elem as usize) - 1).collect(); + + let mut perm = (0..m).collect::>(); + + for (index, &elem) in ipiv.iter().enumerate() { + perm.swap(index, elem); + } + + perm } pub fn get_p< @@ -166,7 +227,7 @@ impl< >( &self, mut arr: Array, - ) -> Array { + ) { let m = self.arr.shape()[0]; assert_eq!( arr.shape(), @@ -178,11 +239,12 @@ impl< arr.shape()[1] ); + let perm = self.get_perm(); + arr.set_zero(); for col in 0..m { - arr[[self.ipiv[col] as usize, col]] = ::one(); + arr[[perm[col], col]] = ::one(); } - arr } } @@ -202,30 +264,129 @@ impl< #[cfg(test)] mod test { + use rlst_common::{assert_array_relative_eq, tools::PrettyPrint}; + use crate::rlst_dynamic_array2; use super::*; + use crate::array::empty_array; #[test] fn test_lu_thick() { - let dim = [3, 5]; - let l_dim = [3, 3]; - let r_dim = [3, 5]; - let p_dim = [3, 3]; + let dim = [8, 20]; + let mut arr = rlst_dynamic_array2!(f64, dim); + + arr.fill_from_seed_normally_distributed(0); + let mut arr2 = rlst_dynamic_array2!(f64, dim); + arr2.fill_from(arr.view()); + + let lu = arr2.into_lu().unwrap(); + + for (index, elem) in lu.get_perm().iter().enumerate() { + println!("Ipiv[{}] = {}", index, elem) + } + + let mut l_mat = empty_array::(); + let mut r_mat = empty_array::(); + let mut p_mat = empty_array::(); + + lu.get_l_resize(l_mat.view_mut()); + lu.get_r_resize(r_mat.view_mut()); + lu.get_p_resize(p_mat.view_mut()); + + let res = crate::array::empty_array::(); + + let res = res.mult_into_resize( + 1.0, + empty_array().mult_into_resize(1.0, p_mat, l_mat, 0.0), + r_mat, + 0.0, + ); + + assert_array_relative_eq!(res, arr, 1E-12) + } + + #[test] + fn test_lu_square() { + let dim = [12, 12]; + // let l_dim = [3, 3]; + // let r_dim = [3, 5]; + // let p_dim = [3, 3]; let mut arr = rlst_dynamic_array2!(f64, dim); arr.fill_from_seed_normally_distributed(0); let mut arr2 = rlst_dynamic_array2!(f64, dim); - arr2.fill_from(arr); + arr2.fill_from(arr.view()); let lu = arr2.into_lu().unwrap(); - let mut l_mat = rlst_dynamic_array2!(f64, l_dim); - let mut r_mat = rlst_dynamic_array2!(f64, r_dim); - let mut p_mat = rlst_dynamic_array2!(f64, p_dim); + for (index, elem) in lu.get_perm().iter().enumerate() { + println!("Ipiv[{}] = {}", index, elem) + } + + let mut l_mat = empty_array::(); + let mut r_mat = empty_array::(); + let mut p_mat = empty_array::(); + + lu.get_l_resize(l_mat.view_mut()); + lu.get_r_resize(r_mat.view_mut()); + lu.get_p_resize(p_mat.view_mut()); + + l_mat.pretty_print(); + r_mat.pretty_print(); + p_mat.pretty_print(); + + let res = crate::array::empty_array::(); + + let res = res.mult_into_resize( + 1.0, + empty_array().mult_into_resize(1.0, p_mat, l_mat, 0.0), + r_mat, + 0.0, + ); + + assert_array_relative_eq!(res, arr, 1E-12) + } + + #[test] + fn test_lu_thin() { + let dim = [12, 8]; + // let l_dim = [3, 3]; + // let r_dim = [3, 5]; + // let p_dim = [3, 3]; + let mut arr = rlst_dynamic_array2!(f64, dim); + + arr.fill_from_seed_normally_distributed(0); + let mut arr2 = rlst_dynamic_array2!(f64, dim); + arr2.fill_from(arr.view()); + + let lu = arr2.into_lu().unwrap(); + + for (index, elem) in lu.get_perm().iter().enumerate() { + println!("Ipiv[{}] = {}", index, elem) + } + + let mut l_mat = empty_array::(); + let mut r_mat = empty_array::(); + let mut p_mat = empty_array::(); + + lu.get_l_resize(l_mat.view_mut()); + lu.get_r_resize(r_mat.view_mut()); + lu.get_p_resize(p_mat.view_mut()); + + l_mat.pretty_print(); + r_mat.pretty_print(); + p_mat.pretty_print(); + + let res = crate::array::empty_array::(); + + let res = res.mult_into_resize( + 1.0, + empty_array().mult_into_resize(1.0, p_mat, l_mat, 0.0), + r_mat, + 0.0, + ); - lu.get_l(l_mat.view_mut()); - lu.get_r(r_mat.view_mut()); - lu.get_p(p_mat.view_mut()); + assert_array_relative_eq!(res, arr, 1E-12) } }