one api并行矩阵乘法

题目

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

并行矩阵乘法是通过将矩阵的乘法操作分解为多个子任务，并同时在多个处理单元上执行这些子任务来提高计算性能的一种方法。其原理可以简要描述如下：
任务分解： 将大矩阵的乘法任务分解为多个小矩阵块的乘法操作。例如，将一个大矩阵A分解为小矩阵块A_ij。
并行计算： 将这些小矩阵块的乘法操作分配给不同的处理单元（例如，多个 CPU 核心或 GPU），使它们可以同时执行。每个处理单元负责处理一个或多个小矩阵块的乘法。
结果合并： 将各个处理单元计算得到的小矩阵块的乘法结果合并以得到最终的矩阵乘法结果。这通常涉及到对结果矩阵的相应元素进行累加操作。

题解

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

1. 分配内存：在主机端分配内存空间用于存储输⼊矩阵和输出矩阵，同时在GPU端分配内存空间用于存储相应
   的输入和输出数据。
2. 数据传输：将输入矩阵数据从主机端内存传输到GPU端内存中。
3. 核函数调用：在SYCL中，矩阵乘法的计算通常会在GPU上使用核函数来实现并行计算。核函数
   会分配线程块和线程来处理不同的数据块。
4. 并行计算：在核函数中，每个线程负责计算输出矩阵的⼀个单独的元素。为了最大限度地利用
   GPU的并行计算能力，通常会使用⼆维线程块和线程网格的方式来处理矩阵的乘法计算。
5. 数据传输：计算完成后，将输出矩阵数据从GPU端内存传输回主机端内存中，以便进⼀步处理或
   分析。
   在并行计算矩阵乘法时，可以利用线程块和线程的层次结构来优化计算。通过合理划分矩阵数据并利用共享内
   存来减少全局内存访问的次数，可以⼤幅提高计算效率。此外，还可以利用GPU上的多个计算单元并执行行矩
   阵乘法，进⼀步提高计算速度。

multiplyMatrices用于执行矩阵乘法的核心逻辑。该函数使用了SYCL库，创建了一个队列对象，指定了使用CPU作为执行设备。该函数还创建了三个缓冲区对象，分别对应对象的A，B，C属性，用于在主机和设备之间传输数据。该函数使用了submit方法，向队列提交了一个任务，该任务包含了一个处理器对象和一个并行执行的内核函数。该内核函数使用了parallel_for方法，根据输出矩阵的行数和列数，为每个元素分配了一个索引。该内核函数还使用了get_access方法，获取了缓冲区对象的访问权限，用于读取和写入数据。该内核函数的主要逻辑是，对于每个输出矩阵的元素，计算其对应的两个输入矩阵的行和列的点积，并将结果存储在输出矩阵的相应位置。

代码

#include <CL/sycl.hpp>
#include <iostream>
#include <fstream>
#include <vector>

const double epsilon = 1e-5;
const char* inputFilename = "test.txt";
const char* outputFilename = "output.txt";

class MatrixMultiplier {
public:
    MatrixMultiplier(int m, int k, int n) : M(m), K(k), N(n), A(M, std::vector<double>(K)), B(K, std::vector<double>(N)), C(M, std::vector<double>(N)) {}

    void readMatrices(std::ifstream& input) {
        readMatrix(input, A);
        readMatrix(input, B);
    }

    void multiplyMatrices() {
        sycl::queue queue(sycl::cpu_selector{});

        sycl::buffer<double, 2> buffer_A(A.data(), sycl::range<2>(M, K));
        sycl::buffer<double, 2> buffer_B(B.data(), sycl::range<2>(K, N));
        sycl::buffer<double, 2> buffer_C(C.data(), sycl::range<2>(M, N));

        queue.submit([&](sycl::handler& cgh) {
            auto access_A = buffer_A.get_access<sycl::access::mode::read>(cgh);
            auto access_B = buffer_B.get_access<sycl::access::mode::read>(cgh);
            auto access_C = buffer_C.get_access<sycl::access::mode::write>(cgh);

            cgh.parallel_for<class matrix_multiplication>(sycl::range<2>(M, N), [=](sycl::item<2> index) {
                double sum = 0.0;
                for (int k = 0; k < K; ++k) {
                    sum += access_A[index[0]][k] * access_B[k][index[1]];
                }
                access_C[index[0]][index[1]] = sum;
            });
        });
    }

    void writeResult(std::ofstream& output) {
        for (int i = 0; i < M; ++i) {
            for (int j = 0; j < N; ++j) {
                output << C[i][j] << ' ';
            }
            output << '\n';
        }
    }

private:
    int M, K, N;
    std::vector<std::vector<double>> A, B, C;

    void readMatrix(std::ifstream& input, std::vector<std::vector<double>>& matrix) {
        for (int i = 0; i < M; ++i) {
            for (int j = 0; j < N; ++j) {
                input >> matrix[i][j];
            }
        }
    }
};

int main() {
    std::ifstream inputFile(inputFilename);
    std::ofstream outputFile(outputFilename);

    if (!inputFile.is_open() || !outputFile.is_open())
        return -1;

    int M, K, N;
    inputFile >> M >> K >> N;

    MatrixMultiplier matrixMultiplier(M, K, N);
    matrixMultiplier.readMatrices(inputFile);
    matrixMultiplier.multiplyMatrices();
    matrixMultiplier.writeResult(outputFile);

    inputFile.close();
    outputFile.close();

    return 0;
}