# include <cmath>
# include <cstdlib>
# include <iomanip>
# include <iostream>

using namespace std;

# include "rref2.hpp"

//****************************************************************************80

bool is_rref ( int m, int n, double *A )

//****************************************************************************80
//
//  Purpose:
//
//    is_rref() determines if a matrix is in reduced row echelon format.
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    15 December 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int m, n: the number of rows and columns in the matrix.'
//
//    double A(m,n): a matrix.
//
//  Output:
//
//    bool is_ref: True if A is in reduced row echelon form.
//
{
  int c;
  int i;
  int j;
  int p;
  int r;
  bool value;

  c = -1;

  for ( r = 0; r < m; r++ )
  {
//
//  Increment the first legal column for next pivot.
//
    c = c + 1;
//
//  Search for pivot p in this row.
//  If none, set p = n.
//
    p = n;
    for ( j = 0; j < n; j++ )
    {
      if ( A[r+j*m] != 0.0 )
      {
        p = j;
        break;
      }
    }
//
//  If p == n, go to next row.
//
    if ( p == n )
    {
      continue;
    }
//
//  If p is too early, fail.
//
    if ( p < c )
    {
      value = false;
      return value;
    }
//
//  Accept p as new c.
//
    c = p;
//
//  If A(r,c) is not 1, fail
//
    if ( A[r+c*m] != 1.0 )
    {
      value = false;
      return value;
    }
//
//  If A(r,c) is not the only nonzero in column c, fail
//
    for ( i = 0; i < m; i++ )
    {
      if ( i != r )
      {
        if ( A[i+c*m] != 0.0 )
        {
          value = false;
          return value;
        }
      }
    }
  }

  value = true;

  return value;
}
//****************************************************************************80

double r8_epsilon ( )

//****************************************************************************80
//
//  Purpose:
//
//    r8_epsilon() returns the R8 roundoff unit.
//
//  Discussion:
//
//    This is a number R which is a power of 2 with the property that,
//    to the precision of the computer's arithmetic,
//      1 < 1 + R
//    but
//      1 = ( 1 + R / 2 )
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    01 September 2012
//
//  Author:
//
//    John Burkardt
//
//  Output:
//
//    double R8_EPSILON, the R8 round-off unit.
//
{
  const double value = 2.220446049250313E-016;

  return value;
}
//****************************************************************************80

void rref_compute ( int m, int n, double *A, double *A_RREF, int *a_cols )

//****************************************************************************80
//
//  Purpose:
//
//    rref2_compute() computes the reduced row echelon form of a matrix.
//
//  Discussion:
//
//    A rectangular matrix is in row reduced echelon form if:
//
//    * The leading nonzero entry in each row has the value 1.
//
//    * All entries are zero above and below the leading nonzero entry 
//      in each row.
//
//    * The leading nonzero in each row occurs in a column to
//      the right of the leading nonzero in the previous row.
//
//    * Rows which are entirely zero occur last.
//
//  Example:
//
//    M = 4, N = 7
//
//    Matrix A:
//
//     1    3    0    2    6    3    1
//    -2   -6    0   -2   -8    3    1
//     3    9    0    0    6    6    2
//    -1   -3    0    1    0    9    3
//
//    RREF(A):
//
//     1    3    0    0    2    0    0
//     0    0    0    1    2    0    0
//     0    0    0    0    0    1   1/3
//     0    0    0    0    0    0    0
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    13 December 2024
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Charles Cullen,
//    An Introduction to Numerical Linear Algebra,
//    PWS Publishing Company, 1994,
//    ISBN: 978-0534936903,
//    LC: QA185.D37.C85.
//
//  Input:
//
//    int m, n: the number of rows and columns in the matrix A.
//
//    double A(M,N), the matrix to be analyzed. 
//
//  Output:
//
//    double A_RREF(M,N), the reduced row echelon form of the matrix.
//
//    int a_cols(n): the columns for which a pivot entry was found.
//
{
  int c;
  int i;
  int j;
  int p;
  double pval;
  int r;
  int rank;
  double temp;
  double tol;

  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < m; i++ )
    {
      A_RREF[i+j*m] = A[i+j*m];
    }
  }
//
//  Initialize rank to 0.
//
  rank = 0;
  for ( j = 0; j < n; j++ )
  {
    a_cols[j] = -1;
  }
