Fix eigenvalue calculation for diagonal matrices

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 ctmp;
-                let stmp = -x_j / r;
+                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;

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));
+        }
     }
 );