function x = lrsolve(A, b)
% Given a matrix  A  and a column vector  b, the function
%   constructs matrices  L, R such that
%       L is lower triangular
%       R is upper triangular
%       A = L * R  (up to roundoff)
%       diag(L) = [1; 1; ...; 1]
%       L, R have minimal generic size
%   and returns an approximate solution x to A x = b
%   using  L, R  by backsubstitution 
%
%   Author: 
%   Date:
    
% Check if A is square
assert(all(size(A) == size(A')));
n = size(A, 2);
    
% Make default vector x and matrices L, R
L = eye(size(A));
R = A;
x = zeros(n, 1);
  
% ...
% TODO: Compute  L, R and c = L\b  explicitly 
% ...

% Check integrity of the LR decomposition
assert(TestLR(A, L, R));
% and that indeed  c = L\b
assert(norm(L*c - b, 'inf') <= 1e-10 * norm(b,'inf'));

% ...
% TODO: Compute x=R\c by backsubstitution
% ...

% Check integrity of the solution
if (norm(A*x-b, 'inf') > 1e-9 * norm(b,'inf'))
    warning('Solution tolerance not met');
end
end

function ok = TestLR(A, L, R)
    ok = false;
    if (~all(all(L == tril(L))))
        warning('L must be lower triangular');
    elseif (max(abs(diag(L) - 1)))
        warning('L(i,i) must be all ones');
    elseif (~all(all(R == triu(R))))
        warning('R must be upper triangular');
    elseif ((norm(A-L*R, 'inf') > 1e-8 * norm(A, 'inf')))
        warning('L*R must approximate A');
    else
        ok = true;
    end
end
