using LinearAlgebra

# perform Gaussian elimination of A without row swaps, returning U,
# while printing a message for each elimination step.
function print_gauss(A)
    m = size(A,1) # number of rows
    U = copy!(similar(A, typeof(inv(A[1,1]))), A)
    for j = 1:m   # loop over m columns
        for i = j+1:m   # loop over rows below the pivot row j
            # subtract a multiple of the pivot row (j)
            # from the current row (i) to cancel U[i,j] = Uᵢⱼ:
            ℓᵢⱼ = U[i,j]/U[j,j]
            println("subtracting $ℓᵢⱼ × (row $j) from (row $i)")
            U[i,:] = U[i,:] - U[j,:] * ℓᵢⱼ
            U[i,j] = 0 # store exact zero to compensate for roundoff errors
        end
    end
    return U
end

A = [4  -2  -7  -4  -8
     9  -6  -6  -1  -5
    -2  -9   3  -5   2
     9   7  -9   5  -8
    -1   6  -3   9   6]

print_gauss(A)

L, U = lu(A, NoPivot())
U

L

L[3,1] * U[1,:] + L[3,2] * U[2,:] + U[3,:]

L

 

L = [ 1  0  0 
     +2  1  0
     -3 +7  1 ]
U = [ 1  2  0
      0  1  1
      0  0 -8 ]
L * U

inv(L)

A = [4  -2  -7  -4  -8
     9  -6  -6  -1  -5
    -2  -9   3  -5   2
     9   7  -9   5  -8
    -1   6  -3   9   6] # our random 5×5 matrix from before
b = rand(-9:9, 5)

_, U_and_c = lu([A b], NoPivot()) # eliminate augmented matrix (without row swaps)
U = UpperTriangular(U_and_c[:, 1:end-1]) # all but last column is U
c = U_and_c[:, end] # last column is c

[U\c  A\b] # print them side by side

L, U = lu(A, NoPivot()) # Gaussian elimination without row swaps
c = L \ b # solve Lc = b for c

c = similar(b, Float64)
for i = 1:length(b)
    print("c[$i] = b[$i]")
    c[i] = b[i]
    for j = 1:i-1
        print("- $(L[i,j]) * c[$j]")
        c[i] = c[i] - L[i,j] * c[j]
    end
    println(" = ", c[i])
end
c

b = [5, 15, -4]
L = [ 1  0  0 
     +2  1  0
     -3 +7  1 ]
U = [ 1  2  0
      0  1  1
      0  0 -8 ]
A = L * U

c = L \ b

U \ c

A \ b

B = [b 2b 3b] # the solutions should be just x = (1,2,3), 2x, and 3x

A \ B

A^-1 * B

LU = lu(A)

[LU\b  A\b] # print them side by side

LU \ 2b

LU \ B

A = [ 1 3  1
      1 3 -1
      3 11 6 ]

b = [9, 1, 35]

A \ b

P2 = [1 0 0
      0 0 1
      0 1 0]
E1 = [1 0 0
     -1 1 0
     -3 0 1]
P2 * E1

L = [1 0 0
     3 1 0
     1 0 1]
U = [1 3 1
     0 2 3
     0 0 -2]
L*U

# matrix from earlier example, does not *require* row swaps

A = [ 1  2   0
      2  5   1
     -3  1  -1]

L, U, p = lu(A)
L

U

p

A[p,:] # A with the rows in order p = PA

L*U

# construct a permutation matrix P from the permutation vector p
permutation_matrix(p) = I(length(p))[p,:] # just reorder the rows of I (returned by I(n))

P = permutation_matrix(p)

P * A

A \ b

b[p] # permute b (equivalent to Pb)

c = L \ b[p] # solve Lc = Pb = b[p]
x = U \ c # solve Ux = c

LU = lu(A)

LU \ b # applies permutation, then forward-sub for L, then back-sub for U