Метапрограммирование в Cuda

template <unsigned int blockSize>
__global__ void reduce6(int *g_idata, int *g_odata, unsigned int n)
{
  extern __shared__ int sdata[];
  unsigned int tid = threadIdx.x;
  unsigned int i = blockIdx.x*(blockSize*2) + tid;
  unsigned int gridSize = blockSize*2*gridDim.x;
  sdata[tid] = 0;
  do { sdata[tid] += g_idata[i] + g_idata[i+blockSize]; i += gridSize; } while (i < n);
  __syncthreads();
  if (blockSize >= 512) { if (tid < 256) { sdata[tid] += sdata[tid + 256]; } __syncthreads(); }
  if (blockSize >= 256) { if (tid < 128) { sdata[tid] += sdata[tid + 128]; } __syncthreads(); }
  if (blockSize >= 128) { if (tid < 64) { sdata[tid] += sdata[tid + 64]; } __syncthreads(); }
  if (tid < 32) {
  if (blockSize >= 64) sdata[tid] += sdata[tid + 32];
  if (blockSize >= 32) sdata[tid] += sdata[tid + 16];
  if (blockSize >= 16) sdata[tid] += sdata[tid + 8];
  if (blockSize >= 8) sdata[tid] += sdata[tid + 4];
  if (blockSize >= 4) sdata[tid] += sdata[tid + 2];
  if (blockSize >= 2) sdata[tid] += sdata[tid + 1];
  }
  if (tid == 0) g_odata[blockIdx.x] = sdata[0];
}
 

Идея этой программы (кстати программа с багами) очень богатая, но не доведена до конца.
Тут сразу направшиваются для вычисления ряда параметров 512, 256, 128... приемы метапрограммирования.
Это я сразу осуществил в виде функции Exp2Log2. Оказалось что эта идея работает, тогда я обнаглел до конца и решил в Cuda сделать своего рода рекурсию (в Cuda делать рекурсию вообще-то для функций работающий в девайсе нельзя, на этот счет есть даже специальный строгий пассаж в NVIDIA_CUDA_Programming_Guide_1.1.pdf, но как выяснилось если сильно хочется -  то можно). Получилось так

//Ôóíêöèÿ f = 2^log2(N) äëÿ âû÷èñëåíèÿ êîíñòàíòíûõ ïàðàìåòðîâ âî âðåìÿ êîìïèëÿöèè
template<int N>
struct Exp2Log2{ enum { value = 2*Exp2Log2<N/2>::value };};
template<>
struct Exp2Log2<1> { enum { value = 1 };};

//Ôóíêöèÿ RecursivSum, ñóììèðóåò "ïî îëèìïèéñêîé ñèñòåìå" 
//âåêòîð sdata ðàçìåðîì blockSize
template<int blockSize, int blockSizeAligned>
__device__ void  RecursivSum(__shared__ myType *sdata)
{  
  if (threadIdx.x < blockSize-blockSizeAligned)
    sdata[threadIdx.x] += sdata[threadIdx.x + blockSizeAligned];
  __syncthreads();
  //ñîçäàåòñÿ íîâàÿ instance ôóíêöèè, ñ ïàðàìåòðàìè blockSize, blockSizeAligned/2:
  if (blockSizeAligned >= 2)  //óñëîâèå ïðîäîëæåíèÿ ðåêóðñèè
    RecursivSum<blockSize, blockSizeAligned/2>(sdata);  
};

//ñóììèðóåò âåêòîð g_idata ðàçìåðîì n, ðåçóëüòàò çàïèñûâàåò â g_odata
template <unsigned int blockSize>
__global__ void VectorSum(const myType *g_idata, myType *g_odata, const unsigned int n)
{
  __shared__ int sdata[blockSize];  
  sdata[threadIdx.x] = 0;
  for (unsigned int i = threadIdx.x; i < n; i += blockSize)
    sdata[threadIdx.x] += g_idata[i];  
  __syncthreads();
  RecursivSum<blockSize, Exp2Log2<blockSize>::value>(sdata);  
  if (threadIdx.x == 0) g_odata[0] = sdata[0];
}
 

Вызывать VectorNorm из основной программы нужно так

   dim3 threads(BLOCK_SIZE);
    dim3 grid(n / threads.x);  
  VectorSum <BLOCK_SIZE> <<< grid, threads >>>(g_idata, g_odata, n);
 

