Реорганизация вычислений

 

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

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

Рассмотрим специфические команды архитектуры. Выполнение следующих действий требует 4 тактов мультипроцессора:

· операции сложения, умножения и madd для типа float;

· целочисленное сложение;

· побитовые операции, сравнение, min, max;

· операции преобразования типов.

Функции 1/x, 1/sqrt(x) и __logf(x) выполняются за 16 тактов.

Целочисленное умножение (32-битовое) также требует 16 тактов, однако во многих случаях вместо него можно воспользоваться функциями __mul24(x,y) и __umul24(x,y), которые занимают всего 4 такта, но осуществляют 24-битовое умножение (т.е. если числа не очень большие - старший байт равен нулю - то эти функции точно заменяют умножение).

Крайне дорогостоящими операциями, которые желательно избегать, оказывается целочисленное деление (x/y) и целочисленный остаток от деления (x%y). В ряде случаев их можно заменять побитовыми операциями, в случае когда идет деление на константу (являющуюся степенью двух) или взятие остатка от деления на константу (вида 2ⁿ-1), то компилятор сам заменяет операцию на побитовую.

Квадратный корень реализуется как суперпозиция двух операций 1/x и 1/sqrt(x) и поэтому занимает 32 такта.

Деление на значение типа float требует 36 тактов, однако есть функция __fdividef(x,y), осуществляющая деление за 20 тактов (но возвращающая ноль при очень больших значениях y).

Функции __sinf(x), __cosf(x) и __expf(x) требуют 32 тактов, функции sinf(x), cosf(x), sincos(x) и tanf(x) гораздо более дорогостоящи.

Кроме того компилятор может вставлять команды для приведения типа в следующих случаях:

· на вход функции, принимающий int (или usnigned int) передается char или short;

· использование констант типа double внутри выражения, использующего float;

· передаче float-значения на вход функции, требующей double (например sin).

Поэтому рекомендуется явно указывать тип float, например вместо 1.3 писать 1.3f и использовать float-аналоги функций (например __sinf или sinf вместо sin).

Обращение к памяти занимает 4 такта, при обращении к глобальной памяти следует также добавить еще от 400 до 600 тактов.

Вызов __syncthreads() требует 4 тактов, если не нужно ждать ни одной нити в warp'е.

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

Как альтернативу данному методу можно пользоваться опцией компилятора nvcc --use_fast_math | -use_fast_math (использовать библиотеку быстрых, но не точных вычислений), которая автоматически при выполнении всех математических операций заменяет их на более быстрые аналоги, что дает прирост производительности, но и все связанные с этим недостатки. Мы предлагаем пользоваться данным методом при планировании ожидаемого выигрыша от реорганизации вычислений.

 

1.6 Исследование производительности

 

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

На рисунке ниже приведен результат преобразования на основе GPGPU, как можно заметить мы достигли некоторого шумоподавления обусловленного применением вейвлета Добеши 8.

Рисунок 12 – Вейвлет фильтрация сейсмического сигнала

Для исследования производительности немного провести серию тестов на разных объемах входных данных. Заметим, что трудоемкость работы метода не зависит от характера входных данных, а только от их объема. Чтобы сориентироваться в достигнутом на данный момент быстродействии вейвлет преобразования будем тестировать одновременно и последовательный метод вейвлет преобразования для CPU, взятый из промышленной свободно распространяемой библиотеки gsl. В таблицах 3 и 4 приведены конфигурации тестовой системы.

В приложении В отражен запуск теста вейвлет преобразования для разных объемов входных данных. В таблице ниже и на графике приведены времена преобразований для CPU и для GPU.

Таблица 2 — Время на Вейвлет фильтрацию

Количество отсчетов	CPU преобразование, мс	CUDA преобразование, мс	CUDA в сравнении с CPU, раз
65536	12	67	0,18
131072	26	67	0,39
262144	81	61	1,33
524288	159	73	2,18
1048576	319	85	3,75
2097152	677	97	6,98
4194304	1401	133	10,53
8388608	2855	209	13,66
16777216	5964	355	16,80
33554432	11979	638	18,78
67108864	31802	1350	23,56

 

Как и ожидалось, время на процессоре растет линейно при росте объемов входных данных, а на GPGPU близка к константе.

Рисунок 13 — Время вейвлет преобразования сигнала

