Solving Linear Equations

Serial Algorithm

First, we construct a linear equations system $Ax=b$.

n = 5
A = rand(-10.0:10.0, (n, n))

x = rand(-10.0:10.0, n)
b = A * x

The code in the following cell converts the problem $Ax=b$ to the upper triangular equation system $Ux=y$.

function convert_to_upper_triangular(A,b)
    # Upper Triangularization: convert Ax=b to Ux=y
    for k in 1:n
        for j in k+1:n
            # Divide by pivot
            A[k,j] = A[k,j] / A[k,k]
        end
        b[k] = b[k] / A[k,k]
        A[k,k] = 1
        # Substract lower rows
        for i in k+1:n 
            for j in k+1:n
                A[i,j]=A[i,j] - A[i,k] * A[k,j]
            end
            b[i] = b[i] - A[i,k] * b[k]
            A[i,k] = 0
        end
    end
    return A, b #U,y
end

The function in the following cell solves the upper triangular equation system using backwards substitution.

function solve_upper_triangular(U,y)
    n = size(U,1)
    for step in reverse(1:n)
        if U[step,step] == 0
            if y[step] != 0
                return "No solution"
            else
                return "Infinity solutions"
            end
        else
        # Backwards substitution
            y[step] = y[step] / U[step,step]
        end
        for row in reverse(1:step-1)
            y[row] -= U[row,step] * y[step]
        end
    end
    return y 
end

U,y = convert_to_upper_triangular(A,b)
sol = solve_upper_triangular(U,y)

We can test if the obtained solution is correct using @test:

using Test
@test sol ≈ x