Эта функция работает довольно быстро, у меня вектор размером 4м она просуммировала за 0.338 мсек.
Естественно ее несложно переписать, что я и сделал, также для вычисления скалярного произведения векторов, нормы векторов и т.п. Но конечно основную ценность эта идея имеет при обработке сложных структур данных, типа деревьев. Там можно будет вместо того, чтоб расписывать каждый узел дерева просто сделать одну __device__ функцию и с помощью template генерировать при компиляции ее нужные instance.

    CUT_SAFE_CALL(cutCreateTimer(&timer));
    CUT_SAFE_CALL(cutStartTimer(timer));
  //......................
    CUT_SAFE_CALL(cutStopTimer(timer));
    CUT_SAFE_CALL(cutDeleteTimer(timer));

2. Из общего времени вызова  VectorSum в 0.338 мсек примерно 0.16 мсек стоит вызов ПУСТОЙ функции, т.е. функции которая вообще ничего не делает.
3. Изменение размерности задачи (т.е размера вектора) в СОТНИ раз время вызова функции VectorSum практически НЕ МЕНЯЕТ. Т.е. время собственно вычислений ничтожно в общем времени вызова функции  VectorSum.

Отсюда следуют важные выводы: следует по возможности ограничить число вызовов глобальных функций, в частности, в итерационных численных методах итерации нужно проводить в самой глобальной функции, может быть путем вызова девайс функций, а не проводить проверку сходимости итераций в хост функции.
Это между прочим объясняет плохие результыты использования cublac при итерационном решении СЛАУ в http://www.gamedev.ru/code/forum/?id=79845

 
__device__ float4 operator+(const float4 &a, const float4 &b) 
{return make_float4(a.x+b.x, a.y+b.y, a.z+b.z, a.w+b.w);}

__device__ float4 operator-(const float4 &a, const float4 &b) 
{return make_float4(a.x-b.x, a.y-b.y, a.z-b.z, a.w-b.w);}

__device__ float operator*(const float4 &a, const float4 &b) 
{return a.x*b.x+a.y*b.y+a.z*b.z+a.w*b.w;}

- В Cuda можно в структурах использовать функции-члены, включая конструкторы и деструкторы. Мне удавалось создавать структуры в шеред памяти с девайс функциями-членами.
- Однако в Cuda нельзя использовать наследование, идентификаторы public/private/protected и template <class T> или template <typename T>.

Все предыдущие утверждения относились глобал и девайс функциям и структурам в глобал и шеред памяти.

- Cuda программа в хосте может использовать все возможности c++ без ограничений, во всяком случае никакого примера противоречащего этому утверждению мне пока не удалось найти.

Последняя правка: 9 мая 2008 17:37

Ну целыми то типами функции параметризуются?

Да, это было бы интересно. Наверное кто-нибудь сделает FFT с шаблонной рекурсией (кстати, хороший термин). Но меня пока интересуют решения Слау и программа для обработки системы частиц. Для первого необходима шаблонная рекурсия при редукции данных, а для второго гораздо важнеее описание частиц как структур со встроенными функциями и, чисто для удобства и наглядности, с перегруженными операторами.

Да, это как раз и используется при шаблонной рекурсии.

>Наверное кто-нибудь сделает FFT с шаблонной рекурсией
Собственно, я и делал, только под DSP. Работало всего в 1.5 раза медленнее круто заоптимизированной и написанной на ассемблере FFT от Analog Devices.

Последняя правка: 10 мая 2008 21:51

http://ray-tracing.ru/articles162.html

мож кого заинтересует

Последняя правка: 11 мая 2008 12:41

>>#include <cutil.h>
и не должно быть. это же только темплейты
это мой Common Source Library, я его включаю в проект с CUDA и использую.

конечно надо было примеры с CUDA сделать, это я понимаю, но мне недосуг сейчас только.
любой код из библиотеки должен съедать nvcc - вот что я имею ввиду

>>v31 = m*v3; на v31 = v3*m;
а вот этого не должно быть
не мог бы ты скинуть сюда сообщение об ошибке и какой у тебя компилятор?

>>Кроме того, ты может быть не в курсе, но в CUDA недопустимы выражения типа template<int n,int m,class T>
что значит, недопустимы? это ты документацию читал что-ли?
а я утверждаю, что допустимы.
хочешь, сделаю проект с CUDA -  специально, где это работает

PS: я не понты кидаю, я подумал мож полезно кому будет.
моя либа как раз и реализует метапрограммирование. и компилится nvcc.

правка:
да, да, мой класс матриц прекрасно компилится под CUDA если ты про
template<int n,int m,class T>
class RectMatrix итд.
правда там надо еще #ifdef на конструкторы поставить, но эт. мелочи

Последняя правка: 11 мая 2008 22:56