From 91b62c586d61adffd168124e38e1890cc87c2e5c Mon Sep 17 00:00:00 2001
From: Andreas Longva <andreas.b.longva@gmail.com>
Date: Mon, 13 Feb 2017 16:12:39 +0100
Subject: [PATCH] Cholesky: validate diagonal entries

---
 src/matrix/decomposition/cholesky.rs | 39 ++++++++++++++++++++++++++--
 1 file changed, 37 insertions(+), 2 deletions(-)

diff --git a/src/matrix/decomposition/cholesky.rs b/src/matrix/decomposition/cholesky.rs
index 60929d9..1566181 100644
--- a/src/matrix/decomposition/cholesky.rs
+++ b/src/matrix/decomposition/cholesky.rs
@@ -51,8 +51,16 @@ impl<T> Cholesky<T> where T: Float {
                 }
             }
 
-            // TODO: Check for zero divisor
-            let divisor = a[[j, j]].sqrt();
+            let diagonal = a[[j, j]];
+            if diagonal.abs() < T::epsilon() {
+                return Err(Error::new(ErrorKind::DecompFailure,
+                    "Matrix is singular to working precision."));
+            } else if diagonal < T::zero() {
+                return Err(Error::new(ErrorKind::DecompFailure,
+                    "Diagonal entries of matrix are not all positive."));
+            }
+
+            let divisor = diagonal.sqrt();
             for k in j .. n {
                 a[[k, j]] = a[[k, j]] / divisor;
             }
@@ -198,6 +206,33 @@ mod tests {
         }
     }
 
+    #[test]
+    fn cholesky_singular_fails() {
+        {
+            let x = matrix![0.0];
+            assert!(Cholesky::decompose(x).is_err());
+        }
+
+        {
+            let x = matrix![0.0, 0.0;
+                            0.0, 1.0];
+            assert!(Cholesky::decompose(x).is_err());
+        }
+
+        {
+            let x = matrix![1.0, 0.0;
+                            0.0, 0.0];
+            assert!(Cholesky::decompose(x).is_err());
+        }
+
+        {
+            let x = matrix![1.0,   3.0,   5.0;
+                            3.0,   9.0,  15.0;
+                            5.0,  15.0,  65.0];
+            assert!(Cholesky::decompose(x).is_err());
+        }
+    }
+
     quickcheck! {
         fn property_cholesky_of_identity_is_identity(n: usize) -> TestResult {
             if n > 30 {