Julia speed

General Usage Performance
1

Hello !
I teach computational Physics in a french university. I learn to my students to code in Modern Fortran and Python. I know that many people claim that Julia may be as fast as C or Fortran. So, I try a little test : to compute pi with the Madhava-Leibniz formula. I wrote the same code in Fortran, Python and Julia.
Here are my codes :
Fortran :
program calcPiMadhava
use, intrinsic :: iso_fortran_env, only : wp => real64
implicit none
real(wp), parameter :: pi = 4.0_wp*atan(1.0_wp)
integer :: k = 0, nb
real(wp) :: S1 = 0.0_wp, S2 = 0.0_wp, eps
logical :: continuer = .True.
print *, "Number of significant digits expected ? "
read *, nb

eps = (0.1)**(nb+1)

do while (continuer)
S1 = S1 + 4.0_wp*((-1.0_wp)**k)/(2.0_wp*k + 1.0_wp)
if (abs(S1-S2) < eps) then
continuer = .False.
else
S2 = S1
k = k+1
end if
end do

print ‘("value of pi calculated with Madhava formula : ", f20.18)’, S1
print *, "pi = ", pi
print ‘("error = ", E40.30)’, abs(S1-pi)
print ‘("number of iterations : ", i10.1)’, k

end program calcPiMadhava

Python :

from numpy import pi
S1 = 0.0
S2 = 0.0
k = 0
continuer = True
nb = float(input(“Number of significant digits expected ?”))
eps = 0.1**(nb+1)

while continuer:
S1 = S1 + 4.0*((-1.0)**k)/(2.0*k + 1.0)
if abs(S1-S2) < eps:
continuer = False
else:
S2 = S1
k = k+1

