function is_rref ( m, n, A )

c*********************************************************************72
c
cc is_rref() determines if a matrix is in reduced row echelon format.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    03 April 2025
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    A(m,n): a matrix.
c
c  Output:
c
c    is_ref: True if A is in reduced row echelon form.
c
      implicit none

      integer m
      integer n

      double precision A(m,n)
      integer c
      integer i
      logical is_rref
      integer j
      integer p
      integer r

      c = 0

      do r = 1, m
c
c  Increment the first legal column for next pivot.
c
        c = c + 1
c
c  Search for pivot p in this row.
c  If none, set p = n + 1.
c
        p = n + 1
        do j = 1, n
          if ( A(r,j) /= 0.0 ) then
            p = j
            exit
          end if
        end do
c
c  If p == n + 1, go to next row.
c
        if ( p == n + 1 ) then
          c = n
          cycle
        end if
c
c  If p is too early, fail.
c
        if ( p < c ) then
          is_rref = .false.
          return
        end if
c
c  Accept p as new c.
c
        c = p
c
c  If A(r,c) is not 1, fail
c
        if ( A(r,c) /= 1.0 ) then
          is_rref = .false.
          return
        end if
c
c  If A(r,c) is not the only nonzero in column c, fail
c
        do i = 1, m
          if ( i /= r ) then
            if ( A(i,c) /= 0.0D+00 ) then
              is_rref = .false.
              return
            end if
          end if
        end do

      end do

      is_rref = .true.

      return
      end
      subroutine rref_compute ( m, n, A, A_RREF, a_cols )

c*********************************************************************72
c
cc rref_compute() computes the reduced row echelon form of a matrix.
c
c  Discussion:
c
c    A rectangular matrix is in row reduced echelon form if:
c
c    * The leading nonzero entry in each row has the value 1.
c
c    * All entries are zero above and below the leading nonzero entry
c      in each row.
c
c    * The leading nonzero in each row occurs in a column to
c      the right of the leading nonzero in the previous row.
c
c    * Rows which are entirely zero occur last.
c
c  Example:
c
c    M = 4, N = 7
c
c    Matrix A:
c
c     1    3    0    2    6    3    1
c    -2   -6    0   -2   -8    3    1
c     3    9    0    0    6    6    2
c    -1   -3    0    1    0    9    3
c
c    RREF(A):
c
c     1    3    0    0    2    0    0
c     0    0    0    1    2    0    0
c     0    0    0    0    0    1   1/3
c     0    0    0    0    0    0    0
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    03 April 2025
c
c  Author:
c
c    Original Python version by ChatGPT.
c    This version by John Burkardt.
c
c  Reference:
c
c    Charles Cullen,
c    An Introduction to Numerical Linear Algebra,
c    PWS Publishing Company, 1994,
c    ISBN: 978-0534936903,
c    LC: QA185.D37.C85.
c
c  Input:
c
c    integer m, n: the number of rows and columns in the matrix A.
c
c    real A(M,N), the matrix to be analyzed.
c
c  Output:
c
c    real A_RREF(M,N), the reduced row echelon form of the matrix.
c
c    integer a_cols(n): the columns for which a pivot entry was found.
c
      implicit none

      integer m
      integer n

      double precision A(m,n)
      integer a_cols(n)
      double precision A_RREF(m,n)
      integer c
      integer i
      integer j
      integer p
      double precision pval
      integer r
      integer rank
      double precision temp
      double precision tol
c
c  Work on a copy of the original matrix.
c
      A_RREF(1:m,1:n) = A(1:m,1:n)
c
c  Initialize rank to 0.
c
      rank = 0
c
c  Initialize indices of independent columns to 0.
c
      a_cols(1:n) = 0
c
c  Set a tolerance for the size of a pivot value.
c
      tol = sqrt ( epsilon ( tol ) )
c
c  Start the pivot search in row R = 1.
c
      r = 1
c
c  Seek a pivot for column C.
c
      do c = 1, n
c
c  Exit if we have run out of rows to examine.
c
        if ( m < r ) then
          exit
        end if
c
c  Find row index P of maximum element in subvector A(R:M,C).
c
        pval = -1.0
        p = -1
        do i = r, m
          if ( pval <= abs ( A_RREF(i,c) ) ) then
            p = i
            pval = abs ( A_RREF(i,c) )
          end if
        end do
c
c  A pivot was not found.  A(r:m,c) was all tiny.
c
        if ( pval <= tol ) then
          A_RREF(p,c) = 0.0
          cycle
        end if
c
c  A pivot was found.
c
        rank = rank + 1
        a_cols(rank) = c
c
c  Swap rows R and P.
c
        do j = c, n
          temp        = A_RREF(r,j)
          A_RREF(r,j) = A_RREF(p,j)
          A_RREF(p,j) = temp
        end do
c
c  Normalize row R so that A(r,c) = 1.
c
        A_RREF(r,c:n) = A_RREF(r,c:n) / A_RREF(r,c)
c
c  Eliminate nonzeros in column C using multiples of row R, except for A(R,C).
c
        do i = 1, m
          if ( i /= r ) then
            A_RREF(i,c:n) = A_RREF(i,c:n) - A_RREF(i,c) * A_RREF(r,c:n)
          end if
        end do
