Рис. 2.3: импульсы пульсара J0528+2200 (предпоследний луч внизу), период 3.746 секунды (наблюдается не каждый импульс) и мерой дисперсии DM=50.937. В боксе справа заметно вызванное дисперсией смещение импульса по частотам, а также очевидно влияние явления фарадеевского вращения (поток пульсара вдоль частот меняется). Включен режим автоматической калибровки и определение наличия дисперсионного сдвига корреляцией данных в разных частотных интервалах (бокс в верхней части рисунка). Именно в таком режиме были сгенерированы все 83096 рисунков в работе (Samodurov et al, 2016).

Для пульсаров с дисперсионным запаздыванием DM<15 импульсы зарегистрированы для каждого второго пульсара указанной выборки из ATNF, с дисперсионным запаздыванием DM<50 - для каждого третьего, для пульсаров 50<DM<100 - для каждого 25-го, для DM>100 в импульсах не найдено ни одного пульсара. Очевидно влияние уширения импульсов пульсара рассеянием на эффективность выделения импульсов из базы данных.

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

Пути усовершенствования методики поиска импульсных явлений теперь довольно ясны (переход на тонкие данные с ?t=12.5 мс, увеличение временных отрезков для анализа).

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



ГЛАВА 3. ПРАКТИЧЕСКАЯ РЕАЛИЗАЦИЯ ОСНОВНЫХ АЛГОРИТМОВ ОБРАБОТКИ

§1 Постановка задачи

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

? Исследование возможности применения технологий GPGPU в области астрономических вычислений

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

? Сравнение эффективности с другими способами анализа, в том числе высокопроизводительными

? Приобщение научного сообщества к теме исследования

? Повышение компетентности в области вычислений с применением графических ускорителей

Исходя из целей работы были сформированы следующие задачи:

? Выбор релевантной технологии GPGPU

? Реализация вычисления спектрограмм Фурье, их агрегации и вычисления сопутствующих статистических метрик при помощи выбранной технологии в виде проекта

? Замер длительности вычислений

? Сравнение количественных результатов вычислений на графических ускорителях и полученных при использовании других технологий

? Выступление с промежуточными результатами на XXXV всероссийской конференции «Актуальные проблемы внегалактической астрономии» в г. Пущино с целью приобщить к теме научное сообщество



§2 Выбор технологии

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

? Nvidia CUDA - одна из самых простых и эффективных в использовании. Однако, технология применима только к видеокартам Nvidia, оставляя неохваченными видеокарты компании AMD. Эффективность можно также объяснить отсутствием задачи портабельности и возможностью производить полную оптимизацию под конкретное железо.

? Microsoft DirectCompute - вычислительные шейдеры из стека DirectX. Шейдеры самые сложные в использовании, но в данном случае это нивелируется интегрированностью DirectX в оконную систему Windows. Таким образом, использование данной технологии проще, чем вычислительных шейдеров OpenGL. Однако вместе с данной технологией приходит ограничение операционной системы - она работает только в среде Windows.

? OpenGL compute shaders - вычислительные шейдеры из стека технологий Khronos Group. Самая сложная в использовании технология, поскольку подразумевает поддержку контекста графического конвейера OpenGL и ручное управление всеми ресурсами. Вместе с DirectCompute являются самыми новыми и неустоявшимися технологиями и в основном нацелены не на полноценную обработку огромных массивов данных, а на дополнительные расчёты в графических приложениях, когда расчёты недостаточно серьёзные, чтобы инициализировать и поддерживать контекст OpenCL.

? OpenCL - технология существует наряду с CUDA и также, как и вычислительные шейдеры OpenGL относится к стеку технологий консорциума Khronos Group. Технология сложнее в использовании, чем CUDA, однако полностью кроссплатформенная, что является её главным плюсом. Одна из самых устоявшихся технологий.

Решение рассчитано на использование на компьютерах, схожих по мощности с домашними ПК, а они могут иметь абсолютно различные видеокарты. Следует отметить, что в большинстве бюджетных сборок используются именно видеокарты AMD, а порой и вовсе не используются (в таких случаях используется «встроенная» видеокарта - часть мощности ЦП уходит на рендер. Тем не менее, все технологии, кроме CUDA, работают). Исходя из этого можно отбросить технологию CUDA, ограничивающую решение до использования только на видеокартах Nvidia. Менее критично ограничение по операционной системе, выдвигаемое вычислительными шейдерами DirectX, тем не менее, учитывая неполную состоятельность технологии, данную технологию стоит также отбросить. По причине несостоятельности и отличного от решаемой задачи предназначения можно также отбросить и технологию вычислительных шейдеров OpenGL. Таким образом, выбранная технология GPGPU - OpenCL, являющаяся полностью кроссплатформенной и состоявшейся.

графический ускоритель фурье bitonic

§3 Особенности работы с графическими ускорителями

Прежде чем разрабатывать конкретную архитектуру, необходимо разобраться в особенностях работы графических ускорителей. Видеокарта в расчётах выступает в роли своего рода удалённого сервера - все данные, необходимые для расчётов копируются из оперативной памяти компьютера в оперативную память графического ускорителя. Это создаёт нагрузку на ЦП, которая к тому же не приносит никакой пользы и занимает время, поэтому возникает основной принцип - необходимо минимизировать объём данных, передаваемых и считываемых с видеокарты, а также количество данных операций. В идеале - необходимо один раз передать данные на видеокарту, полностью их обработать и считать конечный результат, что конечно не представляется возможным ввиду ограничения общего объёма памяти на видеокарте и максимального объёма единовременной инициализации буффера на графическом ускорителе - даже на самых мощных видеокартах максимальный размер буфера в глобальной памяти не превышает 2-3 Гбайта. Таким образом, необходимо найти баланс между производительностью и техническими ограничениями.

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

