c# - LCP 的投影 Gauss-Seidel

Question

我正在寻找用于解决线性互补问题的 Projected Gauss-Seidel 算法的 C# 实现。到目前为止，我在Bullet库中找到了用 C++ 编写的那个，但不幸的是，它经过了高度优化（因此很难将其翻译成 C#）。

在类似的问题中，有人提议查看.NET 的数值库。它们都只包含求解线性方程组的算法。

编辑：即使我找到了一个，它似乎并不完整，所以问题仍然存在。

Answer 1

您在没有投影的情况下实现了 Gauss Seidel。对于投影的 Gauss Seidel，您需要在下限和上限内投影（或钳制）解：

public static double[] Solve (double[,] matrix, double[] right,
                              double relaxation, int iterations, 
                              double[] lo, double[] hi)
{
    // Validation omitted
    var x = right;
    double delta;

    // Gauss-Seidel with Successive OverRelaxation Solver
    for (int k = 0; k < iterations; ++k) {
        for (int i = 0; i < right.Length; ++i) {
            delta = 0.0f;

            for (int j = 0; j < i; ++j)
                delta += matrix [i, j] * x [j];
            for (int j = i + 1; j < right.Length; ++j)
                delta += matrix [i, j] * x [j];

            delta = (right [i] - delta) / matrix [i, i];
            x [i] += relaxation * (delta - x [i]);
    // Project the solution within the lower and higher limits
            if (x[i]<lo[i])
                x[i]=lo[i];
            if (x[i]>hi[i])
                x[i]=hi[i];
        }
    }
    return x;
}

这是一个小的修改。这是一个要点，展示了如何从 Bullet 物理库中提取 A 矩阵和 b 向量并使用投影的 Gauss Seidel 解决它：https ://gist.github.com/erwincoumans/6666160

Answer 2

经过一周的搜索，我终于找到了这本出版物（俄文，基于 Kenny Erleben 的作品）。在那里描述了一个投影的 Gauss-Seidel 算法，然后用SOR和终止条件进行了扩展。所有这些都带有 C++ 中的示例，我用于这个基本的 C# 实现：

public static double[] Solve (double[,] matrix, double[] right,
                              double relaxation, int iterations)
{
    // Validation omitted
    var x = right;
    double delta;

    // Gauss-Seidel with Successive OverRelaxation Solver
    for (int k = 0; k < iterations; ++k) {
        for (int i = 0; i < right.Length; ++i) {
            delta = 0.0f;

            for (int j = 0; j < i; ++j)
                delta += matrix [i, j] * x [j];
            for (int j = i + 1; j < right.Length; ++j)
                delta += matrix [i, j] * x [j];

            delta = (right [i] - delta) / matrix [i, i];
            x [i] += relaxation * (delta - x [i]);
        }
    }

    return x;
}