LU 分解与部分旋转 Matlab

Question

我正在尝试通过部分旋转来实现我自己的 LU 分解。我的代码在下面，显然工作正常，但对于某些矩阵，与内置的相比，它给出了不同的结果[L, U, P] = lu(A)matlab 中函数

谁能看出哪里错了？

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

这是我测试过的两个矩阵。第一个是正确的，而第二个有一些颠倒的元素。

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];

更新

我已经检查了我的代码并纠正了一些错误，但是部分旋转仍然缺少一些东西。在第一列中，最后两行总是反转（与 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

Answer 1

我忘记了，如果矩阵 PI 中有交换，则还必须交换矩阵 L。因此，只需在交换 P 后添加下一行，一切都会很好地工作。

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