% generate a random matrix of size dimension x dimension
dimension = 1500;
A = rand(dimension);

% iterate a number of times so the timing is meaningful
n_iterations = 100;

fprintf('\n======================\n');
fprintf('Use backslash %d times\n\n',n_iterations);
tic                         % start timer
for i=1:n_iterations
    b = rand(dimension,1);  % generate a random rhs
    x = A\b;                % LU factorize and then back/forward sub.

    if mod(i,10)==0
        fprintf('finished step %d of %d\n',i,n_iterations);
    end
end
toc                         % end timer

fprintf('\nLU factorize once then back/forward sub. %d times\n\n',n_iterations);
tic                         % start timer
[L,U] = lu(A);              % lu factorize once
for i=1:n_iterations

    b = rand(dimension,1);  % generate a random rhs
    y = L\b;                % forward substitution to find y
    x = U\y;                % back substitution to find x
    
    if mod(i,10)==0
        fprintf('finished step %d of %d\n',i,n_iterations);
    end
end
toc                         % end timer

fprintf('\n======================\n');
fprintf('\nCompare computation time for inv(A) and lu(A)\n');
fprintf('\nTime for inv(A):\n');
tic;
A_inv = inv(A); % compute A^-1
toc;
fprintf('\nTime for lu(A):\n');
tic;
[L,U] = lu(A);  % compute LU factorization
toc;

fprintf('\n======================\n');
