Повышение эффективности Compact/Scatter в CUDA

Question

Резюме:

Любые идеи о том, как улучшить базовую операцию разброса в CUDA? Особенно, если знать, что он будет использоваться только для сжатия большего массива в меньший? или почему не работают приведенные ниже методы векторизации операций памяти и разделяемой памяти? Я чувствую, что может быть что-то фундаментальное, что я пропускаю, и любая помощь будет оценена.

РЕДАКТИРОВАТЬ 09/09/15: Таким образом, я нашел это сообщение в блоге " Параллельно для всех " "Оптимизированная фильтрация с использованием агрегированных по деформации атомик". Я предполагал, что атомы будут медленнее для этой цели, однако я был неправ - тем более, что я не думаю, что меня волнует поддержание порядка элементов в массиве во время моего моделирования. Я должен подумать об этом еще немного, а затем реализовать это, чтобы увидеть, что происходит!

РЕДАКТИРОВАТЬ 01/04/16: я понял, что я никогда не писал о своих результатах. К сожалению, в этом посте "Параллельно для всех" они сравнили глобальный атомарный метод для компактного метода с префиксной суммой Thrust, который на самом деле довольно медленный. Device:: IF CUB намного быстрее, чем Thrust - как и версия с префиксной суммой, которую я написал, используя собственный код CUB Device::Scan +. Глобальный атомарный метод, основанный на деформации, все еще быстрее, примерно на 5-10%, но нигде в 3-4 раза быстрее, на что я надеялся, основываясь на результатах в блоге. Я по-прежнему использую метод prefix-sum, так как при поддержании порядка элементов нет необходимости, я предпочитаю согласованность результатов prefix-sum, а преимущество от атомики не очень велико. Я все еще пробую различные методы для улучшения компактности, но пока только незначительные улучшения (2%) в лучшем случае для значительного увеличения сложности кода.

---

Подробности:

Я пишу симулятор в CUDA, где я сжимаю элементы, которые меня больше не интересуют, симулируя каждые 40-60 временных шагов. Из профилирования кажется, что операция сжатия занимает большую часть времени при сжатии - больше, чем ядро ​​фильтра или сумма префикса. Прямо сейчас я использую довольно простую функцию разброса:

    __global__ void scatter_arrays(float * new_freq, const float * const freq, const int * const flag, const int * const scan_Index, const int freq_Index){
            int myID =  blockIdx.x*blockDim.x + threadIdx.x;
            for(int id = myID; id < freq_Index; id+= blockDim.x*gridDim.x){
                 if(flag[id]){
                    new_freq[scan_Index[id]] = freq[id];
                 }
             } 
    }

freq_Index - количество элементов в старом массиве. Массив флагов является результатом фильтра. Scan_ID является результатом суммы префикса в массиве флагов.

Попытки улучшить его состоят в том, чтобы сначала прочитать отмеченные частоты в разделяемой памяти, а затем записать из разделяемой памяти в глобальную память - идея состоит в том, что записи в глобальную память будут более объединены между деформациями (например, вместо потока 0). запись в позицию 0 и поток 128 запись в позицию 1, поток 0 записывает в 0, а поток 1 записывает в 1). Я также попытался векторизовать операции чтения и записи - вместо чтения и записи чисел с плавающей точкой / целых я читал / записывал файлы float4/int4 из глобальных массивов, когда это возможно, поэтому по четыре числа за раз. Я подумал, что это может ускорить рассеяние, если меньше операций памяти будет передавать больший объем памяти. Код "кухонной раковины" с загрузкой / хранением векторизованной памяти и общей памятью приведен ниже:

    const int compact_threads = 256;
    __global__ void scatter_arrays2(float * new_freq, const float * const freq, const int * const flag, const int * const scan_Index, const int freq_Index){
        int gID =  blockIdx.x*blockDim.x + threadIdx.x; //global ID
        int tID = threadIdx.x; //thread ID within block
        __shared__ float row[4*compact_threads];
        __shared__ int start_index[1];
        __shared__ int end_index[1];
        float4 myResult;
        int st_index;
        int4 myFlag;
        int4 index;
        for(int id = gID; id < freq_Index/4; id+= blockDim.x*gridDim.x){
            if(tID == 0){
                index = reinterpret_cast<const int4*>(scan_Index)[id];
                myFlag = reinterpret_cast<const int4*>(flag)[id];
                start_index[0] = index.x;
                st_index = index.x;
                myResult = reinterpret_cast<const float4*>(freq)[id];
                if(myFlag.x){ row[0] = myResult.x; }
                if(myFlag.y){ row[index.y-st_index] = myResult.y; }
                if(myFlag.z){ row[index.z-st_index] = myResult.z; }
                if(myFlag.w){ row[index.w-st_index] = myResult.w; }
            }
            __syncthreads();
            if(tID > 0){
                myFlag = reinterpret_cast<const int4*>(flag)[id];
                st_index = start_index[0];
                index = reinterpret_cast<const int4*>(scan_Index)[id];
                myResult = reinterpret_cast<const float4*>(freq)[id];
                if(myFlag.x){ row[index.x-st_index] = myResult.x; }
                if(myFlag.y){ row[index.y-st_index] = myResult.y; }
                if(myFlag.z){ row[index.z-st_index] = myResult.z; }
                if(myFlag.w){ row[index.w-st_index] = myResult.w; }
                if(tID == blockDim.x -1 || gID == mutations_Index/4 - 1){ end_index[0] = index.w + myFlag.w; }
            }
            __syncthreads();
            int count = end_index[0] - st_index;

            int rem = st_index & 0x3; //equivalent to modulo 4
            int offset = 0;
            if(rem){ offset = 4 - rem; }

            if(tID < offset && tID < count){
                new_mutations_freq[population*new_array_Length+st_index+tID] = row[tID];
            }

            int tempID = 4*tID+offset;
            if((tempID+3) < count){
                reinterpret_cast<float4*>(new_freq)[tID] = make_float4(row[tempID],row[tempID+1],row[tempID+2],row[tempID+3]);
            }

            tempID = tID + offset + (count-offset)/4*4;
            if(tempID < count){ new_freq[st_index+tempID] = row[tempID]; }
        }
        int id = gID + freq_Index/4 * 4; 
        if(id < freq_Index){
            if(flag[id]){
                new_freq[scan_Index[id]] = freq[id];
            }
        }
    }

