Stationary Iterative Method Report
注： 由于 GitHub 目前不支持 LaTeX 语法，并且对内联 HTML 的支持也很不友好，所以可能页面阅读不是很方便，建议阅读PDF版本。
雅可比法（Jacobi Method）是一种解对角元素几乎都是各行和各列的绝对值最大的值的线性方程组的算法。求解出每个对角元素并插入近似值，不断迭代直至收敛。雅克比迭代法的计算公式简单，每迭代一次只需计算一次矩阵和向量的乘法，且计算过程中原始矩阵
Gauss-Seidel 方法也称为 Liebmann 方法或连续位移方法，与雅可比法非常相似，同样是基于矩阵分解的原理。虽然它可以应用于对角线上具有非零元素的任何矩阵，但只能在矩阵是对角线主导的或对称的和正定的情况下，保证收敛。
D. M. Young 于 20 世纪 70 年代提出逐次超松弛 (Successive Over Relaxation) 迭代法，简称 SOR 方法，是一种经典的迭代算法。它是为了解决大规模系统的线性等式提出来的，在 G-S 法基础上为提高收敛速度，采用加权平均而得到的新算法。
针对形如
我们研究了很久给出的 Equation193.stiff 和 Equation4800.stiff 两个文件，但仍没有理解文件中子矩阵的分解方式，因此我们使用 Matlab 来生成实验数据。
其中，方阵
然后随机生成解向量
将所得结果保存在 txt 文件中，作为程序的输入。除此之外，程序还要去输入矩阵大小、计算精度、输出文件名，对于超松弛迭代法，还需输入超松弛常数。
- 迭代次数
- 计算时间
- 误差
- 结果（保存在 txt 文件里）
-
首先将方程组中的系数矩阵
$A$ 分解成下三角阵$L$ 、对角阵$D$ 、上三角阵型$U$ 三部分，即$A=L+D+U$ 。 -
然后确定迭代格式：$X^{(k+1)}=B_{J} X^{(k)}+f_{J}$ $$ \left{\begin{array}{l}{{x}{1}^{({k}+1)}=\left(-{a}{12} {x}{2}^{({k})} {-} {a}{13} {x}{3}^{({k})} {-} \ldots-{a}{1 {n}} {x}{{n}}^{({k})}+{b}{1}\right) / {a}{11}} \ {{x}{2}^{({k}+1)}=\left(-{a}{21} {x}{1}^{({k})}-{a}{23} {x}{3}^{({k})} {-} \ldots {a}{2 {n}} {x}{{n}}^{({k})}+{b}{2}\right) / {a}{22}} \ {\cdots} \ {{x}{{n}}^{({k}+1)}=\left(-{a}{{n} 1} {x}{1}^{({k})}-{a}{{n} 2} {x}{2}^{({k})}{-} \ldots-{a}{{n}, {n}-1} {x}{{n}}^{({k})}+{b}{{n}}\right) / {a}{{nn}}}\end{array}\right. $$ 也即 $$ x{i}^{(k+1)}=\frac{1}{a_{i i}} \sum_{j=1}^{n} a_{i j} x_{j}^{(k)}, i=1,2, \ldots, n $$
-
最后选取初始迭代向量
$x^{(0)}$ ，开始逐次迭代。
与 Jacobi 迭代法相似，只是将 Jacobi 迭代法的迭代公式改为：
$$
{c}{x_{i}^{(k+1)}=\left(b_{i}-\sum_{j=1}^{i-1} a_{i j} x_{j}^{(k+1)}-\sum_{j=i+1}^{n} a_{i j} x_{j}^{(k)}\right) / a_{i i}} \ {(i=1,2, \cdots, n ; k=0,1,2, \cdots,n)}
$$
这里假定
设已求得
-
用 Gauss-Seidel(GS) 迭代求得 $$ x_{i}^{\prime}=\frac{1}{a_{i i}}\left(-\sum_{j=1}^{i-1} a_{i j} x_{j}^{(k)}-\sum_{j=i+1}^{n} a_{i j} x_{j}^{(k-1)}+b_{i}\right) $$
-
计算 $x_i' $ 与第
$k-1$ 次迭代值$x_i^{(k-1)}$ 的加权平均$\omega$ 作为第$k$ 次迭代值 $$ x_{i}^{(k)}=(1-\omega) x_{i}^{(k-1)}+\omega x_{i}^{\prime} $$ 或者可以整理成 $$ x_{i}^{(k)}=x_{i}^{(k-1)}+\frac{\omega\left(-\sum_{j=1}^{i-1} a_{i j} x_{j}^{(k)}-\sum_{j=i}^{n} a_{i j} x_{j}^{(k-1)}+b_{i}\right)}{a_{i i}}, i=(1,2, \ldots, n) $$ 其中，参数$\omega$ 称为松弛因子，$0<\omega <2$。当$\omega >1$ 时，上式称为逐次超松弛迭代法；当$\omega =1$ 时，上式为Gauss-Seidel 迭代法；当$0<\omega <1$ 时，上式称为低松弛迭代法。
- CPU: Intel i7-7700HQ
- GPU: NVIDIA GeForce GTX 1050
- Microsoft Visual Studio Community 2017 15.9.10
- Nvidia Cuda 10.1.120
- Visual C++ 2017
我们使用了 C++ 实现了以上三种算法的 CPU 版本，对于不同的迭代法，我们做了关于系数矩阵维度、目标精度对其收敛速度影响的实验。由于这三种算法在性能上远弱于 GPU 版本，没有太多探讨的价值，因此这里就不展示详细的结论和相关图表，只对实验结果进行总结。
高斯-赛德尔迭代、超松驰迭代基本收敛规律与雅可比迭代相同。在目标精度相同的条件下，总体来说，系数矩阵维度越大，收敛至目标精度的迭代次数越多，与此同时，系数矩阵的数据分布等自身特点也会很大程度影响到迭代次数。
相比于雅可比迭代法，高斯-赛德尔迭代和超松驰迭代收敛速度更快，迭代次数更少。但超松驰常数
以上三种算法中，由于 Gauss-Seidel 迭代法和 SOR 迭代法中每一次迭代中都需要用到上一次迭代的结果，因此可挖掘的并行性不多（实际上它们对于雅可比法的改进思路中已经包含了对其并行性的挖掘），我们着重进行了对 Jacobi 法的迭代实现。
雅可比迭代法中可并行计算的部分包括每轮迭代中
-
算法设计：将网格设置为一维网格，其中包含若干线程块，每个线程块也设置为一维形式，其中包含若干线程，这些线程中的每个线程负责计算每轮迭代中的
$x$ 中的一个元素，即以步长为$1$ 的循环遍历方式去计算 A 中的一行乘上一个$x$ （一列）。
-
实验结果：
Fig. 3.3-1: CUDA 基本版本与 CPU 版本的加速比
从以上结果中可以看出，相对于 CPU 版本，该算法对于大小为 4096 的矩阵能够加速到 8 倍左右。
-
算法分析：
-
优点：从 CPU 迭代算法中
$x$ 中元素串行计算，变成为这些元素的并行计算，从而通过基本的并行加速减少计算时间； -
缺点：该算法中每个线程块中的线程会多次访问位于 device 全局内存的数据，即上一轮迭代的
$x$ ，$A$ 的一整行与$x$ 一整列的计算仅由一个线程负责，单个线程的负担过重，并行性不够高，效率太低。
-
优点：从 CPU 迭代算法中
-
算法设计：针对上面方法存在的问题，我们将一个线程中每次步长为
$1$ 的循环遍历计算，改为步长为线程块大小去遍历计算，在每轮循环中使用循环展开；同时为每个线程块设置一个大小与线程块大小等同的共享内存数组，然后每轮循环遍历中一个线程块中的每个线程都会并行地加载一个共享内存数组中的一个元素，从而每轮循环只需访问全局内存中的$x$ 一次，该轮循环中剩下需要访问$x$ 的时候只需从共享内存中读取。
Fig. 3.4-1: 同一线程块内的线程使用共享内存的示意图（墨绿色部分为使用共享内存的部分） -
实验结果：
Fig. 3.4-2: 共享内存版本与基本版本的加速比
Fig. 3.4-3: 共享内存版本与 CPU 版本的加速比
从图中可以看出，当维度达到
-
算法分析：
- 优点：循环展开中的语句是没有数据相关的，从而增加了并发执行的机会；将从全局内存的读取大部分替换为从读取速度更快的块内共享内存，从而减少内存数据的读取消耗。
-
缺点：下一个
$x$ 元素的计算还是仅由一个线程来完成，所以其并行性还是不够高。
-
算法设计：针对上面的问题，将网格还是设置为一维的，但其长度设置为与系数矩阵的行数等同，每个线程块的大小设置为
$32$ （与 GPU 中流处理器 SM 执行基本单位线程束的大小一致），一个线程块中的线程负责$A$ 中的一行乘上一个$x$ （一列），每个线程对A中一行里间隔$32$ 个进行一次乘积，并累加，每个线程都并行的做乘积累加到各自的终点，然后对这$32$ 个线程之间通过 Warp Shuffle 来进行规约求和。
Fig. 3.5-1: 规约思想示意图
-
实验结果：
Fig. 3.5-2: 线程束洗牌与共享内存版本的加速比
Fig. 3.5-3: 线程束洗牌版本与 CPU 版本的加速比
-
算法分析：
-
优点：一个
$x$ 元素的计算是由$32$ 个线程共同负责，同时$A$ 的一行和$x$ （上一个）的数据加载也是由这$32$ 个线程共同负责，增加了并行性；同时使用规约的方式来对$32$ 个线程的累加和进行进一步求和，减少了加法的次数；使用了 Warp Shuffle 来让线程直接读取同一线程束的其他线程的寄存器值，只有很低的延迟，同时也不消耗额外的内存资源来执行数据交换。
-
优点：一个
- 当前维度下的改进策略
- 判断最大误差的步骤并行化
- 使用动态并行策略，将是否迭代的条件判断放到核函数里，让核函数调用核函数，从而使得整个算法全在核函数内执行（难度很大）
- 系数矩阵A的维度达到几万甚至十几万更高时
- GPU的显存已经无法放下这些数据，需要采用分批计算的策略，这个时候就可以使用多GPU或者分布式GPU进行并行加速计算。
Software Engineering, Fall 2019 & Software Test, Spring 2020.
It is a collaborative and interdisciplinary project drawing on expertise from School of Software Engineering & College of Civil Engineering, Tongji University.
