用于LCP的预计Gauss-Seidel

Question

我正在寻找用于解决线性互补问题的Projected Gauss-Seidel算法的C#实现.到目前为止,我已经在Bullet库中找到了用C++编写的那个,但遗憾的是它经过了高度优化(因此很难将其转换为C#).

在类似的问题中,人们建议看一下.NET的数值库.它们都只包含用于求解线性方程组的算法.

编辑:即使我找到一个,它似乎并不完整,所以问题仍然是开放的.

Answer 1

经过一周的搜索,我终于找到了这本出版物(俄文版,基于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;
}

Answer 2

你没有投影就实现了高斯赛德尔.对于投影高斯赛德尔,需要项目(或钳)上限和下限范围内的解决方案:

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矢量,并使用投影高斯赛德尔解决它:https://gist.github.com/erwincoumans/6666160