algorithmgauss

Variant of LU-decomposition without pivoting


I want to write a variant of an algorithm that takes a n x n matrix A and overwrites A with its LU-decomposition, meaning that the strict lower-triangular part of the returned matrix (without the main diagonal) represents L and the upper-triangular part U.

I know the classical implementation of this algorithm that builds L column-wise and U row-wise:

# classical implementation (originated in gaussian-elimination)
for i = 1, ..., n-1 do 
    for j = i + 1, …, n do
        A_j,i = A_j,i / A_i,i;              # update L
        for k = i + 1, …, n do
            A_j,k = A_j,k - A_j,i * A_i,k;  # update R
        end for
    end for
end for

I am now interested in a variant that builds L row- and U column-wise changing each non-trivial entry exactly once:

# variant
for i = 2, ..., n do 
    for j = 1, …, i - 1 do
        A_i,j = ???;        # update L
        for k = 1, …, j do
            A_j+1,i = ???;  # update R    
        end for
    end for
end for

I guess the second ??? should be A_j+1,i = A_j+1,i - A_j+1,k * A_k,i; but I do not know how to update the entries of L (on which in turn the entries on U depend on)


Solution

  • Building L row-wise and U column-wise will not yield an efficient algorithm. The reason lies in the order of the LU multiplication. Matrix multiplication requires the dot product of row of the first matrix (L in this case) and column of the second (U in this case). Doolittle algorithm caches partial sums (dot products) in this fashion, resulting in the existing loop pattern (L column-wise, etc.)