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);
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);
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)
A = [1.0 -1 2;
1.0 0 3;
-1.0 1 2]
lusteps(A);
A = [1.0 -1 2;
-1.0 1 2
1.0 0 3]
L,U = lusteps(A);
@show L,U
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
lupsteps([1.0e-20 1; 1 1])
A = [1.0 -1 2;
-1.0 1 2
1.0 0 3]
lupsteps(A)
lupsteps([1.0 2 3;
4.0 5 6;
7.0 8 0])
lupsteps(randn(4,4))