LU/LUP decomposition: can some of U's diagonal elements be zero?

519 Views Asked by At

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?