0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
/**
 * @brief  implemenntation of the GEMM cuda in global memory
 * 
 * @param a  the first matrix in m*n in device memory
 * @param b  the second matrix in n*k in device memory
 * @param c  the result matrix in m*k in device memory
 * @param m  matrix size
 * @param n  matrix size
 * @param k  matrix size 
 * @return the result in ptr c
 */
__global__ void gpu_matrix_mult_gm(float *a, float *b, float *c, int m, int n, int k)
{
    int row = blockIdx.y * blockDim.y + threadIdx.y
    int col = blockIdx.x * blockDim.x + threadIdx.x
    float temp = 0
    if (row < m && col < k) // Ensure bounds are within the matrix dimensions
    {
        for (int i = 0 i < n i++)
        {
            temp += a[row * n + i] * b[i * k + col]
        }
        c[row * k + col] = temp
    }
}

0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
/**
 * @brief  implemenntation of the GEMM cuda in shared memory
 *
 * @param a  the first matrix in m*n in device memory
 * @param b  the second matrix in n*k in device memory
 * @param c  the result matrix in m*k in device memory
 * @param m  matrix size
 * @param n  matrix size
 * @param k  matrix size
 * @return the result in ptr c
 */
__global__ void gpu_matrix_mult_sm(float *a, float *b, float *c, int m, int n, int k)
{
    __shared__ float tile_a[BLOCK_SIZE][BLOCK_SIZE]
    __shared__ float tile_b[BLOCK_SIZE][BLOCK_SIZE]
    int row = blockIdx.y * BLOCK_SIZE + threadIdx.y
    int col = blockIdx.x * BLOCK_SIZE + threadIdx.x
    float temp = 0

    int block_num = (n + BLOCK_SIZE - 1) / BLOCK_SIZE // Divide the n dimension by block size
    for (int i = 0 i < block_num ++i)
    {
        // Copy the ith block into shared memory
        if (row < m && i * BLOCK_SIZE + threadIdx.x < n)
        {
            int a_index = row * n + i * BLOCK_SIZE + threadIdx.x
            tile_a[threadIdx.y][threadIdx.x] = a[a_index]
        }
        else
            tile_a[threadIdx.y][threadIdx.x] = 0 // Handle edge case for smaller matrix

        if (col < k && i * BLOCK_SIZE + threadIdx.y < n)
        {
            int b_index = (i * BLOCK_SIZE + threadIdx.y) * k + col
            tile_b[threadIdx.x][threadIdx.y] = b[b_index]
        }
        else
            tile_b[threadIdx.y][threadIdx.x] = 0

        __syncthreads() // Synchronize threads in the block before computation

        // Compute the contribution for c[row][col]
        for (int j = 0 j < BLOCK_SIZE j++)
        {
            temp += tile_a[threadIdx.y][j] * tile_b[j][threadIdx.x]
        }
        __syncthreads() // Ensure all threads finish their computation before moving on
    }
    if (row < m && col < k)
    {
        c[row * k + col] = temp
    }
}

0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
        // Copy the ith block into shared memory
        if (row < m && i * BLOCK_SIZE + threadIdx.x < n)
        {
            int a_index = row * n + i * BLOCK_SIZE + threadIdx.x
            tile_a[threadIdx.y][threadIdx.x] = a[a_index]
        }
        else
            tile_a[threadIdx.y][threadIdx.x] = 0 // Handle edge case for smaller matrix

        if (col < k && i * BLOCK_SIZE + threadIdx.y < n)
        {
            int b_index = (i * BLOCK_SIZE + threadIdx.y) * k + col
            tile_b[threadIdx.x][threadIdx.y] = b[b_index]
        }
        else
            tile_b[threadIdx.x][threadIdx.y] = 0

        __syncthreads() // Synchronize threads in the block before computation

        // Compute the contribution for c[row][col]
        for (int j = 0 j < BLOCK_SIZE j++)
        {
            temp += tile_a[threadIdx.y][j] * tile_b[threadIdx.x][j]
        }

0
1
2
3
4
5
6
//   m*k k*n = m*n
	cublasH     cublasHandle_t handle
    cublasCreate(&handle)
    // Configure cuBLAS operations
    cublasOperation_t opA = CUBLAS_OP_N // No transpose for A
    cublasOperation_t opB = CUBLAS_OP_N // No transpose for B
	cublasSgemm(handle, opA, opB, n, m, k, &alpha, d_B, n, d_A, k, &beta, d_C, n) // 我日，cublas 居然是列优先的

