%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Constructs an LU decomposition of a matrix A
% Input:
% A - matrix
% Output:
% [ L, U, P ] - a list consisting of three matrices L,U,P
% such that LU = PA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [L,U,P] = lu_decomp(A)
% Check that A is a square matrix
if size(A,1) ~= size(A,2)
error( 'Error: Matrix is not square.' );
end
n = size(A,1); % size of A
% Construct matrices L, U, P
U = A; % U is initially coincides with A
L=eye(n); % L is initially a unit matrix
% Vector Perm will save the permutation of rows of A
% Initially it is equal to [1,2,...,n];
Perm=1:n;
% Consider columns from the first to the last but one
for col=1:n-1
% 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
isInvert=0;
for row=col:n
if U(row,col) ~= 0
isInvert=1;
break;
end
end
if isInvert==0
error( 'Matrix is not invertible!' );
end
% Exchange col-th and row-th elements of P if necessary
if row>col
% save a permutation
i=Perm(col);
Perm(col)=Perm(row);
Perm(row)=i;
if col>1
i = L(col, 1:col-1);
L(col,1:col-1) = L(row, 1:col-1);
L(row,1:col-1) = i;
end
% 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
%Perform row operations
for row=[col+1:n]
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
end
% create a permutation matrix from vector Perm
P = zeros(n);
for i=1:n
P(i,Perm(i)) = 1;
end
end