In [2]:
function forward_sweep(L,b)
y = zeros(size(L,1))
for i=1:size(L,1)
r = b[i]
for j=1:i-1
r -= L[i,j]*y[j]
end
y[i] = r/L[i,i]
end
return y
end

function backward_sweep(U,y)
x = zeros(size(U,1))
for i=size(U,1):-1:1 # reverse order
r = y[i]
for j=i+1:size(U,2)
r -= U[i,j]*x[j]
end
x[i] = r/U[i,i]
end
return x
end

L = [1.0  0 0
1.0  1 0
-1.0  0 1]
U = [1.0 -1 2
0.0 1  1
0.0 0  4]
A = L*U

b = [3.0; 2.0; 1.0]
@show A \ b
y = forward_sweep(L,b)
@show x = backward_sweep(U,y)
@show norm(x - A\b);

A \ b = [-1.0,-2.0,1.0]
x = backward_sweep(U,y) = [-1.0,-2.0,1.0]
norm(x - A \ b) = 0.0

In [3]:
function lusteps(A)
A = deepcopy(A)
A0 = copy(A)
println()
println("LU factorization of")
println("A = ")
println(A)
println()

# Begin real stuff...
n = size(A,1)
@assert n == size(A,2)
L = eye(n)
for k=1:n-1 # outer step of Gaussian Elimination
for i = k+1:n # for each row
L[i,k] = A[i,k]/A[k,k]
for j = k:n # for each column
A[i,j] = A[i,j] - L[i,k]*A[k,j]
end
end
# end real stuff!

@printf("after step %i\n", k)
println("L = ")
println(L)
println("L^{-1} A = ")
println(L\A0)
println()
#@show L*A, A0, L*A - A0
end

# One more real thing!
U = triu(A)
return L, U
end
#lusteps([3.0 6 9; 2 5 -2; 1 3 -1])
lusteps(A);

LoadError: UndefVarError: A not defined

In [4]:
A = [3.0 6 9; 2 5 -2; 1 3 -1];
b = [39; 3; 2.0]

L,U = lusteps(A)
x = backward_sweep(U,forward_sweep(L,b))
@show x
@show norm(x - A\b)

LU factorization of
A =
[3.0 6.0 9.0
2.0 5.0 -2.0
1.0 3.0 -1.0]

after step 1
L =
[1.0 0.0 0.0
0.6666666666666666 1.0 0.0
0.3333333333333333 0.0 1.0]
L^{-1} A =
[3.0 6.0 9.0
0.0 1.0 -8.0
0.0 1.0 -4.0]

after step 2
L =
[1.0 0.0 0.0
0.6666666666666666 1.0 0.0
0.3333333333333333 1.0 1.0]
L^{-1} A =
[3.0 6.0 9.0
0.0 1.0 -8.0
0.0 0.0 4.0]

x = [2.0,1.0,3.0]
norm(x - A \ b) = 0.0

Out[4]:
0.0
In [ ]:
A = [1.0 -1 2;
1.0 0  3;
-1.0 1  2]

lusteps(A);

In [6]:
A = [1.0 -1 2;
-1.0 1  2
1.0 0  3]
L,U = lusteps(A);
@show L,U

LU factorization of
A =
[1.0 -1.0 2.0
-1.0 1.0 2.0
1.0 0.0 3.0]

after step 1
L =
[1.0 0.0 0.0
-1.0 1.0 0.0
1.0 0.0 1.0]
L^{-1} A =
[1.0 -1.0 2.0
0.0 0.0 4.0
0.0 1.0 1.0]

after step 2
L =
[1.0 0.0 0.0
-1.0 1.0 0.0
1.0 Inf 1.0]
L^{-1} A =
[1.0 -1.0 2.0
0.0 0.0 4.0
NaN NaN -Inf]

