## Here's partial pivoting
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!
    α = abs(A[1,1])
    newrow = 1
    for j=2:n
      if abs(A[j,1]) > α
        α = abs(A[j,1])
        newrow = j
        break
      end
    end
    if α <= 10*eps(1.0) # need a better tolerance here... but this will do
      error("the system is singular, with a size = $(size(A)) system")
    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)

UndefVarError: `norm` not defined

Stacktrace:
 [1] macro expansion
   @ show.jl:1181 [inlined]
 [2] top-level scope
   @ In[1]:43

## Let's test this on a singular system
using StableRNGs
n = 5
r = 2
A = rand(n,n) 
A[:,1] .= 0 # make it singular
b = ones(5) 

xmy = solve1_pivot!(copy(A),b) # this changes A.

the system is singular, with a size = (5, 5) system

Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] solve1_pivot!(A::Matrix{Float64}, b::Vector{Float64})
   @ Main ./In[1]:21
 [3] top-level scope
   @ In[2]:9

## 
# This fails on the first step step because that
# column is zero. But the rank isn't 0! It's n-1.

## Let's modify this for complete pivoting 
function solve1_pivot_complete!(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
    if A[1] <= 10*eps(1.0)
      error("the system is singular, with a size = $(size(A)) system")
    end
    return [b[1]/A[1]]
  else
    # let's find the equation and variable to eliminate
    eqvarpair  = argmax(abs.(A)) # I don't like this as it makes a copy
    newrow = eqvarpair[1]
    newvar = eqvarpair[2]
    if abs(A[newrow,newvar]) <= 10*eps(1.0) # need a better tolerance here... but this will do
      error("the system is singular, with a size = $(size(A)) system")
    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
    if newvar != 1
      tmp = A[:,1]
      A[:,1] .= A[:,newvar]
      A[:,newvar] .= tmp
    end
    D = A[2:end,2:end]
    c = A[1,2:end]
    d = A[2:end,1]
    α = A[1,1]
    y = solve1_pivot_complete!(D-d*c'/α, b[2:end]-b[1]/α*d)
    γ = (b[1] - c'*y)/α
    x = pushfirst!(y,γ)
    # swap variables back
    x[1], x[newvar] = x[newvar], x[1]
    return x 
  end
end

A = [0 1.0; 1.0 1.0]
b = [5; 6.0]
xtrue = A\b
xmy = solve1_pivot_complete!(A,b) # this changes A.
@show norm(xtrue-xmy)

UndefVarError: `norm` not defined

Stacktrace:
 [1] macro expansion
   @ show.jl:1181 [inlined]
 [2] top-level scope
   @ In[4]:48

## Difference between partial and complete pivoting
A = rand(100,100)
x = ones(100)
b = A*x

x1 = solve1_pivot!(copy(A),copy(b))
x2 = solve1_pivot_complete!(copy(A),copy(b))
z = A\b 

norm(x1-ones(100)), norm(x2-ones(100)), norm(z-ones(100))

UndefVarError: `norm` not defined

Stacktrace:
 [1] top-level scope
   @ In[5]:10

##
A = rand(5,5)
A[:,1] .= 0
b = ones(5)
x2 = solve1_pivot_complete!(copy(A),copy(b))
# This now fails on the last step when we have a 1,1
# system, which suggests the rank is n-1. 

the system is singular, with a size = (1, 1) system

Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] solve1_pivot_complete!(A::Matrix{Float64}, b::Vector{Float64})
   @ Main ./In[4]:8
 [3] solve1_pivot_complete!(A::Matrix{Float64}, b::Vector{Float64}) (repeats 4 times)
   @ Main ./In[4]:35
 [4] top-level scope
   @ In[6]:5

## This now fails on the