8 Special programming patterns

8.1  Matrix/vector expressions

8.2  Loop parallelization/serialization

1. All variables that are being used inside the loop must have a static type (i.e. explicitly typed as vec, mat, cube, see Section 2.2↑) or a static type can be inferred from the context (type inference or through explicit/implicit specialization, see Chapter 6↑). Practically speaking: only types that can be used inside __kernel__ or __device__ functions are allowed.
2. For the best performance, for slicing operations A[a..b,2], the dimensions a and b must be constant and known at compile time (either specified explicitly, or obtained through constant propagation). When these constants are not known, a nested kernel function will be generated. Depending on the context, this may or may not involve dynamic kernel memory (see Section 8.3↓)
3. The for-loop nest must be perfect (no code in between subsequent for statements, no dependencies between the loop boundaries, no break, no continue) for example:
   for m=0..size(x,0)-1
   for n=5..size(y,0)-1
   endfor
   endfor
   This is an example of a imperfect loop:
   for m=0..size(x,0)-1
   for n=m..size(y,0)-m
   endfor
   endfor
   (Note: certain types of imperfect loops can also be handled using the imperfect loop function transform, this requires placing {!function_transform enable="imperfectloop"} inside the for-loop)
4. Only types that can be used inside kernel or device functions are allowed (e.g., types based on class and mutable class, but not dynamic class).
5. When host (i.e. non kernel/device) functions are called from a loop, the loop is not eligible for parallelization/serialization. In particular, only a limited number of built-in functions are supported. Functions that interact with/take variables with unspecified types (such as print, load, save, ...) are not supported.
6. Data dependencies/conflicts between different iterations are detected and not allowed. In case a dependency is detected, the loop can be serialized. In this case, the for-loop will be natively compiled (in C++) and executed on the CPU in single-threaded mode.
7. Advanced kernel function features such as shared memory and thread synchronization (see Subsection 2.4.4↑) are not supported, because these functions often require low-level access to the block position (blkpos) and dimensions (blkdim).

• {!parallel for}: forces the loop to be parallelized, despite dependencies (i.e., potential data races) detected between the variables. In the latter case, a warning message is shown.
• {!parallel for; dim=n}: forces the following n loops to be parallelized jointly. Note that this operation requires that the for-loop openings are placed directly after each other, without any intermediate statements (apart from comments). Explicitly specifying the dimension allows the outer n loops to be parallelized even in cases that there are some additional inner loops. When the dims parameter is omitted, the compiler attempt to select a maximal value for dims depending on the inner for-loops.
• {!parallel for; multi_device=true}: causes the for-loop to be parallelized over multiple devices (at least, for runtime engines that support multiple devices, see Chapter 11↓).
• {!serial for}: serializes the for-loop, for execution on devices that support serial execution (e.g., the CPU).
• {!interpreted for}: forces the loop to be interpreted. This comes at a computational cost, but can still be useful for e.g., debugging purposes

8.2.1 While-loop serialization

8.2.2 Example: gamma correction

8.3  Dynamic kernel memory

• When the algorithm requires to process several medium-sized (16x16) to large-sized (e.g. 128x128) matrices. For example, a kernel function that performs matrix operations for every pixel in the image. The size of the matrices may or may not be known in advance. [K]  [K] Note that to avoid use of dynamic kernel memory in some cases, types can be annotated with size information (see Section 3.3↑).
• Efficient handling of multivariate functions that are applied to (non-overlapping or overlapping) image blocks.
• When the algorithm works with dynamic data structures such as linked lists, trees, it is also often necessary to allocate “nodes” on the fly.
• To use some sort of “/scratch” memory that does not fit into the GPU shared memory (note: the GPU shared memory is 32K, but this needs to be shared between all threads - for 1024 threads this is 32 bytes private memory per thread). Dynamic memory does not have such a stringent limitation. Moreover, dynamic memory is not shared and disposed either 1) immediately when the memory is not needed anymore or 2) when a GPU/CPU thread exists. Correspondingly, when 1024 threads would use 32K each, this will require less than 32MB, because the threads are logically in parallel, but not physically.

• Memory can be allocated and released on the fly (for example, scratch pads within a kernel that exceed the shared memory size)
• Functions can be written in a generic way so that they can handle any vector/matrix dimensions

• Despite the internal parallel memory allocator being quite efficient, code that uses dynamic kernel memory is generally 10x-25x slower than code that does not rely on dynamic kernel memory. This is due to 1) the execution cost of allocation and disposal operations (which are performed behind the scenes) and 2) the storage of the data in the global memory which is slower than, e.g., register memory.
• The restrictions of the dynamic memory subsystem may cause runtime errors to be generated, for example when insufficient dynamic kernel memory is available. To alleviate this problem, it is possible to increase the amount of dynamic kernel memory by a configuration setting (in Redshift: Program Settings/runtime/kernel dynamic memory reserve), although this memory is then not available for other purposes.

