[Haskell-cafe] Purely funcional LU decomposition

[Rafael Gustavo da Cunha Pereira Pinto]

          Hello folks


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

         At first I thought it would be necessary to use a mutable or
monadic version of Array, but then I figured out it is a purely interactive
process.

         I am releasing this code fragment as LGPL.

{-
    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published
by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.

-}

import Data.Array.IArray

type Dim=(Int,Int)

lu::Array Dim Double -> (Array Dim Double,Array Dim Double)
lu a = (aa l,aa u)
      where
          (l,u)=lu' a [] []
          aa = accumArray (+) 0 (bounds a)
          lu'::(Floating e) => Array Dim e -> [(Dim,e)] -> [(Dim,e)] ->
([(Dim,e)],[(Dim,e)])
          lu' a l u=if (ui==li) then (
((ui,uj),1.0):l,((ui,uj),a!(ui,uj)):u) else (lu' an (l++ln) (u++un))
            where
              k=li
              ((li,lj),(ui,uj))=bounds a
              lik i=(a!(i,k)/a!(k,k))
              un=[((k,j),a!(k,j)) | j<-[lj..uj]]
              ln=((lj,lj),1.0):[((i,k),lik i) | i<-[li+1..ui]]
              an=array ((li+1,lj+1),(ui,uj)) [((i,j),e_an i j) |
i<-[li+1..ui],j<-[lj+1..uj]]
              e_an i j=a!(i,j)-(lik i)*a!(k,j)



          Still needed:

1) Partial pivoting
2) Profiling... Lots!
3) Parallelize


-- 
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/20090203/54e6c578/attachment.htm