c
c  Move to the next row.
c
        r = r + 1

      end do

      return
      end
      subroutine rref_determinant ( n, A, a_det )

c*********************************************************************72
c
cc rref_determinant() uses reduced row echelon form (RREF) to compute a determinant.
c
c  Discussion:
c
c    The procedure will fail if A is not square.
c
c    This is simply a demonstration of how the RREF can be used to compute
c    the determinant.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    03 April 2025
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer n: the number of rows and columns in the matrix A.
c
c    real A(N,N), a square matrix.
c
c  Output:
c
c    real A_DET: the determinant.
c
      implicit none

      integer n

      double precision A(n,n)
      integer a_cols(n)
      double precision a_det
      double precision A_RREF(n,n)
      integer i
c
c  Do an RREF on A.
c
      call rref_compute ( n, n, A, A_RREF, a_cols )
c
c  Compute product of diagonal entries.
c
      a_det = 1.0
      do i = 1, n
        a_det = a_det * A_RREF(i,i)
      end do

      return
      end
      subroutine rref_inverse ( n, A, A_INV )

c*********************************************************************72
c
cc rref_inverse() uses reduced row echelon form (RREF) to compute an inverse.
c
c  Discussion:
c
c    The procedure will fail if A is not square, or not invertible.
c
c    This is simply a demonstration of how RREF can be used to compute
c    the inverse.  But:
c    * there are better ways to compute the inverse, including
c      B = inv ( A )
c    * the inverse matrix is usually not the appropriate tool for solving
c      linear systems.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    03 April 2025
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer n: the number of rows and columns in the matrix A.
c
c    real A(N,N), a square, invertible matrix.
c
c  Output:
c
c    real A_INV(N,N), the inverse of A, computed using the rref_compute().
c
      implicit none

      integer n

      double precision A(n,n)
      integer a_cols(n)
      double precision A_INV(n,n)
      double precision A_RREF(n,n)
      double precision AI(n,2*n)
      double precision AI_RREF(n,2*n)
      integer ai_cols(2*n)
      integer i
      integer j
c
c  First do an RREF on A alone.
c  The second argument is the independent columns found in the input matrix.
c
      call rref_compute ( n, n, A, A_RREF, a_cols )

      do i = 1, n
        if ( a_cols(i) == 0 ) then
          write ( *, '(a)' ) ''
          write ( *, '(a)' ) 'rref_inverse(): Warning!'
          write ( *, '(a)' ) '  The input matrix seems to be singular.'
          write ( *, '(a)' ) '  The inverse could not be computed.'
          stop ( 0 )
        end if
      end do
c
c  If A has N independent columns, append the identity and repeat the RREF.
c
      AI(1:n,1:n) = A(1:n,1:n)
      do i = 1, n
        do j = n + 1, 2 * n
          AI(i,j) = 0.0
        end do
        AI(i,i+n) = 1.0
      end do

      call rref_compute ( n, n + n, AI, AI_RREF, ai_cols )
c
c  The last N columns of AI_RREF are the inverse.
c
      A_INV = AI_RREF(1:n,n+1:2*n)

      return
      end
      subroutine rref_rank ( m, n, A, a_rank )

c*********************************************************************72
c
cc rref_rank() returns the rank of a matrix, using the reduced row echelon form.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    03 April 2025
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer m, n: the number of rows and columns in the matrix A.
c
c    real A(M,N), the matrix to be analyzed.
c
c  Output:
c
c    integer A_RANK, the estimated rank of A.
c    0 <= A_RANK <= min ( M, N ).
c
      implicit none

      integer m
      integer n

      double precision A(m,n)
      integer a_cols(n)
      integer a_rank
      double precision A_RREF(m,n)
      integer i

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

      a_rank = 0
      do i = 1, n
        if ( a_cols(i) /= 0 ) then
          a_rank = a_rank + 1
        end if
      end do

      return
      end
      subroutine rref_solve ( m, n, A, nb, b, x )

c*********************************************************************72
c
cc rref_solve() uses reduced row echelon form (RREF) to solve a linear system.
c
c  Discussion:
c
c    A linear system A*x=b is given, where
c    A is an MxN1 matrix, possibly singular.
c    b is an Mx1 right hand side, or MxN2 collection of right hand sides.
c    x is the desired N1x1 solution, or N1xN2 collection of solutions.
c
c    The right hand sides are assumed to be consistent, that is,
c    each right hand side is assumed to be a linear combination of
c    columns of A.  If this is not so, the procedure will still produce
c    a result x, but it will not satisfy the equation A*x=b.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    03 April 2025
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    real A(M,N), a matrix.
c
c    real B(M,NB), one or more right hand sides consistent with A.
c
c  Output:
c
c    real X(N,NB), one or more solution vectors.
c
      implicit none

      integer m
      integer n
      integer nb

      double precision A(m,n)
      double precision AI(m,n+nb)
      double precision AI_RREF(m,n+nb)
      integer ai_cols(n+nb)
      double precision b(m,nb)
      double precision x(n,nb)

      AI(1:m,1:n) = A(1:m,1:n)
      AI(1:m,n+1:n+nb) = b(1:m,1:nb)

      call rref_compute ( m, n + nb, AI, AI_RREF, ai_cols )

      x(1:n,1:nb) = AI_RREF(1:n,n+1:n+nb)

      return
      end