CUDA编程- 矩阵乘法

青衫客36

已于 2023-10-20 13:26:50 修改

阅读量311

点赞数

文章标签：矩阵 CUDA

于 2023-10-20 12:39:09 首次发布

本文链接：https://blog.csdn.net/weixin_43844521/article/details/133943835

版权

矩阵乘法是一个经典的并行计算示例，因为每个输出元素的计算都可以独立于其他元素进行。以下是一个简单的CUDA程序，用于计算两个矩阵的乘法：

#include <iostream>
#include <cuda_runtime.h>

const int WIDTH = 16;  // 假设矩阵的宽度和高度都是16

/*
    每个线程负责计算结果矩阵C中的一个元素。
    为了计算这个元素，它会取矩阵A的一个行向量和矩阵B的一个列向量，
    进行点乘，然后将结果存入结果矩阵C的对应位置。
*/
__global__ void matrixMultiply(float* A, float* B, float* C, int width) {
    /*
        计算当前线程应该处理的元素在结果矩阵中的行和列。
        blockIdx.x 和 blockIdx.y: 当前块的ID。在一个2D的网格中，每个块都有一个(x, y)的ID。
        blockDim.x 和 blockDim.y: 每个块中的线程数。这是在调用核函数时通过<<<blocksPerGrid, threadsPerBlock>>>指定的。
        threadIdx.x 和 threadIdx.y: 当前线程在其所属的块中的ID。
    */
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

    float value = 0;
    if (row < width && col < width) {
        /*
            计算矩阵A的第row行和矩阵B的第col列的点乘
            这个循环将两个矩阵的对应元素相乘，然后累加这些乘积来计算结果矩阵C的对应位置。
        */
        for (int k = 0; k < width; ++k) {
            value += A[row * width + k] * B[k * width + col];
        }
        C[row * width + col] = value;
    }
}

int main() {
    float* h_A, * h_B, * h_C;
    float* d_A, * d_B, * d_C;

    int size = WIDTH * WIDTH * sizeof(float);

    // 分配主机内存
    h_A = new float[WIDTH * WIDTH];
    h_B = new float[WIDTH * WIDTH];
    h_C = new float[WIDTH * WIDTH];

    // 初始化矩阵数据
    for (int i = 0; i < WIDTH * WIDTH; ++i) {
        h_A[i] = rand() % 10;
        h_B[i] = rand() % 10;
    }

    std::cout << "matrix A:" << std::endl;
    for (int i = 0; i < WIDTH; ++i) {
        for (int j = 0; j < WIDTH; ++j) {
            std::cout << h_A[i * WIDTH + j] << " ";
        }
        std::cout << std::endl;
    }

    std::cout << "matrix B:" << std::endl;
    for (int i = 0; i < WIDTH; ++i) {
        for (int j = 0; j < WIDTH; ++j) {
            std::cout << h_B[i * WIDTH + j] << " ";
        }
        std::cout << std::endl;
    }

    // 分配设备内存
    cudaMalloc(&d_A, size);
    cudaMalloc(&d_B, size);
    cudaMalloc(&d_C, size);

    // 复制数据到设备
    cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice);
    cudaMemcpy(d_B, h_B, size, cudaMemcpyHostToDevice);

    // 执行核函数
    dim3 dimBlock(WIDTH, WIDTH);
    dim3 dimGrid(1, 1);
    matrixMultiply << <dimGrid, dimBlock >> > (d_A, d_B, d_C, WIDTH);

    // 复制结果回主机
    cudaMemcpy(h_C, d_C, size, cudaMemcpyDeviceToHost);

    // 打印结果
    std::cout << "Result matrix C:" << std::endl;
    for (int i = 0; i < WIDTH; ++i) {
        for (int j = 0; j < WIDTH; ++j) {
            std::cout << h_C[i * WIDTH + j] << " ";
        }
        std::cout << std::endl;
    }

    // 清理
    delete[] h_A;
    delete[] h_B;
    delete[] h_C;
    cudaFree(d_A);
    cudaFree(d_B);
    cudaFree(d_C);

    return 0;
}

