Reduced row echelon form
| This page uses content from Wikipedia. The original article was at Rref#Pseudocode. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance) |
You are encouraged to solve this task according to the task description, using any language you may know.
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). Built-in functions or this pseudocode (from Wikipedia) may be used:
function ToReducedRowEchelonForm(Matrix M) is
lead := 0
rowCount := the number of rows in M
columnCount := the number of columns in M
for 0 ≤ r < rowCount do
if columnCount ≤ lead then
stop
end if
i = r
while M[i, lead] = 0 do
i = i + 1
if rowCount = i then
i = r
lead = lead + 1
if columnCount = lead then
stop
end if
end if
end while
Swap rows i and r
Divide row r by M[r, lead]
for 0 ≤ i < rowCount do
if i ≠ r do
Subtract M[i, lead] multiplied by row r from row i
end if
end for
lead = lead + 1
end for
end function
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
D
Be warned that this is definitely not the most efficient implementation, as it was coded from a thought process. <lang d> import std.stdio; real[][]example = [
[1.0,2,-1,-4],
[2,3,-1,-11],
[-2,0,-3,22]
];
real getFirstNZ(real[]row) {
int i = getOrder(row);
if (i == row.length) return 0.0;
return row[i];
}
int getOrder(real[]row) {
foreach(i,ele;row) {
if (ele != 0) return i;
}
return row.length;
}
real[][]rref(real[][]input) {
real[][]ret;
// copy the input matrix
foreach(row;input) {
real[]currow = row.dup;
// fix leading digit to be positive to make things easier,
// by inverting the signs of all elements of the row if neccesary
if (getFirstNZ(currow) < 0) {
currow[] *= -1;
}
// normalize to the first element
currow[] /= getFirstNZ(currow);
ret ~= currow;
}
ret.sort.reverse;
// get the matrix into reduced row echelon form
for(int i = 0;i<ret.length;i++) {
int order = getOrder(ret[i]);
if (order == ret[i].length) {
// this row has no non-zero digits in it, so we're done,
// since any subsequent rows will also have all 0s, because of the sorting done earlier
break;
}
foreach(j,ref row;ret) {
// prevent us from erasing the current row
if (i == j) {
continue;
}
// see if we need to modify the current row
if (row[order] != 0.0) {
// cancel rows with awesome vector ops!
real[]tmp = ret[i].dup;
tmp[]*= -row[order];
row[]+=tmp[];
// ensure row is renormalized
int corder = getOrder(row);
if (corder < row.length && row[corder] != 1) {
row[]/=row[corder];
}
}
}
}
return ret;
}
void main() {
auto result = rref(example);
writefln(result);
} </lang>
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>
Mathematica
<lang Mathematica>
RowReduce[{{1, 2, -1, -4}, {2, 3, -1, -11}, {-2, 0, -3, 22}}]
</lang> gives back: <lang Mathematica>
{{1, 0, 0, -8}, {0, 1, 0, 1}, {0, 0, 1, -2}}
</lang>
OCaml
<lang ocaml>let swap_rows m i j =
let tmp = m.(i) in m.(i) <- m.(j); m.(j) <- tmp;
let rref m =
try
let lead = ref 0
and rows = Array.length m
and cols = Array.length m.(0) in
for r = 0 to pred rows do
if cols <= !lead then
raise Exit;
let i = ref r in
while m.(!i).(!lead) = 0 do
incr i;
if rows = !i then begin
i := r;
incr lead;
if cols = !lead then
raise Exit;
end
done;
swap_rows m !i r;
let lv = m.(r).(!lead) in
m.(r) <- Array.map (fun v -> v / lv) m.(r);
for i = 0 to pred rows do
if i <> r then
let lv = m.(i).(!lead) in
m.(i) <- Array.mapi (fun i iv -> iv - lv * m.(r).(i)) m.(i);
done;
incr lead;
done
with Exit -> ()
let () =
let m =
[| [| 1; 2; -1; -4 |];
[| 2; 3; -1; -11 |];
[| -2; 0; -3; 22 |]; |]
in
rref m;
Array.iter (fun row ->
Array.iter (fun v ->
Printf.printf " %d" v
) row;
print_newline()
) m</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>
R
<lang R>rref <- function(m) {
pivot <- 1
norow <- nrow(m)
nocolumn <- ncol(m)
for(r in 1:norow) {
if ( nocolumn <= pivot ) break;
i <- r
while( m[i,pivot] == 0 ) {
i <- i + 1
if ( norow == i ) {
i <- r
pivot <- pivot + 1
if ( nocolumn == pivot ) return(m)
}
}
trow <- m[i, ]
m[i, ] <- m[r, ]
m[r, ] <- trow
m[r, ] <- m[r, ] / m[r, pivot]
for(i in 1:norow) {
if ( i != r )
m[i, ] <- m[i, ] - m[r, ] * m[i, pivot]
}
pivot <- pivot + 1
}
return(m)
}
m <- matrix(c(1, 2, -1, -4,
2, 3, -1, -11,
-2, 0, -3, 22), 3, 4, byrow=TRUE)
print(m) print(rref(m))</lang>
Ruby
<lang ruby>require 'rational'
- returns an 2-D array where each element is a 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
rary
end
def convert_to_rational(ary)
new = []
ary.each_index do |row|
new << ary[row].collect {|elem| Rational(elem)}
end
new
end
- type should be one of :to_s, :to_i, :to_f, :to_r
def convert_to(ary, type)
new = []
ary.each_index do |row|
new << ary[row].collect {|elem| elem.send(type)}
end
new
end
def print_matrix(m)
max = m[0].collect {-1}
m.each {|row| row.each_index {|i| max[i] = [max[i], row[i].to_s.length].max}}
m.each {|row| row.each_index {|i| print "%#{max[i]}s " % row[i].to_s}; puts ""}
end
mtx = [
[ 1, 2, -1, -4], [ 2, 3, -1,-11], [-2, 0, -3, 22]
] print_matrix convert_to(reduced_row_echelon_form(mtx), :to_s) puts ""
mtx = [
[ 1, 2, 3, 7], [-4, 7,-2, 7], [ 3, 3, 0, 7]
] reduced = reduced_row_echelon_form(mtx) print_matrix convert_to(reduced, :to_r) print_matrix convert_to(reduced, :to_f) </lang>
1 0 0 -8 0 1 0 1 0 0 1 -2 1 0 0 2/3 0 1 0 5/3 0 0 1 1 1.0 0.0 0.0 0.666666666666667 0.0 1.0 0.0 1.66666666666667 0.0 0.0 1.0 1.0
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
TI-89 BASIC
rref([1,2,–1,–4; 2,3,–1,–11; –2,0,–3,22])
Output (in prettyprint mode):
Matrices can also be stored in variables, and entered interactively using the Data/Matrix Editor.
Ursala
The most convenient representation for a matrix in Ursala is as a list of lists. Several auxiliary functions are defined to make this task more manageable. The pivot function reorders the rows to position the first column entry with maximum magnitude in the first row. The descending function is a second order function abstracting the pattern of recursion down the major diagonal of a matrix. The reflect function allows the code for the first phase in the reduction to be reused during the upward traversal by appropriately permuting the rows and columns. The row_reduce function adds a multiple of the top row to each subsequent row so as to cancel the first column. These are all combined in the main rref function.
<lang Ursala>#import std
- import flo
pivot = -<x fleq+ abs~~bh descending = ~&a^&+ ^|ahPathS2fattS2RpC/~& reflect = ~&lxPrTSx+ *iiD ~&l-~brS+ zipp0 row_reduce = ^C/vid*hhiD *htD minus^*p/~&r times^*D/vid@bh ~&l rref = reflect+ (descending row_reduce)+ reflect+ descending row_reduce+ pivot
- show+
test =
printf/*=*'%8.4f' rref <
<1.,2.,-1.,-4.>, <2.,3.,-1.,-11.>, <-2.,0.,-3.,22.>></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
An alternative and more efficient solution is to use the msolve library function as shown, which interfaces with the lapack library if available. This solution is applicable only if the input is a non-singular augmented square matrix. <lang Ursala>
- import lin
rref = @ySzSX msolve; ^plrNCTS\~& ~&iiDlSzyCK9+ :/1.+ 0.!*t</lang>