0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
// m n, depth分别是矩阵的行，列，通道数
float *padMatrix(const float *matrix, float *padded_matrix, int m, int n, int depth, int padding)
{
    int rows = m
    int cols = n
    int paddedCols = n + 2 * padding
    int paddedRows = m + 2 * padding
    int paddedBlcok = paddedRows * paddedCols
    int originblock = rows * cols
    for (int k = 0 k < depth k++)
    {
        for (int i = 0 i < rows ++i)
        {
            for (int j = 0 j < cols ++j)
            {
                padded_matrix[k * paddedBlcok + (i + padding) * paddedCols + j + padding] = matrix[k * originblock + i * n + j]
            }
        }
    }
    return padded_matrix
}   
   memset(temp, 0, sizeof(float) * paddedRows * paddedCols * input.matrix_nums)  //需要把原来的内存设置为0，否则padding的时候可能会出错
    padMatrix(h_input, temp, m, n, input.matrix_nums, padding)

0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
//输入m,n是原本的行列数
__global__ void conv2d_cal(float *input, float *output, int output_block, int m, int n, const Kernel kernelSet, int stride, int padding)
{
    int row = blockIdx.y * blockDim.y + threadIdx.y
    int col = blockIdx.x * blockDim.x + threadIdx.x
    int paddedm = m + 2 * padding
    int paddedn = n + 2 * padding

    int conved_n = (n - kernelSet.cols + 2 * padding) / stride + 1  //计算卷积后的列数
    
    // extern __shared__ float temp[] // 使用动态共享内存，适配 kernelSet.numKernels 注意会出现竞争！！！
    float temp

    if (row % stride == 0 && col % stride == 0 && row < paddedm - kernelSet.rows + 1 && col < paddedn - kernelSet.cols + 1)
    {
        for (int k = 0 k < kernelSet.numKernels k++)
        {
            for (int i = 0 i < kernelSet.rows i++)
            {
                for (int j = 0 j < kernelSet.cols j++)
                {
                    int input_row = row + i
                    int input_col = col + j
                    int kernel_idx = k * kernelSet.rows * kernelSet.cols + i * kernelSet.cols + j
                    temp += input[input_row * paddedn + input_col] * kernelSet.deviceKernels[kernel_idx]
                }
            }
            // 写入输出           
        }
        output[row * conved_n + col] = temp
    }
}

0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
row: 0, col: 2, i: 0, j:0, input:0.000000, kernel:1.000000, temp: 0.000000
row: 0, col: 2, i: 0, j:1, input:0.000000, kernel:2.000000, temp: 0.000000
row: 0, col: 2, i: 0, j:2, input:0.000000, kernel:3.000000, temp: 0.000000
row: 0, col: 2, i: 1, j:0, input:2.000000, kernel:4.000000, temp: 8.000000
row: 0, col: 2, i: 1, j:1, input:3.000000, kernel:5.000000, temp: 23.000000
row: 0, col: 2, i: 1, j:2, input:4.000000, kernel:6.000000, temp: 47.000000
row: 0, col: 2, i: 2, j:0, input:7.000000, kernel:7.000000, temp: 96.000000
row: 0, col: 2, i: 2, j:1, input:8.000000, kernel:8.000000, temp: 160.000000
row: 0, col: 2, i: 2, j:2, input:9.000000, kernel:9.000000, temp: 241.000000
Custom Kernel Time for sliding conv: 4.198948 ms
Printing matrix of printing 0th matrix
384.000000 606.000000 723.000000 570.000000 312.000000 
318.000000 513.000000 648.000000 603.000000 354.000000 
483.000000 738.000000 873.000000 648.000000 339.000000 
318.000000 513.000000 648.000000 603.000000 354.000000 
150.000000 228.000000 291.000000 264.000000 150.000000

0
1
2
3
4
5
6
7
8
9
10
11
12
13
row: 0, col: 3, i: 0, j:0, input:0.000000, kernel:1.000000, temp: 0.000000
row: 0, col: 3, i: 0, j:1, input:0.000000, kernel:2.000000, temp: 0.000000
row: 0, col: 3, i: 0, j:2, input:0.000000, kernel:3.000000, temp: 0.000000
row: 0, col: 3, i: 1, j:0, input:3.000000, kernel:4.000000, temp: 12.000000
row: 0, col: 3, i: 1, j:1, input:4.000000, kernel:5.000000, temp: 32.000000
row: 0, col: 3, i: 1, j:2, input:5.000000, kernel:6.000000, temp: 62.000000
row: 0, col: 3, i: 2, j:0, input:8.000000, kernel:7.000000, temp: 118.000000
row: 0, col: 3, i: 2, j:1, input:9.000000, kernel:8.000000, temp: 190.000000
row: 0, col: 3, i: 2, j:2, input:0.000000, kernel:9.000000, temp: 190.000000
Custom Kernel Time for sliding conv: 4.153343 ms
Printing matrix of printing 0th matrix
384.000000 570.000000 
318.000000 603.000000 

