%
% QR 分解稳定性演示
%

clear all;

%Name_Problem='Random Matrix';
Name_Problem='Hilbert';

if strcmp(Name_Problem,'Hilbert')
    NN=2:2:20;
else
    NN = 50:50:500;
end

for k = 1:length(NN)
    n = NN(k);    
    if strcmp(Name_Problem,'Hilbert')
        H = hilb(n);
    else
        H = randn(n);
    end

    fprintf('***** n=%d *****\n',n)
    
    [Q1,R1] = qr_GS(H);
    
    Err_GS(1,k) = norm(Q1*R1 - H);
    Err_GS(2,k) = norm(Q1'*Q1 - eye(n));
    fprintf('qr_GS:  ||Q*R-H||=%.2e, ||Q''*Q-I||=%.2e\n', Err_GS(1,k),Err_GS(2,k))
    
    [Q2,R2] = qr_MGS(H);
    Err_MGS(1,k) = norm(Q2*R2 - H);
    Err_MGS(2,k) = norm(Q2'*Q2 - eye(n));
    fprintf('qr_MGS: ||Q*R-H||=%.2e, ||Q''*Q-I||=%.2e\n', Err_MGS(1,k),Err_MGS(2,k)) 

    [Q3,R3] = qr_HH(H);
    Err_HH(1,k) = norm(Q3*R3 - H);
    Err_HH(2,k) = norm(Q3'*Q3 - eye(n));
    fprintf('qr_HH:  ||Q*R-H||=%.2e, ||Q''*Q-I||=%.2e\n', Err_HH(1,k),Err_HH(2,k))

    [Q,R]   = qr(H);  % build-in qr
    Err(1) = norm(Q*R - H);
    Err(2) = norm(Q'*Q - eye(n));
    fprintf('qr:     ||Q*R-H||=%.2e, ||Q''*Q-I||=%.2e\n\n', Err(1),Err(2))
end

figure
semilogy(NN,Err_GS(2,:),'o-', NN,Err_MGS(2,:),'+-', NN,Err_HH(2,:),'s-')
legend('GS', 'MGS', 'Householder')


function [Q,R] = qr_GS(A)
n = size(A,1);
Q = zeros(n); R = zeros(n);
R(1,1) = norm(A(:,1));
Q(:,1) = A(:,1)/R(1,1);
for k = 2 : n
    q = A(:,k);
    qtmp = q;
    for i = 1 : k-1
        R(i,k) = qtmp'*Q(:,i);
        q = q - R(i,k)*Q(:,i);
    end
    R(k,k) = norm(q);
    Q(:,k) = q/R(k,k);
end
end

function [Q,R] = qr_MGS(A)
n = size(A,1);
Q = zeros(n); R = zeros(n);
R(1,1) = norm(A(:,1));
Q(:,1) = A(:,1)/R(1,1);
for k = 2 : n
    q = A(:,k); 
    for i = 1 : k-1
        R(i,k) = q'*Q(:,i);
        q = q - R(i,k)*Q(:,i);
    end
    R(k,k) = norm(q);
    Q(:,k) = q/R(k,k);
end
end

function [Q,R] = qr_HH(A)
n = size(A,1);
Q = eye(n); 
for k = 1 : n
    x = A(k:n,k); 
    [beta,v] = House(x);
    A(k:n,k:n) = A(k:n,k:n) - (beta*v)*(v'*A(k:n,k:n));
    Q(:,k:n) = Q(:,k:n) - Q(:,k:n)*(beta*v)*v';
end
R = A;
end
