# include <cstdlib>
# include <cstring>
# include <ctime>
# include <iomanip>
# include <iostream>

using namespace std;

# include "rref2.hpp"

int main ( );
void is_rref_test ( );
void rref_compute_test ( );
void rref_columns_test ( );
void rref_determinant_test ( );
void rref_inverse_test ( );
void rref_rank_test ( );
void rref_solve_test ( );
void i4vec_print ( int n, int a[], string title );
double *r8mat_mm_new ( int n1, int n2, int n3, double a[], double b[] );
void r8mat_print ( int m, int n, double a[], string title );
void r8mat_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi,
  int jhi, string title );
void timestamp ( );

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

int main ( )

//****************************************************************************80
//
//  Purpose:
//
//    rref2_test() tests rref2().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    15 December 2024
//
//  Author:
//
//    John Burkardt
//
{
  timestamp ( );
  cout << "\n";
  cout << "rref2_test():\n";
  cout << "  C++ version\n";
  cout << "  Test rref(), which analyzes matrices using the\n";
  cout << "  reduced row echelon form (RREF)\n";

  is_rref_test ( );
  rref_columns_test ( );
  rref_compute_test ( );
  rref_determinant_test ( );
  rref_inverse_test ( );
  rref_rank_test ( );
  rref_solve_test ( );
//
//  Terminate.
//
  cout << "\n";
  cout << "rref2_test():\n";
  cout << "  Normal end of execution.\n";
  cout << "\n";
  timestamp ( );

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

void is_rref_test ( )

//****************************************************************************80
//
//  Purpose:
//
//    is_rref_test() tests is_rref().
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    15 December 2024
//
//  Author:
//
//    John Burkardt
//
{
  int m = 4;
  int n = 5;
//
//  Zero rows must come last.
//
  double A0[4*5] = {
    1, 0, 0, 0, 
    0, 0, 0, 0, 
    0, 1, 0, 0, 
    9, 0, 0, 1, 
    4, 8, 0, 0  };
//
//  First nonzero must be to right of first nonzero in previous row.
//
  double A1[4*5] = {
    1, 0, 0, 0, 
    0, 0, 0, 0, 
    0, 0, 1, 0, 
    9, 1, 0, 1, 
    4, 0, 8, 0  };
//
//  First nonzero must be a 1.
//
  double A2[4*5] = {
    1, 0, 0, 0, 
    0, 1, 0, 0, 
    0, 0, 3, 0, 
    9, 2, 0, 0, 
    4, 8, 0, 0  };
//
//  First nonzero must only nonzero in its column.
//
  double A3[4*5] = {
    1, 0, 0, 0, 
    0, 1, 0, 0, 
    3, 0, 1, 0, 
    9, 2, 0, 0, 
    4, 8, 0, 0 };
//
//  RREF example.
//
  double A4[4*5] = {
    1, 0, 0, 0, 
    0, 1, 0, 0, 
    3, 2, 0, 0, 
    0, 0, 1, 0, 
    4, 8, 0, 0 };

  cout << "\n";
  cout << "is_rref_test():\n";
  cout << "  is_rref() reports if a matrix is in reduced row echelon format.\n";

  r8mat_print ( m, n, A0, "  Matrix A0:" );
  cout << "  is_rref(A0) = " << is_rref ( m, n, A0 ) << "\n";

  r8mat_print ( m, n, A1, "  Matrix A1:" );
  cout << "  is_rref(A1) = " << is_rref ( m, n, A1 ) << "\n";

  r8mat_print ( m, n, A2, "  Matrix A2:" );
  cout << "  is_rref(A2) = " << is_rref ( m, n, A2 ) << "\n";

  r8mat_print ( m, n, A3, "  Matrix A3:" );
  cout << "  is_rref(A3) = " << is_rref ( m, n, A3 ) << "\n";

  r8mat_print ( m, n, A4, "  Matrix A4:" );
  cout << "  is_rref(A4) = " << is_rref ( m, n, A4 ) << "\n";

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

void rref_columns_test ( )

//****************************************************************************80
//
//  Purpose:
//
//    rref_columns_test() tests rref_columns().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    13 December 2024
//
//  Author:
//
//    John Burkardt
//
{
  double A[7*4] = {
    1,  2,  3,  4,  5,  6,  7, 
    2,  4,  6,  8, 10, 12, 14, 
    3,  9,  0,  0,  6,  6,  2, 
    1,  3,  0,  2,  6,  3,  1 };
  int *a_cols;
  double *A_COLUMNS;
  double *A_RREF;
  int i;
  int j;
  int j2;
  int j3;
  int m = 7;
  int n = 4;
  int n2;

  cout << "\n";
  cout << "rref_columns_test():\n";
  cout << "  rref_columns() uses the reduced row echelon form (RREF)\n";
  cout << "  of a matrix to find the linearly independent columns.\n";

  r8mat_print ( m, n, A, "  Matrix A:" );

  A_RREF = new double[m*n];
  a_cols = new int[n];
  rref_compute ( m, n, A, A_RREF, a_cols );

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

  cout << "\n";
  cout << "  Number of independent columns is " << n2 << "\n";

  i4vec_print ( n2, a_cols, "  Independent column indices" );

  A_COLUMNS = new double[m*n2];

  j2 = 0;
  for ( j = 0; j < n; j++ )
  {
    if ( a_cols[j] != -1 )
    {
      for ( i = 0; i < m; i++ )
      {
        j3 = a_cols[j];
        A_COLUMNS[i+j2*m] = A[i+j3*m];
      }
      j2 = j2 + 1;
    }
  }

  r8mat_print ( m, n2, A_COLUMNS, "  Independent columns:" );

  delete [] A_RREF;
  delete [] a_cols;
  delete [] A_COLUMNS;

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

void rref_compute_test ( )

//****************************************************************************80
//
//  Purpose:
//
//    ref_compute_test() tests rref_compute().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    13 December 2024
//
//  Author:
//
//    John Burkardt
//
{
  double A[4*7] = {
    1, -2,  3, -1, 
    3, -6,  9, -3, 
    0,  0,  0,  0, 
    2, -2,  0,  1, 
    6, -8,  6,  0, 
    3,  3,  6,  9, 
    1,  1,  2,  3 };
  int *a_cols;
  double *A_RREF;
  int m = 4;
  int n = 7;

  cout << "\n";
  cout << "rref_compute_test():\n";
  cout << "  rref_compute() is a user-written code to compute the\n";
  cout << "  reduced row echelon form (RREF) of a matrix.\n";

  r8mat_print ( m, n, A, "  Matrix A:" );

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

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

  r8mat_print ( m, n, A_RREF, "  rref_compute(A):" );

  i4vec_print ( n, a_cols, "  Column indices" );

  delete [] A_RREF;
  delete [] a_cols;

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

void rref_determinant_test ( )

//****************************************************************************80
//
//  Purpose:
//
//    rref_determinant_test() tests rref_determinant().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    13 December 2024
//
//  Author:
//
//    John Burkardt
//
{
  double A[4*4] = {
    5,  7,  6,  5, 
    7, 10,  8,  7, 
    6,  8, 10,  9, 
    5,  7,  9, 10 };
  double a_det;
  int n = 4;

  cout << "\n";
  cout << "rref_determinant_test():\n";
  cout << "  rref_determinant() uses the reduced row echelon form \n";
  cout << "  of a square matrix to compute the determinant.\n";

  r8mat_print ( n, n, A, "  matrix A:" );

  a_det = rref_determinant ( n, A );
  cout << "\n";
  cout << "  Estimated determinant of A = " << a_det << "\n";

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

void rref_inverse_test ( )

//****************************************************************************80
//
//  Purpose:
//
//    rref_inverse_test() tests rref_inverse().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    13 December 2024
//
//  Author:
//
//    John Burkardt
//
{
  double A[4*4] = {
    5,  7,  6,  5, 
    7, 10,  8,  7, 
    6,  8, 10,  9, 
    5,  7,  9, 10 };
  double *A_INV;
  int n = 4;
  double *P;

  cout << "\n";
  cout << "rref_inverse_test():\n";
  cout << "  rref_inverse() uses the reduced row echelon form\n";
  cout << "  of a square matrix to compute its inverse.\n";

  r8mat_print ( n, n, A, "  matrix A:" );

  A_INV = rref_inverse ( n, A );

  r8mat_print ( n, n, A_INV, "  Estimated inverse A_inv:" );

  P = r8mat_mm_new ( n, n, n, A_INV, A );

  r8mat_print ( n, n, P, "  Product A_inv * A:" );

  delete [] A_INV;
  delete [] P;

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

void rref_rank_test ( )

//****************************************************************************80
//
//  Purpose:
//
//    rref_rank_test() tests rref_rank().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    13 December 2024
//
//  Author:
//
//    John Burkardt
//
{
  double A[4*7] = {
    1, -2,  3, -1, 
    3, -6,  9, -3, 
    0,  0,  0,  0, 
    2, -2,  0,  1, 
    6, -8,  6,  0, 
    3,  3,  6,  9, 
    1,  1,  2,  3 };
  int a_rank;
  int m = 4;
  int n = 7;

  cout << "\n";
  cout << "rref_rank_test():\n";
  cout << "  rref_rank() uses the reduced row echelon form\n";
  cout << "  of a matrix to estimate its rank.\n";

  r8mat_print ( m, n, A, "  matrix A:" );

  a_rank = rref_rank ( m, n, A );
  cout << "\n";
  cout << "  A has rank " << a_rank << "\n";

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

void rref_solve_test ( )

//****************************************************************************80
//
//  Purpose:
//
//    rref_solve_test() tests rref_solve().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    13 December 2024
//
//  Author:
//
//    John Burkardt
//
{
  double A[4*4] = {
    5,  7,  6,  5, 
    7, 10,  8,  7, 
    6,  8, 10,  9, 
    5,  7,  9, 10 };
  double *b;
  double *b2;
  int m = 4;
  int n = 4;
  int nb = 1;
  double x[4*1] = { 1, 2, 3, 4 };
  double *x2;

  cout << "\n";
  cout << "rref_solve_test():\n";
  cout << "  rref_solve() uses the reduced row echelon form\n";
  cout << "  of a square matrix to solve a linear system.\n";

  r8mat_print ( m, n, A, "  matrix A:" );

  b = r8mat_mm_new ( m, n, nb, A, x );
  r8mat_print ( m, nb, b,  "  Right hand side b:" );

  x2 = rref_solve ( m, n, A, nb, b );

  r8mat_print ( n, nb, x2,  "  Estimated solution:" );

  b2 = r8mat_mm_new ( m, n, nb, A, x2 );
  r8mat_print ( m, nb, b2,  "  Product A * x:" );

  delete [] b;
  delete [] b2;
  delete [] x2;

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

void i4vec_print ( int n, int a[], string title )

//****************************************************************************80
//
//  Purpose:
//
//    i4vec_print() prints an I4VEC.
//
//  Discussion:
//
//    An I4VEC is a vector of I4's.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    30 November 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N, the number of components of the vector.
//
//    int A[N], the vector to be printed.
//
//    string TITLE, a title.
//
{
  int i;

  cout << "\n";
  cout << title << "\n";
  cout << "\n";
  for ( i = 0; i < n; i++ )
  {
    cout << "  " << setw(8) << i
         << ": " << setw(8) << a[i]  << "\n";
  }
  return;
}
//****************************************************************************80

double *r8mat_mm_new ( int n1, int n2, int n3, double a[], double b[] )

//****************************************************************************80
//
//  Purpose:
//
//    r8mat_mm_new() multiplies two matrices.
//
//  Discussion:
//
//    An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
//    in column-major order.
//
//    For this routine, the result is returned as the function value.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    18 October 2005
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N1, N2, N3, the order of the matrices.
//
//    double A[N1*N2], double B[N2*N3], the matrices to multiply.
//
//  Output:
//
//    double R8MAT_MM_NEW[N1*N3], the product matrix C = A * B.
//
{
  double *c;
  int i;
  int j;
  int k;

  c = new double[n1*n3];

  for ( i = 0; i < n1; i++ )
  {
    for ( j = 0; j < n3; j++ )
    {
      c[i+j*n1] = 0.0;
      for ( k = 0; k < n2; k++ )
      {
        c[i+j*n1] = c[i+j*n1] + a[i+k*n1] * b[k+j*n2];
      }
    }
  }

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

void r8mat_print ( int m, int n, double a[], string title )

//****************************************************************************80
//
//  Purpose:
//
//    r8mat_print() prints an R8MAT.
//
//  Discussion:
//
//    An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
//    in column-major order.
//
//    Entry A(I,J) is stored as A[I+J*M]
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    10 September 2009
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int M, the number of rows in A.
//
//    int N, the number of columns in A.
//
//    double A[M*N], the M by N matrix.
//
//    string TITLE, a title.
//
{
  r8mat_print_some ( m, n, a, 1, 1, m, n, title );

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

void r8mat_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi,
  int jhi, string title )

//****************************************************************************80
//
//  Purpose:
//
//    r8mat_print_some() prints some of an R8MAT.
//
//  Discussion:
//
//    An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
//    in column-major order.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    26 June 2013
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int M, the number of rows of the matrix.
//    M must be positive.
//
//    int N, the number of columns of the matrix.
//    N must be positive.
//
//    double A[M*N], the matrix.
//
//    int ILO, JLO, IHI, JHI, designate the first row and
//    column, and the last row and column to be printed.
//
//    string TITLE, a title.
//
{
# define INCX 5

  int i;
  int i2hi;
  int i2lo;
  int j;
  int j2hi;
  int j2lo;

  cout << "\n";
  cout << title << "\n";

  if ( m <= 0 || n <= 0 )
  {
    cout << "\n";
    cout << "  (None)\n";
    return;
  }
//
//  Print the columns of the matrix, in strips of 5.
//
  for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX )
  {
    j2hi = j2lo + INCX - 1;
    if ( n < j2hi )
    {
      j2hi = n;
    }
    if ( jhi < j2hi )
    {
      j2hi = jhi;
    }
    cout << "\n";
//
//  For each column J in the current range...
//
//  Write the header.
//
    cout << "  Col:    ";
    for ( j = j2lo; j <= j2hi; j++ )
    {
      cout << setw(7) << j - 1 << "       ";
    }
    cout << "\n";
    cout << "  Row\n";
    cout << "\n";
//
//  Determine the range of the rows in this strip.
//
    if ( 1 < ilo )
    {
      i2lo = ilo;
    }
    else
    {
      i2lo = 1;
    }
    if ( ihi < m )
    {
      i2hi = ihi;
    }
    else
    {
      i2hi = m;
    }

    for ( i = i2lo; i <= i2hi; i++ )
    {
//
//  Print out (up to) 5 entries in row I, that lie in the current strip.
//
      cout << setw(5) << i - 1 << ": ";
      for ( j = j2lo; j <= j2hi; j++ )
      {
        cout << setw(12) << a[i-1+(j-1)*m] << "  ";
      }
      cout << "\n";
    }
  }

  return;
# undef INCX
}
//****************************************************************************80

void timestamp ( )

//****************************************************************************80
//
//  Purpose:
//
//    timestamp() prints the current YMDHMS date as a time stamp.
//
//  Example:
//
//    31 May 2001 09:45:54 AM
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    19 March 2018
//
//  Author:
//
//    John Burkardt
//
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct std::tm *tm_ptr;
  std::time_t now;

  now = std::time ( NULL );
  tm_ptr = std::localtime ( &now );

  std::strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm_ptr );

  std::cout << time_buffer << "\n";

  return;
# undef TIME_SIZE
}