Classing for example ISPC as low level invalidates this comment for me; ISPC among others provides a great interface to SIMD model for CPUs that simply isn’t available in Julia.
That being said, packages such as SIMD.jl [0], and LoopVectorization.jl [1] are making fantastic progress, to the point that LoopVectorization forms the basis of a legitimate BLAS contender, in pure Julia [2]. It's not totally there yet, but it's close enough that real work is being done in LV at OpenBLAS-like speeds.
As an aside, I find it incredible that these kinds of extensions can be built in packages thanks to the fact that Julia's compiler is extensible enough to allow for direct manipulation of the LLVM intrinsics being emitted by user code.
[0] https://github.com/eschnett/SIMD.jl [1] https://github.com/chriselrod/LoopVectorization.jl [2] https://github.com/MasonProtter/Gaius.jl
> find it incredible that these kinds of extensions can be built in packages thanks to the fact that Julia's compiler is extensible
Come on, jeez.. Julia’s compiler is a Lisp-based LLVM driver: of course it can do these things.
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. :)
sorry for the flippant remark then! it's great to be in discover mode, enjoy ;)
Indeed: ISPC provides constructs to provide structure-of-arrays, and other low-level memory layout details. This is a good thing: these details have significant implications on the speed of your program.
Nonetheless, any language which wrangles with manual details of memory layout, or new/delete based memory allocation, is inevitably going to be classified as low level in my books.
> ISPC among others provides a great interface to SIMD model for CPUs that simply isn’t available in Julia.
If Julia can compile into GPU-assembly (which is innately SIMD), I'm sure an AVX-based compile could work eventually.
They may have to target AVX512 (since most GPU-assembly requires per-thread exec-masks), but the general concept is being proven as Julia can now compile down into PTX or AMDGPU assembly.
Julia's compile down to GPU-assembly / SIMD code is not supported in the "general case", only in select circumstances. But that's still an incredible boon for a high-level language.
Eventually is nice but it’s also easy to imagine that it never gets done. ISPC works today.
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.
My comment was more geared to the CPU side of things since I write code whose users don’t usually have GPUs available.
There are also problem sizes which don’t fit the GPU’s more stark memory system separation: tens of CPU cores with tens of MBs of cache can move past GPUs in terms of memory bandwidth for such cases so banking on CUDA just doesn’t work for everyone.
I guess we aren’t solving the same problems: memory allocation is trivial is my domain, mapping complex nested control flow in SIMD is the hard part.
> Julia can now compile down into PTX or AMDGPU assembly.
Sure but Julia-as-syntax is nothing special now, Numba does this for a Python as well.
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.
https://www.youtube.com/watch?v=Gh578e98qAk
> Suppose I want to write a library to find compression polynomials for reed-solomon encoding system (see mary wootters' talks on youtube) for my storage product. I need a LU decomposition algorithm that operates on GF256
Surely you use isa-l[1].
[1] https://github.com/intel/isa-l
> Now I want to accelerate the search by running the search on GPUs.
GPUs are float oriented so I don't think you'll get the performance you hope for out of 8 bit integer operations. If you have interesting results to share I'd like to read them.
You've never seen GPU Hashcat, or GPU Bitcoin / Ethereum mining?
Vega now has a huge focus on INT8 operations. NVidia can perform int operations in parallel with float operations (superscalar GPU cores)
I've heard of it, but afaik it hasn't been profitable to mine using a GPU for a long time due to the competition for hashes, power consumption, and rate they mine at. This is in contrast to ASICs which can mine even faster for less capex and opex.
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)
This is an interesting remark I have to sleep on. I usually don't see this as possible, since I have a bunch of arrays which are allocated once, then a bunch of nested but static control flow. The only way I've found to go beyond single thread performance is with "whole program" SIMD which only worked in the ISPC/OpenCL/CUDA programming models.
But its a different paradigm to do things: something to try if your standard if/else stuff isn't working out.
EDIT: If you can make due with just stack-allocation (or queue-allocation), if you have singular-sized tasks with predictable sizes, if all the information fits inside of __shared__ memory, and if thread-divergence is normally a problem... this methodology will probably work.
Oh, and for a hint:
myIdx = prefixSum(__activemask());
myTask = stack[stackTail + myIdx];
if(myIdx == 0){
stackTail -= __popc(__activemask()); // Lulz, I made a bug. Whatever, I think you understand...
}
__syncThreads(); // Barrier is important