% C1

clc
nExamples = 1000;
disp('sn = sigma_min(A)')
disp('snU = sigma_min(U) where PA = LU')
disp('snEst = sigma_min estimator applied to U')
disp(' ')
disp('nExample = 1000 = Number of examples per choice of n and sigma_min(A)...')
disp(' ')
disp(' ')
disp(' sn n sn/snU snEst/snU snEst/Snu')
disp(' (Average) (Average) (Maximum)')
disp('--------------------------------------------------------------------')

for d = [2 4 8 12]
sn = 10^(-d);
for n = [10 20 50 100 200]
R1 = 0;
R2 = 0;
M1 = 0;
for k=1:nExamples
A = randn(n,n);
[U,S,V] = svd(A);
S(n,n) = sn;
A = U*S*V';
[L,U,P] = lu(A);
snEst = SigmaMinEst(U);
snU = min(svd(U));
R1 = R1 + sn/snU;
R2 = R2 + snEst/snU;
M1 = max(M1,snEst/snU);
end
R1 = R1/nExamples;
R2 = R2/nExamples;
fprintf('10-%1d %3d %10.3f %10.6f %10.6f\n',d,n,R1,R2,M1)
end
disp(' ')
end