Обобщенное вычисление скользящего окна на графическом процессоре - PullRequest
7 голосов
/ 05 октября 2011

Вот некоторый код Python, который реализует вычисление скользящего окна для двух трехмерных матриц, X и Y.

import numpy

def sliding_dot( X,Y ) :

    assert X.ndim == Y.ndim == 3
    iw,ih,id = X.shape
    fw,fh,fd = Y.shape

    assert id == fd
    assert fw < iw and fh < ih

    ow,oh = iw-fw+1,ih-fh+1
    out = numpy.zeros( [ow,oh] )

    for x in xrange(ow) :
        for y in xrange(oh) :
            window = X[x:x+fw,y:y+fh,:]
            out[x,y] = numpy.dot( window.flatten(),Y.flatten() )

    return out

#################    

A_dims = (640,480,32)
B_dims = (6,6,32)

A = numpy.random.rand(*A_dims)
B = numpy.random.rand(*B_dims)

sliding_dot(A,B)

В общем, Y всегда намного меньше X по первому и второму измерениям, ноони равны в третьем измерении.

Обратите внимание, что мы могли бы заменить numpy.dot () на любую функцию Y и окна.Это немного отличается от свертки тем, что Y скользит только по первому и второму измерениям X. Я ищу эффективную стратегию для эффективной реализации этого вида вычисления скользящего окна с использованием CUDA.Кто-нибудь хочет предложить мне какое-то направление?Приветствия!

Обновление : вы можете посмотреть, как я работаю в процессе оптимизации, с помощью других пользователей в моем ответе ниже.

Ответы [ 4 ]

4 голосов
/ 11 октября 2011

Попытка спроектировать «обобщенную» реализацию, которая может вместить практически любую операцию, которая вам может понадобиться, станет огромным компромиссом в архитектуре, такой как CUDA. Для вашего примера продукта с точечной точкой, который является типичной операцией сокращения, это довольно полезная реализация:

__constant__ int ldaX[3];
__constant__ int ldaY[3];
__constant__ int dimX[3];
__constant__ int dimY[3];

template<typename real,int blocksize>
__global__ void sliding_k(const real *X, const real *Y, real *out)
{
    __shared__ volatile real buffer[blocksize];

    int tid = threadIdx.x;
    int gid = blockIdx.x * gridDim.y + blockIdx.y;

    real value = (real)0;
    int xpos = (blockIdx.y * ldaX[2]) + (blockIdx.x * ldaX[1]);
    int ypos = 0;
    for(int i=0; i<dimY[0]; i++) {
        for(int jk=tid; jk<ldaY[1]; jk+=blocksize) {
            value += X[xpos+jk] * Y[ypos+jk];
        }
        xpos += ldaX[1];
        ypos += ldaY[1];
    }

    buffer[tid] = value;
    __syncthreads();

# pragma unroll
    for(int i=(tid+32); ((tid<32)&&(i<blocksize)); i+=32)
        buffer[tid] += buffer[i];

    if (tid < 16) buffer[tid] += buffer[tid + 16];
    if (tid < 8)  buffer[tid] += buffer[tid + 8];
    if (tid < 4)  buffer[tid] += buffer[tid + 4];
    if (tid < 2)  buffer[tid] += buffer[tid + 2];
    if (tid == 0) out[gid] = buffer[0] + buffer[1];
}

Вы можете заменить любой тип оператора сокращения на операцию сложения / сложения с умножением с плавающей запятой, которую использует скалярное произведение, и код должен работать нормально. Каждое окно расчета выполняется одним блоком. Параллельной работы достаточно, чтобы при этом размере окна оправдать блок на окно. Это обеспечивает объединенный доступ к глобальной памяти, а на картах Fermi - большое количество обращений к кэш-памяти первого уровня.

Здесь у меня есть только одно допущение в коде, согласно которому третье измерение исходного массива и массива окон равны. Это позволяет внутренним двум циклам «объединяться» в одну операцию из-за общей схемы памяти, которую они разделяют. Запустив тестовый жгут в Python, используя улучшенную версию вашего ссылочного кода, с кодом хоста, написанным на PyCUDA, я получаю следующее:

In [15]: %timeit -n3 -r3 out2=sliding_cuda(A,B)
3 loops, best of 3: 49.8 ms per loop

In [16]: %timeit -n3 -r3 out=sliding_dot(A,B)
3 loops, best of 3: 2.18 s per loop

In [17]: (numpy.abs(out2-out)/numpy.abs(out)).max()
Out[17]: 4.2921323635558404e-15

