haskell - 为什么使用 Repa 进行矩阵乘法比使用 hmatrix 更快？

Question

有趣的是，Data.Array.Repa 实际上比 hmatrix 快，这是出乎意料的，因为 hmatrix 是使用 LAPACK 实现的。这是因为 Repa 使用了未装箱的类型吗？

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.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)

运行时间：20.714s

Answer 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 $ a@@>(100,100)
    print $ b@@>(100,100)
    print $ (a <> b) @@> (900,900)

$ ghc matrixproduct.hs -O

$ 时间 ./matrixproduct

1.0
1.0
1000.0

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