1
\$\begingroup\$

Classic pointer jumping algorithm summarizing values from an array adapted to run on a GPU instead of a PRAM.

Read about openCL and Aparapi yesterday for the first time, so this is my first trial: please let me know if I use Aparapi correctly and if something can be improved.

I've tested the code on my integrated intel GPU (maxWorkWorkGroupSize=256, maxComputeUnits=48) with random arrays of size up to 16M elements and it works correctly, but slower than sequential adding on the CPU: for-loop on the CPU takes about 20ms on average, while my pointer-jumping implementation on the GPU about 80ms.
My guess is, that pointer-jumping is a too much memory-bound algorithm to be effective on a GPU, but as I'm new to GPU programming, I may be completely wrong here...
The other reason may be that intel has only 16 barrier registers (and only 64kB local memory shared among running work-groups), so only up to 16 work-groups can run in parallel.

Java language level is set to 11 and aparapi is version 3.0.0.

// Copyright (c) Piotr Morgwai Kotarbinski, Licensed under the Apache License, Version 2.0
package pl.morgwai.sample;

import com.aparapi.Kernel;
import com.aparapi.Range;



/**
 * Performs pointer-jumping on an array using a GPU with Aparapi. By default sums values from the
 * array, but subclasses may override {@link #accumulateValue(int, int)} to do something else.
 * <p>Usage:</p>
 * <pre>
 * double[] myArray = getDoubleArray();
 * double sum = PointerJumpingKernel.calculateSum(myArray);</pre>
 */
public class PointerJumpingKernel extends Kernel implements AutoCloseable {



    public static double calculateSum(double[] input) {
        try (var kernel = new PointerJumpingKernel()) {
            return kernel.accumulateArray(input);
        }
    }

    protected PointerJumpingKernel() {}



    double[] input;  // input values to accumulate
    double[] results; // results from all work-groups

    /**
     * Group's local copy of its slice of the {@code input}.
     */
    @Local protected double[] localSlice;
    @Local int[] next; // next[i] is a pointer to the next value in localSlice not yet accumulated
            // by the processing element with localId == i, initially next[i] == i+1 for all i



    /**
     * Triggers GPU execution of {@link #run() pointer-jumping} on distinct slices of {@code input}
     * in separate work-groups, after that recursively accumulates results from all work-groups.
     * Recursion stops when everything is accumulated into a single value.
     * @return value accumulated from the whole input.
     */
    protected double accumulateArray(double[] input) {
        this.input = input;
        int groupSize = Math.min(input.length, getTargetDevice().getMaxWorkGroupSize());
        int numberOfGroups = input.length / groupSize;
        if (groupSize * numberOfGroups < input.length) numberOfGroups++;  // padding the last group
        int paddedSize = groupSize * numberOfGroups;
        results = new double[numberOfGroups];
        localSlice = new double[groupSize];
        next = new int[groupSize];
        execute(Range.create(paddedSize, groupSize));
        if (numberOfGroups == 1) return results[0];
        return accumulateArray(results);
    }



    /**
     * Pointer-jumping procedure executed by a single processing element. Each work-group first
     * copies its slice of the {@code input} to {@link #localSlice} array in group's local memory.
     * Next, the main pointer-jumping loop is performed on the {@link #localSlice}.
     * Finally, the 1st processing element writes group's accumulated result to {@link #results}
     * array.
     */
    @Override
    public final void run() {
        int i = getLocalId();
        int globalIndex = getGlobalId();
        int groupSize = getLocalSize();
        int acivityIndicator = i;// Divided by 2 at each step of the main loop until odd.
                // When odd, the given processing element stays idle (just checks-in at the barrier)

        // copy group's slice into local memory and initialize next pointers
        if (globalIndex < input.length - 1) {
            next[i] = i + 1;
            localSlice[i] = input[globalIndex];
        } else {
            next[i] = getGlobalSize(); // padding in the last group: point beyond the array
            if (globalIndex == input.length - 1) localSlice[i] = input[globalIndex];
        }
        localBarrier();

        // main pointer-jumping loop
        while (next[0] < groupSize) { // run until the whole group is accumulated at index 0
            if ( (acivityIndicator & 1) == 0 && next[i] < groupSize) {
                accumulateValue(next[i], i);
                next[i] = next[next[i]];
                acivityIndicator >>= 1;
            }
            localBarrier();
        }
        if (i == 0) results[getGroupId()] = localSlice[0];
    }

