MatMult.hs

{-# OPTIONS -cpp #-}
{-# LANGUAGE BangPatterns #-}

{- Matrix multiplication using a torus (gentleman algorithm) -- FR10 --    -}

{-
RL/JB ParCo2005: eliminate result communication (Maybe-Type)

JB PhD2008: adapt for simple PhD skeleton tests

JB MSR07/2008: modified to use all available toroid skeletons.

JB MSR07/2008: derived a straight-forward GpH program using identical
               helpers and strategies

JB optimised prodEscalar

JB for ghc-6.9: replaced Control.Parallel.Strategies by a workaround
   (reexporting what should work)

-}
        
module Main(main) where

import System.Environment
import Data.List hiding (foldl', foldl1')

import ListAux
import Control.DeepSeq

import Control.Monad.Par

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

type Vector = [Int]
type Matrix = [Vector]

-- main computation, different versions:
mult :: Int -> Matrix -> Matrix -> Int -> [[Maybe Matrix]]
mult 0 m1 m2 _ = 
#ifdef OUTPUT
        [[Just $ multMatricesTr m1 (transpose m2)]]
#else
        rnf (multMatricesTr m1 (transpose m2)) `seq` [[Nothing]]
#endif
mult v m1 m2 c = results
 where results ::  [[Maybe Matrix]]
#ifdef OUTPUT
       results = [[Just computed]]
#else
       results = (rnf computed `seq` [[Nothing]])
#endif
       computed = multMatricesTr m1 m2Tr
--       strats' = strats ++ repeat undef
       m2Tr = transpose m2

prMM' :: (Matrix,Matrix) -> Matrix
prMM' (c,mt) = [[prVV f c | c <- mt]|f <-c]
prVV :: Vector -> Vector -> Int
prVV f c = sum (zipWith (*) f c)
 
shiftRight c [] = []
shiftRight c (xs:xss) = (xs2++xs1):shiftRight (c-1) xss
 where (xs1,xs2) = splitAt c xs

shiftDown c xss = transpose (shiftRight c (transpose xss))

join2 :: Matrix -> Matrix -> Matrix
join2 xs ys = zipWith (++) xs ys
join :: [Matrix] -> Matrix
join xss = foldr join2 (repeat []) xss
       
splitIntoClusters :: Int -> Matrix -> [[Matrix]]
splitIntoClusters c m | c < 1 = splitIntoClusters 1 m
splitIntoClusters c m1 = mss
  where bh  = kPartition (length m1) c
        bhsplit [] [] = []
        bhsplit [] _  = error "some elements left over"
	bhsplit (t:ts) xs = hs : (bhsplit ts rest)
	  	  where (hs,rest) = splitAt t xs   
        ms  = bhsplit bh m1 -- blocks of rows 
        mss = map (colsplit bh) ms
        colsplit [] _  = []
        colsplit (t:ts) rs
         | head rs == [] = []
         | otherwise = (cab:colsplit ts resto)
          where  (cab,resto) = unzip (map (splitAt t) rs)
        
--        mss = map (repartir (length m1 `div` c)) ms
--        repartir c xs
--         | head xs == [] = []
--         | otherwise = (cab:repartir c resto)
--          where  (cab,resto) = unzip (map (splitAt c) xs)

-- helper for splitIntoClusters (formerly bresenham)
kPartition :: Int -> Int -> [Int]
kPartition n k = zipWith (+) ((replicate (n `mod` k) 1) ++ repeat 0)
                             (replicate k (n `div` k))

          
mult' :: Int -> Int -> ((Matrix,Matrix),[Matrix],[Matrix]) -> (Maybe Matrix,[Matrix],[Matrix])
mult' nc nr ((sm1,sm2),sm1s,sm2s) 
#ifdef OUTPUT
    =  (Just result,toRight,toDown)
#else
    =  (rnf result `seq` Nothing ,toRight,toDown)
#endif
  where toRight  = take (nc-1) (sm1:sm1s)
        toDown   = take (nr-1) (sm2':sm2s)
        sm2'     = transpose sm2
        sms      = zipWith multMatricesTr (sm1:sm1s) (sm2':sm2s)
        result = foldl1' addMatrices sms  -- foldr1: not enough demand??

        
addMatrices :: Matrix -> Matrix -> Matrix
addMatrices m1 m2 = zipWith addVectors m1 m2
  where addVectors :: Vector -> Vector -> Vector
        addVectors v1 v2 = zipWith (+) v1 v2

-- Assumes the second matrix has already been transposed        
multMatricesTr :: Matrix -> Matrix -> Matrix
multMatricesTr m1 m2 =
  runPar $ parMap (\row -> [prodEscalar2 row col | col <- m2]) m1

-- JB 2008: a lot faster, directly consuming lists, and tail-recursive (optimised with -O2)
prodEscalar2JB :: Vector -> Vector -> Int
prodEscalar2JB v1 v2 = addProd v1 v2 0
    where addProd :: Vector -> Vector -> Int -> Int
          addProd (v:vs) (w:ws) acc = addProd vs ws (acc + v*w)
          addProd [] [] n = n
          addProd _  _  _ = error "addProd: length does not match"

-- JB 2008: identical when using ghc-6.8.3, avoids bug in ghc-HEAD. Version suggested by SM
prodEscalar2 :: Vector -> Vector -> Int
prodEscalar2 v1 v2 = addProd v1 v2 0
addProd :: Vector -> Vector -> Int -> Int
addProd (v:vs) (w:ws) !acc = addProd vs ws (acc + v*w)
addProd _ _ !n = n


prodEscalar :: Vector -> Vector -> Int
prodEscalar v1 v2 = sum (zipWith (*) v1 v2)


------- foldl, strict in head element
foldl1' :: NFData a => (a->a->a) -> [a] -> a
foldl1' f (x:xs) = foldl' f x xs

foldl'           ::  NFData a => (a -> b -> a) -> a -> [b] -> a
foldl' f a []     = a
foldl' f a (x:xs) = -- whnf, not enough( (foldl' f) $! (f a x)) xs
		    let first = f a x 
		    in rnf first `seq` foldl' f first xs


usage :: String -> String
usage name = "Cannon's algorithm: Usage:\n\t "++
	     name ++ " <matrix size> <version> <blocksPerRow> \n" ++
             "Version selects from " -- ++ show (zip [0..] names)

main = do 
       args <- getArgs
       let l = length args
       if l == 0 then do n <- getProgName
                         putStrLn (usage n)
                         putStrLn "\n *** defaults: size 100, seq. computation ***"
                 else return () --putStrLn "Cannon's algorithm"

       let    size  = if null args then 100 else read (head args)
	      opt   = if length args < 2 then 0 else read (args!!1)
	      chunk = if length args < 3 then 1 
                                         else read (args!!2)
	      a = "Matrices of size " ++ show size ++ 
--                  " with skeleton " ++ ((names++repeat "UNDEF")!!opt) ++
                  " using chunk parameter " ++ show chunk ++ "\n"
	      res = mult opt (mA size) (mB size) chunk
	      b = multMatricesTr (mA size) (transpose (mB size))
       -- putStrLn a
#ifdef OUTPUT
       putStrLn "Output wanted, checking result for correctness..."
       let computed = map (map fromJust) res
           computed' = concat (map join computed)
       printMat computed'
       if (b == computed')  
                 then putStrLn "Correct!"
                 else do putStrLn "WRONG RESULT! Should be"
                         printMat b
#else
       -- putStrLn "No Output, matrix stays distributed."
       putStrLn (show res)
#endif        


         
m1 size = replicate size [1..size]
m2 size = listToListList size [1..size*size]
mA size = if size <= 4000 then m1 size else listToListList size (concat (take 20 (repeat [1..(size*size `div` 20)])))
mB size = if size <= 4000 then m1 size else listToListList size (concat (take 20 (repeat [0,2.. ((size*size) `div` 20)-2])))
listToListList c m 
 | length m <= c = [m]
 | otherwise = c1 : listToListList c resto
  where (c1,resto) = splitAt c m


printMat :: Matrix -> IO ()
printMat m = putStrLn ("Matrix: " ++ (show (length (head m))) 
                       ++ " x " ++ (show $ length m) ++ "\n" 
                       ++ (showMat m))

-- instance Show a => Show (Matrix a) where
showMat m_ = "<<" ++ unlines (map (concatMap (\ x -> show x ++ " ")) m_) ++ ">>"

fromJust :: Maybe a -> a
fromJust (Just x) = x
fromJust Nothing  = error "fromJust"