%
%   This matlab script
%   is designed to demonstrate the validity of the bound.
%
%   norm(x1 -x)/norm(x) .le.  (cond(A)/cos(theta))*norm(b1 - b)/norm(b)
%
%   where Ax =  b ,    Ax1 =  b1.
%             2             2
%
%   i.e. x and x1 solve least squares problem for rhs b and rhs b1
%   respectively
%
%   Here the angle theta is increased  while Pb (Orth. projection of b)
%   remains constant.   
%
%        b1 =  b + tau*w    , where tau = eps*norm(b)
%
%   and norm(w) = 1
%
%   D.C. Sorensen  10 Oct 01
%
 format short e
 m = 100;
 n = 30;
 d = 20;
 k = 10;
 k1 = 5;

 Xout = zeros(k1,5);
 A = randn(m,n);
 [U,S,V] = svd(A,0);
%
%  prepare rhs.  b  
%  
    b = U(:,1:d)*randn(d,1);
%
%
 z = randn(m,1);
 z = z - U*(U'*z);
 z = z - U*(U'*z);
 z = z/norm(z);

 b0 = b/norm(b);
%
%   Note bo is Pb for all of the constructed right hand sides
%
%   while z is a unit vector orthogonal to Range(A)
%
 S0 = S;
 fac = 10;
 for itr = 1:k, 
    
    b = b0 + fac*z;
%
%   construct RHS b  
%
%   Note  Pz = 0     so  Pb = b0
%
%   and  norm(b) = sqrt(1 + fac^2)
%
%   since b0'*z = 0  and both are unit vectors.
%
    fac = fac*10;

    tau = norm(b)*eps;
    w = randn(m,1); w = w/norm(w);

    s = diag(S0);
    b1 = b + tau*w;
    for j = 1:k1,
%
%       increase condition by factor of 10
%
        s(d+1:n) = s(d+1:n)/10; 
        S = diag(s);
%
%       form new matrix and solve for x1 and x
%
        A1 = U*S*V';
        x1 = A1\b1;
        x = A1\b;
        condA1 = cond(A1);
        xerr = norm(x1 - x)/norm(x);
        berr1 = norm(b1 - b)/norm(b);
        berr2 = norm(b1 - b)/norm(b0);
        Xout(j,:) =  [ xerr condA1*berr2 condA1*norm(b)/norm(b0) condA1*berr1  condA1 ];
    end

    disp('   ')
    disp('  estimate1 =  cond(A) norm(b1 - b)/norm(b)' )
    disp('  estimate2 =  (cond(A)/cos(theta)) norm(b1 - b)/norm(b)' )


    disp('   ')
    disp('     x-err       estimate2   cond(A1)/    estimate1    cond(A1)   ')
    disp('                             cos(theta)                           ')
    disp(Xout)
    disp('   ')
    disp('    norm(b)    1/cos(theta)   ')
    disp([ norm(b)  norm(b)/norm(b0)])
 pause
 end