The Accelerating Adoption of Julia

Producer-consumer queues say otherwise. As does Async. There's also message-passing. There's also thread-pools.

Task-based programming doesn't wait for the tasks to finish. So without the "join", you don't have an easy synchronization point to rely upon, and have to use explicit locks for synchronization.

------

Producer - consumer can't work with fork/join, because there's no "join" at all! The producer keeps producing, and the consumers just keep consuming.

Async is... a tangled mess. I'm pretty sure its just the modern "unstructured goto" style of spaghetti that has been rediscovered. Yeah, its super-efficient, but no, its not really easy to program.

Task-based parallelism is probably my next favorite style. Fork-join is super easy, but the join causes a significant waiting period and underutilization of the processor. Task-based provides more flexibility (at the cost of a little bit more complexity).

------

There's a ton of other models of parallelism. But fork-join seems to be the methodology that more and more work is being focused on.

If you need to synchronize, wait for the next fork/join cycle instead of doing explicit mutexes. The "join" provides your synchronization point.

If you have a complicated set of highly-synchronous calculations, then you don't fork at all. You simply use the "Single Thread" to perform all those calculations (and therefore negate the need of cross-thread synchronization).

Fork-join has low-utilization, but its very easy to program. In some cases, fork-join remains efficient (ex: spawn a thread per pixel on the screen), because all pixels do NOT need to synchronize with each other (or call mutexes).

-----

If a sync point within the children are needed, you usually make due with a barrier instead of a mutex.

> What does 'Single Instruction Multiple Data' have to do with threads?

SIMD processors, such as GPUs, emulate a thread in their SIMD units. Its called a "CUDA Thread". Its not a true thread in the sense of CPU-threads, but the performance by-and-large scales as if it were real threads. (With exception of the "thread-divergence problem")

Ultimately, the fork-join model translates trivially to SIMD. Any practitioner of CUDA, OpenCL, ROCm, or ISPC can prove that to you easily.

I don't know either.

However, GPUs basically do SMT behind the scenes. Because GPUs do not reorder instructions like a CPU they stall on memory accesses but instead of waiting they simply switch to another thread. This means you need a lot of threads for good performance.

I didn't take it as such; there are legitimate shortcomings to any tool, I just wanted to provide pointers to other readers that the devs are aware of it, and that there is ongoing development to address it. :)

> Come on, jeez.. Julia’s compiler is a Lisp-based LLVM driver: of course it can do these things.

As someone who, before Julia, was firmly entrenched in C/C++/Python land, I suppose I am discovering many of these "obvious" things for the first time. :)

Then suppose I rewrite the LU decomposition algorithm, using python. Now I want to accelerate the search by running the search on GPUs. I have to re-rewrite the code from scratch. Each GF256 encoding has to have rejiggered operators, and so I need to rewrite custom GPU kernels, then figure out how to resequence the operations (* looks different for each GF256 encoding), etc.

This is all SUPER easy in julia.

Those statements are just confusing to me. In most GPU-code, you have crazy amounts of parallelism and therefore don't really care what order those statements execute in.

As such, if you can allocate memory, and map those if/else statements into a collection of queues and/or stacks (depending on if you want a breadth-first search, or depth-first search pattern).

In effect: you use memory allocation to solve the complex control flow issue. Maybe its more obvious with code:

Instead of doing:

    if(baz()){
      foo();
    } else {
      bar(); // Thread divergence!!
    }

    if(baz()){
      pushIntoFooQueue();
    } else {
      pushIntoBarQueue(); 
      // Thread divergence, but not much of a penalty
    }

    while(fooIsNotEmpty()){ // No thread divergence at all
        task = SIMDPopFoo();
        task.execute();
    }

    while(barIsNotEmpty()){
        task = SIMDPopBar();
        task.execute();
    }

Bonus points: Not only is thread-divergence negated, but you also achieve effective load-balancing across your workgroup. If Thread#0 spawns 20 items for FooQueue, after pushing/popping from the queue, those 20 Foo-tasks will be assigned to 20 different threads.

-----

From there on out, you're just pushing / popping different parts of your code to various queues and/or stacks. But this is only really a valid solution if you have a decent memory allocator that knows where and how to clean up these queues / stacks (especially if you have nodes starting to point to each other for dependency management)

I haven't really solved this problem "in general", but is clear that the queues should be sorted into topological order, and that any tasks that depend onto each other need to be run in different iterations. It really depends on how much you're willing to spend on organizing this execution information.

In any case: the memory allocation issue is one-and-the-same with the thread-divergence / complex instruction flow issue to me. You need to create a data-structure that organizes the instruction flow, and any complex data-structure will need memory management as soon as you start linking things together. The above uses a queue or stack, but things can get more complicated.

