CUDA编程

__global__ void kernel3D3D(float *input, int dataNum)
{
    // thread在block中位置计算：
    int threadInBlock = threadIdx.x + threadIdx.y*blockDim.x + threadIdx.z*blockDim.x*blockDim.y;
    // block在整个grid中的位置计算：
    int blockInGrid = blockIdx.x + blockIdx.y*gridDim.x + blockIdx.z*gridDim.x*gridDim.y;
    // 一个block有多少个线程计算：
    int oneBlockSize = blockDim.x*blockDim.y*blockDim.z;
    // 位置索引：
    int idx = threadInBlock + oneBlockSize*blockInGrid;
}

 限制：必须返回void类型。
 不能是类的成员函数。
 通常用于并行计算任务的主体逻辑。

 可以有返回值
 用途：实现GPU内部的辅助函数或工具函数。
 限制：无法直接从CPU调用。

#include <stdio.h>
#include <assert.h>
#include <iostream>
#define M 64
#define K 256
#define N 128
#define MAX_ERR 1e-4

__global__ void matrix_vector_multiplication(float* vector_result, float *matrix_a, float *matrix_b, int m_row, int n_col)
{
    //__shared__ float temp[M*N]; 
    //extern __shared__ float temp[]; //matrix_vector_multiplication<<<M,N,sizeof(float) * (M * N)>>>

    // blockIdx.x => 0代表矩阵的第一行与矩阵B相乘，1代表矩阵的第二行与矩阵B相乘
    // blockDim.x => 输出维度，对应矩阵B的col
    // threadIdx.x => 当前block内某个线程的ID，就是矩阵B中某列的索引

    // Unique tid which can index each single element in the output matrix
    int tid = blockIdx.x * blockDim.x + threadIdx.x;

    // 上面这一生成索引会有tid超出矩阵边界的风险，因此需要用下面的条件语句加限制
    int size_of_the_matrix = m_row*n_col;
    if(tid<size_of_the_matrix)
    {
        for(int i=0;i<K;i++)
        {
            vector_result[tid]+=matrix_a[blockIdx.x*K+i]*matrix_b[i*N+threadIdx.x];
        }
    }

    //__syncthreads(); // synchronize all threads
}


int main(void)
{
    float *x, *y;
    float *r;
    // Allocate Unified Memory – accessible from CPU or GPU
    cudaMallocManaged(&x, M*K*sizeof(float));
    cudaMallocManaged(&y, K*N*sizeof(float));
    cudaMallocManaged(&r, M*N*sizeof(float));

    // initialize x and y arrays on the host
    for(int i=0; i<M; i++){
        for(int j=0; j<K; j++){
            x[i*K+j] = float(j);
        }
    }

    for(int i=0; i<K; i++){
        for(int j=0; j<N; j++){
            y[i*N+j] = float(i);
        }
     }

    // Run kernel on 1M elements on the GPU
    matrix_vector_multiplication<<<M,N>>>(r,x,y,M,N);

    // Wait for GPU to finish before accessing on host
    cudaDeviceSynchronize();

    float maxError = 0.0f;
    for (int i=0; i <M; i++){
        for (int j=0; j<N;j++){
            maxError = fmax(maxError, fabs(r[i*N+j]-5559680.0f));
            if(maxError!=0.0f)
            {
                std::cout << "Max error: " << maxError << std::endl;
            }
        }
     }

    // Free memory
    cudaFree(x);
    cudaFree(y);
    cudaFree(r);
    return 0;
}