first.cu 关于CUDA的第一个程序，备注理解架构

// includes, system
#include <stdio.h>
#include <stdlib.h>

#define RADIUS 1
#define BLOCKDIM 16
#define N 13

__global__ static void set_global_idx(int n, int *d_a)
{
    int idx = blockIdx.x*blockDim.x + threadIdx.x;
    if (idx < n)
        d_a[idx] = idx;
   
}

__global__ static void set_block_idx(int n, int *d_a)
{
    int idx = blockIdx.x*blockDim.x + threadIdx.x;
    if (idx < n)
        d_a[idx] = blockIdx.x;
   
}

__global__ static void set_thread_idx(int n, int *d_a)
{
    int idx = blockIdx.x*blockDim.x + threadIdx.x;
    if (idx < n)
        d_a[idx] = threadIdx.x;
   
}

__global__ static void stencil_naive(int n, int *in, int *out)
{
    int idx = blockIdx.x*blockDim.x + threadIdx.x;
    int i;
    int value = 0;
    if (idx < n) {
        for (i=-RADIUS; i<=RADIUS; ++i)
        {
            value += in[idx+i+RADIUS];
        }
    }
    out[idx] = value;
    
}

__global__ static void stencil(int n, int *in, int *out)
{
    __shared__ int shared[BLOCKDIM + 2*RADIUS];
    int idx = blockIdx.x*blockDim.x + threadIdx.x;
    int idx_local = threadIdx.x + RADIUS;
    shared[idx_local] = in[idx + RADIUS];
    int i;
    if (idx < n) {
        if( threadIdx.x < RADIUS)
        {
            shared[idx_local - RADIUS] = in[idx ];
            shared[idx_local + BLOCKDIM] = in[idx + BLOCKDIM + RADIUS];
        }
    }
    __syncthreads();
    int value = 0;
    for (i=-RADIUS; i<=RADIUS; i++)
    {
        value += shared[idx_local + i];
    }
    out[idx] = value;
    
}

__global__ static void matrix_multiply_naive(int n, int *in_a, int *in_b, int *out)
{
    int tmp = gridDim.x*blockDim.x;
    int row = threadIdx.x + blockIdx.x*blockDim.x;
    int col = threadIdx.y + blockIdx.y*blockDim.y;
    int idx = row*tmp + col;
    int sum = 0;
    if (idx < n)
    {
        for (int i=0; i<N; i++)
        {
            sum += in_a[row*N+i]*in_b[i*N+col];
        }
    }
    out[row*N+col] = sum;

    /*int row = threadIdx.x + blockIdx.x*blockDim.x;
    int col = threadIdx.y + blockIdx.y*blockDim.y;
    int idx = row*gridDim.x*blockDim.x + col;
    if (idx < n)
    {
        out[idx] = idx;
    }*/
    


}


