Размещено на http://www.allbest.ru/

• Содержание
• Введение
• 1. Основные понятия
• 1.1 Сплайны
• 1.2 B-сплайны
• 1.3 Кубические сплайны
• 1.4 Сплайн-поверхность Смоляка
• 2. Инструментарий
• 2.1 Технология CUDA
• 2.1.1 Модель программирования CUDA
• 2.1.2 Модель памяти CUDA
• 2.1.3 CUDA и язык C
• 2.1.4 Библиотека cuBLAS
• 2.1.5 Библиотека cuSPARSE
• 2.1.6 Библиотека cuSOLVER
• 2.2 Инструмент gSOAP
• 2.2.1 Пример использования
• 2.3 MATLAB
• 2.3.1 Синтаксис языка
• 3. Реализация
• 3.1 Веб-служба
• 3.1.1 Решение СЛАУ на графических ускорителях
• 3.1.2 Создание веб-службы с помощью gSOAP
• 3.2 Клиент на MATLAB
• 4. Тестирование
• 4.1 Эллиптический параболоид
• 4.2 Трансцендентная поверхность 1
• 4.3 Трансцендентная поверхность 2
• Заключение
• Список литературы
• ПРИЛОЖЕНИЕ А
• ПРИЛОЖЕНИЕ Б
• ПРИЛОЖЕНИЕ В


Введение

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

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

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

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

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

Поскольку работа на горнодобывающих предприятиях никогда не останавливается и ведется круглосуточно весь срок службы карьера, все расчеты должны производится в режиме реального времени. Необходимость вычислений в режиме реального времени накладывает определенные требования на мощность вычислительного узла. В данной работе предлагается использовать графические ускорители и распараллеливание программы для проведения расчетов в реальном времени. Мы будем использовать графические ускорители от компании Nvidia, лидера в этой области, и их технологию CUDA, предназначенную для проведения массированных расчетов на их графических ускорителях (см. раздел 2.1 "Технология CUDA").

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

Рассматриваемые в следующей главе варианты аппроксимирующих поверхностей, кроме сплайн-поверхности Смоляка, строятся лишь по точкам, размещенным на регулярной сетке с фиксированным шагом. Таким образом, данные поверхности не подходят для моделирования изменения поверхности местности проводимых горных работ, строящегося по нерегулярной сетке.

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

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



1. Основные понятия

1.1 Сплайны

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

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

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

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

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

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

(1.1.1)

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

(1.1.2)

Система из 4-х уравнений



(1.1.3)

позволяет однозначно определить четыре коэффициента полинома. Для полинома 5-й степени следует добавить условие равенства 2-й производной на концах отрезка и т. д. Из сказанного должно быть ясно, почему сплайны строят в основном из полиномов нечётных степеней (с чётным количеством коэффициентов).

Для полиномов чётных степеней при сборке системы (1.1.3):

· остаётся неопределённой производная в одном из концов отрезка;

· и условие равенства производных (гладкости кривой) не будет выполняться,

поэтому для полинома 2-й степени невозможно достичь равенства 1-й производной в точках стыка, а для 4-й степени -- 2-й производной и т. д. Для построения сплайнов с чётными степенями искусственно добавляют дополнительные условия, чтобы сформировать систему уравнений, подобную (1.1.3). Если производные полинома сплайна определяются также, как соответствующие производные интерполируемой функции, сплайн называется является эрмитовым.

(1.1.4)

Существуют локальные методы построения сплайнов Бесселя и Акими, B -- сплайны. В основном, когда речь идет о сплайнах, то имеют в виду сплайны, построенные из алгебраических полиномов. Именно к ним относится приведённое выше определение. Именно эти сплайны являются наиболее изученными. Однако, сплайн может состоять из фрагментов функций любого класса. В рассмотрено построение таких сплайнов и исследуются их свойства. Автор не даёт общего определения построенных сплайнов. Очевидно, что для любых классов функций, из которых состоит сплайн, приведённое в начале статьи определение не совсем подходит. Например, если сплайн состоит из отрезков экспоненты, то понятие дефекта сплайна теряет смысл. Хотя количество непрерывных производных останется важной характеристикой.

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

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

1.2 B-сплайны

Сплайны с локальным носителем. (B - сплайны). В последнее время в вычислительной практике широкое распространение получили B - сплайны (от английского слова bell -- колокол), сосредоточенные на конечном носителе. Они используются как для интерполяции функций, так и в качестве базисных функций при построении методов типа конечных элементов.

