0

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.

Blake
  • 149
  • 1
  • 8
  • What is partial pivoting? Only searching the first non-zero column of the remaining lower matrix block instead of the whole block? – Lutz Lehmann Dec 12 '16 at 19:23
  • Yes, that is right. But in this algorithm it is implemented in different way (I think). – Blake Dec 12 '16 at 20:30

1 Answers1

1

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