-- -*- haskell -*-
-- Time-stamp: <Sat Stardate: Stardate: [-27]1960.21 hwloidl>
--
-- Perform a factorial computation, a product over a range, using a 
-- parallel divide-and-conquer paradigm, with a threshhold
-- Haskell98 version.
-- Tested on ghc 7.6.3 (with -threaded extensions)
--
-- To install all necessary library (only needs to be done once):
-- > cabal update
-- > cabal install deepseq
-- > cabal install parallel
--
-- Compile: ghc -O2 -rtsopts -threaded -o parfact_thr parfact.hs
-- or:      ghc -cpp -XBangPatterns -i/home/hwloidl/DISTS/packages/deepseq-1.1.0.0 -i/home/hwloidl/DISTS/packages/parallel-3.2.0.3 -O2 -rtsopts -threaded -o parfact_thr parfact.hs
-- Seq Run: ./parfact_thr 900k 100 +RTS -N1    
-- Par Run: ./parfact_thr 900k 100 +RTS -N7 -sstderr 
-- NOTE: second argument is ignored at the moment: chunksize is calculated based on the number of capabilities
--
-- Batch job of measurements + visualisation:
-- > for ((i=1;i<9;i++)) ; do ./parfact_thr 400k 12 +RTS -N${i} -sstderr 1>>LL 2>&1 ; done
-- > cat LL | sed -ne '/TASKS/H;/Total   time/H;${G;p}' | sed -e 's/^.*-N\([0-9]*\).*$/\1/' | sed -e 's/[.0-9]*.*(\(.*\)$/\1/'  | sed -e 's/s elapsed.*$//' | sed -e '/[.]/a\X' | sed ':a;N;$!ba;s/\n/ /g' | sed -e 's/X/\n/g' | sed -e '1d' > _rt.dat
-- > echo "set term x11; plot '_rt.dat' with lines ; pause 10" | gnuplot

-----------------------------------------------------------------------------

module Main(main) where

import System.Environment(getArgs)
import Control.Parallel 
import Control.Parallel.Strategies
import Control.DeepSeq
import GHC.Conc (numCapabilities)
import Data.Time.Clock (NominalDiffTime, diffUTCTime, getCurrentTime)


-- main function
main = do args <- getArgs       -- read command-line arguments
          let 
            n = readIntegerK (args!!0)  -- size of the interval
            t = readIntegerK (args!!1)  -- threshold/chunksize
            z = n `div` (fromIntegral numCapabilities) 
            res = pfact n                -- result, using naice dnc version
            res_dnc = parfact_dnc 1 n t  -- result, using dnc version
            res_dp = parfact_dp 1 n z    -- result, using data par version
          putStrLn ("Workers: "++(show numCapabilities))
          putStrLn ("Calculated chunk size: "++(show z))
          putStrLn ("Running ...")
          t0 <- getCurrentTime -- take time before doing computation
          -- putStrLn (res_dnc `deepseq` "done (divide-and conquer version)") -- force it only, ignore result
          putStrLn (res_dp `deepseq` "done (data-parallel) ") -- force it only, ignore result
          t1 <- getCurrentTime -- take time after having finished computation
          -- putStrLn ("parallel factorial of " ++ (show n) ++                     " with threshold " ++ (show t) ++ " = " ++ (show res_dnc)) 
          -- putStrLn ("parallel factorial of " ++ (show n) ++                     " with chunk size " ++ (show z) ++ " = " ++ (show res_dp)) 
          putStrLn ("Elapsed time: "++(show $ diffUTCTime t1 t0)++" secs")
          t0 <- getCurrentTime
          putStrLn (res_seq `deepseq` "done (sequential)") -- force it only, ignore result
          -- putStrLn (res_dnc `deepseq` "done (dnc naive) ") -- force it only, ignore result
          t1 <- getCurrentTime
          putStrLn ("Elapsed time: "++(show $ diffUTCTime t1 t0)++" secs")
          -- check the result:
          if (res_dp == res_seq) -- parfact_seq 1 n)
            then putStrLn "++ Result OK"
            else putStrLn "** Result WRONG"

-- read large integer values using shorthand such as 90k
readIntegerK :: String -> Integer
readIntegerK str = case (last str) of
                     'K' -> 1000 * (read (init str) :: Integer)
                     'k' -> 1000 * (read (init str) :: Integer)
                     'M' -> 1000000 * (read (init str) :: Integer)
                     'm' -> 1000000 * (read (init str) :: Integer)
                     'G' -> 1000000000 * (read (init str) :: Integer)
                     'g' -> 1000000000 * (read (init str) :: Integer)
                     _  -> (read str :: Integer)

-- wrapper function
{-
parfact_dnc :: Integer -> Integer -> Integer
parfact_dnc n t = parfact_dnc 1 n t
-}

-- concise version:
-- parfact_seq m n = product [m..n]
-- explicit recursion with tail recursion:
parfact_seq m n = parfact_seq' 1 m n
                  where parfact_seq' acc m n | m>n = acc
                                             | otherwise = parfact_seq' (acc*m) (m+1) n

-- parallel worker function, with thresholding
parfact_dnc :: Integer -> Integer -> Integer -> Integer
parfact_dnc m n t | (n-m) <= t = parfact_seq m n         -- seq version if interval size <= t
             	  | otherwise  = left `par` right `pseq` -- par d&c version
             	                 (left * right)
             	                 where mid = (m + n) `div` 2
             	                       left = parfact_dnc m mid t 
             	                       right = parfact_dnc (mid+1) n t

-- parallel divide-and-conquer, without thresholding (from slides)
pfact :: Integer -> Integer
pfact n = pfact' 1 n

pfact' :: Integer -> Integer -> Integer
pfact' m n  | m == n    = m
            | otherwise = left `par` right `pseq` (left * right)
			  where mid = (m + n) `div` 2
			  	left = pfact' m mid
			  	right = pfact' (mid+1) n

-- data-parallel
parfact_dp :: Integer -> Integer -> Integer -> Integer
parfact_dp m n t | m==n       = m -- (n-m) <= t = parfact_seq m n               -- seq version if interval size <= t
             	 | otherwise  = left `par` right `pseq`       -- data par version
             	                (left * right)
             	                where next = (min (m+t) n)
             	                      left =  parfact_seq m (next-1)
             	                      right = parfact_dp next n t