基于oneAPI的C++/SYCL的并⾏矩阵乘法

一、题目描述

描述

        编写⼀个基于oneAPI 的 C++/SYCL 程序来执行矩阵乘法操作。需要考虑大尺寸矩阵的乘法操作以及不同线程之间的数据依赖关系。通常在实现矩阵乘法时，可以使用块矩阵乘法以及共享内存来提高计算效率。

分析

        利用基于SYCL 的编程模型在 GPU 上实现矩阵乘法的计算，步骤如下：

        1. 分配内存：在主机端分配内存空间用于存储输⼊矩阵和输出矩阵，同时在 GPU 端分配内存空间用于存储相应的输入和输出数据。

        2. 数据传输：将输入矩阵数据从主机端内存传输到 GPU 端内存中。

        3. 核函数调用：在 SYCL 中，矩阵乘法的计算通常会在 GPU 上使用核函数来实现并行计算。核函数 会分配线程块和线程来处理不同的数据块。

        4. 并行计算：在核函数中，每个线程负责计算输出矩阵的⼀个单独的元素。为了最大限度地利用 GPU的并行计算能力，通常会使用⼆维线程块和线程网格的方式来处理矩阵的乘法计算。

        5. 数据传输：计算完成后，将输出矩阵数据从 GPU 端内存传输回主机端内存中，以便进⼀步处理或 分析。

        在并行计算矩阵乘法时，可以利用线程块和线程的层次结构来优化计算。通过合理划分矩阵数据并利用共享内存来减少全局内存访问的次数，可以⼤幅提高计算效率。此外，还可以利用GPU 上的多个计算单元并执行行矩阵乘法，进⼀步提高计算速度。

二、解题思路

        设相乘矩阵的维度分别为M*K 和K*N，故结果矩阵的维度为M*N，设立设置M*N个进程并行执行，第i, j个线程的任务则是计算第一个矩阵的第i行与第二个矩阵的第j列进行矩阵相乘运算所得结果，并将结果存储。

三、代码

#include<CL/sycl.hpp>
#include <iostream>
#include <fstream>
#include <string>
#include <fstream>
#include <iomanip>

// Matrix dimensions
constexpr size_t M = 44;
constexpr size_t N = 50;
constexpr size_t K = 96;

using namespace std;
using namespace sycl;

// Helper function to initialize matrices with random values
void initializeMatrix(float* matrix, size_t rows, size_t cols, string filename) {
    ifstream infile(filename);
    // 逐行读取文件内容
    string line;
    int cnt = 3;
    while(cnt--) {
        getline(infile, line);
    }
    int index = 0;
    while (getline(infile, line)) {
        istringstream iss(line);
        double value;
        while (iss >> value) {
            matrix[index++] = value;
        }
    }
    infile.close();
}

int main() {
    // Allocate host memory for matrices
    float *matrixA = new float[M * K];
    float *matrixB = new float[K * N];
    float *result = new float[M * N];

    string filename1 = "problem-1-AxB.txt";
    string filename2 = "problem-1-AxB_1.txt";

    // Initialize matrices with random values
    initializeMatrix(matrixA, M, K, filename1);
    initializeMatrix(matrixB, K, N, filename2);

    sycl::queue q;

    // Allocate device memory for matrices
    sycl::buffer<float, 2> bufferA(matrixA, sycl::range<2>{M, K});
    sycl::buffer<float, 2> bufferB(matrixB, sycl::range<2>{K, N});
    sycl::buffer<float, 2> bufferResult(result, sycl::range<2>{M, N});

    // Submit a SYCL command group for matrix multiplication
    q.submit([&](sycl::handler &h) {
        // Accessors to matrices
        auto accessorA = bufferA.get_access<sycl::access::mode::read>(h);
        auto accessorB = bufferB.get_access<sycl::access::mode::read>(h);
        auto accessorResult = bufferResult.get_access<sycl::access::mode::write>(h);

        // Define a range representing the work items in a 2D grid
        sycl::range<2> globalRange{M, N};

        // Execute the kernel
        h.parallel_for<class MatrixMultiply>(globalRange, [=](sycl::id<2> idx) {
            float sum = 0.0f;
            for (size_t k = 0; k < K; ++k) {
                sum += accessorA[idx[0]][k] * accessorB[k][idx[1]];
            }
            accessorResult[idx] = sum;
        });
    }).wait(); // Wait for the kernel to finish

    // Transfer results back to host

    std::ofstream outputFile("problem-1-result.txt");
    for (size_t i = 0; i < M; ++i) {
        for (size_t j = 0; j < N; ++j) {
            outputFile << std::fixed << std::setprecision(2) << result[i * N + j];
            outputFile << " ";
        }
        outputFile << "\n";
    }
    outputFile.close();
    return 0;
}