Обратное преобразование.

Прямое преобразование

Выполняется сверка целевого сигнала с низкочастотным и высокочастотными фильтрами с одновременным прореживанием, что при последовательном выполнении дает трудоемкость кратную N/2.

Но так как данную операцию на исходной последовательности данных нужно проводить Log₂(N) раз, каждый раз уменьшая размерность на 2, то получаем трудоемкость:

 

Трешолдинг

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

Обратное преобразование.

Обратное преобразование производится аналогично прямому и имеет ту же трудоемкость:

При использовании любых последовательных средств вычислений мы сталкиваемся с жесткой зависимостью от объема входных данных, выраженную в трудоемкости:

 

1.3 Фильтрация на GPGPU

 

GPGPU – сопроцессор к центральному процессору для массово параллельных вычислений. Теоретически при использовании данного устройства для алгоритмов, допускающих мелкозернистый параллелизм, будет отсутствовать сильная связь с объемом входных данных. То есть при любом объеме мы будем получать результат за ограниченное время.

Перечислим возможности GPGPU как массивно параллельного устройства вычисления:

· множество процессоров виртуально может быть представлено в любой конфигурации с тремя с измерениями;

· каждый процессор имеет уникальный идентификатор;

· процессоры разделены на сильно связные компоненты – потоки в блоках, и слабосвязные – блоки в сетке;

· для оптимизации работы реализована иерархия памяти, технологии группового чтения данных, кэшированного чтения данных;

· одновременно можно запустить на выполнение множество тысяч потоков.

 

Представим каждую из трех частей в параллельном виде, доступном для запуска на массивно параллельном вычислителе.

Прямое вейвлет преобразование.

При прямом вейвлет преобразовании выполняется свертка и для каждого значения ряда сигнала вычисляется новое значение в соответствии с массивами коэффициентов.

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

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

Что и выполнено ниже:

__global__ void

 cudaDWTStepForward(float *a, const size_t n, const size_t nmod, const size_t stride, const float * const h1, const float* const g1, const size_t nc, float*scratch)

