{-# LANGUAGE ParallelArrays, ParallelListComp #-}
{-# OPTIONS -fvectorise #-}

module DPH_mmult_wrapper (matMult_wrapper, Matrix_wrapper) where

import qualified Prelude
import Data.Array.Parallel
import Data.Array.Parallel.Prelude.Double as D
import Data.Array.Parallel.Prelude.Int as I

type MMultType = Double
type Matrix = [:[:MMultType:]:]
type MVector = [:MMultType:]
type Matrix_wrapper = PArray (PArray MMultType)

{-# NOINLINE matMult_wrapper #-}
matMult_wrapper :: Matrix_wrapper -> Matrix_wrapper -> Matrix_wrapper
matMult_wrapper mA mB = toPArrayP (mapP toPArrayP (matMult (fromNestedPArrayP mA) (fromNestedPArrayP mB)))

matMult :: Matrix -> Matrix -> Matrix
matMult mA mB = mapP (\row -> mapP (\col -> dotp row col) mB) mA
--matMult mA mB = mapP (\row -> mapP (\col -> dotp row col) (transposeP mB)) mA

dotp :: MVector -> MVector -> MMultType
dotp row col = D.sumP (zipWithP (D.*) row col)

{-
transposeP :: Matrix -> Matrix
transposeP m = 
    let
        h = lengthP m
        w = lengthP (m !: 0)
        rh = I.enumFromToP 0 (h I.- 1)
        rw = I.enumFromToP 0 (w I.- 1)
    in
        if h I.== 0 then [: :]
        else mapP (\y -> mapP (\x -> m !: x !: y) rh) rw
-}