Proiettato Gauss-Seidel per LCP

Question

Sto cercando l'implementazione C# dell'algoritmo Proed Gauss-Seidel per la risoluzione dello linear complementarity problem. Finora ho trovato quello scritto in C++ nella libreria Bullet, ma sfortunatamente è altamente ottimizzato (quindi sarebbe difficile tradurlo in C#).

Nella domanda similar si è proposto di dare un'occhiata allo numerical libraries for .NET. Tutti contengono solo algoritmi per la risoluzione di systems of linear equations.

Modifica: Anche se ne ho trovato uno, non sembra essere completo, quindi la domanda è ancora aperta.

Answer 1

Hai implementato Gauss Seidel senza proiezione. Per proiettata Gauss Seidel, è necessario proiettare (o bloccare) la soluzione entro i limiti inferiore e superiore:

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

Si tratta di una piccola modifica. Ecco un sommario che mostra come estrarre il vettore A e il vettore b dalla Bullet Physics Library e risolverlo usando Gauss Seidel: https://gist.github.com/erwincoumans/6666160

Answer 2

Dopo una settimana di ricerche ho finalmente trovato la pubblicazione this (in russo, basata sul lavoro di Kenny Erleben). Un algoritmo proiettato Gauss-Seidel è descritto qui e quindi esteso con SOR e condizioni di terminazione. Tutto ciò con esempi in C++, che ho usato per questa implementazione C# di base:

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