Определение. B - сплайном, или базисным сплайном степени N - 1 дефекта 1 относительно узлов называется функция

(1.2.1)

Пусть т.е. рассматривается случай равномерной сетки.

Рассмотрим несколько частных случаев В - сплайнов.

Первый случай. N = 2. В этом случае сплайн строится наиболее просто.

(1.2.2)

или

(1.2.3)



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

Рисунок 1.1 - График функции "крышка".

Рассмотрим случай B - сплайна 2 - го порядка, задаваемого формулой

При t < t_{k - 2}, t > t_{k - 2}, . Построенный сплайн обладает следующими свойствами:

· S't (t_{k - 2}) = S't(t_{k + 2}) = 0 ;

· S(t_{k - 1}) = S(t_{k + 1}) = 1 ;

· S(t_{k - 2}) = S(t_{k + 2}) = 0.

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

b_{i - 1}S_{i - 1} + b_iS_i + b_{i + 1}S_{i + 1}= f_i, (1.2.4)

где b -- коэффициенты интерполяции, S -- B - сплайн, индекс указывает на точку носителя, в которой сплайн достигает своего максимума. Система таких соотношений, естественно, дополняется граничными условиями. Известно [1], что получившаяся система для определения коэффициентов разложения будет иметь трехдиагональную матрицу с диагональным преобладанием при выполнении ограничения на длины соседних шагов: они должны различаться не более чем в раза.

Второй случай. N = 4 (кубический B - сплайн ) имеет вид

(1.2.5)

или, после несложных упрощений:

(1.2.6)

Рисунок 2 - Рассматриваемый сплайн.

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

Для последнего класса задач также используются функции Бернштейна:

(1.2.7)

Функции Бернштейна иногда записывают в форме рекуррентного соотношения:

(1.2.8)

Такие рекуррентные последовательности применяются с целью уменьшения ошибок округления.

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

Подробнее о функциях Бернштейна в [2].

1.3 Кубические сплайны

Некоторая функция f(x) задана на отрезке , разбитом на части



, (1.3.1).

Кубическим сплайном дефекта 1 называется функция , которая:

· на каждом отрезке является многочленом степени не выше третьей;

· имеет непрерывные первую и вторую производные на всём отрезке ;

· в точках выполняется равенство , т. е. сплайн интерполирует функцию f в точках .

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

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

(1.3.2).

Теорема: Для любой функции и любого разбиения отрезка существует ровно один естественный сплайн S(x), удовлетворяющий перечисленным выше условиям.

Эта теорема является следствием более общей теоремы Шёнберга-Уитни об условиях существования интерполяционного сплайна.

Обычно для сплайна выбирают кубический полином

,

определенный на интервале из .

При этом вся кривая представляет собой набор таких кубических полиномов (рис 1.4), с определенным образом подобранными коэффициентами - параметр сплайна [3, 4].

Рисунок 3 - Схема методов сплайнов.

Коэффициенты на каждом интервале определяются из условий сопряжения в узлах:

. (1.3.2)

Кроме того, на границе при и ставятся условия

. (1.3.3)

Будем искать кубический полином в виде

.

Из условия имеем



(1.3.4)

Вычислим производные:

.

и потребуем их непрерывности при :

. (1.3.5)

Общее число неизвестных коэффициентов, очевидно, равно , число уравнений (1.3.4) и (1.3.5) равно . Недостающие два уравнения получаем из условия (1.3.3) при и :

.

Выражение из (1.3.5) , подставляя это выражение в (1.3.4) и исключая , получим:

.



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

. (1.3.6)

С краевыми условиями:

(1.3.7)

Условие эквивалентно условию и уравнению .

Разностное уравнение (1.3.4) с условиями (1.3.5) можно решить методом прогонки, представив в виде системы линейных алгебраических уравнений вида , где вектор соответствует вектору , вектор F поэлементно равен правой части уравнения (1.3.6), а матрица A имеет следующий вид:

,

где и

Метод прогонки, основан на предположении, что искомые неизвестные связаны рекуррентным соотношением:

.

Используя это соотношение, выразим и через и подставим в i-e уравнение:

,

где - правая часть i-го уравнения. Это соотношение будет выполняться независимо от решения, если потребовать:

;

.

Отсюда следует:

;

.

Из первого уравнения получим:

;

.

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

.

Сплайновая интерполяция хороша тем, что требует знания в узлах только значений функции, но не ее производных [5, 6, 7].