//
//  Search for a pivot for each column, starting with column C = 0.
//
  c = 0;
//
//  Start the search for a pivot in row R = 0.
//
  r = 0;

  tol = sqrt ( r8_epsilon ( ) );

  while ( true )
  {
//
//  Find row index P of maximum element in subvector A(R:M,C).
//
    pval = -1.0;
    p = -1;
    for ( i = r; i < m; i++ )
    {
      if ( pval <= fabs ( A_RREF[i+c*m] ) )
      {
        p = i;
        pval = fabs ( A_RREF[i+c*m] );
      }
    }
//
//  If A(r:m,c) was all zero, there is no pivot for this column.
//  Here, we use a tolerance TOL instead of checking for an exact zero.
//
    if ( pval <= tol )
    {
      A_RREF[p+c*m] = 0.0;
//
//  Increment the column index for the next search.
//
      c = c + 1;
//
//  If there are no more columns to search, terminate the algorithm.
//
      if ( n <= c )
      {
        break;
      }
//
//  Restart loop with the next column index.
//
      continue;
    }
//
//  A pivot was found.  Increase rank by 1.
//
    a_cols[rank] = c;
    rank = rank + 1;
//
//  Swap rows R and P.
//
    for ( j = c; j < n; j++ )
    {
      temp          = A_RREF[r+j*m];
      A_RREF[r+j*m] = A_RREF[p+j*m];
      A_RREF[p+j*m] = temp;
    }
//
//  Normalize row R so that A(r,c) = 1.
//
    for ( j = c + 1; j < n; j++ )
    {
      A_RREF[r+j*m] = A_RREF[r+j*m] / A_RREF[r+c*m];
    }
    A_RREF[r+c*m] = 1.0;
//
//  Eliminate nonzeros in column C using multiples of row R, except for A(R,C).
//
    for ( i = 0; i < m; i++ )
    {
      if ( i != r )
      {
        temp = A_RREF[i+c*m];
        for ( j = c; j < n; j++ )
        {
          A_RREF[i+j*m] = A_RREF[i+j*m] - temp * A_RREF[r+j*m];
        }
      }
    }
//
//  Increment C, but terminate if we reached the value C=N.
//
    if ( c + 1 < n )
    {
      c = c + 1;
    }
    else
    {
      break;
    }
//
//  Increment R, but terminate if we reached the value R=M.
//
    if ( r + 1 < m )
    {
      r = r + 1;
    }
    else
    {
      break;
    }

  }

  return;
}
//****************************************************************************80

double rref_determinant ( int n, double *A )

//****************************************************************************80
//
//  Purpose:
//
//    rref_determinant() uses reduced row echelon form (RREF) to compute a determinant.
//
//  Discussion:
//
//    The procedure will fail if A is not square.
//
//    This is simply a demonstration of how the RREF can be used to compute
//    the determinant.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    13 December 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int n: the number of rows and columns in the matrix A.
//
//    double A(N,N), a square matrix.
//
//  Output:
//
//    double rref_determinant: the determinant.
//
{
  int *a_cols;
  double a_det;
  double *A_RREF;
  int i;
//
//  Do an RREF on A.
//
  A_RREF = new double[n*n];
  a_cols = new int[n];
  rref_compute ( n, n, A, A_RREF, a_cols );
//
//  Compute product of diagonal entries.
//
  a_det = 1.0;
  for ( i = 0; i < n; i++ )
  {
    a_det = a_det * A_RREF[i+i*n];
  }

  delete [] A_RREF;
  delete [] a_cols;

  return a_det;
}
//****************************************************************************80

double *rref_inverse ( int n, double *A )

