function [x,A,l] = lusolve(A,b)
% x = lusolve(A,b) solve the linear system Ax=b using Gaussian elimination
% with scaled partial pivoting.
%
% [x,LU,l] = lusolve(A,b) in addition returns a matrix LU which contains 
% the LU-decomposition of A together with a vector l containing the row
% exchanges that have taken place.
%
% INPUT: A (Coefficient matrix, assumed invertible.)
%        b (Right hand side.)
%
% OUTPUT: x (Solution vector.)
%         A (The resulting matrix from the forward elimination procedure.)
%         l (Index vector. If the rows of A are viewed in the order given
%            in l, the upper triangular part holds the transformed 
%            coefficient matrix and the strictly lower triangular part
%            holds the multipliers used in the elimination procedure.)
%         
%


%% initialization and preparation
% check if the dimensions of the problem somehow match and return an error
% message if not
n = length(b);
dimA = size(A);
if(length(dimA) > 2 || dimA(1)~= n || dimA(2) ~= n)
    error('The dimensions of A and b do not match.')
end

% initialize the solution x as zero vector of the correct size
x = zeros(size(b));

% initialize the index vector l
l = 1:n;

% initialize the vector of ratios needed for determining the pivot element
r = zeros(size(b));

% compute the scale vector
s = max(abs(A),[],2);

%% forward elimination
for k=1:n-1
    % compute the ratios between (relevant) entries in column scale and
    % respective scale.
    r(k:end) = abs(A(l(k:end),k))./s(l(k:end));
    % store the position where the maximum is attained in the variable pos.
    % Note that the position is calculated in the truncated vector
    % containing only the entries from k to n. Thus the actual position is
    % obtained by adding k-1.
    [~, pos] = max(r(k:end));
    % If necessary simulate a row permutation.
    if pos > 1
        l([k pos+k-1]) = l([pos+k-1 k]);
    end
    % compute the multipliers and store them in the matrix A
    A(l(k+1:end),k) = A(l(k+1:end),k)/A(l(k),k);
    
    % subtract multiples of line k from the next lines
    A(l(k+1:end),k+1:end) = A(l(k+1:end),k+1:end)-A(l(k+1:end),k)*A(l(k),k+1:end);
end

%% back substitution
% process the right hand side in the same way as the matrix A
for k=1:n-1
    b(l(k+1:end)) = b(l(k+1:end))-A(l(k+1:end),k)*b(l(k));
end

% finally compute the solution of the equation
for i=n:-1:1
    x(i) = (b(l(i))-A(l(i),i+1:end)*x(i+1:end))/A(l(i),i);
end