Dev Cloud for oneAPI-矩阵乘法

最新推荐文章于 2024-08-20 14:59:10 发布

Attente d'amour

最新推荐文章于 2024-08-20 14:59:10 发布

阅读量51

点赞数

文章标签： oneapi 矩阵

本文链接：https://blog.csdn.net/busicooo/article/details/134744872

版权

实验环境—JupyterLab

本次实验可以在Dev Cloud中使用JupyterLab进行。

需要先进行Dev Cloud的注册。

Dev Cloud注册与使用

点击idzcn.com/devcloud.htm进行注册即可。

选择地区时中国大陆的位置不太好找，大致在这个位置。

随后需要访问Overview | Intel® DevCloud来进行服务激活，新用户可以获得两个月的免费使用期。激活并登录后，在此页面选择Get Started，来到最下方就可以看到JupyterLab的入口。

作业提交与运行

使用%%writefile 创建一个作业脚本（job script）-hello-world-example，默认是在根目录

%%writefile hello-world-example
cd $PBS_O_WORKDIR
echo "* Hello world from compute server `hostname`!"
echo "* The current directory is ${PWD}."
echo "* Compute server's CPU model and number of logical CPUs:"
lscpu | grep 'Model name\\|^CPU(s)'
echo "* Python available to us:"
which python
python --version
echo "* The job can create files, and they will be visible back in the Notebook." > newfile.txt
sleep 10
echo "*Bye"
# Remember to have an empty line at the end of the file; otherwise the last command will not run

qsub提交作业至队列
```
!qsub hello-world-example
```
qstat查询作业状态。
```
!qstat
```

S下方的字母表示状态：

Q：在队列中
R：在运行
E：错误或运行完成

结果输出：stduot和stderr会分别写在两个单独的文件中，命名格式为：
- stdout: [Job Name].o[Job ID]. Example: hello-world-example.o12345
- stderr: [Job Name].e[Job ID]. Example: hello-world-example.e12345

可以通过以下命令查看输出内容：

%ls hello-world-example*
%cat hello-world-example.o*
%cat hello-world-example.e*

简化作业提交流程—cfxmagic

正常运行python文件的流程为：

%%writefile my_job_script
echo "Running myapplication.py"
python myapplication.py

import cfxmagic包后可以直接使用%%qsub来提交作业。

%%qsub
cd $PBS_O_WORKDIR
python PythonDemo.py

这个代码块会把代码当作一个名为STDIN的作业进行提交，并不需要单独写一个运行py文件的脚本。

检测性能

以算术强度（以FLOPs/Byte测量）作为X轴，以性能（以GFLOPs/秒测量）作为Y轴绘制的，两者都以对数尺度表示。

水平线代表您的硬件在给定的时间内可以执行的给定类型的浮点计算的数量。对角线代表给定的内存子系统每秒可以传输的数据字节数。

每个点代表一个循环，一般来说，一个点离最上面的"roofs"越远，就有越大的改进空间。根据阿姆达尔定律，优化占程序总运行时间最大部分的循环将比优化占较小部分的循环带来更大的速度提升。

一个点上方的"roofs"代表阻止其达到更高性能的限制，尽管下方的"roofs"可以有所贡献。每个"roof"代表在不利用特定优化的情况下可以达到的最大性能，这与上一个"roof"有关。例如，“Scalar Add Peak"代表在不利用向量化的情况下可能达到的最大性能，如上一个"roof"所示，即"Vector Add Peak”。

下图中绿色的点D、E表示提升空间比较小，占用时间少，而红色的A、G则表示耗时很长，再结合roof，可以看出B的提升空间最大。

矩阵乘法

主函数

由于要对矩阵乘法不断优化，所以选择写一个主函数，在主函数中记录运行时长，随后只需要为不同的方法实现矩阵乘法的函数即可。

#include <CL/sycl.hpp>
#include <ctime>
#include <chrono>

using namespace sycl;

// 不同的实现方法
void mm_kernel(queue &q, std::vector<float> &matrix_a, std::vector<float> &matrix_b, std::vector<float> &matrix_c, size_t N, size_t M);

// 验证计算结果
bool almost_equal(float a, float b){
    float tolerance = 1e-6;
    float diff = fabs(a - b);
    a = fabs(a);
    b = fabs(b);
    float bigger = (b > a) ? b : a;
    if(diff <= bigger * tolerance) return true;
    return false;
}

int main(int argc, char *argv[]) {
    
    size_t N = 1024;
    size_t M = 16;
    int VERIFY = 0;
    int PRINT_OUTPUT_MATRIX = 0;

    // 定义矩阵，a、b为输入，c为输出，d用于检测
    std::vector<float> matrix_a(N*N);
    std::vector<float> matrix_b(N*N);
    std::vector<float> matrix_c(N*N);
    std::vector<float> matrix_d(N*N);
    
    // 初始化
    float v1 = 2.f;
    float v2 = 3.f;
    for (int i=0; i<N; i++)
        for (int j=0; j<N; j++){
            matrix_a[i*N+j] = v1++;
            matrix_b[i*N+j] = v2++;
            matrix_c[i*N+j] = 0.f;
            matrix_d[i*N+j] = 0.f;
    }
    
    queue q(property::queue::enable_profiling{});
    std::cout << "Offload Device        : " << q.get_device().get_info<info::device::name>() << "\n";
    std::cout << "max_work_group_size   : " << q.get_device().get_info<info::device::max_work_group_size>() << "\n";
    
    // 记录起始时间
    auto start = std::chrono::high_resolution_clock::now().time_since_epoch().count();
    
    // 调用不同方式的实现方法
    mm_kernel(q, matrix_a, matrix_b, matrix_c, N, M);
    
    // 计算运行时间
    auto duration = std::chrono::high_resolution_clock::now().time_since_epoch().count() - start;
    std::cout << "Compute Duration      : " << duration / 1e+9 << " seconds\n";
    
    // 输出结果
    if (PRINT_OUTPUT_MATRIX){
        for (int i=0; i<N; i++){
            for (int j=0; j<N; j++){
                std::cout << matrix_c[i*N+j] << " ";
            }
            std::cout << "\n";
        }
    } else {
        std::cout << " [0][0] = " << matrix_c[0] << "\n";
    }
    
    // 与普通的乘法进行比较
    if (VERIFY){
        int fail = 0;
        for(int i=0; i<N; i++){
            for (int j = 0; j < N; j++) {
                for(int k=0; k<N; k++){
                    matrix_d[i*N+j] += matrix_a[i*N+k] * matrix_b[k*N+j];
                }
                if(!almost_equal(matrix_c[i*N+j], matrix_d[i*N+j])) fail = 1;
            }
        }
        if(fail == 1){
            std::cout << "FAIL\n";
        } else {
            std::cout << "PASS\n";
        }
    }
    return 0;
}