Рабочие группы представляют собой объединение рабочих элементов. Главным плюсом их использования является быстрая общая память, которая как правило являет собой около 32 Кбайт. При этом, стоит помнить, что общее количество рабочих элементов должно быть кратно размеру рабочей группы. Так могут появиться ненужные запуски ядра, поэтому прежде, чем выполнять какую-либо работу, всегда необходимо проверять ID ядра.

В контексте OpenCL ядра пишутся на языке C стандарта C99 с расширением, позволяющим использовать функции, присущие ядрам, а также давать специальные сигнатуры функциям (__kernel) и параметрам (__global, __local, __private). Факт, что это язык C, позволит впоследствии использовать оптимизационные директивы препроцессора при компиляции ядер.

§4 Разработка архитектуры проекта

Первая задача заключается в применении анализа Фурье на сырые данные окном в 200 секунд и сдвигом в 1 минуту. Так получится 57 спектрограмм для канала за 1 час наблюдений (3 уйдут из-за нехватки данных в конце). Далее спектрограммы по всем каналам необходимо агрегировать - усреднить и высчитать мощность. После этого необходимо вывести данные.

Исходя из разобранных выше особенностей работы графических ускорителей и задачи мы имеем следующую архитектуру:

? В оперативной памяти 2 буфера:

_ Буфер с сырыми данными (H_raw)

_ Буфер с результатом (H_res)

? В графической памяти 3 буфера:

_ Буфер с сырыми данными (копия) (D_raw)

_ Буфер с результатом преобразований Фурье, ещё не агрегированный (D_fft)

_ Буфер с результатом (копия) (D_res)

Внутренняя структура буферов представлена на Рисунке 3.1.

Рисунок 3.1. Виртуальная структура буферов: сырых данных (сверху), результатов преобразований Фурье (посередине) и агрегированный буфер с результатами преобразований Фурье (снизу).

Последовательность действий следующая:

1. Считать данные из файла в оперативную память (H_raw)

2. Скопировать данные из H_raw на видеокарту в D_raw

3. Произвести преобразование Фурье, задав целевым буфер D_fft

4. Агрегировать данные из D_fft в буфер D_res

5. Скопировать данные из D_res в оперативную память в буфер H_res

6. Записать результат из H_res в файл

Наглядное отображение данной логики представлено на Рисунке 3.2.

Рисунок 3.2. Наглядная структура работы программы (кроме метрик).

Вторая задача заключается в вычислении статистических метрик для каждого канала. Их необходимо искать для каждого фрагмента в 10-15 секунд для каждого из каналов и часть из них необходимо вычислять без учёта 33% максимальных значений и 1% минимальных значений для ликвидации выбросов. Метрики включают:

1. Максимальное значение

2. Индекс максимального значения

3. Минимальное значение

4. Индекс минимального значения

5. Среднее медианное после после урезания

6. Математическое ожидание для всего окна

7. Математическое ожидание после урезания

8. Дисперсия для всего окна

9. Дисперсия для окна после урезания

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

Для вывода результатов необходим новый буфер (D_met) и его копия в оперативной памяти (H_met), а в качестве исходных данных будет использован созданный в прошлом пункте буфер D_raw.

Что касается промежуточных буферов, то они не представляются необходимыми. Даже на “быстрой” моде количество элементов в “окне” будет примерно равно 800-1600, что даже с учётом удвоения ввиду сохранения их ID при сортировке, займёт около 8 Кбайт локальной памяти, а это не способно превысить лимит локальной памяти. Так представляется возможным расчёт всех метрик с использованием только локальной памяти, что приведёт к сильному увеличению быстродействия, но и потребует завершения всех вычислений за один запуск ядра, чтобы не потерять локальную память.

Таким образом, последовательность действий следующая:

1. Скопировать “окно” из глобальной памяти видеокарты в локальный массив, параллельно запоминая их ID в смежный массив.

2. Отсортировать массив, параллельно меняя порядок элементов и в массиве с ID.

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

4. Рассчитать математическое ожидание до и после отброса части элементов. Это будет необходимо при расчёте дисперсий.

5. Рассчитать дисперсии.

6. Вывести метрики в глобальный буфер D_met, считать его в H_met и вывести в файл.

§5 Реализация алгоритмов Фурье

В качестве среды разработки была выбрана Visual Studio 2015/2017 благодаря её удобству. В качестве API библиотеки OpenCL использовался Intel ocl SDK. Библиотека для проведения преобразований Фурье - clFFT, так как является официально поддерживаемой компанией AMD open-source библиотекой. В качестве тестовых данных был взят 1 час наблюдений в «быстрой» моде из 2014 года весом в 2 Гб. Он представляет собой 280 тыс. точек по 32 каналам и 48 лучам. Точки имеют тип float (4 байта). Для отправки на видеокарту был выбран наиболее удобный формат в 280 тыс. точек по 32 каналам одного из 48 лучей файла. Так данные будут наиболее удобны для агрегирования.

Первое, с чем возникли проблемы - считывание данных из файла. Для данной задачи уже был готовый код в файле open_science_data.h. Однако, считывание файла в 2 Гбайта занимало порядка 30 секунд независимо от скорости работы жёсткого диска, что показалось довольно странным. При более подробном рассмотрении была найдена причина , которая заключалась в том, что считывание последовательно производилось маленьким окном в 4 байта (Рисунок 3.3, пункт 1). В качестве решения считывание было перенесено на два уровня цикла выше, а размер буфера увеличен (Рисунок 3.3, пункт 2). Это уменьшило время считывания файла с диска с 30 до ~7 секунд. Дальнейшие улучшения могут быть произведены при переосмыслении считывания и ином порядке переменных при записи в массив.



Рисунок 3.3. Код для считывания данных из файлов .ptr и .pthr. Цифры: 1 - то, как считывалось до оптимизации, 2 - новая логика.

