Shootout/Spectral
< Shootout
Jump to navigation
Jump to search
The printable version is no longer supported and may have rendering errors. Please update your browser bookmarks and please use the default browser print function instead.
A Shootout Entry for the spectral norm test
Timings
Debian Linux x86, n=2,500
||Entry || Time|| ||dons || 9.407|| ||old || 10.564||
Proposed 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))