Почему матричное умножение быстрее с Repa, чем с hmatrix?

Интересно, что Data.Array.Repa на самом деле быстрее hmatrix, что неожиданно, поскольку hmatrix реализован с использованием LAPACK. Это потому, что Repa использует unboxed-тип?

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)

время работы с 1 ядром: 7.011s
время работы с 2 ядром: 3,975 с

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)

Продолжительность: 20.714s

Ответы

Ответ 1

Возможно, вы используете неоптимизированную библиотеку LAPACK. На моем компьютере, используя libatlas-base, время работы ~ 0.4 с.

$cat matrixproduct.hs

import Numeric.LinearAlgebra

main = do
    let a = (1000><1000) $ replicate (1000*1000) (1::Double)
        b = konst 1 (1000,1000)
    print $ [email protected]@>(100,100)
    print $ [email protected]@>(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