LCP的投影Gauss

Question

我正在寻找投影Gauss算法的C#实现，用于求解线性互补问题。到目前为止，我已经在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;
}

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

Answer 2

经过一周的搜索，我终于找到了这的出版物(俄文，基于Kenny的作品)。描述了投影Gauss算法，并利用沙尔和终止条件对算法进行了扩展.所有这些都是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;
}