Как оптимизировать это ядро ​​CUDA

Question

Я профилировал мою модель, и кажется, что это ядро ​​составляет около 2/3 от моей общей продолжительности. Я искал предложения по его оптимизации. Код выглядит следующим образом.

__global__ void calcFlux(double* concs, double* fluxes, double* dt) 
{ 
    int idx = blockIdx.x * blockDim.x + threadIdx.x; 
    fluxes[idx]=knowles_flux(idx, concs); 
    //fluxes[idx]=flux(idx, concs); 
} 

__device__ double knowles_flux(int r, double *conc) 
{ 
    double frag_term = 0; 
    double flux = 0; 
    if (r == ((maxlength)-1)) 
    { 
     //Calculation type : "Max" 
     flux = -km*(r)*conc[r]+2*(ka)*conc[r-1]*conc[0]; 
    } 
    else if (r > ((nc)-1)) 
    { 
     //Calculation type : "F" 
     //arrSum3(conc, &frag_term, r+1, maxlength-1); 
     for (int s = r+1; s < (maxlength); s++) 
     { 
      frag_term += conc[s]; 
     } 
     flux = -(km)*(r)*conc[r] + 2*(km)*frag_term - 2*(ka)*conc[r]*conc[0] + 2*(ka)*conc[r-1]*conc[0]; 
    } 
    else if (r == ((nc)-1)) 
    { 
     //Calculation type : "N" 
     //arrSum3(conc, &frag_term, r+1, maxlength-1); 
     for (int s = r+1; s < (maxlength); s++) 
     { 
      frag_term += conc[s]; 
     } 
     flux = (kn)*pow(conc[0],(nc)) + 2*(km)*frag_term - 2*(ka)*conc[r]*conc[0]; 
    } 
    else if (r < ((nc)-1)) 
    { 
    //Calculation type : "O" 
     flux = 0; 
    } 
    return flux; 
} 

Просто чтобы дать вам представление о том, почему цикл является проблемой, это ядро ​​запускается на массиве около MAXLENGTH = 9000 элементов. В наших целях теперь nc находится в диапазоне 2-6. Вот иллюстрация того, как это ядро ​​обрабатывает входящий массив (conc). Для этого массива необходимо использовать пять разных типов вычислений для разных групп элементов.

Array element : 0 1 2 3 4 5 6 7 8 9 ... 8955 8956 8957 8958 8959 8960 
Type of calc : M O O O O O N F F F ... F F F F F Max 

Потенциальные проблемы я пытался бороться с прямо сейчас ветвью дивергенция от четверки, если-то еще и для цикла.

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

Чтобы справиться с циклом for, вы заметите, что есть прокомментированная функция arrSum3, которую я написал, основываясь на моем ранее (и, вероятно, плохо) написанном параллельном ядре сокращения. Использование его вместо цикла for резко увеличило время выполнения. Я чувствую, что есть умный способ добиться того, что я пытаюсь сделать с циклом for, но я просто не настолько умный, и мой советник устал от того, что я «теряю время», думая об этом.

Цените любую помощь.

EDIT

Полный код находится здесь: https://stackoverflow.com/q/21170233/1218689

Answer 1

Предполагая SGN() и абс() не являются производными от "если" с и "еще" s

