Много параллельных применений последовательного преобразования в репа

Question

В Repa я хотел бы применить определенный dлинейное преобразование параллельно через самое внутреннее измерение моего массива, т. е. по всем векторам "столбцов".

В общем, такое преобразование может быть выражено в виде матрицы Mи каждая запись M*v это просто скалярное произведение соответствующего ряда M с v, Так что я мог бы просто использовать traverse с функцией, которая вычисляет соответствующее скалярное произведение. Это стоит d^2 Работа.

Тем не менее, мой M является особенным: он допускает последовательный алгоритм линейной работы. Например, M может быть нижней треугольной матрицы с 1по всему нижнему треугольнику. затем M*v это просто вектор частичных сумм v (он же "скан"). Эти суммы могут быть вычислены последовательно очевидным образом, но нужно (i-1)Первая запись результата для вычисления iВступление эффективно. (У меня несколько таких M, все из которых могут быть рассчитаны так или иначе в линейном последовательном времени.)

Я не вижу очевидного способа использования traverse (или любые другие функции Repa), чтобы воспользоваться этим свойством M, Это можно сделать? Это будет довольно расточительно использовать d^2алгоритм работы (даже с большим параллелизмом), когда есть такой быстрый алгоритм линейной работы.

(Я видел некоторые старые сообщения SO (например, здесь), задающие подобные вопросы, но ничего, что вполне соответствует моей ситуации.)

ОБНОВИТЬ

По запросу здесь приведен иллюстративный код для M который вычисляет частичные суммы (как описано выше). Как я и ожидал, время выполнения (работа) растет суперлинейно в dвторой аргумент экстента массива (ext). Это связано с тем, что mulM' только указывает, как вычислить iая запись вывода, независимая от всех остальных записей. Хотя в общем размере массива есть алгоритм линейного времени, я не знаю, как это выразить в Repa.

Интересно, если я удалю строку, которая определяет манифест array' от mainто время выполнения масштабируется только линейно в общем размере массива! Поэтому, когда массивы задерживаются "на всем пути вниз", слияние / оптимизация должны каким-то образом извлекать алгоритм линейной работы, но без какой-либо явной помощи от меня. Это удивительно, но и не очень полезно для меня, потому что на самом деле мне нужно позвонить 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"

Answer 1

Вот возможная модифицированная версия вашего кода с использованием отложенных массивов:

        {-# 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 для получения дополнительной информации о его возможностях и как лучше всего выразить свои вычисления.