[Haskell-cafe] Purely funcional LU decomposition
Rafael Gustavo da Cunha Pereira Pinto
RafaelGCPP.Linux at gmail.com
Wed Feb 4 15:57:35 EST 2009
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
More information about the Haskell-Cafe
mailing list