Transposition de la matrice CUDA avec la mémoire partagée

Question

Je dois implémenter une fonction de transposition matricielle sur un GPU utilisant la mémoire partagée. Je l'ai fait d'une manière simple sans mémoire partagée qui fonctionne bien et aussi une tentative avec SM. Mais malheureusement, le calcul n'est pas correct et je ne peux pas comprendre pourquoi. Un exemple de travail complet peut être trouvé here et au bas de cette question.

EDIT 1

je sais en outre que le premier index du résultat où j'ai une valeur incorrecte est index 32 (de la matrice méplat, de sorte que matr[0][32] dans un mode à deux dimensions).

S'il y a plus d'informations, je les préviendrai volontiers.

Un court extrait du code entier qui ressemble à la fonction de travail est non énumérés ci-dessous:

__global__ void notSoNaivaTransKernel(float *matrB, float *matrA, const int width, 
    const int height, const int nreps) 
{ 
    __shared__ float tile[TILE_DIM][TILE_DIM + 1]; 
    int blockIdx_y = blockIdx.x; 
    int blockIdx_x = (blockIdx.x + blockIdx.y) % gridDim.x; 
    int xIndex = blockIdx_x * TILE_DIM + threadIdx.x; 
    int yIndex = blockIdx_y * TILE_DIM + threadIdx.y; 
    int index_in = xIndex + (yIndex)* width; 

    xIndex = blockIdx_y * TILE_DIM + threadIdx.x; 
    yIndex = blockIdx_x * TILE_DIM + threadIdx.y; 
    int index_out = xIndex + (yIndex)* height; 

    int r, i; 
#pragma unroll 
    for (r = 0; r < nreps; r++) 
    { 
#pragma unroll 
     for (i = 0; i < TILE_DIM; i += BLOCK_ROWS) 
      tile[threadIdx.y + i][threadIdx.x] = matrA[index_in + i * width]; 

     __syncthreads(); 

#pragma unroll 
     for (i = 0; i < TILE_DIM; i += BLOCK_ROWS) 
      if (index_in + i * width < width * height) 
       matrB[index_out + i * height] = tile[threadIdx.x][threadIdx.y + i]; 
    } 
} 

La sortie ressemble à ceci:

Avg. CPU Transpose Time: 0.106048 ms, Bandwidth: 3.771873 GB/s 

Avg. GPU Naive Trans Time: 0.009871 ms, bandwidth: 40.520836 GB/s 
    Correct: 50000, Wrong: 0 

Avg. GPU Trans with SM Time: 0.007598 ms, bandwidth: 52.643482 GB/s 
    Correct: 12352, Wrong: 37648 

Voici le plein exemple de travail. J'ai enlevé la plupart du code inutile de sorte qu'il est moins bourré:

#include "cuda_runtime.h" 
#include "device_functions.h" 
#include "device_launch_parameters.h" 

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

#define TILE_DIM 32 
#define BLOCK_ROWS 8 
#define BLOCK_COLS 32 

cudaError_t matrMagicCuda(float *matrB, float *matrA, const int width, const int height, const int nreps, const int operation); 
void cpuMatrTrans(float *matrB, float *matrA, const int width, const int height, const int nreps); 
__global__ void naiveTransKernel(float *matrB, float *matrA, const int width, const int height, const int nreps); 
__global__ void notSoNaivaTransKernel(float *matrB, float *matrA, const int width, const int height, const int nreps); 