1.4 Сплайн-поверхность Смоляка

Сплайн-поверхности были предложены в 1970 г. отечественным математиком С.А.Смоляком [8, 9] а затем, в 1972 г., американскими инженерами Хардером и Демарэ [10] как обобщение модели упругой пластинки бесконечной протяженности, которая деформируется лишь изгибом, причем ее отклонения от исходного состояния задаются в конечном числе независимых точек, где приложены точечные нагрузки. Равновесие бесконечной идеально тонкой пластинки определяется из условия минимума ее свободной энергии.

Мы не будем здесь подробно описывать теорию построения сплайн-поверхностей, которую можно найти в прекрасном изложении В.О.Ашкеназы [11], но приведем взятые из этой публикации минимально необходимые ключевые формулы, чтобы можно было ориентироваться в дальнейшем изложении.

Уравнение сглаживающего сплайна является решением вариационной задачи:



(1.4.1)

где - функционал гладкости сплайн-поверхности , а

(1.4.2)

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

В случае двумерного признакового пространства Z={z₁,z₂}, решение вариационной задачи (1.4.1) ищется в виде:

(1.4.3)

где r_i - это евклидово расстояние от точки i до точки z₁, z₂:

, (1.4.4)

а на коэффициенты c_i накладываются ограничения:

(1.4.5)



Рассмотрев уравнения (1.4.3, 1.4.4, 1.4.5), можно заметить, что на больших расстояниях от всех исходных точек Z_i уравнение D^m-сплайна вырождается в полином:

(1.4.6)

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

(1.4.7)

Решение этой системы существует, если N > 3. Решение является единственным, если среди точек имеются хотя бы три точки, не лежащие на одной прямой.

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

(1.4.8)



Эта система решается любым методом исключения Гаусса с выбором главного элемента (что необходимо, поскольку главная диагональ матрицы системы имеет нули). Для решения данной системы потребуется примерно 2(N + 3)³ арифметических операций и (N + 3)² / 2 ячеек памяти ЭВМ, что в несколько ограничивает максимальное число точек интерполяции N при данном подходе.

Данная сплайн-поверхность полностью подходит для наших целей: строится по нерегулярной сетке точек и является интерполяционной. Нам нужно реализовать алгоритм решения СЛАУ (1.4.8) для нахождения коэффициентов сплайна на графическом ускорителе с использованием технологии CUDA. После решения коэффициенты сплайна будут передаваться клиенту для дальнейших работ со сплайном описанных во введении данной работы.



2. Инструментарий

2.1 Технология CUDA

CUDA (англ. Compute Unified Device Architecture) -- программно-аппаратная архитектура параллельных вычислений, которая позволяет существенно увеличить вычислительную производительность благодаря использованию графических процессоров фирмы Nvidia.

Вычисления с GPU-ускорением заключаются в использовании графического процессора (GPU) совместно с CPU для ускорения приложений в области науки, аналитики, проектирования и на предприятиях. Вычисления с GPU-ускорением были изобретены компанией NVIDIA в 2007 году, теперь GPU обеспечивают мощностью энергоэффективные центры обработки данных в правительственных лабораториях, университетах и на предприятиях всех размеров по всему миру. GPU ускоряют приложения на различных платформах, начиная от автомобилей, мобильных телефонов и планшетов и заканчивая беспилотными летательными аппаратами и роботами.

Вычисления с GPU-ускорением предлагают беспрецедентную производительность приложений благодаря тому, что GPU обрабатывает части приложения, требующие большой вычислительной мощности, при этом остальная часть приложения выполняется на CPU. С точки зрения пользователя, приложение просто работает значительно быстрее.



Рисунок 2.1.1 - Схема взаимодействия GPU и CPU.

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

Рисунок 2.1.2 - GPU состоит из тысяч ядер для эффективной обработки параллельных задач.



Платформа параллельных вычислений CUDA обеспечивает набор расширений для языков C и С++, позволяющих выражать как параллелизм данных, так и параллелизм задач на уровне мелких и крупных структурных единиц. Программист может выбрать средства разработки: языки высокого уровня, такие как C, C++, Fortran или же открытые стандарты, такие как директивы OpenACC. Платформа параллельных вычислений CUDA используется на сегодняшний день в тысячах GPU-ускоренных приложений и тысячах опубликованных научных статьях.

Рисунок 2.1.3 - Сравнение однопоточного и параллельного кода на C.