(L,U) = (
[1.0 0.0 0.0
-1.0 1.0 0.0
1.0 Inf 1.0],

[1.0 -1.0 2.0
0.0 0.0 4.0
0.0 0.0 -Inf])

Out[6]:
(
3x3 Array{Float64,2}:
1.0    0.0  0.0
-1.0    1.0  0.0
1.0  Inf    1.0,

3x3 Array{Float64,2}:
1.0  -1.0     2.0
0.0   0.0     4.0
0.0   0.0  -Inf  )
In [1]:
function swaprows!(A,i,r,cols=:)
temp = A[i,cols]
A[i,cols] = A[r,cols]
A[r,cols] = temp
end
function lupsteps(A)
A = deepcopy(A)
A0 = copy(A)
println()
println("Pivoted LU factorization of")
println("A = ")
println(A)
println()

# Begin real stuff...
n = size(A,1)
@assert n == size(A,2)
L = eye(n)
P = eye(n)
for k=1:n-1 # outer step of Gaussian Elimination
# find Pivot
r = indmax(abs(A[k:n,k]))
@show r
if r > 1 # then we need to pivot
swaprows!(A,k,k+r-1)
swaprows!(P,k,k+r-1)
swaprows!(L,k,k+r-1,1:k-1)
end
for i = k+1:n # for each row
L[i,k] = A[i,k]/A[k,k]
for j = k:n # for each column
A[i,j] = A[i,j] - L[i,k]*A[k,j]
end
end
# end real stuff!

@printf("after step %i\n", k)
println("L = ")
println(L)
println("L^{-1} P A = ")
println(L\(P*A0))
println()
#@show L*A, A0, L*A - A0
end

# One more real thing!
U = triu(A)
return L, U, P
end

Out[1]:
lupsteps (generic function with 1 method)
In [ ]:
In [5]:
lupsteps([1.0e-20 1; 1 1])

Pivoted LU factorization of
A =
[1.0e-20 1.0
1.0 1.0]

r = 2
after step 1
L =
[1.0 0.0
1.0e-20 1.0]
L^{-1} P A =
[1.0 1.0
0.0 1.0]


Out[5]:
(
2x2 Array{Float64,2}:
1.0      0.0
1.0e-20  1.0,

2x2 Array{Float64,2}:
1.0  1.0
0.0  1.0,

2x2 Array{Float64,2}:
0.0  1.0
1.0  0.0)
In [ ]:
In [8]:
A = [1.0 -1 2;
-1.0 1  2
1.0 0  3]
lupsteps(A)

Pivoted LU factorization of
A =
[1.0 -1.0 2.0
-1.0 1.0 2.0
1.0 0.0 3.0]

r = 1
after step 1
L =
[1.0 0.0 0.0
-1.0 1.0 0.0
1.0 0.0 1.0]
L^{-1} P A =
[1.0 -1.0 2.0
0.0 0.0 4.0
0.0 1.0 1.0]

r = 2
after step 2
L =
[1.0 0.0 0.0
1.0 1.0 0.0
-1.0 0.0 1.0]
L^{-1} P A =
[1.0 -1.0 2.0
0.0 1.0 1.0
0.0 0.0 4.0]


Out[8]:
(
3x3 Array{Float64,2}:
1.0  0.0  0.0
1.0  1.0  0.0
-1.0  0.0  1.0,

3x3 Array{Float64,2}:
1.0  -1.0  2.0
0.0   1.0  1.0
0.0   0.0  4.0,

3x3 Array{Float64,2}:
1.0  0.0  0.0
0.0  0.0  1.0
0.0  1.0  0.0)
In [9]:
lupsteps([1.0 2 3;
4.0 5 6;
7.0 8 0])

Pivoted LU factorization of
A =
[1.0 2.0 3.0
4.0 5.0 6.0
7.0 8.0 0.0]

