并行类研究：双共轭梯度稳定法

• 下载 VB.NET 2012 源代码 - 14.3 KB

摘要

在本文中，.NET 并行类被用于提升双共轭梯度稳定算法的性能，该算法是一种求解线性方程组的迭代方法。本文开发了一个具有适当运算符、方法和属性的 Matrix 类。开发了一个能够在独立线程上执行算法的求解器类，并对大小为 1000、1500、2000、2500、3000、4000 和 6000 的一系列矩阵进行了测试和讨论。在这项研究中，实施并行类可以将 CPU 时间平均减少 41%，最大减少 65%。

引言

.NET Framework 并行类为开发人员提供了一个简单的结构化循环，该循环可以在多个核心上并行执行一个代码块。这个类最初随 .NET Framework 4.0 引入，从那时起，它已被应用于多项研究。在本文中，并行类在一种广为人知的求解线性方程组的算法中实现，即双共轭梯度稳定法。这项工作通过拆分独立的计算并将它们在并行循环的迭代中执行来完成。对于任何这类方程，输入数据包括 A、X 和 b 的矩阵。X 和 b 的矩阵大小为 n，而 A 的矩阵大小为 n^2。因此，A 的矩阵需要大量的计算，并且在算法中出现的任何地方都可以将其视为并行计算的目标。

双共轭梯度稳定法

双共轭梯度稳定法，简称 BiConGradStab，是一种先进的迭代求解线性方程组的方法。这种方法和同系列的其他方法，如共轭梯度法，由于在计算中实现大小为 n 的向量而不是大小为 n^2 的矩阵，因此非常适合内存管理。此外，这种方法不限于特定类型或大小的矩阵，可以应用于任何稀疏矩阵 A。下面说明了 BiConGradStab 方法的算法。

1. r₀ = b − Ax₀
2. 选择一个任意向量 r̂₀ 使得 (r̂₀, r₀) ≠ 0，例如，r̂₀ = r₀
3. ρ₀ = α = ω₀ = 1
4. v₀ = p₀ = 0
5. 对于 i = 1, 2, 3, …
   1. ρ_i = (r̂₀, r_i−1)
   2. β = (ρ_i/ρ_i−1)(α/ω_i−1)
   3. p_i = r_i−1 + β(p_i−1 − ω_i−1v_i−1)
   4. v_i = Ap_i
   5. α = ρ_i/(r̂₀, v_i)
   6. s = r_i−1 − αv_i
   7. t = As
   8. ω_i = (t, s)/(t, t)
   9. x_i = x_i−1 + αp_i + ω_is
   10. 如果 x_i 足够精确则退出
   11. r_i = s − ω_it

如上所示，涉及大小为 n^2 的矩阵的步骤 4 和 7 中，矩阵可以被分割成更小的矩阵，然后在并行进程中单独计算。

矩阵类

为了执行清晰易懂的矩阵计算，有必要开发一个符合算法和计算要求的矩阵类。下面列出了所开发的矩阵类的一些属性、方法和运算符。

属性	描述
n (整型)	返回 i 元素的数量
m (整型)	返回 j 元素的数量
Values (双精度型(,))	(默认属性) 获取或设置 i 和 j 处元素的值

 

方法	描述
New ()	重载：1. void，2. 获取矩阵大小，3. 返回单位矩阵，4. 获取值。
T()	返回调用矩阵的转置矩阵
Split()	将当前矩阵拆分为相同 n 的矩阵并返回它们的数组
Join()	连接相同 n 的矩阵并返回一个矩阵
CopyTo()	将当前矩阵中的元素复制到另一个矩阵中

 

运算符	描述
-/+	减去/加上适当大小的矩阵
*	乘以适当大小的矩阵
/	将一个矩阵中 (0,0) 处的元素除以另一个矩阵中相同位置的元素，并返回一个 1*1 大小的结果矩阵。

求解器类

Solver 类包含执行 BiConGradStab 计算的方法和事件。SolverThread 在单独的 SolverThreadClass 类中执行。Solve() 方法需要矩阵 A、b、X 的初始猜测、并行进程数 nParallel 和最终迭代 Iteration。

Public Sub Solve(A As Matrix, X As Matrix, b As Matrix, Parallel As Integer, Iterations As Integer)
        A.CopyTo(Me.A)
        X.CopyTo(Me.X)
        b.CopyTo(Me.b)

        SolverThread = New SolverThreadClass
        AddHandler SolverThread.Solved, AddressOf _Solved
        Dim T As New Threading.Thread(New Threading.ParameterizedThreadStart_
                                     (AddressOf SolverThread.SolverThread))
        T.Start(New Object() {_A, _X, _b, Parallel, Iterations})
    End Sub 