使用KML库

void gemm(queue &exec_queue, transpose transa, transpose transb, 
          std::int64_t m, std::int64_t n, std::int64_t k, T alpha, 
          buffer<T, 1> &a, std::int64_t lda, buffer<T, 1> &b, std::int64_t ldb, 
          T beta, buffer<T, 1> &c, std::int64_t ldc)

#include <CL/sycl.hpp>
#include "oneapi/mkl/blas.hpp"  // oneMKL DPC++ 接口

using namespace sycl;

void mm_kernel(queue &q, std::vector<float> &matrix_a, std::vector<float> &matrix_b, std::vector<float> &matrix_c, size_t N, size_t M) {
    std::cout << "Configuration         : MATRIX_SIZE= " << N << "x" << N << "\n";
    
    // 创建矩阵缓冲
    buffer a(matrix_a);
    buffer b(matrix_b);
    buffer c(matrix_c);

    // 乘法引子
    float alpha = 1.f, beta = 1.f;

    // 参与乘法运算的矩阵的转置
    oneapi::mkl::transpose transA = oneapi::mkl::transpose::nontrans;
    oneapi::mkl::transpose transB = oneapi::mkl::transpose::nontrans;

    // 调用函数
    oneapi::mkl::blas::gemm(q, transA, transB, N, N, N, alpha, b, N, a, N, beta, c, N);
    c.get_access<access::mode::read>();
}

parallel_for

此时只使用了最小单位work-item。

void mm_kernel(queue &q, std::vector<float> &matrix_a, std::vector<float> &matrix_b, std::vector<float> &matrix_c, size_t N, size_t M) {
    std::cout << "Configuration         : MATRIX_SIZE= " << N << "x" << N << "\n";
    
    // 创建缓冲
    buffer a(matrix_a);
    buffer b(matrix_b);
    buffer c(matrix_c);

    // 提交指令
    auto e = q.submit([&](handler &h){
        // 获取缓冲数据
        auto A = a.get_access<access::mode::read>(h);
        auto B = b.get_access<access::mode::read>(h);
        auto C = c.get_access<access::mode::write>(h);

        // 遍历c的所有坐标
        h.parallel_for(range<2>{N,N}, [=](item<2> item){
            const int i = item.get_id(0);
            const int j = item.get_id(1);
            for (int k = 0; k < N; k++) {
                C[i*N+j] += A[i*N+k] * B[k*N+j];
            }
        });
    });
    c.get_access<access::mode::read>();
    
    // 统计运行时间
    auto kernel_duration = (e.get_profiling_info<info::event_profiling::command_end>() - e.get_profiling_info<info::event_profiling::command_start>());
    std::cout << "Kernel Execution Time : " << kernel_duration / 1e+9 << " seconds\n";
}

ND-Range提高并行性

OneAPI可以通过将数据划分为work-Item与work-Group来提高数据的并行性，如下图所示，这样划分之后就可以以work-group为单位进行并行计算，整体的数据则被称为ND-Range。

因此只需要将上述并行代码的for循环进行以下修改即可：

// 定义ND-Range和work-group的大小
range<2> global_size(N,N);
range<2> work_group_size(M,M);

// 以group为单位进行矩阵乘法
h.parallel_for(nd_range<2>{global_size, work_group_size}, [=](nd_item<2> item){
    const int i = item.get_global_id(0);
    const int j = item.get_global_id(1);
    for (int k = 0; k < N; k++) {
        C[i*N+j] += A[i*N+k] * B[k*N+j];
    }
});

运算的示意图如下：

使用局部内存（Local Memory）

工作组开始计算时，其本地内存的内容是未初始化的，并且在工作组执行完毕后不会被保留。对于许多GPU或者计算设备，本地都有专门的内存空间，因此在这些设备上，通过本地内存进行通信应该比通过全局内存进行通信的性能更好。

在ND-Range中的for循环中需要多次写入全局内存，所以可以考虑使用一个临时变量记录A和B矩阵work-group相乘的结果，在for循环执行结束后再写入C矩阵，这样就可以减少对全局内存的读写，也可以提高运行速度。

// 定义ND-Range和work-group的大小
range<2> global_size(N,N);
range<2> work_group_size(M,M);

// 以group为单位进行矩阵乘法
h.parallel_for(nd_range<2>{global_size, work_group_size}, [=](nd_item<2> item){
    const int i = item.get_global_id(0);
    const int j = item.get_global_id(1);
    float temp = 0.f;
    for (int k = 0; k < N; k++) {
        temp += A[i*N+k] * B[k*N+j];
    }
    C[i*N+j] = c;
});