next up previous
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_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:

# 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 up previous
Next: Help File Up: Maple Implementation Previous: Maple Implementation



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