CUDA矩阵相乘

这几天研究了一下CUDA,发现其并行的思想和普通的CPU多线程思想不太一致,但还是挺不错。主要是将任务划分成一个个block,然后每个block里面再划分成细的线程。然后每个线程做自己做的
事情。这种并行思想很适用于像矩阵运算这些元素与元素之间的运算并不耦合得很厉害,但整体数据很大的情况,这只是我对CUDA的初步感觉。
矩阵相乘的CPU程序如下:

//C = A*B
void MatrixMulCPU(float* _C,const float *_A,const float *_B,int _wa,int _ha,int _wb)
{
    
float sum = 0;
    
for (int i = 0; i < _ha; ++i)
    {
        
for (int j = 0; j < _wb; ++j)
        {
            sum 
= 0;
            
for (int k = 0; k < _wa; ++k)
            {
                sum 
+= (float)_A[i*_wa+k]*(float)_B[k*_wb+ j];
            }
            _C[i
*_wb+j] = (float)sum;
        }
    }
}

从上面可以看出,C(i,j) = sum { A(i,k)*B(k,j) } 0<=k < _wa;耦合程度很小,所以我们可以通过划分区域的方法,让每个线程负责一个区域。
怎么划分呢?首先最初的想法是让每一个线程计算一个C(i,j),那么估算一下,应该需要height_c*width_c,也就是ha*wb个线程。进一步,我们将矩阵按一个大方格Grid划分,如果一个
方格Grid大小是16*16,那么矩阵80*48的可以表示为5(*16) * 3(*16),即16*16个大格子(block),每一个格子内,自然就是(height_c/16) *(width_c/16)个线程了。
好了,划分完后,内核代码如下:
计算版本0:
__global__ void matrix_kernel_0(float* _C,const float* _A,const float *_B,int _wa,int _wb)
{
    
float sum = 0;
    
//找出该线程所在的行列
    int row = blockIdx.y*blockDim.y + threadIdx.y;
    
int col = blockIdx.x*blockDim.x + threadIdx.x;

    
//线程Thread(row,col)负责计算C(row,col)
    for (int i = 0; i < _wa; ++i)
    {
        sum 
+= _A[row*_wa + i]*_B[i*_wb + col];
    }
    _C[row
*_wb + col] = sum;
}

另外一种思路,我们不让每一个线程完整计算一个C(i,j),通过C(i,j) = sum { A(i,k)*B(k,j) }发现,我们还可以再细度划分:
Csub(i,j) = sum{A(i,ksub+offsetA)*B(ksub+offsetB,j)}  0<=ksub < blockSize
C(i,j) = sum{Csub(i,j)}
就是把矩阵分成n*n个大的子块,然后每一个block负责计算子块i 和 子块j的子乘积,计算完毕后加起来则可。这里主要使用了共享显存作优化。

计算版本1:
__global__ void matrix_kernel_1(float* _C,const float* _A,const float *_B,int _wa,int _wb)
{
    
int bx = blockIdx.x;
    
int by = blockIdx.y;
    
int tx = threadIdx.x;
    
int ty = threadIdx.y;

    
//该block要处理的A
    int aBegin = _wa*(by*BLOCK_SIZE);//A(0,by)
    int aEnd = aBegin + _wa - 1;
    
int aStep = BLOCK_SIZE;//offsetA

    
int bBegin = BLOCK_SIZE*bx;//B(bx,0)
    int bStep = BLOCK_SIZE*_wb;//offsetB
    
    
float cSub = 0;
    
for (int a = aBegin,b = bBegin; a <= aEnd; a += aStep,b += bStep)
    {
        __shared__ 
float As[BLOCK_SIZE][BLOCK_SIZE];
        __shared__ 
float Bs[BLOCK_SIZE][BLOCK_SIZE];
        
//每个线程负责一个元素拷贝
        As[ty][tx] = _A[a + _wa*ty + tx];
        Bs[ty][tx] 
= _B[b + _wb*ty + tx];

        __syncthreads();
        
        
//每个线程负责计算一个子块i 和 子块j的子乘积
        for (int k = 0; k < BLOCK_SIZE; ++k)
        {
            cSub 
+= As[ty][k]*Bs[k][tx];
        }

        __syncthreads();
    }

    
//全局地址,向全局寄存器写回去
    
//一个线程负责一个元素,一个block负责一个子块
    int cIndex = (by*BLOCK_SIZE + ty)*_wb + (bx*BLOCK_SIZE + tx);
    _C[cIndex] 
= cSub;
}


最后写一个面向Host的接口函数:

void matrixMulGPU(float* _C,const float *_A,const float *_B,int _wa,int _ha,int _wb)
{
    
float* d_a = myNewOnGPU<float>(_wa*_ha);
    
float* d_b = myNewOnGPU<float>(_wb*_wa);
    
float* d_c = myNewOnGPU<float>(_wb*_ha);
    copyFromCPUToGPU(_A,d_a,_wa
*_ha);
    copyFromCPUToGPU(_B,d_b,_wb
*_wa);
    dim3 threads(BLOCK_SIZE,BLOCK_SIZE);
    dim3 blocks(WC
/BLOCK_SIZE,HC/BLOCK_SIZE);
    matrix_kernel_0
<<<blocks,threads>>>(d_c,d_a,d_b,_wa,_wb);
    cudaThreadSynchronize();
    copyFromGPUToCPU(d_c,_C,_wb
*_ha);

    myDeleteOnGPU(d_a);
    myDeleteOnGPU(d_b);
    myDeleteOnGPU(d_c);
}


