[project @ 1999-10-06 11:10:40 by simonmar]
authorsimonmar <unknown>
Wed, 6 Oct 1999 11:10:41 +0000 (11:10 +0000)
committersimonmar <unknown>
Wed, 6 Oct 1999 11:10:41 +0000 (11:10 +0000)
Add new subdirectory ghc/ for programs that use ghc extensions.  The
old GHC_ONLY directory used to be populated with micro-tests, all of
which appear to have been moved to ghc/tests.

ghc/matrix is Pascal Serrarens' conversion of a Clean benchmark, that
reputedly runs 4 times slower than 0.29 and 40 times slower than the
clean implementation.  We have some ideas why this might be so:
inlining and specialisation isn't working very well for array
operations, so the bounds checking isn't being compiled down into
decent code.  The changes I just made to the bounds representation in
the various array datatypes improve things by 20% or so.

ghc/Makefile [new file with mode: 0644]
ghc/matrix/Main.hs [new file with mode: 0644]
ghc/matrix/Makefile [new file with mode: 0644]
ghc/matrix/Matrix.hs [new file with mode: 0644]
ghc/matrix/matrix.stdout [new file with mode: 0644]

diff --git a/ghc/Makefile b/ghc/Makefile
new file mode 100644 (file)
index 0000000..ad19d7a
--- /dev/null
@@ -0,0 +1,7 @@
+TOP = ..
+include $(TOP)/mk/boilerplate.mk
+
+SUBDIRS = matrix
+
+include $(TOP)/mk/target.mk
+
diff --git a/ghc/matrix/Main.hs b/ghc/matrix/Main.hs
new file mode 100644 (file)
index 0000000..f9ab7b1
--- /dev/null
@@ -0,0 +1,54 @@
+module Main where
+
+import ST
+import MutableArray
+import Matrix
+
+eps :: Double
+eps = 1.0E-12
+
+size :: Int
+size = 25
+
+size2 :: Int
+size2 = size * size
+
+main = print (runST (do
+   let a = fplMatrix size
+   x0 <- createVector (take size2 (repeat 1))
+   b  <- createVector (take size2 (repeat (0.1)))
+   x  <- cg a x0 b
+   readDoubleArray x 0
+   )) 
+
+cg :: Matrix -> Vector s -> Vector s -> ST s (Vector s)
+cg a x b = do
+   d      <- newVector x
+   h      <- newVector x
+   g      <- newVector x
+   g      <- amx g a x
+   g      <- xmy g b
+   delta0 <- innerProd g g
+   if delta0 <= eps
+      then
+         return x
+      else do
+         d <- negx g d
+         cgLoop a x g delta0 d h
+
+cgLoop :: Matrix -> Vector s -> Vector s -> Double -> Vector s -> Vector s ->
+   ST s (Vector s)
+cgLoop a x g delta0 d h = do
+   h      <- _scc_ "amx" amx h a d
+   ip_dh  <- innerProd d h
+   let tau = delta0 / ip_dh
+   x1     <- _scc_ "saxpy1" saxpy tau d x
+   g      <- _scc_ "saxpy2" saxpy tau h g
+   delta1 <- innerProd g g
+   if delta1 <= eps
+      then
+         return x
+      else do
+         let beta = delta1 / delta0
+         d       <- _scc_ "saxmy" saxmy beta d g
+         cgLoop a x1 g delta1 d h 
diff --git a/ghc/matrix/Makefile b/ghc/matrix/Makefile
new file mode 100644 (file)
index 0000000..b9798ec
--- /dev/null
@@ -0,0 +1,7 @@
+TOP = ../..
+include $(TOP)/mk/boilerplate.mk
+
+SRC_HC_OPTS += -fglasgow-exts
+
+include $(TOP)/mk/target.mk
+
diff --git a/ghc/matrix/Matrix.hs b/ghc/matrix/Matrix.hs
new file mode 100644 (file)
index 0000000..b055960
--- /dev/null
@@ -0,0 +1,159 @@
+module Matrix (
+   Matrix(..),
+   Vector(..),
+   fplMatrix,
+   createVector,
+   newVector,
+   amx,
+   saxpy,
+   saxmy,
+   xmy,
+   negx,
+   innerProd
+) where
+
+import MutableArray
+import ByteArray
+import ST
+
+type Vector s = MutableByteArray s Int
+type Diagonal = ByteArray Int
+type Matrix = (Int, Diagonal, Diagonal, Diagonal)
+
+fplMatrix :: Int -> Matrix
+fplMatrix size = (size, d0, d1, d2)
+   where
+      n = size * size
+      d0 = al n (\i -> 4)
+      d1 = al (n - 1) (\i -> if ((i + 1) `mod` size == 0) then 0 else (-1))
+      d2 = al (n - size) (\i -> -1)
+
+      al n f = runST (do
+         a <- newDoubleArray (0 :: Int, n - 1)
+         a <- al_ 0 a
+         freezeDoubleArray a)
+         where
+            al_ i a
+               | i >= n    = return a
+               | otherwise = do
+                  writeDoubleArray a i (f i)
+                  al_ (i + 1) a
+
+createVector :: [Double] -> ST s (Vector s)
+createVector xs = do 
+   a <- newDoubleArray (0 :: Int, length xs - 1)
+   createVector_ xs 0 a
+   where
+      createVector_ []     i a = return a
+      createVector_ (x:xs) i a = do
+         writeDoubleArray a i x
+         createVector_ xs (i + 1) a
+
+newVector :: Vector s -> ST s (Vector s)
+newVector v = newDoubleArray (p, q)
+   where
+      (p, q) = boundsOfMutableByteArray v
+
+saxmy :: Double -> Vector s -> Vector s -> ST s (Vector s)
+saxmy a x y =
+   saxmy' 0 n x y
+      where 
+         n = div (sizeofMutableByteArray x) 8
+
+         saxmy' :: Int -> Int -> Vector s -> Vector s -> ST s (Vector s)
+         saxmy' i n x y
+            | i >= n    = return x
+            | otherwise = do
+               xe <- readDoubleArray x i
+               ye <- readDoubleArray y i
+               writeDoubleArray x i (a * xe - ye)
+               saxmy' (i + 1) n x y
+
+saxpy :: Double -> Vector s -> Vector s -> ST s (Vector s)
+saxpy a x y =
+   saxpy' 0 n x y
+      where
+         n = div (sizeofMutableByteArray x) 8
+
+         saxpy' :: Int -> Int -> Vector s -> Vector s -> ST s (Vector s)
+         saxpy' i n x y
+            | i >= n    = return y
+            | otherwise = do
+               xe <- readDoubleArray x i
+               ye <- readDoubleArray y i
+               writeDoubleArray y i (a * xe + ye)
+               saxpy' (i + 1) n x y
+
+xmy :: Vector s -> Vector s -> ST s (Vector s)
+xmy x y =
+   xmy' 0 n x y
+      where
+         n = div (sizeofMutableByteArray x) 8
+
+         xmy' :: Int -> Int -> Vector s -> Vector s -> ST s (Vector s)
+         xmy' i n x y
+            | i >= n    = return x
+            | otherwise = do
+               xe <- readDoubleArray x i
+               ye <- readDoubleArray y i
+               writeDoubleArray x i (xe - ye)
+               xmy' (i + 1) n x y 
+
+negx :: Vector s -> Vector s -> ST s (Vector s)
+negx x u =
+   negx' 0 n x u
+      where
+         n = div (sizeofMutableByteArray x) 8
+
+         negx' :: Int -> Int -> Vector s -> Vector s -> ST s (Vector s)
+         negx' i n x u
+            | i >= n    = return u
+            | otherwise = do
+               xe <- readDoubleArray x i
+               writeDoubleArray u i (negate xe)
+               negx' (i + 1) n x u
+
+innerProd :: Vector s -> Vector s -> ST s Double
+innerProd x y =
+   innerProd' 0 0 n x y
+      where
+         n = div (sizeofMutableByteArray x) 8
+
+         innerProd' :: Double -> Int -> Int -> Vector s -> Vector s -> ST s Double
+         innerProd' r i n x y
+            | i >= n    = return r
+            | otherwise = do
+                xe <- readDoubleArray x i
+                ye <- readDoubleArray y i
+                innerProd' (r + xe * ye) (i + 1) n x y
+
+amx :: Vector s -> Matrix -> Vector s -> ST s (Vector s)
+amx u (offset, d0, d1, d2) v = do
+   u <- mul0 0 n d0 v u
+   u <- mul1 0 1 n d1 v u
+   u <- mul1 0 offset n d2 v u
+   return u
+      where
+         n = div (sizeofMutableByteArray v) 8
+
+--         mul0 :: Int -> Int -> Diagonal -> Vector s -> Vector s -> ST s (Vector s)
+         mul0 i (n :: Int) d0 v u 
+            | (i :: Int) >= n    = return u
+            | otherwise = do
+               ve    <- readDoubleArray v i
+               let de = indexDoubleArray d0 i
+               writeDoubleArray u i (de * ve)
+               mul0 (i + 1) n d0 v u
+
+--         mul1 :: Int -> Int -> Int -> Diagonal -> Vector s -> Vector s -> ST s (Vector s)
+         mul1 i1 i2 (n :: Int) d v u
+            | (i2 :: Int) >= n   = return u
+            | otherwise = do
+               let de = indexDoubleArray d i1
+               e1    <- readDoubleArray u i1
+               e2    <- readDoubleArray u i2
+               ve1   <- readDoubleArray v i1
+               ve2   <- readDoubleArray v i2
+               writeDoubleArray u i1 (e1 + de * ve2)
+               writeDoubleArray u i2 (e2 + de * ve1)
+               mul1 (i1 + 1) (i2 + 1) n d v u
diff --git a/ghc/matrix/matrix.stdout b/ghc/matrix/matrix.stdout
new file mode 100644 (file)
index 0000000..02c602f
--- /dev/null
@@ -0,0 +1 @@
+0.18918007884682553