Таблица 3 ‑ Аппаратная часть тестовой платформы

Центральный процессор:	Intel Core 2 Duo E2200 2.2 ГГц 2 Мб Cache; частота шины 800 МГц;
Оперативная память:	2048 Мб DDR3 800 МГц;
Чипсет:	Intel G33;
Видеоадаптер:	Дискретный стоковый адаптер NVidia Corporation; NVIDIA GeForce 9800 GT 1024 Мб GDDR-3;

 

Таблица 4 ‑ Программная часть тестовой платформы

Операционная система	Ubuntu Linux 10.04 x86 32бит;
Драйвер	NVIDIA Driver for Linux 195.35.08 ;
Версия CUDA SDK	CUDA SDK 3.0 for Linux
Версия CUDA toolkit	CUDA Toolkit 3.0 for Linux

 

Заключение

 

В результате выполнения данной работы были достигнуты поставленные цели.

· Разработан модуль вейвлет преобразования на основе GPGPU.

· Исследована производительность модуля, время выполнения на GPGPU NVIDIA 9800 GT в сравнении с CPU Intel Core 2 DUO E2200 в 23 раза меньше.

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

· Модуль внедрен в действующий программный комплекс по анализу сейсмических сигналов Seismo Detector.

· Исследована практическая эффективность Cuda технологии для решения многопроходных задач, допускающих мелкозернистый параллелизм.

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

 

Список использованных источников

 

1. Хайретдинов М.С., Клименко С.М. Программная система автоматизированной локации и визуализации сейсмических источников. Труды III Межд. конф. "Мониторинг ядерных испытаний и их последствий", Боровое, 2004.–Вестник НЯЦ РК, С.70–76

2. NVIDIA CUDA Programming Guide / NVIDIA Corporation, 2008 – 111 с.

3. Luebke, D., GPGPU: General Purpose Computation On Graphics Hardware / Luebke D., Harris M., Kruger J., и др. – SIGGRAPH, 2005. – 277 c

4. Tarditi, D.Accelerator: Using Data Parallelism to Program GPUs for General-Purpose Uses / Tarditi D., Puri S., Oglesby J – Microsoft Research, 2006 – 11 с.

Приложение А. Исходный код Вейвлет преобразования на GPGPU

dwt_kernel_float.cu — ядра GPGPU

 

#ifndef _DWT_KERNEL_H_

#define _DWT_KERNEL_H_

 

#include <stdio.h>

#include "dwt_float.h"

 

__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)

{

__shared__ float sh_h1[NC_SZ]; // h1 в общей памяти

__shared__ float sh_g1[NC_SZ]; // h2 в общей памяти

 

const size_t tidx = threadIdx.x;

if(tidx < nc){ // загружаем коэффициенты в общую память их у нас nc

sh_h1[tidx] = h1[tidx]; // выполняя все потоки блока мы заполняем общую память для блока

sh_g1[tidx] = g1[tidx]; //

} 

 

const size_t idbx = (blockDim.x - nc) * blockIdx.x; // индекс данных соответствующих началу блока

 

__shared__ float sh_ai[2*BL_SZ]; // включает nc/2 элементов перед + (blockDim.x - nc/2) нормальных элементов.

 

const int a_idx = idbx + tidx; //индекс текущего элемента данных

 

const size_t ni = 2*a_idx + nmod;

if( ni >= n) { // имитируем заполнение нулями несуществующих предыдущих элементов

       sh_ai[2*tidx] = 0.0f;

       sh_ai[2*tidx + 1] = 0.0f;

} else {     // заполняем нужными значениями

       sh_ai[2*tidx] = a[ni];

       sh_ai[2*tidx + 1] = a[ni + 1];

}

 

__syncthreads(); // синхронизация потоков для нормальности общей регистровой памяти памяти

 

const size_t nh = n>>1;

if(a_idx < nh && tidx < blockDim.x - nc){ //

   float h = 0, g = 0;

 

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

     {

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

//       const float ani = a[stride * jf];

       const float ani = sh_ai[2*tidx + k];

       h += sh_h1[k] * ani;

       g += sh_g1[k] * ani;

     }

         

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

   scratch[a_idx + n/2] = g;

}

}

 

// Сделать шаг обратной фильтрации - сейчас прямая

__global__ void

 cudaDWTStepBack(float *a, const size_t n, const size_t nmod, const size_t stride, const float * const h2, const float* const g2, const size_t nc, float*scratch)