--------

Not that Julia, Python, ISPC or anything really... solves this problem. But memory management is very useful for this "style" (I'm mostly doing ref-counted C++ with a custom allocator myself. But such trees or graphs of links can grow into the CPU-side and end up using the default memory allocator)

It's a different language for a different computer. But the fundamentals are still there, laying the groundwork for the future.

Sure, ispc can win on latency. But the raw compute girth of GPU SIMD cannot be underestimated.

In any case, Julia has proven itself capable of compiling down to a SIMD instruction set and achieving nearly the full performance of those GPUs. Even if it's not a computer you prefer, the programming model and technology demonstrated here is clear.

https://github.com/fonsp/Pluto.jl is a super-nice way of just starting, without too much setup cost. I think it will soon be self-contained as to what packages & versions are needed, too.

A bigger barrier here is that we had to know how to apply propagation of errors, so although we could use a tool like this to do the calculation we still had to generate a set of equations for the lower bound, best estimate, and upper bound. There's a seductive argument to be made for just plugging those equations into Excel instead of using them to validate that the Julia functions are working correctly.

http://mesh-spreadsheet.com

If the language is single dispatch though, animal1.f(animal2) will accurately call the Cat or Dog method for the first argument (before the dot), but the second argument will work like above, it will call the Animal version and not the Cat or Dog ones. If you want a method to call Cat or Dog method based on runtime information you'd have to redispatch on runtime within f(Animal), which is basically a pattern known as Visitor Pattern [1].

In Julia (and multiple dispatch languages) defining f(Cat, Dog) and f(Animal, Animal) means that even when I call f(animal1, animal2) the compiler will convert to the runtime values and if the first animal is Cat and the second is Dog it will call the first method (as a form of specialization) and otherwise it will go for f(Animal, Animal) (which is the more general method). And the Julia compiler is really good at inferring runtime types at compile time, so even though it's runtime values it's still static dispatch (that's why it's also fast).

[1] https://en.wikipedia.org/wiki/Visitor_pattern

You can find examples where operator overloading and multiple dispatch actually behaves differently and gives different answers, but the comparison works as a first order approximation.

Remember, I am presenting what I think are counter-examples to "is strictly superior to Python in everyway" (excluding "existing projects"). I still think that's true even if Python weren't popular at all.

As an example of the influence, the ":" at the end of lines which start a block in Python is not required, in the language sense. However, ABC testing showed it made the language more useful.

As an example of improving Python for learners, one of the reasons for Python's change from 1/2 == 0 -> 1/2 == 0.5 was feedback from the Alice developers, where they found that students did not understand van Rossum's decision to follow C's semantics.

Python started being used in schools when it wasn't popular. One of the presenters at the Python conference in 1998 or so was a high school AP comp sci teacher. At that time C++ and Java were the AP languages. He reported better success teaching students Python first and then teaching C++ or Java, than spending all of the time in one of the latter languages.

I would do this exactly too, if I had to do number crunching.

These days, Nim is a better Cython but with less dynamic temptations and more powerful metaprogramming (and, yes, a much smaller ecosystem..maybe not that much smaller than Python in the late 90s, though). { Not that this is all Nim is...It's actually a really good everything-language that's tricky to summarize in just a few words. }

1. It's very easy to inspect generated code of your kernels where you really need SIMD. For instance:

    julia> function my_kernel(xs)
             total = zero(eltype(xs))
             for x in xs
               @fastmath total += x
             end
             total
           end
    julia> @code_native debuginfo=:none my_kernel(rand(10))
    ...
    L96:
     vaddpd 8(%rcx,%rax,8), %ymm0, %ymm0
     vaddpd 40(%rcx,%rax,8), %ymm1, %ymm1
     vaddpd 72(%rcx,%rax,8), %ymm2, %ymm2
     vaddpd 104(%rcx,%rax,8), %ymm3, %ymm3
     addq $16, %rax
     cmpq %rax, %rdi
     jne L96
    ...
    

3. You can make LLVM explain why it did not vectorize certain loops just like in clang.

For performance you virtually always need to write 'vectorized' code, which is often difficult and memory heavy.

The daily churn of handling input parsing, along with contorting code into vectorized form takes up an extraordinary amount of time when writing Matlab code.

Matlab is certainly not 'good enough' for me! It's just what I'm stuck with at work.

That said, in terms of composability you can jit over the closure to achieve a lot of what you might want, e.g.

    def make_loop(f):
        @jit
        def fn(x):
            for i in range(x.shape[0]):
                 x[i] = f(x[i])
        return fn

For introspection and simplicity, I think, in high performance, you simply have to choose two of fast, simple and generic. Julia clearly chooses fast and generic.