{

const size_t idx = blockDim.x * blockIdx.x + tidx; // индекс элемента данных в первой половине = РазмерБлока*номерБлока + индексПотокаВБлоке

 

const size_t n1 = n - 1;

const size_t nh = n >> 1;

 

if(idx < nh){

   float h = 0, g = 0;

   const size_t ni = 2*idx + nmod;

 

 

   for (size_t k = 0; k < nc; k++) // работа с коэффициентами

     {

       const size_t jf = n1 & (ni + k);

       const float ani = a[stride * jf];

       h += h1[k] * ani;

       g += g1[k] * ani;

     }

         

   scratch[idx] += h;  // Преобразованные данные

   scratch[idx + nh] += g;

}

Также необходимо произвести расчет числа процессоров под задачу прямого преобразования исходя из формулы:

Число процессоров = числу входных элементов.

Вызов функции на выполнение:

// параметры конфигурации исполнения

const size_t block_size = BL_SZ;

dim3 threads( block_size, 1);

const size_t grid_size = (n/2) / (block_size - nc) + (((n/2)%(block_size - nc))>0); // n/2 потому проходится лишь первая половина

dim3 grid(grid_size, 1);

 

// запуск ядра на видео-карте

cudaDWTStepForward<<< grid, threads >>>(d_a, n, nmod/*, stride*/, d_h1, d_g1, nc, d_scratch);

 

Теоретическая трудоемкость данного кода будет приближаться к O(Log2(N)), так как необходимо последовательно выполнить данную функцию для прямого преобразования в соответствии с алгоритмом, описанном выше.

 

Пороговая обработка

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

 

__global__ void

 cudaDenoiseKernel(float *data, const size_t n,const float threshold)

{

const size_t idx = blockDim.x * blockIdx.x + threadIdx.x; // индекс элемента данных в первой половине = РазмерБлока*номерБлока + индексПотокаВБлоке

if(idx < n){

float d = data[idx];

const float fd = fabsf(d);

float t = fd - threshold;

 t = (t + fabsf(t)) / 2.f;

 

 // Signum

 if (d != 0.f) d = d / fd * t;

 data[idx] = d;

}

}

Для успешного выполнения обработки необходимо жестко контролировать нахождение алгоритма внутри допустимого множества вычисляемых данных. Контроль необходимо осуществлять как при проектировании ядра, так и во время вычисления количества и конфигурации потоков для запуска ядер:

 

// параметры конфигурации исполнения

const size_t block_size = 16;

dim3 threads( block_size, 1);

const size_t grid_size = len / (block_size) + ((len%(block_size))>0); // n/2 потому проходится лишь первая половина

dim3 grid(grid_size, 1);

 

// запуск ядра на видео-карте

cudaDenoiseKernel<<< grid, threads >>>(d_a, len, threshold);

 

Обратное вейвлет преобразование.

Обратное вейвлет преобразование мало отличается от прямого, разве что формулой вычисления.

Также для каждого входного элемента можно создать отдельный поток.

Рассмотренные выше методы распараллеливания довольно тривиальны и очевидны, но, к сожалению, использование только лишь здравого смысла как показывает практика - не дает роста производительности по сравнению с реализацией на универсальном вычислительном устройстве. А дает лишь набор проблем при работе с модулем и неточности результатов.

 

1.4 Проблемы распараллеливания

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

Приведем проблемы, которые необходимо преодолеть для получения эффективной программы.

· Проблема точности вычислений: GPGPU имеют короткую историю и фактически начали активно развиваться лишь с 2007-го года: в настоящее время поддержка плавающих чисел двойной точности всё еще является предметом обсуждения и встроена не во все имеющиеся GPGPU; также следует отметить потерю производительности при использовании плавающих чисел двойной точности.

· Проблема балансировки загрузки памяти, процессоров, шин памяти для эффективного использования устройства – что в конечном итоге сказывается на производительности.

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

· После проведенных работ по оптимизации поддержка кода становится невозможной.

 

 

1.5 Методы оптимизации

 

Работа с памятью

 

Глобальная память

Глобальная память (GPU DRAM) не кэшируется и поэтому одно обращение к ней может занимать до 600 тактов. Именно оптимизация работы с глобальной памятью может дать огромный прирост производительности.

Обращение к глобальной памяти осуществляется через чтение/запись 32/64/128-битовых слов.

__device__ type dev [32];             // 32 items of type in global memorytype        data = dev [threadIdx.x]; // reading one item from global memory

 

Для того, чтобы обращение в глобальную память было скомпилировано в одну команду, важно, чтобы sizeof(type) был равен 4/8/16 байтам и расположение данных было выровнено по sizeof(type) в памяти.

__device__ float data [256];float3      v; . . . . v = data [i];

 

Рассмотрим пример выше - размер float3 равен 12 байтам и, поэтому, многие элементы массива data не будут надлежащим образом выровнены. Все это ведет к неэффективному обращению к глобальной памяти. Избежать этого можно несколькими способами.

Простейшим является замена float3 на float4 - хотя мы при этом будем использовать немного больше памяти, но размер каждого элемента массива будет равен 16 байтам и каждый элемент будет выровнен в памяти по 16 байтам. Это позволит скомпилировать каждое обращение к элементу такого массива в одну команду.

Точно такая же ситуация возникает при использовании структур. Если размер структуры превышает 16 байт, то будет произведено несколько обращений к памяти. Для наибольшей эффективности желательно, чтобы размер структуры был кратен 4/8/16 байтам и все экземпляры данной структуры были выровнены в памяти. Так приводимый ниже пример использует опцию __align__ для обеспечения выравнивания в памяти по 16 байтам (что позволит минимизировать число операций чтения/записи).

struct __align__(16){ float a; float b; float c;};

 

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

Крайне важной особенностью GPU является возможность объединения (coalescing global memory accesses) нескольких запросов к глобальной памяти в одну операцию над блоком (транзакцию). Подобное объединение запросов в один запрос чтения/записи одного блока длиной 32/64/128 байт может происходить в пределах одного half-warp'а. Для того, чтобы такое объединение произошло, должны быть выполнены специальные требования на то, как отдельные нити half-warp'а обращаются к памяти.

Адрес блока должен быть выровнен в памяти. Также должны быть выполнены дополнительные требования, зависящие от Compute Capability GPU.

Для GPU с Compute Capability 1.0 и 1.1 объединения запросов в одну транзакцию будет происходит при выполнении следующих условий:

· нити должны обращаться либо к 32-битовым словам, давая при этом в результате один 64-байтовый блок (транзакцию), либо к 64-битовым словам, давая при этом один 128-байтовый блок (транзакцию);

· все 16 слов должны лежать в пределах данного блока (64 или 128 байт)

· нити должны обращаться к словам последовательно, т.е. k-ая нить должна обращаться к k-му слову (при этом допускается, что какие-то нити вообще не производят обращения к соответствующим им словам).

Если нити half-warp'а не удовлетворяют какому-либо из данных условий, то каждое обращение к памяти происходит как отдельная транзакция. На следующих рисунках приводятся типичные паттерны обращения, дающие объединения и не дающие объединения.

Рисунок 7 - Паттерны обращения к памяти, дающие объединение для Compute Capability 1.0 и 1.1.

На рисунке 7 приведены типичные паттерны обращения к памяти, приводящие к объединению запросов в одну транзакцию. Слева у нас выполнены все условия, справа - просто для части нитей пропущено обращение к соответствующим словам (что равно позволяет добавить фиктивные обращения и свести к случаю слева).

Рисунок 8 - Паттерны обращения к памяти, не дающие объединение для Compute Capability 1.0 и 1.1.

На рисунке 8 слева для нитей 4 и 5 нарушен порядок обращения к словам, а справа нарушено условие выравнивание - хотя слова, к которым идет обращения и образуют непрерывный блок из 64 байт, но начало этого блока (по адресу 132) не кратно его размеру.

Для GPU с Compute Capability 1.2 и выше объединения запросов в один будет происходить, если слова, к которым идет обращение нитей, лежат в одном сегменте размера 32 байта (если все нити обращаются к 8-битовым словам), 64 байта (если все нити обращаются к 16-битовым словам) и 128 байт (если все нити обращаются к 32-битовым или 64-битовым словам). Получающийся сегмент (блок) должен быть выровнен по 32/64/128 байтам.

Обратите внимание, что в этом случае порядок, в котором нити обращаются к словам, не играет никакой роли и ситуация на рисунке 8 слева приведет к объединению всех запросов в одну транзакцию.

Если идет обращения к n соответствующим сегментам, то происходит группировка запросов в n транзакций (только для GPU с Compute Capability 1.2 и выше).

Еще одним моментом работы с глобальной памятью является то, что гораздо эффективнее с точки зрения объединения запросов к памяти является использование не массивов структур, а структуры массивов (отдельных компонент исходной структуры).

struct A __align__(16){ float a; float b; uint c;}; A array [1024]; . . .A a = array [threadIdx.x];

 

В приведенном выше примере чтение структуры каждой нитью не даст объединения обращений к памяти и на доступ к каждому элементу массива A понадобится отдельная транзакция.

Вместо одного массива A можно сделать три массива его компонент.

float a [1024];float b [1024];uint c [1024];. . .float fa = a [threadIdx.x];float fb = b [threadIdx.x];uint uc = c [threadIdx.x];

 

В результате такого разбиения каждый из трех запросов к очередной компоненте исходной структуры приведет к объединению и в результате нам понадобится всего по три транзакции на half-warp (вместо 16 ранее).

Еще одним важным моментом при работе с глобальной памятью является то, что при большом числе нитей, выполняемых на мультипроцессоре (мультипроцессор поддерживает до 768 нитей), время ожидания warp'ом доступа к памяти может быть использовано для выполнения других warp'ов. Чередование вычислений с обращениями к памяти позволяет более оптимально использовать ресурсы GPU.

Так первому warp'у нужен доступ к памяти. Управление передается другому warp'у, выполняется одна команда для него, далее управление опять передается следующему warp'у и т.д. Для того, чтобы избежать "простоя" мультипроцессора, достаточно обеспечить большое количество warp'ов, которые смогут выполняться в то время, когда первый warp ждет.

Для оптимального доступа к памяти и регистрам желательно, чтобы количество нитей в блоке было кратным 64 и было не менее 192.

Работы с shared-памятью

Для повышения пропускной способности вся shared-память разбита на 16 банков. Каждый банк работает независимо и если идет обращение к набору слов, лежащих в разных банках, то все эти запросы будут обработаны одновременно.

Если же сразу несколько нитей обращаются к одному и тому же банку памяти, то происходит конфликт банков (bank conflict) и эти запросы будут выполняться последовательно один за другим (т.е. в несколько раз медленнее).

Банки памяти организованы следующим образом - последовательно идущие 32-битовые слова попадают в последовательно идущие банки. GPU при обращении к shared-памяти независимо обрабатывает запросы от нитей из первого и от второго half-warp'ов. Таким образом, нити, принадлежащие разным half-warp'ам не могут вызывать конфликт банков - конфликтовать могут только нити в пределах каждого half-warp'а между собой.

Рассмотрим типичный пример, когда адрес для обращения к памяти линейно зависит от номера нити:

__shared__ float buf [128];float       v = buf [baseIndex + threadIdx.x * s];

 

При такой схеме обращения к shared-памяти конфликтов не будет только тогда, когда s нечетно. Однако обратите внимание, что если идет доступ к элементам меньшего размера, то ситуация изменяется:

__shared__ char buf [128];char       v = buf [baseIndex + threadIdx.x];

 

Как легко заметить buf[0], buf[1], buf[2] и buf[3] лежат в одном и том же банке. Поэтому в приведенном выше случае будет иметь место 4-кратный конфликт. Однако его легко можно избежать, поменяв схему обращения к памяти:

__shared__ char buf [128];char       v = buf [baseIndex + threadIdx.x * 4];

 

На следующем рисунке приведены типичные способы обращения нитей к shared-памяти, когда конфликтов банков не происходит - фактически это означает что существует взаимно-однозначное соответствие между нитями half-warp'а и банками памяти.

Рисунок 9 - Паттерны безконфликтного доступа к банкам shared-памяти.

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

Рисунок 10 - Паттерны доступа, приводящие к конфликтам.

На рисунке 10 приведены два способа обращения к shared-памяти, приводящие к возникновению конфликтов.

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

При этом сама загрузка (когда каждая нить загружает по одному элементу) обычно приводит к объединению запросов в глобальную память. Однако при обработке полученной небольшой матрицы в shared-памяти следует иметь в виду, что если ее размер кратен 16, то каждый ее столбец размещается в одном банке. Поэтому если каждая нить по очереди перебирает элементы, сначала из первого столбца, затем из второго и т.д., то мы сразу получаем конфликт банков, причем порядка 16.

Самым простым способом избежать этого является дополнение матрицы еще одним столбцом (т.е. от матрицы из 16 столбцов мы переходим к матрице из 17 столбцов). За счет этого расположение столбцов в памяти сдвигается, и каждый столбец теперь равномерно распределяется по всем банкам памяти.

Рисунок 11 - Дополнение матрицы.

Вообще подобный прием - добавление к каждой строке матрицы нескольких дополнительных байт довольно распространен в CUDA - за счет него при доступе к глобальной памяти мы легко гарантируем выравнивание для начала строки, а при доступе к shared-памяти - равномерное распределение столбца по банкам памяти.

Локальная память

Когда нитям блока не хватает имеющихся на мультипроцессоре регистров, то для размещения локальных переменных используется локальная память GPU. Локальная память размещена в DRAM и не кэшируется, поэтому доступ к ней такой же медленный как и к глобальной памяти. Однако компилятор выполняет запросы отдельных нитей к локальной памяти, объединяя в одну транзакцию (coalescing).