并行 LU 分解

实验2 LU 分解 （OMP 与 MPI 实现）

191240046 孙博文

实验原理

在线性代数中，LU 分解是矩阵分解的一种，可以将一个矩阵分解为一个下三角矩阵和一个上三角矩阵的乘积（有时是它们和一个置换矩阵的乘积）。LU分解主要应用在数值分析中，用来解线性方程、求反矩阵或计算行列式。

在 LU 分解的过程中主要的计算是利用主行 i 对其余各行 就 (j > i) 做初等行变化，各行计算之间没有数据相关关系。因此可以对矩阵按行划分来实现并行计算。考虑到计算过程中处理器负载的均衡，对矩阵采用行交叉划分；假设处理器个数为 p，矩阵的阶数为 n，则每个处理器处理的行数为 $m = \lceil n / p\rceil$。

OpenMP 是基于线程的编程模型，设计基于多线程的 OpenMP 的 LU 分解算法，其主要思想为：外层设置一个列循环，在每次循环中开设 THREAD_NUMS 个线程，每个线程处理的矩阵 A 的行为上述的 m，一次循环过后则完成对应列的变换，这样在 n 次循环过后便可完成矩阵 A 的 LU 分解。即 L 为 A 中下三角部分的元素，其对角线上元素为 1，其他为 0。U 为 A 中 的上三角部分。

对于 MPI 而言，则是将矩阵数据由 0 进程按行交叉划分发送到各个进程中。在计算过程中，主行数据所在的进程将主行的数据广播到其他各进程。

实现思路

核心伪代码：

Begin
  对所有处理器，同时执行如下算法：
  for i = 0 to m - 1 do
    for j = 0 to p - 1 do
      if (my_rank=j) then
        v = i * p + j   // 主行
        for k = v to n - 1 do
          f[k] = a[k, v]
        end for
        将 f 广播到所有处理器
      else
        v = i * p + j
        接收广播的主行到 f
      end if
      u = i * p + my_rank
      if (my_rank <= j) then
        for k = u + my_rank to n - 1 step p do
          a[k, v] = a[k, v] / f[v]
          for w = v + 1 to n - 1 do
            a[k, w] = a[k, w] - f[w] * a[k, v]
          end for
        end for
      else
        for k = u to n - 1 step p do
          a[k, v] = a[k, v] / f[v]
          for w = v + 1 to n - 1 do
            a[k, w] = a[k, w] - f[w] * a[k, v]
          end for
        end for
      end if
    end for
  end for
End     

对于 MPI 来说，仅仅是使用 Bcast 而已；而对于 OMP 编程，除了对第二层循环添加 #pragma omp parallel 以外，还要在 j 循环的末尾加上 #pragma omp barrier，保证当前主行的运算全部结束后，才能选出下一个主行。

实验结果

OMP

在 LU 分解运算开始之前用 omp_get_wtime() 记录下开始的时间，在运算完成后再用当前时间减去开始时间得到运算过程持续的时间。

输入的矩阵阶数 $n = 1000$，是一个下三角矩阵和上三角矩阵的乘积。期望 LU 分解能将其还原。

$$ A = \begin{aligned} \begin{bmatrix} 1 & 0 & 0 & ... & 0\\ 2 & 1 & 0 & ... & 0\\ 3 & 2 & 1 & ... & 0\\ ... & ... & ... & ... & ...\\ n & n - 1 & n - 2 & ... & 1\\ \end{bmatrix} \begin{bmatrix} 1 & 2 & 3 & ... & n\\ 0 & 1 & 2 & ... & n - 1\\ 0 & 0 & 1 & ... & n - 2\\ ... & ... & ... & ... & ... 0 & 0 & 0 & ... & 1\\ \end{bmatrix} \end{aligned} $$

编译时开启优化选项 -O2，得到的可执行程序在 6 核心 12 线程，3.4GHz 的平台上，得到计算过程的运行时间结果如下：

线程数	1	2	4	8	16
1	0.233030	0.150938	0.093188	0.069713	0.108846
2	0.226725	0.146629	0.097024	0.066818	0.113264
3	0.244155	0.170408	0.094506	0.069337	0.151547
4	0.259125	0.161575	0.096396	0.070211	0.186803
5	0.215795	0.149947	0.099653	0.072364	0.166726
AVG	0.235766	0.155899	0.096153	0.069689	0.145437
加速比	1.00	1.51	2.45	3.38	1.62
效率	1.00	0.75	0.61	0.42	0.10

MPI

输入与上述相同，结果为：

线程数	1	2	4	8	16
1	0.252182	0.177158	0.065233	0.066492	0.086463
2	0.240585	0.173630	0.077885	0.064615	0.086954
3	0.251120	0.156510	0.083122	0.062328	0.101553
4	0.223808	0.168837	0.077668	0.067501	0.087048
5	0.289730	0.164849	0.085826	0.066900	0.138214
AVG	0.251485	0.168197	0.077947	0.065567	0.100046
加速比	1.00	1.49	3.23	3.84	2.51
效率	1.00	0.75	0.81	0.48	0.16