Reduced row echelon form: Difference between revisions
(→{{header|Euphoria}}: Euphoria example added) |
(Shorter D code) |
||
Line 700:
=={{header|D}}==
<lang d>import std.stdio, std.algorithm, std.array, std.conv;
const ncols =
if (ncols <= lead) return;▼
▲ foreach (r; 0 .. rowCount) {
▲ return;
int i = r;
while (M[i][lead] == 0)
i = r;
}
swap(M[i], M[r]);
M[r][] /= M[r][lead];
foreach (j; 0 ..
if (j != r)
M[j][] -= M[r][] * M[j][lead];
Line 742 ⟶ 725:
void main() {
auto
toReducedRowEchelonForm(
writeln("[", array(map!text(A)).join("\n "), "]");
}</lang>
Output:
<pre>[[1, 0, 0, -8]
[0, 1, 0, 1]
[0, 0, 1, -2]]</pre>
=={{header|Euphoria}}==
|
Revision as of 19:56, 23 July 2011
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
Ada
matrices.ads: <lang Ada>generic
type Element_Type is private; Zero : Element_Type; with function "-" (Left, Right : in Element_Type) return Element_Type is <>; with function "*" (Left, Right : in Element_Type) return Element_Type is <>; with function "/" (Left, Right : in Element_Type) return Element_Type is <>;
package Matrices is
type Matrix is array (Positive range <>, Positive range <>) of Element_Type; function Reduced_Row_Echelon_form (Source : Matrix) return Matrix;
end Matrices;</lang>
matrices.adb: <lang Ada>package body Matrices is
procedure Swap_Rows (From : in out Matrix; First, Second : in Positive) is Temporary : Element_Type; begin for Col in From'Range (2) loop Temporary := From (First, Col); From (First, Col) := From (Second, Col); From (Second, Col) := Temporary; end loop; end Swap_Rows;
procedure Divide_Row (From : in out Matrix; Row : in Positive; Divisor : in Element_Type) is begin for Col in From'Range (2) loop From (Row, Col) := From (Row, Col) / Divisor; end loop; end Divide_Row;
procedure Subtract_Rows (From : in out Matrix; Subtrahend, Minuend : in Positive; Factor : in Element_Type) is begin for Col in From'Range (2) loop From (Minuend, Col) := From (Minuend, Col) - From (Subtrahend, Col) * Factor; end loop; end Subtract_Rows;
function Reduced_Row_Echelon_form (Source : Matrix) return Matrix is Result : Matrix := Source; Lead : Positive := Result'First (2); I : Positive; begin Rows : for Row in Result'Range (1) loop exit Rows when Lead > Result'Last (2); I := Row; while Result (I, Lead) = Zero loop I := I + 1; if I = Result'Last (1) then I := Row; Lead := Lead + 1; exit Rows when Lead = Result'Last (2); end if; end loop; if I /= Row then Swap_Rows (From => Result, First => I, Second => Row); end if; Divide_Row (From => Result, Row => Row, Divisor => Result (Row, Lead)); for Other_Row in Result'Range (1) loop if Other_Row /= Row then Subtract_Rows (From => Result, Subtrahend => Row, Minuend => Other_Row, Factor => Result (Other_Row, Lead)); end if; end loop; Lead := Lead + 1; end loop Rows; return Result; end Reduced_Row_Echelon_form;
end Matrices;</lang>
Example use: main.adb: <lang Ada>with Matrices; with Ada.Text_IO; procedure Main is
package Float_IO is new Ada.Text_IO.Float_IO (Float); package Float_Matrices is new Matrices ( Element_Type => Float, Zero => 0.0); procedure Print_Matrix (Matrix : in Float_Matrices.Matrix) is begin for Row in Matrix'Range (1) loop for Col in Matrix'Range (2) loop Float_IO.Put (Matrix (Row, Col), 0, 0, 0); Ada.Text_IO.Put (' '); end loop; Ada.Text_IO.New_Line; end loop; end Print_Matrix; My_Matrix : Float_Matrices.Matrix := ((1.0, 2.0, -1.0, -4.0), (2.0, 3.0, -1.0, -11.0), (-2.0, 0.0, -3.0, 22.0)); Reduced : Float_Matrices.Matrix := Float_Matrices.Reduced_Row_Echelon_form (My_Matrix);
begin
Print_Matrix (My_Matrix); Ada.Text_IO.Put_Line ("reduced to:"); Print_Matrix (Reduced);
end Main;</lang>
Output:
1.0 2.0 -1.0 -4.0 2.0 3.0 -1.0 -11.0 -2.0 0.0 -3.0 22.0 reduced to: 1.0 0.0 0.0 -8.0 -0.0 1.0 0.0 1.0 -0.0 -0.0 1.0 -2.0
ALGOL 68
<lang algol68>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
<lang c>#include <stdio.h>
- define TALLOC(n,typ) malloc(n*sizeof(typ))
- define EL_Type int
typedef struct sMtx {
int dim_x, dim_y; EL_Type *m_stor; EL_Type **mtx;
} *Matrix, sMatrix;
typedef struct sRvec {
int dim_x; EL_Type *m_stor;
} *RowVec, sRowVec;
Matrix NewMatrix( int x_dim, int y_dim ) {
int n; Matrix m; m = TALLOC( 1, sMatrix); n = x_dim * y_dim; m->dim_x = x_dim; m->dim_y = y_dim; m->m_stor = TALLOC(n, EL_Type); m->mtx = TALLOC(m->dim_y, EL_Type *); for(n=0; n<y_dim; n++) { m->mtx[n] = m->m_stor+n*x_dim; } return m;
}
void MtxSetRow(Matrix m, int irow, EL_Type *v) {
int ix; EL_Type *mr; mr = m->mtx[irow]; for(ix=0; ix<m->dim_x; ix++) mr[ix] = v[ix];
}
Matrix InitMatrix( int x_dim, int y_dim, EL_Type **v) {
Matrix m; int iy; m = NewMatrix(x_dim, y_dim); for (iy=0; iy<y_dim; iy++) MtxSetRow(m, iy, v[iy]); return m;
}
void MtxDisplay( Matrix m ) {
int iy, ix; const char *sc; for (iy=0; iy<m->dim_y; iy++) { printf(" "); sc = " "; for (ix=0; ix<m->dim_x; ix++) { printf("%s %3d", sc, m->mtx[iy][ix]); sc = ","; } printf("\n"); } printf("\n");
}
void MtxMulAndAddRows(Matrix m, int ixrdest, int ixrsrc, EL_Type mplr) {
int ix; EL_Type *drow, *srow; drow = m->mtx[ixrdest]; srow = m->mtx[ixrsrc]; for (ix=0; ix<m->dim_x; ix++) drow[ix] += mplr * srow[ix];
// printf("Mul row %d by %d and add to row %d\n", ixrsrc, mplr, ixrdest); // MtxDisplay(m); }
void MtxSwapRows( Matrix m, int rix1, int rix2) {
EL_Type *r1, *r2, temp; int ix; if (rix1 == rix2) return; r1 = m->mtx[rix1]; r2 = m->mtx[rix2]; for (ix=0; ix<m->dim_x; ix++) temp = r1[ix]; r1[ix]=r2[ix]; r2[ix]=temp;
// printf("Swap rows %d and %d\n", rix1, rix2); // MtxDisplay(m); }
void MtxNormalizeRow( Matrix m, int rix, int lead) {
int ix; EL_Type *drow; EL_Type lv; drow = m->mtx[rix]; lv = drow[lead]; for (ix=0; ix<m->dim_x; ix++) drow[ix] /= lv;
// printf("Normalize row %d\n", rix); // MtxDisplay(m); }
- define MtxGet( m, rix, cix ) m->mtx[rix][cix]
void MtxToReducedREForm(Matrix m) {
int lead; int rix, iix; EL_Type lv; int rowCount = m->dim_y;
lead = 0; for (rix=0; rix<rowCount; rix++) { if (lead >= m->dim_x) return; iix = rix; while (0 == MtxGet(m, iix,lead)) { iix++; if (iix == rowCount) { iix = rix; lead++; if (lead == m->dim_x) return; } } MtxSwapRows(m, iix, rix ); MtxNormalizeRow(m, rix, lead ); for (iix=0; iix<rowCount; iix++) { if ( iix != rix ) { lv = MtxGet(m, iix, lead ); MtxMulAndAddRows(m,iix, rix, -lv) ; } } lead++; }
}
int main() {
Matrix m1; static EL_Type r1[] = {1,2,-1,-4}; static EL_Type r2[] = {2,3,-1,-11}; static EL_Type r3[] = {-2,0,-3,22}; static EL_Type *im[] = { r1, r2, r3 };
m1 = InitMatrix( 4,3, im ); printf("Initial\n"); MtxDisplay(m1); MtxToReducedREForm(m1); printf("Reduced R-E form\n"); MtxDisplay(m1); return 0;
}</lang>
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
C#
<lang csharp>using System;
namespace rref {
class Program { static void Main(string[] args) { int[,] matrix = new int[3, 4]{ { 1, 2, -1, -4 }, { 2, 3, -1, -11 }, { -2, 0, -3, 22 } }; matrix = rref(matrix); }
private static int[,] rref(int[,] matrix) { int lead = 0, rowCount = matrix.GetLength(0), columnCount = matrix.GetLength(1); for (int r = 0; r < rowCount; r++) { if (columnCount <= lead) break; int i = r; while (matrix[i, lead] == 0) { i++; if (i == rowCount) { i = r; lead++; if (columnCount == lead) { lead--; break; } } } for (int j = 0; j < columnCount; j++) { int temp = matrix[r, j]; matrix[r, j] = matrix[i, j]; matrix[i, j] = temp; } int div = matrix[r, lead]; for (int j = 0; j < columnCount; j++) matrix[r, j] /= div; for (int j = 0; j < rowCount; j++) { if (j != r) { int sub = matrix[j, lead]; for (int k = 0; k < columnCount; k++) matrix[j, k] -= (sub * matrix[r, k]); } } lead++; } return matrix; } }
}</lang>
Common Lisp
Direct implementation of the pseudo-code given.
<lang lisp>(defun convert-to-row-echelon-form (matrix)
(let* ((dimensions (array-dimensions matrix))
(row-count (first dimensions)) (column-count (second dimensions)) (lead 0))
(labels ((find-pivot (start lead)
(let ((i start)) (loop :while (zerop (aref matrix i lead)) :do (progn (incf i) (when (= i row-count) (setf i start) (incf lead) (when (= lead column-count) (return-from convert-to-row-echelon-form matrix)))) :finally (return (values i lead))))) (swap-rows (r1 r2) (loop :for c :upfrom 0 :below column-count :do (rotatef (aref matrix r1 c) (aref matrix r2 c)))) (divide-row (r value) (loop :for c :upfrom 0 :below column-count :do (setf (aref matrix r c) (/ (aref matrix r c) value)))))
(loop
:for r :upfrom 0 :below row-count :when (<= column-count lead) :do (return matrix) :do (multiple-value-bind (i nlead) (find-pivot r lead) (setf lead nlead) (swap-rows i r) (divide-row r (aref matrix r lead)) (loop :for i :upfrom 0 :below row-count :when (/= i r) :do (let ((scale (aref matrix i lead))) (loop :for c :upfrom 0 :below column-count :do (decf (aref matrix i c) (* scale (aref matrix r c)))))) (incf lead)) :finally (return matrix)))))</lang>
D
<lang d>import std.stdio, std.algorithm, std.array, std.conv;
void toReducedRowEchelonForm(T)(T[][] M) pure {
const nrows = M.length; if (!nrows) return; const ncols = M[0].length; int lead; foreach (r; 0 .. nrows) { if (ncols <= lead) return; int i = r; while (M[i][lead] == 0) if (nrows == ++i) { i = r; if (ncols == ++lead) return; } swap(M[i], M[r]); M[r][] /= M[r][lead]; foreach (j; 0 .. nrows) if (j != r) M[j][] -= M[r][] * M[j][lead]; lead++; }
}
void main() {
auto A = [[ 1, 2, -1, -4], [ 2, 3, -1, -11], [-2, 0, -3, 22]]; toReducedRowEchelonForm(A); writeln("[", array(map!text(A)).join("\n "), "]");
}</lang> Output:
[[1, 0, 0, -8] [0, 1, 0, 1] [0, 0, 1, -2]]
Euphoria
<lang euphoria>function ToReducedRowEchelonForm(sequence M)
integer lead,rowCount,columnCount,i sequence temp lead = 1 rowCount = length(M) columnCount = length(M[1]) for r = 1 to rowCount do if columnCount <= lead then exit end if i = r while M[i][lead] = 0 do i += 1 if rowCount = i then i = r lead += 1 if columnCount = lead then exit end if end if end while temp = M[i] M[i] = M[r] M[r] = temp M[r] /= M[r][lead] for j = 1 to rowCount do if j != r then M[j] -= M[j][lead]*M[r] end if end for lead += 1 end for return M
end function
? ToReducedRowEchelonForm(
{ { 1, 2, -1, -4 }, { 2, 3, -1, -11 }, { -2, 0, -3, 22 } })</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>
Go
From WP pseudocode: <lang go>package main
import "fmt"
type matrix [][]float64
func (m matrix) print() {
for _, r := range m { fmt.Println(r) } fmt.Println("")
}
func main() {
m := matrix{ { 1, 2, -1, -4}, { 2, 3, -1, -11}, {-2, 0, -3, 22}, } m.print() rref(m) m.print()
}
func rref(m matrix) {
lead := 0 rowCount := len(m) columnCount := len(m[0]) for r := 0; r < rowCount; r++ { if lead >= columnCount { return } i := r for m[i][lead] == 0 { i++ if rowCount == i { i = r lead++ if columnCount == lead { return } } } m[i], m[r] = m[r], m[i] f := 1 / m[r][lead] for j, _ := range m[r] { m[r][j] *= f } for i = 0; i < rowCount; i++ { if i != r { f = m[i][lead] for j, e := range m[r] { m[i][j] -= e * f } } } lead++ }
}</lang> Output (not so pretty, sorry)
[1 2 -1 -4] [2 3 -1 -11] [-2 0 -3 22] [1 0 0 -8] [-0 1 0 1] [-0 -0 1 -2]
Groovy
This solution implements the transformation to reduced row echelon form with optional pivoting. Options are provided for both partial pivoting and scaled partial pivoting. The default option is no pivoting at all. <lang groovy>enum Pivoting {
NONE({ i, it -> 1 }), PARTIAL({ i, it -> - (it[i].abs()) }), SCALED({ i, it -> - it[i].abs()/(it.inject(0) { sum, elt -> sum + elt.abs() } ) }); public final Closure comparer private Pivoting(Closure c) { comparer = c }
}
def isReducibleMatrix = { matrix ->
def m = matrix.size() m > 1 && matrix[0].size() > m && matrix[1..<m].every { row -> row.size() == matrix[0].size() }
}
def reducedRowEchelonForm = { matrix, Pivoting pivoting = Pivoting.NONE ->
assert isReducibleMatrix(matrix) def m = matrix.size() def n = matrix[0].size() (0..<m).each { i -> matrix[i..<m].sort(pivoting.comparer.curry(i)) matrix[i][i..<n] = matrix[i][i..<n].collect { it/matrix[i][i] } ((0..<i) + ((i+1)..<m)).each { k -> (i..<n).reverse().each { j -> matrix[k][j] -= matrix[i][j]*matrix[k][i] } } } matrix
}</lang>
This test first demonstrates the test case provided, and then demonstrates another test case designed to show the dangers of not using pivoting on an otherwise solvable matrix. Both test cases exercise all three pivoting options. <lang groovy>def matrixCopy = { matrix -> matrix.collect { row -> row.collect { it } } }
println "Tests for matrix A:" def a = [
[1, 2, -1, -4], [2, 3, -1, -11], [-2, 0, -3, 22]
] a.each { println it } println()
println "pivoting == Pivoting.NONE" reducedRowEchelonForm(matrixCopy(a)).each { println it } println() println "pivoting == Pivoting.PARTIAL" reducedRowEchelonForm(matrixCopy(a), Pivoting.PARTIAL).each { println it } println() println "pivoting == Pivoting.SCALED" reducedRowEchelonForm(matrixCopy(a), Pivoting.SCALED).each { println it } println()
println "Tests for matrix B (divides by 0 without pivoting):"
def b = [
[1, 2, -1, -4], [2, 4, -1, -11], [-2, 0, -6, 24]
] b.each { println it } println()
println "pivoting == Pivoting.NONE" try {
reducedRowEchelonForm(matrixCopy(b)).each { println it } println()
} catch (e) {
println "KABOOM! ${e.message}" println()
}
println "pivoting == Pivoting.PARTIAL" reducedRowEchelonForm(matrixCopy(b), Pivoting.PARTIAL).each { println it } println() println "pivoting == Pivoting.SCALED" reducedRowEchelonForm(matrixCopy(b), Pivoting.SCALED).each { println it } println()</lang>
Output:
Tests for matrix A: [1, 2, -1, -4] [2, 3, -1, -11] [-2, 0, -3, 22] pivoting == Pivoting.NONE [1, 0, 0, -8] [0, 1, 0, 1] [0, 0, 1, -2] pivoting == Pivoting.PARTIAL [1, 0.0, 0E-11, -7.9999999997000000000150] [0, 1, 0E-10, 0.999999999700000000010] [0, 0.0, 1, -2.00000000030] pivoting == Pivoting.SCALED [1, 0, 0, -8] [0, 1, 0, 1] [0, 0, 1, -2] Tests for matrix B (divides by 0 without pivoting): [1, 2, -1, -4] [2, 4, -1, -11] [-2, 0, -6, 24] pivoting == Pivoting.NONE KABOOM! Division undefined pivoting == Pivoting.PARTIAL [1, 0, 0.00, -3.00] [0, 1, 0.00, -2.00] [0, 0, 1, -3] pivoting == Pivoting.SCALED [1, 0, 0, -3] [0, 1, 0, -2] [0, 0, 1, -3]
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>
J
The reduced row echelon form of a matrix can be obtained using the gauss_jordan
verb from the linear.ijs script, available as part of the math/misc
addon.
<lang j> ]mymatrix=: _4]\ 1 2 _1 _4 2 3 _1 _11 _2 0 _3 22
1 2 _1 _4 2 3 _1 _11
_2 0 _3 22
require 'math/misc/linear' gauss_jordan mymatrix
1 0 0 _8 0 1 0 1 0 0 1 _2</lang>
Java
<lang java>import java.utils.Arrays;
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);
System.out.println(Arrays.deepToString(mtx)); }
}</lang>
JavaScript
for the print()
function.
Extends the Matrix class defined at Matrix Transpose#JavaScript <lang javascript>// modifies the matrix in-place Matrix.prototype.toReducedRowEchelonForm = function() {
var lead = 0; for (var r = 0; r < this.rows(); r++) { if (this.columns() <= lead) { return; } var i = r; while (this.mtx[i][lead] == 0) { i++; if (this.rows() == i) { i = r; lead++; if (this.columns() == lead) { return; } } }
var tmp = this.mtx[i]; this.mtx[i] = this.mtx[r]; this.mtx[r] = tmp;
var val = this.mtx[r][lead]; for (var j = 0; j < this.columns(); j++) { this.mtx[r][j] /= val; }
for (var i = 0; i < this.rows(); i++) { if (i == r) continue; val = this.mtx[i][lead]; for (var j = 0; j < this.columns(); j++) { this.mtx[i][j] -= val * this.mtx[r][j]; } } lead++; } return this;
}
var m = new Matrix([
[ 1, 2, -1, -4], [ 2, 3, -1,-11], [-2, 0, -3, 22]
]); print(m.toReducedRowEchelonForm()); print();
m = new Matrix([
[ 1, 2, 3, 7], [-4, 7,-2, 7], [ 3, 3, 0, 7]
]); print(m.toReducedRowEchelonForm());</lang> output
1,0,0,-8 0,1,0,1 0,0,1,-2 1,0,0,0.6666666666666663 0,1,0,1.666666666666667 0,0,1,1
Lua
<lang lua>function ToReducedRowEchelonForm ( M )
local lead = 1 local n_rows, n_cols = #M, #M[1]
for r = 1, n_rows do if n_cols <= lead then break end local i = r while M[i][lead] == 0 do i = i + 1 if n_rows == i then i = r lead = lead + 1 if n_cols == lead then break end end end M[i], M[r] = M[r], M[i]
local m = M[r][lead] for k = 1, n_cols do M[r][k] = M[r][k] / m end for i = 1, n_rows do if i ~= r then local m = M[i][lead] for k = 1, n_cols do M[i][k] = M[i][k] - m * M[r][k] end end end lead = lead + 1 end
end
M = { { 1, 2, -1, -4 },
{ 2, 3, -1, -11 }, { -2, 0, -3, 22 } }
res = ToReducedRowEchelonForm( M )
for i = 1, #M do
for j = 1, #M[1] do io.write( M[i][j], " " ) end io.write( "\n" )
end</lang> Output:
1 0 0 -8 0 1 0 1 0 0 1 -2
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>
MATLAB
<lang MATLAB>rref([1, 2, -1, -4; 2, 3, -1, -11; -2, 0, -3, 22])</lang>
Maxima
Maxima does not have a built-in function for reduced row echelon form, only echelon form.
<lang maxima>(%i28) echelon(matrix([ 1, 2, -1, -4],
[ 2, 3, -1, -11], [-2, 0, -3, 22]));
(%o28) matrix([1,0,3/2,-11],[0,1,-4/3,11/3],[0,0,1,-2])</lang>
However, rref can be programmed as a Maxima batch file[1]:
<lang maxima>request_rational_matrix(m, pos, fn) :=
if every('identity, map(lambda([s], every('ratnump,s)), args(m))) then true else print("Some entries in the matrix are not rational numbers. The result might be wrong.");
rowswap(m,i,j) := block([n, p, r],
require_matrix(m, "first", "rowswap"), require_integer(i, "second", "rowswap"), require_integer(j, "third", "rowswap"), n : length(m), if (i < 1) or (i > n) or (j < 1) or (j > n) then error("Array index out of bounds"), p : copymatrix(m), r : p[i], p[i] : p[j], p[j] : r, p);
addrow(m,i,j,k) := block([n,p],
require_matrix(m, "first", "addrow"), require_integer(i, "second", "addrow"), require_integer(j, "third", "addrow"), require_rational(k, "fourth", "addrow"), n : length(m), if (i < 1) or (i > n) or (j < 1) or (j > n) then error("Array index out of bounds"), p : copymatrix(m), p [i] : p[i] + k * p[j], p);
rowmul(m,i,k) := block([n,p],
require_matrix(m, "first", "addrow"), require_integer(i, "second", "addrow"), require_rational(k, "fourth", "addrow"), n : length(m), if (i < 1) or (i > n) then error("Array index out of bounds"), p : copymatrix(m), p [i] : k * p[i], p);
rref(m):= block([p,nr,nc,i,j,k,pivot,pivot_row,debug],
debug : 0, request_rational_matrix(m," ","rref"), nc: length(first(m)), nr: length(m), if nc = 0 or nr = 0 then error ("The argument to 'rref' must be a matrix with one or more rows and columns"),
p:copymatrix(m), ci : 1, cj : 1, while (ci<=nr) and (cj<=nc) do ( if (debug = 1) then (
disp(p), print("curseur en ligne ",ci," et colonne ",cj)),
pivot_row : 0, pivot : 0, for k : ci thru nr do ( if ( abs(p[k,cj]) > pivot ) then ( pivot_row : k, pivot : abs(p[k,cj]))), if (debug = 1) then
print("colonne ",cj," : pivot trouve ligne ", pivot_row,", valeur : ",pivot),
if (pivot = 0) then (cj : cj +1) else ( p : rowswap(p,ci,pivot_row),
if (debug = 1) then print (".. Echange : ",p),
p : rowmul(p,ci,1/p[ci,cj]),
if (debug = 1) then print (".. Normalisation : ",p),
for k : 1 thru nr do ( if not (k=ci) then (p : addrow (p,k,ci,-p[k,cj]))), ci : ci+1, cj : cj+1)), p
);</lang>
If this Maxima batch file is placed in one of the paths listed in the variable file_search_maxima
under the name rref.mac
, then it can be invoked as follows:
<lang maxima>load(rref)$ m:matrix(
[1,2,-1,-4], [2,3,-1,-11], [-2,0,-3,22]
); rref(m);</lang>
yielding the output: <lang maxima>matrix([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>
PARI/GP
PARI has a built-in matrix type, but no commands for row-echelon form. A dimension-limited one can be constructed from the built-in matsolve
command:
<lang parigp>rref(M)={
my(d=matsize(M)); if(d[1]+1 != d[2], error("Bad size in rref"), d=d[1]); concat(matid(d), matsolve(matrix(d,d,x,y,M[x,y]), M[,d+1]))
};</lang> Example: <lang parigp>rref([1,2,-1,-4;2,3,-1,-11;-2,0,-3,22])</lang> Output:
%1 = [1 0 0 -8] [0 1 0 1] [0 0 1 -2]
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>
Perl 6
<lang perl6>sub rref (@m is rw) {
@m or return; my ($lead, $rows, $cols) = 0, +@m, +@m[0]; for ^$rows -> $r { $lead < $cols or return @m; my $i = $r; until @m[$i][$lead] { ++$i == $rows or next; $i = $r; ++$lead == $cols and return @m; }
@m[$i, $r] = @m[$r, $i];
my $lv = @m[$r][$lead]; @m[$r] »/=» $lv;
for ^$rows -> $n { next if $n == $r; @m[$n] »-=» @m[$r] »*» @m[$n][$lead]; } ++$lead; } @m;
}
sub rat_or_int ($num is rw) {
return $num unless $num ~~ Rat; return $num.Int if $num.denominator == 1; return $num.perl;
}
sub say_it ($message, @array) {
say "\n$message"; $_».&rat_or_int.fmt(" %5s").say for @array;
}
my @M = (
[ # base test case [ 1, 2, -1, -4 ], [ 2, 3, -1, -11 ], [ -2, 0, -3, 22 ], ], [ # mix of number styles [ 3, 0, -3, 1 ], [ .5, 3/2, -3, -2 ], [ .2, 4/5, -1.6, .3 ], ], [ # degenerate case [ 1, 2, 3, 4, 3, 1], [ 2, 4, 6, 2, 6, 2], [ 3, 6, 18, 9, 9, -6], [ 4, 8, 12, 10, 12, 4], [ 5, 10, 24, 11, 15, -4], ], [ # larger matrix [1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, -1, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0], [0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, -1, 0, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0], [0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, -1, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0], [0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, -1, 0, 0, 0, 0, 0], [0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0], [0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, -1, 0], [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, -1, 0, 0], [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1], [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, -1, 0, 0, 0], ]
);
for @M -> @matrix {
say_it( 'Original Matrix', @matrix ); say_it( 'Reduced Row Echelon Form Matrix', rref(@matrix) ); say "\n";
}</lang>
Perl 6 handles rational numbers internally as a ratio of two integers to maintain precision. For some situations it is useful to return the ratio rather than the floating point result.
Output:
Original Matrix 1 2 -1 -4 2 3 -1 -11 -2 0 -3 22 Reduced Row Echelon Form Matrix 1 0 0 -8 0 1 0 1 0 0 1 -2 Original Matrix 3 0 -3 1 1/2 3/2 -3 -2 1/5 4/5 -8/5 3/10 Reduced Row Echelon Form Matrix 1 0 0 -41/2 0 1 0 -217/6 0 0 1 -125/6 Original Matrix 1 2 3 4 3 1 2 4 6 2 6 2 3 6 18 9 9 -6 4 8 12 10 12 4 5 10 24 11 15 -4 Reduced Row Echelon Form Matrix 1 2 0 0 3 4 0 0 1 0 0 -1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Original Matrix 1 0 0 0 0 0 1 0 0 0 0 -1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 -1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 -1 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 -1 0 0 0 0 1 0 0 0 0 0 0 1 0 0 -1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 -1 0 0 0 1 0 0 0 1 0 0 0 0 0 -1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 -1 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 -1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 -1 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 -1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 -1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 1 0 0 0 0 1 0 0 0 -1 0 0 0 Reduced Row Echelon Form Matrix 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 17/39 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4/13 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 20/39 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 28/39 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 19/39 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 8/39 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 11/39 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1/3 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 20/39 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 25/39 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 28/39 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 10/13 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 20/39 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 32/39
Re-implemented without the pseudocode, expressed as elementary matrix row operations. See http://unapologetic.wordpress.com/2009/08/27/elementary-row-and-column-operations/ and http://unapologetic.wordpress.com/2009/09/03/reduced-row-echelon-form/
First, a procedural version: <lang perl6>sub swap_rows ( @M, $r1, $r2 ) { @M[ $r1, $r2 ] = @M[ $r2, $r1 ] }; sub scale_row ( @M, $scale, $r ) { @M[$r] = @M[$r] X* $scale }; sub shear_row ( @M, $scale, $r1, $r2 ) { @M[$r1] = @M[$r1] Z+ ( @M[$r2] X* $scale ) }; sub reduce_row ( @M, $r, $c ) { scale_row( @M, 1/@M[$r][$c], $r ) }; sub clear_column ( @M, $r, $c ) {
for @M.keys.grep( * != $r ) -> $row_num { shear_row( @M, -1*@M[$row_num][$c], $row_num, $r ); }
}
my @M = (
[< 1 2 -1 -4 >], [< 2 3 -1 -11 >], [< -2 0 -3 22 >],
);
my $column_count = +@( @M[0] );
my $current_col = 0; while all( @M».[$current_col] ) == 0 {
$current_col++; return if $current_col == $column_count; # Matrix was all-zeros.
}
for @M.keys -> $current_row {
reduce_row( @M, $current_row, $current_col ); clear_column( @M, $current_row, $current_col ); $current_col++; return if $current_col == $column_count;
}
say @($_)».fmt(' %4g') for @M;</lang>
And the same code, recast into OO. Also, scale and shear are recast as unscale and unshear, which fit the problem better. <lang perl6>class Matrix is Array {
method unscale_row ( @M: $scale, $row ) { @M[$row] = @M[$row] X/ $scale; } method unshear_row ( @M: $scale, $r1, $r2 ) { @M[$r1] = @M[$r1] Z- ( @M[$r2] X* $scale ); } method reduce_row ( @M: $row, $col ) { @M.unscale_row( @M[$row][$col], $row ); } method clear_column ( @M: $row, $col ) { for @M.keys.grep( * != $row ) -> $scanning_row { @M.unshear_row( @M[$scanning_row][$col], $scanning_row, $row ); } } method reduced_row_echelon_form ( @M: ) { my $column_count = +@( @M[0] );
my $current_col = 0; # Skip past all-zero columns. while all( @M».[$current_col] ) == 0 { $current_col++; return if $current_col == $column_count; # Matrix was all-zeros. }
for @M.keys -> $current_row { @M.reduce_row( $current_row, $current_col ); @M.clear_column( $current_row, $current_col ); $current_col++; return if $current_col == $column_count; } }
}
my $M = Matrix.new.push(
[< 1 2 -1 -4 >], [< 2 3 -1 -11 >], [< -2 0 -3 22 >],
);
$M.reduced_row_echelon_form;
say @($_)».fmt(' %4g') for @($M);</lang>
Note that both versions can be simplified using Z+=, Z-=, X*=, and X/= to scale and shear. Currently, Rakudo has a bug related to Xop= and Zop=.
Note that the negative zeros in the output are innocuous, and also occur in the Perl 5 version.
PHP
<lang php><?php
function rref($matrix) {
$lead = 0; $rowCount = count($matrix); if ($rowCount == 0) return $matrix; $columnCount = 0; if (isset($matrix[0])) { $columnCount = count($matrix[0]); } for ($r = 0; $r < $rowCount; $r++) { if ($lead >= $columnCount) break; { $i = $r; while ($matrix[$i][$lead] == 0) { $i++; if ($i == $rowCount) { $i = $r; $lead++; if ($lead == $columnCount) return $matrix; } } $temp = $matrix[$r]; $matrix[$r] = $matrix[$i]; $matrix[$i] = $temp; } { $lv = $matrix[$r][$lead]; for ($j = 0; $j < $columnCount; $j++) { $matrix[$r][$j] = $matrix[$r][$j] / $lv; } } for ($i = 0; $i < $rowCount; $i++) { if ($i != $r) { $lv = $matrix[$i][$lead]; for ($j = 0; $j < $columnCount; $j++) { $matrix[$i][$j] -= $lv * $matrix[$r][$j]; } } } $lead++; } return $matrix;
} ?></lang>
PicoLisp
<lang PicoLisp>(de reducedRowEchelonForm (Mat)
(let (Lead 1 Cols (length (car Mat))) (for (X Mat X (cdr X)) (NIL (loop (T (seek '((R) (n0 (get R 1 Lead))) X) @ ) (T (> (inc 'Lead) Cols)) ) ) (xchg @ X) (let D (get X 1 Lead) (map '((R) (set R (/ (car R) D))) (car X) ) ) (for Y Mat (unless (== Y (car X)) (let N (- (get Y Lead)) (map '((Dst Src) (inc Dst (* N (car Src))) ) Y (car X) ) ) ) ) (T (> (inc 'Lead) Cols)) ) ) Mat )</lang>
Output:
(reducedRowEchelonForm '(( 1 2 -1 -4) ( 2 3 -1 -11) (-2 0 -3 22)) ) -> ((1 0 0 -8) (0 1 0 1) (0 0 1 -2))
Python
This closely follows the pseudocode given. If you want accurate results, use the fractions package.
<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>
REXX
Reduced Row Rchelon Form on a matrix, with optimization added. <lang rexx> /*REXX program to do Reduced Row Echelon Form (RREF) on a matrix. */
cols=0 /*maximum columns in any row. */ maxW=0 /*maximum width of any element. */ @.= /*matrix to be constructed. */ mat.=
mat.1=' 1 2 -1 -4 ' mat.2=' 2 3 -1 -11 ' mat.3=' -2 0 -3 22 '
do r=1 until mat.r== /*build @.row.col from mat.n */ _=mat.r do c=1 until _= parse var _ @.r.c _ maxW=max(maxW,length(@.r.c)) end cols=max(cols,c) end
rows=r-1 /*adjust the row count. */ maxW=maxW+1 /*bump the max width, better view*/ call showMat 'original matrix' /*show the original matrix. */
/*───────────────────────────────────Reduced Row Echelon Form on matrix.*/ !=1 /*set the pointer to one. */
do r=1 for rows while cols>! /*start to do the heavy lifting. */ j=r do while @.j.!==0 j=j+1 if j==rows then do j=r !=!+1; if cols==! then leave r end end /*while*/
do w=1 for cols while j\==r /*swap rows J,R (but not if same)*/ _.j.w=@.j.w @.j.w=@.r.w @.r.w=_.j.w end /*w*/
?=@.r.!
do d=1 for cols while ?\=1 /*divide row J by @.r.p--unless 1*/ @.r.d=@.r.d/? end do k=1 for rows /*sub (row K) *@.r.s from row K */ if k==r then iterate /*skip if row k is the same as R */ ?=@.k.! do s=1 for cols while ?\=0 /*but not if @.r.! is 0*/ @.k.s=@.k.s-?*@.r.s end end /*k*/ !=!+1 end /*r*/
call showMat 'matrix RREF' /*show reduced row echelon form. */ exit
/*─────────────────────────────────────SHOWMAT subroutine───────────────*/
showMat: parse arg title; say
say center(title,3+(cols+1)*maxW,'─'); say /*build a pretty title.*/
do r=1 for rows _= do c=1 for cols if @.r.c== then do; say; say '*** error! ***'; say say "matrix element isn't defined:" say 'row' row", column" c'.'; say exit 13 end _=_ right(@.r.c,maxW) end /*c*/ say _ end /*r*/
say return </lang> Output:
────original matrix──── 1 2 -1 -4 2 3 -1 -11 -2 0 -3 22 ──────matrix RREF────── 1 0 0 -8 0 1 0 1 0 0 1 -2
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
Sage
<lang sage> sage: m = matrix(ZZ, [[1,2,-1,-4],[2,3,-1,-11],[-2,0,-3,22]]) sage: m.rref() [ 1 0 0 -8] [ 0 1 0 1] [ 0 0 1 -2] </lang>
Scheme
<lang scheme>(define (reduced-row-echelon-form matrix)
(define (clean-down matrix from-row column) (cons (car matrix) (if (zero? from-row) (map (lambda (row) (map - row (map (lambda (element) (/ (* element (list-ref row column)) (list-ref (car matrix) column))) (car matrix)))) (cdr matrix)) (clean-down (cdr matrix) (- from-row 1) column)))) (define (clean-up matrix until-row column) (if (zero? until-row) matrix (cons (map - (car matrix) (map (lambda (element) (/ (* element (list-ref (car matrix) column)) (list-ref (list-ref matrix until-row) column))) (list-ref matrix until-row))) (clean-up (cdr matrix) (- until-row 1) column)))) (define (normalise matrix row with-column) (if (zero? row) (cons (map (lambda (element) (/ element (list-ref (car matrix) with-column))) (car matrix)) (cdr matrix)) (cons (car matrix) (normalise (cdr matrix) (- row 1) with-column)))) (define (repeat procedure matrix indices) (if (null? indices) matrix (repeat procedure (procedure matrix (car indices) (car indices)) (cdr indices)))) (define (iota start stop) (if (> start stop) (list) (cons start (iota (+ start 1) stop)))) (let ((indices (iota 0 (- (length matrix) 1)))) (repeat normalise (repeat clean-up (repeat clean-down matrix indices) indices) indices)))</lang>
Example: <lang scheme>(define matrix
(list (list 1 2 -1 -4) (list 2 3 -1 -11) (list -2 0 -3 22)))
(display (reduced-row-echelon-form matrix)) (newline)</lang> Output: <lang>((1 0 0 -8) (0 1 0 1) (0 0 1 -2))</lang>
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
<lang ti89b>rref([1,2,–1,–4; 2,3,–1,–11; –2,0,–3,22])</lang>
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>
References
- ↑ a version of rref by Antoine Chambert-Loir/France, taken from UT Math: Maxima Mailing List