Re: [Haskell-cafe] Purely funcional LU decomposition

On Wed, Feb 4, 2009 at 17:09, Dan Piponi wrote:

2009/2/3 Rafael Gustavo da Cunha Pereira Pinto

...

: 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.