This file demonstrates the A=LU decomposition both with and without partial pivoting.
The determinant of the A matrix can be found from the U matrix.
Any section/theorem/algorithm references are for "Matrix Computations" by Golub and Loan.with( LinearAlgebra ):The input matrix which we are computing the LU decomposition and the determinant for.n := 3;
A := RandomMatrix(n,n);We use the native Maple implementation to verify that our implementation is correct.
Note that the determinant is found by multiplying the diagonal values of the U matrix. (Theorem 3.2.1)
Note also that row switching by the LU decomposition is not yet taken into account with the determinant sign.
The evalf statement forces the pivoting to use the largest magnitude value in each column (thanks to Robert Israel, comp.soft-sys.math.maple)(mP, mL, mU) := LUDecomposition( evalf(A), output=['P','L','U'] );
mDet := Determinant( A );
mDetUpp := mul( i, i=Diagonal( mU ) ); LU Decomposition without Partial Pivoting.
Section 3.2.4 "Upper Triangularizing"
Section 3.2.7 "Where is L?"g_vec := Vector(n):
L := IdentityMatrix(n,n,compact=false):
U := Matrix( A ):
for k from 1 to n do
# Without partial pivoting, we cannot handle case where leading element on diagonal is zero
if U[k,k] = 0 then break end if:
# Find multipliers needed to make leading value of remaining rows zero (use to subtract in next step)
g_vec[k+1..n] := U[k+1..n,k] / U[k,k]:
# Ignoring rows already computed, subtract an outter product from rest of matrix
U[k+1..n,1..n] := U[k+1..n,1..n] - g_vec[k+1..n] . U[k,1..n]:
# Each column of L is made from the multipliers for the associated row of U
L[k+1..n,k] := g_vec[k+1..n];
end do:
(L,U);
det := mul( i, i=Diagonal( U ) );Outer Product Gaussian Elimination without Partial Pivoting
Section 3.2.6 "Some Practical Details"
Note that both the L and U matrices are packed in the overwritten input matrixoA := Matrix( A ):
for k from 1 to (n-1) do
# Without partial pivoting, we cannot handle case where leading element on diagonal is zero
if oA[k,k] = 0 then break end if:
# Stash first L column multipliers
oA[k+1..n,k] := oA[k+1..n,k] / oA[k,k];
# Apply multipliers to submatrix
oA[k+1..n,k+1..n] := oA[k+1..n,k+1..n] - oA[k+1..n,k] . oA[k,k+1..n];
end do:
oA;Gauss Elimination with Partial Pivoting
Section 3.4.3 "Partial Pivoting Details"gA := Matrix( A ):
p := Vector( n-1 ):
for k from 1 to (n-1) do
# Find partial pivot (note that we pulled vector norm out of loop as invariant)
cNorm := VectorNorm( gA[k..n,k], infinity );
for mu from k to n do
if abs(gA[mu,k]) = cNorm then break end if:
end do;
# Swap row with largest pivot to top of current submatrix
tempRow := gA[mu,k..n];
gA[mu,k..n] := gA[k,k..n];
gA[k,k..n] := tempRow;
# Update permutation vector with swap we just made
p[k] := mu;
# Use outer product gaussian elimination on remaining submatrix
if gA[k,k] <> 0 then
# Stash first L column multipliers
gA[k+1..n,k] := gA[k+1..n,k] / gA[k,k];
# Apply multipliers to submatrix
gA[k+1..n,k+1..n] := gA[k+1..n,k+1..n] - gA[k+1..n,k] . gA[k,k+1..n];
end if;
end do:
gA;
p;