【Haskell】再帰でpivot選択付きLU分解 ついに完成!

import Data.Array

-- MatrixとVectorの定義
type Matrix a = Array (Int, Int) a
type Vector a = Array Int a

-- リストからMatrixをつくる
makeMatrix :: [[a]] -> Matrix a
makeMatrix xs = listArray ((1, 1), (m, n)) $ concat xs
    where m = length xs
          n = length $ head xs

-- リストからVectorをつくる
makeVector :: [a] -> Vector a
makeVector xs = listArray (1, n) xs
    where n = length xs

-- Matrixの行を入れ替える
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

-- Matrixのpivotが何行目にあるか
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]]

-- pivot選択
pivoting :: Ord a => Matrix a -> Matrix a
pivoting a = swapRows 1 (pivot a) a

-- 再帰的なLU分解(次のMatrix, L_Matrix, U_Matrix, pivotの行のリスト)
next :: (Matrix Double, Matrix Double, Matrix Double, [(Int, Int)]) -> (Matrix Double, Matrix Double, Matrix Double, [(Int, 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, k), 1.0) | k <- [1..max]] ++ [((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' ++ [(max, max)]
             |otherwise = p' ++ [(1 + (max-n), (pivot a) + (max-n))]

-- LU分解(L_Matrix, U_Matrix, pivotの行のリスト)
ludecomp :: Matrix Double -> (Matrix Double, Matrix Double, [(Int, Int)])
ludecomp a = (l, u, p)
    where ((1, 1), (m, n)) = bounds a
          x = listArray ((1, 1), (n, n)) $ repeat 0.0
          (a', l, u, p) = next (a, x, x, [])

-- 連立方程式(ax = b)をLU分解で解く
solve :: Matrix Double -> Vector Double -> Vector Double
solve a b = x
    where (l, u, p) = ludecomp a
          b' = swap p b
          y = lowerSolve l b'
          x = upperSolve u y

-- pivotの行のリストを用いてVectorの要素を入れ替える
swap :: [(Int, Int)] -> Vector a -> Vector a
swap [] v = v
swap (n:ns) v = swap ns (v // [(i, v!j), (j, v!i)])
    where (i, j) = n

-- 下三角行列lの連立方程式(ly = b)を解く
lowerSolve :: Matrix Double -> Vector Double -> Vector Double
lowerSolve l b = y
    where (1, n) = bounds b
          y = array (1, n) ([(1, b!1 / l!(1, 1))] ++ [(k, (b!k - sum [l!(k, i) * y!i | i <- [1..k-1]]) / l!(k, k)) | k <- [2..n]])

-- 上三角行列uの連立方程式(ux = y)を解く
upperSolve :: Matrix Double -> Vector Double -> Vector Double
upperSolve u y = x
    where (1, n) = bounds y
          x = array (1, n) ([(n, y!n / u!(n, n))] ++ [(k, (y!k - sum [u!(k, i) * x!i | i <- [k+1..n]]) / u!(k, k)) | k <- reverse [1..n-1]])