gistfile1.hs

module Main where

omitElementAt :: Int -> [a] -> [a]
omitElementAt i list = let (left, right) = splitAt i list in left ++ (tail right)

jacobiProductSum :: Int -> [Double] -> [Double] -> Double
jacobiProductSum i a x = sum [ aj * xj | (aj, xj) <- aFiltered `zip` xFiltered ]
    where
        aFiltered = omitElementAt i a
        xFiltered = omitElementAt i x 

jacobiXi :: Int -> [[Double]] -> [Double] -> [Double] -> Double
jacobiXi i a b x = (bi - sigma) / aii
    where
        bi    = b !! i
        ai    = a !! i
        aii   = ai !! i
        sigma = jacobiProductSum i ai x

jacobiIteration :: [[Double]] -> [Double] -> [Double] -> [Double]
jacobiIteration a b x = [ jacobiXi i a b x | i <- take (length x) [0..] ]

jacobiIterate :: Int -> [[Double]] -> [Double] -> [Double] -> [Double]
jacobiIterate 1 a b x = jacobiIteration a b x
jacobiIterate n a b x = jacobiIterate (n - 1) a b (jacobiIteration a b x)

-- Wikipedia example
-- answers:
-- x1 = [5, 1.143]
-- x2 = [4.929, -1.713]
-- x25 = [7.111, -3.222]
a  = [ [2, 1]
     , [5, 7]
     ]
b  = [11, 13]
x0 = [1, 1]

main = do
    putStrLn $ show $ jacobiIterate 1 a b x0
    putStrLn $ show $ jacobiIterate 2 a b x0
    putStrLn $ show $ jacobiIterate 25 a b x0