CODEWITHSUNDEEP
numerical-methodsmatrix-decomposition
Blake
Blake

Reputation: 157

How does LU decomposition with partial pivoting work?

I am using the method in which initially the elements on the main diagonal of L are set to ones (think that is Doolittle’s method, but not sure because I have seen it named differently). I know there tons of documentations, papers and books but I could not find simple examples that were not using Gauss elimination for finding L.

Upvotes: 0

Views: 2147

Answers (1)

user1911226
user1911226

Reputation:

Partial Pivoting, as compared to full pivoting, uses row interchanging only as compared to full pivoting which also pivots columns. The primary purpose of partial pivoting as shown below in the picture and the code is to swap the rows to find the maximum u there as to avoid dividing by a very small one in that for loop which would cause a large condition number.

If you try a naive implementation of the LU decomposition and some ill-conditioned matrix like some arbitrary diagonally dominant matrix it should explode.

enter image description here

function [L,U,P] = my_lu_piv(A)
n = size(A,1);
I = eye(n);
O = zeros(n);
L = I;
U = O;
P = I;
    function change_rows(k,p)
        x = P(k,:); P(k,:) = P(p,:); P(p,:) = x;
        x = A(k,:); A(k,:) = A(p,:); A(p,:) = x;
        x = v(k); v(k) = v(p); v(p) = x;
    end

    function change_L(k,p)
        x = L(k,1:k-1); L(k,1:k-1) = L(p,1:k-1);
        L(p,1:k-1) = x;
    end
for k = 1:n
    if k == 1, v(k:n) = A(k:n,k);
    else
        z = L(1:k-1,1:k -1)\ A(1:k-1,k);
        U(1:k-1,k) = z;
        v(k:n) = A(k:n,k)-L(k:n,1:k-1)*z;
    end
    if k<n
        x = v(k:n); p = (k-1)+find(abs(x) == max(abs(x))); % find index p
        change_rows(k,p);
        L(k+1:n,k) = v(k+1:n)/v(k);
        if k > 1, change_L(k,p); end
    end
    U(k,k) = v(k);
end
end

Upvotes: 1

Related Questions