при работе на Phenom II 3 ГГц с GTX470 с использованием 64 потоковых блоков в 2D сетке 635x475 - т.е. Ускорение примерно в 50 раз, включая загрузку модулей, настройку и передачу памяти с использованием распределяемой памяти хоста. Само ядро ​​примерно в 100 раз быстрее, чем Python, без учета передачи памяти и накладных расходов на установку. Обратите внимание, что это версия с двойной точностью - Python по умолчанию использует арифметику с плавающей запятой двойной точности.

1 голос
/ 06 октября 2011

Хорошо, вот некоторые соображения:

Вы выполняете ~ 640 * 480 итераций numpy.dot, который сам обрабатывает 6 * 6 * 32 элементов.Распараллеливание точечного продукта едва ли того стоит: 192 параллельных потока недостаточно для GPU, а сокращение на CUDA - дополнительные проблемы.Итак, IMO, лучший способ распараллелить вашу задачу - назначить по одному элементу выходного массива каждому потоку.

Теперь о памяти: выходной массив будет в глобальной памяти, выбор невелик.Для входных данных A выглядит неплохо для текстурной памяти, поскольку смежные потоки обращаются к смежным элементам.В качестве альтернативы, вы можете вручную «кэшировать» его в разделяемой памяти, но в этом случае это выглядит не слишком выгодно по сравнению с простым использованием текстуры.Для B разделяемая память не годится, так как это может привести к конфликтам банков, поскольку при вычислении точечного продукта все потоки в полусфере имеют доступ к одному и тому же элементу B (вы можете начать суммирование с разных элементов в разных потоках, ноэто (опять же) не выглядит многообещающе).Таким образом, выбор является либо текстурой, либо постоянным.Я голосую за константу, так как (а) постоянная память подходит для данных, к которым имеют доступ все потоки на устройстве, (б) вы не будете загрязнять кеш текстур.

Выше приведены только мои догадки, ичтобы реально добиться хорошей производительности, вам лучше попробовать разные варианты ...

Обновление, касающееся вашей наивной реализации

for (int Yi = 0; Yi < Ydims[0]; Yi++ )

Здесь вы получаете доступ к глобальной памяти накаждая итерация.Это огромный убийца производительности.Поскольку у вас есть 3 измерения, вам лучше заменить int *Ydims на int3 Ydims (то же самое для Xdims и outdims).

out[out_indx] += X[X_indx]*Y[Y_indx];

Опять же, очень плохая идея.Создайте переменную регистра и выполните все операции с ней.Запись в глобальный массив только один раз в конце ядра.

Эти оптимизации - первое, что вы должны сделать.Во-вторых, вам нужно сделать X и Y 3D-текстуры, чтобы доступ к ним был кэширован.Полагаю, после этого CUDA превзойдет ЦП.

Для дальнейшей оптимизации вам лучше прочитать Руководство по рекомендациям CUDA C .Это должно быть прочитано, и вы получите гораздо лучшее представление о том, как написать эффективный код для графического процессора (сейчас ваша реализация слишком проста)

0 голосов
/ 10 октября 2011

Вы можете попытаться отделить ваши чтения от ваших сумм из ваших магазинов.

Таким образом, каждое ядро ​​должно иметь 3 раздела:

  1. Чтение из памяти текстур, сохранение в общей памяти для всего блока

    __shared blockX[ Ydims.z ][ Ydims.y ][ Ydims.x ];
    __shared blockY[ Ydims.z ][ Ydims.y ][ Ydims.x ];
    // NOTE: MAKE EACH THREAD LOAD k ELEMENTs * 2 rather than each thread loading Ydims.X*Y*Z elements
    blockX[k][yj][yi] = ...
    blockY[k][yj][yi] = ...
    __syncthreads(); // <-- critical -- all threads in block must finish
    // reading from shared memory before any may use the values.
    
  2. #pragma Разверните ваши for петли.
    Это значительно увеличит ваш ILP и будет иметь гораздо меньше ветвлений для ваших постоянных размеров петли

  3. Убедитесь, что ваш доступ к общей памяти ограничен надлежащим образом, в противном случае конфликты в банках приведут к снижению вашей производительности.

0 голосов
/ 10 октября 2011

v0.1 - Наивная реализация

Вот моя первая, наивная попытка сделать эту работу:

