42、高斯消元法的并行实现与执行时间分析

高斯消元法的并行实现与执行时间分析

1. 高斯消元法的数据存储结构

高斯消元法在并行计算中有多种数据存储方式。对于行循环实现，有一种替代的数据结构，它由一个一维指针数组和 n 个长度为 n + 1 的一维数组组成。每个长度为 n + 1 的数组包含矩阵 a 的一行和向量 b 对应的元素，指针数组的元素指向这些行数组。

这种存储方案有诸多优点，它不仅便于行的交换，还适用于分布式存储。在分布式内存中，每个处理器 Pq 存储完整的指针数组，但仅存储行索引 i 满足 i mod p = q 的行，其他指针为 NULL 指针。例如，当 n = 8 时，这种存储方式能有效组织数据，避免将 b 的元素复制到公共缓冲区，并且计算 a 和 b 的新值只需一个包含 n + 1 次迭代的循环。

2. 棋盘分布的并行实现

棋盘分布的并行实现使用参数化数据分布，其分布向量为 ((p1, b1), (p2, b2)) 。这里存在一个 p1 × p2 的虚拟处理器网格，包含 p1 行、p2 列和 p1 · p2 = p 个处理器，b1 和 b2 分别是数据块的行数和列数。通过函数 G : P → N²，每个处理器被映射到处理器网格中的唯一位置。

基于此，定义了 p1 个行组 Rq 和 p2 个列组 Cq ，它们是对所有处理器集合的划分。矩阵 A 的第 i 行分布在一个行组 Ro(i) 中，第 j 列分布在一个列组 Co(j) 中。

以一个 12×12 的矩阵（n = 12）、4 个处理器 {P1, P2, P3, P4} 和分布向量 ((2, 2), (2, 3)) 为例，虚拟处理器网格大小为 2 × 2，数据块大小为 2 × 3。此时有两个行组 R1 和 R2 ，两个列组 C1 和 C2 。矩阵 A 的分布情况表明，第 5 行分布在行组 R1 中，第 7 列分布在列组 C1 中。

使用棋盘分布计算 A(k) 的具体步骤如下：
1. 确定局部主元元素 ：由于列 k 分布在列组 Co(k) 的处理器中，只有这些处理器在其局部的列 k 元素中，确定行索引 ≥ k 且绝对值最大的元素。
2. 确定全局主元元素 ：列组 Co(k) 中的处理器在组内执行单累积操作，将各自的局部主元元素提供出来。通过最大操作，计算出全局主元元素及其行索引，并由拥有元素 a(k)kk 的根处理器将信息发送给其他处理器。
3. 交换主元行 ：主元行 r 分布在行组 Ro(r) 中，行 k 分布在行组 Ro(k) 中。若 Ro(r) = Ro(k) ，行组 Ro(k) 的处理器在其拥有的列内局部交换行 k 和行 r 的元素；若 Ro(r) ≠ Ro(k) ，行组 Ro(k) 的每个处理器将行 k 的部分元素发送给行组 Ro(r) 中属于同一列组的对应处理器。
4. 分布主元行 ：每个行组 Ro(r) 中的处理器在其列组内执行面向组的单广播操作，将主元行的部分元素发送给其他处理器。
5. 计算消元因子 ：列组 Co(k) 的处理器根据公式计算列 k 中元素的消元因子 lik 。
5a. 分布消元因子 ：列组 Co(k) 的每个处理器在其行组 Ro(i) 内执行面向组的单广播操作，将消元因子 lik 广播到该组的所有处理器。
6. 计算矩阵元素 ：每个处理器根据公式，使用其拥有的 A(k) 和 b(k) 的元素，计算 A(k + 1) 和 b(k + 1) 的元素。

3. 回代计算结果向量

为计算结果向量 x 的 n 个元素，回代过程分 n 个连续步骤进行，每个步骤包含以下计算：
1. 行组 Ro(k) 的每个处理器计算行 k 中局部元素对和式 (\sum_{j = k + 1}^{n} a_{kj}^{(n)} x_j) 的贡献。
2. 行组 Ro(k) 的处理器通过以存储元素 a(n)kk 的处理器 Pq 为根的面向组的单累积操作，确定整个和式的值。
3. 处理器 Pq 根据公式计算 xk 的值。
4. 处理器 Pq 通过单广播操作将 xk 的值发送给其他处理器。

4. 棋盘分布高斯消元法的伪代码实现

