Ich hatte die Idee einer Warp-basierten Parallelreduktion, da per Definition alle Threads eines Warps synchron sind.

Die Idee war also, dass die Eingabedaten um den Faktor 64 reduziert werden können (jeder Thread reduziert zwei Elemente), ohne dass eine Synchronisierung erforderlich ist.

Wie bei der ursprünglichen Implementierung von Mark Harris wird die Reduzierung auf Blockebene angewendet, und die Daten befinden sich im gemeinsam genutzten Speicher.http://gpgpu.org/static/sc2007/SC07_CUDA_5_Optimization_Harris.pdf

Ich habe einen Kernel erstellt, um seine Version und meine Warp-basierte Version zu testen.
Der Kernel selbst speichert BLOCK_SIZE-Elemente völlig identisch im gemeinsamen Speicher und gibt das Ergebnis an seinem eindeutigen Blockindex in einem Ausgabearray aus.

Der Algorithmus selbst funktioniert gut. Getestet mit einer ganzen Reihe von Einsen, um das "Zählen" zu testen.

Funktionskörper der Implementierungen:

/**
 * Performs a parallel reduction with operator add 
 * on the given array and writes the result with the thread 0
 * to the given target value
 *
 * @param inValues T* Input float array, length must be a multiple of 2 and equal to blockDim.x
 * @param targetValue float 
 */
__device__ void reductionAddBlockThread_f(float* inValues,
    float &outTargetVar)
{
    // code of the below functions
}

1. Umsetzung seiner Version:

if (blockDim.x >= 1024 && threadIdx.x < 512)
    inValues[threadIdx.x] += inValues[threadIdx.x + 512];
__syncthreads();
if (blockDim.x >= 512 && threadIdx.x < 256)
    inValues[threadIdx.x] += inValues[threadIdx.x + 256];
__syncthreads();
if (blockDim.x >= 256 && threadIdx.x < 128)
    inValues[threadIdx.x] += inValues[threadIdx.x + 128];
__syncthreads();
if (blockDim.x >= 128 && threadIdx.x < 64)
    inValues[threadIdx.x] += inValues[threadIdx.x + 64];
__syncthreads();

//unroll last warp no sync needed
if (threadIdx.x < 32)
{
    if (blockDim.x >= 64) inValues[threadIdx.x] += inValues[threadIdx.x + 32];
    if (blockDim.x >= 32) inValues[threadIdx.x] += inValues[threadIdx.x + 16];
    if (blockDim.x >= 16) inValues[threadIdx.x] += inValues[threadIdx.x + 8];
    if (blockDim.x >= 8) inValues[threadIdx.x] += inValues[threadIdx.x + 4];
    if (blockDim.x >= 4) inValues[threadIdx.x] += inValues[threadIdx.x + 2];
    if (blockDim.x >= 2) inValues[threadIdx.x] += inValues[threadIdx.x + 1];

    //set final value
    if (threadIdx.x == 0)
        outTargetVar = inValues[0];
}

Ressourcen:

4 Sync-Threads verwendet
12 if-Anweisungen verwendet
11 Lese-, Add- und Schreibvorgänge
1 abschließender Schreibvorgang
5 Nutzung registrieren

Performance:

Durchschnitt aus fünf Testläufen: ~ 19,54 ms

2. Warp-basierter Ansatz: (Gleicher Funktionskörper wie oben)

/*
 * Perform first warp based reduction by factor of 64
 *
 * 32 Threads per Warp -> LOG2(32) = 5
 *
 * 1024 Threads / 32 Threads per Warp = 32 warps
 * 2 elements compared per thread -> 32 * 2 = 64 elements per warp
 *
 * 1024 Threads/elements divided by 64 = 16
 * 
 * Only half the warps/threads are active
 */
if (threadIdx.x < blockDim.x >> 1)
{
    const unsigned int warpId = threadIdx.x >> 5;
    // alternative threadIdx.x & 31
    const unsigned int threadWarpId = threadIdx.x - (warpId << 5);
    const unsigned int threadWarpOffset = (warpId << 6) + threadWarpId;

    inValues[threadWarpOffset] += inValues[threadWarpOffset + 32];
    inValues[threadWarpOffset] += inValues[threadWarpOffset + 16];
    inValues[threadWarpOffset] += inValues[threadWarpOffset + 8];
    inValues[threadWarpOffset] += inValues[threadWarpOffset + 4];
    inValues[threadWarpOffset] += inValues[threadWarpOffset + 2];
    inValues[threadWarpOffset] += inValues[threadWarpOffset + 1];
}

// synchronize all warps - the local warp result is stored
// at the index of the warp equals the first thread of the warp
__syncthreads();

// use first warp to reduce the 16 warp results to the final one
if (threadIdx.x < 8)
{
    // get first element of a warp
    const unsigned int warpIdx = threadIdx.x << 6;

    if (blockDim.x >= 1024) inValues[warpIdx] += inValues[warpIdx + 512];
    if (blockDim.x >= 512) inValues[warpIdx] += inValues[warpIdx + 256];
    if (blockDim.x >= 256) inValues[warpIdx] += inValues[warpIdx + 128];
    if (blockDim.x >= 128) inValues[warpIdx] += inValues[warpIdx + 64];

    //set final value
    if (threadIdx.x == 0)
        outTargetVar = inValues[0];
}

Ressourcen:

1 Sync-Thread verwendet
7 if-Anweisungen
10 Lese- und Schreiboperationen
1 abschließender Schreibvorgang
5 Nutzung registrieren

5 Bitverschiebungen
1 hinzufügen
1 Unter

Performance:

Durchschnitt aus fünf Testläufen: ~ 20,82 ms

Testen Sie beide Kernel mehrmals auf aGeforce 8800 GT 512 mb mit256 mb von Float-Werten. Und laufender Kernel mit256 Threads pro Block (100% Belegung).

Die auf Warp basierende Version ist ~1.28 Millisekunden langsamer.

Wenn zukünftige Karten größere Blockgrößen erlauben, würde der Warp-basierte Ansatz immer noch keine weitere Sync-Anweisung benötigen, da das Maximum 4096 ist, was auf 64 reduziert wird, was durch das endgültige Warp auf 1 reduziert wird

Warum geht es nicht schneller? Oder wo liegt der Fehler in der Idee, Kernel?

Aus Ressourcennutzung sollte der Warp-Ansatz voraus sein?

Edit1: Der Kernel wurde dahingehend korrigiert, dass nur die Hälfte der Threads aktiv ist, was nicht zu unzulässigen Lesevorgängen führt. Es wurden neue Leistungsdaten hinzugefügt