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")
  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], :)
    #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)
    U = L_temp*U
  P = P'
  if !isequal(P*A,L*U)
    error("LU failed\n")

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?

  • 1
    $\begingroup$ It's possible if $A$ is singular, otherwise it can't happen. However it may happen in the process that a pivot is zero even if the matrix is nonsingular, and then you need to take care of this: usually, you swap rows or columns, and it's common to always swap to get the maximum possible pivot (in absolute value). $\endgroup$ – Jean-Claude Arbaut Jun 4 '18 at 14:07
  • $\begingroup$ I've checked and A really was singular. Thank you. $\endgroup$ – olegst Jun 5 '18 at 21:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.