% Produce figures 2.1, 2.2 and 2.3 in text
% Written by James Demmel, Sept 16, 1995
%                Modified, Aug  18, 1996
%
% 
% Matrix dimensions to try
if (1==0)
nn = [(5:5:100)];
% Number of random matrices of each dimension
iter = 5*ones(size(nn));
%
% eps = 2^(-52) (built into Matlab)
% macheps = eps/2 (because we round correctly)
macheps = eps / 2;
%
% Quantities computed and plotted for each test matrix:
%
% cnd1 = condition number norm(A)*norm(inv(A))
% cnd2 = Hager's condition number estimate 
% cnd3 = relative condition number, using x from GEPP
% cnd4 = relative condition number, using x from GECP
%
% gpp = pivot growth factor for Gaussian Elimination
%       with Partial Pivoting (GEPP) 
%     = max_{ij} |U(i,j)| / max_{ij} |A(i,j)|
% bndpp1, bndpp2, bndpp3 store backward errors bounds for GEPP
% bndpp1, the coarsest estimate,   = macheps*3*n^3*gpp
% bndpp2, a better estimate,       = macheps*3*n*|| |L| |U| || / || A ||
% bndpp3, the true backward error  = macheps*|| A*x-b || / ( ||A|| ||x|| )
% errpp = norm(computed x - true x)/norm(computed x)
% bndpp4 = bndpp3*cnd2 = computed bound on errpp
% bndpp5 = componentwise relative backward error 
% bndpp6 = bndpp5*cnd3 = another computed bound on errpp
% bndpp7 = norm(abs(inv(A))*(abs(r)))/norm(computed x)
%        = another computed bound on errpp
% bndpp8 = norm(abs(inv(A))*
%              (abs(r) + (n+1)*eps*(abs(A)*abs(x)+abs(b))))/norm(computed x)
%        = another computed bound on errpp
% errppitref = norm(computed x from 1 step of working precision 
%                   iterative refinement  - true x)/norm(computed x)
% bndpp9 = norm(abs(inv(A))*(abs(r)))/norm(computed x)
%        = computed bound on errppitref
% bndpp10= norm(abs(inv(A))*
%              (abs(r) + (n+1)*eps*(abs(A)*abs(x)+abs(b))))/norm(computed x)
%        = computed bound on errppitref
% bndpp11= componentwise relative backward error in solution after 
%          1 step of iterative refinement
% normrpp     = norm(residual of xpp)
% normrppitref= norm(residual of xppitref)
%   (all norms are the infinity-norm)
% 
%
% gcp, bndcp1, bndcp2, bndcp3, bndcp4, bndcp5, bndcp6
%      bndcp7, bndcp8, bndcp9, bndcp10, bndcp11, normrcp, normrcpiterf
% are the corresponding quantities for 
% Gaussian Elimination with Complete Pivoting (GECP)
%
% Initialize variables
%
gpp    = 0;
bndpp1 = 0;
bndpp2 = 0;
bndpp3 = 0;
bndpp4 = 0;
bndpp5 = 0;
bndpp6 = 0;
bndpp7 = 0;
bndpp8 = 0;
bndpp9 = 0;
bndpp10= 0;
bndpp11= 0;
errpp  = 0;
normrpp= 0;
errppitref   = 0;
normrppitref = 0;
gcp    = 0;
bndcp1 = 0;
bndcp2 = 0;
bndcp3 = 0;
bndcp4 = 0;
bndcp5 = 0;
bndcp6 = 0;
bndcp7 = 0;
bndcp8 = 0;
bndcp9 = 0;
bndcp10= 0;
bndcp11= 0;
errcp  = 0;
normrcp= 0;
errcpitref   = 0;
normrcpitref = 0;
cnd1 = 0;
cnd2 = 0;
cnd3 = 0;
cnd4 = 0;
dim = 0;
%
% Begin loop over test matrices
%
j=0;
k=0;
%
% For each dimension
%
for n = nn;
   k = k+1;
%
%  For each random sample of that dimension
%
   for i=1:iter(k);
      j = j+1;
%
%     Generate a random matrix
%     This may be modified to illustrate different points.
%     The following choice illustrates the utility of one
%     step of iterative refinement in the same precision
%     as the original data, and is used in Example 2.5
%
      a = eye(n) + .0000001*randn(n,n);
      dd = diag((10*ones(1,n)).^(14*(0:n-1)/(n-1)));
      a = dd*a;
