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を実行したら、しかし：。
= [Matrix C] [Matrix A] X [Matrix B]（代わりに、この記事にいくつかの巨大な行列を貼り付けるしかし、あなたはそれを通じてその半分の方法を書くために開始さ見ることができますcのマトリックスでそれらにリンクされています私のグラフィックスカードには32x32マトリックスの1024スレッドの制限があります。また、100x100マトリックスを実行すると、マトリックスCはすべて0です。

mem_size_Xをsizeof（float）* size_X、size_Xをブロックのサイズは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(); 

カーネルコード：私はあなたのearlier, almost identical questionにあなたに言ったよう

// 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などの場合にしか使用できません。have done with dynamically allocated shared memoryはこれを変更しません。

これを確認するには、カーネルを実行し、実行したかどうかをチェックし、その出力をホスト上で行われた単純な参照計算と比較する完全なコンパイル可能な再現ケースから始めましょう。このコードは、公開されたカーネルと、起動パラメータの計算のコアです。 stdinからサイズを読み込み、ケースを実行します。結果が特定の許容値を超えて異なる場合は、アサートエラーが発生します。ここでは、コードがある、それはCUDA 3.0でコンパイル以降および任意のCUDA互換GPU上で実行する必要があります。

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

だから今、のサイズが異なるため、このコードを実行してみましょう。 GPU上のコードが間違っていないことを確認するために、境界外のメモリアクセスを検出できるcuda-memcheckユーティリティを使用して実行します。次のテストはblocksize=16を使用して、計算能力を1.2カードとCUDA 3.2をOS X 10.6のマシン上で行われたのすべて：

$ 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] 

さんが行列blocksize最初

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

未満の場合を試してみましょう

ここでは、無効な構成引数エラーでカーネルを実行できませんでした。どうして？このため：0グリッドサイズm,n < blocksizeになり

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

。

次はのは、このケース16には、ブロックサイズの最小ラウンドの複数を試してみましょう：エラーなしで実行される、または失敗を主張

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

を。

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 

、我々は境界メモリアクセス検出と期待されている打ち上げ失敗エラーの出る：今度は17にサイズを大きくしてみましょう。

$ 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 

境界エラーのうちが発生しているなぜあなたは従わない場合でも、あなたは、少なくともある今64、96、および128を試すことができますこのコードは本当にブロックサイズの倍数である行列に対してのみ正しく動作するということを受け入れることを望んでいますか？