## 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)

norm(xtrue - xmy) = 4.768226942950618e-15

4.768226942950618e-15

##
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)

4×4 Array{Float64,2}:
 0.0284233     2.22045e-16    0.0        0.0     
 0.0        -148.923        -35.2488   -67.3763  
 0.0         -59.3247       -12.3848   -25.8155  
 0.0          -0.606229      -1.20104    0.209151

##
# 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

4×4 Array{Float64,2}:
 0.0  -2.22045e-16  0.0           0.0        
 0.0  -2.69784e-14  3.55271e-15  -2.13163e-14
 0.0  -1.42109e-14  2.77556e-15   1.33227e-15
 0.0   0.0          0.0          -2.22045e-16