r = 3
after step 1
L =
[1.0 0.0 0.0
0.5714285714285714 1.0 0.0
0.14285714285714285 0.0 1.0]
L^{-1} P A =
[7.0 8.0 0.0
0.0 0.4285714285714288 6.0
0.0 0.8571428571428572 3.0]

r = 2
after step 2
L =
[1.0 0.0 0.0
0.14285714285714285 1.0 0.0
0.5714285714285714 0.5000000000000002 1.0]
L^{-1} P A =
[7.0 8.0 0.0
0.0 0.8571428571428572 3.0
0.0 0.0 4.499999999999999]


Out[9]:
(
3x3 Array{Float64,2}:
1.0       0.0  0.0
0.142857  1.0  0.0
0.571429  0.5  1.0,

3x3 Array{Float64,2}:
7.0  8.0       0.0
0.0  0.857143  3.0
0.0  0.0       4.5,

3x3 Array{Float64,2}:
0.0  0.0  1.0
1.0  0.0  0.0
0.0  1.0  0.0)
In [10]:
lupsteps(randn(4,4))

Pivoted LU factorization of
A =
[1.3536006267991767 -0.6203480402843149 1.5010001572343763 1.610451834998707
0.16405135363831408 -0.16784458699816776 -0.3265101827949094 -0.5244591110414828
0.44477684695433733 -1.1647885807349385 -0.28029146783200054 -0.5358697932291413
-0.5729518937189758 1.46590400214269 2.0053370063498295 -1.1892766968680581]

r = 1
after step 1
L =
[1.0 0.0 0.0 0.0
0.1211962748763215 1.0 0.0 0.0
0.3285879439979939 0.0 1.0 0.0
-0.4232798673223282 0.0 0.0 1.0]
L^{-1} P A =
[1.3536006267991767 -0.6203480402843149 1.5010001572343763 1.610451834998707
0.0 -0.09266071538888257 -0.5084258104404886 -0.7196398743110625
0.0 -0.9609496936147308 -0.7735020234383099 -1.065044850599163
0.0 1.2033231659574788 2.64068015375479 -0.5076048578208056]

r = 3
after step 2
L =
[1.0 0.0 0.0 0.0
-0.4232798673223282 1.0 0.0 0.0
0.3285879439979939 -0.7985799000637601 1.0 0.0
0.1211962748763215 -0.07700401522242187 0.0 1.0]
L^{-1} P A =
[1.3536006267991767 -0.6203480402843149 1.5010001572343763 1.610451834998707
0.0 1.2033231659574788 2.64068015375479 -0.5076048578208056
0.0 0.0 1.3352920698475448 -1.4704078872295812
0.0 0.0 -0.3050828356832075 -0.758727486509671]

r = 1
after step 3
L =
[1.0 0.0 0.0 0.0
-0.4232798673223282 1.0 0.0 0.0
0.3285879439979939 -0.7985799000637601 1.0 0.0
0.1211962748763215 -0.07700401522242187 -0.22847648284022232 1.0]
L^{-1} P A =
[1.3536006267991767 -0.6203480402843149 1.5010001572343763 1.610451834998707
0.0 1.2033231659574788 2.64068015375479 -0.5076048578208056
0.0 0.0 1.3352920698475448 -1.4704078872295812
0.0 0.0 0.0 -1.094681108924408]


Out[10]:
(
4x4 Array{Float64,2}:
1.0        0.0        0.0       0.0
-0.42328    1.0        0.0       0.0
0.328588  -0.79858    1.0       0.0
0.121196  -0.077004  -0.228476  1.0,

4x4 Array{Float64,2}:
1.3536  -0.620348  1.501     1.61045
0.0      1.20332   2.64068  -0.507605
0.0      0.0       1.33529  -1.47041
0.0      0.0       0.0      -1.09468 ,

4x4 Array{Float64,2}:
1.0  0.0  0.0  0.0
0.0  0.0  0.0  1.0
0.0  0.0  1.0  0.0
0.0  1.0  0.0  0.0)