int main() 
{ 
    int i, width, height, nreps, size, wrong, correct; 
    double cpuTime, cpuBandwidth; 
    cudaError_t cudaStatus; 

    float *matrA, *matrATC, *matrATG, *matrAC; 

    srand(time(NULL)); 

    nreps = 10000; 
    width = 500; 
    height = 100; 
    size = width * height; 

    matrA = (float*)malloc(size * sizeof(float)); // matrix A 
    matrAC = (float*)malloc(size * sizeof(float)); // matrix A copied 
    matrATC = (float*)malloc(size * sizeof(float)); // matrix A transposed by CPU 
    matrATG = (float*)malloc(size * sizeof(float)); // matrix A transposed by GPU 

    for (i = 0; i < size; i++) 
    { 
     matrA[i] = (float)i; 
    } 

    auto start = std::chrono::high_resolution_clock::now(); 

    //CPU Transpose 
    cpuMatrTrans(matrATC, matrA, width, height, nreps); 

    auto end = std::chrono::high_resolution_clock::now(); 

    std::chrono::duration<double> diff = end - start; 
    cpuTime = (diff.count() * 1000)/nreps; 
    cpuBandwidth = (sizeof(float) * size * 2)/(cpuTime * 1000000);//scaling from ms to s and B to GB doen implicitly, shortened in fraction, times two for read and write 
    printf("Avg. CPU Transpose Time: %f ms, Bandwidth: %f GB/s\n\n", cpuTime, cpuBandwidth); 

    correct = 0; 
    wrong = 0; 

    //Naive transpose 
    matrMagicCuda(matrATG, matrA, width, height, nreps, 1); 

    //Check if calc was correct 
    for (i = 0; i < size; i++) 
    { 
     if (matrATC[i] != matrATG[i]) 
     { 
      /*printf("ERROR - %d - ATC:%f - ATG:%f\n\n", i, matrATC[i], matrATG[i]); 
      return;*/ 
      wrong++; 
     } 
     else 
     { 
      correct++; 
     } 
    } 

    printf("\tCorrect: %d, Wrong: %d\n\n", correct, wrong); 
    correct = 0; 
    wrong = 0; 

    //Transpose with shared memory 
    matrMagicCuda(matrATG, matrA, width, height, nreps, 2); 

    //Check if calc was correct 
    for (i = 0; i < size; i++) 
    { 
     if (matrATC[i] != matrATG[i]) 
     { 
      /*printf("ERROR - %d - ATC:%f - ATG:%f\n\n", i, matrATC[i], matrATG[i]); 
      return;*/ 
      wrong++; 
     } 
     else 
     { 
      correct++; 
     } 
    } 

    //printf("\tTranspose with SM on GPU was executed correctly.\n\n"); 
    printf("\tCorrect: %d, Wrong: %d\n\n", correct, wrong); 
    correct = 0; 
    wrong = 0; 

    // cudaDeviceReset must be called before exiting in order for profiling and 
    // tracing tools such as Nsight and Visual Profiler to show complete traces. 
    cudaStatus = cudaDeviceReset(); 
    if (cudaStatus != cudaSuccess) 
    { 
     fprintf(stderr, "cudaDeviceReset failed!\n"); 
     return 1; 
    } 

    return 0; 
} 

