Blackscholes ported from Python to Julia for use with ParallelAccelerator

blackscholes.jl

using ParallelAccelerator

const RISKFREE = 0.02
const VOLATILITY = 0.30

#ParallelAccelerator.ParallelIR.PIRSetFuseLimit(0)
#ParallelAccelerator.CGen.setvectorizationlevel(ParallelAccelerator.CGen.VECDISABLE)

@acc begin

@inline function cnd(d)
    A1 = 0.31938153
    A2 = -0.356563782
    A3 = 1.781477937
    A4 = -1.821255978
    A5 = 1.330274429
    RSQRT2PI = 0.39894228040143267793994605993438
    K = 1.0 ./ (1.0 .+ 0.2316419 .* abs(d))
    ret_val = (RSQRT2PI .* exp(-0.5 .* d .* d) .*
               (K .* (A1 .+ K .* (A2 .+ K .* (A3 .+ K .* (A4 .+ K .* A5))))))
    m = d .> 0
    ret_val[m] = 1.0 .- ret_val[m]
    return ret_val
end

function black_scholes(stockPrice, optionStrike, optionYears,
                  Riskfree, Volatility)
    S = stockPrice
    X = optionStrike
    T = optionYears
    R = Riskfree
    V = Volatility
    sqrtT = sqrt(T)
    d1 = (log(S ./ X) .+ (R .+ 0.5 .* V .* V) .* T) ./ (V .* sqrtT)
    d2 = d1 .- V .* sqrtT
    cndd1 = cnd(d1)
    cndd2 = cnd(d2)

    expRT = exp(- R .* T)
    return (S .* cndd1 .- X .* expRT .* cndd2), (X .* expRT .* (1.0 - cndd2) .- S .* (1.0 .- cndd1))
end

end

function randfloat(N, lo, hi)
    rand(N) .* (hi - lo) .+ lo
end

function main()
    OPT_N = 4000000
    iterations = 10
    if length(ARGS) >= 2
        iterations = int(ARGS[1])
    end
    srand(0)
    callResult = zeros(OPT_N)
    putResult = -ones(OPT_N)
    stockPrice = randfloat(OPT_N, 5.0, 30.0)
    optionStrike = randfloat(OPT_N, 1.0, 100.0)
    optionYears = randfloat(OPT_N, 0.25, 10.0)

    black_scholes(Float64[], Float64[],
                      Float64[], Float64[], Float64[])
    tic()
    for i in 1:iterations
        callResult, putResult = black_scholes(stockPrice, optionStrike,
                      optionYears, RISKFREE, VOLATILITY)
    end
    t = toq()
    println("Checksum: ", sum(callResult), " ", sum(putResult))
    println("Time: ", (t/ iterations * 1000), " msec")
end

main()