[Haskell-cafe] Gauss Elimination -> More Clean2Haskell

[Philippos Apolinarius]

I am keeping with my project of translating programs from Clean to Haskell. As far as arrays go, I don't understand well how to use them in Haskell. Therefore, I will appreciate if somebody can find time to check the program below, and make suggestions to improve it. My Haskell program is about 4 times slower than the Clean version. It would be great if one could reduce the execution time by half, approaching to the speed of Clean and Scheme. Here are the constraints:

1 --- The program must be implemented using arrays.  Update must be done in place, in order to minimize use of the garbage collector. I have used Data.Array.IO, but I guess that  Data.Array.ST is better. Is it easy to rewrite the program in order to use DataArray.ST?

2 --  I liked very much the forM_ monad. I wonder if there is an accumulating monad that play the role of a program like leftSide.

3 -- Clean program almost double its speed, if one uses strictness annotations. Is it possible to use similar anotations for Haskell?

Here is how I compiled the program:

ghc -O2 gel.hs --make

In order to run the program with 2000 equations,  type

gel.exe 2000 +RTS -sstderr

The program will create a linear system with 2000 equations so that all elements of the solution is equal to 1. It prints 20 elements of the solution.

Here is the program:

{- File: gel.hs
Compilation: ghc -O2 gel.hs --make
Run: time gel.exe 2000
-}
import Control.Monad
import Data.Array.IO
import System.IO
import System.Random
import System (getArgs)

prtSol i n1 arr | i < 1= return ()
prtSol i n1 arr= do
    b <- readArray arr (i, n1)
    putStr ((show b)++" ")
    prtSol (i-1) n1 arr

fillArray xs s (i, n) (j, m) arr |  i > n= return ()
fillArray xs s (i,n) (j, m) arr | i==n && j>m= do
  writeArray arr (i, j) s
  return ()
fillArray xs s (i, n) (j, m) arr | j > m  = do
   writeArray arr (i, j) s
   fillArray xs 0.0 (i+1, n) (1, m) arr
fillArray (val:xs) s (i, n) (j, m) arr= do
   writeArray arr (i, j) val
   fillArray xs (s+val) (i, n) (j+1, m) arr

sLU arr n= sIJ 2 1 2 n arr

sIJ i j k n arr | i > n = return ()
sIJ i j k n arr | k > n = sIJ (i+1) i (i+1) n arr
sIJ i j k n arr = do
  im <- pmax (j+1) j
  swap j im 1
  a <- readArray arr (k, j)
  forM_ [j..n+1] $  \l -> do
      ajj <- readArray arr (j, j)
      ajl <- readArray arr (j, l)
      akl <- readArray arr (k, l) 
      writeArray arr (k, l) (akl-a*(ajl/ajj))
  sIJ i j (k+1) n arr where
     pmax line imax | line > n = return imax
     pmax line imax = do
       alj <- readArray arr (line, j)
       aij <- readArray arr (imax, j)
       if (abs alj)> (abs aij) 
          then pmax (line+1) line
          else pmax (line+1) imax
     swap r s q | q>n+1 = return ()
     swap r s q | r==s = return ()
     swap r s q = do
        arq <- readArray arr (r,q)
        asq <- readArray arr (s,q)
        writeArray arr (s,q) arq
        writeArray arr (r,q) asq
        swap r s (q+1)
     
  
leftSide acc i j n arr | j>n= return acc
leftSide acc i j n arr = do
   v <- readArray arr (j, n+1)
   a <- readArray arr (i, j)
   leftSide (acc-v*a) i (j+1) n arr

solvit i n arr | i<1= return ()
solvit i n arr= do
   a <- readArray arr (i, i)
   acc <- readArray arr (i, n+1)
   v <- leftSide acc i (i+1) n arr
   writeArray arr (i, n+1) (v/a)
   solvit (i-1) n arr

rnList :: (Double, Double) -> IO [Double]
rnList r=getStdGen>>=(\x->return(randomRs r x))

dims [input] = (read input, read input)
dims _ = (1000, 1000)

main = do
     xs <- rnList (1.0,1000.0)
     args <- getArgs
     let (n, m)= dims args
     arr <- newArray_ ((1,1),(n,m+1)) :: 
                IO (IOUArray (Int, Int) Double)
     fillArray xs 0.0 (1,n) (1,m) arr
     sLU arr n
     solvit n n arr
     prtSol (min 20 n) (n+1) arr
     print "Done"




      __________________________________________________________________
Be smarter than spam. See how smart SpamGuard is at giving junk email the boot with the All-new Yahoo! Mail.  Click on Options in Mail and switch to New Mail today or register for free at http://mail.yahoo.ca
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://www.haskell.org/pipermail/haskell-cafe/attachments/20091103/0494694a/attachment.html