{ // 128 сильно зависит от nc и количество обрабатываемых данных в одном блоке при слиянии

 

__shared__ float sh_h2[NC_SZ]; // nc коэффициентов

__shared__ float sh_g2[NC_SZ]; // nc коэффициентов

 

const size_t tidx = threadIdx.x; // индекс потока

// if(tidx < 2*nc){ // проверить вероятно лишнее

   

if(tidx < nc){ // загружаем коэффициенты в общую память их у нас nc

sh_h2[tidx] = h2[tidx];

sh_g2[tidx] = g2[tidx];

}

const size_t nc_div_2 = nc/2; 

     

const size_t idbx = (blockDim.x - nc_div_2) * blockIdx.x; // индекс данных соответствующих началу блока

 

const int a_idx = idbx + (tidx - nc_div_2); //индекс текущего элемента данных

   

__shared__ float sh_ai[BL_SZ]; // включает nc/2 элементов перед + (blockDim.x - nc/2) нормальных элементов.

__shared__ float sh_ai1[BL_SZ]; // то же самое, но из второй половины

 

const size_t nh = n>>1;

if(a_idx < nh){ // работаем только с нужным диапозоном

if( a_idx < 0) { // имитируем заполнение нулями несуществующих предыдущих элементов

       sh_ai[tidx] = 0.0f;

       sh_ai1[tidx] = 0.0f;

} else {     // заполняем нужными значениями

       sh_ai[tidx] = a[a_idx];

       sh_ai1[tidx] = a[a_idx + nh];

}

}

   

__syncthreads(); // синхронизация потоков для нормальности общей регистровой памяти

 

// Здесь выполняем работу по слиянию

if(a_idx < nh && tidx >= nc_div_2){ // работаем с потоками больше nc/2, и в нужном диапозоне данных

float ai;

float ai1;

   float s = 0.0f;

float s2 = 0.0f;

 

   for (size_t k = 0, k_div_2 = 0; k < nc; k+=2, k_div_2++) { // вычисление четного и нечетного членов при слиянии

       ai = sh_ai[tidx - k_div_2];

       ai1 = sh_ai1[tidx - k_div_2];

 

   s += ai*sh_h2[k] + ai1*sh_g2[k]; // для четных элементов

   s2 += ai*sh_h2[k + 1] + ai1*sh_g2[k + 1]; // для нечетных элементов

}

         

const size_t ni = (n - 1)&(2*a_idx + nmod);

scratch[ni] = s;

scratch[ni + 1] = s2;

}

// }

}

__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;

// tmp[idx] = t;

}

}

 

#endif

 

dwt_float.cu — промежуточный слой между CPU и GPU слоями

 

// Utilities and system includes

#include <shrUtils.h>

#include "cutil_inline.h"

 

#include "dwt_float.h" // Загрузка объявлений функций ядер и их параметров

 

void initGPUDevice()

{

cudaSetDevice(cutGetMaxGflopsDeviceId());

}

 

void exitGPUDevice()

{

cudaThreadExit();

}

 

// перед запуском нужно еще инициализировать адаптер где-то

void cudaDStepForward(float *a, float *d_a, const size_t n, const size_t nmod, const size_t stride, const float * const d_h1, const float* const d_g1, const size_t nc, float*d_scratch)

{

 

// Выделение памяти на видео-карте под данные и коэффициенты

cutilSafeCall(cudaMemcpy(d_a, a, sizeof(*a)*n,

                         cudaMemcpyHostToDevice) );

 

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

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);

 

// Сохранение обработанных данных в память хоста

cutilSafeCall(cudaMemcpy( a, d_scratch, sizeof(*a)*n,

                         cudaMemcpyDeviceToHost) );

 

 

// Выйти из всех потоков

cudaThreadExit();

}

 

// перед запуском нужно еще инициализировать адаптер где-то

void cudaDStepInverse(float *a, float *d_a, const size_t n, const size_t nmod, const size_t stride, const float * const d_h2, const float* const d_g2, const size_t nc, float*d_scratch)

