Merge pull request #164 from dshizzle/master
Fix QR algorithm for diagonal matrices
This commit is contained in:
commit
aa1a4ceb58
|
@ -145,9 +145,23 @@ pub fn eigen_qr<N, V, VS, M>(m: &M, eps: &N, niter: usize) -> (M, V)
|
|||
for k in 0 .. curdim {
|
||||
let x_i = unsafe { eigenvalues.unsafe_at((k, k)) };
|
||||
let x_j = unsafe { eigenvalues.unsafe_at((k + 1, k)) };
|
||||
let r = (x_i * x_i + x_j * x_j).sqrt();
|
||||
let ctmp = x_i / r;
|
||||
let stmp = -x_j / r;
|
||||
|
||||
let ctmp;
|
||||
let stmp;
|
||||
|
||||
if x_j.abs() < *eps {
|
||||
ctmp = N::one();
|
||||
stmp = N::zero();
|
||||
}
|
||||
else if x_i.abs() < *eps {
|
||||
ctmp = N::zero();
|
||||
stmp = -N::one();
|
||||
}
|
||||
else {
|
||||
let r = x_i.hypot(x_j);
|
||||
ctmp = x_i / r;
|
||||
stmp = -x_j / r;
|
||||
}
|
||||
|
||||
c[k] = ctmp;
|
||||
s[k] = stmp;
|
||||
|
|
15
tests/mat.rs
15
tests/mat.rs
|
@ -104,6 +104,21 @@ macro_rules! test_eigen_qr_impl(
|
|||
|
||||
assert!(na::approx_eq_eps(&randmat, &recomp, &1.0e-2));
|
||||
}
|
||||
|
||||
for _ in (0usize .. 10000) {
|
||||
let randmat : $t = random();
|
||||
// Take only diagonal part
|
||||
let randmat: $t = Diag::from_diag(&randmat.diag());
|
||||
let (eigenvectors, eigenvalues) = na::eigen_qr(&randmat, &1e-13, 100);
|
||||
|
||||
let diag: $t = Diag::from_diag(&eigenvalues);
|
||||
let recomp = eigenvectors * diag * na::transpose(&eigenvectors);
|
||||
println!("eigenvalues: {:?}", eigenvalues);
|
||||
println!(" mat: {:?}", randmat);
|
||||
println!("recomp: {:?}", recomp);
|
||||
|
||||
assert!(na::approx_eq_eps(&randmat, &recomp, &1.0e-2));
|
||||
}
|
||||
}
|
||||
);
|
||||
|
||||
|
|
Loading…
Reference in New Issue