Направление вычислений эволюционирует от "централизованной обработки данных" на центральном процессоре до "совместной обработки" на CPU и GPU. Для реализации новой вычислительной парадигмы компания NVIDIA изобрела архитектуру параллельных вычислений CUDA, на данный момент представленную в графических процессорах GeForce, ION, Quadro и Tesla и обеспечивающую необходимую базу разработчикам ПО.

Говоря о потребительском рынке, стоит отметить, что почти все основные приложения для работы с видео уже оборудованы, либо будут оснащены поддержкой CUDA-ускорения, включая продукты от Elemental Technologies, MotionDSP и LoiLo.

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

На финансовом рынке компании Numerix и CompatibL анонсировали поддержку CUDA в новом приложении анализа риска контрагентов и достигли ускорения работы в 18 раз. Numerix используется почти 400 финансовыми институтами.

Показателем роста применения CUDA является также рост использования графических процессоров Tesla в GPU вычислениях. На данный момент более 700 GPU кластеров установлены по всему миру в компаниях из списка Fortune 500, таких как Schlumberger и Chevron в энергетическом секторе, а также BNP Paribas в секторе банковских услуг.

Благодаря системам Microsoft Windows 7 и Apple Snow Leopard, вычисления на GPU займут свои позиции в секторе массовых решений. В этих новых операционных системах GPU предстанет не только графическим процессором, но также и универсальным процессором для параллельных вычислений, работающим с любым приложением.

2.1.1 Модель программирования CUDA

Вычисления в CUDA производятся по принципу SIMD, где каждый процессор параллельно работает над разными данными, но выполняют одни инструкции. Графический процессор является вспомогательным устройством, помогающим основному центральному процессору. Графический процессор имеет свою собственную память и выполняет большее число потоков одновременно. Каждый поток на графическом процессоре выполняет свой экземпляр функции ядра (kernel).

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

Программирование в CUDA предполагает группирование потоков в блоки (thread block) -- в одномерные или двумерные сетки потоков, взаимодействующих между собой при помощи разделяемой памяти и точек синхронизации. В каждый момент времени исполняется одна сетка (grid) блоков потоков (см. рисунок 2.1.4), которая исполняет ядро (kernel). Блоки могут содержать в себе обновременно 512 потоков и могут быть одно-, двух- или трехмерными.

Рисунок 2.1.4 - Схема организации потоков на CUDA.



Блоки потоков, в свою очередь, объединяются в группы (warp). Размер группы- 32 потока, что является минимальным объём данных, который может обрабатываться на графическом ускорителе. Для удобства, CUDA позволяет работать и с блоками, состоящими из от 64 до 512 потоков.

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

2.1.2 Модель памяти CUDA

Адресация в CUDA позволяет использовать побайтовую адресацию, а также поддерживает как gather, так и scatter. Каждому потоку доступно большое число регистров, до 1024 штук. Работа с регистрами является очень быстрой. В регистрах можно хранить 32-битные целые или числа с плавающей точкой.

Имеется следующие типы памяти, доступные потокам:



Рисунок 2.1.5 - Схема доступа потоков к различным типам памяти CUDA.

Глобальная память -- самый медленный тип памяти, из доступных GPU. Глобальные переменные можно выделить с помощью спецификатора __global__, а так же динамически, с помощью функций из семейства cudMallocXXX. Глобальная память в основном служит для хранения больших объемов данных, поступивших на device с host'а, данное перемещение осуществляется с использованием функций cudaMemcpyXXX. В алгоритмах, требующих высокой производительности, количество операций с глобальной памятью необходимо свести к минимуму.

Локальная память -- это область в глобальной памяти, выделенная компилятором для хранения локальных значений потоков.

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

Память констант -- область памяти объемом 64 килобайта (то же -- для нынешних GPU), доступная только для чтения всеми мультипроцессорами. Она кэшируется по 8 килобайт на каждый мультипроцессор. Довольно медленная -- задержка в несколько сот тактов при отсутствии нужных данных в кэше.

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

Естественно, что глобальная, локальная, текстурная и память констант -- это физически одна и та же память, известная как локальная видеопамять видеокарты. Их отличия в различных алгоритмах кэширования и моделях доступа. Центральному процессору доступны только следующие типы памяти: глобальную, константную и текстурную.

Рисунок 2.1.6 - Схема зон видимости различных типов памяти CUDA.



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

Типичный, но не обязательный шаблон решения задач:

· задача разбивается на подзадачи;

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

· каждый блок обрабатывается блоком потоков;

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

· над данными в разделяемой памяти проводятся соответствующие вычисления;

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

2.1.3 CUDA и язык C

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

В состав CUDA входят runtime библиотеки:

· общая часть, предоставляющая встроенные векторные типы и подмножества вызовов RTL, поддерживаемые на CPU и GPU;

· CPU-компонента, для управления одним или несколькими GPU;

· GPU-компонента, предоставляющая специфические функции для GPU.

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

Для примера приведен кода для сложения векторов, представленный в CUDA (листинг 1).

Листинг 2.1 - пример кода на CUDA.

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

Этапы построения приложения CUDA изображены на рисунке 9.

Рисунок 2.1.7 - Стадии компиляции CUDA-приложения.

Для компиляции CUDA программ имеется специальный компилятор NVCC, который является оболочкой над другими инструментами: cudacc, g++, cl и др. NVCC генерирует машинный код для центрального процессора, который компилируется вместе с остальными частями приложения, написанными на чистом Си, и объектный код PTX для графического ускорителя. Выполнение программ CUDA требует наличия библиотек CUDA runtime library (cudart) и CUDA core library (cuda).

Язык для CUDA имеет ряд дополнительных расширений для языка C, которые необходимы для написания кода для графического ускортеля:

1. Спецификаторы функций, которые показывают, как и откуда буду выполняться функции.

2. Спецификаторы переменных, которые служат для указания типа используемой памяти GPU.

3. Спецификаторы запуска ядра GPU.

4. Встроенные переменные для идентификации нитей, блоков и др. параметров при исполнении кода в ядре GPU .

5. Дополнительные типы переменных.

Спецификаторы функций, определяющие как и откуда будут вызываться функции в CUDA существует три:

1. __host__ -- функция исполняется и вызывается с центрального процессора.

2. __global__ -- исполняется на графическом ускорителе, вызвывается с центрального процессора.

3. __device__ -- функция исполняется и вызывается с графического процесора.

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

myKernelFunc<<<gridSize, blockSize, sharedMemSize, cudaStream>>>(float* param1,float* param2)

, где

1. gridSize - размерность выделенной сетки блоков (dim3),

2. blockSize - размер выделенного блока (dim3)

3. sharedMemSize - размер дополнительной памяти, выделяемой при запуске ядра,

4. cudaStream - переменная cudaStream_t, задающая поток, в котором будет произведен вызов.

myKernelFunc - функция ядра (спецификатор __global__). Некоторые переменные при вызове ядра можно опускать, например sharedMemSize и cudaStream. Также в CUDA существует несколько встроенных переменных:

5. gridDim - размерность выделенной сетки, имеет тип dim3.

6. blockDim - размерность выделенного при текущем вызове блока, так же имеет тип dim3.

7. blockIdx - индекс текущего блока в вычислении на графическом ускорителе, имеет тип uint3.

8. threadIdx - индекс текущей нити в вычислении на графическом ускорителе, имеет тип uint3.

9. warpSize - размер warp'а, имеет тип int

2.1.4 Библиотека cuBLAS

cuBLAS -- реализация интерфейса программирования приложений для создания библиотек, выполняющих основные операции линейной алгебры BLAS (Basic Linear Algebra Subprograms) для CUDA. Он позволяет получить доступ к вычислительным ресурсам графических процессоров NVIDIA. Библиотека является самодостаточной на уровне API, то есть, прямого взаимодействия с драйвером CUDA не происходит. cuBLAS прикрепляется к одному GPU и автоматически не распараллеливается между несколькими GPU.

Основные функции библиотеки cuBLAS: создание матриц и векторных объектов в пространстве памяти GPU, заполнение их данными, вызов последовательных функций cuBLAS, и загрузка результатов из области памяти GPU обратно к хосту. Чтобы достичь этого, cuBLAS предоставляет вспомогательные функции для создания и уничтожения объектов в памяти GPU, и для записи данных и извлечения информации из этих объектов.

Для максимальной совместимости с существующими средами Fortran, cuBLAS использует хранения в столбцах и индексирование с 1. Так как C и C++ используют построчное хранение, приложения не могут использовать родные семантики массивов для двумерных массивов. Вместо этого, макросы или встроенные функции должны быть определены для того чтобы использовать матрицы при помощи одномерных массивов. Для Fortran'а код портирован на C механическим способом, при котором сохраняется индексирование с 1. В этом случае индекс массива из матричного элемента в строке i и в столбце j могут быть вычислены с помощью следующего макроса:

