M-Labs
/
nalgebra
3
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'dev' 

# Conflicts:
#	src/linalg/cholesky.rs
This commit is contained in:
sebcrozet 2020-06-07 10:36:34 +02:00
parent 8fbd2b6d7e 2198b0e6b4
commit d5cbe56332
25 changed files with 1215 additions and 129 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
CHANGELOG.md
View File

@ -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.

2
Cargo.toml
View File

@ -1,6 +1,6 @@
[package]
name = "nalgebra"
version = "0.21.0"
version = "0.21.1"
authors = [ "Sébastien Crozet <developer@crozet.re>" ]
description = "Linear algebra library with transformations and statically-sized or dynamically-sized matrices."

24
nalgebra-glm/src/constructors.rs
View File

@ -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<N: Scalar>(x: N, y: N, z: N, w: N) -> TVec4<N> {
TVec4::new(x, y, z, w)
}
/// Create a new 2x2 matrix.
#[rustfmt::skip]
pub fn mat2<N: Scalar>(m11: N, m12: N,
m21: N, m22: N) -> TMat2<N> {
TMat::<N, U2, U2>::new(
@ -45,6 +44,7 @@ pub fn mat2<N: Scalar>(m11: N, m12: N,
}
/// Create a new 2x2 matrix.
#[rustfmt::skip]
pub fn mat2x2<N: Scalar>(m11: N, m12: N,
m21: N, m22: N) -> TMat2<N> {
TMat::<N, U2, U2>::new(
@ -54,6 +54,7 @@ pub fn mat2x2<N: Scalar>(m11: N, m12: N,
}
/// Create a new 2x3 matrix.
#[rustfmt::skip]
pub fn mat2x3<N: Scalar>(m11: N, m12: N, m13: N,
m21: N, m22: N, m23: N) -> TMat2x3<N> {
TMat::<N, U2, U3>::new(
@ -63,6 +64,7 @@ pub fn mat2x3<N: Scalar>(m11: N, m12: N, m13: N,
}
/// Create a new 2x4 matrix.
#[rustfmt::skip]
pub fn mat2x4<N: Scalar>(m11: N, m12: N, m13: N, m14: N,
m21: N, m22: N, m23: N, m24: N) -> TMat2x4<N> {
TMat::<N, U2, U4>::new(
@ -72,6 +74,7 @@ pub fn mat2x4<N: Scalar>(m11: N, m12: N, m13: N, m14: N,
}
/// Create a new 3x3 matrix.
#[rustfmt::skip]
pub fn mat3<N: Scalar>(m11: N, m12: N, m13: N,
m21: N, m22: N, m23: N,
m31: N, m32: N, m33: N) -> TMat3<N> {
@ -83,6 +86,7 @@ pub fn mat3<N: Scalar>(m11: N, m12: N, m13: N,
}
/// Create a new 3x2 matrix.
#[rustfmt::skip]
pub fn mat3x2<N: Scalar>(m11: N, m12: N,
m21: N, m22: N,
m31: N, m32: N) -> TMat3x2<N> {
@ -94,6 +98,7 @@ pub fn mat3x2<N: Scalar>(m11: N, m12: N,
}
/// Create a new 3x3 matrix.
#[rustfmt::skip]
pub fn mat3x3<N: Scalar>(m11: N, m12: N, m13: N,
m21: N, m22: N, m23: N,
m31: N, m32: N, m33: N) -> TMat3<N> {
@ -105,6 +110,7 @@ pub fn mat3x3<N: Scalar>(m11: N, m12: N, m13: N,
}
/// Create a new 3x4 matrix.
#[rustfmt::skip]
pub fn mat3x4<N: Scalar>(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<N> {
@ -116,6 +122,7 @@ pub fn mat3x4<N: Scalar>(m11: N, m12: N, m13: N, m14: N,
}
/// Create a new 4x2 matrix.
#[rustfmt::skip]
pub fn mat4x2<N: Scalar>(m11: N, m12: N,
m21: N, m22: N,
m31: N, m32: N,
@ -129,6 +136,7 @@ pub fn mat4x2<N: Scalar>(m11: N, m12: N,
}
/// Create a new 4x3 matrix.
#[rustfmt::skip]
pub fn mat4x3<N: Scalar>(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<N: Scalar>(m11: N, m12: N, m13: N,
}
/// Create a new 4x4 matrix.
#[rustfmt::skip]
pub fn mat4x4<N: Scalar>(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<N: Scalar>(m11: N, m12: N, m13: N, m14: N,
}
/// Create a new 4x4 matrix.
#[rustfmt::skip]
pub fn mat4<N: Scalar>(m11: N, m12: N, m13: N, m14: N,
m21: N, m22: N, m23: N, m24: N,
m31: N, m32: N, m33: N, m34: N,

15
src/base/array_storage.rs
View File

@ -48,6 +48,21 @@ where
/// Renamed to [ArrayStorage].
pub type MatrixArray<N, R, C> = ArrayStorage<N, R, C>;
impl<N, R, C> Default for ArrayStorage<N, R, C>
where
R: DimName,
C: DimName,
R::Value: Mul<C::Value>,
Prod<R::Value, C::Value>: ArrayLength<N>,
N: Default,
{
fn default() -> Self {
ArrayStorage {
data: Default::default(),
}
}
}
impl<N, R, C> Hash for ArrayStorage<N, R, C>
where
N: Hash,

12
src/base/construction_slice.rs
View File

@ -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<N, ..., S>
=> R: DimName, => C: DimName; // Type parameters for impl<N, ..., S>
R::name(), C::name(); // Arguments for `_generic` constructors.
); // Arguments for non-generic constructors.
=> R: DimName, => C: DimName; // Type parameters for impl<N, ..., S>
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<N, ..., S>
=> R: DimName, => C: DimName; // Type parameters for impl<N, ..., S>
R::name(), C::name(); // Arguments for `_generic` constructors.
); // Arguments for non-generic constructors.
=> R: DimName, => C: DimName; // Type parameters for impl<N, ..., S>
R::name(), C::name(); // Arguments for `_generic` constructors.
); // Arguments for non-generic constructors.
impl_constructors_mut!(R, Dynamic;
=> R: DimName;

15
src/base/matrix.rs
View File

@ -94,6 +94,21 @@ impl<N: Scalar, R: Dim, C: Dim, S: fmt::Debug> fmt::Debug for Matrix<N, R, C, S>
}
}
impl<N, R, C, S> Default for Matrix<N, R, C, S>
where
N: Scalar,
R: Dim,
C: Dim,
S: Default,
{
fn default() -> Self {
Matrix {
data: Default::default(),
_phantoms: PhantomData,
}
}
}
#[cfg(feature = "serde-serialize")]
impl<N, R, C, S> Serialize for Matrix<N, R, C, S>
where

14
src/geometry/quaternion.rs
View File

@ -33,6 +33,14 @@ pub struct Quaternion<N: Scalar + SimdValue> {
pub coords: Vector4<N>,
}
impl<N: RealField> Default for Quaternion<N> {
fn default() -> Self {
Quaternion {
coords: Vector4::zeros(),
}
}
}
#[cfg(feature = "abomonation-serialize")]
impl<N: SimdRealField> Abomonation for Quaternion<N>
where
@ -1536,6 +1544,12 @@ where
}
}
impl<N: RealField> Default for UnitQuaternion<N> {
fn default() -> Self {
Self::identity()
}
}
impl<N: RealField + fmt::Display> fmt::Display for UnitQuaternion<N> {
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
if let Some(axis) = self.axis() {

172
src/linalg/cholesky.rs
View File

@ -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<N, D>: Deserialize<'de>"))
)]
#[derive(Clone, Debug)]
pub struct Cholesky<N: ComplexField, D: Dim>
pub struct Cholesky<N: SimdComplexField, D: Dim>
where
DefaultAllocator: Allocator<N, D, D>,
{
chol: MatrixN<N, D>,
}
impl<N: ComplexField, D: Dim> Copy for Cholesky<N, D>
impl<N: SimdComplexField, D: Dim> Copy for Cholesky<N, D>
where
DefaultAllocator: Allocator<N, D, D>,
MatrixN<N, D>: Copy,
{
}
impl<N: SimdComplexField, D: Dim> Cholesky<N, D>
where
DefaultAllocator: Allocator<N, D, D>,
{
/// 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<N, D>) -> 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<N, D> {
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<N, D> {
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<N, D> {
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<N, D> {
&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<R2: Dim, C2: Dim, S2>(&self, b: &mut Matrix<N, R2, C2, S2>)
where
S2: StorageMut<N, R2, C2>,
ShapeConstraint: SameNumberOfRows<R2, D>,
{
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<R2: Dim, C2: Dim, S2>(&self, b: &Matrix<N, R2, C2, S2>) -> MatrixMN<N, R2, C2>
where
S2: Storage<N, R2, C2>,
DefaultAllocator: Allocator<N, R2, C2>,
ShapeConstraint: SameNumberOfRows<R2, D>,
{
let mut res = b.clone_owned();
self.solve_mut(&mut res);
res
}
/// Computes the inverse of the decomposed matrix.
pub fn inverse(&self) -> MatrixN<N, D> {
let shape = self.chol.data.shape();
let mut res = MatrixN::identity_generic(shape.0, shape.1);
self.solve_mut(&mut res);
res
}
}
impl<N: ComplexField, D: Dim> Cholesky<N, D>
where
DefaultAllocator: Allocator<N, D, D>,
@ -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<N, D> {
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<N, D> {
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<N, D> {
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<N, D> {
&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<R2: Dim, C2: Dim, S2>(&self, b: &mut Matrix<N, R2, C2, S2>)
where
S2: StorageMut<N, R2, C2>,
ShapeConstraint: SameNumberOfRows<R2, D>,
{
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<R2: Dim, C2: Dim, S2>(&self, b: &Matrix<N, R2, C2, S2>) -> MatrixMN<N, R2, C2>
where
S2: Storage<N, R2, C2>,
DefaultAllocator: Allocator<N, R2, C2>,
ShapeConstraint: SameNumberOfRows<R2, D>,
{
let mut res = b.clone_owned();
self.solve_mut(&mut res);
res
}
/// Computes the inverse of the decomposed matrix.
pub fn inverse(&self) -> MatrixN<N, D> {
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]

492
src/linalg/exp.rs Normal file
View File

@ -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<N, D>
where
N: RealField,
D: DimMin<D>,
DefaultAllocator: Allocator<N, D, D> + Allocator<(usize, usize), DimMinimum<D, D>>,
{
use_exact_norm: bool,
ident: MatrixN<N, D>,
a: MatrixN<N, D>,
a2: Option<MatrixN<N, D>>,
a4: Option<MatrixN<N, D>>,
a6: Option<MatrixN<N, D>>,
a8: Option<MatrixN<N, D>>,
a10: Option<MatrixN<N, D>>,
d4_exact: Option<N>,
d6_exact: Option<N>,
d8_exact: Option<N>,
d10_exact: Option<N>,
d4_approx: Option<N>,
d6_approx: Option<N>,
d8_approx: Option<N>,
d10_approx: Option<N>,
}
impl<N, D> ExpmPadeHelper<N, D>
where
N: RealField,
D: DimMin<D>,
DefaultAllocator: Allocator<N, D, D> + Allocator<(usize, usize), DimMinimum<D, D>>,
{
fn new(a: MatrixN<N, D>, use_exact_norm: bool) -> Self {
let (nrows, ncols) = a.data.shape();
ExpmPadeHelper {
use_exact_norm,
ident: MatrixN::<N, D>::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<N, D>, MatrixN<N, D>) {
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<N, D>, MatrixN<N, D>) {
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<N, D>, MatrixN<N, D>) {
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<N, D>, MatrixN<N, D>) {
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<N, D>, MatrixN<N, D>) {
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::<f64, N>(2.0_f64.powf(-s));
self.calc_a2();
self.calc_a4();
self.calc_a6();
let mb2 = self.a2.as_ref().unwrap() * convert::<f64, N>(2.0_f64.powf(-2.0 * s));
let mb4 = self.a4.as_ref().unwrap() * convert::<f64, N>(2.0.powf(-4.0 * s));
let mb6 = self.a6.as_ref().unwrap() * convert::<f64, N>(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<N, D>(a: &MatrixN<N, D>, p: u64) -> N
where
N: RealField,
D: Dim,
DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
{
let nrows = a.data.shape().0;
let mut v = crate::VectorN::<N, D>::repeat_generic(nrows, U1, convert(1.0));
let m = a.transpose();
for _ in 0..p {
v = &m * v;
}
v.max()
}
fn ell<N, D>(a: &MatrixN<N, D>, m: u64) -> u64
where
N: RealField,
D: Dim,
DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
{
// 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<N, D>(u: MatrixN<N, D>, v: MatrixN<N, D>) -> MatrixN<N, D>
where
N: ComplexField,
D: DimMin<D, Output = D>,
DefaultAllocator: Allocator<N, D, D> + Allocator<(usize, usize), DimMinimum<D, D>>,
{
let p = &u + &v;
let q = &v - &u;
q.lu().solve(&p).unwrap()
}
fn one_norm<N, D>(m: &MatrixN<N, D>) -> N
where
N: RealField,
D: Dim,
DefaultAllocator: Allocator<N, D, D>,
{
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<N: RealField, D> MatrixN<N, D>
where
D: DimMin<D, Output = D>,
DefaultAllocator:
Allocator<N, D, D> + Allocator<(usize, usize), DimMinimum<D, D>> + Allocator<N, D>,
{
/// 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::<f64, N>(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);
}
}

2
src/linalg/mod.rs
View File

@ -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::*;

334
src/linalg/solve.rs
View File

@ -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<N: ComplexField, D: Dim, S: Storage<N, D, D>> SquareMatrix<N, D, S> {
true
}
}
/*
*
* SIMD-compatible unchecked versions.
*
*/
impl<N: SimdComplexField, D: Dim, S: Storage<N, D, D>> SquareMatrix<N, D, S> {
/// 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<R2: Dim, C2: Dim, S2>(
&self,
b: &Matrix<N, R2, C2, S2>,
) -> MatrixMN<N, R2, C2>
where
S2: Storage<N, R2, C2>,
DefaultAllocator: Allocator<N, R2, C2>,
ShapeConstraint: SameNumberOfRows<R2, D>,
{
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<R2: Dim, C2: Dim, S2>(
&self,
b: &Matrix<N, R2, C2, S2>,
) -> MatrixMN<N, R2, C2>
where
S2: Storage<N, R2, C2>,
DefaultAllocator: Allocator<N, R2, C2>,
ShapeConstraint: SameNumberOfRows<R2, D>,
{
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<R2: Dim, C2: Dim, S2>(
&self,
b: &mut Matrix<N, R2, C2, S2>,
) where
S2: StorageMut<N, R2, C2>,
ShapeConstraint: SameNumberOfRows<R2, D>,
{
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<R2: Dim, S2>(&self, b: &mut Vector<N, R2, S2>)
where
S2: StorageMut<N, R2, U1>,
ShapeConstraint: SameNumberOfRows<R2, D>,
{
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<R2: Dim, C2: Dim, S2>(
&self,
b: &mut Matrix<N, R2, C2, S2>,
diag: N,
) where
S2: StorageMut<N, R2, C2>,
ShapeConstraint: SameNumberOfRows<R2, D>,
{
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<R2: Dim, C2: Dim, S2>(
&self,
b: &mut Matrix<N, R2, C2, S2>,
) where
S2: StorageMut<N, R2, C2>,
ShapeConstraint: SameNumberOfRows<R2, D>,
{
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<R2: Dim, S2>(&self, b: &mut Vector<N, R2, S2>)
where
S2: StorageMut<N, R2, U1>,
ShapeConstraint: SameNumberOfRows<R2, D>,
{
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<R2: Dim, C2: Dim, S2>(
&self,
b: &Matrix<N, R2, C2, S2>,
) -> MatrixMN<N, R2, C2>
where
S2: Storage<N, R2, C2>,
DefaultAllocator: Allocator<N, R2, C2>,
ShapeConstraint: SameNumberOfRows<R2, D>,
{
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<R2: Dim, C2: Dim, S2>(
&self,
b: &Matrix<N, R2, C2, S2>,
) -> MatrixMN<N, R2, C2>
where
S2: Storage<N, R2, C2>,
DefaultAllocator: Allocator<N, R2, C2>,
ShapeConstraint: SameNumberOfRows<R2, D>,
{
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<R2: Dim, C2: Dim, S2>(
&self,
b: &mut Matrix<N, R2, C2, S2>,
) where
S2: StorageMut<N, R2, C2>,
ShapeConstraint: SameNumberOfRows<R2, D>,
{
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<R2: Dim, C2: Dim, S2>(
&self,
b: &mut Matrix<N, R2, C2, S2>,
) where
S2: StorageMut<N, R2, C2>,
ShapeConstraint: SameNumberOfRows<R2, D>,
{
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<R2: Dim, C2: Dim, S2>(
&self,
b: &Matrix<N, R2, C2, S2>,
) -> MatrixMN<N, R2, C2>
where
S2: Storage<N, R2, C2>,
DefaultAllocator: Allocator<N, R2, C2>,
ShapeConstraint: SameNumberOfRows<R2, D>,
{
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<R2: Dim, C2: Dim, S2>(
&self,
b: &Matrix<N, R2, C2, S2>,
) -> MatrixMN<N, R2, C2>
where
S2: Storage<N, R2, C2>,
DefaultAllocator: Allocator<N, R2, C2>,
ShapeConstraint: SameNumberOfRows<R2, D>,
{
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<R2: Dim, C2: Dim, S2>(
&self,
b: &mut Matrix<N, R2, C2, S2>,
) where
S2: StorageMut<N, R2, C2>,
ShapeConstraint: SameNumberOfRows<R2, D>,
{
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<R2: Dim, C2: Dim, S2>(
&self,
b: &mut Matrix<N, R2, C2, S2>,
) where
S2: StorageMut<N, R2, C2>,
ShapeConstraint: SameNumberOfRows<R2, D>,
{
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<R2: Dim, S2>(
&self,
b: &mut Vector<N, R2, S2>,
conjugate: impl Fn(N) -> N,
dot: impl Fn(
&DVectorSlice<N, S::RStride, S::CStride>,
&DVectorSlice<N, S2::RStride, S2::CStride>,
) -> N,
) where
S2: StorageMut<N, R2, U1>,
ShapeConstraint: SameNumberOfRows<R2, D>,
{
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<R2: Dim, S2>(
&self,
b: &mut Vector<N, R2, S2>,
conjugate: impl Fn(N) -> N,
dot: impl Fn(
&DVectorSlice<N, S::RStride, S::CStride>,
&DVectorSlice<N, S2::RStride, S2::CStride>,
) -> N,
) where
S2: StorageMut<N, R2, U1>,
ShapeConstraint: SameNumberOfRows<R2, D>,
{
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;
}
}
}
}

4
src/linalg/svd.rs
View File

@ -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

2
tests/core/conversion.rs
View File

@ -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,

24
tests/core/edition.rs
View File

@ -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,

2
tests/core/matrix.rs
View File

@ -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::*;

27
tests/core/matrix_slice.rs
View File

@ -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,

2
tests/core/mod.rs
View File

@ -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;

3
tests/linalg/eigen.rs
View File

@ -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,

129
tests/linalg/exp.rs Normal file
View File

@ -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));
}
}

12
tests/linalg/full_piv_lu.rs
View File

@ -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<f64>);
}
/*
#[test]
fn swap_rows() {

9
tests/linalg/inverse.rs
View File

@ -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

6
tests/linalg/lu.rs
View File

@ -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) => {

3
tests/linalg/mod.rs
View File

@ -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;

11
tests/linalg/schur.rs
View File

@ -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,

19
tests/linalg/svd.rs
View File

@ -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<f64>);
}
// 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::<f64>::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));
}
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)
);
}