这个程序首先初始化两个随机矩阵，然后在GPU上执行矩阵乘法，并将结果返回给主机。在matrixMultiply核函数中，每个线程计算输出矩阵中的一个元素。

运行结果如下：

matrix A:
1 4 9 8 2 5 1 1 5 7 1 2 2 1 8 7
1 9 7 5 3 2 3 1 3 7 2 7 3 9 6 0
8 0 4 6 0 0 6 3 9 4 6 6 1 4 6 7
8 9 3 9 4 7 3 1 4 0 7 7 6 1 5 5
1 0 1 6 7 7 7 3 5 9 1 2 6 3 0 2
0 4 8 5 2 1 4 6 8 8 1 3 8 0 7 8
8 7 0 3 9 5 0 9 8 4 4 9 5 3 3 4
8 8 0 8 9 9 2 8 8 0 8 5 6 2 2 8
2 5 9 4 1 2 1 7 2 5 0 3 1 1 5 2
2 8 7 2 9 1 6 8 5 4 1 5 8 9 9 7
9 6 3 7 6 8 5 4 6 6 3 5 2 2 3 1
6 1 6 9 5 3 4 2 3 4 5 6 4 8 7 8
2 7 4 4 2 3 7 8 3 3 2 2 3 2 7 3
3 1 3 5 5 2 2 4 7 3 0 5 1 6 5 0
5 9 6 8 1 8 1 8 9 0 0 0 0 6 3 4
2 5 7 5 9 3 2 3 7 5 4 2 0 1 1 3
matrix B:
7 0 4 8 4 5 7 1 2 6 4 3 2 6 5 6
8 2 9 4 1 3 4 1 8 4 7 9 1 8 5 2
6 2 8 5 9 0 1 8 3 4 0 6 8 9 3 8
2 1 5 8 0 6 6 5 2 9 3 2 0 6 7 4
2 0 4 0 1 7 7 3 9 8 8 6 0 8 1 5
9 7 3 1 0 6 0 1 9 4 4 9 8 8 7 3
3 1 7 4 6 1 9 6 3 8 9 2 5 3 7 3
0 8 6 1 8 7 2 2 9 7 1 8 1 4 5 6
6 3 2 6 8 1 9 6 9 2 0 9 8 9 3 5
5 6 7 2 4 9 9 2 5 9 2 0 3 3 6 9
7 9 5 6 5 2 4 1 1 3 5 3 8 5 6 8
2 4 2 1 0 6 9 9 4 0 9 3 2 8 5 1
5 4 4 4 9 0 7 7 0 8 5 0 7 9 5 4
2 1 2 3 1 2 6 0 0 4 0 3 0 7 0 4
4 5 4 2 6 3 4 4 2 7 6 7 5 7 8 9
6 7 4 7 8 0 2 6 1 5 0 6 6 1 7 0
Result matrix C:
323 250 326 271 301 225 301 287 264 368 197 331 293 405 353 334
302 210 344 236 239 253 390 264 272 355 280 309 225 470 310 339
317 260 306 340 351 222 420 307 235 360 240 309 304 407 370 349
413 281 369 352 290 278 411 311 318 409 368 388 320 530 430 338
257 215 278 202 240 266 374 236 292 389 243 239 237 366 299 285
332 304 375 294 429 229 391 360 306 431 235 352 343 448 391 371
362 310 351 282 318 356 461 292 427 415 343 429 272 537 379 364
445 361 403 376 361 351 460 319 448 482 373 494 352 584 461 381
221 195 283 182 249 200 224 212 239 283 159 272 190 322 257 277
372 312 442 307 440 295 488 374 384 517 360 458 327 588 414 429
371 251 363 297 270 337 432 266 380 426 328 374 276 495 388 367
350 281 358 347 339 280 438 335 256 453 297 343 306 497 405 388
270 241 334 226 300 221 310 239 287 367 259 330 237 376 338 305
195 155 215 182 202 216 317 210 246 279 190 261 168 359 226 265
338 235 330 296 282 242 294 225 363 349 182 436 255 459 327 304
280 198 302 229 248 236 316 232 321 323 217 321 232 390 264 303

