LDU Decomposition By MATLAB

(LDU Decomposition만을 진행하며, 예외 처리 등은 없는 코드입니다. 직접 작성해 오류나 비효율적인 부분이 존재할 수 있습니다.)

 

function [L,D,U] = decomLDU(A)
% Decompose the given matrix A as LDU form
    i = 1;
    j = 1;
    [R, C] = size(A); % square matrix A
    if (R == C)
        N = R; % since A is N by N. Maybe redundant.
        L = eye(N);
        D = eye(N);
            
        while i <= N && j <= N
              L(i,j) = A(i,j);
              D(i,j) = A(i,j);
              A(i,j:N) = A(i,j:N)./A(i,j);
          for k = [i+1:N]
              L(k,j) = A(k,j);
              A(k,j:N) = A(k,j:N) - A(k,j).*A(i,j:N);
          end
          i = i + 1;
          j = j + 1;
        end
    U = A;
    L = L/D;
    else
        L = NaN;
    end
end