Loops in different languages

[jcandy]

These results are unexpected. Am I doing something wrong in Julia?

Results are in in order of fastest to slowest (fortran and C++ compiled with -O3):

PYTHON
answer = 0.06344748845410242
time = 1.243s user + 0.792s sys

import numpy as np

n = 16000
m = np.zeros([n,n])
i = np.arange(n)+1
j = np.arange(n)+1

m[:,:] = np.log(
    1+abs(
        np.outer(np.sin(i[:]),np.cos(j[:]))
    )
)

print(m[n//2-1,n//2-1])

FORTRAN
answer = 6.3447488454102419E-002
time = 2.919s user + 0.441s sys

program ftest
  
  real, dimension(:,:), allocatable :: m
  integer :: n

  n=16000
  allocate(m(n,n))

  do i=1,n
     do j=1,n
        m(i,j) = log(1+abs(sin(i*1.0)*cos(j*1.0)))
     enddo
  enddo

  print *,m(n/2,n/2)
  
end program ftest

C++ 
answer = 0.0634474884541024
time = 3.922s user + 0.535s sys

#include <iostream>
#include <cmath>
using namespace std;

int main()
{
  const int n = 16000;
  auto *m = new double [n][n];
    
  for(int i=0; i<n; i++)   
    {
      for(int j=0; j<n; j++) 
	{
	  m[i][j] = log(1.0+abs(sin(i+1.0)*cos(j+1.0)));
	}
    }

  cout.precision(16);
  cout << fixed << m[n/2-1][n/2-1] << endl;

}

JULIA
answer = 0.06344748845410242
time = 21.794s user + 1.308s sys

n=16000
m = Array{Float64}(undef,n,n)

for i in 1:n
   for j in i:n
      m[i,j] = log(1+abs(sin(i)*cos(j)))
   end
end

println(m[div(n,2),div(n,2)])

[orso82]

Yes, you are working in the global scope. See the first point of https://docs.julialang.org/en/v1/manual/performance-tips/

try this:

function f(n::Int)
    m = Array{Float64}(undef,n,n)

    for i in 1:n
       for j in 1:n
          @inbounds m[i,j] = log(1+abs(sin(i)*cos(j)))
       end
    end

    return m
end

n=16000
@time m=f(n)
println(m[div(n,2),div(n,2)])

on my laptop:

  5.776092 seconds (2 allocations: 1.907 GiB, 0.65% gc time)
0.06344748845410242

NOTE: The first time you run it it will take longer to run because it compiles the function

This variation

function f1(n::Int)
    m = Array{Float64}(undef, n, n)
    sines = sin.(1:n)
    cosines = cos.(1:n)
    
    @inbounds for j in 1:n
        for i in 1:n
            m[i, j] = log(1.0 + abs(sines[i] * cosines[j]))
        end
    end
    
    return m
end

brings it to

  1.254190 seconds (6 allocations: 1.908 GiB, 2.84% gc time)
0.06344748845410242

This is actually faster mainly because of the reduced calls to sin and cos and the contiguous access in memory to the matrix m

[bclyons12]

@jcandy One thing to note is that Julia (and Fortran) are column major, so they should be accessed with the leftmost index as the inner loop, as Orso does here. C++ is the opposite.

[jcandy]

Edit: Switching loop order increases julia time by about 50%, but only by about 10% in both fortran and C++ perhaps because of the optimizer.

[orso82]

This vectorized implementation is actually also ok, although it uses more memory

function f(n::Int)    
    return log.(1.0 .+ abs.(sin.(1:n) * cos.((1:n)')))
end

n=16000
@time m=f(n)
println(m[div(n,2),div(n,2)])

  1.422594 seconds (10 allocations: 3.815 GiB, 1.53% gc time)
0.06344748845410242

[orso82]

I asked ChatGPT why the two have such different speed:

---

There are a few reasons why the first function you provided (which uses broadcasting in Julia) is faster than the second (which uses nested loops):

Vectorized Computations: In the first function, you are making use of Julia's broadcasting capabilities (. operator). Julia's broadcasting mechanism is well-optimized and allows for operations on arrays without the explicit use of loops. Broadcasting effectively transforms the operations into highly efficient looped code behind the scenes. It leverages the memory layout and cache efficiently.

Memory Allocation: The first function directly computes the result without any preallocation, whereas in the second function, you first allocate an array of size n x n and then fill it in using a nested loop. This preallocation and then assignment can be slower than directly computing and placing values, especially when n is large.

Inherent Parallelism: Broadcasting in Julia is designed to benefit from inherent parallelism. Under the hood, the compiler may opt to parallelize certain broadcasting operations if it detects that it can provide a speed-up, which can result in faster computations.

Function Calls Overhead: In the loop-based version, you're making multiple function calls (sin, cos, abs, and log) in each iteration of the nested loops. Even though these function calls are fast, the overhead adds up when you have n^2 iterations. In the broadcasting version, the operations are more "batched" which can reduce this overhead.

That being said, it's worth noting that modern CPUs are pretty good at executing tight loops, so the difference may not always be substantial, especially for small values of n. But as n grows, the benefits of the first function's approach should become more pronounced.

[orso82]

In your example the numpy implementation is actually fast because it uses BLAS underneath. The Python implementation with for loops takes 17 minutes to run.

import numpy as np
import time

def f(n):
    m = np.zeros([n,n])
    i = np.arange(n)+1
    j = np.arange(n)+1

    for i in np.arange(n)+1:
       for j in np.arange(n)+1:
          m[i-1,j-1] = np.log(1+np.abs(np.sin(i)*np.cos(j)))

    return m

n=16000
start = time.time()
m = f(n)
print(m[n//2-1,n//2-1])
end = time.time()
print(end - start)

[jcandy]

@orso82

Thank you. Timings improve with Julia inside a function (not global scope). Of course defining the vectors (as I did for numpy) is cheating for the purposes of this test and will also speed up fortran and C++.

JULIA function
answer = 0.06344748845410242
time = 6.480s user + 0.982s sys

function f1(n::Int)
   m = Array{Float64}(undef,n,n)

   for i in 1:n
      for j in 1:n
         m[i,j] = log(1+abs(sin(i)*cos(j)))
      end
   end
   return m
end

n = 16000
m = f1(n)
println(m[div(n,2),div(n,2)])

[bclyons12]

Did you mean to write j in i:n instead of j in 1:n? I don't think your other code snippets do that.

[jcandy]

Good catch! I am blind. I will edit ...

[bclyons12]

I was also getting this function to be very slow. After some research, I was directed to the LoopVectorization.@avx:

This library provides the @avx macro, which may be used to prefix a for loop or broadcast statement. It then tries to vectorize the loop to improve runtime performance.The macro assumes that loop iterations can be reordered. It also currently supports simple nested loops, where loop bounds of inner loops are constant across iterations of the outer loop, and only a single loop at each level of loop nest. These limitations should be removed in a future version.

Anyway, I changed the code to:

using LoopVectorization
function f4(n::Int)
   m = Array{Float64}(undef,n,n)
   @avx for j in 1:n
      for i in 1:n
         m[i,j] = log(1+abs(sin(i)*cos(j)))
      end
   end
   return m
end
n = 16000
@time m = f4(n)
println(m[div(n,2),div(n,2)]);

And now get:

1.488121 seconds (2 allocations: 1.907 GiB, 15.25% gc time)
0.06344748845410242

[jcandy]

Does this use multithreading? It seems like the startup is using multiple threads even though the evaluation isn't (?)

For comparison, I get user+sys = 2.7s+1.0s (fortran is 2.9s+0.36s)

[bclyons12]

Perhaps. When I run on one thread I get:

@time m = f4(n);
  2.743712 seconds (2 allocations: 1.907 GiB, 1.59% gc time)

[bclyons12]

Digging into this, no, it was not using multithreading. Running it many times on 1 thread gives an average time of 1.44 seconds on 1 thread:

@benchmark f4($n)
BenchmarkTools.Trial: 20 samples with 1 evaluation.
 Range (min … max):  1.305 s …    1.878 s  ┊ GC (min … max): 0.03% … 6.75%
 Time  (median):     1.429 s               ┊ GC (median):    8.09%
 Time  (mean ± σ):   1.440 s ± 132.487 ms  ┊ GC (mean ± σ):  7.37% ± 2.58%

     ▄   ▁    █                                               
  ▆▆▆█▁▆▁█▁▁▁▁█▁▆▆▁▁▁▁▁▁▁▆▆▆▆▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▆ ▁
  1.3 s          Histogram: frequency by time         1.88 s <

 Memory estimate: 1.91 GiB, allocs estimate: 2.

[jcandy]

The multithreading/parallelism happens right at startup, which I suppose means the compilation is parallel

[bclyons12]

This seems like weird behavior to me:

using BenchmarkTools
g(i, j) = log(1+abs(sin(i)*cos(j)))
@btime g(14798, 2166)
i = 14798
j = 2166
@btime g($i,$j)

 1.140 ns (0 allocations: 0 bytes)
  50.681 ns (0 allocations: 0 bytes)

Why is it so much slower when using variables?

[jcandy]

For comparison I would use the unix timer. Not sure if I trust the julia time for apples-to-apples comparison

[bclyons12]

I'm not sure how to use that.

[jcandy]

$ time <yourprogram>

[orso82]

you cannot use the unix time with Julia, since the first execution is always slower, because it compiles the code. Subsequent calls will be faster.

[jcandy]

Yes, I understand there is a startup cost, which you can measure if you run with n=16 instead of 16k.

[jguterl]

The following function gives way faster performances(on apple hardware, for Intel looks for MKL.jl or IntelVectorMath.jl) function test_vec_(v_n::Vector{Float64}, m::Array{Float64,2}) vsin = similar(v_n) vcos = similar(v_n) AppleAccelerate.sin!(vsin, v_n) AppleAccelerate.cos!(vcos, v_n) m .= 1.0 .+ abs.(vsin .* vcos') AppleAccelerate.log!(m, m) nothing end The math functions are not vectorized efficiently in Julia by default.

[bclyons12]

This won't give an apples-to-apples (no pun intended) comparison because it only calls sin and cos N times instead of N^2 times. It would be good to separately test how this compares to the equivalent Fortran code though, and the equivalent Julia code with @avx

[jcandy]

Note that we can call the fortran from Julia:

FORTRAN:

module ffunc

contains

  double precision function f(n)
    double precision, dimension(:,:), allocatable :: m
    integer :: n

    allocate(m(n,n))
    
    do i=1,n
       do j=1,n
          m(i,j) = log(1+abs(sin(i*1d0)*cos(j*1d0)))
       enddo
    enddo

    f = m(n/2,n/2)

    deallocate(m)
 
  end function f
     
end module ffunc

SHARED OBJECT:

$ gfortran -o ffunc.so ffunc.f90 -shared -fPIC

JULIA:

m0 = Int32[16000]
result = @ccall "./ffunc.so".__ffunc_MOD_f(m0::Ptr{Int32})::Float64
println(result)

[jguterl]

Possible vectorization with avx instructions: using BenchmarkTools using LoopVectorization n = 16000 mr = rand(n, n) mi = ones(n) function test_sin(n, m::Matrix{Float64}) for j in axes(m, 2) for i in axes(m, 1) m[i, j] = log(1.0 + abs(sin(i) * cos(j))) end end end function test_sin_avx(n, m::Matrix{Float64}) @avx for j in axes(m, 2) for i in axes(m, 1) m[i, j] = log(1.0 + abs(sin(i) * cos(j))) end end end function test_sin_simd(n, m::Matrix{Float64}) @simd for j in axes(m, 2) for i in axes(m, 1) m[i, j] = log(1.0 + abs(sin(i) * cos(j))) end end end #AppleAccelerate.sin!(mr, mr) @Btime test_sin(n, mr) @Btime test_sin_avx(n, mr) @Btime test_sin_simd(n, mr)

…

------------------------------------- julia> @Btime test_sin(n, mr) 4.799 s (0 allocations: 0 bytes) julia> @Btime test_sin_avx(n, mr) 1.248 s (0 allocations: 0 bytes) julia> @Btime test_sin_simd(n, mr) 4.856 s (0 allocations: 0 bytes)

On Thu, Aug 17, 2023 at 4:22 PM Jeff Candy ***@***.***> wrote: Note that we can call the fortran from Julia: FORTRAN: module ffunc contains double precision function f(n) double precision, dimension(:,:), allocatable :: m integer :: n allocate(m(n,n)) do i=1,n do j=1,n m(i,j) = log(1+abs(sin(i*1d0)*cos(j*1d0))) enddo enddo f = m(n/2,n/2) deallocate(m) end function f end module ffunc SHARED OBJECT: $ gfortran -o ffunc.so ffunc.f90 -shared -fPIC JULIA: m0 = Int32[16000] result = @CCall "./ffunc.so".__ffunc_MOD_f(m0::Ptr{Int32})::Float64 println(result) — Reply to this email directly, view it on GitHub <#396 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AEESMZEDIMBFG5OMYCPOAS3XV2RTZANCNFSM6AAAAAA3TRC3NE> . You are receiving this because you commented.Message ID: ***@***.*** com>

[jcandy]

Thanks for the helpful feedback! Now I am thinking seriously about calling fortran from Julia.