I'm learning LUP decomposition. So far I've wrote Doolittle implementation in GNU Octave:
function [L,U,P] = mylu(A)
[r,c] = size(A)
if (r != c)
error("not square\n")
endif
N = r
P = eye(N)
U = A
L = eye(N)
#iterate over columns
for n = 1:N-1
if U(n,n) == 0
U([n n+1], :) = U([n+1 n], :)
P([n n+1], :) = P([n+1 n], :)
endif
#calculate L
L_temp = eye(N)
# down the column
for i = n+1:N
L_temp(i,n) = - U(i,n)/U(n,n)
L(i,n) = -L_temp(i,n)
endfor
U = L_temp*U
endfor
P = P'
if !isequal(P*A,L*U)
error("LU failed\n")
endif
endfunction
which seems to be working, at least P*A really equals L*U.
I go from 1-st to (N-1)-th column eliminating elements under diagonal.
However, the last diagonal element of U is always zero. Is it inherent to the algorithm or am I doing something wrong? Can any of $u_{n,n}$ (diagonal elements of U) be zero? If not, how should I handle it?