Difference between revisions of "Shootout/Spectral"
< Shootout
Jump to navigation
Jump to search
DonStewart (talk | contribs) (moved) |
DonStewart (talk | contribs) |
||
Line 10: | Line 10: | ||
||old || 10.564|| |
||old || 10.564|| |
||
− | == |
+ | == Current entry == |
+ | |||
+ | Better gcc flags. Bang patterns. |
||
+ | [https://alioth.debian.org/tracker/index.php?func=detail&aid=304453&group_id=30402&atid=411646 Submitted] |
||
+ | |||
+ | <haskell> |
||
+ | {-# OPTIONS -O -fglasgow-exts -fbang-patterns -funbox-strict-fields -fexcess-precision -optc-O2 -optc-march=pentium4 -optc-mfpmath=sse -optc-msse2 #-} |
||
+ | -- |
||
+ | -- The Great Computer Language Shootout |
||
+ | -- http:--shootout.alioth.debian.org/ |
||
+ | -- |
||
+ | -- Contributed by Don Stewart |
||
+ | -- |
||
+ | |||
+ | import Monad |
||
+ | import System |
||
+ | import Text.Printf |
||
+ | import Data.Array.IO |
||
+ | import Data.Array.Base |
||
+ | |||
+ | main = getArgs >>= approximate . read . head >>= printf "%.9f\n" |
||
+ | |||
+ | approximate n = do |
||
+ | u <- newArray (0,n-1) 1 :: IO (IOUArray Int Double) |
||
+ | v <- newArray (0,n-1) 0 :: IO (IOUArray Int Double) |
||
+ | sequence_ $ replicate 10 $ multiplyAtAv n u v >> multiplyAtAv n v u |
||
+ | |||
+ | let loop !vbv !vv !i |
||
+ | | i >= n = return (vbv,vv) |
||
+ | | otherwise = do ui <- unsafeRead u i |
||
+ | vi <- unsafeRead v i |
||
+ | loop (vbv + ui * vi) (vv + vi * vi) (i+1) |
||
+ | (vbv,vv) <- loop 0 0 0 |
||
+ | return $! sqrt (vbv/vv) |
||
+ | |||
+ | -- return element i,j of infinite matrix A |
||
+ | a i j = 1 / fromIntegral (x*(x+1) `div` 2 + i + 1) where x = i+j |
||
+ | |||
+ | -- multiply vector v by matrix A */ |
||
+ | multiplyAv !n !v !av = loop 0 |
||
+ | where loop i = when (i < n) $ do avi <- loop' i 0 0 |
||
+ | unsafeWrite av i avi >> loop (i+1) |
||
+ | loop' !i !j !av | j >= n = return av |
||
+ | | otherwise = do vj <- unsafeRead v j |
||
+ | loop' i (j+1) (av + a i j * vj) |
||
+ | |||
+ | -- multiply vector v by matrix A transposed |
||
+ | multiplyAtv !n !v !atv = loop 0 |
||
+ | where loop i = when (i < n) $ do atvi <- loop' i 0 0 |
||
+ | unsafeWrite atv i atvi >> loop (i+1) |
||
+ | loop' !i !j !atvi |
||
+ | | j >= n = return atvi |
||
+ | | otherwise = do vj <- unsafeRead v j |
||
+ | loop' i (j+1) (atvi + a j i * vj) |
||
+ | |||
+ | -- multiply vector v by matrix A and then by matrix A transposed */ |
||
+ | multiplyAtAv !n !v !atav = do |
||
+ | u <- newArray (0,n-1) 0 :: IO (IOUArray Int Double) |
||
+ | multiplyAv n v u >> multiplyAtv n u atav |
||
+ | </haskell> |
||
+ | |||
+ | == Old entry == |
||
GHC unboxed this better. Careful attention payed to the unboxing. |
GHC unboxed this better. Careful attention payed to the unboxing. |
Revision as of 04:21, 8 February 2007
A Shootout Entry for the spectral norm test
Timings
Debian Linux x86, n=2,500
||Entry || Time|| ||dons || 9.407|| ||old || 10.564||
Current entry
Better gcc flags. Bang patterns. Submitted
{-# OPTIONS -O -fglasgow-exts -fbang-patterns -funbox-strict-fields -fexcess-precision -optc-O2 -optc-march=pentium4 -optc-mfpmath=sse -optc-msse2 #-}
--
-- The Great Computer Language Shootout
-- http:--shootout.alioth.debian.org/
--
-- Contributed by Don Stewart
--
import Monad
import System
import Text.Printf
import Data.Array.IO
import Data.Array.Base
main = getArgs >>= approximate . read . head >>= printf "%.9f\n"
approximate n = do
u <- newArray (0,n-1) 1 :: IO (IOUArray Int Double)
v <- newArray (0,n-1) 0 :: IO (IOUArray Int Double)
sequence_ $ replicate 10 $ multiplyAtAv n u v >> multiplyAtAv n v u
let loop !vbv !vv !i
| i >= n = return (vbv,vv)
| otherwise = do ui <- unsafeRead u i
vi <- unsafeRead v i
loop (vbv + ui * vi) (vv + vi * vi) (i+1)
(vbv,vv) <- loop 0 0 0
return $! sqrt (vbv/vv)
-- return element i,j of infinite matrix A
a i j = 1 / fromIntegral (x*(x+1) `div` 2 + i + 1) where x = i+j
-- multiply vector v by matrix A */
multiplyAv !n !v !av = loop 0
where loop i = when (i < n) $ do avi <- loop' i 0 0
unsafeWrite av i avi >> loop (i+1)
loop' !i !j !av | j >= n = return av
| otherwise = do vj <- unsafeRead v j
loop' i (j+1) (av + a i j * vj)
-- multiply vector v by matrix A transposed
multiplyAtv !n !v !atv = loop 0
where loop i = when (i < n) $ do atvi <- loop' i 0 0
unsafeWrite atv i atvi >> loop (i+1)
loop' !i !j !atvi
| j >= n = return atvi
| otherwise = do vj <- unsafeRead v j
loop' i (j+1) (atvi + a j i * vj)
-- multiply vector v by matrix A and then by matrix A transposed */
multiplyAtAv !n !v !atav = do
u <- newArray (0,n-1) 0 :: IO (IOUArray Int Double)
multiplyAv n v u >> multiplyAtv n u atav
Old entry
GHC unboxed this better. Careful attention payed to the unboxing. -O2 -optc-O -optc-ffast-math -fexcess-precision
{-# OPTIONS -optc-O #-}
--
-- The Great Computer Language Shootout
-- http:--shootout.alioth.debian.org/
--
-- Contributed by Don Stewart
--
-- gcc miscompiles this program with -O3
--
import Monad
import System
import Numeric
import Data.Array.IO
import Data.Array.Base
main = getArgs >>= approximate . read . head >>= putStrLn . \s -> showFFloat (Just 9) s []
approximate n = do
u <- newArray (0,n-1) 1 :: IO (IOUArray Int Double)
v <- newArray_ (0,n-1) :: IO (IOUArray Int Double)
sequence_ $ replicate 10 $ multiplyAtAv n u v >> multiplyAtAv n v u
loop 0 0 0 n u v
loop vbv vv i n u v | vbv `seq` vv `seq` i `seq` n `seq` u `seq` v `seq` False = undefined
loop vbv vv i n u v | i >= n = return $! sqrt (vbv/vv)
| otherwise = do ui <- unsafeRead u i
vi <- unsafeRead v i
loop (vbv + ui * vi) (vv + vi * vi) (i+1) n u v
-- return element i,j of infinite matrix A
a i j | i `seq` j `seq` False = undefined
a i j = 1 / fromIntegral (x*(x+1) `div` 2 + i + 1) where x = i+j
-- multiply vector v by matrix A */
multiplyAv n v av | n `seq` v `seq` av `seq` False = undefined
multiplyAv n v av = loop 0
where loop i = when (i < n) $ loop' i 0 0 >>= unsafeWrite av i >> loop (i+1)
loop' i j av | i `seq` j `seq` av `seq` False = undefined
loop' i j av | j >= n = return av
| otherwise = do vj <- v `unsafeRead` j
loop' i (j+1) (av + a i j * vj)
-- multiply vector v by matrix A transposed
multiplyAtv n v atv | n `seq` v `seq` atv `seq` False = undefined
multiplyAtv n v atv = loop 0
where loop i = when (i < n) $ loop' i 0 0 >>= unsafeWrite atv i >> loop (i+1)
loop' i j atvi | j `seq` atvi `seq` False = undefined
loop' i j atvi | j >= n = return atvi
| otherwise = do vj <- v `unsafeRead` j
loop' i (j+1) (atvi + a j i * vj)
-- multiply vector v by matrix A and then by matrix A transposed */
multiplyAtAv n v atav | n `seq` v `seq` atav `seq` False = undefined
multiplyAtAv n v atav = do u <- newArray_ (0,n-1) :: IO (IOUArray Int Double)
multiplyAv n v u >> multiplyAtv n u atav
Current entry
-- The Great Computer Language Shootout
-- http:--shootout.alioth.debian.org/
--
-- Original C contributed by Sebastien Loisel
-- Conversion to C++ by Jon Harrop
-- Conversion to Haskell by Einar Karttunen
import Control.Monad.ST
import Data.Array.Base
import Data.Array.ST
import Numeric
import System
eval_A :: Int -> Int -> Double
eval_A i j = 1 / fromIntegral ((i+j)*(i+j+1) `div` 2 + i + 1)
plusAt :: STUArray s Int Double -> Int -> Double -> ST s ()
plusAt a i v = do o <- unsafeRead a i
unsafeWrite a i (v+o)
eval_A_Times_u :: STUArray s Int Double -> STUArray s Int Double -> ST s ()
eval_A_Times_u u au = outer (snd $ bounds u)
where outer 0 = unsafeWrite au 0 0 >> inner 0 (snd $ bounds u)
outer i = unsafeWrite au i 0 >> inner i (snd $ bounds u) >> outer (i-1)
inner i 0 = unsafeRead u 0 >>= \uj -> plusAt au i (eval_A i 0 * uj)
inner i j = unsafeRead u j >>= \uj -> plusAt au i (eval_A i j * uj) >> inner i (j-1)
eval_At_Times_u :: STUArray s Int Double -> STUArray s Int Double -> ST s ()
eval_At_Times_u u au = outer (snd $ bounds u)
where outer 0 = unsafeWrite au 0 0 >> inner 0 (snd $ bounds u)
outer i = unsafeWrite au i 0 >> inner i (snd $ bounds u) >> outer (i-1)
inner i 0 = unsafeRead u 0 >>= \uj -> plusAt au i (eval_A 0 i * uj)
inner i j = unsafeRead u j >>= \uj -> plusAt au i (eval_A j i * uj) >> inner i (j-1)
eval_AtA_Times_u u v = do w <- newArray (bounds u) 0
eval_A_Times_u u w >> eval_At_Times_u w v
main = do
n <- getArgs >>= return.read.head
let (vBv,vv) = runST (do u <- newArray (0,n-1) 1
v <- newArray (0,n-1) 0
sequence_ $ replicate 10 (eval_AtA_Times_u u v >> eval_AtA_Times_u v u)
vLoop u v n (0, 0))
putStrLn $ showFFloat (Just 9) (sqrt (vBv/vv)) ""
vLoop :: STUArray s Int Double -> STUArray s Int Double -> Int -> (Double,Double) -> ST s (Double,Double)
vLoop u v 0 a = return a
vLoop u v (i+1) (vBv,vv) = vLoop u v i =<< op
where op = do ui <- unsafeRead u i
vi <- unsafeRead v i
return (vBv+(ui*vi),vv+(vi*vi))