# [Haskell-cafe] Re: Purely funcional LU decomposition

Neal Alexander wqeqweuqy at hotmail.com
Wed Feb 4 02:15:49 EST 2009

```Array is no good man! Quad Tree matrices perform much nicer from what
I've seen.

I wrote some matrix stuff based on D. Stott Parker's "Randomized
Gaussian elimination" papers (http://www.cs.ucla.edu/~stott/ge/). He
presents some recursive block based methods of solving without pivoting.

I want to upload the full library to Hackage, but it'd be nice to have
some people look through it before - mainly because i never took linear
algebra heh. QuickCheck seems to be holding things together though.

Check these definitions (where the input matrix has been pre-transformed
with the butterfly if needed -- revertRBT . invert . applyRBT):

gauss (M.Scalar v)         = (M.Scalar 1, M.Scalar v)
gauss (M.Matrix (a,b,c,d)) = (M.Matrix (l, M.Scalar 0, x, l'),
M.Matrix (u, y, M.Scalar 0, u'))

where
(l,u)   = gauss a
(l',u') = gauss \$ d - (x * y)

x  = c * (invert u)
y  = (invert l) * b

invert (M.Scalar v)         = M.Scalar (1 / v)
invert (M.Matrix (a,b,c,d)) = M.Matrix (w,x,y,z)

where

w   = aI - (aIb * y)
x   = negate \$ aIb * z
y   = negate \$ z * caI
z   = invert \$ d - (caI * b)

-----------

aI  = invert a
caI = c * aI
aIb = aI * b

benchmarked against (in order of speed):

1. Mesa GL math in C called through FFI.
2. Finger tree based port.
3. Mutable array based port.

2 and 3 were scrapped right away.

Heres some of the benchmark vs the C version.

------------------------------------
Matrix Inversion - 999999999x
------------------------------------

ContextSwitches - 224017
First level fills = 0
Second level fills = 0
ETime(   0:00:41.328 ) UTime(   0:00:40.875 ) KTime(   0:00:00.000 )
ITime(   0:00:41.484 )

---------- C ---------------

ContextSwitches - 640016
First level fills = 0
Second level fills = 0
ETime(   0:01:56.109 ) UTime(   0:01:55.359 ) KTime(   0:00:00.000 )
ITime(   0:01:54.328 )

------------------------------------
Matrix Inversion - 9999999x
------------------------------------

ContextSwitches - 6204
First level fills = 0
Second level fills = 0
ETime(   0:00:00.687 ) UTime(   0:00:00.421 ) KTime(   0:00:00.000 )
ITime(   0:00:01.031 )

---------- C ---------------

ContextSwitches - 7396
First level fills = 0
Second level fills = 0
ETime(   0:00:01.171 ) UTime(   0:00:01.171 ) KTime(   0:00:00.000 )
ITime(   0:00:01.250 )

-------------------------------------
Matrix multiplication - 999999999x
-------------------------------------

ContextSwitches - 234644
First level fills = 0
Second level fills = 0
ETime(   0:00:42.093 ) UTime(   0:00:41.609 ) KTime(   0:00:00.000 )
ITime(   0:00:42.187 )

---------- C ---------------

ContextSwitches - 228596
First level fills = 0
Second level fills = 0
ETime(   0:00:41.609 ) UTime(   0:00:41.375 ) KTime(   0:00:00.000 )
ITime(   0:00:41.515 )

```