// Big_simeq.java  solve simultaneous equations AX=Y 
// solve real linear equations for X where Y = A * X
// method: Gauss-Jordan elimination using maximum pivot
// usage:  Big_simeq(A,Y,X); or  Big_simeq(n,A,Y,X)
//    Translated to java by : Jon Squire , 26 March 2003
//    First written by Jon Squire December 1959 for IBM 650, translated to
//    other languages  e.g. Fortran converted to Ada converted to C
//    then converted to java, converted to BigDecimal

import java.math.BigDecimal;

public class Big_simeq
{
  int nbits = 332; // default bits presision, 100 digits

  Big_simeq(final BigDecimal A[][], final BigDecimal Y[], BigDecimal X[],
            int nbit)
  {
    nbits = nbit;
    new Big_simeq(A, Y, X);
  }

  Big_simeq(final BigDecimal A[][], final BigDecimal Y[], BigDecimal X[])
  {
    int n=A.length;
    BigDecimal B[][] = new BigDecimal[n][n+1];  // working matrix
    if(A[0].length!=n || Y.length!=n || X.length!=n)
    {
      System.out.println("Error in Big_simeq, inconsistent array sizes.");
    }
    // build working data structure
    for(int i=0; i<n; i++)
    {
      for(int j=0; j<n; j++) B[i][j] = A[i][j];
      B[i][n] = Y[i];
    }
    solve(n, B, X);
  }

  Big_simeq(int n, final BigDecimal A[][], final BigDecimal Y[], BigDecimal X[])
  {
    BigDecimal B[][] = new BigDecimal[n][n+1];  // working matrix
    solve(n, B, X);
  }

  Big_simeq(int n, final BigDecimal AA[], final BigDecimal Y[], BigDecimal X[])
  {
    BigDecimal B[][] = new BigDecimal[n][n+1];  // working matrix

    for(int i=0; i<n; i++)
    {
      for(int j=0; j<n; j++) B[i][j] = AA[i*n+j];
      B[i][n] = Y[i];
    }
    solve(n, B, X);
  }

  Big_simeq(int n, final BigDecimal B[][], BigDecimal X[])
  {
    solve(n, B, X);
  }


  Big_simeq(int n, final double A[][], double U[])
  {
    BigDecimal B[][] = new BigDecimal[n][n+1];  // working matrix
    BigDecimal X[] = new BigDecimal[n];

    // A to B
    for(int i=0; i<n; i++)
      for(int j=0; j<=n; j++)
      {
	  B[i][j] = new BigDecimal(A[i][j]);
          B[i][j] = B[i][j].setScale(nbits, BigDecimal.ROUND_DOWN);
      }
    solve(n, B, X);
    // X to U
    for(int i=0; i<n; i++)
      U[i] = new Double(X[i].toString());
  }

  private void solve(final int n, final BigDecimal B[][], BigDecimal X[])
  {
    int hold , I_pivot;             // pivot indicies
    BigDecimal pivot;                   // pivot element value
    BigDecimal abs_pivot;
    int row[] = new int[n];           // row interchange indicies

    // set up row interchange vectors
    for(int k=0; k<n; k++)
    {
      row[k] = k;
    }
    //  begin main reduction loop
    for(int k=0; k<n; k++)
    {
      // find largest element for pivot
      pivot = B[row[k]][k];
      abs_pivot = pivot.abs();
      I_pivot = k;
      for(int i=k; i<n; i++)
      {
	if((B[row[i]][k].abs()).compareTo(abs_pivot)>0)
        {
          I_pivot = i;
          pivot = B[row[i]][k];
          abs_pivot = pivot.abs();
        }
      }
      // have pivot, interchange row indicies
      hold = row[k];
      row[k] = row[I_pivot];
      row[I_pivot] = hold;
      // check for near singular
      if(abs_pivot.compareTo(new BigDecimal(1.0E-99))<0)
      {
        for(int j=k+1; j<n+1; j++)
        {
	    B[row[k]][j] = new BigDecimal(0.0);
        }
        System.out.println("redundant row (singular) "+row[k]);
      } // singular, delete row
      else
      {
        // reduce about pivot
        for(int j=k+1; j<n+1; j++)
        {
	  B[row[k]][j] = B[row[k]][j].divide(B[row[k]][k],BigDecimal.ROUND_DOWN);
	  B[row[k]][j] = B[row[k]][j].setScale(nbits,BigDecimal.ROUND_DOWN);
        }
        //  inner reduction loop
        for(int i=0; i<n; i++)
        {
          if( i != k)
          {
            for(int j=k+1; j<n+1; j++)
            {
	      B[row[i]][j] = B[row[i]][j].subtract(
                             B[row[i]][k].multiply(B[row[k]][j]));
              B[row[i]][j] = B[row[i]][j].setScale(nbits,BigDecimal.ROUND_DOWN);
            }
          }
        }
      }
      //  finished inner reduction
    }
    //  end main reduction loop
    //  build  X  for return, unscrambling rows
    for(int i=0; i<n; i++)
    {
      X[i] = B[row[i]][n].setScale(nbits,BigDecimal.ROUND_DOWN);
    }
  } // end solve
} // end class Big_simeq