调用的主函数如下:
#include <stdio.h>
#include 
<cuda_runtime.h>
#include 
<cutil.h>
#include 
<cutil_inline.h>
#include 
<stdlib.h>
#include 
<time.h>
#include 
<math.h>
#include 
<string.h>
#include 
<Windows.h>
#include 
"CUDACommon.h"
#include 
"MatrixMulCPU.h"
#include 
"MatrixMulGPU.h"

void randomInit(float* _data,int _size)
{
    
for (int i = 0; i < _size; ++i)
    {
        _data[i] 
= rand()/(float)RAND_MAX;
    }
}

bool checkError(const float* _A,const float* _B,int _size)
{
    
for (int i = 0 ; i < _size; ++i)
    {
        
if (fabs(_A[i] - _B[i]) > 1.0e-3)
        {
            printf(
"%f \t %f\n",_A[i],_B[i]);
            
return false;
        }
    }
    
return true;
}

int main(int argc, char* argv[])
{
    srand(
13);
    
if(!InitCUDA()) {
        
return 0;
    }

    
float* A = myNewOnCPU<float>(WA*HA);
    
float* B = myNewOnCPU<float>(WB*HB);
    randomInit(A,WA
*HA);
    randomInit(B,WB
*HB);
    
float* C = myNewOnCPU<float>(WC*HC);
    memset(C,
0,sizeof(float)*WC*HC);
    
    
float* C2 = myNewOnCPU<float>(WC*HC);
    memset(C2,
0,sizeof(float)*WC*HC);
    
    unsigned 
int tick1 = GetTickCount();
    MatrixMulCPU(C2,A,B,WA,HA,WB);
    printf(
"CPU use Time : %dms\n",GetTickCount() - tick1);
    unsigned 
int timer = 0;
    cutilCheckError(cutCreateTimer(
&timer));
    cutilCheckError(cutStartTimer(timer));
    {
        matrixMulGPU(C,A,B,WA,HA,WB);
    }
    cutilCheckError(cutStopTimer(timer));
    printf(
"GPU use time: %f (ms) \n", cutGetTimerValue(timer));
    cutilCheckError(cutDeleteTimer(timer));

    
if (checkError(C,C2,WC*HC))
    {
        printf(
"Accept\n");
    }
    
else
    {
        printf(
"Worng Answer\n");
    }

    myDeleteOnCPU(A);
    myDeleteOnCPU(B);
    myDeleteOnCPU(C);
    myDeleteOnCPU(C2);

    
return 0;
}

运算结果如下:
版本0:



版本1:


可以看出,GPU并行性能比CPU好很多,而且版本1优于版本0
### 关于CUDA矩阵相乘的示例代码 #### CPU实现 在传统的CPU上执行矩阵相乘的操作相对简单,主要通过三重循环遍历两个输入矩阵并计算其对应位置元素的成绩之和得到结果矩阵。然而,在现代高性能计算环境中,利用GPU加速成为提升效率的关键手段之一。 #### GPU实现概述 为了充分发挥NVIDIA GPU的强大算力,采用CUDA技术编写程序能够显著加快大规模数据处理速度。具体来说,对于矩阵相乘任务而言,可以通过设计合理的线程结构以及有效地管理设备端资源(如共享内存),从而达到优化性能的目的[^1]。 #### 第一步:添加CUDA内存传输 当涉及到主机与设备间的数据交换时,`cudaMemcpy()` 函数扮演着重要角色。它负责将待运算的数据从宿主机复制到显卡上的特定地址空间内;同样地,在完成所有必要的计算之后再把最终成果搬移回来供后续分析使用。这里需要注意的是,针对不同维度大小的对象可能还需要考虑边界条件等问题以确保正确无误地搬运全部所需资料[^3]。 #### 第二步:定义Kernel函数 核心部分在于构建适用于目标硬件架构特点下的核函数(Kernel),即由多个独立工作单元共同协作完成同一项指令集的过程描述符。下面给出了一段简化版的C++风格伪代码表示形式: ```cpp __global__ void matrixMulKernel(float* C, float * A, float * B, int numARows, int numAColumns, int numBRows, int numBColumns){ // 计算当前线程对应的行列索引值 int row = blockIdx.y * blockDim.y + threadIdx.y; int col = blockIdx.x * blockDim.x + threadIdx.x; if (row < numARows && col < numBColumns) { float value = 0; // 执行实际乘积累加操作 for(int k = 0; k < numAColumns ; ++k ) value += A[row*numAColumns+k]*B[k*numBColumns+col]; C[row*numBColumns+col]=value; } } ``` 这段代码展示了最基础版本的矩阵相乘逻辑,其中涉及到了二维网格布局的选择、局部变量声明及其初始化过程等内容[^4]。 #### 第三步:调用Kernel 最后就是安排好启动参数配置细节并向驱动层发出请求让上述编写的kernel得以运行起来。这通常意味着指定合适的block尺寸及grid规模以便充分利用可用流处理器数量的同时兼顾负载均衡方面的要求: ```cpp dim3 threadsPerBlock(16, 16); dim3 blocksPerGrid((numColsA + threadsPerBlock.x - 1)/threadsPerBlock.x, (numRowsA + threadsPerBlock.y - 1)/threadsPerBlock.y); matrixMulKernel<<<blocksPerGrid, threadsPerBlock>>>(d_C,d_A,d_B,numRowsA,numColsA,numRowsB,numColsB); ``` 此处选择了\(16 \times 16\) 的线程块大小作为默认选项,并据此推导出了整个作业所需的总区块数目。当然,根据实际情况调整这些超参往往可以获得更好的实验效果[^2]。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值