Now that I've said my piece about syntax, we can get down to some
nitty gritty details of the language ...
The following are two different ways of factoring the L and U sections
needed in Crout reduction. The first version is a direct translation
of some FORTRANish pseudo-code and the second version is my attempt to
build the array monolithically.
lu_decomp
:: (Fractional b)
=> (Array (Int, Int) b) -> a -> a -> Array (Int, Int) b
lu_decomp_direct a m n = lu
where
ab = ((m,m), (n,n))
lu = lu1//[(n,n-1) := a!(n,n-1), (n,n) := a!(n,n) - a!(n,n-1)*lu!(n-1,n)]
lu1 = foldl f lu0 [m+1..n-1]
lu0 = array ab [(m,m) := a!(m,m), (m,m+1) := a!(m,m+1)/a!(m,m)]
f lu i = lu'
where
li1 = a!(i,i-1)
li2 = a!(i,i) - li1*lu!(i-1,i)
lu' = lu//[(i,i-1) := li1, (i,i) := li2, (i,i+1) := a!(i,i+1)/li2]
lu_decomp a m n = lu
where
ab = ((m,m), (n,n))
lu = array ab (lu_elems lun)
lu_elems = foldl f lu0 [m+1..n-1]
lu0 v = ((m,m) := a!(m,m)):(((m,m+1) := a!(m,m+1)/a!(m,m)):v)
lun = [(n,n-1) := a!(n,n-1), (n,n) := a!(n,n) - a!(n,n-1)*lu!(n-1,n)]
f le i = le'
where
li1 = a!(i,i-1)
li2 = a!(i,i) - li1*lu!(i-1,i)
le' v = le (((i,i-1) := li1):((i,i) := li2):((i,i+1) := a!(i,i+1)/li2):v)
My question is (and it is probably aimed at implementors): are either
likely to be optimised very well? e.g. update in place for the first
version and doing away with the thunks creating the list in the
second?
Alternately can anyone suggest any methods that are likely to be more
efficient? For example, I have a third version that uses a function
in a similar way to one of my Choleski solutions (see an Haskell
library), but there are so many guards I find it hard to believe it
will be particularly quick. Any timings I've done are distorted by
the fact that I only have one compiler (Hbc), exellent though it is.
BTW feel free to (publically) pick holes in the code if you think it
is "wrong" in any way e.g. should it be foldr instead foldl? (I think
not, but I'm never too sure :-)
bevan
