%  Test Gaussian elimination, computing the inverse of A, and Cramer's
%  rule for solving linear systems to see how they work on ill-conditioned
%  problems.  See which algorithms appear to be backward stable.  When
%  the system is ill-conditioned, we cannot expect any of the methods
%  to give an answer that is close to the true solution, but if an algorithm
%  is backward stable, the residual norm ||b-Ax||/(||A|| ||x||)
%  should be on the order of machine precision even for the most ill-conditioned
%  linear systems.

%  Set problem size.
n = 10;

%  Print headings for table.
fprintf(' condno      Gaussian elim          A inverse            Cramers rule \n')
fprintf('           error   residual      error   residual      error   residual\n')
fprintf('-----------------------------------------------------------------------\n')

%  Loop over condition numbers.

for k=0:4:16,
  condno = 10^k;
  A = matgen(n,condno);   % Form matrix with condition number condno.

  xtrue = randn(n,1);     % Set true solution.
  b = A*xtrue;            % Form right-hand side vector.

  x_ge = A\b;             % Solve with Gaussian elimination.
  x_inv = inv(A)*b;       % Solve by computing inv(A).
  x_cramer = zeros(n,1);  % Solve by Cramer's rule. 
  detA = det(A);
  for j=1:n, 
    Aj = A; Aj(:,j) = b; 
    detAj = det(Aj); 
    x_cramer(j) = detAj/detA; 
  end;

%    Compute errors and residuals.
  err_ge = norm(x_ge - xtrue)/norm(xtrue);
  resid_ge = norm(b - A*x_ge)/(norm(A)*norm(x_ge));
  err_inv = norm(x_inv - xtrue)/norm(xtrue);
  resid_inv = norm(b-A*x_inv)/(norm(A)*norm(x_inv));
  err_cramer = norm(x_cramer - xtrue)/norm(xtrue);
  resid_cramer = norm(b-A*x_cramer)/(norm(A)*norm(x_cramer));

%    Print results.
  fprintf('%6.0e    %6.2e  %6.2e    %6.2e  %6.2e    %6.2e  %6.2e\n',...
          condno,err_ge,resid_ge,err_inv,resid_inv,err_cramer,resid_cramer)

end;