#define IDX2F(i,j,ld) ((((j)-1)*(ld))+((i)-1))

Здесь ld -- это размерности матрицы, в случае хранения в столбцах, является количеством строк. Для кода изначально написанного на С и С++, можно было бы использовать индексирование с 0, в этом случае макрос выглядит так:

#define IDX2C(i,j,ld) (((j)*(ld))+(i))

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

#include <stdio.h> #include <stdlib.h> #include <math.h> #include "CUBLAS.h" #define IDX2F(i,j,ld) ((((j)-1)*(ld))+((i)-1)) void modify (float *m, int ldm, int n, int p, int q, float alpha, float beta) { CUBLASSscal (n-p+1, alpha, &m[IDX2F(p,q,ldm)], ldm); CUBLASSscal (ldm-p+1, beta, &m[IDX2F(p,q,ldm)], 1); } #define M 6 #define N 5 int main (void) { int i, j; CUBLASStatus stat; float* devPtrA; float* a = 0; a = (float *)malloc (M * N * sizeof (*a)); if (!a) { printf ("host memory allocation failed"); return EXIT_FAILURE; } for (j = 1; j <= N; j++) { for (i = 1; i <= M; i++) { a[IDX2F(i,j,M)] = (i-1) * M + j; } } CUBLASInit(); stat = CUBLASAlloc (M*N, sizeof(*a), (void**)&devPtrA); if (stat != CUBLAS_STATUS_SUCCESS) { printf ("device memory allocation failed"); return EXIT_FAILURE; } CUBLASSetMatrix (M, N, sizeof(*a), a, M, devPtrA, M); modify (devPtrA, M, N, 2, 3, 16.0f, 12.0f); CUBLASGetMatrix (M, N, sizeof(*a), devPtrA, M, a, M); CUBLASFree (devPtrA); CUBLASShutdown(); for (j = 1; j <= N; j++) { for (i = 1; i <= M; i++) { printf ("%7.0f", a[IDX2F(i,j,M)]); } printf ("\n"); } return EXIT_SUCCESS; }

Листинг 2.1.2 - пример создания матрицы размерности 5 на 6.

2.1.5 Библиотека cuSPARSE

Библиотека линейной алгебры разреженных матриц для NVIDIA CUDA (cuSPARSE) представляет набор базовых подпрограмм линейной алгебры для разреженных матриц и обеспечивает до 8 раз более высокую производительность, чем последняя библиотека MKL. Библиотека cuSPARSE предназначена для вызовов из C или C++, и последняя версия библиотеки содержит треугольный решатель разреженных систем.

Программы библиотеки могут быть отнесены к четырем категориям:

1. уровень 1: операции между вектором в разреженном формате и вектором в плотном формате;

2. уровень 2: операции между матрицей в разреженном формате и вектором в плотном формате;

3. уровень 3: операции между матрицей в разреженном формате и набором векторов в плотном формате;

4. преобразование: операции, которые позволяют преобразование между разными форматами матриц.

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

API cuSPARSE предполагает, что входные и выходные данные размещены на устройстве (графическом ускорителе), до тех пор пока точно не указано иначе в параметре DevHostPtr в названии параметра функции (например, parameter*resultDevHostPtr в функции cusparse<t>doti()).

Разработчик сам отвечает за выделение памяти и копировании данных между памятью графического ускорителя и памятью центрального процессора используя стандартные операции CUDA API, такие как cudaMalloc(), cudaFree(), cudaMemcpy(), и cudaMemcpyAsync().

Заметьте: библиотека cuSPARSE требует устройство с вычислительной возможностью (compute capability) не ниже 2.0.

2.1.6 Библиотека cuSOLVER

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

cuSOLVER является высокоуровневой сборкой, основанной на библиотеках cuBLAS и cuSPARSE. Он объединяет три отдельные библиотеки в одну оболочку, каждая из которых может быть использована независимо или вместе с другими библиотеками.

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

Первая часть cuSolver называется cuSolverDN, и имеет дело с плотной матричной прогонкой и процедурами решения, такими как LU, QR, SVD и LDLT, а также полезными утилитами, такими как матричные и векторные перестановки.

