[解決済み] 部分ピボッティングによるLU分解 Matlab

質問

部分ピボットによるLU分解を独自に実装しようとしています。私のコードは以下のとおりで、どうやらうまくいっているようなのですが、いくつかの行列については、組み込みの [L, U, P] = lu(A) 関数

どこが間違っているのか、どなたかわかりますか？

function [L, U, P] = lu_decomposition_pivot(A)
    n = size(A,1);
    Ak = A;
    L = zeros(n);
    U = zeros(n);
    P = eye(n);
    for k = 1:n-1
        for i = k+1:n
            [~,r] = max(abs(Ak(:,k)));

            Ak([k r],:) = Ak([r k],:);
            P([k r],:) = P([r k],:);

            L(i,k) = Ak(i,k) / Ak(k,k);
            for j = k+1:n
                U(k,j-1) = Ak(k,j-1);
                Ak(i,j) = Ak(i,j) - L(i,k)*Ak(k,j);
            end
        end
    end
    L(1:n+1:end) = 1;
    U(:,end) = Ak(:,end);
return

以下は、私がテストした2つのマトリックスです。最初のものは正しいのですが、2番目のものはいくつかの要素が反転しています。

A = [1 2 0; 2 4 8; 3 -1 2];

A = [0.8443 0.1707 0.3111;
     0.1948 0.2277 0.9234;
     0.2259 0.4357 0.4302];

アップデイト

私のコードをチェックし、いくつかのバグを修正しましたが、部分ピボットについてはまだ何かが欠けています。最初の列では、最後の2行は常に反転しています（matlabのlu()の結果と比較して）。

function [L, U, P] = lu_decomposition_pivot(A)
    n = size(A,1);
    Ak = A;
    L = eye(n);
    U = zeros(n);
    P = eye(n);
    for k = 1:n-1
        [~,r] = max(abs(Ak(k:end,k)));
        r = n-(n-k+1)+r;
        Ak([k r],:) = Ak([r k],:);
        P([k r],:) = P([r k],:);
        for i = k+1:n
            L(i,k) = Ak(i,k) / Ak(k,k);
            for j = 1:n
                U(k,j) = Ak(k,j);
                Ak(i,j) = Ak(i,j) - L(i,k)*Ak(k,j);
            end
        end
    end
    U(:,end) = Ak(:,end);
return

解決方法は？

行列Pを入れ替えた場合、行列Lも入れ替えなければならないことを忘れていました。そこで、Pを入れ替えた後に次の行を追加すれば、すべてがうまくいきます。

L([k r],:) = L([r k],:);