{

// Выделение памяти на видео-карте под данные и коэффициенты

 

cutilSafeCall(cudaMemcpy(d_a, a, sizeof(*a)*n,

                         cudaMemcpyHostToDevice) );

 

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

const size_t block_size = 512;

dim3 threads( block_size, 1);

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

dim3 grid(grid_size, 1);

 

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

cudaDWTStepBack<<< grid, threads >>>(d_a, n, nmod, stride, d_h2, d_g2, nc, d_scratch);

cudaThreadSynchronize(); // синхронизация потоков выполнения на видеокарте

 

// Сохранение обработанных данных в память хоста

cutilSafeCall(cudaMemcpy( a, d_scratch, sizeof(*a)*n, cudaMemcpyDeviceToHost) );

 

// Выйти из всех потоков

cudaThreadExit();

}

 

#define ELEMENT(a,stride,i) ((a)[(stride)*(i)])

void

cudaDwtStepDevice (const gsl_wavelet_float * w, float *a, float *d_a, size_t stride, size_t n,

     gsl_wavelet_direction dir, gsl_wavelet_workspace_float * work, float * d_h, float* d_g, float* d_scratch)

{

size_t nmod;

nmod = w->nc * n ;

nmod -= w->offset;      // center support

 

if (dir == gsl_wavelet_forward)

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

   cudaDStepForward(a, d_a, n, nmod, stride, d_h, d_g, w->nc, d_scratch);

}

else

{  // обратное преобразование

   cudaDStepInverse(a, d_a, n, nmod, stride, d_h, d_g, w->nc, d_scratch);

}

 

}

 

void

cudaDwtTransform (const gsl_wavelet_float * w, float *data, size_t stride, size_t n,

     gsl_wavelet_direction dir, gsl_wavelet_workspace_float * work){

 

float* d_a;

cutilSafeCall(cudaMalloc((void**) &d_a, sizeof(*data)*n));

cutilSafeCall(cudaMemcpy(d_a, data, sizeof(*data)*n,

                         cudaMemcpyHostToDevice) );

   

float* d_scratch;

cutilSafeCall(cudaMalloc((void**) &d_scratch, sizeof(*data)*n));

 

int i;

if (dir == gsl_wavelet_forward)

{

   

float* d_h1;

cutilSafeCall(cudaMalloc((void**) &d_h1, sizeof(*w->h1)*w->nc));

cutilSafeCall(cudaMemcpy(d_h1, w->h1, sizeof(*w->h1)*w->nc,

                         cudaMemcpyHostToDevice) );

float* d_g1;

cutilSafeCall(cudaMalloc((void**) &d_g1, sizeof(*w->g1)*w->nc));

cutilSafeCall(cudaMemcpy(d_g1, w->g1, sizeof(*w->g1)*w->nc,

                         cudaMemcpyHostToDevice) );

 

for (i = n; i >= 2; i >>= 1)

   {

     cudaDwtStepDevice (w, data, d_a, stride, i, dir, work, d_h1, d_g1, d_scratch);

   }

       

cutilSafeCall(cudaFree(d_h1));

cutilSafeCall(cudaFree(d_g1));

 

}

else

{

   

float* d_h2;

cutilSafeCall(cudaMalloc((void**) &d_h2, sizeof(*w->h2)*w->nc));

cutilSafeCall(cudaMemcpy(d_h2, w->h2, sizeof(*w->h2)*w->nc,

                         cudaMemcpyHostToDevice) );

float* d_g2;

cutilSafeCall(cudaMalloc((void**) &d_g2, sizeof(*w->g2)*w->nc));

cutilSafeCall(cudaMemcpy(d_g2, w->g2, sizeof(*w->g2)*w->nc,

                         cudaMemcpyHostToDevice) );

for (i = 2; i <= n; i <<= 1)

   {

     cudaDwtStepDevice (w, data, d_a, stride, i, dir, work, d_h2, d_g2, d_scratch);

   }

       

cutilSafeCall(cudaFree(d_h2));

cutilSafeCall(cudaFree(d_g2));

}

 

cutilSafeCall(cudaMemcpy( work->scratch, d_scratch, sizeof(*d_scratch)*n,

                         cudaMemcpyDeviceToHost) );

cutilSafeCall(cudaMemcpy( data, d_a, sizeof(*d_a)*n,

                         cudaMemcpyDeviceToHost) );

cutilSafeCall(cudaFree(d_scratch));

cutilSafeCall(cudaFree(d_a));

 

}

 

 

void cudaDenoise(float* data, const int len, const float threshold) // Функция самой фильтрации