// Start the main SDK sample here
int main(int argc, char** argv)
{
    printf(" ###########################\n");
    printf(" #test for cuda toolkit ...#\n");
    printf(" ###########################\n");

    int dim = 16;
    int mem_size = dim*sizeof(int);

    int *d_a, *h_a, *d_b;
    h_a = (int*)malloc(mem_size);
    cudaMalloc((void**)&d_a, mem_size);
    cudaMalloc((void**)&d_b, mem_size);

    if ( 0==h_a || 0==d_a)
    {
        printf("could not allocate memory\n");
        return 1;
    }
    cudaMemset(d_a, 0, mem_size);
    cudaMemcpy(h_a, d_a, mem_size, cudaMemcpyDeviceToHost);

    for (int i=0; i<dim; ++i)
    {
        printf("%2d ", h_a[i]);
    }
    printf("\n");

    int dim_block = 3;
    int dim_grid = (dim + dim_block -1)/dim_block;
    set_global_idx<<<dim_grid, dim_block>>> (dim, d_a);
    cudaMemcpy(h_a, d_a, mem_size, cudaMemcpyDeviceToHost);

    for (int i=0; i<dim; ++i)
    {
        printf("%2d ", h_a[i]);
    }
    printf("\n");

    set_block_idx<<<dim_grid, dim_block>>> (dim, d_a);
    cudaMemcpy(h_a, d_a, mem_size, cudaMemcpyDeviceToHost);

    for (int i=0; i<dim; ++i)
    {
        printf("%2d ", h_a[i]);
    }
    printf("\n");

    set_thread_idx<<<dim_grid, dim_block>>> (dim, d_a);
    cudaMemcpy(d_b, d_a, mem_size, cudaMemcpyDeviceToDevice);
    cudaMemcpy(h_a, d_b, mem_size, cudaMemcpyDeviceToHost);

    for (int i=0; i<dim; ++i)
    {
        printf("%2d ", h_a[i]);
    }
    printf("\n");

    free(h_a);
    cudaFree(d_a);
    cudaFree(d_b);

    ///
    printf("\n");
    printf(" #######################\n");
    printf(" #stencil test ......  #\n");
    printf(" #######################\n");
    dim = 80;
    int dim_pro = dim + RADIUS*2;
    int mem_size_pro = sizeof(int)*dim_pro;
    mem_size = sizeof(int)*dim;
    h_a = (int*)malloc(mem_size_pro);
    memset(h_a, 0, mem_size_pro);
    int *h_b;
    cudaMalloc((void**)&d_a, mem_size_pro);
    cudaMalloc((void**)&d_b, mem_size);
    for (int i=0; i<dim; i++)
    {
        h_a[i+RADIUS] = i;
    }
    h_b = &(h_a[RADIUS]);
    printf(" original vector  ... \n");
    for (int i=0; i<dim; ++i)
    {
        if (i%10 == 0)
            printf("\n");
        printf("%5d ", h_b[i]);
    }
    printf("\n\n\n");
    cudaMemcpy(d_a, h_a, mem_size_pro, cudaMemcpyHostToDevice);
    dim_grid = (dim + BLOCKDIM -1)/BLOCKDIM;
    //cudaMemset(d_b, 0, mem_size);
    stencil<<<dim_grid, BLOCKDIM>>> (dim, d_a, d_b);
    //stencil_naive<<<dim_grid, BLOCKDIM>>> (dim, d_a, d_b);
    cudaMemcpy(h_b, d_b, mem_size, cudaMemcpyDeviceToHost);

    
    printf(" after stencil   ... \n");
    for (int i=0; i<dim; ++i)
    {
        if (i%10 == 0)
            printf("\n");
        printf("%5d ", h_b[i]);
    }
    printf("\n");

    free(h_a);
    cudaFree(d_a);
    cudaFree(d_b);

    ///
    printf("\n");
    printf(" ##########################\n");
    printf(" #matrix multiply ......  #\n");
    printf(" ##########################\n");

    int row = N;
    int col = N;
    dim = row*col;
    mem_size = dim*sizeof(int);
    h_a = (int*)malloc(mem_size);
    int *d_c;
    cudaMalloc((void**)&d_a, mem_size);
    cudaMalloc((void**)&d_b, mem_size);
    cudaMalloc((void**)&d_c, mem_size);
    printf(" matrix A & B : \n");
    for (int i=0; i<dim; ++i)
    {
        h_a[i] = 1;
        if (i%row == 0)
        {
            printf("\n");
        }
        printf("%5d ", h_a[i]);
    }
    cudaMemcpy(d_a, h_a, mem_size, cudaMemcpyHostToDevice);
    cudaMemcpy(d_b, h_a, mem_size, cudaMemcpyHostToDevice);
    dim_block = 2;
    dim3 block(dim_block, dim_block);
    dim3 grid((row+dim_block-1)/dim_block, (col+dim_block-1)/dim_block);
    cudaMemset(d_c, 0, mem_size);
    matrix_multiply_naive<<<grid, block>>> (dim, d_a, d_b, d_c);

    //memset(h_a, 0,mem_size);
    cudaMemcpy(h_a, d_c, mem_size, cudaMemcpyDeviceToHost);
    printf("\n matrix C : \n");
    for (int i=0; i<dim; ++i)
    {
        if (i%row == 0)
        {
            printf("\n");
        }
        printf("%5d ", h_a[i]);
    }
    
    
    return 0;
}

——————————————————————————————————————————————————————————————————

#include <stdio.h>

#include <stdlib.h>
#define RADIUS 3
#define BLOCK_SIZE 16


__global__ void stencil(int* in, int* out)

{
  __shared__ int shared[BLOCK_SIZE + 2 * RADIUS];
  int globIdx = blockIdx.x * blockDim.x + threadIdx.x;
  int locIdx  = threadIdx.x + RADIUS;
  shared[locIdx] = in[globIdx];
 
  if (threadIdx.x < RADIUS)
  {
    shared[locIdx – RADIUS] = in[globIdx – RADIUS];
    shared[locIdx + BLOCK_DIMX] = in[globIdx + BLOCK_SIZE];
  }
  __syncthreads();
 
  int value = 0;
  for (offset = - RADIUS; offset <= RADIUS; offset++)
    value += shared[locIdx + offset];
  out[globIdx] = value;

}




void checkCUDAError(const char *msg)
{
    cudaError_t err = cudaGetLastError();
    if( cudaSuccess != err)
    {
        fprintf(stderr, "Cuda error: %s: %s.\n", msg, cudaGetErrorString( err) );
        exit(-1);
    }                         
}

int main( int argc, char** argv)
{
    int *h_a;
    int *d_a;
    int numBlocks = 16;
    int numThreadsPerBlock = 16;

    size_t memSize = numBlocks * numThreadsPerBlock * sizeof(int);
    h_a = (int *) malloc(memSize);
    cudaMalloc((void**)&d_a, memSize );

    if( 0==h_a || 0==d_a )
    {
        printf("couldn't allocate memory\n");
        return 1;
    }
    
    dim3 dimGrid(16,16);
    dim3 dimBlock(16);
    stencil<<< dimGrid,dimBlock>>>(RADIUS);

   cudaThreadSynchronize();

   checkCUDAError("kernel execution");

   cudaMemcpy(h_a, d_a, memSize, cudaMemcpyDeviceToHost);

   checkCUDAError("cudaMemcpy");

   for (int i = 0; i < 16 ; i++)
    {
        for (int j = 0; j < 16 ; j++)
        {
            assert(h_a[i * numThreadsPerBlock + j] == 1000 * i + j);
        }
    }


    cudaFree(d_a);
    free(h_a);
    printf("Correct!\n");
    return 0;
}