注意：在实际应用中，可能需要使用更优化的算法，例如瓦片（tiling）技术，以更有效地利用GPU的共享内存和减少全局内存的访问。

关键代码分析1：矩阵乘法核函数matrixMultiply的调用

	dim3 dimBlock(WIDTH, WIDTH);
    dim3 dimGrid(1, 1);
    matrixMultiply << <dimGrid, dimBlock >> > (d_A, d_B, d_C, WIDTH);

在这里，我们为矩阵乘法定义了线程块的大小和网格的大小，并使用这些参数来调用核函数。

dim3 dimBlock(WIDTH, WIDTH);

dim3: CUDA提供的一个数据类型，用于定义三维的大小。尽管名为三维，但我们可以只使用其中的一部分（例如，在这里我们只使用了两维）。
dimBlock: 定义了每个线程块的大小。在此例中，每个线程块的大小是WIDTH x WIDTH，意味着每个线程块有WIDTH行和WIDTH列的线程。

dim3 dimGrid(1, 1);

dimGrid: 定义了整个网格的大小。在这里，我们的网格大小为1 x 1，意味着我们在整个网格中只有一个线程块。这样选择可能是因为WIDTH x WIDTH的线程块大小已经足够大，可以覆盖整个矩阵。

matrixMultiply << <dimGrid, dimBlock >> > (d_A, d_B, d_C, WIDTH);

这是核函数matrixMultiply的调用。<<<dimGrid, dimBlock>>>是CUDA语法，用于指定线程的网格和块的大小。在这里，核函数的网格大小为1 x 1，线程块的大小为WIDTH x WIDTH。
d_A, d_B, 和 d_C 是在GPU上的矩阵数据。WIDTH是矩阵的大小。

总之，这段代码的目的是在GPU上启动WIDTH x WIDTH个线程，所有这些线程都在一个线程块中，并执行矩阵乘法操作。每个线程都会计算结果矩阵的一个元素。

关键代码分析2：矩阵乘法的实现部分

这段代码是矩阵乘法的核心部分，用于计算结果矩阵中的每个元素。

int row = blockIdx.y * blockDim.y + threadIdx.y;
int col = blockIdx.x * blockDim.x + threadIdx.x;

blockIdx.x 和 blockIdx.y: 当前块的ID。在一个2D的网格中，每个块都有一个(x, y)的ID。
blockDim.x 和 blockDim.y: 每个块中的线程数。这是在调用核函数时通过<<<blocksPerGrid, threadsPerBlock>>>指定的。
threadIdx.x 和 threadIdx.y: 当前线程在其所属的块中的ID。

这两行代码的目的是为了计算当前线程应该处理的元素在结果矩阵中的行和列。

float value = 0;

初始化累加器value，用于存储结果矩阵的一个元素的计算结果。

if (row < width && col < width) {

这个条件语句确保当前线程对应的行和列在矩阵的范围内。这是因为，基于性能的原因，我们可能会启动更多的线程块和线程，超过了所需的数量，所以需要检查以避免数组越界。

    for (int k = 0; k < width; ++k) {
        value += A[row * width + k] * B[k * width + col];
    }

这是矩阵乘法的核心部分。我们正在计算矩阵A的第row行和矩阵B的第col列的点乘。这个循环将两个矩阵的对应元素相乘，然后累加这些乘积来计算结果矩阵C的对应位置。

    C[row * width + col] = value;
}

将计算得到的value存入结果矩阵C的对应位置。

这段代码实质上是：每个线程负责计算结果矩阵C中的一个元素。为了计算这个元素，它会取矩阵A的一个行向量和矩阵B的一个列向量，进行点乘，然后将结果存入结果矩阵C的对应位置。

几点说明：

（1）在CUDA中，blockIdx.x 和 blockIdx.y 是内置变量，它们用于表示当前线程块在整个网格（grid）中的索引。和线程块（block）的维度一样，网格也可以具有最多三个维度，即x、y和z。

blockIdx.x：表示当前线程块在网格的x维度上的索引。
blockIdx.y：表示当前线程块在网格的y维度上的索引。

