diff --git a/.github/workflows/nalgebra-ci-build.yml b/.github/workflows/nalgebra-ci-build.yml index fd3ec273..3a56df13 100644 --- a/.github/workflows/nalgebra-ci-build.yml +++ b/.github/workflows/nalgebra-ci-build.yml @@ -36,14 +36,20 @@ jobs: run: cargo build; - name: Build --features serde-serialize run: cargo build --features serde-serialize - - name: Build --all-features - run: cargo build --all-features; - - name: Build nalgebra-glm - run: cargo build -p nalgebra-glm --all-features; - name: Build nalgebra-lapack run: cd nalgebra-lapack; cargo build; - name: Build nalgebra-sparse run: cd nalgebra-sparse; cargo build; + # Run this on it’s own job because it alone takes a lot of time. + # So it’s best to let it run in parallel to the other jobs. + build-nalgebra-all-features: + runs-on: ubuntu-latest + steps: + # Needed because the --all-features build which enables cuda support. + - uses: Jimver/cuda-toolkit@v0.2.4 + - uses: actions/checkout@v2 + - run: cargo build --all-features; + - run: cargo build -p nalgebra-glm --all-features; test-nalgebra: runs-on: ubuntu-latest # env: @@ -65,10 +71,10 @@ jobs: - name: test nalgebra-sparse # Manifest-path is necessary because cargo otherwise won't correctly forward features # We increase number of proptest cases to hopefully catch more potential bugs - run: PROPTEST_CASES=10000 cargo test --manifest-path=nalgebra-sparse/Cargo.toml --features compare,proptest-support + run: PROPTEST_CASES=10000 cargo test --manifest-path=nalgebra-sparse/Cargo.toml --features compare,proptest-support,io - name: test nalgebra-sparse (slow tests) # Unfortunately, the "slow-tests" take so much time that we need to run them with --release - run: PROPTEST_CASES=10000 cargo test --release --manifest-path=nalgebra-sparse/Cargo.toml --features compare,proptest-support,slow-tests slow + run: PROPTEST_CASES=10000 cargo test --release --manifest-path=nalgebra-sparse/Cargo.toml --features compare,proptest-support,io,slow-tests slow test-nalgebra-macros: runs-on: ubuntu-latest steps: @@ -106,3 +112,20 @@ jobs: run: xargo build --verbose --no-default-features --features alloc --target=x86_64-unknown-linux-gnu; - name: build thumbv7em-none-eabihf run: xargo build --verbose --no-default-features --target=thumbv7em-none-eabihf; + - name: build x86_64-unknown-linux-gnu nalgebra-glm + run: xargo build --verbose --no-default-features -p nalgebra-glm --target=x86_64-unknown-linux-gnu; + - name: build thumbv7em-none-eabihf nalgebra-glm + run: xargo build --verbose --no-default-features -p nalgebra-glm --target=thumbv7em-none-eabihf; + build-cuda: + runs-on: ubuntu-latest + steps: + - uses: Jimver/cuda-toolkit@v0.2.4 + - name: Install nightly-2021-12-04 + uses: actions-rs/toolchain@v1 + with: + toolchain: nightly-2021-12-04 + override: true + - uses: actions/checkout@v2 + - run: rustup target add nvptx64-nvidia-cuda + - run: cargo build --no-default-features --features cuda + - run: cargo build --no-default-features --features cuda --target=nvptx64-nvidia-cuda \ No newline at end of file diff --git a/CHANGELOG.md b/CHANGELOG.md index a55a6a5f..88006c2e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,52 @@ documented here. This project adheres to [Semantic Versioning](https://semver.org/). +## [0.30.0] (02 Jan. 2022) + +### Breaking changes +- The `Dim` trait is now marked as unsafe. +- The `Matrix::pow` and `Matrix::pow_mut` methods only allow positive integer exponents now. To compute negative + exponents, the user is free to invert the matrix before calling `pow` with the exponent’s absolute value. + +### Modified +- Use more concise debug impls for matrices and geometric transformation types. +- The singular values computed by the SVD are now sorted in increasing order by default. Use `SVD::new_unordered` + instead to reproduce the older behavior without the sorting overhead. +- The `UnitDualQuaternion::sclerp` method will no longer panic when given two equal rotations. +- The `Matrix::select_rows` and `Matrix::select_columns` methods no longer require the matrix components to implement + the trait `Zero`. +- The `Matrix::pow` and `Matrix::pow_mut` methods will now also work with integer matrices. + +### Added +- Added the conversion trait `From>` and method `from_vec_storage` for `RowDVector`. +- Added implementation of `From` and `Into` for converting between `nalgebra` types and types from + `glam 0.18`. These can be enabled by enabling the `convert-glam018` cargo features. +- Added the methods `Matrix::product`, `::row_product`, `::row_product_tr`, and `::column_product` to compute the + product of the components, rows, or columns, of a single matrix or vector. +- The `Default` trait is now implemented for most geometric types: `Point`, `Isometry`, `Rotation`, `Similarity`, + `Transform`, `UnitComplex`, and `UnitQuaternion`. +- Added the `Scale` geometric type for representing non-uniform scaling. +- Added `Cholesky::new_with_substitute` that will replace diagonal elements by a given constant whenever `Cholesky` + meets a non-definite-positiveness. +- Re-added the conversion from a vector/matrix slice to a static array. +- Added the `cuda` feature that enables the support of [rust-cuda](https://github.com/Rust-GPU/Rust-CUDA) for using + `nalgebra` features with CUDA kernels written in Rust. +- Added special-cases implementations for the 2x2 and 3x3 SVDs for better accuracy and performances. +- Added the methods `Matrix::polar`, `Matrix::try_polar`, and `SVD::to_polar` to compute the polar decomposition of + a matrix, based on its SVD. +- `nalgebra-sparse`: provide constructors for unsorted but otherwise valid data using the CSR format. +- `nalgebra-sparse`: added reading MatrixMarked data files to a sparse `CooMatrix`. + +### Fixed +- Fixed a potential unsoundness with `matrix.get(i)` and `matrix.get_mut(i)` where `i` is an `usize`, and `matrix` + is a matrix slice with non-default strides. +- Fixed potential unsoundness with `vector.perp` where `vector` isn’t actually a 2D vector as expected. +- Fixed linkage issue with `nalgebra-lapack`: the user of `nalgebra-lapack` no longer have to add + `extern crate lapack-src` to their `main.rs`. +- Fixed the `no-std` build of `nalgebra-glm`. +- Fix the `pow` and `pow_mut` functions (the result was incorrect for some exponent values). + + ## [0.29.0] ### Breaking changes - We updated to the version 0.6 of `simba`. This means that the trait bounds `T: na::RealField`, `na::ComplexField`, @@ -242,21 +288,21 @@ for details about this switch and its benefits. geometric types. ### Modified * Use of traits like `alga::general::{RealField, ComplexField}` have now been replaced by - `simba::scalar::{RealField, ComplexField}`. + `simba::scalar::{RealField, ComplexField}`. * The implementation of traits from the __alga__ crate (and well as the dependency to _alga__) are now omitted unless the `alga` cargo feature is activated. ### Removed * The `Neg` unary operator is no longer implemented for `UnitComplex` and `UnitQuaternion`. This caused hard-to-track errors when we mistakenly write, e.g., `-q * v` instead of `-(q * v)`. * The `na::convert_unchecked` is no longer marked as unsafe. - + ## [0.20.0] ### Added * `cholesky.rank_one_update(...)` which performs a rank-one update on the cholesky decomposition of a matrix. * `From<&Matrix>` is now implemented for matrix slices. * `.try_set_magnitude(...)` which sets the magnitude of a vector, while keeping its direction. * Implementations of `From` and `Into` for the conversion between matrix slices and standard (`&[N]` `&mut [N]`) slices. - + ### Modified * We started some major changes in order to allow non-Copy types to be used as scalar types inside of matrices/vectors. @@ -264,16 +310,16 @@ for details about this switch and its benefits. ### Added * `.remove_rows_at` and `remove_columns_at` which removes a set of rows or columns (specified by indices) from a matrix. * Several formatting traits have been implemented for all matrices/vectors: `LowerExp`, `UpperExp`, `Octal`, `LowerHex`, - `UpperHex`, `Binary`, `Pointer`. + `UpperHex`, `Binary`, `Pointer`. * `UnitQuaternion::quaternions_mean(...)` which computes the mean rotation of a set of unit quaternions. This implements - the algorithm from _Oshman, Yaakov, and Avishy Carmi, "Attitude estimation from vector observations using a genetic-algorithm-embedded quaternion particle filter." + the algorithm from _Oshman, Yaakov, and Avishy Carmi, "Attitude estimation from vector observations using a genetic-algorithm-embedded quaternion particle filter." ### Modified * It is now possible to get the `min/max` element of unsigned integer matrices. ### Added to nalgebra-glm * Some infinite and reversed perspectives: `::infinite_perspective_rh_no`, `::infinite_perspective_rh_zo`, - `::reversed_perspective_rh_zo`, and `::reversed_infinite_perspective_rh_zo`. + `::reversed_perspective_rh_zo`, and `::reversed_infinite_perspective_rh_zo`. ## [0.18.0] This release adds full complex number support to nalgebra. This includes all common vector/matrix operations as well @@ -287,12 +333,12 @@ as matrix decomposition. This excludes geometric type (like `Isometry`, `Rotatio * Add `.left_div, .right_div` for quaternions. * Add `.renormalize` to `Unit<...>` and `Rotation3` to correct potential drift due to repeated operations. Those drifts could cause them not to be pure rotations anymore. - + #### Convolution * `.convolve_full(kernel)` returns the convolution of `self` by `kernel`. * `.convolve_valid(kernel)` returns the convolution of `self` by `kernel` after removal of all the elements relying on zero-padding. * `.convolve_same(kernel)` returns the convolution of `self` by `kernel` with a result of the same size as `self`. - + #### Complex number support * Add the `::from_matrix` constructor too all rotation types to extract a rotation from a raw matrix. * Add the `::from_matrix_eps` constructor too all rotation types to extract a rotation from a raw matrix. This takes @@ -329,22 +375,22 @@ Note that all the other BLAS operation will continue to work for all fields, inc * Add `.slerp` and `.try_slerp` to unit vectors. * Add `.lerp` to vectors. * Add `.to_projective` and `.as_projective` to `Perspective3` and `Orthographic3` in order to - use them as `Projective3` structures. + use them as `Projective3` structures. * Add `From/Into` impls to allow the conversion of any transformation type to a matrix. * Add `Into` impls to convert a matrix slice into an owned matrix. * Add `Point*::from_slice` to create a point from a slice. * Add `.map_with_location` to matrices to apply a map which passes the component indices to the user-defined closure alongside - the component itself. + the component itself. * Add impl `From` for `Point`. * Add impl `From` for `Quaternion`. * Add impl `From` for `Translation`. * Add the `::from_vec` constructor to construct a matrix from a `Vec` (a `DMatrix` will reuse the original `Vec` - as-is for its storage). + as-is for its storage). * Add `.to_homogeneous` to square matrices (and with dimensions higher than 1x1). This will increase their number of row - and columns by 1. The new column and row are filled with 0, except for the diagonal element which is set to 1. + and columns by 1. The new column and row are filled with 0, except for the diagonal element which is set to 1. * Implement `Extend` for matrices with a dynamic storage. The provided `Vec` is assumed to represent a column-major - matrix with the same number of rows as the one being extended. This will effectively append new columns on the right of - the matrix being extended. + matrix with the same number of rows as the one being extended. This will effectively append new columns on the right of + the matrix being extended. * Implement `Extend` for vectors with a dynamic storage. This will concatenate the vector with the given `Vec`. * Implement `Extend>` for matrices with dynamic storage. This will concatenate the columns of both matrices. * Implement `Into` for the `MatrixVec` storage. @@ -353,10 +399,10 @@ Note that all the other BLAS operation will continue to work for all fields, inc ### Modified * The orthographic projection no longer require that `bottom < top`, that `left < right`, and that `znear < zfar`. The - only restriction now ith that they must not be equal (in which case the projection would be singular). + only restriction now ith that they must not be equal (in which case the projection would be singular). * The `Point::from_coordinates` methods is deprecated. Use `Point::from` instead. * The `.transform_point` and `.transform_vector` methods are now inherent methods for matrices so that the user does not have to - explicitly import the `Transform` trait from the alga crate. + explicitly import the `Transform` trait from the alga crate. * Renamed the matrix storage types: `MatrixArray` -> `ArrayStorage` and `MatrixVec` -> `VecStorage`. * Renamed `.unwrap()` to `.into_inner()` for geometric types that wrap another type. This is for the case of `Unit`, `Transform`, `Orthographic3`, `Perspective3`, `Rotation`. @@ -552,8 +598,8 @@ overview of all the added/modified features. This version is a major rewrite of the library. Major changes are: * Algebraic traits are now defined by the [alga](https://crates.io/crates/alga) crate. - All other mathematical traits, except `Axpy` have been removed from - **nalgebra**. + All other mathematical traits, except `Axpy` have been removed from + **nalgebra**. * Methods are now preferred to free functions because they do not require any trait to be used anymore. * Most algebraic entities can be parametrized by type-level integers @@ -746,11 +792,11 @@ you [there](https://users.nphysics.org)! ### Modified * Vectors are now multipliable with isometries. This will result into a pure rotation (this is how - vectors differ from point semantically: they design directions, so they are not translatable). + vectors differ from point semantically: they design directions, so they are not translatable). * `{Isometry3, Rotation3}::look_at` reimplemented and renamed to `::look_at_rh` and `::look_at_lh` to agree - with the computer graphics community (in particular, the GLM library). Use the `::look_at_rh` - variant to build a view matrix that - may be successfully used with `Persp` and `Ortho`. + with the computer graphics community (in particular, the GLM library). Use the `::look_at_rh` + variant to build a view matrix that + may be successfully used with `Persp` and `Ortho`. * The old `{Isometry3, Rotation3}::look_at` implementations are now called `::new_observer_frame`. * Rename every `fov` on `Persp` to `fovy`. * Fixed the perspective and orthographic projection matrices. diff --git a/Cargo.toml b/Cargo.toml index 2b4a7487..8294398e 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "nalgebra" -version = "0.29.0" +version = "0.30.0" authors = [ "Sébastien Crozet " ] description = "General-purpose linear algebra library with transformations and statically-sized or dynamically-sized matrices." @@ -32,6 +32,7 @@ compare = [ "matrixcompare-core" ] libm = [ "simba/libm" ] libm-force = [ "simba/libm_force" ] macros = [ "nalgebra-macros" ] +cuda = [ "cust", "simba/cuda" ] # Conversion convert-mint = [ "mint" ] @@ -41,6 +42,7 @@ convert-glam014 = [ "glam014" ] convert-glam015 = [ "glam015" ] convert-glam016 = [ "glam016" ] convert-glam017 = [ "glam017" ] +convert-glam018 = [ "glam018" ] # Serialization ## To use serde in a #[no-std] environment, enable the @@ -72,7 +74,7 @@ num-traits = { version = "0.2", default-features = false } num-complex = { version = "0.4", default-features = false } num-rational = { version = "0.4", default-features = false } approx = { version = "0.5", default-features = false } -simba = { version = "0.6", default-features = false } +simba = { version = "0.7", default-features = false } alga = { version = "0.9", default-features = false, optional = true } rand_distr = { version = "0.4", default-features = false, optional = true } matrixmultiply = { version = "0.3", optional = true } @@ -85,12 +87,16 @@ pest = { version = "2", optional = true } pest_derive = { version = "2", optional = true } bytemuck = { version = "1.5", optional = true } matrixcompare-core = { version = "0.1", optional = true } -proptest = { version = "1", optional = true, default-features = false, features = ["std"] } +proptest = { version = "1", optional = true, default-features = false, features = ["std"] } glam013 = { package = "glam", version = "0.13", optional = true } glam014 = { package = "glam", version = "0.14", optional = true } glam015 = { package = "glam", version = "0.15", optional = true } glam016 = { package = "glam", version = "0.16", optional = true } glam017 = { package = "glam", version = "0.17", optional = true } +glam018 = { package = "glam", version = "0.18", optional = true } + +[target.'cfg(not(target_os = "cuda"))'.dependencies] +cust = { version = "0.2", optional = true } [dev-dependencies] @@ -127,3 +133,4 @@ lto = true [package.metadata.docs.rs] # Enable certain features when building docs for docs.rs features = [ "proptest-support", "compare", "macros", "rand" ] + diff --git a/Makefile b/Makefile deleted file mode 100644 index 0af9f3f2..00000000 --- a/Makefile +++ /dev/null @@ -1,12 +0,0 @@ -all: - cargo test --features "debug arbitrary serde-serialize abomonation-serialize compare" - # cargo check --features "debug arbitrary serde-serialize" - -doc: - cargo doc --no-deps --features "debug arbitrary serde-serialize abomonation" - -bench: - cargo bench - -test: - cargo test --features "debug arbitrary serde-serialize abomonation-serialize compare" diff --git a/README.md b/README.md index bed6183d..fa1e0904 100644 --- a/README.md +++ b/README.md @@ -9,7 +9,7 @@ Build status - crates.io + crates.io @@ -30,10 +30,18 @@ ----- -## Platinum sponsors -Rapier is supported by: +## Acknowledgements +nalgebra is supported by our **platinum** sponsors:

- +

+ +And our gold sponsors: + +

+ + + +

