module Cholesky (Arr, cholesky) where

import Data.Array.IArray
import Data.Array.MArray
import Data.Array.Unboxed
import Data.Array.ST

type Idx = (Int,Int)
type Arr = UArray Idx Double



get :: Arr -> Arr -> Idx -> Double
get a l (i,j) | i == j = sqrt $ a!(j,j) - dot
              | i  > j = (a!(i,j) - dot) / l!(j,j)
              | otherwise = 0
  where dot = sum [l!(i,k) * l!(j,k) | k <- [0..j-1]]




cholesky :: Arr -> Arr
cholesky a = let n = maxBnd a
             in runSTUArray $ do
               l <- thaw a
               mapM_ (update a l) [(i,j) | i <- [0..n], j <- [0..n]]
               return l
  where maxBnd = fst . snd . bounds
        update a l i = unsafeFreeze l >>= \l' -> writeArray l i (get a l' i)