nalgebra/tests/linalg/cholesky.rs

127 lines
5.3 KiB
Rust
Raw Normal View History

#![cfg(all(feature = "arbitrary", feature = "debug"))]
macro_rules! gen_tests(
($module: ident, $scalar: ty) => {
mod $module {
use na::debug::RandomSDP;
use na::dimension::{U4, Dynamic};
use na::{DMatrix, DVector, Matrix4x3, Vector4};
use rand::random;
#[allow(unused_imports)]
2019-03-23 21:29:07 +08:00
use crate::core::helper::{RandScalar, RandComplex};
use std::cmp;
quickcheck! {
fn cholesky(n: usize) -> bool {
let m = RandomSDP::new(Dynamic::new(n.max(1).min(50)), || random::<$scalar>().0).unwrap();
let l = m.clone().cholesky().unwrap().unpack();
relative_eq!(m, &l * l.adjoint(), epsilon = 1.0e-7)
}
fn cholesky_static(_m: RandomSDP<f64, U4>) -> bool {
let m = RandomSDP::new(U4, || random::<$scalar>().0).unwrap();
let chol = m.cholesky().unwrap();
let l = chol.unpack();
if !relative_eq!(m, &l * l.adjoint(), epsilon = 1.0e-7) {
false
}
else {
true
}
}
fn cholesky_solve(n: usize, nb: usize) -> bool {
let n = n.max(1).min(50);
let m = RandomSDP::new(Dynamic::new(n), || random::<$scalar>().0).unwrap();
let nb = cmp::min(nb, 50); // To avoid slowing down the test too much.
let chol = m.clone().cholesky().unwrap();
let b1 = DVector::<$scalar>::new_random(n).map(|e| e.0);
let b2 = DMatrix::<$scalar>::new_random(n, nb).map(|e| e.0);
let sol1 = chol.solve(&b1);
let sol2 = chol.solve(&b2);
relative_eq!(&m * &sol1, b1, epsilon = 1.0e-7) &&
relative_eq!(&m * &sol2, b2, epsilon = 1.0e-7)
}
fn cholesky_solve_static(_n: usize) -> bool {
let m = RandomSDP::new(U4, || random::<$scalar>().0).unwrap();
let chol = m.clone().cholesky().unwrap();
let b1 = Vector4::<$scalar>::new_random().map(|e| e.0);
let b2 = Matrix4x3::<$scalar>::new_random().map(|e| e.0);
let sol1 = chol.solve(&b1);
let sol2 = chol.solve(&b2);
relative_eq!(m * sol1, b1, epsilon = 1.0e-7) &&
relative_eq!(m * sol2, b2, epsilon = 1.0e-7)
}
fn cholesky_inverse(n: usize) -> bool {
let m = RandomSDP::new(Dynamic::new(n.max(1).min(50)), || random::<$scalar>().0).unwrap();
let m1 = m.clone().cholesky().unwrap().inverse();
let id1 = &m * &m1;
let id2 = &m1 * &m;
id1.is_identity(1.0e-7) && id2.is_identity(1.0e-7)
}
fn cholesky_inverse_static(_n: usize) -> bool {
let m = RandomSDP::new(U4, || random::<$scalar>().0).unwrap();
let m1 = m.clone().cholesky().unwrap().inverse();
let id1 = &m * &m1;
let id2 = &m1 * &m;
id1.is_identity(1.0e-7) && id2.is_identity(1.0e-7)
}
2019-11-02 22:11:14 +08:00
fn cholesky_rank_one_update(_n: usize) -> bool {
2019-11-03 01:27:01 +08:00
let mut m = RandomSDP::new(U4, || random::<$scalar>().0).unwrap();
let x = Vector4::<$scalar>::new_random().map(|e| e.0);
2019-11-02 23:45:30 +08:00
2019-11-03 02:04:07 +08:00
// this is dirty but $scalar is not a scalar type (its a Rand) in this file
2019-11-02 23:45:30 +08:00
let zero = random::<$scalar>().0 * 0.;
let one = zero + 1.;
let sigma = random::<f64>(); // needs to be a real
let sigma_scalar = zero + sigma;
2019-11-02 22:56:59 +08:00
2019-11-03 02:04:07 +08:00
// updates cholesky decomposition and reconstructs m updated
2019-11-02 22:56:59 +08:00
let mut chol = m.clone().cholesky().unwrap();
chol.rank_one_update(&x, sigma);
let m_chol_updated = chol.l() * chol.l().adjoint();
2019-11-02 22:11:14 +08:00
// updates m manually
2019-11-03 01:27:01 +08:00
m.gerc(sigma_scalar, &x, &x, one); // m += sigma * x * x.adjoint()
2019-11-02 22:11:14 +08:00
2019-11-02 22:56:59 +08:00
relative_eq!(m, m_chol_updated, epsilon = 1.0e-7)
2019-11-02 22:11:14 +08:00
}
2019-11-03 21:33:35 +08:00
fn cholesky_remove_column(n: usize) -> bool {
let n = n.max(1).min(5);
let j = random::<usize>() % n;
let m = RandomSDP::new(Dynamic::new(n), || random::<$scalar>().0).unwrap();
// remove column from cholesky decomposition and rebuild m
let chol = m.clone().cholesky().unwrap().remove_column(j);
let m_chol_updated = chol.l() * chol.l().adjoint();
// remove column from m
let m_updated = m.remove_column(j).remove_row(j);
println!("n={} j={}", n, j);
println!("chol:{}", m_chol_updated);
println!("m up:{}", m_updated);
relative_eq!(m_updated, m_chol_updated, epsilon = 1.0e-7)
}
}
}
}
);
gen_tests!(complex, RandComplex<f64>);
gen_tests!(f64, RandScalar<f64>);