haskell REPA中顺序变换的许多并行应用

of1yzvn4  于 2023-10-19  发布在  其他
关注(0)|答案(1)|浏览(117)

在Repa中,我想在数组的最里面的维度上并行地应用某个d维线性变换,即在所有“列”向量上。
一般来说,这样的变换可以表示为矩阵MM*v的每一项只是M的适当行与v的点积。所以我可以使用traverse和一个函数来计算适当的点积。这花费了d^2的工作。
不过,我的M很特别:它允许线性工作顺序算法。例如,M可能是一个下三角形矩阵,1 s贯穿整个下三角形。那么M*v就是v的部分和的向量(也就是这些和可以以明显的方式顺序计算,但是需要结果的第(i-1)项来有效地计算第i项。(我有几个这样的M,所有这些都可以在线性顺序时间中以一种或另一种方式计算。
我没有看到任何明显的方法来使用traverse(或任何其他Repa函数)来利用M的这个属性。能做到吗?当有这样一个快速的线性运算算法可用时,使用d^2-work算法(即使具有丰富的并行性)也是相当浪费的。
(我看过一些老的SO帖子(例如,here)问类似的问题,但没有一个与我的情况完全匹配。

更新

根据要求,这里有一些说明性的代码,用于计算部分和的M(如上所述)。正如我所料,运行时(工作)在d中超线性增长,这是数组范围的第二个参数(ext)。这是因为mulM'只指定了如何计算输出的第i项,与所有其他项无关。尽管数组的总大小有一个线性时间算法,但我不知道如何在Repa中表达它。
有趣的是,如果我从main中删除定义清单array'的那一行,那么运行时只能线性地扩展数组的总大小!因此,当阵列被延迟“一路向下”时,融合/优化必须以某种方式提取线性功算法,但没有我的任何明确帮助。这很神奇,但对我来说也不是很有用,因为在现实中,我需要在清单数组上调用mulM

{-# LANGUAGE TypeOperators, ScopedTypeVariables, FlexibleContexts #-}

module Main where

import Data.Array.Repa as R

-- multiplication by M across innermost dimension
mulM arr = traverse arr id mulM'
    where mulM' _ idx@(i' :. i) =
              sumAllS $ extract (Z:.0) (Z:.(i+1)) $ slice arr (i' :. All)

ext = Z :. (1000000::Int) :. (10::Int) -- super-linear runtime in 2nd arg
--ext = Z :. (10::Int) :. (1000000::Int) -- takes forever

array = fromFunction ext (\(Z:.j:.i) -> j+i)

main :: IO ()
main = do
  -- apply mulM to a manifest array
  array' :: Array U DIM2 Int <- computeP $ array
  ans :: Array U DIM2 Int <- computeP $ mulM array'
  print "done"
dgtucam1

dgtucam11#

下面是你的代码的一个可能的修改版本,使用延迟数组:

{-# LANGUAGE TypeOperators #-}

  import Data.Array.Repa as R

  mulM :: (Num a, Source r a) => Array r DIM2 a -> Array D DIM2 a
  mulM arr = traverse arr id mulM'
  where
  mulM' _ idx@(i' :. i) =
  sumAllS $ extract (Z:.0) (Z:.(i+1)) $ slice arr (i' :. All)

  ext :: DIM2
  ext = Z :. (1000000::Int) :. (10::Int)

  array :: Array D DIM2 Int
  array = fromFunction ext (\(Z:.j:.i) -> j+i)

  main :: IO ()
  main = do
  let delayedArray :: Array D DIM2 Int
  delayedArray = delay array
  result = computeUnboxedS $ mulM delayedArray
  print "done"

请注意,这是一个简化的示例,您可能需要进一步调整和试验Repa的功能,以完全捕获特定矩阵M所需的线性运算算法。请务必参考Repa的文档,以了解有关其功能以及如何最好地表达计算的更多信息。

相关问题