0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
    int padded_m = h_input.rows + 2 * padding  //padding
    int padded_n = h_input.cols + 2 * padding
    int convable_m = (h_input.rows + 2 * padding - kernelSet.rows) / stride + 1  //after convolution
    int convable_n = (h_input.cols + 2 * padding - kernelSet.cols) / stride + 1

    int converted_n = convable_m * convable_n  // converted to a flat matrix
    int converted_m = kernelSet.rows * kernelSet.cols

    int block = converted_m * converted_n
    int padded_block = padded_m * padded_n

    float *temp
    // printf("hello\n")
    cudaMallocHost((void **)&temp, sizeof(float) * padded_block * h_input.matrix_nums)  //create a temp to get the padded matrix in host mem

    padMatrix(h_input.h_input, temp, h_input.rows, h_input.cols, h_input.matrix_nums, padding)


    float *h_temp
    cudaMallocHost((void **)&h_temp, sizeof(float) * h_input.matrix_nums * block)  //for converted matrix
    for (int i = 0 i +kernelSet.rows-1 < padded_m i+=stride) // the next 2 iterations are for sliding the conv window
    {
        for (int j = 0  j+kernelSet.cols-1<padded_n j+=stride)
        {
            for (int a = 0 a < kernelSet.rows a++) // the naxt 2 iterations are for inside the conv window
            {
                for (int b = 0 b < kernelSet.cols b++)
                {
                    for (int k = 0 k < h_input.matrix_nums k++) // the k iteration for the input matrixs(could be more than 1)
                    {
                        int li = k * block + (a * kernelSet.cols + b) * converted_n + (i/stride) * convable_n + (j/stride)
                        int ri = k * padded_block + (i + a) * padded_m + j + b
                        h_temp[li] = temp[ri]
                    }
                }
            }
        }
    }

0
1
2
3
4
5
6
7
8
9
10
11
12
void gpu_gemm_cublas(const float *d_a, const float *d_b, float *d_result, int m, int n, int k)
{
    //   m*n x n*k =m *k
    cublasHandle_t handle
    cublasCreate(&handle)
    // Configure cuBLAS operations
    const float alpha = 1.0, beta = 0.0
    cublasOperation_t opA = CUBLAS_OP_N                                           // No transpose for A
    cublasOperation_t opB = CUBLAS_OP_N                                           // No transpose for B
    cublasSgemm(handle, opA, opB, k, m, n, &alpha, d_b, k, d_a, n, &beta, d_result, k) // 我日，cublas 居然是列优先的
    // cudaMemcpy(h_output.h_input, d_result, sizeof(float)*converted_n, cudaMemcpyDeviceToHost)
}

0
1
2
3
4
    int blockx = 256, blocky = 1  //和warp对应，还是需要blocksize是32的倍数， 但这个时候，设置为矮宽的形状
    dim3 blockDim(blockx, blocky)
    int gridx = (converted_n + blockx - 1) / blockx  //计算网格
    int gridy = 1  //由于左边向量的高度还是1，所以考虑块的排布也都是扁平的
    dim3 gridDim(gridx, gridy)

0
<<<BLock, Grid>> kernel_function()//我在这里搞反了顺序,Grid应该在前面

0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
stride :1
Custom Kernel Time for img2col conv: 450.193542 ms
Printing matrix of printing 0th matrix
384.000000 606.000000 723.000000 570.000000 312.000000 
318.000000 513.000000 648.000000 603.000000 354.000000 
483.000000 738.000000 873.000000 648.000000 339.000000 
318.000000 513.000000 648.000000 603.000000 354.000000 
150.000000 228.000000 291.000000 264.000000 150.000000
stride:2
Custom Kernel Time for img2col conv: 446.803680 ms
Printing matrix of printing 0th matrix
384.000000 723.000000 312.000000 
483.000000 873.000000 339.000000 
150.000000 291.000000 150.000000 
stride:3
Custom Kernel Time for img2col conv: 424.988342 ms
Printing matrix of printing 0th matrix
384.000000 570.000000 
318.000000 603.000000