% ピボット選択型LU分解
% n次の正方行列AのLU分解PA = LUを計算するピボット選択型アルゴリズム

A = [ 0, 50, 20; ...
      2,  1, 10; ...
     10,  0,  5 ];
n = 3;

L = eye(n,n);
U = zeros(n,n);
row = 1:n;
for k = 1:n
    % ピボット選択
    [value,p] = max(abs(A(k:n,k)));
    p = p + (k-1);
    % ピボット選択に伴う行列AとLの更新
    A([k,p],k:n) = A([p,k],k:n);
    L([k,p],1:k-1) = L([p,k],1:k-1);
    row([k,p]) = row([p,k]);
    % 行列UとLの計算
    U(k,k:n) = A(k,k:n);
    L(k+1:n,k) = A(k+1:n,k)/U(k,k);
    % 行列Aの更新
    A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - L(k+1:n,k)*U(k,k+1:n);
end

% 置換行列Pの計算
P = zeros(n,n);
for k = 1:n
    P(k,row(k)) = 1;
end

P
L
U

fprintf("P'*L*U\n");
P'*L*U