以下是用 C 语言表示的 SPMD 程序的伪代码，使用 MPI 操作实现棋盘分布的高斯消元法：

double * gauss double cyclic (double **a, double *b)
{
    double *x, *buf, *elim buf;
    int i,j,k,r,q, ql, size, buf size, elim size, psz;
    struct { double val; int pvtline; } z,y;
    MPI Status status;
    x = (double *) malloc(n * sizeof(double));
    buf = (double *) malloc((n+1) * sizeof(double));
    elim buf = (double *) malloc((n+1) * sizeof(double));
    for (k=0; k<n-1; k++) {
        if (member(me, Co(k))) {
            r = max col loc(a,k);
            z.pvtline = r; z.val = fabs(a[r][k]);
            MPI Reduce(&z,&y,1,MPI DOUBLE INT,MPI MAXLOC,
                       grp leader(Co(k)),comm(Co(k)));
        }
        MPI Bcast(&y,1,MPI DOUBLE INT,grp leader(Co(k)),MPI COMM WORLD);
        r = y.pvtline;
        if(Ro(k) == Ro(r)){
            /*pivot row and row k are in the same row group */
            if (member(me, Ro(k))) {
                if (r != k)
                    exchange row loc(a,b,r,k);
                copy row loc(a,b,k,buf); }
        }
        else /* pivot row and row k are in different row groups */
            if (member(me, Ro(k))) {
                copy row loc(a,b,k,buf);
                q = compute partner(Ro(r),me);
                psz = compute size(n,k,Ro(k));
                MPI Send(buf+k,psz,MPI DOUBLE,q,tag,MPI COMM WORLD); }
            else
                if (member(me,Ro(r))) {
                    /* executing processor contains a part of the pivot row */
                    q = compute partner(Ro(k),me);
                    psz = compute size(n,k,Ro(r));
                    MPI Recv(buf+k,psz,MPI DOUBLE,q,tag,MPI COMM WORLD,&status);
                    exchange row buf(a,b,r,buf,k);
                }
        for (q=0; q<p; q++) /* broadcast of pivot row */
            if (member(q,Ro(r)) && member(me,Cop(q))) {
                ql = rank(q,Cop(q)); buf size = compute size(n,k,Ro(k));
                MPI Bcast(buf+k,buf size,MPI DOUBLE,ql,comm(Cop(q)));}
        if ((Ro(k) != Ro(r)) && (member(me,Ro(k)))
            copy row loc(a,b,buf,k);
        if
            (member(me,Co(k))) elim buf = compute elim fact loc(a,b,k,buf);
        for (q=0; q<p; q++) /* broadcast of elimination factors */
            if (member(q,Co(k)) && member(me,Rop(q))) {
                ql = rank(q,Rop(q)); elim size = compute size(n,k,Co(k));
                MPI Bcast(elim buf,elim size,MPI DOUBLE,ql,comm(Rop(q))); }
        compute local entries(a,b,k,elim buf,buf); }
    backward substitution(a,b,x);
    return x;
}

这个伪代码中的计算与行循环分布的伪代码类似，但增加了一些用于组织处理器组计算的函数，如 Co(k) 和 Ro(k) 分别表示拥有列 k 或行 k 的列组或行组，member(me, G) 用于判断处理器 me 是否属于组 G 等。

5. 并行执行时间的分析

分析高斯消元法的并行执行时间，需使用与并行机器特性相关的函数来表示计算和通信时间。对于棋盘分布的高斯消元法，忽略主元和回代过程，主要考虑第 4、5、5a 和 6 阶段。由于每个阶段后有屏障同步，可分别考虑这四个 SPMD 计算阶段。

通信阶段使用描述集体通信操作时间的公式，以单广播操作为例，p 个处理器和消息大小为 m 的单广播通信时间记为 Tsb(p, m) 。假设不相交处理器集上的独立通信操作可并行进行，p 和 m 的值取决于数据分布、行组和列组的大小以及高斯消元法的步骤 k ，随着 k 的增加，消息会变小。

计算阶段的时间也与步骤数 k 有关，随着 k 的增加，需要计算的消元因子或矩阵元素减少，算术运算次数也随之减少。每个算术运算的时间记为 top ，由于处理器执行 SPMD 程序，执行算术运算最多的处理器决定了该阶段的计算时间。

下面是各阶段通信和计算时间的具体分析：
| 阶段 | 描述 | 时间计算 |
| ---- | ---- | ---- |
| 4 | 广播主元行 k | 每个处理器 q ∈ Ro(k) 发送行 k 中列索引 ≥ k 的元素，通信时间为 (\max_{q \in Ro(k)} Tsb(#Cop(q), N_{row\geq k}^q)) ，其中 (N_{row\geq k}^q := #{(k, j) \in I_q | j \geq k}) |
| 5 | 计算消元因子 | 列组 Co(k) 中的处理器计算列 k 元素的消元因子，计算时间为 (\max_{q \in Co(k)} N_{col> k}^q \cdot top) ，其中 (N_{col> k}^q := #{(i, k) \in I_q | i > k}) |
| 5a | 分布消元因子 | 列组 Co(k) 中的处理器将消元因子发送到其行组 Rop(q) ，通信时间为 (\max_{q \in Co(k)} Tsb(#Rop(q), N_{col> k}^q)) |
| 6 | 计算矩阵元素 | 每个处理器重新计算其拥有的矩阵元素，计算时间为 (\max_{q \in P} N_{col> k}^q \cdot N_{row> k}^q \cdot 2top) |

总的并行执行时间为：
[T(n, p) = \sum_{k = 1}^{n - 1} \left(\max_{q \in Ro(k)} Tsb(#Cop(q), N_{row\geq k}^q) + \max_{q \in Co(k)} N_{col> k}^q \cdot top + \max_{q \in Co(k)} Tsb(#Rop(q), N_{col> k}^q) + \max_{q \in P} N_{col> k}^q \cdot N_{row> k}^q \cdot 2top\right)]

通过估计消息大小和算术运算次数，可将并行执行时间表示为数据分布参数 (p1, b1)、(p2, b2)、问题规模 n 和步骤数 k 的函数。考虑更大的数据块（超级块），其由 p1 × p2 个大小为 b1 × b2 的连续块组成，有 ( \lfloor \frac{n}{p1b1} \rfloor) 个超级块在行方向，( \lfloor \frac{n}{p2b2} \rfloor) 个在列方向。

处理器 q 在行 k 中列索引 ≥ k 的元素数量可估计为：
[N_{row\geq k}^q \leq \lfloor \frac{n - k + 1}{p2b2} \rfloor b2 \leq (\frac{n - k + 1}{p2b2} + 1)b2 = \frac{n - k + 1}{p2} + b2]

处理器 q 在列 k 中行索引 > k 的元素数量可估计为：
[N_{col> k}^q \leq \lfloor \frac{n - k}{p1b1} \rfloor b1 \leq (\frac{n - k}{p1b1} + 1)b1 = \frac{n - k}{p1} + b1]

由此，并行执行时间可近似为：
[T(n, p) \approx \sum_{k = 1}^{n - 1} \left(Tsb(p1, \frac{n - k + 1}{p2} + b2) + (\frac{n - k}{p1} + b1) \cdot top + Tsb(p2, \frac{n - k}{p1} + b1) + (\frac{n - k}{p1} + b1)(\frac{n - k}{p2} + b2) \cdot 2top\right)]

假设单广播操作的通信时间为 (Tsb(p, m) = \log p \cdot (\tau + m \cdot tc)) ，通过求和公式计算各阶段的通信和计算时间，最终得到总并行执行时间的表达式。由于执行时间中块大小 bi （1 ≤ bi ≤ n/pi）作为因子存在，当 b1 = b2 = 1 时可实现最小执行时间。忽略一些与 p1 和 p2 选择无关或相对较小的项后，执行时间可表示为：
[TS(p1) = \frac{n(n - 1)}{2} \left(\frac{p1 \log p1}{p} + \frac{\log p - \log p1}{p1}\right) tc + \frac{n(n - 1)}{2} \left(\frac{1}{p1} + \frac{p1}{p}\right) 2top]

通过上述分析，我们可以深入了解高斯消元法在棋盘分布下的并行实现及其执行时间特性，为优化并行计算提供理论依据。

mermaid 格式流程图如下：

graph TD
    A[开始] --> B[确定局部主元元素]
    B --> C[确定全局主元元素]
    C --> D{Ro(r) == Ro(k)?}
    D -- 是 --> E[局部交换行 k 和行 r]
    D -- 否 --> F[发送行 k 部分元素]
    E --> G[广播主元行]
    F --> G
    G --> H[计算消元因子]
    H --> I[广播消元因子]
    I --> J[计算矩阵元素]
    J --> K[回代计算结果向量]
    K --> L[结束]

综上所述，高斯消元法的并行实现通过合理的数据存储和通信策略，结合对执行时间的精确分析，能够在并行计算环境中高效地求解线性方程组。不同的存储结构和分布方式适用于不同的场景，而对执行时间的分析有助于我们选择合适的参数，进一步优化算法性能。

高斯消元法的并行实现与执行时间分析

6. 各阶段时间计算的详细推导

在前面分析各阶段通信和计算时间时，我们给出了大致的公式，下面详细推导各阶段时间的具体计算过程。

6.1 阶段 4：广播主元行

阶段 4 中，每个处理器 (q \in Ro(k)) 发送行 (k) 中列索引 (\geq k) 的元素。通信时间为 (\max_{q \in Ro(k)} Tsb(#Cop(q), N_{row\geq k}^q)) ，其中 (N_{row\geq k}^q := #{(k, j) \in I_q | j \geq k}) 。假设单广播操作的通信时间为 (Tsb(p, m) = \log p \cdot (\tau + m \cdot tc)) 。

[
\begin{align }
\sum_{k = 1}^{n - 1}Tsb(p1, \frac{n - k + 1}{p2} + b2)&=\sum_{k = 1}^{n - 1}\log p1\left((\frac{n - k + 1}{p2} + b2)tc + \tau\right)\
&=\log p1\left(\left(\frac{\sum_{k = 1}^{n - 1}(n - k + 1)}{p2}\right)tc+(n - 1)b2tc+(n - 1)\tau\right)
\end{align }
]

根据求和公式 (\sum_{k = 1}^{n - 1}(n - k + 1)=\sum_{k = 2}^{n} k=\frac{n(n + 1)}{2}-1) ，可得：

[
\sum_{k = 1}^{n - 1}Tsb(p1, \frac{n - k + 1}{p2} + b2)=\log p1\left(\left(\frac{n(n + 1)}{2}-1\right)\frac{1}{p2}tc+(n - 1)b2tc+(n - 1)\tau\right)
]

6.2 阶段 5：计算消元因子

阶段 5 中，列组 (Co(k)) 中的处理器计算列 (k) 元素的消元因子，计算时间为 (\sum_{k = 1}^{n - 1}(\frac{n - k}{p1} + b1) \cdot top) 。

根据求和公式 (\sum_{k = 1}^{n - 1}(n - k)=\frac{n(n - 1)}{2}) ，可得：

[
\sum_{k = 1}^{n - 1}(\frac{n - k}{p1} + b1) \cdot top=\left(\frac{n(n - 1)}{2p1}+(n - 1)b1\right) \cdot top
]

6.3 阶段 5a：分布消元因子

阶段 5a 中，列组 (Co(k)) 中的处理器将消元因子发送到其行组 (Rop(q)) ，通信时间为 (\sum_{k = 1}^{n - 1}Tsb(p2, \frac{n - k}{p1} + b1)) 。

[
\begin{align }
\sum_{k = 1}^{n - 1}Tsb(p2, \frac{n - k}{p1} + b1)&=\sum_{k = 1}^{n - 1}\log p2\left((\frac{n - k}{p1} + b1)tc + \tau\right)\
&=\log p2\left(\left(\frac{\sum_{k = 1}^{n - 1}(n - k)}{p1}\right)tc+(n - 1)b1tc+(n - 1)\tau\right)
\end{align }
]

根据求和公式 (\sum_{k = 1}^{n - 1}(n - k)=\frac{n(n - 1)}{2}) ，可得：

[
\sum_{k = 1}^{n - 1}Tsb(p2, \frac{n - k}{p1} + b1)=\log p2\left(\left(\frac{n(n - 1)}{2}\right)\frac{1}{p1}tc+(n - 1)b1tc+(n - 1)\tau\right)
]

6.4 阶段 6：计算矩阵元素

阶段 6 中，每个处理器重新计算其拥有的矩阵元素，计算时间涉及到 (\sum_{k = 1}^{n - 1}(\frac{n - k}{p1} + b1)(\frac{n - k}{p2} + b2) \cdot 2top) 。

利用求和公式 (\sum_{k = 1}^{n - 1}\frac{n - k}{p1}\cdot\frac{n - k}{p2}=\frac{1}{p}\sum_{k = 1}^{n - 1}k^2=\frac{1}{p}\frac{n(n - 1)(2n - 1)}{6}) ，可进一步化简该阶段的计算时间。

7. 参数选择对算法性能的影响

从总并行执行时间的表达式 (TS(p1) = \frac{n(n - 1)}{2} \left(\frac{p1 \log p1}{p} + \frac{\log p - \log p1}{p1}\right) tc + \frac{n(n - 1)}{2} \left(\frac{1}{p1} + \frac{p1}{p}\right) 2top) 可以看出，参数 (p1) 、 (p2) 、 (b1) 和 (b2) 对算法性能有重要影响。

• 块大小 (b1) 和 (b2) ： 由于执行时间中块大小 (b_i) （ (1 \leq b_i \leq \frac{n}{p_i}) ）作为因子存在，当 (b1 = b2 = 1) 时可实现最小执行时间。这是因为较小的块大小可以减少不必要的数据传输和存储开销。
• 处理器网格参数 (p1) 和 (p2) ： 虽然在分析中忽略了一些与 (p1) 和 (p2) 选择无关或相对较小的项，但 (p1) 和 (p2) 的选择仍然会影响通信和计算的平衡。不同的 (p1) 和 (p2) 组合会导致不同的行组和列组划分，从而影响数据分布和通信模式。例如，当 (p1) 较大时，行组数量增多，可能会增加行之间通信的开销；当 (p2) 较大时，列组数量增多，可能会增加列之间通信的开销。

为了选择合适的参数，可以通过实验的方法，在不同的问题规模 (n) 和处理器数量 (p) 下，测试不同参数组合的执行时间，找到使算法性能最优的参数。

8. 高斯消元法并行实现的优化建议

基于前面的分析，以下是一些优化高斯消元法并行实现的建议：

• 数据预分配和缓存优化： 在算法开始前，预先分配好所需的内存空间，如 (x) 、 (buf) 和 (elim_buf) 等。同时，合理利用处理器的缓存，减少内存访问延迟。例如，将经常访问的数据块尽量存储在缓存中，提高数据读取速度。
• 通信优化： 减少不必要的通信操作，尽量使通信和计算重叠进行。例如，在处理器进行计算的同时，让其他处理器进行数据传输，提高整体效率。另外，选择合适的通信模式，如使用非阻塞通信，避免处理器在等待通信完成时处于空闲状态。
• 负载均衡： 确保各个处理器的计算负载均衡，避免出现某些处理器计算任务过重，而其他处理器闲置的情况。可以通过合理的数据分布和任务分配来实现负载均衡，例如根据处理器的性能和资源情况，动态调整每个处理器的计算任务。

9. 总结

本文详细介绍了高斯消元法的并行实现，包括行循环实现的数据存储结构和棋盘分布的并行实现方法。通过对棋盘分布的高斯消元法进行执行时间分析，我们得到了各阶段通信和计算时间的具体公式，并详细推导了各阶段时间的计算过程。

不同的数据存储结构和分布方式适用于不同的并行计算场景，而行循环和棋盘分布是两种常见且有效的方式。对并行执行时间的分析有助于我们理解算法的性能瓶颈，通过合理选择参数，如块大小 (b1) 和 (b2) 、处理器网格参数 (p1) 和 (p2) ，可以优化算法性能。

同时，我们还给出了一些优化高斯消元法并行实现的建议，包括数据预分配和缓存优化、通信优化以及负载均衡等。通过这些优化措施，可以进一步提高高斯消元法在并行计算环境中的效率，更高效地求解线性方程组。

为了更直观地展示各阶段的关系，下面给出另一个 mermaid 格式流程图：

graph LR
    A[通信时间分析] --> B[阶段 4 广播主元行]
    A --> C[阶段 5a 分布消元因子]
    D[计算时间分析] --> E[阶段 5 计算消元因子]
    D --> F[阶段 6 计算矩阵元素]
    B --> G[总并行执行时间计算]
    C --> G
    E --> G
    F --> G
    G --> H[参数选择与优化]

优化方面	具体建议
数据预分配和缓存优化	预先分配内存空间，合理利用处理器缓存
通信优化	减少不必要通信，使通信和计算重叠，使用非阻塞通信
负载均衡	合理数据分布和任务分配，确保处理器计算负载均衡

通过以上对高斯消元法并行实现和执行时间的全面分析，我们可以更好地掌握该算法在并行计算中的应用，为解决大规模线性方程组问题提供有力支持。