Speed up the Cholesky factorization with a fold.

[numerical-analysis.git] / src / Linear / Matrix.hs

diff --git a/src/Linear/Matrix.hs b/src/Linear/Matrix.hs
index a06151db14262838394f473f83841f4766c62266..4c2c7f3e65c7addad3bf7a89f2cae6299b5c9cbf 100644 (file)
--- a/src/Linear/Matrix.hs
+++ b/src/Linear/Matrix.hs
@@ -348,15 +348,21 @@ identity_matrix =
 --   >>> frobenius_norm (r - (transpose expected)) < 1e-12
 --   True
 --
-cholesky :: forall m n a. (Algebraic.C a, Arity m, Arity n)
-         => (Mat m n a) -> (Mat m n a)
-cholesky m = construct r
+cholesky :: forall m a. (Algebraic.C a, Arity m)
+         => (Mat m m a) -> (Mat m m a)
+cholesky m = ifoldl2 f zero m
   where
-    r :: Int -> Int -> a
-    r i j | i == j = sqrt(m !!! (i,j) - sum [(r k i)^2 | k <- [0..i-1]])
-          | i < j =
-              (((m !!! (i,j)) - sum [(r k i) NP.* (r k j) | k <- [0..i-1]]))/(r i i)
-          | otherwise = 0
+    f :: Int -> Int -> (Mat m m a) -> a -> (Mat m m a)
+    f i j cur_R _ = set_idx cur_R (i,j) (r cur_R i j)
+
+    r :: (Mat m m a) -> Int -> Int -> a
+    r cur_R i j
+      | i == j = sqrt(m !!! (i,j) - sum [(cur_R !!! (k,i))^2 | k <- [0..i-1]])
+      | i < j = (((m !!! (i,j))
+                  - sum [(cur_R !!! (k,i)) NP.* (cur_R !!! (k,j))
+                             | k <- [0..i-1]]))/(cur_R !!! (i,i))
+      | otherwise = 0
+
 
 
 -- | Returns True if the given matrix is upper-triangular, and False