cudaError_t matrMagicCuda(float *matrB, float *matrA, const int width, const int height, const int nreps, const int operation) 
{ 
    float elapsed = 0; 
    float *dev_matrA = 0; 
    float *dev_matrB = 0; 
    cudaError_t cudaStatus; 
    dim3 dim_grid, dim_block; 
    double gpuBandwidth; 

    int size = width * height; 

    dim_block.x = TILE_DIM; 
    dim_block.y = BLOCK_ROWS; 
    dim_block.z = 1; 

    dim_grid.x = (width + TILE_DIM - 1)/TILE_DIM; 
    dim_grid.y = (height + TILE_DIM - 1)/TILE_DIM; 
    dim_grid.z = 1; 

    // Choose which GPU to run on, change this on a multi-GPU system. 
    cudaStatus = cudaSetDevice(0); 
    if (cudaStatus != cudaSuccess) 
    { 
     fprintf(stderr, "cudaSetDevice failed! Do you have a CUDA-capable GPU installed?"); 
     goto Error; 
    } 

    // Allocate GPU buffers for three matrix 
    cudaStatus = cudaMalloc((void**)&dev_matrA, size * sizeof(float)); 
    if (cudaStatus != cudaSuccess) 
    { 
     fprintf(stderr, "cudaMalloc failed!"); 
     goto Error; 
    } 

    cudaStatus = cudaMalloc((void**)&dev_matrB, size * sizeof(float)); 
    if (cudaStatus != cudaSuccess) 
    { 
     fprintf(stderr, "cudaMalloc failed!"); 
     goto Error; 
    } 

    // Copy input matrix from host memory to GPU buffers. 
    cudaStatus = cudaMemcpy(dev_matrA, matrA, size * sizeof(float), cudaMemcpyHostToDevice); 
    if (cudaStatus != cudaSuccess) 
    { 
     fprintf(stderr, "cudaMemcpy failed!"); 
     goto Error; 
    } 

    cudaEvent_t start, stop; 
    cudaEventCreate(&start); 
    cudaEventCreate(&stop); 

    switch (operation) 
    { 
     case(1): 
     { 
      cudaEventRecord(start); 
      // Launch a kernel on the GPU with one thread for each element. 
      naiveTransKernel << <dim_grid, dim_block >> >(dev_matrB, dev_matrA, width, height, nreps); 

      cudaEventRecord(stop); 
      cudaEventSynchronize(stop); 

      cudaEventElapsedTime(&elapsed, start, stop); 
      cudaEventDestroy(start); 
      cudaEventDestroy(stop); 

      elapsed /= nreps; 

      gpuBandwidth = (sizeof(float) * size * 2)/(elapsed * 1000000);//scaling from ms to s and B to GB doen implicitly, shortened in fraction, times two for read and write 
      printf("Avg. GPU Naive Trans Time: %f ms, bandwidth: %f GB/s\n", elapsed, gpuBandwidth); 

      break; 
     } 

     case(2): 
     { 
      cudaEventRecord(start); 
      // Launch a kernel on the GPU with one thread for each element. 
      notSoNaivaTransKernel << <dim_grid, dim_block >> >(dev_matrB, dev_matrA, width, height, nreps); 

      cudaEventRecord(stop); 
      cudaEventSynchronize(stop); 

      cudaEventElapsedTime(&elapsed, start, stop); 
      cudaEventDestroy(start); 
      cudaEventDestroy(stop); 

      elapsed /= nreps; 

      gpuBandwidth = (sizeof(float) * size * 2)/(elapsed * 1000000);//scaling from ms to s and B to GB doen implicitly, shortened in fraction, times two for read and write 
      printf("Avg. GPU Trans with SM Time: %f ms, bandwidth: %f GB/s\n", elapsed, gpuBandwidth); 

      break; 
     } 

    default: 
     printf("No matching opcode was found.\n"); 
    } 

    // Check for any errors launching the kernel 
    cudaStatus = cudaGetLastError(); 
    if (cudaStatus != cudaSuccess) 
    { 
     fprintf(stderr, "Kernel launch failed: %s\n", cudaGetErrorString(cudaStatus)); 
     goto Error; 
    } 

    // cudaDeviceSynchronize waits for the kernel to finish, and returns 
    // any errors encountered during the launch. 
    cudaStatus = cudaDeviceSynchronize(); 
    if (cudaStatus != cudaSuccess) 
    { 
     fprintf(stderr, "cudaDeviceSynchronize returned error code %d after launching Kernel!\n", cudaStatus); 
     goto Error; 
    } 

    // Copy output matrix from GPU buffer to host memory. 
    cudaStatus = cudaMemcpy(matrB, dev_matrB, size * sizeof(float), cudaMemcpyDeviceToHost); 
    if (cudaStatus != cudaSuccess) 
    { 
     fprintf(stderr, "cudaMemcpy failed!"); 
     goto Error; 
    } 

Error: 
    cudaFree(dev_matrB); 
    cudaFree(dev_matrA); 

    return cudaStatus; 
} 

void cpuMatrTrans(float *matrB, float *matrA, const int width, const int height, const int nreps) 
{ 
    int i, j, r; 

#pragma unroll 
    for (r = 0; r < nreps; r++) 
#pragma unroll 
     for (i = 0; i < height; i++) 
#pragma unroll 
      for (j = 0; j < width; j++) 
       matrB[j * height + i] = matrA[i * width + j]; 
} 

__global__ void naiveTransKernel(float *matrB, float *matrA, const int width, const int height, const int nreps) 
{ 
    int i, r; 
    int row = blockIdx.x * TILE_DIM + threadIdx.x; 
    int col = blockIdx.y * TILE_DIM + threadIdx.y; 
    int index_in = row + width * col; 
    int index_out = col + height * row; 

#pragma unroll 
    for (r = 0; r < nreps; r++) 
#pragma unroll 
     for (i = 0; i < TILE_DIM; i += BLOCK_ROWS) 
      if (index_in + i * width < width * height) 
       matrB[index_out + i] = matrA[index_in + i * width]; 
} 

