## A simple algorithm to solve a linear system
using LinearAlgebra
function solve1(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
D = A[2:end,2:end]
c = A[1,2:end]
d = A[2:end,1]
α = A[1,1]
y = solve1(D-d*c'/α, b[2:end]-b[1]/α*d)
γ = (b[1] - c'*y)/α
return pushfirst!(y,γ)
end
end
A = randn(4,4)
b = randn(4)
xtrue = A\b
xmy = solve1(A,b)
@show norm(xtrue-xmy)
##
function myreduce1(A)
n = size(A,1)
D = A[2:end,2:end]
c = A[1,2:end]
d = A[2:end,1]
α = A[1,1]
L = Matrix(1.0I,n,n)
L[2:end,1] = -d/α
U = Matrix(1.0I,n,n)
U[1,2:end] = -c/α
return L*A*U
end
myreduce1(A)
##
# This function saves the elimination structure as a set of matrices.
function myreduce_all(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)
for i=1:n-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
end
L,U,d = myreduce_all(A)
L*Diagonal(d)*U - A