%
%     Generate a random right hand side
%
      x = randn(n,1);
      b = a*x;
%
%     Solve using Gaussian Elimination with Partial Pivoting (GEPP)
%
      [lpp,upp,ppp]=lu(a);
      xpp = upp\(lpp\(ppp*b));
%
%     Solve using GEPP with 1 step of iterative refinement
%
      rppt = a*xpp-b;
      dpp = upp\(lpp\(ppp*rppt));
      xppitref = xpp - dpp;
%
%     Compute pivot growth for GEPP
%
      gpp(j) = max(max(abs(upp)))/max(max(abs(a)));
%
%     Compute backward error bounds for GEPP
%
      bndpp1(j) = 3*n^3*gpp(j)*macheps;
      bndpp2(j) = 3*n*norm(abs(lpp)*abs(upp),'inf')/norm(a,'inf')*macheps;
      bndpp3(j) = norm(a*xpp-b,'inf')/(norm(a,'inf')*norm(xpp,'inf'));
      bndpp5(j) = max(abs(a*xpp-b)./(abs(a)*abs(xpp)+abs(b)));
      bndpp11(j)= max(abs(a*xppitref-b)./(abs(a)*abs(xppitref)+abs(b)));
%
%     Compute error for GEPP
%
      errpp(j) = norm(xpp-x,'inf')/norm(xpp,'inf');
%
%     Compute error for GEPP with 1 step of iterative refinement
%
      errppitref(j) = norm(xppitref-x,'inf')/norm(xpp,'inf');
