GPU Offloading¶

This document describes how LFortran decides where a parallel loop runs when a GPU backend is selected with --gpu=metal, --gpu=cuda, --gpu=cuda_cpu or --gpu=hip. It is a design document: where the current implementation differs, the differences are listed at the end.

Goals¶

Correctness first. An offloaded loop computes exactly what the same loop computes on the CPU.
Everything is offloaded. Every do concurrent loop runs on the device, no matter how slow the lowering is. A loop that the compiler cannot offload yet is a compiler bug to fix, never a reason to quietly use the CPU.
No silent CPU execution. A loop runs on the CPU only for a short, explicitly documented list of constructs, and only when the user asks for it.
Performance decisions come later. Choosing the CPU because it is faster for a given loop is a separate, later feature, built on top of being able to offload everything correctly.

Outcomes¶

Every parallel loop compiled with a GPU backend ends in exactly one of these outcomes:

Situation	Without `--gpu-allow-cpu-fallback`	With `--gpu-allow-cpu-fallback`
The loop uses nothing from the unsupported list	offloaded	offloaded
The loop uses something from the unsupported list	error, reported by the unsupported-construct check	the loop runs on the CPU, with a warning
The offloading pipeline fails on the loop	error (a compiler bug)	error (a compiler bug)

--gpu-allow-cpu-fallback therefore softens only the unsupported list. It never hides a failure of the offloading pipeline.

The unsupported list¶

The unsupported list names the constructs that LFortran does not offload on a given device, either not yet or never. It is intentionally short, and every entry states why it is there. Current proposal:

Construct	Devices	Why	Plan
`real(8)`	Metal	the Metal Shading Language has no 64-bit float	emulate later
`real(10)`, `real(16)`	all	no device has these widths	emulate, or never
input/output statements, including internal-file I/O	all	no Fortran units or formats on a device	never, or a host-side buffer later

error stop and stop are not on the list: every device runs them, through a status buffer the host reads after the launch (see Stopping on the device). That covers a stop in a procedure the loop calls (purity is not enforced yet, see the open questions) and in an offloaded !$omp parallel do loop. A stop directly in a do concurrent body is invalid Fortran -- it is an image control statement, which the standard does not allow there and gfortran rejects -- and is to become a semantic error on every target, not an entry of this list; until then it runs like any other stop.

Everything else is not on the list, even when it is not implemented yet. For example integer(8) on Metal (the Metal Shading Language has 64-bit integers), complex, character, run-time sized locals or derived types with allocatable components are all offloaded, or are an error until the pipeline supports them.

Adding an entry to the list is a design decision, not a way to make a failing loop compile.

Device identifiers¶

Valid Fortran names such as half, kernel, int and device need not be valid Metal identifiers. The unique_symbols pass gives Metal device data and derived types collision-free names in their defining scopes, updating component references and binding scalar record component arguments by symbol. It preserves types, layout, host ABI and kernel registration names, without enabling otherwise ineligible global mangling. Device helper declarations and calls use their existing emission names. CUDA (including cuda_cpu emulation) uses the same symbol-based renaming for C/C++ keywords such as int, float and class, without applying Metal-only naming restrictions. Ordinary CPU compilation is unchanged.

Host launch expressions retain host component definitions, not aliases retargeted to device copies. In particular, workspace extents such as size(records(selected)%component, dim) are bound to the actual host array's derived type before device identifiers are renamed.

Scalar logarithms¶

Real log in a kernel or a device-only procedure retains its intrinsic identity through the intrinsic-function pass and lowers to the device's native log operation. This includes transitive helper calls. Host procedures retain their normal runtime lowering; Metal still requires real(real32).

Optional captures¶

Optional captures do not suppress a GPU launch: nonoptional work still runs when a captured argument is absent. The optional-argument ASR pass supplies the kernel's explicit presence flag, including when forwarding to device procedures. Host buffer addresses and shape queries use a safe view associated with the actual argument only when it is present; absent arrays use an empty view. Decomposed allocatable components use the root argument's presence before associating their host views or reading their bounds. Scalar ABI placeholders likewise never read absent storage. These placeholders are not Fortran default values: guarded device accesses retain the original absence semantics. Views preserve proven contiguity without a host array copy; potentially strided arguments still use the normal copy-in/copy-out path.

GPU synchronization recognizes the address-taken backing storage in the generated ASR associations, preserving host read/write and scope-exit synchronization for captured non-TARGET arguments. These body-local aliases do not change declared TARGET attributes or the separately compiled dummy ABI.

The unsupported-construct check¶

The check decides from the program as written, before any lowering, whether a loop uses something on the unsupported list. It does not try to lower the loop and it does not depend on how far lowering would get.

It runs once per parallel loop, when the execution target is chosen (the parallel_dispatch pass, see exec_target):

If the loop passes, it is assigned to the device. From then on the offloading pipeline is committed to it.
If the loop fails, it is an error at the offending construct (the I/O statement, the real(8) variable), or, with --gpu-allow-cpu-fallback, the loop is assigned to the host with a warning naming that construct.

What the check has to look at¶

A do concurrent loop is constrained by the Fortran standard, which keeps the check small:

A procedure referenced in a do concurrent construct must be pure. A pure procedure cannot do external input/output or stop. So external I/O and stop can only appear directly in the loop body.
A pure procedure can still do internal-file I/O (write to a character variable), so that has to be looked for in the called procedures too.
gfortran (-std=f2018) rejects stop in a do concurrent body as an image control statement, while print in the body and error stop in a pure procedure are accepted. (To be confirmed against the standard text.)

Types, on the other hand, can come from anywhere the loop reaches: the loop body, the called procedures, their locals and the data passed to the device.

For an !$omp parallel do loop offloaded with --gpu-offload-omp-loops the standard gives no such guarantee: the body and its callees can do anything, so the check has to walk every procedure the loop reaches.

Since purity is not enforced yet (see the open questions), the check does not rely on it: for every loop it walks the loop body and every procedure the loop reaches transitively, looking at their statements, the declared types of their dummy arguments, results and locals, and the type of every expression they evaluate. A compile-time constant is not evaluated on the device, so what is below it is not looked at.

For derived types the check uses this rule:

A component reference s%m is judged by the type of m alone: only that component reaches the device. A loop that reads only real(4) components of a type that also has a real(8) component is not rejected for it.
A derived-type value used as a whole (passed to a procedure, assigned, copied), and a derived-type variable declared in a procedure the loop reaches or in a block of the loop, is judged by all of its components, recursively: its whole layout reaches the device.

The check lives in src/libasr/pass/gpu_unsupported_check.cpp, and each failing loop is reported once, at the first construct found, with the loop also pointed at when the construct is in a called procedure. Every failing loop of a file is reported before the compilation stops.

--show-gpu-kernel-source assigns a failing loop to the host, as --gpu-allow-cpu-fallback does, so that the kernels of the other loops are still shown. A failure of the offloading pipeline is an error there too.

The offloading pipeline¶

Once a loop is assigned to the device, the offloading pipeline extracts the kernel, lowers it through the shared ASR passes, lays out its arguments and generates the device source. It never falls back to the CPU:

There is no CPU alternative kept alongside the kernel, and nothing to undo.
Outlining copies nested block and associate scopes with their bodies. Calls inside nested parallel regions name those copied scopes and the kernel's copies of their helper procedures, not the host's symbols.
A construct the pipeline cannot lower yet is a compiler error at the loop ("not yet implemented"). It is fixed like any other compiler bug: reduce it to a minimal reproducer, fix the pipeline, add an integration test.
The pipeline does not repeat the unsupported-construct check. Where it needs a fact the check guarantees (for example that no real(8) reaches a Metal kernel), it asserts it.
An associate name bound to a variable -- a scalar, an element, a component or a section -- is an alias of it, in the loop body and in a procedure the loop calls, so what is written through the name reaches the variable: the selector takes the place of a scalar name (and of every name in the loop body, also inside the block and associate constructs in it), and an array name of a procedure points at the selector. So the kernel is handed the variable once, however many names the loop reads and writes it through (associate (w => self%w) with self also used in the loop). A selector designates the object it designates on entry to the construct, so a subscript of it that the construct may change (associate (e => w(j)); j = 3; e = 1) is read into a variable of its own on entry.
An associate name of a construct the loop is in, bound to an expression (associate (n => size(a))), is the value the expression had when the construct was entered, which the kernel is handed as a scalar: the expression is not evaluated again, and what it reads is not handed over.
A call of a function that returns an array, in an array statement of the loop body (a(1:n, p) = x(p)%values() * 2, sum(x(p)%values()), an operand of matmul), is made once, before the statement, and its result assigned to an array of the iteration, which the statement then reads element by element; so is the call an associate name of a construct in the loop is bound to (associate (y => x(p)%values())), on entry. The array, like any allocatable of device code assigned the result of a function, takes the shape the function gives its result, over the actual arguments of the call.
An array of a size known at compile time with unspecified locality that an iteration defines whole before it references it (t = ..., or do k = 1, m; t(k) = ... over its whole declared extent, before anything reads t), on every path through the body that references it, is a temporary of the iteration, as Fortran 2018, 11.1.7.5 allows: each thread gets its own copy, a local array, as if t were declared in a block around the body. This holds whether every iteration, only some, or none take such a path, since none reads a value it did not define. It is only done when no statement that can run after the loop -- the rest of the enclosing statements, the condition of an enclosing do while, or wherever an exit or cycle goes -- may read the array before defining it again, so the array keeps the value it had before the loop and nothing can tell. Locality is a correctness decision, independent of array size or stack capacity. After array temporaries, inlining and materialization, gpu_kernel_storage materializes updating aggregate arguments and places the actual lowered arrays, before gpu_memory_space specializes device calls. Small fixed-size arrays use Thread storage within the device stack budget (64 KiB with Metal, 512 KiB with CUDA), reserving 8 KiB for device compiler frames and register spills and accounting for the remaining out-of-line device routines, including their by-value aggregate parameters. Non-spillable aggregate locals and argument copies consume that budget too. Larger arrays use Global storage, recorded as per-thread workspaces in Function.gpu by gpu_kernel_finalize. Existing launch sizing and indexing give every iteration a separate slice; spilling never turns a private array back into a shared value. Every other array of unspecified locality -- one whose elements the iterations define one each (y(j) = ...), one partly defined, one an iteration reads before defining it, a dummy argument, an array that may be read after the loop, or one a contained procedure, a pointer, an associate around the loop or a locality clause names -- is one array that all threads share. Runtime-sized eligible arrays use the same workspace machinery as automatic BLOCK locals: the host evaluates launch-resolvable extents, including the existing replay over thread indices, and allocates a slice per iteration. Privatized automatic arrays capture their existing bounds and extents on the host; they do not re-evaluate specification expressions whose operands may have changed since entry to the procedure. If an extent cannot be rebuilt on the host, the workspace pre-flight reports a compile-time error; the compiler never silently leaves that array shared. Workspace memory therefore grows with the launch size.
A scalar with unspecified locality that the loop assigns is copied back to the host after the launch only when a statement that can run after the loop may read the value it is left with, by the same reading as for the arrays above: a later t = 0, or a later loop every iteration of which assigns t first, does not. It is kept when an iteration may read it before assigning it, which needs the value it had before the loop, and whenever the reading cannot follow every access to it (a pointer, an equivalence, a contained procedure or an associate around the loop). A loop of one iteration leaves the value it assigns defined (Fortran 2018, 11.1.7.5), so a later read keeps it.
An array an iteration needs whose size is known only at run time -- a local of a block with run-time extents, an allocatable local, or a temporary the compiler introduces, such as the result of a function or the copy an overlapping assignment needs -- lives in a buffer the host allocates before the launch, when the host can size it (the fast path), and on the device heap otherwise (see Allocating on the device). Every thread uses a slice of a host-sized buffer, one slice per iteration, which only the device holds: what a slice holds on entry is undefined, as for any local, and nothing is copied in or back. When the size changes from one iteration to the next (allocate(w(m(k))), or x(k)%values() whose result has the size of a component of the iteration's element), the host works it out for every iteration before the launch, from the loop limits and the values the kernel gives the loop indices, and gives every slice room for the largest. Each thread uses the leading part its own size needs. The size may read the loop indices, the arguments of the loop and the shapes of arrays. A size that is the same in every iteration may also call a function without side effects of the translation unit, such as the one an array intrinsic (maxval(n), sum(n)) is lowered to, which the host calls before the launch as the device does. The copy an overlapping assignment needs (b(1:n(l),l) = ... b(1:n(l+1),l+1) ..., anywhere in the body) is sized by the whole of each dimension of the assigned array whose extent changes inside the iteration, and uses its leading part. A size may not read a value the loop writes, which the host does not know before the launch: such an array is allocated on the device heap instead when it is allocatable, and is an error otherwise ("workspace ... cannot be sized on the host", or "a temporary array cannot be sized on the host" for a temporary the compiler introduces).
A procedure the loop calls runs on the device as a procedure, wherever the call is: in an expression, an associate selector, a condition, or another procedure. A device without variable-length arrays or a heap (Metal, CUDA) cannot declare a local array of such a procedure whose size is known only at run time -- an allocatable local, an automatic array, or a temporary the compiler introduces, such as the result of a function it calls. Each such local becomes an argument of the procedure, and every call passes an array of its own for it, with the extents of the local written over the actual arguments of that call. Up the chain of calls this ends in the kernel, where the array is a per-thread workspace like any other (see the previous point). The extents may read the arguments of the procedure and named constants, and are written as they are on entry to the procedure: a name it binds once (an associate name, a scalar assigned once and changed through no other name) stands for what it is bound to, the extent of a section for the extents of its ranges, and a function without side effects is called with arguments of that kind (allocate(a(maxval(n))) with associate (n => self%nodes_)). The value of an intent(in) argument is fixed for the call unless it is a pointer, whose target the procedure may change. An extent that is evaluated after entry -- of a local of a block, or of the array a call is given for its callee -- may not read a dummy argument the procedure can change before (n = n + 4), whose value the callers do not know. A procedure with such a local that calls itself, directly or not, would need one array per level of the recursion, and is an error, and so is an extent that reads anything else ("workspace ... cannot be sized on the host"), except for an allocatable local, which the procedure allocates on the device heap itself. A call that is the whole right-hand side of an assignment is spliced into the kernel instead: the body of the procedure replaces the call, with the actual arguments in place of the dummy arguments, which more often makes the extents constants. An allocatable dummy array that the procedure neither allocates, deallocates, asks about with allocated, nor hands on to an allocatable dummy, is handed over as a plain array, with the lower bounds of the actual argument as arguments of their own. Inside device code an allocatable local is given the shape of the allocate, reallocation or move that gives it one, also when that shape is read from another such local, and when several of them give it the same shape, as shared lowering does when it allocates a temporary and then reallocates it to the shape it already has. That shape has to be written as it is on entry, as an extent is, so that it is one value at every allocation; any other allocatable of device code is an error. A reallocation by an assignment gives the shape of the array the assigned value takes its shape from: r = 2*c(2:) gives r the shape of c(2:), also when r is the result of a function. A rank-one allocatable local whose extent is only known where it is allocated -- inside a loop of the procedure, as the temporaries of v(1:n(l)) * a(1:n(l), l-1) in do l = 2, m -- points instead to the leading part of a buffer, which is a run-time sized local like those above, of an upper bound of every such extent known on entry: the extent of a section whose subscripts read the loop index is bounded by the extent of its array in that dimension, as a section never has more elements there. The local itself keeps the exact extent of each allocation.
present(x) of an optional dummy argument of the procedure the loop is in is evaluated by the host before the launch, and the kernel reads the result like any other scalar of the host: the loop runs whether the argument is present or not. An allocatable the loop only uses when it is allocated (if (present(cost)) c(i) = ..., with c allocated only when cost is present) is handed to the kernel as no storage when it is not allocated. A loop that names an optional argument other than in present is launched only when the argument is present, for now.
An internal procedure the loop calls reads the variables and named constants of its host by host association, which device code does not have. The kernel's copy of the procedure is given them instead: a named constant becomes a named constant of the procedure, and a variable a dummy argument, which every call hands the variable. A procedure that calls one that is handed a variable is handed it too, to pass it on.
A derived-type value that device code creates -- a local of a procedure, or the result of a function such as a structure constructor, associate (t => tensor_t(x)) -- has no storage for an allocatable array component on the device either. Each such component is stored in an array of its own, of the shape the procedure gives the component: the extents of its allocate, of the array an assignment reallocates it to, of the same component of a value assigned to the whole variable, or the shape a procedure gives the component of its intent(out) or intent(inout) argument when the value is passed for it, in terms of the actual arguments of that call. The shape is written as it is on entry to the procedure, over its arguments: a name the procedure binds once (an associate name, a scalar assigned once) stands for what it is bound to, and the extent of a section for the extents of its ranges. The array is then a run-time sized local like those above, given to the callers up to the kernel and sized on the host; a procedure the value is passed to is handed the array as the component, as it is handed the device buffer of a component of a kernel argument. An assignment of the whole value copies the component into that array. When a value of a type cannot be stored this way -- a component given several different shapes, or a shape read from a value the procedure computes, an array of values of the type, or a value held in a component of another value -- every value of that type device code creates is on the device heap instead (a heap type, see Allocating on the device). Component buffers do not necessarily share the address space of the derived-type value: a fixed-size result temporary can have thread-local buffers, while a host buffer or a lifted workspace is in global memory. After storage placement, gpu_memory_space reads the canonical component-to-storage Associate bindings. A device specialization that needs non-global component buffers takes them, their extents and their lower bounds as explicit ASR arguments, with ordinary array memory-space types. Its dummy components are bound to those arguments with Associate, so forwarding them to another routine preserves each buffer's space. Code generation emits those arguments rather than adding a second, implicitly global component interface.
An argument a procedure defines (intent(out) or intent(inout)), a scalar or a derived-type value, is passed by reference, to a value in the memory of the calling thread, as the procedure's locals are. An actual argument that lives in device memory instead -- an element of a kernel argument (call total(i, s1(i))), of an array dummy in device memory, or a derived-type kernel argument -- is copied into a local of the caller around the call; the allocatable components of a derived-type value are not copied, the procedure is handed their device buffers.
An array of device code is an address and, for each dimension, a lower bound, an extent and a stride, as a Fortran descriptor is, and every one of them reaches the code that indexes the array. An array that does not start at one -- an array argument of the kernel, an allocatable component of a derived-type argument, the component of an element of an array of a derived type, an allocatable dummy array -- is counted from its lower bounds: the host passes those of the component of an argument, and of each element's component in the buffer of its sizes, a routine gets those of an allocatable dummy as dummies of their own, a(lb:), and a name bound to any of these has them too. A name bound to a section -- associate (z => w(1:n, 1:m, l)), or the temporary shared lowering binds to a section it hands to a routine -- is the address of the section's first element, its extents, and the strides of the dimensions of the array it ranges over, times its steps; so is a descriptor dummy of a routine, which every call hands the strides of its actual argument. An assumed-shape dummy that some call hands such a name, or another descriptor dummy, stays a descriptor dummy -- that of the routine matmul(w(1:n, 1:m, l), x) is lowered to, or of a procedure of the program handed z -- while the shared passes turn the others into explicit-shape ones. An explicit-shape dummy is contiguous, so a section that is not contiguous is gathered into a contiguous temporary before such a call and scattered back after it, by the kernel and by a routine the kernel calls alike; the temporary is a run-time sized local as above.
A derived-type value reaches the kernel as a copy of the host's bytes, laid out by the device exactly as the host lays it out. An allocatable array component is handed over in device buffers of its own. So is one of a component that is a derived-type value itself (self%input_map_%slope_): the component (self%input_map_) becomes an argument of its own of every kernel and device procedure that uses it, read or passed on to a procedure, and every call and launch passes it beside the value it is a component of (gpu_struct_component_args.cpp), at any depth. An allocatable scalar component, a character(len=:), allocatable one included, is carried as the host's bytes (the address of the value, and its length for a character), which the device keeps in place and copies back unchanged, but never reads: a loop that reads only the other components of the value, also when it passes the whole value to a procedure, runs on the device. A loop that names such a component (s%name, allocated(s%r)), or assigns a whole value of a type holding one, which would have to allocate it, is an error for now ("the allocatable scalar component r on the device, which has no buffer of its own").
When a loop writes an allocatable array component of an element of an array of derived type that the host passes to the kernel, and the value written is not on the device heap, the storage the kernel writes into has to exist before the launch: the kernel writes into the buffers of the kernel argument, which cannot be reallocated during the launch. The write may be t(i) = f(...), where f gives the component of its result a size, t(i) = tt(...), or t(i)%v = .... That follows the rule every assignment follows: an allocatable scalar is always allocated or reallocated automatically, but an allocatable array, including an allocatable array component such as t(i)%v, only with --realloc-lhs-arrays.
- With the option, the host allocates t(i)%v before the launch, or reallocates it when it is allocated with another size.
- Without it, nothing is allocated, and the component has to be allocated with the right size already. With bounds checking on (the default, without --fast), the launch checks this with the same run-time check an assignment gets ("Array 't(i)%v' is not allocated", "Array shape mismatch in assignment"), pointing at the write. Allocation status and size are known only at run time, so the check cannot be made at compile time. Without bounds checking the kernel writes storage that does not exist.
- Only the elements the loop writes are allocated or checked. The offload pass builds this host code from the loop as the source wrote it, before it turns the loop into a kernel. The host code runs the loop's iterations again, but only what decides the writes: the if and select case tests around a write, the block and associate constructs it is in, and the scalar assignments at the top of the body that the tests, subscripts and sizes read. At each write, also when the loop writes the component in more than one place (in both branches of an if, in several cases of a select case, or twice in an iteration), it allocates or checks the element the write picks, with the size the write gives it: the extents of the array assigned to the component, or those of the allocate or array assignment in f (in the body a submodule gives it, for a separate module procedure), also inside a block or associate construct of f, in terms of the actual arguments of the call. A type-bound f called through an object of a type that is not polymorphic (net%infer(x)) is the procedure its binding names, with the object as its passed-object argument. An element the loop does not write is left as it is. The limits of the loops are evaluated once, before this host code runs, and the launch reads the same values.
- The host does not guess what it cannot work out before the loop runs:
  - The size. It may come from a call (also one whose result f assigns whole, r = g(x)), from a value the loop itself writes (s(i) = ... earlier in the iteration, so the value before the loop would be a wrong size), from a bound or an element of a dummy array of f, or from a dummy argument f changes (n = n + 1 before allocate(r%v(n)), where n is a value dummy, also through an associate name: associate (m => n); m = m + 1; end associate), whose value the actual argument does not tell. Or f may give the component one of several different sizes (if (k > 1) then; allocate(r%v(5)); else; allocate(r%v(1)); end if; when every way through f gives it the same size, that size is used). Or two writes that can both run in an iteration give it different sizes (t(i) = f(2); t(i) = f(3)), since a kernel cannot reallocate between them; writes in different branches of one if or select case each give their own size. In all these cases the host does not allocate the component. With the option, with or without bounds checking, the type of t becomes a heap type instead: the value written is made on the device heap and stored into t(i) after the launch, with the size the loop gave it (see Allocating on the device). Without the option, with bounds checking on, the launch stops if the component is not allocated, with the usual "Array ... is not allocated"; a wrong preallocated size is not detected, and without bounds checking nothing is checked and the kernel writes storage that does not exist.
  - Which elements the loop writes. The element may be picked by a call (t(g(i))) or by a value the loop writes. A test may call a procedure or read a value the loop writes. The write may be inside another construct (a loop, a where), or in a loop with a cycle, exit or return that can end an iteration early, or a pointer association. Then no allocated component is changed and nothing is checked, for any write of that component. With the option, when every write of the component gives it the same extents, and the host can evaluate them before the loop runs and they do not depend on the element, every element whose component is not allocated is allocated with them, so an element the loop does not write can end up allocated. Without the option, a written component that is not allocated stays unallocated, even with bounds checking on.
- A do concurrent mask (do concurrent (i = 1:n, mask(i))) is currently ignored on every backend (#12778), and so is a cycle in an offloaded loop (#12856), so such a loop writes elements it should not.
A kernel cannot end the program, so a stop or error stop the kernel reaches, in the loop body or in a procedure it calls, is split between the device and the host (see the next section).

Stopping on the device¶

Every launch of a kernel that reaches a stop or error stop is given a small status buffer, an extra argument of the kernel, which the host clears once, on entry to the procedure that launches it. The offload pass numbers each such statement (its site) and rewrites it, in ASR, before any backend sees it:

In the kernel, the statement records its stop code, when it has an integer one, in the element of the buffer that belongs to its site, then its site in the first element, and returns: that thread runs nothing more. The other threads carry on; the program ends after the launch anyway.
A procedure the kernel calls that holds such a statement, or calls one that does, gets two hidden scalar arguments, the site and the stop code. The statement sets them and returns, and every call to such a procedure is followed by a test that returns from the caller in turn when they are set. The kernel records them in the buffer as if it had reached the statement itself. A call in the condition of an if, do while, select case or in the limits of a do loop is tested at the start of each body as well, so no body runs with the value a stopped call returned. (The rest of the statement that made the call does run with that value.)
When the host next waits for the launch (see Asynchronous launches and data residency) it reads the buffer, and if a site is recorded it runs that same statement on the host, with the recorded stop code. The program therefore prints what it prints on the CPU (ERROR STOP 42, ERROR STOP evaluate: unknown activation) and ends with the same status.
A thread records its stop only while no site is recorded: a site an earlier launch recorded, which the host has not read yet, is where the program ended, and the first launch that stopped is the one reported.

The code is written before the site, and only by a thread that reached that site, so whichever thread's site the host reads, the code next to it was written by a thread that reached the same site: no atomic operations are needed, and no backend has anything to lower but assignments and return. When several threads stop, which one is reported is not specified, as the iterations of a do concurrent have no order.

A check that no call can fail is removed before any of this, so it costs the kernel nothing and does not make it a kernel that can stop (gpu_device_guard_fold.cpp). Each kernel has copies of the procedures it calls, so all the calls of a copy are known: when every call passes an intent(in) integer dummy a value in a range (modulo(k - 1, 4) is in 0..3), the tests of the copy are judged with the dummy in that range, and an if whose test the ranges of the integers it compares decide is replaced by the branch it takes. The error stop of

pure elemental integer function byte_of(word, part)
    integer, intent(in) :: word, part
    if (part < 0 .or. part > 3) error stop "part must be in 0..3"
    byte_of = ibits(word, 8*part, 8)
end function

goes when every call is byte_of(w((k-1)/4+1), modulo(k-1, 4)), and is kept, and stops the program, when a call can pass 4. The ranges come from constants, modulo, mod, iand, ibits, min, max, integer arithmetic and the indices of counted loops.

A character stop code has to be a constant, which the host statement repeats. One computed at run time (error stop "bad " // name) is an error for now ("a character stop code that is not a constant"): it would have to be formatted on the device or copied back. Run-time checks in the kernel (bounds checking) would use the same buffer once the device code checks them; it does not yet.

Allocating on the device¶

Every allocation of device code -- an allocate, the reallocation of an assignment to an allocatable, a temporary the compiler introduces, a run-time sized local, the components of a derived-type value a procedure creates -- goes to exactly one of two places, decided in ASR by the offload passes, never by a code generator:

A buffer the host sizes before the launch (the fast path): the per-thread workspaces, the storage of a derived-type local's components, and the host preallocation of a component of a kernel argument the loop writes, as described above. It is used whenever the host can work out the size before the launch by those rules.
The device heap otherwise: memory the kernel takes at run time from arenas the host allocates for the launch.

Which one is used is decided per array and per derived type:

An allocatable array of device code whose extents cannot be written over the kernel's parameters and its thread index, as a workspace needs, is a heap array: its allocate reads a value the device computes, such as a count the iteration works out (allocate(w(count(x(:, i) > 0.5)))). It is allocated where the allocate is. An array whose extents read the shape of a heap array or of a component of a heap value -- an explicit-shape local or function result such as real :: r(size(self%values_)), or the temporary the compiler makes for such a result -- is a heap array too, allocated where its scope starts (the temporary: just before the call). An assignment of a whole heap array to another one that is not used afterwards (as the compiler's temporaries are) moves its descriptor instead of copying the elements. (Any other array of an explicit shape the host cannot size is an error, and so is a heap array passed for an allocatable or pointer dummy argument, which the procedure would allocate itself: "the array passed for the allocatable dummy argument 'r' of 'f', on the device heap is not supported yet".)
A derived type with an allocatable array component, directly or in a component, is a heap type of a kernel when its device code creates a value of it that host-sized storage cannot hold: an array of such values (the result of an elemental constructor, pair_t(a(:), b(:))), a value held in a component or in an element of an allocatable component of another value, or a local whose component's shape cannot be told on entry to its procedure -- or when, with --realloc-lhs-arrays, the loop gives a component of an element of an array argument a size the host cannot tell before the launch (see the host preallocation above). A type with a component of a heap type is a heap type too. Every value of a heap type that device code creates -- a local, a temporary, a function result, a dummy argument of a device procedure, an element of such an array -- is then on the heap, also one the host could have sized, so that a type has one representation in device code (apart from the kernel arguments, see below). Every kernel has its own copy of each derived type its device code holds, which is what the decision is about; the host's type of an array argument is the kernel's copy the launch passes it for. Two types of one name, from two modules, are two types.

Scalar reductions used as the sole actual of a side-effect-free scalar function, such as rounded(sum(a*b)), can use the same elementwise accumulator as a direct s = sum(a*b): the shared array-temporary pass assigns the typed reduction to a scalar at the original evaluation site. The product then needs neither a workspace nor a device-heap allocation. This normalization excludes array-valued results, calls within the reduction, volatile/asynchronous reads, competing function evaluations and where bodies; sections must already support in-place elementwise reads. Other arrays and workspaces are unaffected.

The arenas¶

A launch whose kernel allocates on the heap gets one arena per element type the heap holds -- each intrinsic type, and each image type (below) -- and a small integer(4) heap state buffer holding the next free element of every arena. They are arguments the offload passes add to the kernel in ASR, and a device procedure that allocates, or calls one that does, receives them as hidden dummy arguments, as it receives the hidden arguments of a stop. An arena is a buffer only the device holds (a scratch buffer of the GpuKernelLaunch in ASR): the runtime allocates it on the device for the launch, and the host neither allocates memory for it nor copies anything in. Only a launch whose kernel stores heap values into an argument, which the host copies out of the arenas after the launch (see below), reads the arenas back, and then only the elements the kernel took from each, the leading part the heap state counts.

Allocating n elements of an arena is a relaxed atomic add of n to the arena's next free element, which returns the old value (GpuAtomicFetchAdd in ASR: atomic_fetch_add_explicit on Metal, atomicAdd on CUDA). The elements are the n that follow it. That is the only atomic operation, and the only operation a code generator has to lower for the heap: everything else is ordinary ASR.
When the arena does not have n elements left, the allocation fails: the procedure runs an error stop with the message "the device heap is exhausted ...", which is lowered like every other error stop device code reaches (see the previous section). The thread stops, and the host ends the program after the launch with that message.
deallocate -- explicit, at the end of a scope, of an intent(out) argument, or before a reallocation -- only marks the array unallocated. The memory is reclaimed when the launch ends, so the heap has to hold everything one launch allocates, over all its iterations.
The size of the heap, in bytes, is read from the environment variable LFORTRAN_GPU_HEAP_SIZE at every launch that uses it (a number, with an optional K, M or G suffix), and is 256 MiB when it is not set. It is split evenly between the arenas of the kernel. Positions in an arena are integer(4), so an arena holds at most 2^31-1 elements.

Heap arrays and heap values¶

A heap array is described by a descriptor: the position in its arena of the element before its first (-1 when the array is not allocated), and a lower bound and an extent per dimension. Its elements are in array element order. The descriptor of a local array of device code is a set of integer(4) locals.
For every heap type the offload passes declare an image type. Each allocatable array component becomes its descriptor (members of the image type), a component of a heap type becomes a component of its image type, and every other component is kept. Every value of a heap type that device code creates has the image type. The elements of an allocatable component whose type is a heap type are images, in the arena of their image type.
A reference to a heap array is written over its descriptor: an element is an element of the arena, at the position worked out from the subscripts, the bounds and the extents; size, lbound and ubound read the descriptor, allocated tests its position. An array passed whole to a procedure is passed as the section of the arena that holds its elements, which sequence association hands to the explicit-shape dummy the shared passes give every array dummy of device code.
An assignment to an allocatable array on the heap (with --realloc-lhs-arrays, or an allocate) allocates it again when its shape changes. An assignment of a whole value of a heap type, which copies its allocatable components in Fortran, allocates them again and copies their elements, recursively. A new image -- a local, an intent(out) dummy on entry, an element of a newly allocated array of images -- starts with every allocatable component unallocated.
A value of an image type is passed to a procedure as any derived-type value is: by reference to a dummy the procedure defines, by value otherwise. A type-bound procedure called on an image is called directly, the binding of its declared type: a heap type is never polymorphic in device code.

Heap values and kernel arguments¶

A kernel argument of a heap type, or an element of one, keeps the host's representation in the kernel: its components are in the device buffers the host hands over, and it is never on the heap. Where it meets heap values:

Used where a heap value is needed -- the actual argument of a dummy of a device procedure, the value assigned to a heap value -- it is copied into a new image on the heap first, its components read from those buffers.
Assigned a heap value -- by an assignment, or through the actual argument of a dummy the procedure defines, which is then a local image of the caller copied in before the call and out after it -- it is not written in the kernel. The kernel stores the image in a copy-out record: an array of images with one element per element of the kernel argument, which the host allocates before the launch, and a flag per element. After the launch, and after a recorded stop would have ended the program, the host assigns every flagged element its image, as a Fortran assignment of the whole value does whatever --realloc-lhs-arrays says: each allocatable component of the element is deallocated, allocated again with the bounds of the image's descriptor, and its elements copied from the arena, recursively through components of heap types, by a host procedure the offload passes generate for every heap type. A component that the image does not have allocated is left deallocated.
Read after the iteration stored an image into it -- only that iteration can, in a do concurrent -- an element of such an argument is that image, which the record holds; before, it is copied in from the buffers the host handed over. (outs(i) = f(ins(i)) after ins(i) = g(...) in one iteration reads the new ins(i).)

Asynchronous launches and data residency¶

A launch does not wait for its kernel, and the data a kernel writes stays on the device until the host needs it. The host waits only where a host statement needs what a kernel wrote, so a host loop that launches kernels over and over (a training loop) runs without waiting for the device in every iteration. The decisions are made in ASR, by the passes below; the runtimes only act on what the compiler tells them.

What the runtime does¶

Every buffer argument of a launch is handed to the runtime with how the kernel uses it, which gpu_kernel_finalize records in the kernel's layout (gpu_kernel_argument.written): a buffer is read-only when every reference of the kernel to it only reads it -- an operand, a subscript, the value of an assignment, an inquiry about its shape, an actual argument of an intent(in) dummy, the target of a pointer the kernel binds (for the section an elemental function is called with) that nothing writes through -- and written otherwise.

The host bytes of a read-only buffer are copied when the launch is made. The CUDA runtime keeps the copy of a range of at least 64 KiB in device memory until the next lfortran_gpu_wait or lfortran_gpu_sync, and a later launch that reads the same range uses it without copying again: the host does not change what a launch since the last wait reads before it waits (the gpu_launch_sync pass places the waits). The launches of a time loop that only read an array, between which the host touches nothing they use, copy it once. On an MI210, where a copy runs at about 25 GB/s, that makes a matrix-vector product 27 times faster (24 to 634 GB/s) and a sparse one 26 times. A wait right after a GpuKeep keeps the copies of the arrays it names (see Where the host waits), so a time loop whose host reads the result of a reduction at every iteration copies the arrays it reads once too. A kept copy that no launch used since the wait before the last one is freed.
A written buffer is resident: the runtime keeps a device buffer for its host range in a table. When the device does not hold newer values than the host (the entry is not dirty), the host bytes are copied in; after the launch the entry is dirty. A later launch handed the same range -- to read or to write -- uses the same device buffer, without copying anything.
A written buffer whose launch records a whole-array overwrite before read can use LFORTRAN_GPU_WRITE_DISCARD when its host coverage condition holds. CUDA and Metal skip copying it in unless another buffer argument overlaps it, in which case it remains an ordinary write. This does not establish gpu_kernel_argument.noalias; that proof is independent. The accumulators of a reduction are always recorded as overwritten: the kernel sets each one to the identity before it combines into it. Strict integer guards are narrowed without overflowing their endpoints, including dynamic integer(8) bounds. A strict guard requiring a value below the smallest integer or above the largest cannot establish an overwrite and retains ordinary write access.
lfortran_gpu_sync waits for every launch, copies every dirty entry back to the host, and marks it clean.
lfortran_gpu_discard forgets the dirty entries of a host range that is about to be freed, without copying anything back; in the unified and managed memory models it first waits for a pending launch that uses the range in place.
A range that overlaps a dirty one without being the same range makes the runtime wait first, so that the bytes it copies are current. The Metal runtime does not wait for a range a kernel only reads that holds the dirty ones whole (and that no range the kernel writes overlaps): it copies the host bytes and then, on the device, the dirty ranges over them (a blit). A program that writes the newest positions of a cache and then reads all of it, as attention over a KV cache does, so keeps the launches going.
The runtime never writes a device buffer that a pending launch may still use: a copy into a resident buffer that is still in use goes into a fresh buffer, and pooled buffers are reused once the launches that used them have completed. A scratch buffer -- a run-time sized local of the kernel (a workspace), an arena of the device heap -- is only on the device.

The CPU emulation (--gpu=cuda_cpu) runs this same CUDA runtime (lfortran_gpu_cuda.c), on an emulation of the CUDA runtime API that it uses (src/libasr/runtime/cuda_cpu/), so it needs no GPU and no CUDA toolkit; the host compiler builds everything. The emulated device memory is separate from host memory: cudaMalloc makes an allocation of its own, filled with 0xff bytes, and only the copies move data, so a missing or late copy shows as it does on a device. In the emulation the runtime also refills every unaliased discard buffer with 0xff bytes (cudaMemsetAsync) before each launch, so a wrong overwrite proof cannot pass on what an earlier launch left in a reused buffer; natively the buffer keeps those bytes. The stream runs every copy and launch to completion before the call that issues it returns, in the order they are issued, which a program that waits for its results before reading them cannot tell apart from a device. A kernel is reached through a thunk that the CUDA code generator emits for it in emulation and registers in its place; the emulated cudaLaunchKernel runs the blocks one after another, and the threads of a block one after another, so a __syncthreads() cannot be honoured and stops the program. The emulated device reports a 64 MiB L2 cache. The runtime builds against the CUDA runtime API of CUDA 12 and of CUDA 13, which differ where the runtime prefetches managed memory (cudaMemPrefetchAsync takes a cudaMemLocation in CUDA 13). The emulation provides the CUDA 13 API, so --gpu=cuda_cpu builds the runtime as a CUDA 13 toolkit does; it provides the CUDA 12 API when CUDART_VERSION is defined below 13000, which the test_gpu_cuda_api_state_cuda12* unit tests do. The CUDA runtime keeps the same table and orders copies and launches on one stream. It keeps scratch buffers, and upload buffers for read-only ranges below 64 KiB, per kernel and argument slot, reusing their capacity and deferring replaced allocations until a synchronization. Preparing arguments can itself synchronize when the deferred queue or residency table fills. After such a synchronization, CUDA prepares all bindings again before dispatch: writable residents must still be present and dirty for the new launch, and reads must use their current bindings. This does not add a synchronization to each launch. Discard buffers participate in this same retry and writable-residency preparation. Same-address writable aliases, including discards downgraded by overlap, use one covering extent so that preparation cannot alternate between their different byte counts.

--gpu=hip runs the same CUDA runtime on AMD GPUs, through ROCm. HIP is CUDA C++ for AMD GPUs, so the backend shares everything with --gpu=cuda up to the toolchain: the ASR passes, the device capabilities and the device code generator, whose output only gains #include <hip/hip_runtime.h>. The one capability that differs is the unrolling of short sequential loops in device code: ROCm's device compiler does not unroll a loop whose trip count it does not know, so HIP unrolls them four times (a matrix-vector product with a thread per row ran 3x faster on an MI300A), while CUDA unrolls none. The HIP driver (hipcc by default, --device-compiler to choose another) compiles it as HIP for --gpu-arch (gfx90a, gfx942, ...; the GPUs of the compiling machine by default) and links the program. The GPU runtime is lfortran_gpu_cuda.c itself, which the host C compiler builds against hip/cuda_runtime.h, an implementation of the part of the CUDA runtime API it calls with HIP, and the headers of the ROCm installation the HIP driver belongs to. The program finds the HIP runtime library of that installation through its run path.

Memory models¶

A GPU with memory of its own -- every discrete GPU, NVIDIA's and AMD's -- runs a kernel on copies of the host's arrays in device memory: the copy model above. A GPU that shares the host's memory and can access what the host allocated, such as the AMD MI300A with XNACK on, can also run a kernel on the host's arrays where they are: the unified model. On an MI300A a kernel runs at the same speed either way, but the copies cost: a copy between host and device memory runs at about 58 GB/s there, while the kernels run at about 3.5 TB/s, so a launch that copies its arrays spends nearly all its time copying.

The unified model is a special case of the copy model rather than a design of its own: the compiler emits the same launches, waits and discards for both, and the CUDA runtime keeps the same residency table, in which a range used in place is an entry whose device memory is the host range itself. Only what the runtime does with an entry changes:

A range a kernel may write is used in place: nothing is allocated, copied in or copied back, and lfortran_gpu_sync only waits. The host does not access such a range before the next sync, as in the copy model.
The first touch of a page by the GPU is a page fault that the driver serves, at about 2 GB/s on an MI300A, while glibc serves a large allocation with a fresh mmap and unmaps it when it is freed: a temporary array a loop allocates at every iteration, or the accumulators of a reduction, would be fresh memory at every launch (a kernel writing a fresh 512 MB array took 0.2-0.4 s instead of 0.2 ms). The runtime therefore makes glibc serve every allocation from its heap and never trim it (mallopt(M_MMAP_MAX, 0) and mallopt(M_TRIM_THRESHOLD, -1)), so freed memory stays mapped, for the host and the GPU, and is reused. LFORTRAN_GPU_KEEP_MEMORY=0 leaves malloc as it is.
A range a kernel only reads of at most 64 KiB is still copied; a larger one is used in place. The host does not change or free it while the launch is pending: the gpu_launch_sync pass puts a GpuWait before a host statement that may (see "Where the host waits").
lfortran_gpu_discard of a range that a pending launch uses in place waits for the launch, since the host is about to free memory the kernel uses.
The runtime gives the driver no advice about a range it uses in place. On an MI300A, hipMemAdvise(hipMemAdviseSetPreferredLocation) for the GPU made such kernels erratic: a STREAM triad ran at 2.7 TB/s in one run and at 3 GB/s in the next, and a matrix-vector product at 1 GB/s in most runs. Without it the triad runs at 3.2 TB/s and the product at 330 GB/s in every run.

A GPU with memory of its own whose driver can share memory with the host page by page -- CUDA managed memory on an NVIDIA GPU with concurrent managed access, such as an H100 -- has a third model, the managed model, which is what nvfortran's -stdpar=gpu uses for allocatable arrays. In the copy model every launch copies the arrays it only reads to the device again (the host may have changed them), and a host access copies back every array a kernel wrote, even to read one element: on an H100 over PCIe a STREAM triad of 1 GB arrays runs at 13 GB/s that way, against 1.9 TB/s in the kernel. In the managed model:

The CUDA runtime replaces the program's allocator (the default allocator of lfortran_intrinsics.c) when the program starts, so that an allocation of at least 1 MiB -- an allocatable array -- is managed memory. Smaller allocations (strings, temporaries), and memory the allocator does not make (stack and static arrays), stay as they are, in the copy model.
A range inside a managed allocation is used in place and has no entry: the driver moves its pages to the device when a kernel touches them and back when the host does, one page at a time. Nothing is copied at a launch, lfortran_gpu_sync only waits, and a host read of one element moves one page.
The first launch that uses an allocation moves all of it to the device (cudaMemPrefetchAsync): 5 times faster than the page faults of the kernel when the host has just filled it. Later launches leave the pages where they are; prefetching ranges that are on the device already costs more than the faults it saves.
As in the unified model, the host does not change a range a pending launch reads before a GpuWait; freeing managed memory waits for every launch.

The runtime selects the model at the first allocation it could make managed, or when it starts. LFORTRAN_GPU_MEMORY=copy, unified or managed selects one; auto, the default, selects the unified model on a GPU that reports that it is integrated with the host and can access its pageable memory (on AMD, an APU with XNACK on), the managed model on an NVIDIA GPU with concurrent managed access, and the copy model otherwise. A separate AMD GPU keeps the copy model: with XNACK off a kernel reads the managed pages the host touched from host memory (a STREAM triad ran at 7 GB/s on an MI210), and with it on its speed varied from 5 to 1,380 GB/s between identical runs.

XNACK is what lets an AMD GPU serve its page faults, and so access memory the host allocated. ROCm turns it on for a process run with HSA_XNACK=1, and reads the variable once, when the HIP runtime starts. The HIP runtime of LFortran sets HSA_XNACK=1 before the program starts unless the user set it, so an MI300A runs in the unified model by default and HSA_XNACK=0 keeps XNACK off (and the copy model). A GPU or kernel that cannot serve the page faults keeps XNACK off whatever the variable says. Programs the program starts inherit the variable. unified on a GPU that cannot access host memory is an error. LFORTRAN_GPU_VERBOSE=1 prints the model selected. The CPU emulation reports a separate device that can access host memory but not managed memory concurrently, so its default is the copy model, and LFORTRAN_GPU_MEMORY= unified or managed runs the other models on the CPU. In the unified model on AMD the host waits for the GPU by sleeping until the driver wakes it (hipDeviceScheduleBlockingSync) rather than by spinning: the CPU and the GPU of an APU share the power of their socket, and on an MI300A a host thread spinning while a kernel streams through memory made the kernel 3-5% slower.

An Apple GPU shares the host's memory, but a kernel can only use memory a Metal buffer covers, which must be whole pages; making a buffer over a range the host allocated (newBufferWithBytesNoCopy) costs about 0.1 ms, while the first launch that uses a new buffer maps its pages for the GPU (about 10 ms for 1.4 GB). So the Metal runtime has the managed model, chosen as on CUDA (LFORTRAN_GPU_MEMORY=copy or managed; auto selects managed on a GPU that shares the host's memory, every Apple GPU, and copy otherwise):

The runtime replaces the program's allocator, so that an allocation of at least 1 MiB takes whole pages. The first launch that uses it makes a buffer over its pages, kept until the allocation is freed, and a range in it is handed to the kernel in place, at its offset in that buffer: nothing is copied in at a launch or back at a sync. Freeing the allocation waits for the launches that use it, and forgets the buffer, so an allocation made later at the same address gets a buffer of its own.
A range a launch only reads of at most 64 KiB is still copied, unless a pending launch writes the allocation in place, and so is a range that overlaps one the same launch writes, so that the kernel reads it as the host holds it; the runtime first waits for pending launches that write the allocation.
lfortran_gpu_wait waits for every launch when one that reads host memory in place may be pending, and returns at once otherwise. lfortran_gpu_discard of a range in a managed allocation that a pending launch uses waits for the launches; of any other range it waits for none.
Smaller allocations, memory the allocator does not make (stack and static arrays), and an allocation the device cannot map keep the copy model.

Decoding a 27B-parameter language model whose weights are allocatable arrays, the copy model copied 28.8 GB of weights into fresh buffers for every token; the managed model copies none of them.

The Metal runtime encodes the launches one after another into one compute encoder and commits its command buffer every 64 launches and at every wait: encoding a dispatch costs the host about 2 us, committing a command buffer several times that, and each command buffer has an overhead of its own on the device. The encoder runs its dispatches concurrently, with a barrier before a dispatch that uses a device buffer that a dispatch since the last barrier writes, or writes one that such a dispatch reads: launches that do not depend on each other overlap, and the others run in launch order. At most three committed command buffers are pending; a host that launches faster than the device runs waits for the oldest, which also bounds how many buffers the pending launches hold. A scratch buffer is kept per kernel and argument, grown when a launch needs more, and used by every launch of the kernel, which writes it, so launches of one kernel that share it never overlap. That holds for a large one too, such as a device-heap arena or the workspaces of a kernel with many threads: a new buffer's pages are mapped for the GPU at its first use, and a projection of a language model that got new workspaces of 60-360 MB at every launch spent more time on that than in its kernel. A launch that runs in chunks encodes every chunk with its one kernel handle, so the chunks share all its scratch buffers: what an earlier chunk leaves in a scratch buffer is still there for the later chunks and for the host to read back. Between launches, what a launch leaves in a scratch buffer is dead.

Where the host waits¶

The gpu_launch_sync pass, which runs just before device_launch_expand, places the GpuSync statements. In each procedure that launches a kernel it follows the statements in order -- through if, select case, block, associate, exit and cycle, and to a fixed point through loops -- and keeps what pending launches may write and which of them may have recorded a stop. It waits:

before a host statement that may read or write, allocate or reallocate storage a pending launch may write;
before a statement that acts outside the program -- input/output, stop, or a call to a procedure that is not pure and whose body may do input/output or stop -- while a launch that may have recorded a stop is pending, so that the stop the kernel recorded ends the program first;
before a launch of a kernel that can stop while a launch of another such kernel is pending, so that the stop reported is the first one;
before return, goto and the end of the procedure while any launch is pending, except in a deferring twin (below);
before any statement it does not know.

A wait in a loop that launches nothing is placed before the loop instead. A call of a subroutine whose first statement is if (c) return, such as a debug hook that returns unless a dump directory is set, tests c first (the gpu_call_guards pass, which runs before the temporaries of actual arguments are made), so its wait goes into the branch that calls it. It does so only when that is not observable: c calls no function and names no dummy or local of the subroutine, the subroutine has no intent(out) dummy, and the actual arguments call no function. After every GpuSync the pass runs the stop replay of each launch that may be pending there (GpuKernelLaunch.stop_replay), in launch order.

A module subroutine with an array or derived-type dummy -- not bind(c), elemental, a separate module procedure or device code, and with no pointer, target, procedure or frame dummy -- gets a deferring twin, <name>__gpu_deferred, which takes over its body, while the subroutine itself only calls the twin and then waits as before. The pass makes the twins first, of the modules of the file and of those loaded from module files alike: the choice depends only on the declaration, which is in the module file, so every file that uses the module makes the same twins. A module compiled with --separate-compilation and without a GPU device (a library built for the host) records that its object defines its procedures and no twins: other files call its procedures themselves. In the default mode the file with the main program compiles the procedures of the modules it uses, with their twins only if it is compiled with --gpu, so it must be compiled with --gpu too when other files that call them are. The calls the pass sees by name are redirected to the twin; a call through a procedure pointer, a dummy procedure or a binding, and one from a file compiled without a GPU device or from C, reaches the subroutine, which waits. The twin may return with launches pending, provided they write and read only what its dummy arguments hold -- not what an intent(in) dummy holds, a module variable, what a pointer points to -- and none may have recorded a stop (whose replay reads a status buffer that is a local of the twin): otherwise it waits at its exits as before, once the launches on its locals are discarded. A call of a twin is followed like a launch: it may leave launches pending that write the actual arguments of the dummies that are not intent(in), and read the others -- or, for a twin of this file, only those that the launches pending at its exits do. A caller does not wait before the call for the launches pending on an array actual argument that only the twin can name (a local or dummy no pointer reaches); the twin starts as if launches its caller made were pending on what its array dummies hold, and waits where it accesses them on the host. Before any other actual argument, a module variable, and when a launch that may have recorded a stop is pending, the caller waits as before.

Storage that a launch only reads the runtime may hand the kernel in place (the unified and managed memory models), and the host must not change it under the kernel; in the copy model the runtime keeps its copy in device memory until the next wait, and the host must not change it before that either. The pass keeps what pending launches read as well, and puts a GpuWait -- wait for the launches that read host memory in place, copying nothing back -- before a host statement that may change or free such storage. Any statement that accesses the storage counts, except where it only reads it: where it evaluates an expression, such as the value of an assignment, or hands it to an intent(in) dummy of a deferring twin that cannot change it through a pointer (a weight matrix handed to every layer's projection while the previous projection still reads it in place). A launch itself never waits, so launches that only share what they read, such as the iterations of a time loop, run back to back.

A wait ends what the copy model keeps of what launches read, unless the host leaves it unchanged until launches read it again: a time loop whose host reads the result of a reduction at every iteration waits there, and the arrays the reduction reads would be copied again at every iteration. So, with every wait placed, the pass follows the statements once more, and puts a GpuKeep before a wait that names the arrays a launch since the previous wait read (or that the previous wait kept) which, on every path from the wait to the next one, no host statement may change, and no launch since the last GpuSync writes. Here a read is told apart from a write: an assignment changes only its target, an expression only what the procedures it calls may change, and a procedure contained in this one changes a variable of this one only if it, or a contained procedure it calls, names it. Only whole contiguous arrays that are locals or dummy arguments of the procedure are named, not an assumed-shape dummy, which may be a section with gaps, nor a pointer. Locals of a block end with it, and a procedure's arrays may change once it returns, so a procedure other than the main program in which a copy is still kept when it returns puts a GpuWait before the return (or at its end), unless the keeping wait is the statement right before. On an MI210 that makes a time loop of dot products whose host reads the sum at every iteration 36 times faster (23 to 830 GB/s).

Storage that a pending launch writes and the host frees is not waited for: a local of a block at the end of the block, or where exit or cycle leaves it, and an array the host deallocates, get a GpuDiscard instead, since nothing can read those values any more. A local array that a pending launch reads gets a GpuDiscard where it ends too, at the end of its block or when a deferring twin returns: the runtime waits there only for a launch that uses that storage in place, not for every launch that reads host memory in place. A twin such as linear(x, w, y), whose second launch reads a local workspace the first one writes while the weights w are read in place, so returns without waiting for its launches.

What a statement accesses is worked out over access paths: a variable and a chain of components (workspace%dcdw), whatever the subscripts. A pointer the procedure associates with one designator only -- an associate name, or a pointer the compiler creates for a section -- stands for it. An array a launch is handed through a local pointer the procedure associates with a few designators, such as the pointer the lowering of an optional argument associates with the argument or with a placeholder, is each of them. Two paths overlap when one is a prefix of the other, except that a copy of only the bytes of a derived-type value itself (a polymorphic value the host copies for a launch) does not reach the data of its allocatable components. Asking about the shape of an array (size, lbound, allocated) does not access its elements, and neither does disassociating a pointer (nullify) or storing the address of a variable in the frame of a contained procedure. A call accesses its actual arguments, and may access every variable that a pointer may point to or that another procedure can reach (a module variable, a variable of a procedure this one is contained in, a target, a variable a pointer that another procedure can name may point to: not one only a local pointer of this procedure points to, such as the pointer the compiler makes for a section it hands to a procedure); a call to a procedure contained in this one may access anything. Of a derived-type value it is handed, a call whose callee's body is known accesses only the components the body names. A call of a procedure whose body is known, launches nothing and has no pointer, target or procedure dummy accesses what that body accesses: the same parts of its actual arguments, the variables of its host procedure it names, what another procedure can reach only when the body names or calls something that can, and it acts outside the program when the body may, stop included. So a call that only asks for the shapes of its arguments (call check_shapes(x, y)) does not wait for the launch that writes them.

The data the host copies for a launch itself -- a contiguous copy of a strided array, the components of an array of derived type, the arrays of a packed launch, the plain components of a polymorphic value -- are host accesses of the launch, and device_launch_expand waits after the launch itself before it copies any of them back.

Whole-array assignments between launches¶

In a host loop that launches kernels, a whole-array assignment on the data the kernels use -- clearing an accumulator (dcdw = 0.), or reducing what a kernel wrote (dcdb = sum(reduce_dcdb, dim=3)) -- would make the host wait for the device and copy the data back in every iteration, only to hand it to the next launch again. Such an assignment is a loop over the elements of its target whose iterations are independent, so the gpu_array_statements pass, which runs right after parallel_dispatch, makes it a parallel loop and sends it to the device too:

dcdb = sum(reduce_dcdb, dim=3)
! becomes
do concurrent (i1 = 1:size(dcdb, 1), i2 = 1:size(dcdb, 2))
    dcdb(lbound(dcdb, 1) + i1 - 1, lbound(dcdb, 2) + i2 - 1) = &
        sum(reduce_dcdb(lbound(reduce_dcdb, 1) + i1 - 1, &
                        lbound(reduce_dcdb, 2) + i2 - 1, :))
end do

An assignment is sent to the device when all of this holds:

It is in a sequential host loop (do, do while) whose body has, at any depth, a parallel loop the dispatch sent to the device, and its target or an array it reads is one that such a loop names, or one that another assignment sent to the device names.
Its target is a whole array of an intrinsic type -- a variable, a component, or an associate name bound to a whole array -- and the assignment cannot change its shape: not an allocatable it may reallocate (with --realloc-lhs-arrays), not a pointer other than an associate name, and not a defined assignment.
Its value is computed element by element: scalars, arrays of the target's rank, intrinsic operators, casts and elemental intrinsic functions of them, and sum or product with a constant dim of an array of one more dimension. A scalar is evaluated by every iteration instead of once, so it may not call a function. The value may read the target only element for element (x = 2*x); a value that may read other elements of it (x = x(n:1:-1), x = x(1)) keeps the assignment on the host.
The loop passes the unsupported-construct check: an assignment the device cannot run stays on the host, without an error or a warning that the source did not have.

The host still evaluates the bounds of the loop (size, lbound), which read only the shape of the arrays.

An assignment whose every element is assigned again before anything reads it (dcdw = 0. followed by dcdw = sum(reduce_dcdw, dim=4)) is removed before this, by the dead_array_store pass, so it costs no launch at all.

Fast math in device code¶

With --fast (and without --no-fast-math) the CUDA and HIP device code computes 1/sqrt(x) as rsqrt(x), one operation where a correctly rounded quotient of a correctly rounded root takes two long ones: an n-body force loop ran 16% faster on an H100. The result can differ from the exact one in the last bit, as host code under --fast can.

Packed data: unrolling by the divisor¶

A kernel that decodes packed data reads, for each index d, a word w((d - 1)/4 + 1) and a part modulo(d - 1, 4) of it, often through a small pure function that checks the part, as the decode of a quantized matrix-vector product does:

do d = (g - 1)*64 + 1, g*64
    acc = acc + real(packed_byte(w((d - 1)/4 + 1, n), modulo(d - 1, 4)))*x(d)
end do

The gpu_residue_unroll pass unrolls such a loop (a do or a one-index do concurrent of device code) by the constant divisor, in groups that start where the numerator d - 1 is a multiple of it. In each of the four copies the part is a constant and the quotient the same, so the word is read once per group, and a call of a small pure function with a constant argument is inlined so that its check folds away: no reloaded word, no modulo, no call and no test for a stop remain in the loop. A check that still fails, such as packed_byte(w, 4), still stops the program. The iterations run in their original order, so the results are the same bits.

The groups need the numerator to be non-negative and at least four iterations; what the bounds of the loops around it do not prove is tested before the loop, which runs as written otherwise, and the iterations before the first group and after the last run as written too. On an Apple M4 Max, a 5120x17408 q8 product whose kernel runs the loop above took 1.02 ms of GPU time instead of 2.73 ms.

Integer powers in device code¶

A real power with a constant integer exponent from 0 to 8 in device code, such as x**2, is written as the products the host computes it with (gpu_integer_powers): x*x, x*(x*x), (x*x)*(x*x) and so on, with x evaluated once, into a variable of its own unless it is one already. The device's pow takes a logarithm and an exponential under fast math, which is slower and not always exact (a product of powers of two came out inexact on an Apple M1), where x*x is the correctly rounded square the host computes. A power whose base can only be evaluated where it is written (in the test of a do while, or an operand of .and. that need not be evaluated) is kept as it is, unless its base is a variable or an element.

Lanes: an iteration on a SIMD-group¶

A loop with few iterations -- a training step over a mini-batch of one, a loop over the layers of a network -- leaves most of the GPU idle, and runs each iteration as one long chain of dependent memory accesses on one thread: a GPU thread does not run ahead of a load the way a CPU core does, so every element of every array statement inside the iteration costs a memory latency. The gpu_vector_lanes pass, right after gpu_kernel_finalize, rewrites such a kernel so that each iteration runs on the lanes of a SIMD-group, which share the loops inside it (32 lanes on Metal and on an NVIDIA GPU; a device's capabilities say how many, see GpuDeviceCapabilities::lanes_per_iteration, and which loops they share):

A loop whose iterations are independent -- it assigns every shared array at one subscript list holding its variable, reads such an array only at that subscript list (or at one provably different element: a(i, l) and a(k, l - 1)), and every scalar it assigns is assigned first in each iteration and dead after it -- is shared out: lane k takes iterations k, k + lanes, ..., and the lanes wait for each other after it (GpuLaneBarrier). Adjacent loops over the same iterations are joined first, and a loop whose body is one more such loop shares out the pairs of iterations of both.
A loop that accumulates into private scalars -- the inner product of a matrix-vector product, the forces on one body -- by statements s = s op e (+, -, *, max, min, one operation per scalar), and otherwise only assigns scalars of its own in every iteration and integers it adds constants to (the second index of an array statement whose sections start at different subscripts), in statements, blocks or loops of its body, is shared out as well where the lanes can combine values (lane_reductions: an NVIDIA GPU, whose warps run them, and Metal, for values of 4 bytes): the lanes but the first start every accumulator at the identity, and after the loop every lane gets them combined (GpuLaneReduce: shuffles within the warp on CUDA, simd_sum, simd_product, simd_min and simd_max on Metal). A loop that reads a matrix along a row (a(i, j) over j) is left to one lane per iteration, since the threads of the launch already read a column of it together, and so are a loop of fewer than 64 iterations known at compile time and a loop whose bounds are read from an array, as the rows of a compressed sparse matrix are, which usually hold a few elements. An independent loop that holds such a loop of a length not known at compile time has that sum shared out rather than its own iterations, which may be fewer than the lanes (a long sum split in a few parts). A launch with one lane per iteration runs the loop as it was, in the other branch of a test of the lanes. On an H100 the forces of an n-body loop of 32768 bodies run 1.5 times faster this way and the column sums of a 12288 by 12288 matrix 2.6 times faster.
A statement with an effect the lanes must not repeat -- one that writes shared memory outside such a loop, calls a procedure that is not pure, or takes elements of the device heap -- runs on the first lane, and the others take the scalars it assigned from it (GpuLaneBroadcast).
Everything else runs on every lane, with the same values.

The pass sees the loops of the kernel's own body, not those of the procedures it calls, so a call in a kernel of a small pure function that runs a loop, such as sum_squares(x(:, m)) returning sum(x**2), is inlined first (gpu_inline_loop_calls, using the inliner of the loop transformations, GpuInliner): its loop is then one of the kernel's, and the lanes share it and combine its sum. An RMS norm of a 5120-element row computes its sum of squares this way.

A loop shared out between the lanes steps by the lanes, so each lane runs a loop of its own with one load after another. On Metal such a loop is unrolled sixteen times like a sequential loop of device code (gpu_sequential_unroll, its copies at group + k*lanes), so that each lane has sixteen independent loads in flight: the lanes of a 5120-element sum of squares then wait for memory once per sixteen elements, not once per element.

A lane that stops in a shared-out loop records the stop and returns, as a thread does, and the other lanes return after the loop. A kernel is left as it is when a statement in the way cannot be run on one lane (it returns, or assigns a private array), and when it writes an element of a kernel argument that every iteration writes -- an array declared outside the loop that each iteration uses as a temporary, which the launch shares between all threads (a known limitation, which lanes would only make show more often).

The host chooses the lanes per launch: the SIMD-group width for a launch of at most lanes_iterations_max iterations (1024 on Metal, 65536 on an NVIDIA GPU, which keeps about 230,000 threads resident), one otherwise, where the iterations fill the GPU anyway and the lanes would only repeat each other's work. A kernel whose shared-out loops all read along the first subscript of the arrays they name, at least one along a column the iteration takes (a later subscript depends on it), and that repeats no such column walk on every lane, gets lanes up to lanes_column_iterations_max iterations (huge(0)/32 on Metal): one lane per iteration would read 32 columns at once there, a few elements at a time, while the lanes read one column together. Every row of a quantized matrix-vector product, and every row of a gated-delta state update, is such an iteration. On an NVIDIA GPU only reduction loops are shared out (lane_independent_loops is off), and a kernel with workspaces or that takes memory from the device heap keeps one lane per iteration. The kernel takes the lanes as its last parameter (gpu_kernel_layout.lanes); its thread index then numbers the iterations, the host starts lanes threads per iteration, and a workspace slice belongs to an iteration, shared by its lanes.

Deliberate relaxations¶

Two relaxations are deliberate. A pure procedure the host calls while a launch that may have recorded a stop is pending is not waited for, although it may itself error stop: when both stop, the program reports the procedure's stop. And an error of a kernel (a fault, the watchdog) is reported where the host next waits, and not at all when the program ends with a stop before that.

Profiling the Metal runtime¶

The Metal runtime reports where a program's GPU time goes when the environment variable LFORTRAN_GPU_PROFILE is set:

LFORTRAN_GPU_PROFILE=1 counts dispatches, committed command buffers, host syncs and waits, the bytes copied to and from the device, the device buffers allocated, the time spent compiling kernels, and the GPU time of the command buffers (their sum, and their union). For each kernel it counts the launches, the bytes they copy to the device, and the host time they take in lfortran_gpu_launch until they are encoded: the copies of their arguments, the buffers allocated for them, and any sync they need.
LFORTRAN_GPU_PROFILE=2 also times each kernel on the GPU. Every launch is then committed in a command buffer of its own, and the GPU time of that command buffer (GPUEndTime - GPUStartTime) is the kernel's. That runs the launches one after another, without the overlap a concurrent encoder allows, and adds a commit to each launch; for a kernel that takes a few microseconds the time is an upper bound.

At exit, or at lfortran_gpu_shutdown, the runtime waits for the launches still pending and prints the totals and a table of the kernels, sorted by GPU time, to stderr, or to the file LFORTRAN_GPU_PROFILE_FILE names. A kernel is named by its entry point, __lfortran_gpu_kernel_<n>_<scopes>: the number of the kernel in its file and the names of the program units and procedures the loop is in, each prefixed by its length (for example __lfortran_gpu_kernel_0_6solver6update for the first kernel of the file, in the procedure update of the module solver), followed by _<id> with --separate-compilation. A procedure with a binding label, bind(c, name="<label>"), is named by its label, prefixed by b and its length (__lfortran_gpu_kernel_0_b8solve_c1 for name="solve_c1"), since only the label of such a procedure is global. Scopes whose names together are longer than 128 characters are replaced by h and a hash of them, __lfortran_gpu_kernel_<n>_h<digits>. The totals are also available to C code through lfortran_gpu_metal_telemetry_snapshot (lfortran_gpu_runtime.h). To get the cost of one step of a program, such as one generated token, run it for two numbers of steps and subtract the reports.

LFORTRAN_GPU_PROFILE=2 LFORTRAN_GPU_PROFILE_FILE=profile.txt ./a.out

Testing¶

Every entry of the unsupported list has a test for the error, and a test that the same program runs correctly on the CPU with --gpu-allow-cpu-fallback.
Every pipeline fix comes with an integration test that runs on the GPU backends (metal, cuda labels) and checks the results.
One cuda label covers every mode of the CUDA device code, and every cuda test must pass in all of them: integration_tests -b cuda_cpu runs the cuda tests in the CPU emulation (--gpu=cuda_cpu, what CI runs), integration_tests -b cuda runs the same tests natively on an NVIDIA GPU (--gpu=cuda), and integration_tests -b hip on an AMD GPU (--gpu=hip). There are no cuda_cpu and hip labels. The memory models are tested by running a backend with LFORTRAN_GPU_MEMORY=copy and unified: the CPU emulation in either, an AMD APU with XNACK on in either, and a discrete GPU in the copy model.
Compiling is not enough: an offloaded loop can compile and still compute the wrong values. Integration tests therefore check their results in the program, and third-party codes whose test suites check values (Formal, Fiats) are run with --gpu=metal in CI.

The test_gpu_cuda_api_state* unit tests compile the production CUDA runtime against the CPU emulation of the CUDA runtime API that --gpu=cuda_cpu uses, with kernels of their own, and check how the runtime uses the API through the statistics the emulation keeps (cuda_cpu/cuda_cpu_emulator.h). They exercise deferred-queue saturation for scratch and read-only buffers, and residency-table eviction, and check readback, allocation lifetime, and that no memory is freed while a copy or launch issued since the last synchronization may still use it. Discard variants exercise both flush paths without uploading the output. Alias cases also check that a discard shares the covering writable resident, falls back to one upload, and preserves values read before an overlapping write. These are API-state tests, not native GPU execution or performance evidence:

pixi run -e llvm11 ctest -R '^test_gpu_cuda_api_state' -V > cuda-api-state.log 2>&1

The same sequences also run against the real CUDA runtime in src/tests/gpu_cuda_scratch_order.cu, which calls the CUDA API directly and is therefore not an integration test. It needs NVCC and an NVIDIA device: 64 kernel caches with 30 small scratch buffers reach the 4096-entry deferred queue boundary before dispatch 201, whose output and scratch readback must both be 201. The script builds it and runs the scratch sequence, its variant with a discarded output, and the mixed write/discard alias cases. Run it on a GPU machine even when the host API-state tests pass:

src/tests/gpu_cuda_scratch_order.sh > cuda-scratch-order.log 2>&1

gpu_metal_507 checks a full overwrite through a target array that reads its old elements through a pointer alias, on the ordinary, CPU-emulated, Metal and native CUDA paths. The CUDA API-state tests do not replace any of these backend tests.

gpu_discard_guard_extremes checks strict and inclusive guards at both integer(8) extrema, constant and dynamic bounds, reversed comparisons, conjoined guards, and narrower integer(4) bounds. Its CUDA_CPU checks reject poisoned NaNs as well as finite wrong values; it also runs with LLVM, GFortran, and native CUDA.

Open questions¶

Purity is not enforced yet. LFortran currently accepts an impure procedure (for example one that prints) called from a do concurrent loop, and stop in its body. The check relies on the standard's constraints, so these need to become semantic errors first.
What counts as using a type. The check uses the rule above. The pipeline does not always follow it yet: a derived type whose components are read only through allocatable array components is handed to the kernel one component at a time, but one with other components read, or used whole (passed to a procedure), is handed over whole, its allocatable components included, and the loop reads those through it, so an unused real(8) component there is still an error on Metal (reported by the pipeline, not the check). Likewise an unused component of a kind the device does not interpret -- a pointer, or a character that is not allocatable -- still keeps a value handed over whole off the device, while an unused allocatable scalar component no longer does: the device keeps the host's bytes for it.
Procedures without a visible body. A bind(c) interface is how a kernel calls a function the device provides (for example sqrt), but a host-only C routine looks the same. Either keep an allowlist of device-provided functions, or put every other call to a procedure with no body on the list.
Separate compilation. A called procedure can live in another file or submodule. The check reads the body of a separate module procedure where the pipeline does -- the submodule of the translation unit, or the module file of a submodule compiled separately -- so it looks at the body the kernel would be built from. Reading a module file at every compilation that calls into it costs time; recording the facts the check needs (which types a procedure uses, whether it does input/output) in the module file instead would avoid that.
Waiting at the end of a procedure. A procedure that launches kernels waits for them before it returns, so a procedure that launches a kernel and is called in a loop waits once per call. Keeping its launches pending across the return needs what the callers do with the data, which this analysis, one procedure at a time, does not see.

Differences from the current implementation¶

--gpu=hip has the device capabilities of --gpu=cuda (the CUDA entry of gpu_device_capabilities): its 512 KiB of stack per thread, and tiles of 32 iterations, the width of a warp rather than of an AMD wavefront (64). A capability entry of its own, with the stack an AMD GPU gives a thread, needs the passes that ask whether gpu_offload_cuda is set rather than what the device can do to ask the capabilities instead.
A copy from pageable host memory (malloc) to the device waits for the launches before it on the stream, in CUDA and in HIP, except for small copies: the host waits for the device at every launch that copies an array. Copying through pinned staging buffers would leave those launches asynchronous.
The pipeline still has its own checks for device limits (the width of the data that reaches a kernel, input/output statements), spread over several passes and the device code generator. They report a compile error rather than assert, because the check does not guarantee everything they look at: symbols a rewrite introduces, widths that are not on the list (integer(8) on Metal), and a derived type with an unused real(8) component that is handed to the kernel whole (see the open questions).
Purity of procedures called from do concurrent is not enforced.
A run-time sized local of a procedure the loop calls whose extent reads the index of a do loop of the procedure is sized only when the extent is that of a section (bounded by the extent of its array), and only for a rank-one allocatable: allocate(w(n(l))) in do l = 1, m, or a local of higher rank, cannot be sized yet, even when the host could work out the largest value before the launch by going over the loop.
A heap value can be stored into an element of an array argument of the loop, but not into a scalar argument or a section, and a loop that also writes a component of an element of an array it stores heap values into is an error. A heap array of rank greater than one cannot be passed whole, a section of a heap array whose elements are not contiguous is not lowered yet, and neither is an array of a heap type of rank greater than one, a component of a heap type that is an array of values of a heap type of a fixed size, an allocatable scalar, pointer or polymorphic component of a heap type, an allocatable component of a heap type in a value the host hands over, or an array of a heap type the host hands over passed to a procedure. A value the host hands over is copied into an image even when its allocatable component is not allocated (the device sees it as allocated with no elements).