% script file: testerrors.m
%
% test two inequalities:
%    ||b - A*xx||/(||A||*||xx||) <= rho*beta^(-t)
% and
%    ||x - xx||/||xx|| <= rho*cond(A)*beta^(-t)
% where xx is the solution computed by the Gaussian elimination
% with partial pivoting, x is the exact solution, beta is the
% base (=2), and t is the precision (=53).
% The first inequality shows that the relative residual (left side)
% can be used as an estimate for the backward error (right side).
% The second inequality shows that the actual (forward) error
% is bounded by the the product of cond(A) and the backward error.
%
% Dependency: decomp.m, solve.m

n = input('order of the matrix: ');

% Generate an n-by-n matrix whose elements are random integers between
% -10 and 10, so they are exact floating-point numbers
A = ceil(10*(2*rand(n,n) - ones(n,n)));
% make A singular
A(:,3) = 2*A(:,2) - A(:,1);

nrmA = norm(A, 1);

cnd = input('desired condition number (at least 1): ');
k = ceil(log(cnd/nrmA)/log(2));
% perturbation, 2^(-k), an exact floating-point number
pert = 1.0;
if (k > 0)
    for i=1:k
        pert = pert/2;
    end
end
% perturb A
A = A + pert*eye(n,n);

% construct the right side b so that the solution is ones(n,1)
b = A(:,1);
for j=2:n
    b = b + A(:,j);
end

AA = A;		% save A

[A, rcondA, pvt] = decomp(A);
if (1.0 + rcondA == 1.0)
    display('singularity is detected');
    return
end

x = b;
x = solve(A, x, pvt);

nrmx = norm(x, 1);
res = norm(b - AA*x, 1);
errx = norm(ones(n,1) - x, 1);

% compute rho
rho = res/(eps*nrmA*nrmx);
rho,

% compare the two sides of the second inequality
lside2 = errx/nrmx;
rside2 = rho*eps/rcondA;
lside2, rside2,