Julia Community 🟣

Discussion on: Replacing legacy code with Julia

Collapse
 
omendezmorales profile image
Orlando Méndez

NIce article Petr... and the question I guess every Julia enthusiast will ask next: is the Julia Code on par or faster with the legacy Fortran one?

Thanks,

Collapse
 
lmiq profile image
Leandro Martínez • Edited

I can provide a half-way comparison, for the case of the dense matrix. The corresponding Fortran code is (in more modern notation, which is quite similar to Julia's):

subroutine dense(M,F)
    implicit none
    integer, parameter :: dp = kind(1.d0)
    integer :: M, j, i, r
    real(dp) :: s, t
    real(dp), dimension(M,M) :: F
    do j = 2, M
        do i = 1, j-1
            F(i, j) = F(i, j) - dot_product(F(1:i-1, i), F(1:i-1, j))
        end do
        s = F(j, j)
        do r = 1, j-1
            t = F(r, j) / F(r, r)
            s = s - t * F(r, j)
            F(r, j) = t
        end do
        F(j, j) = s
    end do
end subroutine dense
Enter fullscreen mode Exit fullscreen mode

You can compile it with:

gfortran -O3 -shared -o factor_jl.so factor_jl.f90
Enter fullscreen mode Exit fullscreen mode

Then in Julia, you can do:

julia> using LinearAlgebra: dot

julia> @views function dense_ldlt!(F)
           M = size(F, 1)
           for j in 2:M
               for i in 1:j-1
                   F[i, j] -= dot(F[1:i-1, i], F[1:i-1, j])
               end
               s = F[j, j]
               for r in 1:j-1
                   t = F[r, j] / F[r, r]
                   s -= t * F[r, j]
                   F[r, j] = t
               end
               F[j, j] = s
           end
           F
       end
dense_ldlt! (generic function with 1 method)
Enter fullscreen mode Exit fullscreen mode

And the benchmark:

julia> a = rand(100,100);

julia> @btime dense_ldlt!(c) setup=(c=copy($a)) evals=1;
  108.731 μs (0 allocations: 0 bytes)

julia> f(a) = ccall((:dense_,"./factor_jl.so"),Nothing,(Ref{Int},Ref{Float64},),size(a,1),a)
f (generic function with 1 method)

julia> @btime f(c) setup=(c=copy($a)) evals=1
  100.116 μs (0 allocations: 0 bytes)

Enter fullscreen mode Exit fullscreen mode

The second one is the call to the fortran routine, which turns to be slightly faster. Note that I had to add @views to the Julia function, otherwise the slices used in the dot operation, by default, allocate new arrays (and the function becomes much slower).

Collapse
 
petrkryslucsd profile image
Petr Krysl

Did you use julia -O3 as well?

Thread Thread
 
lmiq profile image
Leandro Martínez

No, but I don't see any difference in this case. (nor adding @inbounds makes any measurable difference for me).