import Data.Array
type Matrix a = Array (Int, Int) a
makeMatrix :: [[a]] -> Matrix a
makeMatrix xs = listArray ((1, 1), (m, n)) $ concat xs
where m = length xs
n = length $ head xs
swapRows :: Int -> Int -> Matrix a -> Matrix a
swapRows i j a = a // ([((i, k), a!(j, k)) | k <- [1..n]] ++
[((j, k), a!(i, k)) | k <- [1..n]])
where ((1, 1), (m, n)) = bounds a
pivot :: Ord a => Matrix a -> Int
pivot a = head [i | i <- [1..m], a!(i, 1) == max]
where ((1, 1), (m, n)) = bounds a
max = maximum [a!(i, 1) | i <- [1..m]]
pivoting :: Ord a => Matrix a -> Matrix a
pivoting a = swapRows 1 (pivot a) a
next :: (Matrix Double, Matrix Double, Matrix Double, [Int]) -> (Matrix Double, Matrix Double, Matrix Double, [Int])
next (a, x', y', p') |length p' == max = (a, x', y', p')
|otherwise = next (b, x, y, p)
where ((1, 1), (m, n)) = bounds a
((1, 1), (maxm, max)) = bounds x'
a' = pivoting a
b = array ((1, 1), (n-1, n-1)) [((i, j), a'!(i+1, j+1) - (a'!(i+1, 1) / a'!(1, 1)) * a'!(1, j+1)) | i <- [1..n-1], j <- [1..n-1]]
x |p' == [] = x' // [((k, 1), a'!(k, 1) / a'!(1, 1)) | k <- [2..n]]
|n == 1 = x'
|otherwise = x' // [((k+(max-n), max-n+1), a'!(k, 1) / a'!(1, 1)) | k <- [2..n]]
y |p' == [] = y' // [((1, l), a'!(1, l)) | l <- [1..n]]
|n == 1 = y' // [((max, max), a!(1, 1))]
|otherwise = y' // [((max-n+1, q+(max-n)), a'!(1, q)) | q <- [1..n]]
p |n == 1 = p' ++ [pivot a]
|otherwise = p' ++ [(pivot a) + (max-n)]
ludecomp :: Matrix Double -> (Matrix Double, Matrix Double, [Int])
ludecomp a = (l, u, p)
where ((1, 1), (m, n)) = bounds a
x = listArray ((1, 1), (n, n)) $ repeat 0.0
y = x
(a', l, u, p) = next (a, x, y, [])
