<br> Hello folks<br><br><br> 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...<br><br> 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. <br>
<br> I am releasing this code fragment as LGPL. <br><br><span style="font-family: courier new,monospace;">{-</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> This program is free software: you can redistribute it and/or modify</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> it under the terms of the GNU Lesser General Public License as published by</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> the Free Software Foundation, either version 3 of the License, or</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> (at your option) any later version.</span><br style="font-family: courier new,monospace;"><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> This program is distributed in the hope that it will be useful,</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> but WITHOUT ANY WARRANTY; without even the implied warranty of</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> GNU General Public License for more details.</span><br style="font-family: courier new,monospace;"><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> You should have received a copy of the GNU General Public License</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> along with this program. If not, see <<a href="http://www.gnu.org/licenses/">http://www.gnu.org/licenses/</a>>.</span><br style="font-family: courier new,monospace;">
<br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">-}</span><br style="font-family: courier new,monospace;"><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">import Data.Array.IArray</span><br style="font-family: courier new,monospace;">
<br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">type Dim=(Int,Int)</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"></span><span style="font-family: courier new,monospace;"> </span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">lu::Array Dim Double -> (Array Dim Double,Array Dim Double)</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">lu a = (aa l,aa u) </span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> where</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> (l,u)=lu' a [] [] </span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> aa = accumArray (+) 0 (bounds a)</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> lu'::(Floating e) => Array Dim e -> [(Dim,e)] -> [(Dim,e)] -> ([(Dim,e)],[(Dim,e)])</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> 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))</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> where</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> k=li</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> ((li,lj),(ui,uj))=bounds a</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> lik i=(a!(i,k)/a!(k,k))</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> un=[((k,j),a!(k,j)) | j<-[lj..uj]]</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> ln=((lj,lj),1.0):[((i,k),lik i) | i<-[li+1..ui]]</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> an=array ((li+1,lj+1),(ui,uj)) [((i,j),e_an i j) | i<-[li+1..ui],j<-[lj+1..uj]]</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> e_an i j=a!(i,j)-(lik i)*a!(k,j)</span><br style="font-family: courier new,monospace;"><br style="font-family: courier new,monospace;"><br><br> Still needed:<br>
<br>1) Partial pivoting<br>2) Profiling... Lots!<br>3) Parallelize<br><br clear="all"><br>-- <br>Rafael Gustavo da Cunha Pereira Pinto<br>Electronic Engineer, MSc.<br>