Algorithmus :

Ich schreibe ein Programm mit CUDA und das Problem ist folgendes:

Zwei Matrizen A (n * 128) und B (m * 128)

Ich nehme die erste Reihe von A und berechne nacheinander den Abstand zwischen diesem Vektor und allen Reihen von B.

Ich schreibe das Ergebnis jeder Distanz in eine Zeile einer Matrix C, so dass das Element C (i, j) von C die Distanz zwischen Zeile i von A und Zeile j von B enthält.

und ich fahre mit der nächsten Reihe von A fort.

Ich habe es so implementiert: Ich habe ein Raster aus (n * m) Blöcken und 128 Threads pro Block. (1 * 128).

FRAGE: Das Programm wird erfolgreich mit den erwarteten Ergebnissen ausgeführt, die Zeitausführung ist jedoch nur etwa fünf- bis zehnmal schneller als bei der One-Threaded-CPU-Version. Also würde ich gerne wissenSo erhöhen Sie die Arbeit pro Thread vor der Reduzierung, um die Leistung zu steigern.

Kernel Code (Original: Nicht optimiert)

 __global__ void EuclideanDistances( float *A, float *B , float *C , int n , int m)
{
    // SIZE is equal to 128
__shared__ float accumResult[SIZE];
float sA;
float sB;

    // MAPPING
int bx = blockIdx.x;  // n
int by = blockIdx.y;  // m
int ty = threadIdx.y; // 128
int tx = threadIdx.x; // 1


sA = A [bx * SIZE + ty];
sB = B [by * SIZE + ty];
__syncthreads();


accumResult[ty] = (sA - sB) * (sA - sB);
__syncthreads();


// Parallel tree-reduction
for (int stride = SIZE/2 ; stride > 0 ; stride >>= 1)
    if (ty < stride)
    {
        accumResult[ty] += accumResult [stride + ty];
          __syncthreads();
    }

    // Writing results to output matrix
if ((threadIdx.y == 0))
    C [bx * m + by] = accumResult[ty];
       __syncthreads();
}

AKTUALISIEREN

Jetzt verwende ich ein anderes Mapping: Anstatt ein Raster vonn durchm Blöcke und ein Block von128 Threads, ich erhöhe die Anzahl der Threads innerhalb eines Blocks, um die Anzahl der Blöcke zu verringern.

Neues Mapping:

Block aus128 durch8 Threads (insgesamt 1024 Threads, das ist die maximale Größe)

Gitter vonn/8 durchm/8 Blöcke

Leider gibt es falsche Ergebnisse).

Optimierter Kernel-Code (zu aktualisieren)

__global__ void EuclideanDistances( float *A, float *B , float *C, int n , int m)
{
    __shared__ float accumResult[SIZE][8];
__shared__ float sA[SIZE][8];
__shared__ float sB[SIZE][8];

int bx = blockIdx.x;  // n / 8
int by = blockIdx.y;  // m / 8
int tx = threadIdx.x; // 8
int ty = threadIdx.y; // 128
int i = bx * tx * SIZE + ty;
int j = by * tx * SIZE + ty;

sA[ty][tx] = A [i];
sB[ty][tx] = B[j];
__syncthreads();


accumResult[ty][tx] = (sA[ty][tx] - sB[ty][tx]) * (sA[ty][tx] - sB[ty][tx]);
__syncthreads();

// Reduction
for (int stride = SIZE/2 ; stride > 0 ; stride>>=1)
    if (ty < stride)
    {
        accumResult[ty][tx] += accumResult [stride + ty][tx];
        __syncthreads();
    }

    C[bx *  m + by] = accumResult[0][tx];
}

HOST CODE (Zuweisungen + Kernelaufrufe)

    int main()
{
     int m = 20000; //MatrixA size : m * SIZE
     int n = 4000;  //MatrixB size : n * SIZE

     srand((unsigned)time(0));

     // Host Allocations
     float *matrixA = (float *) malloc (n * SIZE * sizeof(float));
     for(int i=0; i < n * SIZE; i++)
         matrixA[i] = (float) (rand()%100)+1;

     float *matrixB = (float *) malloc (m * SIZE * sizeof(float));
     for(int i=0; i < m * SIZE; i++)
         matrixB[i] = (float) (rand()%100)+1;

     float *results_kernel1 = (float *) malloc (n * m * sizeof(float));
     float *results_kernel2 = (float *) malloc (n * m * sizeof(float));


     //Device Allocation
     float *d_matrixA;
     float *d_matrixB;
     cudaMalloc((void **)&d_matrixA, n * SIZE * sizeof(float));
     cudaMalloc((void **)&d_matrixB, m * SIZE * sizeof(float));
     cudaMemcpy(d_matrixA , matrixA , n * SIZE * sizeof(float) , cudaMemcpyHostToDevice);
     cudaMemcpy(d_matrixB , matrixB , m * SIZE * sizeof(float) , cudaMemcpyHostToDevice);

     float *d_results_kernel1;
     float *d_results_kernel2;
     cudaMalloc((void **)&d_results_kernel1 , n * m * sizeof(float));
     cudaMalloc((void **)&d_results_kernel2 , n * m * sizeof(float));

     dim3 threads1 (1 , 128);
     dim3 blocks1  (n , m);
     EuclideanDistances1 <<<blocks1 , threads1>>> (d_matrixA , d_matrixB , d_results_kernel1 , n , m);
     cudaDeviceSynchronize();
     cudaMemcpy(results_kernel1 , d_results_kernel1 , n * m *sizeof(float) , cudaMemcpyDeviceToHost);
     cudaFree(d_results_kernel1);

     dim3 threads2 (8 , 128);   // 1024 threads per block (maximum)
     dim3 blocks2  (ceil((float)n/8) , ceil((float)m/8));
     EuclideanDistances2 <<<blocks2 , threads2>>> (d_matrixA , d_matrixB , d_results_kernel2 , n , m);
     cudaDeviceSynchronize();
     cudaMemcpy(results_kernel2 , d_results_kernel2 , n * m *sizeof(float) , cudaMemcpyDeviceToHost);
     cudaFree(d_results_kernel2);

     // Visualising and comparing results
     for (int i = 0 ; i < 50 ; i++)
         std::cout << "kernel1 : " << results_kernel1[i] << "  |  kernel2 : " << results_kernel2[i] << std::endl;

     free(matrixA);
     free(matrixB);
     free(results_kernel1);
     free(results_kernel2);

     return 0;
}

PS: Ich habe CUDA 6.0 mit einer NVIDIA GTX 650 (Rechenkapazität 3.0)