Так необходима портабельность между производителями видеокарт, а также подразумевается использование старых видеокарт, необходимо задать использование стандарта OpenCL 1.2 который поддерживается видеокартами NVidia. Самые новые стандарты (2.0 и 2.1) ими не поддерживается.Это возможно сделать при помощи директивы препроцессора .

Применение преобразований Фурье не представляется трудным ввиду огромного количества open-source библиотек. В качестве таковой выбрана библиотека clFFT.

Однако, несмотря на кажущуюся простоту, в процессе реализации возникла проблема со стыками. Связано это с тем, что библиотека clFFT реализует лишь одну систему batch-ей, которая позволяет задать отступы и количество необходимых преобразований Фурье только в одних экземплярах, а так как за раз мы хотим отправить как можно больше данных на видеокарту, то они будут храниться в виде одномерного массива (каналы идут подряд) и при указании смещения, как у “окна”, на стыках 32 каналов “окно” анализа произведёт лишние преобразования Фурье. (см. Рисунок 3.4) При решении этой проблемы возникает два логичных решения:

? Допустив копирование данных, перекопировать только нужные первичные данные в новый буфер, пройдясь по ним “окном”. Так количество передаваемых данных увеличится в 3 раза, так как одна и та же точка участвует в преобразования по нескольку раз, а также будут затрачены дополнительные силы на само копирование.

? Ничего не делать на этапе преобразований фурье, допустив произведение 5% большего количества преобразований Фурье. Затем их придётся игнорировать на этапе преобразований.

Рисунок 3.4. Наглядно показанная работа преобразований Фурье и проблема с произведением лишних вычислений на стыках каналов.

Было решено реализовать оба варианта и экспериментально проверить, что предпочтительнее: увеличение нагрузки на канал передачи данных между графическим ускорителем или на его ядра процессора. Результаты можно увидеть в Таблице 3.1.

Таблица 3.1

	Увеличение буфера	Увеличение расчётов
Суммарное время RAM-GPU-RAM файла в 2 Гб	4.3 сек	0.75 сек



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

§6 Реализация алгоритмов агрегирующей статистики данных

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

Следующая задача - поиск метрик. Для неё необходимо реализовать специальное ядро, которое будет производить вычисления за 1 раз. Для удобства был выбран размер “окна” в 1024 элемента, как ближайшая степень двойки. Это сильно упростит задачу при сортировке массива и последующих вычислениях.

Первая и самая сложная задача - сортировка массива. Самый эффективный алгоритм для параллельной сортировки - Bitonic. В непараллельном виде он имеет время выполнения равное , что больше, чем остальных массивов, однако в параллельном виде способен выдавать , что сильно меньше стандартных алгоритмов QuickSort и MergeSort.

Ввиду необходимости производить сортировку локально, внутри ядра и для известного размера “окна”, было принято решение реализовать алгоритм заново в виде функции. Принцип его работы для массивов размером степени двойки можно наглядно увидеть на Рисунке 3.5. Стрелкой показано направление сравнения, то есть в возрастающем или убывающем виде впоследствии будут сравниваемые элементы. Изначально выбирается размер окна и производятся сравнения в симметричных относительно друг друга размерах, то есть первое окно в возрастающем порядке, второе в убывающем, третье опять в возрастающем и т.д. Размер окна постоянно растёт в 2 раза, пока не достигнет длины всей последовательности. После применения “окна” размером N, применяются все прежде использованные “окна”, при этом их направление сравнения будет равно направлению главного текущего “окна”. На рисунке главные окна обозначены буквами A, B, D, G и K. Внутри самого окна сравнение проходят элементы, расстояние между которыми равно половине размера окна.

Рисунок 3.5. Наглядная работа алгоритма сортировки Bitonic для 32 элементов. Данный способ валиден только для массивов с количеством элементов равным степеням двойки (Sergio Loff, 2013)

Реализацию алгоритма Bitonic для 1024 элементов можно увидеть на Рисунке 3.6. При реализации алгоритма на GPU было необходимо учитывать, что за каждую итерацию производится сравнений, что равно 512. Таким образом для наибольшей эффективности количество рабочих элементов в рабочей группе должно быть равно 512. Для каждой итерации вычисляется ID элемента и его пары, после чего производится сравнение и замена элементов в основном и вспомогательном массиве.

Между итерациями рабочим элементам необходима синхронизация, для чего следует использовать функцию (81 строка кода на Рисунке 3.6). Данная функция позволяет остановить дальнейшее выполнение всех локальных рабочих элементов, пока все остальные не доберутся до барьера. Однако, данная функция несёт в себе опасность ошибки: даже те рабочие элементы, которые не должны выполнять никакой работы всё равно должны выполнить данную функцию, иначе программа сломается.

На 67-й строке кода на Рисунке 3.6 можно заметить оптимизационную директиву препроцессора. (www.khronos.org) Её возможно было использовать в связи с тем, что нам известны все размеры в цикле на 68 строке, поэтому данная директива позволяет компилятору “развернуть” внешний цикл в прямолинейную программу, что избавляет от инициализации и инкремента ненужных переменных, за счёт чего ускоряется выполнение кода. Применение данной директивы снизило время прохождения всего ядра на текущей аппаратуре примерно c 0.27 до 0.15 секунд, что явило собой ускорение на ~40%.

Рисунок 3.6. Реализация алгоритма Bitonic для массива в 1024 элемента в OpenCL

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

Их вычисления производятся по формулам: для математического ожидания и для дисперсии. Данные вычисления производятся примерно одинаково, поэтому целесообразно разобрать только самую показательную функцию, которой является расчёт урезанной дисперсии. Любой расчёт производится в 2 этапа. На первом этапе вспомогательный массив заполняется значениями с определёнными коэффициентами без суммы, то есть в случае с дисперсией: . Для переменных по урезанному ряду также происходит обнуление крайних элементов: 0. Между этапами происходит синхронизация локальной памяти, а на втором этапе происходит суммирование, которое для удобства вынесено в отдельную функцию sumReduce.

