function x = solveLrPivot(A, P, b)
    len = length(b)
    z = zeros(len,1)
    
    % solve Lz = b
    for i = 1:len
        z(i) = b(i)
        
        for j = 1:(i-1)
            z(i) = z(i) - z(j) * A(P(i), j)
        end

        % no division necessary as L is unipotent
    end

    x = zeros(len,1)

    % solve Rx = z
    for i = len:-1:1
        x(i) = z(i)

        for j = len:-1:(i+1)
            x(i) = x(i) - x(j) * A(P(i), j)
        end
        
        x(i) = x(i) / A(P(i), i)
    end

    x_ret = x
end