• Use of array size constraints (see Section 3.3↑). The resulting variables have the advantages that the memory is stored in the registers or on the stack (which is highly efficient).
• Use of size-parametrized generic types (see Section 6.7↑).
• Memory preallocation: allocate one large matrix (using zeros, uninit functions) outside the kernel function and pass the matrix to the kernel function

8.3.1 Examples

Kernel version

Loop version

8.3.2 Memory models

1. Concurrent memory model (default)
   In the concurrent memory model, the computation device (e.g. GPU) autonomously manages a separate memory heap that is reserved for dynamic objects. The size of the heap can be configured in Quasar and is typically 32MB.\begin_inset Separator latexpar\end_inset
   The concurrent memory model is extremely efficient when all threads (e.g.  ≥  512) request dynamic memory at the same time. The memory allocation is done by a specialized parallel allocation algorithm that significantly differs from traditional sequential allocators.
   For efficiency, there are some internal limitations on the size of the allocated blocks:
   • The total amount of kernel dynamic memory is by default 128 MB, but can be configured (see Redshift/program settings/Runtime/Memory management)
   • The maximum size is by default 64 KiB (65536 bytes), but can be configured (see Redshift/program settings/Runtime/Memory management).
   • The minimum size is by default 64 bytes and scales linearly with the setting of maximum size. When the allocation request size is smaller than this minimum size, the request size is rounded up to the minimum size.
   • All Quasar modules need to use the same kernel dynamic memory settings: this is necessary so that dynamic memory blocks can be accessed correctly (allocated, deallocated) in between different modules
   Also note that the maximum size also incorporates the management overhead (reference counting, block sizes etc.), so in case you intend to allocate a 128 × 128 matrix in 32-bit precision float mode (64K) it is required to select 128 KiB.
   The concurrent memory model is primarily designed to enable a large number of concurrent allocations of relatively small memory blocks (for example to perform a local matrix calculation). It is generally not memory-efficient to allocate large memory blocks from a kernel function in the concurrent memory model. This is because often thousands of threads may be launched on the GPU, and correspondingly, thousands of dynamically allocated memory blocks may be active for a certain local calculation. However, the temporary memory blocks may easily consume several megabytes of GPU memory, and this may restrict the amount of GPU memory available for other purposes. For example, consider 2800 concurrent threads using each 128KB of data. This corresponds to 377MB of GPU memory!
   Therefore, for large dynamic memory allocations, consider using the cooperative memory model.
2. Cooperative memory model\begin_inset Separator latexpar\end_inset
   In the cooperative memory model, the kernel function requests memory directly to the Quasar allocator. This way, there are no limitations on the size of the allocated memory. Also, the allocated memory is automatically garbage collected. To enable the cooperative memory model from a kernel function, use:{!kernel memorymodel="cooperative"}inside the kernel function.
   Because the GPU cannot launch callbacks to the CPU, this memory model requires the kernel function to be executed on the CPU.
   Advantages:
   • The maximum block size and the total amount of allocated memory only depend on the available system resources.
   Limitations:
   • Can only be used from CPU kernel functions.
   • The Quasar memory allocator uses locking (to limited extend), so simultaneous memory allocations on all processor cores may be expensive.
     

8.3.3 Features

• Device functions can also use dynamic memory. The functions may even return objects that are dynamically allocated.
• The following built-in functions are supported and can be used from within kernel and device functions:\begin_inset Separator latexpar\end_inset
  zeros, czeros, ones, uninit, eye,
  copy, reshape, repmat, shuffledims,
  seq, linspace, real, imag, complex,
  mathematical functions matrix/matrix
  multiplication matrix/vector multiplication

  Note: when the above functions are applied to arrays with size constraints (see Section 3.3↑), no dynamic kernel memory will be used.

8.3.4 Performance considerations