Функция sumReduce (см. Рисунок 3.7) работает по принципу уменьшения длины ряда в 2 раза при каждой итерации. Рабочий элемент берёт элемент симметричный данному и суммирует их. При этом количество рабочих элементов, реально выполняющих работу, уменьшается вместе с длиной ряда и для каждой итерации их количество равно половине текущей длины ряда. Последняя итерация происходит с использованием всего 1 рабочего элемента. Между каждой итерацией происходит локальная синхронизация (строка 98).

Так как изначально известен размер массива, можно, как и в случае с Bitonic попытаться использовать “развёртку” цикла. Тем не менее, это заметно не повлияло на время работы ядра.



Рисунок 3.7. Реализация алгоритма свёрточного суммирования в OpenCL для массивов длиной, равной степеням двойки,

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

§7 Анализ производительности

Расчёты (без метрик) были на мощном ПК и на ноутбуке с менее мощным графическим ускорителем - Nvidia Geforce GTX 950M. Имеет смысл привести сравнение. (см. Таблица 3.2)

Таблица 3.2

Модель видеокарты	GTX 950M	GTX 1080TI
Количество процессоров CUDA	640	3584
Мощность	1.4 TFlops	12 TFlops
Время расчета тестового файла (1.8 Гб)	4.37 сек	0.75 сек



Из приведенной таблицы можно увидеть прямую зависимость скорости расчётов от количества ядер графического ускорителя. Время различается примерно в 6 раз и примерно в столько раз больше CUDA-процессоров в более мощной видеокарте. В то же прямой зависимости от мощности вычислений с плавающей запятой выявлено не было, так как мощность различается примерно в 10 раз, а время всего в 6. Можно предположить, что отставание по мощности компенсируется не таким сильным отставанием по скорости копирования буферов в памяти, так как пропускного канала между оперативной памятью и графическим ускорителем как раз хватает для мгновенной передачи одним пакетом, как на более, так и на менее мощной видеокарте.

Следует также сравнить время выполнения с анализом, использующим другие способы анализа: (см. Таблица 3.3)

Таблица 3.3

Платформа	GPU (OpenCL)	сервер (24 ядра)	CPU
Стандартные данные (1 час, 46 Мб)	0.3 секунды	9 секунд	...
Быстрые данные (1 час, 1.7 Гб)	0.8 секунды	45 секунд	> 10 минут

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

После введения поиска метрик время расчётов кода было увеличено до 1.0-1.1 секунды с 0.8 секунд, однако при помощи оптимизационных директив препроцессора опустилось до 0.9 секунд.

Время считывания файла после оптимизации снизилось с 30 секунд до 7. Запись файла же зависит от вида записи. В случае бинарной записи, она занимает суммарно порядка 0.1-0.2 секунд, однако в случае записи текстом она может доходить до 1.5-2 минут. Итого: суммарное время выполнения приложения для одного файла с выводом в бинарный файл составляет порядка 8 секунд.

§8 Подведение итогов

По результатам данного исследования был произведен доклад в рамках XXXV Всероссийской конференции “Актуальные проблемы внегалактической астрономии”:

Павлов С.Ю., Тюрин В.А., Самодуров В.А.: Обработка данных круглосуточного обзора БСА ФИАН (110 МГц) методами высокопроизводительных вычислений. // XXXV Всероссийская конференция "АКТУАЛЬНЫЕ ПРОБЛЕМЫ ВНЕГАЛАКТИЧЕСКОЙ АСТРОНОМИИ", тезисы докладов. 24-27 апреля 2018 г., Пущино, с. 43

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

Как итог, в рамках работы были выполнены следующие задачи:

? Проанализированы основные технологии вычислений при помощи графических ускорителей.

? Выбрана наиболее подходящая технология для анализа имеющегося массива астрономических данных в конкретных условиях - OpenCL.

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

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

? С применением разработанной архитектуры было выполнено решение с использованием технологии OpenCL и языков C\С++.

? Во время создания решения были экспериментально проверены некоторые особенности работы с графическими ускорителями.

? Решение было оптимизировано с использованием директив препроцессора.

? Решение было проверено на тестовых данных на двух разных видеокартах Nvidia - GTX 1080TI и GTX 950M, различающихся по мощности в 5-6 раз.

? Был проведен сравнительный анализ времени выполнения двух разных решений проблемы стыков и выбрано наилучшее решение.

? Был проведен сравнительный анализ времени выполнения на двух разных по мощности графических ускорителях.

? Был проведен сравнительный анализ времени выполнения анализа на графических ускорителях и анализа с использованием других способов обработки данных.



ЗАКЛЮЧЕНИЕ

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

В процессе выполнения дипломной работы были созданы решения для типовых задач, которые значительно ускоряют вычисления (на 2-3 порядка) при анализе астрономических данных. Это: фурье анализ поступающих данных и потоковый анализ данных на предмет нахождения среднемедианных значений, дисперсий, стандартных отклонений и т.п. для отрезков данных в несколько секунд.

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

Дополнительно предложенные решения является полностью независимыми от операционной системы и производителя железа, могут быть скомпилированы под любую из основных операционных систем и работать на любом графическом ускорителе, поддерживающем стандарт OpenCL и имеющем в себе около 500 мегабайт памяти.Выступление на XXXV всероссийской радиоастрономической конференции с промежуточными результатами показало, что научное сообщество заинтересовано в дальнейшем развитии темы данной работы и способно найти практическое применение созданному решению.



СПИСОК ЛИТЕРАТУРЫ