%
%     Solve using Gaussian Elimination with Complete Pivoting (GECP)
%
      [pcprow,lcp,ucp,pcpcol,err]=gecp(a);
      xcp = pcpcol*(ucp\(lcp\(pcprow'*b)));
%
%     Solve using GECP with 1 step of iterative refinement
%
      rcpitreft = a*xcp-b;
      dcp = pcpcol*(ucp\(lcp\(pcprow'*rcpitreft)));
      xcpitref = xcp - dcp;
%
%     Compute pivot growth for GECP
%
      gcp(j) = max(max(abs(ucp)))/max(max(abs(a)));
%
%     Compute backward error bounds for GECP
%
      bndcp1(j) = 3*n^3*gcp(j)*macheps;
      bndcp2(j) = 3*n*norm(abs(lcp)*abs(ucp),'inf')/norm(a,'inf')*macheps;
      bndcp3(j) = norm(a*xcp-b,'inf')/(norm(a,'inf')*norm(xcp,'inf'));
      bndcp5(j) = max(abs(a*xcp-b)./(abs(a)*abs(xcp)+abs(b)));
      bndcp11(j)= max(abs(a*xcpitref-b)./(abs(a)*abs(xcpitref)+abs(b)));
%
%     Compute error for GECP
%
      errcp(j) = norm(xcp-x,'inf')/norm(xcp,'inf');
%
%     Compute error for GECP with 1 step of iterative refinement
%
      errcpitref(j) = norm(xcpitref-x,'inf')/norm(xcp,'inf');
%
%     Compute condition numbers
%
      inverse_a = pcpcol*(ucp\(lcp\(pcprow'*eye(n))));
      cnd1(j) = norm(a,'inf')*norm(inverse_a,'inf');
      [cnd2(j), v] = condest(a');
      cnd3(j) = norm(abs(inverse_a)*(abs(a)*abs(xpp)+abs(b)),'inf') ...
                /norm(xpp,'inf');
      cnd4(j) = norm(abs(inverse_a)*(abs(a)*abs(xcp)+abs(b)),'inf') ...
                /norm(xcp,'inf');
%
%     Compute forward error bounds for GEPP 
%
      rpp       = a*xpp - b;
      normrpp(j)= norm(rpp,'inf');
      bndpp4(j) = bndpp3(j) * cnd2(j);
      bndpp6(j) = bndpp5(j) * cnd3(j);
      bndpp7(j) = norm(abs(inverse_a)*abs(rpp),'inf')/norm(xpp,'inf');
      bndpp8(j) = norm(abs(inverse_a)*(abs(rpp)+ ...
                  (n+1)*eps*(abs(a)*abs(xpp)+abs(b))),'inf')/ ...
                  norm(xpp,'inf');
      rppitref  = a*xppitref - b;
      normrppitref(j)= norm(rppitref,'inf');
      bndpp9(j) = norm(abs(inverse_a)*abs(rppitref),'inf')/norm(xpp,'inf');
      bndpp10(j)= norm(abs(inverse_a)*(abs(rppitref)+ ...
                  (n+1)*eps*(abs(a)*abs(xppitref)+abs(b))),'inf')/ ...
                  norm(xpp,'inf');
%
%     Compute forward error bounds for GECP 
%
      rcp = a*xcp-b;
      normrcp(j)= norm(rcp,'inf');
      bndcp4(j) = bndcp3(j) * cnd2(j);
      bndcp6(j) = bndcp5(j) * cnd3(j);
      bndcp7(j) = norm(abs(inverse_a)*abs(rcp),'inf')/norm(xcp,'inf');
      bndcp8(j) = norm(abs(inverse_a)*(abs(rcp)+ ...
                  (n+1)*eps*(abs(a)*abs(xcp)+abs(b))),'inf')/ ...
                  norm(xcp,'inf');
      rcpitref  = a*xcpitref - b;
      normrcpitref(j)= norm(rcpitref,'inf');
      bndcp9(j) = norm(abs(inverse_a)*abs(rcpitref),'inf')/norm(xcp,'inf');
      bndcp10(j)= norm(abs(inverse_a)*(abs(rcpitref)+ ...
                  (n+1)*eps*(abs(a)*abs(xcpitref)+abs(b))),'inf')/ ...
                  norm(xcp,'inf');
%
%     Save dimension for plotting
%
      dim(j) = n;
    end
    disp(['just finished n = ',int2str(n)])
end
end
%
% Plot results
%
% Plot Pivot Growth, gpp and gcp
% Figure 2.1 in textbook
%
disp('Figure 2.1 in textbook')
figure(1), hold off, clg
plot(dim,gpp,'ow',dim,gcp,'+w'),
grid
title('Pivot Growth Factors, GEPP = o, GECP = +')
xlabel('Matrix Dimension')
hold on
% print -deps PivotGrowth1.eps
disp('hit return to continue'), pause
%
% Plot Backward Error Bounds for GEPP
% Figure 2.2 in textbook
%
disp('Figure 2.2 in textbook')
figure(2), hold off, clg
subplot(211)
semilogy(dim,max(bndpp1,1e-18),'xw')
axis([0,max(nn),1e-18,1e-8])
hold on
semilogy(dim,max(bndpp2,1e-18),'+w')
semilogy(dim,max(bndpp3,1e-18),'ow')
grid
plot([0,100],macheps*[1,1],'-w')
% Label plot
title('Backward error in Gaussian elimination with partial pivoting')
xlabel('Matrix dimension')
% text(35,2e-9,'x = 3 n^3 g_pp')
% text(35,1e-11,'+ = 3 n norm(abs(L_pp)*abs(U_pp))/norm(A)')
% text(35,1e-15,'o = norm(A*x_pp-b)/norm(A)*norm(x_pp)')
text(min(nn),1e-9,'macheps = 2^(-53) = 1.1e-16')
%
% Plot Backward Error Bounds for GECP
%
subplot(212)
semilogy(dim,max(bndcp1,1e-18),'xw')
axis([0,max(nn),1e-18,1e-8])
hold on
semilogy(dim,max(bndcp2,1e-18),'+w')
semilogy(dim,max(bndcp3,1e-18),'ow')
plot([0,100],macheps*[1,1],'-w')
grid
% Label plot
title('Backward error in Gaussian elimination with complete pivoting')
xlabel('Matrix dimension')
% text(35,2e-9,'x = 3 n^3 g_cp')
% text(35,1e-11,'+ = 3 n norm(abs(L_cp)*abs(U_cp))/norm(A)')
% text(35,1e-15,'o = norm(A*x_cp-b)/norm(A)*norm(x_cp)')
text(min(nn),1e-9,'macheps = 2^(-53) = 1.1e-16')
% print -deps PivotGrowth2.eps
disp('hit return to continue'), pause
%
% Plot condition number
%
figure(3), hold off, clg
semilogy(dim,cnd1,'+w')
grid
title('Condition number of A')
xlabel('Matrix dimension')
disp('hit return to continue'), pause
%
% Plot ratio of Hager's estimate to true condition number, cnd2./cnd1
%
figure(4), hold off, clg
plot(dim,cnd2./cnd1,'+w')
grid
title('Ratio of Hager''s Estimate to True Condition Number')
xlabel('Matrix dimension')
disp('hit return to continue'), pause
%
% Plot true error bound vs error bound, 
%    bndpp4 vs errpp and bndcp4 vs errcp
% Figure 2.3 or top left of 2.4, in textbook
%
disp('Figure 2.3 or top left of 2.4, in textbook')
lmaxis = floor(log10( ...
         max([bndpp4,bndcp4,bndpp6,bndcp6,bndpp7,bndcp7, ...
              bndpp8,bndcp8,bndpp9,bndcp9,bndpp10,bndcp10])));
maxis = 10^lmaxis;
figure(5), hold off, clg
loglog(errpp,bndpp4,'ow')
% axis([1e-16,1e-10,1e-16,1e-10]);
% lmaxis = floor(log10(max([bndpp4,bndcp4]))); maxis = 10^lmaxis;
axis([1e-16,maxis,1e-16,maxis]);
axis('square');
hold on;
loglog(errcp,bndcp4,'+w')
% loglog([1e-16,1e-10],[1e-16,1e-10],'w-')
loglog([1e-16,maxis],[1e-16,maxis],'w-')
% for i=1:5
for i=1:lmaxis+15,
%  loglog([1e-16,1e-10/10^i],[1e-16*10^i,1e-10],'--w')
   loglog([1e-16,maxis/10^i],[1e-16*10^i,maxis],'--w')
end
title('True Error vs. Error Bound (2.13), o = GEPP, + = GECP')
xlabel('True Error'), ylabel('Error Bound')
% print -deps PivotGrowth3.eps
disp('hit return to continue'), pause
%
% Plot componentwise backward errors, bndpp5, bndcp5
% Figure 2.4, bottom, in textbook
%
disp('Figure 2.4, bottom, in textbook')
figure(6), hold off, clg
semilogy(dim,bndpp5,'ow',dim,bndcp5,'+w'), grid
title('Componentwise relative backward error, o = GEPP, + = GECP')
xlabel('Matrix dimension')
% print -deps PivotGrowth7.eps
disp('hit return to continue'), pause
%
% Plot true error bound vs error bound, Skeel-style, 
%   bndpp6 vs errpp and bndcp6 vs errcp
%
figure(7), hold off, clg
loglog(errpp,bndpp6,'ow')
% axis([1e-16,1e-10,1e-16,1e-10]);
% lmaxis = floor(log10(max([bndpp6,bndcp6]))); maxis = 10^lmaxis;
axis([1e-16,maxis,1e-16,maxis]);
axis('square');
hold on;
loglog(errcp,bndcp6,'+w')
% loglog([1e-16,1e-10],[1e-16,1e-10],'w-')
loglog([1e-16,maxis],[1e-16,maxis],'w-')
% for i=1:5
for i=1:lmaxis+15,
%  loglog([1e-16,1e-10/10^i],[1e-16*10^i,1e-10],'--w')
   loglog([1e-16,maxis/10^i],[1e-16*10^i,maxis],'--w')
end
title('True Error vs. Error Bound (2.7), o = GEPP, + = GECP')
xlabel('True Error'), ylabel('Error Bound')
% print -deps PivotGrowth3a.eps
disp('hit return to continue'), pause
%
% Plot true error bound vs error bound, LAPACK-style, 
%   bndpp7 vs errpp and bndcp7 vs errcp
% Figure 2.4, top right, in textbook
%
disp('Figure 2.4, top right, in textbook')
figure(8), hold off, clg
loglog(errpp,bndpp7,'ow')
% axis([1e-16,1e-10,1e-16,1e-10]);
% lmaxis = floor(log10(max([bndpp7,bndcp7]))); maxis = 10^lmaxis;
axis([1e-16,maxis,1e-16,maxis]);
axis('square');
hold on;
loglog(errcp,bndcp7,'+w')
% loglog([1e-16,1e-10],[1e-16,1e-10],'w-')
loglog([1e-16,maxis],[1e-16,maxis],'w-')
% for i=1:5
for i=1:lmaxis+15,
%  loglog([1e-16,1e-10/10^i],[1e-16*10^i,1e-10],'--w')
   loglog([1e-16,maxis/10^i],[1e-16*10^i,maxis],'--w')
end
title('True Error vs. Error Bound (2.14), o = GEPP, + = GECP')
xlabel('True Error'), ylabel('Error Bound')
% print -deps PivotGrowth4.eps
disp('hit return to continue'), pause
%
% Plot true error bound vs error bound, guaranteed LAPACK-style, 
%   bndpp8 vs errpp and bndcp8 vs errcp
%
figure(9), hold off, clg
loglog(errpp,bndpp8,'ow'), grid
% axis([1e-16,1e-10,1e-16,1e-10]);
% lmaxis = floor(log10(max([bndpp8,bndcp8]))); maxis = 10^lmaxis;
axis([1e-16,maxis,1e-16,maxis]);
axis('square');
hold on;
loglog(errcp,bndcp8,'+w'), grid
% loglog([1e-16,1e-10],[1e-16,1e-10],'w-')
loglog([1e-16,maxis],[1e-16,maxis],'w-')
% for i=1:5
for i=1:lmaxis+15,
%  loglog([1e-16,1e-10/10^i],[1e-16*10^i,1e-10],'--w')
   loglog([1e-16,maxis/10^i],[1e-16*10^i,maxis],'--w')
end
title('True Error vs. Error Bound (2.14), |r| increased, o = GEPP, + = GECP')
xlabel('True Error'), ylabel('Error Bound')
% print -deps PivotGrowth5.eps
disp('hit return to continue'), pause
%
% Plot true error bound vs error bound, LAPACK-style, 
%   after 1 step of iterative refinement
%   bndpp9 vs errppitref and bndcp9 vs errcpitref
%
figure(10), hold off, clg
loglog(errppitref,bndpp9,'ow')
% axis([1e-16,1e-10,1e-16,1e-10]);
% lmaxis = floor(log10(max([bndpp9,bndcp9]))); maxis = 10^lmaxis;
axis([1e-16,maxis,1e-16,maxis]);
axis('square');
hold on;
loglog(errcpitref,bndcp9,'+w')
% loglog([1e-16,1e-10],[1e-16,1e-10],'w-')
loglog([1e-16,maxis],[1e-16,maxis],'w-')
% for i=1:5
for i=1:lmaxis+15,
%  loglog([1e-16,1e-10/10^i],[1e-16*10^i,1e-10],'--w')
   loglog([1e-16,maxis/10^i],[1e-16*10^i,maxis],'--w')
end
title('True Error after Iter Ref vs. Error Bound (2.14), o = GEPP, + = GECP')
xlabel('True Error'), ylabel('Error Bound')
% print -deps PivotGrowth6.eps
disp('hit return to continue'), pause
%
% Plot true error bound vs error bound, guaranteed LAPACK-style, 
%   after 1 step of iterative refinement
%   bndpp10 vs errppitref and bndcp10 vs errcpitref
%
figure(11), hold off, clg
loglog(errppitref,bndpp10,'ow')
% axis([1e-16,1e-10,1e-16,1e-10]);
% lmaxis = floor(log10(max([bndpp10,bndcp10]))); maxis = 10^lmaxis;
axis([1e-16,maxis,1e-16,maxis]);
axis('square');
hold on;
loglog(errcpitref,bndcp10,'+w')
% loglog([1e-16,1e-10],[1e-16,1e-10],'w-')
loglog([1e-16,maxis],[1e-16,maxis],'w-')
% for i=1:5
for i=1:lmaxis+15,
%  loglog([1e-16,1e-10/10^i],[1e-16*10^i,1e-10],'--w')
   loglog([1e-16,maxis/10^i],[1e-16*10^i,maxis],'--w')
end
title('True Error after Iter Ref vs. Error Bound (2.14), |r| increased, o = GEPP, + = GECP')
xlabel('True Error'), ylabel('Error Bound')
disp('hit return to continue'), pause
%
% Plot componentwise backward errors after iterative refinement, 
%     bndpp11, bndcp11
%
figure(12), hold off, clg
semilogy(dim,bndpp11,'ow',dim,bndcp11,'+w'), grid
title('Componentwise backward error after Iter Ref, o = GEPP, + = GECP')
xlabel('Matrix dimension')