Loading

求解三角系统-从理论到代码优化

理论

在求解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\)

从而使得求解 \(Ax=b\) ,变成了求解 \(L(Ux)=b\) :

  • 正向替换为: \(Ly=b\)
  • 反向替换为: \(Ux=y\)

直接求解 \(A\) 的复杂度为 \(O(n^3)\) ,但是求解三角系统的复杂度为 \(O(n^2)\)

推导正向替换

由公式 \(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

原始代码

根据上文中的公式,可以很快写出C语言代码如下:

void solve_triangular_system_swapped( int n, double **L, double *b, double *y ){
    for(int i=0;i<n;++i){
        y[i]=b[i];
        for(int j=0;j<i;++j){
            y[i]=y[i]-y[j]*L[i][j];
        }
    }
}

Pipeline 优化

三种典型的Pipeline:

  1. 加速多条指令。
  2. 加速单条指令。
  3. 接续指令流水(指将计算机指令处理过程拆分为多个步骤,并通过多个硬件处理单元并行执行来加快指令执行速度)。

第三种Pipline,通常情况下每个指令长度不一。
image

对算法进行Pipeline分析

传统指令Pipeline

image

第三种Pipeline

使 \(y_1\) 在下一个流水中可用:

image

因此我们需要重写算法和公式

\(Ly=b\) 前五步为例:

  1. \(y:=b\)

  2. \(\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}\)

  3. \(\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}\)

  4. \(\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<n; i++) y[i] = b[i];
    for(int i=1; i<n; i++)
    {
        for(int j=i; j<n; j++)
            y[j] = y[j] - L[j][i-1]*y[i-1];
    }
}

优化结果

VC++编译优化参数为\Od

image

总结

  1. 使用Pipline可以在矩阵规模较小时,带来性能上的提升。
  2. 矩阵规模增大时(本机为2000),由于代码非缓存友好代码,性能将比原始代码要差。

Unroll

使用循环展开对程序的性能进行优化也是常用的手段。循环展开层数取决于CPU逻辑寄存器的数量。因为Pipline在处理大矩阵时,性能较差,其次在\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
image

结论

  1. 当矩阵规模小时,Unroll可以带来可观的性能提升。
  2. 当矩阵规模增大时,当前程序主要受制于访存(由Roofline模型亦可得到相同结论),可以优化空间有限。
  3. 当前算法无法使用SIMD指令来获得加速增益。

并行优化

使用多个线程来提升性能:

void solve_triangular_system_swapped_omp( int n, double **L, double *b, double *y )
{
    for(int i=0; i<n; i++) y[i] = b[i];
    for(int i=1; i<n; i++)
    {
        #pragma omp parallel shared(L,y) private(j)
        {
            #pragma omp for
            for(int j=i; j<n; j++)
                y[j] = y[j] - L[j][i-1]*y[i-1];
        }
    }
}
posted @ 2022-03-17 08:04  PcDack  阅读(184)  评论(0编辑  收藏  举报