· Брауде С. Я., Мень А. В., Содин. Л.Г. Радиотелескоп декаметрового диапазона волн УТР-2 // Антенны. - М.: Связь, 1978, Т. 21, № 2. -- С. 83-131.

· Р.Д.Дагкесаманский, В.А.Самодуров, К.А.Лапаев: «Обзор неба на частоте 102.5 МГц: радиоисточники в интервалах склонений 27.5 - 33.5 и 67.5 - 70.5 градуса.» // А.Ж., 2000, т. 77, N 1, с. 38-46

· А. А. Коноваленко и др.: Астрофизические исследования с помощью малоразмерных низкочастотных радиотелескопов нового поколения // Радиофизика и радиоастрономия. - 2016. - Т. 21, № 2. - С. 83-131.

· Родин А.Е., Орешко В.В., Самодуров В.А.: Поиск источников периодического излучения на радиотелескопе БСА ФИАН // Астрономический Журнал, 2017, 94, т. 1, с. 35-43.

· Рудницкий Г.М.: Конспект лекций по курсу "радиоастрономия". Учебное пособие для студентов. Электронный ресурс, 2001, http://heritage.sai.msu.ru/ucheb/Rudnickij/

· Орешко В.В., 2014: Радиотелескопы ПРАО - состояние и перспективы. http://www.prao.ru/conf/rrc2014/docs/22092014/06_Oreshko.pdf

· Самодуров В.А. и др.: Поиск радио-откликов на гравитационно-волновые события LIGO/VIRGO в данных обзора БСА ФИАН на 110 Мгц. // Radio Physics and Radio Astronomy, 2017, Т. 22. № 4. с. 284-293

· Тюльбашев С.А. и др.: Обнаружение новых пульсаров на частоте 111 МГц // Астрономический журнал, 2016, Том 93, Выпуск №2, стр. 177-190

· Тюльбашев С.А. и др.: Обнаружение пяти новых источников типа RRAT на частоте 111 МГц // Астрономический журнал, 2018, Том 95, Выпуск №1, стр. 68-77

· Astronet.ru: Прозрачность земной атмосферы http://www.astronet.ru/db/msg/1188575

· Bagchi, Manjari: A Unified Model for Repeating and Non-repeating Fast Radio Bursts // The Astrophysical Journal Letters, 2017, Volume 838, Issue 2, article id. L16, 5 pp.

· Chandra, Poonam; Frail, Dale A.: A Radio-selected Sample of Gamma-Ray Burst Afterglows // The Astrophysical Journal, 2012, Volume 746, Issue 2, article id. 156, 28 pp.

· Coenen, T. et al.: The LOFAR pilot surveys for pulsars and fast radio transients. Astronomy & Astrophysics, 2014, vol. 570, id.A60, 16 pp

· Large Synoptic Survey Telescope, LSST - www.lsst.org

· Lingam, Manasvi; Loeb, Abraham: Fast Radio Bursts from Extragalactic Light Sails // The Astrophysical Journal Letters, 2017, Volume 837, Issue 2, article id. L23, 5 pp.

· Katz, J. I.: Fast radio bursts -- A brief review: Some questions, fewer answers // Modern Physics Letters A, 2016, Volume 31, Issue 14, id. 1630013

· Lorimer, D. R. et al; Crawford, F.,: A Bright Millisecond Radio Burst of Extragalactic Origin // Science, 2007, Vol. 318, Is. 5851. -- P. 777-780

· Hewish A., Bell S. J., Pilkington J. D. H., Scott P.F., Collins R. A.: Observation of a Rapidly Pulsating Radio Source // Nature, 1968, Vol. 217, pp. 709

· Khronos group - www.khronos.org

· Manchester R.N., Hobbs G.B., Teoh A. & Hobbs M. : "The ATNF Pulsar Catalogue" // Astronomical Journal, 2005, 129, 1993-2006, http://www.atnf.csiro.au/research/pulsar/psrcat for updated versions.

· Petroff, E. et al: FRBCAT: The Fast Radio Burst Catalogue // Publications of the Astronomical Society of Australia, 2016, Volume 33, id.e045, 7 pp., on-line: http://frbcat.org/

· Pozanenko A.S. et al: GRB 170817A Associated with GW170817: Multi-frequency Observations and Modeling of Prompt Gamma-Ray Emission // Astrophysical Journal Letters, 2018. Vol. 852. No. 2. P. 30-48.

· Samodurov V.A. et al: The daily 110 MHz sky survey (BSA FIAN): online database, science goals data processing by distributed computing. // Труды ХVII международной конференции «Аналитика и управление данными в областях с интенсивным использованием данных» (“Data Analytics and Management in Data Intensive Domains” (DAMDID)). Обнинск: НИЯУ МИФИ, 2015. P. 127-128.

· Samodurov V. A. et al: Observations of Transient Phenomena in BSA Radio Survey at 110 MHz // XVIII International Conference, DAMDID/RCDL 2016, Ershovo, Moscow, Russia, October 11 -14, 2016, Revised Selected Papers, Springer International Publishing, 2017, P. 130-141, http://doi.org/10.1007/978-3-319-57135-5

· Scholz, P. et al: The Repeating Fast Radio Burst FRB 121102: Multi-wavelength Observations and Additional Bursts // The Astrophysical Journal, 2016, Volume 833, Issue 2, article id. 177, 17 pp.

· SKA TELESCOPE. SQUARE KILOMETRE ARRAY https://www.skatelescope.org/

· Sйrgio Loff: Sorting algorithms on a GPU [Электронный ресурс] // Performance guidelines, 2013 - URL: http://performanceguidelines.blogspot.ru/2013/08/sorting-algorithms-on-gpu.html (Дата обращения: 2.5.2018).

· Taylor, G.B. et al: First Light for the First Station of the Long Wavelength Array. J. Astron. Instrum., vol. 1, p. 1-56 (2012)

