//============================================================================
//
//                              Program/Module 
//                                   from
//                 Selfverifying Solvers for Dense Systems of 
//                     Linear Equations Realized in C-XSC
//
//                     Carlos Holbig and Walter Kraemer, 2003.
//
//       Program developed in C-XSC by Bernardo Frederes Kramer Alcalde,
//           Paulo Sergio Morandi Junior and Carlos Amaral Holbig
//
//============================================================================
//---------------------------------------------------------------------------
// File: lss_aprx (implementation)
// Purpose: Compute an approximate inverse of the square matrix A
// Global functions:
//   MINV(): Computes the approximate inverse of a square matrix A using the
//           Gauss-Jordan algorithm.
//---------------------------------------------------------------------------

#include <lss_aprx.hpp>

using namespace cxsc;
using namespace std;

const real eps = 1e-15;

//---------------------------------------------------------------------------
// MINV computes approximate inverse of the matrix A by use of the
// Gauss-Jordan algorithm.
// A   : input matrix and outpout of approximate inverse
// err : Error indicator = 0 everything ok, = 1 matrix singular
//---------------------------------------------------------------------------
void MINV( rmatrix& A, int& error )
{
  int       i, j, k, au1, au2, ao1, ao2;  
  bool      ok;
  real      h;
  rvector   aux( Lb(A,ROW),Ub(A,ROW) );
  intvector v( Lb(A,ROW),Ub(A,ROW) );
  int       l; 

  au1 = Lb(A,ROW); ao1 = Ub(A,ROW);
  au2 = Lb(A,COL); ao2 = Ub(A,COL);

  ok = ( (ao1 - au1) == (ao2 - au2) )? true : false; // square matrix ?

  // Vector that will control what's happening with the column
  for(i = au1; i <= ao1; i++) v[i] = i;

  i = au1-1;  // Row index
  j = au2-1;  // Column index

  while ( (j < ao2) && ok )
  {
    i++; j++;  // i and j run simultaniously

	// pivot search
	if ( j < ao2)
    {
      l = i;
      for( k = i+1; k <= ao1; k++)
        if( abs(A[k][j] ) > abs( A[l][j] ) )
          l = k;
      
	  // row interchange
	  if(i != l)
      {
        aux = A[l]; A[l] = A[i]; A[i] = aux;
        k = v[l]; v[l] = v[i]; v[i] = k;
      }
    }

    // transformation
    if( A[i][j] == 0 ) ok = false;
      else
      {
        h = 1.0 / A[i][j];
        Col(A,j) = h * rvector( Col(A,j) );
        A[i][j] = h;
        for( k = au2; k <= ao2; k++)
          if( k != j)
          {
            for( l = au1; l <= ao1; l++)
              if( l != i )
              {
                h = A[l][k] - A[l][j] * A[i][k];
                if( h == 0.0) A[l][k] *= eps;
                  else A[l][k] = h;
              }
            A[i][k] *= -A[i][j];  // A[i][k] = -A[i][k] * A[i][j]
          }
      }
  }

  // column interchange
  if(ok)
    for( k = au1; k <= ao1; k++)
    {
      l = au2 + k - au1;
      while( v[k] != k )
      {
        i = v[k]; v[k] = v[i]; v[i] = i;
        j = au2 + i - au1;
        aux = rvector( Col(A,j) ); Col(A,j) = Col(A,l); Col(A,l) = aux;
      }
    }
  error = !ok; 

} // end MINV