% This script generates a a matrix A of order 100 with condition number
% 10^sig.  It also generates the right hand side of a linear system Ax=b
% and an approximate solution xt that is accurate to about 16-sig
% decimal digits.  The linear system is then solved with an underlying
% precision of rho decimal digits and performs four steps of iterative
% refiniment with the residual evaluated in tau-digit decimal floating
% point arithmetic.

global FlapRlevel

% Set the basic rounding level.

SetFlapRlevel(rho)

% Generate a system with log(cond(A)) = sig.

n = 100;
s = logspace(0,-sig, n);
[U, trash] = qr(randn(n));
[V, trash] = qr(randn(n));
A = U*diag(s)*V';
A = flap(A);
b = A*ones(n,1);
xt = A.d\b.d;

% Solve the system with iterative refinement.
% The residuals are evaluated to tau digits.

x = A\b;
fprintf('tau = %2d\n', fix(tau));
fprintf('%8.1e\n', norm(x.d - xt)/norm(xt));
for i=1:4
   SetFlapRlevel(tau);
      r = b - A*x;
   SetFlapRlevel(rho);
   x = x + A\r;
fprintf('%8.1e\n', norm(x.d - xt)/norm(xt));
end