Далее, cuSolverSP предоставляет новый набор редких процедур, основанных на разреженной QR факторизации. Не все матрицы имеют хорошую разреженности для параллелизма в факторизации, поэтому библиотека cuSolverSP также обеспечивает путь к центральному процессору для обработки этих последовательных матриц. Для тех матриц с обильной параллельностью, решение на графических ускорителях будет поставлять более высокую производительность. Библиотека предназначена для вызовов из C и C ++.

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

Путь GPU библиотеки cuSolver принимает данные уже находящиеся в памяти устройства. Это ответственность разработчика выделить память и копировать данные между памятью GPU и памяти CPU с использованием стандартных API CUDA во время выполнения процедуры, таких как cudaMalloc(), cudaFree(), cudaMemcpy(), и cudaMemcpyAsync().

2.2 Инструмент gSOAP

gSOAP - это инструмент разработки SOAP и REST XML веб-служб на языках C и C++, а также шаблонных привязок данных XML и C/C++. Данное средство анализирует WSDL и XML схемы (по отдельности или вместе) и переносит типы XML схемы и протокола сообщений SOAP/REST XML в простой и легкий в использовании эффективный C и C++ код. Также поддерживает преобразование C и C++ программ в XML веб-сервисы с помощью автогенерируемой XML сериализации кода и спецификации WSDL. Или же, вы можете просто использовать его для автоматического преобразования XML в/из C и C++ данных. gSOAP поддерживает возможность создавать классический ANSI C или C++ с или без STL.

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

gSOAP предлагает прозрачное решение для привязки данных C/C++ и XML через автоматически создаваемые коды. Это экономит время на разработку подобной системы. Дополнительно, использование привязки данных значительно упрощает использование XML в приложении путем автоматического сопоставления XML и C/C++ типов данных. Разработчики приложений больше не должны подгонять логику приложения к специфичным библиотекам.

Данное средство подходит к WS-I Basic Profile 1.0a, 1.1, и 1.2, прошедшим тесты W3C Databinding, и может оперировать с другими SOAP/REST сервисами, такими как WCF, .NET, и Axis.

Инструмент также реализует привязку данных ANSI C и C++ с XML через использование технологий компилирования

Эти технологии сопоставляют схемы XML с определениями C/C++ и наоборот. Существует два основных преимущества использования данного подхода. Во-первых, строгая типизация обеспечивает безопасное использование типов данных в вашем приложении, а также гарантирует, что данные прошли проверку XML. Во-вторых, схема, основанная на компиляторе специфичная на схеме парсинг и валидация является более эффективной, чем большинство техник парсинга.

Кроме того, gSOAP создает WSDL и XML схемы (XSD). Это сопоставляет существующие типы данных и функции C/C++ в WSDL, таким образом поддерживая и упрощая перенос старого кода с веб-сервисом. Переносимость кода была достигнута для многих платформ, таких как встраиваемые системы и для программ реального времени.

2.2.1 Пример использования

Установщик gSOAP включает два инструмента компилятора для разработки вашего приложения:

· Парсер WSDL "wsdl2h". Этот инструмент конвертирует WSDL и XSD файлы в документированный C/C++ определения.

· Компилятор заглушек и каркаса "soapcpp2". Данный компилятор создает RPC код и XML сериализаторы из документированных C/C++ определений.

Для начала нужно создать заголовочный файл C и описать в нем те функции и типы данных, которые будут использоваться создаваемой веб-службой. Этот файл должен быть описан в определённом формате, чтобы gSOAP смог его распознать.

Затем мы запускаем компилятор "soapcpp2" и подаем ему наш заголовочный файл, чтобы он создал нам исходный код нашего будущего клиента. Компилятор заглушки и каркаса "soapcpp2" создает прокси-класс (и RPC заглушку) для вашего клиентского приложения, как показано ниже.

Рисунок 2.2.1 - Схема создания клиента на gSOAP.



Библиотека gSOAP предоставляет транспортный уровень с HTTP стеком поверх стека TCP/IP с поддержкой SSL шифрования и прикреплением DIME/MIME.

Чтобы создать веб-службу, запустите обработчик "wsdl2h" и подайте ему ваш WSDL файл, который сформирует вам заголовочный файл. Затем этот файл, так же, как и для клиента, обработайте компилятором "soapcpp2".

Рисунок 2.2.2 - Схема создания веб-сервиса на gSOAP.

Компилятор "soapcpp2" генерирует каркас веб-службы на C/C++. Вы также можете взять старый код C/C++ и разработать веб-сервис просто добавив код веб-службы на C/C++ и типы данных в заголовочный файл. Компилятор "soapcpp2" сгенерирует исходный код вашего проекта и создаст WSDL для вашей веб-службы.