{

float* d_a;

cutilSafeCall(cudaMalloc((void**) &d_a, sizeof(*data)*len));

cutilSafeCall(cudaMemcpy(d_a, data, sizeof(*data)*len,

                         cudaMemcpyHostToDevice) );

 

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

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);

cudaThreadSynchronize(); // синхронизация потоков выполнения на видеокарте

     

 

cutilSafeCall(cudaMemcpy( data, d_a, sizeof(*data)*len, cudaMemcpyDeviceToHost) );

 

cutilSafeCall(cudaFree(d_a));

}

 

dwt_float.h — заголовочный файл с описание ядер и их параметров

 

extern "C" void initGPUDevice();

extern "C" void exitGPUDevice();

 

__global__ void

 cudaDWTStep(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);

 

__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);

 

// Сделать шаг обратной фильтрации - сейчас прямая

 __global__ void

 cudaDWTStepBack(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);

 

#include "wavelet/gsl_wavelet.h"

 

extern "C" void

cudaDwtTransform (const gsl_wavelet_float * w, float *a, size_t stride, size_t n,

     gsl_wavelet_direction dir, gsl_wavelet_workspace_float * work);

 

__global__ void

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

 

extern "C" void cudaDenoise(float* data, int len, float threshold);

 

#define NC_SZ 128

#define BL_SZ 256

 

wavelet_denoise.cpp – общий ход Вейвлет преобразования

 

int cuda_wavelet_denoise_auto(double* data, int len) // must work on GPGPU

{

float* floatdata = (float*) malloc(len * sizeof(float));

 

/* Объявление переменных */

double* threshold_rescales;     /* Массив масштабов порогов для каждого из уровней разложения */

gsl_wavelet *wavelet;           /* Структура представления вейвлета */

gsl_wavelet_workspace *work;    /* Служебная структура для промежуточных данных */

gsl_wavelet_float waveletfloat;       /* Структура представления вейвлета */

gsl_wavelet_workspace_float workfloat;      /* Служебная структура для промежуточных данных */

int decomp_level;         /* Уровень дискретного вейвлет-разложения */

double threshold;         /* Порог для пороговой обработки коэффициентов */

double thresh_curr;       /* Значение порога, смасштабированное для конкретного уровня разложения */

double len_curr;          /* Длина отрезка массива коэеффициентов для конкретного уровня разложения */

int k;                    /* Счётчик цикла */

int i;

 

/* Инициализация переменных */

 

decomp_level = get_power_of_2(len); /* Вычисление уровня вейвлет-разложения */

if (decomp_level == -1) return -2; /* TODO: код ошибки: "массив данных имеет длину не равную степени 2-ки" */

 

 

threshold_rescales = (double*) malloc(decomp_level * sizeof(double));

if (threshold_rescales == 0) return -1;     //TODO: код ошибки: "Не достаточно памяти"

wavelet = gsl_wavelet_alloc (gsl_wavelet_daubechies, 8);

work = gsl_wavelet_workspace_alloc (len);

 

waveletfloat.nc = wavelet->nc;

waveletfloat.type = wavelet->type;

waveletfloat.offset = wavelet->offset;

waveletfloat.g1 = (float*)malloc(wavelet->nc*sizeof(float));

waveletfloat.g2 = (float*)malloc(wavelet->nc*sizeof(float));

waveletfloat.h1 = (float*)malloc(wavelet->nc*sizeof(float));

waveletfloat.h2 = (float*)malloc(wavelet->nc*sizeof(float));

 

doubleCopyToFloat(wavelet->g1, (float*)waveletfloat.g1, wavelet->nc);

doubleCopyToFloat(wavelet->g2, (float*)waveletfloat.g2, wavelet->nc);

doubleCopyToFloat(wavelet->h1, (float*)waveletfloat.h1, wavelet->nc);

doubleCopyToFloat(wavelet->h2, (float*)waveletfloat.h2, wavelet->nc);

 

workfloat.n = work->n;

workfloat.scratch = (float*)malloc(sizeof(float)*work->n);

doubleCopyToFloat(work->scratch, workfloat.scratch, work->n);

 

initGPUDevice();

 

QTime FullFilterTime;

FullFilterTime.start();

 

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

doubleCopyToFloat(data, floatdata, len);

cudaWaveletTransform_forward (&waveletfloat, floatdata, 1, len, &workfloat);

 

/* ----------------------------------------- */

/* Начало фильтрации */

noise_estimation_float(floatdata, threshold_rescales, decomp_level); /* Вычисление масштабных коэфф-в для каждого уровня разложения */

threshold = select_threshold(len); /* Оценка значения базового порога */

 

i = len / 2;

len_curr = len / 2;

for (k = decomp_level - 1; k >= 0; k--)

{

       thresh_curr = threshold * threshold_rescales[decomp_level - 1 - k];

 

       cudaDenoise( &(floatdata[i]), len_curr, (float)thresh_curr);

 

       len_curr /= 2;

       i -= len_curr;

}

 

/* Конец фильтрации */

/* ----------------------------------------- */

 

/* Обратное вейвлет-преобразование (восстановление сигнала) */

cudaWaveletTransform_inverse (&waveletfloat, floatdata, 1, len, &workfloat);

floatCopyToDouble(floatdata, data, len);

 

exitGPUDevice();

 

qDebug("Cuda FullFilterTime: %d ms", FullFilterTime.elapsed());

 

free((void*)waveletfloat.g1);

free((void*)waveletfloat.g2);

free((void*)waveletfloat.h1);

free((void*)waveletfloat.h2);

free((void*)workfloat.scratch);

 

free (threshold_rescales);

gsl_wavelet_free (wavelet);

gsl_wavelet_workspace_free (work);

 

return 0;

}

 

