In [1]:
## LU Decompositions
function lusimple(Aref) # julia always passes by reference, so make sure to make a copy
A = copy(Aref)
n = size(A,1)
L = eye(n,n)

for i = 1:n-1
alpha = A[i,i]
b = A[i,(i+1):end]
c = A[(i+1):end,i]
L[(i+1):end,i] = c/alpha
A[(i+1):end,(i+1):end] = A[(i+1):end,(i+1):end] - c*b'/alpha
end
U = triu(A)
return L,U
end

A = [2.0 1 -1; -3 -1 2; -2 1 2]
L,U = lusimple(A)
display(L)
display(U)

3×3 Array{Float64,2}:
1.0  0.0  0.0
-1.5  1.0  0.0
-1.0  4.0  1.0
3×3 Array{Float64,2}:
2.0  1.0  -1.0
0.0  0.5   0.5
0.0  0.0  -1.0
In [2]:
## Next we want to see what happens for a few cases

function lusimple_show_work(Aorig)

A = copy(Aorig)
A = convert(Array{Float64,2},A)
n = size(A,1)
L = eye(n,n)

for i = 1:n-1
@printf("After %i steps ... \n", i-1)
@printf("L = \n")
display(L)
@printf("L^{-1} A = \n");
display(L\Aorig)

sleep(2)

alpha = A[i,i]
b = A[i,(i+1):end]
c = A[(i+1):end,i]
L[(i+1):end,i] = c/alpha
A[(i+1):end,(i+1):end] = A[(i+1):end,(i+1):end] - c*b'/alpha

end
U = triu(A)
return L,U
end

A = [1.0 2 4; -2 -4 6; 3 2 1]
display(rank(A))

L,U = lusimple_show_work(A)
@printf("Done!")
display(L)
display(U)

3
After 0
3×3 Array{Float64,2}:
1.0  0.0  0.0
0.0  1.0  0.0
0.0  0.0  1.0
In [3]:
##
A = [2.0 1 -1; -3 -1 2; -2 1 2];
L,U = lusimple(A)
display(L)
display(U)

3×3 Array{Float64,2}:
1.0  0.0  0.0
-1.5  1.0  0.0
-1.0  4.0  1.0
3×3 Array{Float64,2}:
2.0  1.0  -1.0
0.0  0.5   0.5
0.0  0.0  -1.0
In [4]:
##
A = [1e-20 1; 1 1]
L,U = lusimple(A)

Out[4]:
([1.0 0.0; 1.0e20 1.0], [1.0e-20 1.0; 0.0 -1.0e20])
In [5]:
##
A = [1 1; 1e-20 1]
L,U = lusimple(A)

Out[5]:
([1.0 0.0; 1.0e-20 1.0], [1.0 1.0; 0.0 1.0])
In [6]:
## A more advanced LU code that doesn't allocate any additional storage
# but still no pivoting. This still needs a bit of work :(
# don't worry if you don't understand this, I'm just providing it as a
# "get started" code in case you want to see what you'd do to make the
# Julia LU code start to go faster.
using SugarBLAS
function lubetter!(A) # this will modify A in place, and store L in the lower-triangle
n = size(A,1)
for i = 1:n-1
alpha = A[i,i]
b = @view A[i,(i+1):end] # this doesn't make a copy
c = @view A[(i+1):end,i]
B = @view A[(i+1):end,(i+1):end]
@blas! c = c/alpha
SugarBLAS.@ger! B -= 1.0*c*b' # c changed in place!
end
return A
end
function lubetter(A)
B = lubetter!(copy(A))
L = tril(B,1) + eye(B) # could improve this step!
U = triu!(B,0)
return L, U
end

# Precompile
lubetter(randn(10,10));
lusimple(randn(10,10));

# Time
@time lubetter(randn(1024,1024));
@time lubetter!(randn(1024,1024)); # if you didn't need L, U explicitly!
@time lusimple(randn(1024,1024));