% This example illustrates the fact that the standard backward
% rounding-error bound for the Cholesky decomposition can be much larger
% than the backward error itself.  It generates a symmetric positive
% definite matrix of order n with norm 1 and condition number 10^t.  It
% then computes the Cholesky factor R in single precision and the
% backward error A - R'*R in double precision.  The classical bound is
% printed out for comparison (see Higham, Accuracy and Stability of
% Numerical Algorithms, 2nd edition, p. 198).

% Generate A.

n = 100;
t = 1;
D = diag(logspace(0,-t,n));
[Q,trash] = qr(randn(n));
A = Q'*(D*Q);

% Change A to a single precision flap and compute
% its Cholesky decomposition.

SetFlapRlevel(15);
A = flap(A);
R = chold(A);

% Print out the norm of the backward error and the bound
% from the rounding-error analysis.

BackwardError = norm(A.d - R.d'*R.d)
BackwardErrorBound = 4*n*(3*n + 1)*5.96e-8