Merge pull request #736 from rustsim/cholesky_aosoa
This commit is contained in:
commit
a8f746f911
|
@ -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,
|
||||
|
|
|
@ -57,7 +57,9 @@ where
|
|||
N: Default,
|
||||
{
|
||||
fn default() -> Self {
|
||||
ArrayStorage { data: Default::default() }
|
||||
ArrayStorage {
|
||||
data: Default::default(),
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
@ -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;
|
||||
|
|
|
@ -36,7 +36,7 @@ pub struct Quaternion<N: Scalar + SimdValue> {
|
|||
impl<N: RealField> Default for Quaternion<N> {
|
||||
fn default() -> Self {
|
||||
Quaternion {
|
||||
coords: Vector4::zeros()
|
||||
coords: Vector4::zeros(),
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
|
@ -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,21 +24,123 @@ 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: ComplexField, D: DimSub<Dynamic>> Cholesky<N, D>
|
||||
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]
|
||||
|
|
|
@ -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;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
|
@ -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,
|
||||
|
|
|
@ -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,
|
||||
|
|
|
@ -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;
|
||||
|
|
|
@ -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,
|
||||
|
|
|
@ -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() {
|
||||
|
|
|
@ -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
|
||||
|
|
|
@ -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) => {
|
||||
|
|
|
@ -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,
|
||||
|
|
|
@ -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)
|
||||
);
|
||||
}
|
||||
|
|
Loading…
Reference in New Issue