private double[] lsolve(double[][] A,double[] b){
  int N=b.length;
  for (int p=0; p < N; p++) {
    int max=p;
    for (int i=p + 1; i < N; i++) {
      if (Math.abs(A[i][p]) > Math.abs(A[max][p])) {
        max=i;
      }
    }
    double[] temp=A[p];
    A[p]=A[max];
    A[max]=temp;
    double t=b[p];
    b[p]=b[max];
    b[max]=t;
    if (Math.abs(A[p][p]) <= Settings.EPS) {
      throw new RuntimeException("Matrix is singular or nearly singular");
    }
    for (int i=p + 1; i < N; i++) {
      double alpha=A[i][p] / A[p][p];
      b[i]-=alpha * b[p];
      for (int j=p; j < N; j++) {
        A[i][j]-=alpha * A[p][j];
      }
    }
  }
  double[] x=new double[N];
  for (int i=N - 1; i >= 0; i--) {
    double sum=0.0;
    for (int j=i + 1; j < N; j++) {
      sum+=A[i][j] * x[j];
    }
    x[i]=(b[i] - sum) / A[i][i];
  }
  return x;
}

private double [ ] lsolve ( double [ ] [ ] A , double [ ] b ) { int N = b . length ; for ( int p = 0 ; p < N ; p ++ ) { int max = p ; for ( int i = p + 1 ; i < N ; i ++ ) { if ( Math . abs ( A [ i ] [ p ] ) > Math . abs ( A [ max ] [ p ] ) ) { max = i ; } } double [ ] temp = A [ p ] ; A [ p ] = A [ max ] ; A [ max ] = temp ; double t = b [ p ] ; b [ p ] = b [ max ] ; b [ max ] = t ; if ( Math . abs ( A [ p ] [ p ] ) <= Settings . EPS ) { throw new RuntimeException ( "Matrix is singular or nearly singular" ) ; } for ( int i = p + 1 ; i < N ; i ++ ) { double alpha = A [ i ] [ p ] / A [ p ] [ p ] ; b [ i ] -= alpha * b [ p ] ; for ( int j = p ; j < N ; j ++ ) { A [ i ] [ j ] -= alpha * A [ p ] [ j ] ; } } } double [ ] x = new double [ N ] ; for ( int i = N - 1 ; i >= 0 ; i -- ) { double sum = 0.0 ; for ( int j = i + 1 ; j < N ; j ++ ) { sum += A [ i ] [ j ] * x [ j ] ; } x [ i ] = ( b [ i ] - sum ) / A [ i ] [ i ] ; } return x ; }
gaussian elimination with partial pivoting .