__global__ void notSoNaivaTransKernel(float *matrB, float *matrA, const int width, const int height, const int nreps) 
{ 
    __shared__ float tile[TILE_DIM][TILE_DIM + 1]; 
    int blockIdx_y = blockIdx.x; 
    int blockIdx_x = (blockIdx.x + blockIdx.y) % gridDim.x; 
    int xIndex = blockIdx_x * TILE_DIM + threadIdx.x; 
    int yIndex = blockIdx_y * TILE_DIM + threadIdx.y; 
    int index_in = xIndex + (yIndex)* width; 

    xIndex = blockIdx_y * TILE_DIM + threadIdx.x; 
    yIndex = blockIdx_x * TILE_DIM + threadIdx.y; 
    int index_out = xIndex + (yIndex)* height; 

    int r, i; 
#pragma unroll 
    for (r = 0; r < nreps; r++) 
    { 
#pragma unroll 
     for (i = 0; i < TILE_DIM; i += BLOCK_ROWS) 
      tile[threadIdx.y + i][threadIdx.x] = matrA[index_in + i * width]; 

     __syncthreads(); 

#pragma unroll 
     for (i = 0; i < TILE_DIM; i += BLOCK_ROWS) 
      if (index_in + i * width < width * height) 
       matrB[index_out + i * height] = tile[threadIdx.x][threadIdx.y + i]; 
    } 
} 

Answer 1

Il y avait un certain nombre de problèmes avec ce code. Je ne suis pas sûr de pouvoir tous les couvrir.

Le problème le plus important était peut-être le manque de (et le manque de compréhension) d'un bon contrôle de fil 2D. Votre algorithme crée une grille de threads plus grande dans les deux dimensions que la taille du problème. Cela crée des threads logiques en dehors de de la dimension de votre matrice, dans les deux dimensions.

Vous avez essayé de créer un contrôle de fil 2D comme ceci:

 if (index_in + i * width < width * height) 

Cela ne fonctionnera pas. Supposons que j'ai une matrice 3x3 et une grille de fils 4x4. Le thread à (3,0) est clairement hors-limites pour votre matrice, mais passerait votre vérification de thread 2D.

Dans cette situation, une vérification de thread appropriée doit tester chaque dimension séparément, pas en tant que produit. Notez que cette erreur logique existe également dans votre noyau de transposition "naïf", et vous pouvez confirmer que si vous exécutez votre code avec cuda-memcheck. Il indiquera les erreurs d'accès hors-limites dans ce noyau, même s'il semble fonctionner correctement.

Il y avait également divers autres problèmes. La plupart d'entre elles concernaient l'indexation dans votre noyau de mémoire partagée. Il n'est pas clair pour moi que vous avez compris la manipulation d'indexation nécessaire pour un shared memory transpose. Il y a deux indexation distincts transpose que nous devons faire dans ce cas:

1. Transposer le bloc (tuiles) indices
2. Transposer les indices de fil

Transposition des indices de fil se fait lors de la lecture/écriture de mémoire partagée.Vous avez cela correctement expliqué avec une inversion de l'utilisation de threadIdx.x et threadIdx.y pour la lecture/écriture de la mémoire partagée. Mais aussi près que je peux dire, votre génération d'index pour l'inversion des indices de blocs (dont l'inversion est utilisée pendant la lecture/écriture à global mémoire) a été brisée. C'était l'autre problème majeur à régler.

Le code suivant a ces questions et d'autres problèmes résolus, et semble fonctionner correctement pour moi:

$ cat t33.cu 
#include "cuda_runtime.h" 
#include "device_functions.h" 
#include "device_launch_parameters.h" 

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

#define TILE_DIM 32 
#define BLOCK_ROWS 8 
#define BLOCK_COLS 32 

cudaError_t matrMagicCuda(float *matrB, float *matrA, const int width, const int height, const int nreps, const int operation); 
void cpuMatrTrans(float *matrB, float *matrA, const int width, const int height, const int nreps); 
__global__ void naiveTransKernel(float *matrB, float *matrA, const int width, const int height, const int nreps); 
__global__ void notSoNaivaTransKernel(float *matrB, float *matrA, const int width, const int height, const int nreps); 

