Добавочная параллельная сумма для большого набора данных с использованием OpenCL

Question

У меня есть трехмерный вектор F, который в настоящее время заполняется последовательно следующим образом:

// for each particle
for(int p = 0; p < pmax; p++) {
    // and each dimension (n x m x o)
    for(int i = 0; i < n; i++) {
        for(int j = 0; j < m; j++) {
            for(int k = 0; k < o; k++) {
                // do some calulations depending on p, i, j and k and add onto F
                F[i][j][k] += p * i * j * k;
            }   
        }
    }
}

Значение F вычисляется постепенно для каждой частицы.

Теперь я хотел чтобы ускорить процесс с использованием OpenCL. Моей первой попыткой было распараллелить внутренние циклы с помощью 3D-Range-Kernel, который вызывается для каждой частицы. Это действительно хорошо работает, но даже медленнее на GPU, чем на CPU (без применения каких-либо оптимизаций - например, F копируется с хоста на GPU и обратно в каждую итерацию).

Во второй попытке я попытался распараллелить внешний l oop, что (я думаю) является лучшей идеей, чем первая попытка, поскольку ожидается, что pmax будет намного больше, чем m * n * o. Из того, что я прочитал, я думаю, что проблему можно решить с помощью параллельного уменьшения суммы, чего можно привести множество примеров. Однако для этого подхода мне нужно иметь копию F для каждого потока (или рабочего элемента), который не помещается в памяти (получая CL_OUT_OF_RESOURCES здесь).

Мой вопрос сейчас : возможно ли даже распараллелить такую ​​инкрементную сумму, не сохраняя F в памяти несколько раз для нескольких частиц? Если да, как бы выглядел хороший подход? Стоит ли вместо этого придерживаться своей первой попытки и попытаться оптимизировать ее?

Обратите внимание, что я новичок в OpenCL и не очень разбираюсь в методах распараллеливания в целом. Буду признателен за любые подсказки или ссылки на полезные лекции или примеры, спасибо!

Кстати, мои сеансы поиска в Google по этой теме c только приведут меня к вычислению суммы префикса.

Answer 1

Вы можете начать с простого трехмерного ядра, например:

import numpy as np
import pyopencl as cl
import time

m=32
n=32
o=32
pmax=16

ctx = cl.create_some_context()
queue = cl.CommandQueue(ctx)

t1=time.time()
F=np.zeros((m,n,o)).astype(np.int32).flatten(order='F')

mf = cl.mem_flags
F_buf = cl.Buffer(ctx, mf.WRITE_ONLY, F.nbytes)

prg = cl.Program(ctx, """
    __kernel void calc(int pmax, __global int *F) {

    int i = get_global_id(0);
    int j = get_global_id(1);
    int k = get_global_id(2);
    int m = get_global_size(0);
    int n = get_global_size(1);
    int o = get_global_size(2);

    int tmp = i * j * k;
    int sum = 0;
    for(int p = 0; p < pmax; p++) 
        sum += p * tmp;
    F[i*n*o+j*o+k] = sum;
    }
    """).build()

prg.calc(queue, (m,n,o), None, np.int32(pmax), F_buf)
cl.enqueue_copy(queue, F, F_buf)

F=F.reshape((m,n,o), order='F')
F=np.transpose(F, (2, 1, 0))
t2=time.time()

t3=time.time()    
Ftest=np.zeros((m,n,o)).astype(np.int32)
for p in range(pmax):
    for i in range(m):
        for j in range(n):
            for k in range(o):
                Ftest[i][j][k] += p*i*j*k

t4=time.time()

print "OpenCL time:", (t2-t1)
print "CPU time:", (t4-t3)

Результаты теста:

$ python test.py 
Choose platform:
[0] <pyopencl.Platform 'Intel(R) OpenCL HD Graphics' at 0x557007fed680>
[1] <pyopencl.Platform 'Portable Computing Language' at 0x7fab67ff0020>
Choice [0]:
Set the environment variable PYOPENCL_CTX='' to avoid being asked again.
OpenCL time: 0.0124819278717
CPU time: 1.03352808952

$ python test.py 
Choose platform:
[0] <pyopencl.Platform 'Intel(R) OpenCL HD Graphics' at 0x55a2650505a0>
[1] <pyopencl.Platform 'Portable Computing Language' at 0x7fd80775d020>
Choice [0]:1
Set the environment variable PYOPENCL_CTX='1' to avoid being asked again.
OpenCL time: 0.0148649215698
CPU time: 1.11784911156

В зависимости от размера матрицы производительность может отличаться. Также это не значит, что это будет быстрее, чем ваша текущая реализация процессора. В OpenCL производительность зависит от многих факторов, таких как, сколько данных должно быть передано на устройство, достаточно ли вычислений, чтобы сделать его понятным, как осуществляется доступ к данным в ядре: в каком порядке по рабочим элементам ), какой тип памяти - глобальный, локальный, регистры - и т. д.