CUDA Matrix Multiplication записывает в неправильную область памяти

Question

Идея моей простой программы, которую я пытался написать, заключается в том, чтобы получить информацию от пользователя, чтобы увидеть, как много матрицы умножается.

dd@cuda-Linux:~/Desktop/multi$ ./program
What is the rowSize of a? 33
What is the colSize of a? 33
What is the rowSize of b? 33
What is the colSize of b? 33
Would you like to write the results to a file?(y or n)
y
Creating the random numbers now
Writing Matrix A to file now...
Writing Matrix B to file now...
Starting it on the device
Writing Matrix C to file now...
Finish

Однако проблема заключается в моих вычислениях потока.Я могу перейти к матрице 32x32, и она будет работать нормально и даст мне правильные результаты.Однако, когда я запускаю 33x33, я получаю результаты, подобные следующим:
[Матрица A] x [Матрица B] = [Матрица C] (связанныйим вместо того, чтобы вставлять несколько огромных матриц в этот пост. Но с матрицей c вы можете видеть, что на полпути она начинает писать неправильные числа. Моя видеокарта имеет ограничение в 1024 потока, что является матрицей 32x32. Также, когда я идудля запуска матрицы 100x100 Матрица C - все 0s.

Пусть mem_size_X будет sizeof (float) * size_X, а size_X - это высота * ширина матрицы. Сейчас высота и ширина должны быть одинаковыми, поэтому 32x32Также "block_size" - это просто высота. Таким образом, для матрицы 32x32 размер блока соответствует 32.
Код хоста (запуск):

    float* deviceMatrixA;
    float* deviceMatrixB;
    cudaMalloc((void**) &deviceMatrixA, mem_size_A);//allocate mem_size_x on the device. 
    cudaMalloc((void**) &deviceMatrixB, mem_size_B);


    cudaMemcpy(deviceMatrixA, a.elements, mem_size_A, cudaMemcpyHostToDevice);
    cudaMemcpy(deviceMatrixB, b.elements, mem_size_B, cudaMemcpyHostToDevice);



    int size_C = c.rowSize * c.colSize;
    int mem_size_C = sizeof(float) * size_C;
    c.elements = (float*) malloc(mem_size_C);


    float* deviceMatrixC;
    cudaMalloc((void**) &deviceMatrixC, mem_size_C);


    dim3 threads(block_size, block_size);
    dim3 grid(c.colSize / threads.x, c.rowSize / threads.y);



    matrixMul<<< grid, threads,2*block_size*block_size*sizeof(float)>>>(deviceMatrixC, deviceMatrixA, deviceMatrixB, a.colSize, b.colSize, block_size);//sizeof(float)*block_size*block_size
    cudaThreadSynchronize();

Код ядра:

// CUDA Kernel
__global__ void matrixMul( float* C, float* A, float* B, int wA, int wB,size_t block_size)
{
    int bx = blockIdx.x;
    int by = blockIdx.y;
    int tx = threadIdx.x;
    int ty = threadIdx.y;

    int aBegin = wA * block_size * by;
    int aEnd   = aBegin + wA - 1;
    int aStep  = block_size;

    int bBegin = block_size * bx;

    int bStep  = block_size * wB;
    float Csub=0;

    for (int a = aBegin, b = bBegin; a <= aEnd; a += aStep, b += bStep) 
    {
        extern __shared__ float As[];
        extern __shared__ float Bs[];
        extern __shared__ float smem[];

        smem[ty*block_size+tx] = A[a + wA * ty + tx];

        smem[block_size*block_size+ty*block_size+tx]  = B[b + wB * ty + tx];

        __syncthreads();

        for (int k = 0; k < block_size; ++k)
            Csub += smem[ty*block_size+k] * smem[block_size*block_size+k*block_size+tx] ;

        __syncthreads();
    }

    int c = wB * block_size * by + block_size * bx;
    C[c + wB * ty + tx] = Csub;


}

Спасибо

Answer 1

Как я уже говорил вам на вашем ранее почти идентичном вопросе , этот код умножения матриц предназначен только для вычислений на матрицах, размеры которых кратны block_size.Если вы выберете block_size = 32, то его можно использовать только для 32x32, 64x64, 96x96, 128x128 и т. Д. Ничто из того, что вы сделали с динамически распределенной общей памятью , не изменит этого.

Чтобы убедиться, что это так, давайте начнем с полного, скомпилируемого репро-случая, который запустит ваше ядро, проверит, выполнено ли оно, и сравнит его вывод с простыми ссылочными вычислениями, выполненными на хосте.Этот код - ваше опубликованное ядро, а также ядро ​​вычислений параметров запуска.Он будет читать размер из стандартного ввода, а затем запустить дело.Если результаты отличаются более чем на определенное отклонение, возникнет ошибка подтверждения.Вот код, который должен скомпилироваться на CUDA 3.0 или более поздней версии и работать на любом CUDA-совместимом графическом процессоре:

#include <assert.h>
#include <cstdio>
#include <cstdlib>
#include <cmath>

inline void GPUassert(cudaError_t code, char * file, int line, bool Abort=true)
{
    if (code != 0) {
        fprintf(stderr, "GPUassert: %s %s %d\n", cudaGetErrorString(code),file,line);
        if (Abort) exit(code);
    }       
}

#define GPUerrchk(ans) { GPUassert((ans), __FILE__, __LINE__); }

__global__ void matrixMul( float* C, float* A, float* B, int wA, int wB, size_t block_size)
{
    int bx = blockIdx.x;
    int by = blockIdx.y;
    int tx = threadIdx.x;
    int ty = threadIdx.y;

    int aBegin = wA * block_size * by;
    int aEnd   = aBegin + wA - 1;
    int aStep  = block_size;
    int bBegin = block_size * bx;
    int bStep  = block_size * wB;

    float Csub=0.f;
    for (int a = aBegin, b = bBegin; a <= aEnd; a += aStep, b += bStep) 
    {
        extern __shared__ float smem[];

        smem[ty*block_size+tx] = A[a + wA * ty + tx];
        smem[block_size*block_size+ty*block_size+tx]  = B[b + wB * ty + tx];

        __syncthreads();

        for (int k = 0; k < block_size; ++k)
            Csub += smem[ty*block_size+k] * smem[block_size*block_size+k*block_size+tx] ;

        __syncthreads();
    }

    int c = wB * block_size * by + block_size * bx;
    C[c + wB * ty + tx] = Csub;
}

inline float frand(){
    return (float)rand()/(float)RAND_MAX;
}

void matmul(float *C, const float *A, const float *B,  int wA, int wB)
{
    for(int k=0; k<wB; k++) {
        for(int j=0; j<wB; j++) {
            float dotp = 0.f;
            for(int i=0; i<wA; i++) {
                dotp += A[j*wA+i] * B[i*wB+k];
            }
            C[j*wB+k] = dotp;
        }
    }
}

int main(int argc, char ** argv)
{
    int val = 128;

    if ( argc == 2 ) {
        val = atoi(argv[1]);
    }

    int m = val, n = val, mn = m*n;
    size_t sz = size_t(mn) * sizeof(float);

    srand(time(NULL));

    float * A = new float[mn], * B = new float[mn], * C= new float[mn];
    float * A_, * B_, * C_;

    for(int i=0; i<mn; i++) {
        A[i] = frand(); B[i] = frand();
    }

    GPUerrchk( cudaMalloc((void **)&A_, sz) );
    GPUerrchk( cudaMalloc((void **)&B_, sz) );
    GPUerrchk( cudaMalloc((void **)&C_, sz) );

    GPUerrchk( cudaMemcpy(A_, A, sz, cudaMemcpyHostToDevice) );
    GPUerrchk( cudaMemcpy(B_, B, sz, cudaMemcpyHostToDevice) );

    // Launch configuration
    // Note that the input matrice sizes *must* be a round
    // multiple of blocksize for this code to work correctly.
    const int blocksize=16;
    const int shmsz = size_t(2*blocksize*blocksize) * sizeof(float);
    dim3 block=dim3(blocksize,blocksize), grid = dim3(m/block.x,m/block.y); 

    matrixMul<<<grid,block,shmsz>>>(C_,A_,B_,m,n,blocksize);
    GPUerrchk( cudaPeekAtLastError() );

    GPUerrchk( cudaMemcpy(C, C_, sz, cudaMemcpyDeviceToHost) );

    // Verfication on host
    float * Cref = new float[mn];
    matmul(Cref,A,B,m,n); 
    const float tol = 5e-5f;
    for(int i=0; i<mn; i++) {
        assert(fabs(C[i]-Cref[i])/C[i] < tol);
    }

    GPUerrchk( cudaThreadExit() ); // CUDA 3.2 compatible

    return 0;
}

Итак, теперь давайте запустим этот код для разных размеров.Чтобы убедиться, что код на графическом процессоре не сделал ничего плохого, я запустим его с помощью утилиты cuda-memcheck, которая может обнаруживать доступ к памяти вне пределов.Все следующие тесты были выполнены на компьютере с OS X 10.6 с картой Compute 1.2 и CUDA 3.2 с использованием blocksize=16:

$ nvcc -arch=sm_12 -Xcompiler="-Wall" -Xptxas="-v" -o matmul2 matmul2.cu
ptxas info    : Compiling entry function '_Z9matrixMulPfS_S_iim' for 'sm_12'
ptxas info    : Used 16 registers, 32+16 bytes smem, 4 bytes cmem[1]

Давайте рассмотрим случай, когда матрицы меньше blocksizefirst

$ cuda-memcheck ./matmul2 4
========= CUDA-MEMCHECK
GPUassert: invalid configuration argument matmul2.cu 101
========= ERROR SUMMARY: 0 errors

Здесь мы не смогли запустить ядро ​​с ошибкой неверного аргумента конфигурации.Зачем?Из-за этого:

    dim3 block=dim3(blocksize,blocksize), grid = dim3(m/block.x,m/block.y); 

, что приводит к размеру сетки 0, когда m,n < blocksize.

Далее давайте попробуем наименьшее круглое кратное размера блока, в данном случае 16:

$ cuda-memcheck ./matmul2 16
========= CUDA-MEMCHECK
========= ERROR SUMMARY: 0 errors

, который работает без ошибок или подтверждает ошибку.Давайте теперь увеличим размер до 17:

cuda-memcheck ./matmul2 17
========= CUDA-MEMCHECK
GPUassert: unspecified launch failure matmul2.cu 103
========= Invalid __global__ read of size 4
=========     at 0x000001f8 in matrixMul
=========     by thread (0,2,0) in block (0,0)
=========     Address 0x001009c8 is out of bounds
=========
========= ERROR SUMMARY: 1 error

, и мы выйдем за границы обнаруженного доступа к памяти и ошибки при запуске, что ожидается.Теперь давайте попробуем 64, 96 и 128:

$ cuda-memcheck ./matmul2 64
========= CUDA-MEMCHECK
========= ERROR SUMMARY: 0 errors

$ cuda-memcheck ./matmul2 96
========= CUDA-MEMCHECK
========= ERROR SUMMARY: 0 errors

$ cuda-memcheck ./matmul2 128
========= CUDA-MEMCHECK
========= ERROR SUMMARY: 0 errors

и, наконец, попробуем 129:

$ cuda-memcheck ./matmul2 129
========= CUDA-MEMCHECK
GPUassert: unspecified launch failure matmul2.cu 103
========= Invalid __global__ read of size 4
=========     at 0x000001f8 in matrixMul
=========     by thread (0,1,0) in block (0,0)
=========     Address 0x00120904 is out of bounds
=========
========= ERROR SUMMARY: 1 error

Даже если вы не понимаете, почему возникают ошибки вне границ,Вы хотя бы согласны с тем, что этот код действительно работает правильно только для матриц, кратных размеру блока?