Rendimiento matemático de Haskell en operación de adición múltiple

Question

Estoy escribiendo un juego en Haskell, y mi pase actual en la interfaz de usuario implica una gran cantidad de generación de geometría de procedimiento. Actualmente estoy enfocado en la identificación de realización de una operación en particular (pseudocódigo C-ish):

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

Es decir, un pantano-estándar de multiplicar-complemento de cuatro flotadores, el tipo de cosas maduro para la optimización SIMD.

El resultado va a un búfer de vértices OpenGL, por lo que tiene que ser arrojado a una matriz plana C eventualmente. Por la misma razón, los cálculos probablemente deberían hacerse en C 'float'.

He buscado una biblioteca o una solución idiomática nativa para hacer este tipo de cosas rápidamente en Haskell, pero cada solución que he encontrado parece rondar el 2% del rendimiento (es decir, 50x) más lento) en comparación con C de GCC con los indicadores de la derecha. Por supuesto, comencé con Haskell hace un par de semanas, por lo que mi experiencia es limitada, y es por eso que me dirijo a ustedes. ¿Alguno de ustedes puede ofrecer sugerencias para una implementación más rápida de Haskell, o sugerencias para la documentación sobre cómo escribir código Haskell de alto rendimiento?

En primer lugar, la solución Haskell más reciente (relojes de aproximadamente 12 segundos). Probé las cosas de patrones de explosión de this SO post, pero no hizo una diferencia AFAICT. Al reemplazar 'multAdd' por '(\ i v -> v * 4)', el tiempo de ejecución se redujo a 1.9 segundos, por lo que las cosas a nivel de bit (y los consiguientes desafíos para la optimización automática) no parecen tener demasiada culpa.

{-# 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) 

Esto es lo que tengo en C. El código aquí tiene algunos #ifdefs la que le impide ser compilado recto-para arriba; desplácese hacia abajo para el controlador de prueba.

#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); 
} 

Este script compilará y ejecutará las pruebas con una serie de combinaciones de banderas gcc. El mejor rendimiento lo obtuvo cmath-64-native-O3-restrict-vector-nocopy en mi sistema, con 0.22 segundos.

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

romana Leschinkskiy responde:

En realidad, el núcleo se ve bien en su mayoría a mí. El uso de insafeIndex en lugar de (!) hace que el programa sea más del doble que rápido (see my answer above). El programa a continuación es mucho más rápido, aunque (y limpiador, IMO). Sospecho que la diferencia restante entre esto y el programa C se debe a la sutileza general de GHC cuando se trata de punto flotante . La CABEZA produce los mejores resultados con la NCG y -msse2

En primer lugar, definir un nuevo tipo de datos 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 

asegurarnos de poder almacenarla en una matriz

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 

Valores y métodos en este tipo:

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 

Con GHC 6.12.1, -O2 -fasm:

• 1,752

con la cabeza GHC (junio 26), -O2 -fasm -msse2

• 1,708

Esto se ve como la forma más idiomática de escribir una matriz Vec4, y obtiene el mejor rendimiento ce (11 veces más rápido que tu original). (Y esto podría convertirse en un punto de referencia para el backend LLVM de GHC)

Answer 2

Bueno, esto es mejor. 3.5s en lugar de 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 [] 

¿Qué es mejor de lo que era:

$ 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 

multisuma compila bien:

 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) 

Sin embargo, estás haciendo 4-elemento a la vez se multiplica en el código C , así que tendremos que hacer eso directamente, en lugar de simularlo mediante el bucle y el enmascaramiento . GCC probablemente esté desenrollando el ciclo también.

Para obtener un rendimiento idéntico, necesitamos el vector multiplicar (un poco difícil, posiblemente a través del servidor de LLVM) y desenrollar el bucle (posiblemente fusionándolo). Voy a ceder a Roman aquí para ver si hay otras cosas obvias.

Una idea podría ser utilizar realmente un Vector Vec4, en lugar de aplanarlo.