英特尔oneAPI—高斯消元算法并行化
英特尔oneAPI简介
Intel oneAPI是一个跨行业、开放、基于标准的统一的编程模型，旨在提供一个适用于各类计算架构的统一编程模型和应用程序接口。也就是说，应用程序的开发者只需要开发一次代码，就可以让代码在跨平台的异构系统上执行，底层的硬件架构可以是CPU、GPU、FPGA、神经网络处理器等。由此可见，使用oneAPI编写的程序既可以利用加速器提高程序性能，又具有可移植性。
一个oneAPI运行环境由一个主机和一系列设备组成。主机通常是一个多核CPU，而设备是一个或多个GPU、FPGA，或是其他加速器。主机的处理器也可以进行并行计算。
oneAPI为一系列的数据并行加速器提供了一个通用的开发者接口（见下图）。
英特尔DevCloud简介
英特尔DevCloud是一个可以在线开发oneAPI程序的平台，DevCloud除了预装了oneAPI开发套件之外，还提供了有关oneAPI的教程，并且免费提供GPU、FPGA等加速器资源供我们使用，因此我们可以很方便地在DevCloud上学习oneAPI知识，并测试我们自己开发的oneAPI程序。
问题描述
高斯消元法
高斯消元法(Gaussian elimination)是求解线性方阵组的一种算法，它也可用来求矩阵的秩，以及求可逆方阵的逆矩阵。它通过逐步消除未知数来将原始线性系统转化为另一个更简单的等价的系统。它的实质是通过初等行变化(Elementary row operations)，将线性方程组的增广矩阵转化为行阶梯矩阵(row echelon form)。
算法分析
高斯消去的计算模式如上图所示，在第k步时，首先将第k行除以A[k,k]，此时A[k,k] = 1，然后再将第k行到第n行减去第一行乘以A[i,k]，此时位于A[k,k]下方的所有元素全部为0。总共进行n步操作后，矩阵A就变成了对角线元素全为1的上三角矩阵。
算法伪代码如下：
procedure LU (A)
beginfor k := 1 to n dofor j := k+1 to n doA[k, j ] := A[k, j]/A[k, k];endfor;A[k, k] := 1.0;for i := k + 1 to n dofor j := k + 1 to n doA[i, j ] := A[i, j ] − A[i, k]×A[k, j];endfor;A[i, k] := 0;endfor;endfor;
end LU算法设计
串行算法
串行算法的代码与伪代码类似，具体代码如下：
for (int k = 0; k < n; k++) {for (int j = k + 1; j < n; j++) {m[k][j] = m[k][j] / m[k][k];}m[k][k] = 1;for (int i = k + 1; i < n; i++) {for (int j = k + 1; j < n; j++) {m[i][j] = m[i][j] - m[i][k] * m[k][j];}m[i][k] = 0;}
}并行算法
观察串行算法，可以发现在第k步中包含两个部分，第一个部分是一个1重循环，作用是将第k行每个元素都除以A[k,k]，第二个部分是一个2重循环，作用是将右下角的k*k大小的矩阵的第一列全部消成0。我们可以对这两个部分进行并行化处理。
基于buffer实现并行算法
首先我们创建一个任务队列，用于给加速器提交任务：
queue q;我们还可以在创建队列时，将队列绑定到不同设备上，例如：
- queue(cpu_selector{});
- queue(gpu_selector{});
可以由此指定我们的并行代码要在什么设备上执行，我们不需要了解GPU/FPGA等加速器具体的编程语言，就可以将计算任务以统一的编程方式卸载到这些加速器上。
然后创建一个n*n大小的矩阵：
const int n = 1024;
buffer<float, 2> buf(range(n, n));针对第一个部分，我们可以使用parallel_for语句将其卸载到加速器上并行计算，具体代码如下：
q.submit([&](handler& h) {accessor m{ buf, h, read_write };h.parallel_for(range(n - k), [=](auto idx) {int j = k + idx;m[k][j] = m[k][j] / m[k][k];});
});针对第二个部分，我们仍然可以使用parallel_for，不过需要注意的是，这里使用的不再是一维的数据划分，而是二维的数据划分，具体代码如下：
q.submit([&](handler& h) {accessor m{ buf, h, read_write };h.parallel_for(range(n - (k + 1), n - (k + 1)), [=](auto idx) {int i = k + 1 + idx.get_id(0);int j = k + 1 + idx.get_id(1);m[i][j] = m[i][j] - m[i][k] * m[k][j];});
});除此之外，还有一点需要注意，由于A[i,k]元素在此部分中会被所有线程使用，因此不能在此阶段置0，需要等这部分全部执行完成之后再置0。因此需要补充第三个阶段，将右下角k*k大小的矩阵的第一列除A[k,k]以外的元素置0，具体代码如下：
q.submit([&](handler& h) {accessor m{ buf, h, read_write };h.parallel_for(range(n - (k + 1)), [=](auto idx) {int i = k + 1 + idx;m[i][k] = 0;});
});在最外层n重循环执行完之后，其实只是将计算任务下发到了加速器， 这些计算任务不一定真的执行完了，因此我们需要手动调用q.wait()，等待任务队列中所有任务都执行完毕，该算法再结束。
完整代码如下：
void gauss_oneapi(buffer<float, 2>& buf, queue& q) {device my_device = q.get_device();std::cout << "Device: " << my_device.get_info<info::device::name>() << std::endl;int n = buf.get_range()[0];for (int k = 0; k < n; k++) {q.submit([&](handler& h) {accessor m{ buf, h, read_write };h.parallel_for(range(n - k), [=](auto idx) {int j = k + idx;m[k][j] = m[k][j] / m[k][k];});});q.submit([&](handler& h) {accessor m{ buf, h, read_write };h.parallel_for(range(n - (k + 1), n - (k + 1)), [=](auto idx) {int i = k + 1 + idx.get_id(0);int j = k + 1 + idx.get_id(1);m[i][j] = m[i][j] - m[i][k] * m[k][j];});});q.submit([&](handler& h) {accessor m{ buf, h, read_write };h.parallel_for(range(n - (k + 1)), [=](auto idx) {int i = k + 1 + idx;m[i][k] = 0;});});}q.wait();
}基于USM实现并行算法
USM的全称是Unified Shared Memory，程序员可以不管数据是处在主机的内存当中还是加速器的内存当中，仍然使用C/C++中的指针来访问数据，它给程序员提供了一个统一的内存视图（如下图所示），有助于程序员更快地将现有C/C++程序移植到DPC++。
但是和buffer+accessor的访问方式相比，如果对数据依赖解决不正确，可能会存在竞争，导致计算结果错误。
在申请矩阵内存时，我们可以使用malloc_shared函数：
float* buf = (float *)malloc_shared(n * n * sizeof(float), q);在访问内存时，则不需要使用accessor，直接访问buf即可。需要注意由于不再使用accessor，因此oneAPI不能推断出数据依赖，所以需要程序员手动指出数据依赖，其中有3种方式：
- q.wait()
- in_order队列属性
- h.depends_on(e)方法
由于本算法中所有kernel都是依次执行的，因此在创建队列时指定in_order属性比较方便。
queue q{property::queue::in_order()};具体算法如下：
void gauss_oneapi(float* m, int n, queue& q) {device my_device = q.get_device();std::cout << "Device: " << my_device.get_info<info::device::name>() << std::endl;for (int k = 0; k < n; k++) {q.submit([&](handler& h) {h.parallel_for(range(n - k), [=](auto idx) {int j = k + idx;m[k*n+j] = m[k * n + j] / m[k * n + k];});});q.submit([&](handler& h) {h.parallel_for(range(n - (k + 1), n - (k + 1)), [=](auto idx) {int i = k + 1 + idx.get_id(0);int j = k + 1 + idx.get_id(1);m[i * n + j] = m[i * n + j] - m[i * n + k] * m[k * n + j];});});q.submit([&](handler& h) {h.parallel_for(range(n - (k + 1)), [=](auto idx) {int i = k + 1 + idx;m[i * n + k] = 0;});});}q.wait();
}
实验设计
串行算法与基于buffer实现的并行算法的性能对比
为了对比串行算法与并行算法的性能差异，设置了不同规模的矩阵，分别测量在不同问题规模下，串行算法与并行算法的性能表现。
基于buffer实现的并行算法与基于USM实现的并行算法的性能对比
为了对比基于buffer实现的并行算法与基于USM实现的并行算法之间的性能差异，也应测量在不同问题规模下各算法的性能表现。
时间测量方法
时间测量
我们使用C++11标准中的std::chrono::high_resolution_clock测量时间，而不使用平台特定的API(如Windows平台的QueryPerformanceCounter、Linux平台的clock等)，这样可以增加程序的可移植性。
除此之外，还应使用多次测量求平均的方式保证运行时间测量的准确性。
auto start = std::chrono::high_resolution_clock::now();// your codeauto end = std::chrono::high_resolution_clock::now();预热
为了防止初始化操作、cache等其他因素影响算法执行时间的测量，因此在开始计时之前，先让算法运行一次，第一次不计入运行时间，然后再运行k次，以这k次运行的平均时间作为算法执行时间。
结果分析
测试环境
性能测试部分全部在DevCloud平台上完成。
GPU型号：Intel(R) UHD Graphics P630
CPU型号：Intel(R) Xeon(R) E-2176G CPU @ 3.70GHz
串行算法与基于buffer实现的并行算法的性能对比
| 问题规模n | 串行算法 | GPU并行算法 |
| 16 | 0.470433 | 1.42439 |
| 32 | 0.812061 | 2.292 |
| 64 | 6.41777 | 4.03059 |
| 128 | 50.598 | 5.60625 |
| 256 | 406.665 | 7.42548 |
| 512 | 3236.67 | 75.7885 |
| 1024 | 25885.7 | 590.462 |
当问题规模n非常小的情况下，串行算法的效率是最高的。因为如果我们使用GPU做并行化，那么就需要将矩阵从主机的内存复制到GPU的内存中，然后再让GPU开始计算，等GPU计算完之后，还要把矩阵再从GPU的内存复制到主机的内存中。在问题规模n非常小的情况下，这些额外的开销占比是非常大的，因此在此时GPU版本算法要比串行算法的执行时间更长。
当问题规模n比较大时，GPU并行算法耗时更短，比较符合预期效果。
基于buffer实现的并行算法与基于USM实现的并行算法的性能对比
| 问题规模n | GPU并行算法（buffer） | GPU并行算法（USM） |
| 16 | 1.42439 | 1.59059 |
| 32 | 2.292 | 2.02017 |
| 64 | 4.03059 | 2.9549 |
| 128 | 5.60625 | 4.81008 |
| 256 | 7.42548 | 5.91026 |
| 512 | 75.7885 | 89.6216 |
| 1024 | 590.462 | 616.103 |
由上表可以看出，基于buffer的GPU并行算法与基于USM的GPU并行算法之间性能差异不大，因此程序员可以使用USM更快地将现有C/C++程序移植到DPC++，并且相比于使用buffer+accessor，不会带来很多额外的性能开销。
代码
zouxianyu / oneAPI gaussian elimination · GitCode
参考资料
Get Started | Intel® DevCloud
高斯消元法详解_ityanger的博客-CSDN博客_高斯消元法
英特尔oneAPI—高斯消元算法并行化相关推荐
- Rocksdb Ribbon Filter : 结合 XOR-filter 以及 高斯消元算法 实现的 高效filter
文章目录 前言 XOR-filter 实现原理 xor filter 的构造原理 xor filter 构造总结 XOR-filter 和 ADD-filter对比 XOR-filter 在计算上的优 ...
- 数论 - 高斯消元算法
1.高斯消元 (1)定义 高斯消元法是求解线性方阵组的一种算法,它也可用来求矩阵的秩,以及求可逆方阵的逆矩阵.它通过逐步消除未知数来将原始线性系统转化为另一个更简单的等价的系统.它的实质是通过初等行变 ...
- 【C++】高斯消元算法
矩阵初等行变换法则 任一行可以与另一行进行加减. 任一行可以乘或除以一个非零常数(除其实就是乘一个倒数). 任两行可以交换位置. 线性方程组 形如 a1,1x1+a1,2x2+⋯+a1,nxn=b1a ...
- [算法] 高斯消元详解
原文链接(强烈建议安利) 0.前置知识 知道如何解三元一次方程组 有手,有脑子 1.答案的表示与存储 先解一个方程组: 2x+3y+5z=31 x-4y -z=-6 4x+2y-5z=9 我们把这个方 ...
- 计算机编程 高斯消元,高斯-若尔当消元法
高斯-若尔当消元法(英语:Gauss-Jordan Elimination),或译为高斯-约旦消元法,简称G-J消元法,是数学中的一个算法,是高斯消元法的另一个版本.它在线性代数中用来找出线性方程组的 ...
- 第三十四章 数论——高斯消元解线性方程组
第三十四章 数论--高斯消元解线性方程组 一.高斯消元 1.线性方程组 2.高斯消元步骤 (1)数学知识铺垫 增广矩阵和阶梯矩阵 初等变换 (2)高斯消元步骤 二.代码模板 1.问题: 2.代码 一. ...
- 基于高斯消元的BATS码的改进译码算法
基于高斯消元的BATS码的改进译码算法 摘要 BATS 码的编译码原理: 改进的高斯消元的结论: 系数矩阵不满秩下可以解码部分包的情况: 摘要 BATS码:(Batched Sparse Codes, ...
- 【算法设计与分析基础】高斯消元
参考<算法导论>第七部分 算法问题选编中的第28章 矩阵运算. #include<iostream> using namespace std;/*** 一些矩阵例子*//* d ...
- {算法}高斯消元不高斯
高斯消元蛋白,由开挂数学家高斯命名.它可以在 O(n2) O(n^2)的时间内解出n元1次方程组. 功能这么强大,其实原理就是常用的"加减消元法". 首先,还是要将方程标准化,都( ...
最新文章
- Oracle 12c(12.1.0.5) oem agent silent install(静默安装agent)
- 嵌入式C语言代码的调试技巧
- 怎么用latex写ppt呢？
- Openstack 与VMware 不同CPU迁移原理
- Java的io类的使用场景
- 下取整函数的含义_取整函数解读
- 大数据之_SCALA工作笔记001---Centos7.3安装scala
- stm32中的延时函数
- 毕设题目：Matlab风电功率预测
- IDEA打开README.md文件时卡死
- 一款黑科技神器：uTools
- vpa函数python_Biopython序列比对
- iPhone屏幕数据
- 马斯克的“星链”会不会威胁中国太空安全？肯定会！
- 3维图像处理的新星--Open3D（实操过程持续更新ing....
- 三分钟带你了解DCMM
- 医院挂号系统代码_智慧医院中心是怎样做的？分诊叫号系统如何正确使用！
- 君康人寿2019年排名_君康人寿保险靠谱吗？
- 计算机教室是使用计划书,教师个人计算机学习计划书_老师计算机学习计划范文...
- 解决backtrack5连接不上ssh问题