print(f"value of pi calculated with Madhava formula : {S2}“)
print(f"Number of iterations : {k}”)
print(f"π = {pi}“)
print(f"error = {abs(S2-pi)}”)

Julia :

S1 = 0.0
S2 = 0.0
k = 0
continuer = true
print(“Number of significant digits expected ?”)
nb = parse(Int,readline())
eps = 0.1^(nb+1)

while continuer
global continuer, S1, S2, k, eps
S1 = S1 + 4.0*((-1.0)^k)/(2.0*k + 1.0)
if abs(S1-S2) < eps
continuer = false
else
S2 = S1
k = k+1
end
end

println("Value of pi calculated with Madhava formula : ", S2)
println("Number of iterations : ", k)
println("π = ", BigFloat(pi))
println("error = ", abs(S2-π))

If execute this 3 codes under Bash with the instructions time and ask 7 significant digits, it takes 17 s to Fortran to perform the calculation, 2.46 minutes to Python and 3.23 minutes to Julia.
So, Julia is the slowest …
I’m new to Julia, and my code is certainly perfectible, but I’d like to know how to obtain fast code in Julia.
At this time, I learn to my students to prototype codes under Python and translate them in Fortran if it is too slow. If Julia was as fast as Fortran, I would learn them Julia … And this is the reason of my question… Is Julia really fast for finite differences, FDTD, solving time dependent Schrodinger equation …
Many thanks in advance for your answers and help !
PE

2

Welcome to the community!

Please read: Performance Tips · The Julia Language

Key point: Put your code into a function.

Additional remark: If you share code, please enclose it with triple backticks, like ``` and indent it
for readability.

12 Likes
3

I just did what @ufechner7 suggested and used @btime to take out the compilation time…

using BenchmarkTools
function main(nb)
    S1 = 0.0
    S2 = 0.0
    k = 0
    continuer = true
    eps = 0.1^(nb+1)

    while continuer
#      global continuer, S1, S2, k, eps
        S1 = S1 + 4.0*((-1.0)^k)/(2.0*k + 1.0)
        if abs(S1-S2) < eps
            continuer = false
        else
            S2 = S1
            k = k+1
        end
    end

    println("Value of pi calculated with Madhava formula : ", S2)
    println("Number of iterations : ", k)
    println("π = ", BigFloat(pi))
    println("error = ", abs(S2-π))
end

print("Number of significant digits expected ?")
nb = parse(Int, readline())

@btime main($nb)

I responded 7 to the initial question. Just reporting the last output:

Value of pi calculated with Madhava formula : 3.1415926585894076
Number of iterations : 199999997
π = 3.141592653589793238462643383279502884197169399375105820974944592307816406286198
error = 4.999614500178495e-9
  12.927 s (53 allocations: 6.31 KiB)
4

Many thanks for your answers ! I will look and try what you told me.
Sorry for such basics questions but it’s been a long time that I’ve been trying to learn Julia and see if I should at the end teach it to my students.
But I have still a lot of things to lear ! :slight_smile:
Many thanks again for your help !

1 Like
5

Same as above, but measuring everything:

script pi.jl 
function compute_pi(S1, S2, k, eps)
    continuer = true
    while continuer
        S1 = S1 + 4.0*((-1.0)^k)/(2.0*k + 1.0)
        if abs(S1-S2) < eps
            continuer = false
        else
            S2 = S1
            k = k+1
        end
    end
    return S2, k
end

S1 = 0.0
S2 = 0.0
k = 0
print("Number of significant digits expected ?")
nb = parse(Int,readline())
eps = 0.1^(nb+1)

S2, k = compute_pi(S1, S2, k, eps)

println("Value of pi calculated with Madhava formula : ", S2)
println("Number of iterations : ", k)
println("π = ", BigFloat(pi))
println("error = ", abs(S2-π))

Running:

% time julia pi.jl
Number of significant digits expected ?7
Value of pi calculated with Madhava formula : 3.1415926585894076
Number of iterations : 199999997
π = 3.141592653589793238462643383279502884197169399375105820974944592307816406286198
error = 4.999614500178495e-9

real	0m11,298s
user	0m10,658s
sys	0m0,855s

The key thing here is that Julia compiles the functions to specialized machine code the first time they are executed. The compilation time here is negligible. If one writes the code outside functions, however, the code is not specialized and runs interpreted, which is slow.

4 Likes
6

Many Thanks for your answer ! :slight_smile:

7

Thanks ! :slight_smile:

8

Many thanks for your explanation ! :slight_smile:

9

@peDurand welcome! Please ask more questions - the community here is very friendly.

Also please check out Pluto notebooks. They might be useful in your teaching.

7 Likes
10

Slightly simpler:

function main(nb)
    S1 = 0.0
    S2 = 1.0
    k = 0
    eps = 0.1^(nb+1)

    while abs(S1-S2) >= eps
        S2  = S1
        S1 += 4 * ((-1)^k)/(2k + 1)
        k  += 1
    end

    println("Value of pi calculated with Madhava formula : ", S2)
    println("Number of iterations : ", k)
    println("π = ", BigFloat(pi))
    println("error = ", abs(S2-π))
end
2 Likes
11

And a version that is much, much faster:

function main(nb)
    S1 = 0.0
    S2 = 1.0
    k = 0
    eps = 0.1^(nb+1)
    fac = 1

    while abs(S1-S2) >= eps
        S2  = S1
        S1 += 4 * fac/(2k + 1)
        k  += 1
        fac *= -1
    end

    println("Value of pi calculated with Madhava formula : ", S2)
    println("Number of iterations : ", k)
    println("π = ", BigFloat(pi))
    println("error = ", abs(S2-π))
end

julia> @time main(7)
Value of pi calculated with Madhava formula : 3.1415926585894076
Number of iterations : 199999998
π = 3.141592653589793238462643383279502884197169399375105820974944592307816406286198
error = 4.999614500178495e-9
  0.751861 seconds (52 allocations: 6.013 KiB)

The power operator in Julia is slow in this case… Perhaps a fast path for (-1)^n if n is an integer should be added?

16 Likes
12

Not adding much over the user above me, but another technique one can use with these alternating sums is to simply unroll them by pairs. Then one doesn’t have to deal with alternating signs on variables and can simply hard-code them.

I also go further in precomputing the number of iterations (maybe with a small excess, but no big deal with this many – also this stopping condition is less prone to roundoff error). Then I transform this into the equivalent sum and use Julia’s sum function.

Unroll 2x 
function calcpi_1(nb)
    S1 = 0.0
    S2 = 1.0
    k = 0
    eps = 0.1^(nb+1)

    while abs(S1-S2) >= eps
        # S2  = S1
        S1 += 4 / (2k + 1)
        k  += 1
        S2  = S1
        S1 -= 4 / (2k + 1)
        k  += 1
    end

    println("Value of pi calculated with Madhava formula : ", S2)
    println("Number of iterations : ", k)
    println("π = ", BigFloat(pi))
    println("error = ", abs(S2-π))
    return S2
end
Precompute number of iterations 
function calcpi_2(nb)
    S1 = 0.0
    S2 = 1.0
    # 4*halfkmax + 3 > 4 * 10^(nb+1)
    halfkmax = 10^(nb+1) # 4 / (4halfkmax + 3) < 10^-(nb+1)
    for halfk in 0:halfkmax
        S2 = S1 + 4 / (4halfk + 1)
        S1 = S2 - 4 / (4halfk + 3)
    end

    println("Value of pi calculated with Madhava formula : ", S2)
    println("Number of iterations : ", 2+2halfkmax)
    println("π = ", BigFloat(pi))
    println("error = ", abs(S2-π))
    return S2
end
Simplify and convert to sum 
function calcpi_3(nb)
    # 4*halfkmax + 3 > 4 * 10^(nb+1)
    halfkmax = 10^(nb+1) # 4 / (4halfkmax + 3) < 10^-(nb+1)

    # 4 / (4halfk + 1) - 4 / (4halfk - 1) == -8 / (16halfk^2 - 1)
    S2 = 4.0 + sum(halfk -> -8 / ((4.0*halfk)^2 - 1), 1:halfkmax)

    println("Value of pi calculated with Madhava formula : ", S2)
    println("Number of iterations : ", 2+2halfkmax)
    println("π = ", BigFloat(pi))
    println("error = ", abs(S2-π))
    return S2
end

These are not one-to-one comparisons with your implementations in other languages (be sure to compare like algorithms). Here’s a performance comparison:

julia> @time main(7) # by ufechner7
Value of pi calculated with Madhava formula : 3.1415926585894076
Number of iterations : 199999998
π = 3.141592653589793238462643383279502884197169399375105820974944592307816406286198
error = 4.999614500178495e-9
  0.237817 seconds (52 allocations: 6.013 KiB)

julia> @time calcpi_1(7); # unroll
Value of pi calculated with Madhava formula : 3.1415926585894076
Number of iterations : 199999998
π = 3.141592653589793238462643383279502884197169399375105820974944592307816406286198
error = 4.999614500178495e-9
  0.203474 seconds (52 allocations: 6.013 KiB)

julia> @time calcpi_2(7); # precompute number of iterations
Value of pi calculated with Madhava formula : 3.1415926585894076
Number of iterations : 200000002
π = 3.141592653589793238462643383279502884197169399375105820974944592307816406286198
error = 4.999614500178495e-9
  0.191706 seconds (52 allocations: 6.013 KiB)

julia> @time calcpi_3(7); # simplify series and use sum
Value of pi calculated with Madhava formula : 3.141592658589793
Number of iterations : 200000002
π = 3.141592653589793238462643383279502884197169399375105820974944592307816406286198
error = 4.999999969612645e-9
  0.065315 seconds (53 allocations: 6.310 KiB)

I suspect that the sum version is using vectorization (aka SIMD) to further accelerate this computation. The result is a 3x speedup but a slight change to the final result.

2 Likes
13

Vectorization I don’t know. But sum uses Kahan pairwise summation (thanks for the correction @mikmoore) if possible and thus gives more accurate results than a naive, linear summation.

14

This code is as fast as my Fortran code on my computer !
I am very impressed ! Many thanks ! :slight_smile:

2 Likes
15

I have to learn how to work under Julia REPL. Usually, I prefer using an editor, Geany is my favourite, and executing codes under Geany or bash. But the results of this way of using Julia REPL is amazing ! I still have a lot to learn. Thank you very much to have shown me the power of Julia ! Now, I know that Julia is a serious option for my next courses of computational physics ! :slight_smile:
Thanks a lot !

4 Likes
16

Thank you very much for your answer. I have to admit that these things are a little bit complicated for me. I am just beginning to learn Julia, but I will try to understand. Thanks a lot !

1 Like
17

No call to Base Julia’s sum currently uses Kahan summation, but it does use pairwise summation in many cases. Pairwise summation does give accuracy benefits (error growth O(\log{N})) over serial (O(N)) but is not as accurate as Kahan summation (O(1)).

But for iterators (for now) and some other complicated cases, it will fall back to serial summation.

I strongly suspect that the sum in this case is using pairwise. I don’t see why it would need to fall back to serial, although I’d have to check the generated code to be sure.

3 Likes
18

Thanks a lot to all of you who took the time to explain me how to obtain a fast code under Julia. That’s very kind of you ! I learn a lot with your explanations.
For the past years, I’ve try several time to work with Julia. But I had problems with installations of some packages, my codes where slow … So, I was doubting of the power of Julia … I think that your help will be decisive in motivating me to seriously learn Julia in order to use it for my courses. !
Many many thanks ! :slight_smile:

4 Likes
19

It’s not totally obvious that it would be, but cospi can be used for this calculation (although returning a float rather than integer, which isn’t a problem in this case since the result is immediately converted to a float) and is efficient:

cospi(x::Integer) = isodd(x) ? -one(float(x)) : one(float(x))
3 Likes
20

just

f(n)=iseven(n) ? 1 : -1
is much faster than (-1)^n

4 Likes