Nice blogpost comparing Julia w/ 6 other languages to solve a Logistic DE

New to Julia
21

I was wondering about this too. The post says Fortran is used only in legacy projects but I seem to know a lot of people who use it on a daily basis. Could you elaborate on what else you think was wrong?

1 Like
22

I wonder if using packages is allowed, @turbo brings another 4X speedup over the optimized version on my computer. An extra 4X is also possible with @tturbo. IMO, any fair benchmark should allow packages and tools of the language especially if written in the same language and not in another low-level alternative.

using LoopVectorization

function julia_logistic_turbo!(N,dt,u) 
    # Parameters
    u0 = 1e-5
    # Right hand side function
    f(U)= U * (1.0-U)
    # Discretization
    u[1] = u0
    @turbo for k = 1:N-1
        u[k+1] = u[k] + dt*f(u[k])
    end
end

Timing:

@btime julia_logistic_turbo!(N, 1.0f0, u) setup=(N=10000; u=zeros(Float32,N))
# The last optimized version on the website
@btime julia_logistic!(N, 1.0f0, u) setup=(N=10000; u=zeros(Float32,N))

  9.600 μs (0 allocations: 0 bytes)
  35.000 μs (0 allocations: 0 bytes)
1 Like
23

I think that the test problem is too simple to be a meaningful performance comparison and the comments about both C++ and Fortran are just silly. Actually, the Python example isn’t very meaningful either - any real computation people are doing in Python will either be vectorized operations calling out to numpy or a JIT-compiled kernel using Numba or an extension module written in C.

Re. C++: His mocking tone when people pointed out to him that the performance was probably because he’s using it with a notebook kernel is off putting. I’m not sure what that infrastructure is like, but I bet that in the original case he was timing JIT compilation of the C++ code rather than runtime - the notebook kernel for C++ I’ve used was LLVM-based JIT of each cell.

Re. Fortran: I’m not a very competent Fortran developer but I’ve worked a bit in a modern Fortran code base, and it can be a very nice language for scientific computing with all the modern features that have been added. Native support for multidimensional arrays makes it very nice to write array-heavy scientific code in, much like Julia.

Of course these statically typed languages are more work to develop a code in than Julia, but they have their place for some things for now. For me one of the best parts of Julia is metaprogramming, but this can make it impossible to create a static executable. My experience has also been that when performance really matters, I have to spend significant time profiling and optimizing Julia code to get the same performance as what I’d get the first time writing it in C++.

2 Likes
24

There is an option to compile with -fcheck=bounds.

My own C++ code uses llvm::SmallVectors, which while poorly documented basically acts like a drop in replacement for std::vector. It bounds checks by default in debug builds:
https://llvm.org/doxygen/SmallVector_8h_source.html#l00273
Compiling with -DNDEBUG should disable the asserts.

llvm is a bit heavy of a dependency, so it’s probably not worth depending on it just for their SmallVector data type. Perhaps boost’s version is equally good.
C++ is definitely also hurt here by their bad package management situation.

FWIW, I like loops and indexing (instead of patterns like operating on iterators and using higher order functions, which seems to be the preferred C++ style), and have gotten tons of bounds errors/asserts.
Having bounds checks is a huge time saver when debugging.

Unfortunately, even #include <debug/vector> does not bounds check by default. =/

2 Likes
25

Agree with preferring loops and indexing but higher-order functions working on iterators is a very powerful pattern for writing libraries.

It’s very simple to make a debug vector type if you want one:

template <class T>
class my_vector : public std::vector<T>
{
public:
    using std::vector<T>::vector; // Has all of vector's constructors
    
    const T& operator[](size_t i) const
    {
#ifndef NDEBUG
        return std::vector<T>::at(i); // vector::at is the bounds checked version
#else
        return std::vector<T>::operator[](i);
#endif
    }
// Add the non-const overload of operator[] here
};

Might contain syntax errors. Fixed in godbolt

2 Likes