%
%   This matlab script
%   is designed to demonstrate the validity of the bound.
%
%       norm(x1 -x)/norm(x) .le. cond(A)*norm(b1 - b)/norm(b)
%
%   where Ax = b , Ax1 = b1.
%
%   Also demonstrates small residual does not imply accurate answer.
%   You MUST know condition number.
%
%   D.C. Sorensen  9/27/99
%
%   modified 10 Oct 01 (DCS)

 format short e
 n = 100;
 m = 90;
 k = 13;
 Xout = zeros(k,5);
 A = randn(n);
 [U,S,V] = svd(A);
%
%  prepare rhs.  b  
%  
    b = U(:,1:m)*randn(m,1);
%
%
 z = randn(n,1);
 tau = 100*eps;
 x = A\b;
 kk = cond(A);
 s = diag(S);
 b1 = b + tau*z;
 for j = 1:k,
%
%    increase condition by factor of 10
%
     s(m+1:100) = s(m+1:100)/10; 
     S = diag(s);
%
%    form new matrix and solve for x1 and x
%
     A1 = U*S*V';
     x1 = A1\b1;
     x = A1\b;
     resid = norm(b - A1*x1);
     condA1 = cond(A1);
     xerr = norm(x1 - x)/norm(x);
     berr = norm(b1 - b)/norm(b);
     Xout(j,:) =  [ xerr  condA1*berr  round(-log10(xerr))  condA1 resid];
 end

 disp('   ')
 disp('     x-err       estimate      digits      cond(A1)       resid ')
 disp(Xout)