%
%  risoluzione di un sistema lineare fortemente mal-condizionato
%
stde = 0*1e-6;  % deviazione standard del rumore additivo aggiunto
n=64;
[U, R] = qr(rand(n)); 
[V, R] = qr(rand(n)); 
if 0
  S = diag(exp(1e1*[0:-1:-n+1]));
else
  S = diag(exp(0.6*[0:-1:-n+1]));
end
figure(1), semilogy(log10(diag(S)));
A = U*S*V';
K2_A = cond(A);
disp(['numero di condizionamento di A = ' num2str(K2_A,'%g')]);
campa=1.1*hamming(26)'; 
x_vero = [zeros(5,1); 0.75*ones(15,1); zeros(5,1); 0.2*ones(2,1); zeros(7,1); campa-campa(1); zeros(4,1)];
b = A*x_vero + stde*randn(n,1);
x_LU = A \ b;
n2resLU = norm(A*x_LU - b);
figure(3); clf, plot(x_vero,'b-'); hold on, pause, plot(x_LU,'m'); pause,
%
n2res_hist = zeros(n,1);
limsup_beta = 2 / S(1)^2;
beta = limsup_beta*0.48;
ATA = A'*A;
x_LW = zeros(n,1);
n_iter = 200;
for k=1:n_iter
  disp(['k=' num2str(k)]);
  x_LW = x_LW + beta * (A'*b - ATA*x_LW);
  n2xTSVD = norm(x_LW);
  res = A*x_LW - b;
  n2res = norm(res);
  n2res_hist(k) = n2res;
  %
  %figure(3); plot(x_TSVD,'m'); %pause,
  figure(4); clf, plot(x_vero,'b-'); hold on, plot(x_LU,'m-'); plot(x_LW,'r'); axis([1 n -.5 1.5]); %pause,
end

figure(10); clf, hold on,
            semilogy([1 n_iter],log10([stde stde]),'r--');
            semilogy(log10(n2res_hist),'b-');
            axis([1 n_iter log10(1e-18) log10(1e1)]);

            