同样地，还可能会看到：

blockIdx.z：表示当前线程块在网格的z维度上的索引。

假设我们创建了一个二维网格，其中包含多个线程块：

dim3 blocksInGrid(10, 5);  // 这表示整个网格的尺寸为10x5的线程块

在这种情况下，blockIdx.x的值会在0到9之间变化，而blockIdx.y的值会在0到4之间变化。

在核函数（kernel）中，使用这些值可以帮助我们确定当前线程块在整个网格中的位置。这对于索引输入和输出数据非常有用，尤其是当数据结构是多维的时候。

简而言之，当我们有一个多维的问题空间或数据结构时，blockIdx 和 blockDim 变量的x、y和z维度提供了一个直观的方式来组织和索引我们的线程和线程块。

（2）在CUDA中，线程块（block）和线程（thread）都可以具有多达三个维度，即x、y和z。这些维度提供了一种更灵活的方式来组织和索引线程，从而更好地匹配某些应用程序的数据结构和问题领域。

blockDim.x：表示线程块的x维度大小，也就是该方向上的线程数量。
blockDim.y：表示线程块的y维度大小，也就是该方向上的线程数量。

同理，还可能会看到：

blockDim.z：表示线程块的z维度大小。

例如，当我们创建一个二维线程块的网格（grid）时，可能会这样定义：

dim3 threadsPerBlock(16, 16);  // 这表示每个线程块包含16x16的线程

在此情况下，blockDim.x的值为16，blockDim.y的值也为16。

当在核函数（kernel）内部使用这些值时，它们可以帮助我们确定当前线程在线程块内的位置。例如，threadIdx.x会在0到blockDim.x-1之间变化，threadIdx.y会在0到blockDim.y-1之间变化。

选择适当的线程块大小和维度可以影响核函数的性能，因为它影响到共享内存的使用、线程块中的线程数量以及可能的寄存器使用情况。所以，选择最佳的线程块维度和大小通常需要根据特定应用程序的需要进行实验调整。

（3）在CUDA中，threadIdx.x 和 threadIdx.y 是预定义的变量，表示当前线程在其所属线程块（block）中的索引。和网格（grid）的维度一样，线程块也可以有最多三个维度，即x、y和z。

threadIdx.x：表示当前线程在线程块的x维度上的索引。
threadIdx.y：表示当前线程在线程块的y维度上的索引。

另外，还可能会遇到：

threadIdx.z：表示当前线程在线程块的z维度上的索引。

如果我们创建了一个二维线程块，比如：

dim3 threadsInBlock(16, 16);  // 这表示线程块的尺寸是16x16的线程

在这种情况下，threadIdx.x 的值会在0到15之间变化，而 threadIdx.y 的值也会在0到15之间变化。

在核函数（kernel）内部，我们通常使用 threadIdx 和 blockIdx 的组合来确定当前线程在整个数据集或问题空间中的位置。这特别适用于多维数据结构，例如矩阵或图像，其中可能会使用两个或三个维度。

总之，通过threadIdx.x, threadIdx.y, 和 threadIdx.z 我们可以确定当前线程在它所属的线程块内的位置。而这些位置信息，结合网格和线程块的维度，是在CUDA并行编程中确定每个线程应执行的数据部分的关键。

青衫客36

关注

0
点赞
踩
3

收藏

觉得还不错? 一键收藏
打赏
0
评论
CUDA编程- 矩阵乘法

为了计算这个元素，它会取矩阵A的一个行向量和矩阵B的一个列向量，进行点乘，然后将结果存入结果矩阵C的对应位置。这是因为，基于性能的原因，我们可能会启动更多的线程块和线程，超过了所需的数量，所以需要检查以避免数组越界。注意：在实际应用中，可能需要使用更优化的算法，例如瓦片（tiling）技术，以更有效地利用GPU的共享内存和减少全局内存的访问。在这里，我们为矩阵乘法定义了线程块的大小和网格的大小，并使用这些参数来调用核函数。这两行代码的目的是为了计算当前线程应该处理的元素在结果矩阵中的行和列。
复制链接

扫一扫