# Julia's assignment behavior differs from Fortran?

**URL:** <https://discourse.julialang.org/t/julias-assignment-behavior-differs-from-fortran/50389>\
**Category:** General Usage\
**Tags:** fortran\
**Created:** [November 18, 2020, 1:24pm UTC](https://discourse.julialang.org/t/julias-assignment-behavior-differs-from-fortran/50389 "2020-11-18T13:24:15Z")\
**Posts on this page:** 9\
**Page:** 3

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [November 18, 2020, 7:30pm UTC](https://discourse.julialang.org/t/julias-assignment-behavior-differs-from-fortran/50389/47 "2020-11-18T19:30:47Z")

</div>

`=` is not an equality test, but an assignment. Confusing, perhaps. But entirely consistent across many languages.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [November 18, 2020, 10:34pm UTC](https://discourse.julialang.org/t/julias-assignment-behavior-differs-from-fortran/50389/48 "2020-11-18T22:34:06Z")

</div>

When I learned that the syntax `x[i] = v` is converted to `setindex!(x, i, v)` in Julia, I didn’t need to think too much about what is mutable or not, it is pretty clear that something is being set to a new value. The syntax `x = v` is a change of name as others have already pointed out.

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [November 18, 2020, 10:57pm UTC](https://discourse.julialang.org/t/julias-assignment-behavior-differs-from-fortran/50389/49 "2020-11-18T22:57:04Z")

</div>

> [@lmiq](#):
>
> However, the values occupying these types of memory positions do not have an address in the usual sense. You cannot change the value associated to that position in memory because that value in that position in memory is not persistent, that is, it will be discarded as soon as possible.

FWIW, Fortran compilers normally don’t have a problem putting all the arrays you’re mutating on the stack. For example, see the gfortran flag [-fstack-arrays](https://gcc.gnu.org/onlinedocs/gfortran/Code-Gen-Options.html#Code-Gen-Options).  
In Julia, `StaticArrays.MArrays` will also generally be placed on the stack, so long as they do not escape the function creating them or get passed as an argument to a function that isn’t inlined.

A feature request I have in Julia (maybe I’ll make a PR) is support for `alloca`, so it’s easier to stack allocate memory that’s freed when the function returns, but that you can pass as an argument to non-inlined functions. `llvmcall` doesn’t work for this, because the memory is immediately freed when the llvmcall returns, not when the Julia function calling llvmcall returns.

```julia
julia> using StaticArrays, BenchmarkTools

julia> function f()
         s = 0.
         for i in 1:1000
           x = SVector{3,Float64}(i, sqrt(i), i^2)
           for j in 1:3
             s = s + x[j]
           end
         end
         s
       end
f (generic function with 1 method)

julia> function g()
         s = 0.
         for i in 1:1000
           x = MVector{3,Float64}(undef)
           x[1] = i
           x[2] = sqrt(i)
           x[3] = i^2
           for j in 1:3
             s = s + x[j]
           end
         end
         s
       end
g (generic function with 1 method)

julia> @btime f()
  2.605 μs (0 allocations: 0 bytes)
3.343550974558874e8

julia> @btime g()
  2.605 μs (0 allocations: 0 bytes)
3.343550974558874e8

```

However, `SVector`s are easier for the compiler; using the `MVector` can prevent some optimizations; see what happens when we add `@fastmath`:

```julia
julia> @fastmath function f()
         s = 0.
         for i in 1:1000
           x = SVector{3,Float64}(i, sqrt(i), i^2)
           for j in 1:3
             s = s + x[j]
           end
         end
         s
       end
f (generic function with 1 method)

julia> @fastmath function g()
         s = 0.
         for i in 1:1000
           x = MVector{3,Float64}(undef)
           x[1] = i
           x[2] = sqrt(i)
           x[3] = i^2
           for j in 1:3
             s = s + x[j]
           end
         end
         s
       end
g (generic function with 1 method)

julia> @btime f()
  744.343 ns (0 allocations: 0 bytes)
3.3435509745588744e8

julia> @btime g()
  2.604 μs (0 allocations: 0 bytes)
3.343550974558874e8

```

Suddenly `f` (using `SVector`s) is much faster, but the speed of `g()` is almost unchanged.

But my point here is that, just because the `MVector` is mutable and `getindex` and `setindex!` are [defined using loads and stores to/from their pointers](https://github.com/JuliaArrays/StaticArrays.jl/blob/8033332b1b905cb3a10998b0d42e0f90034f3f3c/src/MArray.jl#L80) doesn’t mean it has to be heap allocated.

EDIT:  
The simplest way to speed up `g` is to use `@inbounds`:

```julia
julia> @fastmath function h()
         s = 0.
         @inbounds for i in 1:1000
           x = MVector{3,Float64}(undef)
           x[1] = i
           x[2] = sqrt(i)
           x[3] = i^2
           for j in 1:3
             s = s + x[j]
           end
         end
         s
       end
h (generic function with 1 method)

julia> @btime h()
  744.317 ns (0 allocations: 0 bytes)
3.3435509745588744e8

```

A less simple way is to manually unroll the inner loop:

```julia
julia> @fastmath function g()
         s = 0.
         for i in 1:1000
           x = MVector{3,Float64}(undef)
           x[1] = i
           x[2] = sqrt(i)
           x[3] = i^2
           Base.Cartesian.@nexprs 3 j -> begin
             s = s + x[j]
           end
         end
         s
       end
g (generic function with 1 method)

julia> @btime g()
  744.324 ns (0 allocations: 0 bytes)
3.3435509745588744e8

```

I guess it’s close to a limit where the compiler gives up optimizing it.

EDIT:  
It just occurred to me that AVX512 is probably needed for the `@fastmath` versions of this function to be fast, because AVX512 is needed for SIMD `Int64 -> Float64` conversion.  
If you don’t have AVX512, I’d recommend

```julia
julia> @fastmath function h2()
         s = 0.; i = 0.0
         @inbounds for _ in 1:1000
           x = MVector{3,Float64}(undef);
           i += 1
           x[1] = i
           x[2] = sqrt(i)
           x[3] = i^2
           for j in 1:3
             s = s + x[j]
           end
         end
         s
       end
h2 (generic function with 1 method)

julia> @btime h2()
  711.417 ns (0 allocations: 0 bytes)
3.3435509745588744e8

```

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [November 19, 2020, 12:39pm UTC](https://discourse.julialang.org/t/julias-assignment-behavior-differs-from-fortran/50389/50 "2020-11-19T12:39:57Z")

</div>

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):

```julia
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:

```julia
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:

```julia
y = ones(3,n)

```

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

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

```

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

```julia
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:

```fortran
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
 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**
>
> ```julia
> 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
> 
>                                                  
> 
> ```

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [November 19, 2020, 3:56pm UTC](https://discourse.julialang.org/t/julias-assignment-behavior-differs-from-fortran/50389/51 "2020-11-19T15:56:37Z")

</div>

On computers with AVX, you can improve performance by padding out `P` to be 256 bits, i.e. make it 4D but just ignore the last.  
Additionally, gfortran generally needs the compiler flag `-fno-semantic-interposition` for performance. Without it, it will not inline functions.

I also added `@inbounds` to the Julia version to make it equivalent to the Fortran.

> **Code**
>
> ```fortran
> module mutate
> 
> use ISO_C_BINDING
> 
> contains
> subroutine update_position(x)
> implicit none
> real(C_double), dimension(:), intent(inout) :: x
> x = 2*x
> end subroutine update_position
> 
> subroutine f1(n,x,m,list) BIND(C, name="f1")
> implicit none
> integer(C_long) :: i, n, m
> integer(C_long) :: list(m)
> real(C_double) :: x(3,n)
> do i = 1, m
> call update_position(x(:,list(i)))
> end do
> end subroutine f1
> 
> subroutine f2(n,x,m,list) BIND(C, name="f2")
> implicit none
> integer(C_long) :: i, n, m
> integer(C_long) :: list(m)
> real(C_double) :: x(4,n)
> do i = 1, m
> call update_position(x(:,list(i)))
> end do
> end subroutine f2
> 
> end module mutate
> 
> ```
> 
> I compiled with  
> `gfortran -Ofast -march=native -Wall -fno-semantic-interposition -shared -fPIC mutate.f90 -o libmutate.so`
> 
> ```julia
> using Test, BenchmarkTools
> 
> struct P{T}
> x :: T
> y :: T
> z :: T
> end
> update_position( x :: P ) = P( 2*x.x, 2*x.y, 2*x.z )
> 
> struct Ppad{T}
> x :: T
> y :: T
> z :: T
> p :: T
> end
> function Ppad(_x, _y, _z)
> x, y, z = promote(_x, _y, _z)
> Ppad(x, y, z, zero(x))
> end
> update_position( x :: Ppad ) = Ppad( 2x.x, 2x.y, 2x.z, 2x.p )
> 
> function f1!(x,list)
> @inbounds for i in list
> x[i] = update_position(x[i])
> end
> end
> 
> const LIBMUTATE = "libmutate.so"
> f2_3!(x,list) = ccall((:f1, LIBMUTATE), Nothing,
> (Ref{Clong}, Ref{Float64}, Ref{Clong}, Ref{Clong}),
> size(x,2), x, size(list,1), list)     
> f2_4!(x,list) = ccall((:f2, LIBMUTATE), Nothing,
> (Ref{Clong}, Ref{Float64}, Ref{Clong}, Ref{Clong}),
> size(x,2), x, size(list,1), list)     
> f2!(x, list) = size(x,1) == 3 ? f2_3!(x, list) : f2_4!(x, list)
> 
> n = 1_000_000
> m = 1_000
> list = rand(1:n,m);
> 
> setup = ( x = [P(1.,1.,1.) for i in 1:n] );
> print(" Julia with P: "); @btime f1!(x,$list) setup = setup evals = 1
> setup = ( x = [Ppad(1.,1.,1.) for i in 1:n] );
> print(" Julia with Ppad: "); @btime f1!(x,$list) setup = setup evals = 1
> 
> setup = ( y = ones(3,n) );
> print(" Fortran: "); @btime f2!(y,$list) setup=setup evals=1
> setup = ( y = ones(4,n) );
> print(" Fortran: "); @btime f2!(y,$list) setup=setup evals=1
> 
> ```

This yields:

```julia
julia> setup = ( x = [P(1.,1.,1.) for i in 1:n] );

julia> print(" Julia with P: "); @btime f1!(x,$list) setup = setup evals = 1
 Julia with P: 978.000 ns (0 allocations: 0 bytes)

julia> setup = ( x = [Ppad(1.,1.,1.) for i in 1:n] );

julia> print(" Julia with Ppad: "); @btime f1!(x,$list) setup = setup evals = 1
 Julia with Ppad: 752.000 ns (0 allocations: 0 bytes)

julia> setup = ( y = ones(3,n) );

julia> print(" Fortran: "); @btime f2!(y,$list) setup=setup evals=1
 Fortran: 983.000 ns (0 allocations: 0 bytes)

julia> setup = ( y = ones(4,n) );

julia> print(" Fortran: "); @btime f2!(y,$list) setup=setup evals=1
 Fortran: 852.000 ns (0 allocations: 0 bytes)

```

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [November 19, 2020, 4:26pm UTC](https://discourse.julialang.org/t/julias-assignment-behavior-differs-from-fortran/50389/52 "2020-11-19T16:26:00Z")

</div>

Interesting, that padding stuff is advanced stuff for me 🙂

Concerning the compiler flags, I get the best performance (the same as that with your flags) with

`-O3 -march=native`

Another thing: the `-Ofast` option, as I checked, included `fastmath`. Probably one should compare with the loop in Julia using `@fastmath`, but I didn’t see any important difference here.

By the way, thank you very much for the example on how to call Fortran routines from Julia. I struggled a lot to find good examples of that, the manual is quite short in examples showing both codes.

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [November 19, 2020, 5:34pm UTC](https://discourse.julialang.org/t/julias-assignment-behavior-differs-from-fortran/50389/53 "2020-11-19T17:34:46Z")

</div>

> [@lmiq](#):
>
> Concerning the compiler flags, I get the best performance (the same as that with your flags) with

Does this mean you get better performance without `-fno-semantic-interposition`?  
I guess it depends on compiler options. [Godbolt](https://godbolt.org/z/obsGe7) shows the functions inlining regardless of `-fno-semantic-interposition`, but when I compile locally (checking asm with `-S -masm=intel`), the functions do not inline unless I add `-fno-semantic-interposition`.

Fast math makes no difference there; associativity isn’t needed for optimizations here.  
I just didn’t bother editing it when jumping up in my history to the last `gfortran` invocation.

> [@lmiq](#):
>
> Interesting, that padding stuff is advanced stuff for me 🙂

SIMD (Single Instruction Multiple Data) lets the CPU operate on multiple numbers per instruction.  
With `P`,

```julia
julia> xp = [Ppad(1.,1.,1.) for i in 1:n];

julia> x = [P(1.,1.,1.) for i in 1:n];

julia> @code_llvm debuginfo=:none f1!(x, list)

julia> @code_llvm debuginfo=:none f1!(xp, list)

```

This yields, without padding:

```llvm
L9: ; preds = %L9, %L9.lr.ph
  %value_phi14 = phi i64 [0, %L9.lr.ph], [%21, %L9]
  %14 = getelementptr inbounds i64, i64* %12, i64 %value_phi14
  %15 = load i64, i64* %14, align 8
  %16 = add i64 %15, -1
  %.elt = getelementptr inbounds [3 x double], [3 x double]* %13, i64 %16, i64 0
  %17 = bitcast double* %.elt to <2 x double>*
  %18 = load <2 x double>, <2 x double>* %17, align 8
  %.elt7 = getelementptr inbounds [3 x double], [3 x double]* %13, i64 %16, i64 2
  %.unpack8 = load double, double* %.elt7, align 8
  %19 = fmul <2 x double> %18, <double 2.000000e+00, double 2.000000e+00>
  %20 = fmul double %.unpack8, 2.000000e+00
  store <2 x double> %19, <2 x double>* %17, align 8
  store double %20, double* %.elt7, align 8
  %21 = add nuw nsw i64 %value_phi14, 1
  %exitcond.not = icmp eq i64 %21, %8
  br i1 %exitcond.not, label %L28, label %L9

```

It performs 2 loads:

```llvm
  %18 = load <2 x double>, <2 x double>* %17, align 8
  %.unpack8 = load double, double* %.elt7, align 8

```

Two multiplies (these later get turned into adds):

```llvm
  %19 = fmul <2 x double> %18, <double 2.000000e+00, double 2.000000e+00>
  %20 = fmul double %.unpack8, 2.000000e+00

```

and two stores:

```llvm
  store <2 x double> %19, <2 x double>* %17, align 8
  store double %20, double* %.elt7, align 8

```

One of these is of `<2 x double>`, and the other a single double.

With padding:

```llvm
L9: ; preds = %L9, %L9.lr.ph
  %value_phi19 = phi i64 [0, %L9.lr.ph], [%20, %L9]
  %14 = getelementptr inbounds i64, i64* %12, i64 %value_phi19
  %15 = load i64, i64* %14, align 8
  %16 = add i64 %15, -1
  %.elt = getelementptr inbounds [4 x double], [4 x double]* %13, i64 %16, i64 0
  %17 = bitcast double* %.elt to <4 x double>*
  %18 = load <4 x double>, <4 x double>* %17, align 8
  %19 = fmul <4 x double> %18, <double 2.000000e+00, double 2.000000e+00, double 2.000000e+00, double 2.000000e+00>
  store <4 x double> %19, <4 x double>* %17, align 8
  %20 = add nuw nsw i64 %value_phi19, 1
  %exitcond.not = icmp eq i64 %20, %8
  br i1 %exitcond.not, label %L31, label %L9

```

We now have just a single load:

```llvm
  %18 = load <4 x double>, <4 x double>* %17, align 8

```

A single multiply:

```llvm
  %19 = fmul <4 x double> %18, <double 2.000000e+00, double 2.000000e+00, double 2.000000e+00, double 2.000000e+00>

```

and a single store:

```llvm
  store <4 x double> %19, <4 x double>* %17, align 8

```

So on CPUs with the AVX instruction set (which provides 256 bit vectors, letting them operate on `<4 x double>`s at a time), padding to add an extra operation will actually decrease the number of instructions that must be performed: instead of `2 + 1 = 3`, we can do `4` at a time.  
Of course, there are a lot of other instructions aside from just these, so the actual performance difference is smaller.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [November 19, 2020, 5:54pm UTC](https://discourse.julialang.org/t/julias-assignment-behavior-differs-from-fortran/50389/54 "2020-11-19T17:54:45Z")

</div>

> [@Elrod](#):
>
> Does this mean you get better performance without `-fno-semantic-interposition` ?  
> I guess it depends on compiler options. [Godbolt](https://godbolt.org/z/obsGe7)

No, sorry, I get roughly the same performance. I could not see any difference. I didn’t know about that tool, very interesting.

ps. By chance do you teach all this stuff somewhere? The way you approach these problems is a complete new world to me.

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [November 20, 2020, 6:14am UTC](https://discourse.julialang.org/t/julias-assignment-behavior-differs-from-fortran/50389/55 "2020-11-20T06:14:55Z")

</div>

> [@lmiq](#):
>
> ps. By chance do you teach all this stuff somewhere? The way you approach these problems is a completely new world to me.

No, but I’ll work on some material eventually. I’ll probably write some for LoopVectorization’s eventual v1.0 release (which is months away).

[Previous page](https://discourse.julialang.org/t/julias-assignment-behavior-differs-from-fortran/50389.md?page=2)