Why is matrix multiplication faster with Repa than with hmatrix?

Question

It's interesting that Data.Array.Repa is actually faster than hmatrix, which is unexpected since hmatrix is implemented using LAPACK. Is this because Repa uses the unboxed type?

import Data.Array.Repa
import Data.Array.Repa.Algorithms.Matrix

main = do
    let
        a = fromListUnboxed (Z:.1000:.1000::DIM2) $ replicate (1000*1000) 1.0 :: Array U DIM2 Double
        b = fromListUnboxed (Z:.1000:.1000::DIM2) $ replicate (1000*1000) 1.0 :: Array U DIM2 Double
    m <- (a `mmultP` b)
    print $ m!(Z:.900:.900)

running time with 1 core: 7.011s
running time with 2 core: 3.975s

import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.LAPACK

main = do
    let
        a = (1000><1000) $ replicate (1000*1000) 1.0
        b = (1000><1000) $ replicate (1000*1000) 1.0
    print $ (a `multiplyR` b) @@> (900,900)

Running time: 20.714s

Answer 1

Perhaps you are using a non-optimized LAPACK library. In my computer, using libatlas-base, the running time is ~0.4s.

$ cat matrixproduct.hs

import Numeric.LinearAlgebra

main = do
    let a = (1000><1000) $ replicate (1000*1000) (1::Double)
        b = konst 1 (1000,1000)
    print $ a@@>(100,100)
    print $ b@@>(100,100)
    print $ (a <> b) @@> (900,900)

$ ghc matrixproduct.hs -O

$ time ./matrixproduct

1.0
1.0
1000.0

real    0m0.331s
user    0m0.512s
sys     0m0.016s