# L-19 MCS 471 Wed 5 Oct 2022 : gsqr.jl

# This script provides a simple implementation of the
# Gram-Schmidt method for the QR factorization.

using Printf
Base.show(io::IO, f::Float64) = @printf(io, "%.3e", f)

using LinearAlgebra

"""
    gsqr(A::Array{Float64,2})

Gram-Schmidt orthogonalization for a QR decomposition:

    Q, R = gsqr(A) returns a QR decomposition of A,
computed with 64-bit floating-point arithmetic.

ON ENTRY :

    A         an n-by-k Float64 matrix, n >= k.

ON RETURN :

    Q         an orthogonal matrix: Q'Q = I, whose column span 
              is the same as the column span of A;

    R         an upper triangular matrix: Q*R = A.

EXAMPLE :

    a = rand(4,3);    # a random 4-by-3 matrix
    q, r = gsqr(a);
    a - q*r           # check residual
    transpose(q)*q    # check orthogonality
"""
function gsqr(A::Array{Float64,2})
    k = size(A,2)
    Q = deepcopy(A)
    R = zeros(k,k)
    for j = 1:k
        p = A[:,j]
        for i = 1:j-1
            R[i,j] = transpose(Q[:,i])*A[:,j]
            p = p - Q[:,i]*R[i,j]
        end
        Q[:,j] = p/norm(p,2);
        R[j,j] = transpose(Q[:,j])*A[:,j]
    end
    return Q, R
end

"""
    gsqr(A::Array{BigFloat,2})

Gram-Schmidt orthogonalization for a QR decomposition:

    Q, R = gsqr(A) returns a QR decomposition of A,
computed with multiprecision floating-point arithmetic.

ON ENTRY :

    A         an n-by-k BigFloat matrix, n >= k.

ON RETURN :

    Q         an orthogonal matrix: Q'Q = I, whose column span 
              is the same as the column span of A;
    R         an upper triangular matrix: Q*R = A.

EXAMPLE :

    a = rand(4,3);            # a random 4-by-3 matrix
    A = Array{BigFloat,2}(a)  # convert to BigFloat
    q, r = gsqr(a);
    a - q*r                   # check residual
    transpose(q)*q            # check orthogonality
"""
function gsqr(A::Array{BigFloat,2})
    k = size(A,2)
    Q = deepcopy(A)
    R = Array{BigFloat,2}(zeros(k,k))
    for j = 1:k
        p = A[:,j]
        for i = 1:j-1
            R[i,j] = transpose(Q[:,i])*A[:,j]
            p = p - Q[:,i]*R[i,j]
        end
        Q[:,j] = p/norm(p,2);
        R[j,j] = transpose(Q[:,j])*A[:,j]
    end
    return Q, R
end

"""
Generates a random 4-by-3 matrix
and tests the gsqr function
with 64-bit floating-point numbers.
"""
function testFloat64()
    println("Testing the 64-bit version ...")
    a = rand(4, 3)
    q, r = gsqr(a)
    res = norm(a - q*r)
    stres = @sprintf("%.2e", res)
    println("The residual : $stres")
    println("Testing orthogonality, Q^T*Q is")
    show(stdout, "text/plain", transpose(q)*q); println("");
end

"""
Generates a random 4-by-3 matrix
and tests the gsqr function
with 64-bit floating-point numbers.
"""
function testBigFloat()
    println("Testing the BigFloat version ...")
    a = rand(4, 3)
    A = Array{BigFloat,2}(a)
    q, r = gsqr(A)
    res = norm(A - q*r)
    stres = @sprintf("%.2e", res)
    println("The residual : $stres")
    println("Testing orthogonality, Q^T*Q is")
    show(stdout, "text/plain", transpose(q)*q); println("");
end

"""
Runs the tests on gsqr.
"""
function main()
    testFloat64()
    testBigFloat()
end

main()