Очевидно, это становится немного сложнее.:) Хотя вышеприведенное ядро ​​кажется стабильным, когда в массиве сотни тысяч элементов, я заметил состояние гонки, когда массив насчитывает десятки миллионов. Я все еще пытаюсь отследить ошибку.

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

Answer 1

Я обнаружил, что алгоритм, упомянутый в этом постере (аналогичный алгоритм, также обсуждаемый в этой статье), работает довольно хорошо, особенно для сжатия больших массивов. Он использует меньше памяти, чтобы сделать это, и немного быстрее, чем мой предыдущий метод (5-10%). Я внес несколько изменений в алгоритм плаката: 1) устранение окончательного уменьшения тасования основы в фазе 1, можно просто суммировать элементы по мере их вычисления, 2) дать функции возможность работать не только с массивами, измеренными как кратное 1024 + добавление циклов с шагом сетки и 3) разрешение каждому потоку загружать свои регистры одновременно в фазе 3, а не по одному за раз. Я также использую CUB вместо суммы Thrust for Inclusive для ускорения сканирования. Я могу сделать больше твиков, но пока это хорошо.

//kernel phase 1
int myID =  blockIdx.x*blockDim.x + threadIdx.x;
//padded_length is nearest multiple of 1024 > true_length
for(int id = myID; id < (padded_length >> 5); id+= blockDim.x*gridDim.x){
    int lnID = threadIdx.x % warp_size;
    int warpID = id >> 5;

    unsigned int mask;
    unsigned int cnt=0;//;//

    for(int j = 0; j < 32; j++){
        int index = (warpID<<10)+(j<<5)+lnID;

        bool pred;
        if(index > true_length) pred = false;
        else pred = predicate(input[index]);
        mask = __ballot(pred); 

        if(lnID == 0) {
            flag[(warpID<<5)+j] = mask;
            cnt += __popc(mask);
        }
    }

    if(lnID == 0) counter[warpID] = cnt; //store sum
}

//kernel phase 2 -> CUB Inclusive sum transforms counter array to scan_Index array

//kernel phase 3
int myID =  blockIdx.x*blockDim.x + threadIdx.x;

for(int id = myID; id < (padded_length >> 5); id+= blockDim.x*gridDim.x){
    int lnID = threadIdx.x % warp_size;
    int warpID = id >> 5;

    unsigned int predmask;
    unsigned int cnt;

    predmask = flag[(warpID<<5)+lnID];
    cnt = __popc(predmask);

    //parallel prefix sum
#pragma unroll
    for(int offset = 1; offset < 32; offset<<=1){
        unsigned int n = __shfl_up(cnt, offset);
        if(lnID >= offset) cnt += n;
    }

    unsigned int global_index = 0;
    if(warpID > 0) global_index = scan_Index[warpID - 1];

    for(int i = 0; i < 32; i++){
        unsigned int mask = __shfl(predmask, i); //broadcast from thread i
        unsigned int sub_group_index = 0;
        if(i > 0) sub_group_index = __shfl(cnt, i-1);
        if(mask & (1 << lnID)){
            compacted_array[global_index + sub_group_index + __popc(mask & ((1 << lnID) - 1))] = input[(warpID<<10)+(i<<5)+lnID]; 
        }
    }
}

}

РЕДАКТИРОВАТЬ: есть более новая статья подмножеством авторов плаката, где они исследуют более быстрое изменение компакта, чем написанное выше. Тем не менее, их новая версия не сохраняет порядок, поэтому для меня это бесполезно, и я не реализовал ее для тестирования. Тем не менее, если ваш проект не зависит от порядка объектов, их более новая компактная версия, вероятно, может ускорить ваш алгоритм.