//****************************************************************************80
//
//  Purpose:
//
//    rref_inverse() uses reduced row echelon form (RREF) to compute an inverse.
//
//  Discussion:
//
//    The procedure will fail if A is not square, or not invertible.
//
//    This is simply a demonstration of how RREF can be used to compute
//    the inverse.  But:
//    * there are better ways to compute the inverse, including
//      B = inv ( A )
//    * the inverse matrix is usually not the appropriate tool for solving
//      linear systems.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    13 December 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int n: the number of rows and columns in the matrix A.
//
//    double A(N,N), a square, invertible matrix.
//
//  Output:
//
//    double rref_inverse[N*N]: the inverse of A, computed using rref_compute().
//
{
  int *a_cols;
  double *A_INV;
  double *A_RREF;
  double *AI;
  double *AI_RREF;
  int *ai_cols;
  int i;
  int j;
//
//  First do an RREF on A alone.
//  The second argument is the independent columns found in the input matrix.
//
  A_RREF = new double[n*n];
  a_cols = new int[n];
  rref_compute ( n, n, A, A_RREF, a_cols );

  for ( i = 0; i < n; i++ )
  {
    if ( a_cols[i] == -1 )
    {
      cout << "\n";
      cout << "rref_inverse(): Warning\n";
      cout << "  The input matrix seems to be singular.\n";
      cout << "  The inverse could not be computed.\n";
      exit ( 1 );
    }
  }
//
//  If A has N independent columns, append the identity and repeat the RREF.
//
  AI = new double[n*2*n];

  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < n; i++ )
    {
      AI[i+j*n] = A[i+j*n];
      AI[i+(j+n)*n] = 0.0;
    }
    AI[j+(j+n)*n] = 1.0;
  }

  AI_RREF = new double[n*2*n];
  ai_cols = new int[2*n];
  rref_compute ( n, n + n, AI, AI_RREF, ai_cols );
//
//  The last N columns of AI_RREF are the inverse.
//
  A_INV = new double[n*n];
  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < n; i++ )
    {
      A_INV[i+j*n] = AI_RREF[i+(j+n)*n];
    }
  }

  delete [] A_RREF ;
  delete [] a_cols ;
  delete [] AI ;
  delete [] AI_RREF;
  delete [] ai_cols;

  return A_INV;
}
//****************************************************************************80

int rref_rank ( int m, int n, double *A )

//****************************************************************************80
//
//  Purpose:
//
//    rref_rank() returns the rank of a matrix, using the reduced row echelon form.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    13 December 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int m, n: the number of rows and columns in the matrix A.
//
//    double A(M,N), the matrix to be analyzed.
//
//  Output:
//
//    int rref_rank: the estimated rank of A.
//    0 <= rref_rank <= min ( M, N ).
//
{
  int *a_cols;
  int a_rank;
  double *A_RREF;
  int i;

  A_RREF = new double[m*n];
  a_cols = new int[n];

  rref_compute ( m, n, A, A_RREF, a_cols );

  a_rank = 0;
  for ( i = 0; i < n; i++ )
  {
    if ( 0 <= a_cols[i] )
    {
      a_rank = a_rank + 1;
    }
  }

  delete [] A_RREF;
  delete [] a_cols;

  return a_rank;
}
//****************************************************************************80

double *rref_solve ( int m, int n, double *A, int nb, double *b )

//****************************************************************************80
//
//  Purpose:
//
//    rref_solve() uses reduced row echelon form (RREF) to solve a linear system.
//
//  Discussion:
//
//    A linear system A*x=b is given, where
//    A is an MxN1 matrix, possibly singular.
//    b is an Mx1 right hand side, or MxN2 collection of right hand sides.
//    x is the desired N1x1 solution, or N1xN2 collection of solutions.
//
//    The right hand sides are assumed to be consistent, that is,
//    each right hand side is assumed to be a linear combination of 
//    columns of A.  If this is not so, the procedure will still produce
//    a result x, but it will not satisfy the equation A*x=b. 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    13 December 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    double A(M,N), a matrix.
//
//    double B(M,NB), one or more right hand sides consistent with A.
//
//  Output:
//
//    double rref_solve(N,NB), one or more solution vectors.
//
{
  double *AI;
  double *AI_RREF;
  int *ai_cols;
  int i;
  int j;
  double *x;

  AI = new double[m*(n+nb)];

  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < m; i++ )
    {
      AI[i+j*m] = A[i+j*m];
    }
  }
  for ( j = 0; j < nb; j++ )
  {
    for ( i = 0; i < m; i++ )
    {
      AI[i+(j+n)*m] = b[i+j*m];
    }
  }

  AI_RREF = new double[m*(n+nb)];
  ai_cols = new int[n+nb];
  rref_compute ( m, n + nb, AI, AI_RREF, ai_cols );

  x = new double[n*nb];

  for ( j = 0; j < nb; j++ )
  {
    for ( i = 0; i < m; i++ )
    {
      x[i+j*m] = AI_RREF[i+(j+n)*m];
    }
  }

  delete [] AI;
  delete [] AI_RREF;
  delete [] ai_cols;

  return x;
}