• Global memory access: code relying on dynamic memory may be slow (for linear filters on GPU: 4x-8x slower), not because of the allocation algorithms, but because of the global memory accesses. However, it all depends on what you want to do: for example, for non-overlapping block-based processing (e.g., blocks of a fixed size), the dynamic kernel memory is an excellent choice.
• Static vs. dynamic allocation: when the size of the matrices is known in advanced, static allocation (e.g. outside the kernel function may be used as well). The dynamic allocation approach relieves the programmer from writing code to pre-allocate memory and calculating the size as a function of the size of the data dimensions. The cost of calling the functions uninit, zeros is negligible to the global memory access times (one memory allocation is comparable to 4-8 memory accesses on average - 16-32 bytes is still small compared to the typical sizes of allocated memory blocks). Because dynamic memory is disposed whenever possible when a particular threads exists, the maximum amount of dynamic memory that is in use at any time is much smaller than the amount of memory required for pre-allocation.
• Use vecX, matXxY, cubeXxYxZ, ... types (with size constraint) whenever your algorithm allows it (see Section 3.3↑). This completely avoids using global memory, by using the registers instead. Once a vector of length 17 is created, the vector is allocated as dynamic kernel memory.
• Avoid writing code that leads to thread divergence: in CUDA, instructions execute in warps of (typically) 32 threads. Within a warp, instructions of different threads are executed at the same time (called implicit SIMD). Control flow instructions (if, match, repeat, while) can negatively affect the performance by causing threads of the same warp to diverge; that is, to follow different execution paths. Then,the different execution paths must be serialized, because all of the threads of a warp share a program counter. Consequently, the total number of instructions executed for this warp is increased. When all the different execution paths have completed, the threads converge back to the same execution path.
• To obtain best performance in cases where the control flow depends on the position (i.e., pos or blkpos), the controlling condition should be written so as to minimize the number of divergent warps.

8.4  Map and Reduce pattern

• The following operators are supported: addition (+=), subtraction (-=), multiplication (*=), minimum (__=), maximum (^^=), bitwise AND (&=), bitwise OR (|=) and bitwise XOR (~=).
• The compiler allows maximal flexibility in defining the body of the loop: i.e., it is possible to use control structures (if, loop, while, ...). Also, other intermediate values can be calculated, hence it is not a prerequisite that the loop only consists of summation statements such as in the above example.
• The reduction may also perform a dimension reduction (for example, summing a matrix along the rows). Currently, the Quasar compiler supports dimension reductions for which the number of dimensions is reduced with 1 at the time (dimension reduction in multiple dimensions at the time may be supported in tfuture versions).
• There is an upper limit on the number of accumulators (due to the size limit of the shared memory). For 32-bit floating point, up to 32 accumulators and for 64-bit floating point, up to 32 accumulators are supported. When the upper limit is exceeded, the generated code will still work, but the block size will silently be reduced (see also Section 9.9↓). This, together with the impact on the occupancy (due to high number of registers being used) might lead to a performance degradation.

8.5  Cumulative maps (prefix sum)

• The following operators are supported: addition (+=), subtraction (-=), multiplication (*=), minimum (__=), maximum (^^=), bitwise AND (&=), bitwise OR (|=) and bitwise XOR (~=).
• Again, the compiler allows maximal flexibility in defining the body of the loop: i.e., it is possible to use control structures (if, loop, while, ...).
• One restriction is that the dependency should be im[m,n]  ←  im[m-k,n] where k is always 1. Higher order cumulative patterns (e.g., as in higher-order IIR filters) are not supported yet. In this case, it is recommended to decompose the filter into a cascade of first order IIR filters. Each filter of the cascade can then by individually parallelized using the cumulative map pattern.

8.6  Meta functions

• $eval(.): compile-time evaluation of expressions
• $str(.): conversion of an expression to string
• $subs(a=b,.): substitution of a variable by another variable or expression
• $check(.): checks the satisfiability of a given condition (the result is either valid, satisfiable or unsatisfiable), based on the information that the compiler has at this point.
• $assump(.): returns an expression with the assertions of a given variable
• $simplify(.): simplifies boolean expressions (based on the information of the compiler, for example constant values etc.)
• $args[in](.): returns an expression with the input arguments of a given function.
• $args[out](.): returns an expression with the input arguments of a given function.
• $nops(.): returns the number of operands in the expression
• $op(.,n): returns the n-th operand of the expression
• $ubound(.): calculates an upper bound for the given expression
• $specialize(func,.): performs function specialization (see Section 6.5↑)
• $inline(lambda)(...): performs inlining of lambda expressions/functions
• $ftype(x) with x="__host__"/"__device__"/"__kernel__": determines whether we are inside a host, device or kernel function.
• $typerecon(x,y) : reconstructs the type of the specified function with a given set of (specialized) input parameters. For example, given a function f = x -> x, $typerecon(f,type(x,scalar)) will return [scalar->scalar]. This meta function is mainly used internally in the compiler.
• $target(.): indicates whether we are compiling for the given target platform (see Subsection 17.2.6↓).

• Most of these functions (and in particular $eval, $check, $specialize, $typerecon and $inline) are only provided for testing and should not be used from user-code.
• The function $ftype is useful in combination with reductions with where clause (Subsection 4.9.4↑), to express that the reduction may only be applied in a device/kernel or host function (also see [functions](Functions-in-Quasar)). For example:\begin_inset Separator latexpar\end_inset
  reduction x -> log(x) = x - 1 where abs(x - 1) < 1e-1 && $ftype("__device__")

Example: copying the type and assumptions from one variable to another