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;
}
}
