** 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

> 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

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`; quitThe 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.

Wed Jan 31 11:33:59 EST 1996