using LinearAlgebra

U = qr(randn(80,80)).Q # U = random 80×80 orthogonal matrix
V = qr(randn(80,80)).Q # V = random 80×80 orthogonal matrix
Σ = Diagonal(2.0 .^ (-1:-1:-80)) # Σ = diagonal matrix with entries 2^-1, 2^-2, ..., 2^-80
A = U * Σ * V'

# Classical Gram–Schmidt (Trefethen algorithm 7.1), implemented in the simplest way
# (We could make it faster by unrolling loops to avoid temporaries arrays etc.)
function clgs(A)
    m,n = size(A)
    T = float(eltype(A))
    Q = similar(A, T)
    R = zeros(T,n,n)
    @views for j = 1:n
        aⱼ = A[:,j]
        vⱼ = copy(aⱼ) # use copy so that modifying vⱼ doesn't change aⱼ
        for i = 1:j-1
            qᵢ = Q[:,i]
            R[i,j] = qᵢ'aⱼ
            vⱼ .-= R[i,j] .* qᵢ
        end
        R[j,j] = norm(vⱼ)
        Q[:,j] .= vⱼ ./ R[j,j]
    end
    return Q, R
end

# Modified Gram–Schmidt (Trefethen algorithm 8.1)
function mgs(A)
    m,n = size(A)
    T = float(eltype(A))
    Q = similar(A, T)
    R = zeros(T,n,n)
    @views for j = 1:n
        aⱼ = A[:,j]
        vⱼ = copy(aⱼ)
        for i = 1:j-1
            qᵢ = Q[:,i]
            R[i,j] = qᵢ'vⱼ # ⟵ NOTICE: mgs has vⱼ, clgs has aⱼ
            vⱼ .-= R[i,j] .* qᵢ
        end
        R[j,j] = norm(vⱼ)
        Q[:,j] .= vⱼ ./ R[j,j]
    end
    return Q, R
end

QC, RC = clgs(A);
QM, RM = mgs(A);

using PyPlot

n = size(A,2)
semilogy(diag(RC), "bo")
semilogy(diag(RM), "rx")
semilogy(2.0 .^ -(1:n), "k-")
semilogy(fill(sqrt(eps()), n), "b--")
semilogy(fill(eps(), n), "r--")
legend(["classical","modified",L"2^{-j}", L"\sqrt{\epsilon_\mathrm{mach}}", L"\epsilon_\mathrm{mach}"], loc="lower left")
ylabel(L"r_{jj}")
xlabel(L"j")

QCbig, RCbig = clgs(big.(A))

semilogy(Float64.(diag(RCbig)), "o", mfc="none") # convert back to Float64 for plotting, since PyPlot doesn't understand BigFloat
semilogy(diag(RM), "rx")
legend(["classical GS in arbitrary precision","modified GS"], loc="upper right")
ylabel(L"r_{jj}")
xlabel(L"j")

opnorm(QC'*QC - I, 1), opnorm(QM'*QM - I, 1), opnorm(QCbig'*QCbig - I, 1)

QRH = qr(A)
QH, RH = QRH.Q, QRH.R
semilogy(abs.(diag(RH)), "s", mfc="none") # need abs because householder has arbitrary phase factor in R
semilogy(diag(RM), "rx")
legend(["Householder QR","modified GS"], loc="upper right")
ylabel(L"r_{jj}")
xlabel(L"j")

println("L1 NORM of Householder Q'Q - I: ", opnorm(QH'*QH - I, 1))