求解三角系统-从理论到代码优化
理论
在求解n维线性系统 \(Ax=b\) ,我们通常将因子\(A\)分解为两个三角矩阵，即 \(A=LU\) :
- \(L\) 是下三角，其中 \(L=[l_{i,j}]\) , 并满足当\(j>i\) 时，\(l_{i,j}=0\) 和 \(l_{i,i}=1\) 。
- \(U\) 是上三角，其中 \(U=[u_{i,j}]\) ,且当 \(i>j\) 时， \(u_{i,j}=0\) 。
- 正向替换为: \(Ly=b\) 。
- 反向替换为： \(Ux=y\) 。
推导正向替换 由公式 \(Ly=b\) 展开可得：
\[\begin{array}{l} y_1&=b_1 \\ l_{2,1}y_1+y_2&=b_2 \\ l_{3,1}y_1+l_{3,2}y_2+y_3&=b_3 \\ ... \\ l_{n,1}y_1+l_{n,2}y_2+l_{n,3}y_3+...+l_{n,n-1}y_{n-1}+y_n&=b_n \end{array} \]
求出对角元素:
\[\begin{array}{l} y_1&=b_1 \\ y_2&=b_2-l_{2,1}y_1 \\ y_3&=b_3-l_{3,1}y_1-l_{3,2}y_2 \\ ...\\ y_n&=b_n-l_{n,1}y_1-l_{n,2}y_2-...l_{n,n-1}y_{n-1} \end{array} \]
公式和算法 由上面推到过程可得，对于 \(k=1,2,...,n\) ：
\[\begin{array}{r} y_k=b_k-\sum_{i=1}^{k-1}{l_{k,i}y_i} \end{array} \]
伪代码算法实现：
for k from 1 to n do:
y[k]=b[k];
for i from 1 to k-1 do:
y[k]=y[k]-l[k][i]*y[i];
反向替换过程类似代码实现和优化 准备 编译器
- MSVC2019
void solve_triangular_system_swapped( int n, double **L, double *b, double *y ){
for(int i=0;
iPipeline 优化 三种典型的Pipeline：
- 加速多条指令。
- 加速单条指令。
- 接续指令流水（指将计算机指令处理过程拆分为多个步骤，并通过多个硬件处理单元并行执行来加快指令执行速度）。
文章图片
对算法进行Pipeline分析
传统指令Pipeline
文章图片
第三种Pipeline
使 \(y_1\) 在下一个流水中可用：
文章图片
因此我们需要重写算法和公式
以 \(Ly=b\) 前五步为例：
- \(y:=b\)
- \(\begin{array}{l} y_2:=y_2-l_{2,1}*y_1;
\\ y_3:=y_3-l_{3,1}*y_1;
\\ y_4:=y_4-l_{4,1}*y_1;
\\ y_5:=y_5-l_{5,1}*y_1;
\\ \end{array}\)
- \(\begin{array}{l} y_3:=y_3-l_{3,2}*y_2;
\\ y_4:=y_4-l_{4,2}*y_2;
\\ y_5:=y_5-l_{5,2}*y_2;
\end{array}\)
- \(\begin{array}{l} y_4:=y_4-l_{4,3}*y_3;
\\ y_5:=y_5-l_{5,3}*y_3;
\end{array}\)
y:=b;
for i from 2 to n do
for j from i to n do
y_j=y_j-l[j,i-1]*y_{i-1};
由伪代码可写C语言代码：
void solve_triangular_system_swapped( int n, double **L, double *b, double *y )
{
for(int i=0;
i优化结果
VC++编译优化参数为
\Od
文章图片
总结
- 使用Pipline可以在矩阵规模较小时，带来性能上的提升。
- 矩阵规模增大时(本机为2000)，由于代码非缓存友好代码，性能将比原始代码要差。
\O1 优化的情况下原始代码有着和Pipline优化后相仿的性能，因此这里只对原始代码进行Unroll。int j = 0, i = 0;
for (i = 0;
i < n;
++i) {
y[i] = b[i];
for (j = 0;
j < i;
j += 4) {
y[i] = y[i] - y[j] * L[i][j];
y[i] = y[i] - y[j + 1] * L[i][j + 1];
y[i] = y[i] - y[j + 2] * L[i][j + 2];
y[i] = y[i] - y[j + 3] * L[i][j + 3];
}
for (;
j < i;
j++) {
y[i] = y[i] - y[j] * L[i][j];
}
}
优化结果
VC++编译优化参数为
\O1
文章图片
结论
- 当矩阵规模小时，Unroll可以带来可观的性能提升。
- 当矩阵规模增大时，当前程序主要受制于访存（由Roofline模型亦可得到相同结论），可以优化空间有限。
- 当前算法无法使用SIMD指令来获得加速增益。
void solve_triangular_system_swapped_omp( int n, double **L, double *b, double *y )
{
for(int i=0;
i推荐阅读
- win10系统下基于tensorrt的yolov5l网络部署
- 最棒的|最棒的 7 个 Laravel admin 后台管理系统推荐 - 卡拉云
- 全局视角系统学习《推荐系统》，实战中提升竞争力吾爱学习无mi
- Java实战之兼职平台系统的实现
- C语言实现学生宿舍信息管理系统课程设计
- C++实现寝室卫生管理系统
- 【|CSS 高级技巧（精灵图、字体图标、三角）
- 傻瓜笔记|2.7css精灵图 字体图标 三角 用户界面样式 布局技巧 文字溢出省略号
- CSS|CSS高级技巧——精灵图，字体图标，三角形等
- 1024程序员节|CSS精灵图/字体图标/三角/用户界面模式/vertical-align/溢出文字省略号/常见布局技巧/CSS初始化