· Tyul'bashev, S. A.; Chernikov P. A.: Properties of CSS radio sources from 102 MHz interplanetary scintillation observations // Astronomy and Astrophysics, v.373, p.381-393 (2001)

· Van Haarlem, M. P. et al.: LOFAR: The LOw-Frequency ARray // Astronomy & Astrophysics, 2013, Vol. 556, id.A2, 53 pp.

· Wang, F. Y.; Yu, H.: SGR-like behaviour of the repeating FRB 121102

· // Journal of Cosmology and Astroparticle Physics, 2017, Issue 03, article id. 023

· Whiting M., Humphreys B.: Source-Finding for the Australian Square Kilometre Array Pathfinder // Publications of the Astronomical Society of Australia, 2012, v. 29, Issue 3, pp. 371-381



Приложение 1. Исходный код ядер (kernels.cl)

__kernel

void aggregateFFTs(

__global float *fftBuff,

__global float *target,

const int resultBandwidth,

const int totalBandwidth2,

const int totalWidth2,

const int bandCount)

{

int id = get_global_id(0);

if (id < resultBandwidth)

{

float reSum = 0.0f;

float imSum = 0.0f;

float powSum = 0.0f;

float re = 0.0f, im = 0.0f;

for (size_t i = 0; i < totalWidth2; i += totalBandwidth2)

{

re = fftBuff[id * 2 + i] + fftBuff[id * 2 + i];

im = fftBuff[id * 2 + 1 + i] + fftBuff[id * 2 + 1 + i];

powSum += re * re + im * im;

reSum += re;

imSum += im;

}

target[id * 3] = reSum / bandCount;

target[id * 3 + 1] = imSum / bandCount;

target[id * 3 + 2] = powSum / bandCount;

}

}

inline void compare(

__local float *A,

__local float *B,

__local float *idA,

__local float *idB,

bool dir)

{

if ((*A > *B) == dir) {

//#pragma HLS PIPELINE

float t;

t = *A; *A = *B; *B = t;

t = *idA; *idA = *idB; *idB = t;

}

}

inline void BitonicSort1024(

__local float * arr,

__local float * arrHelper,

const unsigned int pairCount,

__private int localId)

{

#pragma unroll

for (int globalGran = 2; globalGran <= 1024; globalGran <<= 1)

{

for (int localGran = globalGran; localGran >= 2; localGran >>= 1)

{

if (localId < pairCount)

{

int halfLocalGran = localGran >> 1;

int arrId = localId % halfLocalGran + (localId / halfLocalGran) * localGran;

int arrPairId = arrId + halfLocalGran;

bool dir = ((arrId / globalGran) % 2) == 0; // true = asc, false = desc

compare(&arr[arrId], &arr[arrPairId], &arrHelper[arrId], &arrHelper[arrPairId], dir);

}

barrier(CLK_LOCAL_MEM_FENCE);

}

}

}

inline void sumReduce(

__local float * arr,

const unsigned int pairCount,

__private int localId)

{

for (int i = pairCount; i > 0; i >>= 1)

{

if (localId < i)

arr[localId] += arr[i + localId];

barrier(CLK_LOCAL_MEM_FENCE);

}

}

inline void mathExpectationCompute(

__local float * arr,

__local float * arrHelper,

__local float * mathExp,

const unsigned int pairCount,

__private int localId)

{

if (localId < pairCount)

{

arrHelper[localId + localId] = arr[localId + localId] / (pairCount << 1);

arrHelper[localId + localId + 1] = arr[localId + localId + 1] / (pairCount << 1);

}

barrier(CLK_LOCAL_MEM_FENCE);

sumReduce(arrHelper, pairCount, localId);

if (localId == 0)

*mathExp = arrHelper[0];

}

inline void mathExpectationShortenedCompute(

__local float * arr,

__local float * arrHelper,

__local float * mathExpShortened,

const unsigned int pairCount,

__private int localId)

{

if (localId < pairCount)

if (localId >= 5 && localId < 343) {

arrHelper[localId + localId] = arr[localId + localId] / (676);

arrHelper[localId + localId + 1] = arr[localId + localId + 1] / (676);

}

else {

arrHelper[localId + localId] = 0.0f;

arrHelper[localId + localId + 1] = 0.0f;

}

barrier(CLK_LOCAL_MEM_FENCE);

sumReduce(arrHelper, pairCount, localId);

if (localId == 0)

*mathExpShortened = arrHelper[0];

}

inline void varianceCompute(

__local float * arr,

__local float * arrHelper,

__local float * variance,

__local float * mathExp,

const unsigned int pairCount,

__private int localId)

{

if (localId < pairCount)

{

float val;

val = arr[localId + localId] - *mathExp; arrHelper[localId + localId] = val * val / (pairCount << 1);

val = arr[localId + localId + 1] - *mathExp; arrHelper[localId + localId + 1] = val * val / (pairCount << 1);

}

barrier(CLK_LOCAL_MEM_FENCE);

sumReduce(arrHelper, pairCount, localId);

if (localId == 0)

*variance = arrHelper[0];

}

inline void varianceShortenedCompute(

__local float * arr,

__local float * arrHelper,

__local float * varianceShortened,

__local float * mathExpShortened,

const unsigned int pairCount,

__private int localId)

{

if (localId < pairCount)

{

if (localId >= 5 && localId < 343)

{

float val;

val = arr[localId + localId] - *mathExpShortened; arrHelper[localId + localId] = val * val / (676);

val = arr[localId + localId + 1] - *mathExpShortened; arrHelper[localId + localId + 1] = val * val / (676);

}

else

{

arrHelper[localId + localId] = 0.0f;

arrHelper[localId + localId + 1] = 0.0f;

}

}

barrier(CLK_LOCAL_MEM_FENCE);

sumReduce(arrHelper, pairCount, localId);

if (localId == 0)

*varianceShortened = arrHelper[0];

}

