# include <math.h>
# include <stdbool.h>
# include <stdio.h>
# include <stdlib.h>

# include "rref2.h"

/******************************************************************************/

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

/******************************************************************************/
/*
  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;
}
/******************************************************************************/

double r8_epsilon ( )

/******************************************************************************/
/*
  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;
}
/******************************************************************************/

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

/******************************************************************************/
/*
  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:

    08 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;
}
/******************************************************************************/

double rref_determinant ( int n, double *A )

/******************************************************************************/
/*
  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:

    08 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 = ( double * ) malloc ( n * n * sizeof ( double ) );
  a_cols = ( int * ) malloc ( n * sizeof ( int ) );
  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];
  }

  free ( A_RREF );
  free ( a_cols );

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

double *rref_inverse ( int n, double *A )

/******************************************************************************/
/*
  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:

    07 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 = ( double * ) malloc ( n * n * sizeof ( double ) );
  a_cols = ( int * ) malloc ( n * sizeof ( int ) );
  rref_compute ( n, n, A, A_RREF, a_cols );

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

  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 = ( double * ) malloc ( n * 2 * n * sizeof ( double ) );
  ai_cols = ( int * ) malloc ( 2 * n * sizeof ( int ) );
  rref_compute ( n, n + n, AI, AI_RREF, ai_cols );
/*
  The last N columns of AI_RREF are the inverse.
*/
  A_INV = ( double * ) malloc ( n * n * sizeof ( double ) );
  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < n; i++ )
    {
      A_INV[i+j*n] = AI_RREF[i+(j+n)*n];
    }
  }

  free ( A_RREF );
  free ( a_cols );
  free ( AI );
  free ( AI_RREF );
  free ( ai_cols );

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

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

/******************************************************************************/
/*
  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:

    08 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 = ( double * ) malloc ( m * n * sizeof ( double ) );
  a_cols = ( int * ) malloc ( n * sizeof ( int ) );

  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;
    }
  }

  free ( A_RREF );
  free ( a_cols );

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

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

/******************************************************************************/
/*
  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:

    08 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 = ( double * ) malloc ( m * ( n + nb ) * sizeof ( double ) );

  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 = ( double * ) malloc ( m * ( n + nb ) * sizeof ( double ) );
  ai_cols = ( int * ) malloc ( ( n + nb ) * sizeof ( int ) );
  rref_compute ( m, n + nb, AI, AI_RREF, ai_cols );

  x = ( double * ) malloc ( n * nb * sizeof ( double ) );

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

  free ( AI );
  free ( AI_RREF );
  free ( ai_cols );

  return x;
}