Wydajność matematyczna Haskell w operacji multiply-add

Question

Piszę grę w Haskell, a moja obecna przepustka w interfejsie użytkownika wiąże się z wieloma procesami generowania geometrii. Obecnie koncentruje się na określeniu wydajności jednej konkretnej operacji (C-owski pseudokod):

Vec4f multiplier, addend; 
Vec4f vecList[]; 
for (int i = 0; i < count; i++) 
    vecList[i] = vecList[i] * multiplier + addend; 

Oznacza to, że bog standard pomnożyć-add czterech pływaków, niby rzecz dojrzałych optymalizacji SIMD.

Rezultatem jest przejście do bufora wierzchołkowego OpenGL, więc ostatecznie musi zostać wrzucony do płaskiej macierzy C. Z tego samego powodu, obliczenia powinny prawdopodobnie zostać wykonane na typach C 'float'.

Szukałem biblioteki lub natywnego rozwiązania idiomatycznego, aby szybko zrobić takie rzeczy w Haskell, ale każde rozwiązanie, które wymyśliłem, zdaje się krążyć wokół 2% wydajności (czyli 50x wolniej) w porównaniu z C z GCC z prawymi flagami. Oczywiście, zacząłem od Haskella kilka tygodni temu, więc moje doświadczenie jest ograniczone - dlatego właśnie przybywam do was. Czy ktokolwiek z was może zaproponować szybszą implementację Haskella lub wskazówki do dokumentacji na temat pisania wysokowydajnego kodu Haskella?

Najpierw najnowsze rozwiązanie Haskell (zegary około 12 sekund). Próbowałem rzeczy z hukiem od this SO post, ale to nie miało znaczenia AFAICT. Zastąpienie "multAdd" przez "(\ i v -> v * 4)" spowolniło czas wykonywania do 1,9 sekundy, więc błędy bitowe (i wynikające z nich wyzwania dla automatycznej optymalizacji) nie wydają się zbyt duże.

