with(LinearAlgebra): ####################################################################### # Constructs an LU decomposition of a matrix A # Input: # A - matrix # Output: # [ P, L, U ] - a list consisting of three matrices L,U,P # such that LU = PA ####################################################################### lu := proc( A :: Matrix( numeric ) ) local N, i : local L, U, Perm, P : local col, row : local tmp, f : # Check that A is a square matrix if ( Dimension(A)[2] <> Dimension(A)[1] ) then error( "Matrix is not square." ) : end if : N := Dimension(A)[1] : # Construct matrices L, U, P # U is initially coincides with A U := Copy(A) ; # L is initially a unit matrix L := Matrix(N, N) : for i from 1 by 1 to N do L(i,i):=1 : end do : # Vector Perm will save the permutation of rows of A # Initially it is equal to [1,2,...,N]; Perm := Vector(1..N, i->i) : # Consider columns from the first to the last but one for col from 1 by 1 to N - 1 do # Find a row that does not have a zero in this column # If all entries on and below the diagonal are 0 # the matrix is not invertible for row from col by 1 to N do if ( U(row,col) <> 0 ) then break : end if : end do : if ( row > N ) then error( "Matrix is not invertible!" ) : end if : # Exchange col-th and row-th elements of P if necessary if ( row > col ) then # save a permutation i := Perm(col) : Perm(col) := Perm(row) : Perm(row) := i : if ( col > 1 ) then i := L(col, 1..col-1) : L(col, 1..col-1) := L(row, 1..col-1) : L(row,1..col-1) := i : end if : # Swap rows in U. # Entries to the left of col have already been annihilated. tmp := U(col,col..N) : U(col,col..N) := U(row,col..N) : U(row,col..N) := tmp : end if : #Perform row operations for row from col + 1 by 1 to N do f := U(row,col) / U(col,col) : U(row,col) := 0 : U(row,col+1..N) := U(row,col+1..N) - f * U(col,col+1..N) : L(row, col) := f : end do : end do : # create a permutation matrix from vector Perm P := Matrix(N, N) : for i from 1 by 1 to N do P(i,Perm(i)) := 1 : end do : return [P, L, U] : end proc :