How to get the correct answer for LU Decomposition with partial pivoting?

I am trying to solve a matrix using LU decomposition with crout's method with partial pivoting. I found this code online at this website.

https://gist.github.com/02015678/353ea94c27e4c5f89c9d

I added few lines to get the my x from Ax=b equation but it seems that the code is not giving me the correct matrices at the end. Sometimes it even gives me error. I can't figure out where in the code is giving a problem.

A = [0   1  -1   1   0    0    0  ;
     1  1   0   0   0    0    0  ;
     0   0     0   -1  0   1   0  ;
     0   1     0    0 -1  -0.3   0  ;
     0   0     0    0  0    0    1  ;
     0   1    0    0  0    0   -1  ;
     0   1    0    0  1    0    0  ]   % given matrix to solve
b = [0; 120; 0; 100; 0; 0; 0]

[n,n]=size(A);
L=eye(n); P=L; U=A;
for k=1:n
    [pivot m]=max(abs(U(k:n,k)));
    m=m+k-1;
    if m~=k
        % interchange rows m and k in U
        temp=U(k,:);
        U(k,:)=U(m,:);
        U(m,:)=temp;
        % interchange rows m and k in P
        temp=P(k,:);
        P(k,:)=P(m,:);
        P(m,:)=temp;
        if k >= 2
            temp=L(k,1:k-1);
            L(k,1:k-1)=L(m,1:k-1);
            L(m,1:k-1)=temp;
        end
    end
    for j=k+1:n
        L(j,k)=U(j,k)/U(k,k);
        U(j,:)=U(j,:)-L(j,k)*U(k,:);
    end
end
L,U

Here is the code I added to get x.

Y=zeros(n,1);
Y(1)=b(1)/L(1,1);
for j=2:n
    Y(j)=(b(j)-L(j,1:j-1)*Y(1:j-1))/L(j,j);
end
Y
X=zeros(n,1);
X(n)=Y(n)/U(n,n);
for j=n-1:-1:1
    X(j)=(Y(j)-U(j,j+1:n)*X(j+1:n))/U(j,j);
end
X

When I solve it with x=A\b, i am getting completely different answers. What should I do?

Answer 1

[Jan]

You forgot to apply the permutation matrix P to b. The results are correct:

P * L * U - A   % == zeros, fine
b = P * b;

Then A\b replies the same as your X.

A simplification of your code:

% Replace:
temp=U(k,:);
U(k,:)=U(m,:);
U(m,:)=temp;

% by:
U([m, k], :) = U([k, m], :);

[Kenny Gee]

thank you!