diff --git a/CHANGELOG.md b/CHANGELOG.md index 7487b236..57f18f25 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,15 @@ documented here. This project adheres to [Semantic Versioning](https://semver.org/). +## [0.22.0] - WIP + +### Added + * `Cholesky::new_unchecked` which build a Cholesky decomposition without checking that its input is + positive-definite. It can be use with SIMD types. + * The `Default` trait is now implemented for matrices, and quaternions. They are all filled with zeros, + except for `UnitQuaternion` which is initialized with the identity. + * Matrix exponential `matrix.exp()`. + ## [0.21.0] In this release, we are no longer relying on traits from the __alga__ crate for our generic code. diff --git a/Cargo.toml b/Cargo.toml index 21a365e3..4adf2ab1 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "nalgebra" -version = "0.21.0" +version = "0.21.1" authors = [ "Sébastien Crozet " ] description = "Linear algebra library with transformations and statically-sized or dynamically-sized matrices." diff --git a/nalgebra-glm/src/constructors.rs b/nalgebra-glm/src/constructors.rs index 949ea9e4..175626a1 100644 --- a/nalgebra-glm/src/constructors.rs +++ b/nalgebra-glm/src/constructors.rs @@ -1,9 +1,8 @@ -#![cfg_attr(rustfmt, rustfmt_skip)] - -use na::{Scalar, RealField, U2, U3, U4}; -use crate::aliases::{TMat, Qua, TVec1, TVec2, TVec3, TVec4, TMat2, TMat2x3, TMat2x4, TMat3, TMat3x2, TMat3x4, - TMat4, TMat4x2, TMat4x3}; - +use crate::aliases::{ + Qua, TMat, TMat2, TMat2x3, TMat2x4, TMat3, TMat3x2, TMat3x4, TMat4, TMat4x2, TMat4x3, TVec1, + TVec2, TVec3, TVec4, +}; +use na::{RealField, Scalar, U2, U3, U4}; /// Creates a new 1D vector. /// @@ -34,8 +33,8 @@ pub fn vec4(x: N, y: N, z: N, w: N) -> TVec4 { TVec4::new(x, y, z, w) } - /// Create a new 2x2 matrix. +#[rustfmt::skip] pub fn mat2(m11: N, m12: N, m21: N, m22: N) -> TMat2 { TMat::::new( @@ -45,6 +44,7 @@ pub fn mat2(m11: N, m12: N, } /// Create a new 2x2 matrix. +#[rustfmt::skip] pub fn mat2x2(m11: N, m12: N, m21: N, m22: N) -> TMat2 { TMat::::new( @@ -54,6 +54,7 @@ pub fn mat2x2(m11: N, m12: N, } /// Create a new 2x3 matrix. +#[rustfmt::skip] pub fn mat2x3(m11: N, m12: N, m13: N, m21: N, m22: N, m23: N) -> TMat2x3 { TMat::::new( @@ -63,6 +64,7 @@ pub fn mat2x3(m11: N, m12: N, m13: N, } /// Create a new 2x4 matrix. +#[rustfmt::skip] pub fn mat2x4(m11: N, m12: N, m13: N, m14: N, m21: N, m22: N, m23: N, m24: N) -> TMat2x4 { TMat::::new( @@ -72,6 +74,7 @@ pub fn mat2x4(m11: N, m12: N, m13: N, m14: N, } /// Create a new 3x3 matrix. +#[rustfmt::skip] pub fn mat3(m11: N, m12: N, m13: N, m21: N, m22: N, m23: N, m31: N, m32: N, m33: N) -> TMat3 { @@ -83,6 +86,7 @@ pub fn mat3(m11: N, m12: N, m13: N, } /// Create a new 3x2 matrix. +#[rustfmt::skip] pub fn mat3x2(m11: N, m12: N, m21: N, m22: N, m31: N, m32: N) -> TMat3x2 { @@ -94,6 +98,7 @@ pub fn mat3x2(m11: N, m12: N, } /// Create a new 3x3 matrix. +#[rustfmt::skip] pub fn mat3x3(m11: N, m12: N, m13: N, m21: N, m22: N, m23: N, m31: N, m32: N, m33: N) -> TMat3 { @@ -105,6 +110,7 @@ pub fn mat3x3(m11: N, m12: N, m13: N, } /// Create a new 3x4 matrix. +#[rustfmt::skip] pub fn mat3x4(m11: N, m12: N, m13: N, m14: N, m21: N, m22: N, m23: N, m24: N, m31: N, m32: N, m33: N, m34: N) -> TMat3x4 { @@ -116,6 +122,7 @@ pub fn mat3x4(m11: N, m12: N, m13: N, m14: N, } /// Create a new 4x2 matrix. +#[rustfmt::skip] pub fn mat4x2(m11: N, m12: N, m21: N, m22: N, m31: N, m32: N, @@ -129,6 +136,7 @@ pub fn mat4x2(m11: N, m12: N, } /// Create a new 4x3 matrix. +#[rustfmt::skip] pub fn mat4x3(m11: N, m12: N, m13: N, m21: N, m22: N, m23: N, m31: N, m32: N, m33: N, @@ -142,6 +150,7 @@ pub fn mat4x3(m11: N, m12: N, m13: N, } /// Create a new 4x4 matrix. +#[rustfmt::skip] pub fn mat4x4(m11: N, m12: N, m13: N, m14: N, m21: N, m22: N, m23: N, m24: N, m31: N, m32: N, m33: N, m34: N, @@ -155,6 +164,7 @@ pub fn mat4x4(m11: N, m12: N, m13: N, m14: N, } /// Create a new 4x4 matrix. +#[rustfmt::skip] pub fn mat4(m11: N, m12: N, m13: N, m14: N, m21: N, m22: N, m23: N, m24: N, m31: N, m32: N, m33: N, m34: N, diff --git a/src/base/array_storage.rs b/src/base/array_storage.rs index 3ff06e89..1743f06b 100644 --- a/src/base/array_storage.rs +++ b/src/base/array_storage.rs @@ -48,6 +48,21 @@ where /// Renamed to [ArrayStorage]. pub type MatrixArray = ArrayStorage; +impl Default for ArrayStorage +where + R: DimName, + C: DimName, + R::Value: Mul, + Prod: ArrayLength, + N: Default, +{ + fn default() -> Self { + ArrayStorage { + data: Default::default(), + } + } +} + impl Hash for ArrayStorage where N: Hash, diff --git a/src/base/construction_slice.rs b/src/base/construction_slice.rs index d63b374b..06b8be9f 100644 --- a/src/base/construction_slice.rs +++ b/src/base/construction_slice.rs @@ -220,9 +220,9 @@ macro_rules! impl_constructors( // FIXME: this is not very pretty. We could find a better call syntax. impl_constructors!(R, C; // Arguments for Matrix - => R: DimName, => C: DimName; // Type parameters for impl - R::name(), C::name(); // Arguments for `_generic` constructors. - ); // Arguments for non-generic constructors. +=> R: DimName, => C: DimName; // Type parameters for impl +R::name(), C::name(); // Arguments for `_generic` constructors. +); // Arguments for non-generic constructors. impl_constructors!(R, Dynamic; => R: DimName; @@ -279,9 +279,9 @@ macro_rules! impl_constructors_mut( // FIXME: this is not very pretty. We could find a better call syntax. impl_constructors_mut!(R, C; // Arguments for Matrix - => R: DimName, => C: DimName; // Type parameters for impl - R::name(), C::name(); // Arguments for `_generic` constructors. - ); // Arguments for non-generic constructors. +=> R: DimName, => C: DimName; // Type parameters for impl +R::name(), C::name(); // Arguments for `_generic` constructors. +); // Arguments for non-generic constructors. impl_constructors_mut!(R, Dynamic; => R: DimName; diff --git a/src/base/matrix.rs b/src/base/matrix.rs index 783f8c19..e4821bf8 100644 --- a/src/base/matrix.rs +++ b/src/base/matrix.rs @@ -94,6 +94,21 @@ impl fmt::Debug for Matrix } } +impl Default for Matrix +where + N: Scalar, + R: Dim, + C: Dim, + S: Default, +{ + fn default() -> Self { + Matrix { + data: Default::default(), + _phantoms: PhantomData, + } + } +} + #[cfg(feature = "serde-serialize")] impl Serialize for Matrix where diff --git a/src/geometry/quaternion.rs b/src/geometry/quaternion.rs index 45a03c78..0c996ede 100755 --- a/src/geometry/quaternion.rs +++ b/src/geometry/quaternion.rs @@ -33,6 +33,14 @@ pub struct Quaternion { pub coords: Vector4, } +impl Default for Quaternion { + fn default() -> Self { + Quaternion { + coords: Vector4::zeros(), + } + } +} + #[cfg(feature = "abomonation-serialize")] impl Abomonation for Quaternion where @@ -1536,6 +1544,12 @@ where } } +impl Default for UnitQuaternion { + fn default() -> Self { + Self::identity() + } +} + impl fmt::Display for UnitQuaternion { fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { if let Some(axis) = self.axis() { diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index e02209b6..d9f33f27 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -3,6 +3,7 @@ use serde::{Deserialize, Serialize}; use num::One; use simba::scalar::ComplexField; +use simba::simd::SimdComplexField; use crate::allocator::Allocator; use crate::base::{DefaultAllocator, Matrix, MatrixMN, MatrixN, SquareMatrix, Vector}; @@ -23,20 +24,122 @@ use crate::storage::{Storage, StorageMut}; MatrixN: Deserialize<'de>")) )] #[derive(Clone, Debug)] -pub struct Cholesky +pub struct Cholesky where DefaultAllocator: Allocator, { chol: MatrixN, } -impl Copy for Cholesky +impl Copy for Cholesky where DefaultAllocator: Allocator, MatrixN: Copy, { } +impl Cholesky +where + DefaultAllocator: Allocator, +{ + /// Computes the Cholesky decomposition of `matrix` without checking that the matrix is definite-positive. + /// + /// If the input matrix is not definite-positive, the decomposition may contain trash values (Inf, NaN, etc.) + pub fn new_unchecked(mut matrix: MatrixN) -> Self { + assert!(matrix.is_square(), "The input matrix must be square."); + + let n = matrix.nrows(); + + for j in 0..n { + for k in 0..j { + let factor = unsafe { -*matrix.get_unchecked((j, k)) }; + + let (mut col_j, col_k) = matrix.columns_range_pair_mut(j, k); + let mut col_j = col_j.rows_range_mut(j..); + let col_k = col_k.rows_range(j..); + col_j.axpy(factor.simd_conjugate(), &col_k, N::one()); + } + + let diag = unsafe { *matrix.get_unchecked((j, j)) }; + let denom = diag.simd_sqrt(); + + unsafe { + *matrix.get_unchecked_mut((j, j)) = denom; + } + + let mut col = matrix.slice_range_mut(j + 1.., j); + col /= denom; + } + + 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) -> MatrixN { + self.chol.fill_upper_triangle(N::zero(), 1); + self.chol + } + + /// Retrieves the lower-triangular factor of the Cholesky decomposition, without zeroing-out + /// its strict upper-triangular part. + /// + /// The values of the strict upper-triangular part are garbage and should be ignored by further + /// computations. + pub fn unpack_dirty(self) -> MatrixN { + self.chol + } + + /// Retrieves the lower-triangular factor of the Cholesky decomposition with its strictly + /// uppen-triangular part filled with zeros. + pub fn l(&self) -> MatrixN { + self.chol.lower_triangle() + } + + /// Retrieves the lower-triangular factor of the Cholesky decomposition, without zeroing-out + /// its strict upper-triangular part. + /// + /// This is an allocation-less version of `self.l()`. The values of the strict upper-triangular + /// part are garbage and should be ignored by further computations. + pub fn l_dirty(&self) -> &MatrixN { + &self.chol + } + + /// Solves the system `self * x = b` where `self` is the decomposed matrix and `x` the unknown. + /// + /// The result is stored on `b`. + pub fn solve_mut(&self, b: &mut Matrix) + where + S2: StorageMut, + ShapeConstraint: SameNumberOfRows, + { + self.chol.solve_lower_triangular_unchecked_mut(b); + self.chol.ad_solve_lower_triangular_unchecked_mut(b); + } + + /// Returns the solution of the system `self * x = b` where `self` is the decomposed matrix and + /// `x` the unknown. + pub fn solve(&self, b: &Matrix) -> MatrixMN + where + S2: Storage, + DefaultAllocator: Allocator, + ShapeConstraint: SameNumberOfRows, + { + let mut res = b.clone_owned(); + self.solve_mut(&mut res); + res + } + + /// Computes the inverse of the decomposed matrix. + pub fn inverse(&self) -> MatrixN { + let shape = self.chol.data.shape(); + let mut res = MatrixN::identity_generic(shape.0, shape.1); + + self.solve_mut(&mut res); + res + } +} + impl Cholesky where DefaultAllocator: Allocator, @@ -82,71 +185,6 @@ where Some(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) -> MatrixN { - self.chol.fill_upper_triangle(N::zero(), 1); - self.chol - } - - /// Retrieves the lower-triangular factor of the Cholesky decomposition, without zeroing-out - /// its strict upper-triangular part. - /// - /// The values of the strict upper-triangular part are garbage and should be ignored by further - /// computations. - pub fn unpack_dirty(self) -> MatrixN { - self.chol - } - - /// Retrieves the lower-triangular factor of the Cholesky decomposition with its strictly - /// uppen-triangular part filled with zeros. - pub fn l(&self) -> MatrixN { - self.chol.lower_triangle() - } - - /// Retrieves the lower-triangular factor of the Cholesky decomposition, without zeroing-out - /// its strict upper-triangular part. - /// - /// This is an allocation-less version of `self.l()`. The values of the strict upper-triangular - /// part are garbage and should be ignored by further computations. - pub fn l_dirty(&self) -> &MatrixN { - &self.chol - } - - /// Solves the system `self * x = b` where `self` is the decomposed matrix and `x` the unknown. - /// - /// The result is stored on `b`. - pub fn solve_mut(&self, b: &mut Matrix) - where - S2: StorageMut, - ShapeConstraint: SameNumberOfRows, - { - let _ = self.chol.solve_lower_triangular_mut(b); - let _ = self.chol.ad_solve_lower_triangular_mut(b); - } - - /// Returns the solution of the system `self * x = b` where `self` is the decomposed matrix and - /// `x` the unknown. - pub fn solve(&self, b: &Matrix) -> MatrixMN - where - S2: Storage, - DefaultAllocator: Allocator, - ShapeConstraint: SameNumberOfRows, - { - let mut res = b.clone_owned(); - self.solve_mut(&mut res); - res - } - - /// Computes the inverse of the decomposed matrix. - pub fn inverse(&self) -> MatrixN { - let shape = self.chol.data.shape(); - let mut res = MatrixN::identity_generic(shape.0, shape.1); - - self.solve_mut(&mut res); - res - } - /// Given the Cholesky decomposition of a matrix `M`, a scalar `sigma` and a vector `v`, /// performs a rank one update such that we end up with the decomposition of `M + sigma * (v * v.adjoint())`. #[inline] diff --git a/src/linalg/exp.rs b/src/linalg/exp.rs new file mode 100644 index 00000000..c11fd757 --- /dev/null +++ b/src/linalg/exp.rs @@ -0,0 +1,492 @@ +//! This module provides the matrix exponent (exp) function to square matrices. +//! +use crate::{ + base::{ + allocator::Allocator, + dimension::{Dim, DimMin, DimMinimum, U1}, + storage::Storage, + DefaultAllocator, + }, + convert, try_convert, ComplexField, MatrixN, RealField, +}; + +// https://github.com/scipy/scipy/blob/c1372d8aa90a73d8a52f135529293ff4edb98fc8/scipy/sparse/linalg/matfuncs.py +struct ExpmPadeHelper +where + N: RealField, + D: DimMin, + DefaultAllocator: Allocator + Allocator<(usize, usize), DimMinimum>, +{ + use_exact_norm: bool, + ident: MatrixN, + + a: MatrixN, + a2: Option>, + a4: Option>, + a6: Option>, + a8: Option>, + a10: Option>, + + d4_exact: Option, + d6_exact: Option, + d8_exact: Option, + d10_exact: Option, + + d4_approx: Option, + d6_approx: Option, + d8_approx: Option, + d10_approx: Option, +} + +impl ExpmPadeHelper +where + N: RealField, + D: DimMin, + DefaultAllocator: Allocator + Allocator<(usize, usize), DimMinimum>, +{ + fn new(a: MatrixN, use_exact_norm: bool) -> Self { + let (nrows, ncols) = a.data.shape(); + ExpmPadeHelper { + use_exact_norm, + ident: MatrixN::::identity_generic(nrows, ncols), + a, + a2: None, + a4: None, + a6: None, + a8: None, + a10: None, + d4_exact: None, + d6_exact: None, + d8_exact: None, + d10_exact: None, + d4_approx: None, + d6_approx: None, + d8_approx: None, + d10_approx: None, + } + } + + fn calc_a2(&mut self) { + if self.a2.is_none() { + self.a2 = Some(&self.a * &self.a); + } + } + + fn calc_a4(&mut self) { + if self.a4.is_none() { + self.calc_a2(); + let a2 = self.a2.as_ref().unwrap(); + self.a4 = Some(a2 * a2); + } + } + + fn calc_a6(&mut self) { + if self.a6.is_none() { + self.calc_a2(); + self.calc_a4(); + let a2 = self.a2.as_ref().unwrap(); + let a4 = self.a4.as_ref().unwrap(); + self.a6 = Some(a4 * a2); + } + } + + fn calc_a8(&mut self) { + if self.a8.is_none() { + self.calc_a2(); + self.calc_a6(); + let a2 = self.a2.as_ref().unwrap(); + let a6 = self.a6.as_ref().unwrap(); + self.a8 = Some(a6 * a2); + } + } + + fn calc_a10(&mut self) { + if self.a10.is_none() { + self.calc_a4(); + self.calc_a6(); + let a4 = self.a4.as_ref().unwrap(); + let a6 = self.a6.as_ref().unwrap(); + self.a10 = Some(a6 * a4); + } + } + + fn d4_tight(&mut self) -> N { + if self.d4_exact.is_none() { + self.calc_a4(); + self.d4_exact = Some(one_norm(self.a4.as_ref().unwrap()).powf(convert(0.25))); + } + self.d4_exact.unwrap() + } + + fn d6_tight(&mut self) -> N { + if self.d6_exact.is_none() { + self.calc_a6(); + self.d6_exact = Some(one_norm(self.a6.as_ref().unwrap()).powf(convert(1.0 / 6.0))); + } + self.d6_exact.unwrap() + } + + fn d8_tight(&mut self) -> N { + if self.d8_exact.is_none() { + self.calc_a8(); + self.d8_exact = Some(one_norm(self.a8.as_ref().unwrap()).powf(convert(1.0 / 8.0))); + } + self.d8_exact.unwrap() + } + + fn d10_tight(&mut self) -> N { + if self.d10_exact.is_none() { + self.calc_a10(); + self.d10_exact = Some(one_norm(self.a10.as_ref().unwrap()).powf(convert(1.0 / 10.0))); + } + self.d10_exact.unwrap() + } + + fn d4_loose(&mut self) -> N { + if self.use_exact_norm { + return self.d4_tight(); + } + + if self.d4_exact.is_some() { + return self.d4_exact.unwrap(); + } + + if self.d4_approx.is_none() { + self.calc_a4(); + self.d4_approx = Some(one_norm(self.a4.as_ref().unwrap()).powf(convert(0.25))); + } + + self.d4_approx.unwrap() + } + + fn d6_loose(&mut self) -> N { + if self.use_exact_norm { + return self.d6_tight(); + } + + if self.d6_exact.is_some() { + return self.d6_exact.unwrap(); + } + + if self.d6_approx.is_none() { + self.calc_a6(); + self.d6_approx = Some(one_norm(self.a6.as_ref().unwrap()).powf(convert(1.0 / 6.0))); + } + + self.d6_approx.unwrap() + } + + fn d8_loose(&mut self) -> N { + if self.use_exact_norm { + return self.d8_tight(); + } + + if self.d8_exact.is_some() { + return self.d8_exact.unwrap(); + } + + if self.d8_approx.is_none() { + self.calc_a8(); + self.d8_approx = Some(one_norm(self.a8.as_ref().unwrap()).powf(convert(1.0 / 8.0))); + } + + self.d8_approx.unwrap() + } + + fn d10_loose(&mut self) -> N { + if self.use_exact_norm { + return self.d10_tight(); + } + + if self.d10_exact.is_some() { + return self.d10_exact.unwrap(); + } + + if self.d10_approx.is_none() { + self.calc_a10(); + self.d10_approx = Some(one_norm(self.a10.as_ref().unwrap()).powf(convert(1.0 / 10.0))); + } + + self.d10_approx.unwrap() + } + + fn pade3(&mut self) -> (MatrixN, MatrixN) { + let b: [N; 4] = [convert(120.0), convert(60.0), convert(12.0), convert(1.0)]; + self.calc_a2(); + let a2 = self.a2.as_ref().unwrap(); + let u = &self.a * (a2 * b[3] + &self.ident * b[1]); + let v = a2 * b[2] + &self.ident * b[0]; + (u, v) + } + + fn pade5(&mut self) -> (MatrixN, MatrixN) { + let b: [N; 6] = [ + convert(30240.0), + convert(15120.0), + convert(3360.0), + convert(420.0), + convert(30.0), + convert(1.0), + ]; + self.calc_a2(); + self.calc_a6(); + let u = &self.a + * (self.a4.as_ref().unwrap() * b[5] + + self.a2.as_ref().unwrap() * b[3] + + &self.ident * b[1]); + let v = self.a4.as_ref().unwrap() * b[4] + + self.a2.as_ref().unwrap() * b[2] + + &self.ident * b[0]; + (u, v) + } + + fn pade7(&mut self) -> (MatrixN, MatrixN) { + let b: [N; 8] = [ + convert(17297280.0), + convert(8648640.0), + convert(1995840.0), + convert(277200.0), + convert(25200.0), + convert(1512.0), + convert(56.0), + convert(1.0), + ]; + self.calc_a2(); + self.calc_a4(); + self.calc_a6(); + let u = &self.a + * (self.a6.as_ref().unwrap() * b[7] + + self.a4.as_ref().unwrap() * b[5] + + self.a2.as_ref().unwrap() * b[3] + + &self.ident * b[1]); + let v = self.a6.as_ref().unwrap() * b[6] + + self.a4.as_ref().unwrap() * b[4] + + self.a2.as_ref().unwrap() * b[2] + + &self.ident * b[0]; + (u, v) + } + + fn pade9(&mut self) -> (MatrixN, MatrixN) { + let b: [N; 10] = [ + convert(17643225600.0), + convert(8821612800.0), + convert(2075673600.0), + convert(302702400.0), + convert(30270240.0), + convert(2162160.0), + convert(110880.0), + convert(3960.0), + convert(90.0), + convert(1.0), + ]; + self.calc_a2(); + self.calc_a4(); + self.calc_a6(); + self.calc_a8(); + let u = &self.a + * (self.a8.as_ref().unwrap() * b[9] + + self.a6.as_ref().unwrap() * b[7] + + self.a4.as_ref().unwrap() * b[5] + + self.a2.as_ref().unwrap() * b[3] + + &self.ident * b[1]); + let v = self.a8.as_ref().unwrap() * b[8] + + self.a6.as_ref().unwrap() * b[6] + + self.a4.as_ref().unwrap() * b[4] + + self.a2.as_ref().unwrap() * b[2] + + &self.ident * b[0]; + (u, v) + } + + fn pade13_scaled(&mut self, s: u64) -> (MatrixN, MatrixN) { + let b: [N; 14] = [ + convert(64764752532480000.0), + convert(32382376266240000.0), + convert(7771770303897600.0), + convert(1187353796428800.0), + convert(129060195264000.0), + convert(10559470521600.0), + convert(670442572800.0), + convert(33522128640.0), + convert(1323241920.0), + convert(40840800.0), + convert(960960.0), + convert(16380.0), + convert(182.0), + convert(1.0), + ]; + let s = s as f64; + + let mb = &self.a * convert::(2.0_f64.powf(-s)); + self.calc_a2(); + self.calc_a4(); + self.calc_a6(); + 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] + &mb4 * b[11] + &mb2 * b[9]); + let u = &mb * (&u2 + &mb6 * b[7] + &mb4 * b[5] + &mb2 * b[3] + &self.ident * b[1]); + let v2 = &mb6 * (&mb6 * b[12] + &mb4 * b[10] + &mb2 * b[8]); + let v = v2 + &mb6 * b[6] + &mb4 * b[4] + &mb2 * b[2] + &self.ident * b[0]; + (u, v) + } +} + +fn factorial(n: u128) -> u128 { + if n == 1 { + return 1; + } + 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: &MatrixN, p: u64) -> N +where + N: RealField, + D: Dim, + DefaultAllocator: Allocator + Allocator, +{ + let nrows = a.data.shape().0; + let mut v = crate::VectorN::::repeat_generic(nrows, U1, convert(1.0)); + let m = a.transpose(); + + for _ in 0..p { + v = &m * v; + } + + v.max() +} + +fn ell(a: &MatrixN, m: u64) -> u64 +where + N: RealField, + D: Dim, + DefaultAllocator: Allocator + Allocator, +{ + // 2m choose m = (2m)!/(m! * (2m-m)!) + + let a_abs_onenorm = onenorm_matrix_power_nonm(&a.abs(), 2 * m + 1); + + if a_abs_onenorm == N::zero() { + 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); + let alpha = a_abs_onenorm / one_norm(a); + let alpha: f64 = try_convert(alpha).unwrap() / abs_c_recip as f64; + + let u = 2_f64.powf(-53.0); + let log2_alpha_div_u = (alpha / u).log2(); + let value = (log2_alpha_div_u / (2.0 * m as f64)).ceil(); + if value > 0.0 { + value as u64 + } else { + 0 + } +} + +fn solve_p_q(u: MatrixN, v: MatrixN) -> MatrixN +where + N: ComplexField, + D: DimMin, + DefaultAllocator: Allocator + Allocator<(usize, usize), DimMinimum>, +{ + let p = &u + &v; + let q = &v - &u; + + q.lu().solve(&p).unwrap() +} + +fn one_norm(m: &MatrixN) -> N +where + N: RealField, + D: Dim, + DefaultAllocator: Allocator, +{ + let mut max = N::zero(); + + for i in 0..m.ncols() { + let col = m.column(i); + max = max.max(col.iter().fold(N::zero(), |a, b| a + b.abs())); + } + + max +} + +impl MatrixN +where + D: DimMin, + DefaultAllocator: + Allocator + Allocator<(usize, usize), DimMinimum> + Allocator, +{ + /// Computes exponential of this matrix + pub fn exp(&self) -> Self { + // Simple case + if self.nrows() == 1 { + return self.map(|v| v.exp()); + } + + let mut h = ExpmPadeHelper::new(self.clone(), true); + + let eta_1 = N::max(h.d4_loose(), h.d6_loose()); + if eta_1 < convert(1.495585217958292e-002) && ell(&h.a, 3) == 0 { + let (u, v) = h.pade3(); + return solve_p_q(u, v); + } + + let eta_2 = N::max(h.d4_tight(), h.d6_loose()); + if eta_2 < convert(2.539398330063230e-001) && ell(&h.a, 5) == 0 { + let (u, v) = h.pade5(); + return solve_p_q(u, v); + } + + let eta_3 = N::max(h.d6_tight(), h.d8_loose()); + if eta_3 < convert(9.504178996162932e-001) && ell(&h.a, 7) == 0 { + let (u, v) = h.pade7(); + return solve_p_q(u, v); + } + if eta_3 < convert(2.097847961257068e+000) && ell(&h.a, 9) == 0 { + let (u, v) = h.pade9(); + return solve_p_q(u, v); + } + + let eta_4 = N::max(h.d8_loose(), h.d10_loose()); + let eta_5 = N::min(eta_3, eta_4); + let theta_13 = convert(4.25); + + let mut s = if eta_5 == N::zero() { + 0 + } else { + let l2 = try_convert((eta_5 / theta_13).log2().ceil()).unwrap(); + + if l2 < 0.0 { + 0 + } else { + l2 as u64 + } + }; + + s += ell(&(&h.a * convert::(2.0_f64.powf(-(s as f64)))), 13); + + let (u, v) = h.pade13_scaled(s); + let mut x = solve_p_q(u, v); + + for _ in 0..s { + x = &x * &x; + } + x + } +} + +#[cfg(test)] +mod tests { + #[test] + 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); + + assert_eq!(super::one_norm(&m), 19.0); + } +} diff --git a/src/linalg/mod.rs b/src/linalg/mod.rs index de1108f7..f96cef0c 100644 --- a/src/linalg/mod.rs +++ b/src/linalg/mod.rs @@ -5,6 +5,7 @@ mod bidiagonal; mod cholesky; mod convolution; mod determinant; +mod exp; mod full_piv_lu; pub mod givens; mod hessenberg; @@ -26,6 +27,7 @@ mod symmetric_tridiagonal; pub use self::bidiagonal::*; pub use self::cholesky::*; pub use self::convolution::*; +pub use self::exp::*; pub use self::full_piv_lu::*; pub use self::hessenberg::*; pub use self::lu::*; diff --git a/src/linalg/solve.rs b/src/linalg/solve.rs index 56db4ade..ac5bff46 100644 --- a/src/linalg/solve.rs +++ b/src/linalg/solve.rs @@ -1,4 +1,5 @@ use simba::scalar::ComplexField; +use simba::simd::SimdComplexField; use crate::base::allocator::Allocator; use crate::base::constraint::{SameNumberOfRows, ShapeConstraint}; @@ -432,3 +433,336 @@ impl> SquareMatrix { true } } + +/* + * + * SIMD-compatible unchecked versions. + * + */ + +impl> SquareMatrix { + /// Computes the solution of the linear system `self . x = b` where `x` is the unknown and only + /// the lower-triangular part of `self` (including the diagonal) is considered not-zero. + #[inline] + pub fn solve_lower_triangular_unchecked( + &self, + b: &Matrix, + ) -> MatrixMN + where + S2: Storage, + DefaultAllocator: Allocator, + ShapeConstraint: SameNumberOfRows, + { + let mut res = b.clone_owned(); + self.solve_lower_triangular_unchecked_mut(&mut res); + res + } + + /// Computes the solution of the linear system `self . x = b` where `x` is the unknown and only + /// the upper-triangular part of `self` (including the diagonal) is considered not-zero. + #[inline] + pub fn solve_upper_triangular_unchecked( + &self, + b: &Matrix, + ) -> MatrixMN + where + S2: Storage, + DefaultAllocator: Allocator, + ShapeConstraint: SameNumberOfRows, + { + let mut res = b.clone_owned(); + self.solve_upper_triangular_unchecked_mut(&mut res); + res + } + + /// Solves the linear system `self . x = b` where `x` is the unknown and only the + /// lower-triangular part of `self` (including the diagonal) is considered not-zero. + pub fn solve_lower_triangular_unchecked_mut( + &self, + b: &mut Matrix, + ) where + S2: StorageMut, + ShapeConstraint: SameNumberOfRows, + { + for i in 0..b.ncols() { + self.solve_lower_triangular_vector_unchecked_mut(&mut b.column_mut(i)); + } + } + + fn solve_lower_triangular_vector_unchecked_mut(&self, b: &mut Vector) + where + S2: StorageMut, + ShapeConstraint: SameNumberOfRows, + { + let dim = self.nrows(); + + for i in 0..dim { + let coeff; + + unsafe { + let diag = *self.get_unchecked((i, i)); + coeff = *b.vget_unchecked(i) / diag; + *b.vget_unchecked_mut(i) = coeff; + } + + b.rows_range_mut(i + 1..) + .axpy(-coeff, &self.slice_range(i + 1.., i), N::one()); + } + } + + // FIXME: add the same but for solving upper-triangular. + /// Solves the linear system `self . x = b` where `x` is the unknown and only the + /// lower-triangular part of `self` is considered not-zero. The diagonal is never read as it is + /// assumed to be equal to `diag`. Returns `false` and does not modify its inputs if `diag` is zero. + pub fn solve_lower_triangular_with_diag_unchecked_mut( + &self, + b: &mut Matrix, + diag: N, + ) where + S2: StorageMut, + ShapeConstraint: SameNumberOfRows, + { + let dim = self.nrows(); + let cols = b.ncols(); + + for k in 0..cols { + let mut bcol = b.column_mut(k); + + for i in 0..dim - 1 { + let coeff = unsafe { *bcol.vget_unchecked(i) } / diag; + bcol.rows_range_mut(i + 1..) + .axpy(-coeff, &self.slice_range(i + 1.., i), N::one()); + } + } + } + + /// Solves the linear system `self . x = b` where `x` is the unknown and only the + /// upper-triangular part of `self` (including the diagonal) is considered not-zero. + pub fn solve_upper_triangular_unchecked_mut( + &self, + b: &mut Matrix, + ) where + S2: StorageMut, + ShapeConstraint: SameNumberOfRows, + { + for i in 0..b.ncols() { + self.solve_upper_triangular_vector_unchecked_mut(&mut b.column_mut(i)) + } + } + + fn solve_upper_triangular_vector_unchecked_mut(&self, b: &mut Vector) + where + S2: StorageMut, + ShapeConstraint: SameNumberOfRows, + { + let dim = self.nrows(); + + for i in (0..dim).rev() { + let coeff; + + unsafe { + let diag = *self.get_unchecked((i, i)); + coeff = *b.vget_unchecked(i) / diag; + *b.vget_unchecked_mut(i) = coeff; + } + + b.rows_range_mut(..i) + .axpy(-coeff, &self.slice_range(..i, i), N::one()); + } + } + + /* + * + * Transpose and adjoint versions + * + */ + /// Computes the solution of the linear system `self.transpose() . x = b` where `x` is the unknown and only + /// the lower-triangular part of `self` (including the diagonal) is considered not-zero. + #[inline] + pub fn tr_solve_lower_triangular_unchecked( + &self, + b: &Matrix, + ) -> MatrixMN + where + S2: Storage, + DefaultAllocator: Allocator, + ShapeConstraint: SameNumberOfRows, + { + let mut res = b.clone_owned(); + self.tr_solve_lower_triangular_unchecked_mut(&mut res); + res + } + + /// Computes the solution of the linear system `self.transpose() . x = b` where `x` is the unknown and only + /// the upper-triangular part of `self` (including the diagonal) is considered not-zero. + #[inline] + pub fn tr_solve_upper_triangular_unchecked( + &self, + b: &Matrix, + ) -> MatrixMN + where + S2: Storage, + DefaultAllocator: Allocator, + ShapeConstraint: SameNumberOfRows, + { + let mut res = b.clone_owned(); + self.tr_solve_upper_triangular_unchecked_mut(&mut res); + res + } + + /// Solves the linear system `self.transpose() . x = b` where `x` is the unknown and only the + /// lower-triangular part of `self` (including the diagonal) is considered not-zero. + pub fn tr_solve_lower_triangular_unchecked_mut( + &self, + b: &mut Matrix, + ) where + S2: StorageMut, + ShapeConstraint: SameNumberOfRows, + { + for i in 0..b.ncols() { + self.xx_solve_lower_triangular_vector_unchecked_mut( + &mut b.column_mut(i), + |e| e, + |a, b| a.dot(b), + ) + } + } + + /// Solves the linear system `self.transpose() . x = b` where `x` is the unknown and only the + /// upper-triangular part of `self` (including the diagonal) is considered not-zero. + pub fn tr_solve_upper_triangular_unchecked_mut( + &self, + b: &mut Matrix, + ) where + S2: StorageMut, + ShapeConstraint: SameNumberOfRows, + { + for i in 0..b.ncols() { + self.xx_solve_upper_triangular_vector_unchecked_mut( + &mut b.column_mut(i), + |e| e, + |a, b| a.dot(b), + ) + } + } + + /// Computes the solution of the linear system `self.adjoint() . x = b` where `x` is the unknown and only + /// the lower-triangular part of `self` (including the diagonal) is considered not-zero. + #[inline] + pub fn ad_solve_lower_triangular_unchecked( + &self, + b: &Matrix, + ) -> MatrixMN + where + S2: Storage, + DefaultAllocator: Allocator, + ShapeConstraint: SameNumberOfRows, + { + let mut res = b.clone_owned(); + self.ad_solve_lower_triangular_unchecked_mut(&mut res); + res + } + + /// Computes the solution of the linear system `self.adjoint() . x = b` where `x` is the unknown and only + /// the upper-triangular part of `self` (including the diagonal) is considered not-zero. + #[inline] + pub fn ad_solve_upper_triangular_unchecked( + &self, + b: &Matrix, + ) -> MatrixMN + where + S2: Storage, + DefaultAllocator: Allocator, + ShapeConstraint: SameNumberOfRows, + { + let mut res = b.clone_owned(); + self.ad_solve_upper_triangular_unchecked_mut(&mut res); + res + } + + /// Solves the linear system `self.adjoint() . x = b` where `x` is the unknown and only the + /// lower-triangular part of `self` (including the diagonal) is considered not-zero. + pub fn ad_solve_lower_triangular_unchecked_mut( + &self, + b: &mut Matrix, + ) where + S2: StorageMut, + ShapeConstraint: SameNumberOfRows, + { + for i in 0..b.ncols() { + self.xx_solve_lower_triangular_vector_unchecked_mut( + &mut b.column_mut(i), + |e| e.simd_conjugate(), + |a, b| a.dotc(b), + ) + } + } + + /// Solves the linear system `self.adjoint() . x = b` where `x` is the unknown and only the + /// upper-triangular part of `self` (including the diagonal) is considered not-zero. + pub fn ad_solve_upper_triangular_unchecked_mut( + &self, + b: &mut Matrix, + ) where + S2: StorageMut, + ShapeConstraint: SameNumberOfRows, + { + for i in 0..b.ncols() { + self.xx_solve_upper_triangular_vector_unchecked_mut( + &mut b.column_mut(i), + |e| e.simd_conjugate(), + |a, b| a.dotc(b), + ) + } + } + + #[inline(always)] + fn xx_solve_lower_triangular_vector_unchecked_mut( + &self, + b: &mut Vector, + conjugate: impl Fn(N) -> N, + dot: impl Fn( + &DVectorSlice, + &DVectorSlice, + ) -> N, + ) where + S2: StorageMut, + ShapeConstraint: SameNumberOfRows, + { + let dim = self.nrows(); + + for i in (0..dim).rev() { + let dot = dot(&self.slice_range(i + 1.., i), &b.slice_range(i + 1.., 0)); + + unsafe { + let b_i = b.vget_unchecked_mut(i); + let diag = conjugate(*self.get_unchecked((i, i))); + *b_i = (*b_i - dot) / diag; + } + } + } + + #[inline(always)] + fn xx_solve_upper_triangular_vector_unchecked_mut( + &self, + b: &mut Vector, + conjugate: impl Fn(N) -> N, + dot: impl Fn( + &DVectorSlice, + &DVectorSlice, + ) -> N, + ) where + S2: StorageMut, + ShapeConstraint: SameNumberOfRows, + { + for i in 0..self.nrows() { + let dot = dot(&self.slice_range(..i, i), &b.slice_range(..i, 0)); + + unsafe { + let b_i = b.vget_unchecked_mut(i); + let diag = conjugate(*self.get_unchecked((i, i))); + *b_i = (*b_i - dot) / diag; + } + } + } +} diff --git a/src/linalg/svd.rs b/src/linalg/svd.rs index 6d74a438..4be0ecdb 100644 --- a/src/linalg/svd.rs +++ b/src/linalg/svd.rs @@ -95,7 +95,7 @@ where /// # 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 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 @@ -626,7 +626,7 @@ where /// # 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 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 diff --git a/tests/core/conversion.rs b/tests/core/conversion.rs index 74aa41e9..b7a8c5f8 100644 --- a/tests/core/conversion.rs +++ b/tests/core/conversion.rs @@ -1,4 +1,4 @@ -#![cfg(feature = "arbitrary")] +#![cfg(all(feature = "arbitrary", feature = "alga"))] use alga::linear::Transformation; use na::{ self, Affine3, Isometry3, Matrix2, Matrix2x3, Matrix2x4, Matrix2x5, Matrix2x6, Matrix3, diff --git a/tests/core/edition.rs b/tests/core/edition.rs index 7c2dee0f..dac25d45 100644 --- a/tests/core/edition.rs +++ b/tests/core/edition.rs @@ -1,13 +1,11 @@ -#![cfg_attr(rustfmt, rustfmt_skip)] - -use na::{Matrix, - DMatrix, - Matrix3, Matrix4, Matrix5, - Matrix4x3, Matrix3x4, Matrix5x3, Matrix3x5, Matrix4x5, Matrix5x4}; +use na::{ + DMatrix, Matrix, Matrix3, Matrix3x4, Matrix3x5, Matrix4, Matrix4x3, Matrix4x5, Matrix5, + Matrix5x3, Matrix5x4, +}; use na::{Dynamic, U2, U3, U5}; - #[test] +#[rustfmt::skip] fn upper_lower_triangular() { let m = Matrix4::new( 11.0, 12.0, 13.0, 14.0, @@ -173,6 +171,7 @@ fn upper_lower_triangular() { } #[test] +#[rustfmt::skip] fn swap_rows() { let mut m = Matrix5x3::new( 11.0, 12.0, 13.0, @@ -194,6 +193,7 @@ fn swap_rows() { } #[test] +#[rustfmt::skip] fn swap_columns() { let mut m = Matrix3x5::new( 11.0, 12.0, 13.0, 14.0, 15.0, @@ -211,6 +211,7 @@ fn swap_columns() { } #[test] +#[rustfmt::skip] fn remove_columns() { let m = Matrix3x5::new( 11, 12, 13, 14, 15, @@ -261,6 +262,7 @@ fn remove_columns() { } #[test] +#[rustfmt::skip] fn remove_columns_at() { let m = DMatrix::from_row_slice(5, 5, &[ 11, 12, 13, 14, 15, @@ -317,8 +319,8 @@ fn remove_columns_at() { assert_eq!(m.remove_columns_at(&[0,3,4]), expected3); } - #[test] +#[rustfmt::skip] fn remove_rows() { let m = Matrix5x3::new( 11, 12, 13, @@ -374,6 +376,7 @@ fn remove_rows() { } #[test] +#[rustfmt::skip] fn remove_rows_at() { let m = DMatrix::from_row_slice(5, 5, &[ 11, 12, 13, 14, 15, @@ -424,8 +427,8 @@ fn remove_rows_at() { assert_eq!(m.remove_rows_at(&[0,3,4]), expected3); } - #[test] +#[rustfmt::skip] fn insert_columns() { let m = Matrix5x3::new( 11, 12, 13, @@ -490,6 +493,7 @@ fn insert_columns() { } #[test] +#[rustfmt::skip] fn insert_columns_to_empty_matrix() { let m1 = DMatrix::repeat(0, 0, 0); let m2 = DMatrix::repeat(3, 0, 0); @@ -502,6 +506,7 @@ fn insert_columns_to_empty_matrix() { } #[test] +#[rustfmt::skip] fn insert_rows() { let m = Matrix3x5::new( 11, 12, 13, 14, 15, @@ -573,6 +578,7 @@ fn insert_rows_to_empty_matrix() { } #[test] +#[rustfmt::skip] fn resize() { let m = Matrix3x5::new( 11, 12, 13, 14, 15, diff --git a/tests/core/matrix.rs b/tests/core/matrix.rs index 9e25db58..6ade807f 100644 --- a/tests/core/matrix.rs +++ b/tests/core/matrix.rs @@ -945,7 +945,7 @@ mod normalization_tests { } } -#[cfg(feature = "arbitrary")] +#[cfg(all(feature = "arbitrary", feature = "alga"))] // FIXME: move this to alga ? mod finite_dim_inner_space_tests { use super::*; diff --git a/tests/core/matrix_slice.rs b/tests/core/matrix_slice.rs index 3a6a1aa5..5fd90b0f 100644 --- a/tests/core/matrix_slice.rs +++ b/tests/core/matrix_slice.rs @@ -1,18 +1,15 @@ #![allow(non_snake_case)] -#![cfg_attr(rustfmt, rustfmt_skip)] +use na::{ + DMatrix, DMatrixSlice, DMatrixSliceMut, Matrix2, Matrix2x3, Matrix2x4, Matrix2x6, Matrix3, + Matrix3x2, Matrix3x4, Matrix4x2, Matrix6x2, MatrixSlice2, MatrixSlice2x3, MatrixSlice2xX, + MatrixSlice3, MatrixSlice3x2, MatrixSliceMut2, MatrixSliceMut2x3, MatrixSliceMut2xX, + MatrixSliceMut3, MatrixSliceMut3x2, MatrixSliceMutXx3, MatrixSliceXx3, RowVector4, Vector3, +}; use na::{U2, U3, U4}; -use na::{DMatrix, - RowVector4, - Vector3, - Matrix2, Matrix3, - Matrix2x3, Matrix3x2, Matrix3x4, Matrix4x2, Matrix2x4, Matrix6x2, Matrix2x6, - MatrixSlice2, MatrixSlice3, MatrixSlice2x3, MatrixSlice3x2, - MatrixSliceXx3, MatrixSlice2xX, DMatrixSlice, - MatrixSliceMut2, MatrixSliceMut3, MatrixSliceMut2x3, MatrixSliceMut3x2, - MatrixSliceMutXx3, MatrixSliceMut2xX, DMatrixSliceMut}; #[test] +#[rustfmt::skip] fn nested_fixed_slices() { let a = Matrix3x4::new(11.0, 12.0, 13.0, 14.0, 21.0, 22.0, 23.0, 24.0, @@ -38,6 +35,7 @@ fn nested_fixed_slices() { } #[test] +#[rustfmt::skip] fn nested_slices() { let a = Matrix3x4::new(11.0, 12.0, 13.0, 14.0, 21.0, 22.0, 23.0, 24.0, @@ -63,6 +61,7 @@ fn nested_slices() { } #[test] +#[rustfmt::skip] fn slice_mut() { let mut a = Matrix3x4::new(11.0, 12.0, 13.0, 14.0, 21.0, 22.0, 23.0, 24.0, @@ -82,6 +81,7 @@ fn slice_mut() { } #[test] +#[rustfmt::skip] fn nested_row_slices() { let a = Matrix6x2::new(11.0, 12.0, 21.0, 22.0, @@ -105,6 +105,7 @@ fn nested_row_slices() { } #[test] +#[rustfmt::skip] fn row_slice_mut() { let mut a = Matrix6x2::new(11.0, 12.0, 21.0, 22.0, @@ -129,6 +130,7 @@ fn row_slice_mut() { } #[test] +#[rustfmt::skip] fn nested_col_slices() { let a = Matrix2x6::new(11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 21.0, 22.0, 23.0, 24.0, 25.0, 26.0); @@ -146,6 +148,7 @@ fn nested_col_slices() { } #[test] +#[rustfmt::skip] fn col_slice_mut() { let mut a = Matrix2x6::new(11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 21.0, 22.0, 23.0, 24.0, 25.0, 26.0); @@ -163,6 +166,7 @@ fn col_slice_mut() { } #[test] +#[rustfmt::skip] fn rows_range_pair() { let a = Matrix3x4::new(11.0, 12.0, 13.0, 14.0, 21.0, 22.0, 23.0, 24.0, @@ -180,6 +184,7 @@ fn rows_range_pair() { } #[test] +#[rustfmt::skip] fn columns_range_pair() { let a = Matrix3x4::new(11.0, 12.0, 13.0, 14.0, 21.0, 22.0, 23.0, 24.0, @@ -198,6 +203,7 @@ fn columns_range_pair() { } #[test] +#[rustfmt::skip] fn new_slice() { let data = [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, @@ -228,6 +234,7 @@ fn new_slice() { } #[test] +#[rustfmt::skip] fn new_slice_mut() { let data = [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, diff --git a/tests/core/mod.rs b/tests/core/mod.rs index 6e897c85..ec1c4e3e 100644 --- a/tests/core/mod.rs +++ b/tests/core/mod.rs @@ -3,9 +3,9 @@ mod abomonation; mod blas; mod conversion; mod edition; +mod empty; mod matrix; mod matrix_slice; -mod empty; #[cfg(feature = "mint")] mod mint; mod serde; diff --git a/tests/linalg/eigen.rs b/tests/linalg/eigen.rs index 8b2f8ed1..1269ed45 100644 --- a/tests/linalg/eigen.rs +++ b/tests/linalg/eigen.rs @@ -1,5 +1,3 @@ -#![cfg_attr(rustfmt, rustfmt_skip)] - use na::DMatrix; #[cfg(feature = "arbitrary")] @@ -67,6 +65,7 @@ mod quickcheck_tests { // Test proposed on the issue #176 of rulinalg. #[test] +#[rustfmt::skip] fn symmetric_eigen_singular_24x24() { let m = DMatrix::from_row_slice( 24, diff --git a/tests/linalg/exp.rs b/tests/linalg/exp.rs new file mode 100644 index 00000000..75122107 --- /dev/null +++ b/tests/linalg/exp.rs @@ -0,0 +1,129 @@ +#[cfg(test)] +mod tests { + //https://github.com/scipy/scipy/blob/c1372d8aa90a73d8a52f135529293ff4edb98fc8/scipy/sparse/linalg/tests/test_matfuncs.py + #[test] + fn exp_static() { + use nalgebra::{Matrix1, Matrix2, Matrix3}; + + { + let m = Matrix1::new(1.0); + + let f = m.exp(); + + assert!(relative_eq!(f, Matrix1::new(1_f64.exp()), epsilon = 1.0e-7)); + } + + { + let m = Matrix2::new(0.0, 1.0, 0.0, 0.0); + + assert!(relative_eq!( + m.exp(), + Matrix2::new(1.0, 1.0, 0.0, 1.0), + epsilon = 1.0e-7 + )); + } + + { + let a: f64 = 1.0; + let b: f64 = 2.0; + let c: f64 = 3.0; + let d: f64 = 4.0; + let m = Matrix2::new(a, b, c, d); + + let delta = ((a - d).powf(2.0) + 4.0 * b * c).sqrt(); + let delta_2 = delta / 2.0; + let ad_2 = (a + d) / 2.0; + let m11 = ad_2.exp() * (delta * delta_2.cosh() + (a - d) * delta_2.sinh()); + let m12 = 2.0 * b * ad_2.exp() * delta_2.sinh(); + let m21 = 2.0 * c * ad_2.exp() * delta_2.sinh(); + let m22 = ad_2.exp() * (delta * delta_2.cosh() + (d - a) * delta_2.sinh()); + + let f = Matrix2::new(m11, m12, m21, m22) / delta; + assert!(relative_eq!(f, m.exp(), epsilon = 1.0e-7)); + } + + { + // https://mathworld.wolfram.com/MatrixExponential.html + use rand::{ + distributions::{Distribution, Uniform}, + thread_rng, + }; + let mut rng = thread_rng(); + let dist = Uniform::new(-10.0, 10.0); + loop { + let a: f64 = dist.sample(&mut rng); + let b: f64 = dist.sample(&mut rng); + let c: f64 = dist.sample(&mut rng); + let d: f64 = dist.sample(&mut rng); + let m = Matrix2::new(a, b, c, d); + + let delta_sq = (a - d).powf(2.0) + 4.0 * b * c; + if delta_sq < 0.0 { + continue; + } + + let delta = delta_sq.sqrt(); + let delta_2 = delta / 2.0; + let ad_2 = (a + d) / 2.0; + let m11 = ad_2.exp() * (delta * delta_2.cosh() + (a - d) * delta_2.sinh()); + let m12 = 2.0 * b * ad_2.exp() * delta_2.sinh(); + let m21 = 2.0 * c * ad_2.exp() * delta_2.sinh(); + let m22 = ad_2.exp() * (delta * delta_2.cosh() + (d - a) * delta_2.sinh()); + + let f = Matrix2::new(m11, m12, m21, m22) / delta; + println!("a: {}", m); + assert!(relative_eq!(f, m.exp(), epsilon = 1.0e-7)); + break; + } + } + + { + let m = Matrix3::new(1.0, 3.0, 0.0, 0.0, 1.0, 5.0, 0.0, 0.0, 2.0); + + let e1 = 1.0_f64.exp(); + let e2 = 2.0_f64.exp(); + + let f = Matrix3::new( + e1, + 3.0 * e1, + 15.0 * (e2 - 2.0 * e1), + 0.0, + e1, + 5.0 * (e2 - e1), + 0.0, + 0.0, + e2, + ); + + assert!(relative_eq!(f, m.exp(), epsilon = 1.0e-7)); + } + } + + #[test] + fn exp_dynamic() { + use nalgebra::DMatrix; + + let m = DMatrix::from_row_slice(3, 3, &[1.0, 3.0, 0.0, 0.0, 1.0, 5.0, 0.0, 0.0, 2.0]); + + let e1 = 1.0_f64.exp(); + let e2 = 2.0_f64.exp(); + + let f = DMatrix::from_row_slice( + 3, + 3, + &[ + e1, + 3.0 * e1, + 15.0 * (e2 - 2.0 * e1), + 0.0, + e1, + 5.0 * (e2 - e1), + 0.0, + 0.0, + e2, + ], + ); + + assert!(relative_eq!(f, m.exp(), epsilon = 1.0e-7)); + } +} diff --git a/tests/linalg/full_piv_lu.rs b/tests/linalg/full_piv_lu.rs index 320fe241..0bb832cd 100644 --- a/tests/linalg/full_piv_lu.rs +++ b/tests/linalg/full_piv_lu.rs @@ -1,8 +1,7 @@ -#![cfg_attr(rustfmt, rustfmt_skip)] - use na::Matrix3; #[test] +#[rustfmt::skip] fn full_piv_lu_simple() { let m = Matrix3::new( 2.0, -1.0, 0.0, @@ -22,11 +21,11 @@ fn full_piv_lu_simple() { } #[test] +#[rustfmt::skip] fn full_piv_lu_simple_with_pivot() { - let m = Matrix3::new( - 0.0, -1.0, 2.0, - -1.0, 2.0, -1.0, - 2.0, -1.0, 0.0); + let m = Matrix3::new(0.0, -1.0, 2.0, + -1.0, 2.0, -1.0, + 2.0, -1.0, 0.0); let lu = m.full_piv_lu(); assert_eq!(lu.determinant(), -4.0); @@ -175,7 +174,6 @@ mod quickcheck_tests { gen_tests!(f64, RandScalar); } - /* #[test] fn swap_rows() { diff --git a/tests/linalg/inverse.rs b/tests/linalg/inverse.rs index f0be4dd7..ab641a79 100644 --- a/tests/linalg/inverse.rs +++ b/tests/linalg/inverse.rs @@ -1,5 +1,3 @@ -#![cfg_attr(rustfmt, rustfmt_skip)] - use na::{Matrix1, Matrix2, Matrix3, Matrix4, Matrix5}; #[test] @@ -11,6 +9,7 @@ fn matrix1_try_inverse() { } #[test] +#[rustfmt::skip] fn matrix2_try_inverse() { let a = Matrix2::new( 5.0, -2.0, -10.0, 1.0); @@ -23,6 +22,7 @@ fn matrix2_try_inverse() { } #[test] +#[rustfmt::skip] fn matrix3_try_inverse() { let a = Matrix3::new(-3.0, 2.0, 0.0, -6.0, 9.0, -2.0, @@ -37,6 +37,7 @@ fn matrix3_try_inverse() { } #[test] +#[rustfmt::skip] fn matrix4_try_inverse_issue_214() { let m1 = Matrix4::new( -0.34727043, 0.00000005397217, -0.000000000000003822135, -0.000000000000003821371, @@ -58,6 +59,7 @@ fn matrix4_try_inverse_issue_214() { } #[test] +#[rustfmt::skip] fn matrix5_try_inverse() { // Dimension 5 is chosen so that the inversion happens by Gaussian elimination. // (at the time of writing dimensions <= 3 are implemented as analytic formulas, but we choose @@ -90,6 +92,7 @@ fn matrix1_try_inverse_scaled_identity() { } #[test] +#[rustfmt::skip] fn matrix2_try_inverse_scaled_identity() { // A perfectly invertible matrix with // very small coefficients @@ -103,6 +106,7 @@ fn matrix2_try_inverse_scaled_identity() { } #[test] +#[rustfmt::skip] fn matrix3_try_inverse_scaled_identity() { // A perfectly invertible matrix with // very small coefficients @@ -118,6 +122,7 @@ fn matrix3_try_inverse_scaled_identity() { } #[test] +#[rustfmt::skip] fn matrix5_try_inverse_scaled_identity() { // A perfectly invertible matrix with // very small coefficients diff --git a/tests/linalg/lu.rs b/tests/linalg/lu.rs index 69c387d2..7fab6b01 100644 --- a/tests/linalg/lu.rs +++ b/tests/linalg/lu.rs @@ -1,8 +1,7 @@ -#![cfg_attr(rustfmt, rustfmt_skip)] - use na::Matrix3; #[test] +#[rustfmt::skip] fn lu_simple() { let m = Matrix3::new( 2.0, -1.0, 0.0, @@ -21,6 +20,7 @@ fn lu_simple() { } #[test] +#[rustfmt::skip] fn lu_simple_with_pivot() { let m = Matrix3::new( 0.0, -1.0, 2.0, @@ -41,7 +41,7 @@ fn lu_simple_with_pivot() { #[cfg(feature = "arbitrary")] mod quickcheck_tests { #[allow(unused_imports)] - use crate::core::helper::{RandScalar, RandComplex}; + use crate::core::helper::{RandComplex, RandScalar}; macro_rules! gen_tests( ($module: ident, $scalar: ty) => { diff --git a/tests/linalg/mod.rs b/tests/linalg/mod.rs index 234cac39..7fc01396 100644 --- a/tests/linalg/mod.rs +++ b/tests/linalg/mod.rs @@ -1,7 +1,9 @@ mod balancing; mod bidiagonal; mod cholesky; +mod convolution; mod eigen; +mod exp; mod full_piv_lu; mod hessenberg; mod inverse; @@ -10,5 +12,4 @@ mod qr; mod schur; mod solve; mod svd; -mod convolution; mod tridiagonal; diff --git a/tests/linalg/schur.rs b/tests/linalg/schur.rs index f9c923a2..2086ce2d 100644 --- a/tests/linalg/schur.rs +++ b/tests/linalg/schur.rs @@ -1,12 +1,11 @@ -#![cfg_attr(rustfmt, rustfmt_skip)] - use na::{DMatrix, Matrix3, Matrix4}; #[test] +#[rustfmt::skip] fn schur_simpl_mat3() { let m = Matrix3::new(-2.0, -4.0, 2.0, - -2.0, 1.0, 2.0, - 4.0, 2.0, 5.0); + -2.0, 1.0, 2.0, + 4.0, 2.0, 5.0); let schur = m.schur(); let (vecs, vals) = schur.unpack(); @@ -83,6 +82,7 @@ mod quickcheck_tests { } #[test] +#[rustfmt::skip] fn schur_static_mat4_fail() { let m = Matrix4::new( 33.32699857679677, 46.794945978960044, -20.792148817005838, 84.73945485997737, @@ -95,6 +95,7 @@ fn schur_static_mat4_fail() { } #[test] +#[rustfmt::skip] fn schur_static_mat4_fail2() { let m = Matrix4::new( 14.623586538485966, 7.646156622760756, -52.11923331576265, -97.50030223503413, @@ -107,6 +108,7 @@ fn schur_static_mat4_fail2() { } #[test] +#[rustfmt::skip] fn schur_static_mat3_fail() { let m = Matrix3::new( -21.58457553143394, -67.3881542667948, -14.619829849784338, @@ -119,6 +121,7 @@ fn schur_static_mat3_fail() { // Test proposed on the issue #176 of rulinalg. #[test] +#[rustfmt::skip] fn schur_singular() { let m = DMatrix::from_row_slice(24, 24, &[ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, diff --git a/tests/linalg/svd.rs b/tests/linalg/svd.rs index ca7bab4c..cd44b61d 100644 --- a/tests/linalg/svd.rs +++ b/tests/linalg/svd.rs @@ -1,4 +1,3 @@ -#![cfg_attr(rustfmt, rustfmt_skip)] use na::{DMatrix, Matrix6}; #[cfg(feature = "arbitrary")] @@ -160,9 +159,9 @@ mod quickcheck_tests { gen_tests!(f64, RandScalar); } - // Test proposed on the issue #176 of rulinalg. #[test] +#[rustfmt::skip] fn svd_singular() { let m = DMatrix::from_row_slice(24, 24, &[ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, @@ -202,6 +201,7 @@ fn svd_singular() { // Same as the previous test but with one additional row. #[test] +#[rustfmt::skip] fn svd_singular_vertical() { let m = DMatrix::from_row_slice(25, 24, &[ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, @@ -241,6 +241,7 @@ fn svd_singular_vertical() { // Same as the previous test but with one additional column. #[test] +#[rustfmt::skip] fn svd_singular_horizontal() { let m = DMatrix::from_row_slice(24, 25, &[ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, @@ -299,6 +300,7 @@ fn svd_identity() { } #[test] +#[rustfmt::skip] fn svd_with_delimited_subproblem() { let mut m = DMatrix::::from_element(10, 10, 0.0); m[(0,0)] = 1.0; m[(0,1)] = 2.0; @@ -334,6 +336,7 @@ fn svd_with_delimited_subproblem() { } #[test] +#[rustfmt::skip] fn svd_fail() { let m = Matrix6::new( 0.9299319121545955, 0.9955870335651049, 0.8824725266413644, 0.28966880207132295, 0.06102723649846409, 0.9311880746048009, @@ -351,6 +354,12 @@ fn svd_fail() { fn svd_err() { let m = DMatrix::from_element(10, 10, 0.0); let svd = m.clone().svd(false, false); - assert_eq!(Err("SVD recomposition: U and V^t have not been computed."), svd.clone().recompose()); - assert_eq!(Err("SVD pseudo inverse: the epsilon must be non-negative."), svd.clone().pseudo_inverse(-1.0)); -} \ No newline at end of file + assert_eq!( + Err("SVD recomposition: U and V^t have not been computed."), + svd.clone().recompose() + ); + assert_eq!( + Err("SVD pseudo inverse: the epsilon must be non-negative."), + svd.clone().pseudo_inverse(-1.0) + ); +}