\ No newline at end of file diff --git a/benches/linalg/svd.rs b/benches/linalg/svd.rs index b84f60d6..c95e44de 100644 --- a/benches/linalg/svd.rs +++ b/benches/linalg/svd.rs @@ -1,4 +1,18 @@ -use na::{Matrix4, SVD}; +use na::{Matrix2, Matrix3, Matrix4, SVD}; + +fn svd_decompose_2x2_f32(bh: &mut criterion::Criterion) { + let m = Matrix2::::new_random(); + bh.bench_function("svd_decompose_2x2", move |bh| { + bh.iter(|| std::hint::black_box(SVD::new_unordered(m.clone(), true, true))) + }); +} + +fn svd_decompose_3x3_f32(bh: &mut criterion::Criterion) { + let m = Matrix3::::new_random(); + bh.bench_function("svd_decompose_3x3", move |bh| { + bh.iter(|| std::hint::black_box(SVD::new_unordered(m.clone(), true, true))) + }); +} fn svd_decompose_4x4(bh: &mut criterion::Criterion) { let m = Matrix4::::new_random(); @@ -114,6 +128,8 @@ fn pseudo_inverse_200x200(bh: &mut criterion::Criterion) { criterion_group!( svd, + svd_decompose_2x2_f32, + svd_decompose_3x3_f32, svd_decompose_4x4, svd_decompose_10x10, svd_decompose_100x100, diff --git a/examples/cargo/Cargo.toml b/examples/cargo/Cargo.toml index 454c88b0..501ae23e 100644 --- a/examples/cargo/Cargo.toml +++ b/examples/cargo/Cargo.toml @@ -4,7 +4,7 @@ version = "0.0.0" authors = [ "You" ] [dependencies] -nalgebra = "0.29.0" +nalgebra = "0.30.0" [[bin]] name = "example" diff --git a/examples/mvp.rs b/examples/mvp.rs index 3ccb2549..add7ac8f 100644 --- a/examples/mvp.rs +++ b/examples/mvp.rs @@ -3,6 +3,7 @@ extern crate nalgebra as na; use na::{Isometry3, Perspective3, Point3, Vector3}; +use std::f32::consts; fn main() { // Our object is translated along the x axis. @@ -15,7 +16,7 @@ fn main() { let view = Isometry3::look_at_rh(&eye, &target, &Vector3::y()); // A perspective projection. - let projection = Perspective3::new(16.0 / 9.0, 3.14 / 2.0, 1.0, 1000.0); + let projection = Perspective3::new(16.0 / 9.0, consts::PI / 2.0, 1.0, 1000.0); // The combination of the model with the view is still an isometry. let model_view = view * model; diff --git a/examples/raw_pointer.rs b/examples/raw_pointer.rs index d8c84917..bb3dccd5 100644 --- a/examples/raw_pointer.rs +++ b/examples/raw_pointer.rs @@ -19,6 +19,7 @@ fn main() { /* Then pass the raw pointers to some graphics API. */ + #[allow(clippy::float_cmp)] unsafe { assert_eq!(*v_pointer, 1.0); assert_eq!(*v_pointer.offset(1), 0.0); diff --git a/examples/screen_to_view_coords.rs b/examples/screen_to_view_coords.rs index b4d3d255..e16baa98 100644 --- a/examples/screen_to_view_coords.rs +++ b/examples/screen_to_view_coords.rs @@ -3,9 +3,10 @@ extern crate nalgebra as na; use na::{Perspective3, Point2, Point3, Unit}; +use std::f32::consts; fn main() { - let projection = Perspective3::new(800.0 / 600.0, 3.14 / 2.0, 1.0, 1000.0); + let projection = Perspective3::new(800.0 / 600.0, consts::PI / 2.0, 1.0, 1000.0); let screen_point = Point2::new(10.0f32, 20.0); // Compute two points in clip-space. diff --git a/examples/transform_conversion.rs b/examples/transform_conversion.rs index 71660a07..f284e4b3 100644 --- a/examples/transform_conversion.rs +++ b/examples/transform_conversion.rs @@ -1,6 +1,7 @@ extern crate nalgebra as na; use na::{Isometry2, Similarity2, Vector2}; +use std::f32::consts; fn main() { // Isometry -> Similarity conversion always succeeds. @@ -8,8 +9,8 @@ fn main() { let _: Similarity2 = na::convert(iso); // Similarity -> Isometry conversion fails if the scaling factor is not 1.0. - let sim_without_scaling = Similarity2::new(Vector2::new(1.0f32, 2.0), 3.14, 1.0); - let sim_with_scaling = Similarity2::new(Vector2::new(1.0f32, 2.0), 3.14, 2.0); + let sim_without_scaling = Similarity2::new(Vector2::new(1.0f32, 2.0), consts::PI, 1.0); + let sim_with_scaling = Similarity2::new(Vector2::new(1.0f32, 2.0), consts::PI, 2.0); let iso_success: Option> = na::try_convert(sim_without_scaling); let iso_fail: Option> = na::try_convert(sim_with_scaling); diff --git a/examples/transform_matrix4.rs b/examples/transform_matrix4.rs index 74649d74..b7228e51 100644 --- a/examples/transform_matrix4.rs +++ b/examples/transform_matrix4.rs @@ -3,6 +3,7 @@ extern crate approx; extern crate nalgebra as na; use na::{Matrix4, Point3, Vector3}; +use std::f32::consts; fn main() { // Create a uniform scaling matrix with scaling factor 2. @@ -28,7 +29,7 @@ fn main() { ); // Create rotation. - let rot = Matrix4::from_scaled_axis(&Vector3::x() * 3.14); + let rot = Matrix4::from_scaled_axis(Vector3::x() * consts::PI); let rot_then_m = m * rot; // Right-multiplication is equivalent to prepending `rot` to `m`. let m_then_rot = rot * m; // Left-multiplication is equivalent to appending `rot` to `m`. diff --git a/examples/transformation_pointer.rs b/examples/transformation_pointer.rs index a39c5568..f2074906 100644 --- a/examples/transformation_pointer.rs +++ b/examples/transformation_pointer.rs @@ -12,6 +12,7 @@ fn main() { /* Then pass the raw pointer to some graphics API. */ + #[allow(clippy::float_cmp)] unsafe { assert_eq!(*iso_pointer, 1.0); assert_eq!(*iso_pointer.offset(5), 1.0); diff --git a/examples/unit_wrapper.rs b/examples/unit_wrapper.rs index 229fad1a..81339ddf 100644 --- a/examples/unit_wrapper.rs +++ b/examples/unit_wrapper.rs @@ -1,3 +1,4 @@ +#![allow(clippy::float_cmp)] extern crate nalgebra as na; use na::{Unit, Vector3}; diff --git a/nalgebra-glm/Cargo.toml b/nalgebra-glm/Cargo.toml index 5c2ae2e1..287bb8c7 100644 --- a/nalgebra-glm/Cargo.toml +++ b/nalgebra-glm/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "nalgebra-glm" -version = "0.15.0" +version = "0.16.0" authors = ["sebcrozet "] description = "A computer-graphics oriented API for nalgebra, inspired by the C++ GLM library." @@ -22,9 +22,20 @@ std = [ "nalgebra/std", "simba/std" ] arbitrary = [ "nalgebra/arbitrary" ] serde-serialize = [ "nalgebra/serde-serialize-no-std" ] abomonation-serialize = [ "nalgebra/abomonation-serialize" ] +cuda = [ "nalgebra/cuda" ] + +# Conversion +convert-mint = [ "nalgebra/mint" ] +convert-bytemuck = [ "nalgebra/bytemuck" ] +convert-glam013 = [ "nalgebra/glam013" ] +convert-glam014 = [ "nalgebra/glam014" ] +convert-glam015 = [ "nalgebra/glam015" ] +convert-glam016 = [ "nalgebra/glam016" ] +convert-glam017 = [ "nalgebra/glam017" ] +convert-glam018 = [ "nalgebra/glam018" ] [dependencies] num-traits = { version = "0.2", default-features = false } approx = { version = "0.5", default-features = false } -simba = { version = "0.6", default-features = false } -nalgebra = { path = "..", version = "0.29", default-features = false } +simba = { version = "0.7", default-features = false } +nalgebra = { path = "..", version = "0.30", default-features = false } diff --git a/nalgebra-glm/src/constructors.rs b/nalgebra-glm/src/constructors.rs index e998dd23..8afdaa20 100644 --- a/nalgebra-glm/src/constructors.rs +++ b/nalgebra-glm/src/constructors.rs @@ -75,6 +75,21 @@ pub fn mat2x4(m11: T, m12: T, m13: T, m14: T, } /// Create a new 3x3 matrix. +/// +/// # Example +/// ``` +/// # use nalgebra_glm::mat3; +/// let m = mat3( +/// 1.0, 2.0, 3.0, +/// 4.0, 5.0, 6.0, +/// 7.0, 8.0, 9.0 +/// ); +/// assert!( +/// m.m11 == 1.0 && m.m12 == 2.0 && m.m13 == 3.0 && +/// m.m21 == 4.0 && m.m22 == 5.0 && m.m23 == 6.0 && +/// m.m31 == 7.0 && m.m32 == 8.0 && m.m33 == 9.0 +/// ); +/// ``` #[rustfmt::skip] pub fn mat3(m11: T, m12: T, m13: T, m21: T, m22: T, m23: T, @@ -105,8 +120,8 @@ pub fn mat3x3(m11: T, m12: T, m13: T, m31: T, m32: T, m33: T) -> TMat3 { TMat::::new( m11, m12, m13, - m31, m32, m33, m21, m22, m23, + m31, m32, m33, ) } diff --git a/nalgebra-glm/src/traits.rs b/nalgebra-glm/src/traits.rs index a09a95f2..ccdb4a02 100644 --- a/nalgebra-glm/src/traits.rs +++ b/nalgebra-glm/src/traits.rs @@ -1,9 +1,9 @@ use approx::AbsDiffEq; use num::{Bounded, Signed}; +use core::cmp::PartialOrd; use na::Scalar; use simba::scalar::{ClosedAdd, ClosedMul, ClosedSub, RealField}; -use std::cmp::PartialOrd; /// A number that can either be an integer or a float. pub trait Number: diff --git a/nalgebra-glm/src/trigonometric.rs b/nalgebra-glm/src/trigonometric.rs index 90227a8d..dce45738 100644 --- a/nalgebra-glm/src/trigonometric.rs +++ b/nalgebra-glm/src/trigonometric.rs @@ -53,7 +53,7 @@ pub fn degrees(radians: &TVec) -> TVec(degrees: &TVec) -> TVec { degrees.map(|e| e * T::pi() / na::convert(180.0)) } diff --git a/nalgebra-lapack/Cargo.toml b/nalgebra-lapack/Cargo.toml index b7f296f0..7cdf9f78 100644 --- a/nalgebra-lapack/Cargo.toml +++ b/nalgebra-lapack/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "nalgebra-lapack" -version = "0.20.0" +version = "0.21.0" authors = [ "Sébastien Crozet ", "Andrew Straw " ] description = "Matrix decompositions using nalgebra matrices and Lapack bindings." @@ -29,17 +29,17 @@ accelerate = ["lapack-src/accelerate"] intel-mkl = ["lapack-src/intel-mkl"] [dependencies] -nalgebra = { version = "0.29", path = ".." } +nalgebra = { version = "0.30", path = ".." } num-traits = "0.2" num-complex = { version = "0.4", default-features = false } -simba = "0.5" +simba = "0.7" serde = { version = "1.0", features = [ "derive" ], optional = true } lapack = { version = "0.19", default-features = false } lapack-src = { version = "0.8", default-features = false } # clippy = "*" [dev-dependencies] -nalgebra = { version = "0.29", features = [ "arbitrary", "rand" ], path = ".." } +nalgebra = { version = "0.30", features = [ "arbitrary", "rand" ], path = ".." } proptest = { version = "1", default-features = false, features = ["std"] } quickcheck = "1" approx = "0.5" diff --git a/nalgebra-lapack/src/eigen.rs b/nalgebra-lapack/src/eigen.rs index f6628bfe..8eab62d8 100644 --- a/nalgebra-lapack/src/eigen.rs +++ b/nalgebra-lapack/src/eigen.rs @@ -294,7 +294,7 @@ where let mut res = Matrix::zeros_generic(nrows, Const::<1>); for i in 0..res.len() { - res[i] = Complex::new(wr[i], wi[i]); + res[i] = Complex::new(wr[i].clone(), wi[i].clone()); } res @@ -306,7 +306,7 @@ where pub fn determinant(&self) -> T { let mut det = T::one(); for e in self.eigenvalues.iter() { - det *= *e; + det *= e.clone(); } det diff --git a/nalgebra-lapack/src/lib.rs b/nalgebra-lapack/src/lib.rs index 84fa03fa..9cf0d73d 100644 --- a/nalgebra-lapack/src/lib.rs +++ b/nalgebra-lapack/src/lib.rs @@ -73,6 +73,9 @@ html_root_url = "https://nalgebra.org/rustdoc" )] +extern crate lapack; +extern crate lapack_src; + extern crate nalgebra as na; extern crate num_traits as num; diff --git a/nalgebra-lapack/src/schur.rs b/nalgebra-lapack/src/schur.rs index 13dfc05e..52f5dced 100644 --- a/nalgebra-lapack/src/schur.rs +++ b/nalgebra-lapack/src/schur.rs @@ -155,7 +155,7 @@ where let mut out = Matrix::zeros_generic(self.t.shape_generic().0, Const::<1>); for i in 0..out.len() { - out[i] = Complex::new(self.re[i], self.im[i]) + out[i] = Complex::new(self.re[i].clone(), self.im[i].clone()) } out diff --git a/nalgebra-lapack/src/symmetric_eigen.rs b/nalgebra-lapack/src/symmetric_eigen.rs index 8cbe63f8..0a4b09ce 100644 --- a/nalgebra-lapack/src/symmetric_eigen.rs +++ b/nalgebra-lapack/src/symmetric_eigen.rs @@ -140,7 +140,7 @@ where pub fn determinant(&self) -> T { let mut det = T::one(); for e in self.eigenvalues.iter() { - det *= *e; + det *= e.clone(); } det @@ -153,7 +153,7 @@ where pub fn recompose(&self) -> OMatrix { let mut u_t = self.eigenvectors.clone(); for i in 0..self.eigenvalues.len() { - let val = self.eigenvalues[i]; + let val = self.eigenvalues[i].clone(); u_t.column_mut(i).mul_assign(val); } u_t.transpose_mut(); diff --git a/nalgebra-lapack/tests/lib.rs b/nalgebra-lapack/tests/lib.rs index 973b6d17..6bf16e30 100644 --- a/nalgebra-lapack/tests/lib.rs +++ b/nalgebra-lapack/tests/lib.rs @@ -6,9 +6,6 @@ compile_error!("Tests must be run with `proptest-support`"); extern crate nalgebra as na; extern crate nalgebra_lapack as nl; -extern crate lapack; -extern crate lapack_src; - mod linalg; #[path = "../../tests/proptest/mod.rs"] mod proptest; diff --git a/nalgebra-macros/Cargo.toml b/nalgebra-macros/Cargo.toml index f77fc32e..bd37b689 100644 --- a/nalgebra-macros/Cargo.toml +++ b/nalgebra-macros/Cargo.toml @@ -21,5 +21,5 @@ quote = "1.0" proc-macro2 = "1.0" [dev-dependencies] -nalgebra = { version = "0.29.0", path = ".." } +nalgebra = { version = "0.30.0", path = ".." } trybuild = "1.0.42" diff --git a/nalgebra-sparse/Cargo.toml b/nalgebra-sparse/Cargo.toml index 76fb408a..6f7a7b4a 100644 --- a/nalgebra-sparse/Cargo.toml +++ b/nalgebra-sparse/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "nalgebra-sparse" -version = "0.5.0" +version = "0.6.0" authors = [ "Andreas Longva", "Sébastien Crozet " ] edition = "2018" description = "Sparse matrix computation based on nalgebra." @@ -16,19 +16,24 @@ license = "Apache-2.0" proptest-support = ["proptest", "nalgebra/proptest-support"] compare = [ "matrixcompare-core" ] +# Enable matrix market I/O +io = [ "pest", "pest_derive" ] + # Enable to enable running some tests that take a lot of time to run slow-tests = [] [dependencies] -nalgebra = { version="0.29", path = "../" } +nalgebra = { version="0.30", path = "../" } num-traits = { version = "0.2", default-features = false } proptest = { version = "1.0", optional = true } matrixcompare-core = { version = "0.1.0", optional = true } +pest = { version = "2", optional = true } +pest_derive = { version = "2", optional = true } [dev-dependencies] itertools = "0.10" matrixcompare = { version = "0.3.0", features = [ "proptest-support" ] } -nalgebra = { version="0.29", path = "../", features = ["compare"] } +nalgebra = { version="0.30", path = "../", features = ["compare"] } [package.metadata.docs.rs] # Enable certain features when building docs for docs.rs diff --git a/nalgebra-sparse/src/csr.rs b/nalgebra-sparse/src/csr.rs index c64be915..4324d18d 100644 --- a/nalgebra-sparse/src/csr.rs +++ b/nalgebra-sparse/src/csr.rs @@ -10,6 +10,7 @@ use crate::{SparseEntry, SparseEntryMut, SparseFormatError, SparseFormatErrorKin use nalgebra::Scalar; use num_traits::One; +use std::iter::FromIterator; use std::slice::{Iter, IterMut}; /// A CSR representation of a sparse matrix. @@ -170,6 +171,77 @@ impl CsrMatrix { Self::try_from_pattern_and_values(pattern, values) } + /// Try to construct a CSR matrix from raw CSR data with unsorted column indices. + /// + /// It is assumed that each row contains unique column indices that are in + /// bounds with respect to the number of columns in the matrix. If this is not the case, + /// an error is returned to indicate the failure. + /// + /// An error is returned if the data given does not conform to the CSR storage format + /// with the exception of having unsorted column indices and values. + /// See the documentation for [CsrMatrix](struct.CsrMatrix.html) for more information. + pub fn try_from_unsorted_csr_data( + num_rows: usize, + num_cols: usize, + row_offsets: Vec, + col_indices: Vec, + values: Vec, + ) -> Result + where + T: Scalar, + { + use SparsityPatternFormatError::*; + let count = col_indices.len(); + let mut p: Vec = (0..count).collect(); + + if col_indices.len() != values.len() { + return Err(SparseFormatError::from_kind_and_msg( + SparseFormatErrorKind::InvalidStructure, + "Number of values and column indices must be the same", + )); + } + + if row_offsets.len() == 0 { + return Err(SparseFormatError::from_kind_and_msg( + SparseFormatErrorKind::InvalidStructure, + "Number of offsets should be greater than 0", + )); + } + + for (index, &offset) in row_offsets[0..row_offsets.len() - 1].iter().enumerate() { + let next_offset = row_offsets[index + 1]; + if next_offset > count { + return Err(SparseFormatError::from_kind_and_msg( + SparseFormatErrorKind::InvalidStructure, + "No row offset should be greater than the number of column indices", + )); + } + if offset > next_offset { + return Err(NonmonotonicOffsets).map_err(pattern_format_error_to_csr_error); + } + p[offset..next_offset].sort_by(|a, b| { + let x = &col_indices[*a]; + let y = &col_indices[*b]; + x.partial_cmp(y).unwrap() + }); + } + + // permute indices + let sorted_col_indices: Vec = + Vec::from_iter((p.iter().map(|i| &col_indices[*i])).cloned()); + + // permute values + let sorted_values: Vec = Vec::from_iter((p.iter().map(|i| &values[*i])).cloned()); + + return Self::try_from_csr_data( + num_rows, + num_cols, + row_offsets, + sorted_col_indices, + sorted_values, + ); + } + /// Try to construct a CSR matrix from a sparsity pattern and associated non-zero values. /// /// Returns an error if the number of values does not match the number of minor indices diff --git a/nalgebra-sparse/src/io/matrix_market.pest b/nalgebra-sparse/src/io/matrix_market.pest new file mode 100644 index 00000000..aeeba406 --- /dev/null +++ b/nalgebra-sparse/src/io/matrix_market.pest @@ -0,0 +1,53 @@ +WHITESPACE = _{ " "|"\t" } + +// + +Sparsity = {^"coordinate" | ^"array"} +DataType = {^"real" | ^"complex" | ^"pattern" | ^"integer" } +StorageScheme = {^"symmetric" | ^"general" | ^"skew-symmetric" | ^"hermitian"} +// Only consider matrices here. +Header = { ^"%%matrixmarket matrix" ~ Sparsity ~ DataType ~ StorageScheme } + +// + +Comments = _{ "%" ~ (!NEWLINE ~ ANY)* } + +// + +Dimension = @{ ASCII_DIGIT+ } +SparseShape = { Dimension ~ Dimension ~ Dimension} +DenseShape = { Dimension ~ Dimension} +Shape = {SparseShape | DenseShape } + +// + +// grammar from https://doc.rust-lang.org/std/primitive.f64.html#grammar + +Sign = {("+" | "-")} +Exp = @{ ^"e" ~ Sign? ~ ASCII_DIGIT+} +Number = @{ ((ASCII_DIGIT+ ~ "." ~ ASCII_DIGIT*) | (ASCII_DIGIT* ~ "." ~ASCII_DIGIT+) | ASCII_DIGIT+ ) ~ Exp? } +Real = @{ Sign? ~ ("inf" | "NaN" | Number) } + + +SparseReal = {Dimension~ Dimension~ Real } +SparseComplex = {Dimension ~ Dimension ~ Real ~ Real} +SparsePattern = {Dimension ~ Dimension} + +DenseReal = {Real} +DenseComplex = {Real ~ Real} + + +Entry = { SparseComplex | SparseReal | SparsePattern | DenseComplex | DenseReal } + +// end of file, a silent way, see https://github.com/pest-parser/pest/issues/304#issuecomment-427198507 +eoi = _{ !ANY } + +Document = { + SOI ~ + NEWLINE* ~ + Header ~ + (NEWLINE ~ Comments)* ~ + (NEWLINE ~ Shape) ~ + (NEWLINE ~ Entry?)* ~ + eoi +} \ No newline at end of file diff --git a/nalgebra-sparse/src/io/matrix_market.rs b/nalgebra-sparse/src/io/matrix_market.rs new file mode 100644 index 00000000..dea284ee --- /dev/null +++ b/nalgebra-sparse/src/io/matrix_market.rs @@ -0,0 +1,1331 @@ +//! Implementation of matrix market io code. +//! +//! See the [website](https://math.nist.gov/MatrixMarket/formats.html) or the [paper](https://www.researchgate.net/publication/2630533_The_Matrix_Market_Exchange_Formats_Initial_Design) for more details about matrix market. +use crate::coo::CooMatrix; +use crate::SparseFormatError; +use crate::SparseFormatErrorKind; +use nalgebra::Complex; +use pest::iterators::Pairs; +use pest::Parser; +use std::cmp::PartialEq; +use std::convert::Infallible; +use std::convert::TryFrom; +use std::fmt; +use std::fmt::Formatter; +use std::fs; +use std::num::ParseIntError; +use std::num::TryFromIntError; +use std::path::Path; +use std::str::FromStr; + +/// A description of the error that occurred during importing a matrix from a matrix market format data. +#[derive(Debug)] +pub struct MatrixMarketError { + error_kind: MatrixMarketErrorKind, + message: String, +} + +/// Errors produced by functions that expect well-formed matrix market format data. +#[non_exhaustive] +#[derive(Copy, Clone, Debug, PartialEq)] +pub enum MatrixMarketErrorKind { + /// Parsing failure. + /// + /// Indicates that the parser failed, for example due to an unexpected string. + /// + /// Examples + /// -------- + /// ```rust + /// # use nalgebra_sparse::io::load_coo_from_matrix_market_str; + /// # use nalgebra_sparse::io::MatrixMarketErrorKind; + /// let str = r#" + /// %%MatrixMarket invalid invalid invalid invalid + /// 1 1 1 + /// 1 1 5 + /// "#; + /// let matrix_result = load_coo_from_matrix_market_str::(str); + /// assert_eq!(matrix_result.is_err(), true); + /// assert_eq!(matrix_result.unwrap_err().kind(), MatrixMarketErrorKind::ParsingError); + /// ``` + ParsingError, + + /// Indicates that the matrix market header is invalid. + /// + /// Examples + /// -------- + /// ```rust + /// # use nalgebra_sparse::io::load_coo_from_matrix_market_str; + /// # use nalgebra_sparse::io::MatrixMarketErrorKind; + /// let str = r#" + /// %%MatrixMarket matrix coordinate real hermitian + /// % a real matrix can't be hermitian + /// 1 1 1 + /// 1 1 5 + /// "#; + /// let matrix_result = load_coo_from_matrix_market_str::(str); + /// assert_eq!(matrix_result.is_err(), true); + /// assert_eq!(matrix_result.unwrap_err().kind(),MatrixMarketErrorKind::InvalidHeader); + /// ``` + InvalidHeader, + + /// Indicates that the number of data entries in the matrix market file does not match the header. + /// + /// Examples + /// -------- + /// ```rust + /// # use nalgebra_sparse::io::load_coo_from_matrix_market_str; + /// # use nalgebra_sparse::io::MatrixMarketErrorKind; + /// let str = r#" + /// %%matrixmarket matrix coordinate real general + /// % it has one more data entry than specified. + /// 3 3 1 + /// 2 2 2 + /// 2 3 2 + /// "#; + /// let matrix_result = load_coo_from_matrix_market_str::(str); + /// assert_eq!(matrix_result.is_err(), true); + /// assert_eq!(matrix_result.unwrap_err().kind(),MatrixMarketErrorKind::EntryMismatch); + /// ``` + EntryMismatch, + + /// Indicates that the scalar type requested is not compatible with the scalar type stored + /// in the matrix market file. + /// + /// Examples + /// -------- + /// ```rust + /// # use nalgebra_sparse::io::load_coo_from_matrix_market_str; + /// # use nalgebra_sparse::io::MatrixMarketErrorKind; + /// let str = r#" + /// %%matrixmarket matrix coordinate real general + /// % it should be loaded with load_coo_from_matrix_market_str::(str) (or f32) + /// 3 3 2 + /// 2 2 2.22 + /// 2 3 2.22 + /// "#; + /// let matrix_result = load_coo_from_matrix_market_str::(str); + /// assert_eq!(matrix_result.is_err(), true); + /// assert_eq!(matrix_result.unwrap_err().kind(),MatrixMarketErrorKind::TypeMismatch); + /// ``` + TypeMismatch, + + /// Indicates that zero has been used as an index in the data. + /// + /// **Note**: The matrix market format uses 1-based indexing. + /// + /// Examples + /// -------- + /// ```rust + /// # use nalgebra_sparse::io::load_coo_from_matrix_market_str; + /// # use nalgebra_sparse::io::MatrixMarketErrorKind; + /// let str = r#" + /// %%matrixmarket matrix coordinate real general + /// 1 1 1 + /// 0 0 10 + /// "#; + /// let matrix_result = load_coo_from_matrix_market_str::(str); + /// assert_eq!(matrix_result.is_err(), true); + /// assert_eq!(matrix_result.unwrap_err().kind(),MatrixMarketErrorKind::ZeroError); + /// ``` + ZeroError, + + /// Indicates [SparseFormatError] while creating the sparse matrix. + /// + /// + /// Examples + /// -------- + /// ```rust + /// # use nalgebra_sparse::io::load_coo_from_matrix_market_str; + /// # use nalgebra_sparse::io::MatrixMarketErrorKind; + /// # use nalgebra_sparse::SparseFormatErrorKind; + /// let str = r#" + /// %%matrixmarket matrix coordinate real general + /// 1 1 1 + /// 4 2 10 + /// "#; + /// let matrix_result = load_coo_from_matrix_market_str::(str); + /// assert_eq!(matrix_result.is_err(), true); + /// assert_eq!(matrix_result.unwrap_err().kind(), + /// MatrixMarketErrorKind::SparseFormatError(SparseFormatErrorKind::IndexOutOfBounds)); + /// ``` + SparseFormatError(SparseFormatErrorKind), + + /// Indicates that a wrong diagonal element has been provided to the matrix. + /// + /// Examples + /// -------- + /// ```rust + /// # use nalgebra_sparse::io::load_coo_from_matrix_market_str; + /// # use nalgebra_sparse::io::MatrixMarketErrorKind; + /// # use nalgebra::Complex; + /// let str = r#" + /// %%matrixmarket matrix coordinate real skew-symmetric + /// % skew-symmetric matrix can't have element on diagonal + /// 5 5 2 + /// 1 1 10 + /// 2 1 5 + /// "#; + /// let matrix_result = load_coo_from_matrix_market_str::(str); + /// assert_eq!(matrix_result.is_err(), true); + /// assert_eq!(matrix_result.unwrap_err().kind(),MatrixMarketErrorKind::DiagonalError); + /// + /// let str = r#" + /// %%matrixmarket matrix coordinate complex hermitian + /// % hermitian matrix diagonal element must be a real number + /// 5 5 2 + /// 1 1 10 2 + /// 2 1 5 2 + /// "#; + /// let matrix_result = load_coo_from_matrix_market_str::>(str); + /// assert_eq!(matrix_result.is_err(), true); + /// assert_eq!(matrix_result.unwrap_err().kind(),MatrixMarketErrorKind::DiagonalError); + /// ``` + /// Here the skew matrix shouldn't have an element on the diagonal. + DiagonalError, + + /// Indicates an [IO error](`std::io::Error`) while reading the data from file. + /// + /// Examples + /// -------- + /// ```rust + /// # use nalgebra_sparse::io::load_coo_from_matrix_market_file; + /// # use nalgebra_sparse::io::MatrixMarketErrorKind; + /// let matrix_result = load_coo_from_matrix_market_file::("matrix.mtx"); + /// assert_eq!(matrix_result.is_err(), true); + /// assert_eq!(matrix_result.unwrap_err().kind(),MatrixMarketErrorKind::IOError(std::io::ErrorKind::NotFound)); + /// ``` + IOError(std::io::ErrorKind), + + /// Indicates that a (skew-)symmetric (or hermitian) matrix is not a lower triangular matrix. + /// + /// Examples + /// -------- + /// ```rust + /// # use nalgebra_sparse::io::load_coo_from_matrix_market_str; + /// # use nalgebra_sparse::io::MatrixMarketErrorKind; + /// let str = r#" + /// %%matrixmarket matrix coordinate integer symmetric + /// 5 5 2 + /// 1 1 10 + /// 2 3 5 + /// "#; + /// let matrix_result = load_coo_from_matrix_market_str::(str); + /// assert_eq!(matrix_result.is_err(), true); + /// assert_eq!(matrix_result.unwrap_err().kind(),MatrixMarketErrorKind::NotLowerTriangle); + /// ``` + NotLowerTriangle, + + /// Indicates that a (skew-)symmetric (or hermitian) matrix is not a square matrix. + /// + /// Examples + /// -------- + /// ```rust + /// # use nalgebra_sparse::io::load_coo_from_matrix_market_str; + /// # use nalgebra_sparse::io::MatrixMarketErrorKind; + /// let str = r#" + /// %%matrixmarket matrix coordinate integer symmetric + /// 5 4 2 + /// 1 1 10 + /// 3 2 5 + /// "#; + /// let matrix_result = load_coo_from_matrix_market_str::(str); + /// assert_eq!(matrix_result.is_err(), true); + /// assert_eq!(matrix_result.unwrap_err().kind(),MatrixMarketErrorKind::NonSquare); + /// ``` + NonSquare, +} + +impl MatrixMarketError { + fn from_kind_and_message(error_type: MatrixMarketErrorKind, message: String) -> Self { + Self { + error_kind: error_type, + message, + } + } + + /// The matrix market error kind. + #[must_use] + pub fn kind(&self) -> MatrixMarketErrorKind { + self.error_kind + } + + /// The underlying error message. + #[must_use] + pub fn message(&self) -> &str { + self.message.as_str() + } +} + +impl fmt::Display for MatrixMarketError { + fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result { + write!(f, "Matrix Market error: ")?; + match self.kind() { + MatrixMarketErrorKind::ParsingError => { + write!(f, "ParsingError,")?; + } + MatrixMarketErrorKind::InvalidHeader => { + write!(f, "InvalidHeader,")?; + } + MatrixMarketErrorKind::EntryMismatch => { + write!(f, "EntryNumUnmatched,")?; + } + MatrixMarketErrorKind::TypeMismatch => { + write!(f, "TypeMismatch,")?; + } + MatrixMarketErrorKind::SparseFormatError(_) => { + write!(f, "SparseFormatError,")?; + } + MatrixMarketErrorKind::ZeroError => { + write!(f, "ZeroError,")?; + } + MatrixMarketErrorKind::IOError(_) => { + write!(f, "IOError,")?; + } + MatrixMarketErrorKind::DiagonalError => { + write!(f, "DiagonalError,")?; + } + MatrixMarketErrorKind::NotLowerTriangle => { + write!(f, "NotLowerTriangle,")?; + } + MatrixMarketErrorKind::NonSquare => { + write!(f, "NotSquareMatrix,")?; + } + } + write!(f, " message: {}", self.message) + } +} + +impl std::error::Error for MatrixMarketError {} + +impl MatrixMarketError { + fn from_pest_error(error: pest::error::Error) -> Self + where + T: fmt::Debug + std::hash::Hash + std::marker::Copy + Ord, + { + Self::from_kind_and_message( + MatrixMarketErrorKind::ParsingError, + format!("Can't parse the data.\n Error: {}", error), + ) + } +} + +impl From for MatrixMarketError { + fn from(err: ParseIntError) -> Self { + Self::from_kind_and_message( + MatrixMarketErrorKind::ParsingError, + format!("Can't parse data as i128.\n Error: {}", err), + ) + } +} + +impl From for MatrixMarketError { + fn from(err: SparseFormatError) -> Self { + Self::from_kind_and_message( + MatrixMarketErrorKind::SparseFormatError(*err.kind()), + format!("{}", &err), + ) + } +} + +impl From for MatrixMarketError { + fn from(err: std::io::Error) -> Self { + Self::from_kind_and_message( + MatrixMarketErrorKind::IOError(err.kind()), + format!("{}", &err), + ) + } +} + +impl From for MatrixMarketError { + fn from(err: TryFromIntError) -> Self { + Self::from_kind_and_message( + MatrixMarketErrorKind::TypeMismatch, + format!( + "Please consider using a larger integer type. Error message: {}", + &err + ), + ) + } +} + +// This is needed when calling `i128::try_from(i: i128)` +// but it won't happen +impl From for MatrixMarketError { + fn from(_err: Infallible) -> Self { + Self::from_kind_and_message( + MatrixMarketErrorKind::TypeMismatch, + format!("This won't happen"), + ) + } +} + +#[derive(Debug, PartialEq)] +enum Sparsity { + Sparse, + Dense, +} +#[derive(Debug, PartialEq)] +enum DataType { + Real, + Complex, + Pattern, + Integer, +} +#[derive(Debug, PartialEq)] +enum StorageScheme { + Symmetric, + General, + Skew, + Hermitian, +} +#[derive(Debug, PartialEq)] +struct Typecode { + sparsity: Sparsity, + datatype: DataType, + storagescheme: StorageScheme, +} + +impl FromStr for Sparsity { + type Err = MatrixMarketError; + /// Assumes that `word` is already lower case. + fn from_str(word: &str) -> Result { + match word { + "coordinate" => Ok(Sparsity::Sparse), + "array" => Ok(Sparsity::Dense), + _ => Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::ParsingError, + format!("keyword {} is unknown", word), + )), + } + } +} + +impl FromStr for DataType { + type Err = MatrixMarketError; + /// Assumes that `word` is already lower case. + fn from_str(word: &str) -> Result { + match word { + "real" => Ok(DataType::Real), + "complex" => Ok(DataType::Complex), + "integer" => Ok(DataType::Integer), + "pattern" => Ok(DataType::Pattern), + _ => Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::ParsingError, + format!("keyword {} is unknown", word), + )), + } + } +} + +impl FromStr for StorageScheme { + type Err = MatrixMarketError; + /// Assumes that `word` is already lower case. + fn from_str(word: &str) -> Result { + match word { + "skew-symmetric" => Ok(StorageScheme::Skew), + "general" => Ok(StorageScheme::General), + "symmetric" => Ok(StorageScheme::Symmetric), + "hermitian" => Ok(StorageScheme::Hermitian), + _ => Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::ParsingError, + format!("keyword {} is unknown", word), + )), + } + } +} + +/// Precheck if it's a valid header. +/// +/// For more details, please check +/// +/// Boisvert, Ronald F., Roldan Pozo, and Karin A. Remington. +/// The matrix market formats: Initial design. +/// Technical report, Applied and Computational Mathematics Division, NIST, 1996. Section 3. +fn typecode_precheck(tc: &Typecode) -> Result<(), MatrixMarketError> { + match tc { + Typecode { + datatype: DataType::Real, + storagescheme: StorageScheme::Hermitian, + .. + } => Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::InvalidHeader, + String::from("Real matrix can't be hermitian."), + )), + Typecode { + datatype: DataType::Integer, + storagescheme: StorageScheme::Hermitian, + .. + } => Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::InvalidHeader, + String::from("Integer matrix can't be hermitian."), + )), + Typecode { + datatype: DataType::Pattern, + storagescheme: StorageScheme::Hermitian, + .. + } => Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::InvalidHeader, + String::from("Pattern matrix can't be hermitian."), + )), + Typecode { + datatype: DataType::Pattern, + storagescheme: StorageScheme::Skew, + .. + } => Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::InvalidHeader, + String::from("Pattern matrix can't be skew-symmetric."), + )), + Typecode { + datatype: DataType::Pattern, + sparsity: Sparsity::Dense, + .. + } => Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::InvalidHeader, + String::from("Dense matrix can't be pattern matrix."), + )), + // precheck success + _ => Ok(()), + } +} + +/// Scalar types supported by the matrix market parser. +mod internal { + use crate::io::MatrixMarketError; + use na::{Complex, Scalar}; + + pub trait SupportedMatrixMarketScalar: Scalar { + /// When the matrix is an integer matrix, it will convert a [i128] number to this type. + fn from_i128(i: i128) -> Result; + /// When matrix is a Real matrix, it will convert a [f64] number to this type. + fn from_f64(f: f64) -> Result; + /// When matrix is a Complex matrix, it will convert a [Complex] number to this type. + fn from_c64(c: Complex) -> Result; + /// When matrix is a Pattern matrix, it will convert a unit type [unit] to this type. + fn from_pattern(p: ()) -> Result; + /// When matrix is a Skew-symmetric matrix, it will convert itself to its negative. + fn negative(self) -> Result; + /// When matrix is a Hermitian matrix, it will convert itself to its conjugate. + fn conjugate(self) -> Result; + } +} + +/// A marker trait for supported matrix market scalars. +/// +/// This is a sealed trait; it cannot be implemented by external crates. This is done in order to prevent leaking +/// some of the implementation details we currently rely on. We may relax this restriction in the future. +pub trait MatrixMarketScalar: internal::SupportedMatrixMarketScalar {} + +/// Implement MatrixMarketScalar for primitive integer types. +macro_rules! mm_int_impl { + ($T:ty) => { + impl MatrixMarketScalar for $T {} + + impl internal::SupportedMatrixMarketScalar for $T { + #[inline] + fn from_i128(i: i128) -> Result { + Ok(Self::try_from(i)?) + } + #[inline] + fn from_f64(_f: f64) -> Result { + Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::TypeMismatch, + format!("Int type can't be parsed from f64"), + )) + } + #[inline] + fn from_c64(_c: Complex) -> Result { + Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::TypeMismatch, + format!("Int type can't be parsed from Complex"), + )) + } + #[inline] + fn from_pattern(_p: ()) -> Result { + Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::TypeMismatch, + format!("Int type can't be parsed from ()"), + )) + } + #[inline] + fn conjugate(self) -> Result { + Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::TypeMismatch, + format!("Int type has no conjugate"), + )) + } + #[inline] + fn negative(self) -> Result { + Ok(-self) + } + } + }; +} +/// Implement MatrixMarketScalar for primitive real types. +macro_rules! mm_real_impl { + ($T:ty) => { + impl MatrixMarketScalar for $T {} + + impl internal::SupportedMatrixMarketScalar for $T { + #[inline] + fn from_i128(_i: i128) -> Result { + Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::TypeMismatch, + format!("real type can't be parsed from i128"), + )) + } + #[inline] + fn from_f64(f: f64) -> Result { + Ok(f as Self) + } + #[inline] + fn from_c64(_c: Complex) -> Result { + Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::TypeMismatch, + format!("real type can't be parsed from Complex"), + )) + } + #[inline] + fn from_pattern(_p: ()) -> Result { + Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::TypeMismatch, + format!("real type can't be parsed from ()"), + )) + } + #[inline] + fn conjugate(self) -> Result { + Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::TypeMismatch, + format!("real type has no conjugate"), + )) + } + #[inline] + fn negative(self) -> Result { + Ok(-self) + } + } + }; +} + +/// Implement MatrixMarketScalar for primitive complex types. +macro_rules! mm_complex_impl { + ($T:ty) => { + impl MatrixMarketScalar for Complex<$T> {} + + impl internal::SupportedMatrixMarketScalar for Complex<$T> { + #[inline] + fn from_i128(_i: i128) -> Result { + Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::TypeMismatch, + format!("Complex type can't be parsed from i128"), + )) + } + #[inline] + fn from_f64(_f: f64) -> Result { + Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::TypeMismatch, + format!("Complex type can't be parsed from f64"), + )) + } + #[inline] + fn from_c64(c: Complex) -> Result { + Ok(Self { + re: c.re as $T, + im: c.im as $T, + }) + } + #[inline] + fn from_pattern(_p: ()) -> Result { + Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::TypeMismatch, + format!("Complex type can't be parsed from ()"), + )) + } + #[inline] + fn conjugate(self) -> Result { + Ok(self.conj()) + } + #[inline] + fn negative(self) -> Result { + Ok(-self) + } + } + }; +} +/// Implement MatrixMarketScalar for primitive unit types. +macro_rules! mm_pattern_impl { + ($T:ty) => { + impl MatrixMarketScalar for $T {} + + impl internal::SupportedMatrixMarketScalar for $T { + #[inline] + fn from_i128(_i: i128) -> Result { + Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::TypeMismatch, + format!("Pattern type can't be parsed from i128"), + )) + } + #[inline] + fn from_f64(_f: f64) -> Result { + Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::TypeMismatch, + format!("Pattern type can't be parsed from f64"), + )) + } + #[inline] + fn from_c64(_c: Complex) -> Result { + Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::TypeMismatch, + format!("Pattern type can't be parsed from Complex"), + )) + } + #[inline] + fn from_pattern(p: ()) -> Result { + Ok(p) + } + + #[inline] + fn conjugate(self) -> Result { + Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::TypeMismatch, + format!("Pattern type has no conjugate"), + )) + } + #[inline] + fn negative(self) -> Result { + Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::TypeMismatch, + format!("Pattern type has no negative"), + )) + } + } + }; +} + +mm_int_impl!(i8); +mm_int_impl!(i16); +mm_int_impl!(i32); +mm_int_impl!(i64); +mm_int_impl!(i128); + +mm_real_impl!(f32); +mm_real_impl!(f64); + +mm_complex_impl!(f32); +mm_complex_impl!(f64); + +mm_pattern_impl!(()); + +#[derive(Parser)] +#[grammar = "io/matrix_market.pest"] +struct MatrixMarketParser; + +/// Parses a Matrix Market file at the given path as a `CooMatrix`. +/// +/// The matrix market format specification does not clarify whether duplicate entries are allowed. Our importer +/// assumes that this is permitted and produces a `CooMatrix` with possibly duplicate entries. +/// +/// **Note**: A current restriction of the importer is that you must use a compatible scalar type when importing. +/// For example, in order to import a matrix stored as `integer` in the matrix market format, you must +/// import it as an integer matrix, otherwise a [TypeMismatch](MatrixMarketErrorKind::TypeMismatch) error +/// will be returned. This restriction may be lifted in the future, and is +/// tracked by issue [#1038](https://github.com/dimforge/nalgebra/issues/1038). +/// +/// Errors +/// -------- +/// +/// See [MatrixMarketErrorKind] for a list of possible error conditions. +/// +/// Examples +/// -------- +/// ```no_run +/// use nalgebra_sparse::io::load_coo_from_matrix_market_file; +/// // Use e.g. `i32` for integer matrices +/// let matrix = load_coo_from_matrix_market_file::("path/to/matrix.mtx").unwrap(); +/// ``` +pub fn load_coo_from_matrix_market_file>( + path: P, +) -> Result, MatrixMarketError> +where + T: MatrixMarketScalar, +{ + let file = fs::read_to_string(path)?; + load_coo_from_matrix_market_str(&file) +} + +/// Parses a Matrix Market file described by the given string as a `CooMatrix`. +/// +/// See [load_coo_from_matrix_market_file] for more information. +/// +/// Errors +/// -------- +/// +/// See [MatrixMarketErrorKind] for a list of possible error conditions. +/// +/// Examples +/// -------- +/// ``` +/// use nalgebra_sparse::io::load_coo_from_matrix_market_str; +/// let str = r#" +/// %%matrixmarket matrix coordinate integer general +/// 5 4 2 +/// 1 1 10 +/// 2 3 5 +/// "#; +/// // Use e.g. `i32` for integer matrices +/// let matrix = load_coo_from_matrix_market_str::(str).unwrap(); +/// ``` +pub fn load_coo_from_matrix_market_str(data: &str) -> Result, MatrixMarketError> +where + T: MatrixMarketScalar, +{ + // unwrap() in this function are guaranteed by parsing the data + let file = MatrixMarketParser::parse(Rule::Document, data) + .map_err(MatrixMarketError::from_pest_error)? + .next() + .unwrap(); + + let mut rows: Vec = Vec::new(); + let mut cols: Vec = Vec::new(); + let mut data: Vec = Vec::new(); + let mut lines = file.into_inner(); + + let header_line = lines.next().unwrap(); + let header_type = parse_header(&mut header_line.into_inner()); + typecode_precheck(&header_type)?; + + let shape_line = lines.next().unwrap(); + // shape here is number of rows, columns, non-zeros + let shape: (usize, usize, usize); + match header_type.sparsity { + Sparsity::Sparse => { + shape = parse_sparse_shape(&mut shape_line.into_inner(), &header_type.storagescheme)?; + } + Sparsity::Dense => { + shape = parse_dense_shape(&mut shape_line.into_inner(), &header_type.storagescheme)?; + } + } + + // used when constructing dense matrix. + // If it's sparse matrix, it has no effect. + let mut current_dense_coordinate: (usize, usize) = (0, 0); + if header_type.storagescheme == StorageScheme::Skew { + // for skew dense matrix, the first element starts from (1,0) + current_dense_coordinate = (1, 0); + } + // count how many entries in the matrix data + let count = lines.clone().count(); + if count != shape.2 { + return Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::EntryMismatch, + format!( + "{} entries required for the matrix, but {} was provided", + shape.2, count, + ), + )); + } + + for data_line in lines { + let entry: (usize, usize, T); + match header_type { + Typecode { + sparsity: Sparsity::Sparse, + datatype: DataType::Real, + .. + } => { + entry = parse_sparse_real::(&mut data_line.into_inner())?; + } + Typecode { + sparsity: Sparsity::Sparse, + datatype: DataType::Integer, + .. + } => { + entry = parse_sparse_int::(&mut data_line.into_inner())?; + } + Typecode { + sparsity: Sparsity::Sparse, + datatype: DataType::Pattern, + .. + } => { + entry = parse_sparse_pattern::(&mut data_line.into_inner())?; + } + Typecode { + sparsity: Sparsity::Sparse, + datatype: DataType::Complex, + .. + } => { + entry = parse_sparse_complex::(&mut data_line.into_inner())?; + } + Typecode { + sparsity: Sparsity::Dense, + datatype: DataType::Complex, + .. + } => { + entry = ( + current_dense_coordinate.0, + current_dense_coordinate.1, + parse_dense_complex::(&mut data_line.into_inner())?, + ); + next_dense_coordinate( + &mut current_dense_coordinate, + shape, + &header_type.storagescheme, + ); + } + Typecode { + sparsity: Sparsity::Dense, + datatype: DataType::Real, + .. + } => { + entry = ( + current_dense_coordinate.0, + current_dense_coordinate.1, + parse_dense_real::(&mut data_line.into_inner())?, + ); + + next_dense_coordinate( + &mut current_dense_coordinate, + shape, + &header_type.storagescheme, + ); + } + Typecode { + sparsity: Sparsity::Dense, + datatype: DataType::Integer, + .. + } => { + entry = ( + current_dense_coordinate.0, + current_dense_coordinate.1, + parse_dense_int::(&mut data_line.into_inner())?, + ); + next_dense_coordinate( + &mut current_dense_coordinate, + shape, + &header_type.storagescheme, + ); + } + _ => { + // it shouldn't happen here, because dense matrix can't be pattern. And it will give InvalidHeader error beforehand. + entry = (1, 1, T::from_i128(1)?) + } + } + + let (r, c, d) = entry; + + match header_type.storagescheme { + StorageScheme::General => { + rows.push(r); + cols.push(c); + data.push(d); + } + StorageScheme::Symmetric => { + check_lower_triangle(r, c)?; + rows.push(r); + cols.push(c); + data.push(d.clone()); + // don't need to add twice if the element in on diagonal + if r != c { + rows.push(c); + cols.push(r); + data.push(d); + } + } + StorageScheme::Skew => { + check_lower_triangle(r, c)?; + rows.push(r); + cols.push(c); + data.push(d.clone()); + // skew-symmetric matrix shouldn't have diagonal element + if r == c { + return Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::DiagonalError, + format!( + "There is a diagonal element in skew matrix, in row(and column) {}", + r + 1 + ), + )); + } + rows.push(c); + cols.push(r); + data.push(d.negative()?); + } + StorageScheme::Hermitian => { + check_lower_triangle(r, c)?; + rows.push(r); + cols.push(c); + data.push(d.clone()); + + if r == c && d != d.clone().conjugate()? { + return Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::DiagonalError, + format!( + "There is a diagonal element in hermitian matrix, which is not a real number, in row(and column) {}", + r + 1 + ), + )); + } + // don't need to add twice if the element in on diagonal + if r != c { + rows.push(c); + cols.push(r); + data.push(d.conjugate()?); + } + } + } + } + Ok(CooMatrix::try_from_triplets( + shape.0, shape.1, rows, cols, data, + )?) +} + +#[inline] +/// do a quick check it the entry is in the lower triangle part of the matrix +fn check_lower_triangle(r: usize, c: usize) -> Result<(), MatrixMarketError> { + if c > r { + return Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::NotLowerTriangle, + format!( + "Entry: row {} col {} should be put into lower triangle", + r, c + ), + )); + } + Ok(()) +} + +#[inline] +/// Parse a pest structure to a Typecode of the matrix. +fn parse_header(inner: &mut Pairs<'_, Rule>) -> Typecode { + // unwrap() in this function are guaranteed by parsing the data + Typecode { + sparsity: inner + .next() + .unwrap() + .as_str() + .to_ascii_lowercase() + .parse::() + .unwrap(), + datatype: inner + .next() + .unwrap() + .as_str() + .to_ascii_lowercase() + .parse::() + .unwrap(), + storagescheme: inner + .next() + .unwrap() + .as_str() + .to_ascii_lowercase() + .parse::() + .unwrap(), + } +} + +// Parse shape starts here------------------------------------------------- + +/// Parse a pest structure to sparse shape information, including 3 int, which are number of rows, cols and non-zeros. +fn parse_sparse_shape( + inner: &mut Pairs<'_, Rule>, + storagescheme: &StorageScheme, +) -> Result<(usize, usize, usize), MatrixMarketError> { + // unwrap() in this function are guaranteed by parsing the data + let shape_inner = inner.next().unwrap(); + if shape_inner.as_rule() != Rule::SparseShape { + return Err(MatrixMarketError::from_kind_and_message(MatrixMarketErrorKind::ParsingError,format!(" + Shape shape line requires 3 int numbers as number of rows, columns and non-zeros, but line {} was provided here. + ",shape_inner.as_str()))); + } + + let mut inner = shape_inner.into_inner(); + + let r = inner.next().unwrap().as_str().parse::().unwrap(); + let c = inner.next().unwrap().as_str().parse::().unwrap(); + let nnz = inner.next().unwrap().as_str().parse::().unwrap(); + + // check for square matrix, when it's not a general matrix + if *storagescheme != StorageScheme::General && r != c { + return Err(MatrixMarketError::from_kind_and_message(MatrixMarketErrorKind::NonSquare, format!("(Skew-)Symmetric or hermitian matrix should be square matrix, but it has dimension {} and {}", r, c))); + } + + Ok((r, c, nnz)) +} + +/// Parse a pest structure to dense shape information, including 2 int, which are number of rows, cols. +fn parse_dense_shape( + inner: &mut Pairs<'_, Rule>, + storagescheme: &StorageScheme, +) -> Result<(usize, usize, usize), MatrixMarketError> { + // unwrap() in this function are guaranteed by parsing the data + let shape_inner = inner.next().unwrap(); + if shape_inner.as_rule() != Rule::DenseShape { + return Err(MatrixMarketError::from_kind_and_message(MatrixMarketErrorKind::ParsingError,format!(" + Shape shape line requires 2 int numbers as number of rows, columns, but line {} was provided here. + ",shape_inner.as_str()))); + } + + let mut inner = shape_inner.into_inner(); + let r = inner.next().unwrap().as_str().parse::().unwrap(); + let c = inner.next().unwrap().as_str().parse::().unwrap(); + + // check for square matrix, when it's not a general matrix + if *storagescheme != StorageScheme::General && r != c { + return Err(MatrixMarketError::from_kind_and_message(MatrixMarketErrorKind::NonSquare, format!("(Skew-)Symmetric or hermitian matrix should be square matrix, but it has dimension {} and {}", r, c))); + } + + let n: usize; + // Calculate the number of entries in the dense matrix + match storagescheme { + StorageScheme::General => { + // general matrix should contain r*c entries + n = r * c; + } + StorageScheme::Symmetric | StorageScheme::Hermitian => { + // it must be square matrix, so r==c is true here + // Symmetric or Hermitian should contain 1+2...+r = r*(r+1)/2 entries + n = r * (r + 1) / 2; + } + StorageScheme::Skew => { + // it must be square matrix, so r==c is true here + // Skew-Symmetric should contain 1+2...+r-1 = r*(r-1)/2 entries + n = r * (r - 1) / 2; + } + } + + Ok((r, c, n)) +} + +// Parse shape ends here------------------------------------------------- + +// Parse entry starts here------------------------------------------------- + +/// Parse a pest structure to sparse real entry, including 2 int, which are number of rows, cols, and a real number as data +fn parse_sparse_real(inner: &mut Pairs<'_, Rule>) -> Result<(usize, usize, T), MatrixMarketError> +where + T: MatrixMarketScalar, +{ + // unwrap() in this function are guaranteed by parsing the data + let entry_inner = inner.next().unwrap(); + if entry_inner.as_rule() != Rule::SparseReal { + return Err(MatrixMarketError::from_kind_and_message(MatrixMarketErrorKind::ParsingError,format!(" + Spare real matrix requires 2 int number as coordinates and 1 real number as data, but line {} was provided. + ",entry_inner.as_str() ))); + } + + let mut inner = entry_inner.into_inner(); + let (r, c) = parse_sparse_coordinate(&mut inner)?; + let d = inner.next().unwrap().as_str().parse::().unwrap(); + Ok((r, c, T::from_f64(d)?)) +} + +/// Parse a pest structure to sparse integer entry, including 2 int, which are number of rows, cols, and a int number as data +fn parse_sparse_int(inner: &mut Pairs<'_, Rule>) -> Result<(usize, usize, T), MatrixMarketError> +where + T: MatrixMarketScalar, +{ + // unwrap() in this function are guaranteed by parsing the data + let entry_inner = inner.next().unwrap(); + // Because integer numbers can also be parsed as float numbers, it will be checked again in `parse::()?` + if entry_inner.as_rule() != Rule::SparseReal { + return Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::ParsingError, + format!( + " + Spare real matrix requires 3 int number as coordinates and data, but line {} was provided. + ", + entry_inner.as_str() + ), + )); + } + let mut inner = entry_inner.into_inner(); + let (r, c) = parse_sparse_coordinate(&mut inner)?; + // Here to guarantee it is an integer number + let d = inner.next().unwrap().as_str().parse::()?; + Ok((r, c, T::from_i128(d)?)) +} + +/// Parse a pest structure to sparse pattern entry, including 2 int, which are number of rows, cols +fn parse_sparse_pattern( + inner: &mut Pairs<'_, Rule>, +) -> Result<(usize, usize, T), MatrixMarketError> +where + T: MatrixMarketScalar, +{ + // unwrap() in this function are guaranteed by parsing the data + let entry_inner = inner.next().unwrap(); + if entry_inner.as_rule() != Rule::SparsePattern { + return Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::ParsingError, + format!( + " + Spare real matrix requires 2 int number as coordinates, but line {} was provided. + ", + entry_inner.as_str() + ), + )); + } + let mut inner = entry_inner.into_inner(); + let (r, c) = parse_sparse_coordinate(&mut inner)?; + Ok((r, c, T::from_pattern(())?)) +} + +/// Parse a pest structure to sparse complex entry, including 2 int, which are number of rows, cols, and 2 real number as complex data +fn parse_sparse_complex( + inner: &mut Pairs<'_, Rule>, +) -> Result<(usize, usize, T), MatrixMarketError> +where + T: MatrixMarketScalar, +{ + // unwrap() in this function are guaranteed by parsing the data + let entry_inner = inner.next().unwrap(); + if entry_inner.as_rule() != Rule::SparseComplex { + return Err(MatrixMarketError::from_kind_and_message(MatrixMarketErrorKind::ParsingError,format!(" + Spare real matrix requires 2 int number as coordinates and 2 real number as complex data, but line {} was provided. + ",entry_inner.as_str() ))); + } + let mut inner = entry_inner.into_inner(); + let (r, c) = parse_sparse_coordinate(&mut inner)?; + let real = inner.next().unwrap().as_str().parse::().unwrap(); + let imag = inner.next().unwrap().as_str().parse::().unwrap(); + let complex = Complex::::new(real, imag); + Ok((r, c, T::from_c64(complex)?)) +} + +/// Parse a pest structure to dense real entry, including a real number as data +fn parse_dense_real(inner: &mut Pairs<'_, Rule>) -> Result +where + T: MatrixMarketScalar, +{ + // unwrap() in this function are guaranteed by parsing the data + let entry_inner = inner.next().unwrap(); + if entry_inner.as_rule() != Rule::DenseReal { + return Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::ParsingError, + format!( + " + Dense real matrix requires 1 real number as data, but line {} was provided. + ", + entry_inner.as_str() + ), + )); + } + let mut inner = entry_inner.into_inner(); + let d = inner.next().unwrap().as_str().parse::().unwrap(); + Ok(T::from_f64(d)?) +} + +/// Parse a pest structure to dense integer entry, including a integer number as data +fn parse_dense_int(inner: &mut Pairs<'_, Rule>) -> Result +where + T: MatrixMarketScalar, +{ + // unwrap() in this function are guaranteed by parsing the data + let entry_inner = inner.next().unwrap(); + // Because integer numbers can also be parsed as float numbers, it will be checked again in `parse::()?` + if entry_inner.as_rule() != Rule::DenseReal { + return Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::ParsingError, + format!( + " + Dense real matrix requires 1 int number as data, but line {} was provided. + ", + entry_inner.as_str() + ), + )); + } + let mut inner = entry_inner.into_inner(); + // Here to guarantee it is an integer number + let d = inner.next().unwrap().as_str().parse::()?; + Ok(T::from_i128(d)?) +} + +/// Parse a pest structure to dense complex entry, including 2 real number as complex data +fn parse_dense_complex(inner: &mut Pairs<'_, Rule>) -> Result +where + T: MatrixMarketScalar, +{ + // unwrap() in this function are guaranteed by parsing the data + let entry_inner = inner.next().unwrap(); + // Note: theoretically, 2 positive integers could also become the complex number, + // but it would be parsed as SparsePattern, because SparsePattern has higher priority. + // But DenseComplex can't have higher priority, + // because, it's more often to deal with "normal" SparsePattern, rather than "unnormal" DenseComplex + if entry_inner.as_rule() != Rule::DenseComplex && entry_inner.as_rule() != Rule::SparsePattern { + return Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::ParsingError, + format!( + " + Dense real matrix requires 2 real number as complex data, but line {} was provided. + ", + entry_inner.as_str() + ), + )); + } + let mut inner = entry_inner.into_inner(); + let real = inner.next().unwrap().as_str().parse::().unwrap(); + let imag = inner.next().unwrap().as_str().parse::().unwrap(); + let complex = Complex::::new(real, imag); + Ok(T::from_c64(complex)?) +} + +// Parse entry ends here------------------------------------------------- + +/// Parse the coordinates information used for sparse matrix +fn parse_sparse_coordinate( + inner: &mut Pairs<'_, Rule>, +) -> Result<(usize, usize), MatrixMarketError> { + // unwrap() in this function are guaranteed by parsing the data + let r = inner.next().unwrap().as_str().parse::().unwrap(); + let c = inner.next().unwrap().as_str().parse::().unwrap(); + if r * c == 0 { + return Err(MatrixMarketError::from_kind_and_message( + MatrixMarketErrorKind::ZeroError, + String::from("The data has to be one-indexed"), + )); + } + // The coordinates in matrix market is one-based, but in CooMatrix is zero-based. + Ok((r - 1, c - 1)) +} + +/// Calculate the next coordinates used for dense matrix +fn next_dense_coordinate( + current_dense_coordinate: &mut (usize, usize), + shape: (usize, usize, usize), + storagescheme: &StorageScheme, +) { + // matrix market is column based format. + // so it follows the order (0,0) -> (1,0) -> ... -> (row, 0) -> (0,1) -> ... ->(row,col) + // current_dense_coordinate is (row, column) + match storagescheme { + StorageScheme::General => { + if current_dense_coordinate.0 < shape.0 - 1 { + current_dense_coordinate.0 += 1 + } else { + // jump to next column, reset row to 1, column add 1 + current_dense_coordinate.0 = 0; + current_dense_coordinate.1 += 1; + } + } + StorageScheme::Symmetric | StorageScheme::Hermitian => { + if current_dense_coordinate.0 < shape.0 - 1 { + current_dense_coordinate.0 += 1 + } else { + // jump to next column, column add 1, then set row equals to current column + // for example (0,0) -> (1,0) -> ... -> (row, 0) -> (1,1) -> ... + current_dense_coordinate.1 += 1; + current_dense_coordinate.0 = current_dense_coordinate.1; + } + } + StorageScheme::Skew => { + if current_dense_coordinate.0 < shape.0 - 1 { + current_dense_coordinate.0 += 1; + } else { + // jump to next column, set row equals to current column, then column add 1 + // skew matrix doesn't have element on diagonal + // for example (1,0) -> (2,0) -> ... -> (row, 0) -> (2,1) -> ... + current_dense_coordinate.1 += 1; + current_dense_coordinate.0 = current_dense_coordinate.1 + 1; + } + } + } +} diff --git a/nalgebra-sparse/src/io/mod.rs b/nalgebra-sparse/src/io/mod.rs new file mode 100644 index 00000000..89b21ffb --- /dev/null +++ b/nalgebra-sparse/src/io/mod.rs @@ -0,0 +1,38 @@ +//! Functionality for importing and exporting sparse matrices to and from files. +//! +//! **Available only when the `io` feature is enabled.** +//! +//! The following formats are currently supported: +//! +//! | Format | Import | Export | +//! | ------------------------------------------------|------------|------------| +//! | [Matrix market](#matrix-market-format) | Yes | No | +//! +//! [Matrix market]: https://math.nist.gov/MatrixMarket/formats.html +//! +//! ## Matrix Market format +//! +//! The Matrix Market format is a simple ASCII-based file format for sparse matrices, and was initially developed for +//! the [NIST Matrix Market](https://math.nist.gov/MatrixMarket/), a repository of example sparse matrices. +//! In later years it has largely been superseded by the +//! [SuiteSparse Matrix Collection](https://sparse.tamu.edu/) (formerly University of Florida Sparse Matrix Collection), +//! which also uses the Matrix Market file format. +//! +//! We currently offer functionality for importing a Matrix market file to an instance of a +//! [CooMatrix](crate::CooMatrix) through the function [load_coo_from_matrix_market_file]. It is also possible to load +//! a matrix stored in the matrix market format with the function [load_coo_from_matrix_market_str]. +//! +//! Export is currently not implemented, but [planned](https://github.com/dimforge/nalgebra/issues/1037). +//! +//! Our implementation is based on the [format description](https://math.nist.gov/MatrixMarket/formats.html) +//! on the Matrix Market website and the +//! [following NIST whitepaper](https://math.nist.gov/MatrixMarket/reports/MMformat.ps): +//! +//! > Boisvert, Ronald F., Roldan Pozo, and Karin A. Remington.
+//! > "*The Matrix Market Exchange Formats: Initial Design.*" (1996). + +pub use self::matrix_market::{ + load_coo_from_matrix_market_file, load_coo_from_matrix_market_str, MatrixMarketError, + MatrixMarketErrorKind, MatrixMarketScalar, +}; +mod matrix_market; diff --git a/nalgebra-sparse/src/lib.rs b/nalgebra-sparse/src/lib.rs index bf845757..edbf83bd 100644 --- a/nalgebra-sparse/src/lib.rs +++ b/nalgebra-sparse/src/lib.rs @@ -19,6 +19,7 @@ //! - Sparsity patterns in CSR and CSC matrices are explicitly represented by the //! [SparsityPattern](pattern::SparsityPattern) type, which encodes the invariants of the //! associated index data structures. +//! - [Matrix market format support](`io`) when the `io` feature is enabled. //! - [proptest strategies](`proptest`) for sparse matrices when the feature //! `proptest-support` is enabled. //! - [matrixcompare support](https://crates.io/crates/matrixcompare) for effortless @@ -142,11 +143,19 @@ )] pub extern crate nalgebra as na; +#[cfg(feature = "io")] +extern crate pest; +#[macro_use] +#[cfg(feature = "io")] +extern crate pest_derive; + pub mod convert; pub mod coo; pub mod csc; pub mod csr; pub mod factorization; +#[cfg(feature = "io")] +pub mod io; pub mod ops; pub mod pattern; diff --git a/nalgebra-sparse/src/ops/mod.rs b/nalgebra-sparse/src/ops/mod.rs index d7e6d432..9a73148c 100644 --- a/nalgebra-sparse/src/ops/mod.rs +++ b/nalgebra-sparse/src/ops/mod.rs @@ -67,7 +67,7 @@ //! As can be seen from the table, only `CSR * Dense` and `CSC * Dense` are supported. //! The other way around, i.e. `Dense * CSR` and `Dense * CSC` are not implemented. //! -//! Additionally, [CsrMatrix](`crate::csr::CsrMatrix`) and [CooMatrix](`crate::coo::CooMatrix`) +//! Additionally, [CsrMatrix](`crate::csr::CsrMatrix`) and [CscMatrix](`crate::csc::CscMatrix`) //! support multiplication with scalars, in addition to division by a scalar. //! Note that only `Matrix * Scalar` works in a generic context, although `Scalar * Matrix` //! has been implemented for many of the built-in arithmetic types. This is due to a fundamental diff --git a/nalgebra-sparse/src/proptest.rs b/nalgebra-sparse/src/proptest.rs index 472c466f..93db7796 100644 --- a/nalgebra-sparse/src/proptest.rs +++ b/nalgebra-sparse/src/proptest.rs @@ -5,11 +5,6 @@ //! The strategies provided here are generally expected to be able to generate the entire range //! of possible outputs given the constraints on dimensions and values. However, there are no //! particular guarantees on the distribution of possible values. - -// Contains some patched code from proptest that we can remove in the (hopefully near) future. -// See docs in file for more details. -mod proptest_patched; - use crate::coo::CooMatrix; use crate::csc::CscMatrix; use crate::csr::CsrMatrix; @@ -31,16 +26,10 @@ fn dense_row_major_coord_strategy( let mut booleans = vec![true; nnz]; booleans.append(&mut vec![false; (nrows * ncols) - nnz]); // Make sure that exactly `nnz` of the booleans are true - - // TODO: We cannot use the below code because of a bug in proptest, see - // https://github.com/AltSysrq/proptest/pull/217 - // so for now we're using a patched version of the Shuffle adapter - // (see also docs in `proptest_patched` - // Just(booleans) - // // Need to shuffle to make sure they are randomly distributed - // .prop_shuffle() - - proptest_patched::Shuffle(Just(booleans)).prop_map(move |booleans| { + Just(booleans) + // Need to shuffle to make sure they are randomly distributed + .prop_shuffle() + .prop_map(move |booleans| { booleans .into_iter() .enumerate() diff --git a/nalgebra-sparse/src/proptest/proptest_patched.rs b/nalgebra-sparse/src/proptest/proptest_patched.rs deleted file mode 100644 index 37c71262..00000000 --- a/nalgebra-sparse/src/proptest/proptest_patched.rs +++ /dev/null @@ -1,146 +0,0 @@ -//! Contains a modified implementation of `proptest::strategy::Shuffle`. -//! -//! The current implementation in `proptest` does not generate all permutations, which is -//! problematic for our proptest generators. The issue has been fixed in -//! https://github.com/AltSysrq/proptest/pull/217 -//! but it has yet to be merged and released. As soon as this fix makes it into a new release, -//! the modified code here can be removed. -//! -/*! - This code has been copied and adapted from - https://github.com/AltSysrq/proptest/blob/master/proptest/src/strategy/shuffle.rs - The original licensing text is: - - //- - // Copyright 2017 Jason Lingle - // - // Licensed under the Apache License, Version 2.0 or the MIT license - // , at your - // option. This file may not be copied, modified, or distributed - // except according to those terms. - -*/ - -use proptest::num; -use proptest::prelude::Rng; -use proptest::strategy::{NewTree, Shuffleable, Strategy, ValueTree}; -use proptest::test_runner::{TestRng, TestRunner}; -use std::cell::Cell; - -#[derive(Clone, Debug)] -#[must_use = "strategies do nothing unless used"] -pub struct Shuffle(pub(super) S); - -impl Strategy for Shuffle -where - S::Value: Shuffleable, -{ - type Tree = ShuffleValueTree; - type Value = S::Value; - - fn new_tree(&self, runner: &mut TestRunner) -> NewTree { - let rng = runner.new_rng(); - - self.0.new_tree(runner).map(|inner| ShuffleValueTree { - inner, - rng, - dist: Cell::new(None), - simplifying_inner: false, - }) - } -} - -#[derive(Clone, Debug)] -pub struct ShuffleValueTree { - inner: V, - rng: TestRng, - dist: Cell>, - simplifying_inner: bool, -} - -impl ShuffleValueTree -where - V::Value: Shuffleable, -{ - fn init_dist(&self, dflt: usize) -> usize { - if self.dist.get().is_none() { - self.dist.set(Some(num::usize::BinarySearch::new(dflt))); - } - - self.dist.get().unwrap().current() - } - - fn force_init_dist(&self) { - if self.dist.get().is_none() { - let _ = self.init_dist(self.current().shuffle_len()); - } - } -} - -impl ValueTree for ShuffleValueTree -where - V::Value: Shuffleable, -{ - type Value = V::Value; - - fn current(&self) -> V::Value { - let mut value = self.inner.current(); - let len = value.shuffle_len(); - // The maximum distance to swap elements. This could be larger than - // `value` if `value` has reduced size during shrinking; that's OK, - // since we only use this to filter swaps. - let max_swap = self.init_dist(len); - - // If empty collection or all swaps will be filtered out, there's - // nothing to shuffle. - if 0 == len || 0 == max_swap { - return value; - } - - let mut rng = self.rng.clone(); - - for start_index in 0..len - 1 { - // Determine the other index to be swapped, then skip the swap if - // it is too far. This ordering is critical, as it ensures that we - // generate the same sequence of random numbers every time. - - // NOTE: The below line is the whole reason for the existence of this adapted code - // We need to be able to swap with the same element, so that some elements remain in - // place rather being swapped - // let end_index = rng.gen_range(start_index + 1..len); - let end_index = rng.gen_range(start_index..len); - if end_index - start_index <= max_swap { - value.shuffle_swap(start_index, end_index); - } - } - - value - } - - fn simplify(&mut self) -> bool { - if self.simplifying_inner { - self.inner.simplify() - } else { - // Ensure that we've initialised `dist` to *something* to give - // consistent non-panicking behaviour even if called in an - // unexpected sequence. - self.force_init_dist(); - if self.dist.get_mut().as_mut().unwrap().simplify() { - true - } else { - self.simplifying_inner = true; - self.inner.simplify() - } - } - } - - fn complicate(&mut self) -> bool { - if self.simplifying_inner { - self.inner.complicate() - } else { - self.force_init_dist(); - self.dist.get_mut().as_mut().unwrap().complicate() - } - } -} diff --git a/nalgebra-sparse/tests/unit.rs b/nalgebra-sparse/tests/unit.rs index 73a95cd7..74d32a40 100644 --- a/nalgebra-sparse/tests/unit.rs +++ b/nalgebra-sparse/tests/unit.rs @@ -1,6 +1,9 @@ //! Unit tests -#[cfg(any(not(feature = "proptest-support"), not(feature = "compare")))] -compile_error!("Tests must be run with features `proptest-support` and `compare`"); +#[cfg(not(all(feature = "proptest-support", feature = "compare", feature = "io",)))] +compile_error!( + "Please enable the `proptest-support`, `compare` and `io` features in order to compile and run the tests. + Example: `cargo test -p nalgebra-sparse --features proptest-support,compare,io`" +); mod unit_tests; diff --git a/nalgebra-sparse/tests/unit_tests/csr.rs b/nalgebra-sparse/tests/unit_tests/csr.rs index dee1ae1e..3ca2f0dc 100644 --- a/nalgebra-sparse/tests/unit_tests/csr.rs +++ b/nalgebra-sparse/tests/unit_tests/csr.rs @@ -5,6 +5,8 @@ use nalgebra_sparse::{SparseEntry, SparseEntryMut, SparseFormatErrorKind}; use proptest::prelude::*; use proptest::sample::subsequence; +use super::test_data_examples::InvalidCsrDataExamples; + use crate::assert_panics; use crate::common::csr_strategy; @@ -171,11 +173,36 @@ fn csr_matrix_valid_data() { } } +#[test] +fn csr_matrix_valid_data_unsorted_column_indices() { + let csr = CsrMatrix::try_from_unsorted_csr_data( + 4, + 5, + vec![0, 3, 5, 8, 11], + vec![4, 1, 3, 3, 1, 2, 3, 0, 3, 4, 1], + vec![5, 1, 4, 7, 4, 2, 3, 1, 8, 9, 6], + ) + .unwrap(); + + let expected_csr = CsrMatrix::try_from_csr_data( + 4, + 5, + vec![0, 3, 5, 8, 11], + vec![1, 3, 4, 1, 3, 0, 2, 3, 1, 3, 4], + vec![1, 4, 5, 4, 7, 1, 2, 3, 6, 8, 9], + ) + .unwrap(); + + assert_eq!(csr, expected_csr); +} + #[test] fn csr_matrix_try_from_invalid_csr_data() { + let invalid_data: InvalidCsrDataExamples = InvalidCsrDataExamples::new(); { // Empty offset array (invalid length) - let matrix = CsrMatrix::try_from_csr_data(0, 0, Vec::new(), Vec::new(), Vec::::new()); + let (offsets, indices, values) = invalid_data.empty_offset_array; + let matrix = CsrMatrix::try_from_csr_data(0, 0, offsets, indices, values); assert_eq!( matrix.unwrap_err().kind(), &SparseFormatErrorKind::InvalidStructure @@ -184,10 +211,8 @@ fn csr_matrix_try_from_invalid_csr_data() { { // Offset array invalid length for arbitrary data - let offsets = vec![0, 3, 5]; - let indices = vec![0, 1, 2, 3, 5]; - let values = vec![0, 1, 2, 3, 4]; - + let (offsets, indices, values) = + invalid_data.offset_array_invalid_length_for_arbitrary_data; let matrix = CsrMatrix::try_from_csr_data(3, 6, offsets, indices, values); assert_eq!( matrix.unwrap_err().kind(), @@ -197,9 +222,7 @@ fn csr_matrix_try_from_invalid_csr_data() { { // Invalid first entry in offsets array - let offsets = vec![1, 2, 2, 5]; - let indices = vec![0, 5, 1, 2, 3]; - let values = vec![0, 1, 2, 3, 4]; + let (offsets, indices, values) = invalid_data.invalid_first_entry_in_offsets_array; let matrix = CsrMatrix::try_from_csr_data(3, 6, offsets, indices, values); assert_eq!( matrix.unwrap_err().kind(), @@ -209,9 +232,7 @@ fn csr_matrix_try_from_invalid_csr_data() { { // Invalid last entry in offsets array - let offsets = vec![0, 2, 2, 4]; - let indices = vec![0, 5, 1, 2, 3]; - let values = vec![0, 1, 2, 3, 4]; + let (offsets, indices, values) = invalid_data.invalid_last_entry_in_offsets_array; let matrix = CsrMatrix::try_from_csr_data(3, 6, offsets, indices, values); assert_eq!( matrix.unwrap_err().kind(), @@ -221,9 +242,7 @@ fn csr_matrix_try_from_invalid_csr_data() { { // Invalid length of offsets array - let offsets = vec![0, 2, 2]; - let indices = vec![0, 5, 1, 2, 3]; - let values = vec![0, 1, 2, 3, 4]; + let (offsets, indices, values) = invalid_data.invalid_length_of_offsets_array; let matrix = CsrMatrix::try_from_csr_data(3, 6, offsets, indices, values); assert_eq!( matrix.unwrap_err().kind(), @@ -233,9 +252,7 @@ fn csr_matrix_try_from_invalid_csr_data() { { // Nonmonotonic offsets - let offsets = vec![0, 3, 2, 5]; - let indices = vec![0, 1, 2, 3, 4]; - let values = vec![0, 1, 2, 3, 4]; + let (offsets, indices, values) = invalid_data.nonmonotonic_offsets; let matrix = CsrMatrix::try_from_csr_data(3, 6, offsets, indices, values); assert_eq!( matrix.unwrap_err().kind(), @@ -245,9 +262,7 @@ fn csr_matrix_try_from_invalid_csr_data() { { // Nonmonotonic minor indices - let offsets = vec![0, 2, 2, 5]; - let indices = vec![0, 2, 3, 1, 4]; - let values = vec![0, 1, 2, 3, 4]; + let (offsets, indices, values) = invalid_data.nonmonotonic_minor_indices; let matrix = CsrMatrix::try_from_csr_data(3, 6, offsets, indices, values); assert_eq!( matrix.unwrap_err().kind(), @@ -257,9 +272,7 @@ fn csr_matrix_try_from_invalid_csr_data() { { // Minor index out of bounds - let offsets = vec![0, 2, 2, 5]; - let indices = vec![0, 6, 1, 2, 3]; - let values = vec![0, 1, 2, 3, 4]; + let (offsets, indices, values) = invalid_data.minor_index_out_of_bounds; let matrix = CsrMatrix::try_from_csr_data(3, 6, offsets, indices, values); assert_eq!( matrix.unwrap_err().kind(), @@ -269,9 +282,7 @@ fn csr_matrix_try_from_invalid_csr_data() { { // Duplicate entry - let offsets = vec![0, 2, 2, 5]; - let indices = vec![0, 5, 2, 2, 3]; - let values = vec![0, 1, 2, 3, 4]; + let (offsets, indices, values) = invalid_data.duplicate_entry; let matrix = CsrMatrix::try_from_csr_data(3, 6, offsets, indices, values); assert_eq!( matrix.unwrap_err().kind(), @@ -280,6 +291,91 @@ fn csr_matrix_try_from_invalid_csr_data() { } } +#[test] +fn csr_matrix_try_from_unsorted_invalid_csr_data() { + let invalid_data: InvalidCsrDataExamples = InvalidCsrDataExamples::new(); + { + // Empty offset array (invalid length) + let (offsets, indices, values) = invalid_data.empty_offset_array; + let matrix = CsrMatrix::try_from_unsorted_csr_data(0, 0, offsets, indices, values); + assert_eq!( + matrix.unwrap_err().kind(), + &SparseFormatErrorKind::InvalidStructure + ); + } + + { + // Offset array invalid length for arbitrary data + let (offsets, indices, values) = + invalid_data.offset_array_invalid_length_for_arbitrary_data; + let matrix = CsrMatrix::try_from_unsorted_csr_data(3, 6, offsets, indices, values); + assert_eq!( + matrix.unwrap_err().kind(), + &SparseFormatErrorKind::InvalidStructure + ); + } + + { + // Invalid first entry in offsets array + let (offsets, indices, values) = invalid_data.invalid_first_entry_in_offsets_array; + let matrix = CsrMatrix::try_from_unsorted_csr_data(3, 6, offsets, indices, values); + assert_eq!( + matrix.unwrap_err().kind(), + &SparseFormatErrorKind::InvalidStructure + ); + } + + { + // Invalid last entry in offsets array + let (offsets, indices, values) = invalid_data.invalid_last_entry_in_offsets_array; + let matrix = CsrMatrix::try_from_unsorted_csr_data(3, 6, offsets, indices, values); + assert_eq!( + matrix.unwrap_err().kind(), + &SparseFormatErrorKind::InvalidStructure + ); + } + + { + // Invalid length of offsets array + let (offsets, indices, values) = invalid_data.invalid_length_of_offsets_array; + let matrix = CsrMatrix::try_from_unsorted_csr_data(3, 6, offsets, indices, values); + assert_eq!( + matrix.unwrap_err().kind(), + &SparseFormatErrorKind::InvalidStructure + ); + } + + { + // Nonmonotonic offsets + let (offsets, indices, values) = invalid_data.nonmonotonic_offsets; + let matrix = CsrMatrix::try_from_unsorted_csr_data(3, 6, offsets, indices, values); + assert_eq!( + matrix.unwrap_err().kind(), + &SparseFormatErrorKind::InvalidStructure + ); + } + + { + // Minor index out of bounds + let (offsets, indices, values) = invalid_data.minor_index_out_of_bounds; + let matrix = CsrMatrix::try_from_unsorted_csr_data(3, 6, offsets, indices, values); + assert_eq!( + matrix.unwrap_err().kind(), + &SparseFormatErrorKind::IndexOutOfBounds + ); + } + + { + // Duplicate entry + let (offsets, indices, values) = invalid_data.duplicate_entry; + let matrix = CsrMatrix::try_from_unsorted_csr_data(3, 6, offsets, indices, values); + assert_eq!( + matrix.unwrap_err().kind(), + &SparseFormatErrorKind::DuplicateEntry + ); + } +} + #[test] fn csr_disassemble_avoids_clone_when_owned() { // Test that disassemble avoids cloning the sparsity pattern when it holds the sole reference diff --git a/nalgebra-sparse/tests/unit_tests/matrix_market.rs b/nalgebra-sparse/tests/unit_tests/matrix_market.rs new file mode 100644 index 00000000..48ff1a78 --- /dev/null +++ b/nalgebra-sparse/tests/unit_tests/matrix_market.rs @@ -0,0 +1,345 @@ +use matrixcompare::assert_matrix_eq; +use nalgebra::dmatrix; +use nalgebra::Complex; +use nalgebra_sparse::io::load_coo_from_matrix_market_str; +use nalgebra_sparse::CooMatrix; + +#[test] +#[rustfmt::skip] +fn test_matrixmarket_sparse_real_general_empty() { + // Test several valid zero-shapes of a sparse matrix + let shapes = vec![ (0, 0), (1, 0), (0, 1) ]; + let strings: Vec = shapes + .iter() + .map(|(m, n)| format!("%%MatrixMarket matrix coordinate real general\n {} {} 0", m, n)) + .collect(); + + for (shape,string) in shapes.iter().zip(strings.iter()) { + let sparse_mat = load_coo_from_matrix_market_str::(string).unwrap(); + assert_eq!(sparse_mat.nrows(), shape.0); + assert_eq!(sparse_mat.ncols(), shape.1); + assert_eq!(sparse_mat.nnz(), 0); + } +} + +#[test] +#[rustfmt::skip] +fn test_matrixmarket_dense_real_general_empty() { + // Test several valid zero-shapes of a dense matrix + let shapes = vec![ (0, 0), (1, 0), (0, 1) ]; + let strings: Vec = shapes + .iter() + .map(|(m, n)| format!("%%MatrixMarket matrix array real general\n {} {}", m, n)) + .collect(); + + for (shape,string) in shapes.iter().zip(strings.iter()) { + let sparse_mat = load_coo_from_matrix_market_str::(string).unwrap(); + assert_eq!(sparse_mat.nrows(), shape.0); + assert_eq!(sparse_mat.ncols(), shape.1); + assert_eq!(sparse_mat.nnz(), 0); + } +} + +#[test] +#[rustfmt::skip] +fn test_matrixmarket_sparse_real_general() { + let file_str = r#" +%%MatrixMarket matrix CoOrdinate real general +% This is also an example of free-format features. +%================================================================================= +% +% This ASCII file represents a sparse MxN matrix with L +% nonzeros in the following Matrix Market format: +% +% +----------------------------------------------+ +% |%%MatrixMarket matrix coordinate real general | <--- header line +% |% | <--+ +% |% comments | |-- 0 or more comment lines +% |% | <--+ +% | M T L | <--- rows, columns, entries +% | I1 J1 A(I1, J1) | <--+ +% | I2 J2 A(I2, J2) | | +% | I3 J3 A(I3, J3) | |-- L lines +% | . . . | | +% | IL JL A(IL, JL) | <--+ +% +----------------------------------------------+ +% +% Indices are 1-based, i.e. A(1,1) is the first element. +% +%================================================================================= + 5 5 8 + 1 1 1 + + 2 2 1.050e+01 + 3 3 1.500e-02 + + + + + + + + + 1 4 6.000e+00 + + 4 2 2.505e+02 + +4 4 -2.800e+02 +4 5 3.332e+01 + 5 5 1.200e+01 +"#; + let sparse_mat = load_coo_from_matrix_market_str::(file_str).unwrap(); + let expected = dmatrix![ + 1.0, 0.0, 0.0, 6.0, 0.0; + 0.0, 10.5, 0.0, 0.0, 0.0; + 0.0, 0.0, 0.015, 0.0, 0.0; + 0.0, 250.5, 0.0, -280.0, 33.32; + 0.0, 0.0, 0.0, 0.0, 12.0; + ]; + assert_matrix_eq!(sparse_mat, expected); +} + +#[test] +#[rustfmt::skip] +fn test_matrixmarket_sparse_int_symmetric() { + let file_str = r#" +%%MatrixMarket matrix coordinate integer symmetric +% + 5 5 9 + 1 1 11 + 2 2 22 + 3 2 23 + 3 3 33 + 4 2 24 + 4 4 44 + 5 1 -15 + 5 3 35 + 5 5 55 +"#; + let sparse_mat = load_coo_from_matrix_market_str::(file_str).unwrap(); + let expected = dmatrix![ + 11, 0, 0, 0, -15; + 0, 22, 23, 24, 0; + 0, 23, 33, 0, 35; + 0, 24, 0, 44, 0; + -15, 0, 35, 0, 55; + ]; + assert_matrix_eq!(sparse_mat, expected); +} + +#[test] +#[rustfmt::skip] +fn test_matrixmarket_sparse_complex_hermitian() { + let file_str = r#" +%%MatrixMarket matrix coordinate complex hermitian +% + 5 5 7 + 1 1 1.0 0.0 + 2 2 10.5 0.0 + 4 2 250.5 22.22 + 3 3 0.015 0.0 + 4 4 -2.8e2 0.0 + 5 5 12.0 0.0 + 5 4 0.0 33.32 + +"#; + let sparse_mat = load_coo_from_matrix_market_str::>(file_str).unwrap(); + let expected = dmatrix![ + Complex::{re:1.0,im:0.0}, Complex::{re:0.0,im:0.0}, Complex::{re:0.0,im:0.0}, Complex::{re:0.0,im:0.0},Complex::{re:0.0,im:0.0}; + Complex::{re:0.0,im:0.0}, Complex::{re:10.5,im:0.0}, Complex::{re:0.0,im:0.0}, Complex::{re:250.5,im:-22.22},Complex::{re:0.0,im:0.0}; + Complex::{re:0.0,im:0.0}, Complex::{re:0.0,im:0.0}, Complex::{re:0.015,im:0.0}, Complex::{re:0.0,im:0.0},Complex::{re:0.0,im:0.0}; + Complex::{re:0.0,im:0.0}, Complex::{re:250.5,im:22.22}, Complex::{re:0.0,im:0.0}, Complex::{re:-280.0,im:0.0},Complex::{re:0.0,im:-33.32}; + Complex::{re:0.0,im:0.0}, Complex::{re:0.0,im:0.0}, Complex::{re:0.0,im:0.0}, Complex::{re:0.0,im:33.32},Complex::{re:12.0,im:0.0}; + ]; + assert_matrix_eq!(sparse_mat, expected); +} + +#[test] +#[rustfmt::skip] +fn test_matrixmarket_sparse_real_skew() { + let file_str = r#" +%%MatrixMarket matrix coordinate real skew-symmetric +% + 5 5 4 + 3 2 -23.0 + 4 2 -24.0 + 5 1 -15.0 + 5 3 -35.0 +"#; + let sparse_mat = load_coo_from_matrix_market_str::(file_str).unwrap(); + let expected = dmatrix![ + 0.0, 0.0, 0.0, 0.0, 15.0; + 0.0, 0.0, 23.0, 24.0, 0.0; + 0.0, -23.0, 0.0, 0.0, 35.0; + 0.0, -24.0, 0.0, 0.0, 0.0; + -15.0, 0.0, -35.0, 0.0, 0.0; + ]; + assert_matrix_eq!(sparse_mat, expected); +} + +#[test] +#[rustfmt::skip] +fn test_matrixmarket_sparse_pattern_general() { + let file_str = r#" +%%MatrixMarket matrix coordinate pattern general +% + 5 5 10 + 1 1 + 1 5 + 2 3 + 2 4 + 3 2 + 3 5 + 4 1 + 5 2 + 5 4 + 5 5 +"#; + let pattern_matrix = load_coo_from_matrix_market_str::<()>(file_str).unwrap(); + let nrows = pattern_matrix.nrows(); + let ncols = pattern_matrix.ncols(); + let (row_idx, col_idx, val) = pattern_matrix.disassemble(); + let values = vec![1; val.len()]; + let sparse_mat = CooMatrix::try_from_triplets(nrows, ncols, row_idx, col_idx, values).unwrap(); + let expected = dmatrix![ + 1, 0, 0, 0, 1; + 0, 0, 1, 1, 0; + 0, 1, 0, 0, 1; + 1, 0, 0, 0, 0; + 0, 1, 0, 1, 1; + ]; + assert_matrix_eq!(sparse_mat, expected); +} + +#[test] +#[rustfmt::skip] +fn test_matrixmarket_dense_real_general() { + let file_str = r#" +%%MatrixMarket matrix array real general +% +4 3 +1.0 +2.0 +3.0 +4.0 +5.0 +6.0 +7.0 +8.0 +9.0 +10.0 +11.0 +12.0 + +"#; + let sparse_mat = load_coo_from_matrix_market_str::(file_str).unwrap(); + let expected = dmatrix![ + 1.0, 5.0, 9.0; + 2.0, 6.0, 10.0; + 3.0, 7.0, 11.0; + 4.0, 8.0, 12.0; + ]; + assert_matrix_eq!(sparse_mat, expected); +} + +#[test] +#[rustfmt::skip] +fn test_matrixmarket_dense_real_symmetric() { + let file_str = r#" +%%MatrixMarket matrix array real symmetric +% +4 4 +1.0 +2.0 +3.0 +4.0 +5.0 +6.0 +7.0 +8.0 +9.0 +10.0 + +"#; + let sparse_mat = load_coo_from_matrix_market_str::(file_str).unwrap(); + let expected = dmatrix![ + 1.0, 2.0, 3.0, 4.0; + 2.0, 5.0, 6.0, 7.0; + 3.0, 6.0, 8.0, 9.0; + 4.0, 7.0, 9.0, 10.0; + ]; + assert_matrix_eq!(sparse_mat, expected); +} + +#[test] +#[rustfmt::skip] +fn test_matrixmarket_dense_complex_hermitian() { + let file_str = r#" +%%MatrixMarket matrix array complex hermitian +% +4 4 +1.0 0.0 +2.0 2.0 +3.0 3.0 +4.0 4.0 +5.0 0.0 +6.0 6.0 +7.0 7.0 +8.0 0.0 +9.0 9.0 +10.0 0.0 + +"#; + let sparse_mat = load_coo_from_matrix_market_str::>(file_str).unwrap(); + let expected = dmatrix![ + Complex::{re:1.0,im:0.0}, Complex::{re:2.0,im:-2.0} ,Complex::{re:3.0,im:-3.0} ,Complex::{re:4.0,im:-4.0}; + Complex::{re:2.0,im:2.0}, Complex::{re:5.0,im:0.0} ,Complex::{re:6.0,im:-6.0} ,Complex::{re:7.0,im:-7.0}; + Complex::{re:3.0,im:3.0}, Complex::{re:6.0,im:6.0} ,Complex::{re:8.0,im:0.0} ,Complex::{re:9.0,im:-9.0}; + Complex::{re:4.0,im:4.0}, Complex::{re:7.0,im:7.0} ,Complex::{re:9.0,im:9.0} ,Complex::{re:10.0,im:0.0}; + ]; + assert_matrix_eq!(sparse_mat, expected); +} + +#[test] +#[rustfmt::skip] +fn test_matrixmarket_dense_int_skew() { + let file_str = r#" +%%MatrixMarket matrix array integer skew-symmetric +% +4 4 +1 +2 +3 +4 +5 +6 +"#; + let sparse_mat = load_coo_from_matrix_market_str::(file_str).unwrap(); + let expected = dmatrix![ + 0,-1,-2,-3; + 1, 0,-4,-5; + 2, 4, 0,-6; + 3, 5, 6, 0; + ]; + assert_matrix_eq!(sparse_mat, expected); +} + +#[test] +#[rustfmt::skip] +fn test_matrixmarket_dense_complex_general() { + let file_str = r#" +%%MatrixMarket matrix array complex general +% +2 2 +1 0 +1 0 +1 0 +1 0 +"#; + let sparse_mat = load_coo_from_matrix_market_str::>(file_str).unwrap(); + let expected = dmatrix![ + Complex::{re:1.0,im:0.0},Complex::{re:1.0,im:0.0}; + Complex::{re:1.0,im:0.0},Complex::{re:1.0,im:0.0}; + ]; + assert_matrix_eq!(sparse_mat, expected); +} diff --git a/nalgebra-sparse/tests/unit_tests/mod.rs b/nalgebra-sparse/tests/unit_tests/mod.rs index ee2166dc..7090a493 100644 --- a/nalgebra-sparse/tests/unit_tests/mod.rs +++ b/nalgebra-sparse/tests/unit_tests/mod.rs @@ -3,6 +3,8 @@ mod convert_serial; mod coo; mod csc; mod csr; +mod matrix_market; mod ops; mod pattern; mod proptest; +mod test_data_examples; diff --git a/nalgebra-sparse/tests/unit_tests/test_data_examples.rs b/nalgebra-sparse/tests/unit_tests/test_data_examples.rs new file mode 100644 index 00000000..20721087 --- /dev/null +++ b/nalgebra-sparse/tests/unit_tests/test_data_examples.rs @@ -0,0 +1,44 @@ +/// Examples of *invalid* raw CSR data `(offsets, indices, values)`. +pub struct InvalidCsrDataExamples { + pub empty_offset_array: (Vec, Vec, Vec), + pub offset_array_invalid_length_for_arbitrary_data: (Vec, Vec, Vec), + pub invalid_first_entry_in_offsets_array: (Vec, Vec, Vec), + pub invalid_last_entry_in_offsets_array: (Vec, Vec, Vec), + pub invalid_length_of_offsets_array: (Vec, Vec, Vec), + pub nonmonotonic_offsets: (Vec, Vec, Vec), + pub nonmonotonic_minor_indices: (Vec, Vec, Vec), + pub minor_index_out_of_bounds: (Vec, Vec, Vec), + pub duplicate_entry: (Vec, Vec, Vec), +} + +impl InvalidCsrDataExamples { + pub fn new() -> Self { + let empty_offset_array = (Vec::::new(), Vec::::new(), Vec::::new()); + let offset_array_invalid_length_for_arbitrary_data = + (vec![0, 3, 5], vec![0, 1, 2, 3, 5], vec![0, 1, 2, 3, 4]); + let invalid_first_entry_in_offsets_array = + (vec![1, 2, 2, 5], vec![0, 5, 1, 2, 3], vec![0, 1, 2, 3, 4]); + let invalid_last_entry_in_offsets_array = + (vec![0, 2, 2, 4], vec![0, 5, 1, 2, 3], vec![0, 1, 2, 3, 4]); + let invalid_length_of_offsets_array = + (vec![0, 2, 2], vec![0, 5, 1, 2, 3], vec![0, 1, 2, 3, 4]); + let nonmonotonic_offsets = (vec![0, 3, 2, 5], vec![0, 1, 2, 3, 4], vec![0, 1, 2, 3, 4]); + let nonmonotonic_minor_indices = + (vec![0, 2, 2, 5], vec![0, 2, 3, 1, 4], vec![0, 1, 2, 3, 4]); + let minor_index_out_of_bounds = + (vec![0, 2, 2, 5], vec![0, 6, 1, 2, 3], vec![0, 1, 2, 3, 4]); + let duplicate_entry = (vec![0, 2, 2, 5], vec![0, 5, 2, 2, 3], vec![0, 1, 2, 3, 4]); + + return Self { + empty_offset_array, + offset_array_invalid_length_for_arbitrary_data, + invalid_first_entry_in_offsets_array, + invalid_last_entry_in_offsets_array, + invalid_length_of_offsets_array, + nonmonotonic_minor_indices, + nonmonotonic_offsets, + minor_index_out_of_bounds, + duplicate_entry, + }; + } +} diff --git a/src/base/array_storage.rs b/src/base/array_storage.rs index 3fc88ade..b46d442f 100644 --- a/src/base/array_storage.rs +++ b/src/base/array_storage.rs @@ -32,6 +32,10 @@ use std::mem; /// A array-based statically sized matrix data storage. #[repr(transparent)] #[derive(Copy, Clone, PartialEq, Eq, Hash)] +#[cfg_attr( + all(not(target_os = "cuda"), feature = "cuda"), + derive(cust::DeviceCopy) +)] pub struct ArrayStorage(pub [[T; R]; C]); impl ArrayStorage { diff --git a/src/base/conversion.rs b/src/base/conversion.rs index 46747f0e..dd71186f 100644 --- a/src/base/conversion.rs +++ b/src/base/conversion.rs @@ -20,10 +20,10 @@ use crate::base::{ MatrixSliceMut, OMatrix, Scalar, }; #[cfg(any(feature = "std", feature = "alloc"))] -use crate::base::{DVector, VecStorage}; +use crate::base::{DVector, RowDVector, VecStorage}; use crate::base::{SliceStorage, SliceStorageMut}; use crate::constraint::DimEq; -use crate::{IsNotStaticOne, RowSVector, SMatrix, SVector}; +use crate::{IsNotStaticOne, RowSVector, SMatrix, SVector, VectorSlice, VectorSliceMut}; use std::mem::MaybeUninit; // TODO: too bad this won't work for slice conversions. @@ -125,6 +125,24 @@ impl From> for [T; D] { } } +impl<'a, T: Scalar, RStride: Dim, CStride: Dim, const D: usize> + From, RStride, CStride>> for [T; D] +{ + #[inline] + fn from(vec: VectorSlice<'a, T, Const, RStride, CStride>) -> Self { + vec.into_owned().into() + } +} + +impl<'a, T: Scalar, RStride: Dim, CStride: Dim, const D: usize> + From, RStride, CStride>> for [T; D] +{ + #[inline] + fn from(vec: VectorSliceMut<'a, T, Const, RStride, CStride>) -> Self { + vec.into_owned().into() + } +} + impl From<[T; D]> for RowSVector where Const: IsNotStaticOne, @@ -197,8 +215,26 @@ impl From<[[T; R]; C]> for SMatrix From> for [[T; R]; C] { #[inline] - fn from(vec: SMatrix) -> Self { - vec.data.0 + fn from(mat: SMatrix) -> Self { + mat.data.0 + } +} + +impl<'a, T: Scalar, RStride: Dim, CStride: Dim, const R: usize, const C: usize> + From, Const, RStride, CStride>> for [[T; R]; C] +{ + #[inline] + fn from(mat: MatrixSlice<'a, T, Const, Const, RStride, CStride>) -> Self { + mat.into_owned().into() + } +} + +impl<'a, T: Scalar, RStride: Dim, CStride: Dim, const R: usize, const C: usize> + From, Const, RStride, CStride>> for [[T; R]; C] +{ + #[inline] + fn from(mat: MatrixSliceMut<'a, T, Const, Const, RStride, CStride>) -> Self { + mat.into_owned().into() } } @@ -453,6 +489,14 @@ impl<'a, T: Scalar> From> for DVector { } } +#[cfg(any(feature = "std", feature = "alloc"))] +impl<'a, T: Scalar> From> for RowDVector { + #[inline] + fn from(vec: Vec) -> Self { + Self::from_vec(vec) + } +} + impl<'a, T: Scalar + Copy, R: Dim, C: Dim, S: RawStorage + IsContiguous> From<&'a Matrix> for &'a [T] { diff --git a/src/base/default_allocator.rs b/src/base/default_allocator.rs index b676b5e3..09197bbd 100644 --- a/src/base/default_allocator.rs +++ b/src/base/default_allocator.rs @@ -195,7 +195,7 @@ where unsafe fn reallocate_copy( rto: Const, cto: Const, - mut buf: >::Buffer, + buf: >::Buffer, ) -> ArrayStorage, RTO, CTO> { let mut res = , Const>>::allocate_uninit(rto, cto); @@ -226,7 +226,7 @@ where unsafe fn reallocate_copy( rto: Dynamic, cto: CTo, - mut buf: ArrayStorage, + buf: ArrayStorage, ) -> VecStorage, Dynamic, CTo> { let mut res = >::allocate_uninit(rto, cto); @@ -257,7 +257,7 @@ where unsafe fn reallocate_copy( rto: RTo, cto: Dynamic, - mut buf: ArrayStorage, + buf: ArrayStorage, ) -> VecStorage, RTo, Dynamic> { let mut res = >::allocate_uninit(rto, cto); diff --git a/src/base/dimension.rs b/src/base/dimension.rs index 8573dd59..86006f3d 100644 --- a/src/base/dimension.rs +++ b/src/base/dimension.rs @@ -13,6 +13,10 @@ use serde::{Deserialize, Deserializer, Serialize, Serializer}; /// Dim of dynamically-sized algebraic entities. #[derive(Clone, Copy, Eq, PartialEq, Debug)] +#[cfg_attr( + all(not(target_os = "cuda"), feature = "cuda"), + derive(cust::DeviceCopy) +)] pub struct Dynamic { value: usize, } @@ -55,7 +59,7 @@ impl IsNotStaticOne for Dynamic {} /// Trait implemented by any type that can be used as a dimension. This includes type-level /// integers and `Dynamic` (for dimensions not known at compile-time). -pub trait Dim: Any + Debug + Copy + PartialEq + Send + Sync { +pub unsafe trait Dim: Any + Debug + Copy + PartialEq + Send + Sync { #[inline(always)] fn is() -> bool { TypeId::of::() == TypeId::of::() @@ -74,7 +78,7 @@ pub trait Dim: Any + Debug + Copy + PartialEq + Send + Sync { fn from_usize(dim: usize) -> Self; } -impl Dim for Dynamic { +unsafe impl Dim for Dynamic { #[inline] fn try_to_usize() -> Option { None @@ -197,6 +201,10 @@ dim_ops!( ); #[derive(Debug, Copy, Clone, PartialEq, Eq, Hash)] +#[cfg_attr( + all(not(target_os = "cuda"), feature = "cuda"), + derive(cust::DeviceCopy) +)] pub struct Const; /// Trait implemented exclusively by type-level integers. @@ -270,7 +278,7 @@ pub trait ToTypenum { type Typenum: Unsigned; } -impl Dim for Const { +unsafe impl Dim for Const { fn try_to_usize() -> Option { Some(T) } diff --git a/src/base/edition.rs b/src/base/edition.rs index 9569294e..760b3950 100644 --- a/src/base/edition.rs +++ b/src/base/edition.rs @@ -14,7 +14,7 @@ use crate::base::{DefaultAllocator, Matrix, OMatrix, RowVector, Scalar, Vector}; use crate::{Storage, UninitMatrix}; use std::mem::MaybeUninit; -/// # Rows and columns extraction +/// # Triangular matrix extraction impl> Matrix { /// Extracts the upper triangular part of this matrix (including the diagonal). #[inline] @@ -41,7 +41,10 @@ impl> Matrix { res } +} +/// # Rows and columns extraction +impl> Matrix { /// Creates a new matrix by extracting the given set of rows from `self`. #[cfg(any(feature = "std", feature = "alloc"))] #[must_use] @@ -95,9 +98,7 @@ impl> Matrix { for (destination, source) in icols.enumerate() { // NOTE: this is basically a copy_frow but wrapping the values insnide of MaybeUninit. res.column_mut(destination) - .zip_apply(&self.column(*source), |out, e| { - *out = MaybeUninit::new(e.clone()) - }); + .zip_apply(&self.column(*source), |out, e| *out = MaybeUninit::new(e)); } // Safety: res is now fully initialized. @@ -1094,7 +1095,7 @@ unsafe fn compress_rows( if new_nrows == 0 || ncols == 0 { // The output matrix is empty, drop everything. - ptr::drop_in_place(data.as_mut()); + ptr::drop_in_place(data); return; } diff --git a/src/base/indexing.rs b/src/base/indexing.rs index 93f41ed3..2c691bd1 100644 --- a/src/base/indexing.rs +++ b/src/base/indexing.rs @@ -1,4 +1,5 @@ //! Indexing +#![allow(clippy::reversed_empty_ranges)] use crate::base::storage::{RawStorage, RawStorageMut}; use crate::base::{ @@ -43,7 +44,7 @@ impl DimRange for usize { #[test] fn dimrange_usize() { - assert_eq!(DimRange::contained_by(&0, Const::<0>), false); + assert!(!DimRange::contained_by(&0, Const::<0>)); assert!(DimRange::contained_by(&0, Const::<1>)); } @@ -68,8 +69,8 @@ impl DimRange for ops::Range { #[test] fn dimrange_range_usize() { - assert_eq!(DimRange::contained_by(&(0..0), Const::<0>), false); - assert_eq!(DimRange::contained_by(&(0..1), Const::<0>), false); + assert!(!DimRange::contained_by(&(0..0), Const::<0>)); + assert!(!DimRange::contained_by(&(0..1), Const::<0>)); assert!(DimRange::contained_by(&(0..1), Const::<1>)); assert!(DimRange::contained_by( &((usize::MAX - 1)..usize::MAX), @@ -110,8 +111,8 @@ impl DimRange for ops::RangeFrom { #[test] fn dimrange_rangefrom_usize() { - assert_eq!(DimRange::contained_by(&(0..), Const::<0>), false); - assert_eq!(DimRange::contained_by(&(0..), Const::<0>), false); + assert!(!DimRange::contained_by(&(0..), Const::<0>)); + assert!(!DimRange::contained_by(&(0..), Const::<0>)); assert!(DimRange::contained_by(&(0..), Const::<1>)); assert!(DimRange::contained_by( &((usize::MAX - 1)..), @@ -204,16 +205,16 @@ impl DimRange for ops::RangeInclusive { #[test] fn dimrange_rangeinclusive_usize() { - assert_eq!(DimRange::contained_by(&(0..=0), Const::<0>), false); + assert!(!DimRange::contained_by(&(0..=0), Const::<0>)); assert!(DimRange::contained_by(&(0..=0), Const::<1>)); - assert_eq!( - DimRange::contained_by(&(usize::MAX..=usize::MAX), Dynamic::new(usize::MAX)), - false - ); - assert_eq!( - DimRange::contained_by(&((usize::MAX - 1)..=usize::MAX), Dynamic::new(usize::MAX)), - false - ); + assert!(!DimRange::contained_by( + &(usize::MAX..=usize::MAX), + Dynamic::new(usize::MAX) + )); + assert!(!DimRange::contained_by( + &((usize::MAX - 1)..=usize::MAX), + Dynamic::new(usize::MAX) + )); assert!(DimRange::contained_by( &((usize::MAX - 1)..=(usize::MAX - 1)), Dynamic::new(usize::MAX) @@ -255,7 +256,7 @@ impl DimRange for ops::RangeTo { #[test] fn dimrange_rangeto_usize() { assert!(DimRange::contained_by(&(..0), Const::<0>)); - assert_eq!(DimRange::contained_by(&(..1), Const::<0>), false); + assert!(!DimRange::contained_by(&(..1), Const::<0>)); assert!(DimRange::contained_by(&(..0), Const::<1>)); assert!(DimRange::contained_by( &(..(usize::MAX - 1)), @@ -292,13 +293,13 @@ impl DimRange for ops::RangeToInclusive { #[test] fn dimrange_rangetoinclusive_usize() { - assert_eq!(DimRange::contained_by(&(..=0), Const::<0>), false); - assert_eq!(DimRange::contained_by(&(..=1), Const::<0>), false); + assert!(!DimRange::contained_by(&(..=0), Const::<0>)); + assert!(!DimRange::contained_by(&(..=1), Const::<0>)); assert!(DimRange::contained_by(&(..=0), Const::<1>)); - assert_eq!( - DimRange::contained_by(&(..=(usize::MAX)), Dynamic::new(usize::MAX)), - false - ); + assert!(!DimRange::contained_by( + &(..=(usize::MAX)), + Dynamic::new(usize::MAX) + )); assert!(DimRange::contained_by( &(..=(usize::MAX - 1)), Dynamic::new(usize::MAX) @@ -566,7 +567,10 @@ where #[doc(hidden)] #[inline(always)] unsafe fn get_unchecked(self, matrix: &'a Matrix) -> Self::Output { - matrix.data.get_unchecked_linear(self) + let nrows = matrix.shape().0; + let row = self % nrows; + let col = self / nrows; + matrix.data.get_unchecked(row, col) } } @@ -585,7 +589,10 @@ where where S: RawStorageMut, { - matrix.data.get_unchecked_linear_mut(self) + let nrows = matrix.shape().0; + let row = self % nrows; + let col = self / nrows; + matrix.data.get_unchecked_mut(row, col) } } diff --git a/src/base/matrix.rs b/src/base/matrix.rs index 4dccc439..652eace1 100644 --- a/src/base/matrix.rs +++ b/src/base/matrix.rs @@ -32,7 +32,7 @@ use crate::{ArrayStorage, SMatrix, SimdComplexField, Storage, UninitMatrix}; use crate::storage::IsContiguous; use crate::uninit::{Init, InitStatus, Uninit}; #[cfg(any(feature = "std", feature = "alloc"))] -use crate::{DMatrix, DVector, Dynamic, VecStorage}; +use crate::{DMatrix, DVector, Dynamic, RowDVector, VecStorage}; use std::mem::MaybeUninit; /// A square matrix. @@ -92,6 +92,7 @@ pub type MatrixCross = /// - [Interpolation `lerp`, `slerp`…](#interpolation) /// - [BLAS functions `gemv`, `gemm`, `syger`…](#blas-functions) /// - [Swizzling `xx`, `yxz`…](#swizzling) +/// - [Triangular matrix extraction `upper_triangle`, `lower_triangle`](#triangular-matrix-extraction) /// /// #### Statistics /// - [Common operations `row_sum`, `column_mean`, `variance`…](#common-statistics-operations) @@ -154,6 +155,10 @@ pub type MatrixCross = /// some concrete types for `T` and a compatible data storage type `S`). #[repr(C)] #[derive(Clone, Copy)] +#[cfg_attr( + all(not(target_os = "cuda"), feature = "cuda"), + derive(cust::DeviceCopy) +)] pub struct Matrix { /// The data storage that contains all the matrix components. Disappointed? /// @@ -188,17 +193,14 @@ pub struct Matrix { // Note that it would probably make sense to just have // the type `Matrix`, and have `T, R, C` be associated-types // of the `RawStorage` trait. However, because we don't have - // specialization, this is not bossible because these `T, R, C` + // specialization, this is not possible because these `T, R, C` // allows us to desambiguate a lot of configurations. _phantoms: PhantomData<(T, R, C)>, } impl fmt::Debug for Matrix { fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> Result<(), fmt::Error> { - formatter - .debug_struct("Matrix") - .field("data", &self.data) - .finish() + self.data.fmt(formatter) } } @@ -411,6 +413,21 @@ impl DVector { } } +// TODO: Consider removing/deprecating `from_vec_storage` once we are able to make +// `from_data` const fn compatible +#[cfg(any(feature = "std", feature = "alloc"))] +impl RowDVector { + /// Creates a new heap-allocated matrix from the given [`VecStorage`]. + /// + /// This method exists primarily as a workaround for the fact that `from_data` can not + /// work in `const fn` contexts. + pub const fn from_vec_storage(storage: VecStorage) -> Self { + // This is sound because the dimensions of the matrix and the storage are guaranteed + // to be the same + unsafe { Self::from_data_statically_unchecked(storage) } + } +} + impl UninitMatrix where DefaultAllocator: Allocator, @@ -681,7 +698,7 @@ impl> Matrix { #[inline] fn transpose_to_uninit( &self, - status: Status, + _status: Status, out: &mut Matrix, ) where Status: InitStatus, @@ -1377,7 +1394,7 @@ impl> Matrix( &self, - status: Status, + _status: Status, out: &mut Matrix, ) where Status: InitStatus, @@ -1777,7 +1794,7 @@ where assert!(self.shape() == other.shape()); self.iter() .zip(other.iter()) - .all(|(a, b)| a.ulps_eq(b, epsilon.clone(), max_ulps.clone())) + .all(|(a, b)| a.ulps_eq(b, epsilon.clone(), max_ulps)) } } @@ -2021,16 +2038,26 @@ impl + SameNumberOfColumns, { - assert!( - self.shape() == (2, 1), - "2D perpendicular product requires (2, 1) vector but found {:?}", - self.shape() + let shape = self.shape(); + assert_eq!( + shape, + b.shape(), + "2D vector perpendicular product dimension mismatch." + ); + assert_eq!( + shape, + (2, 1), + "2D perpendicular product requires (2, 1) vectors {:?}", + shape ); - unsafe { - self.get_unchecked((0, 0)).clone() * b.get_unchecked((1, 0)).clone() - - self.get_unchecked((1, 0)).clone() * b.get_unchecked((0, 0)).clone() - } + // SAFETY: assertion above ensures correct shape + let ax = unsafe { self.get_unchecked((0, 0)).clone() }; + let ay = unsafe { self.get_unchecked((1, 0)).clone() }; + let bx = unsafe { b.get_unchecked((0, 0)).clone() }; + let by = unsafe { b.get_unchecked((1, 0)).clone() }; + + ax * by - ay * bx } // TODO: use specialization instead of an assertion. @@ -2051,17 +2078,14 @@ impl::from_usize(3); - let ncols = SameShapeC::::from_usize(1); - let mut res = Matrix::uninit(nrows, ncols); + let mut res = Matrix::uninit(Dim::from_usize(3), Dim::from_usize(1)); let ax = self.get_unchecked((0, 0)); let ay = self.get_unchecked((1, 0)); @@ -2083,10 +2107,7 @@ impl::from_usize(1); - let ncols = SameShapeC::::from_usize(3); - let mut res = Matrix::uninit(nrows, ncols); + let mut res = Matrix::uninit(Dim::from_usize(1), Dim::from_usize(3)); let ax = self.get_unchecked((0, 0)); let ay = self.get_unchecked((0, 1)); diff --git a/src/base/matrix_simba.rs b/src/base/matrix_simba.rs index 5c259207..6284ab56 100644 --- a/src/base/matrix_simba.rs +++ b/src/base/matrix_simba.rs @@ -42,14 +42,14 @@ where #[inline] fn replace(&mut self, i: usize, val: Self::Element) { - self.zip_apply(&val, |mut a, b| { + self.zip_apply(&val, |a, b| { a.replace(i, b); }) } #[inline] unsafe fn replace_unchecked(&mut self, i: usize, val: Self::Element) { - self.zip_apply(&val, |mut a, b| { + self.zip_apply(&val, |a, b| { a.replace_unchecked(i, b); }) } diff --git a/src/base/matrix_slice.rs b/src/base/matrix_slice.rs index 261d41e2..5fbd4b01 100644 --- a/src/base/matrix_slice.rs +++ b/src/base/matrix_slice.rs @@ -73,7 +73,6 @@ macro_rules! slice_storage_impl( S: $Storage, RStride: Dim, CStride: Dim { - $T::from_raw_parts(storage.$get_addr(start.0, start.1), shape, strides) } } diff --git a/src/base/min_max.rs b/src/base/min_max.rs index 0876fe67..324da0a0 100644 --- a/src/base/min_max.rs +++ b/src/base/min_max.rs @@ -60,7 +60,7 @@ impl> Matrix { T: SimdPartialOrd + Zero, { self.fold_with( - |e| e.map(|e| e.clone()).unwrap_or_else(T::zero), + |e| e.cloned().unwrap_or_else(T::zero), |a, b| a.simd_max(b.clone()), ) } @@ -123,7 +123,7 @@ impl> Matrix { T: SimdPartialOrd + Zero, { self.fold_with( - |e| e.map(|e| e.clone()).unwrap_or_else(T::zero), + |e| e.cloned().unwrap_or_else(T::zero), |a, b| a.simd_min(b.clone()), ) } diff --git a/src/base/norm.rs b/src/base/norm.rs index 3968885b..11f5661f 100644 --- a/src/base/norm.rs +++ b/src/base/norm.rs @@ -434,12 +434,7 @@ impl> Matrix { { let n = self.norm(); let le = n.clone().simd_le(min_norm); - self.apply(|e| { - *e = e - .clone() - .simd_unscale(n.clone()) - .select(le.clone(), e.clone()) - }); + self.apply(|e| *e = e.clone().simd_unscale(n.clone()).select(le, e.clone())); SimdOption::new(n, le) } diff --git a/src/base/ops.rs b/src/base/ops.rs index 5608119e..6b6d0c45 100644 --- a/src/base/ops.rs +++ b/src/base/ops.rs @@ -146,7 +146,7 @@ macro_rules! componentwise_binop_impl( #[inline] fn $method_to_statically_unchecked_uninit(&self, - status: Status, + _status: Status, rhs: &Matrix, out: &mut Matrix) where Status: InitStatus, @@ -699,7 +699,7 @@ where #[inline(always)] fn xx_mul_to_uninit( &self, - status: Status, + _status: Status, rhs: &Matrix, out: &mut Matrix, dot: impl Fn( diff --git a/src/base/properties.rs b/src/base/properties.rs index 7536a4a5..914c0b9e 100644 --- a/src/base/properties.rs +++ b/src/base/properties.rs @@ -7,10 +7,10 @@ use simba::scalar::{ClosedAdd, ClosedMul, ComplexField, RealField}; use crate::base::allocator::Allocator; use crate::base::dimension::{Dim, DimMin}; use crate::base::storage::Storage; -use crate::base::{DefaultAllocator, Matrix, Scalar, SquareMatrix}; +use crate::base::{DefaultAllocator, Matrix, SquareMatrix}; use crate::RawStorage; -impl> Matrix { +impl> Matrix { /// The total number of elements of this matrix. /// /// # Examples: @@ -63,50 +63,18 @@ impl> Matrix { T::Epsilon: Clone, { let (nrows, ncols) = self.shape(); - let d; - - if nrows > ncols { - d = ncols; - - for i in d..nrows { - for j in 0..ncols { - if !relative_eq!(self[(i, j)], T::zero(), epsilon = eps.clone()) { - return false; - } - } - } - } else { - // nrows <= ncols - d = nrows; + for j in 0..ncols { for i in 0..nrows { - for j in d..ncols { - if !relative_eq!(self[(i, j)], T::zero(), epsilon = eps.clone()) { - return false; - } - } - } - } - - // Off-diagonal elements of the sub-square matrix. - for i in 1..d { - for j in 0..i { - // TODO: use unsafe indexing. - if !relative_eq!(self[(i, j)], T::zero(), epsilon = eps.clone()) - || !relative_eq!(self[(j, i)], T::zero(), epsilon = eps.clone()) + let el = unsafe { self.get_unchecked((i, j)) }; + if (i == j && !relative_eq!(*el, T::one(), epsilon = eps.clone())) + || (i != j && !relative_eq!(*el, T::zero(), epsilon = eps.clone())) { return false; } } } - // Diagonal elements of the sub-square matrix. - for i in 0..d { - if !relative_eq!(self[(i, i)], T::one(), epsilon = eps.clone()) { - return false; - } - } - true } } diff --git a/src/base/statistics.rs b/src/base/statistics.rs index fc623c29..320cd12f 100644 --- a/src/base/statistics.rs +++ b/src/base/statistics.rs @@ -1,8 +1,8 @@ use crate::allocator::Allocator; use crate::storage::RawStorage; use crate::{Const, DefaultAllocator, Dim, Matrix, OVector, RowOVector, Scalar, VectorSlice, U1}; -use num::Zero; -use simba::scalar::{ClosedAdd, Field, SupersetOf}; +use num::{One, Zero}; +use simba::scalar::{ClosedAdd, ClosedMul, Field, SupersetOf}; use std::mem::MaybeUninit; /// # Folding on columns and rows @@ -123,7 +123,9 @@ impl> Matrix { /// 4.0, 5.0, 6.0); /// assert_eq!(m.row_sum(), RowVector3::new(5.0, 7.0, 9.0)); /// - /// let mint = Matrix3x2::new(1,2,3,4,5,6); + /// let mint = Matrix3x2::new(1, 2, + /// 3, 4, + /// 5, 6); /// assert_eq!(mint.row_sum(), RowVector2::new(9,12)); /// ``` #[inline] @@ -148,8 +150,10 @@ impl> Matrix { /// 4.0, 5.0, 6.0); /// assert_eq!(m.row_sum_tr(), Vector3::new(5.0, 7.0, 9.0)); /// - /// let mint = Matrix3x2::new(1,2,3,4,5,6); - /// assert_eq!(mint.row_sum_tr(), Vector2::new(9,12)); + /// let mint = Matrix3x2::new(1, 2, + /// 3, 4, + /// 5, 6); + /// assert_eq!(mint.row_sum_tr(), Vector2::new(9, 12)); /// ``` #[inline] #[must_use] @@ -173,8 +177,10 @@ impl> Matrix { /// 4.0, 5.0, 6.0); /// assert_eq!(m.column_sum(), Vector2::new(6.0, 15.0)); /// - /// let mint = Matrix3x2::new(1,2,3,4,5,6); - /// assert_eq!(mint.column_sum(), Vector3::new(3,7,11)); + /// let mint = Matrix3x2::new(1, 2, + /// 3, 4, + /// 5, 6); + /// assert_eq!(mint.column_sum(), Vector3::new(3, 7, 11)); /// ``` #[inline] #[must_use] @@ -189,6 +195,120 @@ impl> Matrix { }) } + /* + * + * Product computation. + * + */ + /// The product of all the elements of this matrix. + /// + /// # Example + /// + /// ``` + /// # use nalgebra::Matrix2x3; + /// + /// let m = Matrix2x3::new(1.0, 2.0, 3.0, + /// 4.0, 5.0, 6.0); + /// assert_eq!(m.product(), 720.0); + /// ``` + #[inline] + #[must_use] + pub fn product(&self) -> T + where + T: ClosedMul + One, + { + self.iter().cloned().fold(T::one(), |a, b| a * b) + } + + /// The product of all the rows of this matrix. + /// + /// Use `.row_sum_tr` if you need the result in a column vector instead. + /// + /// # Example + /// + /// ``` + /// # use nalgebra::{Matrix2x3, Matrix3x2}; + /// # use nalgebra::{RowVector2, RowVector3}; + /// + /// let m = Matrix2x3::new(1.0, 2.0, 3.0, + /// 4.0, 5.0, 6.0); + /// assert_eq!(m.row_product(), RowVector3::new(4.0, 10.0, 18.0)); + /// + /// let mint = Matrix3x2::new(1, 2, + /// 3, 4, + /// 5, 6); + /// assert_eq!(mint.row_product(), RowVector2::new(15, 48)); + /// ``` + #[inline] + #[must_use] + pub fn row_product(&self) -> RowOVector + where + T: ClosedMul + One, + DefaultAllocator: Allocator, + { + self.compress_rows(|col| col.product()) + } + + /// The product of all the rows of this matrix. The result is transposed and returned as a column vector. + /// + /// # Example + /// + /// ``` + /// # use nalgebra::{Matrix2x3, Matrix3x2}; + /// # use nalgebra::{Vector2, Vector3}; + /// + /// let m = Matrix2x3::new(1.0, 2.0, 3.0, + /// 4.0, 5.0, 6.0); + /// assert_eq!(m.row_product_tr(), Vector3::new(4.0, 10.0, 18.0)); + /// + /// let mint = Matrix3x2::new(1, 2, + /// 3, 4, + /// 5, 6); + /// assert_eq!(mint.row_product_tr(), Vector2::new(15, 48)); + /// ``` + #[inline] + #[must_use] + pub fn row_product_tr(&self) -> OVector + where + T: ClosedMul + One, + DefaultAllocator: Allocator, + { + self.compress_rows_tr(|col| col.product()) + } + + /// The product of all the columns of this matrix. + /// + /// # Example + /// + /// ``` + /// # use nalgebra::{Matrix2x3, Matrix3x2}; + /// # use nalgebra::{Vector2, Vector3}; + /// + /// let m = Matrix2x3::new(1.0, 2.0, 3.0, + /// 4.0, 5.0, 6.0); + /// assert_eq!(m.column_product(), Vector2::new(6.0, 120.0)); + /// + /// let mint = Matrix3x2::new(1, 2, + /// 3, 4, + /// 5, 6); + /// assert_eq!(mint.column_product(), Vector3::new(2, 12, 30)); + /// ``` + #[inline] + #[must_use] + pub fn column_product(&self) -> OVector + where + T: ClosedMul + One, + DefaultAllocator: Allocator, + { + let nrows = self.shape_generic().0; + self.compress_columns( + OVector::repeat_generic(nrows, Const::<1>, T::one()), + |out, col| { + out.component_mul_assign(&col); + }, + ) + } + /* * * Variance computation. diff --git a/src/base/uninit.rs b/src/base/uninit.rs index 92d246df..ad2759eb 100644 --- a/src/base/uninit.rs +++ b/src/base/uninit.rs @@ -66,11 +66,11 @@ unsafe impl InitStatus for Uninit { #[inline(always)] unsafe fn assume_init_ref(t: &MaybeUninit) -> &T { - std::mem::transmute(t.as_ptr()) // TODO: use t.assume_init_ref() + &*t.as_ptr() // TODO: use t.assume_init_ref() } #[inline(always)] unsafe fn assume_init_mut(t: &mut MaybeUninit) -> &mut T { - std::mem::transmute(t.as_mut_ptr()) // TODO: use t.assume_init_mut() + &mut *t.as_mut_ptr() // TODO: use t.assume_init_mut() } } diff --git a/src/base/unit.rs b/src/base/unit.rs index cd32b44b..60281b8f 100644 --- a/src/base/unit.rs +++ b/src/base/unit.rs @@ -1,3 +1,4 @@ +use std::fmt; #[cfg(feature = "abomonation-serialize")] use std::io::{Result as IOResult, Write}; use std::ops::Deref; @@ -24,11 +25,21 @@ use crate::{Dim, Matrix, OMatrix, RealField, Scalar, SimdComplexField, SimdRealF /// and [`UnitQuaternion`](crate::UnitQuaternion); both built on top of `Unit`. If you are interested /// in their documentation, read their dedicated pages directly. #[repr(transparent)] -#[derive(Clone, Hash, Debug, Copy)] +#[derive(Clone, Hash, Copy)] +// #[cfg_attr( +// all(not(target_os = "cuda"), feature = "cuda"), +// derive(cust::DeviceCopy) +// )] pub struct Unit { pub(crate) value: T, } +impl fmt::Debug for Unit { + fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> Result<(), fmt::Error> { + self.value.fmt(formatter) + } +} + #[cfg(feature = "bytemuck")] unsafe impl bytemuck::Zeroable for Unit where T: bytemuck::Zeroable {} @@ -111,6 +122,17 @@ mod rkyv_impl { } } +#[cfg(all(not(target_os = "cuda"), feature = "cuda"))] +unsafe impl cust::memory::DeviceCopy + for Unit> +where + T: Scalar, + R: Dim, + C: Dim, + S: RawStorage + Copy, +{ +} + impl PartialEq for Unit> where T: Scalar + PartialEq, diff --git a/src/geometry/dual_quaternion.rs b/src/geometry/dual_quaternion.rs index 11ff46d4..4280668a 100644 --- a/src/geometry/dual_quaternion.rs +++ b/src/geometry/dual_quaternion.rs @@ -39,6 +39,10 @@ use simba::scalar::{ClosedNeg, RealField}; /// See #[repr(C)] #[derive(Debug, Copy, Clone)] +#[cfg_attr( + all(not(target_os = "cuda"), feature = "cuda"), + derive(cust::DeviceCopy) +)] pub struct DualQuaternion { /// The real component of the quaternion pub real: Quaternion, @@ -351,13 +355,14 @@ impl> UlpsEq for DualQuaternion { #[inline] fn ulps_eq(&self, other: &Self, epsilon: Self::Epsilon, max_ulps: u32) -> bool { - self.clone().to_vector().ulps_eq(&other.clone().to_vector(), epsilon.clone(), max_ulps.clone()) || + self.clone().to_vector().ulps_eq(&other.clone().to_vector(), epsilon.clone(), max_ulps) || // Account for the double-covering of S², i.e. q = -q. - self.clone().to_vector().iter().zip(other.clone().to_vector().iter()).all(|(a, b)| a.ulps_eq(&-b.clone(), epsilon.clone(), max_ulps.clone())) + self.clone().to_vector().iter().zip(other.clone().to_vector().iter()).all(|(a, b)| a.ulps_eq(&-b.clone(), epsilon.clone(), max_ulps)) } } -/// A unit quaternions. May be used to represent a rotation followed by a translation. +/// A unit dual quaternion. May be used to represent a rotation followed by a +/// translation. pub type UnitDualQuaternion = Unit>; impl PartialEq for UnitDualQuaternion { @@ -593,8 +598,9 @@ where /// Screw linear interpolation between two unit quaternions. This creates a /// smooth arc from one dual-quaternion to another. /// - /// Panics if the angle between both quaternion is 180 degrees (in which case the interpolation - /// is not well-defined). Use `.try_sclerp` instead to avoid the panic. + /// Panics if the angle between both quaternion is 180 degrees (in which + /// case the interpolation is not well-defined). Use `.try_sclerp` + /// instead to avoid the panic. /// /// # Example /// ``` @@ -627,15 +633,16 @@ where .expect("DualQuaternion sclerp: ambiguous configuration.") } - /// Computes the screw-linear interpolation between two unit quaternions or returns `None` - /// if both quaternions are approximately 180 degrees apart (in which case the interpolation is - /// not well-defined). + /// Computes the screw-linear interpolation between two unit quaternions or + /// returns `None` if both quaternions are approximately 180 degrees + /// apart (in which case the interpolation is not well-defined). /// /// # Arguments /// * `self`: the first quaternion to interpolate from. /// * `other`: the second quaternion to interpolate toward. /// * `t`: the interpolation parameter. Should be between 0 and 1. - /// * `epsilon`: the value below which the sinus of the angle separating both quaternion + /// * `epsilon`: the value below which the sinus of the angle separating + /// both quaternion /// must be to return `None`. #[inline] #[must_use] @@ -650,6 +657,10 @@ where // interpolation. let other = { let dot_product = self.as_ref().real.coords.dot(&other.as_ref().real.coords); + if relative_eq!(dot_product, T::zero(), epsilon = epsilon.clone()) { + return None; + } + if dot_product < T::zero() { -other.clone() } else { @@ -660,13 +671,21 @@ where let difference = self.as_ref().conjugate() * other.as_ref(); let norm_squared = difference.real.vector().norm_squared(); if relative_eq!(norm_squared, T::zero(), epsilon = epsilon) { - return None; + return Some(Self::from_parts( + self.translation() + .vector + .lerp(&other.translation().vector, t) + .into(), + self.rotation(), + )); } - let inverse_norm_squared = T::one() / norm_squared; + let scalar: T = difference.real.scalar(); + let mut angle = two.clone() * scalar.acos(); + + let inverse_norm_squared: T = T::one() / norm_squared; let inverse_norm = inverse_norm_squared.sqrt(); - let mut angle = two.clone() * difference.real.scalar().acos(); let mut pitch = -two * difference.dual.scalar() * inverse_norm.clone(); let direction = difference.real.vector() * inverse_norm.clone(); let moment = (difference.dual.vector() @@ -678,6 +697,7 @@ where let sin = (half.clone() * angle.clone()).sin(); let cos = (half.clone() * angle).cos(); + let real = Quaternion::from_parts(cos.clone(), direction.clone() * sin.clone()); let dual = Quaternion::from_parts( -pitch.clone() * half.clone() * sin.clone(), diff --git a/src/geometry/isometry.rs b/src/geometry/isometry.rs index 4492c6c1..f020c0e9 100755 --- a/src/geometry/isometry.rs +++ b/src/geometry/isometry.rs @@ -54,7 +54,11 @@ use crate::geometry::{AbstractRotation, Point, Translation}; /// * [Conversion to a matrix `to_matrix`…](#conversion-to-a-matrix) /// #[repr(C)] -#[derive(Debug)] +#[derive(Debug, Copy, Clone)] +#[cfg_attr( + all(not(target_os = "cuda"), feature = "cuda"), + derive(cust::DeviceCopy) +)] #[cfg_attr(feature = "serde-serialize-no-std", derive(Serialize, Deserialize))] #[cfg_attr( feature = "serde-serialize-no-std", @@ -170,20 +174,6 @@ where } } -impl Copy for Isometry where - Owned>: Copy -{ -} - -impl Clone for Isometry { - #[inline] - fn clone(&self) -> Self { - Self { - rotation: self.rotation.clone(), - translation: self.translation.clone(), - } - } -} /// # From the translation and rotation parts impl, const D: usize> Isometry { /// Creates a new isometry from its rotational and translational parts. @@ -629,7 +619,7 @@ where #[inline] fn ulps_eq(&self, other: &Self, epsilon: Self::Epsilon, max_ulps: u32) -> bool { self.translation - .ulps_eq(&other.translation, epsilon.clone(), max_ulps.clone()) + .ulps_eq(&other.translation, epsilon.clone(), max_ulps) && self.rotation.ulps_eq(&other.rotation, epsilon, max_ulps) } } diff --git a/src/geometry/isometry_construction.rs b/src/geometry/isometry_construction.rs index 9b855599..cce2932d 100644 --- a/src/geometry/isometry_construction.rs +++ b/src/geometry/isometry_construction.rs @@ -21,6 +21,15 @@ use crate::{ UnitQuaternion, }; +impl, const D: usize> Default for Isometry +where + T::Element: SimdRealField, +{ + fn default() -> Self { + Self::identity() + } +} + impl, const D: usize> Isometry where T::Element: SimdRealField, diff --git a/src/geometry/mod.rs b/src/geometry/mod.rs index 37ca57f9..e7ff7304 100644 --- a/src/geometry/mod.rs +++ b/src/geometry/mod.rs @@ -48,6 +48,14 @@ mod translation_coordinates; mod translation_ops; mod translation_simba; +mod scale; +mod scale_alias; +mod scale_construction; +mod scale_conversion; +mod scale_coordinates; +mod scale_ops; +mod scale_simba; + mod isometry; mod isometry_alias; mod isometry_construction; @@ -95,6 +103,9 @@ pub use self::unit_complex::*; pub use self::translation::*; pub use self::translation_alias::*; +pub use self::scale::*; +pub use self::scale_alias::*; + pub use self::isometry::*; pub use self::isometry_alias::*; diff --git a/src/geometry/orthographic.rs b/src/geometry/orthographic.rs index 731b46a1..18a7852d 100644 --- a/src/geometry/orthographic.rs +++ b/src/geometry/orthographic.rs @@ -19,19 +19,15 @@ use crate::geometry::{Point3, Projective3}; /// A 3D orthographic projection stored as a homogeneous 4x4 matrix. #[repr(C)] +#[cfg_attr( + all(not(target_os = "cuda"), feature = "cuda"), + derive(cust::DeviceCopy) +)] +#[derive(Copy, Clone)] pub struct Orthographic3 { matrix: Matrix4, } -impl Copy for Orthographic3 {} - -impl Clone for Orthographic3 { - #[inline] - fn clone(&self) -> Self { - Self::from_matrix_unchecked(self.matrix.clone()) - } -} - impl fmt::Debug for Orthographic3 { fn fmt(&self, f: &mut fmt::Formatter<'_>) -> Result<(), fmt::Error> { self.matrix.fmt(f) @@ -175,7 +171,7 @@ impl Orthographic3 { ); let half: T = crate::convert(0.5); - let width = zfar.clone() * (vfov.clone() * half.clone()).tan(); + let width = zfar.clone() * (vfov * half.clone()).tan(); let height = width.clone() / aspect; Self::new( diff --git a/src/geometry/perspective.rs b/src/geometry/perspective.rs index 34af6f0b..59b7f9f2 100644 --- a/src/geometry/perspective.rs +++ b/src/geometry/perspective.rs @@ -20,19 +20,15 @@ use crate::geometry::{Point3, Projective3}; /// A 3D perspective projection stored as a homogeneous 4x4 matrix. #[repr(C)] +#[cfg_attr( + all(not(target_os = "cuda"), feature = "cuda"), + derive(cust::DeviceCopy) +)] +#[derive(Copy, Clone)] pub struct Perspective3 { matrix: Matrix4, } -impl Copy for Perspective3 {} - -impl Clone for Perspective3 { - #[inline] - fn clone(&self) -> Self { - Self::from_matrix_unchecked(self.matrix.clone()) - } -} - impl fmt::Debug for Perspective3 { fn fmt(&self, f: &mut fmt::Formatter<'_>) -> Result<(), fmt::Error> { self.matrix.fmt(f) diff --git a/src/geometry/point.rs b/src/geometry/point.rs index 69022671..0953a9ce 100644 --- a/src/geometry/point.rs +++ b/src/geometry/point.rs @@ -40,7 +40,7 @@ use std::mem::MaybeUninit; /// may have some other methods, e.g., `isometry.inverse_transform_point(&point)`. See the documentation /// of said transformations for details. #[repr(C)] -#[derive(Debug, Clone)] +#[derive(Clone)] pub struct OPoint where DefaultAllocator: Allocator, @@ -49,6 +49,15 @@ where pub coords: OVector, } +impl fmt::Debug for OPoint +where + DefaultAllocator: Allocator, +{ + fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> Result<(), fmt::Error> { + self.coords.as_slice().fmt(formatter) + } +} + impl hash::Hash for OPoint where DefaultAllocator: Allocator, @@ -65,6 +74,15 @@ where { } +#[cfg(all(not(target_os = "cuda"), feature = "cuda"))] +unsafe impl cust::memory::DeviceCopy + for OPoint +where + DefaultAllocator: Allocator, + OVector: cust::memory::DeviceCopy, +{ +} + #[cfg(feature = "bytemuck")] unsafe impl bytemuck::Zeroable for OPoint where diff --git a/src/geometry/point_construction.rs b/src/geometry/point_construction.rs index e4e729aa..2136080a 100644 --- a/src/geometry/point_construction.rs +++ b/src/geometry/point_construction.rs @@ -19,6 +19,15 @@ use simba::scalar::{ClosedDiv, SupersetOf}; use crate::geometry::Point; +impl Default for OPoint +where + DefaultAllocator: Allocator, +{ + fn default() -> Self { + Self::origin() + } +} + /// # Other construction methods impl OPoint where diff --git a/src/geometry/quaternion.rs b/src/geometry/quaternion.rs index 0c2c01c7..cae06251 100755 --- a/src/geometry/quaternion.rs +++ b/src/geometry/quaternion.rs @@ -27,12 +27,22 @@ use crate::geometry::{Point3, Rotation}; /// A quaternion. See the type alias `UnitQuaternion = Unit` for a quaternion /// that may be used as a rotation. #[repr(C)] -#[derive(Debug, Copy, Clone)] +#[derive(Copy, Clone)] +#[cfg_attr( + all(not(target_os = "cuda"), feature = "cuda"), + derive(cust::DeviceCopy) +)] pub struct Quaternion { /// This quaternion as a 4D vector of coordinates in the `[ x, y, z, w ]` storage order. pub coords: Vector4, } +impl fmt::Debug for Quaternion { + fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> Result<(), fmt::Error> { + self.coords.as_slice().fmt(formatter) + } +} + impl Hash for Quaternion { fn hash(&self, state: &mut H) { self.coords.hash(state) @@ -1039,9 +1049,9 @@ impl> UlpsEq for Quaternion { #[inline] fn ulps_eq(&self, other: &Self, epsilon: Self::Epsilon, max_ulps: u32) -> bool { - self.as_vector().ulps_eq(other.as_vector(), epsilon.clone(), max_ulps.clone()) || + self.as_vector().ulps_eq(other.as_vector(), epsilon.clone(), max_ulps) || // Account for the double-covering of S², i.e. q = -q. - self.as_vector().iter().zip(other.as_vector().iter()).all(|(a, b)| a.ulps_eq(&-b.clone(), epsilon.clone(), max_ulps.clone())) + self.as_vector().iter().zip(other.as_vector().iter()).all(|(a, b)| a.ulps_eq(&-b.clone(), epsilon.clone(), max_ulps)) } } @@ -1058,6 +1068,9 @@ impl fmt::Display for Quaternion { /// A unit quaternions. May be used to represent a rotation. pub type UnitQuaternion = Unit>; +#[cfg(all(not(target_os = "cuda"), feature = "cuda"))] +unsafe impl cust::memory::DeviceCopy for UnitQuaternion {} + impl PartialEq for UnitQuaternion { #[inline] fn eq(&self, rhs: &Self) -> bool { @@ -1492,18 +1505,18 @@ where let wk = w.clone() * k.clone() * crate::convert(2.0f64); let wj = w.clone() * j.clone() * crate::convert(2.0f64); let ik = i.clone() * k.clone() * crate::convert(2.0f64); - let jk = j.clone() * k.clone() * crate::convert(2.0f64); - let wi = w.clone() * i.clone() * crate::convert(2.0f64); + let jk = j * k * crate::convert(2.0f64); + let wi = w * i * crate::convert(2.0f64); Rotation::from_matrix_unchecked(Matrix3::new( ww.clone() + ii.clone() - jj.clone() - kk.clone(), ij.clone() - wk.clone(), wj.clone() + ik.clone(), - wk.clone() + ij.clone(), + wk + ij, ww.clone() - ii.clone() + jj.clone() - kk.clone(), jk.clone() - wi.clone(), - ik.clone() - wj.clone(), - wi.clone() + jk.clone(), + ik - wj, + wi + jk, ww - ii - jj + kk, )) } diff --git a/src/geometry/rotation.rs b/src/geometry/rotation.rs index 3ac3ca57..2d5e5d24 100755 --- a/src/geometry/rotation.rs +++ b/src/geometry/rotation.rs @@ -54,11 +54,21 @@ use crate::geometry::Point; /// * [Conversion to a matrix `matrix`, `to_homogeneous`…](#conversion-to-a-matrix) /// #[repr(C)] -#[derive(Debug)] +#[cfg_attr( + all(not(target_os = "cuda"), feature = "cuda"), + derive(cust::DeviceCopy) +)] +#[derive(Copy, Clone)] pub struct Rotation { matrix: SMatrix, } +impl fmt::Debug for Rotation { + fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> Result<(), fmt::Error> { + self.matrix.fmt(formatter) + } +} + impl hash::Hash for Rotation where , Const>>::Buffer: hash::Hash, @@ -68,21 +78,6 @@ where } } -impl Copy for Rotation where - , Const>>::Buffer: Copy -{ -} - -impl Clone for Rotation -where - , Const>>::Buffer: Clone, -{ - #[inline] - fn clone(&self) -> Self { - Self::from_matrix_unchecked(self.matrix.clone()) - } -} - #[cfg(feature = "bytemuck")] unsafe impl bytemuck::Zeroable for Rotation where diff --git a/src/geometry/rotation_construction.rs b/src/geometry/rotation_construction.rs index c18c6a58..84e42693 100644 --- a/src/geometry/rotation_construction.rs +++ b/src/geometry/rotation_construction.rs @@ -6,6 +6,15 @@ use crate::base::{SMatrix, Scalar}; use crate::geometry::Rotation; +impl Default for Rotation +where + T: Scalar + Zero + One, +{ + fn default() -> Self { + Self::identity() + } +} + /// # Identity impl Rotation where @@ -15,9 +24,16 @@ where /// /// # Example /// ``` - /// # use nalgebra::Quaternion; - /// let rot1 = Quaternion::identity(); - /// let rot2 = Quaternion::new(1.0, 2.0, 3.0, 4.0); + /// # use nalgebra::{Rotation2, Rotation3}; + /// # use nalgebra::Vector3; + /// let rot1 = Rotation2::identity(); + /// let rot2 = Rotation2::new(std::f32::consts::FRAC_PI_2); + /// + /// assert_eq!(rot1 * rot2, rot2); + /// assert_eq!(rot2 * rot1, rot2); + /// + /// let rot1 = Rotation3::identity(); + /// let rot2 = Rotation3::from_axis_angle(&Vector3::z_axis(), std::f32::consts::FRAC_PI_2); /// /// assert_eq!(rot1 * rot2, rot2); /// assert_eq!(rot2 * rot1, rot2); diff --git a/src/geometry/rotation_specialization.rs b/src/geometry/rotation_specialization.rs index c24514ba..41405c87 100644 --- a/src/geometry/rotation_specialization.rs +++ b/src/geometry/rotation_specialization.rs @@ -60,7 +60,7 @@ impl Rotation2 { impl Rotation2 { /// Builds a rotation from a basis assumed to be orthonormal. /// - /// In order to get a valid unit-quaternion, the input must be an + /// In order to get a valid rotation matrix, the input must be an /// orthonormal basis, i.e., all vectors are normalized, and the are /// all orthogonal to each other. These invariants are not checked /// by this method. @@ -204,7 +204,7 @@ impl Rotation2 { *self = Self::from_matrix_eps(self.matrix(), T::default_epsilon(), 0, c.into()) } - /// Raise the quaternion to a given floating power, i.e., returns the rotation with the angle + /// Raise the rotation to a given floating power, i.e., returns the rotation with the angle /// of `self` multiplied by `n`. /// /// # Example @@ -660,7 +660,7 @@ where other * self.inverse() } - /// Raise the quaternion to a given floating power, i.e., returns the rotation with the same + /// Raise the rotation to a given floating power, i.e., returns the rotation with the same /// axis as `self` and an angle equal to `self.angle()` multiplied by `n`. /// /// # Example @@ -692,7 +692,7 @@ where /// Builds a rotation from a basis assumed to be orthonormal. /// - /// In order to get a valid unit-quaternion, the input must be an + /// In order to get a valid rotation matrix, the input must be an /// orthonormal basis, i.e., all vectors are normalized, and the are /// all orthogonal to each other. These invariants are not checked /// by this method. @@ -846,7 +846,7 @@ impl Rotation3 { } } - /// The rotation axis and angle in ]0, pi] of this unit quaternion. + /// The rotation axis and angle in ]0, pi] of this rotation matrix. /// /// Returns `None` if the angle is zero. /// diff --git a/src/geometry/scale.rs b/src/geometry/scale.rs new file mode 100755 index 00000000..b1d278d6 --- /dev/null +++ b/src/geometry/scale.rs @@ -0,0 +1,443 @@ +use approx::{AbsDiffEq, RelativeEq, UlpsEq}; +use num::{One, Zero}; +use std::fmt; +use std::hash; +#[cfg(feature = "abomonation-serialize")] +use std::io::{Result as IOResult, Write}; + +#[cfg(feature = "serde-serialize-no-std")] +use serde::{Deserialize, Deserializer, Serialize, Serializer}; + +#[cfg(feature = "abomonation-serialize")] +use abomonation::Abomonation; + +use crate::base::allocator::Allocator; +use crate::base::dimension::{DimNameAdd, DimNameSum, U1}; +use crate::base::storage::Owned; +use crate::base::{Const, DefaultAllocator, OMatrix, OVector, SVector, Scalar}; +use crate::ClosedDiv; +use crate::ClosedMul; + +use crate::geometry::Point; + +/// A scale which supports non-uniform scaling. +#[repr(C)] +#[cfg_attr( + all(not(target_os = "cuda"), feature = "cuda"), + derive(cust::DeviceCopy) +)] +#[derive(Copy, Clone)] +pub struct Scale { + /// The scale coordinates, i.e., how much is multiplied to a point's coordinates when it is + /// scaled. + pub vector: SVector, +} + +impl fmt::Debug for Scale { + fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> Result<(), fmt::Error> { + self.vector.as_slice().fmt(formatter) + } +} + +impl hash::Hash for Scale +where + Owned>: hash::Hash, +{ + fn hash(&self, state: &mut H) { + self.vector.hash(state) + } +} + +#[cfg(feature = "bytemuck")] +unsafe impl bytemuck::Zeroable for Scale +where + T: Scalar + bytemuck::Zeroable, + SVector: bytemuck::Zeroable, +{ +} + +#[cfg(feature = "bytemuck")] +unsafe impl bytemuck::Pod for Scale +where + T: Scalar + bytemuck::Pod, + SVector: bytemuck::Pod, +{ +} + +#[cfg(feature = "abomonation-serialize")] +impl Abomonation for Scale +where + T: Scalar, + SVector: Abomonation, +{ + unsafe fn entomb(&self, writer: &mut W) -> IOResult<()> { + self.vector.entomb(writer) + } + + fn extent(&self) -> usize { + self.vector.extent() + } + + unsafe fn exhume<'a, 'b>(&'a mut self, bytes: &'b mut [u8]) -> Option<&'b mut [u8]> { + self.vector.exhume(bytes) + } +} + +#[cfg(feature = "serde-serialize-no-std")] +impl Serialize for Scale +where + Owned>: Serialize, +{ + fn serialize(&self, serializer: S) -> Result + where + S: Serializer, + { + self.vector.serialize(serializer) + } +} + +#[cfg(feature = "serde-serialize-no-std")] +impl<'a, T: Scalar, const D: usize> Deserialize<'a> for Scale +where + Owned>: Deserialize<'a>, +{ + fn deserialize(deserializer: Des) -> Result + where + Des: Deserializer<'a>, + { + let matrix = SVector::::deserialize(deserializer)?; + + Ok(Scale::from(matrix)) + } +} + +#[cfg(feature = "rkyv-serialize-no-std")] +mod rkyv_impl { + use super::Scale; + use crate::base::SVector; + use rkyv::{offset_of, project_struct, Archive, Deserialize, Fallible, Serialize}; + + impl Archive for Scale { + type Archived = Scale; + type Resolver = as Archive>::Resolver; + + fn resolve( + &self, + pos: usize, + resolver: Self::Resolver, + out: &mut core::mem::MaybeUninit, + ) { + self.vector.resolve( + pos + offset_of!(Self::Archived, vector), + resolver, + project_struct!(out: Self::Archived => vector), + ); + } + } + + impl, S: Fallible + ?Sized, const D: usize> Serialize for Scale { + fn serialize(&self, serializer: &mut S) -> Result { + self.vector.serialize(serializer) + } + } + + impl Deserialize, _D> + for Scale + where + T::Archived: Deserialize, + { + fn deserialize(&self, deserializer: &mut _D) -> Result, _D::Error> { + Ok(Scale { + vector: self.vector.deserialize(deserializer)?, + }) + } + } +} + +impl Scale { + /// Inverts `self`. + /// + /// # Example + /// ``` + /// # use nalgebra::{Scale2, Scale3}; + /// let t = Scale3::new(1.0, 2.0, 3.0); + /// assert_eq!(t * t.try_inverse().unwrap(), Scale3::identity()); + /// assert_eq!(t.try_inverse().unwrap() * t, Scale3::identity()); + /// + /// // Work in all dimensions. + /// let t = Scale2::new(1.0, 2.0); + /// assert_eq!(t * t.try_inverse().unwrap(), Scale2::identity()); + /// assert_eq!(t.try_inverse().unwrap() * t, Scale2::identity()); + /// + /// // Returns None if any coordinate is 0. + /// let t = Scale2::new(0.0, 2.0); + /// assert_eq!(t.try_inverse(), None); + /// ``` + #[inline] + #[must_use = "Did you mean to use try_inverse_mut()?"] + pub fn try_inverse(&self) -> Option> + where + T: ClosedDiv + One + Zero, + { + for i in 0..D { + if self.vector[i] == T::zero() { + return None; + } + } + return Some(self.vector.map(|e| T::one() / e).into()); + } + + /// Inverts `self`. + /// + /// # Example + /// ``` + /// # use nalgebra::{Scale2, Scale3}; + /// + /// unsafe { + /// let t = Scale3::new(1.0, 2.0, 3.0); + /// assert_eq!(t * t.inverse_unchecked(), Scale3::identity()); + /// assert_eq!(t.inverse_unchecked() * t, Scale3::identity()); + /// + /// // Work in all dimensions. + /// let t = Scale2::new(1.0, 2.0); + /// assert_eq!(t * t.inverse_unchecked(), Scale2::identity()); + /// assert_eq!(t.inverse_unchecked() * t, Scale2::identity()); + /// } + /// ``` + #[inline] + #[must_use] + pub unsafe fn inverse_unchecked(&self) -> Scale + where + T: ClosedDiv + One, + { + return self.vector.map(|e| T::one() / e).into(); + } + + /// Inverts `self`. + /// + /// # Example + /// ``` + /// # use nalgebra::{Scale2, Scale3}; + /// let t = Scale3::new(1.0, 2.0, 3.0); + /// assert_eq!(t * t.pseudo_inverse(), Scale3::identity()); + /// assert_eq!(t.pseudo_inverse() * t, Scale3::identity()); + /// + /// // Work in all dimensions. + /// let t = Scale2::new(1.0, 2.0); + /// assert_eq!(t * t.pseudo_inverse(), Scale2::identity()); + /// assert_eq!(t.pseudo_inverse() * t, Scale2::identity()); + /// + /// // Inverts only non-zero coordinates. + /// let t = Scale2::new(0.0, 2.0); + /// assert_eq!(t * t.pseudo_inverse(), Scale2::new(0.0, 1.0)); + /// assert_eq!(t.pseudo_inverse() * t, Scale2::new(0.0, 1.0)); + /// ``` + #[inline] + #[must_use] + pub fn pseudo_inverse(&self) -> Scale + where + T: ClosedDiv + One + Zero, + { + return self + .vector + .map(|e| { + if e != T::zero() { + T::one() / e + } else { + T::zero() + } + }) + .into(); + } + + /// Converts this Scale into its equivalent homogeneous transformation matrix. + /// + /// # Example + /// ``` + /// # use nalgebra::{Scale2, Scale3, Matrix3, Matrix4}; + /// let t = Scale3::new(10.0, 20.0, 30.0); + /// let expected = Matrix4::new(10.0, 0.0, 0.0, 0.0, + /// 0.0, 20.0, 0.0, 0.0, + /// 0.0, 0.0, 30.0, 0.0, + /// 0.0, 0.0, 0.0, 1.0); + /// assert_eq!(t.to_homogeneous(), expected); + /// + /// let t = Scale2::new(10.0, 20.0); + /// let expected = Matrix3::new(10.0, 0.0, 0.0, + /// 0.0, 20.0, 0.0, + /// 0.0, 0.0, 1.0); + /// assert_eq!(t.to_homogeneous(), expected); + /// ``` + #[inline] + #[must_use] + pub fn to_homogeneous(&self) -> OMatrix, U1>, DimNameSum, U1>> + where + T: Zero + One + Clone, + Const: DimNameAdd, + DefaultAllocator: Allocator, U1>, DimNameSum, U1>> + + Allocator, U1>, U1>, + { + // TODO: use self.vector.push() instead. We can’t right now because + // that would require the DimAdd bound (but here we use DimNameAdd). + // This should be fixable once Rust gets a more complete support of + // const-generics. + let mut v = OVector::from_element(T::one()); + for i in 0..D { + v[i] = self.vector[i].clone(); + } + return OMatrix::from_diagonal(&v); + } + + /// Inverts `self` in-place. + /// + /// # Example + /// ``` + /// # use nalgebra::{Scale2, Scale3}; + /// let t = Scale3::new(1.0, 2.0, 3.0); + /// let mut inv_t = Scale3::new(1.0, 2.0, 3.0); + /// assert!(inv_t.try_inverse_mut()); + /// assert_eq!(t * inv_t, Scale3::identity()); + /// assert_eq!(inv_t * t, Scale3::identity()); + /// + /// // Work in all dimensions. + /// let t = Scale2::new(1.0, 2.0); + /// let mut inv_t = Scale2::new(1.0, 2.0); + /// assert!(inv_t.try_inverse_mut()); + /// assert_eq!(t * inv_t, Scale2::identity()); + /// assert_eq!(inv_t * t, Scale2::identity()); + /// + /// // Does not perform any operation if a coordinate is 0. + /// let mut t = Scale2::new(0.0, 2.0); + /// assert!(!t.try_inverse_mut()); + /// ``` + #[inline] + pub fn try_inverse_mut(&mut self) -> bool + where + T: ClosedDiv + One + Zero, + { + if let Some(v) = self.try_inverse() { + self.vector = v.vector; + true + } else { + false + } + } +} + +impl Scale { + /// Translate the given point. + /// + /// This is the same as the multiplication `self * pt`. + /// + /// # Example + /// ``` + /// # use nalgebra::{Scale3, Point3}; + /// let t = Scale3::new(1.0, 2.0, 3.0); + /// let transformed_point = t.transform_point(&Point3::new(4.0, 5.0, 6.0)); + /// assert_eq!(transformed_point, Point3::new(4.0, 10.0, 18.0)); + /// ``` + #[inline] + #[must_use] + pub fn transform_point(&self, pt: &Point) -> Point { + self * pt + } +} + +impl Scale { + /// Translate the given point by the inverse of this Scale. + /// + /// # Example + /// ``` + /// # use nalgebra::{Scale3, Point3}; + /// let t = Scale3::new(1.0, 2.0, 3.0); + /// let transformed_point = t.try_inverse_transform_point(&Point3::new(4.0, 6.0, 6.0)).unwrap(); + /// assert_eq!(transformed_point, Point3::new(4.0, 3.0, 2.0)); + /// + /// // Returns None if the inverse doesn't exist. + /// let t = Scale3::new(1.0, 0.0, 3.0); + /// let transformed_point = t.try_inverse_transform_point(&Point3::new(4.0, 6.0, 6.0)); + /// assert_eq!(transformed_point, None); + /// ``` + #[inline] + #[must_use] + pub fn try_inverse_transform_point(&self, pt: &Point) -> Option> { + self.try_inverse().map(|s| s * pt) + } +} + +impl Eq for Scale {} + +impl PartialEq for Scale { + #[inline] + fn eq(&self, right: &Scale) -> bool { + self.vector == right.vector + } +} + +impl AbsDiffEq for Scale +where + T::Epsilon: Clone, +{ + type Epsilon = T::Epsilon; + + #[inline] + fn default_epsilon() -> Self::Epsilon { + T::default_epsilon() + } + + #[inline] + fn abs_diff_eq(&self, other: &Self, epsilon: Self::Epsilon) -> bool { + self.vector.abs_diff_eq(&other.vector, epsilon) + } +} + +impl RelativeEq for Scale +where + T::Epsilon: Clone, +{ + #[inline] + fn default_max_relative() -> Self::Epsilon { + T::default_max_relative() + } + + #[inline] + fn relative_eq( + &self, + other: &Self, + epsilon: Self::Epsilon, + max_relative: Self::Epsilon, + ) -> bool { + self.vector + .relative_eq(&other.vector, epsilon, max_relative) + } +} + +impl UlpsEq for Scale +where + T::Epsilon: Clone, +{ + #[inline] + fn default_max_ulps() -> u32 { + T::default_max_ulps() + } + + #[inline] + fn ulps_eq(&self, other: &Self, epsilon: Self::Epsilon, max_ulps: u32) -> bool { + self.vector.ulps_eq(&other.vector, epsilon, max_ulps) + } +} + +/* + * + * Display + * + */ +impl fmt::Display for Scale { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + let precision = f.precision().unwrap_or(3); + + writeln!(f, "Scale {{")?; + write!(f, "{:.*}", precision, self.vector)?; + writeln!(f, "}}") + } +} diff --git a/src/geometry/scale_alias.rs b/src/geometry/scale_alias.rs new file mode 100644 index 00000000..87932442 --- /dev/null +++ b/src/geometry/scale_alias.rs @@ -0,0 +1,19 @@ +use crate::geometry::Scale; + +/// A 1-dimensional scale. +pub type Scale1 = Scale; + +/// A 2-dimensional scale. +pub type Scale2 = Scale; + +/// A 3-dimensional scale. +pub type Scale3 = Scale; + +/// A 4-dimensional scale. +pub type Scale4 = Scale; + +/// A 5-dimensional scale. +pub type Scale5 = Scale; + +/// A 6-dimensional scale. +pub type Scale6 = Scale; diff --git a/src/geometry/scale_construction.rs b/src/geometry/scale_construction.rs new file mode 100644 index 00000000..02cce69c --- /dev/null +++ b/src/geometry/scale_construction.rs @@ -0,0 +1,123 @@ +#[cfg(feature = "arbitrary")] +use crate::base::storage::Owned; +#[cfg(feature = "arbitrary")] +use quickcheck::{Arbitrary, Gen}; + +use num::One; +#[cfg(feature = "rand-no-std")] +use rand::{ + distributions::{Distribution, Standard}, + Rng, +}; + +use simba::scalar::{ClosedMul, SupersetOf}; + +use crate::base::{SVector, Scalar}; +use crate::geometry::Scale; + +impl Scale { + /// Creates a new identity scale. + /// + /// # Example + /// ``` + /// # use nalgebra::{Point2, Point3, Scale2, Scale3}; + /// let t = Scale2::identity(); + /// let p = Point2::new(1.0, 2.0); + /// assert_eq!(t * p, p); + /// + /// // Works in all dimensions. + /// let t = Scale3::identity(); + /// let p = Point3::new(1.0, 2.0, 3.0); + /// assert_eq!(t * p, p); + /// ``` + #[inline] + pub fn identity() -> Scale + where + T: One, + { + Scale::from(SVector::from_element(T::one())) + } + + /// Cast the components of `self` to another type. + /// + /// # Example + /// ``` + /// # use nalgebra::Scale2; + /// let tra = Scale2::new(1.0f64, 2.0); + /// let tra2 = tra.cast::(); + /// assert_eq!(tra2, Scale2::new(1.0f32, 2.0)); + /// ``` + pub fn cast(self) -> Scale + where + Scale: SupersetOf, + { + crate::convert(self) + } +} + +impl One for Scale { + #[inline] + fn one() -> Self { + Self::identity() + } +} + +#[cfg(feature = "rand-no-std")] +impl Distribution> for Standard +where + Standard: Distribution, +{ + /// Generate an arbitrary random variate for testing purposes. + #[inline] + fn sample(&self, rng: &mut G) -> Scale { + Scale::from(rng.gen::>()) + } +} + +#[cfg(feature = "arbitrary")] +impl Arbitrary for Scale +where + Owned>: Send, +{ + #[inline] + fn arbitrary(rng: &mut Gen) -> Self { + let v: SVector = Arbitrary::arbitrary(rng); + Self::from(v) + } +} + +/* + * + * Small Scale construction from components. + * + */ +macro_rules! componentwise_constructors_impl( + ($($doc: expr; $D: expr, $($args: ident:$irow: expr),*);* $(;)*) => {$( + impl Scale + { + #[doc = "Initializes this Scale from its components."] + #[doc = "# Example\n```"] + #[doc = $doc] + #[doc = "```"] + #[inline] + pub const fn new($($args: T),*) -> Self { + Self { vector: SVector::::new($($args),*) } + } + } + )*} +); + +componentwise_constructors_impl!( + "# use nalgebra::Scale1;\nlet t = Scale1::new(1.0);\nassert!(t.vector.x == 1.0);"; + 1, x:0; + "# use nalgebra::Scale2;\nlet t = Scale2::new(1.0, 2.0);\nassert!(t.vector.x == 1.0 && t.vector.y == 2.0);"; + 2, x:0, y:1; + "# use nalgebra::Scale3;\nlet t = Scale3::new(1.0, 2.0, 3.0);\nassert!(t.vector.x == 1.0 && t.vector.y == 2.0 && t.vector.z == 3.0);"; + 3, x:0, y:1, z:2; + "# use nalgebra::Scale4;\nlet t = Scale4::new(1.0, 2.0, 3.0, 4.0);\nassert!(t.vector.x == 1.0 && t.vector.y == 2.0 && t.vector.z == 3.0 && t.vector.w == 4.0);"; + 4, x:0, y:1, z:2, w:3; + "# use nalgebra::Scale5;\nlet t = Scale5::new(1.0, 2.0, 3.0, 4.0, 5.0);\nassert!(t.vector.x == 1.0 && t.vector.y == 2.0 && t.vector.z == 3.0 && t.vector.w == 4.0 && t.vector.a == 5.0);"; + 5, x:0, y:1, z:2, w:3, a:4; + "# use nalgebra::Scale6;\nlet t = Scale6::new(1.0, 2.0, 3.0, 4.0, 5.0, 6.0);\nassert!(t.vector.x == 1.0 && t.vector.y == 2.0 && t.vector.z == 3.0 && t.vector.w == 4.0 && t.vector.a == 5.0 && t.vector.b == 6.0);"; + 6, x:0, y:1, z:2, w:3, a:4, b:5; +); diff --git a/src/geometry/scale_conversion.rs b/src/geometry/scale_conversion.rs new file mode 100644 index 00000000..2dc670a1 --- /dev/null +++ b/src/geometry/scale_conversion.rs @@ -0,0 +1,233 @@ +use num::{One, Zero}; + +use simba::scalar::{RealField, SubsetOf, SupersetOf}; +use simba::simd::PrimitiveSimdValue; + +use crate::base::allocator::Allocator; +use crate::base::dimension::{DimNameAdd, DimNameSum, U1}; +use crate::base::{Const, DefaultAllocator, OMatrix, OVector, SVector, Scalar}; + +use crate::geometry::{Scale, SuperTCategoryOf, TAffine, Transform}; +use crate::Point; + +/* + * This file provides the following conversions: + * ============================================= + * + * Scale -> Scale + * Scale -> Transform + * Scale -> Matrix (homogeneous) + */ + +impl SubsetOf> for Scale +where + T1: Scalar, + T2: Scalar + SupersetOf, +{ + #[inline] + fn to_superset(&self) -> Scale { + Scale::from(self.vector.to_superset()) + } + + #[inline] + fn is_in_subset(rot: &Scale) -> bool { + crate::is_convertible::<_, SVector>(&rot.vector) + } + + #[inline] + fn from_superset_unchecked(rot: &Scale) -> Self { + Scale { + vector: rot.vector.to_subset_unchecked(), + } + } +} + +impl SubsetOf> for Scale +where + T1: RealField, + T2: RealField + SupersetOf, + C: SuperTCategoryOf, + Const: DimNameAdd, + DefaultAllocator: Allocator, U1>, DimNameSum, U1>> + + Allocator, U1>, U1> + + Allocator, U1>, DimNameSum, U1>>, +{ + #[inline] + fn to_superset(&self) -> Transform { + Transform::from_matrix_unchecked(self.to_homogeneous().to_superset()) + } + + #[inline] + fn is_in_subset(t: &Transform) -> bool { + >::is_in_subset(t.matrix()) + } + + #[inline] + fn from_superset_unchecked(t: &Transform) -> Self { + Self::from_superset_unchecked(t.matrix()) + } +} + +impl + SubsetOf, U1>, DimNameSum, U1>>> for Scale +where + T1: RealField, + T2: RealField + SupersetOf, + Const: DimNameAdd, + DefaultAllocator: Allocator, U1>, DimNameSum, U1>> + + Allocator, U1>, U1> + + Allocator, U1>, DimNameSum, U1>>, +{ + #[inline] + fn to_superset(&self) -> OMatrix, U1>, DimNameSum, U1>> { + self.to_homogeneous().to_superset() + } + + #[inline] + fn is_in_subset(m: &OMatrix, U1>, DimNameSum, U1>>) -> bool { + if m[(D, D)] != T2::one() { + return false; + } + for i in 0..D + 1 { + for j in 0..D + 1 { + if i != j && m[(i, j)] != T2::zero() { + return false; + } + } + } + true + } + + #[inline] + fn from_superset_unchecked( + m: &OMatrix, U1>, DimNameSum, U1>>, + ) -> Self { + let v = m.fixed_slice::(0, 0).diagonal(); + Self { + vector: crate::convert_unchecked(v), + } + } +} + +impl From> + for OMatrix, U1>, DimNameSum, U1>> +where + Const: DimNameAdd, + DefaultAllocator: Allocator, U1>, DimNameSum, U1>> + + Allocator, U1>, U1> + + Allocator>, +{ + #[inline] + fn from(t: Scale) -> Self { + t.to_homogeneous() + } +} + +impl From>> for Scale { + #[inline] + fn from(vector: OVector>) -> Self { + Scale { vector } + } +} + +impl From<[T; D]> for Scale { + #[inline] + fn from(coords: [T; D]) -> Self { + Scale { + vector: coords.into(), + } + } +} + +impl From> for Scale { + #[inline] + fn from(pt: Point) -> Self { + Scale { vector: pt.coords } + } +} + +impl From> for [T; D] { + #[inline] + fn from(t: Scale) -> Self { + t.vector.into() + } +} + +impl From<[Scale; 2]> for Scale +where + T: From<[::Element; 2]>, + T::Element: Scalar, +{ + #[inline] + fn from(arr: [Scale; 2]) -> Self { + Self::from(OVector::from([ + arr[0].vector.clone(), + arr[1].vector.clone(), + ])) + } +} + +impl From<[Scale; 4]> for Scale +where + T: From<[::Element; 4]>, + T::Element: Scalar, +{ + #[inline] + fn from(arr: [Scale; 4]) -> Self { + Self::from(OVector::from([ + arr[0].vector.clone(), + arr[1].vector.clone(), + arr[2].vector.clone(), + arr[3].vector.clone(), + ])) + } +} + +impl From<[Scale; 8]> for Scale +where + T: From<[::Element; 8]>, + T::Element: Scalar, +{ + #[inline] + fn from(arr: [Scale; 8]) -> Self { + Self::from(OVector::from([ + arr[0].vector.clone(), + arr[1].vector.clone(), + arr[2].vector.clone(), + arr[3].vector.clone(), + arr[4].vector.clone(), + arr[5].vector.clone(), + arr[6].vector.clone(), + arr[7].vector.clone(), + ])) + } +} + +impl From<[Scale; 16]> + for Scale +where + T: From<[::Element; 16]>, + T::Element: Scalar, +{ + #[inline] + fn from(arr: [Scale; 16]) -> Self { + Self::from(OVector::from([ + arr[0].vector.clone(), + arr[1].vector.clone(), + arr[2].vector.clone(), + arr[3].vector.clone(), + arr[4].vector.clone(), + arr[5].vector.clone(), + arr[6].vector.clone(), + arr[7].vector.clone(), + arr[8].vector.clone(), + arr[9].vector.clone(), + arr[10].vector.clone(), + arr[11].vector.clone(), + arr[12].vector.clone(), + arr[13].vector.clone(), + arr[14].vector.clone(), + arr[15].vector.clone(), + ])) + } +} diff --git a/src/geometry/scale_coordinates.rs b/src/geometry/scale_coordinates.rs new file mode 100644 index 00000000..5158c62d --- /dev/null +++ b/src/geometry/scale_coordinates.rs @@ -0,0 +1,39 @@ +use std::ops::{Deref, DerefMut}; + +use crate::base::coordinates::{X, XY, XYZ, XYZW, XYZWA, XYZWAB}; +use crate::base::Scalar; + +use crate::geometry::Scale; + +/* + * + * Give coordinates to Scale{1 .. 6} + * + */ + +macro_rules! deref_impl( + ($D: expr, $Target: ident $(, $comps: ident)*) => { + impl Deref for Scale { + type Target = $Target; + + #[inline] + fn deref(&self) -> &Self::Target { + self.vector.deref() + } + } + + impl DerefMut for Scale { + #[inline] + fn deref_mut(&mut self) -> &mut Self::Target { + self.vector.deref_mut() + } + } + } +); + +deref_impl!(1, X, x); +deref_impl!(2, XY, x, y); +deref_impl!(3, XYZ, x, y, z); +deref_impl!(4, XYZW, x, y, z, w); +deref_impl!(5, XYZWA, x, y, z, w, a); +deref_impl!(6, XYZWAB, x, y, z, w, a, b); diff --git a/src/geometry/scale_ops.rs b/src/geometry/scale_ops.rs new file mode 100644 index 00000000..c056a301 --- /dev/null +++ b/src/geometry/scale_ops.rs @@ -0,0 +1,125 @@ +use std::ops::{Mul, MulAssign}; + +use simba::scalar::ClosedMul; + +use crate::base::constraint::{SameNumberOfColumns, SameNumberOfRows, ShapeConstraint}; +use crate::base::dimension::U1; +use crate::base::{Const, SVector, Scalar}; + +use crate::geometry::{Point, Scale}; + +// Scale × Scale +add_sub_impl!(Mul, mul, ClosedMul; + (Const, U1), (Const, U1) -> (Const, U1) + const D; for; where; + self: &'a Scale, right: &'b Scale, Output = Scale; + Scale::from(self.vector.component_mul(&right.vector)); + 'a, 'b); + +add_sub_impl!(Mul, mul, ClosedMul; + (Const, U1), (Const, U1) -> (Const, U1) + const D; for; where; + self: &'a Scale, right: Scale, Output = Scale; + Scale::from(self.vector.component_mul(&right.vector)); + 'a); + +add_sub_impl!(Mul, mul, ClosedMul; + (Const, U1), (Const, U1) -> (Const, U1) + const D; for; where; + self: Scale, right: &'b Scale, Output = Scale; + Scale::from(self.vector.component_mul(&right.vector)); + 'b); + +add_sub_impl!(Mul, mul, ClosedMul; + (Const, U1), (Const, U1) -> (Const, U1) + const D; for; where; + self: Scale, right: Scale, Output = Scale; + Scale::from(self.vector.component_mul(&right.vector)); ); + +// Scale × scalar +add_sub_impl!(Mul, mul, ClosedMul; + (Const, U1), (Const, U1) -> (Const, U1) + const D; for; where; + self: &'a Scale, right: T, Output = Scale; + Scale::from(&self.vector * right); + 'a); + +add_sub_impl!(Mul, mul, ClosedMul; + (Const, U1), (Const, U1) -> (Const, U1) + const D; for; where; + self: Scale, right: T, Output = Scale; + Scale::from(self.vector * right); ); + +// Scale × Point +add_sub_impl!(Mul, mul, ClosedMul; + (Const, U1), (Const, U1) -> (Const, U1) + const D; for; where; + self: &'a Scale, right: &'b Point, Output = Point; + Point::from(self.vector.component_mul(&right.coords)); + 'a, 'b); + +add_sub_impl!(Mul, mul, ClosedMul; + (Const, U1), (Const, U1) -> (Const, U1) + const D; for; where; + self: &'a Scale, right: Point, Output = Point; + Point::from(self.vector.component_mul(&right.coords)); + 'a); + +add_sub_impl!(Mul, mul, ClosedMul; + (Const, U1), (Const, U1) -> (Const, U1) + const D; for; where; + self: Scale, right: &'b Point, Output = Point; + Point::from(self.vector.component_mul(&right.coords)); + 'b); + +add_sub_impl!(Mul, mul, ClosedMul; + (Const, U1), (Const, U1) -> (Const, U1) + const D; for; where; + self: Scale, right: Point, Output = Point; + Point::from(self.vector.component_mul(&right.coords)); ); + +// Scale * Vector +add_sub_impl!(Mul, mul, ClosedMul; + (Const, U1), (Const, U1) -> (Const, U1) + const D; for; where; + self: &'a Scale, right: &'b SVector, Output = SVector; + SVector::from(self.vector.component_mul(&right)); + 'a, 'b); + +add_sub_impl!(Mul, mul, ClosedMul; + (Const, U1), (Const, U1) -> (Const, U1) + const D; for; where; + self: &'a Scale, right: SVector, Output = SVector; + SVector::from(self.vector.component_mul(&right)); + 'a); + +add_sub_impl!(Mul, mul, ClosedMul; + (Const, U1), (Const, U1) -> (Const, U1) + const D; for; where; + self: Scale, right: &'b SVector, Output = SVector; + SVector::from(self.vector.component_mul(&right)); + 'b); + +add_sub_impl!(Mul, mul, ClosedMul; + (Const, U1), (Const, U1) -> (Const, U1) + const D; for; where; + self: Scale, right: SVector, Output = SVector; + SVector::from(self.vector.component_mul(&right)); ); + +// Scale *= Scale +add_sub_assign_impl!(MulAssign, mul_assign, ClosedMul; + const D; + self: Scale, right: &'b Scale; + self.vector.component_mul_assign(&right.vector); + 'b); + +add_sub_assign_impl!(MulAssign, mul_assign, ClosedMul; + const D; + self: Scale, right: Scale; + self.vector.component_mul_assign(&right.vector); ); + +// Scale ×= scalar +add_sub_assign_impl!(MulAssign, mul_assign, ClosedMul; + const D; + self: Scale, right: T; + self.vector *= right; ); diff --git a/src/geometry/scale_simba.rs b/src/geometry/scale_simba.rs new file mode 100755 index 00000000..cb42b715 --- /dev/null +++ b/src/geometry/scale_simba.rs @@ -0,0 +1,49 @@ +use simba::simd::SimdValue; + +use crate::base::OVector; +use crate::Scalar; + +use crate::geometry::Scale; + +impl SimdValue for Scale +where + T::Element: Scalar, +{ + type Element = Scale; + type SimdBool = T::SimdBool; + + #[inline] + fn lanes() -> usize { + T::lanes() + } + + #[inline] + fn splat(val: Self::Element) -> Self { + OVector::splat(val.vector).into() + } + + #[inline] + fn extract(&self, i: usize) -> Self::Element { + self.vector.extract(i).into() + } + + #[inline] + unsafe fn extract_unchecked(&self, i: usize) -> Self::Element { + self.vector.extract_unchecked(i).into() + } + + #[inline] + fn replace(&mut self, i: usize, val: Self::Element) { + self.vector.replace(i, val.vector) + } + + #[inline] + unsafe fn replace_unchecked(&mut self, i: usize, val: Self::Element) { + self.vector.replace_unchecked(i, val.vector) + } + + #[inline] + fn select(self, cond: Self::SimdBool, other: Self) -> Self { + self.vector.select(cond, other.vector).into() + } +} diff --git a/src/geometry/similarity.rs b/src/geometry/similarity.rs index 4cff61ce..fd6b1d36 100755 --- a/src/geometry/similarity.rs +++ b/src/geometry/similarity.rs @@ -23,7 +23,11 @@ use crate::geometry::{AbstractRotation, Isometry, Point, Translation}; /// A similarity, i.e., an uniform scaling, followed by a rotation, followed by a translation. #[repr(C)] -#[derive(Debug)] +#[derive(Debug, Copy, Clone)] +#[cfg_attr( + all(not(target_os = "cuda"), feature = "cuda"), + derive(cust::DeviceCopy) +)] #[cfg_attr(feature = "serde-serialize-no-std", derive(Serialize, Deserialize))] #[cfg_attr( feature = "serde-serialize-no-std", @@ -73,22 +77,6 @@ where } } -impl + Copy, const D: usize> Copy - for Similarity -where - Owned>: Copy, -{ -} - -impl + Clone, const D: usize> Clone - for Similarity -{ - #[inline] - fn clone(&self) -> Self { - Similarity::from_isometry(self.isometry.clone(), self.scaling.clone()) - } -} - impl Similarity where R: AbstractRotation, @@ -415,7 +403,7 @@ where #[inline] fn ulps_eq(&self, other: &Self, epsilon: Self::Epsilon, max_ulps: u32) -> bool { self.isometry - .ulps_eq(&other.isometry, epsilon.clone(), max_ulps.clone()) + .ulps_eq(&other.isometry, epsilon.clone(), max_ulps) && self.scaling.ulps_eq(&other.scaling, epsilon, max_ulps) } } diff --git a/src/geometry/similarity_construction.rs b/src/geometry/similarity_construction.rs index feb5719b..8d1d38b8 100644 --- a/src/geometry/similarity_construction.rs +++ b/src/geometry/similarity_construction.rs @@ -20,6 +20,16 @@ use crate::{ Translation, UnitComplex, UnitQuaternion, }; +impl Default for Similarity +where + T::Element: SimdRealField, + R: AbstractRotation, +{ + fn default() -> Self { + Self::identity() + } +} + impl Similarity where T::Element: SimdRealField, diff --git a/src/geometry/transform.rs b/src/geometry/transform.rs index f9dbeb51..b0b5cced 100755 --- a/src/geometry/transform.rs +++ b/src/geometry/transform.rs @@ -1,6 +1,6 @@ use approx::{AbsDiffEq, RelativeEq, UlpsEq}; use std::any::Any; -use std::fmt::Debug; +use std::fmt::{self, Debug}; use std::hash; use std::marker::PhantomData; @@ -60,14 +60,26 @@ where /// Tag representing the most general (not necessarily inversible) `Transform` type. #[derive(Debug, Copy, Clone, Hash, PartialEq, Eq)] +#[cfg_attr( + all(not(target_os = "cuda"), feature = "cuda"), + derive(cust::DeviceCopy) +)] pub enum TGeneral {} /// Tag representing the most general inversible `Transform` type. #[derive(Debug, Copy, Clone, Hash, PartialEq, Eq)] +#[cfg_attr( + all(not(target_os = "cuda"), feature = "cuda"), + derive(cust::DeviceCopy) +)] pub enum TProjective {} /// Tag representing an affine `Transform`. Its bottom-row is equal to `(0, 0 ... 0, 1)`. #[derive(Debug, Copy, Clone, Hash, PartialEq, Eq)] +#[cfg_attr( + all(not(target_os = "cuda"), feature = "cuda"), + derive(cust::DeviceCopy) +)] pub enum TAffine {} impl TCategory for TGeneral { @@ -157,7 +169,6 @@ super_tcategory_impl!( /// It is stored as a matrix with dimensions `(D + 1, D + 1)`, e.g., it stores a 4x4 matrix for a /// 3D transformation. #[repr(C)] -#[derive(Debug)] pub struct Transform where Const: DimNameAdd, @@ -167,6 +178,16 @@ where _phantom: PhantomData, } +impl Debug for Transform +where + Const: DimNameAdd, + DefaultAllocator: Allocator, U1>, DimNameSum, U1>>, +{ + fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> Result<(), fmt::Error> { + self.matrix.fmt(formatter) + } +} + impl hash::Hash for Transform where Const: DimNameAdd, @@ -186,6 +207,16 @@ where { } +#[cfg(all(not(target_os = "cuda"), feature = "cuda"))] +unsafe impl + cust::memory::DeviceCopy for Transform +where + Const: DimNameAdd, + DefaultAllocator: Allocator, U1>, DimNameSum, U1>>, + Owned, U1>, DimNameSum, U1>>: cust::memory::DeviceCopy, +{ +} + impl Clone for Transform where Const: DimNameAdd, diff --git a/src/geometry/transform_construction.rs b/src/geometry/transform_construction.rs index e9601864..45d61fcf 100644 --- a/src/geometry/transform_construction.rs +++ b/src/geometry/transform_construction.rs @@ -8,6 +8,16 @@ use crate::base::{Const, DefaultAllocator, OMatrix}; use crate::geometry::{TCategory, Transform}; +impl Default for Transform +where + Const: DimNameAdd, + DefaultAllocator: Allocator, U1>, DimNameSum, U1>>, +{ + fn default() -> Self { + Self::identity() + } +} + impl Transform where Const: DimNameAdd, diff --git a/src/geometry/translation.rs b/src/geometry/translation.rs index 8a64b97a..1ce8cba3 100755 --- a/src/geometry/translation.rs +++ b/src/geometry/translation.rs @@ -22,13 +22,23 @@ use crate::geometry::Point; /// A translation. #[repr(C)] -#[derive(Debug)] +#[cfg_attr( + all(not(target_os = "cuda"), feature = "cuda"), + derive(cust::DeviceCopy) +)] +#[derive(Copy, Clone)] pub struct Translation { /// The translation coordinates, i.e., how much is added to a point's coordinates when it is /// translated. pub vector: SVector, } +impl fmt::Debug for Translation { + fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> Result<(), fmt::Error> { + self.vector.as_slice().fmt(formatter) + } +} + impl hash::Hash for Translation where Owned>: hash::Hash, @@ -38,18 +48,6 @@ where } } -impl Copy for Translation {} - -impl Clone for Translation -where - Owned>: Clone, -{ - #[inline] - fn clone(&self) -> Self { - Translation::from(self.vector.clone()) - } -} - #[cfg(feature = "bytemuck")] unsafe impl bytemuck::Zeroable for Translation where diff --git a/src/geometry/translation_construction.rs b/src/geometry/translation_construction.rs index 5371b648..128937bf 100644 --- a/src/geometry/translation_construction.rs +++ b/src/geometry/translation_construction.rs @@ -15,6 +15,12 @@ use simba::scalar::{ClosedAdd, SupersetOf}; use crate::base::{SVector, Scalar}; use crate::geometry::Translation; +impl Default for Translation { + fn default() -> Self { + Self::identity() + } +} + impl Translation { /// Creates a new identity translation. /// diff --git a/src/geometry/unit_complex.rs b/src/geometry/unit_complex.rs index 2c621674..48405dd4 100755 --- a/src/geometry/unit_complex.rs +++ b/src/geometry/unit_complex.rs @@ -31,6 +31,9 @@ use std::cmp::{Eq, PartialEq}; /// * [Conversion to a matrix `to_rotation_matrix`, `to_homogeneous`…](#conversion-to-a-matrix) pub type UnitComplex = Unit>; +#[cfg(all(not(target_os = "cuda"), feature = "cuda"))] +unsafe impl cust::memory::DeviceCopy for UnitComplex {} + impl PartialEq for UnitComplex { #[inline] fn eq(&self, rhs: &Self) -> bool { @@ -458,8 +461,7 @@ impl UlpsEq for UnitComplex { #[inline] fn ulps_eq(&self, other: &Self, epsilon: Self::Epsilon, max_ulps: u32) -> bool { - self.re - .ulps_eq(&other.re, epsilon.clone(), max_ulps.clone()) + self.re.ulps_eq(&other.re, epsilon.clone(), max_ulps) && self.im.ulps_eq(&other.im, epsilon, max_ulps) } } diff --git a/src/geometry/unit_complex_construction.rs b/src/geometry/unit_complex_construction.rs index 0bf0188c..ebf4e81d 100644 --- a/src/geometry/unit_complex_construction.rs +++ b/src/geometry/unit_complex_construction.rs @@ -17,6 +17,15 @@ use crate::geometry::{Rotation2, UnitComplex}; use simba::scalar::{RealField, SupersetOf}; use simba::simd::SimdRealField; +impl Default for UnitComplex +where + T::Element: SimdRealField, +{ + fn default() -> Self { + Self::identity() + } +} + /// # Identity impl UnitComplex where diff --git a/src/lib.rs b/src/lib.rs index 5ce5cb18..28701cfa 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -71,10 +71,11 @@ an optimized set of tools for computer graphics and physics. Those features incl * Insertion and removal of rows of columns of a matrix. */ -#![allow(unused_variables, unused_mut)] #![deny( missing_docs, nonstandard_style, + unused_variables, + unused_mut, unused_parens, unused_qualifications, unused_results, diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index 51da364f..f61a4e63 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -74,6 +74,14 @@ where Cholesky { chol: matrix } } + /// Uses the given matrix as-is without any checks or modifications as the + /// Cholesky decomposition. + /// + /// It is up to the user to ensure all invariants hold. + pub fn pack_dirty(matrix: OMatrix) -> Self { + Cholesky { chol: matrix } + } + /// Retrieves the lower-triangular factor of the Cholesky decomposition with its strictly /// upper-triangular part filled with zeros. pub fn unpack(mut self) -> OMatrix { @@ -163,7 +171,32 @@ where /// /// Returns `None` if the input matrix is not definite-positive. The input matrix is assumed /// to be symmetric and only the lower-triangular part is read. - pub fn new(mut matrix: OMatrix) -> Option { + pub fn new(matrix: OMatrix) -> Option { + Self::new_internal(matrix, None) + } + + /// Attempts to approximate the Cholesky decomposition of `matrix` by + /// replacing non-positive values on the diagonals during the decomposition + /// with the given `substitute`. + /// + /// [`try_sqrt`](ComplexField::try_sqrt) will be applied to the `substitute` + /// when it has to be used. + /// + /// If your input matrix results only in positive values on the diagonals + /// during the decomposition, `substitute` is unused and the result is just + /// the same as if you used [`new`](Cholesky::new). + /// + /// This method allows to compensate for matrices with very small or even + /// negative values due to numerical errors but necessarily results in only + /// an approximation: it is basically a hack. If you don't specifically need + /// Cholesky, it may be better to consider alternatives like the + /// [`LU`](crate::linalg::LU) decomposition/factorization. + pub fn new_with_substitute(matrix: OMatrix, substitute: T) -> Option { + Self::new_internal(matrix, Some(substitute)) + } + + /// Common implementation for `new` and `new_with_substitute`. + fn new_internal(mut matrix: OMatrix, substitute: Option) -> Option { assert!(matrix.is_square(), "The input matrix must be square."); let n = matrix.nrows(); @@ -179,17 +212,25 @@ where col_j.axpy(factor.conjugate(), &col_k, T::one()); } - let diag = unsafe { matrix.get_unchecked((j, j)).clone() }; - if !diag.is_zero() { - if let Some(denom) = diag.try_sqrt() { - unsafe { - *matrix.get_unchecked_mut((j, j)) = denom.clone(); - } - - let mut col = matrix.slice_range_mut(j + 1.., j); - col /= denom; - continue; + let sqrt_denom = |v: T| { + if v.is_zero() { + return None; } + v.try_sqrt() + }; + + let diag = unsafe { matrix.get_unchecked((j, j)).clone() }; + + if let Some(denom) = + sqrt_denom(diag).or_else(|| substitute.clone().and_then(sqrt_denom)) + { + unsafe { + *matrix.get_unchecked_mut((j, j)) = denom.clone(); + } + + let mut col = matrix.slice_range_mut(j + 1.., j); + col /= denom; + continue; } // The diagonal element is either zero or its square root could not diff --git a/src/linalg/decomposition.rs b/src/linalg/decomposition.rs index 5b56b65a..c72babf3 100644 --- a/src/linalg/decomposition.rs +++ b/src/linalg/decomposition.rs @@ -1,8 +1,8 @@ use crate::storage::Storage; use crate::{ Allocator, Bidiagonal, Cholesky, ColPivQR, ComplexField, DefaultAllocator, Dim, DimDiff, - DimMin, DimMinimum, DimSub, FullPivLU, Hessenberg, Matrix, RealField, Schur, SymmetricEigen, - SymmetricTridiagonal, LU, QR, SVD, U1, UDU, + DimMin, DimMinimum, DimSub, FullPivLU, Hessenberg, Matrix, OMatrix, RealField, Schur, + SymmetricEigen, SymmetricTridiagonal, LU, QR, SVD, U1, UDU, }; /// # Rectangular matrix decomposition @@ -17,6 +17,7 @@ use crate::{ /// | LU with partial pivoting | `P⁻¹ * L * U` | `L` is lower-triangular with a diagonal filled with `1` and `U` is upper-triangular. `P` is a permutation matrix. | /// | LU with full pivoting | `P⁻¹ * L * U * Q⁻¹` | `L` is lower-triangular with a diagonal filled with `1` and `U` is upper-triangular. `P` and `Q` are permutation matrices. | /// | SVD | `U * Σ * Vᵀ` | `U` and `V` are two orthogonal matrices and `Σ` is a diagonal matrix containing the singular values. | +/// | Polar (Left Polar) | `P' * U` | `U` is semi-unitary/unitary and `P'` is a positive semi-definite Hermitian Matrix impl> Matrix { /// Computes the bidiagonalization using householder reflections. pub fn bidiagonalize(self) -> Bidiagonal @@ -74,7 +75,31 @@ impl> Matrix { } /// Computes the Singular Value Decomposition using implicit shift. + /// The singular values are guaranteed to be sorted in descending order. + /// If this order is not required consider using `svd_unordered`. pub fn svd(self, compute_u: bool, compute_v: bool) -> SVD + where + R: DimMin, + DimMinimum: DimSub, // for Bidiagonal. + DefaultAllocator: Allocator + + Allocator + + Allocator + + Allocator, U1>> + + Allocator, C> + + Allocator> + + Allocator> + + Allocator> + + Allocator, U1>> + + Allocator<(usize, usize), DimMinimum> + + Allocator<(T::RealField, usize), DimMinimum>, + { + SVD::new(self.into_owned(), compute_u, compute_v) + } + + /// Computes the Singular Value Decomposition using implicit shift. + /// The singular values are not guaranteed to be sorted in any particular order. + /// If a descending order is required, consider using `svd` instead. + pub fn svd_unordered(self, compute_u: bool, compute_v: bool) -> SVD where R: DimMin, DimMinimum: DimSub, // for Bidiagonal. @@ -88,10 +113,12 @@ impl> Matrix { + Allocator> + Allocator, U1>>, { - SVD::new(self.into_owned(), compute_u, compute_v) + SVD::new_unordered(self.into_owned(), compute_u, compute_v) } /// Attempts to compute the Singular Value Decomposition of `matrix` using implicit shift. + /// The singular values are guaranteed to be sorted in descending order. + /// If this order is not required consider using `try_svd_unordered`. /// /// # Arguments /// @@ -119,10 +146,103 @@ impl> Matrix { + Allocator> + Allocator> + Allocator> - + Allocator, U1>>, + + Allocator, U1>> + + Allocator<(usize, usize), DimMinimum> + + Allocator<(T::RealField, usize), DimMinimum>, { SVD::try_new(self.into_owned(), compute_u, compute_v, eps, max_niter) } + + /// Attempts to compute the Singular Value Decomposition of `matrix` using implicit shift. + /// The singular values are not guaranteed to be sorted in any particular order. + /// If a descending order is required, consider using `try_svd` instead. + /// + /// # Arguments + /// + /// * `compute_u` − set this to `true` to enable the computation of left-singular vectors. + /// * `compute_v` − set this to `true` to enable the computation of right-singular vectors. + /// * `eps` − tolerance used to determine when a value converged to 0. + /// * `max_niter` − maximum total number of iterations performed by the algorithm. If this + /// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm + /// continues indefinitely until convergence. + pub fn try_svd_unordered( + self, + compute_u: bool, + compute_v: bool, + eps: T::RealField, + max_niter: usize, + ) -> Option> + where + R: DimMin, + DimMinimum: DimSub, // for Bidiagonal. + DefaultAllocator: Allocator + + Allocator + + Allocator + + Allocator, U1>> + + Allocator, C> + + Allocator> + + Allocator> + + Allocator> + + Allocator, U1>>, + { + SVD::try_new_unordered(self.into_owned(), compute_u, compute_v, eps, max_niter) + } + + /// Computes the Polar Decomposition of a `matrix` (indirectly uses SVD). + pub fn polar(self) -> (OMatrix, OMatrix) + where + R: DimMin, + DimMinimum: DimSub, // for Bidiagonal. + DefaultAllocator: Allocator + + Allocator, R> + + Allocator> + + Allocator + + Allocator, DimMinimum> + + Allocator + + Allocator + + Allocator, U1>> + + Allocator, C> + + Allocator> + + Allocator> + + Allocator> + + Allocator, U1>>, + { + SVD::new_unordered(self.into_owned(), true, true) + .to_polar() + .unwrap() + } + + /// Attempts to compute the Polar Decomposition of a `matrix` (indirectly uses SVD). + /// + /// # Arguments + /// + /// * `eps` − tolerance used to determine when a value converged to 0 when computing the SVD. + /// * `max_niter` − maximum total number of iterations performed by the SVD computation algorithm. + pub fn try_polar( + self, + eps: T::RealField, + max_niter: usize, + ) -> Option<(OMatrix, OMatrix)> + where + R: DimMin, + DimMinimum: DimSub, // for Bidiagonal. + DefaultAllocator: Allocator + + Allocator, R> + + Allocator> + + Allocator + + Allocator, DimMinimum> + + Allocator + + Allocator + + Allocator, U1>> + + Allocator, C> + + Allocator> + + Allocator> + + Allocator> + + Allocator, U1>>, + { + SVD::try_new_unordered(self.into_owned(), true, true, eps, max_niter) + .and_then(|svd| svd.to_polar()) + } } /// # Square matrix decomposition diff --git a/src/linalg/determinant.rs b/src/linalg/determinant.rs index 7b5d6b2c..fe1ae82f 100644 --- a/src/linalg/determinant.rs +++ b/src/linalg/determinant.rs @@ -49,7 +49,7 @@ impl, S: Storage> SquareMatri let m33 = self.get_unchecked((2, 2)).clone(); let minor_m12_m23 = m22.clone() * m33.clone() - m32.clone() * m23.clone(); - let minor_m11_m23 = m21.clone() * m33.clone() - m31.clone() * m23.clone(); + let minor_m11_m23 = m21.clone() * m33 - m31.clone() * m23; let minor_m11_m22 = m21 * m32 - m31 * m22; m11 * minor_m12_m23 - m12 * minor_m11_m23 + m13 * minor_m11_m22 diff --git a/src/linalg/exp.rs b/src/linalg/exp.rs index 835730da..d6f61956 100644 --- a/src/linalg/exp.rs +++ b/src/linalg/exp.rs @@ -11,6 +11,47 @@ use crate::{ use crate::num::Zero; +/// Precomputed factorials for integers in range `0..=34`. +/// Note: `35!` does not fit into 128 bits. +// TODO: find a better place for this array? +const FACTORIAL: [u128; 35] = [ + 1, + 1, + 2, + 6, + 24, + 120, + 720, + 5040, + 40320, + 362880, + 3628800, + 39916800, + 479001600, + 6227020800, + 87178291200, + 1307674368000, + 20922789888000, + 355687428096000, + 6402373705728000, + 121645100408832000, + 2432902008176640000, + 51090942171709440000, + 1124000727777607680000, + 25852016738884976640000, + 620448401733239439360000, + 15511210043330985984000000, + 403291461126605635584000000, + 10888869450418352160768000000, + 304888344611713860501504000000, + 8841761993739701954543616000000, + 265252859812191058636308480000000, + 8222838654177922817725562880000000, + 263130836933693530167218012160000000, + 8683317618811886495518194401280000000, + 295232799039604140847618609643520000000, +]; + // https://github.com/scipy/scipy/blob/c1372d8aa90a73d8a52f135529293ff4edb98fc8/scipy/sparse/linalg/matfuncs.py struct ExpmPadeHelper where @@ -321,8 +362,8 @@ where self.calc_a2(); self.calc_a4(); self.calc_a6(); - let mb2 = self.a2.as_ref().unwrap() * convert::(2.0_f64.powf(-2.0 * s.clone())); - let mb4 = self.a4.as_ref().unwrap() * convert::(2.0.powf(-4.0 * s.clone())); + let mb2 = self.a2.as_ref().unwrap() * convert::(2.0_f64.powf(-2.0 * s)); + let mb4 = self.a4.as_ref().unwrap() * convert::(2.0.powf(-4.0 * s)); let mb6 = self.a6.as_ref().unwrap() * convert::(2.0.powf(-6.0 * s)); let u2 = &mb6 * (&mb6 * b[13].clone() + &mb4 * b[11].clone() + &mb2 * b[9].clone()); @@ -342,15 +383,17 @@ where } } -fn factorial(n: u128) -> u128 { - if n == 1 { - return 1; +/// Compute `n!` +#[inline(always)] +fn factorial(n: usize) -> u128 { + match FACTORIAL.get(n) { + Some(f) => *f, + None => panic!("{}! is greater than u128::MAX", n), } - n * factorial(n - 1) } /// Compute the 1-norm of a non-negative integer power of a non-negative matrix. -fn onenorm_matrix_power_nonm(a: &OMatrix, p: u64) -> T +fn onenorm_matrix_power_nonm(a: &OMatrix, p: usize) -> T where T: RealField, D: Dim, @@ -367,7 +410,7 @@ where v.max() } -fn ell(a: &OMatrix, m: u64) -> u64 +fn ell(a: &OMatrix, m: usize) -> u64 where T: ComplexField, D: Dim, @@ -376,8 +419,6 @@ where + Allocator + Allocator, { - // 2m choose m = (2m)!/(m! * (2m-m)!) - let a_abs = a.map(|x| x.abs()); let a_abs_onenorm = onenorm_matrix_power_nonm(&a_abs, 2 * m + 1); @@ -386,9 +427,11 @@ where return 0; } - let choose_2m_m = - factorial(2 * m as u128) / (factorial(m as u128) * factorial(2 * m as u128 - m as u128)); - let abs_c_recip = choose_2m_m * factorial(2 * m as u128 + 1); + // 2m choose m = (2m)!/(m! * (2m-m)!) = (2m)!/((m!)^2) + let m_factorial = factorial(m); + let choose_2m_m = factorial(2 * m) / (m_factorial * m_factorial); + + let abs_c_recip = choose_2m_m * factorial(2 * m + 1); let alpha = a_abs_onenorm / one_norm(a); let alpha: f64 = try_convert(alpha).unwrap() / abs_c_recip as f64; @@ -510,6 +553,7 @@ where #[cfg(test)] mod tests { #[test] + #[allow(clippy::float_cmp)] fn one_norm() { use crate::Matrix3; let m = Matrix3::new(-3.0, 5.0, 7.0, 2.0, 6.0, 4.0, 0.0, 2.0, 8.0); diff --git a/src/linalg/givens.rs b/src/linalg/givens.rs index c719deb6..30e2bed3 100644 --- a/src/linalg/givens.rs +++ b/src/linalg/givens.rs @@ -47,7 +47,7 @@ impl GivensRotation { if denom > eps { let norm = sign0.scale(denom.clone()); let c = mod0 / denom; - let s = s.clone() / norm.clone(); + let s = s / norm.clone(); Some((Self { c, s }, norm)) } else { None diff --git a/src/linalg/lu.rs b/src/linalg/lu.rs index b0fa065d..01ae46f0 100644 --- a/src/linalg/lu.rs +++ b/src/linalg/lu.rs @@ -317,7 +317,7 @@ where pub fn is_invertible(&self) -> bool { assert!( self.lu.is_square(), - "QR: unable to test the invertibility of a non-square matrix." + "LU: unable to test the invertibility of a non-square matrix." ); for i in 0..self.lu.nrows() { diff --git a/src/linalg/mod.rs b/src/linalg/mod.rs index 0c272494..0209f9b1 100644 --- a/src/linalg/mod.rs +++ b/src/linalg/mod.rs @@ -24,6 +24,8 @@ mod qr; mod schur; mod solve; mod svd; +mod svd2; +mod svd3; mod symmetric_eigen; mod symmetric_tridiagonal; mod udu; diff --git a/src/linalg/pow.rs b/src/linalg/pow.rs index df513643..e8c174d3 100644 --- a/src/linalg/pow.rs +++ b/src/linalg/pow.rs @@ -1,83 +1,71 @@ //! This module provides the matrix exponential (pow) function to square matrices. -use std::ops::DivAssign; - use crate::{ allocator::Allocator, storage::{Storage, StorageMut}, - DefaultAllocator, DimMin, Matrix, OMatrix, + DefaultAllocator, DimMin, Matrix, OMatrix, Scalar, }; -use num::PrimInt; -use simba::scalar::ComplexField; +use num::{One, Zero}; +use simba::scalar::{ClosedAdd, ClosedMul}; -impl Matrix +impl Matrix where + T: Scalar + Zero + One + ClosedAdd + ClosedMul, D: DimMin, S: StorageMut, DefaultAllocator: Allocator + Allocator, { - /// Attempts to raise this matrix to an integral power `e` in-place. If this - /// matrix is non-invertible and `e` is negative, it leaves this matrix - /// untouched and returns `false`. Otherwise, it returns `true` and - /// overwrites this matrix with the result. - pub fn pow_mut(&mut self, mut e: I) -> bool { - let zero = I::zero(); - + /// Raises this matrix to an integral power `exp` in-place. + pub fn pow_mut(&mut self, mut exp: u32) { // A matrix raised to the zeroth power is just the identity. - if e == zero { + if exp == 0 { self.fill_with_identity(); - return true; - } + } else if exp > 1 { + // We use the buffer to hold the result of multiplier^2, thus avoiding + // extra allocations. + let mut x = self.clone_owned(); + let mut workspace = self.clone_owned(); - // If e is negative, we compute the inverse matrix, then raise it to the - // power of -e. - if e < zero && !self.try_inverse_mut() { - return false; - } - - let one = I::one(); - let two = I::from(2u8).unwrap(); - - // We use the buffer to hold the result of multiplier ^ 2, thus avoiding - // extra allocations. - let mut multiplier = self.clone_owned(); - let mut buf = self.clone_owned(); - - // Exponentiation by squares. - loop { - if e % two == one { - self.mul_to(&multiplier, &mut buf); - self.copy_from(&buf); + if exp % 2 == 0 { + self.fill_with_identity(); + } else { + // Avoid an useless multiplication by the identity + // if the exponent is odd. + exp -= 1; } - e /= two; - multiplier.mul_to(&multiplier, &mut buf); - multiplier.copy_from(&buf); + // Exponentiation by squares. + loop { + if exp % 2 == 1 { + self.mul_to(&x, &mut workspace); + self.copy_from(&workspace); + } - if e == zero { - return true; + exp /= 2; + + if exp == 0 { + break; + } + + x.mul_to(&x, &mut workspace); + x.copy_from(&workspace); } } } } -impl> Matrix +impl> Matrix where + T: Scalar + Zero + One + ClosedAdd + ClosedMul, D: DimMin, S: StorageMut, DefaultAllocator: Allocator + Allocator, { - /// Attempts to raise this matrix to an integral power `e`. If this matrix - /// is non-invertible and `e` is negative, it returns `None`. Otherwise, it - /// returns the result as a new matrix. Uses exponentiation by squares. + /// Raise this matrix to an integral power `exp`. #[must_use] - pub fn pow(&self, e: I) -> Option> { - let mut clone = self.clone_owned(); - - if clone.pow_mut(e) { - Some(clone) - } else { - None - } + pub fn pow(&self, exp: u32) -> OMatrix { + let mut result = self.clone_owned(); + result.pow_mut(exp); + result } } diff --git a/src/linalg/qr.rs b/src/linalg/qr.rs index 5839f270..1b06e34b 100644 --- a/src/linalg/qr.rs +++ b/src/linalg/qr.rs @@ -147,6 +147,11 @@ where &self.qr } + #[must_use] + pub(crate) fn diag_internal(&self) -> &OVector> { + &self.diag + } + /// Multiplies the provided matrix by the transpose of the `Q` matrix of this decomposition. pub fn q_tr_mul(&self, rhs: &mut Matrix) // TODO: do we need a static constraint on the number of rows of rhs? diff --git a/src/linalg/svd.rs b/src/linalg/svd.rs index 5f1b0112..3f945a65 100644 --- a/src/linalg/svd.rs +++ b/src/linalg/svd.rs @@ -1,5 +1,6 @@ #[cfg(feature = "serde-serialize-no-std")] use serde::{Deserialize, Serialize}; +use std::any::TypeId; use approx::AbsDiffEq; use num::{One, Zero}; @@ -9,6 +10,7 @@ use crate::base::{DefaultAllocator, Matrix, Matrix2x3, OMatrix, OVector, Vector2 use crate::constraint::{SameNumberOfRows, ShapeConstraint}; use crate::dimension::{Dim, DimDiff, DimMin, DimMinimum, DimSub, U1}; use crate::storage::Storage; +use crate::{Matrix2, Matrix3, RawStorage, U2, U3}; use simba::scalar::{ComplexField, RealField}; use crate::linalg::givens::GivensRotation; @@ -78,9 +80,21 @@ where + Allocator> + Allocator, U1>>, { + fn use_special_always_ordered_svd2() -> bool { + TypeId::of::>() == TypeId::of::>() + && TypeId::of::() == TypeId::of::>() + } + + fn use_special_always_ordered_svd3() -> bool { + TypeId::of::>() == TypeId::of::>() + && TypeId::of::() == TypeId::of::>() + } + /// Computes the Singular Value Decomposition of `matrix` using implicit shift. - pub fn new(matrix: OMatrix, compute_u: bool, compute_v: bool) -> Self { - Self::try_new( + /// The singular values are not guaranteed to be sorted in any particular order. + /// If a descending order is required, consider using `new` instead. + pub fn new_unordered(matrix: OMatrix, compute_u: bool, compute_v: bool) -> Self { + Self::try_new_unordered( matrix, compute_u, compute_v, @@ -91,6 +105,8 @@ where } /// Attempts to compute the Singular Value Decomposition of `matrix` using implicit shift. + /// The singular values are not guaranteed to be sorted in any particular order. + /// If a descending order is required, consider using `try_new` instead. /// /// # Arguments /// @@ -100,7 +116,7 @@ where /// * `max_niter` − maximum total number of iterations performed by the algorithm. If this /// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm /// continues indefinitely until convergence. - pub fn try_new( + pub fn try_new_unordered( mut matrix: OMatrix, compute_u: bool, compute_v: bool, @@ -113,6 +129,21 @@ where ); let (nrows, ncols) = matrix.shape_generic(); let min_nrows_ncols = nrows.min(ncols); + + if Self::use_special_always_ordered_svd2() { + // SAFETY: the reference transmutes are OK since we checked that the types match exactly. + let matrix: &Matrix2 = unsafe { std::mem::transmute(&matrix) }; + let result = super::svd2::svd_ordered2(matrix, compute_u, compute_v); + let typed_result: &Self = unsafe { std::mem::transmute(&result) }; + return Some(typed_result.clone()); + } else if Self::use_special_always_ordered_svd3() { + // SAFETY: the reference transmutes are OK since we checked that the types match exactly. + let matrix: &Matrix3 = unsafe { std::mem::transmute(&matrix) }; + let result = super::svd3::svd_ordered3(matrix, compute_u, compute_v, eps, max_niter); + let typed_result: &Self = unsafe { std::mem::transmute(&result) }; + return Some(typed_result.clone()); + } + let dim = min_nrows_ncols.value(); let m_amax = matrix.camax(); @@ -610,6 +641,144 @@ where } } } + + /// converts SVD results to Polar decomposition form of the original Matrix: `A = P' * U`. + /// + /// The polar decomposition used here is Left Polar Decomposition (or Reverse Polar Decomposition) + /// Returns None if the singular vectors of the SVD haven't been calculated + pub fn to_polar(&self) -> Option<(OMatrix, OMatrix)> + where + DefaultAllocator: Allocator //result + + Allocator, R> // adjoint + + Allocator> // mapped vals + + Allocator // result + + Allocator, DimMinimum>, // square matrix + { + match (&self.u, &self.v_t) { + (Some(u), Some(v_t)) => Some(( + u * OMatrix::from_diagonal(&self.singular_values.map(|e| T::from_real(e))) + * u.adjoint(), + u * v_t, + )), + _ => None, + } + } +} + +impl, C: Dim> SVD +where + DimMinimum: DimSub, // for Bidiagonal. + DefaultAllocator: Allocator + + Allocator + + Allocator + + Allocator, U1>> + + Allocator, C> + + Allocator> + + Allocator> + + Allocator> + + Allocator, U1>> + + Allocator<(usize, usize), DimMinimum> // for sorted singular values + + Allocator<(T::RealField, usize), DimMinimum>, // for sorted singular values +{ + /// Computes the Singular Value Decomposition of `matrix` using implicit shift. + /// The singular values are guaranteed to be sorted in descending order. + /// If this order is not required consider using `new_unordered`. + pub fn new(matrix: OMatrix, compute_u: bool, compute_v: bool) -> Self { + let mut svd = Self::new_unordered(matrix, compute_u, compute_v); + + if !Self::use_special_always_ordered_svd3() && !Self::use_special_always_ordered_svd2() { + svd.sort_by_singular_values(); + } + + svd + } + + /// Attempts to compute the Singular Value Decomposition of `matrix` using implicit shift. + /// The singular values are guaranteed to be sorted in descending order. + /// If this order is not required consider using `try_new_unordered`. + /// + /// # Arguments + /// + /// * `compute_u` − set this to `true` to enable the computation of left-singular vectors. + /// * `compute_v` − set this to `true` to enable the computation of right-singular vectors. + /// * `eps` − tolerance used to determine when a value converged to 0. + /// * `max_niter` − maximum total number of iterations performed by the algorithm. If this + /// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm + /// continues indefinitely until convergence. + pub fn try_new( + matrix: OMatrix, + compute_u: bool, + compute_v: bool, + eps: T::RealField, + max_niter: usize, + ) -> Option { + Self::try_new_unordered(matrix, compute_u, compute_v, eps, max_niter).map(|mut svd| { + if !Self::use_special_always_ordered_svd3() && !Self::use_special_always_ordered_svd2() + { + svd.sort_by_singular_values(); + } + + svd + }) + } + + /// Sort the estimated components of the SVD by its singular values in descending order. + /// Such an ordering is often implicitly required when the decompositions are used for estimation or fitting purposes. + /// Using this function is only required if `new_unordered` or `try_new_unorderd` were used and the specific sorting is required afterward. + pub fn sort_by_singular_values(&mut self) { + const VALUE_PROCESSED: usize = usize::MAX; + + // Collect the singular values with their original index, ... + let mut singular_values = self.singular_values.map_with_location(|r, _, e| (e, r)); + assert_ne!( + singular_values.data.shape().0.value(), + VALUE_PROCESSED, + "Too many singular values" + ); + + // ... sort the singular values, ... + singular_values + .as_mut_slice() + .sort_unstable_by(|(a, _), (b, _)| b.partial_cmp(a).expect("Singular value was NaN")); + + // ... and store them. + self.singular_values + .zip_apply(&singular_values, |value, (new_value, _)| { + value.clone_from(&new_value) + }); + + // Calculate required permutations given the sorted indices. + // We need to identify all circles to calculate the required swaps. + let mut permutations = + crate::PermutationSequence::identity_generic(singular_values.data.shape().0); + + for i in 0..singular_values.len() { + let mut index_1 = i; + let mut index_2 = singular_values[i].1; + + // Check whether the value was already visited ... + while index_2 != VALUE_PROCESSED // ... or a "double swap" must be avoided. + && singular_values[index_2].1 != VALUE_PROCESSED + { + // Add the permutation ... + permutations.append_permutation(index_1, index_2); + // ... and mark the value as visited. + singular_values[index_1].1 = VALUE_PROCESSED; + + index_1 = index_2; + index_2 = singular_values[index_1].1; + } + } + + // Permute the optional components + if let Some(u) = self.u.as_mut() { + permutations.permute_columns(u); + } + + if let Some(v_t) = self.v_t.as_mut() { + permutations.permute_rows(v_t); + } + } } impl, C: Dim, S: Storage> Matrix @@ -626,9 +795,11 @@ where + Allocator, U1>>, { /// Computes the singular values of this matrix. + /// The singular values are not guaranteed to be sorted in any particular order. + /// If a descending order is required, consider using `singular_values` instead. #[must_use] - pub fn singular_values(&self) -> OVector> { - SVD::new(self.clone_owned(), false, false).singular_values + pub fn singular_values_unordered(&self) -> OVector> { + SVD::new_unordered(self.clone_owned(), false, false).singular_values } /// Computes the rank of this matrix. @@ -636,7 +807,7 @@ where /// All singular values below `eps` are considered equal to 0. #[must_use] pub fn rank(&self, eps: T::RealField) -> usize { - let svd = SVD::new(self.clone_owned(), false, false); + let svd = SVD::new_unordered(self.clone_owned(), false, false); svd.rank(eps) } @@ -647,7 +818,31 @@ where where DefaultAllocator: Allocator, { - SVD::new(self.clone_owned(), true, true).pseudo_inverse(eps) + SVD::new_unordered(self.clone_owned(), true, true).pseudo_inverse(eps) + } +} + +impl, C: Dim, S: Storage> Matrix +where + DimMinimum: DimSub, + DefaultAllocator: Allocator + + Allocator + + Allocator + + Allocator, U1>> + + Allocator, C> + + Allocator> + + Allocator> + + Allocator> + + Allocator, U1>> + + Allocator<(usize, usize), DimMinimum> + + Allocator<(T::RealField, usize), DimMinimum>, +{ + /// Computes the singular values of this matrix. + /// The singular values are guaranteed to be sorted in descending order. + /// If this order is not required consider using `singular_values_unordered`. + #[must_use] + pub fn singular_values(&self) -> OVector> { + SVD::new(self.clone_owned(), false, false).singular_values } } diff --git a/src/linalg/svd2.rs b/src/linalg/svd2.rs new file mode 100644 index 00000000..3bbd9a1b --- /dev/null +++ b/src/linalg/svd2.rs @@ -0,0 +1,50 @@ +use crate::{Matrix2, RealField, Vector2, SVD, U2}; + +// Implementation of the 2D SVD from https://ieeexplore.ieee.org/document/486688 +// See also https://scicomp.stackexchange.com/questions/8899/robust-algorithm-for-2-times-2-svd +pub fn svd_ordered2( + m: &Matrix2, + compute_u: bool, + compute_v: bool, +) -> SVD { + let half: T = crate::convert(0.5); + let one: T = crate::convert(1.0); + + let e = (m.m11.clone() + m.m22.clone()) * half.clone(); + let f = (m.m11.clone() - m.m22.clone()) * half.clone(); + let g = (m.m21.clone() + m.m12.clone()) * half.clone(); + let h = (m.m21.clone() - m.m12.clone()) * half.clone(); + let q = (e.clone() * e.clone() + h.clone() * h.clone()).sqrt(); + let r = (f.clone() * f.clone() + g.clone() * g.clone()).sqrt(); + + // Note that the singular values are always sorted because sx >= sy + // because q >= 0 and r >= 0. + let sx = q.clone() + r.clone(); + let sy = q - r; + let sy_sign = if sy < T::zero() { -one.clone() } else { one }; + let singular_values = Vector2::new(sx, sy * sy_sign.clone()); + + if compute_u || compute_v { + let a1 = g.atan2(f); + let a2 = h.atan2(e); + let theta = (a2.clone() - a1.clone()) * half.clone(); + let phi = (a2 + a1) * half; + let (st, ct) = theta.sin_cos(); + let (sp, cp) = phi.sin_cos(); + + let u = Matrix2::new(cp.clone(), -sp.clone(), sp, cp); + let v_t = Matrix2::new(ct.clone(), -st.clone(), st * sy_sign.clone(), ct * sy_sign); + + SVD { + u: if compute_u { Some(u) } else { None }, + singular_values, + v_t: if compute_v { Some(v_t) } else { None }, + } + } else { + SVD { + u: None, + singular_values, + v_t: None, + } + } +} diff --git a/src/linalg/svd3.rs b/src/linalg/svd3.rs new file mode 100644 index 00000000..b36e0889 --- /dev/null +++ b/src/linalg/svd3.rs @@ -0,0 +1,55 @@ +use crate::{Matrix3, SVD, U3}; +use simba::scalar::RealField; + +// For the 3x3 case, on the GPU, it is much more efficient to compute the SVD +// using an eigendecomposition followed by a QR decomposition. +// +// This is based on the paper "Computing the Singular Value Decomposition of 3 x 3 matrices with +// minimal branching and elementary floating point operations" from McAdams, et al. +pub fn svd_ordered3( + m: &Matrix3, + compute_u: bool, + compute_v: bool, + eps: T, + niter: usize, +) -> Option> { + let s = m.tr_mul(&m); + let mut v = s.try_symmetric_eigen(eps, niter)?.eigenvectors; + let mut b = m * &v; + + // Sort singular values. This is a necessary step to ensure that + // the QR decompositions R matrix ends up diagonal. + let mut rho0 = b.column(0).norm_squared(); + let mut rho1 = b.column(1).norm_squared(); + let mut rho2 = b.column(2).norm_squared(); + + if rho0 < rho1 { + b.swap_columns(0, 1); + b.column_mut(1).neg_mut(); + v.swap_columns(0, 1); + v.column_mut(1).neg_mut(); + std::mem::swap(&mut rho0, &mut rho1); + } + if rho0 < rho2 { + b.swap_columns(0, 2); + b.column_mut(2).neg_mut(); + v.swap_columns(0, 2); + v.column_mut(2).neg_mut(); + std::mem::swap(&mut rho0, &mut rho2); + } + if rho1 < rho2 { + b.swap_columns(1, 2); + b.column_mut(2).neg_mut(); + v.swap_columns(1, 2); + v.column_mut(2).neg_mut(); + std::mem::swap(&mut rho0, &mut rho2); + } + + let qr = b.qr(); + + Some(SVD { + u: if compute_u { Some(qr.q()) } else { None }, + singular_values: qr.diag_internal().map(|e| e.abs()), + v_t: if compute_v { Some(v.transpose()) } else { None }, + }) +} diff --git a/src/third_party/glam/mod.rs b/src/third_party/glam/mod.rs index 9d458db6..b1006e91 100644 --- a/src/third_party/glam/mod.rs +++ b/src/third_party/glam/mod.rs @@ -8,3 +8,5 @@ mod v015; mod v016; #[cfg(feature = "glam017")] mod v017; +#[cfg(feature = "glam018")] +mod v018; diff --git a/src/third_party/glam/v018/mod.rs b/src/third_party/glam/v018/mod.rs new file mode 100644 index 00000000..7124fd82 --- /dev/null +++ b/src/third_party/glam/v018/mod.rs @@ -0,0 +1,18 @@ +#[path = "../common/glam_isometry.rs"] +mod glam_isometry; +#[path = "../common/glam_matrix.rs"] +mod glam_matrix; +#[path = "../common/glam_point.rs"] +mod glam_point; +#[path = "../common/glam_quaternion.rs"] +mod glam_quaternion; +#[path = "../common/glam_rotation.rs"] +mod glam_rotation; +#[path = "../common/glam_similarity.rs"] +mod glam_similarity; +#[path = "../common/glam_translation.rs"] +mod glam_translation; +#[path = "../common/glam_unit_complex.rs"] +mod glam_unit_complex; + +pub(self) use glam018 as glam; diff --git a/tests/core/matrix.rs b/tests/core/matrix.rs index 4a35fb20..8a545e97 100644 --- a/tests/core/matrix.rs +++ b/tests/core/matrix.rs @@ -80,8 +80,8 @@ fn iter() { #[test] fn debug_output_corresponds_to_data_container() { let m = Matrix2::new(1.0, 2.0, 3.0, 4.0); - let output_stable = "Matrix { data: [[1, 3], [2, 4]] }"; // Current output on the stable channel. - let output_nightly = "Matrix { data: [[1.0, 3.0], [2.0, 4.0]] }"; // Current output on the nightly channel. + let output_stable = "[[1, 3], [2, 4]]"; // Current output on the stable channel. + let output_nightly = "[[1.0, 3.0], [2.0, 4.0]]"; // Current output on the nightly channel. let current_output = format!("{:?}", m); dbg!(output_stable); dbg!(output_nightly); diff --git a/tests/geometry/dual_quaternion.rs b/tests/geometry/dual_quaternion.rs index 6cc975a5..1926fee9 100644 --- a/tests/geometry/dual_quaternion.rs +++ b/tests/geometry/dual_quaternion.rs @@ -1,7 +1,7 @@ #![cfg(feature = "proptest-support")] #![allow(non_snake_case)] -use na::{DualQuaternion, Point3, UnitDualQuaternion, Vector3}; +use na::{DualQuaternion, Point3, Unit, UnitDualQuaternion, UnitQuaternion, Vector3}; use crate::proptest::*; use proptest::{prop_assert, proptest}; @@ -74,6 +74,98 @@ proptest!( prop_assert!(relative_eq!((dq * t) * p, dq * (t * p), epsilon = 1.0e-7)); } + #[cfg_attr(rustfmt, rustfmt_skip)] + #[test] + fn sclerp_is_defined_for_identical_orientations( + dq in unit_dual_quaternion(), + s in -1.0f64..2.0f64, + t in translation3(), + ) { + // Should not panic. + prop_assert!(relative_eq!(dq.sclerp(&dq, 0.0), dq, epsilon = 1.0e-7)); + prop_assert!(relative_eq!(dq.sclerp(&dq, 0.5), dq, epsilon = 1.0e-7)); + prop_assert!(relative_eq!(dq.sclerp(&dq, 1.0), dq, epsilon = 1.0e-7)); + prop_assert!(relative_eq!(dq.sclerp(&dq, s), dq, epsilon = 1.0e-7)); + + let unit = UnitDualQuaternion::identity(); + prop_assert!(relative_eq!(unit.sclerp(&unit, 0.0), unit, epsilon = 1.0e-7)); + prop_assert!(relative_eq!(unit.sclerp(&unit, 0.5), unit, epsilon = 1.0e-7)); + prop_assert!(relative_eq!(unit.sclerp(&unit, 1.0), unit, epsilon = 1.0e-7)); + prop_assert!(relative_eq!(unit.sclerp(&unit, s), unit, epsilon = 1.0e-7)); + + // ScLERPing two unit dual quaternions with nearly equal rotation + // components should result in a unit dual quaternion with a rotation + // component nearly equal to either input. + let dq2 = t * dq; + prop_assert!(relative_eq!(dq.sclerp(&dq2, 0.0).real, dq.real, epsilon = 1.0e-7)); + prop_assert!(relative_eq!(dq.sclerp(&dq2, 0.5).real, dq.real, epsilon = 1.0e-7)); + prop_assert!(relative_eq!(dq.sclerp(&dq2, 1.0).real, dq.real, epsilon = 1.0e-7)); + prop_assert!(relative_eq!(dq.sclerp(&dq2, s).real, dq.real, epsilon = 1.0e-7)); + + // ScLERPing two unit dual quaternions with nearly equal rotation + // components should result in a unit dual quaternion with a translation + // component which is nearly equal to linearly interpolating the + // translation components of the inputs. + prop_assert!(relative_eq!( + dq.sclerp(&dq2, s).translation().vector, + dq.translation().vector.lerp(&dq2.translation().vector, s), + epsilon = 1.0e-7 + )); + + let unit2 = t * unit; + prop_assert!(relative_eq!(unit.sclerp(&unit2, 0.0).real, unit.real, epsilon = 1.0e-7)); + prop_assert!(relative_eq!(unit.sclerp(&unit2, 0.5).real, unit.real, epsilon = 1.0e-7)); + prop_assert!(relative_eq!(unit.sclerp(&unit2, 1.0).real, unit.real, epsilon = 1.0e-7)); + prop_assert!(relative_eq!(unit.sclerp(&unit2, s).real, unit.real, epsilon = 1.0e-7)); + + prop_assert!(relative_eq!( + unit.sclerp(&unit2, s).translation().vector, + unit.translation().vector.lerp(&unit2.translation().vector, s), + epsilon = 1.0e-7 + )); + } + + #[cfg_attr(rustfmt, rustfmt_skip)] + #[test] + fn sclerp_is_not_defined_for_opposite_orientations( + dq in unit_dual_quaternion(), + s in 0.1f64..0.9f64, + t in translation3(), + t2 in translation3(), + v in vector3(), + ) { + let iso = dq.to_isometry(); + let rot = iso.rotation; + if let Some((axis, angle)) = rot.axis_angle() { + let flipped = UnitQuaternion::from_axis_angle(&axis, angle + std::f64::consts::PI); + let dqf = flipped * rot.inverse() * dq.clone(); + prop_assert!(dq.try_sclerp(&dqf, 0.5, 1.0e-7).is_none()); + prop_assert!(dq.try_sclerp(&dqf, s, 1.0e-7).is_none()); + } + + let dq2 = t * dq; + let iso2 = dq2.to_isometry(); + let rot2 = iso2.rotation; + if let Some((axis, angle)) = rot2.axis_angle() { + let flipped = UnitQuaternion::from_axis_angle(&axis, angle + std::f64::consts::PI); + let dq3f = t2 * flipped * rot.inverse() * dq.clone(); + prop_assert!(dq2.try_sclerp(&dq3f, 0.5, 1.0e-7).is_none()); + prop_assert!(dq2.try_sclerp(&dq3f, s, 1.0e-7).is_none()); + } + + if let Some(axis) = Unit::try_new(v, 1.0e-7) { + let unit = UnitDualQuaternion::identity(); + let flip = UnitQuaternion::from_axis_angle(&axis, std::f64::consts::PI); + let unitf = flip * unit; + prop_assert!(unit.try_sclerp(&unitf, 0.5, 1.0e-7).is_none()); + prop_assert!(unit.try_sclerp(&unitf, s, 1.0e-7).is_none()); + + let unit2f = t * unit * flip; + prop_assert!(unit.try_sclerp(&unit2f, 0.5, 1.0e-7).is_none()); + prop_assert!(unit.try_sclerp(&unit2f, s, 1.0e-7).is_none()); + } + } + #[cfg_attr(rustfmt, rustfmt_skip)] #[test] fn all_op_exist( diff --git a/tests/lib.rs b/tests/lib.rs index 20d38d7a..11c1e23f 100644 --- a/tests/lib.rs +++ b/tests/lib.rs @@ -33,3 +33,10 @@ mod proptest; //#[cfg(all(feature = "debug", feature = "compare", feature = "rand"))] //#[cfg(feature = "sparse")] //mod sparse; + +mod utils { + /// Checks if a slice is sorted in descending order. + pub fn is_sorted_descending(slice: &[T]) -> bool { + slice.windows(2).all(|elts| elts[0] >= elts[1]) + } +} diff --git a/tests/linalg/cholesky.rs b/tests/linalg/cholesky.rs index 6fd83912..891e54ca 100644 --- a/tests/linalg/cholesky.rs +++ b/tests/linalg/cholesky.rs @@ -1,5 +1,16 @@ #![cfg(all(feature = "proptest-support", feature = "debug"))] +#[test] +// #[rustfmt::skip] +fn cholesky_with_substitute() { + // Make a tiny covariance matrix with a small covariance value. + let m = na::Matrix2::new(1.0, f64::NAN, 1.0, 1e-32); + // Show that the cholesky fails for our matrix. We then try again with a substitute. + assert!(na::Cholesky::new(m).is_none()); + // ...and show that we get some result this time around. + assert!(na::Cholesky::new_with_substitute(m, 1e-8).is_some()); +} + macro_rules! gen_tests( ($module: ident, $scalar: ty) => { mod $module { diff --git a/tests/linalg/mod.rs b/tests/linalg/mod.rs index 9c252bfd..d9bd6cd9 100644 --- a/tests/linalg/mod.rs +++ b/tests/linalg/mod.rs @@ -9,6 +9,7 @@ mod full_piv_lu; mod hessenberg; mod inverse; mod lu; +mod pow; mod qr; mod schur; mod solve; diff --git a/tests/linalg/pow.rs b/tests/linalg/pow.rs new file mode 100644 index 00000000..22a1a0c0 --- /dev/null +++ b/tests/linalg/pow.rs @@ -0,0 +1,49 @@ +#[cfg(feature = "proptest-support")] +mod proptest_tests { + macro_rules! gen_tests( + ($module: ident, $scalar: expr, $scalar_type: ty) => { + mod $module { + use na::DMatrix; + #[allow(unused_imports)] + use crate::core::helper::{RandScalar, RandComplex}; + use std::cmp; + + use crate::proptest::*; + use proptest::{prop_assert, proptest}; + + proptest! { + #[test] + fn pow(n in PROPTEST_MATRIX_DIM, p in 0u32..=4) { + let n = cmp::max(1, cmp::min(n, 10)); + let m = DMatrix::<$scalar_type>::new_random(n, n).map(|e| e.0); + let m_pow = m.pow(p); + let mut expected = m.clone(); + expected.fill_with_identity(); + + for _ in 0..p { + expected = &m * &expected; + } + + prop_assert!(relative_eq!(m_pow, expected, epsilon = 1.0e-5)) + } + + #[test] + fn pow_static_square_4x4(m in matrix4_($scalar), p in 0u32..=4) { + let mut expected = m.clone(); + let m_pow = m.pow(p); + expected.fill_with_identity(); + + for _ in 0..p { + expected = &m * &expected; + } + + prop_assert!(relative_eq!(m_pow, expected, epsilon = 1.0e-5)) + } + } + } + } + ); + + gen_tests!(complex, complex_f64(), RandComplex); + gen_tests!(f64, PROPTEST_F64, RandScalar); +} diff --git a/tests/linalg/svd.rs b/tests/linalg/svd.rs index 80aa6a20..deb3d38d 100644 --- a/tests/linalg/svd.rs +++ b/tests/linalg/svd.rs @@ -1,3 +1,4 @@ +use crate::utils::is_sorted_descending; use na::{DMatrix, Matrix6}; #[cfg(feature = "proptest-support")] @@ -14,6 +15,7 @@ mod proptest_tests { use crate::core::helper::{RandScalar, RandComplex}; use crate::proptest::*; use proptest::{prop_assert, proptest}; + use crate::utils::is_sorted_descending; proptest! { #[test] @@ -26,6 +28,7 @@ mod proptest_tests { prop_assert!(s.iter().all(|e| *e >= 0.0)); prop_assert!(relative_eq!(&u * ds * &v_t, recomp_m, epsilon = 1.0e-5)); prop_assert!(relative_eq!(m, recomp_m, epsilon = 1.0e-5)); + prop_assert!(is_sorted_descending(s.as_slice())); } #[test] @@ -38,6 +41,7 @@ mod proptest_tests { prop_assert!(relative_eq!(m, &u * ds * &v_t, epsilon = 1.0e-5)); prop_assert!(u.is_orthogonal(1.0e-5)); prop_assert!(v_t.is_orthogonal(1.0e-5)); + prop_assert!(is_sorted_descending(s.as_slice())); } #[test] @@ -50,6 +54,7 @@ mod proptest_tests { prop_assert!(relative_eq!(m, &u * ds * &v_t, epsilon = 1.0e-5)); prop_assert!(u.is_orthogonal(1.0e-5)); prop_assert!(v_t.is_orthogonal(1.0e-5)); + prop_assert!(is_sorted_descending(s.as_slice())); } #[test] @@ -61,6 +66,7 @@ mod proptest_tests { prop_assert!(s.iter().all(|e| *e >= 0.0)); prop_assert!(relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5)); + prop_assert!(is_sorted_descending(s.as_slice())); } #[test] @@ -71,6 +77,7 @@ mod proptest_tests { prop_assert!(s.iter().all(|e| *e >= 0.0)); prop_assert!(relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5)); + prop_assert!(is_sorted_descending(s.as_slice())); } #[test] @@ -83,6 +90,7 @@ mod proptest_tests { prop_assert!(relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5)); prop_assert!(u.is_orthogonal(1.0e-5)); prop_assert!(v_t.is_orthogonal(1.0e-5)); + prop_assert!(is_sorted_descending(s.as_slice())); } #[test] @@ -95,6 +103,20 @@ mod proptest_tests { prop_assert!(relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5)); prop_assert!(u.is_orthogonal(1.0e-5)); prop_assert!(v_t.is_orthogonal(1.0e-5)); + prop_assert!(is_sorted_descending(s.as_slice())); + } + + #[test] + fn svd_static_square_3x3(m in matrix3_($scalar)) { + let svd = m.svd(true, true); + let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap()); + let ds = Matrix3::from_diagonal(&s.map(|e| ComplexField::from_real(e))); + + prop_assert!(s.iter().all(|e| *e >= 0.0)); + prop_assert!(relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5)); + prop_assert!(u.is_orthogonal(1.0e-5)); + prop_assert!(v_t.is_orthogonal(1.0e-5)); + prop_assert!(is_sorted_descending(s.as_slice())); } #[test] @@ -131,6 +153,25 @@ mod proptest_tests { prop_assert!(relative_eq!(&m * &sol2, b2, epsilon = 1.0e-6)); } } + + #[test] + fn svd_polar_decomposition(m in dmatrix_($scalar)) { + let svd = m.clone().svd_unordered(true, true); + let (p, u) = svd.to_polar().unwrap(); + + assert_relative_eq!(m, &p* &u, epsilon = 1.0e-5); + // semi-unitary check + assert!(u.is_orthogonal(1.0e-5) || u.transpose().is_orthogonal(1.0e-5)); + // hermitian check + assert_relative_eq!(p, p.adjoint(), epsilon = 1.0e-5); + + /* + * Same thing, but using the method instead of calling the SVD explicitly. + */ + let (p2, u2) = m.clone().polar(); + assert_eq!(p, p2); + assert_eq!(u, u2); + } } } } @@ -175,6 +216,7 @@ fn svd_singular() { let ds = DMatrix::from_diagonal(&s); assert!(s.iter().all(|e| *e >= 0.0)); + assert!(is_sorted_descending(s.as_slice())); assert!(u.is_orthogonal(1.0e-5)); assert!(v_t.is_orthogonal(1.0e-5)); assert_relative_eq!(m, &u * ds * &v_t, epsilon = 1.0e-5); @@ -217,6 +259,7 @@ fn svd_singular_vertical() { let ds = DMatrix::from_diagonal(&s); assert!(s.iter().all(|e| *e >= 0.0)); + assert!(is_sorted_descending(s.as_slice())); assert_relative_eq!(m, &u * ds * &v_t, epsilon = 1.0e-5); } @@ -255,6 +298,7 @@ fn svd_singular_horizontal() { let ds = DMatrix::from_diagonal(&s); assert!(s.iter().all(|e| *e >= 0.0)); + assert!(is_sorted_descending(s.as_slice())); assert_relative_eq!(m, &u * ds * &v_t, epsilon = 1.0e-5); } @@ -326,11 +370,41 @@ fn svd_fail() { 0.07311092531259344, 0.5579247949052946, 0.14518764691585773, 0.03502980663114896, 0.7991329455957719, 0.4929930019965745, 0.12293810556077789, 0.6617084679545999, 0.9002240700227326, 0.027153062135304884, 0.3630189466989524, 0.18207502727558866, 0.843196731466686, 0.08951878746549924, 0.7533450877576973, 0.009558876499740077, 0.9429679490873482, 0.9355764454129878); + + // Check unordered ... + let svd = m.clone().svd_unordered(true, true); + let recomp = svd.recompose().unwrap(); + assert_relative_eq!(m, recomp, epsilon = 1.0e-5); + + // ... and ordered SVD. let svd = m.clone().svd(true, true); let recomp = svd.recompose().unwrap(); assert_relative_eq!(m, recomp, epsilon = 1.0e-5); } +#[test] +#[rustfmt::skip] +fn svd3_fail() { + // NOTE: this matrix fails the special case done for 3x3 SVDs. + // It was found on an actual application using SVD as part of the minimization of a + // quadratic error function. + let m = nalgebra::matrix![ + 0.0, 1.0, 0.0; + 0.0, 1.7320508075688772, 0.0; + 0.0, 0.0, 0.0 + ]; + + // Check unordered ... + let svd = m.svd_unordered(true, true); + let recomp = svd.recompose().unwrap(); + assert_relative_eq!(m, recomp, epsilon = 1.0e-5); + + // ... and ordered SVD. + let svd = m.svd(true, true); + let recomp = svd.recompose().unwrap(); + assert_relative_eq!(m, recomp, epsilon = 1.0e-5); +} + #[test] fn svd_err() { let m = DMatrix::from_element(10, 10, 0.0); @@ -344,3 +418,45 @@ fn svd_err() { svd.clone().pseudo_inverse(-1.0) ); } + +#[test] +#[rustfmt::skip] +fn svd_sorted() { + let reference = nalgebra::matrix![ + 1.0, 2.0, 3.0, 4.0; + 5.0, 6.0, 7.0, 8.0; + 9.0, 10.0, 11.0, 12.0 + ]; + + let mut svd = nalgebra::SVD { + singular_values: nalgebra::matrix![1.72261225; 2.54368356e+01; 5.14037515e-16], + u: Some(nalgebra::matrix![ + -0.88915331, -0.20673589, 0.40824829; + -0.25438183, -0.51828874, -0.81649658; + 0.38038964, -0.82984158, 0.40824829 + ]), + v_t: Some(nalgebra::matrix![ + 0.73286619, 0.28984978, -0.15316664, -0.59618305; + -0.40361757, -0.46474413, -0.52587069, -0.58699725; + 0.44527162, -0.83143156, 0.32704826, 0.05911168 + ]), + }; + + assert_relative_eq!( + svd.recompose().expect("valid SVD"), + reference, + epsilon = 1.0e-5 + ); + + svd.sort_by_singular_values(); + + // Ensure successful sorting + assert_relative_eq!(svd.singular_values.x, 2.54368356e+01, epsilon = 1.0e-5); + + // Ensure that the sorted components represent the same decomposition + assert_relative_eq!( + svd.recompose().expect("valid SVD"), + reference, + epsilon = 1.0e-5 + ); +}