%  Solve a random linear system using the LU-decomposition
%  in a decimal floating point system using a mantissa length 
%  of m
%
%  Matthias Heinkenschloss
%  Department of Computational and Applied Mathematics
%  Rice University
%  Feb 12, 2001
%

clear;
% 
fprintf(1, 'Solve a random linear system using the LU-decomposition \n' )
fprintf(1, 'in a decimal floating point system using a mantissa length of m \n')

% Use floating point arithmetic
arith = 'r';

% Set the dimension of the system and compute the matrix
% A and the exact solution x using the random number generator.

n    = 5;
A_ex = 10*rand(n);
x_ex = 10*rand(n,1);

% Compute the corresponding right hand side
b_ex = A_ex * x_ex;

fprintf(1, 'The condition number of A is  %12.6e \n', cond(A_ex) )


% compute the solution of the linear system using 
% the LU-decomposition with pivoting and various 
% mantissa length

% print header for absolute and relative errors
fprintf(1, '       Absolute      Relative  \n')
fprintf(1, '  m    Error         Error  \n')

for m = 1:20
    A = A_ex;
    b = b_ex;
    [A, ipivt, iflag] = lu_pp_fl( A, m, arith );
    if( iflag == 0 ) 
        [b, iflag] = lu_pp_sl_fl( A, b, ipivt, m, arith );
    end

%   if a solution is computed, compute the relative and
%   absolute error between computed and exact solution.

    if( iflag > 0 ) 
%       matrix is detected to be singular. No solution 
%       is computed, set error equal to infinity
        abserr(m) = inf;
        relerr(m) = inf; 
    else
        abserr(m) = norm(x_ex - b);
        relerr(m) = abserr(m) / norm(x_ex);
    end
    
%   print absolute and relative errors
    fprintf(1, '%3i  %12.6e  %12.6e \n', m, abserr(m), relerr(m) );

end
%   plot absolute and relative errors
semilogy((1:20),abserr,'r'); hold on;
semilogy((1:20),relerr,'b'); hold off;
title( 'Absolute (red) and relative (blue) error in solution' );
xlabel(' Mantissa length m ');