From df26fee2d83c1e12c045da1bfc5327d47251f4b1 Mon Sep 17 00:00:00 2001 From: Timo Betcke Date: Thu, 14 Mar 2024 17:34:56 +0000 Subject: [PATCH] Mac accelerate (#64) --- .github/workflows/run-tests.yml | 3 + .github/workflows/run-weekly-tests.yml | 3 + .gitmodules | 6 - Cargo.toml | 4 - blis-src/Cargo.toml | 21 - blis-src/blis | 1 - blis-src/build.rs | 50 -- blis-src/src/lib.rs | 1 - blis/Cargo.toml | 20 - blis/build.rs | 29 -- blis/src/interface.rs | 4 - blis/src/interface/gemm.rs | 128 ----- blis/src/interface/types.rs | 16 - blis/src/lib.rs | 5 - blis/src/raw.rs | 7 - blis/src/wrapper.h | 1 - dense/Cargo.toml | 11 +- dense/src/array/empty_axis.rs | 37 -- dense/src/array/iterators.rs | 114 +---- dense/src/array/slice.rs | 75 --- dense/src/gemm.rs | 146 +++++- dense/src/gemm/blis.rs | 102 ---- dense/src/layout.rs | 21 - dense/src/lib.rs | 3 - dense/src/linalg/inverse.rs | 47 +- dense/src/linalg/lu.rs | 141 ------ dense/src/linalg/pseudo_inverse.rs | 77 --- dense/src/linalg/qr.rs | 109 ----- dense/src/linalg/svd.rs | 136 ------ dense/src/matrix_multiply.rs | 134 ------ dense/src/traits.rs | 9 + gfortran-src/Cargo.toml | 18 - gfortran-src/build.rs | 6 - gfortran-src/src/lib.rs | 1 - lapack/Cargo.toml | 41 -- lapack/src/geqp3.rs | 98 ---- lapack/src/getrf.rs | 32 -- lapack/src/getrs.rs | 58 --- lapack/src/lib.rs | 20 - lapack/src/ormqr.rs | 253 ---------- lapack/src/unmqr.rs | 50 -- netlib-lapack-src/Cargo.toml | 22 - netlib-lapack-src/build.rs | 18 - netlib-lapack-src/lapack | 1 - netlib-lapack-src/src/lib.rs | 3 - operator/Cargo.toml | 1 - operator/src/interface.rs | 5 + operator/src/interface/array_vector_space.rs | 7 - .../src/interface/dense_matrix_operator.rs | 30 -- .../src/operations/conjugate_gradients.rs | 45 -- .../src/operations/modified_gram_schmidt.rs | 66 --- operator/src/operator.rs | 43 -- rlst/Cargo.toml | 25 +- rlst/examples/array_operations.rs | 2 +- rlst/src/dense.rs | 7 - rlst/src/lib.rs | 12 +- rlst/src/prelude.rs | 60 +++ {blis/src/interface => rlst/src}/threading.rs | 15 +- rlst/tests/array_operations.rs | 342 ++++++++++++- rlst/tests/linalg.rs | 449 ++++++++++++++++++ rlst/tests/operator.rs | 137 ++++++ rlst/tests/sparse.rs | 251 ++++++++++ sparse/Cargo.toml | 8 +- sparse/src/index_layout.rs | 2 + .../index_layout/default_mpi_index_layout.rs | 29 -- .../default_serial_index_layout.rs | 21 - sparse/src/lib.rs | 11 +- sparse/src/sparse.rs | 4 + sparse/src/sparse/csc_mat.rs | 76 --- sparse/src/sparse/csr_mat.rs | 77 --- sparse/src/sparse/umfpack.rs | 112 +---- sparse/src/traits.rs | 2 - sparse/src/traits/indexable_matrix.rs | 149 ------ sparse/src/traits/indexable_vector.rs | 109 ----- sparse/src/traits/matrix_traits.rs | 71 --- suitesparse-src/Cargo.toml | 4 - suitesparse-src/build.rs | 26 +- umfpack/Cargo.toml | 3 +- umfpack/examples/umfpack_simple.rs | 61 --- umfpack/src/lib.rs | 2 - 80 files changed, 1485 insertions(+), 2861 deletions(-) delete mode 100644 blis-src/Cargo.toml delete mode 160000 blis-src/blis delete mode 100644 blis-src/build.rs delete mode 100644 blis-src/src/lib.rs delete mode 100644 blis/Cargo.toml delete mode 100644 blis/build.rs delete mode 100644 blis/src/interface.rs delete mode 100644 blis/src/interface/gemm.rs delete mode 100644 blis/src/interface/types.rs delete mode 100644 blis/src/lib.rs delete mode 100644 blis/src/raw.rs delete mode 100644 blis/src/wrapper.h delete mode 100644 dense/src/gemm/blis.rs delete mode 100644 gfortran-src/Cargo.toml delete mode 100644 gfortran-src/build.rs delete mode 100644 gfortran-src/src/lib.rs delete mode 100644 lapack/Cargo.toml delete mode 100644 lapack/src/geqp3.rs delete mode 100644 lapack/src/getrf.rs delete mode 100644 lapack/src/getrs.rs delete mode 100644 lapack/src/lib.rs delete mode 100644 lapack/src/ormqr.rs delete mode 100644 lapack/src/unmqr.rs delete mode 100644 netlib-lapack-src/Cargo.toml delete mode 100644 netlib-lapack-src/build.rs delete mode 160000 netlib-lapack-src/lapack delete mode 100644 netlib-lapack-src/src/lib.rs delete mode 100644 rlst/src/dense.rs create mode 100644 rlst/src/prelude.rs rename {blis/src/interface => rlst/src}/threading.rs (62%) create mode 100644 rlst/tests/linalg.rs create mode 100644 rlst/tests/operator.rs create mode 100644 rlst/tests/sparse.rs delete mode 100644 sparse/src/traits/indexable_matrix.rs delete mode 100644 sparse/src/traits/indexable_vector.rs delete mode 100644 sparse/src/traits/matrix_traits.rs delete mode 100644 umfpack/examples/umfpack_simple.rs diff --git a/.github/workflows/run-tests.yml b/.github/workflows/run-tests.yml index 40b6d6d2..ea8e9b4c 100644 --- a/.github/workflows/run-tests.yml +++ b/.github/workflows/run-tests.yml @@ -37,6 +37,9 @@ jobs: - name: Install Clang run: sudo apt-get install -y libclang-dev + - name: Install Blas + run: sudo apt-get install -y libopenblas-dev + - name: Install CMake run: sudo apt-get install -y cmake diff --git a/.github/workflows/run-weekly-tests.yml b/.github/workflows/run-weekly-tests.yml index 070b6ca7..3ba06a0c 100644 --- a/.github/workflows/run-weekly-tests.yml +++ b/.github/workflows/run-weekly-tests.yml @@ -33,6 +33,9 @@ jobs: - name: Install Clang run: sudo apt-get install -y libclang-dev + - name: Install Blas + run: sudo apt-get install -y libopenblas-dev + - name: Install CMake run: sudo apt-get install -y cmake diff --git a/.gitmodules b/.gitmodules index 03cd710c..dfec6993 100644 --- a/.gitmodules +++ b/.gitmodules @@ -2,10 +2,4 @@ path = suitesparse-src/suitesparse url = https://github.com/DrTimothyAldenDavis/SuiteSparse.git -[submodule "blis-src/blis"] - path = blis-src/blis - url = https://github.com/flame/blis.git -[submodule "netlib-lapack-src/lapack"] - path = netlib-lapack-src/lapack - url = https://github.com/Reference-LAPACK/lapack.git diff --git a/Cargo.toml b/Cargo.toml index d34a6c9d..8aa7f204 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -7,11 +7,7 @@ members = [ "sparse", "io", "suitesparse-src", - "blis-src", "umfpack", - "gfortran-src", - "netlib-lapack-src", - "blis", "rlst", "proc-macro", ] diff --git a/blis-src/Cargo.toml b/blis-src/Cargo.toml deleted file mode 100644 index ff0f5998..00000000 --- a/blis-src/Cargo.toml +++ /dev/null @@ -1,21 +0,0 @@ -[features] -strict = [] - -[package] -name = "rlst-blis-src" -version = "0.1.0" -edition = "2021" -links = "blis" - -# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html - - -[lib] -name = "rlst_blis_src" - -[build-dependencies] -bindgen = "0.65.*" -cmake = "0.1.*" -autotools = "0.2" -cc = "1.*" - diff --git a/blis-src/blis b/blis-src/blis deleted file mode 160000 index d6395548..00000000 --- a/blis-src/blis +++ /dev/null @@ -1 +0,0 @@ -Subproject commit d639554894b6252a86bd3164921bce6fbb9e3b5e diff --git a/blis-src/build.rs b/blis-src/build.rs deleted file mode 100644 index c1fc9f63..00000000 --- a/blis-src/build.rs +++ /dev/null @@ -1,50 +0,0 @@ -use std::env; -use std::fs; -use std::path::PathBuf; -use std::process::Command; - -fn main() { - // Build BLIS - - let dst = PathBuf::from(env::var("OUT_DIR").unwrap()); - let build = dst.join("build"); - let _ = fs::create_dir(&build); - - let configure = PathBuf::from("blis") - .join("configure") - .canonicalize() - .unwrap(); - let mut config_command = Command::new("sh"); - config_command - .args(["-c", "exec \"$0\" \"$@\""]) - .arg(configure); - config_command.arg(format!("--prefix={}", dst.display())); - config_command.arg("--enable-threading=pthreads"); - config_command.args(["--enable-cblas", "auto"]); - let status = match config_command.current_dir(&build).status() { - Ok(status) => status, - Err(e) => panic!("Could not execute configure command with error {}", e), - }; - if !status.success() { - panic!("Configure command failed with error {}", status); - } - - let make_flags = env::var_os("CARGO_MAKEFLAGS").unwrap(); - let mut build_command = Command::new("sh"); - build_command - .args(["-c", "exec \"$0\" \"$@\""]) - .args(["make", "install"]); - build_command.env("MAKEFLAGS", make_flags); - - let status = match build_command.current_dir(&build).status() { - Ok(status) => status, - Err(e) => panic!("Could not execute build command with error {}", e), - }; - if !status.success() { - panic!("Build command failed with error {}", status); - } - - println!("cargo:root={}", dst.display()); - println!("cargo:rustc-link-search={}", dst.join("lib").display()); - println!("cargo:rustc-link-lib=static=blis"); -} diff --git a/blis-src/src/lib.rs b/blis-src/src/lib.rs deleted file mode 100644 index f05efe3b..00000000 --- a/blis-src/src/lib.rs +++ /dev/null @@ -1 +0,0 @@ -#![cfg_attr(feature = "strict", deny(warnings))] diff --git a/blis/Cargo.toml b/blis/Cargo.toml deleted file mode 100644 index 2522b8bf..00000000 --- a/blis/Cargo.toml +++ /dev/null @@ -1,20 +0,0 @@ -[features] -strict = [] - -[package] -name = "rlst-blis" -version = "0.1.0" -edition = "2021" - -# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html - -[dependencies] -rlst-blis-src = { path = "../blis-src" } -num_cpus = "1" -num = "0.4" - -[lib] -name = "rlst_blis" - -[build-dependencies] -bindgen = "0.65.*" diff --git a/blis/build.rs b/blis/build.rs deleted file mode 100644 index b1cbd1a5..00000000 --- a/blis/build.rs +++ /dev/null @@ -1,29 +0,0 @@ -use std::path::PathBuf; - -fn main() { - let root = PathBuf::from(std::env::var("DEP_BLIS_ROOT").unwrap()); - - let bindings = bindgen::Builder::default() - // The input header we would like to generate - // bindings for. - .header("src/wrapper.h") - // Add an include path - .clang_arg(format!("-I{}", root.join("include/blis").display())) - .allowlist_function("bli.*") - .allowlist_type("BLIS.*") - .allowlist_var("BLIS.*") - // Tell cargo to invalidate the built crate whenever any of the - // included header files changed. - .parse_callbacks(Box::new(bindgen::CargoCallbacks)) - // Finish the builder and generate the bindings. - .generate() - // Unwrap the Result and panic on failure. - .expect("Unable to generate bindings"); - - // Write the bindings to the $OUT_DIR/bindings.rs file. - let out_path = PathBuf::from(std::env::var("OUT_DIR").unwrap()); - //println!("Out path: {}", out_path.to_str().unwrap()); - bindings - .write_to_file(out_path.join("bindings.rs")) - .expect("Couldn't write bindings!"); -} diff --git a/blis/src/interface.rs b/blis/src/interface.rs deleted file mode 100644 index 93dd3a34..00000000 --- a/blis/src/interface.rs +++ /dev/null @@ -1,4 +0,0 @@ -//! Public interface to Blis routines. - -// pub mod gemm; -pub mod threading; diff --git a/blis/src/interface/gemm.rs b/blis/src/interface/gemm.rs deleted file mode 100644 index b52f6f66..00000000 --- a/blis/src/interface/gemm.rs +++ /dev/null @@ -1,128 +0,0 @@ -//! Matrix product in Blis. -use super::types::TransMode; -use crate::interface::assert_data_size; -use crate::raw; -use cauchy::{c32, c64}; - -/// Transposition Mode. -#[derive(Clone, Copy, PartialEq)] -#[repr(u32)] - -pub enum TransMode { - /// Complex conjugate of matrix. - ConjNoTrans, - /// No modification of matrix. - NoTrans, - /// Transposition of matrix. - Trans, - /// Conjugate transpose of matrix. - ConjTrans, -} - -/// Performs the matrix multiplication `C = alpha * A * B + beta * C`. -/// # Arguments -/// -/// - `transa` - Transposition mode of A. -/// - `transb` - Transposition mode of B. -/// - `m` - Row dimension of C. -/// - `k` - Column dimension of C. -/// - `k` - Column dimension of A. Same as row dimension of B. -/// - `a` - Reference to data of A. -/// - `alpha` - Scalar `alpha parameter. -/// - `rsa` - Row stride of A. -/// - `csb` - Row stride of A. -/// - `b` - Reference to data of B. -/// - `rsb` - Row stride of B. -/// - `csb` - Column stride of B. -/// - `beta` - Scalar `beta` parameter. -/// - `c` - Reference to data of C. -/// - `rsc` - Row stride of C. -/// - `csc` - Column stride of C. -pub trait Gemm: Sized { - #[allow(clippy::too_many_arguments)] - fn gemm( - transa: TransMode, - transb: TransMode, - m: usize, - n: usize, - k: usize, - alpha: Self, - a: &[Self], - rsa: usize, - csa: usize, - b: &[Self], - rsb: usize, - csb: usize, - beta: Self, - c: &mut [Self], - rsc: usize, - csc: usize, - ); -} - -macro_rules! impl_gemm { - ($scalar:ty, $blis_gemm:ident, $blis_type:ty) => { - impl Gemm for $scalar { - fn gemm( - transa: TransMode, - transb: TransMode, - m: usize, - n: usize, - k: usize, - alpha: Self, - a: &[Self], - rsa: usize, - csa: usize, - b: &[Self], - rsb: usize, - csb: usize, - beta: Self, - c: &mut [Self], - rsc: usize, - csc: usize, - ) { - match transa { - TransMode::NoTrans => assert_data_size(a, [rsa, csa], [m, k]), - TransMode::ConjNoTrans => assert_data_size(a, [rsa, csa], [m, k]), - TransMode::Trans => assert_data_size(a, [rsa, csa], [k, m]), - TransMode::ConjTrans => assert_data_size(a, [rsa, csa], [k, m]), - } - - match transb { - TransMode::NoTrans => assert_data_size(b, [rsb, csb], [k, n]), - TransMode::ConjNoTrans => assert_data_size(b, [rsb, csb], [k, n]), - TransMode::Trans => assert_data_size(b, [rsb, csb], [n, k]), - TransMode::ConjTrans => assert_data_size(b, [rsb, csb], [n, k]), - } - - assert_data_size(c, [rsc, csc], [m, n]); - - unsafe { - raw::$blis_gemm( - transa as u32, - transb as u32, - m as i64, - n as i64, - k as i64, - &alpha as *const _ as *const $blis_type, - a.as_ptr() as *const _ as *const $blis_type, - rsa as i64, - csa as i64, - b.as_ptr() as *const _ as *const $blis_type, - rsb as i64, - csb as i64, - &beta as *const _ as *const $blis_type, - c.as_mut_ptr() as *mut _ as *mut $blis_type, - rsc as i64, - csc as i64, - ) - }; - } - } - }; -} - -impl_gemm!(f32, bli_sgemm, f32); -impl_gemm!(f64, bli_dgemm, f64); -impl_gemm!(c32, bli_cgemm, raw::scomplex); -impl_gemm!(c64, bli_zgemm, raw::dcomplex); diff --git a/blis/src/interface/types.rs b/blis/src/interface/types.rs deleted file mode 100644 index ec2dabdc..00000000 --- a/blis/src/interface/types.rs +++ /dev/null @@ -1,16 +0,0 @@ -// //! Interface to Blis types -// use crate::raw; - -// /// Transposition Mode. -// #[derive(Clone, Copy, PartialEq)] -// #[repr(u32)] -// pub enum TransMode { -// /// Complex conjugate of matrix. -// ConjNoTrans = raw::trans_t_BLIS_CONJ_NO_TRANSPOSE, -// /// No modification of matrix. -// NoTrans = raw::trans_t_BLIS_NO_TRANSPOSE, -// /// Transposition of matrix. -// Trans = raw::trans_t_BLIS_TRANSPOSE, -// /// Conjugate transpose of matrix. -// ConjTrans = raw::trans_t_BLIS_CONJ_TRANSPOSE, -// } diff --git a/blis/src/lib.rs b/blis/src/lib.rs deleted file mode 100644 index 0b8a2b97..00000000 --- a/blis/src/lib.rs +++ /dev/null @@ -1,5 +0,0 @@ -//! Interface to BLIS -#![cfg_attr(feature = "strict", deny(warnings))] - -pub mod interface; -pub mod raw; diff --git a/blis/src/raw.rs b/blis/src/raw.rs deleted file mode 100644 index ab5174be..00000000 --- a/blis/src/raw.rs +++ /dev/null @@ -1,7 +0,0 @@ -//! Raw interface to Blis routines. - -#![allow(non_upper_case_globals)] -#![allow(non_camel_case_types)] -#![allow(non_snake_case)] - -include!(concat!(env!("OUT_DIR"), "/bindings.rs")); diff --git a/blis/src/wrapper.h b/blis/src/wrapper.h deleted file mode 100644 index 17caa1d0..00000000 --- a/blis/src/wrapper.h +++ /dev/null @@ -1 +0,0 @@ -#include "blis.h" diff --git a/dense/Cargo.toml b/dense/Cargo.toml index dd970518..59546886 100644 --- a/dense/Cargo.toml +++ b/dense/Cargo.toml @@ -1,6 +1,3 @@ -[features] -strict = [] - [package] name = "rlst-dense" version = "0.0.1" @@ -22,15 +19,13 @@ num = {version = "0.4", features = ["serde", "rand"] } rand = "0.8" itertools = "0.12" rand_distr = "0.4" -rlst-blis = { path = "../blis" } approx = { version = "0.5", features = ["num-complex"] } -paste = "1" rand_chacha = "0.3" -rlst-blis-src = { path = "../blis-src" } -rlst-netlib-lapack-src = { path = "../netlib-lapack-src" } +paste = "1" lapack = "0.19.*" thiserror = "1.0.40" serde = "1.*" +blas = "0.22.*" [dev-dependencies] @@ -38,3 +33,5 @@ criterion = { version = "0.3", features = ["html_reports"] } [package.metadata.docs.rs] rustdoc-args = ["--html-in-header", "katex-header.html"] + +[features] diff --git a/dense/src/array/empty_axis.rs b/dense/src/array/empty_axis.rs index 1737b37e..679bf1cf 100644 --- a/dense/src/array/empty_axis.rs +++ b/dense/src/array/empty_axis.rs @@ -232,40 +232,3 @@ where orig } - -#[cfg(test)] -mod test { - use crate::traits::{Shape, Stride}; - - use crate::{array::empty_axis::AxisPosition, rlst_dynamic_array3}; - - #[test] - fn test_insert_axis_back() { - let shape = [3, 7, 6]; - let mut arr = rlst_dynamic_array3!(f64, shape); - arr.fill_from_seed_equally_distributed(0); - - let new_arr = arr.view().insert_empty_axis(AxisPosition::Back); - - assert_eq!(new_arr.shape(), [3, 7, 6, 1]); - - assert_eq!(new_arr[[1, 2, 5, 0]], arr[[1, 2, 5]]); - - assert_eq!(new_arr.stride(), [1, 3, 21, 126]) - } - - #[test] - fn test_insert_axis_front() { - let shape = [3, 7, 6]; - let mut arr = rlst_dynamic_array3!(f64, shape); - arr.fill_from_seed_equally_distributed(0); - - let new_arr = arr.view().insert_empty_axis(AxisPosition::Front); - - assert_eq!(new_arr.shape(), [1, 3, 7, 6]); - - assert_eq!(new_arr[[0, 1, 2, 5]], arr[[1, 2, 5]]); - - assert_eq!(new_arr.stride(), [1, 1, 3, 21]) - } -} diff --git a/dense/src/array/iterators.rs b/dense/src/array/iterators.rs index 60bf827a..8ffb893a 100644 --- a/dense/src/array/iterators.rs +++ b/dense/src/array/iterators.rs @@ -2,6 +2,7 @@ use crate::array::*; use crate::layout::convert_1d_nd_from_shape; +use crate::traits::AsMultiIndex; use crate::types::RlstScalar; use super::slice::ArraySlice; @@ -56,12 +57,6 @@ impl, const NDIM: usize> Iterator } } -/// Helper trait that returns from an enumeration iterator a new iterator -/// that converts the 1d index into a multi-index. -pub trait AsMultiIndex, const NDIM: usize> { - fn multi_index(self, shape: [usize; NDIM]) -> MultiIndexIterator; -} - impl, const NDIM: usize> AsMultiIndex for I { fn multi_index(self, shape: [usize; NDIM]) -> MultiIndexIterator { MultiIndexIterator:: { shape, iter: self } @@ -369,110 +364,3 @@ impl< } } } - -#[cfg(test)] -mod test { - - use approx::assert_relative_eq; - - use crate::traits::*; - - #[test] - fn test_iter() { - let mut arr = crate::rlst_dynamic_array3![f64, [1, 3, 2]]; - - for (index, item) in arr.iter_mut().enumerate() { - *item = index as f64; - } - - assert_eq!(arr[[0, 0, 0]], 0.0); - assert_eq!(arr[[0, 1, 0]], 1.0); - assert_eq!(arr[[0, 2, 0]], 2.0); - assert_eq!(arr[[0, 0, 1]], 3.0); - assert_eq!(arr[[0, 1, 1]], 4.0); - assert_eq!(arr[[0, 2, 1]], 5.0); - } - - #[test] - fn test_row_iter() { - let shape = [2, 3]; - let mut arr = crate::rlst_dynamic_array2![f64, shape]; - - arr.fill_from_seed_equally_distributed(0); - - let mut row_count = 0; - for (row_index, row) in arr.row_iter().enumerate() { - for col_index in 0..shape[1] { - assert_eq!(row[[col_index]], arr[[row_index, col_index]]); - } - row_count += 1; - } - assert_eq!(row_count, shape[0]); - } - - #[test] - fn test_row_iter_mut() { - let shape = [2, 3]; - let mut arr = crate::rlst_dynamic_array2![f64, shape]; - let mut arr2 = crate::rlst_dynamic_array2![f64, shape]; - - arr.fill_from_seed_equally_distributed(0); - arr2.fill_from(arr.view()); - - let mut row_count = 0; - for (row_index, mut row) in arr.row_iter_mut().enumerate() { - for col_index in 0..shape[1] { - row[[col_index]] *= 2.0; - assert_relative_eq!( - row[[col_index]], - 2.0 * arr2[[row_index, col_index]], - epsilon = 1E-14 - ); - } - row_count += 1; - } - assert_eq!(row_count, shape[0]); - } - - #[test] - fn test_col_iter() { - let shape = [2, 3]; - let mut arr = crate::rlst_dynamic_array2![f64, shape]; - - arr.fill_from_seed_equally_distributed(0); - - let mut col_count = 0; - for (col_index, col) in arr.col_iter().enumerate() { - for row_index in 0..shape[0] { - assert_eq!(col[[row_index]], arr[[row_index, col_index]]); - } - col_count += 1; - } - - assert_eq!(col_count, shape[1]); - } - - #[test] - fn test_col_iter_mut() { - let shape = [2, 3]; - let mut arr = crate::rlst_dynamic_array2![f64, shape]; - let mut arr2 = crate::rlst_dynamic_array2![f64, shape]; - - arr.fill_from_seed_equally_distributed(0); - arr2.fill_from(arr.view()); - - let mut col_count = 0; - for (col_index, mut col) in arr.col_iter_mut().enumerate() { - for row_index in 0..shape[0] { - col[[row_index]] *= 2.0; - assert_relative_eq!( - col[[row_index]], - 2.0 * arr2[[row_index, col_index]], - epsilon = 1E-14 - ); - } - col_count += 1; - } - assert_eq!(col_count, shape[1]); - } -} diff --git a/dense/src/array/slice.rs b/dense/src/array/slice.rs index 8429cc96..903a3035 100644 --- a/dense/src/array/slice.rs +++ b/dense/src/array/slice.rs @@ -312,78 +312,3 @@ fn compute_raw_range( 1 + convert_nd_raw(end_multi_index, stride), ) } - -#[cfg(test)] -mod test { - - use crate::traits::*; - use crate::{layout::convert_nd_raw, rlst_dynamic_array3}; - - #[test] - fn test_create_slice() { - let shape = [3, 7, 6]; - let mut arr = rlst_dynamic_array3!(f64, shape); - - arr.fill_from_seed_equally_distributed(0); - - let slice = arr.view().slice(1, 2); - - assert_eq!(slice[[1, 5]], arr[[1, 2, 5]]); - - assert_eq!(slice.shape(), [3, 6]); - - let stride_expected = [arr.stride()[0], arr.stride()[2]]; - let stride_actual = slice.stride(); - - assert_eq!(stride_expected, stride_actual); - - let orig_data = arr.data(); - let slice_data = slice.data(); - - let orig_raw_index = convert_nd_raw([1, 2, 5], arr.stride()); - let slice_raw_index = convert_nd_raw([1, 5], slice.stride()); - - assert_eq!(orig_data[orig_raw_index], slice_data[slice_raw_index]); - assert_eq!(slice_data[slice_raw_index], slice[[1, 5]]); - - let last_raw_index = convert_nd_raw([2, 5], slice.stride()); - assert_eq!(slice_data[last_raw_index], slice[[2, 5]]); - } - - #[test] - fn test_multiple_slices() { - let shape = [3, 7, 6]; - let mut arr = rlst_dynamic_array3!(f64, shape); - arr.fill_from_seed_equally_distributed(0); - - let mut slice = arr.view_mut().slice(1, 3).slice(1, 1); - - slice[[2]] = 5.0; - - assert_eq!(slice.shape(), [3]); - assert_eq!(arr[[2, 3, 1]], 5.0); - } - - #[test] - fn test_slice_of_subview() { - let shape = [3, 7, 6]; - let mut arr = rlst_dynamic_array3!(f64, shape); - arr.fill_from_seed_equally_distributed(0); - let mut arr2 = rlst_dynamic_array3!(f64, shape); - arr2.fill_from(arr.view()); - - let slice = arr.into_subview([1, 2, 4], [2, 3, 2]).slice(1, 2); - - assert_eq!(slice.shape(), [2, 2]); - - assert_eq!(slice[[1, 0]], arr2[[2, 4, 4]]); - - let slice_data = slice.data(); - let arr_data = arr2.data(); - - let slice_index = convert_nd_raw([1, 0], slice.stride()); - let array_index = convert_nd_raw([2, 4, 4], arr2.stride()); - - assert_eq!(slice_data[slice_index], arr_data[array_index]); - } -} diff --git a/dense/src/gemm.rs b/dense/src/gemm.rs index 418ac28b..87bc97ce 100644 --- a/dense/src/gemm.rs +++ b/dense/src/gemm.rs @@ -1,5 +1,6 @@ //! Gemm trait for matrix multiplication -use crate::types::TransMode; +use crate::types::{c32, c64, TransMode}; +use blas::{cgemm, dgemm, sgemm, zgemm}; pub trait Gemm: Sized { #[allow(clippy::too_many_arguments)] @@ -23,4 +24,145 @@ pub trait Gemm: Sized { ); } -pub mod blis; +/// Compute expected size of a data slice from stride and shape. +fn get_expected_data_size(stride: [usize; 2], shape: [usize; 2]) -> usize { + if shape[0] == 0 || shape[1] == 0 { + return 0; + } + + 1 + (shape[0] - 1) * stride[0] + (shape[1] - 1) * stride[1] +} + +/// Panic if expected data size is not identical to actual data size. +fn assert_data_size(nelems: usize, stride: [usize; 2], shape: [usize; 2]) { + let expected = get_expected_data_size(stride, shape); + + assert_eq!( + expected, nelems, + "Wrong size for data slice. Actual size {}. Expected size {}.", + nelems, expected + ); +} + +macro_rules! impl_gemm { + ($scalar:ty, $gemm:ident) => { + impl Gemm for $scalar { + fn gemm( + transa: TransMode, + transb: TransMode, + m: usize, + n: usize, + k: usize, + alpha: Self, + a: &[Self], + rsa: usize, + csa: usize, + b: &[Self], + rsb: usize, + csb: usize, + beta: Self, + c: &mut [Self], + rsc: usize, + csc: usize, + ) { + match transa { + TransMode::NoTrans => assert_data_size(a.len(), [rsa, csa], [m, k]), + TransMode::ConjNoTrans => assert_data_size(a.len(), [rsa, csa], [m, k]), + TransMode::Trans => assert_data_size(a.len(), [rsa, csa], [k, m]), + TransMode::ConjTrans => assert_data_size(a.len(), [rsa, csa], [k, m]), + } + + match transb { + TransMode::NoTrans => assert_data_size(b.len(), [rsb, csb], [k, n]), + TransMode::ConjNoTrans => assert_data_size(b.len(), [rsb, csb], [k, n]), + TransMode::Trans => assert_data_size(b.len(), [rsb, csb], [n, k]), + TransMode::ConjTrans => assert_data_size(b.len(), [rsb, csb], [n, k]), + } + + assert_data_size(c.len(), [rsc, csc], [m, n]); + + assert_eq!( + rsa, 1, + "Input matrix for dgemm must have rsa=1. Actual {}", + rsa + ); + assert_eq!( + rsa, 1, + "Input matrix for dgemm must have rsb=1. Actual {}", + rsb + ); + assert_eq!( + rsa, 1, + "Input matrix for dgemm must have rsc=1. Actual {}", + rsc + ); + + let lda = csa as i32; + let ldb = csb as i32; + let ldc = csc as i32; + + let transa = match transa { + TransMode::NoTrans => b'N', + TransMode::ConjNoTrans => { + panic!("TransMode::ConjNoTrans not supported for gemm implementation.") + } + TransMode::Trans => b'T', + TransMode::ConjTrans => b'C', + }; + + let transb = match transb { + TransMode::NoTrans => b'N', + TransMode::ConjNoTrans => { + panic!("TransMode::ConjNoTrans not supported for gemm implementation.") + } + TransMode::Trans => b'T', + TransMode::ConjTrans => b'C', + }; + + unsafe { + $gemm( + transa as u8, + transb as u8, + m as i32, + n as i32, + k as i32, + alpha, + a, + lda, + b, + ldb, + beta, + c, + ldc, + ) + } + + // unsafe { + // raw::$blis_gemm( + // convert_trans_mode(transa), + // convert_trans_mode(transb), + // m as i64, + // n as i64, + // k as i64, + // &alpha as *const _ as *const $blis_type, + // a.as_ptr() as *const _ as *const $blis_type, + // rsa as i64, + // csa as i64, + // b.as_ptr() as *const _ as *const $blis_type, + // rsb as i64, + // csb as i64, + // &beta as *const _ as *const $blis_type, + // c.as_mut_ptr() as *mut _ as *mut $blis_type, + // rsc as i64, + // csc as i64, + // ) + // }; + } + } + }; +} + +impl_gemm!(f32, sgemm); +impl_gemm!(f64, dgemm); +impl_gemm!(c32, cgemm); +impl_gemm!(c64, zgemm); diff --git a/dense/src/gemm/blis.rs b/dense/src/gemm/blis.rs deleted file mode 100644 index 9a9aaa6e..00000000 --- a/dense/src/gemm/blis.rs +++ /dev/null @@ -1,102 +0,0 @@ -//! Blis implementation of Gemm - -extern crate rlst_blis_src; -use crate::gemm::Gemm; -use crate::types::{c32, c64, TransMode}; -use rlst_blis::raw; - -/// Compute expected size of a data slice from stride and shape. -fn get_expected_data_size(stride: [usize; 2], shape: [usize; 2]) -> usize { - if shape[0] == 0 || shape[1] == 0 { - return 0; - } - - 1 + (shape[0] - 1) * stride[0] + (shape[1] - 1) * stride[1] -} - -/// Panic if expected data size is not identical to actual data size. -fn assert_data_size(nelems: usize, stride: [usize; 2], shape: [usize; 2]) { - let expected = get_expected_data_size(stride, shape); - - assert_eq!( - expected, nelems, - "Wrong size for data slice. Actual size {}. Expected size {}.", - nelems, expected - ); -} - -fn convert_trans_mode(trans: TransMode) -> u32 { - match trans { - TransMode::ConjNoTrans => raw::trans_t_BLIS_CONJ_NO_TRANSPOSE, - TransMode::NoTrans => raw::trans_t_BLIS_NO_TRANSPOSE, - TransMode::Trans => raw::trans_t_BLIS_TRANSPOSE, - TransMode::ConjTrans => raw::trans_t_BLIS_CONJ_TRANSPOSE, - } -} - -macro_rules! impl_gemm { - ($scalar:ty, $blis_gemm:ident, $blis_type:ty) => { - impl Gemm for $scalar { - fn gemm( - transa: TransMode, - transb: TransMode, - m: usize, - n: usize, - k: usize, - alpha: Self, - a: &[Self], - rsa: usize, - csa: usize, - b: &[Self], - rsb: usize, - csb: usize, - beta: Self, - c: &mut [Self], - rsc: usize, - csc: usize, - ) { - match transa { - TransMode::NoTrans => assert_data_size(a.len(), [rsa, csa], [m, k]), - TransMode::ConjNoTrans => assert_data_size(a.len(), [rsa, csa], [m, k]), - TransMode::Trans => assert_data_size(a.len(), [rsa, csa], [k, m]), - TransMode::ConjTrans => assert_data_size(a.len(), [rsa, csa], [k, m]), - } - - match transb { - TransMode::NoTrans => assert_data_size(b.len(), [rsb, csb], [k, n]), - TransMode::ConjNoTrans => assert_data_size(b.len(), [rsb, csb], [k, n]), - TransMode::Trans => assert_data_size(b.len(), [rsb, csb], [n, k]), - TransMode::ConjTrans => assert_data_size(b.len(), [rsb, csb], [n, k]), - } - - assert_data_size(c.len(), [rsc, csc], [m, n]); - - unsafe { - raw::$blis_gemm( - convert_trans_mode(transa), - convert_trans_mode(transb), - m as i64, - n as i64, - k as i64, - &alpha as *const _ as *const $blis_type, - a.as_ptr() as *const _ as *const $blis_type, - rsa as i64, - csa as i64, - b.as_ptr() as *const _ as *const $blis_type, - rsb as i64, - csb as i64, - &beta as *const _ as *const $blis_type, - c.as_mut_ptr() as *mut _ as *mut $blis_type, - rsc as i64, - csc as i64, - ) - }; - } - } - }; -} - -impl_gemm!(f32, bli_sgemm, f32); -impl_gemm!(f64, bli_dgemm, f64); -impl_gemm!(c32, bli_cgemm, raw::scomplex); -impl_gemm!(c64, bli_zgemm, raw::dcomplex); diff --git a/dense/src/layout.rs b/dense/src/layout.rs index 1959d89b..e83da536 100644 --- a/dense/src/layout.rs +++ b/dense/src/layout.rs @@ -60,24 +60,3 @@ pub fn convert_1d_nd_from_shape( } res } - -#[cfg(test)] -mod test { - use super::*; - - #[test] - fn test_convert_1d_nd() { - let multi_index: [usize; 3] = [2, 3, 7]; - let shape: [usize; 3] = [3, 4, 8]; - let stride = stride_from_shape(shape); - - let index_1d = convert_nd_raw(multi_index, stride); - let actual_nd = convert_1d_nd_from_shape(index_1d, shape); - - println!("{}, {:#?}", index_1d, actual_nd); - - for (&expected, actual) in multi_index.iter().zip(actual_nd) { - assert_eq!(expected, actual) - } - } -} diff --git a/dense/src/lib.rs b/dense/src/lib.rs index e7c67e8e..f3c836a4 100644 --- a/dense/src/lib.rs +++ b/dense/src/lib.rs @@ -2,9 +2,6 @@ //! #![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/inverse.rs b/dense/src/linalg/inverse.rs index 32d6d32b..eb7a292f 100644 --- a/dense/src/linalg/inverse.rs +++ b/dense/src/linalg/inverse.rs @@ -18,7 +18,7 @@ use super::assert_lapack_stride; /// /// The following command computes the inverse of an array `a`. The content /// of `a` is replaced by the inverse. -/// ``` +/// ```ignore /// # use rlst_dense::rlst_dynamic_array2; /// # use rlst_dense::linalg::inverse::MatrixInverse; /// # let mut a = rlst_dynamic_array2!(f64, [3, 3]); @@ -89,48 +89,3 @@ impl_inverse!(f64, dgetrf, dgetri); impl_inverse!(f32, sgetrf, sgetri); impl_inverse!(c32, cgetrf, cgetri); impl_inverse!(c64, zgetrf, zgetri); - -#[cfg(test)] -mod test { - - use super::*; - - use crate::assert_array_abs_diff_eq; - use paste::paste; - - use crate::array::empty_array; - use crate::rlst_dynamic_array2; - - macro_rules! impl_inverse_tests { - ($scalar:ty, $tol:expr) => { - paste! { - - #[test] - fn []() { - let n = 4; - - let mut a = rlst_dynamic_array2!($scalar, [n, n]); - let mut b = rlst_dynamic_array2!($scalar, [n, n]); - - let mut ident = rlst_dynamic_array2!($scalar, [n, n]); - ident.set_identity(); - - a.fill_from_seed_equally_distributed(0); - b.fill_from(a.view()); - - a.view_mut().into_inverse_alloc().unwrap(); - - let actual = empty_array::<$scalar, 2>().simple_mult_into_resize(a.view(), b.view()); - - assert_array_abs_diff_eq!(actual, ident, $tol); - } - - } - }; - } - - impl_inverse_tests!(f64, 1E-12); - impl_inverse_tests!(f32, 5E-6); - impl_inverse_tests!(c32, 5E-6); - impl_inverse_tests!(c64, 1E-12); -} diff --git a/dense/src/linalg/lu.rs b/dense/src/linalg/lu.rs index 591838ff..52961738 100644 --- a/dense/src/linalg/lu.rs +++ b/dense/src/linalg/lu.rs @@ -455,144 +455,3 @@ impl_lu!(f64, dgetrf, dgetrs); impl_lu!(f32, sgetrf, sgetrs); impl_lu!(c64, zgetrf, zgetrs); impl_lu!(c32, cgetrf, cgetrs); - -#[cfg(test)] -mod test { - - use crate::assert_array_relative_eq; - use approx; - - use crate::rlst_dynamic_array1; - use crate::rlst_dynamic_array2; - - use super::*; - use crate::array::empty_array; - - use paste::paste; - - macro_rules! impl_lu_tests { - - ($scalar:ty, $tol:expr) => { - paste! { - #[test] - fn []() { - let dim = [8, 20]; - let mut arr = rlst_dynamic_array2!($scalar, dim); - - arr.fill_from_seed_normally_distributed(0); - let mut arr2 = rlst_dynamic_array2!($scalar, dim); - arr2.fill_from(arr.view()); - - let lu = LuDecomposition::<$scalar,_>::new(arr2).unwrap(); - - let mut l_mat = empty_array::<$scalar, 2>(); - let mut u_mat = empty_array::<$scalar, 2>(); - let mut p_mat = empty_array::<$scalar, 2>(); - - lu.get_l_resize(l_mat.view_mut()); - lu.get_u_resize(u_mat.view_mut()); - lu.get_p_resize(p_mat.view_mut()); - - let res = crate::array::empty_array::<$scalar, 2>(); - - let res = - res.simple_mult_into_resize(empty_array().simple_mult_into_resize(p_mat, l_mat), u_mat); - - assert_array_relative_eq!(res, arr, $tol) - } - - #[test] - fn []() { - let dim = [12, 12]; - let mut arr = rlst_dynamic_array2!($scalar, dim); - - arr.fill_from_seed_normally_distributed(0); - let mut arr2 = rlst_dynamic_array2!($scalar, dim); - arr2.fill_from(arr.view()); - - let lu = LuDecomposition::<$scalar, _>::new(arr2).unwrap(); - - let mut l_mat = empty_array::<$scalar, 2>(); - let mut u_mat = empty_array::<$scalar, 2>(); - let mut p_mat = empty_array::<$scalar, 2>(); - - lu.get_l_resize(l_mat.view_mut()); - lu.get_u_resize(u_mat.view_mut()); - lu.get_p_resize(p_mat.view_mut()); - - let res = crate::array::empty_array::<$scalar, 2>(); - - let res = - res.simple_mult_into_resize(empty_array().simple_mult_into_resize(p_mat, l_mat), u_mat); - - assert_array_relative_eq!(res, arr, $tol) - } - - #[test] - fn []() { - let dim = [12, 12]; - let mut arr = rlst_dynamic_array2!($scalar, dim); - arr.fill_from_seed_equally_distributed(0); - let mut x_actual = rlst_dynamic_array1!($scalar, [dim[0]]); - let mut rhs = rlst_dynamic_array1!($scalar, [dim[0]]); - x_actual.fill_from_seed_equally_distributed(1); - rhs.view_mut().simple_mult_into_resize(arr.view(), x_actual.view()); - - let lu = LuDecomposition::<$scalar,_>::new(arr).unwrap(); - lu.solve_vec(LuTrans::NoTrans, rhs.view_mut()).unwrap(); - - assert_array_relative_eq!(x_actual, rhs, $tol) - } - - - - #[test] - fn []() { - let dim = [12, 8]; - let mut arr = rlst_dynamic_array2!($scalar, dim); - - arr.fill_from_seed_normally_distributed(0); - let mut arr2 = rlst_dynamic_array2!($scalar, dim); - arr2.fill_from(arr.view()); - - let lu = LuDecomposition::<$scalar,_>::new(arr2).unwrap(); - - let mut l_mat = empty_array::<$scalar, 2>(); - let mut u_mat = empty_array::<$scalar, 2>(); - let mut p_mat = empty_array::<$scalar, 2>(); - - lu.get_l_resize(l_mat.view_mut()); - lu.get_u_resize(u_mat.view_mut()); - lu.get_p_resize(p_mat.view_mut()); - - let res = crate::array::empty_array::<$scalar, 2>(); - - let res = - res.simple_mult_into_resize(empty_array().simple_mult_into_resize(p_mat, l_mat), u_mat); - - assert_array_relative_eq!(res, arr, $tol) - } - - #[test] - fn []() { - let dim = [2, 2]; - let mut arr = rlst_dynamic_array2!($scalar, dim); - arr[[0, 1]] = $scalar::from_real(3.0); - arr[[1, 0]] = $scalar::from_real(2.0); - - let det = LuDecomposition::<$scalar, _>::new(arr).unwrap().det(); - - approx::assert_relative_eq!(det, $scalar::from_real(-6.0), epsilon=$tol); - } - - - - } - }; - } - - impl_lu_tests!(f64, 1E-12); - impl_lu_tests!(f32, 1E-5); - impl_lu_tests!(c64, 1E-12); - impl_lu_tests!(c32, 1E-5); -} diff --git a/dense/src/linalg/pseudo_inverse.rs b/dense/src/linalg/pseudo_inverse.rs index 5ba1f687..ca48e379 100644 --- a/dense/src/linalg/pseudo_inverse.rs +++ b/dense/src/linalg/pseudo_inverse.rs @@ -152,80 +152,3 @@ impl_pinv!(f64); impl_pinv!(f32); impl_pinv!(c32); impl_pinv!(c64); - -#[cfg(test)] -mod test { - use super::*; - - use crate::array::empty_array; - use crate::assert_array_abs_diff_eq; - use paste::paste; - - use crate::rlst_dynamic_array2; - - macro_rules! impl_pinv_tests { - ($scalar:ty, $tol:expr) => { - paste! { - - #[test] - fn []() { - let shape = [5, 10]; - let tol = 0.0; - - let mut mat = rlst_dynamic_array2!($scalar, shape); - let mut mat2 = rlst_dynamic_array2!($scalar, shape); - - let mut pinv = rlst_dynamic_array2!($scalar, [shape[1], shape[0]]); - let mut ident = rlst_dynamic_array2!($scalar, [shape[0], shape[0]]); - ident.set_identity(); - - mat.fill_from_seed_equally_distributed(0); - mat2.fill_from(mat.view()); - - mat2.into_pseudo_inverse_alloc(pinv.view_mut(), tol) - .unwrap(); - - let actual = if shape[0] >= shape[1] { - empty_array::<$scalar, 2>().simple_mult_into_resize(pinv, mat) - } else { - empty_array::<$scalar, 2>().simple_mult_into_resize(mat, pinv) - }; - - assert_array_abs_diff_eq!(actual, ident, $tol); - } - - #[test] - fn []() { - let shape = [10, 5]; - let tol = 0.0; - - let mut mat = rlst_dynamic_array2!($scalar, shape); - let mut mat2 = rlst_dynamic_array2!($scalar, shape); - - let mut pinv = rlst_dynamic_array2!($scalar, [shape[1], shape[0]]); - let mut ident = rlst_dynamic_array2!($scalar, [shape[1], shape[1]]); - ident.set_identity(); - - mat.fill_from_seed_equally_distributed(0); - mat2.fill_from(mat.view()); - - mat2.into_pseudo_inverse_alloc(pinv.view_mut(), tol) - .unwrap(); - - let actual = if shape[0] >= shape[1] { - empty_array::<$scalar, 2>().simple_mult_into_resize(pinv, mat) - } else { - empty_array::<$scalar, 2>().simple_mult_into_resize(mat, pinv) - }; - - assert_array_abs_diff_eq!(actual, ident, $tol); - } - } - }; - } - - impl_pinv_tests!(f32, 1E-5); - impl_pinv_tests!(f64, 1E-12); - impl_pinv_tests!(c32, 1E-5); - impl_pinv_tests!(c64, 1E-12); -} diff --git a/dense/src/linalg/qr.rs b/dense/src/linalg/qr.rs index 3e39df43..ee58deeb 100644 --- a/dense/src/linalg/qr.rs +++ b/dense/src/linalg/qr.rs @@ -656,112 +656,3 @@ implement_qr_real!(f64, dgeqp3, dormqr); implement_qr_real!(f32, sgeqp3, sormqr); implement_qr_complex!(c64, zgeqp3, zunmqr); implement_qr_complex!(c32, cgeqp3, cunmqr); - -#[cfg(test)] -mod test { - - use crate::{assert_array_abs_diff_eq, assert_array_relative_eq, traits::*}; - - use crate::array::empty_array; - use crate::rlst_dynamic_array2; - use paste::paste; - - use super::*; - - macro_rules! implement_qr_tests { - ($scalar:ty, $tol:expr) => { - paste! { - - #[test] - pub fn []() { - let shape = [8, 5]; - let mut mat = rlst_dynamic_array2!($scalar, shape); - let mut mat2 = rlst_dynamic_array2!($scalar, shape); - - mat.fill_from_seed_equally_distributed(0); - mat2.fill_from(mat.view()); - - let mut r_mat = rlst_dynamic_array2!($scalar, [5, 5]); - let mut q_mat = rlst_dynamic_array2!($scalar, [8, 5]); - let mut p_mat = rlst_dynamic_array2!($scalar, [5, 5]); - let mut p_trans = rlst_dynamic_array2!($scalar, [5, 5]); - let actual = rlst_dynamic_array2!($scalar, [8, 5]); - let mut ident = rlst_dynamic_array2!($scalar, [5, 5]); - ident.set_identity(); - - let qr = QrDecomposition::<$scalar,_>::new(mat).unwrap(); - - let _ = qr.get_r(r_mat.view_mut()); - let _ = qr.get_q_alloc(q_mat.view_mut()); - let _ = qr.get_p(p_mat.view_mut()); - - p_trans.fill_from(p_mat.transpose()); - - let actual = empty_array::<$scalar, 2>() - .simple_mult_into_resize(actual.simple_mult_into(q_mat.view(), r_mat.view()), p_trans); - - assert_array_relative_eq!(actual, mat2, $tol); - - let qtq = empty_array::<$scalar, 2>().mult_into_resize( - crate::types::TransMode::ConjTrans, - crate::types::TransMode::NoTrans, - 1.0.into(), - q_mat.view(), - q_mat.view(), - 1.0.into(), - ); - - assert_array_abs_diff_eq!(qtq, ident, $tol); - } - - #[test] - pub fn []() { - let shape = [5, 8]; - let mut mat = rlst_dynamic_array2!($scalar, shape); - let mut mat2 = rlst_dynamic_array2!($scalar, shape); - - mat.fill_from_seed_equally_distributed(0); - mat2.fill_from(mat.view()); - - let mut r_mat = rlst_dynamic_array2!($scalar, [5, 8]); - let mut q_mat = rlst_dynamic_array2!($scalar, [5, 5]); - let mut p_mat = rlst_dynamic_array2!($scalar, [8, 8]); - let mut p_trans = rlst_dynamic_array2!($scalar, [8, 8]); - let actual = rlst_dynamic_array2!($scalar, [5, 8]); - let mut ident = rlst_dynamic_array2!($scalar, [5, 5]); - ident.set_identity(); - - let qr = QrDecomposition::<$scalar, _>::new(mat).unwrap(); - - let _ = qr.get_r(r_mat.view_mut()); - let _ = qr.get_q_alloc(q_mat.view_mut()); - let _ = qr.get_p(p_mat.view_mut()); - - p_trans.fill_from(p_mat.transpose()); - - let actual = empty_array::<$scalar, 2>() - .simple_mult_into_resize(actual.simple_mult_into(q_mat.view(), r_mat), p_trans); - - assert_array_relative_eq!(actual, mat2, $tol); - - let qtq = empty_array::<$scalar, 2>().mult_into_resize( - crate::types::TransMode::ConjTrans, - crate::types::TransMode::NoTrans, - 1.0.into(), - q_mat.view(), - q_mat.view(), - 1.0.into(), - ); - - assert_array_abs_diff_eq!(qtq, ident, $tol); - } - - } - }; - } - - implement_qr_tests!(f32, 1E-6); - implement_qr_tests!(f64, 1E-12); - implement_qr_tests!(c32, 1E-6); - implement_qr_tests!(c64, 1E-12); -} diff --git a/dense/src/linalg/svd.rs b/dense/src/linalg/svd.rs index c0f5d597..3b714a70 100644 --- a/dense/src/linalg/svd.rs +++ b/dense/src/linalg/svd.rs @@ -460,139 +460,3 @@ impl_svd_real!(f64, dgesvd); impl_svd_real!(f32, sgesvd); impl_svd_complex!(c32, cgesvd); impl_svd_complex!(c64, zgesvd); - -#[cfg(test)] -mod test { - - use super::*; - - use crate::assert_array_relative_eq; - use approx::assert_relative_eq; - use paste::paste; - - use crate::array::empty_array; - use crate::{rlst_dynamic_array1, rlst_dynamic_array2}; - - use crate::linalg::qr::*; - - macro_rules! impl_tests { - ($scalar:ty, $tol:expr) => { - paste! { - - #[test] - fn []() { - [](5, 10, $tol); - [](10, 5, $tol); - } - - #[test] - fn []() { - [](10, 5, SvdMode::Reduced, $tol); - [](5, 10, SvdMode::Reduced, $tol); - [](10, 5, SvdMode::Full, $tol); - [](5, 10, SvdMode::Full, $tol); - } - - fn [](m: usize, n: usize, tol: <$scalar as RlstScalar>::Real) { - let k = std::cmp::min(m, n); - let mut mat = rlst_dynamic_array2!($scalar, [m, m]); - let mut q = rlst_dynamic_array2!($scalar, [m, m]); - let mut sigma = rlst_dynamic_array2!($scalar, [m, n]); - - mat.fill_from_seed_equally_distributed(0); - let qr = QrDecomposition::<$scalar,_>::new(mat).unwrap(); - qr.get_q_alloc(q.view_mut()).unwrap(); - - for index in 0..k { - sigma[[index, index]] = ((k - index) as <$scalar as RlstScalar>::Real).into(); - } - - let a = empty_array::<$scalar, 2>().simple_mult_into_resize(q.view(), sigma.view()); - - let mut singvals = rlst_dynamic_array1!(<$scalar as RlstScalar>::Real, [k]); - - a.into_singular_values_alloc(singvals.data_mut()).unwrap(); - - for index in 0..k { - assert_relative_eq!(singvals[[index]], sigma[[index, index]].re(), epsilon = tol); - } - } - - fn [](m: usize, n: usize, mode: SvdMode, tol: <$scalar as RlstScalar>::Real) { - let k = std::cmp::min(m, n); - - let mut mat_u; - let mut u; - let mut mat_vt; - let mut vt; - let mut sigma; - - match mode { - SvdMode::Full => { - mat_u = rlst_dynamic_array2!($scalar, [m, m]); - u = rlst_dynamic_array2!($scalar, [m, m]); - mat_vt = rlst_dynamic_array2!($scalar, [n, n]); - vt = rlst_dynamic_array2!($scalar, [n, n]); - sigma = rlst_dynamic_array2!($scalar, [m, n]); - } - SvdMode::Reduced => { - mat_u = rlst_dynamic_array2!($scalar, [m, k]); - u = rlst_dynamic_array2!($scalar, [m, k]); - mat_vt = rlst_dynamic_array2!($scalar, [k, n]); - vt = rlst_dynamic_array2!($scalar, [k, n]); - sigma = rlst_dynamic_array2!($scalar, [k, k]); - } - } - - mat_u.fill_from_seed_equally_distributed(0); - mat_vt.fill_from_seed_equally_distributed(1); - - let qr = QrDecomposition::<$scalar,_>::new(mat_u).unwrap(); - qr.get_q_alloc(u.view_mut()).unwrap(); - - let qr = QrDecomposition::<$scalar,_>::new(mat_vt).unwrap(); - qr.get_q_alloc(vt.view_mut()).unwrap(); - - for index in 0..k { - sigma[[index, index]] = ((k - index) as <$scalar as RlstScalar>::Real).into(); - } - - let a = empty_array::<$scalar, 2>().simple_mult_into_resize( - empty_array::<$scalar, 2>().simple_mult_into_resize(u.view(), sigma.view()), - vt.view(), - ); - - let mut expected = rlst_dynamic_array2!($scalar, a.shape()); - expected.fill_from(a.view()); - - u.set_zero(); - vt.set_zero(); - sigma.set_zero(); - - let mut singvals = rlst_dynamic_array1!(<$scalar as RlstScalar>::Real, [k]); - - a.into_svd_alloc(u.view_mut(), vt.view_mut(), singvals.data_mut(), mode) - .unwrap(); - - for index in 0..k { - sigma[[index, index]] = singvals[[index]].into(); - } - - let actual = empty_array::<$scalar, 2>().simple_mult_into_resize( - empty_array::<$scalar, 2>().simple_mult_into_resize(u, sigma), - vt, - ); - - assert_array_relative_eq!(expected, actual, tol); - } - - - } - }; - } - - impl_tests!(f32, 1E-5); - impl_tests!(f64, 1E-12); - impl_tests!(c32, 1E-5); - impl_tests!(c64, 1E-12); -} diff --git a/dense/src/matrix_multiply.rs b/dense/src/matrix_multiply.rs index bef0d127..e4868064 100644 --- a/dense/src/matrix_multiply.rs +++ b/dense/src/matrix_multiply.rs @@ -70,137 +70,3 @@ pub fn matrix_multiply< csc, ) } - -#[cfg(test)] -mod test { - - use crate::assert_array_relative_eq; - use crate::types::{c32, c64}; - - use super::*; - use crate::array::Array; - use crate::base_array::BaseArray; - use crate::data_container::VectorContainer; - use crate::rlst_dynamic_array2; - use paste::paste; - - macro_rules! mat_mul_test_impl { - ($ScalarType:ty, $eps:expr) => { - paste! { - fn [](transa: TransMode, transb: TransMode, shape_a: [usize; 2], shape_b: [usize; 2], shape_c: [usize; 2]) { - - let mut mat_a = rlst_dynamic_array2!($ScalarType, shape_a); - let mut mat_b = rlst_dynamic_array2!($ScalarType, shape_b); - let mut mat_c = rlst_dynamic_array2!($ScalarType, shape_c); - let mut expected = rlst_dynamic_array2!($ScalarType, shape_c); - - mat_a.fill_from_seed_equally_distributed(0); - mat_b.fill_from_seed_equally_distributed(1); - //mat_c.fill_from_seed_equally_distributed(2); - - expected.fill_from(mat_c.view_mut()); - - matrix_multiply( - transa, - transb, - <$ScalarType as RlstScalar>::from_real(1.), - &mat_a, - &mat_b, - <$ScalarType as RlstScalar>::from_real(1.), - &mut mat_c, - ); - matrix_multiply_compare( - transa, - transb, - <$ScalarType as RlstScalar>::from_real(1.), - &mat_a, - &mat_b, - <$ScalarType as RlstScalar>::from_real(1.), - &mut expected, - ); - - assert_array_relative_eq!(mat_c, expected, $eps); - } - - #[test] - fn []() { - - [](TransMode::NoTrans, TransMode::NoTrans, [3, 5], [5, 7], [3, 7]); - [](TransMode::ConjNoTrans, TransMode::ConjNoTrans, [3, 5], [5, 7], [3, 7]); - [](TransMode::ConjNoTrans, TransMode::NoTrans, [3, 5], [5, 7], [3, 7]); - [](TransMode::NoTrans, TransMode::ConjNoTrans, [3, 5], [5, 7], [3, 7]); - - [](TransMode::NoTrans, TransMode::Trans, [3, 5], [7, 5], [3, 7]); - [](TransMode::Trans, TransMode::NoTrans, [2, 1], [2, 1], [1, 1]); - [](TransMode::Trans, TransMode::Trans, [5, 3], [7, 5], [3, 7]); - - [](TransMode::NoTrans, TransMode::ConjTrans, [3, 5], [7, 5], [3, 7]); - [](TransMode::ConjTrans, TransMode::NoTrans, [5, 3], [5, 7], [3, 7]); - [](TransMode::ConjTrans, TransMode::ConjTrans, [5, 3], [7, 5], [3, 7]); - - - } - - } - }; - } - - fn matrix_multiply_compare( - transa: TransMode, - transb: TransMode, - alpha: Item, - mat_a: &Array, 2>, 2>, - mat_b: &Array, 2>, 2>, - beta: Item, - mat_c: &mut Array, 2>, 2>, - ) { - let a_shape = match transa { - TransMode::NoTrans => mat_a.shape(), - TransMode::ConjNoTrans => mat_a.shape(), - TransMode::Trans => [mat_a.shape()[1], mat_a.shape()[0]], - TransMode::ConjTrans => [mat_a.shape()[1], mat_a.shape()[0]], - }; - - let b_shape = match transb { - TransMode::NoTrans => mat_b.shape(), - TransMode::ConjNoTrans => mat_b.shape(), - TransMode::Trans => [mat_b.shape()[1], mat_b.shape()[0]], - TransMode::ConjTrans => [mat_b.shape()[1], mat_b.shape()[0]], - }; - - let mut a_actual = rlst_dynamic_array2!(Item, a_shape); - let mut b_actual = rlst_dynamic_array2!(Item, b_shape); - - match transa { - TransMode::NoTrans => a_actual.fill_from(mat_a.view()), - TransMode::ConjNoTrans => a_actual.fill_from(mat_a.view().conj()), - TransMode::Trans => a_actual.fill_from(mat_a.view().transpose()), - TransMode::ConjTrans => a_actual.fill_from(mat_a.view().conj().transpose()), - } - - match transb { - TransMode::NoTrans => b_actual.fill_from(mat_b.view()), - TransMode::ConjNoTrans => b_actual.fill_from(mat_b.view().conj()), - TransMode::Trans => b_actual.fill_from(mat_b.view().transpose()), - TransMode::ConjTrans => b_actual.fill_from(mat_b.view().conj().transpose()), - } - - let m = mat_c.shape()[0]; - let n = mat_c.shape()[1]; - let k = a_actual.shape()[1]; - - for row in 0..m { - for col in 0..n { - mat_c[[row, col]] *= beta; - for index in 0..k { - mat_c[[row, col]] += alpha * a_actual[[row, index]] * b_actual[[index, col]]; - } - } - } - } - - mat_mul_test_impl!(f64, 1E-14); - mat_mul_test_impl!(f32, 1E-5); - mat_mul_test_impl!(c32, 1E-5); - mat_mul_test_impl!(c64, 1E-14); -} diff --git a/dense/src/traits.rs b/dense/src/traits.rs index 848224dd..e75c49e3 100644 --- a/dense/src/traits.rs +++ b/dense/src/traits.rs @@ -127,3 +127,12 @@ pub trait AijIterator { fn iter_aij(&self) -> Self::Iter<'_>; } + +/// Helper trait that returns from an enumeration iterator a new iterator +/// that converts the 1d index into a multi-index. +pub trait AsMultiIndex, const NDIM: usize> { + fn multi_index( + self, + shape: [usize; NDIM], + ) -> crate::array::iterators::MultiIndexIterator; +} diff --git a/gfortran-src/Cargo.toml b/gfortran-src/Cargo.toml deleted file mode 100644 index c04eeb93..00000000 --- a/gfortran-src/Cargo.toml +++ /dev/null @@ -1,18 +0,0 @@ -[features] -strict = [] - -[package] -name = "rlst-gfortran-src" -version = "0.1.0" -edition = "2021" -links = "gfortran" - -# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html - - -[lib] -name = "rlst_gfortran_src" - -[build-dependencies] - - diff --git a/gfortran-src/build.rs b/gfortran-src/build.rs deleted file mode 100644 index b954e736..00000000 --- a/gfortran-src/build.rs +++ /dev/null @@ -1,6 +0,0 @@ -fn main() { - if cfg!(target_os = "macos") { - println!("cargo:rustc-link-search=/opt/homebrew/opt/gfortran/lib/gcc/current"); - } - println!("cargo:rustc-link-lib=dylib=gfortran"); -} diff --git a/gfortran-src/src/lib.rs b/gfortran-src/src/lib.rs deleted file mode 100644 index f05efe3b..00000000 --- a/gfortran-src/src/lib.rs +++ /dev/null @@ -1 +0,0 @@ -#![cfg_attr(feature = "strict", deny(warnings))] diff --git a/lapack/Cargo.toml b/lapack/Cargo.toml deleted file mode 100644 index 4e1af9f0..00000000 --- a/lapack/Cargo.toml +++ /dev/null @@ -1,41 +0,0 @@ -[features] -strict = [] - -[package] -name = "rlst-lapack" -version = "0.0.1" -authors = ["The linalg-rs contributors"] -edition = "2021" -description = "Safe Lapack wrappers." -license = "MIT + Apache 2.0" -homepage = "https://github.com/linalg-rs" -repository = "https://github.com/linalg-rs/rlst" -readme = "README.md" -keywords = ["numerics"] -categories = ["mathematics", "science"] - -[lib] -name = "rlst_lapack" - -[dependencies] -num = "0.4" -rand = "0.8" -itertools = "0.12" -rand_distr = "0.4" -rlst-blis = { path = "../blis"} -approx = { version = "0.5", features=["num-complex"] } -rlst-common = {path = "../common"} -paste = "1" -rand_chacha = "0.3" -rlst-netlib-lapack-src = {path = "../netlib-lapack-src"} -lapack = "0.19.*" - - -[dev-dependencies] -criterion = { version = "0.3", features = ["html_reports"] } - -[package.metadata.docs.rs] -rustdoc-args = [ "--html-in-header", "katex-header.html" ] - - - diff --git a/lapack/src/geqp3.rs b/lapack/src/geqp3.rs deleted file mode 100644 index 31bd6fb5..00000000 --- a/lapack/src/geqp3.rs +++ /dev/null @@ -1,98 +0,0 @@ -//! Interface to dgeqp3 - -use lapack::{cgeqp3, dgeqp3, sgeqp3, zgeqp3}; -use num::Zero; -use rlst_common::types::*; - -pub trait Geqp3: Scalar { - fn geqp3(m: i32, n: i32, a: &mut [Self], lda: i32, jpvt: &mut [i32], tau: &mut [Self]) -> i32; -} - -macro_rules! impl_geqp3_real { - ($scalar:ty, $geqp3:expr) => { - impl Geqp3 for $scalar { - fn geqp3( - m: i32, - n: i32, - a: &mut [Self], - lda: i32, - jpvt: &mut [i32], - tau: &mut [Self], - ) -> i32 { - assert!(m >= 0); - assert!(n >= 0); - assert!(lda >= 0); - assert_eq!(a.len() as i32, lda * n); - assert_eq!(jpvt.len() as i32, n); - assert_eq!(tau.len() as i32, std::cmp::min(m, n)); - - let mut info = 0; - let lwork = -1; - let mut work_query = [::zero()]; - unsafe { $geqp3(m, n, a, lda, jpvt, tau, &mut work_query, lwork, &mut info) }; - assert_eq!(info, 0); - let lwork = work_query[0].re() as i32; - - let mut work = vec![::zero(); lwork as usize]; - unsafe { $geqp3(m, n, a, lda, jpvt, tau, &mut work, lwork, &mut info) }; - info - } - } - }; -} - -macro_rules! impl_geqp3_complex { - ($scalar:ty, $geqp3:expr) => { - impl Geqp3 for $scalar { - fn geqp3( - m: i32, - n: i32, - a: &mut [Self], - lda: i32, - jpvt: &mut [i32], - tau: &mut [Self], - ) -> i32 { - assert!(m >= 0); - assert!(n >= 0); - assert!(lda >= 0); - assert_eq!(a.len() as i32, lda * n); - assert_eq!(jpvt.len() as i32, n); - assert_eq!(tau.len() as i32, std::cmp::min(m, n)); - let mut rwork = vec![<::Real as Zero>::zero(); 2 * n as usize]; - - let mut info = 0; - let lwork = -1; - let mut work_query = [::zero()]; - unsafe { - $geqp3( - m, - n, - a, - lda, - jpvt, - tau, - &mut work_query, - lwork, - &mut rwork, - &mut info, - ) - }; - assert_eq!(info, 0); - let lwork = work_query[0].re() as i32; - - let mut work = vec![::zero(); lwork as usize]; - unsafe { - $geqp3( - m, n, a, lda, jpvt, tau, &mut work, lwork, &mut rwork, &mut info, - ) - }; - info - } - } - }; -} - -impl_geqp3_real!(f32, sgeqp3); -impl_geqp3_real!(f64, dgeqp3); -impl_geqp3_complex!(c32, cgeqp3); -impl_geqp3_complex!(c64, zgeqp3); diff --git a/lapack/src/getrf.rs b/lapack/src/getrf.rs deleted file mode 100644 index 6b71da8b..00000000 --- a/lapack/src/getrf.rs +++ /dev/null @@ -1,32 +0,0 @@ -//! Getrf -use lapack::{cgetrf, dgetrf, sgetrf, zgetrf}; -use rlst_common::types::*; - -pub trait Getrf: Scalar { - fn getrf(m: i32, n: i32, a: &mut [Self], lda: i32, ipiv: &mut [i32]) -> i32; -} - -macro_rules! impl_getrf { - ($scalar:ty, $getrf:expr) => { - impl Getrf for $scalar { - fn getrf(m: i32, n: i32, a: &mut [Self], lda: i32, ipiv: &mut [i32]) -> i32 { - assert!(m >= 0); - assert!(n >= 0); - assert!(lda >= std::cmp::max(m, 1)); - assert_eq!(a.len() as i32, lda * n); - assert_eq!(ipiv.len() as i32, std::cmp::min(m, n)); - - let mut info = 0; - - unsafe { $getrf(m, n, a, lda, ipiv, &mut info) } - - info - } - } - }; -} - -impl_getrf!(f32, sgetrf); -impl_getrf!(f64, dgetrf); -impl_getrf!(c32, cgetrf); -impl_getrf!(c64, zgetrf); diff --git a/lapack/src/getrs.rs b/lapack/src/getrs.rs deleted file mode 100644 index d2f38d20..00000000 --- a/lapack/src/getrs.rs +++ /dev/null @@ -1,58 +0,0 @@ -//! Getrf -use lapack::{cgetrs, dgetrs, sgetrs, zgetrs}; -use rlst_common::types::*; - -pub trait Getrs: Scalar { - fn getrs( - trans: u8, - n: i32, - nrhs: i32, - a: &[Self], - lda: i32, - ipiv: &[i32], - b: &mut [Self], - ldb: i32, - ) -> i32; -} - -macro_rules! impl_getrs { - ($scalar:ty, $getrs:expr) => { - impl Getrs for $scalar { - fn getrs( - trans: u8, - n: i32, - nrhs: i32, - a: &[Self], - lda: i32, - ipiv: &[i32], - b: &mut [Self], - ldb: i32, - ) -> i32 { - assert!(n >= 0); - assert!(nrhs >= 0); - assert!(lda >= std::cmp::max(1, n)); - assert_eq!(a.len() as i32, lda * n); - assert_eq!(ipiv.len() as i32, n); - assert_eq!(b.len() as i32, ldb * nrhs); - assert!(ldb >= std::cmp::max(1, n)); - assert!(trans == b'N' || trans == b'T' || trans == b'C'); - - for &elem in ipiv { - assert!(elem >= 1); - assert!(elem <= n); - } - - let mut info = 0; - - unsafe { $getrs(trans, n, nrhs, a, lda, ipiv, b, ldb, &mut info) } - - info - } - } - }; -} - -impl_getrs!(f64, dgetrs); -impl_getrs!(f32, sgetrs); -impl_getrs!(c64, zgetrs); -impl_getrs!(c32, cgetrs); diff --git a/lapack/src/lib.rs b/lapack/src/lib.rs deleted file mode 100644 index 15cd6f7b..00000000 --- a/lapack/src/lib.rs +++ /dev/null @@ -1,20 +0,0 @@ -//! Safe wrapper for Lapack - -#![allow(clippy::too_many_arguments)] - -pub mod geqp3; -pub mod getrf; -pub mod getrs; -pub mod ormqr; -pub mod unmqr; - -pub use geqp3::Geqp3; -pub use getrf::Getrf; -pub use getrs::Getrs; -pub use ormqr::Ormqr; -pub use unmqr::Unmqr; - -// // Collective Lapack wrapper trait -// pub trait Lapack: Scalar + Getrf + Getrs + Unmqr + Ormqr {} - -// impl Lapack for T {} diff --git a/lapack/src/ormqr.rs b/lapack/src/ormqr.rs deleted file mode 100644 index 6874120e..00000000 --- a/lapack/src/ormqr.rs +++ /dev/null @@ -1,253 +0,0 @@ -//! Implementation of ormqr - -use lapack::cunmqr; -use lapack::dormqr; -use lapack::sormqr; -use lapack::zunmqr; -use num::Zero; -use rlst_common::types::*; - -pub trait Ormqr: Scalar { - fn ormqr( - side: u8, - trans: u8, - m: i32, - n: i32, - k: i32, - a: &mut [Self], - lda: i32, - tau: &mut [Self], - c: &mut [Self], - ldc: i32, - work: Option<&mut [Self]>, - ) -> i32; - - fn ormqr_query_work(side: u8, trans: u8, m: i32, n: i32, k: i32) -> i32; -} - -macro_rules! impl_ormqr_complex { - ($scalar:ty, $ormqr:expr) => { - impl Ormqr for $scalar { - fn ormqr( - side: u8, - trans: u8, - m: i32, - n: i32, - k: i32, - a: &mut [Self], - lda: i32, - tau: &mut [Self], - c: &mut [Self], - ldc: i32, - work: Option<&mut [Self]>, - ) -> i32 { - assert!(side == b'L' || side == b'R'); - assert!(trans == b'C' || trans == b'N'); - assert!(m >= 0); - assert!(n >= 0); - - assert!(if side == b'L' { - k >= 0 && k <= m - } else { - k >= 0 && k <= n - }); - assert!(if side == b'L' { - lda >= std::cmp::max(1, m) - } else { - lda >= std::cmp::max(1, n) - }); - assert_eq!(a.len() as i32, lda * k); - assert_eq!(tau.len() as i32, k); - assert!(ldc >= std::cmp::max(1, m)); - assert_eq!(c.len() as i32, ldc * n); - let mut my_work = Vec::::new(); - let work = if let Some(work) = work { - assert!(if side == b'L' { - work.len() as i32 >= std::cmp::max(1, n) - } else { - work.len() as i32 >= std::cmp::max(1, m) - }); - work - } else { - let len = ::ormqr_query_work(side, trans, m, n, k) as usize; - my_work.resize(len, ::zero()); - &mut my_work - }; - - let mut info = 0; - unsafe { - $ormqr( - side, - trans, - m, - n, - k, - a, - lda, - tau, - c, - ldc, - work, - work.len() as i32, - &mut info, - ) - }; - info - } - - fn ormqr_query_work(side: u8, trans: u8, m: i32, n: i32, k: i32) -> i32 { - let a = [::zero()]; - let tau = [::zero()]; - let mut c = [::zero()]; - let mut work_query = [::zero()]; - let mut info = 0; - - assert!(side == b'L' || side == b'R'); - assert!(trans == b'T' || trans == b'N'); - assert!(m >= 0); - assert!(n >= 0); - assert!(if side == b'L' { - k >= 0 && k <= m - } else { - k >= 0 && k <= n - }); - let lda = if side == b'L' { m } else { n }; - unsafe { - $ormqr( - side, - trans, - m, - n, - k, - &a, - lda, - &tau, - &mut c, - m, - &mut work_query, - -1, - &mut info, - ); - } - assert_eq!(info, 0); - work_query[0].re() as i32 - } - } - }; -} - -macro_rules! impl_ormqr_real { - ($scalar:ty, $ormqr:expr) => { - impl Ormqr for $scalar { - fn ormqr( - side: u8, - trans: u8, - m: i32, - n: i32, - k: i32, - a: &mut [Self], - lda: i32, - tau: &mut [Self], - c: &mut [Self], - ldc: i32, - work: Option<&mut [Self]>, - ) -> i32 { - assert!(side == b'L' || side == b'R'); - assert!(trans == b'T' || trans == b'N'); - assert!(m >= 0); - assert!(n >= 0); - - assert!(if side == b'L' { - k >= 0 && k <= m - } else { - k >= 0 && k <= n - }); - assert!(if side == b'L' { - lda >= std::cmp::max(1, m) - } else { - lda >= std::cmp::max(1, n) - }); - assert_eq!(a.len() as i32, lda * k); - assert_eq!(tau.len() as i32, k); - assert!(ldc >= std::cmp::max(1, m)); - assert_eq!(c.len() as i32, ldc * n); - let mut my_work = Vec::::new(); - let work = if let Some(work) = work { - assert!(if side == b'L' { - work.len() as i32 >= std::cmp::max(1, n) - } else { - work.len() as i32 >= std::cmp::max(1, m) - }); - work - } else { - let len = ::ormqr_query_work(side, trans, m, n, k) as usize; - my_work.resize(len, ::zero()); - &mut my_work - }; - - let mut info = 0; - unsafe { - $ormqr( - side, - trans, - m, - n, - k, - a, - lda, - tau, - c, - ldc, - work, - work.len() as i32, - &mut info, - ) - }; - info - } - - fn ormqr_query_work(side: u8, trans: u8, m: i32, n: i32, k: i32) -> i32 { - let a = [::zero()]; - let tau = [::zero()]; - let mut c = [::zero()]; - let mut work_query = [::zero()]; - let mut info = 0; - - assert!(side == b'L' || side == b'R'); - assert!(trans == b'T' || trans == b'N'); - assert!(m >= 0); - assert!(n >= 0); - assert!(if side == b'L' { - k >= 0 && k <= m - } else { - k >= 0 && k <= n - }); - let lda = if side == b'L' { m } else { n }; - unsafe { - $ormqr( - side, - trans, - m, - n, - k, - &a, - lda, - &tau, - &mut c, - m, - &mut work_query, - -1, - &mut info, - ); - } - assert_eq!(info, 0); - work_query[0].re() as i32 - } - } - }; -} - -impl_ormqr_real!(f32, sormqr); -impl_ormqr_real!(f64, dormqr); -impl_ormqr_complex!(c32, cunmqr); -impl_ormqr_complex!(c64, zunmqr); diff --git a/lapack/src/unmqr.rs b/lapack/src/unmqr.rs deleted file mode 100644 index 56f9d2b0..00000000 --- a/lapack/src/unmqr.rs +++ /dev/null @@ -1,50 +0,0 @@ -//! Implementation of ormqr - -use lapack::cunmqr; -use lapack::dormqr; -use lapack::sormqr; -use lapack::zunmqr; -use num::Zero; -use rlst_common::types::*; - -use crate::Ormqr; - -pub trait Unmqr: Scalar { - fn unmqr( - side: u8, - trans: u8, - m: i32, - n: i32, - k: i32, - a: &mut [Self], - lda: i32, - tau: &mut [Self], - c: &mut [Self], - ldc: i32, - work: Option<&mut [Self]>, - ) -> i32; - - fn unmqr_query_work(side: u8, trans: u8, m: i32, n: i32, k: i32) -> i32; -} - -impl Unmqr for T { - fn unmqr( - side: u8, - trans: u8, - m: i32, - n: i32, - k: i32, - a: &mut [Self], - lda: i32, - tau: &mut [Self], - c: &mut [Self], - ldc: i32, - work: Option<&mut [Self]>, - ) -> i32 { - ::ormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work) - } - - fn unmqr_query_work(side: u8, trans: u8, m: i32, n: i32, k: i32) -> i32 { - ::ormqr_query_work(side, trans, m, n, k) - } -} diff --git a/netlib-lapack-src/Cargo.toml b/netlib-lapack-src/Cargo.toml deleted file mode 100644 index 3d0a68bd..00000000 --- a/netlib-lapack-src/Cargo.toml +++ /dev/null @@ -1,22 +0,0 @@ -[features] -strict = [] - -[package] -name = "rlst-netlib-lapack-src" -version = "0.1.0" -edition = "2021" -links = "netliblapack" - -# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html - -[dependencies] -rlst-blis-src = { path = "../blis-src" } -rlst-gfortran-src = { path = "../gfortran-src"} - -[lib] -name = "rlst_netlib_lapack_src" - -[build-dependencies] -cmake = "0.1.*" - - diff --git a/netlib-lapack-src/build.rs b/netlib-lapack-src/build.rs deleted file mode 100644 index b0f41a77..00000000 --- a/netlib-lapack-src/build.rs +++ /dev/null @@ -1,18 +0,0 @@ -fn main() { - // // Build Lapack static - - let blis_root = std::path::PathBuf::from(std::env::var("DEP_BLIS_ROOT").unwrap()); - - let dst = cmake::Config::new("lapack") - .define("BUILD_SHARED_LIBS", "OFF") - .define("LAPACKE", "ON") - .define("BLA_VENDOR", "FLAME") - .define("CMAKE_PREFIX_PATH", blis_root.display().to_string()) - .define("USE_OPTIMIZED_BLAS", "ON") - .define("CMAKE_POSITION_INDEPENDENT_CODE", "ON") - .build(); - - println!("cargo:rustc-link-search={}", dst.join("lib").display()); - println!("cargo:rustc-link-lib=static=lapack"); - println!("cargo:rustc-link-lib=static=lapacke"); -} diff --git a/netlib-lapack-src/lapack b/netlib-lapack-src/lapack deleted file mode 160000 index ae2f14fe..00000000 --- a/netlib-lapack-src/lapack +++ /dev/null @@ -1 +0,0 @@ -Subproject commit ae2f14fe1f0dd94b2dacc1c35107f0ea5b5c6f07 diff --git a/netlib-lapack-src/src/lib.rs b/netlib-lapack-src/src/lib.rs deleted file mode 100644 index ed9dda9f..00000000 --- a/netlib-lapack-src/src/lib.rs +++ /dev/null @@ -1,3 +0,0 @@ -#![cfg_attr(feature = "strict", deny(warnings))] - -extern crate rlst_gfortran_src; diff --git a/operator/Cargo.toml b/operator/Cargo.toml index 8981d11e..467520bf 100644 --- a/operator/Cargo.toml +++ b/operator/Cargo.toml @@ -11,7 +11,6 @@ edition = "2021" [dependencies] rlst-sparse = {path = "../sparse"} rlst-dense = {path = "../dense"} -rlst-blis = {path = "../blis"} num = "0.4" thiserror = "1.0.40" diff --git a/operator/src/interface.rs b/operator/src/interface.rs index 8f3a81db..20c730f5 100644 --- a/operator/src/interface.rs +++ b/operator/src/interface.rs @@ -3,3 +3,8 @@ pub mod array_vector_space; pub mod dense_matrix_operator; pub mod sparse_operator; + +pub use array_vector_space::{ArrayVectorSpace, ArrayVectorSpaceElement}; +pub use dense_matrix_operator::DenseMatrixOperator; +pub use sparse_operator::CscMatrixOperator; +pub use sparse_operator::CsrMatrixOperator; diff --git a/operator/src/interface/array_vector_space.rs b/operator/src/interface/array_vector_space.rs index 07fc8414..375b3ffd 100644 --- a/operator/src/interface/array_vector_space.rs +++ b/operator/src/interface/array_vector_space.rs @@ -98,10 +98,3 @@ impl Element for ArrayVectorSpaceElement { self.view_mut().scale_in_place(alpha); } } - -#[cfg(test)] -mod test { - - #[test] - fn test_vec() {} -} diff --git a/operator/src/interface/dense_matrix_operator.rs b/operator/src/interface/dense_matrix_operator.rs index a05692ba..f2ef23b8 100644 --- a/operator/src/interface/dense_matrix_operator.rs +++ b/operator/src/interface/dense_matrix_operator.rs @@ -95,33 +95,3 @@ impl< Ok(()) } } - -#[cfg(test)] -mod test { - - use rlst_dense::rlst_dynamic_array2; - - use crate::interface::array_vector_space::ArrayVectorSpace; - use crate::AsApply; - use crate::Element; - use crate::LinearSpace; - use crate::OperatorBase; - - use super::DenseMatrixOperator; - - #[test] - fn test_mat() { - let mut mat = rlst_dynamic_array2!(f64, [3, 4]); - let domain = ArrayVectorSpace::new(4); - let range = ArrayVectorSpace::new(3); - mat.fill_from_seed_equally_distributed(0); - - let op = DenseMatrixOperator::new(mat, &domain, &range); - let mut x = op.domain().zero(); - let mut y = op.range().zero(); - - x.view_mut().fill_from_seed_equally_distributed(0); - - op.apply_extended(1.0, &x, 0.0, &mut y).unwrap(); - } -} diff --git a/operator/src/operations/conjugate_gradients.rs b/operator/src/operations/conjugate_gradients.rs index 3f3c9699..66ed0b6b 100644 --- a/operator/src/operations/conjugate_gradients.rs +++ b/operator/src/operations/conjugate_gradients.rs @@ -125,48 +125,3 @@ impl<'a, Space: InnerProductSpace, Op: AsApply> (self.x, rel_res) } } - -#[cfg(test)] -mod test { - - use crate::interface::{ - array_vector_space::ArrayVectorSpace, dense_matrix_operator::DenseMatrixOperator, - }; - - use super::*; - use rand::prelude::*; - use rlst_dense::rlst_dynamic_array2; - - #[test] - fn test_cg() { - let dim = 10; - let tol = 1E-5; - - let space = ArrayVectorSpace::::new(dim); - let mut residuals = Vec::::new(); - - let mut rng = rand::thread_rng(); - - let mut mat = rlst_dynamic_array2!(f64, [dim, dim]); - - for index in 0..dim { - mat[[index, index]] = rng.gen_range(1.0..=2.0); - } - - let op = DenseMatrixOperator::new(mat.view(), &space, &space); - - let mut rhs = space.zero(); - rhs.view_mut().fill_from_equally_distributed(&mut rng); - - let cg = (CgIteration::new(&op, &rhs)) - .set_callable(|_, res| { - let res_norm = space.norm(res); - residuals.push(res_norm); - }) - .set_tol(tol) - .print_debug(); - - let (_sol, res) = cg.run(); - assert!(res < tol); - } -} diff --git a/operator/src/operations/modified_gram_schmidt.rs b/operator/src/operations/modified_gram_schmidt.rs index eceb6ad2..96040d34 100644 --- a/operator/src/operations/modified_gram_schmidt.rs +++ b/operator/src/operations/modified_gram_schmidt.rs @@ -39,69 +39,3 @@ impl ModifiedGramSchmidt { } } } - -#[cfg(test)] -mod test { - - use approx::{assert_abs_diff_eq, assert_relative_eq}; - use num::Zero; - use rlst_dense::rlst_dynamic_array2; - use rlst_dense::types::c64; - use rlst_dense::{rlst_dynamic_array1, traits::RawAccess}; - - use crate::space::frame::VectorFrame; - use crate::{interface::array_vector_space::ArrayVectorSpace, LinearSpace}; - - use super::*; - - #[test] - pub fn test_gram_schmidt() { - let space = ArrayVectorSpace::::new(5); - let mut vec1 = space.zero(); - let mut vec2 = space.zero(); - let mut vec3 = space.zero(); - - vec1.view_mut().fill_from_seed_equally_distributed(0); - vec2.view_mut().fill_from_seed_equally_distributed(1); - vec3.view_mut().fill_from_seed_equally_distributed(2); - - let mut frame = VectorFrame::default(); - - let mut original = VectorFrame::default(); - - frame.push(vec1); - frame.push(vec2); - frame.push(vec3); - - for elem in frame.iter() { - original.push(space.new_from(elem)); - } - - let mut r_mat = rlst_dynamic_array2!(c64, [3, 3]); - - ModifiedGramSchmidt::orthogonalize(&space, &mut frame, &mut r_mat); - - // Check orthogonality - for index1 in 0..3 { - for index2 in 0..3 { - let inner = space.inner(frame.get(index1).unwrap(), frame.get(index2).unwrap()); - if index1 == index2 { - assert_relative_eq!(inner, c64::one(), epsilon = 1E-12); - } else { - assert_abs_diff_eq!(inner, c64::zero(), epsilon = 1E-12); - } - } - } - - // Check that r is correct. - for (index, col) in r_mat.col_iter().enumerate() { - let mut actual = space.zero(); - let expected = original.get(index).unwrap(); - let mut coeffs = rlst_dynamic_array1!(c64, [frame.len()]); - coeffs.fill_from(col.view()); - frame.evaluate(coeffs.data(), &mut actual); - let rel_diff = (actual.view() - expected.view()).norm_2() / expected.view().norm_2(); - assert_abs_diff_eq!(rel_diff, f64::zero(), epsilon = 1E-12); - } - } -} diff --git a/operator/src/operator.rs b/operator/src/operator.rs index 530591d6..3f31a93c 100644 --- a/operator/src/operator.rs +++ b/operator/src/operator.rs @@ -262,46 +262,3 @@ impl, Op2: AsApply i64; +} +#[cfg(feature = "blis")] /// Get the current number of threads used by Blis. pub fn get_num_threads() -> usize { - let threads = unsafe { raw::bli_thread_get_num_threads() }; + let threads = unsafe { bli_thread_get_num_threads() }; threads as usize } +#[cfg(feature = "blis")] /// Set threads to a given number of threads. pub fn set_num_threads(nthreads: usize) { unsafe { bli_thread_set_num_threads(nthreads as i64) }; } +#[cfg(feature = "blis")] /// Set threads to the number of logical cpus. pub fn enable_threading() { let num_cpus = num_cpus::get(); @@ -22,6 +28,7 @@ pub fn enable_threading() { unsafe { bli_thread_set_num_threads(num_cpus as i64) }; } +#[cfg(feature = "blis")] /// Set number of threads to 1. pub fn disable_threading() { unsafe { bli_thread_set_num_threads(1) }; diff --git a/rlst/tests/array_operations.rs b/rlst/tests/array_operations.rs index 3674831c..b5e024ee 100644 --- a/rlst/tests/array_operations.rs +++ b/rlst/tests/array_operations.rs @@ -1,10 +1,12 @@ //! Tests of array algebray operations use approx::assert_relative_eq; -use rlst::rlst_dynamic_array3; -use rlst_dense::types::*; -use rlst_dense::{array::iterators::AsMultiIndex, layout::convert_1d_nd_from_shape}; -use rlst_dense::{assert_array_relative_eq, traits::*}; +use paste::paste; +use rlst::dense; +use rlst::dense::assert_array_relative_eq; +use rlst::dense::layout::convert_1d_nd_from_shape; +use rlst::prelude::rlst_dynamic_array3; +use rlst::prelude::*; #[test] fn test_addition() { @@ -216,3 +218,335 @@ fn test_to_complex() { assert_array_relative_eq!(res_chunked, expected, 1E-14); assert_array_relative_eq!(res, expected, 1E-14); } + +#[test] +fn test_iter() { + let mut arr = crate::rlst_dynamic_array3![f64, [1, 3, 2]]; + + for (index, item) in arr.iter_mut().enumerate() { + *item = index as f64; + } + + assert_eq!(arr[[0, 0, 0]], 0.0); + assert_eq!(arr[[0, 1, 0]], 1.0); + assert_eq!(arr[[0, 2, 0]], 2.0); + assert_eq!(arr[[0, 0, 1]], 3.0); + assert_eq!(arr[[0, 1, 1]], 4.0); + assert_eq!(arr[[0, 2, 1]], 5.0); +} + +#[test] +fn test_row_iter() { + let shape = [2, 3]; + let mut arr = crate::rlst_dynamic_array2![f64, shape]; + + arr.fill_from_seed_equally_distributed(0); + + let mut row_count = 0; + for (row_index, row) in arr.row_iter().enumerate() { + for col_index in 0..shape[1] { + assert_eq!(row[[col_index]], arr[[row_index, col_index]]); + } + row_count += 1; + } + assert_eq!(row_count, shape[0]); +} + +#[test] +fn test_row_iter_mut() { + let shape = [2, 3]; + let mut arr = crate::rlst_dynamic_array2![f64, shape]; + let mut arr2 = crate::rlst_dynamic_array2![f64, shape]; + + arr.fill_from_seed_equally_distributed(0); + arr2.fill_from(arr.view()); + + let mut row_count = 0; + for (row_index, mut row) in arr.row_iter_mut().enumerate() { + for col_index in 0..shape[1] { + row[[col_index]] *= 2.0; + assert_relative_eq!( + row[[col_index]], + 2.0 * arr2[[row_index, col_index]], + epsilon = 1E-14 + ); + } + row_count += 1; + } + assert_eq!(row_count, shape[0]); +} + +#[test] +fn test_col_iter() { + let shape = [2, 3]; + let mut arr = crate::rlst_dynamic_array2![f64, shape]; + + arr.fill_from_seed_equally_distributed(0); + + let mut col_count = 0; + for (col_index, col) in arr.col_iter().enumerate() { + for row_index in 0..shape[0] { + assert_eq!(col[[row_index]], arr[[row_index, col_index]]); + } + col_count += 1; + } + + assert_eq!(col_count, shape[1]); +} + +#[test] +fn test_col_iter_mut() { + let shape = [2, 3]; + let mut arr = crate::rlst_dynamic_array2![f64, shape]; + let mut arr2 = crate::rlst_dynamic_array2![f64, shape]; + + arr.fill_from_seed_equally_distributed(0); + arr2.fill_from(arr.view()); + + let mut col_count = 0; + for (col_index, mut col) in arr.col_iter_mut().enumerate() { + for row_index in 0..shape[0] { + col[[row_index]] *= 2.0; + assert_relative_eq!( + col[[row_index]], + 2.0 * arr2[[row_index, col_index]], + epsilon = 1E-14 + ); + } + col_count += 1; + } + assert_eq!(col_count, shape[1]); +} + +#[test] +fn test_convert_1d_nd() { + let multi_index: [usize; 3] = [2, 3, 7]; + let shape: [usize; 3] = [3, 4, 8]; + let stride = dense::layout::stride_from_shape(shape); + + let index_1d = dense::layout::convert_nd_raw(multi_index, stride); + let actual_nd = dense::layout::convert_1d_nd_from_shape(index_1d, shape); + + println!("{}, {:#?}", index_1d, actual_nd); + + for (&expected, actual) in multi_index.iter().zip(actual_nd) { + assert_eq!(expected, actual) + } +} + +#[test] +fn test_insert_axis_back() { + let shape = [3, 7, 6]; + let mut arr = rlst_dynamic_array3!(f64, shape); + arr.fill_from_seed_equally_distributed(0); + + let new_arr = arr.view().insert_empty_axis(AxisPosition::Back); + + assert_eq!(new_arr.shape(), [3, 7, 6, 1]); + + assert_eq!(new_arr[[1, 2, 5, 0]], arr[[1, 2, 5]]); + + assert_eq!(new_arr.stride(), [1, 3, 21, 126]) +} + +#[test] +fn test_insert_axis_front() { + let shape = [3, 7, 6]; + let mut arr = rlst_dynamic_array3!(f64, shape); + arr.fill_from_seed_equally_distributed(0); + + let new_arr = arr.view().insert_empty_axis(AxisPosition::Front); + + assert_eq!(new_arr.shape(), [1, 3, 7, 6]); + + assert_eq!(new_arr[[0, 1, 2, 5]], arr[[1, 2, 5]]); + + assert_eq!(new_arr.stride(), [1, 1, 3, 21]) +} + +#[test] +fn test_create_slice() { + let shape = [3, 7, 6]; + let mut arr = rlst_dynamic_array3!(f64, shape); + + arr.fill_from_seed_equally_distributed(0); + + let slice = arr.view().slice(1, 2); + + assert_eq!(slice[[1, 5]], arr[[1, 2, 5]]); + + assert_eq!(slice.shape(), [3, 6]); + + let stride_expected = [arr.stride()[0], arr.stride()[2]]; + let stride_actual = slice.stride(); + + assert_eq!(stride_expected, stride_actual); + + let orig_data = arr.data(); + let slice_data = slice.data(); + + let orig_raw_index = dense::layout::convert_nd_raw([1, 2, 5], arr.stride()); + let slice_raw_index = dense::layout::convert_nd_raw([1, 5], slice.stride()); + + assert_eq!(orig_data[orig_raw_index], slice_data[slice_raw_index]); + assert_eq!(slice_data[slice_raw_index], slice[[1, 5]]); + + let last_raw_index = dense::layout::convert_nd_raw([2, 5], slice.stride()); + assert_eq!(slice_data[last_raw_index], slice[[2, 5]]); +} + +#[test] +fn test_multiple_slices() { + let shape = [3, 7, 6]; + let mut arr = rlst_dynamic_array3!(f64, shape); + arr.fill_from_seed_equally_distributed(0); + + let mut slice = arr.view_mut().slice(1, 3).slice(1, 1); + + slice[[2]] = 5.0; + + assert_eq!(slice.shape(), [3]); + assert_eq!(arr[[2, 3, 1]], 5.0); +} + +#[test] +fn test_slice_of_subview() { + let shape = [3, 7, 6]; + let mut arr = rlst_dynamic_array3!(f64, shape); + arr.fill_from_seed_equally_distributed(0); + let mut arr2 = rlst_dynamic_array3!(f64, shape); + arr2.fill_from(arr.view()); + + let slice = arr.into_subview([1, 2, 4], [2, 3, 2]).slice(1, 2); + + assert_eq!(slice.shape(), [2, 2]); + + assert_eq!(slice[[1, 0]], arr2[[2, 4, 4]]); + + let slice_data = slice.data(); + let arr_data = arr2.data(); + + let slice_index = dense::layout::convert_nd_raw([1, 0], slice.stride()); + let array_index = dense::layout::convert_nd_raw([2, 4, 4], arr2.stride()); + + assert_eq!(slice_data[slice_index], arr_data[array_index]); +} + +macro_rules! mat_mul_test_impl { + ($ScalarType:ty, $eps:expr) => { + paste! { + fn [](transa: TransMode, transb: TransMode, shape_a: [usize; 2], shape_b: [usize; 2], shape_c: [usize; 2]) { + + use rlst::dense::matrix_multiply::matrix_multiply; + + let mut mat_a = rlst_dynamic_array2!($ScalarType, shape_a); + let mut mat_b = rlst_dynamic_array2!($ScalarType, shape_b); + let mut mat_c = rlst_dynamic_array2!($ScalarType, shape_c); + let mut expected = rlst_dynamic_array2!($ScalarType, shape_c); + + mat_a.fill_from_seed_equally_distributed(0); + mat_b.fill_from_seed_equally_distributed(1); + //mat_c.fill_from_seed_equally_distributed(2); + + expected.fill_from(mat_c.view_mut()); + + matrix_multiply( + transa, + transb, + <$ScalarType as RlstScalar>::from_real(1.), + &mat_a, + &mat_b, + <$ScalarType as RlstScalar>::from_real(1.), + &mut mat_c, + ); + matrix_multiply_compare( + transa, + transb, + <$ScalarType as RlstScalar>::from_real(1.), + &mat_a, + &mat_b, + <$ScalarType as RlstScalar>::from_real(1.), + &mut expected, + ); + + assert_array_relative_eq!(mat_c, expected, $eps); + } + + #[test] + fn []() { + + [](TransMode::NoTrans, TransMode::NoTrans, [3, 5], [5, 7], [3, 7]); + + [](TransMode::NoTrans, TransMode::Trans, [3, 5], [7, 5], [3, 7]); + [](TransMode::Trans, TransMode::NoTrans, [2, 1], [2, 1], [1, 1]); + [](TransMode::Trans, TransMode::Trans, [5, 3], [7, 5], [3, 7]); + + [](TransMode::NoTrans, TransMode::ConjTrans, [3, 5], [7, 5], [3, 7]); + [](TransMode::ConjTrans, TransMode::NoTrans, [5, 3], [5, 7], [3, 7]); + [](TransMode::ConjTrans, TransMode::ConjTrans, [5, 3], [7, 5], [3, 7]); + + + } + + } + }; + } + +fn matrix_multiply_compare( + transa: TransMode, + transb: TransMode, + alpha: Item, + mat_a: &Array, 2>, 2>, + mat_b: &Array, 2>, 2>, + beta: Item, + mat_c: &mut Array, 2>, 2>, +) { + let a_shape = match transa { + TransMode::NoTrans => mat_a.shape(), + TransMode::ConjNoTrans => mat_a.shape(), + TransMode::Trans => [mat_a.shape()[1], mat_a.shape()[0]], + TransMode::ConjTrans => [mat_a.shape()[1], mat_a.shape()[0]], + }; + + let b_shape = match transb { + TransMode::NoTrans => mat_b.shape(), + TransMode::ConjNoTrans => mat_b.shape(), + TransMode::Trans => [mat_b.shape()[1], mat_b.shape()[0]], + TransMode::ConjTrans => [mat_b.shape()[1], mat_b.shape()[0]], + }; + + let mut a_actual = rlst_dynamic_array2!(Item, a_shape); + let mut b_actual = rlst_dynamic_array2!(Item, b_shape); + + match transa { + TransMode::NoTrans => a_actual.fill_from(mat_a.view()), + TransMode::ConjNoTrans => a_actual.fill_from(mat_a.view().conj()), + TransMode::Trans => a_actual.fill_from(mat_a.view().transpose()), + TransMode::ConjTrans => a_actual.fill_from(mat_a.view().conj().transpose()), + } + + match transb { + TransMode::NoTrans => b_actual.fill_from(mat_b.view()), + TransMode::ConjNoTrans => b_actual.fill_from(mat_b.view().conj()), + TransMode::Trans => b_actual.fill_from(mat_b.view().transpose()), + TransMode::ConjTrans => b_actual.fill_from(mat_b.view().conj().transpose()), + } + + let m = mat_c.shape()[0]; + let n = mat_c.shape()[1]; + let k = a_actual.shape()[1]; + + for row in 0..m { + for col in 0..n { + mat_c[[row, col]] *= beta; + for index in 0..k { + mat_c[[row, col]] += alpha * a_actual[[row, index]] * b_actual[[index, col]]; + } + } + } +} + +mat_mul_test_impl!(f64, 1E-14); +mat_mul_test_impl!(f32, 1E-5); +mat_mul_test_impl!(c32, 1E-5); +mat_mul_test_impl!(c64, 1E-14); diff --git a/rlst/tests/linalg.rs b/rlst/tests/linalg.rs new file mode 100644 index 00000000..a1dbfb78 --- /dev/null +++ b/rlst/tests/linalg.rs @@ -0,0 +1,449 @@ +//! Tests of array algebray operations + +use paste::paste; +use rlst::dense::assert_array_abs_diff_eq; +use rlst::dense::assert_array_relative_eq; +use rlst::prelude::*; + +macro_rules! impl_inverse_tests { + ($scalar:ty, $tol:expr) => { + paste! { + + #[test] + fn []() { + let n = 4; + + let mut a = rlst_dynamic_array2!($scalar, [n, n]); + let mut b = rlst_dynamic_array2!($scalar, [n, n]); + + let mut ident = rlst_dynamic_array2!($scalar, [n, n]); + ident.set_identity(); + + a.fill_from_seed_equally_distributed(0); + b.fill_from(a.view()); + + a.view_mut().into_inverse_alloc().unwrap(); + + let actual = empty_array::<$scalar, 2>().simple_mult_into_resize(a.view(), b.view()); + + assert_array_abs_diff_eq!(actual, ident, $tol); + } + + } + }; + } + +impl_inverse_tests!(f64, 1E-12); +impl_inverse_tests!(f32, 5E-6); +impl_inverse_tests!(c32, 5E-6); +impl_inverse_tests!(c64, 1E-12); + +macro_rules! impl_lu_tests { + + ($scalar:ty, $tol:expr) => { + paste! { + #[test] + fn []() { + let dim = [8, 20]; + let mut arr = rlst_dynamic_array2!($scalar, dim); + + arr.fill_from_seed_normally_distributed(0); + let mut arr2 = rlst_dynamic_array2!($scalar, dim); + arr2.fill_from(arr.view()); + + let lu = LuDecomposition::<$scalar,_>::new(arr2).unwrap(); + + let mut l_mat = empty_array::<$scalar, 2>(); + let mut u_mat = empty_array::<$scalar, 2>(); + let mut p_mat = empty_array::<$scalar, 2>(); + + lu.get_l_resize(l_mat.view_mut()); + lu.get_u_resize(u_mat.view_mut()); + lu.get_p_resize(p_mat.view_mut()); + + let res = empty_array::<$scalar, 2>(); + + let res = + res.simple_mult_into_resize(empty_array().simple_mult_into_resize(p_mat, l_mat), u_mat); + + assert_array_relative_eq!(res, arr, $tol) + } + + #[test] + fn []() { + let dim = [12, 12]; + let mut arr = rlst_dynamic_array2!($scalar, dim); + + arr.fill_from_seed_normally_distributed(0); + let mut arr2 = rlst_dynamic_array2!($scalar, dim); + arr2.fill_from(arr.view()); + + let lu = LuDecomposition::<$scalar, _>::new(arr2).unwrap(); + + let mut l_mat = empty_array::<$scalar, 2>(); + let mut u_mat = empty_array::<$scalar, 2>(); + let mut p_mat = empty_array::<$scalar, 2>(); + + lu.get_l_resize(l_mat.view_mut()); + lu.get_u_resize(u_mat.view_mut()); + lu.get_p_resize(p_mat.view_mut()); + + let res = empty_array::<$scalar, 2>(); + + let res = + res.simple_mult_into_resize(empty_array().simple_mult_into_resize(p_mat, l_mat), u_mat); + + assert_array_relative_eq!(res, arr, $tol) + } + + #[test] + fn []() { + let dim = [12, 12]; + let mut arr = rlst_dynamic_array2!($scalar, dim); + arr.fill_from_seed_equally_distributed(0); + let mut x_actual = rlst_dynamic_array1!($scalar, [dim[0]]); + let mut rhs = rlst_dynamic_array1!($scalar, [dim[0]]); + x_actual.fill_from_seed_equally_distributed(1); + rhs.view_mut().simple_mult_into_resize(arr.view(), x_actual.view()); + + let lu = LuDecomposition::<$scalar,_>::new(arr).unwrap(); + lu.solve_vec(LuTrans::NoTrans, rhs.view_mut()).unwrap(); + + assert_array_relative_eq!(x_actual, rhs, $tol) + } + + + + #[test] + fn []() { + let dim = [12, 8]; + let mut arr = rlst_dynamic_array2!($scalar, dim); + + arr.fill_from_seed_normally_distributed(0); + let mut arr2 = rlst_dynamic_array2!($scalar, dim); + arr2.fill_from(arr.view()); + + let lu = LuDecomposition::<$scalar,_>::new(arr2).unwrap(); + + let mut l_mat = empty_array::<$scalar, 2>(); + let mut u_mat = empty_array::<$scalar, 2>(); + let mut p_mat = empty_array::<$scalar, 2>(); + + lu.get_l_resize(l_mat.view_mut()); + lu.get_u_resize(u_mat.view_mut()); + lu.get_p_resize(p_mat.view_mut()); + + let res = empty_array::<$scalar, 2>(); + + let res = + res.simple_mult_into_resize(empty_array().simple_mult_into_resize(p_mat, l_mat), u_mat); + + assert_array_relative_eq!(res, arr, $tol) + } + + #[test] + fn []() { + let dim = [2, 2]; + let mut arr = rlst_dynamic_array2!($scalar, dim); + arr[[0, 1]] = $scalar::from_real(3.0); + arr[[1, 0]] = $scalar::from_real(2.0); + + let det = LuDecomposition::<$scalar, _>::new(arr).unwrap().det(); + + approx::assert_relative_eq!(det, $scalar::from_real(-6.0), epsilon=$tol); + } + + + + } + }; + } + +impl_lu_tests!(f64, 1E-12); +impl_lu_tests!(f32, 1E-5); +impl_lu_tests!(c64, 1E-12); +impl_lu_tests!(c32, 1E-5); + +macro_rules! impl_pinv_tests { + ($scalar:ty, $tol:expr) => { + paste! { + + #[test] + fn []() { + let shape = [5, 10]; + let tol = 0.0; + + let mut mat = rlst_dynamic_array2!($scalar, shape); + let mut mat2 = rlst_dynamic_array2!($scalar, shape); + + let mut pinv = rlst_dynamic_array2!($scalar, [shape[1], shape[0]]); + let mut ident = rlst_dynamic_array2!($scalar, [shape[0], shape[0]]); + ident.set_identity(); + + mat.fill_from_seed_equally_distributed(0); + mat2.fill_from(mat.view()); + + mat2.into_pseudo_inverse_alloc(pinv.view_mut(), tol) + .unwrap(); + + let actual = if shape[0] >= shape[1] { + empty_array::<$scalar, 2>().simple_mult_into_resize(pinv, mat) + } else { + empty_array::<$scalar, 2>().simple_mult_into_resize(mat, pinv) + }; + + assert_array_abs_diff_eq!(actual, ident, $tol); + } + + #[test] + fn []() { + let shape = [10, 5]; + let tol = 0.0; + + let mut mat = rlst_dynamic_array2!($scalar, shape); + let mut mat2 = rlst_dynamic_array2!($scalar, shape); + + let mut pinv = rlst_dynamic_array2!($scalar, [shape[1], shape[0]]); + let mut ident = rlst_dynamic_array2!($scalar, [shape[1], shape[1]]); + ident.set_identity(); + + mat.fill_from_seed_equally_distributed(0); + mat2.fill_from(mat.view()); + + mat2.into_pseudo_inverse_alloc(pinv.view_mut(), tol) + .unwrap(); + + let actual = if shape[0] >= shape[1] { + empty_array::<$scalar, 2>().simple_mult_into_resize(pinv, mat) + } else { + empty_array::<$scalar, 2>().simple_mult_into_resize(mat, pinv) + }; + + assert_array_abs_diff_eq!(actual, ident, $tol); + } + } + }; +} + +impl_pinv_tests!(f32, 1E-5); +impl_pinv_tests!(f64, 1E-12); +impl_pinv_tests!(c32, 1E-5); +impl_pinv_tests!(c64, 1E-12); + +macro_rules! implement_qr_tests { + ($scalar:ty, $tol:expr) => { + paste! { + + #[test] + pub fn []() { + let shape = [8, 5]; + let mut mat = rlst_dynamic_array2!($scalar, shape); + let mut mat2 = rlst_dynamic_array2!($scalar, shape); + + mat.fill_from_seed_equally_distributed(0); + mat2.fill_from(mat.view()); + + let mut r_mat = rlst_dynamic_array2!($scalar, [5, 5]); + let mut q_mat = rlst_dynamic_array2!($scalar, [8, 5]); + let mut p_mat = rlst_dynamic_array2!($scalar, [5, 5]); + let mut p_trans = rlst_dynamic_array2!($scalar, [5, 5]); + let actual = rlst_dynamic_array2!($scalar, [8, 5]); + let mut ident = rlst_dynamic_array2!($scalar, [5, 5]); + ident.set_identity(); + + let qr = QrDecomposition::<$scalar,_>::new(mat).unwrap(); + + let _ = qr.get_r(r_mat.view_mut()); + let _ = qr.get_q_alloc(q_mat.view_mut()); + let _ = qr.get_p(p_mat.view_mut()); + + p_trans.fill_from(p_mat.transpose()); + + let actual = empty_array::<$scalar, 2>() + .simple_mult_into_resize(actual.simple_mult_into(q_mat.view(), r_mat.view()), p_trans); + + assert_array_relative_eq!(actual, mat2, $tol); + + let qtq = empty_array::<$scalar, 2>().mult_into_resize( + TransMode::ConjTrans, + TransMode::NoTrans, + 1.0.into(), + q_mat.view(), + q_mat.view(), + 1.0.into(), + ); + + assert_array_abs_diff_eq!(qtq, ident, $tol); + } + + #[test] + pub fn []() { + let shape = [5, 8]; + let mut mat = rlst_dynamic_array2!($scalar, shape); + let mut mat2 = rlst_dynamic_array2!($scalar, shape); + + mat.fill_from_seed_equally_distributed(0); + mat2.fill_from(mat.view()); + + let mut r_mat = rlst_dynamic_array2!($scalar, [5, 8]); + let mut q_mat = rlst_dynamic_array2!($scalar, [5, 5]); + let mut p_mat = rlst_dynamic_array2!($scalar, [8, 8]); + let mut p_trans = rlst_dynamic_array2!($scalar, [8, 8]); + let actual = rlst_dynamic_array2!($scalar, [5, 8]); + let mut ident = rlst_dynamic_array2!($scalar, [5, 5]); + ident.set_identity(); + + let qr = QrDecomposition::<$scalar, _>::new(mat).unwrap(); + + let _ = qr.get_r(r_mat.view_mut()); + let _ = qr.get_q_alloc(q_mat.view_mut()); + let _ = qr.get_p(p_mat.view_mut()); + + p_trans.fill_from(p_mat.transpose()); + + let actual = empty_array::<$scalar, 2>() + .simple_mult_into_resize(actual.simple_mult_into(q_mat.view(), r_mat), p_trans); + + assert_array_relative_eq!(actual, mat2, $tol); + + let qtq = empty_array::<$scalar, 2>().mult_into_resize( + TransMode::ConjTrans, + TransMode::NoTrans, + 1.0.into(), + q_mat.view(), + q_mat.view(), + 1.0.into(), + ); + + assert_array_abs_diff_eq!(qtq, ident, $tol); + } + + } + }; + } + +implement_qr_tests!(f32, 1E-6); +implement_qr_tests!(f64, 1E-12); +implement_qr_tests!(c32, 1E-6); +implement_qr_tests!(c64, 1E-12); + +macro_rules! impl_tests { + ($scalar:ty, $tol:expr) => { + paste! { + + #[test] + fn []() { + [](5, 10, $tol); + [](10, 5, $tol); + } + + #[test] + fn []() { + [](10, 5, SvdMode::Reduced, $tol); + [](5, 10, SvdMode::Reduced, $tol); + [](10, 5, SvdMode::Full, $tol); + [](5, 10, SvdMode::Full, $tol); + } + + fn [](m: usize, n: usize, tol: <$scalar as RlstScalar>::Real) { + let k = std::cmp::min(m, n); + let mut mat = rlst_dynamic_array2!($scalar, [m, m]); + let mut q = rlst_dynamic_array2!($scalar, [m, m]); + let mut sigma = rlst_dynamic_array2!($scalar, [m, n]); + + mat.fill_from_seed_equally_distributed(0); + let qr = QrDecomposition::<$scalar,_>::new(mat).unwrap(); + qr.get_q_alloc(q.view_mut()).unwrap(); + + for index in 0..k { + sigma[[index, index]] = ((k - index) as <$scalar as RlstScalar>::Real).into(); + } + + let a = empty_array::<$scalar, 2>().simple_mult_into_resize(q.view(), sigma.view()); + + let mut singvals = rlst_dynamic_array1!(<$scalar as RlstScalar>::Real, [k]); + + a.into_singular_values_alloc(singvals.data_mut()).unwrap(); + + for index in 0..k { + approx::assert_relative_eq!(singvals[[index]], sigma[[index, index]].re(), epsilon = tol); + } + } + + fn [](m: usize, n: usize, mode: SvdMode, tol: <$scalar as RlstScalar>::Real) { + let k = std::cmp::min(m, n); + + let mut mat_u; + let mut u; + let mut mat_vt; + let mut vt; + let mut sigma; + + match mode { + SvdMode::Full => { + mat_u = rlst_dynamic_array2!($scalar, [m, m]); + u = rlst_dynamic_array2!($scalar, [m, m]); + mat_vt = rlst_dynamic_array2!($scalar, [n, n]); + vt = rlst_dynamic_array2!($scalar, [n, n]); + sigma = rlst_dynamic_array2!($scalar, [m, n]); + } + SvdMode::Reduced => { + mat_u = rlst_dynamic_array2!($scalar, [m, k]); + u = rlst_dynamic_array2!($scalar, [m, k]); + mat_vt = rlst_dynamic_array2!($scalar, [k, n]); + vt = rlst_dynamic_array2!($scalar, [k, n]); + sigma = rlst_dynamic_array2!($scalar, [k, k]); + } + } + + mat_u.fill_from_seed_equally_distributed(0); + mat_vt.fill_from_seed_equally_distributed(1); + + let qr = QrDecomposition::<$scalar,_>::new(mat_u).unwrap(); + qr.get_q_alloc(u.view_mut()).unwrap(); + + let qr = QrDecomposition::<$scalar,_>::new(mat_vt).unwrap(); + qr.get_q_alloc(vt.view_mut()).unwrap(); + + for index in 0..k { + sigma[[index, index]] = ((k - index) as <$scalar as RlstScalar>::Real).into(); + } + + let a = empty_array::<$scalar, 2>().simple_mult_into_resize( + empty_array::<$scalar, 2>().simple_mult_into_resize(u.view(), sigma.view()), + vt.view(), + ); + + let mut expected = rlst_dynamic_array2!($scalar, a.shape()); + expected.fill_from(a.view()); + + u.set_zero(); + vt.set_zero(); + sigma.set_zero(); + + let mut singvals = rlst_dynamic_array1!(<$scalar as RlstScalar>::Real, [k]); + + a.into_svd_alloc(u.view_mut(), vt.view_mut(), singvals.data_mut(), mode) + .unwrap(); + + for index in 0..k { + sigma[[index, index]] = singvals[[index]].into(); + } + + let actual = empty_array::<$scalar, 2>().simple_mult_into_resize( + empty_array::<$scalar, 2>().simple_mult_into_resize(u, sigma), + vt, + ); + + assert_array_relative_eq!(expected, actual, tol); + } + + + } + }; + } + +impl_tests!(f32, 1E-5); +impl_tests!(f64, 1E-12); +impl_tests!(c32, 1E-5); +impl_tests!(c64, 1E-12); diff --git a/rlst/tests/operator.rs b/rlst/tests/operator.rs new file mode 100644 index 00000000..a5b6eb19 --- /dev/null +++ b/rlst/tests/operator.rs @@ -0,0 +1,137 @@ +//! Tests for the operator interface + +use num::traits::{One, Zero}; +use rand::Rng; +use rlst::dense; +use rlst::prelude::*; + +#[test] +fn test_dense_matrix_operator() { + let mut mat = rlst_dynamic_array2!(f64, [3, 4]); + let domain = ArrayVectorSpace::new(4); + let range = ArrayVectorSpace::new(3); + mat.fill_from_seed_equally_distributed(0); + + let op = DenseMatrixOperator::new(mat, &domain, &range); + let mut x = op.domain().zero(); + let mut y = op.range().zero(); + + x.view_mut().fill_from_seed_equally_distributed(0); + + op.apply_extended(1.0, &x, 0.0, &mut y).unwrap(); +} + +#[test] +pub fn test_gram_schmidt() { + let space = ArrayVectorSpace::::new(5); + let mut vec1 = space.zero(); + let mut vec2 = space.zero(); + let mut vec3 = space.zero(); + + vec1.view_mut().fill_from_seed_equally_distributed(0); + vec2.view_mut().fill_from_seed_equally_distributed(1); + vec3.view_mut().fill_from_seed_equally_distributed(2); + + let mut frame = VectorFrame::default(); + + let mut original = VectorFrame::default(); + + frame.push(vec1); + frame.push(vec2); + frame.push(vec3); + + for elem in frame.iter() { + original.push(space.new_from(elem)); + } + + let mut r_mat = rlst_dynamic_array2!(c64, [3, 3]); + + ModifiedGramSchmidt::orthogonalize(&space, &mut frame, &mut r_mat); + + // Check orthogonality + for index1 in 0..3 { + for index2 in 0..3 { + let inner = space.inner(frame.get(index1).unwrap(), frame.get(index2).unwrap()); + if index1 == index2 { + approx::assert_relative_eq!(inner, c64::one(), epsilon = 1E-12); + } else { + approx::assert_abs_diff_eq!(inner, c64::zero(), epsilon = 1E-12); + } + } + } + + // Check that r is correct. + for (index, col) in r_mat.col_iter().enumerate() { + let mut actual = space.zero(); + let expected = original.get(index).unwrap(); + let mut coeffs = rlst_dynamic_array1!(c64, [frame.len()]); + coeffs.fill_from(col.view()); + frame.evaluate(coeffs.data(), &mut actual); + let rel_diff = (actual.view() - expected.view()).norm_2() / expected.view().norm_2(); + approx::assert_abs_diff_eq!(rel_diff, f64::zero(), epsilon = 1E-12); + } +} + +#[test] +fn test_cg() { + let dim = 10; + let tol = 1E-5; + + let space = ArrayVectorSpace::::new(dim); + let mut residuals = Vec::::new(); + + let mut rng = rand::thread_rng(); + + let mut mat = rlst_dynamic_array2!(f64, [dim, dim]); + + for index in 0..dim { + mat[[index, index]] = rng.gen_range(1.0..=2.0); + } + + let op = DenseMatrixOperator::new(mat.view(), &space, &space); + + let mut rhs = space.zero(); + rhs.view_mut().fill_from_equally_distributed(&mut rng); + + let cg = (CgIteration::new(&op, &rhs)) + .set_callable(|_, res| { + let res_norm = space.norm(res); + residuals.push(res_norm); + }) + .set_tol(tol) + .print_debug(); + + let (_sol, res) = cg.run(); + assert!(res < tol); +} + +#[test] +fn test_operator_algebra() { + let mut mat1 = rlst_dynamic_array2!(f64, [4, 3]); + let mut mat2 = rlst_dynamic_array2!(f64, [4, 3]); + + let domain = ArrayVectorSpace::new(3); + let range = ArrayVectorSpace::new(4); + + mat1.fill_from_seed_equally_distributed(0); + mat2.fill_from_seed_equally_distributed(1); + + let op1 = DenseMatrixOperator::new(mat1, &domain, &range); + let op2 = DenseMatrixOperator::new(mat2, &domain, &range); + + let mut x = domain.zero(); + let mut y = range.zero(); + let mut y_expected = range.zero(); + x.view_mut().fill_from_seed_equally_distributed(2); + y.view_mut().fill_from_seed_equally_distributed(3); + y_expected.view_mut().fill_from(y.view()); + + op2.apply_extended(2.0, &x, 3.5, &mut y_expected).unwrap(); + op1.apply_extended(10.0, &x, 1.0, &mut y_expected).unwrap(); + + let sum = op1.scale(5.0).sum(op2.as_ref_obj()); + + sum.apply_extended(2.0, &x, 3.5, &mut y).unwrap(); + + dense::assert_array_relative_eq!(y.view(), y_expected.view(), 1E-12); +} diff --git a/rlst/tests/sparse.rs b/rlst/tests/sparse.rs new file mode 100644 index 00000000..192fbe84 --- /dev/null +++ b/rlst/tests/sparse.rs @@ -0,0 +1,251 @@ +use rlst::dense; +use rlst::prelude::*; + +#[test] +fn test_csc_from_aij() { + // Test the matrix [[1, 2], [3, 4]] + let rows = vec![0, 0, 1, 1, 0]; + let cols = vec![0, 1, 0, 1, 1]; + let data = vec![1.0, 2.0, 3.0, 4.0, 6.0]; + + let csc = CscMatrix::from_aij([2, 2], &rows, &cols, &data).unwrap(); + + assert_eq!(csc.data().len(), 4); + assert_eq!(csc.indices().len(), 4); + assert_eq!(csc.indptr().len(), 3); + assert_eq!(csc.data()[2], 8.0); + + // Test the matrix [[0, 2.0, 0.0], [0, 0, 0], [0, 0, 0]] + let rows = vec![0, 2, 0]; + let cols = vec![1, 2, 1]; + let data = vec![2.0, 0.0, 3.0]; + + let csc = CscMatrix::from_aij([3, 3], &rows, &cols, &data).unwrap(); + + assert_eq!(csc.indptr()[0], 0); + assert_eq!(csc.indptr()[1], 0); + assert_eq!(csc.indptr()[2], 1); + assert_eq!(csc.indptr()[3], 1); + assert_eq!(csc.data()[0], 5.0); +} + +#[test] +fn test_csc_matmul() { + // Test the matrix [[1, 2], [3, 4]] + let rows = vec![0, 0, 1, 1]; + let cols = vec![0, 1, 0, 1]; + let data = vec![1.0, 2.0, 3.0, 4.0]; + + let csc = CscMatrix::from_aij([2, 2], &rows, &cols, &data).unwrap(); + + // Execute 2 * [1, 2] + 3 * A*x with x = [3, 4]; + // Expected result is [35, 79]. + + let x = vec![3.0, 4.0]; + let mut res = vec![1.0, 2.0]; + + csc.matmul(3.0, &x, 2.0, &mut res); + + assert_eq!(res[0], 35.0); + assert_eq!(res[1], 79.0); +} + +#[test] +fn test_csc_aij_iterator() { + let rows = vec![2, 3, 4, 4, 6]; + let cols = vec![1, 1, 3, 3, 4]; + let data = vec![1.0, 2.0, 3.0, 4.0, 5.0]; + + // This matrix has a zero col at the beginning one in between and several zero cols at the end. + let csr = CscMatrix::from_aij([10, 20], &rows, &cols, &data).unwrap(); + + let aij_data: Vec<(usize, usize, f64)> = csr.iter_aij().collect(); + + assert_eq!(aij_data.len(), 4); + + assert_eq!(aij_data[0], (2, 1, 1.0)); + assert_eq!(aij_data[1], (3, 1, 2.0)); + assert_eq!(aij_data[2], (4, 3, 7.0)); + assert_eq!(aij_data[3], (6, 4, 5.0)); +} +#[test] +fn test_csr_from_aij() { + // Test the matrix [[1, 2], [3, 4]] + let rows = vec![0, 0, 1, 1, 0]; + let cols = vec![0, 1, 0, 1, 1]; + let data = vec![1.0, 2.0, 3.0, 4.0, 6.0]; + + let csr = CsrMatrix::from_aij([2, 2], &rows, &cols, &data).unwrap(); + + assert_eq!(csr.data().len(), 4); + assert_eq!(csr.indices().len(), 4); + assert_eq!(csr.indptr().len(), 3); + assert_eq!(csr.data()[1], 8.0); + + //Test the matrix [[0, 0, 0], [2.0, 0, 0], [0, 0, 0]] + let rows = vec![1, 2, 1]; + let cols = vec![0, 2, 0]; + let data = vec![2.0, 0.0, 3.0]; + + let csr = CsrMatrix::from_aij([3, 3], &rows, &cols, &data).unwrap(); + + assert_eq!(csr.indptr()[0], 0); + assert_eq!(csr.indptr()[1], 0); + assert_eq!(csr.indptr()[2], 1); + assert_eq!(csr.indptr()[3], 1); + assert_eq!(csr.data()[0], 5.0); +} + +#[test] +fn test_csr_matmul() { + // Test the matrix [[1, 2], [3, 4]] + let rows = vec![0, 0, 1, 1]; + let cols = vec![0, 1, 0, 1]; + let data = vec![1.0, 2.0, 3.0, 4.0]; + + let csr = CsrMatrix::from_aij([2, 2], &rows, &cols, &data).unwrap(); + + // Execute 2 * [1, 2] + 3 * A*x with x = [3, 4]; + // Expected result is [35, 79]. + + let x = vec![3.0, 4.0]; + let mut res = vec![1.0, 2.0]; + + csr.matmul(3.0, &x, 2.0, &mut res); + + assert_eq!(res[0], 35.0); + assert_eq!(res[1], 79.0); +} + +#[test] +fn test_csr_aij_iterator() { + let rows = vec![2, 3, 4, 4, 6]; + let cols = vec![0, 1, 0, 2, 1]; + let data = vec![1.0, 2.0, 3.0, 4.0, 5.0]; + + // This matrix has a zero row at the beginning one in between and several zero rows at the end. + let csr = CsrMatrix::from_aij([10, 3], &rows, &cols, &data).unwrap(); + + let aij_data: Vec<(usize, usize, f64)> = csr.iter_aij().collect(); + + assert_eq!(aij_data.len(), 5); + + assert_eq!(aij_data[0], (2, 0, 1.0)); + assert_eq!(aij_data[1], (3, 1, 2.0)); + assert_eq!(aij_data[2], (4, 0, 3.0)); + assert_eq!(aij_data[3], (4, 2, 4.0)); + assert_eq!(aij_data[4], (6, 1, 5.0)); +} + +#[test] +fn test_distributed_index_set() { + use mpi; + + let universe = mpi::initialize().unwrap(); + let world = universe.world(); + + let index_layout = DefaultMpiIndexLayout::new(14, &world); + + // Test that the range is correct on rank 0 + assert_eq!(index_layout.index_range(0).unwrap(), (0, 14)); + + // Test that the number of global indices is correct. + assert_eq!(index_layout.number_of_global_indices(), 14); + + // Test that map works + + assert_eq!(index_layout.local2global(2).unwrap(), 2); + + // Return the correct process for an index + + assert_eq!(index_layout.rank_from_index(5).unwrap(), 0); +} + +#[test] +fn test_local_index_layout() { + let index_layout = DefaultSerialIndexLayout::new(14); + + // Test that the range is correct on rank 0 + assert_eq!(index_layout.index_range(0).unwrap(), (0, 14)); + + // Test that the number of global indices is correct. + assert_eq!(index_layout.number_of_global_indices(), 14); + + // Test that map works + + assert_eq!(index_layout.local2global(2).unwrap(), 2); +} + +#[cfg(feature = "umfpack")] +#[test] +fn test_csc_umfpack_f64() { + let n = 5; + + let mut mat = rlst_dynamic_array2!(f64, [n, n]); + let mut x_exact = rlst_dynamic_array1!(f64, [n]); + let mut x_actual = rlst_dynamic_array1!(f64, [n]); + + mat.fill_from_seed_equally_distributed(0); + x_exact.fill_from_seed_equally_distributed(1); + + let rhs = empty_array::().simple_mult_into_resize(mat.view(), x_exact.view()); + + let mut rows = Vec::::with_capacity(n * n); + let mut cols = Vec::::with_capacity(n * n); + let mut data = Vec::::with_capacity(n * n); + + for col_index in 0..n { + for row_index in 0..n { + rows.push(row_index); + cols.push(col_index); + data.push(mat[[row_index, col_index]]); + } + } + + let sparse_mat = CscMatrix::from_aij([n, n], &rows, &cols, &data).unwrap(); + + sparse_mat + .into_lu() + .unwrap() + .solve(rhs.view(), x_actual.view_mut(), TransMode::NoTrans) + .unwrap(); + + dense::assert_array_relative_eq!(x_actual, x_exact, 1E-12); +} + +#[cfg(feature = "umfpack")] +#[test] +fn test_csc_umfpack_c64() { + let n = 5; + + let mut mat = rlst_dynamic_array2!(c64, [n, n]); + let mut x_exact = rlst_dynamic_array1!(c64, [n]); + let mut x_actual = rlst_dynamic_array1!(c64, [n]); + + mat.fill_from_seed_equally_distributed(0); + x_exact.fill_from_seed_equally_distributed(1); + + let rhs = empty_array::().simple_mult_into_resize(mat.view(), x_exact.view()); + + let mut rows = Vec::::with_capacity(n * n); + let mut cols = Vec::::with_capacity(n * n); + let mut data = Vec::::with_capacity(n * n); + + for col_index in 0..n { + for row_index in 0..n { + rows.push(row_index); + cols.push(col_index); + data.push(mat[[row_index, col_index]]); + } + } + + let sparse_mat = CscMatrix::from_aij([n, n], &rows, &cols, &data).unwrap(); + + sparse_mat + .into_lu() + .unwrap() + .solve(rhs.view(), x_actual.view_mut(), TransMode::NoTrans) + .unwrap(); + + dense::assert_array_relative_eq!(x_actual, x_exact, 1E-12); +} diff --git a/sparse/Cargo.toml b/sparse/Cargo.toml index 7a53da12..ce1b6227 100644 --- a/sparse/Cargo.toml +++ b/sparse/Cargo.toml @@ -14,7 +14,7 @@ rand = "0.8" rand_distr = "0.4" approx = { version = "0.5", features=["num-complex"] } itertools = "0.12" -rlst-umfpack = { path = "../umfpack"} +rlst-umfpack = { path = "../umfpack", optional = true} @@ -27,6 +27,8 @@ rand_chacha = "0.3" name = "rlst_sparse" [features] -default = ["mpi","mpi-sys"] -mpi = ["dep:mpi"] +default = [] +mpi = ["dep:mpi", "mpi-sys"] +umfpack = ["dep:rlst-umfpack"] + strict = [] diff --git a/sparse/src/index_layout.rs b/sparse/src/index_layout.rs index 6f7d4158..6269e413 100644 --- a/sparse/src/index_layout.rs +++ b/sparse/src/index_layout.rs @@ -1,7 +1,9 @@ #[cfg(feature = "mpi")] pub mod default_mpi_index_layout; + pub mod default_serial_index_layout; #[cfg(feature = "mpi")] pub use default_mpi_index_layout::*; + pub use default_serial_index_layout::*; diff --git a/sparse/src/index_layout/default_mpi_index_layout.rs b/sparse/src/index_layout/default_mpi_index_layout.rs index 99e4915c..81d73435 100644 --- a/sparse/src/index_layout/default_mpi_index_layout.rs +++ b/sparse/src/index_layout/default_mpi_index_layout.rs @@ -134,32 +134,3 @@ impl<'a, C: Communicator> IndexLayout for DefaultMpiIndexLayout<'a, C> { None } } - -#[cfg(test)] -mod test { - - use super::*; - use mpi; - - #[test] - fn test_distributed_index_set() { - let universe = mpi::initialize().unwrap(); - let world = universe.world(); - - let index_layout = DefaultMpiIndexLayout::new(14, &world); - - // Test that the range is correct on rank 0 - assert_eq!(index_layout.index_range(0).unwrap(), (0, 14)); - - // Test that the number of global indices is correct. - assert_eq!(index_layout.number_of_global_indices(), 14); - - // Test that map works - - assert_eq!(index_layout.local2global(2).unwrap(), 2); - - // Return the correct process for an index - - assert_eq!(index_layout.rank_from_index(5).unwrap(), 0); - } -} diff --git a/sparse/src/index_layout/default_serial_index_layout.rs b/sparse/src/index_layout/default_serial_index_layout.rs index 430a275e..716153b7 100644 --- a/sparse/src/index_layout/default_serial_index_layout.rs +++ b/sparse/src/index_layout/default_serial_index_layout.rs @@ -56,24 +56,3 @@ impl IndexLayout for DefaultSerialIndexLayout { } } } - -#[cfg(test)] -mod test { - - use super::*; - - #[test] - fn test_local_index_layout() { - let index_layout = DefaultSerialIndexLayout::new(14); - - // Test that the range is correct on rank 0 - assert_eq!(index_layout.index_range(0).unwrap(), (0, 14)); - - // Test that the number of global indices is correct. - assert_eq!(index_layout.number_of_global_indices(), 14); - - // Test that map works - - assert_eq!(index_layout.local2global(2).unwrap(), 2); - } -} diff --git a/sparse/src/lib.rs b/sparse/src/lib.rs index bf81b398..08e962ea 100644 --- a/sparse/src/lib.rs +++ b/sparse/src/lib.rs @@ -1,15 +1,12 @@ #![cfg_attr(feature = "strict", deny(warnings))] +#[cfg(feature = "mpi")] pub mod distributed_vector; + +#[cfg(feature = "mpi")] pub mod ghost_communicator; + pub mod index_layout; -// pub mod operator_interface; pub mod sparse; pub mod tools; pub mod traits; - -// #[cfg(feature = "mpi")] -// pub use distributed_vector::*; - -// #[cfg(test)] -// mod tests {} diff --git a/sparse/src/sparse.rs b/sparse/src/sparse.rs index 791e0096..aaf39d0a 100644 --- a/sparse/src/sparse.rs +++ b/sparse/src/sparse.rs @@ -1,7 +1,11 @@ pub mod csc_mat; pub mod csr_mat; + +#[cfg(feature = "mpi")] pub mod mpi_csr_mat; + pub mod tools; +#[cfg(feature = "umfpack")] pub mod umfpack; #[derive(Copy, Clone)] diff --git a/sparse/src/sparse/csc_mat.rs b/sparse/src/sparse/csc_mat.rs index b08ec2d1..f9985798 100644 --- a/sparse/src/sparse/csc_mat.rs +++ b/sparse/src/sparse/csc_mat.rs @@ -198,79 +198,3 @@ impl rlst_dense::traits::Shape<2> for CscMatrix { self.shape } } - -#[cfg(test)] -mod test { - - use rlst_dense::traits::AijIterator; - - use super::*; - - #[test] - fn test_csc_from_aij() { - // Test the matrix [[1, 2], [3, 4]] - let rows = vec![0, 0, 1, 1, 0]; - let cols = vec![0, 1, 0, 1, 1]; - let data = vec![1.0, 2.0, 3.0, 4.0, 6.0]; - - let csc = CscMatrix::from_aij([2, 2], &rows, &cols, &data).unwrap(); - - assert_eq!(csc.data().len(), 4); - assert_eq!(csc.indices().len(), 4); - assert_eq!(csc.indptr().len(), 3); - assert_eq!(csc.data()[2], 8.0); - - // Test the matrix [[0, 2.0, 0.0], [0, 0, 0], [0, 0, 0]] - let rows = vec![0, 2, 0]; - let cols = vec![1, 2, 1]; - let data = vec![2.0, 0.0, 3.0]; - - let csc = CscMatrix::from_aij([3, 3], &rows, &cols, &data).unwrap(); - - assert_eq!(csc.indptr()[0], 0); - assert_eq!(csc.indptr()[1], 0); - assert_eq!(csc.indptr()[2], 1); - assert_eq!(csc.indptr()[3], 1); - assert_eq!(csc.data()[0], 5.0); - } - - #[test] - fn test_csc_matmul() { - // Test the matrix [[1, 2], [3, 4]] - let rows = vec![0, 0, 1, 1]; - let cols = vec![0, 1, 0, 1]; - let data = vec![1.0, 2.0, 3.0, 4.0]; - - let csc = CscMatrix::from_aij([2, 2], &rows, &cols, &data).unwrap(); - - // Execute 2 * [1, 2] + 3 * A*x with x = [3, 4]; - // Expected result is [35, 79]. - - let x = vec![3.0, 4.0]; - let mut res = vec![1.0, 2.0]; - - csc.matmul(3.0, &x, 2.0, &mut res); - - assert_eq!(res[0], 35.0); - assert_eq!(res[1], 79.0); - } - - #[test] - fn test_aij_iterator() { - let rows = vec![2, 3, 4, 4, 6]; - let cols = vec![1, 1, 3, 3, 4]; - let data = vec![1.0, 2.0, 3.0, 4.0, 5.0]; - - // This matrix has a zero col at the beginning one in between and several zero cols at the end. - let csr = CscMatrix::from_aij([10, 20], &rows, &cols, &data).unwrap(); - - let aij_data: Vec<(usize, usize, f64)> = csr.iter_aij().collect(); - - assert_eq!(aij_data.len(), 4); - - assert_eq!(aij_data[0], (2, 1, 1.0)); - assert_eq!(aij_data[1], (3, 1, 2.0)); - assert_eq!(aij_data[2], (4, 3, 7.0)); - assert_eq!(aij_data[3], (6, 4, 5.0)); - } -} diff --git a/sparse/src/sparse/csr_mat.rs b/sparse/src/sparse/csr_mat.rs index 5c59c995..0a4172c5 100644 --- a/sparse/src/sparse/csr_mat.rs +++ b/sparse/src/sparse/csr_mat.rs @@ -196,80 +196,3 @@ impl Shape<2> for CsrMatrix { self.shape } } - -#[cfg(test)] -mod test { - - use rlst_dense::traits::AijIterator; - - use super::*; - - #[test] - fn test_csr_from_aij() { - // Test the matrix [[1, 2], [3, 4]] - let rows = vec![0, 0, 1, 1, 0]; - let cols = vec![0, 1, 0, 1, 1]; - let data = vec![1.0, 2.0, 3.0, 4.0, 6.0]; - - let csr = CsrMatrix::from_aij([2, 2], &rows, &cols, &data).unwrap(); - - assert_eq!(csr.data().len(), 4); - assert_eq!(csr.indices().len(), 4); - assert_eq!(csr.indptr().len(), 3); - assert_eq!(csr.data()[1], 8.0); - - //Test the matrix [[0, 0, 0], [2.0, 0, 0], [0, 0, 0]] - let rows = vec![1, 2, 1]; - let cols = vec![0, 2, 0]; - let data = vec![2.0, 0.0, 3.0]; - - let csr = CsrMatrix::from_aij([3, 3], &rows, &cols, &data).unwrap(); - - assert_eq!(csr.indptr()[0], 0); - assert_eq!(csr.indptr()[1], 0); - assert_eq!(csr.indptr()[2], 1); - assert_eq!(csr.indptr()[3], 1); - assert_eq!(csr.data()[0], 5.0); - } - - #[test] - fn test_csr_matmul() { - // Test the matrix [[1, 2], [3, 4]] - let rows = vec![0, 0, 1, 1]; - let cols = vec![0, 1, 0, 1]; - let data = vec![1.0, 2.0, 3.0, 4.0]; - - let csr = CsrMatrix::from_aij([2, 2], &rows, &cols, &data).unwrap(); - - // Execute 2 * [1, 2] + 3 * A*x with x = [3, 4]; - // Expected result is [35, 79]. - - let x = vec![3.0, 4.0]; - let mut res = vec![1.0, 2.0]; - - csr.matmul(3.0, &x, 2.0, &mut res); - - assert_eq!(res[0], 35.0); - assert_eq!(res[1], 79.0); - } - - #[test] - fn test_aij_iterator() { - let rows = vec![2, 3, 4, 4, 6]; - let cols = vec![0, 1, 0, 2, 1]; - let data = vec![1.0, 2.0, 3.0, 4.0, 5.0]; - - // This matrix has a zero row at the beginning one in between and several zero rows at the end. - let csr = CsrMatrix::from_aij([10, 3], &rows, &cols, &data).unwrap(); - - let aij_data: Vec<(usize, usize, f64)> = csr.iter_aij().collect(); - - assert_eq!(aij_data.len(), 5); - - assert_eq!(aij_data[0], (2, 0, 1.0)); - assert_eq!(aij_data[1], (3, 1, 2.0)); - assert_eq!(aij_data[2], (4, 0, 3.0)); - assert_eq!(aij_data[3], (4, 2, 4.0)); - assert_eq!(aij_data[4], (6, 1, 5.0)); - } -} diff --git a/sparse/src/sparse/umfpack.rs b/sparse/src/sparse/umfpack.rs index 8fccc51a..36301f3a 100644 --- a/sparse/src/sparse/umfpack.rs +++ b/sparse/src/sparse/umfpack.rs @@ -6,18 +6,13 @@ use crate::sparse::csc_mat::CscMatrix; use rlst_dense::array::Array; use rlst_dense::traits::{RawAccess, RawAccessMut, Shape, Stride}; use rlst_dense::types::RlstScalar; +use rlst_dense::types::TransMode; use rlst_dense::types::*; use rlst_dense::types::{RlstError, RlstResult}; use rlst_umfpack as umfpack; use super::csr_mat::CsrMatrix; -pub enum TransposeMode { - NoTrans, - Trans, - ConjugateTrans, -} - pub struct UmfpackLu { shape: [usize; 2], indices: Vec, @@ -46,15 +41,16 @@ impl UmfpackLu { &self, rhs: Array, mut x: Array, - trans: TransposeMode, + trans: TransMode, ) -> rlst_dense::types::RlstResult<()> { assert_eq!(rhs.stride()[0], 1); assert_eq!(x.stride()[0], 1); let sys = match trans { - TransposeMode::NoTrans => umfpack::UMFPACK_A, - TransposeMode::Trans => umfpack::UMFPACK_Aat, - TransposeMode::ConjugateTrans => umfpack::UMFPACK_At, + TransMode::NoTrans => umfpack::UMFPACK_A, + TransMode::Trans => umfpack::UMFPACK_Aat, + TransMode::ConjTrans => umfpack::UMFPACK_At, + _ => panic!("Transpose mode not supported for Umfpack"), }; let n = self.shape[0]; @@ -97,15 +93,16 @@ impl UmfpackLu { &self, rhs: Array, mut x: Array, - trans: TransposeMode, + trans: TransMode, ) -> rlst_dense::types::RlstResult<()> { assert_eq!(rhs.stride()[0], 1); assert_eq!(x.stride()[0], 1); let sys = match trans { - TransposeMode::NoTrans => umfpack::UMFPACK_A, - TransposeMode::Trans => umfpack::UMFPACK_Aat, - TransposeMode::ConjugateTrans => umfpack::UMFPACK_At, + TransMode::NoTrans => umfpack::UMFPACK_A, + TransMode::Trans => umfpack::UMFPACK_Aat, + TransMode::ConjTrans => umfpack::UMFPACK_At, + _ => panic!("Transpose mode not supported for Umfpack"), }; let n = self.shape[0]; @@ -274,90 +271,3 @@ impl CscMatrix { Ok(umfpack_lu) } } - -#[cfg(test)] -mod test { - - use rlst_dense::traits::*; - - use rlst_dense::{ - array::empty_array, assert_array_relative_eq, rlst_dynamic_array1, rlst_dynamic_array2, - traits::MultIntoResize, - }; - - use super::*; - - #[test] - fn test_csc_umfpack_f64() { - let n = 5; - - let mut mat = rlst_dynamic_array2!(f64, [n, n]); - let mut x_exact = rlst_dynamic_array1!(f64, [n]); - let mut x_actual = rlst_dynamic_array1!(f64, [n]); - - mat.fill_from_seed_equally_distributed(0); - x_exact.fill_from_seed_equally_distributed(1); - - let rhs = empty_array::().simple_mult_into_resize(mat.view(), x_exact.view()); - - let mut rows = Vec::::with_capacity(n * n); - let mut cols = Vec::::with_capacity(n * n); - let mut data = Vec::::with_capacity(n * n); - - for col_index in 0..n { - for row_index in 0..n { - rows.push(row_index); - cols.push(col_index); - data.push(mat[[row_index, col_index]]); - } - } - - let sparse_mat = - crate::sparse::csc_mat::CscMatrix::from_aij([n, n], &rows, &cols, &data).unwrap(); - - sparse_mat - .into_lu() - .unwrap() - .solve(rhs.view(), x_actual.view_mut(), TransposeMode::NoTrans) - .unwrap(); - - assert_array_relative_eq!(x_actual, x_exact, 1E-12); - } - - #[test] - fn test_csc_umfpack_c64() { - let n = 5; - - let mut mat = rlst_dynamic_array2!(c64, [n, n]); - let mut x_exact = rlst_dynamic_array1!(c64, [n]); - let mut x_actual = rlst_dynamic_array1!(c64, [n]); - - mat.fill_from_seed_equally_distributed(0); - x_exact.fill_from_seed_equally_distributed(1); - - let rhs = empty_array::().simple_mult_into_resize(mat.view(), x_exact.view()); - - let mut rows = Vec::::with_capacity(n * n); - let mut cols = Vec::::with_capacity(n * n); - let mut data = Vec::::with_capacity(n * n); - - for col_index in 0..n { - for row_index in 0..n { - rows.push(row_index); - cols.push(col_index); - data.push(mat[[row_index, col_index]]); - } - } - - let sparse_mat = - crate::sparse::csc_mat::CscMatrix::from_aij([n, n], &rows, &cols, &data).unwrap(); - - sparse_mat - .into_lu() - .unwrap() - .solve(rhs.view(), x_actual.view_mut(), TransposeMode::NoTrans) - .unwrap(); - - assert_array_relative_eq!(x_actual, x_exact, 1E-12); - } -} diff --git a/sparse/src/traits.rs b/sparse/src/traits.rs index 4c04853f..675dad4c 100644 --- a/sparse/src/traits.rs +++ b/sparse/src/traits.rs @@ -1,5 +1,3 @@ // Interface traits pub mod index_layout; -pub mod indexable_matrix; -pub mod indexable_vector; diff --git a/sparse/src/traits/indexable_matrix.rs b/sparse/src/traits/indexable_matrix.rs deleted file mode 100644 index 6fa6f3bd..00000000 --- a/sparse/src/traits/indexable_matrix.rs +++ /dev/null @@ -1,149 +0,0 @@ -// //! Basic trait for matrices - -// use crate::traits::index_layout::IndexLayout; -// use crate::traits::indexable_vector::IndexableVector; -// use rlst_common::types::Scalar; - -// pub enum DenseMatrixLayout { -// RowMajor((usize, usize)), -// ColMajor((usize, usize)), -// } - -// pub trait IndexableMatrix { -// type Item: Scalar; -// type Ind: IndexLayout; - -// type View<'a> -// where -// Self: 'a; -// type ViewMut<'a> -// where -// Self: 'a; - -// fn view<'a>(&'a self) -> Option>; -// fn view_mut<'a>(&'a mut self) -> Option>; - -// fn column_layout(&self) -> &Self::Ind; -// fn row_layout(&self) -> &Self::Ind; - -// fn shape(&self) -> (usize, usize) { -// ( -// self.row_layout().number_of_global_indices(), -// self.column_layout().number_of_global_indices(), -// ) -// } -// } - -// pub trait IndexableDenseMatrixView { -// type T: Scalar; - -// /// Return a reference to the element at position (`row`, `col`). -// unsafe fn get_unchecked(&self, row: usize, col: usize) -> &Self::T; - -// /// Return a reference to the element at position `index` in one-dimensional numbering. -// unsafe fn get1d_unchecked(&self, index: usize) -> &Self::T; - -// /// Return a reference to the element at position (`row`, `col`). -// fn get(&self, row: usize, col: usize) -> Option<&Self::T>; - -// /// Return a reference to the element at position `index` in one-dimensional numbering. -// fn get1d(&self, elem: usize) -> Option<&Self::T>; - -// /// Get a raw data slice. -// fn data(&self) -> &[Self::T]; -// } - -// pub trait IndexableDenseMatrixViewMut { -// type T: Scalar; - -// /// Return a mutable reference to the element at position (`row`, `col`). -// unsafe fn get_unchecked_mut(&mut self, row: usize, col: usize) -> &mut Self::T; - -// /// Return a mutable reference to the element at position `index` in one-dimensional numbering. -// unsafe fn get1d_unchecked_mut(&mut self, index: usize) -> &mut Self::T; - -// /// Return a mutable reference to the element at position (`row`, `col`). -// fn get_mut(&mut self, row: usize, col: usize) -> Option<&mut Self::T>; - -// /// Return a mutable reference to the element at position `index` in one-dimensional numbering. -// fn get1d_mut(&mut self, elem: usize) -> Option<&mut Self::T>; - -// /// Get a mutable raw data slice. -// fn data(&mut self) -> &mut [Self::T]; -// } - -// pub trait CsrMatrixView { -// type T: Scalar; - -// fn indices(&self) -> &[usize]; - -// fn indptr(&self) -> &[usize]; - -// fn data(&self) -> &[Self::T]; -// } - -// /// Compute the Frobenious norm of a matrix. -// pub trait NormFrob: IndexableMatrix { -// fn norm_frob(&self) -> ::Real; -// } - -// /// Compute the 1-Norm of a matrix. -// pub trait Norm1: IndexableMatrix { -// fn norm_1(&self) -> ::Real; -// } - -// /// Compute the 2-Norm of a matrix. -// pub trait Norm2: IndexableMatrix { -// fn norm_2(&self) -> ::Real; -// } - -// /// Compute the supremum norm of a matrix. -// pub trait NormInfty: IndexableMatrix { -// fn norm_infty(&self) -> ::Real; -// } - -// /// Compute the trace of a matrix. -// pub trait Trace: IndexableMatrix { -// fn trace(&self) -> Self::Item; -// } - -// /// Compute the real part of a matrix. -// pub trait Real: IndexableMatrix { -// fn real(&self) -> ::Real; -// } - -// /// Compute the imaginary part of a matrix. -// pub trait Imag: IndexableMatrix { -// fn imag(&self) -> ::Real; -// } - -// /// Swap entries with another matrix. -// pub trait Swap: IndexableMatrix { -// fn swap(&mut self, other: &mut Self) -> rlst_common::types::RlstResult<()>; -// } - -// /// Fill matrix by copying from another matrix. -// pub trait Fill: IndexableMatrix { -// fn fill(&mut self, other: &Self) -> rlst_common::types::RlstResult<()>; -// } - -// /// Multiply entries with a scalar. -// pub trait ScalarMult: IndexableMatrix { -// fn scalar_mult(&mut self, scalar: Self::Item); -// } - -// /// Compute self -> alpha * other + self. -// pub trait MultSumInto: IndexableMatrix { -// fn mult_sum_into( -// &mut self, -// other: &Self, -// scalar: Self::Item, -// ) -> rlst_common::types::RlstResult<()>; -// } - -// /// Compute y-> alpha A x + y -// pub trait MatMul: IndexableMatrix { -// type Vec: IndexableVector; - -// fn matmul(&self, alpha: Self::Item, x: &Self::Vec, beta: Self::Item, y: &mut Self::Vec); -// } diff --git a/sparse/src/traits/indexable_vector.rs b/sparse/src/traits/indexable_vector.rs deleted file mode 100644 index 092ea567..00000000 --- a/sparse/src/traits/indexable_vector.rs +++ /dev/null @@ -1,109 +0,0 @@ -//! An indexable vector is the standard type for n-dimensional containers - -use crate::traits::index_layout::IndexLayout; -use rlst_dense::types::RlstScalar; - -pub trait IndexableVector { - type Item: RlstScalar; - type Ind: IndexLayout; - - type View<'a>: IndexableVectorView - where - Self: 'a; - type ViewMut<'a>: IndexableVectorView - where - Self: 'a; - - fn view(&self) -> Option>; - fn view_mut(&mut self) -> Option>; - - fn index_layout(&self) -> &Self::Ind; -} - -pub trait IndexableVectorView { - type T: RlstScalar; - type Iter<'a>: std::iter::Iterator - where - Self: 'a; - - fn iter(&self) -> Self::Iter<'_>; - - fn get(&self, index: usize) -> Option<&Self::T>; - - /// # Safety - /// `index` must not exceed bounds. - unsafe fn get_unchecked(&self, index: usize) -> &Self::T; - - fn len(&self) -> usize; - - fn is_empty(&self) -> bool { - self.len() == 0 - } - - fn data(&self) -> &[Self::T]; -} - -pub trait IndexableVectorViewMut: IndexableVectorView { - type IterMut<'a>: std::iter::Iterator - where - Self: 'a; - - fn iter_mut(&mut self) -> Self::IterMut<'_>; - - fn get_mut(&mut self, index: usize) -> Option<&mut Self::T>; - - /// # Safety - /// `index` must not exceed bounds. - unsafe fn get_unchecked_mut(&mut self, index: usize) -> &mut Self::T; - - fn data_mut(&mut self) -> &mut [Self::T]; -} - -/// Inner product with another object. -pub trait Inner: IndexableVector { - fn inner(&self, other: &Self) -> rlst_dense::types::RlstResult; -} - -/// Take the sum of the squares of the absolute values of the entries. -pub trait AbsSquareSum: IndexableVector { - fn abs_square_sum(&self) -> ::Real; -} - -/// Return the 1-Norm (Sum of absolute values of the entries). -pub trait Norm1: IndexableVector { - fn norm_1(&self) -> ::Real; -} - -/// Return the 2-Norm (Sqrt of the sum of squares). -pub trait Norm2: IndexableVector { - fn norm_2(&self) -> ::Real; -} - -/// Return the supremum norm (largest absolute value of the entries). -pub trait NormInfty: IndexableVector { - fn norm_infty(&self) -> ::Real; -} - -/// Swap entries with another vector. -pub trait Swap: IndexableVector { - fn swap(&mut self, other: &mut Self) -> rlst_dense::types::RlstResult<()>; -} - -/// Fill vector by copying from another vector. -pub trait Fill: IndexableVector { - fn fill(&mut self, other: &Self) -> rlst_dense::types::RlstResult<()>; -} - -/// Multiply entries with a scalar. -pub trait ScalarMult: IndexableVector { - fn scalar_mult(&mut self, scalar: Self::Item); -} - -/// Compute self -> alpha * other + self. -pub trait MultSumInto: IndexableVector { - fn mult_sum_into( - &mut self, - other: &Self, - scalar: Self::Item, - ) -> rlst_dense::types::RlstResult<()>; -} diff --git a/sparse/src/traits/matrix_traits.rs b/sparse/src/traits/matrix_traits.rs deleted file mode 100644 index 29620877..00000000 --- a/sparse/src/traits/matrix_traits.rs +++ /dev/null @@ -1,71 +0,0 @@ -//! Traits for sparse matrices - -pub use crate::linalg::indexable_matrix::IndexableMatrix; -pub use crate::linalg::indexable_vector::*; -pub use crate::types::{usize, Scalar}; - -/// Compute the Frobenious norm of a matrix. -pub trait NormFrob: IndexableMatrix { - fn norm_frob(&self) -> ::Real; -} - -/// Compute the 1-Norm of a matrix. -pub trait Norm1: IndexableMatrix { - fn norm_1(&self) -> ::Real; -} - -/// Compute the 2-Norm of a matrix. -pub trait Norm2: IndexableMatrix { - fn norm_2(&self) -> ::Real; -} - -/// Compute the supremum norm of a matrix. -pub trait NormInfty: IndexableMatrix { - fn norm_infty(&self) -> ::Real; -} - -/// Compute the trace of a matrix. -pub trait Trace: IndexableMatrix { - fn trace(&self) -> Self::Item; -} - -/// Compute the real part of a matrix. -pub trait Real: IndexableMatrix { - fn real(&self) -> ::Real; -} - -/// Compute the imaginary part of a matrix. -pub trait Imag: IndexableMatrix { - fn imag(&self) -> ::Real; -} - -/// Swap entries with another matrix. -pub trait Swap: IndexableMatrix { - fn swap(&mut self, other: &mut Self) -> crate::types::RlstResult<()>; -} - -/// Fill matrix by copying from another matrix. -pub trait Fill: IndexableMatrix { - fn fill(&mut self, other: &Self) -> crate::types::RlstResult<()>; -} - -/// Multiply entries with a scalar. -pub trait ScalarMult: IndexableMatrix { - fn scalar_mult(&mut self, scalar: Self::Item); -} - -/// Compute self -> alpha * other + self. -pub trait MultSumInto: IndexableMatrix { - fn mult_sum_into( - &mut self, - other: &Self, - scalar: Self::Item, - ) -> crate::types::RlstResult<()>; -} - -/// Compute y-> alpha A x + y -pub trait MatMul: IndexableMatrix { - type Vec: IndexableVector; - - fn matmul(&self, alpha: Self::Item, x: &Self::Vec, beta: Self::Item, y: &mut Self::Vec); -} diff --git a/suitesparse-src/Cargo.toml b/suitesparse-src/Cargo.toml index c625b8ed..0f2268d9 100644 --- a/suitesparse-src/Cargo.toml +++ b/suitesparse-src/Cargo.toml @@ -13,10 +13,6 @@ links = "suitesparse" [lib] name = "rlst_suitesparse_src" -[dependencies] -rlst-blis-src = { path = "../blis-src"} -rlst-netlib-lapack-src = { path = "../netlib-lapack-src"} - [build-dependencies] bindgen = "0.65.*" cmake = "0.1.*" diff --git a/suitesparse-src/build.rs b/suitesparse-src/build.rs index 9caba2ac..419885d1 100644 --- a/suitesparse-src/build.rs +++ b/suitesparse-src/build.rs @@ -1,13 +1,10 @@ use cmake::Config; macro_rules! build_dep { - ($name:literal, $blas_lib:expr, $lapack_lib:expr) => {{ + ($name:literal) => {{ let out_dir = std::env::var("OUT_DIR").unwrap(); Config::new(format!("suitesparse/{}", $name)) .define("CMAKE_PREFIX_PATH", out_dir) - .define("LAPACK_LIBRARIES", $lapack_lib) - .define("BLAS_LIBRARIES", $blas_lib) - .define("BLA_STATIC", "TRUE") .build() }}; } @@ -16,16 +13,17 @@ fn main() { // These are only needed since the build scripts for Suitesparse insist // on also building shared library versions that require resolution of // the Blas/Lapack symbol names. - let blas_lib = std::env::var("DEP_BLIS_ROOT").unwrap() + "/lib/libblis.a"; - let lapack_lib = std::env::var("DEP_NETLIBLAPACK_ROOT").unwrap() + "/lib/liblapack.a"; - - let suitesparse = build_dep!("SuiteSparse_config", &blas_lib, &lapack_lib); - let _amd = build_dep!("AMD", &blas_lib, &lapack_lib); - let _camd = build_dep!("CAMD", &blas_lib, &lapack_lib); - let _colamd = build_dep!("COLAMD", &blas_lib, &lapack_lib); - let _ccolamd = build_dep!("CCOLAMD", &blas_lib, &lapack_lib); - let _cholmod = build_dep!("CHOLMOD", &blas_lib, &lapack_lib); - let _umfpack = build_dep!("UMFPACK", &blas_lib, &lapack_lib); + + let out_dir = std::env::var("OUT_DIR").unwrap(); + println!("cargo:warning=Out Dir {}", out_dir); + + let suitesparse = build_dep!("SuiteSparse_config"); + let _amd = build_dep!("AMD"); + let _camd = build_dep!("CAMD"); + let _colamd = build_dep!("COLAMD"); + let _ccolamd = build_dep!("CCOLAMD"); + let _cholmod = build_dep!("CHOLMOD"); + let _umfpack = build_dep!("UMFPACK"); println!("cargo:root={}", std::env::var("OUT_DIR").unwrap()); diff --git a/umfpack/Cargo.toml b/umfpack/Cargo.toml index 99226b0d..b4d491b2 100644 --- a/umfpack/Cargo.toml +++ b/umfpack/Cargo.toml @@ -10,8 +10,7 @@ edition = "2021" [dependencies] rlst-suitesparse-src = {path = "../suitesparse-src"} -rlst-blis-src = { path = "../blis-src"} -rlst-netlib-lapack-src = { path = "../netlib-lapack-src"} + [lib] name = "rlst_umfpack" diff --git a/umfpack/examples/umfpack_simple.rs b/umfpack/examples/umfpack_simple.rs deleted file mode 100644 index 2e23c69d..00000000 --- a/umfpack/examples/umfpack_simple.rs +++ /dev/null @@ -1,61 +0,0 @@ -//! Translation of umfpack_simple.c from the UMFPACK User Guide -use std::ffi::c_void; - -use rlst_umfpack as umfpack; - -pub fn main() { - let n = 5; - - let ai = [0, 1, 0, 2, 4, 1, 2, 3, 4, 2, 1, 4]; - let ap = [0, 2, 5, 9, 10, 12]; - let ax = [2.0, 3.0, 3.0, -1.0, 4.0, 4.0, -3.0, 1.0, 2.0, 2.0, 6.0, 1.0]; - let b = [8.0, 45.0, -3.0, 3.0, 19.0]; - let mut symbolic = std::ptr::null_mut::(); - let mut numeric = std::ptr::null_mut::(); - let mut x = vec![0.0; 5]; - - let _ = unsafe { - umfpack::umfpack_di_symbolic( - n, - n, - ap.as_ptr(), - ai.as_ptr(), - ax.as_ptr(), - &mut symbolic, - std::ptr::null(), - std::ptr::null_mut(), - ) - }; - - let _ = unsafe { - umfpack::umfpack_di_numeric( - ap.as_ptr(), - ai.as_ptr(), - ax.as_ptr(), - symbolic, - &mut numeric, - std::ptr::null(), - std::ptr::null_mut(), - ) - }; - - unsafe { umfpack::umfpack_di_free_symbolic(&mut symbolic) }; - - let _ = unsafe { - umfpack::umfpack_di_solve( - umfpack::UMFPACK_A as i32, - ap.as_ptr(), - ai.as_ptr(), - ax.as_ptr(), - x.as_mut_ptr(), - b.as_ptr(), - numeric, - std::ptr::null(), - std::ptr::null_mut(), - ) - }; - - unsafe { umfpack::umfpack_di_free_numeric(&mut numeric) }; - - println!("Solution: {:#?}", x); -} diff --git a/umfpack/src/lib.rs b/umfpack/src/lib.rs index e02ecba1..7ebb12e4 100644 --- a/umfpack/src/lib.rs +++ b/umfpack/src/lib.rs @@ -3,8 +3,6 @@ #![allow(non_snake_case)] #![cfg_attr(feature = "strict", deny(warnings))] -extern crate rlst_blis_src; -extern crate rlst_netlib_lapack_src; extern crate rlst_suitesparse_src; include!(concat!(env!("OUT_DIR"), "/bindings.rs"));