+
+-- | Given a positive-definite matrix @m@, computes the
+-- upper-triangular matrix @r@ with (transpose r)*r == m and all
+-- values on the diagonal of @r@ positive.
+--
+-- Examples:
+--
+-- >>> let m1 = fromList [[20,-1], [-1,20]] :: Mat2 Double
+-- >>> cholesky m1
+-- ((4.47213595499958,-0.22360679774997896),(0.0,4.466542286825459))
+-- >>> (transpose (cholesky m1)) `mult` (cholesky m1)
+-- ((20.000000000000004,-1.0),(-1.0,20.0))
+--
+cholesky :: forall a v w.
+ (RealFloat a,
+ Vector v (Vn w a),
+ Vector w a)
+ => (Mat v w a)
+ -> (Mat v w a)
+cholesky m = construct r
+ 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)*(r k j) | k <- [0..i-1]]))/(r i i)
+ | otherwise = 0
+
+-- | Matrix multiplication. Our 'Num' instance doesn't define one, and
+-- we need additional restrictions on the result type anyway.
+--
+-- Examples:
+--
+-- >>> let m1 = fromList [[1,2,3], [4,5,6]] :: Mat Vec2D Vec3D Int
+-- >>> let m2 = fromList [[1,2],[3,4],[5,6]] :: Mat Vec3D Vec2D Int
+-- >>> m1 `mult` m2
+-- ((22,28),(49,64))
+--
+mult :: (Num a,
+ Vector v (Vn w a),
+ Vector w a,
+ Vector w (Vn z a),
+ Vector z a,
+ Vector v (Vn z a))
+ => Mat v w a
+ -> Mat w z a
+ -> Mat v z a
+mult m1 m2 = construct lambda
+ where
+ lambda i j =
+ sum [(m1 !!! (i,k)) * (m2 !!! (k,j)) | k <- [0..(ncols m1)-1] ]