    /**
     * Accumulates value from {@code fromIndex} in {@link #localSlice} into {@code intoIndex}.
     * Subclasses may override this method to do something else than summing.
     * Subclasses should then provide a static method that creates a kernel and calls
     * {@link #accumulateArray(double[])} similarly to {@link #calculateSum(double[])}.
     */
    protected void accumulateValue(int fromIndex, int intoIndex) {
        localSlice[intoIndex] += localSlice[fromIndex];
    }



    @Override
    public final void close() {
        dispose();
    }
}

For reference here is the code I was using to test it:

static Random random = new Random();

public static void runPointerJumpingExample(int size) {
    double[] values = new double[size];
    for (int i = 0; i < size; i++) {
        values[i] = random.nextDouble();
    }

    double val = 0;
    var start = System.nanoTime();
    for (int i = 0; i < size; i++) {
        val += values[i];
    }
    System.out.println(String.format(
            "cpu: %1$15d,  result: %2$20.12f", System.nanoTime() - start, val));

    start = System.nanoTime();
    double result = PointerJumpingKernel.calculateSum(values);
    System.out.println(String.format(
            "gpu: %1$15d,  result: %2$20.12f\n", System.nanoTime() - start, result));

    if (Math.abs(result - val) > 0.00001) {
        throw new RuntimeException("error!"
                + "\nexpected: " + val
                + "\nresult:   " + result);
    }
}

public static void main(String[] args) {
    for (int i = 0; i < 20; i++)
        runPointerJumpingExample(16*1024*1024);
    System.out.println("bye bye!");
}

...and for non-java openCL folks, here is the generated openCL code:

#pragma OPENCL EXTENSION cl_khr_fp64 : enable

typedef struct This_s{
   __local double *localSlice;
   __global double *input;
   int input__javaArrayLength;
   __local int *next;
   __global double *results;
   int passid;
}This;
int get_pass_id(This *this){
   return this->passid;
}
void pl_morgwai_sample_AparapiSample$PointerJumpingProductKernel__accumulateValue(This *this, int fromIndex, int intoIndex){
   this->localSlice[intoIndex]  = this->localSlice[intoIndex] * this->localSlice[fromIndex];
   return;
}
__kernel void run(
   __local double *localSlice, 
   __global double *input, 
   int input__javaArrayLength, 
   __local int *next, 
   __global double *results, 
   int passid
){
   This thisStruct;
   This* this=&thisStruct;
   this->localSlice = localSlice;
   this->input = input;
   this->input__javaArrayLength = input__javaArrayLength;
   this->next = next;
   this->results = results;
   this->passid = passid;
   {
      int i = get_local_id(0);
      int globalIndex = get_global_id(0);
      int groupSize = get_local_size(0);
      int acivityIndicator = i;
      if (globalIndex<(this->input__javaArrayLength - 1)){
         this->next[i]  = i + 1;
         this->localSlice[i]  = this->input[globalIndex];
      } else {
         this->next[i]  = get_global_size(0);
         if (globalIndex==(this->input__javaArrayLength - 1)){
            this->localSlice[i]  = this->input[globalIndex];
         }
      }
      barrier(CLK_LOCAL_MEM_FENCE);
      for (; this->next[0]<groupSize; barrier(CLK_LOCAL_MEM_FENCE)){
         if ((acivityIndicator & 1)==0 && this->next[i]<groupSize){
            pl_morgwai_sample_AparapiSample$PointerJumpingProductKernel__accumulateValue(this, this->next[i], i);
            this->next[i]  = this->next[this->next[i]];
            acivityIndicator = acivityIndicator >> 1;
         }
      }
      if (i==0){
         this->results[get_group_id(0)]  = this->localSlice[0];
      }
      return;
   }
}
\$\endgroup\$
1
\$\begingroup\$

While pointer-jumping may be applied to arrays, it is rather intended for linked-lists. For arrays parallel reduction algorithm is better suited as it takes advantage of constant access to any element and does not need next pointers.

There are few optimization commonly applied to parallel reduction that can probably be applied to the pointer-jumping in the OP.

The most significant is probably "last warp unrolling" which takes advantage of SIMD. When the number of active threads does not exceed hardware SIMD width ("warp size" in CUDA terms, "wavefront width" in AMD: number of processing elements performing instructions synchronously: 32 on Nvidia, 64 on AMD, 8-32 on Intel), then it is possible to omit local barrier. Pointer to the array must be declared volatile though, which is currently not supported by aparapi :-(

See this Nvidia slideshow for details and other possible optimizations.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.