int main() 
{ 
    int i, width, height, nreps, size, wrong, correct; 
    double cpuTime, cpuBandwidth; 
    cudaError_t cudaStatus; 

    float *matrA, *matrATC, *matrATG, *matrAC; 

    srand(time(NULL)); 

    nreps = 10000; 
    width = 500; 
    height = 100; 


    size = width * height; 

    matrA = (float*)malloc(size * sizeof(float)); // matrix A 
    matrAC = (float*)malloc(size * sizeof(float)); // matrix A copied 
    matrATC = (float*)malloc(size * sizeof(float)); // matrix A transposed by CPU 
    matrATG = (float*)malloc(size * sizeof(float)); // matrix A transposed by GPU 

    for (i = 0; i < size; i++) 
    { 
     matrA[i] = (float)i; 
    } 

    auto start = std::chrono::high_resolution_clock::now(); 

    //CPU Transpose 
    cpuMatrTrans(matrATC, matrA, width, height, nreps); 

    auto end = std::chrono::high_resolution_clock::now(); 

    std::chrono::duration<double> diff = end - start; 
    cpuTime = (diff.count() * 1000)/nreps; 
    cpuBandwidth = (sizeof(float) * size * 2)/(cpuTime * 1000000);//scaling from ms to s and B to GB doen implicitly, shortened in fraction, times two for read and write 
    printf("Avg. CPU Transpose Time: %f ms, Bandwidth: %f GB/s\n\n", cpuTime, cpuBandwidth); 

    correct = 0; 
    wrong = 0; 

    //Naive transpose 
    memset(matrATG, 0, size*sizeof(float)); 
    matrMagicCuda(matrATG, matrA, width, height, nreps, 1); 

    //Check if calc was correct 
    for (i = 0; i < size; i++) 
    { 
     if (matrATC[i] != matrATG[i]) 
     { 
      /*printf("ERROR - %d - ATC:%f - ATG:%f\n\n", i, matrATC[i], matrATG[i]); 
      return;*/ 
      wrong++; 
     } 
     else 
     { 
      correct++; 
     } 
    } 

    printf("\tCorrect: %d, Wrong: %d\n\n", correct, wrong); 
    correct = 0; 
    wrong = 0; 

    //Transpose with shared memory 
    memset(matrATG, 0, size*sizeof(float)); 
    matrMagicCuda(matrATG, matrA, width, height, nreps, 2); 

    //Check if calc was correct 
    for (i = 0; i < size; i++) 
    { 
     if (matrATC[i] != matrATG[i]) 
     { 
      /*printf("ERROR - %d - ATC:%f - ATG:%f\n\n", i, matrATC[i], matrATG[i]); 
      return;*/ 
      wrong++; 
     } 
     else 
     { 
      correct++; 
     } 
    } 

    //printf("\tTranspose with SM on GPU was executed correctly.\n\n"); 
    printf("\tCorrect: %d, Wrong: %d\n\n", correct, wrong); 
    correct = 0; 
    wrong = 0; 

    // cudaDeviceReset must be called before exiting in order for profiling and 
    // tracing tools such as Nsight and Visual Profiler to show complete traces. 
    cudaStatus = cudaDeviceReset(); 
    if (cudaStatus != cudaSuccess) 
    { 
     fprintf(stderr, "cudaDeviceReset failed!\n"); 
     return 1; 
    } 

    return 0; 
} 

