/**
 * Solves linear equations with Gauss' method.
 * 
 * @author (Stefan Edelkamp) 
 * @version (2013)
 */
public class Gauss
{
    static final int DIM = 4;
    double M[][] = {
      {1.0,  2.0, -1.0, -2.0}, {1.0,  3.0, -1.0, -2.0},{2.0,  1.0,  1.0,  1.0}, {3.0,  1.0,  2.0,  1.0}};
    double V[] = {-6.0, -4.0, 11.0, 15.0};
    double S[] = new double[DIM];
    /**
     * Constructor for objects of class Gauss
     */
    public Gauss() {
      if (solve(M, V)) {
        for (int k = DIM-1; k>=0; k--) {
          S[k] = V[k];
          for (int i=(k+1); i<DIM; i++)  S[k] -= (M[k][i]*S[i]);
          S[k] = S[k] / M[k][k];
        } 
        System.out.println("Solution:");
        for (int i=0; i<DIM; i++) 
          System.out.print(" " + S[i]);
        System.out.println();
      }
    } 
    /**
     * Gauss method to compute S in M S = V
     * 
     * @param  M   matrix of linear equation
     * @param  V   vector of linear equation
     * @param  S   solution vector 
     * @return     flag denoting solvability status (full rank) 
     */
    public boolean solve(double M[][], double V[]) {
      for (int k=0; k<DIM-1; k++) {
        double max = M[k][k] > 0 ? M[k][k] : -M[k][k];
        int m = k;
        for (int i=k+1; i<DIM; i++) {
          double abs = M[i][k] > 0 ? M[i][k] : -M[i][k];
          if (max < abs) { max = M[i][k]; m = i; }
        }
        if (m != k) {
          for (int i=k; i<DIM; i++) { double t = M[k][i]; M[k][i] = M[m][i]; M[m][i] = t; }
          double t = V[k]; V[k] = V[m]; V[m] = t; 
        }
        if(M[k][k] == 0) return false; 
        for (int j=(k+1); j<DIM; j++) {
          double f = - M[j][k] / M[k][k];
          for (int i=k; i<DIM; i++) M[j][i] = M[j][i] + f*M[k][i];
          V[j] = V[j] + f*V[k]; 
        }
      }
      return true;
    } 
}