Migrating reduction operations to SYCL in a molecular docking application
Tags: sycl, tips, oneapi
Published: 2023-08-10 22:30:31

I completed porting of a molecular docking application from CUDA to SYCL using the Intel&reg; DPC++ Compatibility Tool (Compatibility Tool) in June 2021. Let me share selected techniques that I used without delving into the details of the docking application. If you want to learn how to use this tool to migrate CUDA applications to SYCL, please refer to [1].

The Compatibility Tool adds comments in the code where manual migration may be required. Typically, the manual changes required fall into two categories. First, changes are required for the code to compile and make the code functionally correct. Other changes are necessary to get better performance. Here, I will cover code that uses the operation of 'reduction'. Reductions are frequently used in High Performance Computing and scientific applications and can be performance hotspots. The first example finds the sum of integers and the second finds the minimum of floats and the identifier of the run that corresponds to the minimum.

## Integer Reductions to find the number of evaluations

The  docking application performs integer reductions to keep a running count of the number of score evaluations. This reduction is  implemented as a multi-line macro in CUDA as shown below.

```CUDA
#define REDUCEINTEGERSUM(value, pAccumulator)     
    if (threadIdx.x == 0)     
    {     
        *pAccumulator = 0;     
    }     
    __threadfence();     
    __syncthreads();     
    if (__any_sync(0xffffffff, value != 0))     
    {     
        uint32_t tgx            = threadIdx.x & cData.warpmask;     
        value                  += __shfl_sync(0xffffffff, value, tgx ^ 1);     
        value                  += __shfl_sync(0xffffffff, value, tgx ^ 2);     
        value                  += __shfl_sync(0xffffffff, value, tgx ^ 4);     
        value                  += __shfl_sync(0xffffffff, value, tgx ^ 8);     
        value                  += __shfl_sync(0xffffffff, value, tgx ^ 16);     
        if (tgx == 0)     
        {     
            atomicAdd(pAccumulator, value);     
        }     
    }     
    __threadfence();     
    __syncthreads();     
    value = *pAccumulator;     
    __syncthreads();
```

Let us review what this code is doing:

* The code is called for each work item (thread) in a work group (warp)
* *\*pAccumulator* is where the final sum is stored summing across all work items
* The combination of *__threadfence()* and *__syncthreads()* guarantees memory consistency and synchronizes threads in the warp at the point of the call.
* The *__any_sync()* call executes the block for those non-exited threads for which *'value != 0'*
* The following *__shfl_sync* calls do a tree-wise summing with the final sum available in the first thread in the warp in variable *value*
* The *value* is then added to the Accumulator atomically with *atomicAdd* and finally all threads assign the sum to the *value* variable.

For more details about these CUDA calls please refer to [2].

The Compatibility tool was not able to automatically migrate this code with the following comments.

```SYCL
/*
DPCT1023:40: The DPC++ sub-group does not support mask options for sycl::ext::oneapi::any_of.

DPCT1023:41: The DPC++ sub-group does not support mask options for shuffle.

DPCT1007:39: Migration of this CUDA API is not supported by the Intel(R) DPC++ Compatibility Tool.
*/
```

However, SYCL supports a rich set of functions for performing reductions. In this case, the reduce_over_group() function in SYCL can be used to create the same functionality as the above code as follows.

```SYCL

#define REDUCEINTEGERSUM(value, pAccumulator)     
        int val = sycl::reduce_over_group(item_ct1.get_group(), value, std::plus<>());      
        *pAccumulator = val;     
        item_ct1.barrier(sycl::access::fence_space::local_space);
```

The *sycl::reduce_over_group* is a collective function. The usage of this function simplifies the macro. The function takes the group, the value to be reduced, and the reduction operation which in this case is plus or summation. The function can adapt to varied sizes of work groups in SYCL and will use the best available optimizations available per the compiler and run-time.

## Finding the minimum energy

In another part of the application, a block of CUDA threads perform shuffles to find the minimum of scores *v0* and the corresponding identifier *k0* of the run in the simulation that is the minimum score. The CUDA code calls a macro WARPMINIMUM2 (not shown) which in turn calls another macro WARPMINIMUMEXCHANGE (shown) with *mask* set to 1, 2, 4, 8, and 16.

```CUDA
#define WARPMINIMUMEXCHANGE(tgx, v0, k0, mask)     
    {     
        float v1    = v0;     
        int k1      = k0;     
        int otgx    = tgx ^ mask;     
        float v2    = __shfl_sync(0xffffffff, v0, otgx);     
        int k2      = __shfl_sync(0xffffffff, k0, otgx);     
        int flag    = ((v1 < v2) ^ (tgx > otgx)) && (v1 != v2);     
        k0          = flag ? k1 : k2;     
        v0          = flag ? v1 : v2;     
    }
```

The *__shfl_sync* provides a way of moving a value from one thread to other threads in the warp in one instruction. In this code snippet *__shfl_sync* gets the *v0* or *k0* value from the thread identified by the *otgx* mask and saves it in *v2*, *k2* variables. We then compare *v1* with *v2* to set *flag* and eventually store the minimum in *v0* and the run identifier for this minimum in *k0*.

Compatibility Tool could not completely migrate this code and included this comment as the reason it could not. However, Compatibility Tool correctly replaced the *__shfl_sync* call with SYCL *shuffle* call as shown in the below diff which shows the manual change.

```SYCL
/*
DPCT1023:57: The DPC++ sub-group does not support mask options for shuffle.
*/
```

This comment indicates that the *shuffle* call in SYCL does not use a mask as shown below.

```diff
#define WARPMINIMUMEXCHANGE(tgx, v0, k0, mask)     
        {     
                float v1 = v0;     
                int k1 = k0;     
                int otgx = tgx ^ mask;     
-               float v2 = item_ct1.get_sub_group().shuffle(energy, otgx);     
+               float v2 = item_ct1.get_sub_group().shuffle(v0, otgx);  
-               int k2 = item_ct1.get_sub_group().shuffle(bestID, otgx);     
+               int k2 = item_ct1.get_sub_group().shuffle(k0, otgx);     
                int flag = ((v1 < v2) ^ (tgx > otgx)) && (v1 != v2);     
                k0 = flag ? k1 : k2;     
                v0 = flag ? v1 : v2;     
        }
```

In this case, Compatibility Tool performed incorrect variable substitution for *v0* and *k0* in the shuffle calls using *energy* and *bestID* variables from the caller function. We manually fixed this by replacing *energy* with *v0* and *bestID* with *k0*. This bug has been fixed in recent versions of the Compatibility Tool.

## Summary

In summary, reduction operations in CUDA applications may not be migrated correctly by the Compatibility Tool. Review the comments provided by the tool to understand if manual migration is necessary and what change might be required. A good understanding of the original CUDA code will then help to make manual changes to develop functionally correct code in SYCL.

[1] <https://www.intel.com/content/www/us/en/docs/dpcpp-compatibility-tool/get-started-guide/2023-1/overview.html>

[2] <https://developer.nvidia.com/blog/using-cuda-warp-level-primitives/>