__global__ void sliding_dot(float *out, int *outdims, float *X, int *Xdims, float *Y, int *Ydims )
{
    int i = threadIdx.x + blockDim.x * blockIdx.x;
    int j = threadIdx.y + blockDim.y * blockIdx.y;
    int Y_indx = 0;
    int X_indx = 0;
    if ( i < outdims[0] & j < outdims[1] )
    {
        int out_indx = j + i*outdims[1];
        for (int Yi = 0; Yi < Ydims[0]; Yi++ )
        {
            for (int Yj = 0; Yj < Ydims[1]; Yj++ )
            {
                for (int k = 0; k < Ydims[2]; k++ )
                {
                    Y_indx = k + Yj*    Ydims[2] + Yi*    Ydims[2]*Ydims[1];
                    X_indx = k + (j+Yj)*Xdims[2] + (i+Yi)*Xdims[2]*Xdims[1];
                    out[out_indx] += X[X_indx]*Y[Y_indx];
                }
            }
        }
    }
}

Пока что результаты менее желательны. С размером блока (32,32,1) и размерами сетки p, q выбрано так, что p * 32> = outdims [0] и q * 32> = outdims [1]:

method=[ sliding_dot ] gputime=[ 7013.280 ] cputime=[ 18.000 ] occupancy=[ 0.667 ] 
method=[ sliding_dot ] gputime=[ 6945.184 ] cputime=[ 7.000 ] occupancy=[ 0.667 ] 
method=[ sliding_dot ] gputime=[ 6990.816 ] cputime=[ 6.000 ] occupancy=[ 0.667 ] 
method=[ sliding_dot ] gputime=[ 6931.648 ] cputime=[ 6.000 ] occupancy=[ 0.667 ] 

v0.2 - texture<float,1>

Я надеюсь, что все узнают так же много, как и я! Я последовал советам @ aland и значительно ускорился:

texture<float,1> X;
texture<float,1> Y;

__global__ void dotconv(float *out, int2 outdims, int3 Xdims, int3 Ydims )
{
    int i = threadIdx.x + blockDim.x * blockIdx.x;
    int j = threadIdx.y + blockDim.y * blockIdx.y;

    if ( i < outdims.x & j < outdims.y )
    {
        int out_indx = j + i*outdims.y;
        float total = 0.0f;
        int X_indx = 0;
        int Y_indx = 0;
        for (int Yi=0; Yi<Ydims.x; Yi++ )
        {
            for (int Yj=0; Yj<Ydims.y; Yj++ )
            {
                for (int k=0; k<Ydims.z; k++ )
                {
                    Y_indx = k + Yj*    Ydims.z + Yi*    Ydims.z*Ydims.y;
                    X_indx = k + (j+Yj)*Xdims.z + (i+Yi)*Xdims.z*Xdims.y;
                    total += tex1Dfetch(X,X_indx)*tex1Dfetch(Y,Y_indx);
                }
            }
        }
        out[out_indx] = total;
    }
}

Но мы все еще работаем не так быстро, как процессор:

method=[ dotconv ] gputime=[ 2224.928 ] cputime=[ 24.000 ] occupancy=[ 0.667 ] 
method=[ dotconv ] gputime=[ 2222.592 ] cputime=[ 7.000 ] occupancy=[ 0.667 ] 
method=[ dotconv ] gputime=[ 2225.216 ] cputime=[ 10.000 ] occupancy=[ 0.667 ] 
method=[ dotconv ] gputime=[ 2222.752 ] cputime=[ 10.000 ] occupancy=[ 0.667 ] 

v0.3 - texture<float,3>

texture<float,3,cudaReadModeElementType> X;
texture<float,3,cudaReadModeElementType> Y;

__global__ void dotconv(float *out, int2 outdims, int3 Xdims, int3 Ydims )
{
    int i = threadIdx.x + blockDim.x * blockIdx.x;
    int j = threadIdx.y + blockDim.y * blockIdx.y;
    if ( i < outdims.x & j < outdims.y )
    {
        int out_indx = j + i*outdims.y;
        float total = 0.0f;
        for (int Yi=0; Yi<Ydims.x; Yi++ )
        {
            for (int Yj=0; Yj<Ydims.y; Yj++ )
            {
                for (int k=0; k<Ydims.z; k++ )
                {
                    total += tex3D(X,k,j+Yj,i+Yi) * tex3D(Y,k,Yj,Yi);   
                }
            }
        }
        out[out_indx] = total;
    }
}

Это на самом деле немного медленнее, чем v0.2

method=[ dotconv ] gputime=[ 2403.360 ] cputime=[ 35.000 ] occupancy=[ 0.667 ] 
method=[ dotconv ] gputime=[ 2392.160 ] cputime=[ 15.000 ] occupancy=[ 0.667 ] 
method=[ dotconv ] gputime=[ 2396.448 ] cputime=[ 15.000 ] occupancy=[ 0.667 ] 
method=[ dotconv ] gputime=[ 2398.880 ] cputime=[ 16.000 ] occupancy=[ 0.667 ] 

Спасибо за ваши предложения!

...