__kernel

void getMetrics(

__global float *buff,

__global float *metrics,

__local float *arr,

__local float *arrHelper,

__local float *metLocal,

const unsigned int pairCount,

const unsigned int smallBatchCount,

const unsigned int smallBatchShift,

const unsigned int bigBatchCount,

const unsigned int bigBatchShift)

{

int bigBatchId, smallBatchId;

int localId = get_local_id(0); // size: 512

int wgId = get_group_id(0); // size: total batch count

if (localId < pairCount)

{

bigBatchId = wgId / smallBatchCount;

smallBatchId = wgId % smallBatchCount;

int buffId = localId + localId + smallBatchId * smallBatchShift + bigBatchId * bigBatchShift;

arr[localId + localId] = buff[buffId];

arr[localId + localId + 1] = buff[buffId + 1];

arrHelper[localId + localId] = localId + localId;

arrHelper[localId + localId + 1] = localId + localId + 1;

}

barrier(CLK_LOCAL_MEM_FENCE); // wait for the array to be copied to the local memory

BitonicSort1024(arr, arrHelper, pairCount, localId);

if (localId == 0) metLocal[0] = arr[1023]; // max val

else if (localId == 1) metLocal[1] = arrHelper[1023]; // max id

else if (localId == 2) metLocal[2] = arr[0]; //min val

else if (localId == 3) metLocal[3] = arrHelper[0]; // min id

else if (localId == 4) metLocal[4] = ((arr[347] + arr[348]) / 2); // median of shortened array: median = (347 + 348) /2

barrier(CLK_LOCAL_MEM_FENCE); // we don't need ids info from arrHelper now.

mathExpectationCompute(arr, arrHelper, &metLocal[5], pairCount, localId);

barrier(CLK_LOCAL_MEM_FENCE);

mathExpectationShortenedCompute(arr, arrHelper, &metLocal[6], pairCount, localId);

barrier(CLK_LOCAL_MEM_FENCE);

varianceCompute(arr, arrHelper, &metLocal[7], &metLocal[5], pairCount, localId);

barrier(CLK_LOCAL_MEM_FENCE);

varianceShortenedCompute(arr, arrHelper, &metLocal[8], &metLocal[6], pairCount, localId);

if (localId >= 0 && localId < 7)

{

int metricsId = (bigBatchId * smallBatchCount + smallBatchId) * 9;

metrics[metricsId + localId] = metLocal[localId];

}

else if (localId >= 7 && localId < 9)

{

int metricsId = (bigBatchId * smallBatchCount + smallBatchId) * 9;

metrics[metricsId + localId] = sqrt(metLocal[localId]);

}

}



Приложение 2. Чтение данных из файла (open_science_data.h)

#pragma once

#include <string>

#include <iostream>

#include <fstream>

#include <cstdlib>

#include <map>

using namespace std;

#define entryType pair<string, string>

namespace DataReader

{

float data[RAY_COUNT][BANDS_COUNT+1][TIME_COUNT];

entryType getEntry(string line) { //divides string into key and value and combines them into map entry

int pos = line.find_first_of(' ');

string key(line, 0, pos);

string value(line, pos);

value = string(value, value.find_first_not_of(' '));

return entryType(key, value);

};

int open_science_data(string filename, int type_of_data = -1)

{

//opening input file for both text and binary reading

ifstream inf(filename.c_str(), ifstream::in | ifstream::binary);

//1. reading header;

map<string, string> params; //map for parameters from file header

string line;

entryType entry;

cout << inf.tellg() << " > "; //same as in further loop

getline(inf, line);

entry = getEntry(line);

params.insert(entry);

cout << entry.first << ": " << entry.second << endl << endl;

int numpar = atoi(entry.second.c_str()); //number of parameters

for (int i = 0; i < numpar - 1; i++)

{

// cout«inf.tellg()«" > ";

getline(inf, line); //reading a line with parameter

entry = getEntry(line); //dividing it into key and value

params.insert(entry); //inserting it into map

//cout«entry.first«": "«entry.second«endl«endl;

cout << entry.first << ": " << entry.second << endl;

};

int npointsB = atoi(params.at("npoints").c_str()); //number of points in big file

int nbandsB = atoi(params.at("nbands").c_str()); //number of bands ...

cout << endl << "task 1 complete: header is read;" << endl << endl;

//2. reading data;

float value[33*RAY_COUNT];

for (int t = 0; t < TIME_COUNT; t++)

{

inf.read((char*)value, 4 * (nbandsB + 1)*RAY_COUNT);

for (int r = 0; r < RAY_COUNT; r++)

{

for (int b = 0; b < nbandsB + 1; b++)

{

//inf.read((char*)value, 4);

if (type_of_data == -1) // no checked band : all bands

data[r][b][t] = value[b + (nbandsB + 1)*r];

else if ((b + 1) == type_of_data) //checked band

data[r][0][t] = value[b + (nbandsB + 1)*r];

}

}

}

inf.close();

cout << "test data was read - data[0][1]: " << data[0][1] << " and data[1][1]: " << data[1][1] << endl;

return(0);

};

}



Приложение 3. Управление контекстом OpenCL (ocl_Astronomy.h)

#pragma once

#include <stdio.h>

#include <math.h>

#include <clFFT.h>

#include <iostream>

#define POINT_TYPE float

#define METRICS_LEN 1024

#define METRICS_COUNT 9

namespace cl_Astronomy

