-----------------------------------------------------------------------------
-- |
-- Module      :  Matrix.Cholesky
-- Copyright   :  (c) Matthew Donadio 2003
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- This module contains a routine that solves the system Ax=b, where A
-- is positive definite, using Cholesky decomposition.
--
-----------------------------------------------------------------------------


module Matrix.Cholesky (cholesky) where

import Data.Array
import Data.Complex

-- * Functions

-- Formulas 2.53--2.55 in Kay

cholesky :: (Ix a, Integral a, RealFloat b) => Array (a,a) (Complex b) -- ^ A
	                                    -> Array a (Complex b)     -- ^ b
					    -> Array a (Complex b)     -- ^ x
cholesky :: forall a b.
(Ix a, Integral a, RealFloat b) =>
Array (a, a) (Complex b)
-> Array a (Complex b) -> Array a (Complex b)
cholesky Array (a, a) (Complex b)
a Array a (Complex b)
b = Array a (Complex b)
x
    where y :: Array a (Complex b)
y = forall i e. Ix i => (i, i) -> [(i, e)] -> Array i e
array (a
1,a
n) ((a
1,Array a (Complex b)
bforall i e. Ix i => Array i e -> i -> e
!a
1) forall a. a -> [a] -> [a]
: [ (a
k, Array a (Complex b)
bforall i e. Ix i => Array i e -> i -> e
!a
k forall a. Num a => a -> a -> a
- forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ Array (a, a) (Complex b)
lforall i e. Ix i => Array i e -> i -> e
!(a
k,a
j) forall a. Num a => a -> a -> a
* Array a (Complex b)
yforall i e. Ix i => Array i e -> i -> e
!a
j | a
j <- [a
1..(a
kforall a. Num a => a -> a -> a
-a
1)] ] ) | a
k <- [a
2..a
n] ])
	  x :: Array a (Complex b)
x = forall i e. Ix i => (i, i) -> [(i, e)] -> Array i e
array (a
1,a
n) ((a
n, Array a (Complex b)
yforall i e. Ix i => Array i e -> i -> e
!a
n forall a. Fractional a => a -> a -> a
/ Array a (Complex b)
dforall i e. Ix i => Array i e -> i -> e
!a
n) forall a. a -> [a] -> [a]
: [ (a
k, Array a (Complex b)
yforall i e. Ix i => Array i e -> i -> e
!a
k forall a. Fractional a => a -> a -> a
/ Array a (Complex b)
dforall i e. Ix i => Array i e -> i -> e
!a
k forall a. Num a => a -> a -> a
- forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ (forall a. Num a => Complex a -> Complex a
conjugate (Array (a, a) (Complex b)
lforall i e. Ix i => Array i e -> i -> e
!(a
j,a
k))) forall a. Num a => a -> a -> a
* Array a (Complex b)
xforall i e. Ix i => Array i e -> i -> e
!a
j | a
j <- [(a
kforall a. Num a => a -> a -> a
+a
1)..a
n] ] ) | a
k <- (forall a. [a] -> [a]
reverse [a
1..(a
nforall a. Num a => a -> a -> a
-a
1)]) ])
	  l :: Array (a, a) (Complex b)
l = forall i e. Ix i => (i, i) -> [(i, e)] -> Array i e
array ((a
1,a
1),(a
n,a
n)) [ ((a
i,a
j), a -> a -> Complex b
lij a
i a
j) | a
i <- [a
2..a
n], a
j <- [a
1..(a
iforall a. Num a => a -> a -> a
-a
1)] ]
	  lij :: a -> a -> Complex b
lij a
i a
j | a
jforall a. Eq a => a -> a -> Bool
==a
1      = Array (a, a) (Complex b)
aforall i e. Ix i => Array i e -> i -> e
!(a
i,a
1) forall a. Fractional a => a -> a -> a
/ Array a (Complex b)
dforall i e. Ix i => Array i e -> i -> e
!a
1
		    | Bool
otherwise = Array (a, a) (Complex b)
aforall i e. Ix i => Array i e -> i -> e
!(a
i,a
j) forall a. Fractional a => a -> a -> a
/ Array a (Complex b)
dforall i e. Ix i => Array i e -> i -> e
!a
j forall a. Num a => a -> a -> a
- forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ Array (a, a) (Complex b)
lforall i e. Ix i => Array i e -> i -> e
!(a
i,a
k) forall a. Num a => a -> a -> a
* Array a (Complex b)
dforall i e. Ix i => Array i e -> i -> e
!a
k forall a. Num a => a -> a -> a
* (forall a. Num a => Complex a -> Complex a
conjugate (Array (a, a) (Complex b)
lforall i e. Ix i => Array i e -> i -> e
!(a
j,a
k))) forall a. Fractional a => a -> a -> a
/ Array a (Complex b)
dforall i e. Ix i => Array i e -> i -> e
!a
j | a
k <- [a
1..(a
jforall a. Num a => a -> a -> a
-a
1)] ]
	  d :: Array a (Complex b)
d = forall i e. Ix i => (i, i) -> [(i, e)] -> Array i e
array (a
1,a
n) ((a
1, Array (a, a) (Complex b)
aforall i e. Ix i => Array i e -> i -> e
!(a
1,a
1)) forall a. a -> [a] -> [a]
: [ (a
i, Array (a, a) (Complex b)
aforall i e. Ix i => Array i e -> i -> e
!(a
i,a
i) forall a. Num a => a -> a -> a
- forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ Array a (Complex b)
dforall i e. Ix i => Array i e -> i -> e
!a
k forall a. Num a => a -> a -> a
* ((forall a. Num a => a -> a
abs (Array (a, a) (Complex b)
lforall i e. Ix i => Array i e -> i -> e
!(a
i,a
k)))forall a b. (Num a, Integral b) => a -> b -> a
^(Int
2::Int)) | a
k <- [a
1..(a
iforall a. Num a => a -> a -> a
-a
1)] ] ) | a
i <- [a
2..a
n]])
	  ((a
_,a
_),(a
n,a
_)) = forall i e. Array i e -> (i, i)
bounds Array (a, a) (Complex b)
a