## A simple algorithm to solve a linear system
using LinearAlgebra
function solve1_pivot!(A::Matrix, b::Vector)
m,n = size(A)
@assert(m==n, "the system is not square")
@assert(n==length(b), "vector b has the wrong length")
if n==1
return [b[1]/A[1]]
else
# let's make sure we have an equation
# that we can eliminate!
α = A[1,1]
newrow = 1
if α == 0
for j=2:n
if A[j,1] != 0
newrow = j
break
end
end
if newrow == 1
error("the system is singular")
end
end
# swap rows 1, and newrow
if newrow != 1
tmp = A[1,:]
A[1,:] .= A[newrow,:]
A[newrow,:] .= tmp
b[1], b[newrow] = b[newrow], b[1]
end
D = A[2:end,2:end]
c = A[1,2:end]
d = A[2:end,1]
α = A[1,1]
y = solve1_pivot!(D-d*c'/α, b[2:end]-b[1]/α*d)
γ = (b[1] - c'*y)/α
return pushfirst!(y,γ)
end
end
A = [0 1.0; 1.0 1.0]
b = [5; 6.0]
xtrue = A\b
xmy = solve1_pivot!(A,b) # this changes A.
@show norm(xtrue-xmy)
##
A = [1e-18 1.0; 1.0 1.0]
b = [5; 6.0]
xtrue = A\b
xmy = solve1_pivot!(A,b)
@show norm(xtrue-xmy)
##
##
## A simple algorithm to solve a linear system
using LinearAlgebra
function solve1_pivot2(A::Matrix, b::Vector)
m,n = size(A)
@assert(m==n, "the system is not square")
@assert(n==length(b), "vector b has the wrong length")
if n==1
return [b[1]/A[1]]
else
# let's make sure we have an equation
# that we can eliminate!
# let's try that again, where we pick the
# largest magnitude entry!
maxval = abs(A[1,1])
newrow = 1
for j=2:n
if abs(A[j,1]) > maxval
newrow = j
maxval = abs(A[j,1])
end
end
if maxval < eps(1.0)
error("the system is singular")
end
@show newrow
# swap rows 1, and newrow
if newrow != 1
tmp = A[1,:]
A[1,:] .= A[newrow,:]
A[newrow,:] .= tmp
b[1], b[newrow] = b[newrow], b[1]
end
D = A[2:end,2:end]
c = A[1,2:end]
d = A[2:end,1]
α = A[1,1]
y = solve1_pivot2(D-d*c'/α, b[2:end]-b[1]/α*d)
γ = (b[1] - c'*y)/α
return pushfirst!(y,γ)
end
end
A = [1e-18 1.0; 1.0 1.0]
b = [5; 6.0]
xtrue = A\b
xmy = solve1_pivot2(A,b)
@show norm(xtrue-xmy)
##
function myreduce_all_pivot(A::Matrix)
A = copy(A) # save a copy
n = size(A,1)
L = Matrix(1.0I,n,n)
U = Matrix(1.0I,n,n)
d = zeros(n)
p = collect(1:n)
for i=1:n-1
maxval = abs(A[i,i])
newrow = i
for j=i+1:n
if abs(A[j,i]) > maxval
newrow = j
maxval = abs(A[j,i])
end
end
if maxval < eps(1.0)
error("the system is singular")
end
j = newrow
# swap the ith row/column
tmp = A[i,:]
A[i,:] .= A[j,:]
A[j,:] .= tmp
p[i],p[j] = p[j], p[i]
L[i,1:i-1], L[j,1:i-1] = L[j,1:i-1], L[i,1:i-1]
α = A[i,i]
d[i] = α
U[i,i+1:end] = A[i,i+1:end]/α
L[i+1:end,i] = A[i+1:end,i]/α
A[i+1:end,i+1:end] -= A[i+1:end,i]*A[i,i+1:end]'/α
end
d[n] = A[n,n]
return L,U,d,p
end
A = randn(4,4)
L,U,d,p = myreduce_all_pivot(A)
display(L*Diagonal(d)*U - A[p,:])
##
A = [1e-18 1.0; 1.0 1.0]
b = [5; 6.0]
L,U,d,p = myreduce_all_pivot(A)
L*Diagonal(d)*U - A[p,:]
##
A = randn(4,4)
L,U,d,p = myreduce_all_pivot(A)
norm(L*Diagonal(d)*U - A[p,:])
##
A = [1e-18 1; 1e-18 1]
b = [5; 6]
L,U,d,p = myreduce_all_pivot(A)
##
lu(A)
## We can re-use an LU factorization
A = randn(4,4)
F = lu(A)
F\randn(4)
F\randn(4)
## Your experiment -- play around with how long A\b takes vs. F\b