Make matrix power work for non-owned matrices.

This commit is contained in:
Crozet Sébastien 2021-04-11 14:07:06 +02:00
parent 24d546d3b6
commit cc4427e52b
3 changed files with 29 additions and 16 deletions

View File

@ -140,9 +140,9 @@ where
} }
/// Computes the determinant of the decomposed matrix. /// Computes the determinant of the decomposed matrix.
pub fn determinant(&self) -> N::SimdRealField { pub fn determinant(&self) -> T::SimdRealField {
let dim = self.chol.nrows(); let dim = self.chol.nrows();
let mut prod_diag = N::one(); let mut prod_diag = T::one();
for i in 0..dim { for i in 0..dim {
prod_diag *= unsafe { *self.chol.get_unchecked((i, i)) }; prod_diag *= unsafe { *self.chol.get_unchecked((i, i)) };
} }

View File

@ -2,23 +2,27 @@
use std::ops::DivAssign; use std::ops::DivAssign;
use crate::{allocator::Allocator, DefaultAllocator, DimMin, MatrixN}; use crate::{
allocator::Allocator,
storage::{Storage, StorageMut},
DefaultAllocator, DimMin, Matrix, OMatrix,
};
use num::PrimInt; use num::PrimInt;
use simba::scalar::ComplexField; use simba::scalar::ComplexField;
impl<N: ComplexField, D> MatrixN<N, D> impl<T: ComplexField, D, S> Matrix<T, D, D, S>
where where
D: DimMin<D, Output = D>, D: DimMin<D, Output = D>,
DefaultAllocator: Allocator<N, D, D>, S: StorageMut<T, D, D>,
DefaultAllocator: Allocator<N, D>, DefaultAllocator: Allocator<T, D, D> + Allocator<T, D>,
{ {
/// Attempts to raise this matrix to an integral power `e` in-place. If this /// Attempts to raise this matrix to an integral power `e` in-place. If this
/// matrix is non-invertible and `e` is negative, it leaves this matrix /// matrix is non-invertible and `e` is negative, it leaves this matrix
/// untouched and returns `Err(())`. Otherwise, it returns `Ok(())` and /// untouched and returns `Err(())`. Otherwise, it returns `Ok(())` and
/// overwrites this matrix with the result. /// overwrites this matrix with the result.
#[must_use] #[must_use]
pub fn pow_mut<T: PrimInt + DivAssign>(&mut self, mut e: T) -> Result<(), ()> { pub fn pow_mut<I: PrimInt + DivAssign>(&mut self, mut e: I) -> Result<(), ()> {
let zero = T::zero(); let zero = I::zero();
// A matrix raised to the zeroth power is just the identity. // A matrix raised to the zeroth power is just the identity.
if e == zero { if e == zero {
@ -34,18 +38,19 @@ where
} }
} }
let one = T::one(); let one = I::one();
let two = T::from(2u8).unwrap(); let two = I::from(2u8).unwrap();
// We use the buffer to hold the result of multiplier ^ 2, thus avoiding // We use the buffer to hold the result of multiplier ^ 2, thus avoiding
// extra allocations. // extra allocations.
let mut multiplier = self.clone(); let mut multiplier = self.clone_owned();
let mut buf = self.clone(); let mut buf = self.clone_owned();
// Exponentiation by squares. // Exponentiation by squares.
loop { loop {
if e % two == one { if e % two == one {
*self *= &multiplier; self.mul_to(&multiplier, &mut buf);
self.copy_from(&buf);
} }
e /= two; e /= two;
@ -57,12 +62,20 @@ where
} }
} }
} }
}
impl<T: ComplexField, D, S: Storage<T, D, D>> Matrix<T, D, D, S>
where
D: DimMin<D, Output = D>,
S: StorageMut<T, D, D>,
DefaultAllocator: Allocator<T, D, D> + Allocator<T, D>,
{
/// Attempts to raise this matrix to an integral power `e`. If this matrix /// Attempts to raise this matrix to an integral power `e`. If this matrix
/// is non-invertible and `e` is negative, it returns `None`. Otherwise, it /// is non-invertible and `e` is negative, it returns `None`. Otherwise, it
/// returns the result as a new matrix. Uses exponentiation by squares. /// returns the result as a new matrix. Uses exponentiation by squares.
pub fn pow<T: PrimInt + DivAssign>(&self, e: T) -> Option<Self> { #[must_use]
let mut clone = self.clone(); pub fn pow<I: PrimInt + DivAssign>(&self, e: I) -> Option<OMatrix<T, D, D>> {
let mut clone = self.clone_owned();
match clone.pow_mut(e) { match clone.pow_mut(e) {
Ok(()) => Some(clone), Ok(()) => Some(clone),

View File

@ -92,7 +92,7 @@ macro_rules! gen_tests(
#[test] #[test]
fn cholesky_determinant_static(_n in PROPTEST_MATRIX_DIM) { fn cholesky_determinant_static(_n in PROPTEST_MATRIX_DIM) {
let m = RandomSDP::new(U4, || random::<$scalar>().0).unwrap(); let m = RandomSDP::new(Const::<4>, || random::<$scalar>().0).unwrap();
let lu_det = m.clone().lu().determinant(); let lu_det = m.clone().lu().determinant();
assert_relative_eq!(lu_det.imaginary(), 0., epsilon = 1.0e-7); assert_relative_eq!(lu_det.imaginary(), 0., epsilon = 1.0e-7);
let chol_det = m.cholesky().unwrap().determinant(); let chol_det = m.cholesky().unwrap().determinant();