在 SolverThreadClass 的 Solve() 方法中，矩阵 A 和 b 使用 Split() 方法被拆分为 nParallel 个较小的矩阵。

Dim A() As Matrix = rawA.Split(nParallel)
Dim b() As Matrix = rawB.Split(nParallel) 

在如下所示的主循环中，算法的步骤 4 和 7 在并行块中执行。

 Do While nIteration < Iterations

            NewRho = r_tilde * r
            beta = (NewRho / Rho) * (Alpha / Omega)
            P = r + beta * (P - Omega * V)

            Parallel.For(0, nParallel, Sub(i As Integer)
                                      V_Split(i) = A(i) * P
                                  End Sub)
            V.Join(V_Split)

            Alpha = NewRho / (r_tilde * V)
            s = r - Alpha * V

            Parallel.For(0, nParallel, Sub(i As Integer)
                                      t_Split(i) = A(i) * s
                                  End Sub)
            t.Join(t_Split)

            Omega = (t.T * s) / (t.T * t)
            X += Alpha * P + Omega * s
            r = s - Omega * t
            NewRho.CopyTo(Rho)

            nIteration += 1

 Loop 

Using the Code

为了在您的应用程序中使用此求解器，您需要如下声明它，以便在触发 Solved 事件时能够捕获到它。

Friend WithEvents mySolver As Solver

Solved 事件处理程序可以编写如下

Sub Solved(X As Matrix) Handles mySolver.Solved
Stop
'Any task with X here.
End Sub 

现在在调用过程中，实例化并调用了线性方程组矩阵的 solver 类。给出一个简单的例子

        mySolver = New Solver

        Dim A As New Matrix(3, 3)
        Dim x As New Matrix(1, 3)
        Dim b As New Matrix(1, 3)

        A(0, 0) = 0.1
        A(1, 1) = 0.4
        A(2, 2) = 0.2

        A(1, 0) = 0.45
        A(0, 1) = 0.67
        A(0, 2) = 0.98

        b(0, 0) = 1
        b(0, 1) = 1.47
        b(0, 2) = 1.58

        'Guesses for X
        x(0, 0) = 4
        x(0, 1) = 2
        x(0, 2) = 1

        mySolver.Solve(A, x, b, 1, 5) 'Number of Parallel=1 & Iterations = 5 

结果

CPU：2.8GHz | 内存：2GB | 精度：双精度

矩阵大小：1000*1000

矩阵大小：2000*2000

矩阵大小：4000*4000

矩阵大小：6000*6000

讨论

平均 CPU 时间图很好地反映了每个矩阵的并行类的整体性能。如图所示，随着矩阵大小的相对增加，CPU 时间逐渐减少。对于大小为 6000 的矩阵，可以将其解释为内存访问问题，因为大小为 6000 的矩阵将占用比 4000 大小矩阵多两倍的内存。此外，最佳 CPU 时间图表明，并行类对于 BiConGradStab 的性能通过增加矩阵大小而平稳提高。

结论

1. 在安装了 Windows 8.1 且条件连续的双核 2.8GHz、2GB RAM 的系统上，并行类在多次运行中将 BiConGradStab 的性能提高了 41% 到 65%，在单次运行中将 CPU 时间提高了 -35% 到 75%。在具有不同硬件和软件能力的其他系统上，这可能会有所不同。
2. 高度确定地，人们会尝试对大小为 1000 或更大的矩阵使用并行类，并期望整体性能得到显著改善。
3. 没有对整体或局部性能的精确估计，因为它受到多种因素的影响。但对于更多的并行进程，我们得到了最少的波动，并且并行提供了更可靠的结果。
4. 在最大并行进程下获得了最佳结果，其中并行类上的迭代次数与矩阵大小相当，例如，对于每边为 3000 的矩阵，有 1000 个并行进程。
5. 对于较小的矩阵，并行类不稳定但仍然有效。

建议

1. 对于更严肃的工作，最好在代码的迭代部分添加容差和/或均方误差检查器。
2. 使用预处理器，可以为 BiConGradStab 方法实现更好的收敛速度。尽管如此，本代码不旨在评估此方法的准确性或收敛速度。
3. 可以开发一个内存和 CPU 检查程序，以便用户在迭代进行之前和进行过程中了解内存和 CPU 的状态。
4. 在许多情况下，用户可能需要为了模拟目的多次执行算法。对于这种情况，如果有一个调查方法在模拟任务之前确定最佳运行结构，那将是极好的。
5. 整个工作可以通过两个独立的矩阵操作循环简单地完成，而无需涉及 BiConGradStab 算法的其他部分，这项工作旨在展示并行类的实际应用。人们可以尝试为其他迭代解决方案编写类似的代码，只要并行计算可行。