Приложение Б. Характеристики сопроцессора на основе GPGPU

 

evgeniy@evgeniy-desktop:~/NVIDIA_GPU_Computing_SDK/C/bin/linux/release$ sudo ./deviceQuery

./deviceQuery Starting...

 CUDA Device Query (Runtime API) version (CUDART static linking)

There is 1 device supporting CUDA

Device 0: "GeForce 9800 GT"

CUDA Driver Version:                      3.0

CUDA Runtime Version:                     3.0

CUDA Capability Major revision number:    1

CUDA Capability Minor revision number:    1

Total amount of global memory:            1073020928 bytes

Number of multiprocessors:                14

Number of cores:                          112

Total amount of constant memory:          65536 bytes

Total amount of shared memory per block:  16384 bytes

Total number of registers available per block: 8192

Warp size:                                32

Maximum number of threads per block:      512

Maximum sizes of each dimension of a block: 512 x 512 x 64

Maximum sizes of each dimension of a grid: 65535 x 65535 x 1

Maximum memory pitch:                     2147483647 bytes

Texture alignment:                        256 bytes

Clock rate:                               1.50 GHz

Concurrent copy and execution:            Yes

Run time limit on kernels:                Yes

Integrated:                               No

Support host page-locked memory mapping:  No

Compute mode:                             Default (multiple host threads can use this device simultaneously)

 

deviceQuery, CUDA Driver = CUDART, CUDA Driver Version = 134564823, CUDA Runtime Version = 3.0, NumDevs = 1, Device = GeForce 9800 GT

Приложение В. Запуск теста вейвлет преобразования

 

evgeniy@evgeniy-desktop:~/work/dissertation/seismo_source/seismo_detector$ ./seismodetector

sample count: 65536

CPU FullFilterTime: 12 ms

Cuda FullFilterTime: 67 ms

sample count: 131072

CPU FullFilterTime: 26 ms

Cuda FullFilterTime: 67 ms

sample count: 262144

CPU FullFilterTime: 81 ms

Cuda FullFilterTime: 61 ms

sample count: 524288

CPU FullFilterTime: 159 ms

Cuda FullFilterTime: 73 ms

sample count: 1048576

CPU FullFilterTime: 319 ms

Cuda FullFilterTime: 85 ms

sample count: 2097152

CPU FullFilterTime: 677 ms

Cuda FullFilterTime: 97 ms

sample count: 4194304

CPU FullFilterTime: 1401 ms

Cuda FullFilterTime: 133 ms

sample count: 8388608

CPU FullFilterTime: 2855 ms

Cuda FullFilterTime: 209 ms

sample count: 16777216

CPU FullFilterTime: 5964 ms

Cuda FullFilterTime: 355 ms

sample count: 33554432

CPU FullFilterTime: 11979 ms

Cuda FullFilterTime: 638 ms

sample count: 67108864

CPU FullFilterTime: 31802 ms

Cuda FullFilterTime: 1350 ms