Ticket #5054: TestCase.hs

File TestCase.hs, 4.8 KB (added by arsenm, 3 years ago)

First test case

Line 
1{-# LANGUAGE NoMonomorphismRestriction #-}
2{-# OPTIONS_GHC -W #-}
3
4import Data.Int
5import Data.Packed
6import Data.Packed.ST
7import Control.Applicative
8import Control.Monad
9--import qualified Control.Monad.Parallel as Parallel
10import Control.Monad.ST
11import Foreign.Storable
12import Foreign.Ptr
13import Foreign.Marshal.Utils
14
15import Control.Parallel.Strategies
16--import Data.Vector.Strategies
17
18import Graphics.Plot
19
20inParallel = parMap rwhnf id
21
22zeroMatrix m n = buildMatrix m n (const 0)
23
24twoMatrix m n = buildMatrix m n (const (Value 2))
25
26data ComputeElement = Constant !Double
27                    | Value !Double
28                    deriving (Eq)
29
30-- We don't care about showing if it's constant or not
31instance Show ComputeElement where
32  show (Constant v) = show v
33  show (Value v) = show v
34
35instance Element ComputeElement
36
37isConstant (Constant _) = True
38isConstant _            = False
39
40fromComputeElement (Constant v) = v
41fromComputeElement (Value v) = v
42
43sizeofDouble = sizeOf (undefined :: Double)
44sizeofInt64 = sizeOf (undefined :: Int64)
45
46
47{-
48typedef struct
49{
50    double v;
51    int64_t c;
52} ComputeElement;
53-}
54
55instance Storable ComputeElement where
56  sizeOf _ = sizeofDouble + sizeofInt64
57  alignment _ = 16
58
59  peek p = do v <- peek (castPtr p)
60              c <- peek (castPtr (p `plusPtr` sizeofDouble))
61              return $ if toBool (c :: Int64)
62                         then Constant v
63                         else Value v
64
65  poke p v = do let c :: Int64
66                    c = fromBool (isConstant v)
67                poke (castPtr p) (fromComputeElement v)
68                poke (castPtr p `plusPtr` sizeofDouble) c
69
70
71
72jacobi :: Element a => Int -> Matrix a -> Matrix a
73jacobi n mat = undefined
74  where
75    core = subMatrix (1, 1) (rows mat - 1, cols mat - 1) mat
76
77
78applyComputeElement _ v@(Constant _) = v
79applyComputeElement f (Value v) = Value (f v)
80
81
82writeMatrix' = uncurry . writeMatrix
83readMatrix' = uncurry . readMatrix
84
85
86zeroRho _ _ = 0
87
88--jacobiST :: Storable t => Matrix t -> Matrix ComputeElement
89-- rho :: Double -> Double -> Double
90
91type STComputeMatrix s = STMatrix s ComputeElement
92
93type RelaxationFunction s =  STComputeMatrix s    -- initial matrix
94                          -> STComputeMatrix s -- new matrix
95                          -> Int               -- i
96                          -> Int               -- j
97                          -> ST s Double       -- new element
98
99
100
101
102applyMethod :: RelaxationFunction s -> STComputeMatrix s -> STComputeMatrix s -> Int -> Int -> ST s ()
103applyMethod f mat mat' i j = do
104  c <- readMatrix mat i j
105  u <- f mat mat' i j
106  writeMatrix mat' i j $ if isConstant c
107                           then c
108                           else Value u
109
110{-# INLINE readElement #-}
111readElement mat x y = fromComputeElement <$> readMatrix mat x y
112
113jacobiST :: (Double -> Double -> Double) -> (Double, Double) -> Matrix ComputeElement -> Matrix ComputeElement
114jacobiST rho (rangeX, rangeY) origMat = runST $ do
115  let m = rows origMat
116      n = cols origMat
117
118      dx = rangeX / fromIntegral (m - 1)
119      dy = rangeY / fromIntegral (n - 1)
120      dd = dx * dy
121
122      rs = [1 .. (m - 2)] -- without borders
123      cs = [1 .. (n - 2)]
124
125      evalRho i j = rho (fromIntegral i * dx) (fromIntegral j * dy)
126
127      gaussSeidel f mat mat' i j = do
128        -- Read from old matrix
129        a1 <- readElement mat (i + 1) j
130        a2 <- readElement mat i       (j + 1)
131
132        -- Read from new matrix
133        b1 <- readElement mat' (i - 1) j
134        b2 <- readElement mat' i       (j - 1)
135        let f = evalRho i j
136            u = 0.25 * (a1 + a2 + b1 + b2) + (pi * f * dd)
137        return u
138
139
140      jacobi mat mat' i j = do
141        a <- readElement mat (i + 1) j
142        b <- readElement mat (i - 1) j
143        c <- readElement mat i       (j + 1)
144        d <- readElement mat i       (j - 1)
145
146        let f = evalRho i j
147            u = 0.25 * (a + b + c + d) + (pi * f * dd)
148        return u
149
150      jacobiThings = applyMethod jacobi
151
152      --iterateJacobi mat mat' = sequence_ [jacobiThings mat mat' r c | r <- rs, c <- cs]
153
154      -- faster
155      iterateJacobi mat mat' = sequence_ $ map (uncurry (jacobiThings mat mat')) [(r, c) | r <- rs, c <- cs]
156
157      -- Swap the matrices. Iterations will be an event number, 2 * n
158      iterateNJacobi n mat mat' = replicateM n (iterateJacobi mat mat' >> iterateJacobi mat' mat)
159
160  mat  <- thawMatrix origMat
161  mat' <- thawMatrix origMat
162
163  iterateNJacobi 4000 mat mat'
164
165  freezeMatrix mat'
166
167
168
169
170constLeftBorder v n = fromColumns (border:replicate (n - 1) rest)
171  where border = buildVector n (const (Constant v))
172        rest = buildVector n (const (Value 0))
173
174
175computeElementMatrixToDouble :: Matrix ComputeElement -> Matrix Double
176computeElementMatrixToDouble = liftMatrix (mapVector fromComputeElement)
177
178
179herp = let whee = jacobiST zeroRho (0, 1) (constLeftBorder 100 128)
180       in writeFile "Something.pgm" $ matrixToPGM (computeElementMatrixToDouble whee)
181
182main = herp
183