I am trying to implement my own LU decomposition with partial pivoting. My code is below and apparently is working fine, but for some matrices it gives different results when comparing with the built-in [L, U, P] = lu(A)
function in matlab
Can anyone spot where is it wrong?
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
Here are the two matrices I've tested with. The first one is correct, whereas the second has some elements inverted.
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];
UPDATE
I have checked my code and corrected some bugs, but still there's something missing with the partial pivoting. In the first column the last two rows are always inverted (compared with the result of lu() in matlab)
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
I forgot that If there was a swap in matrix P I had to swap also the matrix L. So just add the next line after after swapping P and everything will work excellent.
L([k r],:) = L([r k],:);