Плохая производительность с транспонированной и совокупной суммой в репе

Question

Я разработал функцию совокупной суммы, как определено ниже в библиотеке Haskell REPA. Тем не менее, я столкнулся с проблемой при объединении этой функции с операцией транспонирования. Все 3 из следующих операций зайдут за секунду:

cumsum $ cumsum $ cumsum x
transpose $ transpose $ transpose x
transpose $ cumsum x

Однако, если я напишу:

cumsum $ transpose x

Производительность ужасно ухудшается. В то время как каждая отдельная операция в изоляции занимает хорошо секунду на изображении 1920x1080, когда они в совокупности теперь занимают более 30 секунд ...

Есть идеи о том, что может вызвать это? Моя интуиция говорит мне, что это как -то связано с задержкой массивов, не принуждая в нужное время и т. Д., Но у меня еще недостаточно опыта, чтобы отследить это.

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

import Data.Array.Repa as Repa

{-# INLINE indexSlice #-}
indexSlice :: (Shape sh, Elt a) => Int -> Array (sh :. Int) a -> (sh :. Int) -> a
indexSlice from arr (z :. ix) = arr `unsafeIndex` (z :. (ix + from))

{-# INLINE sliceRange #-}
sliceRange :: (Slice sh, Shape sh, Elt a) => Int -> Int -> Array (sh :. Int) a -> Array (sh :. Int) a
sliceRange from to arr = fromFunction (z :. (to - from + 1)) $ indexSlice from arr
    where (z :. _) = extent arr

{-# INLINE cumsum' #-}
cumsum' :: (Slice (SliceShape sh), Slice sh, Shape (FullShape sh), Shape (SliceShape sh), Elt a, Num a) =>
     Array (FullShape sh :. Int) a -> t -> (sh :. Int) -> a
cumsum' arr f (sh :. outer) = Repa.sumAll $ sliceRange 0 outer $ Repa.slice arr (sh :. All)

{-# INLINE cumsum #-}
cumsum :: (FullShape sh ~ sh, Slice sh, Slice (SliceShape sh), Shape sh, Shape (SliceShape sh), Elt a, Num a) =>
    Array (sh :. Int) a -> Array (sh :. Int) a
cumsum arr = Repa.force $ unsafeTraverse arr id $ cumsum' arr

Answer 1

С точки зрения библиотеки, способ отладки, чтобы создать обертку для подозрительной операции, а затем посмотрите на основной код, чтобы увидеть, сработал ли Fusion.

-- Main.hs ---------------------------------------------------
import Solver
import Data.Array.Repa.IO.BMP

main 
 = do   Right img       <- readImageFromBMP "whatever.bmp"
        print $ cumsumBMP img

-- Solver.hs --------------------------------------------------
{-# LANGUAGE TypeOperators, FlexibleContexts, TypeFamilies #-}
module Solver (cumsumBMP) where
import Data.Array.Repa  as Repa
import Data.Word

{- all your defs -}

{-# NOINLINE cumsumBMP #-}
cumsumBMP :: Array DIM3 Word8 -> Array DIM3 Word8
cumsumBMP img = cumsum $ transpose img

Я поместил код «решателя» в отдельный модуль, поэтому нам нужно только пройти через базовый код для определений, которые мы заботимся.

Скомпилировать как:

touch Solver.hs ; ghc -O2 --make Main.hs \
 -ddump-simpl -dsuppress-module-prefixes -dsuppress-coercions  > dump

Перейти к определению cumsumBMP и искать letrec ключевое слово. Searching for letrec это быстрый способ найти внутренние петли.

Не слишком далеко, я вижу это: (слегка переформатировано)

case gen_a1tr
of _ {
  GenManifest vec_a1tv ->
    case sh2_a1tc  `cast` ... of _ { :. sh3_a1iu  sh4_a1iv ->
    case ix'_a1t9  `cast` ... of _ { :. sh1'_a1iz sh2'_a1iA ->
    case sh3_a1iu  `cast` ... of _ { :. sh5_X1n0  sh6_X1n2 ->
    case sh1'_a1iz `cast` ... of _ { :. sh1'1_X1n9 sh2'1_X1nb ->
    case sh5_X1n0             of _ { :. sh7_X1n8   sh8_X1na ->
    ...
    case sh2'1_X1nb           of _ { I# y3_X1nO ->
    case sh4_a1iv             of _ { I# y4_X1nP ->
    case sh2'_a1iA            of _ { I# y5_X1nX ->
    ...
    let { x3_a1x6 :: Int# [LclId]
      x3_a1x6 =
        +#
          (*#
             (+#
                (*#
                   y1_a1iM
                   y2_X1nG)
                y3_X1nO)
             y4_X1nP)
          y5_X1nX } in
    case >=#
           x3_a1x6
           0
    of ...

Стихийное бедствие! А x3_a1x6 Переплет явно выполняет некоторую полезную работу (умножения, дополнения и тому подобное), но оно обернуто в длинную серию операций распаковки, которые также выполняются для каждой итерации цикла. Хуже всего то, что он распаковывает длину и ширину (форму) массива на каждой итерации, и эта информация всегда будет одинаковой. GHC действительно должен выплачивать эти выражения случая из петли, но это еще не. Это случай Выпуск № 4081 на GHC TRAC, что, надеюсь, скоро будет исправлено.

Работа должна применить deepSeqArray до входящего массива. Это ставит спрос на свою стоимость на верхнем уровне (за пределами цикла), что позволяет GHC знать, что можно перемещать корпус, совпадает с дальше. Для такой функции, как cumsumBMP, мы также ожидаем, что входящий массив уже будет проявлен, поэтому мы можем добавить явное совпадение случая для этого:

{-# NOINLINE cumsumBMP #-}
cumsumBMP :: Array DIM3 Word8 -> Array DIM3 Word8
cumsumBMP img@(Array _ [Region RangeAll (GenManifest _)])
  = img `deepSeqArray` cumsum $ transpose img

Соберите снова, внутренняя петля теперь выглядит намного лучше:

letrec {
$s$wfoldlM'_loop_s2mW [...]
  :: Int# -> Word# -> Word# [...]
$s$wfoldlM'_loop_s2mW =
  \ (sc_s2mA :: Int#) (sc1_s2mB :: Word#) ->
    case <=# sc_s2mA a_s2ji of _ {
      False -> sc1_s2mB;
      True ->
        $s$wfoldlM'_loop_s2mW
          (+# sc_s2mA 1)
          (narrow8Word#
             (plusWord#
                sc1_s2mB
                (indexWord8Array#
                   rb3_a2gZ
                   (+#
                      rb1_a2gX
                      (+#
                         (*#
                            (+#
                               (*#
                                  wild19_X1zO
                                  ipv1_X1m5)
                               sc_s2mA)
                            ipv2_X1m0)
                         wild20_X1Ct)))))
    }; } in

Это плотный, хвостовой рекурсивный цикл, который использует только примитивные операции. При условии, что вы компилируете с -fllvm -optlo-O3, нет никаких причин, которые не будут работать так быстро, как эквивалентная программа C.

Хотя при управлении им есть небольшая икота:

desire:tmp benl$ ./Main 
Main: Solver.hs:(50,1)-(51,45): Non-exhaustive patterns in function cumsumBMP

Это просто напоминает нам, что нам нужно заставить массив, прежде чем звонить cumsumBMP.

-- Main.hs ---------------------------------------------------
...
import Data.Array.Repa as Repa
main 
 = do   Right img       <- readImageFromBMP "whatever.bmp"
        print $ cumsumBMP $ Repa.force img

В итоге:

1. Вам нужно добавить немного deepSeqArray и сопоставление рисунков GOOP с вашими функциями верхнего уровня, чтобы обойти текущую беспрепятственность в GHC. Это продемонстрировано окончательной версией cumsumBMP функция выше. Если вы хотите, чтобы GHC -HQ в ближайшее время исправил это, добавьте себя в качестве CC в Выпуск № 4081 на GHC TRAC. Анкет Программы REPA будут намного красивее, когда это будет исправлено.
2. Вам не нужно добавлять GOOP к каждой функции. В этом примере мне не нужно было трогать indexSlice и друзья. Общее правило состоит в том, чтобы добавить GOOP к функциям, которые используют force, fold или же sumAll. Анкет Эти функции создают экземпляры фактических циклов, которые работают над данными массива, то есть они конвертируют задержку массива в манифестное значение.
3. Производительность кода REPA определяется так же, как контекст, в котором он используется как фактический код. Если вы передаете свои функции верхнего уровня задерживаемые массивы, они будут работать очень медленно. Есть больше обсуждения этого в Учебник в REPA.
4. Файлы BMP, прочитанные с помощью библиотеки REPA-IO, не предварительно предварительно, поэтому вам необходимо заставить их перед использованием. Это, наверное, неправильный по умолчанию, поэтому я изменю его в следующей версии.