Program to compute the Row Echelon Decomposition of a general input matrix A.
Copyright © Rob Corless and David Jeffrey, March, 1990. Revised for share library June 1991
Updated to use Normalizer, and added RowEchelonSolve and BackSubstitute, July 1992
The Row Echelon Decomposition [36] of A is a factorization A = PLUR, where P is a permutation matrix, L is a unit lower triangular matrix, U is upper triangular, and R is the (unique) Row Echelon Form of A.
Note that if A is m by n, then P, L, and U are all m by
m and
P and L are nonsingular always. U is assumed nonsingular:
the row echelon form R that is computed is only correct if
U is in fact nonsingular.
This means that the determinant of U must be checked. If
for some values of the parameters, then the row echelon
form must be recomputed for those values of the parameters.
Call:
> R := RowEchelon(A,'dt'); # or > R := RowEchelon(A,'dt','rank','P','L','U');The names
'rank', 'P', 'L', and 'U' indicate
optional arguments (though of course if you want U you
must also ask for rank, P, and L).
The determinant of U will be returned in 'dt', and
the rank of A will be returned in 'rank' if present.
If you wish to use RootOf's in your procedure, you may
specify evala@Normal as the normalizer by simply assigning
(before you call RowEchelon)
> Normalizer := x -> evala(Normal(x));In Maple V, algebraic numbers can be handled this way. If you wish to handle algebraic functions, then you can use instead
> Normalizer := x -> normal(simplify(x));
In Maple 5 Release 2, the normalizer evala should work.
Likewise, evalc is recommended as a normalizer if your matrix
has complex elements in it.
This routine uses userinfo. If you wish progress information
on the execution, type
infolevel[RowEchelon] := 1;before starting. If you wish also to see intermediate results, type
SPM_quotinfolevel[RowEchelon] := 2;". Matrices will be printed if
infolevel is 5 or above.
RowEchelon := proc(A,dt)
local d,localdt,rank,P,L,U,R,
i,j,k,m,n,SparseId,
PivotColumn,
pivot,test,multiplier,
Pt,Ut,Utinv,Lt,Ltinv;
Before we continue with the program, we note that Maple provides
some facility for creating sparse matrices, i.e.
matrices whose elements are mostly zero. Sparseness is `fragile'
in Maple, and the result of multiplying two sparse matrices together
is a dense matrix, as is that of multiplying a sparse matrix with a
dense one. Nevertheless using sparse matrices in the following
programs saves some space.
The following short utility program provides a modifiable identity matrix, using the sparse representation to save space.
SparseId := proc(n)
local i,t;
t := array(1..n,1..n,sparse);
for i to n do t[i,i] := 1; od;
op(t)
end;
if not type(dt,name) then
ERROR(`Expecting a variable name in the 2nd argument
to put the determinant of U in: instead got `,dt);
fi;
userinfo(1,'RowEchelon',`Initialization phase.`);
m := linalg[rowdim](A);
n := linalg[coldim](A);
# PivotColumn[i] is the column in which we find the pivot
# taken from row i. In the nonsingular case this will always
# be equal to i, but in the singular case we may pivot on some
# column to the right instead. 1 <= PivotColum[i] <= n.
PivotColumn := array(1..m);
localdt := 1; # For some reason, a local
# must be used to accumulate dt
# Accumulate the factorization only if we will be returning it.
if nargs > 3 then P := SparseId(m); fi;
if nargs > 4 then L := SparseId(m); fi;
if nargs > 5 then U := SparseId(m); fi;
# The Row Echelon Form of A is the same size as A, and is
# initially equal to A.
R := array(1..m,1..n);
for i to m do
PivotColumn[i] := i;
for j to n do
R[i,j] := A[i,j];
od;
od;
rank := m; # Initially, the rank is thought to be m.
# At this point, A = P L U R, though P, L, and U are
# only identity matrices. dt is always det(U), so dt = 1
# to start.
userinfo(1,'RowEchelon',`Initialization complete.`);
userinfo(2,'RowEchelon', `determinant = `,localdt);
if nargs > 2 then userinfo(5,'RowEchelon',`initial rank = `,rank); fi;
if nargs > 3 then userinfo(5,'RowEchelon',`P = `,print(P)); fi;
if nargs > 4 then userinfo(5,'RowEchelon',`L = `,print(L)); fi;
if nargs > 5 then userinfo(5,'RowEchelon',`U = `,print(U)); fi;
userinfo(5,'RowEchelon',`R = `,print(R));
# Forward Elimination:
userinfo(1,'RowEchelon',`Forward Elimination phase.`);
for i to m while i <= rank do
pivot := 0;
Pt := SparseId(m);
userinfo(1,'RowEchelon',`Choose nonzero pivot.`);
while (pivot = 0 and PivotColumn[i] <= n ) do
for j from i to m while pivot = 0 do
test := Normalizer(R[j,PivotColumn[i]]);
if test <> 0 then
pivot := test;
if i < m and i <> j then
Pt[i,j] := 1; Pt[j,i] := 1; Pt[i,i] := 0; Pt[j,j] := 0;
if nargs > 2 then P := evalm(P &* Pt); fi;
if nargs > 3 then L := evalm( Pt &* L &* Pt); fi;
R := evalm( Pt &* R)
fi
fi;
od;
# A = P L U R is a loop invariant for the above while loop,
# though P changes, so L, U, and R must also change.
# Since U is the identity matrix at this point P and Pt
# commute through U leaving no change.
if (pivot = 0) then
for j from i to m do
PivotColumn[j] := PivotColumn[j] + 1;
od;
fi;
od;
userinfo(1,'RowEchelon',`Pivot chosen.`);
if PivotColumn[i] > n then
rank := i-1;
userinfo(1,'RowEchelon',`Zero pivot found; rank reduced.`);
d := 1
else
userinfo(2,'RowEchelon',
`Pivot from row`,i,
`is from column`,PivotColumn[i],
`and is = `,pivot);
d := pivot
fi;
# The heart of the matter: keep track of the determinant of U,
# by keeping track of the product of the pivots. This could be
# done closer to where U is calculated, but it is available now
# so why not?
localdt := localdt * d;
userinfo(2,'RowEchelon', `determinant = `,localdt);
# Define the elementary matrix to be used for
# forward elimination for this column.
Lt := SparseId(m);
Ltinv := SparseId(m);
for k from i+1 to m while PivotColumn[i] <= n do
multiplier := R[k,PivotColumn[i]];
Lt[k,i] := Normalizer(multiplier/d);
Ltinv[k,i] := -Lt[k,i];
od;
if nargs > 4 then L := map(Normalizer,evalm( L &* Lt)); fi;
R := map(Normalizer,evalm( Ltinv &* R));
# Again after this statement A = P L U R
# though now L has changed, so U and R had to.
# Again since U = I, Lt and Ltinv commute through it.
if nargs > 4 then userinfo(5,'RowEchelon',`L = `,print(L)); fi;
od;
# Now do back substitution:
userinfo(1,'RowEchelon',`Back Substitution phase.`);
for i from m by -1 to 1 do
# Define the elementary matrix for back substitution of this
# column.
Ut := SparseId(m);
Utinv := SparseId(m);
if PivotColumn[i] <= n then
d := R[i,PivotColumn[i]];
else
d := 1;
fi;
Ut[i,i] := d;
Utinv[i,i] := 1/d;
for k to i-1 while PivotColumn[i] <= n do
Ut[k,i] := R[k,PivotColumn[i]];
Utinv[k,i] := -Ut[k,i]/d;
od;
if nargs > 5 then U := map(Normalizer,evalm(U &* Ut)); fi;
R := map(Normalizer,evalm(Utinv &* R));
# Again A = P L U R is maintained, though now U is changing so
# R must also change.
userinfo(1,'RowEchelon',`Back substitution of this column complete.`);
if nargs > 5 then userinfo(5,'RowEchelon',`U = `,print(U)); fi;
userinfo(5,'RowEchelon',`R = `,print(R));
od;
userinfo(1,'RowEchelon',`Row Echelon Decomposition completed.`);
# As maintained throughout, A = P L U R but now R is in Row Echelon
# Form, P is a permutation matrix, L is unit lower triangular, and
# U is presumed nonsingular. The determinant of U has been accumulated
# in localdt, and will be returned in dt.
dt := localdt;
if (nargs > 2) then assign(args[3],rank); fi;
if (nargs > 3) then assign(args[4],op(P)); fi;
if (nargs > 4) then assign(args[5],op(L)); fi;
if (nargs > 5) then assign(args[6],op(U)); fi;
op(R);
end:
# BackSubstituteUtility for solution of upper-triangular systems by back-substitution.
Procedure takes an upper-triangular rectangular matrix and possibly a right-hand side and solves the associated linear system. One can call it with a square matrix and a right-hand side, or one can call it with an augmented matrix.
BackSubstitute := proc(U,FreeName,B,N,M)
local i,j,k,s,x,y,localFreeName,n,m,rhszero;
if nargs < 1 then
ERROR(`Expecting an upper triangular matrix`)
fi;
if nargs < 5 then
m := linalg[coldim](U)
else
m := M
fi;
if nargs < 4 then
n := linalg[rowdim](U)
else
n := N
fi;
if nargs < 3 then
y := array(1..n);
for i to n do y[i] := 0 od
else
y := B
fi;
if nargs < 2 then
localFreeName := FreeParameter # Use a global name
else
localFreeName := FreeName
fi;
x := array(1..n);
for i to n do x[i] := cat(localFreeName,n-i+1) od:
# Free parameters to start
# Now start the backsubstitution proper, taking care to avoid
# division by zero diagonal elements.
i := n;
while i > 0 do
j := i; # U is guaranteed upper triangular
while (j <= n) and (Normalizer(U[i,j]) = 0) do j := j + 1 od;
if j > n then
rhszero := true;
for k from n+1 to m do
rhszero := rhszero and (Normalizer(U[i,k])=0);
od;
if not rhszero then
userinfo(1,BackSubstitute,`Inconsistent Equations`);
RETURN(NULL);
fi;
i := i - 1;
else
s := 0;
for k from n+1 to m do s := s - U[i,k]; od;
for k from j+1 to n do
s := s + U[i,k]*x[k] # May involve free parameters
od;
x[j] := Normalizer( (y[i] - s)/U[i,j] );
i := i - 1
fi
od;
eval(x)
end:
# RowEchelonSolveUser-callable linear system solver.
Use the Row Echelon Decomposition, with its proviso, to solve a linear system. Any free parameters are returned using the global names FreeParameter1, FreeParameter2, etc.
RowEchelonSolve := proc(A,b,dt)
local Aug,R;
Aug := linalg[augment](A,b);
R := RowEchelon(Aug,dt);
BackSubstitute(R);
end:
save `help/text/RowEchelon`,RowEchelon,
RowEchelonSolve, BackSubstitute, `echelon.m`;
quit
The last lines save all the required procedures into the Maple-readable
file `echelon.m' for later loading. This statement is required if
you wish to use readlib to load your package, or if you wish
to put your package in the share library.