Julia Community 🟣

Discussion on: Replacing legacy code with Julia

View post
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).