In [1]:
## 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)

norm(xtrue - xmy) = 0.0
Out[1]:
0.0
In [2]:
##
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)

norm(xtrue - xmy) = 889.1784197001251
Out[2]:
889.1784197001251
In [3]:
##
In [4]:
##
## 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)

newrow = 2
norm(xtrue - xmy) = 0.0
Out[4]:
0.0
In [5]:
##
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,:])

4×4 Array{Float64,2}:
 0.0   0.0           0.0           0.0        
 0.0   0.0          -8.67362e-18   0.0        
 0.0  -6.93889e-18   0.0          -1.11022e-16
 0.0   0.0           0.0          -1.11022e-16
In [6]:
##

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,:]
Out[6]:
2×2 Array{Float64,2}:
 0.0  0.0
 0.0  0.0
In [7]:
##
A = randn(4,4)
L,U,d,p = myreduce_all_pivot(A)
norm(L*Diagonal(d)*U - A[p,:])
Out[7]:
4.554929335617108e-16
In [8]:
##
A = [1e-18 1; 1e-18 1]
b = [5; 6]
L,U,d,p = myreduce_all_pivot(A)

the system is singular

Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] myreduce_all_pivot(::Array{Float64,2}) at ./In[5]:19
 [3] top-level scope at In[8]:4
In [9]:
##
lu(A)
Out[9]:
LU{Float64,Array{Float64,2}}
L factor:
2×2 Array{Float64,2}:
 1.0  0.0
 1.0  1.0
U factor:
2×2 Array{Float64,2}:
 1.0e-18  1.0        
 0.0      1.11022e-16
In [10]:
## We can re-use an LU factorization
A = randn(4,4)
F = lu(A)
F\randn(4)
F\randn(4)
Out[10]:
4-element Array{Float64,1}:
 -0.6036838717289073 
 -1.042742001925395  
 -0.03644492583425636
 -0.4285022544556376
In [11]:
## Your experiment -- play around with how long A\b takes vs. F\b