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     = 
| otherwise = (/u!(j,j)) $ a!(i,j) - sum [l!(i,k) * u!(k,j) | k <- [1..j-1]]
g (i,j) | i > j     = 
| otherwise = a!(i,j) - sum [l!(i,k) * u!(k,j) | k <- [1..i-1]]