2.3 MATLAB

MATLAB - мультипарадигменная вычислительная среда и одноимённый язык программирования, используемый в этом пакете. MATLAB используют более 1 000 000 инженерных и научных работников [12], он работает на большинстве современных операционных систем, включая Linux, Mac OS и Microsoft Windows[13].

Проприетарный язык программирования, разработанный компанией MathWorks, MATLAB позволяет проводить манипуляции с матрицами, строить графики функций и данных, реализовывать алгоритмы, создавать пользовательский интерфейс и взаимодействовать с программами, написанными на других языках, включая C, C++, Java, Fortran и Python.

Хотя MATLAB разработан в основном для числовых операций, дополнительная панель управления MuPAD symbolic engine предоставляет доступ к возиожностям символьных операций. Дополнительный пакет Simulink добавляет графическое мульти-доменную симуляцию и основанную на моделях разработку для динамических и встроенных систем.

2.3.1 Синтаксис языка

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

Переменные определяются с помощью оператора присваивания "=". MATLAB является слабо-типизированным языком, потому что типы являются неявно конвертируемыми. Также он является языком с предполагаемой типизацией, потому что переменные могут быть объявлены без указания их типа, кроме тех случаев, когда они рассматриваются как символьные объекты, и затем их тип может быть изменен. Значения переменных могут быть получены из констант, вычислений, включающих эти же переменные, или со значения, возвращаемого функцией. Например:

>> x = 17

x =

17

>> x = 'hat'

x =

hat

>> y = x + 0

y =

104 97 116

>> x = [3*4, pi/2]

x =

12.0000 1.5708

>> y = 3*sin(x)

y =

-1.6097 3.0000

Листинг 2.3.1 - Пример объявления переменных в MATLAB.

Вектора и матрицы. Простой массив объявляется с помощью двоеточий: init:increment:terminator. Например:

>> array = 1:2:9

array =

1 3 5 7 9

Листинг 2.3.2 - Пример объявления массива в MATLAB.

Increment может быть опущен, тогда инкрементация будет с шагом 1.

>> ari = 1:5

ari =

1 2 3 4 5

Листинг 2.3.3 - Пример объявления массива в MATLAB.

Индексация начинается с 1, поскольку это обычная конвенция в математике.

Матрицы могут быть объявлены разделением элементов строки пробелами или запятыми и используя точку с запятой для отделения каждой строки. Список элементов следует заключать в квадратные строки: []. Скобки используются для доступа к элементу или подмассиву.

>> A = [16 3 2 13; 5 10 11 8; 9 6 7 12; 4 15 14 1]

A =

16 3 2 13

5 10 11 8

9 6 7 12

4 15 14 1

>> A(2,3)

ans =

11

Листинг 2.3.4 - Пример объявления матрицы в MATLAB.

Большинство функций MATLAB могут получить доступ к матрице и применять себя к каждому элементу. Например, mod(2 * J, n) умножит каждый элемент "J" на 2 и затем присвоит осток от деления на n. Следующий код, выдержанный из magic.m, создает магический квадрат M нечетными значениями n (функция MATLAB meshgrid здесь используется для создания квадратной матрицы I на J содержащей 1:n).

[J,I] = meshgrid(1:n);

A = mod(I + J - (n + 3) / 2, n);

B = mod(I + 2 * J - 2, n);

M = n * A + B + 1;

Листинг 2.3.5 - Пример применения функции к матрице в MATLAB.

Структуры. MATLAB имеет структурные типы данных. Поскольку все переменные в MATLAB массивы, более подходящее название "массив структур", в которых каждый элемент массива имеет поле с таким же названием. Дополнительно, MATLAB поддерживает поля с динамическими именами. К сожалению, MATLAB JIT не поддерживает структуры MATLAB, поэтому простое объединение различных переменных в структуру выйдет по стоимости.

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

MATLAB также поддерживает ООП.



3. Реализация

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

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

3.1 Веб-служба

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

Веб-службу можно разделить на два основных блока: блок вычисления коэффициентов сплайна Смоляка и блок передачи данных. Блок вычислений будет решать СЛАУ (1.4.8), описанный в разделе 1.4 и находить искомые коэффициенты сплайна Смоляка используя технологию CUDA. Затем эти коэффициенты будут переданы клиенту кодом, созданным с помощью инструмента gSOAP.