Why doesn't the L in my partial pivot LU decomposition switch rows?

Everything else works well besides my L, please help.

A= [1 2 3 1; 4 5 6 2; 7 8 0 4; 0 1 3 1];
GE = GE_partial(A);

L = 4×4

1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1

L = 4×4

1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1

L = 4×4

1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1

%%
function[L,U,P] = GE_partial(A)
[n,n] = size(A); %nxn matrix A
L = eye(n,n); 
P = eye(n,n);
U = A;

for j =1:n-1 % looping over the columns of U
    %Finding the row width with the largest magnitude 
    [~, r] = max(abs(U(j:n,j)));
    % Calling the index of this row k
    k = r + j - 1;
    % Swap rows of j and k for U and P
    U([j,k],:) = U([k,j],:); %Swaping the first and the third row.
    P([j, k],:) = P([k, j],:); % swaping the first ad third row
    % Swap rows j and k in L for coulmns 1 to j-1
    for c = 1:j-1
        L([j,k],c) = L([k,j],c) 
    end
end
end

[Torsten]

Output j and k in the outer loop (for j =1:n-1) and you'll see the reason.

[Jessica Arroyo]

For U and P my rows are swapped. The output for L is the identity.

[Torsten]

j=1, k=3 means: your loop to modify L is not entered. Do you understand ?

Answer 1

[Piyush Patil]

Hello Jessica,

As @Torsten correctly mentioned, you should check the values of variables j and k in the outer “for” loop.

You will notice that –

For j = 1, value of k = 3 but since value of j is 1, so it won’t enter the for loop that modifies the matrix L

For j = 2, value of k = 2. Therefore, both j and k are pointing to same row and hence swapping of rows won’t happen.

Similarly, for j = 3, value of k = 3 and hence swapping won’t take place.

Now, to add on to this discussion, since you are trying to perform LU Decomposition with partial pivoting, you need to modify matrix L and U by computing multipliers and eliminating the elements below pivot element.

You can do this by adding the following logic to the code –

% Compute multipliers and eliminate elements below pivot
for i=j+1:n
    L(i,j) = U(i,j)/U(j,j);
    U(i,:) = U(i,:) - L(i,j)*U(j,:);
end  

Also, you can replace the for loop that switches rows of matrix L by a single line –

L([j, k], 1:j-1) = L([k, j], 1:j-1);

Both will give the same result.

So, finally your entire function should look like –

function[L,U,P] = GE_partial(A)
    [n,n] = size(A); %nxn matrix A
    L = eye(n,n); 
    P = eye(n,n);
    U = A;
    for j =1:n-1 % looping over the columns of U
        %Finding the row width with the largest magnitude 
        [~, r] = max(abs(U(j:n,j)));
        % Calling the index of this row k
        k = r + j - 1;
        % Swap rows of j and k for U, P, and L
        U([j, k],:) = U([k, j],:);
        P([j, k],:) = P([k, j],:);
        L([j, k], 1:j-1) = L([k, j], 1:j-1);
        % Compute multipliers and eliminate elements below pivot
        for i=j+1:n
            L(i,j) = U(i,j)/U(j,j);
            U(i,:) = U(i,:) - L(i,j)*U(j,:);
        end
    end
end