# Julia: Fast as Fortran, Beautiful as Python
**URL:** https://fortran-lang.discourse.group/t/julia-fast-as-fortran-beautiful-as-python/1405
**Category:** Uncategorized
**Created:** [June 20, 2021, 6:21pm UTC](https://fortran-lang.discourse.group/t/julia-fast-as-fortran-beautiful-as-python/1405 "2021-06-20T18:21:13Z")
**Posts on this page:** 20
**Page:** 7
### Author:  [@themos](https://fortran-lang.discourse.group/u/themos)
#### Post date: [June 28, 2021, 11:06am UTC](https://fortran-lang.discourse.group/t/julia-fast-as-fortran-beautiful-as-python/1405/129 "2021-06-28T11:06:41Z")
Hello and welcome.
You seem to be attempting to time the SUM intrinsic. My experience with it is that Intel’s SUM() is about 3 times faster that NAG’s and GNU’s (slightly different versions, I use ifort 2021.2.0, NAG 7.0 and gfortran 8.3.1). Perhaps you need to repeat the timing and clear caches in between.
It is also possible to write Fortran that does it faster still on a multi-core machine (another 3x) and with a more accurate result.
---
### Author:  [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)
#### Post date: [June 28, 2021, 11:36am UTC](https://fortran-lang.discourse.group/t/julia-fast-as-fortran-beautiful-as-python/1405/130 "2021-06-28T11:36:01Z")
> [@themos](#):
>
> Hello and welcome.
>
> You seem to be attempting to time the SUM intrinsic. My experience with it is that Intel’s SUM() is about 3 times faster that NAG’s and GNU’s (slightly different versions, I use ifort 2021.2.0, NAG 7.0 and gfortran 8.3.1). Perhaps you need to repeat the timing and clear caches in between.
Is it possible to use C interoperability to compile source files with two different Fortran compilers to create an executable? This could be useful if you know that one compiler is faster for part of the code or if you have object code for a numerical library such as NAG or IMSL that is compiler-specific (but would calling the library from another compiler violate the license?).
---
### Author:  [@themos](https://fortran-lang.discourse.group/u/themos)
#### Post date: [June 28, 2021, 11:47am UTC](https://fortran-lang.discourse.group/t/julia-fast-as-fortran-beautiful-as-python/1405/131 "2021-06-28T11:47:41Z")
You would need two different runtime libraries to be linked in and I/O, multithreading etc. could be a nightmare. But for some simple things, you could get away with it. NAG Library products are compiler specific but sometimes the documentation says that they can be used with user code compiled with a different compiler. I would not recommend it. If you are unhappy with your vendor’s performance, write them a note!
---
### Author:  [@cmaapic](https://fortran-lang.discourse.group/u/cmaapic)
#### Post date: [June 28, 2021, 12:47pm UTC](https://fortran-lang.discourse.group/t/julia-fast-as-fortran-beautiful-as-python/1405/132 "2021-06-28T12:47:59Z")
Hi
The Fortran and C++ could also be rewritten to use openmp fairly easily.
Jane and I are writing some additional parallel examples at the moment.
Cheers
Ian
---
### Author:  [@cmaapic](https://fortran-lang.discourse.group/u/cmaapic)
#### Post date: [June 28, 2021, 1:25pm UTC](https://fortran-lang.discourse.group/t/julia-fast-as-fortran-beautiful-as-python/1405/133 "2021-06-28T13:25:58Z")
Hi Themos The intention was to look at the run time of a small example that does initialisation and summation in a variety of languages and with a variety of compilers. The intention was not to test the sum intrinsic, but rather to look at the differences between an interpreted language and some commonly used compiled languages.
---
### Author:  [@certik](https://fortran-lang.discourse.group/u/certik)
#### Post date: [June 28, 2021, 10:16pm UTC](https://fortran-lang.discourse.group/t/julia-fast-as-fortran-beautiful-as-python/1405/134 "2021-06-28T22:16:43Z")
@cmaapic thanks for your posts and welcome to the forum! I’ve been reading your compiler comparisons papers for a long time.
---
### Author:  [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)
#### Post date: [July 9, 2021, 5:20pm UTC](https://fortran-lang.discourse.group/t/julia-fast-as-fortran-beautiful-as-python/1405/135 "2021-07-09T17:20:21Z")
Reading [Concise and beautiful algorithms written in Julia | Hacker News](https://news.ycombinator.com/item?id=27783961) and browsing the associated [code](https://github.com/mossr/BeautifulAlgorithms.jl) at GitHub, I wonder how the corresponding modern Fortran algorithms would look. Having adopted a ? ternary, the next step for Fortran is Greek letters and more mathematical symbols. Just kidding.
---
### Author:  [@certik](https://fortran-lang.discourse.group/u/certik)
#### Post date: [July 9, 2021, 5:42pm UTC](https://fortran-lang.discourse.group/t/julia-fast-as-fortran-beautiful-as-python/1405/136 "2021-07-09T17:42:37Z")
It does look beautiful, although some of the unicode does not render correctly in my Firefox, so it becomes unreadable. The bigger issue is just like with Python, I struggle to understand what the code does. Take this example: [GitHub - mossr/BeautifulAlgorithms.jl: Concise and beautiful algorithms written in Julia](https://github.com/mossr/BeautifulAlgorithms.jl#thompson-sampling)
```auto
using Distributions
function thompson_sampling(𝛂, 𝛃, apply; T=100)
for t in 1:T
𝛉 = rand.(Beta.(𝛂, 𝛃))
x = argmax(𝛉)
r = apply(x)
𝛂[x], 𝛃[x] = (𝛂[x] + r, 𝛃[x] + 1 - r)
end
return Beta.(𝛂, 𝛃)
end
```
Due to the lack of types, I have no immediate idea how to use this function. You have to read the code to “guess” that `alpha` and `beta` are probably 1D arrays (?) and `apply` is a callback. What argument does `apply` accept? Well, it’s the return value of `argmax`, so it’s probably an integer. It took me a few minutes looking through the code to figure it out.
What is nice about Fortran is that you have to declare the interface for `apply`, so you would know right away what it accepts. The price that you pay is that it is more verbose.
---
### Author:  [@jacobwilliams](https://fortran-lang.discourse.group/u/jacobwilliams)
#### Post date: [July 9, 2021, 10:24pm UTC](https://fortran-lang.discourse.group/t/julia-fast-as-fortran-beautiful-as-python/1405/137 "2021-07-09T22:24:35Z")
It is also far from obvious what the “.(“ is all about or why there is a “;” after apply. They definitely made some odd syntax choices.
What we need is a more interactive Fortran with some some kind of generics/type inference/multiple dispatch that isn’t incredibly verbose. Then I don’t see the need for Julia anymore. 😀
---
### Author:  [@oscardssmith](https://fortran-lang.discourse.group/u/oscardssmith)
#### Post date: [July 9, 2021, 10:33pm UTC](https://fortran-lang.discourse.group/t/julia-fast-as-fortran-beautiful-as-python/1405/138 "2021-07-09T22:33:54Z")
The use of `.` for broadcasting is far from unique to Julia (it’s also used by numpy and matlab), and IMO it does by far the best job of it. Having an easy way of vectorizing code is way better than having to define scalar and vector (and matrix and tensor) versions of everything.
---
### Author:  [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)
#### Post date: [July 9, 2021, 10:36pm UTC](https://fortran-lang.discourse.group/t/julia-fast-as-fortran-beautiful-as-python/1405/139 "2021-07-09T22:36:58Z")
> [@oscardssmith](#):
>
> Having an easy way of vectorizing code is way better than having to define scalar and vector (and matrix and tensor) versions of everything.
Fortran has intrinsic and user-defined ELEMENTAL functions for that.
---
### Author:  [@oscardssmith](https://fortran-lang.discourse.group/u/oscardssmith)
#### Post date: [July 9, 2021, 10:40pm UTC](https://fortran-lang.discourse.group/t/julia-fast-as-fortran-beautiful-as-python/1405/140 "2021-07-09T22:40:04Z")
From what I’ve seen there’s a key difference which is that you have to decide whether a function is elemental or not when you write it. This means that for example, you can’t broadcast matrix multiplication without re-writing BLAS. (also, you need different names for the elemental and non-elemental versions)
---
### Author:  [@shahmoradi](https://fortran-lang.discourse.group/u/shahmoradi)
#### Post date: [July 10, 2021, 3:02am UTC](https://fortran-lang.discourse.group/t/julia-fast-as-fortran-beautiful-as-python/1405/141 "2021-07-10T03:02:31Z")
So, does Julia vectorize a compiled external library at the request of the user, or does it recompile the library every time to achieve the goal?
---
### Author:  [@oscardssmith](https://fortran-lang.discourse.group/u/oscardssmith)
#### Post date: [July 10, 2021, 3:35am UTC](https://fortran-lang.discourse.group/t/julia-fast-as-fortran-beautiful-as-python/1405/142 "2021-07-10T03:35:14Z")
It’s somewhere in-between
```auto
using SpecialFunctions # for error function, erf
v = rand(10)
println(erf.(log.(v)))
```
will be lowered to `Base.materialize(Base.broadcasted(erf, Base.broadcasted(log, v)))` which basically turns into
```auto
y = similar(v)
for i in 1:10
y[i] = erf(log(v[i]))
end
```
If `log` and `erf` inline (in the particular they don’t since they are too expensive), Julia might make further optimizations like unrolling the loop, vectorizing the functions, or removing bounds checks.
---
### Author:  [@FortranFan](https://fortran-lang.discourse.group/u/FortranFan)
#### Post date: [July 10, 2021, 1:21pm UTC](https://fortran-lang.discourse.group/t/julia-fast-as-fortran-beautiful-as-python/1405/143 "2021-07-10T13:21:32Z")
> [@shahmoradi](#):
>
> So, does Julia vectorize a compiled external library at the request of the user, or does it recompile the library every time to achieve the goal?
Please see this example where the key computation, the sine of a radian, was moved to an external library, `cordic_sine`: [Simple summation 8x slower than in Julia - #44 by FortranFan](https://fortran-lang.discourse.group/t/simple-summation-8x-slower-than-in-julia/1171/44).
Once that was done, note in all the subsequent scenarios, the ratio of compute times with Fortran vs Julia were around unity, meaning no real difference i.e., the specific advantages with an enthusiast-driven `macro` implementation toward **vectorization** and specific “fast” trig functions leading to the original 8x difference became inapplicable.
---
### Author:  [@chriselrod](https://fortran-lang.discourse.group/u/chriselrod)
#### Post date: [July 10, 2021, 11:30pm UTC](https://fortran-lang.discourse.group/t/julia-fast-as-fortran-beautiful-as-python/1405/144 "2021-07-10T23:30:21Z")
> [@Beliavsky](#):
>
> Fortran has intrinsic and user-defined ELEMENTAL functions for that.
Sometimes the same function or operator has a different meaning when acting on matrices than elementwise.
```julia
julia> A = rand(2,2)
2×2 Matrix{Float64}:
0.198773 0.747052
0.0135511 0.777073
julia> exp.(A)
2×2 Matrix{Float64}:
1.21991 2.11077
1.01364 2.1751
julia> exp(A)
2×2 Matrix{Float64}:
1.22747 1.23599
0.0224203 2.18427
julia> A * A
2×2 Matrix{Float64}:
0.0496342 0.729008
0.0132238 0.613966
julia> A .* A
2×2 Matrix{Float64}:
0.0395108 0.558087
0.000183634 0.603843
julia> A ^ 3
2×2 Matrix{Float64}:
0.0197449 0.603572
0.0109485 0.486976
julia> A .^ 3
2×2 Matrix{Float64}:
0.00785369 0.41692
2.48845e-6 0.46923
```
---
### Author:  [@chriselrod](https://fortran-lang.discourse.group/u/chriselrod)
#### Post date: [July 11, 2021, 3:28am UTC](https://fortran-lang.discourse.group/t/julia-fast-as-fortran-beautiful-as-python/1405/146 "2021-07-11T03:28:43Z")
I’m not sure what you mean. `exp.(A)` and `exp(A)` do something different (elementwise vs [matrix exponential](https://en.wikipedia.org/wiki/Matrix_exponential)). The dot is used to syntactically indicate whether the function is applied to the entire argument, or applied elementwise.
Just like having type information listed makes it easier to read code, this is also informative, i.e. we know that `foo.(A)` is elementwise and `foo(A)` probably isn’t (unless someone defined `foo(A)` to explicitly loop over `A`).
---
### Author:  [@chriselrod](https://fortran-lang.discourse.group/u/chriselrod)
#### Post date: [July 11, 2021, 4:44am UTC](https://fortran-lang.discourse.group/t/julia-fast-as-fortran-beautiful-as-python/1405/148 "2021-07-11T04:44:59Z")
I’m lost here by the combative tone your comments are striking.
In Julia, we can define `A::AbstractMatrix{Float64} * B::AbstractMatrix{Float64}` to be `dgemm` and `A::Float64 * B::Float64` to be scalar multiplication, and then automatically `A::AbstractMatrix{Float64} .* B::AbstractMatrix{Float64}` is elementwise multiplication.
Thus we do need to define two different functions, but the “ELEMENTAL” behavior comes for free by just defining the elemental (scalar) version of the function.
Then syntactically, someone reading our code knows whether our function is applying elementwise, or applying to the entire argument.
Similarly, one knows syntactically if it is fusing with other `.`, so one can look at code and tell how many array temporaries are being created.
Then, if the same operator has a different meaning when applied to a certain container type – like a matrix – we are still free to add a definition for that container type, without needing to use other operators/function names.
My argument would only be specious if we had to explicitly define `exp.()` and `.*`.
The `.` is syntax. All we have to define is `exp(::Float64)` to get both `exp(::Float64)` (no dot) and `exp.(::AbstractArray{Float64})` (elementwise on a matrix) or `exp.((A::AbstractArray .- B::AbstractArray) .^ 2)` (elementwise exponential of squared differences, without allocating any temporaries).
---
### Author:  [@implicitall](https://fortran-lang.discourse.group/u/implicitall)
#### Post date: [July 11, 2021, 10:44am UTC](https://fortran-lang.discourse.group/t/julia-fast-as-fortran-beautiful-as-python/1405/150 "2021-07-11T10:44:09Z")
I wanted to see if I could implement what you are saying just to see how far Fortran is away from that. Take this with a grain of salt as a goofy example. You would typically use intrinsics and elemental functions as noted earlier.
```auto
program test
use, intrinsic :: iso_fortran_env
interface operator(.elem.)
procedure elem
end interface
interface exp
procedure exp_mat
procedure exp_elem
end interface
interface operator(.star.)
procedure star_mat
procedure star_elem
end interface
type elem_ptr
real(real64), pointer :: mat(:, :)
end type
real(real64) :: A(2, 2)
A = reshape([0.198773, 0.747052, &
0.0135511, 0.777073], shape=[2, 2], order=[2, 1])
call print_mat(A)
call print_mat(exp(.elem.A))
call print_mat(exp(A))
call print_mat(A.star.A)
call print_mat(A.star..elem.A)
contains
function elem(x) result(y)
real(real64), target, intent(in) :: x(:, :)
type(elem_ptr) :: y
y%mat => x
end function
function exp_mat(x) result(y)
real(real64), intent(in) :: x(:, :)
real(real64) :: y(size(x, 1), size(x, 2))
real(real64) :: x_k(size(x, 1), size(x, 2))
integer :: i, j, k, fact
y = 0.
do concurrent(i=1:size(x, 1), j=1:size(x, 2))
if(i == j) y(i, j) = 1.
enddo
fact = 1
x_k = x
do k = 1, 8
fact = fact*k
y = y + x_k/fact
x_k = matmul(x_k, x)
enddo
end function
function exp_elem(x) result(y)
type(elem_ptr), intent(in) :: x
real(real64) :: y(size(x%mat, 1), size(x%mat, 2))
integer :: i, j
do concurrent(i=1:size(x%mat, 1), j=1:size(x%mat, 2))
y(i, j) = exp(x%mat(i, j))
enddo
end function
function star_mat(a, b) result(c)
real(real64), intent(in) :: a(:, :), b(:, :)
real(real64) :: c(size(a, 1), size(b, 2))
c = matmul(a, b)
end function
function star_elem(a, b) result(c)
real(real64), intent(in) :: a(:, :)
type(elem_ptr), intent(in) :: b
real(real64) :: c(size(a, 1), size(a, 2))
c = a*b%mat
end function
subroutine print_mat(mat)
real(real64), intent(in) :: mat(:, :)
integer(int32) :: i, j, m, n
write(*, '(i0,"x",i0," Matrix{Float64}")') size(mat, 1), size(mat, 2)
do i = 1, size(mat, 1)
write(*, '(22f15.7)', advance='no')(mat(i, j), j=1, size(mat, 2))
write(*, *)
enddo
end subroutine
end program
```
---
### Author:  [@lmiq](https://fortran-lang.discourse.group/u/lmiq)
#### Post date: [July 11, 2021, 2:52pm UTC](https://fortran-lang.discourse.group/t/julia-fast-as-fortran-beautiful-as-python/1405/151 "2021-07-11T14:52:02Z")
> [@kargl](#):
>
> How is this any different than what Julia does (other than syntax)?
Note that, in `exp` case (and in any other), one could do:
```auto
a = 1.0
A = ones(2,2)
exp(a) # scalar
exp(A) # exp of matrix
v = [a, a, a]
exp.(a) # scalar, element wise
exp.(A) # scalar, element wise
V = [A, A, A] # vector of matices
exp.(V) # Matrix exp on each element of V
```
You can do this with essentially any function in Julia, even if it is function defined in a library. In this case what is needed is one `exp` method defined for scalars and other for matrices. (Because of generics, for many other functions only one method is needed at all)
[Previous page](https://fortran-lang.discourse.group/t/julia-fast-as-fortran-beautiful-as-python/1405.md?page=6)
[Next page](https://fortran-lang.discourse.group/t/julia-fast-as-fortran-beautiful-as-python/1405.md?page=8)