From 179a6560ceb6fa6c74afa05036d1b3619effe4a1 Mon Sep 17 00:00:00 2001
From: Daniel D <sersmicro@gmail.com>
Date: Sat, 14 Nov 2015 15:40:35 +0100
Subject: [PATCH] Fix eigenvalue calculation for diagonal matrices

---
 src/linalg/decompositions.rs | 20 +++++++++++++++++---
 tests/mat.rs                 | 15 +++++++++++++++
 2 files changed, 32 insertions(+), 3 deletions(-)

diff --git a/src/linalg/decompositions.rs b/src/linalg/decompositions.rs
index df8a1c0f..d6e73322 100644
--- a/src/linalg/decompositions.rs
+++ b/src/linalg/decompositions.rs
@@ -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;
diff --git a/tests/mat.rs b/tests/mat.rs
index 72c14877..afae18d8 100644
--- a/tests/mat.rs
+++ b/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));
+        }
     }
 );