[Haskell-cafe] Re: Purely funcional LU decomposition

[Neal Alexander]

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

------- haskell ------------

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

------- haskell ------------

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

------- haskell ------------

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 )