function resid_LUppUnstableExamplePerturbed_demo
% 
% Show the residuals of GEPP when applied to perturbations of the well
% known worst case example for GEPP, along with condition numbers of L, U.
% We very different behavior for additive perturbations (which change
% the zero entries of the unperturbed matrix no matter how small the perturbation
% is) and multiplicative perturbations (which leave the zeros intact, and 
% have no effect for pert smaller than about 10^{-16}).
%
% In the MULTIPLICATIVE perturbation case, we see that for effectively zero
% perturbations, that is pert smaller than or equal to 10^{-17} so that 1 + pert  
% rounds to one, the perturbation has no effect at all and so the 
% GEPP residual is big (the size depends on n), the condition number of L is big, 
% and the condition number of U is much bigger still, overflowing to inf for about 
% n = 80. But for small multiplicative perturbations for which 1 + pert does *not* 
% round to one, that is 10^{-16} and bigger, the effect dramatically disappears, 
% because the perturbations result inpivoting, thus breaking the properties of the 
% bad example, with no difference whether pert is 10^{-16} or 10^{-1}, showing how 
% sensitive and special this worst-case example is!
%
% The results in the ADDITIVE perturbation case are very different. Even
% tiny perturbations to the zero entries in the upper triangle of the  matrix, 
% such as 10^{-20}, get slowly amplified during the elimination and eventually 
% they become large enough to impact the diagonal and subdiagonal in U, resulting
% in pivoting which breaks the bad example, so we see essentially a smooth 
% transition from the effect with perturbations ranging from 10^{-30} to 10^{-16},
% since smaller perturbations take longer to grow large enough to impact
% the pivoting.  For the same reason, once the bad example is broken by the
% pivoting, cond(U) is no longer enormous, even for n=100. In fact, now 
% cond(L) and cond(U) have similar magnitudes, unlike in the multiplicative 
% perturbation case where the example does not break.
% However, because the perturbations with size >= 10^{-16} introduce pivoting
% from the beginning of the elimination, the bad example breaks immediately,
% and so there is almost no difference between the effect for 
% pert = 10^{-16} and 10^{-1}.

% In the residual figure, the only difference between resid1 and resid2 is
% that the first comes from x1=A\b, while the second is from 
% [L,U]=A; x2=U\(L\b), which one expects to give the same result, 
% but it is not exactly the same.
%
% Interesting to run this for:
%  n = 30: this is already big enough to see a dramatic effect
%  n = 50: much more dramatic effect
%  n = 70: now, for multiplicative perturbations only, cond(U) is huge
%          for effectively zero perturbations
%  n = 100: now, for multiplicative perturbations only, cond(U) is inf 
%         for effectively zero perturbations and does not appear in the
%         plot, so the gigantic scale of the vertical axis of the bottom
%         right plot drops dramatically.
%
% Michael Overton, March 1, 2017.
%
n = input('enter n ')
% first additive perturbations: show results in Figures 1, 2
resid_LUppUnstableExamplePerturbed(n,1,1);
% now multiplicative perturbations: show results in Figures 3, 4
resid_LUppUnstableExamplePerturbed(n,2,3);