cudaError_t matrMagicCuda(float *matrB, float *matrA, const int width, const int height, const int nreps, const int operation) 
{ 
    float elapsed = 0; 
    float *dev_matrA = 0; 
    float *dev_matrB = 0; 
    cudaError_t cudaStatus; 
    dim3 dim_grid, dim_block; 
    double gpuBandwidth; 

    int size = width * height; 

    dim_block.x = TILE_DIM; 
    dim_block.y = BLOCK_ROWS; 
    dim_block.z = 1; 

    dim_grid.x = (width + TILE_DIM - 1)/TILE_DIM; 
    dim_grid.y = (height + TILE_DIM - 1)/TILE_DIM; 
    dim_grid.z = 1; 

    // Choose which GPU to run on, change this on a multi-GPU system. 
    cudaStatus = cudaSetDevice(0); 
    if (cudaStatus != cudaSuccess) 
    { 
     fprintf(stderr, "cudaSetDevice failed! Do you have a CUDA-capable GPU installed?"); 
     goto Error; 
    } 

    // Allocate GPU buffers for three matrix 
    cudaStatus = cudaMalloc((void**)&dev_matrA, size * sizeof(float)); 
    if (cudaStatus != cudaSuccess) 
    { 
     fprintf(stderr, "cudaMalloc failed!"); 
     goto Error; 
    } 

    cudaStatus = cudaMalloc((void**)&dev_matrB, size * sizeof(float)); 
    if (cudaStatus != cudaSuccess) 
    { 
     fprintf(stderr, "cudaMalloc failed!"); 
     goto Error; 
    } 

    // Copy input matrix from host memory to GPU buffers. 
    cudaStatus = cudaMemcpy(dev_matrA, matrA, size * sizeof(float), cudaMemcpyHostToDevice); 
    if (cudaStatus != cudaSuccess) 
    { 
     fprintf(stderr, "cudaMemcpy failed!"); 
     goto Error; 
    } 
    cudaMemset(dev_matrB, 0, size * sizeof(float)); 
    cudaEvent_t start, stop; 
    cudaEventCreate(&start); 
    cudaEventCreate(&stop); 

    switch (operation) 
    { 
     case(1): 
     { 
      cudaEventRecord(start); 
      // Launch a kernel on the GPU with one thread for each element. 
      naiveTransKernel << <dim_grid, dim_block >> >(dev_matrB, dev_matrA, width, height, nreps); 

      cudaEventRecord(stop); 
      cudaEventSynchronize(stop); 

      cudaEventElapsedTime(&elapsed, start, stop); 
      cudaEventDestroy(start); 
      cudaEventDestroy(stop); 

      elapsed /= nreps; 

      gpuBandwidth = (sizeof(float) * size * 2)/(elapsed * 1000000);//scaling from ms to s and B to GB doen implicitly, shortened in fraction, times two for read and write 
      printf("Avg. GPU Naive Trans Time: %f ms, bandwidth: %f GB/s\n", elapsed, gpuBandwidth); 

      break; 
     } 

     case(2): 
     { 
      cudaEventRecord(start); 
      // Launch a kernel on the GPU with one thread for each element. 
      notSoNaivaTransKernel << <dim_grid, dim_block >> >(dev_matrB, dev_matrA, width, height, nreps); 

      cudaEventRecord(stop); 
      cudaEventSynchronize(stop); 

      cudaEventElapsedTime(&elapsed, start, stop); 
      cudaEventDestroy(start); 
      cudaEventDestroy(stop); 

      elapsed /= nreps; 

      gpuBandwidth = (sizeof(float) * size * 2)/(elapsed * 1000000);//scaling from ms to s and B to GB doen implicitly, shortened in fraction, times two for read and write 
      printf("Avg. GPU Trans with SM Time: %f ms, bandwidth: %f GB/s\n", elapsed, gpuBandwidth); 

      break; 
     } 

    default: 
     printf("No matching opcode was found.\n"); 
    } 

    // Check for any errors launching the kernel 
    cudaStatus = cudaGetLastError(); 
    if (cudaStatus != cudaSuccess) 
    { 
     fprintf(stderr, "Kernel launch failed: %s\n", cudaGetErrorString(cudaStatus)); 
     goto Error; 
    } 

    // cudaDeviceSynchronize waits for the kernel to finish, and returns 
    // any errors encountered during the launch. 
    cudaStatus = cudaDeviceSynchronize(); 
    if (cudaStatus != cudaSuccess) 
    { 
     fprintf(stderr, "cudaDeviceSynchronize returned error code %d after launching Kernel!\n", cudaStatus); 
     goto Error; 
    } 

    // Copy output matrix from GPU buffer to host memory. 
    cudaStatus = cudaMemcpy(matrB, dev_matrB, size * sizeof(float), cudaMemcpyDeviceToHost); 
    if (cudaStatus != cudaSuccess) 
    { 
     fprintf(stderr, "cudaMemcpy failed!"); 
     goto Error; 
    } 

Error: 
    cudaFree(dev_matrB); 
    cudaFree(dev_matrA); 

    return cudaStatus; 
} 

void cpuMatrTrans(float *matrB, float *matrA, const int width, const int height, const int nreps) 
{ 
    int i, j, r; 

#pragma unroll 
    for (r = 0; r < nreps; r++) 
#pragma unroll 
     for (i = 0; i < height; i++) 
#pragma unroll 
      for (j = 0; j < width; j++) 
       matrB[j * height + i] = matrA[i * width + j]; 
} 