{-# LANGUAGE BangPatterns #-} 
{-# OPTIONS_GHC -O2 -fvia-C -optc-O3 -fexcess-precision -optc-march=native #-} 

import Data.Vector.Storable 
import qualified Data.Vector.Storable as V 
import Foreign.C.Types 
import Data.Bits 

repCount = 10000 
arraySize = 20000 

a = fromList $ [0.2::CFloat, 0.1, 0.6, 1.0] 
m = fromList $ [0.99::CFloat, 0.7, 0.8, 0.6] 

multAdd :: Int -> CFloat -> CFloat 
multAdd !i !v = v * (m ! (i .&. 3)) + (a ! (i .&. 3)) 

multList :: Int -> Vector CFloat -> Vector CFloat 
multList !count !src 
    | count <= 0 = src 
    | otherwise  = multList (count-1) $ V.imap multAdd src 

main = do 
    print $ Data.Vector.Storable.sum $ multList repCount $ 
     Data.Vector.Storable.replicate (arraySize*4) (0::CFloat) 

Oto, co mam w C. Kod tutaj ma kilka #ifdefs, które uniemożliwiają go kompilowane prosto; przewiń w dół dla sterownika testowego.

#include <stdio.h> 
#include <stdlib.h> 
#include <time.h> 

typedef float v4fs __attribute__ ((vector_size (16))); 
typedef struct { float x, y, z, w; } Vector4; 

void setv4(v4fs *v, float x, float y, float z, float w) { 
    float *a = (float*) v; 
    a[0] = x; 
    a[1] = y; 
    a[2] = z; 
    a[3] = w; 
} 

float sumv4(v4fs *v) { 
    float *a = (float*) v; 
    return a[0] + a[1] + a[2] + a[3]; 
} 

void vecmult(v4fs *MAYBE_RESTRICT s, v4fs *MAYBE_RESTRICT d, v4fs a, v4fs m) { 
    for (int j = 0; j < N; j++) { 
     d[j] = s[j] * m + a; 
    } 
} 

void scamult(float *MAYBE_RESTRICT s, float *MAYBE_RESTRICT d, 
      Vector4 a, Vector4 m) { 
    for (int j = 0; j < (N*4); j+=4) { 
     d[j+0] = s[j+0] * m.x + a.x; 
     d[j+1] = s[j+1] * m.y + a.y; 
     d[j+2] = s[j+2] * m.z + a.z; 
     d[j+3] = s[j+3] * m.w + a.w; 
    } 
} 

int main() { 
    v4fs a, m; 
    v4fs *s, *d; 

    setv4(&a, 0.2, 0.1, 0.6, 1.0); 
    setv4(&m, 0.99, 0.7, 0.8, 0.6); 

    s = calloc(N, sizeof(v4fs)); 
    d = s; 

    double start = clock(); 
    for (int i = 0; i < M; i++) { 

#ifdef COPY 
     d = malloc(N * sizeof(v4fs)); 
#endif 

#ifdef VECTOR 
     vecmult(s, d, a, m); 
#else 
     Vector4 aa = *(Vector4*)(&a); 
     Vector4 mm = *(Vector4*)(&m); 
     scamult((float*)s, (float*)d, aa, mm); 
#endif 

#ifdef COPY 
     free(s); 
     s = d; 
#endif 
    } 
    double end = clock(); 

    float sum = 0; 
    for (int j = 0; j < N; j++) { 
     sum += sumv4(s+j); 
    } 
    printf("%-50s %2.5f %f\n\n", NAME, 
      (end - start)/(double) CLOCKS_PER_SEC, sum); 
} 

Ten skrypt skompiluje i uruchomi testy z wieloma kombinacjami flag gcc. Najlepsza wydajność została osiągnięta przez mój system w systemie cmath-64-natywnego wektora-ograniczenia O3, wynosząc 0,22 sekundy.

import System.Process 
import GHC.IOBase 

cBase = ("cmath", "gcc mult.c -ggdb --std=c99 -DM=10000 -DN=20000") 
cOptions = [ 
      [("32", "-m32"), ("64", "-m64")], 
      [("generic", ""), ("native", "-march=native -msse4")], 
      [("O1", "-O1"), ("O2", "-O2"), ("O3", "-O3")], 
      [("restrict", "-DMAYBE_RESTRICT=__restrict__"), 
       ("norestrict", "-DMAYBE_RESTRICT=")], 
      [("vector", "-DVECTOR"), ("scalar", "")], 
      [("copy", "-DCOPY"), ("nocopy", "")] 
      ] 

-- Fold over the Cartesian product of the double list. Probably a Prelude function 
-- or two that does this, but hey. The 'perm' referred to permutations until I realized 
-- that this wasn't actually doing permutations. ' 
permfold :: (a -> a -> a) -> a -> [[a]] -> [a] 
permfold f z [] = [z] 
permfold f z (x:xs) = concat $ map (\a -> (permfold f (f z a) xs)) x 

prepCmd :: (String, String) -> (String, String) -> (String, String) 
prepCmd (name, cmd) (namea, cmda) = 
    (name ++ "-" ++ namea, cmd ++ " " ++ cmda) 

runCCmd name compileCmd = do 
    res <- system (compileCmd ++ " -DNAME=\\\"" ++ name ++ "\\\" -o " ++ name) 
    if res == ExitSuccess 
     then do system ("./" ++ name) 
       return() 
     else putStrLn $ name ++ " did not compile" 

main = do 
    mapM_ (uncurry runCCmd) $ permfold prepCmd cBase cOptions 

Answer 1

Roman Leschinkskiy odpowiada:

Faktycznie, rdzeń wygląda przeważnie ok do mnie. Używanie unsafeIndex zamiast (!) powoduje, że program jest dwukrotnie szybszy niż fast (see my answer above). Poniższy program jest znacznie szybszy, chociaż (i czystszy, IMO). Podejrzewam, że pozostała różnica między tym i programem C wynika z ogólnego przyzwyczajenia GHC , jeśli chodzi o pływający punkt . HEAD produkuje najlepsze wyniki z NCG i -msse2

Najpierw zdefiniować nowy typ danych vec4:

{-# LANGUAGE BangPatterns #-} 

import Data.Vector.Storable 
import qualified Data.Vector.Storable as V 
import Foreign 
import Foreign.C.Types 

-- Define a 4 element vector type 
data Vec4 = Vec4 {-# UNPACK #-} !CFloat 
       {-# UNPACK #-} !CFloat 
       {-# UNPACK #-} !CFloat 
       {-# UNPACK #-} !CFloat 

Upewnij możemy przechowywać je w tablicy

instance Storable Vec4 where 
    sizeOf _ = sizeOf (undefined :: CFloat) * 4 
    alignment _ = alignment (undefined :: CFloat) 

    {-# INLINE peek #-} 
    peek p = do 
      a <- peekElemOff q 0 
      b <- peekElemOff q 1 
      c <- peekElemOff q 2 
      d <- peekElemOff q 3 
      return (Vec4 a b c d) 
    where 
     q = castPtr p 
    {-# INLINE poke #-} 
    poke p (Vec4 a b c d) = do 
      pokeElemOff q 0 a 
      pokeElemOff q 1 b 
      pokeElemOff q 2 c 
      pokeElemOff q 3 d 
    where 
     q = castPtr p 

wartości i metody tego typu:

a = Vec4 0.2 0.1 0.6 1.0 
m = Vec4 0.99 0.7 0.8 0.6 

add :: Vec4 -> Vec4 -> Vec4 
{-# INLINE add #-} 
add (Vec4 a b c d) (Vec4 a' b' c' d') = Vec4 (a+a') (b+b') (c+c') (d+d') 

mult :: Vec4 -> Vec4 -> Vec4 
{-# INLINE mult #-} 
mult (Vec4 a b c d) (Vec4 a' b' c' d') = Vec4 (a*a') (b*b') (c*c') (d*d') 

vsum :: Vec4 -> CFloat 
{-# INLINE vsum #-} 
vsum (Vec4 a b c d) = a+b+c+d 

multList :: Int -> Vector Vec4 -> Vector Vec4 
multList !count !src 
    | count <= 0 = src 
    | otherwise  = multList (count-1) $ V.map (\v -> add (mult v m) a) src 

main = do 
    print $ Data.Vector.Storable.sum 
      $ Data.Vector.Storable.map vsum 
      $ multList repCount 
      $ Data.Vector.Storable.replicate arraySize (Vec4 0 0 0 0) 

repCount, arraySize :: Int 
repCount = 10000 
arraySize = 20000 

z ghc 6.12.1, -O2 -fasm:

• 1,752

ghc z łbem (26 czerwca) -O2 -fasm -msse2

• 1,708

To wygląda jak najbardziej idiomatyczny sposób na napisanie tablicy Vec4 i uzyskanie najlepszej wydajności ce (11 razy szybciej niż oryginał). (I może to stać się punktem odniesienia dla zaplecza LLVM GHC)

Answer 2

Cóż, to jest lepsze. 3.5s zamiast 14s.

{-# LANGUAGE BangPatterns #-} 
{- 

-- multiply-add of four floats, 
Vec4f multiplier, addend; 
Vec4f vecList[]; 
for (int i = 0; i < count; i++) 
    vecList[i] = vecList[i] * multiplier + addend; 

-} 

import qualified Data.Vector.Storable as V 
import Data.Vector.Storable (Vector) 
import Data.Bits 

repCount, arraySize :: Int 
repCount = 10000 
arraySize = 20000 

a, m :: Vector Float 
a = V.fromList [0.2, 0.1, 0.6, 1.0] 
m = V.fromList [0.99, 0.7, 0.8, 0.6] 

multAdd :: Int -> Float -> Float 
multAdd i v = v * (m `V.unsafeIndex` (i .&. 3)) + (a `V.unsafeIndex` (i .&. 3)) 

go :: Int -> Vector Float -> Vector Float 
go n s 
    | n <= 0 = s 
    | otherwise = go (n-1) (f s) 
    where 
    f = V.imap multAdd 

main = print . V.sum $ go repCount v 
    where 
    v :: Vector Float 
    v = V.replicate (arraySize * 4) 0 
      --^a flattened Vec4f [] 

Co jest lepsze niż to było:

$ ghc -O2 --make A.hs 
[1 of 1] Compiling Main    (A.hs, A.o) 
Linking A ... 

$ time ./A 
516748.13 
./A 3.58s user 0.01s system 99% cpu 3.593 total 

multAdd kompiluje dobrze:

 case readFloatOffAddr# 
       rb_aVn 
       (word2Int# 
        (and# (int2Word# sc1_s1Yx) __word 3)) 
       realWorld# 
     of _ { (# s25_X1Tb, x4_X1Te #) -> 
     case readFloatOffAddr# 
       rb11_X118 
       (word2Int# 
        (and# (int2Word# sc1_s1Yx) __word 3)) 
       realWorld# 
     of _ { (# s26_X1WO, x5_X20B #) -> 
     case writeFloatOffAddr# 
       @ RealWorld 
       a17_s1Oe 
       sc3_s1Yz 
       (plusFloat# 
        (timesFloat# x3_X1Qz x4_X1Te) x5_X20B) 

Jednak robisz 4-elementowy przy mnoży czas w kodzie C , więc będziemy musieli to zrobić bezpośrednio, zamiast udawać, że zapętliliśmy i maskowanie. GCC prawdopodobnie również rozwija pętlę.

Aby uzyskać taką samą wydajność, potrzebujemy wektora mnożącego się (trochę twardego, prawdopodobnie za pośrednictwem backendu LLVM) i rozwinąć pętlę (ewentualnie ją fusing). Odłożę tu na Romana, żeby sprawdzić, czy są inne oczywiste rzeczy.

Jednym z pomysłów może być użycie Vector Vec4 zamiast spłaszczania go.