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

Almost!

This commit is contained in:
Violeta Hernández 2021-07-17 04:36:14 -05:00
parent a6b8dd6d78
commit 9a528e23b9
24 changed files with 423 additions and 145 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
nalgebra-lapack/src/cholesky.rs
View File

@ -24,17 +24,17 @@ use lapack;
OMatrix<T, D, D>: Deserialize<'de>"))
)]
#[derive(Clone, Debug)]
pub struct Cholesky<T: Scalar, D: Dim>
pub struct Cholesky<T, D: Dim>
where
DefaultAllocator: Allocator<T, D, D>,
{
l: OMatrix<T, D, D>,
}
impl<T: Scalar + Copy, D: Dim> Copy for Cholesky<T, D>
impl<T: Copy, D: Dim> Copy for Cholesky<T, D>
where
DefaultAllocator: Allocator<T, D, D>,
OMatrix<T, D, D>: Copy,
Owned<T, D, D>: Copy,
{
}

92
src/base/blas.rs
View File

@ -329,22 +329,22 @@ where
if !b.is_zero() {
for i in 0..x.len() {
unsafe {
let y = y.get_unchecked_mut(i * rstride1);
*y = a.inlined_clone()
* x.get_unchecked(i * rstride2).inlined_clone()
* c.inlined_clone()
+ b.inlined_clone() * y.inlined_clone();
}
}
} else {
for i in 0..x.len() {
unsafe {
let y = y.get_unchecked_mut(i * rstride1);
*y = a.inlined_clone()
* x.get_unchecked(i * rstride2).inlined_clone()
* c.inlined_clone();
}
}
}
}
@ -788,11 +788,83 @@ where
for j in 1..ncols2 {
let col2 = a.column(j);
let val = unsafe { x.vget_unchecked(j).inlined_clone() };
let val = x.vget_unchecked(j).inlined_clone() ;
init.axcpy(alpha.inlined_clone(), &col2, val, T::one());
}
}
}
#[inline(always)]
fn xxgemv_z<D2: Dim, D3: Dim, SB, SC>(
&mut self,
alpha: T,
a: &SquareMatrix<T, D2, SB>,
x: &Vector<T, D3, SC>,
dot: impl Fn(
&DVectorSlice<T, SB::RStride, SB::CStride>,
&DVectorSlice<T, SC::RStride, SC::CStride>,
) -> T,
) where
T: One,
SB: Storage<T, D2, D2>,
SC: Storage<T, D3>,
ShapeConstraint: DimEq<D, D2> + AreMultipliable<D2, D2, D3, U1>,
{
let dim1 = self.nrows();
let dim2 = a.nrows();
let dim3 = x.nrows();
assert!(
a.is_square(),
"Symmetric cgemv: the input matrix must be square."
);
assert!(
dim2 == dim3 && dim1 == dim2,
"Symmetric cgemv: dimensions mismatch."
);
if dim2 == 0 {
return;
}
// TODO: avoid bound checks.
let col2 = a.column(0);
let val = unsafe { x.vget_unchecked(0).inlined_clone() };
self.axc(alpha.inlined_clone(), &col2, val);
let mut res = unsafe { self.assume_init_mut() };
res[0] += alpha.inlined_clone() * dot(&a.slice_range(1.., 0), &x.rows_range(1..));
for j in 1..dim2 {
let col2 = a.column(j);
let dot = dot(&col2.rows_range(j..), &x.rows_range(j..));
let val;
unsafe {
val = x.vget_unchecked(j).inlined_clone();
*res.vget_unchecked_mut(j) += alpha.inlined_clone() * dot;
}
res.rows_range_mut(j + 1..).axpy(
alpha.inlined_clone() * val,
&col2.rows_range(j + 1..),
T::one(),
);
}
}
pub fn hegemv_z<D2: Dim, D3: Dim, SB, SC>(
&mut self,
alpha: T,
a: &SquareMatrix<T, D2, SB>,
x: &Vector<T, D3, SC>,
) where
T: SimdComplexField,
SB: Storage<T, D2, D2>,
SC: Storage<T, D3>,
ShapeConstraint: DimEq<D, D2> + AreMultipliable<D2, D2, D3, U1>,
{
self.xxgemv_z(alpha, a, x, |a, b| a.dotc(b))
}
}
impl<T, R1: Dim, C1: Dim, S: StorageMut<MaybeUninit<T>, R1, C1>> Matrix<MaybeUninit<T>, R1, C1, S>
@ -850,7 +922,7 @@ where
// matrixmultiply can be used only if the std feature is available.
let nrows1 = self.nrows();
let (nrows2, ncols2) = a.shape();
let (nrows3, ncols3) = b.shape();
let (_, ncols3) = b.shape();
// Threshold determined empirically.
const SMALL_DIM: usize = 5;
@ -1502,9 +1574,9 @@ where
ShapeConstraint: DimEq<D1, R3> + DimEq<C3, D4>,
DefaultAllocator: Allocator<T, R3>,
{
let work = Matrix::new_uninitialized_generic(R3::from_usize(self.shape().0), Const::<1>);
let mut work = Matrix::new_uninitialized_generic(R3::from_usize(self.shape().0), Const::<1>);
work.gemv_z(T::one(), lhs, &mid.column(0));
let work = unsafe { work.assume_init() };
let mut work = unsafe { work.assume_init() };
self.ger(alpha.inlined_clone(), &work, &lhs.column(0), beta);
@ -1552,9 +1624,9 @@ where
DefaultAllocator: Allocator<T, D3>,
{
// TODO: figure out why type inference wasn't doing its job.
let work = Matrix::new_uninitialized_generic(D3::from_usize(self.shape().0), Const::<1>);
let mut work = Matrix::new_uninitialized_generic(D3::from_usize(self.shape().0), Const::<1>);
work.gemv_z::<D3, D3, R4, S3, _>(T::one(), mid, &rhs.column(0));
let work = unsafe { work.assume_init() };
let mut work = unsafe { work.assume_init() };
self.column_mut(0)
.gemv_tr(alpha.inlined_clone(), rhs, &work, beta.inlined_clone());

9
src/base/conversion.rs
View File

@ -109,13 +109,14 @@ impl<T, const D: usize> From<[T; D]> for SVector<T, D> {
}
}
impl<T, const D: usize> From<SVector<T, D>> for [T; D] {
impl<T, const D: usize> From<SVector<T, D>> for [T; D]
where
T: Clone,
{
#[inline]
fn from(vec: SVector<T, D>) -> Self {
// TODO: unfortunately, we must clone because we can move out of an array.
// Counterpoint: this seems to work?
vec.data.0[0]
vec.data.0[0].clone()
}
}

14
src/base/default_allocator.rs
View File

@ -72,7 +72,9 @@ impl<T, const R: usize, const C: usize> Allocator<T, Const<R>, Const<C>> for Def
_: Const<R>,
_: Const<C>,
) -> Owned<MaybeUninit<T>, Const<R>, Const<C>> {
ArrayStorage([[MaybeUninit::uninit(); R]; C])
// SAFETY: An uninitialized `[MaybeUninit<_>; LEN]` is valid.
let array = unsafe { MaybeUninit::uninit().assume_init() };
ArrayStorage(array)
}
#[inline]
@ -126,9 +128,8 @@ impl<T, C: Dim> Allocator<T, Dynamic, C> for DefaultAllocator {
let mut data = ManuallyDrop::new(uninit.data);
// Safety: MaybeUninit<T> has the same alignment and layout as T.
let new_data = unsafe {
Vec::from_raw_parts(data.as_mut_ptr() as *mut T, data.len(), data.capacity())
};
let new_data =
Vec::from_raw_parts(data.as_mut_ptr() as *mut T, data.len(), data.capacity());
VecStorage::new(uninit.nrows, uninit.ncols, new_data)
}
@ -170,9 +171,8 @@ impl<T, R: DimName> Allocator<T, R, Dynamic> for DefaultAllocator {
let mut data = ManuallyDrop::new(uninit.data);
// Safety: MaybeUninit<T> has the same alignment and layout as T.
let new_data = unsafe {
Vec::from_raw_parts(data.as_mut_ptr() as *mut T, data.len(), data.capacity())
};
let new_data =
Vec::from_raw_parts(data.as_mut_ptr() as *mut T, data.len(), data.capacity());
VecStorage::new(uninit.nrows, uninit.ncols, new_data)
}

21
src/base/edition.rs
View File

@ -178,7 +178,7 @@ impl<T, R: Dim, C: Dim, S: StorageMut<T, R, C>> Matrix<T, R, C, S> {
/// Sets all the elements of this matrix to `f()`.
#[inline]
pub fn fill_fn<F: FnMut() -> T>(&mut self, f: F) {
pub fn fill_fn<F: FnMut() -> T>(&mut self, mut f: F) {
for e in self.iter_mut() {
*e = f();
}
@ -942,8 +942,11 @@ impl<T: Clone> OMatrix<T, Dynamic, Dynamic> {
where
DefaultAllocator: Reallocator<T, Dynamic, Dynamic, Dynamic, Dynamic>,
{
let placeholder =
Matrix::new_uninitialized_generic(Dynamic::new(0), Dynamic::new(0)).assume_init();
// BEEEP!!!! BEEEEEEEP!!!
let placeholder = unsafe {
Matrix::new_uninitialized_generic(Dynamic::new(0), Dynamic::new(0)).assume_init()
};
let old = mem::replace(self, placeholder);
let new = old.resize(new_nrows, new_ncols, val);
let _ = mem::replace(self, new);
@ -966,7 +969,8 @@ where
where
DefaultAllocator: Reallocator<T, Dynamic, C, Dynamic, C>,
{
let placeholder = Matrix::from_fn_generic(Dynamic::new(0), self.data.shape().1, |_, _| val);
let placeholder =
Matrix::from_fn_generic(Dynamic::new(0), self.data.shape().1, |_, _| val.clone());
let old = mem::replace(self, placeholder);
let new = old.resize_vertically(new_nrows, val);
let _ = mem::replace(self, new);
@ -989,7 +993,8 @@ where
where
DefaultAllocator: Reallocator<T, R, Dynamic, R, Dynamic>,
{
let placeholder = Matrix::from_fn_generic(self.data.shape().0, Dynamic::new(0), |_, _| val);
let placeholder =
Matrix::from_fn_generic(self.data.shape().0, Dynamic::new(0), |_, _| val.clone());
let old = mem::replace(self, placeholder);
let new = old.resize_horizontally(new_ncols, val);
let _ = mem::replace(self, new);
@ -1059,11 +1064,7 @@ unsafe fn extend_rows<T>(data: &mut [T], nrows: usize, ncols: usize, i: usize, n
/// Extend the number of columns of the `Matrix` with elements from
/// a given iterator.
#[cfg(any(feature = "std", feature = "alloc"))]
impl<T, R, S> Extend<T> for Matrix<T, R, Dynamic, S>
where
R: Dim,
S: Extend<T>,
{
impl<T, R: Dim, S: Extend<T>> Extend<T> for Matrix<T, R, Dynamic, S> {
/// Extend the number of columns of the `Matrix` with elements
/// from the given iterator.
///

8
src/base/matrix.rs
View File

@ -1249,7 +1249,7 @@ impl<T, R: Dim, C: Dim, S: StorageMut<T, R, C>> Matrix<T, R, C, S> {
/// Fills this matrix with the content of another one, after applying a function to
/// the references of the entries of the other matrix. Both must have the same shape.
#[inline]
pub fn copy_from_fn<U, R2: Dim, C2: Dim, SB, F>(&mut self, other: &Matrix<U, R2, C2, SB>, f: F)
pub fn copy_from_fn<U, R2: Dim, C2: Dim, SB, F>(&mut self, other: &Matrix<U, R2, C2, SB>,mut f: F)
where
SB: Storage<U, R2, C2>,
ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>,
@ -1282,7 +1282,7 @@ impl<T, R: Dim, C: Dim, S: StorageMut<T, R, C>> Matrix<T, R, C, S> {
/// Fills this matrix with the content of another one via moves. Both must have the same shape.
#[inline]
pub fn move_from_fn<U, R2: Dim, C2: Dim, SB, F>(&mut self, other: Matrix<U, R2, C2, SB>, f: F)
pub fn move_from_fn<U, R2: Dim, C2: Dim, SB, F>(&mut self, other: Matrix<U, R2, C2, SB>, mut f: F)
where
SB: Storage<U, R2, C2>,
ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>,
@ -1322,7 +1322,7 @@ impl<T, R: Dim, C: Dim, S: StorageMut<T, R, C>> Matrix<T, R, C, S> {
pub fn tr_copy_from_fn<U, R2: Dim, C2: Dim, SB, F>(
&mut self,
other: &Matrix<U, R2, C2, SB>,
f: F,
mut f: F,
) where
SB: Storage<U, R2, C2>,
ShapeConstraint: DimEq<R, C2> + SameNumberOfColumns<C, R2>,
@ -1359,7 +1359,7 @@ impl<T, R: Dim, C: Dim, S: StorageMut<T, R, C>> Matrix<T, R, C, S> {
pub fn tr_move_from_fn<U, R2: Dim, C2: Dim, SB, F>(
&mut self,
other: Matrix<U, R2, C2, SB>,
f: F,
mut f: F,
) where
SB: Storage<U, R2, C2>,
ShapeConstraint: DimEq<R, C2> + SameNumberOfColumns<C, R2>,

3
src/base/statistics.rs
View File

@ -59,11 +59,12 @@ impl<T, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
}
/// Returns a column vector resulting from the folding of `f` on each column of this matrix.
// BEEEEP!!!! Pretty sure there's something fishy here.
#[inline]
#[must_use]
pub fn compress_columns(
&self,
init: OVector<T, R>,
mut init: OVector<T, R>,
f: impl Fn(&mut OVector<T, R>, VectorSlice<T, R, S::RStride, S::CStride>),
) -> OVector<T, R>
where

3
src/geometry/dual_quaternion.rs
View File

@ -278,7 +278,8 @@ impl<'a, T: Deserialize<'a>> Deserialize<'a> for DualQuaternion<T> {
}
impl<T> DualQuaternion<T> {
fn to_vector(self) -> OVector<T, U8> {
// TODO: Cloning shouldn't be necessary.
fn to_vector(self) -> OVector<T, U8>where T:Clone {
(*self.as_ref()).into()
}
}

4
src/geometry/orthographic.rs
View File

@ -28,7 +28,9 @@ impl<T: Copy> Copy for Orthographic3<T> {}
impl<T: Clone> Clone for Orthographic3<T> {
#[inline]
fn clone(&self) -> Self {
Self::from_matrix_unchecked(self.matrix)
Self {
matrix: self.matrix.clone(),
}
}
}

11
src/geometry/point.rs
View File

@ -215,7 +215,7 @@ where
let mut res = OVector::<_, DimNameSum<D, U1>>::new_uninitialized();
for i in 0..D::dim() {
unsafe {
*res.get_unchecked(i) = MaybeUninit::new(self.coords[i].clone());
*res.get_unchecked_mut(i) = MaybeUninit::new(self.coords[i].clone());
}
}
@ -236,15 +236,16 @@ where
// to avoid double-dropping.
for i in 0..D::dim() {
unsafe {
*res.get_unchecked(i) = MaybeUninit::new(self.coords[i]);
*res.get_unchecked_mut(i) = MaybeUninit::new(*self.coords.get_unchecked(i));
}
}
// Fix double drop
res[(D::dim(), 0)] = MaybeUninit::new(T::one());
unsafe { res.assume_init() }
unsafe {
*res.get_unchecked_mut(D::dim()) = MaybeUninit::new(T::one());
res.assume_init()
}
}
/// Creates a new point with the given coordinates.

3
src/geometry/point_conversion.rs
View File

@ -91,7 +91,8 @@ impl<T, const D: usize> From<[T; D]> for Point<T, D> {
}
}
impl<T, const D: usize> From<Point<T, D>> for [T; D] {
impl<T, const D: usize> From<Point<T, D>> for [T; D] where
T: Clone,{
#[inline]
fn from(p: Point<T, D>) -> Self {
p.coords.into()

4
src/geometry/transform.rs
View File

@ -550,8 +550,8 @@ where
Const<D>: DimNameAdd<U1>,
C: SubTCategoryOf<TProjective>,
DefaultAllocator: Allocator<T, DimNameSum<Const<D>, U1>, DimNameSum<Const<D>, U1>>
+ Allocator<T, DimNameSum<Const<D>, U1>>, // + Allocator<T, D, D>
// + Allocator<T, D>
+ Allocator<T, DimNameSum<Const<D>, U1>>,
Owned<T, DimNameSum<Const<D>, U1>, DimNameSum<Const<D>, U1>>: Clone,
{
/// Transform the given point by the inverse of this transformation.
/// This may be cheaper than inverting the transformation and transforming

5
src/geometry/transform_ops.rs
View File

@ -8,7 +8,7 @@ use simba::scalar::{ClosedAdd, ClosedMul, RealField, SubsetOf};
use crate::base::allocator::Allocator;
use crate::base::dimension::{DimNameAdd, DimNameSum, U1};
use crate::base::{Const, DefaultAllocator, OMatrix, SVector, Scalar};
use crate::base::{Const, DefaultAllocator, OMatrix, SVector, Scalar};use crate::storage::Owned;
use crate::geometry::{
Isometry, Point, Rotation, Similarity, SubTCategoryOf, SuperTCategoryOf, TAffine, TCategory,
@ -586,7 +586,8 @@ md_assign_impl_all!(
const D;
for CA, CB;
where Const<D>: DimNameAdd<U1>, CA: SuperTCategoryOf<CB>, CB: SubTCategoryOf<TProjective>,
DefaultAllocator: Allocator<T, DimNameSum<Const<D>, U1>, DimNameSum<Const<D>, U1>>;
DefaultAllocator: Allocator<T, DimNameSum<Const<D>, U1>, DimNameSum<Const<D>, U1>>,
Owned<T, DimNameSum<Const<D>, U1>, DimNameSum<Const<D>, U1>>: Clone;
self: Transform<T, CA, D>, rhs: Transform<T, CB, D>;
[val] => #[allow(clippy::suspicious_op_assign_impl)] { *self *= rhs.inverse() };
[ref] => #[allow(clippy::suspicious_op_assign_impl)] { *self *= rhs.clone().inverse() };

5
src/geometry/translation_conversion.rs
View File

@ -215,7 +215,10 @@ impl<T, const D: usize> From<Point<T, D>> for Translation<T, D> {
}
}
impl<T, const D: usize> From<Translation<T, D>> for [T; D] {
impl<T, const D: usize> From<Translation<T, D>> for [T; D]
where
T: Clone,
{
#[inline]
fn from(t: Translation<T, D>) -> Self {
t.vector.into()

4
src/linalg/bidiagonal.rs
View File

@ -185,6 +185,7 @@ where
);
}
unsafe {
Bidiagonal {
uv: matrix,
diagonal: diagonal.assume_init(),
@ -192,6 +193,7 @@ where
upper_diagonal,
}
}
}
/// Indicates whether this decomposition contains an upper-diagonal matrix.
#[inline]
@ -300,7 +302,7 @@ where
let axis = self.uv.slice_range(i, i + shift..);
let mut axis_packed = axis_packed.rows_range_mut(i + shift..);
axis_packed.tr_copy_init_from(&axis);
let mut axis_packed = unsafe { axis_packed.slice_assume_init() };
let axis_packed = unsafe { axis_packed.slice_assume_init() };
// TODO: sometimes, the axis might have a zero magnitude.
let refl = Reflection::new(Unit::new_unchecked(axis_packed), T::zero());

4
src/linalg/cholesky.rs
View File

@ -307,7 +307,7 @@ where
);
chol.slice_range_mut(j + 1.., j).copy_init_from(&new_colj);
let chol = unsafe { chol.assume_init() };
let mut chol = unsafe { chol.assume_init() };
// update the bottom right corner
let mut bottom_right_corner = chol.slice_range_mut(j + 1.., j + 1..);
@ -348,7 +348,7 @@ where
.copy_init_from(&self.chol.slice_range(j + 1.., ..j));
chol.slice_range_mut(j.., j..)
.copy_init_from(&self.chol.slice_range(j + 1.., j + 1..));
let chol = unsafe { chol.assume_init() };
let mut chol = unsafe { chol.assume_init() };
// updates the bottom right corner
let mut bottom_right_corner = chol.slice_range_mut(j.., j..);

4
src/linalg/permutation_sequence.rs
View File

@ -99,11 +99,11 @@ where
/// Creates a new sequence of D identity permutations.
#[inline]
pub fn identity_generic(dim: D) -> Self {
unsafe {
Self {
len: 0,
ipiv: OVector::new_uninitialized_generic(dim, Const::<1>),
}
}
}

8
src/linalg/pow.rs
View File

@ -40,18 +40,24 @@ where
// We use the buffer to hold the result of multiplier ^ 2, thus avoiding
// extra allocations.
let (nrows, ncols) = self.data.shape();
let mut multiplier = self.clone_owned();
let mut buf = self.clone_owned();
// TODO: ACTUALLY MAKE BUF USEFUL! BEEEEEEEEP!!
// Exponentiation by squares.
loop {
if e % two == one {
let mut buf = Matrix::new_uninitialized_generic(nrows, ncols);
self.mul_to(&multiplier, &mut buf);
let buf = unsafe { buf.assume_init() };
self.copy_from(&buf);
}
e /= two;
let mut buf = Matrix::new_uninitialized_generic(nrows, ncols);
multiplier.mul_to(&multiplier, &mut buf);
let buf = unsafe { buf.assume_init() };
multiplier.copy_from(&buf);
if e == zero {

57
src/linalg/qr.rs
View File

@ -1,3 +1,5 @@
use std::fmt;
use num::Zero;
#[cfg(feature = "serde-serialize-no-std")]
use serde::{Deserialize, Serialize};
@ -6,7 +8,7 @@ use crate::allocator::{Allocator, Reallocator};
use crate::base::{DefaultAllocator, Matrix, OMatrix, OVector, Unit};
use crate::constraint::{SameNumberOfRows, ShapeConstraint};
use crate::dimension::{Const, Dim, DimMin, DimMinimum};
use crate::storage::{Storage, StorageMut};
use crate::storage::{Owned, Storage, StorageMut};
use simba::scalar::ComplexField;
use crate::geometry::Reflection;
@ -28,8 +30,8 @@ use crate::linalg::householder;
OMatrix<T, R, C>: Deserialize<'de>,
OVector<T, DimMinimum<R, C>>: Deserialize<'de>"))
)]
#[derive(Clone, Debug)]
pub struct QR<T: ComplexField, R: DimMin<C>, C: Dim>
pub struct QR<T, R: DimMin<C>, C: Dim>
where
DefaultAllocator: Allocator<T, R, C> + Allocator<T, DimMinimum<R, C>>,
{
@ -37,14 +39,42 @@ where
diag: OVector<T, DimMinimum<R, C>>,
}
impl<T: ComplexField, R: DimMin<C>, C: Dim> Copy for QR<T, R, C>
impl<T: Copy, R: DimMin<C>, C: Dim> Copy for QR<T, R, C>
where
DefaultAllocator: Allocator<T, R, C> + Allocator<T, DimMinimum<R, C>>,
OMatrix<T, R, C>: Copy,
OVector<T, DimMinimum<R, C>>: Copy,
Owned<T, R, C>: Copy,
Owned<T, DimMinimum<R, C>>: Copy,
{
}
impl<T: Clone, R: DimMin<C>, C: Dim> Clone for QR<T, R, C>
where
DefaultAllocator: Allocator<T, R, C> + Allocator<T, DimMinimum<R, C>>,
Owned<T, R, C>: Clone,
Owned<T, DimMinimum<R, C>>: Clone,
{
fn clone(&self) -> Self {
Self {
qr: self.qr.clone(),
diag: self.diag.clone(),
}
}
}
impl<T: fmt::Debug, R: DimMin<C>, C: Dim> fmt::Debug for QR<T, R, C>
where
DefaultAllocator: Allocator<T, R, C> + Allocator<T, DimMinimum<R, C>>,
Owned<T, R, C>: fmt::Debug,
Owned<T, DimMinimum<R, C>>: fmt::Debug,
{
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
f.debug_struct("QR")
.field("qr", &self.qr)
.field("diag", &self.diag)
.finish()
}
}
impl<T: ComplexField, R: DimMin<C>, C: Dim> QR<T, R, C>
where
DefaultAllocator: Allocator<T, R, C> + Allocator<T, R> + Allocator<T, DimMinimum<R, C>>,
@ -54,18 +84,23 @@ where
let (nrows, ncols) = matrix.data.shape();
let min_nrows_ncols = nrows.min(ncols);
let mut diag =
unsafe { crate::unimplemented_or_uninitialized_generic!(min_nrows_ncols, Const::<1>) };
let mut diag = Matrix::new_uninitialized_generic(min_nrows_ncols, Const::<1>);
if min_nrows_ncols.value() == 0 {
return QR { qr: matrix, diag };
return Self {
qr: matrix,
diag: unsafe { diag.assume_init() },
};
}
for i in 0..min_nrows_ncols.value() {
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);
}
QR { qr: matrix, diag }
Self {
qr: matrix,
diag: unsafe { diag.assume_init() },
}
}
/// Retrieves the upper trapezoidal submatrix `R` of this decomposition.

92
src/linalg/schur.rs
View File

@ -1,16 +1,18 @@
#![allow(clippy::suspicious_operation_groupings)]
use std::cmp;
use std::fmt;
use std::mem::MaybeUninit;
#[cfg(feature = "serde-serialize-no-std")]
use serde::{Deserialize, Serialize};
use approx::AbsDiffEq;
use num_complex::Complex as NumComplex;
use simba::scalar::{ComplexField, RealField};
use std::cmp;
use std::mem::MaybeUninit;
use crate::allocator::Allocator;
use crate::base::dimension::{Const, Dim, DimDiff, DimSub, Dynamic, U1, U2};
use crate::base::storage::Storage;
use crate::base::storage::{Owned, Storage};
use crate::base::{DefaultAllocator, OMatrix, OVector, SquareMatrix, Unit, Vector2, Vector3};
use crate::geometry::Reflection;
@ -32,8 +34,7 @@ use crate::linalg::Hessenberg;
serde(bound(deserialize = "DefaultAllocator: Allocator<T, D, D>,
OMatrix<T, D, D>: Deserialize<'de>"))
)]
#[derive(Clone, Debug)]
pub struct Schur<T: ComplexField, D: Dim>
pub struct Schur<T, D: Dim>
where
DefaultAllocator: Allocator<T, D, D>,
{
@ -41,13 +42,39 @@ where
t: OMatrix<T, D, D>,
}
impl<T: ComplexField, D: Dim> Copy for Schur<T, D>
impl<T: Copy, D: Dim> Copy for Schur<T, D>
where
DefaultAllocator: Allocator<T, D, D>,
OMatrix<T, D, D>: Copy,
Owned<T, D, D>: Copy,
{
}
impl<T: Clone, D: Dim> Clone for Schur<T, D>
where
DefaultAllocator: Allocator<T, D, D>,
Owned<T, D, D>: Clone,
{
fn clone(&self) -> Self {
Self {
q: self.q.clone(),
t: self.t.clone(),
}
}
}
impl<T: fmt::Debug, D: Dim> fmt::Debug for Schur<T, D>
where
DefaultAllocator: Allocator<T, D, D>,
Owned<T, D, D>: fmt::Debug,
{
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
f.debug_struct("Schur")
.field("q", &self.q)
.field("t", &self.t)
.finish()
}
}
impl<T: ComplexField, D: Dim> Schur<T, D>
where
D: DimSub<U1>, // For Hessenberg.
@ -73,8 +100,7 @@ where
/// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm
/// continues indefinitely until convergence.
pub fn try_new(m: OMatrix<T, D, D>, eps: T::RealField, max_niter: usize) -> Option<Self> {
let mut work =
unsafe { crate::unimplemented_or_uninitialized_generic!(m.data.shape().0, Const::<1>) };
let mut work = OVector::new_uninitialized_generic(m.data.shape().0, Const::<1>);
Self::do_decompose(m, &mut work, eps, max_niter, true)
.map(|(q, t)| Schur { q: q.unwrap(), t })
@ -82,7 +108,7 @@ where
fn do_decompose(
mut m: OMatrix<T, D, D>,
work: &mut OVector<T, D>,
work: &mut OVector<MaybeUninit<T>, D>,
eps: T::RealField,
max_niter: usize,
compute_q: bool,
@ -271,7 +297,9 @@ where
}
/// Computes the eigenvalues of the decomposed matrix.
fn do_eigenvalues(t: &OMatrix<T, D, D>, out: &mut OVector<T, D>) -> bool {
fn do_eigenvalues(t: &OMatrix<T, D, D>, out: &mut OVector<MaybeUninit<T>, D>) -> bool {
// TODO: check dropping stuff.
let dim = t.nrows();
let mut m = 0;
@ -279,7 +307,7 @@ where
let n = m + 1;
if t[(n, m)].is_zero() {
out[m] = t[(m, m)];
out[m] = MaybeUninit::new(t[(m, m)]);
m += 1;
} else {
// Complex eigenvalue.
@ -288,7 +316,7 @@ where
}
if m == dim - 1 {
out[m] = t[(m, m)];
out[m] = MaybeUninit::new(t[(m, m)]);
}
true
@ -297,11 +325,13 @@ where
/// Computes the complex eigenvalues of the decomposed matrix.
fn do_complex_eigenvalues(
t: &OMatrix<T, D, D>,
out: &mut OVector<NumComplex<T>, D>,
out: &mut OVector<MaybeUninit<NumComplex<T>>, D>,
) where
T: RealField,
DefaultAllocator: Allocator<NumComplex<T>, D>,
{
// TODO: check for dropping behavior.
let dim = t.nrows();
let mut m = 0;
@ -309,7 +339,7 @@ where
let n = m + 1;
if t[(n, m)].is_zero() {
out[m] = NumComplex::new(t[(m, m)], T::zero());
out[m] = MaybeUninit::new(NumComplex::new(t[(m, m)], T::zero()));
m += 1;
} else {
// Solve the 2x2 eigenvalue subproblem.
@ -391,11 +421,9 @@ where
/// Return `None` if some eigenvalues are complex.
#[must_use]
pub fn eigenvalues(&self) -> Option<OVector<T, D>> {
let mut out = unsafe {
crate::unimplemented_or_uninitialized_generic!(self.t.data.shape().0, Const::<1>)
};
let mut out = OVector::new_uninitialized_generic(self.t.data.shape().0, Const::<1>);
if Self::do_eigenvalues(&self.t, &mut out) {
Some(out)
Some(unsafe { out.assume_init() })
} else {
None
}
@ -408,11 +436,9 @@ where
T: RealField,
DefaultAllocator: Allocator<NumComplex<T>, D>,
{
let mut out = unsafe {
crate::unimplemented_or_uninitialized_generic!(self.t.data.shape().0, Const::<1>)
};
let mut out = OVector::new_uninitialized_generic(self.t.data.shape().0, Const::<1>);
Self::do_complex_eigenvalues(&self.t, &mut out);
out
unsafe { out.assume_init() }
}
}
@ -517,14 +543,14 @@ where
/// Computes the eigenvalues of this matrix.
#[must_use]
pub fn eigenvalues(&self) -> Option<OVector<T, D>> {
// TODO: check drop stuff.
assert!(
self.is_square(),
"Unable to compute eigenvalues of a non-square matrix."
);
let mut work = unsafe {
crate::unimplemented_or_uninitialized_generic!(self.data.shape().0, Const::<1>)
};
let mut work = OVector::new_uninitialized_generic(self.data.shape().0, Const::<1>);
// Special case for 2x2 matrices.
if self.nrows() == 2 {
@ -533,9 +559,9 @@ where
let me = self.fixed_slice::<2, 2>(0, 0);
return match compute_2x2_eigvals(&me) {
Some((a, b)) => {
work[0] = a;
work[1] = b;
Some(work)
work[0] = MaybeUninit::new(a);
work[1] = MaybeUninit::new(b);
Some(unsafe { work.assume_init() })
}
None => None,
};
@ -551,7 +577,7 @@ where
)
.unwrap();
if Schur::do_eigenvalues(&schur.1, &mut work) {
Some(work)
Some(unsafe { work.assume_init() })
} else {
None
}
@ -566,7 +592,7 @@ where
DefaultAllocator: Allocator<NumComplex<T>, D>,
{
let dim = self.data.shape().0;
let mut work = unsafe { crate::unimplemented_or_uninitialized_generic!(dim, Const::<1>) };
let mut work = OVector::new_uninitialized_generic(dim, Const::<1>);
let schur = Schur::do_decompose(
self.clone_owned(),
@ -576,8 +602,8 @@ where
false,
)
.unwrap();
let mut eig = unsafe { crate::unimplemented_or_uninitialized_generic!(dim, Const::<1>) };
let mut eig = OVector::new_uninitialized_generic(dim, Const::<1>);
Schur::do_complex_eigenvalues(&schur.1, &mut eig);
eig
unsafe { eig.assume_init() }
}
}

47
src/linalg/svd.rs
View File

@ -1,3 +1,5 @@
use std::fmt;
#[cfg(feature = "serde-serialize-no-std")]
use serde::{Deserialize, Serialize};
@ -8,7 +10,7 @@ use crate::allocator::Allocator;
use crate::base::{DefaultAllocator, Matrix, Matrix2x3, OMatrix, OVector, Vector2};
use crate::constraint::{SameNumberOfRows, ShapeConstraint};
use crate::dimension::{Dim, DimDiff, DimMin, DimMinimum, DimSub, U1};
use crate::storage::Storage;
use crate::storage::{Owned, Storage};
use simba::scalar::{ComplexField, RealField};
use crate::linalg::givens::GivensRotation;
@ -39,7 +41,6 @@ use crate::linalg::Bidiagonal;
OVector<T::RealField, DimMinimum<R, C>>: Deserialize<'de>"
))
)]
#[derive(Clone, Debug)]
pub struct SVD<T: ComplexField, R: DimMin<C>, C: Dim>
where
DefaultAllocator: Allocator<T, DimMinimum<R, C>, C>
@ -59,12 +60,48 @@ where
DefaultAllocator: Allocator<T, DimMinimum<R, C>, C>
+ Allocator<T, R, DimMinimum<R, C>>
+ Allocator<T::RealField, DimMinimum<R, C>>,
OMatrix<T, R, DimMinimum<R, C>>: Copy,
OMatrix<T, DimMinimum<R, C>, C>: Copy,
OVector<T::RealField, DimMinimum<R, C>>: Copy,
Owned<T, R, DimMinimum<R, C>>: Copy,
Owned<T, DimMinimum<R, C>, C>: Copy,
Owned<T::RealField, DimMinimum<R, C>>: Copy,
{
}
impl<T: ComplexField, R: DimMin<C>, C: Dim> Clone for SVD<T, R, C>
where
DefaultAllocator: Allocator<T, DimMinimum<R, C>, C>
+ Allocator<T, R, DimMinimum<R, C>>
+ Allocator<T::RealField, DimMinimum<R, C>>,
Owned<T, R, DimMinimum<R, C>>: Clone,
Owned<T, DimMinimum<R, C>, C>: Clone,
Owned<T::RealField, DimMinimum<R, C>>: Clone,
{
fn clone(&self) -> Self {
Self {
u: self.u.clone(),
v_t: self.v_t.clone(),
singular_values: self.singular_values.clone(),
}
}
}
impl<T: ComplexField, R: DimMin<C>, C: Dim> fmt::Debug for SVD<T, R, C>
where
DefaultAllocator: Allocator<T, DimMinimum<R, C>, C>
+ Allocator<T, R, DimMinimum<R, C>>
+ Allocator<T::RealField, DimMinimum<R, C>>,
Owned<T, R, DimMinimum<R, C>>: fmt::Debug,
Owned<T, DimMinimum<R, C>, C>: fmt::Debug,
Owned<T::RealField, DimMinimum<R, C>>: fmt::Debug,
{
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
f.debug_struct("SVD")
.field("u", &self.u)
.field("v_t", &self.v_t)
.field("singular_values", &self.singular_values)
.finish()
}
}
impl<T: ComplexField, R: DimMin<C>, C: Dim> SVD<T, R, C>
where
DimMinimum<R, C>: DimSub<U1>, // for Bidiagonal.

42
src/linalg/symmetric_eigen.rs
View File

@ -1,3 +1,5 @@
use std::fmt;
#[cfg(feature = "serde-serialize-no-std")]
use serde::{Deserialize, Serialize};
@ -7,7 +9,7 @@ use num::Zero;
use crate::allocator::Allocator;
use crate::base::{DefaultAllocator, Matrix2, OMatrix, OVector, SquareMatrix, Vector2};
use crate::dimension::{Dim, DimDiff, DimSub, U1};
use crate::storage::Storage;
use crate::storage::{Owned, Storage};
use simba::scalar::ComplexField;
use crate::linalg::givens::GivensRotation;
@ -29,7 +31,6 @@ use crate::linalg::SymmetricTridiagonal;
OVector<T::RealField, D>: Deserialize<'de>,
OMatrix<T, D, D>: Deserialize<'de>"))
)]
#[derive(Clone, Debug)]
pub struct SymmetricEigen<T: ComplexField, D: Dim>
where
DefaultAllocator: Allocator<T, D, D> + Allocator<T::RealField, D>,
@ -44,11 +45,39 @@ where
impl<T: ComplexField, D: Dim> Copy for SymmetricEigen<T, D>
where
DefaultAllocator: Allocator<T, D, D> + Allocator<T::RealField, D>,
OMatrix<T, D, D>: Copy,
OVector<T::RealField, D>: Copy,
Owned<T, D, D>: Copy,
Owned<T::RealField, D>: Copy,
{
}
impl<T: ComplexField, D: Dim> Clone for SymmetricEigen<T, D>
where
DefaultAllocator: Allocator<T, D, D> + Allocator<T::RealField, D>,
Owned<T, D, D>: Clone,
Owned<T::RealField, D>: Clone,
{
fn clone(&self) -> Self {
Self {
eigenvectors: self.eigenvectors.clone(),
eigenvalues: self.eigenvalues.clone(),
}
}
}
impl<T: ComplexField, D: Dim> fmt::Debug for SymmetricEigen<T, D>
where
DefaultAllocator: Allocator<T, D, D> + Allocator<T::RealField, D>,
Owned<T, D, D>: fmt::Debug,
Owned<T::RealField, D>: fmt::Debug,
{
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
f.debug_struct("SymmetricEigen")
.field("eigenvectors", &self.eigenvectors)
.field("eigenvalues", &self.eigenvalues)
.finish()
}
}
impl<T: ComplexField, D: Dim> SymmetricEigen<T, D>
where
DefaultAllocator: Allocator<T, D, D> + Allocator<T::RealField, D>,
@ -270,7 +299,10 @@ where
///
/// This is useful if some of the eigenvalues have been manually modified.
#[must_use]
pub fn recompose(&self) -> OMatrix<T, D, D> {
pub fn recompose(&self) -> OMatrix<T, D, D>
where
Owned<T, D, D>: Clone,
{
let mut u_t = self.eigenvectors.clone();
for i in 0..self.eigenvalues.len() {
let val = self.eigenvalues[i];

57
src/linalg/symmetric_tridiagonal.rs
View File

@ -1,10 +1,13 @@
use std::fmt;
use std::mem::MaybeUninit;
#[cfg(feature = "serde-serialize-no-std")]
use serde::{Deserialize, Serialize};
use crate::allocator::Allocator;
use crate::base::{DefaultAllocator, OMatrix, OVector};
use crate::dimension::{Const, DimDiff, DimSub, U1};
use crate::storage::Storage;
use crate::storage::{Owned, Storage};
use simba::scalar::ComplexField;
use crate::linalg::householder;
@ -25,8 +28,7 @@ use crate::linalg::householder;
OMatrix<T, D, D>: Deserialize<'de>,
OVector<T, DimDiff<D, U1>>: Deserialize<'de>"))
)]
#[derive(Clone, Debug)]
pub struct SymmetricTridiagonal<T: ComplexField, D: DimSub<U1>>
pub struct SymmetricTridiagonal<T, D: DimSub<U1>>
where
DefaultAllocator: Allocator<T, D, D> + Allocator<T, DimDiff<D, U1>>,
{
@ -34,14 +36,42 @@ where
off_diagonal: OVector<T, DimDiff<D, U1>>,
}
impl<T: ComplexField, D: DimSub<U1>> Copy for SymmetricTridiagonal<T, D>
impl<T: Copy, D: DimSub<U1>> Copy for SymmetricTridiagonal<T, D>
where
DefaultAllocator: Allocator<T, D, D> + Allocator<T, DimDiff<D, U1>>,
OMatrix<T, D, D>: Copy,
OVector<T, DimDiff<D, U1>>: Copy,
Owned<T, D, D>: Copy,
Owned<T, DimDiff<D, U1>>: Copy,
{
}
impl<T: Clone, D: DimSub<U1>> Clone for SymmetricTridiagonal<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 {
tri: self.tri.clone(),
off_diagonal: self.off_diagonal.clone(),
}
}
}
impl<T: fmt::Debug, D: DimSub<U1>> fmt::Debug for SymmetricTridiagonal<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("SymmetricTridiagonal")
.field("tri", &self.tri)
.field("off_diagonal", &self.off_diagonal)
.finish()
}
}
impl<T: ComplexField, D: DimSub<U1>> SymmetricTridiagonal<T, D>
where
DefaultAllocator: Allocator<T, D, D> + Allocator<T, DimDiff<D, U1>>,
@ -61,24 +91,21 @@ where
"Unable to compute the symmetric tridiagonal decomposition of an empty matrix."
);
let mut off_diagonal = unsafe {
crate::unimplemented_or_uninitialized_generic!(dim.sub(Const::<1>), Const::<1>)
};
let mut p = unsafe {
crate::unimplemented_or_uninitialized_generic!(dim.sub(Const::<1>), Const::<1>)
};
let mut off_diagonal = OVector::new_uninitialized_generic(dim.sub(Const::<1>), Const::<1>);
let mut p = OVector::new_uninitialized_generic(dim.sub(Const::<1>), Const::<1>);
for i in 0..dim.value() - 1 {
let mut m = m.rows_range_mut(i + 1..);
let (mut axis, mut m) = m.columns_range_pair_mut(i, i + 1..);
let (norm, not_zero) = householder::reflection_axis_mut(&mut axis);
off_diagonal[i] = norm;
off_diagonal[i] = MaybeUninit::new(norm);
if not_zero {
let mut p = p.rows_range_mut(i..);
p.hegemv(crate::convert(2.0), &m, &axis, T::zero());
p.hegemv_z(crate::convert(2.0), &m, &axis);
let p = unsafe { p.slice_assume_init() };
let dot = axis.dotc(&p);
m.hegerc(-T::one(), &p, &axis, T::one());
@ -89,7 +116,7 @@ where
Self {
tri: m,
off_diagonal,
off_diagonal: unsafe { off_diagonal.assume_init() },
}
}

41
src/linalg/udu.rs
View File

@ -1,10 +1,12 @@
use std::fmt;
#[cfg(feature = "serde-serialize-no-std")]
use serde::{Deserialize, Serialize};
use crate::allocator::Allocator;
use crate::base::{Const, DefaultAllocator, OMatrix, OVector};
use crate::dimension::Dim;
use crate::storage::Storage;
use crate::storage::{Owned, Storage};
use simba::scalar::RealField;
/// UDU factorization.
@ -19,8 +21,7 @@ use simba::scalar::RealField;
deserialize = "OVector<T, D>: Deserialize<'de>, OMatrix<T, D, D>: Deserialize<'de>"
))
)]
#[derive(Clone, Debug)]
pub struct UDU<T: RealField, D: Dim>
pub struct UDU<T, D: Dim>
where
DefaultAllocator: Allocator<T, D> + Allocator<T, D, D>,
{
@ -30,14 +31,42 @@ where
pub d: OVector<T, D>,
}
impl<T: RealField, D: Dim> Copy for UDU<T, D>
impl<T: Copy, D: Dim> Copy for UDU<T, D>
where
DefaultAllocator: Allocator<T, D> + Allocator<T, D, D>,
OVector<T, D>: Copy,
OMatrix<T, D, D>: Copy,
Owned<T, D>: Copy,
Owned<T, D, D>: Copy,
{
}
impl<T: Clone, D: Dim> Clone for UDU<T, D>
where
DefaultAllocator: Allocator<T, D> + Allocator<T, D, D>,
Owned<T, D>: Clone,
Owned<T, D, D>: Clone,
{
fn clone(&self) -> Self {
Self {
u: self.u.clone(),
d: self.d.clone(),
}
}
}
impl<T: fmt::Debug, D: Dim> fmt::Debug for UDU<T, D>
where
DefaultAllocator: Allocator<T, D> + Allocator<T, D, D>,
Owned<T, D>: fmt::Debug,
Owned<T, D, D>: fmt::Debug,
{
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
f.debug_struct("UDU")
.field("u", &self.u)
.field("d", &self.d)
.finish()
}
}
impl<T: RealField, D: Dim> UDU<T, D>
where
DefaultAllocator: Allocator<T, D> + Allocator<T, D, D>,