On Wed, Feb 4, 2009 at 17:09, Dan Piponi
<dpiponi@gmail.com> wrote:
2009/2/3 Rafael Gustavo da Cunha Pereira Pinto <RafaelGCPP.Linux@gmail.com>:
> After a discussion on whether is possible to compile hmatrix in
> Windows, I decided to go crazy and do a LU decomposition entirely in
> Haskell...
> import Data.Array.IArray
...
> e_an i j=a!(i,j)-(lik i)*a!(k,j)
There are three different representations of arrays in this code:
arrays, lists and a functional one where f represents an array with
elements (f i j). Looks like a candidate for a stream fusion type
thing so the user only writes code using one representation and some
RULES switch back and forth between representations behind the scenes.
--
Dan
Those different representations are derived from my (very) low Haskell handicap! :-D
1) I used lists for the intermediate L and U matrices so I wouldn't have the overhead of calling accumArray and assocs everywhere
2) I used the function that returns an array element just to shorten the line length, and to make sure I was typing it right!
When I started this, I thought I would end with a couple of folds, but then I turned to the recursive version, so I could implement pivoting on matrix a' more easily.
List comprehension also derived from this. I started thinking of a fold for L matrix and then it hit me that it was a list comprehension...
Matt's DSP library has one of the most elegant implementations I saw, but it is harder to implement pivoting, which is really important for my application.
lu :: Array (Int,Int) Double -- ^ A
-> Array (Int,Int) Double -- ^ LU(A)
lu a = a'
where a' = array bnds [ ((i,j), luij i j) | (i,j) <- range bnds ]
luij i j | i>j = (a!(i,j) - sum [ a'!(i,k) * a'!(k,j) | k <- [1 ..(j-1)] ]) / a'!(j,j)
| i<=j = a!(i,j) - sum [ a'!(i,k) * a'!(k,j) | k <- [1 ..(i-1)] ]
bnds = bounds a
--
Rafael Gustavo da Cunha Pereira Pinto
Electronic Engineer, MSc.