# Apply smoothing operator to a, putting result in b.
function relax{T}( b::Vector{T}, a::Vector{T} )
    assert(length(a)==length(b))
    c = T(0.25)
    d = T(0.5)
    n = length(b)
    b[1] = 1                            # Boundary condition
    @inbounds for i=2:n-1
        b[i] = c*a[i-1]+d*a[i]+c*a[i+1]
    end
    b[n] = 0                            # Boundary condition
end

# Apply smoothing operator 1000 times
function flog{T}( b::Vector{T}, a::Vector{T} )
    for t=1:500
        relax(b,a)
        relax(a,b)
    end
end

# Warmup
flog(rand(Float32,100), rand(Float32,100))

# Do some timing trials.
for i=1:12
    b = zeros(Float32,1000)
    if i%3==0
        # Use tiny noise instead of zeros
        a = rand(Float32,1000) .* 1E-9
        input = "nano-noise"
    else
        a = zeros(Float32,1000)
        input = "zeros"
    end
    if i%3==2
        set_zero_subnormals(true)
        rounding = "subnormals are zero"
    else
        set_zero_subnormals(false)
        rounding = "IEEE"
    end
    a[1] = 1
    print(rpad(string(input," with ",rounding),30))
    @time flog(b,a)
end