Reduced row echelon form
Show how to compute the reduced row echelon form (a.k.a. row canonical form) of a matrix. The matrix can be stored in any datatype that is convenient (for most languages, this will probably be a two-dimensional array). Pseudocode found here may be used, or built-in functions may be used.
You are encouraged to solve this task according to the task description, using any language you may know.
For testing purposes, the RREF of this matrix:
1 2 -1 -4 2 3 -1 -11 -2 0 -3 22
is:
1 0 0 -8 0 1 0 1 0 0 1 -2
ALGOL 68
<lang algol>MODE FIELD = REAL; # FIELD can be REAL, LONG REAL etc, or COMPL, FRAC etc # MODE VEC = [0]FIELD; MODE MAT = [0,0]FIELD;
PROC to reduced row echelon form = (REF MAT m)VOID: (
INT lead col := 2 LWB m;
FOR this row FROM LWB m TO UPB m DO IF lead col > 2 UPB m THEN return FI; INT other row := this row; WHILE m[other row,lead col] = 0 DO other row +:= 1; IF other row > UPB m THEN other row := this row; lead col +:= 1; IF lead col > 2 UPB m THEN return FI FI OD; IF this row /= other row THEN VEC swap = m[this row,lead col:]; m[this row,lead col:] := m[other row,lead col:]; m[other row,lead col:] := swap FI; FIELD scale = 1/m[this row,lead col]; IF scale /= 1 THEN m[this row,lead col] := 1; FOR col FROM lead col+1 TO 2 UPB m DO m[this row,col] *:= scale OD FI; FOR other row FROM LWB m TO UPB m DO IF this row /= other row THEN REAL scale = m[other row,lead col]; m[other row,lead col]:=0; FOR col FROM lead col+1 TO 2 UPB m DO m[other row,col] -:= scale*m[this row,col] OD FI OD; lead col +:= 1 OD; return: EMPTY
);
[3,4]FIELD mat := (
( 1, 2, -1, -4), ( 2, 3, -1, -11), (-2, 0, -3, 22)
);
to reduced row echelon form( mat );
FORMAT
real repr = $g(-7,4)$, vec repr = $"("n(2 UPB mat-1)(f(real repr)", ")f(real repr)")"$, mat repr = $"("n(1 UPB mat-1)(f(vec repr)", "lx)f(vec repr)")"$;
printf((mat repr, mat, $l$))</lang> Output:
(( 1.0000, 0.0000, 0.0000, -8.0000), ( 0.0000, 1.0000, 0.0000, 1.0000), ( 0.0000, 0.0000, 1.0000, -2.0000))
C++
Note: This code is written in generic form. While it slightly complicates the code, it allows to use the same code for both built-in arrays and matrix classes. To use it with a matrix class, either program the matrix class to the specifications given in the matrix_traits comment, or specialize matrix_traits for the specific interface of your matrix class.
The test code uses a built-in array for the matrix.
<lang cpp>
- include <algorithm> // for std::swap
- include <cstddef>
- include <cassert>
// Matrix traits: This describes how a matrix is accessed. By // externalizing this information into a traits class, the same code // can be used both with native arrays and matrix classes. To use the // dfault implementation of the traits class, a matrix type has to // provide the following definitions as members: // // * typedef ... index_type; // - The type used for indexing (e.g. size_t) // * typedef ... value_type; // - The element type of the matrix (e.g. double) // * index_type min_row() const; // - returns the minimal allowed row index // * index_type max_row() const; // - returns the maximal allowed row index // * index_type min_column() const; // - returns the minimal allowed column index // * index_type max_column() const; // - returns the maximal allowed column index // * value_type& operator()(index_type i, index_type k) // - returns a reference to the element i,k, where // min_row() <= i <= max_row() // min_column() <= k <= max_column() // * value_type operator()(index_type i, index_type k) const // - returns the value of element i,k // // Note that the functions are all inline and simple, so the compiler // should completely optimize them away. template<typename MatrixType> struct matrix_traits {
typedef typename MatrixType::index_type index_type; typedef typename MatrixType::value_typ value_type; static index_type min_row(MatrixType const& A) { return A.min_row(); } static index_type max_row(MatrixType const& A) { return A.max_row(); } static index_type min_column(MatrixType const& A) { return A.min_column(); } static index_type max_column(MatrixType const& A) { return A.max_column(); } static value_type& element(MatrixType& A, index_type i, index_type k) { return A(i,k); } static value_type element(MatrixType const& A, index_type i, index_type k) { return A(i,k); }
};
// specialization of the matrix traits for built-in two-dimensional // arrays template<typename T, std::size_t rows, std::size_t columns>
struct matrix_traits<T[rows][columns]>
{
typedef std::size_t index_type; typedef T value_type; static index_type min_row(T const (&)[rows][columns]) { return 0; } static index_type max_row(T const (&)[rows][columns]) { return rows-1; } static index_type min_column(T const (&)[rows][columns]) { return 0; } static index_type max_column(T const (&)[rows][columns]) { return columns-1; } static value_type& element(T (&A)[rows][columns], index_type i, index_type k) { return A[i][k]; } static value_type element(T const (&A)[rows][columns], index_type i, index_type k) { return A[i][k]; }
};
// Swap rows i and k of a matrix A // Note that due to the reference, both dimensions are preserved for // built-in arrays template<typename MatrixType>
void swap_rows(MatrixType& A, typename matrix_traits<MatrixType>::index_type i, typename matrix_traits<MatrixType>::index_type k)
{
matrix_traits<MatrixType> mt; typedef typename matrix_traits<MatrixType>::index_type index_type;
// check indices assert(mt.min_row(A) <= i); assert(i <= mt.max_row(A));
assert(mt.min_row(A) <= k); assert(k <= mt.max_row(A));
for (index_type col = mt.min_column(A); col <= mt.max_column(A); ++col) std::swap(mt.element(A, i, col), mt.element(A, k, col));
}
// divide row i of matrix A by v template<typename MatrixType>
void divide_row(MatrixType& A, typename matrix_traits<MatrixType>::index_type i, typename matrix_traits<MatrixType>::value_type v)
{
matrix_traits<MatrixType> mt; typedef typename matrix_traits<MatrixType>::index_type index_type;
assert(mt.min_row(A) <= i); assert(i <= mt.max_row(A));
assert(v != 0);
for (index_type col = mt.min_column(A); col <= mt.max_column(A); ++col) mt.element(A, i, col) /= v;
}
// in matrix A, add v times row k to row i template<typename MatrixType>
void add_multiple_row(MatrixType& A, typename matrix_traits<MatrixType>::index_type i, typename matrix_traits<MatrixType>::index_type k, typename matrix_traits<MatrixType>::value_type v)
{
matrix_traits<MatrixType> mt; typedef typename matrix_traits<MatrixType>::index_type index_type;
assert(mt.min_row(A) <= i); assert(i <= mt.max_row(A));
assert(mt.min_row(A) <= k); assert(k <= mt.max_row(A));
for (index_type col = mt.min_column(A); col <= mt.max_column(A); ++col) mt.element(A, i, col) += v * mt.element(A, k, col);
}
// convert A to reduced row echelon form template<typename MatrixType>
void to_reduced_row_echelon_form(MatrixType& A)
{
matrix_traits<MatrixType> mt; typedef typename matrix_traits<MatrixType>::index_type index_type;
index_type lead = mt.min_row(A);
for (index_type row = mt.min_row(A); row <= mt.max_row(A); ++row) { if (lead > mt.max_column(A)) return; index_type i = row; while (mt.element(A, i, lead) == 0) { ++i; if (i > mt.max_row(A)) { i = row; ++lead; if (lead > mt.max_column(A)) return; } } swap_rows(A, i, row); divide_row(A, row, mt.element(A, row, lead)); for (i = mt.min_row(A); i <= mt.max_row(A); ++i) { if (i != row) add_multiple_row(A, i, row, -mt.element(A, i, lead)); } }
}
// test code
- include <iostream>
int main() {
double M[3][4] = { { 1, 2, -1, -4 }, { 2, 3, -1, -11 }, { -2, 0, -3, 22 } };
to_reduced_row_echelon_form(M); for (int i = 0; i < 3; ++i) { for (int j = 0; j < 4; ++j) std::cout << M[i][j] << '\t'; std::cout << "\n"; }
return EXIT_SUCCESS;
} </lang> Output:
1 0 0 -8 -0 1 0 1 -0 -0 1 -2
Fortran
<lang fortran>module Rref
implicit none
contains
subroutine to_rref(matrix) real, dimension(:,:), intent(inout) :: matrix
integer :: pivot, norow, nocolumn integer :: r, i real, dimension(:), allocatable :: trow
pivot = 1 norow = size(matrix, 1) nocolumn = size(matrix, 2)
allocate(trow(nocolumn)) do r = 1, norow if ( nocolumn <= pivot ) exit i = r do while ( matrix(i, pivot) == 0 ) i = i + 1 if ( norow == i ) then i = r pivot = pivot + 1 if ( nocolumn == pivot ) return end if end do trow = matrix(i, :) matrix(i, :) = matrix(r, :) matrix(r, :) = trow matrix(r, :) = matrix(r, :) / matrix(r, pivot) do i = 1, norow if ( i /= r ) matrix(i, :) = matrix(i, :) - matrix(r, :) * matrix(i, pivot) end do pivot = pivot + 1 end do deallocate(trow) end subroutine to_rref
end module Rref</lang>
<lang fortran>program prg_test
use rref implicit none
real, dimension(3, 4) :: m = reshape( (/ 1, 2, -1, -4, & 2, 3, -1, -11, & -2, 0, -3, 22 /), & (/ 3, 4 /), order = (/ 2, 1 /) ) integer :: i
print *, "Original matrix" do i = 1, size(m,1) print *, m(i, :) end do
call to_rref(m)
print *, "Reduced row echelon form" do i = 1, size(m,1) print *, m(i, :) end do
end program prg_test</lang>
Haskell
This program was produced by translating from the Python and gradually refactoring the result into a more functional style.
<lang haskell>import Data.List (find)
rref :: Fractional a => a -> a rref m = f m 0 [0 .. rows - 1]
where rows = length m cols = length $ head m
f m _ [] = m f m lead (r : rs) | indices == Nothing = m | otherwise = f m' (lead' + 1) rs where indices = find p l p (col, row) = m !! row !! col /= 0 l = [(col, row) | col <- [lead .. cols - 1], row <- [r .. rows - 1]]
Just (lead', i) = indices newRow = map (/ m !! i !! lead') $ m !! i
m' = zipWith g [0..] $ replace r newRow $ replace i (m !! r) m g n row | n == r = row | otherwise = zipWith h newRow row where h = subtract . (* row !! lead')
replace :: Int -> a -> [a] -> [a] {- Replaces the element at the given index. -} replace n e l = a ++ e : b
where (a, _ : b) = splitAt n l</lang>
Java
<lang java>public class RREF {
public static void toRREF(double[][] M) { int rowCount = M.length; if (rowCount == 0) return;
int columnCount = M[0].length;
int lead = 0; for (int r = 0; r < rowCount; r++) { if (lead >= columnCount) break; { int i = r; while (M[i][lead] == 0) { i++; if (i == rowCount) { i = r; lead++; if (lead == columnCount) return; } } double[] temp = M[r]; M[r] = M[i]; M[i] = temp; }
{ double lv = M[r][lead]; for (int j = 0; j < columnCount; j++) M[r][j] /= lv; }
for (int i = 0; i < rowCount; i++) { if (i != r) { double lv = M[i][lead]; for (int j = 0; j < columnCount; j++) M[i][j] -= lv * M[r][j]; } } lead++; } }
public static void main(String[] args) { double[][] mtx = { { 1, 2, -1, -4}, { 2, 3, -1,-11}, {-2, 0, -3, 22}};
toRREF(mtx);
for (double [] row : mtx) { for (double x : row) System.out.print(x + ", "); System.out.println(); } }
}</lang>
Octave
<lang octave>A = [ 1, 2, -1, -4; 2, 3, -1, -11; -2, 0, -3, 22]; refA = rref(A); disp(refA);</lang>
Perl
Note that the function defined here takes an array reference, not an array. <lang perl>sub rref
{our @m; local *m = shift; @m or return; my ($lead, $rows, $cols) = (0, scalar(@m), scalar(@{$m[0]}));
foreach my $r (0 .. $rows - 1) {$lead < $cols or return; my $i = $r;
until ($m[$i][$lead]) {++$i == $rows or next; $i = $r; ++$lead == $cols and return;}
@m[$i, $r] = @m[$r, $i]; my $lv = $m[$r][$lead]; $_ /= $lv foreach @{ $m[$r] };
my @mr = @{ $m[$r] }; foreach my $i (0 .. $rows - 1) {$i == $r and next; ($lv, my $n) = ($m[$i][$lead], -1); $_ -= $lv * $mr[++$n] foreach @{ $m[$i] };}
++$lead;}}</lang>
Python
This closely follows the pseudocode given in the link. <lang python>def ToReducedRowEchelonForm( M):
if not M: return lead = 0 rowCount = len(M) columnCount = len(M[0]) for r in range(rowCount): if lead >= columnCount: return i = r while M[i][lead] == 0: i += 1 if i == rowCount: i = r lead += 1 if columnCount == lead: return M[i],M[r] = M[r],M[i] lv = M[r][lead] M[r] = [ mrx / lv for mrx in M[r]] for i in range(rowCount): if i != r: lv = M[i][lead] M[i] = [ iv - lv*rv for rv,iv in zip(M[r],M[i])] lead += 1
mtx = [
[ 1, 2, -1, -4], [ 2, 3, -1, -11], [-2, 0, -3, 22],]
ToReducedRowEchelonForm( mtx )
for rw in mtx:
print ', '.join( (str(rv) for rv in rw) )
</lang>
Ruby
<lang ruby>require 'rational'
def reduced_row_echelon_form(ary)
lead = 0 rows = ary.size cols = ary[0].size rary = convert_to_rational(ary) # use rational arithmetic catch :done do rows.times do |r| throw :done if cols <= lead i = r while rary[i][lead] == 0 i += 1 if rows == i i = r lead += 1 throw :done if cols == lead end end # swap rows i and r rary[i], rary[r] = rary[r], rary[i] # normalize row r v = rary[r][lead] rary[r].collect! {|x| x /= v} # reduce other rows rows.times do |i| next if i == r v = rary[i][lead] rary[i].each_index {|j| rary[i][j] -= v * rary[r][j]} end lead += 1 end end convert_to_string(rary)
end
def convert_to_rational(ary)
new = [] ary.each_index do |row| new << ary[row].collect {|elem| Rational(elem)} end new
end
def convert_to_string(ary)
new = [] ary.each_index do |row| new << ary[row].collect {|elem| elem.to_s} end new
end
p mtx = [
[ 1, 2, -1, -4], [ 2, 3, -1,-11], [-2, 0, -3, 22]
] p reduced_row_echelon_form(mtx)
p mtx = [
[ 1, 2, 3, 7], [-4, 7,-2, 7], [ 3, 3, 0, 7]
] p reduced_row_echelon_form(mtx)</lang>
[[1, 2, -1, -4], [2, 3, -1, -11], [-2, 0, -3, 22]] [["1", "0", "0", "-8"], ["0", "1", "0", "1"], ["0", "0", "1", "-2"]] [[1, 2, 3, 7], [-4, 7, -2, 7], [3, 3, 0, 7]] [["1", "0", "0", "2/3"], ["0", "1", "0", "5/3"], ["0", "0", "1", "1"]]
Tcl
Using utility procs defined at Matrix Transpose#Tcl <lang tcl>package require Tcl 8.5 namespace path {::tcl::mathop ::tcl::mathfunc}
proc toRREF {m} {
set lead 0 lassign [size $m] rows cols for {set r 0} {$r < $rows} {incr r} { if {$cols <= $lead} { break } set i $r while {[lindex $m $i $lead] == 0} { incr i if {$rows == $i} { set i $r incr lead if {$cols == $lead} { # Tcl can't break out of nested loops return $m } } } # swap rows i and r foreach idx [list $i $r] row [list [lindex $m $r] [lindex $m $i]] { lset m $idx $row } # divide row r by m(r,lead) set val [lindex $m $r $lead] for {set j 0} {$j < $cols} {incr j} { lset m $r $j [/ [double [lindex $m $r $j]] $val] } for {set i 0} {$i < $rows} {incr i} { if {$i != $r} { # subtract m(i,lead) multiplied by row r from row i set val [lindex $m $i $lead] for {set j 0} {$j < $cols} {incr j} { lset m $i $j [- [lindex $m $i $j] [* $val [lindex $m $r $j]]] } } } incr lead } return $m
}
set m {{1 2 -1 -4} {2 3 -1 -11} {-2 0 -3 22}} print_matrix $m print_matrix [toRREF $m]</lang> outputs
1 2 -1 -4 2 3 -1 -11 -2 0 -3 22 1.0 0.0 0.0 -8.0 -0.0 1.0 0.0 1.0 -0.0 -0.0 1.0 -2.0