Julia's assignment behavior differs from Fortran?

General Usage
50

I will post an example of the use of immutable types to compare with a Fortran routine. I do this to learn things myself, so why not share… I think this is an example on how the Julia syntax is nice when we get used to it:

Problem: Given the positions of N particles in 3D, update the positions of a subset of size M of these particles defined in a list.

Using Julia, I will first define an immutable struct to contain the coordinates, and a constructor that updated the particles with the transformation desired (in this case, multiply the coordinates by 2):

struct P{T}
  x :: T
  y :: T
  z :: T
end
update_position( x :: P ) = P( 2*x.x, 2*x.y, 2*x.z )

Let us generate a vector of these particles:

n = 1_000_000
x = [ P(1.,1.,1.) for i in 1:n ]

For the Fortran routine, the natural way will be to create a matrix, so we do:

y = ones(3,n)

We will create a list of indexes of the particles that should have their positions updated, with:

m = 1_000
list = rand(1:n,m)

Then, the Julia function that updates the positions of these m particles is:

function f1!(x,list)
  for i in list
    x[i] = update_position(x[i])
  end
end

The Fortran routines to do the same job are:

subroutine update_position(i,n,x)
  implicit none
  integer :: j, n
  integer*8 i
  double precision :: x(3,n)
  do j = 1, 3
    x(j,i) = 2*x(j,i)
  end do
end subroutine update_position

subroutine f1(n,x,m,list)
  implicit none
  integer :: i, j, n, m
  integer*8 :: list(m)
  double precision :: x(3,n)
  do i = 1, m
    call update_position(list(i),n,x)
  end do
end subroutine f1

I think that it is pretty clear that the Julia code is much simpler and clear. Bugs turn out to be easier to find at the end, even with the more permissive syntax of Julia. Of course we can suggest changes in both codes, but in particular if we want to keep the function that updates the position of a single particle separated from the one that updates the lists, the Fortran code code will not be much simpler than that.

Finally, the Julia function turns out to be even slightly faster than the Fortran routine:

 Julia with P:   1.869 μs (0 allocations: 0 bytes)
 Fortran:   2.292 μs (0 allocations: 0 bytes)

(although I know that benchmarks are always complicated and a disclaimer: if we wanted to update the positions of all particles, which are contiguous in memory, than we should use a matrix of particles in row-order, allowing loop vectorizations. The case here is that the positions to be updated are not contiguous in memory).

The complete code is here:

Summary 
using Test
using BenchmarkTools

# Using user defined immutable type 
struct P{T}
  x :: T
  y :: T
  z :: T
end
update_position( x :: P ) = P( 2*x.x, 2*x.y, 2*x.z )

function f1!(x,list)
  for i in list
    x[i] = update_position(x[i])
  end
end

# Calling Fortran routine:
"""
subroutine update_position(i,n,x)
  implicit none
  integer :: j, n
  integer*8 i
  double precision :: x(3,n)
  do j = 1, 3
    x(j,i) = 2*x(j,i)
  end do
end subroutine update_position

subroutine f1(n,x,m,list)
  implicit none
  integer :: i, j, n, m
  integer*8 :: list(m)
  double precision :: x(3,n)
  do i = 1, m
    call update_position(list(i),n,x)
  end do
end subroutine f1   

Compiled with:
gfortran -O3 mutate.f90 -shared -o mutate.so
"""
f2!(x,list) = ccall((:f1_,"mutate.so"),Nothing,
                    (Ref{Int32},         Ref{Float64}, Ref{Int32},          Ref{Int64}),
                     Int32(size(x,2)),   x,            Int32(size(list,1)), list)     


# Running

n = 1_000_000
m = 1_000
list = rand(1:n,m)

# With julia immutable struct
x = [ P(1.,1.,1.) for i in 1:n ]
f1!(x,list)

# Fortran Matrix
y = ones(3,n) 
f2!(y,list)

# Testing if the fortran code returns the same thing
@test begin
  t = true
  for i in 1:n
    if x[i].x != y[1,i] || x[i].y != y[2,i] || x[i].z != y[3,i] 
      t = false
    end
  end
  t
end

setup = ( x = [ P(1.,1.,1.) for i in 1:n ] )
print(" Julia with P: "); @btime f1!(x,$list) setup = setup evals = 1

setup = ( y = ones(3,n) )
print(" Fortran: "); @btime f2!(y,$list) setup=setup evals=1
3 Likes