multiMat::(Num a)=>[[a]]->[[a]]->[[a]]
multiMat mat1 mat2= foldr(\x y ->x + y) mz ma
                                where
                                        -- a list of matriess
                                        mz = matZero $ len mat1 
                                        ma = zipWith(\x y -> outer x y) (L.transpose mat1) mat2

void multiMat(int m, int n, const double* a, 
              int n1, int h, const double* b, 
              double* const ret){
  for(int k = 0; k < h; k++){
    for(int i = 0; i < m; i++){

      // dot product
      double d = 0.0;
      for(int j = 0; j < n; j++){
        d += a[i * n + j] * b [j * h + k];
      }

      ret[i*h + k]  = d;
    }
  }
}

foreign import ccall "multiMat"
        c_multiMat:: CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> Ptr CDouble -> IO()

f_multiMat :: [[Double]] -> [[Double]] -> IO [[Double]]
f_multiMat cs cx = do
  let fx = fromIntegral
  let (h, w) = dim cs
  let (h', w') = dim cx 
  if w == w' && h == h' then do
        let cs' = map realToFrac $ join cs
        let cx' = map realToFrac $ join cx
        let nBytes = let n = 0 :: CDouble in w * h * sizeOf n
        allocaBytes nBytes $ \pt -> do
          withArrayLen cs' $ \l1 pt1 -> do
                withArrayLen cx' $ \l2 pt2 -> do
                  c_multiMat (fx h) (fx w) pt1 (fx h') (fx w') pt2 pt 
                  arr <- peekArray (w * h) pt
                  return $ fxx w $ map realToFrac arr 
  else do
        error $ "ERROR114, Invalid dimension, w == w' && h == h'" 
  where
        dim x = (h, w) 
          where
                w = case x of
                         [] -> 0
                         vs:_ -> length vs
                h = length x 

multiMatArr :: (Ix a, Ix b, Ix c, Num n) => Array (a, b) n -> Array (b, c) n -> Array (a, c) n
multiMatArr a1 a2 = array resultBnd [ ((i, k), sum [a1 ! (i, j) * a2 ! (j, k) | j <- range (x1, y1)]) 
                                                                           | i <- range (x0, y0), 
                                                                                 k <- range (x1', y1')
                                                                        ]
  where
        ((x0, x1),  (y0,  y1) ) = bounds a1
        ((x0',x1'), (y0', y1')) = bounds a2
        resultBnd | (x1, y1) == (x0', y0') = ((x0, x1'),(y0, y1'))
                          | otherwise = error "ERROR: dim of matrices are NOT compatible"

-- | Perform the inner product of two matrices using provided addition and multiplication functions.
innerProduct :: Num a => (a -> a -> a) -> (b -> c -> a) -> Matrix b -> Matrix c -> Matrix a
innerProduct f g m1 m2
| cols m1 /= rows m2 = error "SIZE ERROR: Incompatible sizes for inner product"
| otherwise = reshape (rows m1) (cols m2) [ getList r c | r <- [1..rows m1], c <- [1..cols m2]]
where
getList r c = foldl1 f (zipWith g (getRow m1 r) (getCol m2 c))
getRow m r = [get m r c | c <- [1..cols m]]
getCol m c = [get m r c | r <- [1..rows m]]

-- | Transpose a matrix (rows become columns, and vice versa).
transpose :: Matrix a -> Matrix a
transpose matrix = reshape (cols matrix) (rows matrix) $ transpose' matrix
where
transpose' m = [get m c r | r <- [1..cols m], c <- [1..rows m]]

mulMat :: Num a => Matrix a -> Matrix a -> Matrix a
mulMat m1  m2 = innerProduct (+) (*) m1 m2