% CP2p9.m
% Computer problem 2.9 on p.102 Heath
% K. Ming Leung, 02/23/04

clear all; format short e;
xExact = [1; 1];    % The exact solution
n = length(xExact); % Size of the problem 

fprintf(' s   Error w/o Pivoting    Error Iterative Refinement    Error w. Pivoting\n');
fprintf('============================================================================\n');

for k = 1:10;
    % Set up the linear system for each k
    s = 10^(-2*k);      % epsilon
    A = [ s 1; 1 1 ];   % Matrix A
    b = [ 1+s; 2];      % The right-hand-side vector b
    [L,U] = LU_factor(A,n);         % LU factorization of A
    % Solve the linear system:
    y = forwardSubstitution(L,b,n);
    x = backSubstitution(U,y,n);        % Solution using LU w/o pivoting
    
    % Perform one step of iterative refinement:
    r = b - A*x;    % Find the residual
    y = forwardSubstitution(L,r,n);
    s1 = backSubstitution(U,y,n);
    x1 = x+s1;      % Solution using 1 step iterative refinement
    
    % Solution using Gaussian elimination with pivoting
    [xPivot,P,L,U] = LU_rowPivot(A,b,n);    % Solution using LU w/ pivoting
    
    % Compute the errors:
    err = norm(x-xExact);           % Error w/o using pivoting 
    err1 = norm(x1-xExact);         % Error after 1 step of iterative refinement
    errPivot = norm(xPivot-xExact); % Error using pivoting
    fprintf('%8.0e %16.4e %20.4e %20.4e\n', s, err, err1, errPivot);
end


