LU分解

LU分解@Haskell
遅延評価って便利．自分で評価の順番を明示的に与えなくても，
よきにはからってくれる．

import Data.Array

type Matrix = Array (Int, Int) Double

-- assumption: u!(i,i) /= 0
lu :: Matrix -> (Matrix, Matrix)
lu a = (l, u)
    where b = bounds a
          l = listArray b.map f.range $ b
          u = listArray b.map g.range $ b
          f (i,j) | i == j    = 1
                  | i < j     = 0
                  | otherwise = (/u!(j,j)) $ a!(i,j) - sum [l!(i,k) * u!(k,j) | k <- [1..j-1]]
          g (i,j) | i > j     = 0
                  | otherwise = a!(i,j) - sum [l!(i,k) * u!(k,j) | k <- [1..i-1]]