// Simeq.scala  solve for X, linear system of eqyations [A] * |X| = |Y|
//       method : Gauss-Jordan elimination using maximum pivot for pivot
//              optional ways to pass arguments:
//   solve(int n, double A[n][n], double Y[], double X[])
//   solveb(double B[n][n+1], X[])  B is AY
//    written by : Jon Squire                          
//    Original Dec 1959 for IBM 650, translated to other languages
//    e.g. Fortran converted to Ada converted to C, Java, Ruby, Scheme, Scala
//    compile: acalac Simeq.scala  use at object level  var S = Simeq
//    in function call  X = S.solve(A, Y)
  
object Simeq
{
  def solve(A:Array[Array[Double]], Y:Array[Double]):Array[Double]=
  {
    var n = A.length
    var m = A(0).length
    var B = Array.ofDim[Double](n,n+1)
    for(i <- 0 until n)
    {
      for(j <- 0 until n)
      {
        B(i)(j) = A(i)(j)
      }
      B(i)(n) = Y(i)
    }
    return solveb(B)
  }

  def solveb(B:Array[Array[Double]]):Array[Double]=
  {
    var n = B.length
    var row = Array.ofDim[Int](n) // row INTERCHANGE INDICIES 
    var hold:Int = 0
    var i_pivot:Int = 0     // pivot INDICIES 
    var pivot:Double = 0.0  // pivot ELEMENT VALUE 
    var abs_pivot:Double = 0.0
    var X = Array.ofDim[Double](n)

    // SET UP row  INTERCHANGE VECTORS 
    for(k <- 0 until n)
    {
      row(k) = k
    }

    // BEGIN MAIN REDUCTION LOOP 
    for(k <- 0 until n)
    {
      // FIND LARGEST ELEMENT FOR pivot 
      pivot = B(row(k))(k)
      abs_pivot = pivot.abs
      i_pivot = k
      for(i <- k+1 until n)
      {
        if(B(row(i))(k).abs > abs_pivot)
        {
          i_pivot = i
          pivot = B(row(i))(k);
          abs_pivot = pivot.abs
        }
      }

      // HAVE pivot, INTERCHANGE row POINTERS 
      hold = row(k)
      row(k) = row(i_pivot)
      row(i_pivot) = hold

      // CHECK FOR NEAR SINGULAR 
      if( abs_pivot < 1.0E-64 )
      {
        for(j <- k+1 until n+1)
        {
          B(row(k))(j) = 0.0
        }
        println("redundant row (singular)"+ row(k))
      } // singular, delete row 
      else
      {
        // REDUCE ABOUT pivot 
        for(j <- k+1 until n+1)
        {
          B(row(k))(j) = B(row(k))(j) / pivot
        }

        // INNER REDUCTION LOOP 
        for(i <- 0 until n)
        {
          if( i != k)
          {
            for(j <- k+1 until n+1)
            {
              B(row(i))(j) = B(row(i))(j) - B(row(i))(k) * B(row(k))(j)
            }
          }
        }
      }
      // FINISHED INNER REDUCTION 
    }

    // END OF MAIN REDUCTION LOOP 
    // BUILD  X  FOR RETURN, UNSCRAMBLING rowS 
    for(i <- 0 until n)
    {
      X(i) = B(row(i))(n)
    }
    return X
  } // end solveb 

  def inverse(A:Array[Array[Double]]):Array[Array[Double]]=
  {
    var n = A.length
    var row = Array.ofDim[Int](n) // row INTERCHANGE INDICIES 
    var col = Array.ofDim[Int](n) // col INTERCHANGE INDICIES 
    var temp = Array.ofDim[Double](n) 
    var hold:Int = 0
    var i_pivot:Int = 0     // pivot INDICIES 
    var j_pivot:Int = 0     // pivot INDICIES 
    var pivot:Double = 0.0  // pivot ELEMENT VALUE 
    var abs_pivot:Double = 0.0
    var AI = Array.ofDim[Double](n,n)

    for(i <- 0 until n)
    {
      for(j <- 0 until n)
      {
        AI(i)(j) = A(i)(j)
      }
    }

    // set up row and column interchange vectors 
    for(k <- 0 until n)
    {
      row(k) = k 
      col(k) = k 
    }

    //               begin main reduction loop 
    for(k <- 0 until n)
    {
      //             find largest element for pivot 
      pivot = AI(row(k))(col(k)) 
      i_pivot = k
      j_pivot = k
      for(i <- k until n)
      {
        for(j <- k until n)
        {
          abs_pivot = pivot.abs 
          if( (AI(row(i))(col(j))).abs > abs_pivot )
          {
            i_pivot = i 
            j_pivot = j 
            pivot = AI(row(i))(col(j)) 
          }
        }
      }
      if(pivot.abs < 1.0E-64)
      {
        println("inverse: MATRIX is SINGULAR !!!")
        return AI
      }

      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)) = 1.0 / pivot 
      for(j <- 0 until n)
      {
        if( j != k )
        {
          AI(row(k))(col(j)) = AI(row(k))(col(j)) * AI(row(k))(col(k))
        }
      }

      //               inner reduction loop 
      for(i <- 0 until n)
      {
        if( k != i )
        {
          for (j <- 0 until n)
          {
            if( k != j )
            {
              AI(row(i))(col(j)) = AI(row(i))(col(j)) - AI(row(i))(col(k)) *
                                   AI(row(k))(col(j)) 
            }
          }
          AI(row(i))(col (k)) = - AI(row(i))(col(k)) * AI(row(k))(col(k)) 
        }
      }
    }
    //      end of main reduction loop 

    //      unscramble rows 
    for (j <- 0 until n)
    {
      for (i <- 0 until n)
      {
        temp(col(i)) = AI(row(i))(j)
      }
      for (i <- 0 until n)
      {
        AI(i)(j) = temp(i) 
      }
    }

    //      unscramble columns 
    for(i <- 0 until n)
    {
      for(j <- 0 until n)
      {
        temp(row(j)) = AI(i)(col(j)) 
      }
      for(j <- 0 until n)
      {
        AI(i)(j) = temp(j) 
      }
    }

    return AI
  } // end inverse 

} // end Simeq.scala