Next: Help File Up: Maple Implementation Previous: Maple Implementation

## Source Code

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_quot`infolevel[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:
```

```# BackSubstitute
```
Utility 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:
```

```# RowEchelonSolve
```
User-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.

Next: Help File Up: Maple Implementation Previous: Maple Implementation

Robert Corless
Wed Jan 31 11:33:59 EST 1996