// Big_inverse.java  find inverse of |A| into |AI|, |A|*|AI|=|I|
//                   |A| unmodified for possible reuse or testing

import java.math.BigDecimal;

public class Big_inverse
{
  public Big_inverse(final BigDecimal A[][], BigDecimal AI[][])
  {
    int n = A.length;
    int row[] = new int[n];
    int col[] = new int[n];
    BigDecimal temp[] = new BigDecimal[n];
    int hold , I_pivot , J_pivot;
    BigDecimal pivot, abs_pivot;
    BigDecimal one = new BigDecimal("1.0");
    BigDecimal zero = new BigDecimal("0.0");

    if(A[0].length!=n || AI.length!=n || AI[0].length!=n)
    {
      System.out.println("Error in inverse, inconsistent array sizes.");
    }
    one = one.setScale(330,BigDecimal.ROUND_DOWN);

    for(int i=0; i<n; i++)
      for(int j=0; j<n; j++)
	AI[i][j] = A[i][j];

    // set up row and column interchange vectors
    for(int k=0; k<n; k++)
    {
      row[k] = k ;
      col[k] = k ;
    }
    // begin main reduction loop
    for(int k=0; k<n; k++)
    {
      // find largest element for pivot
      pivot = AI[row[k]][col[k]] ;
      I_pivot = k;
      J_pivot = k;
      for(int i=k; i<n; i++)
      {
        for(int j=k; j<n; j++)
        {
          abs_pivot = pivot.abs();
          if((AI[row[i]][col[j]].abs()).compareTo(abs_pivot)>0)
          {
            I_pivot = i ;
            J_pivot = j ;
            pivot = AI[row[i]][col[j]] ;
          }
        }
      }
      if((pivot.abs()).compareTo(new BigDecimal(1.0E-300)) <0)
      {
        System.out.println("Matrix is singular !");
        return;
      }
      hold = row[k];
      row[k]= row[I_pivot];
      row[I_pivot] = hold ;
      hold = col[k];
      col[k]= col[J_pivot];
      col[J_pivot] = hold ;
       // reduce about pivot
      AI[row[k]][col[k]] = one.divide(pivot, BigDecimal.ROUND_DOWN);
      AI[row[k]][col[k]] = AI[row[k]][col[k]].
                           setScale(330,BigDecimal.ROUND_DOWN);
      for(int j=0; j<n; j++)
      {
        if(j != k)
        {
	  AI[row[k]][col[j]] = AI[row[k]][col[j]].multiply(AI[row[k]][col[k]]);
          AI[row[k]][col[j]] = AI[row[k]][col[j]].
                               setScale(330,BigDecimal.ROUND_DOWN);
        }
      }
      // inner reduction loop
      for(int i=0; i<n; i++)
      {
        if(k != i)
        {
          for(int j=0; j<n; j++)
          {
            if( k != j )
            {
	      AI[row[i]][col[j]] = AI[row[i]][col[j]].subtract(
				   AI[row[i]][col[k]].multiply(
                                   AI[row[k]][col[j]])) ;
              AI[row[i]][col[j]] = AI[row[i]][col[j]].
                                   setScale(330,BigDecimal.ROUND_DOWN);
            }
          }
          AI[row[i]][col[k]] = (AI[row[i]][col[k]].multiply(
                                AI[row[k]][col[k]])).negate();
          AI[row[i]][col[k]] = AI[row[i]][col[k]].
                               setScale(330,BigDecimal.ROUND_DOWN);
        }
      }
    }
    // end main reduction loop

    // unscramble rows
    for(int j=0; j<n; j++)
    {
      for(int i=0; i<n; i++)
      {
        temp[col[i]] = AI[row[i]][j];
      }
      for(int i=0; i<n; i++)
      {
        AI[i][j] = temp[i] ;
      }
    }
    // unscramble columns
    for(int i=0; i<n; i++)
    {
      for(int j=0; j<n; j++)
      {
        temp[row[j]] = AI[i][col[j]] ;
      }
      for(int j=0; j<n; j++)
      {
        AI[i][j] = temp[j] ;
      }
    }
  } // end Big_inverse

} // end class Big_inverse