diff --git a/CHANGELOG.md b/CHANGELOG.md index bf6b6a0f..c00c01fc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,7 +4,29 @@ documented here. This project adheres to [Semantic Versioning](https://semver.org/). +## [0.31.0] (30 Apr. 2022) + +### Breaking changes +- Switch to `cust` 0.3 (for CUDA support). +- Switch to `rkyv` 0.7 +- Remove support for serialization based on `abomonation`. +- Remove support for conversions between `nalgebra` types and `glam` 0.13. + +### Modified +- The aliases for `Const` types have been simplified to help `rust-analyzer`. + +### Added +- Add `TryFrom` conversion between `UnitVector2/3/4` and `glam`’s `Vec2/3/4`. +- `nalgebra-sparse`: added support for serialization of sparse matrices with `serde`. +- `nalgebra-sparse`: add a CSC matrix constructor from unsorted (but valid) data. +- `nalgebra-lapack`: add generalized eigenvalues/eigenvectors calculation + QZ decomposition. + +### Fixed +- Improve stability of SVD. +- Fix slerp for `UnitComplex`. + ## [0.30.1] (09 Jan. 2022) + ### Added - Add conversion from/to types of `glam` 0.19 and 0.20. diff --git a/Cargo.toml b/Cargo.toml index 8a3fea5c..732676ec 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "nalgebra" -version = "0.30.1" +version = "0.31.0" authors = [ "Sébastien Crozet " ] description = "General-purpose linear algebra library with transformations and statically-sized or dynamically-sized matrices." @@ -37,7 +37,6 @@ cuda = [ "cust_core", "simba/cuda" ] # Conversion convert-mint = [ "mint" ] convert-bytemuck = [ "bytemuck" ] -convert-glam013 = [ "glam013" ] convert-glam014 = [ "glam014" ] convert-glam015 = [ "glam015" ] convert-glam016 = [ "glam016" ] @@ -54,7 +53,7 @@ convert-glam020 = [ "glam020" ] serde-serialize-no-std = [ "serde", "num-complex/serde" ] serde-serialize = [ "serde-serialize-no-std", "serde/std" ] rkyv-serialize-no-std = [ "rkyv" ] -rkyv-serialize = [ "rkyv-serialize-no-std", "rkyv/std" ] +rkyv-serialize = [ "rkyv-serialize-no-std", "rkyv/std", "bytecheck" ] # Randomness ## To use rand in a #[no-std] environment, enable the @@ -80,7 +79,8 @@ 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 } serde = { version = "1.0", default-features = false, features = [ "derive" ], optional = true } -rkyv = { version = "~0.6.4", default-features = false, features = ["const_generics"], optional = true } +rkyv = { version = "~0.7.1", optional = true } +bytecheck = { version = "~0.6.1", optional = true } mint = { version = "0.5", optional = true } quickcheck = { version = "1", optional = true } pest = { version = "2", optional = true } @@ -88,7 +88,6 @@ 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"] } -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 } diff --git a/examples/cargo/Cargo.toml b/examples/cargo/Cargo.toml index 501ae23e..63e70aab 100644 --- a/examples/cargo/Cargo.toml +++ b/examples/cargo/Cargo.toml @@ -4,7 +4,7 @@ version = "0.0.0" authors = [ "You" ] [dependencies] -nalgebra = "0.30.0" +nalgebra = "0.31.0" [[bin]] name = "example" diff --git a/nalgebra-glm/Cargo.toml b/nalgebra-glm/Cargo.toml index f8087581..e700af37 100644 --- a/nalgebra-glm/Cargo.toml +++ b/nalgebra-glm/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "nalgebra-glm" -version = "0.16.0" +version = "0.17.0" authors = ["sebcrozet "] description = "A computer-graphics oriented API for nalgebra, inspired by the C++ GLM library." @@ -26,7 +26,6 @@ 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" ] @@ -37,4 +36,4 @@ convert-glam018 = [ "nalgebra/glam018" ] num-traits = { version = "0.2", default-features = false } approx = { version = "0.5", default-features = false } simba = { version = "0.7", default-features = false } -nalgebra = { path = "..", version = "0.30", default-features = false } +nalgebra = { path = "..", version = "0.31", default-features = false } diff --git a/nalgebra-lapack/Cargo.toml b/nalgebra-lapack/Cargo.toml index 7cdf9f78..91517a8d 100644 --- a/nalgebra-lapack/Cargo.toml +++ b/nalgebra-lapack/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "nalgebra-lapack" -version = "0.21.0" +version = "0.22.0" authors = [ "Sébastien Crozet ", "Andrew Straw " ] description = "Matrix decompositions using nalgebra matrices and Lapack bindings." @@ -29,7 +29,7 @@ accelerate = ["lapack-src/accelerate"] intel-mkl = ["lapack-src/intel-mkl"] [dependencies] -nalgebra = { version = "0.30", path = ".." } +nalgebra = { version = "0.31", path = ".." } num-traits = "0.2" num-complex = { version = "0.4", default-features = false } simba = "0.7" @@ -39,7 +39,7 @@ lapack-src = { version = "0.8", default-features = false } # clippy = "*" [dev-dependencies] -nalgebra = { version = "0.30", features = [ "arbitrary", "rand" ], path = ".." } +nalgebra = { version = "0.31", features = [ "arbitrary", "rand" ], path = ".." } proptest = { version = "1", default-features = false, features = ["std"] } quickcheck = "1" approx = "0.5" diff --git a/nalgebra-lapack/src/generalized_eigenvalues.rs b/nalgebra-lapack/src/generalized_eigenvalues.rs new file mode 100644 index 00000000..5d1e3ace --- /dev/null +++ b/nalgebra-lapack/src/generalized_eigenvalues.rs @@ -0,0 +1,350 @@ +#[cfg(feature = "serde-serialize")] +use serde::{Deserialize, Serialize}; + +use num::Zero; +use num_complex::Complex; + +use simba::scalar::RealField; + +use crate::ComplexHelper; +use na::allocator::Allocator; +use na::dimension::{Const, Dim}; +use na::{DefaultAllocator, Matrix, OMatrix, OVector, Scalar}; + +use lapack; + +/// Generalized eigenvalues and generalized eigenvectors (left and right) of a pair of N*N real square matrices. +/// +/// Each generalized eigenvalue (lambda) satisfies determinant(A - lambda*B) = 0 +/// +/// The right eigenvector v(j) corresponding to the eigenvalue lambda(j) +/// of (A,B) satisfies +/// +/// A * v(j) = lambda(j) * B * v(j). +/// +/// The left eigenvector u(j) corresponding to the eigenvalue lambda(j) +/// of (A,B) satisfies +/// +/// u(j)**H * A = lambda(j) * u(j)**H * B . +/// where u(j)**H is the conjugate-transpose of u(j). +#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))] +#[cfg_attr( + feature = "serde-serialize", + serde( + bound(serialize = "DefaultAllocator: Allocator + Allocator, + OVector: Serialize, + OMatrix: Serialize") + ) +)] +#[cfg_attr( + feature = "serde-serialize", + serde( + bound(deserialize = "DefaultAllocator: Allocator + Allocator, + OVector: Deserialize<'de>, + OMatrix: Deserialize<'de>") + ) +)] +#[derive(Clone, Debug)] +pub struct GeneralizedEigen +where + DefaultAllocator: Allocator + Allocator, +{ + alphar: OVector, + alphai: OVector, + beta: OVector, + vsl: OMatrix, + vsr: OMatrix, +} + +impl Copy for GeneralizedEigen +where + DefaultAllocator: Allocator + Allocator, + OMatrix: Copy, + OVector: Copy, +{ +} + +impl GeneralizedEigen +where + DefaultAllocator: Allocator + Allocator, +{ + /// Attempts to compute the generalized eigenvalues, and left and right associated eigenvectors + /// via the raw returns from LAPACK's dggev and sggev routines + /// + /// Panics if the method did not converge. + pub fn new(a: OMatrix, b: OMatrix) -> Self { + Self::try_new(a, b).expect("Calculation of generalized eigenvalues failed.") + } + + /// Attempts to compute the generalized eigenvalues (and eigenvectors) via the raw returns from LAPACK's + /// dggev and sggev routines + /// + /// Returns `None` if the method did not converge. + pub fn try_new(mut a: OMatrix, mut b: OMatrix) -> Option { + assert!( + a.is_square() && b.is_square(), + "Unable to compute the generalized eigenvalues of non-square matrices." + ); + + assert!( + a.shape_generic() == b.shape_generic(), + "Unable to compute the generalized eigenvalues of two square matrices of different dimensions." + ); + + let (nrows, ncols) = a.shape_generic(); + let n = nrows.value(); + + let mut info = 0; + + let mut alphar = Matrix::zeros_generic(nrows, Const::<1>); + let mut alphai = Matrix::zeros_generic(nrows, Const::<1>); + let mut beta = Matrix::zeros_generic(nrows, Const::<1>); + let mut vsl = Matrix::zeros_generic(nrows, ncols); + let mut vsr = Matrix::zeros_generic(nrows, ncols); + + let lwork = T::xggev_work_size( + b'V', + b'V', + n as i32, + a.as_mut_slice(), + n as i32, + b.as_mut_slice(), + n as i32, + alphar.as_mut_slice(), + alphai.as_mut_slice(), + beta.as_mut_slice(), + vsl.as_mut_slice(), + n as i32, + vsr.as_mut_slice(), + n as i32, + &mut info, + ); + lapack_check!(info); + + let mut work = vec![T::zero(); lwork as usize]; + + T::xggev( + b'V', + b'V', + n as i32, + a.as_mut_slice(), + n as i32, + b.as_mut_slice(), + n as i32, + alphar.as_mut_slice(), + alphai.as_mut_slice(), + beta.as_mut_slice(), + vsl.as_mut_slice(), + n as i32, + vsr.as_mut_slice(), + n as i32, + &mut work, + lwork, + &mut info, + ); + lapack_check!(info); + + Some(GeneralizedEigen { + alphar, + alphai, + beta, + vsl, + vsr, + }) + } + + /// Calculates the generalized eigenvectors (left and right) associated with the generalized eigenvalues + /// + /// Outputs two matrices. + /// The first output matrix contains the left eigenvectors of the generalized eigenvalues + /// as columns. + /// The second matrix contains the right eigenvectors of the generalized eigenvalues + /// as columns. + pub fn eigenvectors(&self) -> (OMatrix, D, D>, OMatrix, D, D>) + where + DefaultAllocator: + Allocator, D, D> + Allocator, D> + Allocator<(Complex, T), D>, + { + /* + How the eigenvectors are built up: + + Since the input entries are all real, the generalized eigenvalues if complex come in pairs + as a consequence of the [complex conjugate root thorem](https://en.wikipedia.org/wiki/Complex_conjugate_root_theorem) + The Lapack routine output reflects this by expecting the user to unpack the real and complex eigenvalues associated + eigenvectors from the real matrix output via the following procedure + + (Note: VL stands for the lapack real matrix output containing the left eigenvectors as columns, + VR stands for the lapack real matrix output containing the right eigenvectors as columns) + + If the j-th and (j+1)-th eigenvalues form a complex conjugate pair, + then + + u(j) = VL(:,j)+i*VL(:,j+1) + u(j+1) = VL(:,j)-i*VL(:,j+1) + + and + + u(j) = VR(:,j)+i*VR(:,j+1) + v(j+1) = VR(:,j)-i*VR(:,j+1). + */ + + let n = self.vsl.shape().0; + + let mut l = self.vsl.map(|x| Complex::new(x, T::RealField::zero())); + + let mut r = self.vsr.map(|x| Complex::new(x, T::RealField::zero())); + + let eigenvalues = self.raw_eigenvalues(); + + let mut c = 0; + + while c < n { + if eigenvalues[c].0.im.abs() != T::RealField::zero() && c + 1 < n { + // taking care of the left eigenvector matrix + l.column_mut(c).zip_apply(&self.vsl.column(c + 1), |r, i| { + *r = Complex::new(r.re.clone(), i.clone()); + }); + l.column_mut(c + 1).zip_apply(&self.vsl.column(c), |i, r| { + *i = Complex::new(r.clone(), -i.re.clone()); + }); + + // taking care of the right eigenvector matrix + r.column_mut(c).zip_apply(&self.vsr.column(c + 1), |r, i| { + *r = Complex::new(r.re.clone(), i.clone()); + }); + r.column_mut(c + 1).zip_apply(&self.vsr.column(c), |i, r| { + *i = Complex::new(r.clone(), -i.re.clone()); + }); + + c += 2; + } else { + c += 1; + } + } + + (l, r) + } + + /// Outputs the unprocessed (almost) version of generalized eigenvalues ((alphar, alphai), beta) + /// straight from LAPACK + #[must_use] + pub fn raw_eigenvalues(&self) -> OVector<(Complex, T), D> + where + DefaultAllocator: Allocator<(Complex, T), D>, + { + let mut out = Matrix::from_element_generic( + self.vsl.shape_generic().0, + Const::<1>, + (Complex::zero(), T::RealField::zero()), + ); + + for i in 0..out.len() { + out[i] = (Complex::new(self.alphar[i], self.alphai[i]), self.beta[i]) + } + + out + } +} + +/* + * + * Lapack functions dispatch. + * + */ +/// Trait implemented by scalars for which Lapack implements the RealField GeneralizedEigen decomposition. +pub trait GeneralizedEigenScalar: Scalar { + #[allow(missing_docs)] + fn xggev( + jobvsl: u8, + jobvsr: u8, + n: i32, + a: &mut [Self], + lda: i32, + b: &mut [Self], + ldb: i32, + alphar: &mut [Self], + alphai: &mut [Self], + beta: &mut [Self], + vsl: &mut [Self], + ldvsl: i32, + vsr: &mut [Self], + ldvsr: i32, + work: &mut [Self], + lwork: i32, + info: &mut i32, + ); + + #[allow(missing_docs)] + fn xggev_work_size( + jobvsl: u8, + jobvsr: u8, + n: i32, + a: &mut [Self], + lda: i32, + b: &mut [Self], + ldb: i32, + alphar: &mut [Self], + alphai: &mut [Self], + beta: &mut [Self], + vsl: &mut [Self], + ldvsl: i32, + vsr: &mut [Self], + ldvsr: i32, + info: &mut i32, + ) -> i32; +} + +macro_rules! generalized_eigen_scalar_impl ( + ($N: ty, $xggev: path) => ( + impl GeneralizedEigenScalar for $N { + #[inline] + fn xggev(jobvsl: u8, + jobvsr: u8, + n: i32, + a: &mut [$N], + lda: i32, + b: &mut [$N], + ldb: i32, + alphar: &mut [$N], + alphai: &mut [$N], + beta : &mut [$N], + vsl: &mut [$N], + ldvsl: i32, + vsr: &mut [$N], + ldvsr: i32, + work: &mut [$N], + lwork: i32, + info: &mut i32) { + unsafe { $xggev(jobvsl, jobvsr, n, a, lda, b, ldb, alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr, work, lwork, info); } + } + + + #[inline] + fn xggev_work_size(jobvsl: u8, + jobvsr: u8, + n: i32, + a: &mut [$N], + lda: i32, + b: &mut [$N], + ldb: i32, + alphar: &mut [$N], + alphai: &mut [$N], + beta : &mut [$N], + vsl: &mut [$N], + ldvsl: i32, + vsr: &mut [$N], + ldvsr: i32, + info: &mut i32) + -> i32 { + let mut work = [ Zero::zero() ]; + let lwork = -1 as i32; + + unsafe { $xggev(jobvsl, jobvsr, n, a, lda, b, ldb, alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr, &mut work, lwork, info); } + ComplexHelper::real_part(work[0]) as i32 + } + } + ) +); + +generalized_eigen_scalar_impl!(f32, lapack::sggev); +generalized_eigen_scalar_impl!(f64, lapack::dggev); diff --git a/nalgebra-lapack/src/lib.rs b/nalgebra-lapack/src/lib.rs index 9cf0d73d..ea2e2b53 100644 --- a/nalgebra-lapack/src/lib.rs +++ b/nalgebra-lapack/src/lib.rs @@ -83,9 +83,11 @@ mod lapack_check; mod cholesky; mod eigen; +mod generalized_eigenvalues; mod hessenberg; mod lu; mod qr; +mod qz; mod schur; mod svd; mod symmetric_eigen; @@ -94,9 +96,11 @@ use num_complex::Complex; pub use self::cholesky::{Cholesky, CholeskyScalar}; pub use self::eigen::Eigen; +pub use self::generalized_eigenvalues::GeneralizedEigen; pub use self::hessenberg::Hessenberg; pub use self::lu::{LUScalar, LU}; pub use self::qr::QR; +pub use self::qz::QZ; pub use self::schur::Schur; pub use self::svd::SVD; pub use self::symmetric_eigen::SymmetricEigen; diff --git a/nalgebra-lapack/src/qz.rs b/nalgebra-lapack/src/qz.rs new file mode 100644 index 00000000..99f3c374 --- /dev/null +++ b/nalgebra-lapack/src/qz.rs @@ -0,0 +1,321 @@ +#[cfg(feature = "serde-serialize")] +use serde::{Deserialize, Serialize}; + +use num::Zero; +use num_complex::Complex; + +use simba::scalar::RealField; + +use crate::ComplexHelper; +use na::allocator::Allocator; +use na::dimension::{Const, Dim}; +use na::{DefaultAllocator, Matrix, OMatrix, OVector, Scalar}; + +use lapack; + +/// QZ decomposition of a pair of N*N square matrices. +/// +/// Retrieves the left and right matrices of Schur Vectors (VSL and VSR) +/// the upper-quasitriangular matrix `S` and upper triangular matrix `T` such that the +/// decomposed input matrix `a` equals `VSL * S * VSL.transpose()` and +/// decomposed input matrix `b` equals `VSL * T * VSL.transpose()`. +#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))] +#[cfg_attr( + feature = "serde-serialize", + serde( + bound(serialize = "DefaultAllocator: Allocator + Allocator, + OVector: Serialize, + OMatrix: Serialize") + ) +)] +#[cfg_attr( + feature = "serde-serialize", + serde( + bound(deserialize = "DefaultAllocator: Allocator + Allocator, + OVector: Deserialize<'de>, + OMatrix: Deserialize<'de>") + ) +)] +#[derive(Clone, Debug)] +pub struct QZ +where + DefaultAllocator: Allocator + Allocator, +{ + alphar: OVector, + alphai: OVector, + beta: OVector, + vsl: OMatrix, + s: OMatrix, + vsr: OMatrix, + t: OMatrix, +} + +impl Copy for QZ +where + DefaultAllocator: Allocator + Allocator, + OMatrix: Copy, + OVector: Copy, +{ +} + +impl QZ +where + DefaultAllocator: Allocator + Allocator, +{ + /// Attempts to compute the QZ decomposition of input real square matrices `a` and `b`. + /// + /// i.e retrieves the left and right matrices of Schur Vectors (VSL and VSR) + /// the upper-quasitriangular matrix `S` and upper triangular matrix `T` such that the + /// decomposed matrix `a` equals `VSL * S * VSL.transpose()` and + /// decomposed matrix `b` equals `VSL * T * VSL.transpose()`. + /// + /// Panics if the method did not converge. + pub fn new(a: OMatrix, b: OMatrix) -> Self { + Self::try_new(a, b).expect("QZ decomposition: convergence failed.") + } + + /// Computes the decomposition of input matrices `a` and `b` into a pair of matrices of Schur vectors + /// , a quasi-upper triangular matrix and an upper-triangular matrix . + /// + /// Returns `None` if the method did not converge. + pub fn try_new(mut a: OMatrix, mut b: OMatrix) -> Option { + assert!( + a.is_square() && b.is_square(), + "Unable to compute the qz decomposition of non-square matrices." + ); + + assert!( + a.shape_generic() == b.shape_generic(), + "Unable to compute the qz decomposition of two square matrices of different dimensions." + ); + + let (nrows, ncols) = a.shape_generic(); + let n = nrows.value(); + + let mut info = 0; + + let mut alphar = Matrix::zeros_generic(nrows, Const::<1>); + let mut alphai = Matrix::zeros_generic(nrows, Const::<1>); + let mut beta = Matrix::zeros_generic(nrows, Const::<1>); + let mut vsl = Matrix::zeros_generic(nrows, ncols); + let mut vsr = Matrix::zeros_generic(nrows, ncols); + // Placeholders: + let mut bwork = [0i32]; + let mut unused = 0; + + let lwork = T::xgges_work_size( + b'V', + b'V', + b'N', + n as i32, + a.as_mut_slice(), + n as i32, + b.as_mut_slice(), + n as i32, + &mut unused, + alphar.as_mut_slice(), + alphai.as_mut_slice(), + beta.as_mut_slice(), + vsl.as_mut_slice(), + n as i32, + vsr.as_mut_slice(), + n as i32, + &mut bwork, + &mut info, + ); + lapack_check!(info); + + let mut work = vec![T::zero(); lwork as usize]; + + T::xgges( + b'V', + b'V', + b'N', + n as i32, + a.as_mut_slice(), + n as i32, + b.as_mut_slice(), + n as i32, + &mut unused, + alphar.as_mut_slice(), + alphai.as_mut_slice(), + beta.as_mut_slice(), + vsl.as_mut_slice(), + n as i32, + vsr.as_mut_slice(), + n as i32, + &mut work, + lwork, + &mut bwork, + &mut info, + ); + lapack_check!(info); + + Some(QZ { + alphar, + alphai, + beta, + vsl, + s: a, + vsr, + t: b, + }) + } + + /// Retrieves the left and right matrices of Schur Vectors (VSL and VSR) + /// the upper-quasitriangular matrix `S` and upper triangular matrix `T` such that the + /// decomposed input matrix `a` equals `VSL * S * VSL.transpose()` and + /// decomposed input matrix `b` equals `VSL * T * VSL.transpose()`. + pub fn unpack( + self, + ) -> ( + OMatrix, + OMatrix, + OMatrix, + OMatrix, + ) { + (self.vsl, self.s, self.t, self.vsr) + } + + /// outputs the unprocessed (almost) version of generalized eigenvalues ((alphar, alpai), beta) + /// straight from LAPACK + #[must_use] + pub fn raw_eigenvalues(&self) -> OVector<(Complex, T), D> + where + DefaultAllocator: Allocator<(Complex, T), D>, + { + let mut out = Matrix::from_element_generic( + self.vsl.shape_generic().0, + Const::<1>, + (Complex::zero(), T::RealField::zero()), + ); + + for i in 0..out.len() { + out[i] = ( + Complex::new(self.alphar[i].clone(), self.alphai[i].clone()), + self.beta[i].clone(), + ) + } + + out + } +} + +/* + * + * Lapack functions dispatch. + * + */ +/// Trait implemented by scalars for which Lapack implements the RealField QZ decomposition. +pub trait QZScalar: Scalar { + #[allow(missing_docs)] + fn xgges( + jobvsl: u8, + jobvsr: u8, + sort: u8, + // select: ??? + n: i32, + a: &mut [Self], + lda: i32, + b: &mut [Self], + ldb: i32, + sdim: &mut i32, + alphar: &mut [Self], + alphai: &mut [Self], + beta: &mut [Self], + vsl: &mut [Self], + ldvsl: i32, + vsr: &mut [Self], + ldvsr: i32, + work: &mut [Self], + lwork: i32, + bwork: &mut [i32], + info: &mut i32, + ); + + #[allow(missing_docs)] + fn xgges_work_size( + jobvsl: u8, + jobvsr: u8, + sort: u8, + // select: ??? + n: i32, + a: &mut [Self], + lda: i32, + b: &mut [Self], + ldb: i32, + sdim: &mut i32, + alphar: &mut [Self], + alphai: &mut [Self], + beta: &mut [Self], + vsl: &mut [Self], + ldvsl: i32, + vsr: &mut [Self], + ldvsr: i32, + bwork: &mut [i32], + info: &mut i32, + ) -> i32; +} + +macro_rules! qz_scalar_impl ( + ($N: ty, $xgges: path) => ( + impl QZScalar for $N { + #[inline] + fn xgges(jobvsl: u8, + jobvsr: u8, + sort: u8, + // select: ??? + n: i32, + a: &mut [$N], + lda: i32, + b: &mut [$N], + ldb: i32, + sdim: &mut i32, + alphar: &mut [$N], + alphai: &mut [$N], + beta : &mut [$N], + vsl: &mut [$N], + ldvsl: i32, + vsr: &mut [$N], + ldvsr: i32, + work: &mut [$N], + lwork: i32, + bwork: &mut [i32], + info: &mut i32) { + unsafe { $xgges(jobvsl, jobvsr, sort, None, n, a, lda, b, ldb, sdim, alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr, work, lwork, bwork, info); } + } + + + #[inline] + fn xgges_work_size(jobvsl: u8, + jobvsr: u8, + sort: u8, + // select: ??? + n: i32, + a: &mut [$N], + lda: i32, + b: &mut [$N], + ldb: i32, + sdim: &mut i32, + alphar: &mut [$N], + alphai: &mut [$N], + beta : &mut [$N], + vsl: &mut [$N], + ldvsl: i32, + vsr: &mut [$N], + ldvsr: i32, + bwork: &mut [i32], + info: &mut i32) + -> i32 { + let mut work = [ Zero::zero() ]; + let lwork = -1 as i32; + + unsafe { $xgges(jobvsl, jobvsr, sort, None, n, a, lda, b, ldb, sdim, alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr, &mut work, lwork, bwork, info); } + ComplexHelper::real_part(work[0]) as i32 + } + } + ) +); + +qz_scalar_impl!(f32, lapack::sgges); +qz_scalar_impl!(f64, lapack::dgges); diff --git a/nalgebra-lapack/tests/linalg/generalized_eigenvalues.rs b/nalgebra-lapack/tests/linalg/generalized_eigenvalues.rs new file mode 100644 index 00000000..b0d9777c --- /dev/null +++ b/nalgebra-lapack/tests/linalg/generalized_eigenvalues.rs @@ -0,0 +1,72 @@ +use na::dimension::Const; +use na::{DMatrix, OMatrix}; +use nl::GeneralizedEigen; +use num_complex::Complex; +use simba::scalar::ComplexField; + +use crate::proptest::*; +use proptest::{prop_assert, prop_compose, proptest}; + +prop_compose! { + fn f64_dynamic_dim_squares() + (n in PROPTEST_MATRIX_DIM) + (a in matrix(PROPTEST_F64,n,n), b in matrix(PROPTEST_F64,n,n)) -> (DMatrix, DMatrix){ + (a,b) +}} + +proptest! { + #[test] + fn ge((a,b) in f64_dynamic_dim_squares()){ + + let a_c = a.clone().map(|x| Complex::new(x, 0.0)); + let b_c = b.clone().map(|x| Complex::new(x, 0.0)); + let n = a.shape_generic().0; + + let ge = GeneralizedEigen::new(a.clone(), b.clone()); + let (vsl,vsr) = ge.clone().eigenvectors(); + + + for (i,(alpha,beta)) in ge.raw_eigenvalues().iter().enumerate() { + let l_a = a_c.clone() * Complex::new(*beta, 0.0); + let l_b = b_c.clone() * *alpha; + + prop_assert!( + relative_eq!( + ((&l_a - &l_b)*vsr.column(i)).map(|x| x.modulus()), + OMatrix::zeros_generic(n, Const::<1>), + epsilon = 1.0e-5)); + + prop_assert!( + relative_eq!( + (vsl.column(i).adjoint()*(&l_a - &l_b)).map(|x| x.modulus()), + OMatrix::zeros_generic(Const::<1>, n), + epsilon = 1.0e-5)) + }; + } + + #[test] + fn ge_static(a in matrix4(), b in matrix4()) { + + let ge = GeneralizedEigen::new(a.clone(), b.clone()); + let a_c =a.clone().map(|x| Complex::new(x, 0.0)); + let b_c = b.clone().map(|x| Complex::new(x, 0.0)); + let (vsl,vsr) = ge.eigenvectors(); + let eigenvalues = ge.raw_eigenvalues(); + + for (i,(alpha,beta)) in eigenvalues.iter().enumerate() { + let l_a = a_c.clone() * Complex::new(*beta, 0.0); + let l_b = b_c.clone() * *alpha; + + prop_assert!( + relative_eq!( + ((&l_a - &l_b)*vsr.column(i)).map(|x| x.modulus()), + OMatrix::zeros_generic(Const::<4>, Const::<1>), + epsilon = 1.0e-5)); + prop_assert!( + relative_eq!((vsl.column(i).adjoint()*(&l_a - &l_b)).map(|x| x.modulus()), + OMatrix::zeros_generic(Const::<1>, Const::<4>), + epsilon = 1.0e-5)) + } + } + +} diff --git a/nalgebra-lapack/tests/linalg/mod.rs b/nalgebra-lapack/tests/linalg/mod.rs index a6742217..251bbe7b 100644 --- a/nalgebra-lapack/tests/linalg/mod.rs +++ b/nalgebra-lapack/tests/linalg/mod.rs @@ -1,6 +1,8 @@ mod cholesky; +mod generalized_eigenvalues; mod lu; mod qr; +mod qz; mod real_eigensystem; mod schur; mod svd; diff --git a/nalgebra-lapack/tests/linalg/qz.rs b/nalgebra-lapack/tests/linalg/qz.rs new file mode 100644 index 00000000..c7a702ca --- /dev/null +++ b/nalgebra-lapack/tests/linalg/qz.rs @@ -0,0 +1,34 @@ +use na::DMatrix; +use nl::QZ; + +use crate::proptest::*; +use proptest::{prop_assert, prop_compose, proptest}; + +prop_compose! { + fn f64_dynamic_dim_squares() + (n in PROPTEST_MATRIX_DIM) + (a in matrix(PROPTEST_F64,n,n), b in matrix(PROPTEST_F64,n,n)) -> (DMatrix, DMatrix){ + (a,b) +}} + +proptest! { + #[test] + fn qz((a,b) in f64_dynamic_dim_squares()) { + + let qz = QZ::new(a.clone(), b.clone()); + let (vsl,s,t,vsr) = qz.clone().unpack(); + + prop_assert!(relative_eq!(&vsl * s * vsr.transpose(), a, epsilon = 1.0e-7)); + prop_assert!(relative_eq!(vsl * t * vsr.transpose(), b, epsilon = 1.0e-7)); + + } + + #[test] + fn qz_static(a in matrix4(), b in matrix4()) { + let qz = QZ::new(a.clone(), b.clone()); + let (vsl,s,t,vsr) = qz.unpack(); + + prop_assert!(relative_eq!(&vsl * s * vsr.transpose(), a, epsilon = 1.0e-7)); + prop_assert!(relative_eq!(vsl * t * vsr.transpose(), b, epsilon = 1.0e-7)); + } +} diff --git a/nalgebra-macros/Cargo.toml b/nalgebra-macros/Cargo.toml index bd37b689..6e35f9e9 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.30.0", path = ".." } +nalgebra = { version = "0.31.0", path = ".." } trybuild = "1.0.42" diff --git a/nalgebra-sparse/Cargo.toml b/nalgebra-sparse/Cargo.toml index 0ed64acb..390f594e 100644 --- a/nalgebra-sparse/Cargo.toml +++ b/nalgebra-sparse/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "nalgebra-sparse" -version = "0.6.0" +version = "0.7.0" authors = [ "Andreas Longva", "Sébastien Crozet " ] edition = "2018" description = "Sparse matrix computation based on nalgebra." @@ -24,7 +24,7 @@ io = [ "pest", "pest_derive" ] slow-tests = [] [dependencies] -nalgebra = { version="0.30", path = "../" } +nalgebra = { version="0.31", path = "../" } num-traits = { version = "0.2", default-features = false } proptest = { version = "1.0", optional = true } matrixcompare-core = { version = "0.1.0", optional = true } @@ -35,8 +35,8 @@ serde = { version = "1.0", default-features = false, features = [ "derive" ], op [dev-dependencies] itertools = "0.10" matrixcompare = { version = "0.3.0", features = [ "proptest-support" ] } -nalgebra = { version="0.30", path = "../", features = ["compare"] } -tempfile = "3" +nalgebra = { version="0.31", path = "../", features = ["compare"] } +tempfile = "3.3" serde_json = "1.0" [package.metadata.docs.rs] diff --git a/nalgebra-sparse/src/ops/impl_std_ops.rs b/nalgebra-sparse/src/ops/impl_std_ops.rs index 107c38ba..91b6574f 100644 --- a/nalgebra-sparse/src/ops/impl_std_ops.rs +++ b/nalgebra-sparse/src/ops/impl_std_ops.rs @@ -3,7 +3,7 @@ use crate::csr::CsrMatrix; use crate::ops::serial::{ spadd_csc_prealloc, spadd_csr_prealloc, spadd_pattern, spmm_csc_dense, spmm_csc_pattern, - spmm_csc_prealloc, spmm_csr_dense, spmm_csr_pattern, spmm_csr_prealloc, + spmm_csc_prealloc_unchecked, spmm_csr_dense, spmm_csr_pattern, spmm_csr_prealloc_unchecked, }; use crate::ops::Op; use nalgebra::allocator::Allocator; @@ -112,9 +112,9 @@ macro_rules! impl_spmm { } } -impl_spmm!(CsrMatrix, spmm_csr_pattern, spmm_csr_prealloc); +impl_spmm!(CsrMatrix, spmm_csr_pattern, spmm_csr_prealloc_unchecked); // Need to switch order of operations for CSC pattern -impl_spmm!(CscMatrix, spmm_csc_pattern, spmm_csc_prealloc); +impl_spmm!(CscMatrix, spmm_csc_pattern, spmm_csc_prealloc_unchecked); /// Implements Scalar * Matrix operations for *concrete* scalar types. The reason this is necessary /// is that we are not able to implement Mul> for all T generically due to orphan rules. diff --git a/nalgebra-sparse/src/ops/serial/cs.rs b/nalgebra-sparse/src/ops/serial/cs.rs index 86484053..cc13c168 100644 --- a/nalgebra-sparse/src/ops/serial/cs.rs +++ b/nalgebra-sparse/src/ops/serial/cs.rs @@ -20,6 +20,51 @@ fn spmm_cs_unexpected_entry() -> OperationError { /// reversed (since transpose(AB) = transpose(B) * transpose(A) and CSC(A) = transpose(CSR(A)). /// /// We assume here that the matrices have already been verified to be dimensionally compatible. +pub fn spmm_cs_prealloc_unchecked( + beta: T, + c: &mut CsMatrix, + alpha: T, + a: &CsMatrix, + b: &CsMatrix, +) -> Result<(), OperationError> +where + T: Scalar + ClosedAdd + ClosedMul + Zero + One, +{ + assert_eq!(c.pattern().major_dim(), a.pattern().major_dim()); + assert_eq!(c.pattern().minor_dim(), b.pattern().minor_dim()); + let some_val = Zero::zero(); + let mut scratchpad_values: Vec = vec![some_val; b.pattern().minor_dim()]; + for i in 0..c.pattern().major_dim() { + let a_lane_i = a.get_lane(i).unwrap(); + + let mut c_lane_i = c.get_lane_mut(i).unwrap(); + + for (&k, a_ik) in a_lane_i.minor_indices().iter().zip(a_lane_i.values()) { + let b_lane_k = b.get_lane(k).unwrap(); + let alpha_aik = alpha.clone() * a_ik.clone(); + for (j, b_kj) in b_lane_k.minor_indices().iter().zip(b_lane_k.values()) { + // use a dense scatter vector to accumulate non-zeros quickly + unsafe { + *scratchpad_values.get_unchecked_mut(*j) += alpha_aik.clone() * b_kj.clone(); + } + } + } + + //Get indices from C pattern and gather from the dense scratchpad_values + let (indices, values) = c_lane_i.indices_and_values_mut(); + values + .iter_mut() + .zip(indices) + .for_each(|(output_ref, index)| unsafe { + *output_ref = beta.clone() * output_ref.clone() + + scratchpad_values.get_unchecked(*index).clone(); + *scratchpad_values.get_unchecked_mut(*index) = Zero::zero(); + }); + } + + Ok(()) +} + pub fn spmm_cs_prealloc( beta: T, c: &mut CsMatrix, diff --git a/nalgebra-sparse/src/ops/serial/csc.rs b/nalgebra-sparse/src/ops/serial/csc.rs index e5c9ae4e..6a691ef9 100644 --- a/nalgebra-sparse/src/ops/serial/csc.rs +++ b/nalgebra-sparse/src/ops/serial/csc.rs @@ -1,5 +1,7 @@ use crate::csc::CscMatrix; -use crate::ops::serial::cs::{spadd_cs_prealloc, spmm_cs_dense, spmm_cs_prealloc}; +use crate::ops::serial::cs::{ + spadd_cs_prealloc, spmm_cs_dense, spmm_cs_prealloc, spmm_cs_prealloc_unchecked, +}; use crate::ops::serial::{OperationError, OperationErrorKind}; use crate::ops::Op; use nalgebra::{ClosedAdd, ClosedMul, DMatrixSlice, DMatrixSliceMut, RealField, Scalar}; @@ -83,35 +85,81 @@ where { assert_compatible_spmm_dims!(c, a, b); - use Op::{NoOp, Transpose}; + use Op::NoOp; match (&a, &b) { (NoOp(ref a), NoOp(ref b)) => { // Note: We have to reverse the order for CSC matrices spmm_cs_prealloc(beta, &mut c.cs, alpha, &b.cs, &a.cs) } - _ => { - // Currently we handle transposition by explicitly precomputing transposed matrices - // and calling the operation again without transposition - let a_ref: &CscMatrix = a.inner_ref(); - let b_ref: &CscMatrix = b.inner_ref(); - let (a, b) = { - use Cow::*; - match (&a, &b) { - (NoOp(_), NoOp(_)) => unreachable!(), - (Transpose(ref a), NoOp(_)) => (Owned(a.transpose()), Borrowed(b_ref)), - (NoOp(_), Transpose(ref b)) => (Borrowed(a_ref), Owned(b.transpose())), - (Transpose(ref a), Transpose(ref b)) => { - (Owned(a.transpose()), Owned(b.transpose())) - } - } - }; - - spmm_csc_prealloc(beta, c, alpha, NoOp(a.as_ref()), NoOp(b.as_ref())) - } + _ => spmm_csc_transposed(beta, c, alpha, a, b, spmm_csc_prealloc), } } +/// Faster sparse-sparse matrix multiplication, `C <- beta * C + alpha * op(A) * op(B)`. +/// This will not return an error even if the patterns don't match. +/// Should be used for situations where pattern creation immediately preceeds multiplication. +/// +/// Panics if the dimensions of the matrices involved are not compatible with the expression. +pub fn spmm_csc_prealloc_unchecked( + beta: T, + c: &mut CscMatrix, + alpha: T, + a: Op<&CscMatrix>, + b: Op<&CscMatrix>, +) -> Result<(), OperationError> +where + T: Scalar + ClosedAdd + ClosedMul + Zero + One, +{ + assert_compatible_spmm_dims!(c, a, b); + + use Op::NoOp; + + match (&a, &b) { + (NoOp(ref a), NoOp(ref b)) => { + // Note: We have to reverse the order for CSC matrices + spmm_cs_prealloc_unchecked(beta, &mut c.cs, alpha, &b.cs, &a.cs) + } + _ => spmm_csc_transposed(beta, c, alpha, a, b, spmm_csc_prealloc_unchecked), + } +} + +fn spmm_csc_transposed( + beta: T, + c: &mut CscMatrix, + alpha: T, + a: Op<&CscMatrix>, + b: Op<&CscMatrix>, + spmm_kernel: F, +) -> Result<(), OperationError> +where + T: Scalar + ClosedAdd + ClosedMul + Zero + One, + F: Fn( + T, + &mut CscMatrix, + T, + Op<&CscMatrix>, + Op<&CscMatrix>, + ) -> Result<(), OperationError>, +{ + use Op::{NoOp, Transpose}; + + // Currently we handle transposition by explicitly precomputing transposed matrices + // and calling the operation again without transposition + let a_ref: &CscMatrix = a.inner_ref(); + let b_ref: &CscMatrix = b.inner_ref(); + let (a, b) = { + use Cow::*; + match (&a, &b) { + (NoOp(_), NoOp(_)) => unreachable!(), + (Transpose(ref a), NoOp(_)) => (Owned(a.transpose()), Borrowed(b_ref)), + (NoOp(_), Transpose(ref b)) => (Borrowed(a_ref), Owned(b.transpose())), + (Transpose(ref a), Transpose(ref b)) => (Owned(a.transpose()), Owned(b.transpose())), + } + }; + spmm_kernel(beta, c, alpha, NoOp(a.as_ref()), NoOp(b.as_ref())) +} + /// Solve the lower triangular system `op(L) X = B`. /// /// Only the lower triangular part of L is read, and the result is stored in B. diff --git a/nalgebra-sparse/src/ops/serial/csr.rs b/nalgebra-sparse/src/ops/serial/csr.rs index fa317bbf..708c81a3 100644 --- a/nalgebra-sparse/src/ops/serial/csr.rs +++ b/nalgebra-sparse/src/ops/serial/csr.rs @@ -1,5 +1,7 @@ use crate::csr::CsrMatrix; -use crate::ops::serial::cs::{spadd_cs_prealloc, spmm_cs_dense, spmm_cs_prealloc}; +use crate::ops::serial::cs::{ + spadd_cs_prealloc, spmm_cs_dense, spmm_cs_prealloc, spmm_cs_prealloc_unchecked, +}; use crate::ops::serial::OperationError; use crate::ops::Op; use nalgebra::{ClosedAdd, ClosedMul, DMatrixSlice, DMatrixSliceMut, Scalar}; @@ -77,30 +79,73 @@ where { assert_compatible_spmm_dims!(c, a, b); - use Op::{NoOp, Transpose}; + use Op::NoOp; match (&a, &b) { (NoOp(ref a), NoOp(ref b)) => spmm_cs_prealloc(beta, &mut c.cs, alpha, &a.cs, &b.cs), - _ => { - // Currently we handle transposition by explicitly precomputing transposed matrices - // and calling the operation again without transposition - // TODO: At least use workspaces to allow control of allocations. Maybe - // consider implementing certain patterns (like A^T * B) explicitly - let a_ref: &CsrMatrix = a.inner_ref(); - let b_ref: &CsrMatrix = b.inner_ref(); - let (a, b) = { - use Cow::*; - match (&a, &b) { - (NoOp(_), NoOp(_)) => unreachable!(), - (Transpose(ref a), NoOp(_)) => (Owned(a.transpose()), Borrowed(b_ref)), - (NoOp(_), Transpose(ref b)) => (Borrowed(a_ref), Owned(b.transpose())), - (Transpose(ref a), Transpose(ref b)) => { - (Owned(a.transpose()), Owned(b.transpose())) - } - } - }; - - spmm_csr_prealloc(beta, c, alpha, NoOp(a.as_ref()), NoOp(b.as_ref())) - } + _ => spmm_csr_transposed(beta, c, alpha, a, b, spmm_csr_prealloc), } } + +/// Faster sparse-sparse matrix multiplication, `C <- beta * C + alpha * op(A) * op(B)`. +/// This will not return an error even if the patterns don't match. +/// Should be used for situations where pattern creation immediately preceeds multiplication. +/// +/// Panics if the dimensions of the matrices involved are not compatible with the expression. +pub fn spmm_csr_prealloc_unchecked( + beta: T, + c: &mut CsrMatrix, + alpha: T, + a: Op<&CsrMatrix>, + b: Op<&CsrMatrix>, +) -> Result<(), OperationError> +where + T: Scalar + ClosedAdd + ClosedMul + Zero + One, +{ + assert_compatible_spmm_dims!(c, a, b); + + use Op::NoOp; + + match (&a, &b) { + (NoOp(ref a), NoOp(ref b)) => { + spmm_cs_prealloc_unchecked(beta, &mut c.cs, alpha, &a.cs, &b.cs) + } + _ => spmm_csr_transposed(beta, c, alpha, a, b, spmm_csr_prealloc_unchecked), + } +} + +fn spmm_csr_transposed( + beta: T, + c: &mut CsrMatrix, + alpha: T, + a: Op<&CsrMatrix>, + b: Op<&CsrMatrix>, + spmm_kernel: F, +) -> Result<(), OperationError> +where + T: Scalar + ClosedAdd + ClosedMul + Zero + One, + F: Fn( + T, + &mut CsrMatrix, + T, + Op<&CsrMatrix>, + Op<&CsrMatrix>, + ) -> Result<(), OperationError>, +{ + use Op::{NoOp, Transpose}; + + // Currently we handle transposition by explicitly precomputing transposed matrices + // and calling the operation again without transposition + let a_ref: &CsrMatrix = a.inner_ref(); + let b_ref: &CsrMatrix = b.inner_ref(); + let (a, b) = { + use Cow::*; + match (&a, &b) { + (NoOp(_), NoOp(_)) => unreachable!(), + (Transpose(ref a), NoOp(_)) => (Owned(a.transpose()), Borrowed(b_ref)), + (NoOp(_), Transpose(ref b)) => (Borrowed(a_ref), Owned(b.transpose())), + (Transpose(ref a), Transpose(ref b)) => (Owned(a.transpose()), Owned(b.transpose())), + } + }; + spmm_kernel(beta, c, alpha, NoOp(a.as_ref()), NoOp(b.as_ref())) +} diff --git a/nalgebra-sparse/tests/unit_tests/ops.rs b/nalgebra-sparse/tests/unit_tests/ops.rs index f2a02fd8..0a335567 100644 --- a/nalgebra-sparse/tests/unit_tests/ops.rs +++ b/nalgebra-sparse/tests/unit_tests/ops.rs @@ -6,7 +6,8 @@ use nalgebra_sparse::csc::CscMatrix; use nalgebra_sparse::csr::CsrMatrix; use nalgebra_sparse::ops::serial::{ spadd_csc_prealloc, spadd_csr_prealloc, spadd_pattern, spmm_csc_dense, spmm_csc_prealloc, - spmm_csr_dense, spmm_csr_pattern, spmm_csr_prealloc, spsolve_csc_lower_triangular, + spmm_csc_prealloc_unchecked, spmm_csr_dense, spmm_csr_pattern, spmm_csr_prealloc, + spmm_csr_prealloc_unchecked, spsolve_csc_lower_triangular, }; use nalgebra_sparse::ops::Op; use nalgebra_sparse::pattern::SparsityPattern; @@ -543,6 +544,29 @@ proptest! { prop_assert_eq!(&c_pattern, c_csr.pattern()); } + #[test] + fn spmm_csr_prealloc_unchecked_test(SpmmCsrArgs { c, beta, alpha, a, b } + in spmm_csr_prealloc_args_strategy() + ) { + // Test that we get the expected result by comparing to an equivalent dense operation + // (here we give in the C matrix, so the sparsity pattern is essentially fixed) + let mut c_sparse = c.clone(); + spmm_csr_prealloc_unchecked(beta, &mut c_sparse, alpha, a.as_ref(), b.as_ref()).unwrap(); + + let mut c_dense = DMatrix::from(&c); + let op_a_dense = match a { + Op::NoOp(ref a) => DMatrix::from(a), + Op::Transpose(ref a) => DMatrix::from(a).transpose(), + }; + let op_b_dense = match b { + Op::NoOp(ref b) => DMatrix::from(b), + Op::Transpose(ref b) => DMatrix::from(b).transpose(), + }; + c_dense = beta * c_dense + alpha * &op_a_dense * op_b_dense; + + prop_assert_eq!(&DMatrix::from(&c_sparse), &c_dense); + } + #[test] fn spmm_csr_prealloc_test(SpmmCsrArgs { c, beta, alpha, a, b } in spmm_csr_prealloc_args_strategy() @@ -705,6 +729,29 @@ proptest! { prop_assert_eq!(&DMatrix::from(&c_sparse), &c_dense); } + #[test] + fn spmm_csc_prealloc_unchecked_test(SpmmCscArgs { c, beta, alpha, a, b } + in spmm_csc_prealloc_args_strategy() + ) { + // Test that we get the expected result by comparing to an equivalent dense operation + // (here we give in the C matrix, so the sparsity pattern is essentially fixed) + let mut c_sparse = c.clone(); + spmm_csc_prealloc_unchecked(beta, &mut c_sparse, alpha, a.as_ref(), b.as_ref()).unwrap(); + + let mut c_dense = DMatrix::from(&c); + let op_a_dense = match a { + Op::NoOp(ref a) => DMatrix::from(a), + Op::Transpose(ref a) => DMatrix::from(a).transpose(), + }; + let op_b_dense = match b { + Op::NoOp(ref b) => DMatrix::from(b), + Op::Transpose(ref b) => DMatrix::from(b).transpose(), + }; + c_dense = beta * c_dense + alpha * &op_a_dense * op_b_dense; + + prop_assert_eq!(&DMatrix::from(&c_sparse), &c_dense); + } + #[test] fn spmm_csc_prealloc_panics_on_dim_mismatch( (alpha, beta, c, a, b) diff --git a/src/base/array_storage.rs b/src/base/array_storage.rs index b6bd236a..3bc71e1a 100644 --- a/src/base/array_storage.rs +++ b/src/base/array_storage.rs @@ -27,6 +27,11 @@ use std::mem; /// A array-based statically sized matrix data storage. #[repr(transparent)] #[derive(Copy, Clone, PartialEq, Eq, Hash)] +#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))] +#[cfg_attr( + feature = "rkyv-serialize-no-std", + derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize) +)] #[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))] pub struct ArrayStorage(pub [[T; R]; C]); @@ -273,45 +278,3 @@ unsafe impl by for ArrayStorage { } - -#[cfg(feature = "rkyv-serialize-no-std")] -mod rkyv_impl { - use super::ArrayStorage; - use rkyv::{offset_of, project_struct, Archive, Deserialize, Fallible, Serialize}; - - impl Archive for ArrayStorage { - type Archived = ArrayStorage; - type Resolver = <[[T; R]; C] as Archive>::Resolver; - - fn resolve( - &self, - pos: usize, - resolver: Self::Resolver, - out: &mut core::mem::MaybeUninit, - ) { - self.0.resolve( - pos + offset_of!(Self::Archived, 0), - resolver, - project_struct!(out: Self::Archived => 0), - ); - } - } - - impl, S: Fallible + ?Sized, const R: usize, const C: usize> Serialize - for ArrayStorage - { - fn serialize(&self, serializer: &mut S) -> Result { - self.0.serialize(serializer) - } - } - - impl - Deserialize, D> for ArrayStorage - where - T::Archived: Deserialize, - { - fn deserialize(&self, deserializer: &mut D) -> Result, D::Error> { - Ok(ArrayStorage(self.0.deserialize(deserializer)?)) - } - } -} diff --git a/src/base/dimension.rs b/src/base/dimension.rs index de51339f..4be97586 100644 --- a/src/base/dimension.rs +++ b/src/base/dimension.rs @@ -13,6 +13,11 @@ use serde::{Deserialize, Deserializer, Serialize, Serializer}; /// Dim of dynamically-sized algebraic entities. #[derive(Clone, Copy, Eq, PartialEq, Debug)] +#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))] +#[cfg_attr( + feature = "rkyv-serialize-no-std", + derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize) +)] #[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))] pub struct Dynamic { value: usize, @@ -198,6 +203,11 @@ dim_ops!( ); #[derive(Debug, Copy, Clone, PartialEq, Eq, Hash)] +#[cfg_attr( + feature = "rkyv-serialize-no-std", + derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize) +)] +#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))] #[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))] pub struct Const; @@ -233,37 +243,6 @@ impl<'de, const D: usize> Deserialize<'de> for Const { } } -#[cfg(feature = "rkyv-serialize-no-std")] -mod rkyv_impl { - use super::Const; - use rkyv::{Archive, Deserialize, Fallible, Serialize}; - - impl Archive for Const { - type Archived = Self; - type Resolver = (); - - fn resolve( - &self, - _: usize, - _: Self::Resolver, - _: &mut core::mem::MaybeUninit, - ) { - } - } - - impl Serialize for Const { - fn serialize(&self, _: &mut S) -> Result { - Ok(()) - } - } - - impl Deserialize for Const { - fn deserialize(&self, _: &mut D) -> Result { - Ok(Const) - } - } -} - pub trait ToConst { type Const: DimName; } diff --git a/src/base/matrix.rs b/src/base/matrix.rs index f12cb3fa..8f8786c1 100644 --- a/src/base/matrix.rs +++ b/src/base/matrix.rs @@ -150,6 +150,11 @@ pub type MatrixCross = /// some concrete types for `T` and a compatible data storage type `S`). #[repr(C)] #[derive(Clone, Copy)] +#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))] +#[cfg_attr( + feature = "rkyv-serialize-no-std", + derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize) +)] #[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))] pub struct Matrix { /// The data storage that contains all the matrix components. Disappointed? @@ -288,53 +293,6 @@ where { } -#[cfg(feature = "rkyv-serialize-no-std")] -mod rkyv_impl { - use super::Matrix; - use core::marker::PhantomData; - use rkyv::{offset_of, project_struct, Archive, Deserialize, Fallible, Serialize}; - - impl Archive for Matrix { - type Archived = Matrix; - type Resolver = S::Resolver; - - fn resolve( - &self, - pos: usize, - resolver: Self::Resolver, - out: &mut core::mem::MaybeUninit, - ) { - self.data.resolve( - pos + offset_of!(Self::Archived, data), - resolver, - project_struct!(out: Self::Archived => data), - ); - } - } - - impl, _S: Fallible + ?Sized> Serialize<_S> - for Matrix - { - fn serialize(&self, serializer: &mut _S) -> Result { - self.data.serialize(serializer) - } - } - - impl - Deserialize, D> - for Matrix - where - S::Archived: Deserialize, - { - fn deserialize(&self, deserializer: &mut D) -> Result, D::Error> { - Ok(Matrix { - data: self.data.deserialize(deserializer)?, - _phantoms: PhantomData, - }) - } - } -} - impl Matrix { /// Creates a new matrix with the given data without statically checking that the matrix /// dimension matches the storage dimension. diff --git a/src/base/unit.rs b/src/base/unit.rs index bb8b56a1..6fc00092 100644 --- a/src/base/unit.rs +++ b/src/base/unit.rs @@ -21,6 +21,11 @@ use crate::{Dim, Matrix, OMatrix, RealField, Scalar, SimdComplexField, SimdRealF /// in their documentation, read their dedicated pages directly. #[repr(transparent)] #[derive(Clone, Hash, Copy)] +#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))] +#[cfg_attr( + feature = "rkyv-serialize-no-std", + derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize) +)] // #[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))] pub struct Unit { pub(crate) value: T, @@ -58,47 +63,6 @@ impl<'de, T: Deserialize<'de>> Deserialize<'de> for Unit { } } -#[cfg(feature = "rkyv-serialize-no-std")] -mod rkyv_impl { - use super::Unit; - use rkyv::{offset_of, project_struct, Archive, Deserialize, Fallible, Serialize}; - - impl Archive for Unit { - type Archived = Unit; - type Resolver = T::Resolver; - - fn resolve( - &self, - pos: usize, - resolver: Self::Resolver, - out: &mut ::core::mem::MaybeUninit, - ) { - self.value.resolve( - pos + offset_of!(Self::Archived, value), - resolver, - project_struct!(out: Self::Archived => value), - ); - } - } - - impl, S: Fallible + ?Sized> Serialize for Unit { - fn serialize(&self, serializer: &mut S) -> Result { - self.value.serialize(serializer) - } - } - - impl Deserialize, D> for Unit - where - T::Archived: Deserialize, - { - fn deserialize(&self, deserializer: &mut D) -> Result, D::Error> { - Ok(Unit { - value: self.value.deserialize(deserializer)?, - }) - } - } -} - #[cfg(feature = "cuda")] unsafe impl cust_core::DeviceCopy for Unit> where diff --git a/src/geometry/dual_quaternion.rs b/src/geometry/dual_quaternion.rs index 719ae13d..6f1b7053 100644 --- a/src/geometry/dual_quaternion.rs +++ b/src/geometry/dual_quaternion.rs @@ -40,6 +40,11 @@ use simba::scalar::{ClosedNeg, RealField}; /// See #[repr(C)] #[derive(Debug, Copy, Clone)] +#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))] +#[cfg_attr( + feature = "rkyv-serialize-no-std", + derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize) +)] #[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))] pub struct DualQuaternion { /// The real component of the quaternion diff --git a/src/geometry/isometry.rs b/src/geometry/isometry.rs index 0179f1ff..92169742 100755 --- a/src/geometry/isometry.rs +++ b/src/geometry/isometry.rs @@ -66,6 +66,11 @@ use crate::geometry::{AbstractRotation, Point, Translation}; Owned>: Deserialize<'de>, T: Scalar")) )] +#[cfg_attr( + feature = "rkyv-serialize-no-std", + derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize) +)] +#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))] pub struct Isometry { /// The pure rotational part of this isometry. pub rotation: R, @@ -73,66 +78,6 @@ pub struct Isometry { pub translation: Translation, } -#[cfg(feature = "rkyv-serialize-no-std")] -mod rkyv_impl { - use super::Isometry; - use crate::{base::Scalar, geometry::Translation}; - use rkyv::{offset_of, project_struct, Archive, Deserialize, Fallible, Serialize}; - - impl Archive for Isometry - where - T::Archived: Scalar, - { - type Archived = Isometry; - type Resolver = (R::Resolver, as Archive>::Resolver); - - fn resolve( - &self, - pos: usize, - resolver: Self::Resolver, - out: &mut core::mem::MaybeUninit, - ) { - self.rotation.resolve( - pos + offset_of!(Self::Archived, rotation), - resolver.0, - project_struct!(out: Self::Archived => rotation), - ); - self.translation.resolve( - pos + offset_of!(Self::Archived, translation), - resolver.1, - project_struct!(out: Self::Archived => translation), - ); - } - } - - impl, R: Serialize, S: Fallible + ?Sized, const D: usize> - Serialize for Isometry - where - T::Archived: Scalar, - { - fn serialize(&self, serializer: &mut S) -> Result { - Ok(( - self.rotation.serialize(serializer)?, - self.translation.serialize(serializer)?, - )) - } - } - - impl - Deserialize, _D> for Isometry - where - T::Archived: Scalar + Deserialize, - R::Archived: Scalar + Deserialize, - { - fn deserialize(&self, deserializer: &mut _D) -> Result, _D::Error> { - Ok(Isometry { - rotation: self.rotation.deserialize(deserializer)?, - translation: self.translation.deserialize(deserializer)?, - }) - } - } -} - impl hash::Hash for Isometry where Owned>: hash::Hash, diff --git a/src/geometry/orthographic.rs b/src/geometry/orthographic.rs index 1119d4e3..7348f676 100644 --- a/src/geometry/orthographic.rs +++ b/src/geometry/orthographic.rs @@ -19,6 +19,11 @@ use crate::geometry::{Point3, Projective3}; /// A 3D orthographic projection stored as a homogeneous 4x4 matrix. #[repr(C)] +#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))] +#[cfg_attr( + feature = "rkyv-serialize-no-std", + derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize) +)] #[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))] #[derive(Copy, Clone)] pub struct Orthographic3 { diff --git a/src/geometry/perspective.rs b/src/geometry/perspective.rs index 8ebab3e4..351960bb 100644 --- a/src/geometry/perspective.rs +++ b/src/geometry/perspective.rs @@ -20,6 +20,11 @@ use crate::geometry::{Point3, Projective3}; /// A 3D perspective projection stored as a homogeneous 4x4 matrix. #[repr(C)] +#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))] +#[cfg_attr( + feature = "rkyv-serialize-no-std", + derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize) +)] #[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))] #[derive(Copy, Clone)] pub struct Perspective3 { diff --git a/src/geometry/point.rs b/src/geometry/point.rs index cdc590fa..306c18e5 100644 --- a/src/geometry/point.rs +++ b/src/geometry/point.rs @@ -36,6 +36,11 @@ use std::mem::MaybeUninit; /// of said transformations for details. #[repr(C)] #[derive(Clone)] +#[cfg_attr( + feature = "rkyv-serialize-no-std", + derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize) +)] +#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))] pub struct OPoint where DefaultAllocator: Allocator, diff --git a/src/geometry/quaternion.rs b/src/geometry/quaternion.rs index 987c9757..f38dca6f 100755 --- a/src/geometry/quaternion.rs +++ b/src/geometry/quaternion.rs @@ -23,6 +23,11 @@ use crate::geometry::{Point3, Rotation}; /// that may be used as a rotation. #[repr(C)] #[derive(Copy, Clone)] +#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))] +#[cfg_attr( + feature = "rkyv-serialize-no-std", + derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize) +)] #[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))] pub struct Quaternion { /// This quaternion as a 4D vector of coordinates in the `[ x, y, z, w ]` storage order. @@ -97,48 +102,6 @@ where } } -#[cfg(feature = "rkyv-serialize-no-std")] -mod rkyv_impl { - use super::Quaternion; - use crate::base::Vector4; - use rkyv::{offset_of, project_struct, Archive, Deserialize, Fallible, Serialize}; - - impl Archive for Quaternion { - type Archived = Quaternion; - type Resolver = as Archive>::Resolver; - - fn resolve( - &self, - pos: usize, - resolver: Self::Resolver, - out: &mut core::mem::MaybeUninit, - ) { - self.coords.resolve( - pos + offset_of!(Self::Archived, coords), - resolver, - project_struct!(out: Self::Archived => coords), - ); - } - } - - impl, S: Fallible + ?Sized> Serialize for Quaternion { - fn serialize(&self, serializer: &mut S) -> Result { - self.coords.serialize(serializer) - } - } - - impl Deserialize, D> for Quaternion - where - T::Archived: Deserialize, - { - fn deserialize(&self, deserializer: &mut D) -> Result, D::Error> { - Ok(Quaternion { - coords: self.coords.deserialize(deserializer)?, - }) - } - } -} - impl Quaternion where T::Element: SimdRealField, diff --git a/src/geometry/rotation.rs b/src/geometry/rotation.rs index 4dbcfb43..2a8bf427 100755 --- a/src/geometry/rotation.rs +++ b/src/geometry/rotation.rs @@ -49,6 +49,11 @@ use crate::geometry::Point; /// * [Conversion to a matrix `matrix`, `to_homogeneous`…](#conversion-to-a-matrix) /// #[repr(C)] +#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))] +#[cfg_attr( + feature = "rkyv-serialize-no-std", + derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize) +)] #[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))] #[derive(Copy, Clone)] pub struct Rotation { diff --git a/src/geometry/scale.rs b/src/geometry/scale.rs index abaeeccc..37da1ef0 100755 --- a/src/geometry/scale.rs +++ b/src/geometry/scale.rs @@ -17,6 +17,11 @@ use crate::geometry::Point; /// A scale which supports non-uniform scaling. #[repr(C)] +#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))] +#[cfg_attr( + feature = "rkyv-serialize-no-std", + derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize) +)] #[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))] #[derive(Copy, Clone)] pub struct Scale { @@ -84,49 +89,6 @@ where } } -#[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`. /// diff --git a/src/geometry/similarity.rs b/src/geometry/similarity.rs index 9658685e..8c38ff1e 100755 --- a/src/geometry/similarity.rs +++ b/src/geometry/similarity.rs @@ -34,6 +34,11 @@ use crate::geometry::{AbstractRotation, Isometry, Point, Translation}; DefaultAllocator: Allocator>, Owned>: Deserialize<'de>")) )] +#[cfg_attr( + feature = "rkyv-serialize-no-std", + derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize) +)] +#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))] pub struct Similarity { /// The part of this similarity that does not include the scaling factor. pub isometry: Isometry, diff --git a/src/geometry/translation.rs b/src/geometry/translation.rs index e1921d0a..bef66f68 100755 --- a/src/geometry/translation.rs +++ b/src/geometry/translation.rs @@ -17,6 +17,11 @@ use crate::geometry::Point; /// A translation. #[repr(C)] +#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))] +#[cfg_attr( + feature = "rkyv-serialize-no-std", + derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize) +)] #[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))] #[derive(Copy, Clone)] pub struct Translation { @@ -84,49 +89,6 @@ where } } -#[cfg(feature = "rkyv-serialize-no-std")] -mod rkyv_impl { - use super::Translation; - use crate::base::SVector; - use rkyv::{offset_of, project_struct, Archive, Deserialize, Fallible, Serialize}; - - impl Archive for Translation { - type Archived = Translation; - 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 Translation { - fn serialize(&self, serializer: &mut S) -> Result { - self.vector.serialize(serializer) - } - } - - impl Deserialize, _D> - for Translation - where - T::Archived: Deserialize, - { - fn deserialize(&self, deserializer: &mut _D) -> Result, _D::Error> { - Ok(Translation { - vector: self.vector.deserialize(deserializer)?, - }) - } - } -} - impl Translation { /// Creates a new translation from the given vector. #[inline] diff --git a/src/geometry/unit_complex.rs b/src/geometry/unit_complex.rs index caf25493..efe0dac2 100755 --- a/src/geometry/unit_complex.rs +++ b/src/geometry/unit_complex.rs @@ -410,7 +410,8 @@ where #[inline] #[must_use] pub fn slerp(&self, other: &Self, t: T) -> Self { - Self::new(self.angle() * (T::one() - t.clone()) + other.angle() * t) + let delta = other / self; + self * Self::new(delta.angle() * t) } } diff --git a/src/lib.rs b/src/lib.rs index 92b28dcb..1ee1a3ba 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -78,12 +78,13 @@ an optimized set of tools for computer graphics and physics. Those features incl unused_mut, unused_parens, unused_qualifications, - unused_results, rust_2018_idioms, rust_2018_compatibility, future_incompatible, missing_copy_implementations )] +#![cfg_attr(feature = "rkyv-serialize-no-std", warn(unused_results))] // TODO: deny this once bytecheck stops generating warnings. +#![cfg_attr(not(feature = "rkyv-serialize-no-std"), deny(unused_results))] #![doc( html_favicon_url = "https://nalgebra.org/img/favicon.ico", html_root_url = "https://docs.rs/nalgebra/0.25.0" diff --git a/src/third_party/glam/common/glam_isometry.rs b/src/third_party/glam/common/glam_isometry.rs index 3a8d4961..7b188c39 100644 --- a/src/third_party/glam/common/glam_isometry.rs +++ b/src/third_party/glam/common/glam_isometry.rs @@ -1,5 +1,5 @@ use super::glam::{DMat3, DMat4, DQuat, DVec2, DVec3, Mat3, Mat4, Quat, Vec2, Vec3}; -use crate::{Isometry2, Isometry3, Matrix3, Matrix4, Translation3, UnitQuaternion, Vector2}; +use crate::{Isometry2, Isometry3, Matrix3, Matrix4, Translation3, UnitQuaternion}; use std::convert::TryFrom; impl From> for Mat3 { @@ -36,18 +36,18 @@ impl From> for (DVec3, DQuat) { } } -impl From> for (Vec3, Quat) { - fn from(iso: Isometry2) -> (Vec3, Quat) { - let tra = Vec3::new(iso.translation.x, iso.translation.y, 0.0); - let rot = Quat::from_axis_angle(Vec3::Z, iso.rotation.angle()); +impl From> for (Vec2, f32) { + fn from(iso: Isometry2) -> (Vec2, f32) { + let tra = Vec2::new(iso.translation.x, iso.translation.y); + let rot = iso.rotation.angle(); (tra, rot) } } -impl From> for (DVec3, DQuat) { - fn from(iso: Isometry2) -> (DVec3, DQuat) { - let tra = DVec3::new(iso.translation.x, iso.translation.y, 0.0); - let rot = DQuat::from_axis_angle(DVec3::Z, iso.rotation.angle()); +impl From> for (DVec2, f64) { + fn from(iso: Isometry2) -> (DVec2, f64) { + let tra = DVec2::new(iso.translation.x, iso.translation.y); + let rot = iso.rotation.angle(); (tra, rot) } } @@ -64,30 +64,6 @@ impl From<(DVec3, DQuat)> for Isometry3 { } } -impl From<(Vec3, Quat)> for Isometry2 { - fn from((tra, rot): (Vec3, Quat)) -> Self { - Isometry2::new([tra.x, tra.y].into(), rot.to_axis_angle().1) - } -} - -impl From<(DVec3, DQuat)> for Isometry2 { - fn from((tra, rot): (DVec3, DQuat)) -> Self { - Isometry2::new([tra.x, tra.y].into(), rot.to_axis_angle().1) - } -} - -impl From<(Vec2, Quat)> for Isometry2 { - fn from((tra, rot): (Vec2, Quat)) -> Self { - Isometry2::new(tra.into(), rot.to_axis_angle().1) - } -} - -impl From<(DVec2, DQuat)> for Isometry2 { - fn from((tra, rot): (DVec2, DQuat)) -> Self { - Isometry2::new(tra.into(), rot.to_axis_angle().1) - } -} - impl From<(Vec2, f32)> for Isometry2 { fn from((tra, rot): (Vec2, f32)) -> Self { Isometry2::new(tra.into(), rot) @@ -112,18 +88,6 @@ impl From for Isometry3 { } } -impl From for Isometry2 { - fn from(rot: Quat) -> Self { - Isometry2::new(Vector2::zeros(), rot.to_axis_angle().1) - } -} - -impl From for Isometry2 { - fn from(rot: DQuat) -> Self { - Isometry2::new(Vector2::zeros(), rot.to_axis_angle().1) - } -} - impl From for Isometry3 { fn from(tra: Vec3) -> Self { Isometry3::from_parts(tra.into(), UnitQuaternion::identity()) @@ -148,18 +112,6 @@ impl From for Isometry2 { } } -impl From for Isometry2 { - fn from(tra: Vec3) -> Self { - Isometry2::new([tra.x, tra.y].into(), 0.0) - } -} - -impl From for Isometry2 { - fn from(tra: DVec3) -> Self { - Isometry2::new([tra.x, tra.y].into(), 0.0) - } -} - impl TryFrom for Isometry2 { type Error = (); diff --git a/src/third_party/glam/common/glam_matrix.rs b/src/third_party/glam/common/glam_matrix.rs index 80f88054..fa9f713f 100644 --- a/src/third_party/glam/common/glam_matrix.rs +++ b/src/third_party/glam/common/glam_matrix.rs @@ -3,7 +3,11 @@ use super::glam::{ Mat4, UVec2, UVec3, UVec4, Vec2, Vec3, Vec3A, Vec4, }; use crate::storage::RawStorage; -use crate::{Matrix, Matrix2, Matrix3, Matrix4, Vector, Vector2, Vector3, Vector4, U2, U3, U4}; +use crate::{ + Matrix, Matrix2, Matrix3, Matrix4, Unit, UnitVector2, UnitVector3, UnitVector4, Vector, + Vector2, Vector3, Vector4, U2, U3, U4, +}; +use std::convert::TryFrom; macro_rules! impl_vec_conversion( ($N: ty, $Vec2: ty, $Vec3: ty, $Vec4: ty) => { @@ -66,6 +70,63 @@ impl_vec_conversion!(i32, IVec2, IVec3, IVec4); impl_vec_conversion!(u32, UVec2, UVec3, UVec4); impl_vec_conversion!(bool, BVec2, BVec3, BVec4); +const ERR: &'static str = "Normalization failed."; + +macro_rules! impl_unit_vec_conversion( + ($N: ty, $Vec2: ty, $Vec3: ty, $Vec4: ty) => { + impl TryFrom<$Vec2> for UnitVector2<$N> { + type Error = &'static str; + #[inline] + fn try_from(e: $Vec2) -> Result { + Unit::try_new(e.into(), 0.0).ok_or(ERR) + } + } + + impl From> for $Vec2 + { + #[inline] + fn from(e: UnitVector2<$N>) -> $Vec2 { + e.into_inner().into() + } + } + + impl TryFrom<$Vec3> for UnitVector3<$N> { + type Error = &'static str; + #[inline] + fn try_from(e: $Vec3) -> Result { + Unit::try_new(e.into(), 0.0).ok_or(ERR) + } + } + + impl From> for $Vec3 + { + #[inline] + fn from(e: UnitVector3<$N>) -> $Vec3 { + e.into_inner().into() + } + } + + impl TryFrom<$Vec4> for UnitVector4<$N> { + type Error = &'static str; + #[inline] + fn try_from(e: $Vec4) -> Result { + Unit::try_new(e.into(), 0.0).ok_or(ERR) + } + } + + impl From> for $Vec4 + { + #[inline] + fn from(e: UnitVector4<$N>) -> $Vec4 { + e.into_inner().into() + } + } + } +); + +impl_unit_vec_conversion!(f32, Vec2, Vec3, Vec4); +impl_unit_vec_conversion!(f64, DVec2, DVec3, DVec4); + impl From for Vector3 { #[inline] fn from(e: Vec3A) -> Vector3 { @@ -83,6 +144,21 @@ where } } +impl TryFrom for UnitVector3 { + type Error = &'static str; + #[inline] + fn try_from(e: Vec3A) -> Result { + Unit::try_new(e.into(), 0.0).ok_or(ERR) + } +} + +impl From> for Vec3A { + #[inline] + fn from(e: UnitVector3) -> Vec3A { + e.into_inner().into() + } +} + impl From for Matrix2 { #[inline] fn from(e: Mat2) -> Matrix2 { diff --git a/src/third_party/glam/mod.rs b/src/third_party/glam/mod.rs index d24ff7e5..ae2c4514 100644 --- a/src/third_party/glam/mod.rs +++ b/src/third_party/glam/mod.rs @@ -1,5 +1,3 @@ -#[cfg(feature = "glam013")] -mod v013; #[cfg(feature = "glam014")] mod v014; #[cfg(feature = "glam015")] diff --git a/src/third_party/glam/v013/mod.rs b/src/third_party/glam/v013/mod.rs deleted file mode 100644 index 4787fb21..00000000 --- a/src/third_party/glam/v013/mod.rs +++ /dev/null @@ -1,18 +0,0 @@ -#[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 glam013 as glam; diff --git a/tests/geometry/rotation.rs b/tests/geometry/rotation.rs index 9a29772e..84bba676 100644 --- a/tests/geometry/rotation.rs +++ b/tests/geometry/rotation.rs @@ -32,7 +32,9 @@ fn quaternion_euler_angles_issue_494() { #[cfg(feature = "proptest-support")] mod proptest_tests { + use approx::AbsDiffEq; use na::{self, Rotation2, Rotation3, Unit}; + use na::{UnitComplex, UnitQuaternion}; use simba::scalar::RealField; use std::f64; @@ -229,5 +231,74 @@ mod proptest_tests { prop_assert_eq!(r, Rotation3::identity()) } } + + // + //In general, `slerp(a,b,t)` should equal `(b/a)^t * a` even though in practice, + //we may not use that formula directly for complex numbers or quaternions + // + + #[test] + fn slerp_powf_agree_2(a in unit_complex(), b in unit_complex(), t in PROPTEST_F64) { + let z1 = a.slerp(&b, t); + let z2 = (b/a).powf(t) * a; + prop_assert!(relative_eq!(z1,z2,epsilon=1e-10)); + } + + #[test] + fn slerp_powf_agree_3(a in unit_quaternion(), b in unit_quaternion(), t in PROPTEST_F64) { + if let Some(z1) = a.try_slerp(&b, t, f64::default_epsilon()) { + let z2 = (b/a).powf(t) * a; + prop_assert!(relative_eq!(z1,z2,epsilon=1e-10)); + } + } + + // + //when not antipodal, slerp should always take the shortest path between two orientations + // + + #[test] + fn slerp_takes_shortest_path_2( + z in unit_complex(), dtheta in -f64::pi()..f64::pi(), t in 0.0..1.0f64 + ) { + + //ambiguous when at ends of angle range, so we don't really care here + if dtheta.abs() != f64::pi() { + + //make two complex numbers separated by an angle between -pi and pi + let (z1, z2) = (z, z * UnitComplex::new(dtheta)); + let z3 = z1.slerp(&z2, t); + + //since the angle is no larger than a half-turn, and t is between 0 and 1, + //the shortest path just corresponds to adding the scaled angle + let a1 = z3.angle(); + let a2 = na::wrap(z1.angle() + dtheta*t, -f64::pi(), f64::pi()); + + prop_assert!(relative_eq!(a1, a2, epsilon=1e-10)); + } + + } + + #[test] + fn slerp_takes_shortest_path_3( + q in unit_quaternion(), dtheta in -f64::pi()..f64::pi(), t in 0.0..1.0f64 + ) { + + //ambiguous when at ends of angle range, so we don't really care here + if let Some(axis) = q.axis() { + + //make two quaternions separated by an angle between -pi and pi + let (q1, q2) = (q, q * UnitQuaternion::from_axis_angle(&axis, dtheta)); + let q3 = q1.slerp(&q2, t); + + //since the angle is no larger than a half-turn, and t is between 0 and 1, + //the shortest path just corresponds to adding the scaled angle + let q4 = q1 * UnitQuaternion::from_axis_angle(&axis, dtheta*t); + prop_assert!(relative_eq!(q3, q4, epsilon=1e-10)); + + } + + } + + } }