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.