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 f6dbec07d316b7209e5e24450647fab27719358f..4c2c7f3e65c7addad3bf7a89f2cae6299b5c9cbf 100644 (file)
--- a/src/Linear/Matrix.hs
+++ b/src/Linear/Matrix.hs
@@ -40,7 +40,7 @@ import qualified Data.Vector.Fixed as V (
   zipWith )
 import Data.Vector.Fixed.Cont ( Arity, arity )
 import Linear.Vector ( Vec, delete )
-import Naturals ( N1, N2, N3, N4, N5, N6, N7, N8, N9, N10, S, Z )
+import Naturals
 import Normed ( Normed(..) )
 
 import NumericPrelude hiding ( (*), abs )
@@ -67,6 +67,8 @@ type Mat2 a = Mat N2 N2 a
 type Mat3 a = Mat N3 N3 a
 type Mat4 a = Mat N4 N4 a
 type Mat5 a = Mat N5 N5 a
+type Mat6 a = Mat N6 N6 a
+type Mat7 a = Mat N7 N7 a
 
 -- * Type synonyms for 1-by-n row "vectors".
 
@@ -93,7 +95,29 @@ type Col6 a = Col N6 a
 type Col7 a = Col N7 a
 type Col8 a = Col N8 a
 type Col9 a = Col N9 a
-type Col10 a = Col N10 a -- We need a big column for Gaussian quadrature.
+type Col10 a = Col N10 a
+type Col11 a = Col N11 a
+type Col12 a = Col N12 a
+type Col13 a = Col N13 a
+type Col14 a = Col N14 a
+type Col15 a = Col N15 a
+type Col16 a = Col N16 a
+type Col17 a = Col N17 a
+type Col18 a = Col N18 a
+type Col19 a = Col N19 a
+type Col20 a = Col N20 a
+type Col21 a = Col N21 a
+type Col22 a = Col N22 a
+type Col23 a = Col N23 a
+type Col24 a = Col N24 a
+type Col25 a = Col N25 a
+type Col26 a = Col N26 a
+type Col27 a = Col N27 a
+type Col28 a = Col N28 a
+type Col29 a = Col N29 a
+type Col30 a = Col N30 a
+type Col31 a = Col N31 a
+type Col32 a = Col N32 a
 
 
 instance (Eq a) => Eq (Mat m n a) where
@@ -324,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