__global__ void naiveTransKernel(float *matrB, float *matrA, const int width, const int height, const int nreps) 
{ 
    int i, r; 
    int col = blockIdx.x * TILE_DIM + threadIdx.x; 
    int row = blockIdx.y * TILE_DIM + threadIdx.y; 
    int index_in = col + width * row; 
    int index_out = row + height * col; 

#pragma unroll 
    for (r = 0; r < nreps; r++) 
#pragma unroll 
     for (i = 0; i < TILE_DIM; i += BLOCK_ROWS) 
      if ((row+i<height) && (col < width)) 
       matrB[index_out + i] = matrA[index_in + i * width]; 
} 

__global__ void notSoNaivaTransKernel(float *matrB, float *matrA, const int width, const int height, const int nreps) 
{ 
    __shared__ float tile[TILE_DIM][TILE_DIM + 1]; 
    int ciIndex = blockIdx.x * TILE_DIM + threadIdx.x; 
    int riIndex = blockIdx.y * TILE_DIM + threadIdx.y; 
    int coIndex = blockIdx.y * TILE_DIM + threadIdx.x; 
    int roIndex = blockIdx.x * TILE_DIM + threadIdx.y; 
    int index_in = ciIndex + (riIndex)* width; 
    int index_out = coIndex + (roIndex)* height; 


    int r, i; 
#pragma unroll 
    for (r = 0; r < nreps; r++) 
    { 
#pragma unroll 
     for (i = 0; i < TILE_DIM; i += BLOCK_ROWS) 
      if ((ciIndex<width) && (riIndex+i < height)) 
       tile[threadIdx.y + i][threadIdx.x] = matrA[index_in + i * width]; 
     __syncthreads(); 

#pragma unroll 
     for (i = 0; i < TILE_DIM; i += BLOCK_ROWS) 
      if ((coIndex<height) && (roIndex+i < width)) 
       matrB[index_out + i*height] = tile[threadIdx.x][threadIdx.y + i]; 
     __syncthreads(); 
    } 
} 
$ nvcc -std=c++11 -arch=sm_61 -o t33 t33.cu 
t33.cu(25): warning: variable "matrAC" was set but never used 

t33.cu(25): warning: variable "matrAC" was set but never used 

$ cuda-memcheck ./t33 
========= CUDA-MEMCHECK 
Avg. CPU Transpose Time: 0.143087 ms, Bandwidth: 2.795509 GB/s 

Avg. GPU Naive Trans Time: 0.028587 ms, bandwidth: 13.992195 GB/s 
     Correct: 50000, Wrong: 0 

Avg. GPU Trans with SM Time: 0.040328 ms, bandwidth: 9.918678 GB/s 
     Correct: 50000, Wrong: 0 

========= ERROR SUMMARY: 0 errors 
$ ./t33 
Avg. CPU Transpose Time: 0.140469 ms, Bandwidth: 2.847594 GB/s 

Avg. GPU Naive Trans Time: 0.003828 ms, bandwidth: 104.505440 GB/s 
     Correct: 50000, Wrong: 0 

Avg. GPU Trans with SM Time: 0.000715 ms, bandwidth: 559.206604 GB/s 
     Correct: 50000, Wrong: 0 

$ 

Note: Le code tente de mesurer la bande passante. Cependant, vous devez savoir que la bande passante mesurée est affectée par la bande passante du cache. Vos tailles de matrice (500 x 100 = 200 Ko chacune pour l'entrée et la sortie) sont ici assez petites pour tenir dans le cache L2 de la plupart des GPU. Ce fait, couplé avec le fait que vous exécutez la même transposition plusieurs fois (nreps) signifie qu'une grande partie du travail fonctionne directement à partir du cache L2. Par conséquent, dans le cas "optimisé" ci-dessus, nous voyons un nombre de bande passante signalé qui dépasse largement la bande passante mémoire disponible du GPU (il s'agit d'un Pascal Titan X, donc ~ 340Go/s de bande passante mémoire principale disponible). Cela est dû au fait que cette mesure inclut certains avantages du cache L2, dont la bande passante est au moins deux fois plus élevée que la bande passante de la mémoire principale. Vous pouvez éliminer cet effet en utilisant une taille de matrice significativement plus grande et/ou en réduisant nreps à 1.