[Haskell-cafe] Purely funcional LU decomposition

[Rafael Gustavo da Cunha Pereira Pinto]

On Wed, Feb 4, 2009 at 17:09, Dan Piponi <dpiponi at gmail.com> wrote:

> 2009/2/3 Rafael Gustavo da Cunha Pereira Pinto <RafaelGCPP.Linux at 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.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://www.haskell.org/pipermail/haskell-cafe/attachments/20090204/5a34d28c/attachment.htm