# Bend: a new GPU-native language
**URL:** https://discourse.julialang.org/t/bend-a-new-gpu-native-language/114440
**Category:** Offtopic
**Tags:** gpu
**Created:** [May 18, 2024, 10:28pm UTC](https://discourse.julialang.org/t/bend-a-new-gpu-native-language/114440 "2024-05-18T22:28:05Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@MilesCranmer](https://discourse.julialang.org/u/MilesCranmer)
#### Post date: [May 18, 2024, 10:28pm UTC](https://discourse.julialang.org/t/bend-a-new-gpu-native-language/114440/1 "2024-05-18T22:28:05Z")
I thought this new language looked cool: [GitHub - HigherOrderCO/Bend: A massively parallel, high-level programming language](https://github.com/HigherOrderCO/Bend). It’s still in the research stages but presents some interesting ideas.
It compiles a Python-like language directly into GPU kernels – more than just array broadcasting – by exploiting parallelism wherever it can:
> [@](#):
>
> ## Parallel Programming in Bend
>
> To write parallel programs in Bend, all you have to do is… **nothing**. Other than not making it _inherently sequential_! For example, the expression:
>
> ```python
> (((1 + 2) + 3) + 4)
>
> ```
>
> Can **not** run in parallel, because `+4` depends on `+3` which depends on `(1+2)`. But the following expression:
>
> ```python
> ((1 + 2) + (3 + 4))
>
> ```
>
> Can run in parallel, because `(1+2)` and `(3+4)` are independent; and it _will_, per Bend’s fundamental pledge:
>
> > Everything that **can** run in parallel, **will** run in parallel.
>
> For a more complete example, consider:
>
> ```python
> # Sorting Network = just rotate trees!
> def sort(d, s, tree):
> switch d:
> case 0:
> return tree
> case _:
> (x,y) = tree
> lft = sort(d-1, 0, x)
> rgt = sort(d-1, 1, y)
> return rots(d, s, lft, rgt)
>
> # Rotates sub-trees (Blue/Green Box)
> def rots(d, s, tree):
> switch d:
> case 0:
> return tree
> case _:
> (x,y) = tree
> return down(d, s, warp(d-1, s, x, y))
>
> (...)
>
> ```
>
> This [file](https://gist.github.com/VictorTaelin/face210ca4bc30d96b2d5980278d3921) implements a [bitonic sorter](https://en.wikipedia.org/wiki/Bitonic_sorter) with _immutable tree rotations_. It is not the kind of algorithm you’d expect to run fast on GPUs. Yet, since it uses a divide-and-conquer approach, which is _inherently parallel_, Bend will run it multi-threaded. Some benchmarks:
>
> - CPU, Apple M3 Max, 1 thread: **12.15 seconds**
>
> - CPU, Apple M3 Max, 16 threads: **0.96 seconds**
>
> - GPU, NVIDIA RTX 4090, 16k threads: **0.21 seconds**
>
> That’s a **57x speedup** by doing nothing. No thread spawning, no explicit management of locks, mutexes. We just asked Bend to run our program on RTX, and it did. Simple as that.
I wonder if Julia could eventually try to do something like this.
They have a paper on it here: [https://raw.githubusercontent.com/HigherOrderCO/HVM/main/paper/PAPER.pdf](https://raw.githubusercontent.com/HigherOrderCO/HVM/main/paper/PAPER.pdf)
The underlying compiler is here: [GitHub - HigherOrderCO/HVM: A massively parallel, optimal functional runtime in Rust](https://github.com/HigherOrderCO/HVM)
> [@](#):
>
> ## Language
>
> HVM is a low-level compile target for high-level languages. It provides a raw
> syntax for wiring interaction nets. For example:
>
> ```javascript
> @main = a
> & @sum ~ (28 (0 a))
>
> @sum = (?(((a a) @sum__C0) b) b)
>
> @sum__C0 = ({c a} ({$([*2] $([+1] d)) $([*2] $([+0] b))} f))
> &! @sum ~ (a (b $(:[+] $(e f))))
> &! @sum ~ (c (d e))
>
> ```
>
> The file above implements a recursive sum. If that looks unreadable to you - don’t worry, it isn’t meant to. [Bend](https://github.com/HigherOrderCO/Bend) is the human-readable language and should be used both by end users and by languages aiming to target the HVM. If you’re looking to learn more about the core
> syntax and tech, though, please check the paper.
So compiling to this target “HVM”, a language can exploit their GPU-ification.
* * *
Edit: changed title so its clearer that HVM is GPU-native, not just CPU threads
---
### Author:  [@jling](https://discourse.julialang.org/u/jling)
#### Post date: [May 19, 2024, 12:20am UTC](https://discourse.julialang.org/t/bend-a-new-gpu-native-language/114440/2 "2024-05-19T00:20:57Z")
> [@MilesCranmer](#):
>
> That’s a **57x speedup** by doing nothing.
I remember someone pointed out that it’s partly due to their baseline being very slow.
---
### Author:  [@MilesCranmer](https://discourse.julialang.org/u/MilesCranmer)
#### Post date: [May 19, 2024, 12:24am UTC](https://discourse.julialang.org/t/bend-a-new-gpu-native-language/114440/3 "2024-05-19T00:24:55Z")
Do you have a link? I read through the HackerNews thread where someone commented on this number being low but it seems like they didn’t realize that 1 GPU core \<\< 1 CPU core
---
### Author:  [@Tarny\_GG\_Channie](https://discourse.julialang.org/u/Tarny_GG_Channie)
#### Post date: [May 19, 2024, 12:51am UTC](https://discourse.julialang.org/t/bend-a-new-gpu-native-language/114440/4 "2024-05-19T00:51:45Z")
I’m doubtful. I don’t think a language that can automatically allocates a thread every time something gets executed would be optimal. Spawning a thread has an overhead. Locking and unlocking has an overhead. If your language is slow to start with on an each thread, then you can naively speed these up with parallelism. The hard part of many parallel algorithms is not sprinkling in a bunch of “spawn”, or “lock”, it’s designing the algorithm to be parallelizable and perhaps lock-free and so on.
---
### Author:  [@jling](https://discourse.julialang.org/u/jling)
#### Post date: [May 19, 2024, 1:13am UTC](https://discourse.julialang.org/t/bend-a-new-gpu-native-language/114440/5 "2024-05-19T01:13:02Z")
The HY comment points out that python single thread is much faster than Bend single thread, and that pypy is much faster than Bend GPU.
I think this doesn’t invalid the idea of something like Bend, but if the base language is fast (Julia) it’s very hard to automatically find better parallelism algorithms without being slower for some problems sizes.
For example, `all(predicate, array)` can be parallelized, but depending on problem size and probability of early termination, how would the language decide if it should run in parallel or not
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 19, 2024, 1:52am UTC](https://discourse.julialang.org/t/bend-a-new-gpu-native-language/114440/6 "2024-05-19T01:52:32Z")
We’ve had quite a few projects on automating different forms of parallelism. The issue is not whether you can parallelize things in a compiler, it’s whether you want to. Solving the scheduling problem for whether a given form of parallelism will actually make things faster or slower is hard on shared memory since threads have a non-trivial spin up cost, and so you need some pretty good cost models. It’s even harder for something that’s multiprocessed or GPU because you have to factor in memory transfer casts.
Whether parallelism is a good idea is not necessarily fully captured via the program either, you need to know the size of data in many cases in order to know if it makes sense. For example, for a large enough matrix then sending it to the GPU, computing there, and sending it back, will beat out even the best CPU. Whether you should even multithread at all is dependent on if it’s too small (you cannot multithread an 8x8 matmul and expect to win). So this means you need dynamic scheduling, and a dynamic scheduler itself has overhead. If you’d want to mask that overhead, then you need a compiler which can figure out what places you may want to have a more streamlined computation (like an inner loop) and pull it out of there, but not all inner loops because then you’d never multithread a linear algebra implementation.
One other factor to keep in mind is that the choice of scheduling is not something with a unique solution. Most frameworks that “auto” parallel in some form make some kind of choice as to how to perform the scheduling in some set manner. For example, `vmap` is an “auto parallelism” construct in many machine learning libraries which, while it’s good for linear algebra, it’s quite obviously a bad idea for kernels not dominated by linear algebra which is how we can show it’s (PyTorch and Jax’s implementation) 20x-100x slower than doing the right thing.
[https://www.sciencedirect.com/science/article/abs/pii/S0045782523007156](https://www.sciencedirect.com/science/article/abs/pii/S0045782523007156)
Not very surprising of course when you consider how vmap chooses to schedule, but of course that shows you that even a billion dollar company cannot beat a normal human taking into considering how a compute should be happening.
Finally, whenever I hear this kind of discussion, I always think back to the NetworkX vs Lightgraphs.jl early discussions. People lamented for awhile that it would be hard to catch up to NetworkX and all of its parallelism goodies, back before Julia had multithreading. But someone did an independent benchmark and…
> **[Benchmark of popular graph/network packages v2](https://www.timlrx.com/blog/benchmark-of-popular-graph-network-packages-v2)**
>
> A revised benchmark of graphs / network computation packages featuring an updated methodology and more comprehensive testing. Find out how Networkx, igraph, graph-tool, Networkit, SNAP and lightgraphs perform
Well that’s the v2 of it, but same results. Basically, NetworkX has a bunch of papers talking about how great its parallel scaling is, but it’s about 30x slower than a single core implementation for most size graphs you can fit onto a shared memory machine.
Parallelism and scaling is a good property to have, but you cannot act like a reinforcement learning algorithm and just optimize scaling at all costs. Remember, the goal is to make codes faster, not more parallel.
---
### Author:  [@Benny](https://discourse.julialang.org/u/Benny)
#### Post date: [May 19, 2024, 2:17am UTC](https://discourse.julialang.org/t/bend-a-new-gpu-native-language/114440/7 "2024-05-19T02:17:25Z")
I know this is kinda ironic in a Julia forum, but I wish people would stop making comparisons to Python when the langauge is not even similar enough to be a subset or a derivative. You don’t annotate types and the function keyword is `def`, but there’s no `for` or `class`, which may be reasonable high-level limitations to allow parallelism). From their incomplete paper, Python and Haskell are _potential_ targets for compilation to the intermediate HVM2, Bend is just the language they could demonstrate for now.
---
### Author:  [@MilesCranmer](https://discourse.julialang.org/u/MilesCranmer)
#### Post date: [May 19, 2024, 2:50am UTC](https://discourse.julialang.org/t/bend-a-new-gpu-native-language/114440/8 "2024-05-19T02:50:36Z")
I think some people in this thread may be missing the key innovation/coolness here and instead projecting it onto existing libraries.
This is not `vmap`. This is a language that, in its entirety, can generate efficient CUDA kernels, parallelizing any operations which _can_ be run in parallel.
The only thing I’ve ever seen close to this is JAX, which via XLA can fuse CUDA kernels, but in JAX you are severely limited in what you can do, as you are basically using Python for meta-programming C++.
But this is a full compiler that can generate massively parallel GPU kernels from high-level code. I haven’t seen anything like this before.
And don’t think of this as a competitor; the “bend” language is just an example implementation, but [HVM2](https://github.com/HigherOrderCO/HVM) is something that Julia could actually use (or maybe a macro).
* * *
Anyways also just to point out, it’s not optimized yet, it seems like a research language which just got off the ground. i.e., you don’t want to naively compare performance numbers with other languages yet. The scaling is what matters.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 19, 2024, 3:24am UTC](https://discourse.julialang.org/t/bend-a-new-gpu-native-language/114440/9 "2024-05-19T03:24:13Z")
That seems to not address anything that was mentioned above. vmap is one of many forms of parallelism. Others like Cilk are already built into Julia in some form. Remember, the precurser to Julia in some sense was StarP, a parallel MATLAB compiler
> **[star-p.pdf](https://math.mit.edu/~edelman/publications/star-p.pdf)**
>
> 1185.62 KB
Early projects such as the ParallelAccelerator subset (made by Intel):
> **[LIPIcs.ECOOP.2017.4.pdf](https://drops.dagstuhl.de/storage/00lipics/lipics-vol074-ecoop2017/LIPIcs.ECOOP.2017.4/LIPIcs.ECOOP.2017.4.pdf)**
>
> 825.97 KB
before it was abandoned due to Tapir.
Out-of-core scheduler projects still exist such as Dagger:
> **[GitHub - JuliaParallel/Dagger.jl: A framework for out-of-core and parallel...](https://github.com/JuliaParallel/Dagger.jl)**
>
> A framework for out-of-core and parallel execution - JuliaParallel/Dagger.jl
It’s not difficult to setup an overlay table so that every function call turns into a `Dagger.@spawn` call and is thus handled by a scheduler. That would effectively give you Bend. The reason why people don’t do that is because you’d just get slower code in most scenarios because of scheduler overhead.
Julia has been taking a step-by-step development to get there, starting from single-core performance and then really focusing on multithreading performance. Then extending the scheduler to support distributed scheduling with an appropriate cost model. You can dig up older documents that outline exactly the plans:
But of course, we’re still in the phase of optimizing Tapir. You’ll notice some familiar names:
> **[Tapir | Proceedings of the 22nd ACM SIGPLAN Symposium on Principles and...](https://dl.acm.org/doi/10.1145/3018743.3018758)**
>
>
> \## Introduction -- What is Tapir
> \[Tapir\](http://cilk.mit.edu/tapir/) is a paral…lel IR extension to LLVM. For the interested I recommend
> perusing the \[Tapir paper\](https://dl.acm.org/citation.cfm?doid=3018743.3018758). The key takeaway is that parallel (non-concurrent) programs, can be effectively model with cilk-style task parallelism and that given the serial-projection property (serial execution is always a valid execution), it is possible to reason about parallelism in the LLVM compiler.
>
> By doing so Tapir solves one primary problem: Traditionally introducing parallelism into a program, inhibits compiler optimisations. This is due to a variety of reasons, but chiefly that most implementations of parallelism choose to do early-outlining of parallel thunks. Causing the optimizer to only see calls into the runtime/program thunks without context. A classical optimisation that is inhibited by this is loop-invariant-code-movement. In Julia we encounter a different problem (#15276) in which using a closure to outline a thunk can cause performance issues.
>
> \### Tapir concepts
> \#### Syncregion
> An opaque token that is used to associate the various parallel IR statements with each other, so that during \`sync\` only synchronizes tasks that it is responsible for. Important for nested parallelism and inlining of functions containing parallel constructs.
>
> \#### Detach
> Think of this as a "function call" to the parallel region.
> \`detach within %syncregion, %label, %reattach\`. The \`label\` points to the basic-block that starts-off the parallel region and the \`reattach\` label points past a \`reattach\` statement and represents the execution on the task that is spawning the parallel region.
>
> \#### Reattach
> This is the "return" of a parallel region. It reattaches the parallel region to the original code and the \`label\` should point to the same basic-block that the \`reattach\` label in \`detach\` is pointing to.
>
> \#### Sync
> Synchronises all tasks with the same \`syncregion\`
>
> \## Goal of this PR
> This is very much ongoing research on how to best integrate the ideas from Tapir and the technology behind it into Julia. I want to lay a foundation on which we can build and experiment in the future. While the full-benefits will only be realised if one uses a Tapir enabled LLVM build, one
> of my goals is to bring the concepts of tapir into the Julia IR and thereby enable us to do optimizations on parallel code in the Julia IR even on a LLVM that doesn't have the Tapir extension. Right now we are in the very early stages of supporting Tapir in Julia.
>
> It is important to note that the semantics of this representation are parallel and not concurrent,
> by this extent this will not and cannot replace Julia Tasks. In order to exemplify this issue see the following Julia task code:
>
> \`\`\`julia
> @sync begin
> ch1 = Channel(0)
> ch2 = Channel(0)
> @async begin
> take!(ch1)
> put!(ch2, 1)
> end
> @async begin
> put!(ch1, 1)
> take!(ch2)
> end
> end
> \`\`\`
> Doing a serial projection of this code leads to a deadlock.
>
> \## User interface
> In \`test/tapir.jl\` I have placed some functions that I have been experimenting with. I do not expect users to directly use \`@syncregion\`, \`@spawn\` and \`@sync\_end\`, but rather I think the prototype implementation of a parallel for loop and \`@sync\`, \`@spawn\`.
>
> \`\`\`julia
> @par for i in 1:10
> ...
> end
>
> function fib(N)
> if N \<= 1
> return N
> end
> x = Ref{Int64}()
> @sync begin # different sync than Tasks
> @spawn begin
> x\[\] = fib(N-2)
> end
> y = fib(N-1)
> end
> return x\[\] + y
> end
> \`\`\`
>
> \## Changes/Current Status
> \- Buildsystem support for Tapir/LLVM
> \- New expr nodes:
> - \`syncregion\`: Obtain a token to synchronize spawned tasks
> - \`spawn\`: Spawn a block in a task
> - \`sync\`: Synchronize all tasks using the same token
> \- New IR nodes:
> - \`detach\`: Detach a parallel region
> - \`reattach\`: Join a parallel region
> \- Codegen support for \`syncregion\`, \`detach\`, \`reattach\`, \`sync\`
>
> \## Examples
>
> \## TODO:
> \- \[x\] loop information
> \- \[\] tests!!!
> \- \[\] \`fib2\`
> \- \[\] early lowering (in codegen) to PARTR
> \- \[\] late lowering as a LLVM pass to PARTR
> \- \[\] runtime support for GC/PTLS
> \- \[\] interpreter
> \- \[x\] cleanup PR
>
> \## Notes
> \### Make.user
> \`\`\`
> LLVM\_VER=svn
> USE\_TAPIR=1
> BUILD\_LLVM\_CLANG=1
> LLVM\_GIT\_VER="WIP-taskinfo"
> LLVM\_GIT\_VER\_CLANG="WIP-csi-tapir-exceptions"
> LLVM\_GIT\_VER\_COMPILER\_RT="WIP-cilksan-bugfixes"
> override CC=gcc-7
> override CXX=g++-7
> \`\`\`
>
> \## Acknowledgments
> Many thanks to T.B. Schardl (@neboat) for the many discussions around Tapir and LLVM.
Probably the most recent step towards this was:
>
>
> \*\*tl;dr\*\* How about adding optimizable task parallel API to Julia?
>
> \## Introdu…ction
>
> \### How to teach parallelism to the Julia compiler
>
> (If you've already seen @vchuravy's PR #31086, maybe this part is redundant.)
>
> How we currently implement the task parallel API in Julia introduces a couple of obstacles for supporting high-performance parallel programs. In particular, the compiler cannot analyze and optimize the child tasks in the context of the surrounding code. This limits the benefit parallel programs can obtain from existing analysis and optimizations like type inference and constant propagations. Furthermore, the notion of tasks in Julia supports complex concurrent communication mechanisms that imposes a hard limitation for the scheduler to implement an efficient scheduling strategy.
>
> \*Tapir\* (\[Schardl et al., 2019\](https://doi.org/10.1145/3365655)) is a parallel IR that can be added to existing IR in the SSA form. They demonstrated that Tapir can be added to LLVM (aka \_Tapir/LLVM\_; a part of \[OpenCilk\](https://cilk.mit.edu/)) relatively "easily" and existing programs written in Cilk can benefit from \*pre-existing\* optimizations such as loop invariant code motion (LICM), common sub-expression elimination (CSE), and others that were \*already written for serial programs\*. In principle, a similar strategy should work for any existing compiler with SSA IR developed for serial programs, including the Julia compiler. That is to say, with Tapir in the Julia IR, we should be able to unlock the optimizations in Julia for parallel programs.
>
> Tapir enables the parallelism in the compiler by limiting its focus to parallel programs with the \*serial-projection property\* (@vchuravy and I call this type of parallelism the \_may-happen in parallel parallelism\_ for clarity). Although this type of parallel programs cannot use unconstrained concurrency communication primitives (e.g., \`take!(channel)\`), it can be used for a vast majority of the parallel programs that are compute-intensive; i.e., the type of programs for which Julia is already optimized/targeting. Furthermore, having a natural way to express this type of computation can be beneficial not only for the compiler but also for the scheduler.
>
> \### A strategy for optimizable task parallel API
>
> This PR implements Tapir in the Julia compiler, in Julia. I've been working with @vchuravy and TB Schardl (@neboat) to extend and complete @vchuravy's PR #31086 that uses OpenCilk (which includes a fork of LLVM and clang) for supporting Tapir at LLVM level (Tapir/LLVM). This project still has to solve many obstacles since extracting out Julia tasks at a late stage of LLVM pass is very hard (you can see my \*\*VERY\*\* work-in-progress fork at https://github.com/cesmix-mit/julia). We observed that many benefits of Tapir can actually be realized at the level of Julia's compilation passes (see below). For example, type inference and constant propagation implemented in the Julia compiler can penetrate through \`@spawn\` code blocks with Tapir. So, I implemented Tapir in pure Julia and made it work without the dependency on OpenCilk. Since this can be done without taking on a dependency on a non-standard extension to the LLVM compilation pipeline, we think this is a good sweet spot for start adding parallelism support to Julia in a way fully integrated into the compiler.
>
> Although there are more works to be done to turn this into a mergeable/releasable (but experimental) state, I think it'd be beneficial to open a PR at this stage because to start discussing:
>
> \* if we want parallelism support integrated into the compiler
> \* what kind of frontend API we want
> \* implementation strategy
>
> I'm very interested in receiving feedback!
>
> \## Proposed API
>
> Here is an example usage of the Tapir API implemented in this PR for the moment. Following the tradition, it computes the Fibonacci number in parallel:
>
> \`\`\`julia
> using Base.Experimental: Tapir
>
> function fib(N)
> if N \<= 1
> return N
> end
> local x1, x2 # task output
> Tapir.@sync begin # -.
> Tapir.@spawn x1 = fib(N - 1) # child task # | syncregion
> x2 = fib(N - 2) # |
> end # -'
> return x1 + x2
> end
> \`\`\`
>
> \* \`Tapir.@sync begin ... end\` denotes the \_syncregion\_ in which Tapir tasks can be spawned.
> \* \`Tapir.@spawn ...\` denotes the code block that \_may\_ run in parallel with respect to other tasks (including the parent task; e.g., \`fib(N - 2)\` in the example above).
> \* The \_task output\_ variables, i.e., the variables that are accessed after the \_syncregion\_ must be declared with \`local\`. In the example above, \`x1\` and \`x2\` are the task output. Declaring with \`local\` is required because, just like the standard \`@sync\`, \`Tapir.@sync\` creates a scope (i.e., it is a \`let\` block). Task output variables can also be initialized before \`Tapir.@sync\`.
> \* The code written with \`Tapir\` is expected to have the \*serial projection property\* (see below for the motivation). In particular, the code must be valid after removing \`Tapir\`-related syntax (i.e., replacing \`Tapir.@sync\` and \`Tapir.@spawn\` with \`let\` blocks).
> \* It is undefined behavior (for now \[^1\]) to assign local variables in one task and read them in another task, since such behavior leads to data races.
> \* \`@goto\` into or out of the tasks is not allowed
>
> One of the important consideration is to keep this API very minimal to let us add more optimizations like OpenCilk-based lowering at a later stage of LLVM. I've designed this API first with my experimental branch of OpenCilk-based Tapir support. My OpenCilk-based implementation was not perfect but I think this is reasonably constrained to fully implement it.
>
> Loop-based API: Although OpenCilk and @vchuravy's original PR support parallel \`for\` loop, I propose to not include it for the initial support in \`Base\`. On one hand, it is straight forward to implement a simple parallel loop framework from just this API. On the other hand, there are a lot of consideration to be made in the design space if we want to have extensible data parallelism and I feel it's beyond the scope of this PR.
>
> \---
>
> \[^1\]: I think we can make it an error by analyzing this in the front end, in principle.
>
> \### Questions
>
> \#### More explicit task output declaration?
>
> Current handling of task output variables may be "too automatic" and makes reasoning about the code harder for Julia programmers. For example, consider following code that has no race for now:
>
> \`\`\`julia
> function ...
> # ... hundreds of lines (no code mentioning x) ...
> Tapir.@sync begin
> Tapir.@spawn begin
> x = f()
> h(x)
> end
> x = g()
> h(x)
> end
> ...
> end
> \`\`\`
>
> After sometimes, it may be re-written as
>
> \`\`\`julia
> function ...
> if ...
> x = ...
> return x
> end
> # ... hundreds of lines (no code mentioning x) ...
> Tapir.@sync begin
> Tapir.@spawn begin
> x = f()
> h(x)
> end
> x = g()
> h(x)
> end
> ...
> end
> \`\`\`
>
> This code now has a race since \`x\` is updated by the child and parent tasks. Although this particular example can be rejected by the compiler, it is impossible to do so in general. So, it may be worth considering adding a declaration of task outputs. Maybe something like
>
> \`\`\`julia
> Tapir.@output a b c
> Tapir.@sync begin
> Tapir.@spawn ...
> ...
> end
> \`\`\`
>
> or
>
> \`\`\`julia
> Tapir.@sync (a, b, c) begin
> Tapir.@spawn ...
> ...
> end
> \`\`\`
>
> The compiler can then compare and verify that its understanding of task output variables and the variables specified by the programmer.
>
> It would be even better if we can support a unified syntax like this for all types of tasks (\`@async\` and \`Threads.@spawn\`).
>
> \#### Bikeshedding the name
>
> On one hand, the name of the module \`Base.Experimental.Tapir\` is not entirely appropriate since Tapir is the concept for IR and not the user-facing API. On the other hand, we cannot come up with a more appropriate name for this. It is mainly because there is no crisp terminology for this type of parallelism (even though Cilk has been showing the importance of this approach with multiple aspects). Another name could be \`Base.Experimental.Parallel\`. But "parallel" can also mean distributed or GPU-based parallelism. Providing \`@psync\` and \`@pspawn\` macros directly from \`Experimental\` is another approach.
>
> \## Performance improvements
>
> \### Demo 1: type inference
>
> The return type of the example \`fib\` above and the simple threaded mapreduce example \`mapfold\` implemented in \`test/tapir.jl\` can be inferred even though they contain \`Tapir.@spawn\`:
>
> \`\`\`julia
> julia\> include("test/tapir.jl");
>
> julia\> @code\_typed fib(3)
> CodeInfo(
> ...
> ) =\> Int64
>
> julia\> @code\_typed mapfold(x -\> isodd(x) ? missing : x, +, 1:3)
> CodeInfo(
> ...
> ) =\> Union{Missing, Int64}
> \`\`\`
>
> As we all know, the improvement in type inference drastically change the performance when there is a tight loop following the parallel code (without a function boundary).
>
> \### Demo 2: constant propagation
>
> Here is another very minimal example for demonstrating performance benefit we observed. It (naively) computes the average on the sliding window in parallel:
>
> \`\`\`julia
> using Base.Experimental: Tapir
>
> @inline function avgfilter!(ys, xs, N)
> @assert axes(ys) == axes(xs)
> for offset in firstindex(xs)-1:lastindex(xs)-N
> y = zero(eltype(xs))
> for k in 1:N
> y += @inbounds xs\[offset + k\]
> end
> @inbounds ys\[offset + 1\] = y / N
> end
> return ys
> end
>
> function demo!(ys1, ys2, xs1, xs2)
> N = 32
> Tapir.@sync begin
> Tapir.@spawn avgfilter!(ys1, xs1, N)
> avgfilter!(ys2, xs2, N)
> end
> return ys1, ys2
> end
> \`\`\`
>
> For comparison, here are the same code using the current task system (\`Threads.@spawn\`) and the sequential version.:
>
> \`\`\`julia
> function demo\_current!(ys1, ys2, xs1, xs2)
> N = 32
> @sync begin
> Threads.@spawn avgfilter!(ys1, xs1, N)
> avgfilter!(ys2, xs2, N)
> end
> return ys1, ys2
> end
>
> function demo\_seq!(ys1, ys2, xs1, xs2)
> N = 32
> avgfilter!(ys1, xs1, N)
> avgfilter!(ys2, xs2, N)
> return ys1, ys2
> end
> \`\`\`
>
> We can then run the benchmarks with
>
> \`\`\`julia
> using BenchmarkTools
> suite = BenchmarkGroup()
> xs1 = randn(2^20)
> ys1 = zero(xs1)
> xs2 = randn(length(xs1))
> ys2 = zero(xs2)
> @assert demo!(zero(ys1), zero(ys2), xs1, xs2) == demo\_current!(ys1, ys2, xs1, xs2)
> @assert demo!(zero(ys1), zero(ys2), xs1, xs2) == demo\_seq!(ys1, ys2, xs1, xs2)
> suite\["tapir"\] = @benchmarkable demo!($ys1, $ys2, $xs1, $xs2)
> suite\["current"\] = @benchmarkable demo\_current!($ys1, $ys2, $xs1, $xs2)
> suite\["seq"\] = @benchmarkable demo\_seq!($ys1, $ys2, $xs1, $xs2)
> results = run(suite, verbose = true)
> \`\`\`
>
> With \`julia\` started with a single thread \[^2\], we get
>
> \`\`\`
> 3-element BenchmarkTools.BenchmarkGroup:
> tags: \[\]
> "tapir" =\> Trial(5.614 ms)
> "seq" =\> Trial(5.588 ms)
> "current" =\> Trial(23.537 ms)
> \`\`\`
>
> i.e., Tapir and sequential programs have identical performance while the performance of the code written with \`Threads.@spawn\` (current) is much worse. This is because Julia can propagate the constant across task boundaries with Tapir. Subsequently, since LLVM can see that the innermost loop has a fixed loop count, the loop can be unrolled and vectorized.
>
> It can be observed by introspecting the generated code:
>
> \`\`\`
> julia\> @code\_typed demo!(ys1, ys2, xs1, xs2)
> CodeInfo(
> 1 ── %1 = $(Expr(:syncregion))
> │ %2 = (Base.Tapir.taskgroup)()::Channel{Any}
> │ %3 = $(Expr(:new\_opaque\_closure, Tuple{}, false, Union{}, Any, opaque closure @0x00007fec1d9ce5e0 in Main, Core.Argument(2), Core.Ar
> gument(4), :(%1)))::Any
> ...
> ) =\> Nothing
>
> julia\> m = Base.unsafe\_pointer\_to\_objref(Base.reinterpret(Ptr{Cvoid}, 0x00007fec1d9ce5e0))
> opaque closure @0x00007fec1d9ce5e0 in Main
>
> julia\> Base.uncompressed\_ir(m)
> CodeInfo(
> ...
> │ ││┌ @ range.jl:740 within \`iterate'
> │ │││┌ @ promotion.jl:409 within \`=='
> │ ││││ %56 = (%51 === 32)::Bool
> │ │││└
> └
> ...
> \`\`\`
>
> i.e., \`N = 32\` is successfully propagated. On the other hand, in the current task system:
>
> \`\`\`
> julia\> @code\_typed demo\_current!(ys1, ys2, xs1, xs2)
> CodeInfo(
> ...
> │ %13 = π (32, Core.Const(32))
> │ %14 = %new(var"#14#15"{Vector{Float64}, Vector{Float64}, Int64}, ys1, xs1, %13)::Core.PartialStruct(var"#14#15"{Vector{Float64}, Vec
> tor{Float64}, Int64}, Any\[Vector{Float64}, Vector{Float64}, Core.Const(32)\])
> \`\`\`
>
> i.e., \`N = 32\` (\`%32\`) is captured as an \`Int64\`.
>
> Indeed, the performance of the sequential parts of the code is crucial for observing the speedup. With \`julia -t2\`, we see:
>
> \`\`\`
> 3-element BenchmarkTools.BenchmarkGroup:
> tags: \[\]
> "tapir" =\> Trial(2.799 ms)
> "seq" =\> Trial(5.583 ms)
> "current" =\> Trial(20.767 ms)
> \`\`\`
>
> I think an important aspect of this example is that even a "little bit" of compiler optimizations enabled on the Julia side can be enough for triggering optimizations on the LLVM side yielding a substantial effect.
>
> \---
>
> \[^2\]: Since this PR is mainly about compiler optimization and not about the scheduler, single-thread performance compared with sequential program is more informative than multi-thread performance.
>
> \### Demo 3: dead code elimination
>
> The optimizations that can be done with forward analysis such as type inference and constant propagation are probably implementable for the current threading system in Julia with reasonable amount of effort. However, optimizations such as dead code elimination (DCE) that require backward analysis may be significantly more challenging given unconstrained concurrency of Julia's \`Task\`. In contrast, enabling Tapir at Julia IR level automatically triggers Julia's DCE (which, in turn, can trigger LLVM's DCE):
>
> \`\`\`julia
> using Base.Experimental: Tapir
>
> @inline function eliminatable\_computation(xs)
> a = typemax(UInt)
> b = 0
> for x in xs
> a = (typemax(a) + a) ÷ ifelse(x == 0, 1, x) # arbitrary computation that LLVM can eliminate
> b += x
> end
> return (a, b)
> end
>
> function demo\_dce()
> local b1, b2
> Tapir.@sync begin
> Tapir.@spawn begin
> a1, b1 = eliminatable\_computation(UInt(1):UInt(33554432))
> end
> a2, b2 = eliminatable\_computation(UInt(3):UInt(33554432))
> end
> return b1 + b2 # a1 and a2 not used
> end
> \`\`\`
>
> In single thread \`julia\`, this takes 30 ms while an equivalent code with current threading system (i.e., replace \`Tapir.\` with \`Threads.\`) takes 250 ms.
>
> Note that Julia only has to eliminate \`b1\` and not the code inside \`eliminatable\_computation\`. The rest of DCE can happen inside LLVM (which does not have to understand Julia's parallelism).
>
> \## Motivation behind the restricted task semantics
>
> As explained in the proposed API, the serial projection property restrict the set of programs expressible with \`Tapir\`. Although we have some degree of interoperability (see below), reasoning and guaranteeing forward progress are easy only when the programmer stick with the existing patterns. In particular, it means that we will not be able to drop \`Threads.@spawn\` or \`@async\` for expressing programs with unconstrained concurrency. However, this restricted semantics is still enough for expressing compute-oriented programs and allows more optimizations in the compiler and the scheduler.
>
> As shown above, enabling Tapir at Julia level already unlocks some appealing set of optimizations for parallel programs. The serial projection property is useful for supporting optimizations that require backward analysis such as DCE in a straightforward manner. Once we manage to implement a proper OpenCilk integration, we expect to see the improvements from much richer set of existing optimizations at the LLVM level. This would have multiplicative effects when more LLVM-side optimizations are enabled (e.g., #36832). Furthermore, there is ongoing research on more cleverly using the parallel IR for going beyond unlocking pre-existing optimizations. For example, fusing arbitrary multiple \`Task\`s to be executed in a single \`Task\` can introduce deadlock. However, this is not the case in Tapir tasks thanks to the serial projection property. This, in turn, let us implement more optimizations inside the compiler such as a task coalescing pass that is aware of a downstream vecotrizer pass. In addition to optimizing user programs, this type of task improves productivity tools such as race detector (ref \[productivity tools provided by OpenCilk\](https://cilk.mit.edu/tools/)).
>
> In addition to the performance improvements on the side of the compiler, may-happen in parallel parallelism can help the parallel task runtime to handle the task scheduling cleverly. For example, having parallel IR makes it easy to implement continuation-stealing as used in Cilk (although it may not be compatible with the depth-first approach in Julia). Another possibility is to use a version of \`schedule\` that \_may fail\_ to reduce contention when there are a lot of small tasks (\[as I discussed here\](https://discourse.julialang.org/t/ann-foldsthreads-jl-a-zoo-of-pluggable-thread-based-data-parallel-execution-mechanisms/54662/5)). This is possible because we don't allow concurrency primitives in Tapir tasks and the task is not guaranteed to be executed in a dedicated \`Task\`. Since Tapir makes the call/task tree analyzable by the compiler, it may also help improving the depth-first scheduler.
>
> \## Composability with existing task system
>
> If we are going to have two notions of parallel tasks, it is important that the two systems can work together seamlessly. Of course, since the Tapir tasks cannot use arbitrary concurrency APIs, we can't support some task APIs inside \`Tapir.@spawn\` (e.g., \`take!(channel)\`). However, Tapir only requires the forward progress guarantee to be independent of other tasks within the same syncregion but not with respect to the code outside. Simplify put, we can invoke concurrency API as long as it does not communicate with other Tapir tasks in the same syncregion. For example, following code is valid:
>
> \`\`\`julia
> @sync begin
> Threads.@spawn begin
> Tapir.@sync begin
> Tapir.@spawn begin
> put!(bounded\_channel, 0) # Thunk 1
> end
> f() # Thunk 2
> end
> end
> Threads.@spawn take!(bounded\_channel) # Thunk 3
> end
> \`\`\`
>
> as long as \`f()\` does not use the \`bounded\_channel\`. That is to say, the author of this code guarantees the forward progress of Thunk 1 and Thunk 2 independent of each other but \_not\_ with respect to Thunk 3. Therefore, the example above is a valid program.
>
> Another class of useful composable idiom is the use of concurrency API that unconditionally guarantees forward progress. For example, \`put!(unbounded\_channel, item)\` can make forward progress independent of other tasks (but this is not true for \`take!\`). This is also true for \`schedule(::Task)\`. Thus, invoking \`Threads.@spawn\` inside \`Tapir.@sync\` is valid if (bot not only if \[^3\]) we do not invoke \`wait(task)\` in the \`Tapir.@sync\`. For example, the following code is valid
>
> \`\`\`julia
> unbuffered\_channel = Channel(0)
> @sync begin
> Tapir.@sync begin
> Tapir.@spawn begin
> t1 = Threads.@spawn put!(unbuffered\_channel, 0)
> # wait(t1)
> end
> t2 = Threads.@spawn take!(unbuffered\_channel)
> # wait(t2)
> end
> end
> \`\`\`
>
> However, it is invalid if \`wait(t1)\` and \`wait(t2)\` are uncommented.
>
> \---
>
> \[^3\]: Note that invoking \`wait\` in a Tapir task can be valid in some cases. For example, if it is known that the set of the tasks spawned by \`Threads.@spawn\` eventually terminate independent of any forward progress in other Tapir task in the same syncregion, it is valid to \`wait\` on these tasks. For example, it is valid:
>
> \`\`\`julia
> unbuffered\_channel = Channel(0)
> Tapir.@sync begin
> Tapir.@spawn begin
> @sync begin
> Threads.@spawn put!(unbuffered\_channel, 0)
> take!(unbuffered\_channel)
> end
> end
> do\_something\_else() # does not touch \`unbuffered\_channel\`
> end
> \`\`\`
>
> However, this is not an example of the idiom of using API that unconditionally guarantees forward progress.
>
> \## Implementation strategy
>
> \### Outlining
>
> At the end of Julia's optimization passes (in \`run\_passes\`) the child tasks are outlined into \*opaque closures\* and wrapped into a \`Task\` (see \`lower\_tapir!(ir::IRCode)\` function). (Aside: @Keno's opaque closure was \*very\* useful for implementing outlining at a late stage! Actually, @vchuravy has been suggesting opaque closure would be the enabler for this type of transformation since the beginning of the project. But it's nice to see how things fit together in a real code.) The outlined tasks are spawned and synced using the helper functions defined in \`Base.Tapir\`.
>
> \### Questions
>
> \* We need to create a new closure at the end of Julia's optimization phase. I'm using \`jl\_new\_code\_info\_uninit\` and \`jl\_make\_opaque\_closure\_method\` for this. Are they allowed to be invoked in this part of the compiler?
> \* In this PR, task outputs are handled by manual ad-hoc reg2mem-like passes (which is not very ideal). Currently (Julia 1.7-DEV), it looks like \[left over slots are rejected by the verifier\](https://github.com/JuliaLang/julia/blob/3129a5bef56bb7216024ae606c02b413b00990e3/base/compiler/ssair/verify.jl#L47-L48). Would it make sense to allow slots in the SSA IR? Can code other than Tapir make use of it? For example, maybe slots can work as an alternative to \`Core.Box\` when combined with opaque closure?
>
> \## Acknowledgment
>
> Many thanks to @vchuravy for pre-reviewing the PR!
>
> \## TODOs
>
> To reviewers: please feel free to add new ones or move the things in wishlist to here
>
> \- \[\] Decide the API
> \- \[x\] Pass tests of existing packages (\[TapirBenchmarks.jl\](https://github.com/cesmix-mit/TapirBenchmarks.jl) and \[FoldsTapir.jl\](https://github.com/JuliaFolds/FoldsTapir.jl))
> \- \[\] Resolve all TODO/ASK comments in the code
>
> \## Wishlist
>
> This is a list of nice-to-have things that are maybe not strictly required.
>
> \- IR verifier for checking the invariance of Tapir (e.g., disallow \`@goto\` into a task)
> \- Proper data flow analysis integration (do we need to use slot in the optimizer?)
> \- Support nested syncregion (or at least error out in the frontend?)
> \- Refine exception handling (esp. in the continuation)
> \- More tests
> \- Documentation (docstrings and \`manual/multi-threading.md\`)
> \- Easy sub-CFG-to-opaque-closure interface? (useful for outlining of exceptions for GPUs)
Basically APIs that let someone opt functions into may-happen parallelism.
The part of this which ended up getting deployed the soonest was ModelingToolkit, since that took a lot of the principles and then applied to do a domain-specific space where more assumptions could be made to simplify the problem. This was the actual first topic of the symbolic-numeric toolset:
[ ](https://www.youtube.com/watch?v=UNkXNZZ3hSw)
> [@MilesCranmer](#):
>
> but [HVM2](https://github.com/HigherOrderCO/HVM) is something that Julia could actually use (or maybe a macro).
The non-trivial part would be the cost model, which would need to be heavily tailored to the language. I don’t know how much you’d gain by pulling the rest of the scheduler along.
---
### Author:  [@MilesCranmer](https://discourse.julialang.org/u/MilesCranmer)
#### Post date: [May 19, 2024, 4:51am UTC](https://discourse.julialang.org/t/bend-a-new-gpu-native-language/114440/10 "2024-05-19T04:51:43Z")
I still don’t see anything related to GPUs though? `Dagger.@spawn` does not give you CUDA kernels… CUDA.jl still involves writing CUDA-like code.
The GPU part is what is so crazy about HVM. You write high-level code – _which looks nothing like a CUDA kernel_ – and get PTX out. It is not limited to vectorized array operations either.
> [@ChrisRackauckas](#):
>
> It’s not difficult to setup an overlay table so that every function call turns into a `Dagger.@spawn` call and is thus handled by a scheduler. That would effectively give you Bend.
I think you might be underestimating things a bit here… Author says they’ve been working on this for 10 years – [x.com](https://x.com/VictorTaelin/status/1791213162525524076). HVM2 is the fourth iteration of their framework. It’s really sounds not so easy to do this.
---
### Author:  [@Benny](https://discourse.julialang.org/u/Benny)
#### Post date: [May 19, 2024, 6:56am UTC](https://discourse.julialang.org/t/bend-a-new-gpu-native-language/114440/11 "2024-05-19T06:56:35Z")
> [@MilesCranmer](#):
>
> But this is a full compiler that can generate massively parallel GPU kernels from high-level code. I haven’t seen anything like this before.
> [@MilesCranmer](#):
>
> You write high-level code – _which looks nothing like a CUDA kernel_ – and get PTX out
What is the line you’re drawing here exactly? CUDA C++ is technically high level code compiled to PTX. They got tutorials of traversing and generating trees too.
---
### Author:  [@MilesCranmer](https://discourse.julialang.org/u/MilesCranmer)
#### Post date: [May 19, 2024, 7:18am UTC](https://discourse.julialang.org/t/bend-a-new-gpu-native-language/114440/12 "2024-05-19T07:18:24Z")
> [@](#):
>
> C++ is technically high level code
But we are on a Julia forum…
I guess maybe you are saying “high level” is a relative concept, but we are speaking relative to CUDA rather than relative to PTX.
Not to mention there are many things which are extremely difficult to express in CUDA (on a GPU in general I suppose) or with a library of array operations. Which is where something like HVM is proposed to help.
---
### Author:  [@xiaoxi](https://discourse.julialang.org/u/xiaoxi)
#### Post date: [May 19, 2024, 7:33am UTC](https://discourse.julialang.org/t/bend-a-new-gpu-native-language/114440/13 "2024-05-19T07:33:11Z")
> [@jling](#):
>
> someone pointed out that it’s partly due to their baseline being very slow.
According to the HVM2 paper:
> In single-thread CPU evaluation, HVM2, is, baseline, still about 5x slower than GHC, and this number can grow to 100x on programs that involve loops and mutable arrays, since HVM2 doesn’t feature these yet.
> [@ChrisRackauckas](#):
>
> The part of this which ended up getting deployed the soonest was ModelingToolkit, since that took a lot of the principles and then applied to do a domain-specific space where more assumptions could be made to simplify the problem. This was the actual first topic of the symbolic-numeric toolset:
It’s a shame that so many interesting projects in the Julia ecosystem don’t get more promotion. Fortunately, Calculus with Julia has been promoted today in Hacker News with great success.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 19, 2024, 1:08pm UTC](https://discourse.julialang.org/t/bend-a-new-gpu-native-language/114440/14 "2024-05-19T13:08:52Z")
> [@MilesCranmer](#):
>
> I think you might be underestimating things a bit here… Author says they’ve been working on this for 10 years – [x.com](https://x.com/VictorTaelin/status/1791213162525524076). HVM2 is the fourth iteration of their framework. It’s really sounds not so easy to do this.
Well Julia v1.0 was also the 4th iteration in some sense, going back to Star-P in 2004, Julia initial release, Julia experimental Distributed and Multithreading, then finally Tapir-based parallelism and constructs now going beyond that. So we’re going on 20 years 😅. I don’t think there’s much point to such counting, but this ain’t the first rodeo around here.
> **[25\_Edelman\_Pres.pdf](https://archive.ll.mit.edu/HPEC/agendas/proc06/Day2/25_Edelman_Pres.pdf)**
>
> 4.62 MB
> **[interactive\_supercomputing2.pdf](https://math.mit.edu/~edelman/publications/interactive_supercomputing2.pdf)**
>
> 81.21 KB
> **[ADA433445.pdf](https://apps.dtic.mil/sti/pdfs/ADA433445.pdf)**
>
> 402.78 KB
> [@MilesCranmer](#):
>
> The GPU part is what is so crazy about HVM. You write high-level code – _which looks nothing like a CUDA kernel_ – and get PTX out. It is not limited to vectorized array operations either.
Julia code can compile down to give custom kernels in CUDA.jl. There’s then abstractions written on top of that with tools like KernelAbstractions.jl
> **[GitHub - JuliaGPU/KernelAbstractions.jl: Heterogeneous programming in Julia](https://github.com/JuliaGPU/KernelAbstractions.jl)**
>
> Heterogeneous programming in Julia. Contribute to JuliaGPU/KernelAbstractions.jl development by creating an account on GitHub.
For example, the matmul kernel is effectively just Julia code with an added piece for describing the global index:
```julia
@kernel function matmul_kernel!(output, a, b)
i, j = @index(Global, NTuple)
# creating a temporary sum variable for matrix multiplication
tmp_sum = zero(eltype(output))
for k = 1:size(a)[2]
tmp_sum += a[i,k] * b[k, j]
end
output[i,j] = tmp_sum
end
```
> [@MilesCranmer](#):
>
> Not to mention there are many things which are extremely difficult to express in CUDA (on a GPU in general I suppose) or with a library of array operations. Which is where something like HVM is proposed to help.
With Bend if I’m not mistaken you still have to figure out how to represent the code using bend and fold constructs. Note this is pretty close to what I was linking to before with Tapir extensions, where Taka’s proposed Tapir extensions were coupled with transducer-type parallelism approaches.
> **[GitHub - JuliaFolds/FoldsTapir.jl: A PoC Tapir-JuliaFolds interface](https://github.com/JuliaFolds/FoldsTapir.jl)**
>
> A PoC Tapir-JuliaFolds interface. Contribute to JuliaFolds/FoldsTapir.jl development by creating an account on GitHub.
>
>
> \*\*tl;dr\*\* How about adding optimizable task parallel API to Julia?
>
> \## Introdu…ction
>
> \### How to teach parallelism to the Julia compiler
>
> (If you've already seen @vchuravy's PR #31086, maybe this part is redundant.)
>
> How we currently implement the task parallel API in Julia introduces a couple of obstacles for supporting high-performance parallel programs. In particular, the compiler cannot analyze and optimize the child tasks in the context of the surrounding code. This limits the benefit parallel programs can obtain from existing analysis and optimizations like type inference and constant propagations. Furthermore, the notion of tasks in Julia supports complex concurrent communication mechanisms that imposes a hard limitation for the scheduler to implement an efficient scheduling strategy.
>
> \*Tapir\* (\[Schardl et al., 2019\](https://doi.org/10.1145/3365655)) is a parallel IR that can be added to existing IR in the SSA form. They demonstrated that Tapir can be added to LLVM (aka \_Tapir/LLVM\_; a part of \[OpenCilk\](https://cilk.mit.edu/)) relatively "easily" and existing programs written in Cilk can benefit from \*pre-existing\* optimizations such as loop invariant code motion (LICM), common sub-expression elimination (CSE), and others that were \*already written for serial programs\*. In principle, a similar strategy should work for any existing compiler with SSA IR developed for serial programs, including the Julia compiler. That is to say, with Tapir in the Julia IR, we should be able to unlock the optimizations in Julia for parallel programs.
>
> Tapir enables the parallelism in the compiler by limiting its focus to parallel programs with the \*serial-projection property\* (@vchuravy and I call this type of parallelism the \_may-happen in parallel parallelism\_ for clarity). Although this type of parallel programs cannot use unconstrained concurrency communication primitives (e.g., \`take!(channel)\`), it can be used for a vast majority of the parallel programs that are compute-intensive; i.e., the type of programs for which Julia is already optimized/targeting. Furthermore, having a natural way to express this type of computation can be beneficial not only for the compiler but also for the scheduler.
>
> \### A strategy for optimizable task parallel API
>
> This PR implements Tapir in the Julia compiler, in Julia. I've been working with @vchuravy and TB Schardl (@neboat) to extend and complete @vchuravy's PR #31086 that uses OpenCilk (which includes a fork of LLVM and clang) for supporting Tapir at LLVM level (Tapir/LLVM). This project still has to solve many obstacles since extracting out Julia tasks at a late stage of LLVM pass is very hard (you can see my \*\*VERY\*\* work-in-progress fork at https://github.com/cesmix-mit/julia). We observed that many benefits of Tapir can actually be realized at the level of Julia's compilation passes (see below). For example, type inference and constant propagation implemented in the Julia compiler can penetrate through \`@spawn\` code blocks with Tapir. So, I implemented Tapir in pure Julia and made it work without the dependency on OpenCilk. Since this can be done without taking on a dependency on a non-standard extension to the LLVM compilation pipeline, we think this is a good sweet spot for start adding parallelism support to Julia in a way fully integrated into the compiler.
>
> Although there are more works to be done to turn this into a mergeable/releasable (but experimental) state, I think it'd be beneficial to open a PR at this stage because to start discussing:
>
> \* if we want parallelism support integrated into the compiler
> \* what kind of frontend API we want
> \* implementation strategy
>
> I'm very interested in receiving feedback!
>
> \## Proposed API
>
> Here is an example usage of the Tapir API implemented in this PR for the moment. Following the tradition, it computes the Fibonacci number in parallel:
>
> \`\`\`julia
> using Base.Experimental: Tapir
>
> function fib(N)
> if N \<= 1
> return N
> end
> local x1, x2 # task output
> Tapir.@sync begin # -.
> Tapir.@spawn x1 = fib(N - 1) # child task # | syncregion
> x2 = fib(N - 2) # |
> end # -'
> return x1 + x2
> end
> \`\`\`
>
> \* \`Tapir.@sync begin ... end\` denotes the \_syncregion\_ in which Tapir tasks can be spawned.
> \* \`Tapir.@spawn ...\` denotes the code block that \_may\_ run in parallel with respect to other tasks (including the parent task; e.g., \`fib(N - 2)\` in the example above).
> \* The \_task output\_ variables, i.e., the variables that are accessed after the \_syncregion\_ must be declared with \`local\`. In the example above, \`x1\` and \`x2\` are the task output. Declaring with \`local\` is required because, just like the standard \`@sync\`, \`Tapir.@sync\` creates a scope (i.e., it is a \`let\` block). Task output variables can also be initialized before \`Tapir.@sync\`.
> \* The code written with \`Tapir\` is expected to have the \*serial projection property\* (see below for the motivation). In particular, the code must be valid after removing \`Tapir\`-related syntax (i.e., replacing \`Tapir.@sync\` and \`Tapir.@spawn\` with \`let\` blocks).
> \* It is undefined behavior (for now \[^1\]) to assign local variables in one task and read them in another task, since such behavior leads to data races.
> \* \`@goto\` into or out of the tasks is not allowed
>
> One of the important consideration is to keep this API very minimal to let us add more optimizations like OpenCilk-based lowering at a later stage of LLVM. I've designed this API first with my experimental branch of OpenCilk-based Tapir support. My OpenCilk-based implementation was not perfect but I think this is reasonably constrained to fully implement it.
>
> Loop-based API: Although OpenCilk and @vchuravy's original PR support parallel \`for\` loop, I propose to not include it for the initial support in \`Base\`. On one hand, it is straight forward to implement a simple parallel loop framework from just this API. On the other hand, there are a lot of consideration to be made in the design space if we want to have extensible data parallelism and I feel it's beyond the scope of this PR.
>
> \---
>
> \[^1\]: I think we can make it an error by analyzing this in the front end, in principle.
>
> \### Questions
>
> \#### More explicit task output declaration?
>
> Current handling of task output variables may be "too automatic" and makes reasoning about the code harder for Julia programmers. For example, consider following code that has no race for now:
>
> \`\`\`julia
> function ...
> # ... hundreds of lines (no code mentioning x) ...
> Tapir.@sync begin
> Tapir.@spawn begin
> x = f()
> h(x)
> end
> x = g()
> h(x)
> end
> ...
> end
> \`\`\`
>
> After sometimes, it may be re-written as
>
> \`\`\`julia
> function ...
> if ...
> x = ...
> return x
> end
> # ... hundreds of lines (no code mentioning x) ...
> Tapir.@sync begin
> Tapir.@spawn begin
> x = f()
> h(x)
> end
> x = g()
> h(x)
> end
> ...
> end
> \`\`\`
>
> This code now has a race since \`x\` is updated by the child and parent tasks. Although this particular example can be rejected by the compiler, it is impossible to do so in general. So, it may be worth considering adding a declaration of task outputs. Maybe something like
>
> \`\`\`julia
> Tapir.@output a b c
> Tapir.@sync begin
> Tapir.@spawn ...
> ...
> end
> \`\`\`
>
> or
>
> \`\`\`julia
> Tapir.@sync (a, b, c) begin
> Tapir.@spawn ...
> ...
> end
> \`\`\`
>
> The compiler can then compare and verify that its understanding of task output variables and the variables specified by the programmer.
>
> It would be even better if we can support a unified syntax like this for all types of tasks (\`@async\` and \`Threads.@spawn\`).
>
> \#### Bikeshedding the name
>
> On one hand, the name of the module \`Base.Experimental.Tapir\` is not entirely appropriate since Tapir is the concept for IR and not the user-facing API. On the other hand, we cannot come up with a more appropriate name for this. It is mainly because there is no crisp terminology for this type of parallelism (even though Cilk has been showing the importance of this approach with multiple aspects). Another name could be \`Base.Experimental.Parallel\`. But "parallel" can also mean distributed or GPU-based parallelism. Providing \`@psync\` and \`@pspawn\` macros directly from \`Experimental\` is another approach.
>
> \## Performance improvements
>
> \### Demo 1: type inference
>
> The return type of the example \`fib\` above and the simple threaded mapreduce example \`mapfold\` implemented in \`test/tapir.jl\` can be inferred even though they contain \`Tapir.@spawn\`:
>
> \`\`\`julia
> julia\> include("test/tapir.jl");
>
> julia\> @code\_typed fib(3)
> CodeInfo(
> ...
> ) =\> Int64
>
> julia\> @code\_typed mapfold(x -\> isodd(x) ? missing : x, +, 1:3)
> CodeInfo(
> ...
> ) =\> Union{Missing, Int64}
> \`\`\`
>
> As we all know, the improvement in type inference drastically change the performance when there is a tight loop following the parallel code (without a function boundary).
>
> \### Demo 2: constant propagation
>
> Here is another very minimal example for demonstrating performance benefit we observed. It (naively) computes the average on the sliding window in parallel:
>
> \`\`\`julia
> using Base.Experimental: Tapir
>
> @inline function avgfilter!(ys, xs, N)
> @assert axes(ys) == axes(xs)
> for offset in firstindex(xs)-1:lastindex(xs)-N
> y = zero(eltype(xs))
> for k in 1:N
> y += @inbounds xs\[offset + k\]
> end
> @inbounds ys\[offset + 1\] = y / N
> end
> return ys
> end
>
> function demo!(ys1, ys2, xs1, xs2)
> N = 32
> Tapir.@sync begin
> Tapir.@spawn avgfilter!(ys1, xs1, N)
> avgfilter!(ys2, xs2, N)
> end
> return ys1, ys2
> end
> \`\`\`
>
> For comparison, here are the same code using the current task system (\`Threads.@spawn\`) and the sequential version.:
>
> \`\`\`julia
> function demo\_current!(ys1, ys2, xs1, xs2)
> N = 32
> @sync begin
> Threads.@spawn avgfilter!(ys1, xs1, N)
> avgfilter!(ys2, xs2, N)
> end
> return ys1, ys2
> end
>
> function demo\_seq!(ys1, ys2, xs1, xs2)
> N = 32
> avgfilter!(ys1, xs1, N)
> avgfilter!(ys2, xs2, N)
> return ys1, ys2
> end
> \`\`\`
>
> We can then run the benchmarks with
>
> \`\`\`julia
> using BenchmarkTools
> suite = BenchmarkGroup()
> xs1 = randn(2^20)
> ys1 = zero(xs1)
> xs2 = randn(length(xs1))
> ys2 = zero(xs2)
> @assert demo!(zero(ys1), zero(ys2), xs1, xs2) == demo\_current!(ys1, ys2, xs1, xs2)
> @assert demo!(zero(ys1), zero(ys2), xs1, xs2) == demo\_seq!(ys1, ys2, xs1, xs2)
> suite\["tapir"\] = @benchmarkable demo!($ys1, $ys2, $xs1, $xs2)
> suite\["current"\] = @benchmarkable demo\_current!($ys1, $ys2, $xs1, $xs2)
> suite\["seq"\] = @benchmarkable demo\_seq!($ys1, $ys2, $xs1, $xs2)
> results = run(suite, verbose = true)
> \`\`\`
>
> With \`julia\` started with a single thread \[^2\], we get
>
> \`\`\`
> 3-element BenchmarkTools.BenchmarkGroup:
> tags: \[\]
> "tapir" =\> Trial(5.614 ms)
> "seq" =\> Trial(5.588 ms)
> "current" =\> Trial(23.537 ms)
> \`\`\`
>
> i.e., Tapir and sequential programs have identical performance while the performance of the code written with \`Threads.@spawn\` (current) is much worse. This is because Julia can propagate the constant across task boundaries with Tapir. Subsequently, since LLVM can see that the innermost loop has a fixed loop count, the loop can be unrolled and vectorized.
>
> It can be observed by introspecting the generated code:
>
> \`\`\`
> julia\> @code\_typed demo!(ys1, ys2, xs1, xs2)
> CodeInfo(
> 1 ── %1 = $(Expr(:syncregion))
> │ %2 = (Base.Tapir.taskgroup)()::Channel{Any}
> │ %3 = $(Expr(:new\_opaque\_closure, Tuple{}, false, Union{}, Any, opaque closure @0x00007fec1d9ce5e0 in Main, Core.Argument(2), Core.Ar
> gument(4), :(%1)))::Any
> ...
> ) =\> Nothing
>
> julia\> m = Base.unsafe\_pointer\_to\_objref(Base.reinterpret(Ptr{Cvoid}, 0x00007fec1d9ce5e0))
> opaque closure @0x00007fec1d9ce5e0 in Main
>
> julia\> Base.uncompressed\_ir(m)
> CodeInfo(
> ...
> │ ││┌ @ range.jl:740 within \`iterate'
> │ │││┌ @ promotion.jl:409 within \`=='
> │ ││││ %56 = (%51 === 32)::Bool
> │ │││└
> └
> ...
> \`\`\`
>
> i.e., \`N = 32\` is successfully propagated. On the other hand, in the current task system:
>
> \`\`\`
> julia\> @code\_typed demo\_current!(ys1, ys2, xs1, xs2)
> CodeInfo(
> ...
> │ %13 = π (32, Core.Const(32))
> │ %14 = %new(var"#14#15"{Vector{Float64}, Vector{Float64}, Int64}, ys1, xs1, %13)::Core.PartialStruct(var"#14#15"{Vector{Float64}, Vec
> tor{Float64}, Int64}, Any\[Vector{Float64}, Vector{Float64}, Core.Const(32)\])
> \`\`\`
>
> i.e., \`N = 32\` (\`%32\`) is captured as an \`Int64\`.
>
> Indeed, the performance of the sequential parts of the code is crucial for observing the speedup. With \`julia -t2\`, we see:
>
> \`\`\`
> 3-element BenchmarkTools.BenchmarkGroup:
> tags: \[\]
> "tapir" =\> Trial(2.799 ms)
> "seq" =\> Trial(5.583 ms)
> "current" =\> Trial(20.767 ms)
> \`\`\`
>
> I think an important aspect of this example is that even a "little bit" of compiler optimizations enabled on the Julia side can be enough for triggering optimizations on the LLVM side yielding a substantial effect.
>
> \---
>
> \[^2\]: Since this PR is mainly about compiler optimization and not about the scheduler, single-thread performance compared with sequential program is more informative than multi-thread performance.
>
> \### Demo 3: dead code elimination
>
> The optimizations that can be done with forward analysis such as type inference and constant propagation are probably implementable for the current threading system in Julia with reasonable amount of effort. However, optimizations such as dead code elimination (DCE) that require backward analysis may be significantly more challenging given unconstrained concurrency of Julia's \`Task\`. In contrast, enabling Tapir at Julia IR level automatically triggers Julia's DCE (which, in turn, can trigger LLVM's DCE):
>
> \`\`\`julia
> using Base.Experimental: Tapir
>
> @inline function eliminatable\_computation(xs)
> a = typemax(UInt)
> b = 0
> for x in xs
> a = (typemax(a) + a) ÷ ifelse(x == 0, 1, x) # arbitrary computation that LLVM can eliminate
> b += x
> end
> return (a, b)
> end
>
> function demo\_dce()
> local b1, b2
> Tapir.@sync begin
> Tapir.@spawn begin
> a1, b1 = eliminatable\_computation(UInt(1):UInt(33554432))
> end
> a2, b2 = eliminatable\_computation(UInt(3):UInt(33554432))
> end
> return b1 + b2 # a1 and a2 not used
> end
> \`\`\`
>
> In single thread \`julia\`, this takes 30 ms while an equivalent code with current threading system (i.e., replace \`Tapir.\` with \`Threads.\`) takes 250 ms.
>
> Note that Julia only has to eliminate \`b1\` and not the code inside \`eliminatable\_computation\`. The rest of DCE can happen inside LLVM (which does not have to understand Julia's parallelism).
>
> \## Motivation behind the restricted task semantics
>
> As explained in the proposed API, the serial projection property restrict the set of programs expressible with \`Tapir\`. Although we have some degree of interoperability (see below), reasoning and guaranteeing forward progress are easy only when the programmer stick with the existing patterns. In particular, it means that we will not be able to drop \`Threads.@spawn\` or \`@async\` for expressing programs with unconstrained concurrency. However, this restricted semantics is still enough for expressing compute-oriented programs and allows more optimizations in the compiler and the scheduler.
>
> As shown above, enabling Tapir at Julia level already unlocks some appealing set of optimizations for parallel programs. The serial projection property is useful for supporting optimizations that require backward analysis such as DCE in a straightforward manner. Once we manage to implement a proper OpenCilk integration, we expect to see the improvements from much richer set of existing optimizations at the LLVM level. This would have multiplicative effects when more LLVM-side optimizations are enabled (e.g., #36832). Furthermore, there is ongoing research on more cleverly using the parallel IR for going beyond unlocking pre-existing optimizations. For example, fusing arbitrary multiple \`Task\`s to be executed in a single \`Task\` can introduce deadlock. However, this is not the case in Tapir tasks thanks to the serial projection property. This, in turn, let us implement more optimizations inside the compiler such as a task coalescing pass that is aware of a downstream vecotrizer pass. In addition to optimizing user programs, this type of task improves productivity tools such as race detector (ref \[productivity tools provided by OpenCilk\](https://cilk.mit.edu/tools/)).
>
> In addition to the performance improvements on the side of the compiler, may-happen in parallel parallelism can help the parallel task runtime to handle the task scheduling cleverly. For example, having parallel IR makes it easy to implement continuation-stealing as used in Cilk (although it may not be compatible with the depth-first approach in Julia). Another possibility is to use a version of \`schedule\` that \_may fail\_ to reduce contention when there are a lot of small tasks (\[as I discussed here\](https://discourse.julialang.org/t/ann-foldsthreads-jl-a-zoo-of-pluggable-thread-based-data-parallel-execution-mechanisms/54662/5)). This is possible because we don't allow concurrency primitives in Tapir tasks and the task is not guaranteed to be executed in a dedicated \`Task\`. Since Tapir makes the call/task tree analyzable by the compiler, it may also help improving the depth-first scheduler.
>
> \## Composability with existing task system
>
> If we are going to have two notions of parallel tasks, it is important that the two systems can work together seamlessly. Of course, since the Tapir tasks cannot use arbitrary concurrency APIs, we can't support some task APIs inside \`Tapir.@spawn\` (e.g., \`take!(channel)\`). However, Tapir only requires the forward progress guarantee to be independent of other tasks within the same syncregion but not with respect to the code outside. Simplify put, we can invoke concurrency API as long as it does not communicate with other Tapir tasks in the same syncregion. For example, following code is valid:
>
> \`\`\`julia
> @sync begin
> Threads.@spawn begin
> Tapir.@sync begin
> Tapir.@spawn begin
> put!(bounded\_channel, 0) # Thunk 1
> end
> f() # Thunk 2
> end
> end
> Threads.@spawn take!(bounded\_channel) # Thunk 3
> end
> \`\`\`
>
> as long as \`f()\` does not use the \`bounded\_channel\`. That is to say, the author of this code guarantees the forward progress of Thunk 1 and Thunk 2 independent of each other but \_not\_ with respect to Thunk 3. Therefore, the example above is a valid program.
>
> Another class of useful composable idiom is the use of concurrency API that unconditionally guarantees forward progress. For example, \`put!(unbounded\_channel, item)\` can make forward progress independent of other tasks (but this is not true for \`take!\`). This is also true for \`schedule(::Task)\`. Thus, invoking \`Threads.@spawn\` inside \`Tapir.@sync\` is valid if (bot not only if \[^3\]) we do not invoke \`wait(task)\` in the \`Tapir.@sync\`. For example, the following code is valid
>
> \`\`\`julia
> unbuffered\_channel = Channel(0)
> @sync begin
> Tapir.@sync begin
> Tapir.@spawn begin
> t1 = Threads.@spawn put!(unbuffered\_channel, 0)
> # wait(t1)
> end
> t2 = Threads.@spawn take!(unbuffered\_channel)
> # wait(t2)
> end
> end
> \`\`\`
>
> However, it is invalid if \`wait(t1)\` and \`wait(t2)\` are uncommented.
>
> \---
>
> \[^3\]: Note that invoking \`wait\` in a Tapir task can be valid in some cases. For example, if it is known that the set of the tasks spawned by \`Threads.@spawn\` eventually terminate independent of any forward progress in other Tapir task in the same syncregion, it is valid to \`wait\` on these tasks. For example, it is valid:
>
> \`\`\`julia
> unbuffered\_channel = Channel(0)
> Tapir.@sync begin
> Tapir.@spawn begin
> @sync begin
> Threads.@spawn put!(unbuffered\_channel, 0)
> take!(unbuffered\_channel)
> end
> end
> do\_something\_else() # does not touch \`unbuffered\_channel\`
> end
> \`\`\`
>
> However, this is not an example of the idiom of using API that unconditionally guarantees forward progress.
>
> \## Implementation strategy
>
> \### Outlining
>
> At the end of Julia's optimization passes (in \`run\_passes\`) the child tasks are outlined into \*opaque closures\* and wrapped into a \`Task\` (see \`lower\_tapir!(ir::IRCode)\` function). (Aside: @Keno's opaque closure was \*very\* useful for implementing outlining at a late stage! Actually, @vchuravy has been suggesting opaque closure would be the enabler for this type of transformation since the beginning of the project. But it's nice to see how things fit together in a real code.) The outlined tasks are spawned and synced using the helper functions defined in \`Base.Tapir\`.
>
> \### Questions
>
> \* We need to create a new closure at the end of Julia's optimization phase. I'm using \`jl\_new\_code\_info\_uninit\` and \`jl\_make\_opaque\_closure\_method\` for this. Are they allowed to be invoked in this part of the compiler?
> \* In this PR, task outputs are handled by manual ad-hoc reg2mem-like passes (which is not very ideal). Currently (Julia 1.7-DEV), it looks like \[left over slots are rejected by the verifier\](https://github.com/JuliaLang/julia/blob/3129a5bef56bb7216024ae606c02b413b00990e3/base/compiler/ssair/verify.jl#L47-L48). Would it make sense to allow slots in the SSA IR? Can code other than Tapir make use of it? For example, maybe slots can work as an alternative to \`Core.Box\` when combined with opaque closure?
>
> \## Acknowledgment
>
> Many thanks to @vchuravy for pre-reviewing the PR!
>
> \## TODOs
>
> To reviewers: please feel free to add new ones or move the things in wishlist to here
>
> \- \[\] Decide the API
> \- \[x\] Pass tests of existing packages (\[TapirBenchmarks.jl\](https://github.com/cesmix-mit/TapirBenchmarks.jl) and \[FoldsTapir.jl\](https://github.com/JuliaFolds/FoldsTapir.jl))
> \- \[\] Resolve all TODO/ASK comments in the code
>
> \## Wishlist
>
> This is a list of nice-to-have things that are maybe not strictly required.
>
> \- IR verifier for checking the invariance of Tapir (e.g., disallow \`@goto\` into a task)
> \- Proper data flow analysis integration (do we need to use slot in the optimizer?)
> \- Support nested syncregion (or at least error out in the frontend?)
> \- Refine exception handling (esp. in the continuation)
> \- More tests
> \- Documentation (docstrings and \`manual/multi-threading.md\`)
> \- Easy sub-CFG-to-opaque-closure interface? (useful for outlining of exceptions for GPUs)
> **[GitHub - JuliaFolds/Transducers.jl: Efficient transducers for Julia](https://github.com/JuliaFolds/Transducers.jl)**
>
> Efficient transducers for Julia. Contribute to JuliaFolds/Transducers.jl development by creating an account on GitHub.
You can see a demonstration of it here:
[https://juliafolds.github.io/data-parallelism/tutorials/quick-introduction/](https://juliafolds.github.io/data-parallelism/tutorials/quick-introduction/)
And it includes a FoldsCUDA.jl for a version with CUDA support. The idea was to integrate the DAG construction into the compiler (that’s the PR) and give a macro so users can opt-in easily (as opposed to being fully parallel, so that the general scheduling problem does not have to be solved).
But there are some issues with the bend/fold approach. It’s not new, and Guy Steele is probably the person you can watch who has had the most on this. His early parallel programming language Fortress is one of the ones that heavily influenced Julia’s designs. He has a nice talk on the limits of such an approach:
[ ](https://vimeo.com/6624203)
And was a keynote at an early JuliaCon:
[ ](https://www.youtube.com/watch?v=EZD3Scuv02g)
This is a particular project that I would be interested in reviving if I ever find the right person and grant funding for it (though I’m not the right person to do the day-to-day reviews on it 😅). I think integrating bend/fold with opt-in may-happen parallelism (i.e. opt into allowing the compiler to choose parallelism or not, and how) would be a nice thing to have, though I’m personally skeptical of the number of codes I have that it could match my parallelism on so I tend to keep this stuff off “the critical path” for now.
---
### Author:  [@MilesCranmer](https://discourse.julialang.org/u/MilesCranmer)
#### Post date: [May 19, 2024, 1:49pm UTC](https://discourse.julialang.org/t/bend-a-new-gpu-native-language/114440/15 "2024-05-19T13:49:39Z")
Maybe the missing piece here: bend _works_. It’s a full programming language that can compile to correct CUDA kernels. It’s not a limited set of constructs that are already easy to write kernels for (FoldsCUDA.jl) or a CUDA wrapper (CUDA.jl), but rather a full GPU-native language that looks like regular code, and it _currently_ generates correct PTX. (and as the saying goes, the last 20% takes 80% of the time)
Of course the ideas for this are not new, the author themselves cites a 1997 paper for the design. But the fact they actually got it working is novel.
I’d love to be mistaken of course, but last I checked I still need to write CUDA kernels for any operations that doesn’t fit nicely into broadcasting or reduction. While CUDA.jl lets you use Julia inside the kernel, it still needs to look like and act like a regular CUDA kernel. And with CUDA kernels there are some things that are hugely complicated to do well, to the point I don’t even bother trying.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 19, 2024, 3:24pm UTC](https://discourse.julialang.org/t/bend-a-new-gpu-native-language/114440/16 "2024-05-19T15:24:16Z")
> [@MilesCranmer](#):
>
> It’s not a limited set of constructs that are already easy to write kernels for ([FoldsCUDA.jl](https://juliahub.com/ui/Packages/General/FoldsCUDA)) or a CUDA wrapper ([CUDA.jl](https://juliahub.com/ui/Packages/General/CUDA)), but rather a full GPU-native language that looks like regular code, and it _currently_ generates correct PTX. (and as the saying goes, the last 20% takes 80% of the time)
Ehh, marketing. Read the actual docs, not the marketing material. Bend takes lineage from those kind of DAG-builder approaches but mixes some of Clojure’s pieces in:
> **[Language: Concurrency and Parallelism](https://clojure-doc.org/articles/language/concurrency_and_parallelism/)**
>
> This guide covers:
Where all data structures are immutable and recursive
> Finally, let’s get straight to the fun part: how do we implement parallel algorithms with Bend? Just kidding. Before we get there, let’s talk about loops. You might have noticed we have avoided them so far. That wasn’t by accident. There is an important aspect on which Bend diverges from Python, and aligns with Haskell: **variables are immutable**. Not “by default”. They just **are**. For example, in Bend, we’re not allowed to write:
>
> Which is immutable. If that sounds annoying, that’s because **it is**. Don’t let anyone tell you otherwise. We are aware of that, and we have many ideas on how to improve this, making Bend feel even more Python-like. For now, we have to live with it. But, wait… if variables are immutable… how do we even do loops?
The solution is of course, transducers. Wait, I mean bend and fold.
> Fold: consuming recursive datatypes
> Bend: generating recursive datatypes
Bend and fold are just transducers. Transducers.jl comes from the same place.
> Transducers.jl provides composable algorithms on “sequence” of inputs. They are called _[transducers](https://clojure.org/reference/transducers)_ , first introduced in Clojure language by Rich Hickey.
> **[GitHub - JuliaFolds/Transducers.jl: Efficient transducers for Julia](https://github.com/JuliaFolds/Transducers.jl)**
>
> Efficient transducers for Julia. Contribute to JuliaFolds/Transducers.jl development by creating an account on GitHub.
JuliaFolds then built constructs on top of those because most people find writing transducer-based programs quite hard to do:
> **[GitHub - JuliaFolds/Folds.jl: A unified interface for sequential, threaded,...](https://github.com/JuliaFolds/Folds.jl)**
>
> A unified interface for sequential, threaded, and distributed fold - JuliaFolds/Folds.jl
So then it becomes a skeleton-based parallel framework on top of a transducer-based immutable programming layer. Which again I mentioned before as one of the research directions that was ongoing (that I would like to revive). One of the things required to make this more performant in Julia (and thus the barriers for now) include optimizations for immutability and improved escape analysis, which are some major projects right now with the new Memory type of v1.11 and to re-building of constructs on that piece.
So… “rather a full GPU-native language that looks like regular code”, if your regular code only consists of immutable recursive data structures manipulated without loops but instead by transducers then yes, its a full language that looks like regular code! I happen to use loops from time to time, and arrays, and heaps, and some other non-recursive structures. So to me, this is not “regular code”
This of course was the same major adoption barrier to Clojure.
> [@MilesCranmer](#):
>
> Maybe the missing piece here: bend _works_.
No it doesn’t. Parallelism works when it speeds things up. It currently doesn’t speed up code.
But there’s also many projects like this. DawnCC was a C compiler that attempted to automate parallel constructs:
[https://homepages.dcc.ufmg.br/~fernando/publications/papers\_pt/Breno16Tools.pdf](https://homepages.dcc.ufmg.br/~fernando/publications/papers_pt/Breno16Tools.pdf)
Bones was another such project:
> **[GitHub - tue-es/bones: Research compiler based on algorithmic skeletons](https://github.com/tue-es/bones)**
>
> Research compiler based on algorithmic skeletons. Contribute to tue-es/bones development by creating an account on GitHub.
It did well on codes that were amenable to a loop analysis but had difficulties speeding things up that were more general. These types of parallelism models have been called “Skeleton-based programming models”. Other languages to look at in this space include:
> **[SkePU 2: Flexible and Type-Safe Skeleton Programming for Heterogeneous...](https://link.springer.com/article/10.1007/s10766-017-0490-5)**
>
> In this article we present SkePU 2, the next generation of the SkePU C++ skeleton programming framework for heterogeneous parallel systems. We critically examine the design and limitations of the SkePU 1 programming interface. We present a new,...
[https://skelcl.github.io/](https://skelcl.github.io/)
[https://muesli.uni-muenster.de](https://muesli.uni-muenster.de)
[https://parallelme.github.io/](https://parallelme.github.io/)
> **[Valente\_2022.pdf](https://run.unl.pt/bitstream/10362/155512/1/Valente_2022.pdf)**
>
> 5.72 MB
There’s a project from NIST that I’m missing too…
Most are skeleton-based because that’s usually a good way to get something working (i.e. fast enough to do some example well) before taking off to be a whole thing.
Probably the thing I would point to that worked the most in this space so far is hiCUDA, which if you read their papers you’ll see very similar language:
> We have designed hiCUDA,
> a high-level directive-based language for CUDA programming. It allows programmers to perform these tedious tasks
> in a simpler manner, and directly to the sequential code.
Its promise was that you could write something pretty close to normal C++ but it would automatically construct the parallel code. The main difference from Bend is that it required you control memory allocations so that you could enforce locality, as this was pretty necessary for the “compiler-guessed CUDA” to get close to the handwritten optimal CUDA. But getting the data handling right is hard, so it was an interesting thing to just opt-out of that part and only focus the automation on the code under the assumption the user will locate the data properly.
> [@MilesCranmer](#):
>
> I’d love to be mistaken of course, but last I checked I still need to write CUDA kernels for any operations that doesn’t fit nicely into broadcasting or reduction. While [CUDA.jl](https://juliahub.com/ui/Packages/General/CUDA) lets you use Julia inside the kernel, it still needs to look like and act like a regular CUDA kernel. And with CUDA kernels there are some things that are hugely complicated to do well, to the point I don’t even bother trying.
There have been many codes to automatically construct DAGs to be alternatively compiled. Those projects have effectively been thrown away because they weren’t actually fast, just like Bend. The pieces that you actually see shared and remaining are pieces with opt-in parallelism because again, anything that was not opt-in was not able to guarantee it wouldn’t slow down normal execution. Removing the opt-in nature is a choice, not a technical hurdle. Just no one has found a good enough solution to the scheduling problem to justify moving to a more automated interface. And clearly Bend hasn’t solved that problem either.
---
### Author:  [@wsmoses](https://discourse.julialang.org/u/wsmoses)
#### Post date: [May 19, 2024, 3:34pm UTC](https://discourse.julialang.org/t/bend-a-new-gpu-native-language/114440/17 "2024-05-19T15:34:45Z")
FWIW I’m experimenting with an XLA-based Julia fuser in [GitHub - EnzymeAD/Reactant.jl](https://github.com/EnzymeAD/Reactant.jl)
---
### Author:  [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)
#### Post date: [May 19, 2024, 3:43pm UTC](https://discourse.julialang.org/t/bend-a-new-gpu-native-language/114440/18 "2024-05-19T15:43:09Z")
The most exciting (to me) approach to high level GPU programming is based on, SOACS [https://futhark-lang.org/](https://futhark-lang.org/)
I suspect that some of the futhark compilation rules could be used in a Julia dedicated package but I did not find time (yet) to explore this.
---
### Author:  [@MilesCranmer](https://discourse.julialang.org/u/MilesCranmer)
#### Post date: [May 19, 2024, 4:27pm UTC](https://discourse.julialang.org/t/bend-a-new-gpu-native-language/114440/19 "2024-05-19T16:27:16Z")
@ChrisRackauckas thanks for all of these links. Of course Transducers.jl and Folds.jl are awesome but they definitely constrain what sorts of code you can write. I think most things possible on GPU via Transducers.jl and Folds.jl would be fairly easy to write as explicit CUDA kernels already (and likely faster when written in explicit CUDA).
Again, the GPU-ification of arbitrary high-level code is the cool part.
Imagine writing a compiler that could be run natively on GPUs. That seems impossibly difficult in CUDA and CUDA wrappers like CUDA/CUDAFolds/etc, but seems possible in a high-level framework like this.
> [@ChrisRackauckas](#):
>
> No it doesn’t. Parallelism works when it speeds things up. It currently doesn’t speed up code.
I would push back on this as Bend seems to already get good scaling:
> [@](#):
>
> - CPU, Apple M3 Max, 1 thread: **12.15 seconds**
> - CPU, Apple M3 Max, 16 threads: **0.96 seconds**
> - GPU, NVIDIA RTX 4090, 16k threads: **0.21 seconds**
It only got released two days ago though. I’m just excited it seems to scale like this at all!
> [@ChrisRackauckas](#):
>
> But there’s also many projects like this
I took a look through these but (a) the GPU-compatible ones look highly-constrained in terms of language features, and (b) the more flexible ones seem incompatible with GPUs. Bones looks the most related but seems like it is missing some language features, and also it has been abandoned for 10 years, sadly.
> [@wsmoses](#):
>
> FWIW I’m experimenting with an XLA-based Julia fuser in [GitHub - EnzymeAD/Reactant.jl](https://github.com/EnzymeAD/Reactant.jl)
Thanks, looks useful!
> [@LaurentPlagne](#):
>
> The most exciting (to me) approach to high level GPU programming is based on, SOACS [https://futhark-lang.org/](https://futhark-lang.org/)
Thanks, also looks cool. Although I will note this language declares itself to be an “array language” (which means the CUDA kernels are likely easy to write by hand in many cases)
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 19, 2024, 4:31pm UTC](https://discourse.julialang.org/t/bend-a-new-gpu-native-language/114440/20 "2024-05-19T16:31:32Z")
> [@MilesCranmer](#):
>
> Of course [Transducers.jl](https://juliahub.com/ui/Packages/General/Transducers) and [Folds.jl](https://juliahub.com/ui/Packages/General/Folds) are awesome but they definitely constrain what sorts of code you can write. I think most things possible on GPU via [Transducers.jl](https://juliahub.com/ui/Packages/General/Transducers) and [Folds.jl](https://juliahub.com/ui/Packages/General/Folds) would be fairly easy to write as explicit CUDA kernels already (and likely faster when written in explicit CUDA).
Why don’t you say Bend also constrains the codes you can write? It has the same constraints as folds and transcuders, right? The difference is that Bend is a language where only fold and transducer constructs exist, while Transducers.jl is a system inside of Julia which requires that you use only folds and transducers to get the parallelism. But given they have the same constraints (which I pointed exactly to the points in the documentation which say this), why would you not say Bend also constrains the way you write code?
> [@MilesCranmer](#):
>
> I would push back on this as Bend seems to already get good scaling:
Scaling isn’t performance. Remember, NetworkX is more scalable than LightGraphs.jl, there just happens to be no size graph that fits on a modern computer for which NetworkX is faster.
[Next page](https://discourse.julialang.org/t/bend-a-new-gpu-native-language/114440.md?page=2)