__device__ double knowles_flux(int r, double *conc) 
{ 
    double frag_term = 0; 
    double flux = 0; 

     //Calculation type : "Max" 
     //no divergence 
     //should prefer 20-30 extra cycles instead of a branching. 
     //may not be good for CPU 
     fluxA = (1-abs(sgn(r-(maxlength-1)))) * (-km*(r)*conc[r]+2*(ka)*conc[r-1]*conc[0]); 
     //is zero if r and maxlength-1 are not equal 

     //always compute this in shared memory so work will be equal for all cores, no divergence 

     // you should divide kernel into several pieces to do a reduction 
     // but if you dont want that, then you can try : 
     for (int s = 0;s<someLimit ; s++) // all count for same number of cycles so no divergence 
     { 
      frag_term += conc[s] * ( abs(sgn(s-maxlength))*sgn(1- sgn(s-maxlength)) )* (  sgn(1+sgn(s-(r+1))) ); 
     } 
     //but you can make easier of this using "add and assign" operation 
     // in local memory (was it __shared in CUDA?) 
     // global conc[] to local concL[] memory(using all cores)(100 cycles) 
     // for(others from zero to upper_limit) 
     // if(localID==0) 
     // { 
     // frag_termL[0]+=concL[s]    // local to local (10 cycles/assign.) 
     // frag_termL[0+others]=frag_termL[0]; // local to local (10 cycles/assign.) 
     // } -----> uses nearly same number of cycles but uses much less energy 
     //using single core (2000 instr. with single core vs 1000 instr. with 2k cores) 
     // in local memory, then copy it to private registers accordingly using all cores 



     //Calculation type : "F" 

     fluxB = ( abs(sgn(r-(nc-1)))*sgn(1+sgn(r-(nc-1))) )*(-(km)*(r)*conc[r] + 2*(km)*frag_term - 2*(ka)*conc[r]*conc[0] + 2*(ka)*conc[r-1]*conc[0]); 
     // is zero if r is not greater than (nc-1) 



     //Calculation type : "N" 


     fluxC = ( 1-abs(sgn(r-(nc-1))) )*((kn)*pow(conc[0],(nc)) + 2*(km)*frag_term - 2*(ka)*conc[r]*conc[0]); 
     //zero if r and nc-1 are not equal 



    flux=fluxA+fluxB+fluxC; //only one of these can be different than zero 

    flux=flux*( -sgn(r-(nc-1))*sgn(1-sgn(r-(nc-1))) ) 
    //zero if r > (nc-1) 

    return flux; 
} 

Хорошо, позвольте мне немного открыть:

if(a>b) x+=y; 

могут быть приняты в качестве

if a-b is negative sgn(a-b) is -1 
then adding 1 to that -1 gives zero ==> satisfies lower part of comparison(a<b) 
x+= (sgn(a-b) +1) = 0 if a<b (not a>b), x unchanged 

if(a-b) is zero, sgn(a-b) is zero 
then we should multiply the upper solution with sgn(a-b) too! 
x+= y*(sgn(a-b) +1)*sgn(a-b) 
means 
x+= y*(0 + 1) * 0 = 0   a==b is satisfied too! 

lets check what happens if a>b 
x+= y*(sgn(a-b) +1)*sgn(a-b) 
x+= y*(1 +1)*1 ==> y*2 is not acceptable, needs another sgn on outherside 

x+= y* sgn((sgn(a-b)+1)*sgn(a-b)) 

x+= y* sgn((1+1)*1) 

x+= y* sgn(2) 

x+= y only when a is greater than b 

, когда есть слишком много

abs(sgn(r-(nc-1)) 

, то вы можете повторно использовать его в качестве

tmp=abs(sgn(r-(nc-1)) 

..... *tmp*(tmp-1) .... 
...... +tmp*zxc[s] ..... 
...... ...... 

уменьшить общий цикл еще больше! Регистрация доступа может быть на уровне терабайт/с, поэтому не должна быть проблемой. Подобно тому, как это делается для глобального доступа:

tmpGlobal= conc[r]; 

...... tmpGlobal * tmp ..... 
.... tmpGlobal +x -y .... 

все частные регистры занимаются вещами в терабайтах в секунду.

Предупреждение: Чтение от conc [-1] не должно вызывать никаких ошибок, если оно умножается на ноль, если реальный адрес conc [0] уже не является нулевым. Но написание является опасным.

Если вам нужно убежать от conc [-1], в любом случае вы можете умножить индекс с некоторым абсолютным значением! См:

tmp=conc[i-1] becomes tmp=conc[abs((i-1))] will always read from positive index, the value will be multiplied by zero later anyway. This was lower bound protection. 
    You can apply a higher bound protection too. Just this adds even more cycles. 

Подумайте об использовании операции векторного воспроизведения в случайном порядке при работе на чистом скалярные значения не достаточно быстро при доступе к конц [г-1] и конц [г + 1]. Операция в случайном порядке между элементами вектора быстрее, чем копирование его через локальный mem на другой ядро ​​/ поток.