reorthogonalization.jl

reortho.jl

import Base.LinAlg.BlasFloat
import Base.LinAlg.BLAS: gemv!, gemv, axpy!

using BenchmarkTools

function classical_gram_schmidt!{T<:BlasFloat}(V::StridedMatrix{T}, w::StridedVector{T})
    
    # orthogonalize
    h = gemv('T', one(T), V, w)
    gemv!('N', -one(T), V, h, one(T), w)

    # Normalize
    nrm = norm(w)
    scale!(w, one(T) / nrm)

    h, nrm
end

function iterative_classical_gram_schmidt!{T<:BlasFloat}(V::StridedMatrix{T}, w::StridedVector{T}; η = one(T) / √2)
    
    # orthogonalize
    h = gemv('T', one(T), V, w)
    gemv!('N', -one(T), V, h, one(T), w)

    # Repeat
    nrm = norm(w)

    if nrm < η * norm(h)
        increment = gemv('T', one(T), V, w)
        gemv!('N', -one(T), V, increment, one(T), w)
        axpy!(one(T), increment, h)
        nrm = norm(w)
    end

    # Normalize
    scale!(w, one(T) / nrm)

    h, nrm
end

function modified_gram_schmidt!{T<:BlasFloat}(V::StridedMatrix{T}, w::StridedVector{T})
    k = size(V, 2)
    h = zeros(T, k)

    # Orthogonalize
    @inbounds for i = 1 : k
        column = @view(V[:, i])
        h[i] = dot(column, w)
        axpy!(-h[i], column, w)
    end

    # Normalize
    nrm = norm(w)
    scale!(w, one(T) / nrm)

    h, nrm
end

function modified_gram_schmidt!{numT <: BlasFloat, vecT<:StridedVector}(V::Vector{vecT}, w::StridedVector{numT})
    k = length(V)
    h = zeros(numT, k)

    # Orthogonalize
    @inbounds for i = 1 : k
        h[i] = dot(V[i], w)
        axpy!(-h[i], V[i], w)
    end

    # Normalize
    nrm = norm(w)
    scale!(w, one(numT) / nrm)

    h, nrm
end

function test(; m = 50)
    srand(1)
    suite = BenchmarkGroup()
    
    for n = [1_000, 10_000, 100_000, 1_000_000]
        key = string(n)
        suite[key] = BenchmarkGroup(["cgs", "icgs", "mgs"])
        
        V, _ = qr(rand(n, m))
        new_vec = rand(n)
        suite[key]["cgs"] = @benchmark classical_gram_schmidt!($V, a) setup=(a = copy($new_vec))
        suite[key]["icgs"] = @benchmark iterative_classical_gram_schmidt!($V, b) setup=(b = copy($new_vec))
        suite[key]["mgs"] = @benchmark modified_gram_schmidt!($V, c) setup=(c = copy($new_vec))
    end

    suite
end

results

n = 1000

"icgs" => Trial(41.635 μs)
"mgs" => Trial(182.874 μs)
"cgs" => Trial(18.854 μs)


n = 10.000

"icgs" => Trial(616.911 μs)
"mgs" => Trial(789.294 μs)
"cgs" => Trial(319.054 μs)


n = 100.000

"icgs" => Trial(10.128 ms)
"mgs" => Trial(7.047 ms)
"cgs" => Trial(5.054 ms)


n = 1.000.000

"icgs" => Trial(106.763 ms)
"mgs" => Trial(117.694 ms)
"cgs" => Trial(53.784 ms)