The Maple session below shows how you use the addrow commands in Maple to step through a modified Gaussian Elimination algorithm that finds an LU decomposition of a matrix. Recall, L is the lower triangular factor and U is the upper triangular factor. Note: two of the addrow steps are shown. The rest of the steps are omitted to save space. At the end, the product of LU is computed and shown equal to A.
NOTE: The ACC would like us to keep the amount of printing to a minimum. Please do not print this file unless absolutely necessary.
--------------------------------------------------------------------------------
> with(linalg);
Warning: new definition for norm
Warning: new definition for trace
[BlockDiagonal, GramSchmidt, JordanBlock, Wronskian, add, addcol, addrow, adj,
adjoint, angle, augment, backsub, band, basis, bezout, blockmatrix,
charmat, charpoly, col, coldim, colspace, colspan, companion, concat, cond,
copyinto, crossprod, curl, definite, delcols, delrows, det, diag, diverge,
dotprod, eigenvals, eigenvects, entermatrix, equal, exponential, extend,
ffgausselim, fibonacci, frobenius, gausselim, gaussjord, genmatrix, grad,
hadamard, hermite, hessian, hilbert, htranspose, ihermite, indexfunc,
innerprod, intbasis, inverse, ismith, iszero, jacobian, jordan, kernel,
laplacian, leastsqrs, linsolve, matrix, minor, minpoly, mulcol, mulrow,
multiply, norm, normalize, nullspace, orthog, permanent, pivot, potential,
randmatrix, randvector, rank, ratform, row, rowdim, rowspace, rowspan,
rref, scalarmul, singularvals, smith, stack, submatrix, subvector,
sumbasis, swapcol, swaprow, sylvester, toeplitz, trace, transpose,
vandermonde, vecpotent, vectdim, vector]
--------------------------------------------------------------------------------
> A:=matrix(4,4,[4,1,-1,0,1,3,-1,0,-1,-1,5,2,0,0,2,4]);
[ 4 1 -1 0 ]
[ ]
[ 1 3 -1 0 ]
A := [ ]
[ -1 -1 5 2 ]
[ ]
[ 0 0 2 4 ]
--------------------------------------------------------------------------------
> L:=matrix(4,4);
L := array(1 .. 4, 1 .. 4, [])
--------------------------------------------------------------------------------
> for i from 1 to 4 do
> for j from 1 to 4 do
> L[i,j]:=0;
> od;
> od;
--------------------------------------------------------------------------------
> evalm(L);
[ 0 0 0 0 ]
[ ]
[ 0 0 0 0 ]
[ ]
[ 0 0 0 0 ]
[ ]
[ 0 0 0 0 ]
--------------------------------------------------------------------------------
> evalm(A);
[ 4 1 -1 0 ]
[ ]
[ 1 3 -1 0 ]
[ ]
[ -1 -1 5 2 ]
[ ]
[ 0 0 2 4 ]
--------------------------------------------------------------------------------
> m:=-A[2,1]/A[1,1];
m := -1/4
--------------------------------------------------------------------------------
> U:=addrow(A,1,2,m);
[ 4 1 -1 0 ]
[ ]
[ 0 11/4 -3/4 0 ]
U := [ ]
[ -1 -1 5 2 ]
[ ]
[ 0 0 2 4 ]
--------------------------------------------------------------------------------
> L[2,1]:=-m;
L[2, 1] := 1/4
--------------------------------------------------------------------------------
> m:=-U[3,1]/U[1,1];
m := 1/4
--------------------------------------------------------------------------------
> U:=addrow(U,1,3,m);
[ 4 1 -1 0 ]
[ ]
[ 0 11/4 -3/4 0 ]
U := [ ]
[ 0 -3/4 19/4 2 ]
[ ]
[ 0 0 2 4 ]
--------------------------------------------------------------------------------
> L[3,1]:=-m;
L[3, 1] := -1/4
--------------------------------------------------------------------------------
> m:=-U[4,3]/U[3,3];
11
m := - ----
25
--------------------------------------------------------------------------------
> U:=addrow(U,3,4,m);
[ 4 1 -1 0 ]
[ ]
[ 0 11/4 -3/4 0 ]
[ ]
[ 50 ]
U := [ 0 0 ---- 2 ]
[ 11 ]
[ ]
[ 78 ]
[ 0 0 0 ---- ]
[ 25 ]
--------------------------------------------------------------------------------
> L[4,3]:=-m;
11
L[4, 3] := ----
25
--------------------------------------------------------------------------------
> evalm(L);
[ 0 0 0 0 ]
[ ]
[ 1/4 0 0 0 ]
[ ]
[ -1/4 -3/11 0 0 ]
[ ]
[ 11 ]
[ 0 0 ---- 0 ]
[ 25 ]
--------------------------------------------------------------------------------
> for i from 1 to 4 do
> L[i,i]:=1;
> od;
L[1, 1] := 1
L[2, 2] := 1
L[3, 3] := 1
L[4, 4] := 1
--------------------------------------------------------------------------------
> evalm(L);
[ 1 0 0 0 ]
[ ]
[ 1/4 1 0 0 ]
[ ]
[ -1/4 -3/11 1 0 ]
[ ]
[ 11 ]
[ 0 0 ---- 1 ]
[ 25 ]
--------------------------------------------------------------------------------
> evalm(L&*U);
[ 4 1 -1 0 ]
[ ]
[ 1 3 -1 0 ]
[ ]
[ -1 -1 5 2 ]
[ ]
[ 0 0 2 4 ]
--------------------------------------------------------------------------------
> evalm(A);
[ 4 1 -1 0 ]
[ ]
[ 1 3 -1 0 ]
[ ]
[ -1 -1 5 2 ]
[ ]
[ 0 0 2 4 ]
--------------------------------------------------------------------------------
>