function [A0, resid1, resid2] = resid_LUppUnstableExamplePerturbed(n,pert_type,fig)
%
% compute residual of solution of Ax=b, where b is random,
% on perturbations of the extremely rare matrix where GEPP 
% (LU factorization with partial pivoting) is unstable (see A&G, p. 109).
% n is the order of the matrix.
% pert_type = 1 means additive perturbations (making changes even to the zeros)
% pert_type = 2 means multiplicative perturbations (leaving zeros intact).
% fig specifies the figure number: plots are done in figures fig,fig+1
% The plot shows that as long as the perturbation is 10^{-16} or bigger, the
% residual is very small, as usual for GEPP, but for zero or tiny 
% perturbations and n >= 20, GEPP is unstable and the residual is not small. 
% Such a situation is VERY RARE.
% Note that the CONDITION NUMBER of A is NOT large, 
% even when there is no perturbation, but the condition numbers of L and U
% are large, as shown in Figure 2, for sufficiently tiny perturbations, 
% again indicating the instability of GEPP on this rare example.
%
fsize = 18; % font size for plots
msize = 10; % marker size for plots
b = randn(n,1); % ones(n,1) is too special: do not get large residuals
A0 = LUppUnstableExample(n);
for j=1:30
    if j == 20
        keyboard
    end
    pert(j) = 10^(-j);
    if pert_type == 1 
        A = A0 + pert(j)*randn(n); % changes the zeros for all levels of pert
    elseif pert_type == 2
        A = A0 .* (1 + pert(j)*randn(n)); % leaves zeros of A0 intact
    else
        error('invalid pert_type')
    end
    x1 = A\b;      % uses partial pivoting (row interchanges)
    [L,U] = lu(A); % same algorithm another way
    x2 = U\(L\b);  % should be the same as x1, though not always identical
    resid1(j) = norm(A*x1 - b); % not much different from relative residual in this case
    resid2(j) = norm(A*x2 - b); 
    condA(j) = cond(A)
    condL(j) = cond(L);
    condU(j) = cond(U);
end
figure(fig)
clf
loglog(pert,resid1,'rx',pert,resid2,'kd','MarkerSize',msize)
legend('resid1','resid2')
xlabel('perturbation made to worst case example','FontSize',fsize)
ylabel('residual','FontSize',fsize)
if pert_type == 1
    title(sprintf('Residual from GEPP with n=%d (add. perts.)',n))
else
    title(sprintf('Residual from GEPP with n=%d (mult. perts.)',n))
end
%
set(gca,'FontSize',fsize)
figure(fig+1)
clf
loglog(pert,condA,'k+',pert,condL,'bo',pert,condU,'m*','MarkerSize',msize)
legend('cond(A)','cond(L)','cond(U)')
xlabel('perturbation made to worst case example','FontSize',fsize)
ylabel('condition number','FontSize',fsize)
if pert_type == 1
    title(sprintf('Cond. Nos. of A and L,U (GEPP), n=%d (add. perts.)',n))
else
        title(sprintf('Cond. Nos. of A and L,U (GEPP), n=%d (mult. perts.)',n))
end
set(gca,'FontSize',fsize)
figure(fig),shg,figure(fig+1),shg

    