%
% QR ·Ö½âÎÈ¶¨ÐÔ²âÊÔ
%
% random matrix with given condition number
%

clearvars;
close all;

M = 256; N = 64;
Cond = N.^[1:7];  % condition number
Kind = 2; % 1 - matrix with cluster singular values and 
          % 2 - matrix with exponential linear singular values 
if Kind == 1
    fig_tit_cond = ': cluster singular values';
elseif Kind == 2 
    fig_tit_cond = ': exponential linear singular values';
else
    fprintf('Error: Kind=%d\n',Kind); return;
end

for k = 1 : length(Cond)
    fprintf('***** M=%d, N=%d, cond=%d^%d *****\n', M,N,N,k)
    H = setmatrix(M,N,Cond(k),Kind);
    fprintf('***** M=%d, N=%d, cond=%.2e *****\n', M,N,cond(H))

    hnorm = norm(H,'fro');
    Inorm = norm(eye(N),'fro');

    [Q_GS,R_GS] = qr_GS(H);
    Err_GS(1,k) = norm(Q_GS*R_GS - H,'fro')/hnorm;
    Err_GS(2,k) = norm(Q_GS'*Q_GS - eye(N),'fro')/Inorm;
    fprintf('qr_GS:  ||Q*R-H||/||H||=%.2e, ||Q''*Q-I||/||I||=%.2e\n', Err_GS(1,k),Err_GS(2,k))

    [Q_MGS,R_MGS] = qr_MGS2(H);
    Err_MGS21(1,k) = norm(Q_MGS*R_MGS - H,'fro')/hnorm;
    Err_MGS21(2,k) = norm(Q_MGS'*Q_MGS - eye(N),'fro')/Inorm;
    fprintf('qr_MGS: ||Q*R-H||/||H||=%.2e, ||Q''*Q-I||/||I||=%.2e\n', ...
        Err_MGS21(1,k), Err_MGS21(2,k))

    [Q_MGS2,R_MGS2] = qr_MGS2(Q_MGS);
    Err_MGS2(1,k) = norm(Q_MGS2*R_MGS2*R_MGS - H,'fro')/hnorm;
    Err_MGS2(2,k) = norm(Q_MGS2'*Q_MGS2 - eye(N),'fro')/Inorm;
    fprintf('qr_MGS2:||Q*R-H||/||H||=%.2e, ||Q''*Q-I||/||I||=%.2e\n', ...
        Err_MGS2(1,k),Err_MGS2(2,k))

    [Q_H,R_H] = qr_HH(H);
    Err_HH(1,k) = norm(Q_H*R_H - H,'fro')/hnorm;
    Err_HH(2,k) = norm(Q_H'*Q_H - eye(N),'fro')/Inorm;
    fprintf('qr_HH:  ||Q*R-H||/||H||=%.2e, ||Q''*Q-I||/||I||=%.2e\n', Err_HH(1,k),Err_HH(2,k))

    [Q_Chol,R_Chol] = qr_chol(H);
    Err_chol(1,k) = norm(Q_Chol*R_Chol - H,'fro')/hnorm;
    Err_chol(2,k) = norm(Q_Chol'*Q_Chol - eye(N),'fro')/Inorm;
    fprintf('qr_chol:||Q*R-H||/||H||=%.2e, ||Q''*Q-I||/||I||=%.2e\n', Err_chol(1,k),Err_chol(2,k))

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

FF=figure;
% orient(FF,'landscape');
% FF.PaperType='A3';
xx = 1:length(Cond);
FF1=semilogy(xx,Err_GS(2,:),'o-', xx,Err_MGS21(2,:),'+-', ...
    xx,Err_MGS2(2,:),'d-', xx,Err_HH(2,:),'s-', ...
    xx,Err_chol(2,:),'v-', 'MarkerSize',8,'LineWidth',1);
legend('GS', 'MGS', 'MGS2', 'Householder','Chol','Location','northwest')
% FF1=semilogy(xx,Err_GS(2,:),'o-', xx,Err_MGS21(2,:),'^-', ...
%     xx,Err_MGS2(2,:),'d-', xx,Err_HH(2,:),'s-', ...
%     'MarkerSize',12,'LineWidth',2);
legend('GS', 'MGS', 'MGS(double)', 'Householder','Cholesky','Location','northwest')
xlabel('Condition number: $64^x$', 'Interpreter','latex');
ylabel('$\frac{\|Q^T Q-I\|_F}{\|I\|_F}$', 'Interpreter','latex');
xticks(xx);
xticklabels(xx);
set(gca,'Fontsize',21)
%set(get(gca,'xlabel'),'fontsize',20)
%set(get(gca,'ylabel'),'fontsize',20)
%set(gcf,'Position',[100,100,600,600])          
fig_tit = ['Random Matrix $256 \times 64$',fig_tit_cond];
title(fig_tit, 'Interpreter','latex')
set(gcf, 'Units', 'normalized', 'Position', [0.1 0.25 0.57 0.65]);
set(gca, 'Units', 'normalized', 'Position', [0.12 0.12 0.82 0.82]);
ylim([1e-16,1e1]);
FF.PaperPositionMode = 'auto';
fig_pos = FF.PaperPosition;
fig.PaperSize = [fig_pos(3) fig_pos(4)];
figname = ['QR_3methods_256_64'];
% print('-fillpage',figname,'-dpdf')
print('-dpng','-r200',figname)



%%%%%%%%%%%%%%%%%%%%%%%%%%
function A = setmatrix(M,N,Kappa,Kind)
% set random matrix with given condition number
%
if Kind == 1 % cluster singular values
    Sigma = (Kappa-1)*rand(N-2,1) + 1;  % singular values
    Sigma = [Kappa; Sigma; 1]/Kappa;    % set max and min singular values
elseif Kind == 2 % exponential linear singular values
    Sigma = 10.^(linspace(0,log10(Kappa),N))/Kappa;
end

[U,S,V] = svd(randn(M,N),'econ');
A = U*diag(Sigma)*V';
end


function [Q,R] = qr_GS(A)
[m,n] = size(A);
Q = zeros(m,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

% column-oriented MGS
function [Q,R] = qr_MGS(A)
[m,n] = size(A);
Q = zeros(m,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

% row-oriented MGS
function [Q,R] = qr_MGS2(A)
[m,n] = size(A);
Q = A; R = zeros(n);
R(1,1) = norm(Q(:,1));
Q(:,1) = Q(:,1)/R(1,1);
for k = 2 : n
    R(k-1,k:n) = Q(:,k-1)'*Q(:,k:n);
    for j = k : n
        Q(:,j) = Q(:,j) - R(k-1,j)*Q(:,k-1);
    end
    R(k,k) = norm(Q(:,k));
    Q(:,k) = Q(:,k)/R(k,k);
end
end

function [Q,R] = qr_chol(A)
n = size(A,2);
B = A'*A;
R = chol(B+10*eps*max(B(:))*eye(n));  % ±ÜÃâ B ÆæÒì
Q = A/R;
end

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