Checkpoint #10

This commit is contained in:
Violeta Hernández 2021-07-17 02:52:57 -05:00
parent 87fe2b30df
commit a6b8dd6d78
16 changed files with 267 additions and 85 deletions

View File

@ -31,6 +31,9 @@ type DefaultUninitBuffer<T, R, C> =
* Allocator. * Allocator.
* *
*/ */
/// A helper struct that controls how the storage for a matrix should be allocated.
///
/// This struct is useless on its own. Instead, it's used in trait
/// An allocator based on `GenericArray` and `VecStorage` for statically-sized and dynamically-sized /// An allocator based on `GenericArray` and `VecStorage` for statically-sized and dynamically-sized
/// matrices respectively. /// matrices respectively.
pub struct DefaultAllocator; pub struct DefaultAllocator;

View File

@ -152,7 +152,7 @@ pub type MatrixCross<T, R1, C1, R2, C2> =
/// dynamically-sized column vector should be represented as a `Matrix<T, Dynamic, U1, S>` (given /// dynamically-sized column vector should be represented as a `Matrix<T, Dynamic, U1, S>` (given
/// some concrete types for `T` and a compatible data storage type `S`). /// some concrete types for `T` and a compatible data storage type `S`).
#[repr(C)] #[repr(C)]
#[derive(Clone, Copy)] #[derive(Clone, Copy, Debug)]
pub struct Matrix<T, R, C, S> { pub struct Matrix<T, R, C, S> {
/// The data storage that contains all the matrix components. Disappointed? /// The data storage that contains all the matrix components. Disappointed?
/// ///
@ -192,15 +192,6 @@ pub struct Matrix<T, R, C, S> {
_phantoms: PhantomData<(T, R, C)>, _phantoms: PhantomData<(T, R, C)>,
} }
impl<T, R: Dim, C: Dim, S: fmt::Debug> fmt::Debug for Matrix<T, R, C, S> {
fn fmt(&self, formatter: &mut fmt::Formatter) -> Result<(), fmt::Error> {
formatter
.debug_struct("Matrix")
.field("data", &self.data)
.finish()
}
}
impl<T, R: Dim, C: Dim, S: Default> Default for Matrix<T, R, C, S> { impl<T, R: Dim, C: Dim, S: Default> Default for Matrix<T, R, C, S> {
fn default() -> Self { fn default() -> Self {
Matrix { Matrix {

View File

@ -228,7 +228,7 @@ impl<T> Unit<T> {
/// Wraps the given reference, assuming it is already normalized. /// Wraps the given reference, assuming it is already normalized.
#[inline] #[inline]
pub fn from_ref_unchecked(value: &T) -> &Self { pub fn from_ref_unchecked(value: &T) -> &Self {
unsafe { &*(value as *const T as *const Self) } unsafe { &*(value as *const _ as *const Self) }
} }
/// Retrieves the underlying value. /// Retrieves the underlying value.
@ -331,7 +331,7 @@ impl<T> Deref for Unit<T> {
#[inline] #[inline]
fn deref(&self) -> &T { fn deref(&self) -> &T {
unsafe { &*(self as *const Self as *const T) } unsafe { &*(self as *const _ as *const T) }
} }
} }

View File

@ -59,14 +59,14 @@ use std::ops::{
impl<T> AsRef<[T; 8]> for DualQuaternion<T> { impl<T> AsRef<[T; 8]> for DualQuaternion<T> {
#[inline] #[inline]
fn as_ref(&self) -> &[T; 8] { fn as_ref(&self) -> &[T; 8] {
unsafe { &*(self as *const Self as *const [T; 8]) } unsafe { &*(self as *const _ as *const [T; 8]) }
} }
} }
impl<T> AsMut<[T; 8]> for DualQuaternion<T> { impl<T> AsMut<[T; 8]> for DualQuaternion<T> {
#[inline] #[inline]
fn as_mut(&mut self) -> &mut [T; 8] { fn as_mut(&mut self) -> &mut [T; 8] {
unsafe { &mut *(self as *mut Self as *mut [T; 8]) } unsafe { &mut *(self as *mut _ as *mut [T; 8]) }
} }
} }

View File

@ -18,26 +18,27 @@ use crate::base::{Matrix4, Vector, Vector3};
use crate::geometry::{Point3, Projective3}; use crate::geometry::{Point3, Projective3};
/// A 3D orthographic projection stored as a homogeneous 4x4 matrix. /// A 3D orthographic projection stored as a homogeneous 4x4 matrix.
#[repr(C)]
pub struct Orthographic3<T> { pub struct Orthographic3<T> {
matrix: Matrix4<T>, matrix: Matrix4<T>,
} }
impl<T: RealField> Copy for Orthographic3<T> {} impl<T: Copy> Copy for Orthographic3<T> {}
impl<T: RealField> Clone for Orthographic3<T> { impl<T: Clone> Clone for Orthographic3<T> {
#[inline] #[inline]
fn clone(&self) -> Self { fn clone(&self) -> Self {
Self::from_matrix_unchecked(self.matrix) Self::from_matrix_unchecked(self.matrix)
} }
} }
impl<T: RealField> fmt::Debug for Orthographic3<T> { impl<T: fmt::Debug> fmt::Debug for Orthographic3<T> {
fn fmt(&self, f: &mut fmt::Formatter) -> Result<(), fmt::Error> { fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
self.matrix.fmt(f) self.matrix.fmt(f)
} }
} }
impl<T: RealField> PartialEq for Orthographic3<T> { impl<T: PartialEq> PartialEq for Orthographic3<T> {
#[inline] #[inline]
fn eq(&self, right: &Self) -> bool { fn eq(&self, right: &Self) -> bool {
self.matrix == right.matrix self.matrix == right.matrix
@ -45,7 +46,7 @@ impl<T: RealField> PartialEq for Orthographic3<T> {
} }
#[cfg(feature = "serde-serialize-no-std")] #[cfg(feature = "serde-serialize-no-std")]
impl<T: RealField + Serialize> Serialize for Orthographic3<T> { impl<T: Serialize> Serialize for Orthographic3<T> {
fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error> fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error>
where where
S: Serializer, S: Serializer,
@ -55,7 +56,7 @@ impl<T: RealField + Serialize> Serialize for Orthographic3<T> {
} }
#[cfg(feature = "serde-serialize-no-std")] #[cfg(feature = "serde-serialize-no-std")]
impl<'a, T: RealField + Deserialize<'a>> Deserialize<'a> for Orthographic3<T> { impl<'a, T: Deserialize<'a>> Deserialize<'a> for Orthographic3<T> {
fn deserialize<Des>(deserializer: Des) -> Result<Self, Des::Error> fn deserialize<Des>(deserializer: Des) -> Result<Self, Des::Error>
where where
Des: Deserializer<'a>, Des: Deserializer<'a>,
@ -66,7 +67,8 @@ impl<'a, T: RealField + Deserialize<'a>> Deserialize<'a> for Orthographic3<T> {
} }
} }
impl<T: RealField> Orthographic3<T> { /// # Basic methods and casts.
impl<T> Orthographic3<T> {
/// Creates a new orthographic projection matrix. /// Creates a new orthographic projection matrix.
/// ///
/// This follows the OpenGL convention, so this will flip the `z` axis. /// This follows the OpenGL convention, so this will flip the `z` axis.
@ -110,8 +112,11 @@ impl<T: RealField> Orthographic3<T> {
/// assert_relative_eq!(proj.project_point(&p8), Point3::new(-1.0, -1.0, -1.0)); /// assert_relative_eq!(proj.project_point(&p8), Point3::new(-1.0, -1.0, -1.0));
/// ``` /// ```
#[inline] #[inline]
pub fn new(left: T, right: T, bottom: T, top: T, znear: T, zfar: T) -> Self { pub fn new(left: T, right: T, bottom: T, top: T, znear: T, zfar: T) -> Self
let matrix = Matrix4::<T>::identity(); where
T: RealField,
{
let matrix = Matrix4::identity();
let mut res = Self::from_matrix_unchecked(matrix); let mut res = Self::from_matrix_unchecked(matrix);
res.set_left_and_right(left, right); res.set_left_and_right(left, right);
@ -145,7 +150,10 @@ impl<T: RealField> Orthographic3<T> {
/// Creates a new orthographic projection matrix from an aspect ratio and the vertical field of view. /// Creates a new orthographic projection matrix from an aspect ratio and the vertical field of view.
#[inline] #[inline]
pub fn from_fov(aspect: T, vfov: T, znear: T, zfar: T) -> Self { pub fn from_fov(aspect: T, vfov: T, znear: T, zfar: T) -> Self
where
T: RealField,
{
assert!( assert!(
znear != zfar, znear != zfar,
"The far plane must not be equal to the near plane." "The far plane must not be equal to the near plane."
@ -188,7 +196,10 @@ impl<T: RealField> Orthographic3<T> {
/// ``` /// ```
#[inline] #[inline]
#[must_use] #[must_use]
pub fn inverse(&self) -> Matrix4<T> { pub fn inverse(&self) -> Matrix4<T>
where
T: RealField,
{
let mut res = self.to_homogeneous(); let mut res = self.to_homogeneous();
let inv_m11 = T::one() / self.matrix[(0, 0)]; let inv_m11 = T::one() / self.matrix[(0, 0)];
@ -257,7 +268,8 @@ impl<T: RealField> Orthographic3<T> {
#[inline] #[inline]
#[must_use] #[must_use]
pub fn as_projective(&self) -> &Projective3<T> { pub fn as_projective(&self) -> &Projective3<T> {
unsafe { &*(self as *const Orthographic3<T> as *const Projective3<T>) } // Safety: Self and Projective3 are both #[repr(C)] of a matrix.
unsafe { &*(self as *const _ as *const Projective3<T>) }
} }
/// This transformation seen as a `Projective3`. /// This transformation seen as a `Projective3`.
@ -301,7 +313,10 @@ impl<T: RealField> Orthographic3<T> {
pub fn unwrap(self) -> Matrix4<T> { pub fn unwrap(self) -> Matrix4<T> {
self.matrix self.matrix
} }
}
/// # Mathematical methods.
impl<T: RealField> Orthographic3<T> {
/// The left offset of the view cuboid. /// The left offset of the view cuboid.
/// ///
/// ``` /// ```

View File

@ -33,7 +33,7 @@ impl<T: RealField> Clone for Perspective3<T> {
} }
impl<T: RealField> fmt::Debug for Perspective3<T> { impl<T: RealField> fmt::Debug for Perspective3<T> {
fn fmt(&self, f: &mut fmt::Formatter) -> Result<(), fmt::Error> { fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
self.matrix.fmt(f) self.matrix.fmt(f)
} }
} }
@ -139,7 +139,7 @@ impl<T: RealField> Perspective3<T> {
#[inline] #[inline]
#[must_use] #[must_use]
pub fn as_projective(&self) -> &Projective3<T> { pub fn as_projective(&self) -> &Projective3<T> {
unsafe { &*(self as *const Perspective3<T> as *const Projective3<T>) } unsafe { &*(self as *const _ as *const Projective3<T>) }
} }
/// This transformation seen as a `Projective3`. /// This transformation seen as a `Projective3`.

View File

@ -12,13 +12,13 @@ impl<T: Scalar + SimdValue> Deref for Quaternion<T> {
#[inline] #[inline]
fn deref(&self) -> &Self::Target { fn deref(&self) -> &Self::Target {
unsafe { &*(self as *const Self as *const Self::Target) } unsafe { &*(self as *const _ as *const Self::Target) }
} }
} }
impl<T: Scalar + SimdValue> DerefMut for Quaternion<T> { impl<T: Scalar + SimdValue> DerefMut for Quaternion<T> {
#[inline] #[inline]
fn deref_mut(&mut self) -> &mut Self::Target { fn deref_mut(&mut self) -> &mut Self::Target {
unsafe { &mut *(self as *mut Self as *mut Self::Target) } unsafe { &mut *(self as *mut _ as *mut Self::Target) }
} }
} }

View File

@ -18,14 +18,14 @@ macro_rules! deref_impl(
#[inline] #[inline]
fn deref(&self) -> &Self::Target { fn deref(&self) -> &Self::Target {
unsafe { &*(self as *const Translation<T, $D> as *const Self::Target) } unsafe { &*(self as *const _ as *const Self::Target) }
} }
} }
impl<T: Scalar> DerefMut for Translation<T, $D> { impl<T: Scalar> DerefMut for Translation<T, $D> {
#[inline] #[inline]
fn deref_mut(&mut self) -> &mut Self::Target { fn deref_mut(&mut self) -> &mut Self::Target {
unsafe { &mut *(self as *mut Translation<T, $D> as *mut Self::Target) } unsafe { &mut *(self as *mut _ as *mut Self::Target) }
} }
} }
} }

View File

@ -7,7 +7,7 @@ use crate::allocator::Allocator;
use crate::base::{DefaultAllocator, Matrix, OMatrix, OVector, Unit}; use crate::base::{DefaultAllocator, Matrix, OMatrix, OVector, Unit};
use crate::dimension::{Const, Dim, DimDiff, DimMin, DimMinimum, DimSub, U1}; use crate::dimension::{Const, Dim, DimDiff, DimMin, DimMinimum, DimSub, U1};
use crate::storage::{Owned, Storage}; use crate::storage::{Owned, Storage};
use crate::{Dynamic, }; use crate::Dynamic;
use simba::scalar::ComplexField; use simba::scalar::ComplexField;
use crate::geometry::Reflection; use crate::geometry::Reflection;

View File

@ -30,7 +30,6 @@ use crate::linalg::{householder, PermutationSequence};
PermutationSequence<DimMinimum<R, C>>: Deserialize<'de>, PermutationSequence<DimMinimum<R, C>>: Deserialize<'de>,
OVector<T, DimMinimum<R, C>>: Deserialize<'de>")) OVector<T, DimMinimum<R, C>>: Deserialize<'de>"))
)] )]
#[derive(Clone, Debug)]
pub struct ColPivQR<T: ComplexField, R: DimMin<C>, C: Dim> pub struct ColPivQR<T: ComplexField, R: DimMin<C>, C: Dim>
where where
DefaultAllocator: Allocator<T, R, C> DefaultAllocator: Allocator<T, R, C>
@ -53,6 +52,24 @@ where
{ {
} }
impl<T: ComplexField, R: DimMin<C>, C: Dim> Clone for ColPivQR<T, R, C>
where
DefaultAllocator: Allocator<T, R, C>
+ Allocator<T, DimMinimum<R, C>>
+ Allocator<(usize, usize), DimMinimum<R, C>>,
OMatrix<T, R, C>: Clone,
PermutationSequence<DimMinimum<R, C>>: Clone,
OVector<T, DimMinimum<R, C>>: Clone,
{
fn clone(&self) -> Self {
Self {
col_piv_qr: self.col_piv_qr.clone(),
p: self.p.clone(),
diag: self.diag.clone(),
}
}
}
impl<T: ComplexField, R: DimMin<C>, C: Dim> ColPivQR<T, R, C> impl<T: ComplexField, R: DimMin<C>, C: Dim> ColPivQR<T, R, C>
where where
DefaultAllocator: Allocator<T, R, C> DefaultAllocator: Allocator<T, R, C>
@ -66,14 +83,13 @@ where
let min_nrows_ncols = nrows.min(ncols); let min_nrows_ncols = nrows.min(ncols);
let mut p = PermutationSequence::identity_generic(min_nrows_ncols); let mut p = PermutationSequence::identity_generic(min_nrows_ncols);
let mut diag = let mut diag = Matrix::new_uninitialized_generic(min_nrows_ncols, Const::<1>);
unsafe { crate::unimplemented_or_uninitialized_generic!(min_nrows_ncols, Const::<1>) };
if min_nrows_ncols.value() == 0 { if min_nrows_ncols.value() == 0 {
return ColPivQR { return ColPivQR {
col_piv_qr: matrix, col_piv_qr: matrix,
p, p,
diag, diag: unsafe { diag.assume_init() },
}; };
} }
@ -83,13 +99,13 @@ where
matrix.swap_columns(i, col_piv); matrix.swap_columns(i, col_piv);
p.append_permutation(i, col_piv); p.append_permutation(i, col_piv);
householder::clear_column_unchecked(&mut matrix, &mut diag[i], i, 0, None); householder::clear_column_unchecked(&mut matrix, diag[i].as_mut_ptr(), i, 0, None);
} }
ColPivQR { ColPivQR {
col_piv_qr: matrix, col_piv_qr: matrix,
p, p,
diag, diag:unsafe{diag.assume_init()},
} }
} }

View File

@ -1,14 +1,11 @@
//! This module provides the matrix exponent (exp) function to square matrices. //! This module provides the matrix exponent (exp) function to square matrices.
//! //!
use crate::{ use crate::{ComplexField, OMatrix, RealField, base::{
base::{
allocator::Allocator, allocator::Allocator,
dimension::{Const, Dim, DimMin, DimMinimum}, dimension::{Const, Dim, DimMin, DimMinimum},
storage::Storage, storage::Storage,
DefaultAllocator, DefaultAllocator,
}, }, convert, storage::Owned, try_convert};
convert, try_convert, ComplexField, OMatrix, RealField,
};
use crate::num::Zero; use crate::num::Zero;
@ -433,6 +430,7 @@ where
+ Allocator<T, D> + Allocator<T, D>
+ Allocator<T::RealField, D> + Allocator<T::RealField, D>
+ Allocator<T::RealField, D, D>, + Allocator<T::RealField, D, D>,
Owned<T, D, D>: Clone,
{ {
/// Computes exponential of this matrix /// Computes exponential of this matrix
#[must_use] #[must_use]

View File

@ -1,3 +1,5 @@
use std::fmt;
#[cfg(feature = "serde-serialize-no-std")] #[cfg(feature = "serde-serialize-no-std")]
use serde::{Deserialize, Serialize}; use serde::{Deserialize, Serialize};
@ -27,8 +29,7 @@ use crate::linalg::PermutationSequence;
OMatrix<T, R, C>: Deserialize<'de>, OMatrix<T, R, C>: Deserialize<'de>,
PermutationSequence<DimMinimum<R, C>>: Deserialize<'de>")) PermutationSequence<DimMinimum<R, C>>: Deserialize<'de>"))
)] )]
#[derive(Clone, Debug)] pub struct FullPivLU<T, R: DimMin<C>, C: Dim>
pub struct FullPivLU<T: ComplexField, R: DimMin<C>, C: Dim>
where where
DefaultAllocator: Allocator<T, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>, DefaultAllocator: Allocator<T, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>,
{ {
@ -40,11 +41,41 @@ where
impl<T: ComplexField, R: DimMin<C>, C: Dim> Copy for FullPivLU<T, R, C> impl<T: ComplexField, R: DimMin<C>, C: Dim> Copy for FullPivLU<T, R, C>
where where
DefaultAllocator: Allocator<T, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>, DefaultAllocator: Allocator<T, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>,
OMatrix<T, R, C>: Copy,
PermutationSequence<DimMinimum<R, C>>: Copy, PermutationSequence<DimMinimum<R, C>>: Copy,
OMatrix<T, R, C>: Copy,
{ {
} }
impl<T: ComplexField, R: DimMin<C>, C: Dim> Clone for FullPivLU<T, R, C>
where
DefaultAllocator: Allocator<T, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>,
PermutationSequence<DimMinimum<R, C>>: Clone,
OMatrix<T, R, C>: Clone,
{
fn clone(&self) -> Self {
Self {
lu: self.lu.clone(),
p: self.p.clone(),
q: self.q.clone(),
}
}
}
impl<T: ComplexField, R: DimMin<C>, C: Dim> fmt::Debug for FullPivLU<T, R, C>
where
DefaultAllocator: Allocator<T, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>,
PermutationSequence<DimMinimum<R, C>>: fmt::Debug,
OMatrix<T, R, C>: fmt::Debug,
{
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
f.debug_struct("FullPivLU")
.field("lu", &self.lu)
.field("p", &self.p)
.field("q", &self.q)
.finish()
}
}
impl<T: ComplexField, R: DimMin<C>, C: Dim> FullPivLU<T, R, C> impl<T: ComplexField, R: DimMin<C>, C: Dim> FullPivLU<T, R, C>
where where
DefaultAllocator: Allocator<T, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>, DefaultAllocator: Allocator<T, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>,

View File

@ -1,10 +1,14 @@
use std::fmt;
use std::mem::MaybeUninit;
#[cfg(feature = "serde-serialize-no-std")] #[cfg(feature = "serde-serialize-no-std")]
use serde::{Deserialize, Serialize}; use serde::{Deserialize, Serialize};
use crate::allocator::Allocator; use crate::allocator::Allocator;
use crate::base::{DefaultAllocator, OMatrix, OVector}; use crate::base::{DefaultAllocator, OMatrix, OVector};
use crate::dimension::{Const, DimDiff, DimSub, U1}; use crate::dimension::{Const, DimDiff, DimSub, U1};
use crate::storage::Storage; use crate::storage::{Owned, Storage};
use crate::Matrix;
use simba::scalar::ComplexField; use simba::scalar::ComplexField;
use crate::linalg::householder; use crate::linalg::householder;
@ -25,7 +29,6 @@ use crate::linalg::householder;
OMatrix<T, D, D>: Deserialize<'de>, OMatrix<T, D, D>: Deserialize<'de>,
OVector<T, DimDiff<D, U1>>: Deserialize<'de>")) OVector<T, DimDiff<D, U1>>: Deserialize<'de>"))
)] )]
#[derive(Clone, Debug)]
pub struct Hessenberg<T: ComplexField, D: DimSub<U1>> pub struct Hessenberg<T: ComplexField, D: DimSub<U1>>
where where
DefaultAllocator: Allocator<T, D, D> + Allocator<T, DimDiff<D, U1>>, DefaultAllocator: Allocator<T, D, D> + Allocator<T, DimDiff<D, U1>>,
@ -37,20 +40,46 @@ where
impl<T: ComplexField, D: DimSub<U1>> Copy for Hessenberg<T, D> impl<T: ComplexField, D: DimSub<U1>> Copy for Hessenberg<T, D>
where where
DefaultAllocator: Allocator<T, D, D> + Allocator<T, DimDiff<D, U1>>, DefaultAllocator: Allocator<T, D, D> + Allocator<T, DimDiff<D, U1>>,
OMatrix<T, D, D>: Copy, Owned<T, D, D>: Copy,
OVector<T, DimDiff<D, U1>>: Copy, Owned<T, DimDiff<D, U1>>: Copy,
{ {
} }
impl<T: ComplexField, D: DimSub<U1>> Clone for Hessenberg<T, D>
where
DefaultAllocator: Allocator<T, D, D> + Allocator<T, DimDiff<D, U1>>,
Owned<T, D, D>: Clone,
Owned<T, DimDiff<D, U1>>: Clone,
{
fn clone(&self) -> Self {
Self {
hess: self.hess.clone(),
subdiag: self.subdiag.clone(),
}
}
}
impl<T: ComplexField, D: DimSub<U1>> fmt::Debug for Hessenberg<T, D>
where
DefaultAllocator: Allocator<T, D, D> + Allocator<T, DimDiff<D, U1>>,
Owned<T, D, D>: fmt::Debug,
Owned<T, DimDiff<D, U1>>: fmt::Debug,
{
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
f.debug_struct("Hessenberg")
.field("hess", &self.hess)
.field("subdiag", &self.subdiag)
.finish()
}
}
impl<T: ComplexField, D: DimSub<U1>> Hessenberg<T, D> impl<T: ComplexField, D: DimSub<U1>> Hessenberg<T, D>
where where
DefaultAllocator: Allocator<T, D, D> + Allocator<T, D> + Allocator<T, DimDiff<D, U1>>, DefaultAllocator: Allocator<T, D, D> + Allocator<T, D> + Allocator<T, DimDiff<D, U1>>,
{ {
/// Computes the Hessenberg decomposition using householder reflections. /// Computes the Hessenberg decomposition using householder reflections.
pub fn new(hess: OMatrix<T, D, D>) -> Self { pub fn new(hess: OMatrix<T, D, D>) -> Self {
let mut work = unsafe { let mut work = OVector::new_uninitialized_generic(hess.data.shape().0, Const::<1>);
crate::unimplemented_or_uninitialized_generic!(hess.data.shape().0, Const::<1>)
};
Self::new_with_workspace(hess, &mut work) Self::new_with_workspace(hess, &mut work)
} }
@ -58,7 +87,10 @@ where
/// ///
/// The workspace containing `D` elements must be provided but its content does not have to be /// The workspace containing `D` elements must be provided but its content does not have to be
/// initialized. /// initialized.
pub fn new_with_workspace(mut hess: OMatrix<T, D, D>, work: &mut OVector<T, D>) -> Self { pub fn new_with_workspace(
mut hess: OMatrix<T, D, D>,
work: &mut OVector<MaybeUninit<T>, D>,
) -> Self {
assert!( assert!(
hess.is_square(), hess.is_square(),
"Cannot compute the hessenberg decomposition of a non-square matrix." "Cannot compute the hessenberg decomposition of a non-square matrix."
@ -76,19 +108,29 @@ where
"Hessenberg: invalid workspace size." "Hessenberg: invalid workspace size."
); );
let mut subdiag = unsafe { let mut subdiag = Matrix::new_uninitialized_generic(dim.sub(Const::<1>), Const::<1>);
crate::unimplemented_or_uninitialized_generic!(dim.sub(Const::<1>), Const::<1>)
};
if dim.value() == 0 { if dim.value() == 0 {
return Hessenberg { hess, subdiag }; return Self {
hess,
subdiag: unsafe { subdiag.assume_init() },
};
} }
for ite in 0..dim.value() - 1 { for ite in 0..dim.value() - 1 {
householder::clear_column_unchecked(&mut hess, &mut subdiag[ite], ite, 1, Some(work)); householder::clear_column_unchecked(
&mut hess,
subdiag[ite].as_mut_ptr(),
ite,
1,
Some(work),
);
} }
Hessenberg { hess, subdiag } Self {
hess,
subdiag: unsafe { subdiag.assume_init() },
}
} }
/// Retrieves `(q, h)` with `q` the orthogonal matrix of this decomposition and `h` the /// Retrieves `(q, h)` with `q` the orthogonal matrix of this decomposition and `h` the
@ -117,7 +159,10 @@ where
/// This is less efficient than `.unpack_h()` as it allocates a new matrix. /// This is less efficient than `.unpack_h()` as it allocates a new matrix.
#[inline] #[inline]
#[must_use] #[must_use]
pub fn h(&self) -> OMatrix<T, D, D> { pub fn h(&self) -> OMatrix<T, D, D>
where
Owned<T, D, D>: Clone,
{
let dim = self.hess.nrows(); let dim = self.hess.nrows();
let mut res = self.hess.clone(); let mut res = self.hess.clone();
res.fill_lower_triangle(T::zero(), 2); res.fill_lower_triangle(T::zero(), 2);

View File

@ -51,7 +51,7 @@ pub fn clear_column_unchecked<T: ComplexField, R: Dim, C: Dim>(
diag_elt: *mut T, diag_elt: *mut T,
icol: usize, icol: usize,
shift: usize, shift: usize,
bilateral: Option<&mut OVector<T, R>>, bilateral: Option<&mut OVector<MaybeUninit<T>, R>>,
) where ) where
DefaultAllocator: Allocator<T, R, C> + Allocator<T, R>, DefaultAllocator: Allocator<T, R, C> + Allocator<T, R>,
{ {
@ -88,11 +88,14 @@ pub fn clear_row_unchecked<T: ComplexField, R: Dim, C: Dim>(
{ {
let (mut top, mut bottom) = matrix.rows_range_pair_mut(irow, irow + 1..); let (mut top, mut bottom) = matrix.rows_range_pair_mut(irow, irow + 1..);
let mut axis = axis_packed.rows_range_mut(irow + shift..); let mut axis = axis_packed.rows_range_mut(irow + shift..);
axis.tr_copy_from(&top.columns_range(irow + shift..)); axis.tr_copy_init_from(&top.columns_range(irow + shift..));
let mut axis = unsafe { axis.assume_init_mut() };
let (reflection_norm, not_zero) = reflection_axis_mut(&mut axis); let (reflection_norm, not_zero) = reflection_axis_mut(&mut axis);
axis.conjugate_mut(); // So that reflect_rows actually cancels the first row. axis.conjugate_mut(); // So that reflect_rows actually cancels the first row.
unsafe{ *diag_elt = reflection_norm;} unsafe {
*diag_elt = reflection_norm;
}
if not_zero { if not_zero {
let refl = Reflection::new(Unit::new_unchecked(axis), T::zero()); let refl = Reflection::new(Unit::new_unchecked(axis), T::zero());

View File

@ -1,3 +1,6 @@
use std::fmt;
use std::mem;
#[cfg(feature = "serde-serialize-no-std")] #[cfg(feature = "serde-serialize-no-std")]
use serde::{Deserialize, Serialize}; use serde::{Deserialize, Serialize};
@ -5,9 +8,8 @@ use crate::allocator::{Allocator, Reallocator};
use crate::base::{DefaultAllocator, Matrix, OMatrix, Scalar}; use crate::base::{DefaultAllocator, Matrix, OMatrix, Scalar};
use crate::constraint::{SameNumberOfRows, ShapeConstraint}; use crate::constraint::{SameNumberOfRows, ShapeConstraint};
use crate::dimension::{Dim, DimMin, DimMinimum}; use crate::dimension::{Dim, DimMin, DimMinimum};
use crate::storage::{Storage, StorageMut}; use crate::storage::{Owned, Storage, StorageMut};
use simba::scalar::{ComplexField, Field}; use simba::scalar::{ComplexField, Field};
use std::mem;
use crate::linalg::PermutationSequence; use crate::linalg::PermutationSequence;
@ -27,8 +29,7 @@ use crate::linalg::PermutationSequence;
OMatrix<T, R, C>: Deserialize<'de>, OMatrix<T, R, C>: Deserialize<'de>,
PermutationSequence<DimMinimum<R, C>>: Deserialize<'de>")) PermutationSequence<DimMinimum<R, C>>: Deserialize<'de>"))
)] )]
#[derive(Clone, Debug)] pub struct LU<T, R: DimMin<C>, C: Dim>
pub struct LU<T: ComplexField, R: DimMin<C>, C: Dim>
where where
DefaultAllocator: Allocator<T, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>, DefaultAllocator: Allocator<T, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>,
{ {
@ -36,14 +37,42 @@ where
p: PermutationSequence<DimMinimum<R, C>>, p: PermutationSequence<DimMinimum<R, C>>,
} }
impl<T: ComplexField, R: DimMin<C>, C: Dim> Copy for LU<T, R, C> impl<T: Copy, R: DimMin<C>, C: Dim> Copy for LU<T, R, C>
where where
DefaultAllocator: Allocator<T, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>, DefaultAllocator: Allocator<T, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>,
OMatrix<T, R, C>: Copy,
PermutationSequence<DimMinimum<R, C>>: Copy, PermutationSequence<DimMinimum<R, C>>: Copy,
Owned<T, R, C>: Copy,
{ {
} }
impl<T: Clone, R: DimMin<C>, C: Dim> Clone for LU<T, R, C>
where
DefaultAllocator: Allocator<T, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>,
PermutationSequence<DimMinimum<R, C>>: Clone,
Owned<T, R, C>: Clone,
{
fn clone(&self) -> Self {
Self {
lu: self.lu.clone(),
p: self.p.clone(),
}
}
}
impl<T: ComplexField, R: DimMin<C>, C: Dim> fmt::Debug for LU<T, R, C>
where
DefaultAllocator: Allocator<T, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>,
PermutationSequence<DimMinimum<R, C>>: fmt::Debug,
Owned<T, R, C>: fmt::Debug,
{
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
f.debug_struct("LU")
.field("lu", &self.lu)
.field("p", &self.p)
.finish()
}
}
/// Performs a LU decomposition to overwrite `out` with the inverse of `matrix`. /// Performs a LU decomposition to overwrite `out` with the inverse of `matrix`.
/// ///
/// If `matrix` is not invertible, `false` is returned and `out` may contain invalid data. /// If `matrix` is not invertible, `false` is returned and `out` may contain invalid data.

View File

@ -1,3 +1,4 @@
use std::fmt;
use std::mem::MaybeUninit; use std::mem::MaybeUninit;
#[cfg(feature = "serde-serialize-no-std")] #[cfg(feature = "serde-serialize-no-std")]
@ -11,7 +12,9 @@ use crate::base::{DefaultAllocator, Matrix, OVector, Scalar};
#[cfg(any(feature = "std", feature = "alloc"))] #[cfg(any(feature = "std", feature = "alloc"))]
use crate::dimension::Dynamic; use crate::dimension::Dynamic;
use crate::dimension::{Dim, DimName}; use crate::dimension::{Dim, DimName};
use crate::storage::StorageMut; use crate::iter::MatrixIter;
use crate::storage::{Owned, StorageMut};
use crate::{Const, U1};
/// A sequence of row or column permutations. /// A sequence of row or column permutations.
#[cfg_attr(feature = "serde-serialize-no-std", derive(Serialize, Deserialize))] #[cfg_attr(feature = "serde-serialize-no-std", derive(Serialize, Deserialize))]
@ -25,7 +28,6 @@ use crate::storage::StorageMut;
serde(bound(deserialize = "DefaultAllocator: Allocator<(usize, usize), D>, serde(bound(deserialize = "DefaultAllocator: Allocator<(usize, usize), D>,
OVector<(usize, usize), D>: Deserialize<'de>")) OVector<(usize, usize), D>: Deserialize<'de>"))
)] )]
#[derive(Clone, Debug)]
pub struct PermutationSequence<D: Dim> pub struct PermutationSequence<D: Dim>
where where
DefaultAllocator: Allocator<(usize, usize), D>, DefaultAllocator: Allocator<(usize, usize), D>,
@ -41,6 +43,32 @@ where
{ {
} }
impl<D: Dim> Clone for PermutationSequence<D>
where
DefaultAllocator: Allocator<(usize, usize), D>,
OVector<MaybeUninit<(usize, usize)>, D>: Clone,
{
fn clone(&self) -> Self {
Self {
len: self.len,
ipiv: self.ipiv.clone(),
}
}
}
impl<D: Dim> fmt::Debug for PermutationSequence<D>
where
DefaultAllocator: Allocator<(usize, usize), D>,
OVector<MaybeUninit<(usize, usize)>, D>: fmt::Debug,
{
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
f.debug_struct("PermutationSequence")
.field("len", &self.len)
.field("ipiv", &self.ipiv)
.finish()
}
}
impl<D: DimName> PermutationSequence<D> impl<D: DimName> PermutationSequence<D>
where where
DefaultAllocator: Allocator<(usize, usize), D>, DefaultAllocator: Allocator<(usize, usize), D>,
@ -74,7 +102,7 @@ where
unsafe { unsafe {
Self { Self {
len: 0, len: 0,
ipiv: OVector::new_uninitialized(dim), ipiv: OVector::new_uninitialized_generic(dim, Const::<1>),
} }
} }
} }
@ -88,7 +116,7 @@ where
self.len < self.ipiv.len(), self.len < self.ipiv.len(),
"Maximum number of permutations exceeded." "Maximum number of permutations exceeded."
); );
self.ipiv[self.len] = (i, i2); self.ipiv[self.len] = MaybeUninit::new((i, i2));
self.len += 1; self.len += 1;
} }
} }
@ -99,8 +127,8 @@ where
where where
S2: StorageMut<T, R2, C2>, S2: StorageMut<T, R2, C2>,
{ {
for i in self.ipiv.rows_range(..self.len).iter().map(MaybeUninit::assume_init) { for perm in self.iter() {
rhs.swap_rows(i.0, i.1) rhs.swap_rows(perm.0, perm.1)
} }
} }
@ -110,8 +138,8 @@ where
where where
S2: StorageMut<T, R2, C2>, S2: StorageMut<T, R2, C2>,
{ {
for i in 0..self.len { for perm in self.iter().rev() {
let (i1, i2) = self.ipiv[self.len - i - 1]; let (i1, i2) = perm;
rhs.swap_rows(i1, i2) rhs.swap_rows(i1, i2)
} }
} }
@ -122,8 +150,8 @@ where
where where
S2: StorageMut<T, R2, C2>, S2: StorageMut<T, R2, C2>,
{ {
for i in self.ipiv.rows_range(..self.len).iter() { for perm in self.iter() {
rhs.swap_columns(i.0, i.1) rhs.swap_columns(perm.0, perm.1)
} }
} }
@ -135,8 +163,8 @@ where
) where ) where
S2: StorageMut<T, R2, C2>, S2: StorageMut<T, R2, C2>,
{ {
for i in 0..self.len { for perm in self.iter().rev() {
let (i1, i2) = self.ipiv[self.len - i - 1]; let (i1, i2) = perm;
rhs.swap_columns(i1, i2) rhs.swap_columns(i1, i2)
} }
} }
@ -163,4 +191,27 @@ where
-T::one() -T::one()
} }
} }
/// Iterates over the permutations that have been initialized.
pub fn iter(
&self,
) -> std::iter::Map<
std::iter::Copied<
std::iter::Take<
MatrixIter<
MaybeUninit<(usize, usize)>,
D,
U1,
Owned<MaybeUninit<(usize, usize)>, D, U1>,
>,
>,
>,
impl FnMut(MaybeUninit<(usize, usize)>) -> (usize, usize),
> {
self.ipiv
.iter()
.take(self.len)
.copied()
.map(|e| unsafe { e.assume_init() })
}
} }