Как развитие данной работы планируется разработка гетерогенных вычислительных методов, которые будут использовать и универсальные процессоры и специализированные параллельно массивные вычислители 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; /* Уровень дискретного вейвлет-разложения */