Projected Gauss-Seidel for LCP

Question

I'm looking for the C# implementation of Projected Gauss-Seidel algorithm for solving the linear complementarity problem. So far I've found the one written in C++ in the Bullet library, but unfortunately it's highly optimized (so it would be hard to translate it into C#).

In the similar question one proposed to take a look of the numerical libraries for .NET. All of them contain only algorithms for solving systems of linear equations.

Edit: Even as I've found one, it doesn't appear to be complete, so the question is still open.

Answer 1

You implemented Gauss Seidel without projection. For projected Gauss Seidel, you need to project (or clamp) the solution within lower and upper limits:

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;
}

It is a minor modification. Here is a gist that shows how to extract the A matrix and b vector from the Bullet Physics Library and solve it using projected Gauss Seidel: https://gist.github.com/erwincoumans/6666160

Answer 2

After a week of searching I've finally found this publication (in Russian, based on Kenny Erleben's work). A projected Gauss-Seidel algorithm is described there and then extended with SOR and termination conditions. All that with examples in C++, which I used for this basic C# implementation:

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;
}