{

char * kernelSource;

char * getKernelSource()

{

FILE *fp;

char *source_str;

size_t source_size, program_size;

fp = fopen("kernels.cl", "rb");

if (!fp) {

printf("Failed to load kernel\n");

return "\0";

}

fseek(fp, 0, SEEK_END);

program_size = ftell(fp);

rewind(fp);

source_str = (char*)malloc(program_size + 1);

source_str[program_size] = '\0';

fread(source_str, sizeof(char), program_size, fp);

fclose(fp);

return source_str;

}

cl_int err;

cl_platform_id platform = 0;

cl_device_id device = 0;

cl_context_properties props[3] = { CL_CONTEXT_PLATFORM, 0, 0 };

cl_context ctx = 0;

cl_command_queue queue = 0;

cl_mem D_raw, D_fft, D_sum, D_met;

cl_event event = NULL;

clfftPlanHandle planHandle;

cl_program program;

cl_kernel aggregateKernel, metricsKernel;

size_t globalSize, localSize; // Addition kernel workgroup size

size_t metricsGlobalSize, metricsLocalSize;

int windowSize = 8; // analysis window size. FFT is applied to the sequence of these numbers

int bandCount = 1;

int timeCount = 36000;

int shiftSize = 1;

int fftsPerBand = 1; // fft productive calls per band analysis

int fftTotalCount = 0; // fft total (including non-poductive) calls per ray

int halfedWindowSize = 1; // 1 window fft analysis result productive width (without duplicate ending)

int resultBandwidth = 1; // productive width of 1 band fft analysis result

int totalBandwidth2 = 1; // productive + non-productive width of 1 band fft analysis result

int totalWidth2 = 1; // total width of 1 ray fft analysis results

int metSmallBatchCount = 1; // small batch executions count

int metSmallBatchSize = 1; // small batch points size

int metSmallBatchShift = 1; // shift between small batch executions

int metBigBatchCount = 1; // bigger batch count. Defines how many times should small batches be executed. (IS NOT TOTAL BATCH COUNT!!!)

int metBigBatchSize = 1; // bigger batch size. Defines the point size of bigger batch. Is total number of points in big batch

int metBigBatchShift = 1; // shift between bigger batch executions.

int metPairCount = 1; // window pair count. This should be 512 (window size is 1024)! Don't change this unless you've changed the kernel logic.

// Buffer length

int rawBufLen = 0,

fftBufLen = 0,

sumBufLen = 0,

metBufLen = 0;

float *H_raw, *H_sum, *H_met; // Host device buffer pointers

// call this #1

void setup_ocl()

{

char platform_name[128];

char device_name[128];

/* Setup OpenCL environment. */

err = clGetPlatformIDs(1, &platform, NULL);

size_t ret_param_size = 0;

err = clGetPlatformInfo(platform, CL_PLATFORM_NAME, sizeof(platform_name), platform_name, &ret_param_size);

printf("Platform found: %s\n", platform_name);

err = clGetDeviceIDs(platform, CL_DEVICE_TYPE_DEFAULT, 1, &device, NULL);

err = clGetDeviceInfo(device, CL_DEVICE_NAME, sizeof(device_name), device_name, &ret_param_size);

printf("Device found on the above platform: %s\n", device_name);

props[1] = (cl_context_properties)platform;

ctx = clCreateContext(props, 1, &device, NULL, NULL, &err);

queue = clCreateCommandQueue(ctx, device, 0, &err);

}

void finalize_ocl()

{

/* Release OpenCL memory objects. */

clReleaseMemObject(D_raw);

clReleaseMemObject(D_fft);

clReleaseMemObject(D_sum);

clReleaseMemObject(D_met);

delete[] H_sum;

delete[] H_met;

/* Release OpenCL working objects. */

clReleaseCommandQueue(queue);

err = clReleaseContext(ctx);

}

// call this #2

void initBuffers(const int &wSize, const int &tCount, const int &bCount, const int &sSize)

{

windowSize = wSize;

bandCount = bCount;

shiftSize = sSize;

timeCount = tCount;

fftsPerBand = (timeCount - windowSize) / shiftSize + 1; // how many shifts + 1 initial call

halfedWindowSize = ceil((windowSize - 1) / 2.0f) + 1; // first element + uprounded half of what's left

resultBandwidth = halfedWindowSize * fftsPerBand; // productive width of 1 band fft analysis result

totalBandwidth2 = (resultBandwidth + OMIT_SHIFTS * halfedWindowSize) * 2; // productive + non-productive width of 1 band fft analysis result multiplied by 2 for real and imaginary

totalWidth2 = totalBandwidth2 * bandCount; // total width of 1 ray fft analysis results for real and imaginary

rawBufLen = timeCount * bandCount;

sumBufLen = fftsPerBand * halfedWindowSize * 3;

fftTotalCount = (fftsPerBand + OMIT_SHIFTS) * bandCount - OMIT_SHIFTS; // productive ffts + omitted. Don't do the last ones

fftBufLen = fftTotalCount * halfedWindowSize * 2;

metSmallBatchCount = timeCount/METRICS_LEN;

metSmallBatchSize = METRICS_LEN;

metSmallBatchShift = METRICS_LEN;

metBigBatchCount = bandCount;

metBigBatchSize = timeCount;

metBigBatchShift = timeCount;

metPairCount = (METRICS_LEN >> 1);

metBufLen = metSmallBatchCount * metBigBatchCount * METRICS_COUNT;

/* Prepare OpenCL memory objects */

D_raw = clCreateBuffer(ctx, CL_MEM_READ_WRITE, rawBufLen * sizeof(POINT_TYPE), NULL, &err);

D_fft = clCreateBuffer(ctx, CL_MEM_READ_WRITE, fftBufLen * sizeof(POINT_TYPE), NULL, &err);

D_sum = clCreateBuffer(ctx, CL_MEM_READ_WRITE, sumBufLen * sizeof(POINT_TYPE), NULL, &err);

D_met = clCreateBuffer(ctx, CL_MEM_READ_WRITE, metBufLen * sizeof(POINT_TYPE), NULL, &err);

...