<br><br><div class="gmail_quote">On Wed, Feb 4, 2009 at 17:09, Dan Piponi <span dir="ltr">&lt;<a href="mailto:dpiponi@gmail.com" target="_blank">dpiponi@gmail.com</a>&gt;</span> wrote:<br><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">

2009/2/3 Rafael Gustavo da Cunha Pereira Pinto &lt;<a href="mailto:RafaelGCPP.Linux@gmail.com" target="_blank">RafaelGCPP.Linux@gmail.com</a>&gt;:<br>
<div>&gt; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;After a discussion on whether is possible to compile hmatrix in<br>
&gt; Windows, I decided to go crazy and do a LU decomposition entirely in<br>
&gt; Haskell...<br>
</div>&gt; import Data.Array.IArray<br>
...<br>
<div>&gt; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; e_an i j=a!(i,j)-(lik i)*a!(k,j)<br>
<br>
</div>There are three different representations of arrays in this code:<br>
arrays, lists and a functional one where f represents an array with<br>
elements (f i j). Looks like a candidate for a stream fusion type<br>
thing so the user only writes code using one representation and some<br>
RULES switch back and forth between representations behind the scenes.<br>
--<br>
<font color="#888888">Dan<br>
</font></blockquote></div><br>Those different representations are derived from my (very) low Haskell handicap! :-D<br><br>1) I used lists for the intermediate L and U matrices so I wouldn&#39;t have the overhead of calling accumArray and assocs everywhere<br>

2) I used the function that returns an array element just to shorten the line length, and to make sure I was typing it right!<br><br>When I started this, I thought I would end with a couple of folds, but then I turned to the recursive version, so I could implement pivoting on matrix a&#39; more easily.<br>

<br>List comprehension also derived from this. I started thinking of a fold for L matrix and then it hit me that it was a list comprehension...<br><br>Matt&#39;s DSP library has one of the most elegant implementations I saw, but it is harder to implement pivoting, which is really important for my application.<br>
<br><pre><font style="font-family: courier new,monospace;" size="2">lu :: Array (Int,Int) Double -- ^ A<br>   -&gt; Array (Int,Int) Double -- ^ LU(A)<br><br>lu a = a&#39;<br>    where a&#39; = array bnds [ ((i,j), luij i j) | (i,j) &lt;- range bnds ]<br>
          luij i j | i&gt;j  = (a!(i,j) - sum [ a&#39;!(i,k) * a&#39;!(k,j) | k &lt;- [1 ..(j-1)] ]) / a&#39;!(j,j)<br>                   | i&lt;=j =  a!(i,j) - sum [ a&#39;!(i,k) * a&#39;!(k,j) | k &lt;- [1 ..(i-1)] ]<br>          bnds = bounds a</font><br>
</pre><br clear="all"><br>-- <br>Rafael Gustavo da Cunha Pereira Pinto<br>Electronic Engineer, MSc.<br>