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

Министерство образования и науки Российской Федерации Федеральное государственное бюджетное образовательное учреждение высшего профессионального образования «Самарский государственный аэрокосмический университет имени академика С. П. Королева (национальный исследовательский университет)»

Факультет информатики Кафедра технической кибернетики

Выпускная квалификационная работа специалиста на тему: Разработка и реализация алгоритма решения систем линейных алгебраических уравнений трёхдиагонального вида на графическом вычислительном устройстве

Самара, 2013

Реферат

Объектом исследования является метод циклической редукции и алгоритм Томаса для решения СЛАУ с матрицей трехдиагонального вида.

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

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

Развитие математики (теория разностных схем, проекционные и вариационные методы и другое), физики (оптика, теория теплопроводности, газовая динамика и другое) и моделирования приводит к необходимости разработки новых способов решения различных задач вычислительной физики. При исследовании некоторых из указанных задач возникает потребность в рассмотрении трехдиагональных систем линейных алгебраических уравнений (СЛАУ).

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

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

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

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

CUDA (англ. Compute Unified Device Architecture) — технология GPGPU (англ. General-Purposecomputing on Graphics Processing Units), позволяющая программистам реализовывать на упрощённом языке программирования Си алгоритмы, выполнимые на графических процессорах ускорителей GeForce восьмого поколения и старше.

Часто при решении больших трехдиагональных систем на GPU, возникает проблема нехватки памяти («разделяемая память»). Поэтому, в рамках данной дипломной работы был разработан параллельный метод решения трехдиагональных систем линейных алгебраических уравнений (СЛАУ), решающий эту проблему.

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

1. Методы решения систем линейных уравнений трехдигонального вида

1.1 Метод прогонки

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

Ax=f (1.1)

с ленточной матрицей A. В этом случае уравнения всегда могут быть упорядочены так, что каждое неизвестное появляется только в нескольких уравнениях, «соседних» с i-м уравнением. Формально мы говорим, что матрица имеет верхнюю ширину ленты q, если для всех , и нижнюю ширину ленты р, если для всех [8]. При q=p=1 такая матрица называется трехдиагональной:

.

В системе (1.1) — вектор неизвестных, — вектор правых частей.

Данную систему можно переписать в виде:

(1.2)

где

Пусть требуется найти решение этой системы уравнений.

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

Рассмотрим алгоритм решения системы вида (1.2). Проведем исключение неизвестных в (1.2) для этого запишем систему в следующем виде:

(1.3)

где

Произведем замену в первом уравнении системы (1.3):

; .

Возьмем первые два уравнения системы (1.3):

Умножим первое уравнение на и сложим со вторым уравнением получим:

.

Разделим полученное уравнение на, тогда будем иметь:

.

Произведем замену в уравнении:

.

Тогда получим:

.

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

(1.4)

где

Система (1.4) имеет аналогичную системе (1.2) структуру и не содержит неизвестное. Если решить эту систему, то потом неизвестное можно будет найти по следующей формуле:

.

К системе (1.4) можно опять применить описанный выше способ исключения неизвестных. Таким образом, на втором шаге будет исключено неизвестное, на третьем и так далее. В результате l-го шага придем к системе для неизвестных, , …, :

(1.5)

где

А так же получим формулы для нахождения ():

(1.6)

где .

Коэффициенты и находятся по следующим формулам:

где

Полагая в системе (1.5), получим систему для и :

Из этой системы найдем и :

Объединяя эти равенства с равенством (1.6) для случая, получим окончательные формулы для нахождения неизвестных:

(1.7)

где ;

; (1.8)

. (1.9)

Формулы (1.7−1.9) описывают метод прогонки. Коэффициенты и называются прогоночными коэффициентами, формулы (1.8, 1.9) описывают прямой ход прогонки, а (1.7) — обратный ход.

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

;

.

Подсчет арифметический действий в (1.7) — (1.9) показывает, что реализация метода прогонки по этим формулам требует арифметических действий[1].

В рассмотренном методе правой прогонки определение производится последовательно справа налево.

Общий последовательный алгоритм правой прогонки состоит из следующей последовательности шагов:

Шаг 1. Задать начальные значения коэффициентов a, b, c и правых частей уравнений f.

Шаг 2. Рассчитать для прямого хода прогонки и .

Шаг 3. Рассчитать значения прогоночных коэффициентов и

Шаг 4. Найти значение .

Шаг 5. Используя формулу найти остальные неизвестные.

Аналогично выписываются формулы левой прогонки .

Рассчитываются прогоночные коэффициенты:

; (1.10)

. (1.11)

И неизвестные:

(1.12)

Общий последовательный алгоритм левой прогонки состоит из шагов аналогичных шагам правой прогонки.

1.2 Метод встречных прогонок

Иногда оказывается удобно комбинировать правую и левую прогонки, получая так называемый метод встречных прогонок. Этот метод целесообразно применять, если надо найти только одного неизвестное или группу идущих подряд неизвестных. В этом методе в области вычисляются прогоночные коэффициенты, по формулам (1.8, 1.9) и в области коэффициенты, по формулам (1.10, 1.11).

Выпишем формулы (1.7) и (1.12) для обратного хода правой и левой прогонок для. Получим систему:

Из которой найдем :

.

Используя найденное, по формуле (1.6) для найдем последовательно, а по формуле (1.12) для вычислим остальные.

Итак, формулы метода встречных прогонок имеют вид:

; (1.8)

. (1.9)

; (1.10)

. (1.11)

для прогоночных коэффициентов и для определения решения.

1.3 Метод циклической редукции

Вернемся к СЛАУ вида (1.1) и рассмотрим метод циклической редукции применительно к указанной системе.

Перепишем систему (1.1) в виде (1.12):

.

Идея метода циклической редукции состоит в последовательном исключении из уравнений системы (1.12) неизвестных сначала с нечетными номерами i, затем из оставшихся уравнений с номерами i, кратными 2, затем 4 и так далее.

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

(1.13)

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

Тогда получаем:

.

Произведем следующие замены:

Тогда получаем:

. (1.14)

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

.

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

.

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

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

(1.15)

где с шагом до и, когда они имеют место.

Диаграмма маршрутизации этого алгоритма изображена на рисунке Рисунок 2 для. Для удобства обозначим вектор .

Рисунок 2 — Диаграмма маршрутизации алгоритма циклической редукции решение линейный уравнение трехдигональный Общий последовательный алгоритм циклической редукции состоит из следующей последовательности шагов:

Шаг 1. Задать начальные значения коэффициентов a, b, c и правых частей уравнений f.

Шаг 2. Пересчитать для каждого шага редукции значения коэффициентов a и c.

Шаг 3. Пересчитать для каждого шага редукции значения правых частей уравнений f.

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

Шаг 5. Пересчитать для каждого шага редукции значения оставшихся неизвестных.

2. Параллельные алгоритмы решения систем линейных алгебраических уравнений трехдиагонального вида

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

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

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

2.1 Параллельный алгоритм встречных прогонок

Рассмотрим двухзадачный алгоритм, основанный на методе встречных прогонок.

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

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

Рисунок 3а — Прямой ход прогонок Рисунок 3б — Обмен данными между задачами Рисунок 3в — Обратный ход прогонок В течение прямого хода прогонок первый процессор находит прогоночные коэффициенты, , второй процессор —,. Далее, первый процессор передает второму пару чисел, , необходимых для запуска обратного хода прогонки на втором процессоре, который в это время производит аналогичную передачу, необходимую первому процессору. На третьем шаге оба процессора одновременно реализуют обратный ход прогонок, находя решение СЛАУ.

Ускорение такого алгоритма можно оценить величиной [4]:

;

где M — размерность СЛАУ, С₁M_а — длительность расчета по последовательному алгоритму, если _а — время производства одной арифметической операции с плавающей точкой, _к=_у+_п — длительность одной коммуникации между процессорами (_у — время установки на обмен, _п — время передачи или приема (считаем их равными) одного пакета по сети). При пакетной модели передачи данных время передачи разного объема информации остается константой до тех пор, пока объем передаваемого массива не превысит размер пакета. Здесь это условие соблюдается (передается всего пара чисел), в дальнейшем также будем полагать его верным. К сожалению, этот подход нельзя развить для большего количества задач (алгоритм не масштабируем [10]), если сеточная область одномерна.

2.2 Метод декомпозиции области

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

Ax=b;

где А — трехдиагональная матрица размера nЧn;

n=pq+p-1 для некоторых целых p и q (p>q).

Алгоритм основан на разбиении множества неизвестных на непересекающиеся подмножества D_i и множество разделителей S_i, как это показано на рисунке 4:

Рисунок 4 — Разбиение множества неизвестных Затем переменные переупорядочиваются путем присвоения последних индексов элементам множества S, что приводит к системе стреловидного вида:

(2.1)

Вводя следующие обозначения A_l=diag (A₁…A_p), B^T=(B₁…B_p), C=(C₁…C_p), x_l=(x₁…x_p), b_l=(b₁…b_p), описанную выше систему можно представить в виде:

(2.2)

(2.3)

Предположив, что матрица A_l невырожденная, умножим уравнение (2.3) на C и вычтем полученное из уравнениия (2.2):

(2.4)

где

Как только будет решена система (2.4) относительно x_s, остальные неизвестные можно найти из систем:

A_lx_i= b_l - Bx_s;

где .

Данные системы можно решать независимо друг от друга. Кроме того, слагаемые CB и Cb также можно находить независимо. Схема алгоритма декомпозиции областей представлена на рисунке 5:

Рисунок 5 — Схема алгоритма декомпозиции области Ускорение алгоритма оценивается величиной:

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

2.3 Параллельный алгоритм правой прогонки

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

Применим блочный подход к разделению данных: пусть каждый поток обрабатывает m=n/p строк матрицы A, т. е. k-й поток обрабатывает строки с номерами. кратно числу потоков, в общем случае изменится только число уравнений в последнем потоке. Ниже, на рисунке 6, представлено разделение данных для трех потоков в случае системы из 12 уравнений.

Рисунок 6 — Разделение данных для трех потоков системы из 12 уравнений

В пределах полосы матрицы, обрабатываемой k-м потоком, можно организовать исключение поддиагональных элементов матрицы (прямой ход метода). Для этого осуществляется вычитание строки i, умноженной на константу, из строки с тем, чтобы результирующий коэффициент при неизвестной в ()-й строке оказался нулевым. Если исключение первым потоком поддиагональных переменных не добавит в матрицу новых коэффициентов, то исключение поддиагональных элементов в остальных потоках приведет к возникновению столбца отличных от нуля коэффициентов: во всех блоках (кроме первого) число ненулевых элементов в строке не изменится, но изменится структура уравнений.

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

Рисунок 7 — Иллюстрация прямого хода параллельной прогонки.

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

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

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

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

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

2.4 Параллельный алгоритм циклической редукции

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

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

Диаграмма маршрутизации для параллельного варианта показана на рисунке 12. Для удобства обозначим вектор .

Рисунок 12 — Диаграмма маршрутизации алгоритма параллельного циклической редукции На каждом уровне уравнения редукции для всех п уравнений решаются параллельно. Трудность возникает тогда, когда требуются данные, которые не попадают в описанный диапазон . Однако можно увидеть, что правильный результат получается, если принять

(2.5)

для и, или .

Добавив эти специальные значения в изначальную систему получим уравнение

для и; (2.6)

которое дает правильные значения при попадании в систему (1.13). Теперь можем рассматривать либо решение исходного конечного множества уравнений, либо, как вариант, решение расширенного бесконечного множества с коэффициентами (2.5). Это равнозначно только суммирующим уравнениям, таким как (2.6), вне диапазона исходной задачи. Любая точка зрения равно справедлива, но последняя больше подходит для параллельного варианта циклической редукции, так как она описывает необходимые значения и . вне первоначально обозначенной задачи. После подсчета (фактически требуются только значения и ). Значение х получается из уравнения (1.15):

.

Члены х в правой части уравнения (1.15) отсутствуют, так как на этом уровне редукции они лежат вне диапазона и по уравнению (2.2) равны нулю.

Число операций в алгоритме параллельной циклической редукции составляет:

;

где ,.

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

3. Разработка алгоритма для реализации на GPU

В предложенном алгоритме предполагается проведение вычислений, на центральном процессорном устройстве (CPU) и на графическом процессорном устройстве (GPU), то есть, на так называемой гетерогенной вычислительной среде. Вычисления, которые можно произвести независимо, будут выполняться на GPU. Участки, которые нельзя распараллелить, будут выполняться на СPU. Запускать параллельные участки кода на графическом устройстве позволяет технология CUDA.

3.1 Основные возможности и особенности технологии CUDA

Технология CUDA — это программно-аппаратная вычислительная архитектура NVIDIA, основанная на расширении языка Си, которая даёт возможность организации доступа к набору инструкций графического ускорителя и управления его памятью при организации параллельных вычислений. CUDA помогает реализовывать алгоритмы, выполнимые на графических процессорах видеоускорителей GeForce восьмого поколения и старше (серии GeForce 8, GeForce 9, GeForce 200), а также Quadro и Tesla.

Хотя трудоёмкость программирования GPU при помощи CUDA довольно велика, она ниже, чем с ранними GPGPU решениями. Такие программы требуют разбиения приложения между несколькими мультипроцессорами подобно MPI программированию, но без разделения данных, которые хранятся в общей видеопамяти. И так как CUDA программирование для каждого мультипроцессора подобно OpenMP программированию, оно требует хорошего понимания организации памяти.

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

Модель программирования в CUDA предполагает группирование потоков. Потоки объединяются в блоки потоков (thread block) — одномерные или двумерные сетки потоков, взаимодействующих между собой при помощи разделяемой памяти и точек синхронизации. Программа (ядро, kernel) исполняется над сеткой (grid) блоков потоков (thread blocks). Эта модель представлена на рисунке 13.

Рисунок 13 — Модель программирования в CUDA

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

Каждый поток имеет доступ к шести типам памяти, что отображено на рисунке 14.

Рисунок 14 — Виды памяти в CUDA

Глобальная память — самый большой объём памяти, доступный для всех мультипроцессоров на видеочипе, размер составляет от 256 мегабайт до 1.5 гигабайт на текущих решениях. Обладает высокой пропускной способностью, более 100 гигабайт/с для топовых решений NVIDIA, но очень большими задержками в несколько сот тактов. Не кэшируется, поддерживает обобщённые инструкции, и обычные указатели на память.

Локальная память — это небольшой объём памяти, к которому имеет доступ только один потоковый процессор. Она относительно медленная — такая же, как и глобальная.

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

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

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

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

Рисунок 15 — Модель памяти в CUDA

3.2 Реализация разработанного алгоритма на GPU

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

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

Разработанный алгоритм^* продемонстрирован схематически на рисунке 16 для случая, когда в системе 8 уравнений.

Рисунок 16 — Схема разработанного алгоритма На входе мы имеем одну трехдиагональную систему линейных уравнений, состоящую из восьми уравнений. На первом этапе данная система разбивается на две подсистемы, состоящих из четырех уравнений, по строению эти системы аналогичны первой. При необходимости разбиение можно продолжать до того момента, когда подсистемы будут помещаться в разделяемую память GPU. На последнем этапе полученные системы решаются с помощью метода прогонки.

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

Первое ядро осуществляет редукцию, разбивая систему до размеров, которые умещаются в разделяемую память GPU. Редукция осуществляется до уровня. При этом должно соответствовать условию:

;

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

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

;

где — количество блоков, количество уравнений в системе, максимально возможное количество нитей в одном блоке, а общее количество выделенных блоков составит:

;

где общее количество выделенных блоков, а количество решаемых систем.

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

.

На рисунке 17 изображена схема алгоритма.

Рисунок 17 — Схема реализации алгоритма на вычислительной системе

Ниже изложен алгоритм реализации разработанного метода с использованием одного ГПУ.

1.Инициализация исходных массивов данных.

2. Выделение глобальной памяти на GPU под массивы коэффициентов a, b, c, f и массив неизвестных x.

3. Передача массивов коэффициентов a, b, c, f на GPU.

4. Вызов ядра. Разделение системы на более мелкие.

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

6.Передача результирующего x вектора на CPU.

Чтобы функции выполнялись не процессором, а устройством, они имеют квалификатор __global__.

Объявление функции ядра, осуществляющего алгоритм редукции, может иметь следующий вид:

__global__ void reduction (float *a, float *b, float *c, float *f) {

int thread = blockIdx. x*N+threadIdx.x+1;

int block = blockIdx. x;

for (int i=0; i

int min = thread-pow (2.0, i);

int max = thread+pow (2.0, i);

if (min

if (max>(block+1)*N) max=0;

float alfa=-a[thread]/b[min];

float betta=-c[thread]/b[max];

float zn0=alfa*a[min];

float zn1=b[thread]+alfa*c[min]+betta*a[max];

float zn2=betta*c[max];

float zn3=f[thread]+alfa*f[min]+betta*f[max];

__threadfence ();

a[thread]=zn0;

b[thread]=zn1;

c[thread]=zn2;

f[thread]=zn3;

}

}

Вызов ядра из основной программы осуществляется следующим образом:

reduction<<>>(dev_a, dev_b, dev_c, dev_f).

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

— выделение памяти на устройстве:

cudaMalloc ((void**) &dev_a, (N*M+1)*sizeof (float));

— освобождения памяти на устройстве:

cudaFree (dev_a);

— пересылка данных с ЦПУ на ГПУ:

cudaMemcpy (dev_a, a,(N*M+1)*sizeof (float), cudaMemcpyHostToDevice);

— пересылка данных с ГПУ на ЦПУ:

cudaMemcpy (a, dev_a,(N*M+1)*sizeof (float), cudaMemcpyDeviceToHost);

4. Исследования алгоритмов

Исследования алгоритмов будем проводить на вычислительной системе с графическим вычислительным устройством GeForce GT 540M, ЦПУ Intel Core i5−2410М 2,3 Ггц, операционной системой Microsoft Windows 7 с установленным драйвером NVIDIA CUDA Version 4.0.

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

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

Так как в разработанном алгоритме решения трехдиагональных систем линейных алгебраических уравнений для разбиения больших систем на малые используется метод циклической редукции, проведем эксперименты для разных степеней разбиения системы. Зафиксируем значения M и N и будем изменять степень разбиения системы k. Результат запишем в виде таблицы и графика (для лучшей наглядности)

4.1 Исследования ускорения алгоритма

Проведем исследования алгоритма, зафиксируем количество систем M, а количество уравнений N будем изменять. В таблице 1 представлены результаты исследования отношения времени выполнения вычислений последовательного алгоритма на CPU (T_CPU) ко времени выполнения вычислений разработанного алгоритма на GPU (T_GPU) от размеров сеточной области. Число систем уравнений равняется 300, 500 и 700 для одних и тех же N .

Таблица 1

Для иллюстрации приведенных в таблице 6 данных на рисунке 18 приведен график зависимости ускорения S алгоритма от .

Рисунок 18

На графике, изображенном на рисунке 18, видно, что при небольшой размерности систем (до 2³) ускорение меньше единицы, то есть разработанный алгоритм для GPU работает медленнее, чем последовательный. Это связано с затратами времени на выделение памяти на GPU и последующим копирование в нее исходных данных копирование превышают выигрыш в скорости выполнения арифметических операций.

На больших же размерах системы алгоритм эффективен, так как ускорение больше единицы. Из графика видно, что ускорение увеличивается с ростом числа уравнений. Но при размерности 2¹⁰ и 2¹³ ускорение начинает колебаться около девяти. Это объясняется ограничением роста производительности вычислительной системы с увеличением количества вычислителей. Согласно этому ограничению, ускорение выполнения программы за счёт распараллеливания её инструкций на множестве вычислителей ограничено временем, необходимым для выполнения её последовательных инструкций (закон Амдала). Кроме того из графиков видно, что чем больше число систем, тем быстрее растет ускорение при увеличении количества уравнений в системе. Это можно объяснить тем, что алгоритм построен таким образом, что количество задействованных блоков в ядре напрямую зависит от числа систем в задаче, а с ростом числа задействованных блоков растет и производительность.

Теперь зафиксируем количество уравнений N, а изменять будем количество систем M. В таблице 1 представлены результаты исследования отношения времени выполнения вычислений по алгоритму на CPU (T_CPU)ко времени выполнения вычислений на GPU (T_GPU) от размеров сеточной области. Число уравнений в системе равняется 512, 1024 и 2048 для одних и тех же M.

Таблица 2

Для иллюстрации приведенных в таблице 2 данных на рисунке 19 приведен график зависимости ускорения S алгоритма от M.

Рисунок 19

Из рисунка 19 видны границы эффективного применения алгоритма. Можно сделать выводы, что при увеличении количества систем отношение времени выполнения вычислений последовательного алгоритма на CPU (T_CPU) ко времени выполнения вычислений разработанного алгоритма на GPU (T_GPU) возрастает с увеличением числа решаемых систем. При количестве систем до 100, ускорение меньше единицы. А значит, параллельный алгоритм работает медленнее, чем последовательный, что связано с неполной загруженностью. На больших же количествах систем алгоритм эффективен, так как ускорение больше единицы. Из графика видно, что при количестве систем больше 350, ускорение начинает колебаться около 9. Это объясняется ограничением роста производительности вычислительной системы.

В разработанном алгоритме для разбиениея систем используется циклическая редукция. В таблице 3 представлены результаты исследования отношения времени выполнения вычислений по алгоритму на CPU (T_CPU) ко времени выполнения вычислений на GPU (T_GPU) в зависимости от степени разбиения систем.

Таблица 3

Для иллюстрации приведенных в таблице 3 данных на рисунке 20 приведен график зависимости ускорения S алгоритма от .

Рисунок 20

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

Заключение

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

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

Для алгоритма, в ходе экспериментов, было установлено, что при размерности систем N в пределах от 2⁵ до 2¹³ ускорение изменяется от 1 до 3,5. Так же было установлено, что при количестве систем M в пределах от 100 до 600 ускорение изменяется от 1 до 3,5. Это говорит о том, что в данных диапазонах размерностей сеточной области лучше воспользоваться параллельным алгоритмом на GPU, чем последовательным на CPU. Также было выяснено, что при N>2¹⁰ M>400 ускорение начинает колебаться около 3,5. Было выяснено, что разбиение системы оказывает негативное влияние на ускорение разработанного алгоритма, однако оно необходимо для использования разделяемой памяти при больших размерностях систем.

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

1 Самарский, А. А. Методы решения сеточных уравнений [Текст]/А.А. Самарский, Е. С. Николаев. — М.: Наука, 1978. — 592 с.

2 Хокни Р., Джессхоуп К. Параллельные ЭВМ. Архитектура, программирование и алгоритмы: Пер. с англ. — М.: Радио и связь, 1986. — 392

3 Головашкин, Д. Л. Методы параллельных вычислений. Часть 2. [Текст]: учебное пособие/Д.Л. Головвкашкин, С. П. Головашкина. — Самар. гос. аэрокосм. ун-т. — Самара, 2003. — 103 с.

4 Ортега, Дж.

Введение

в параллельные и векторные методы решения линейных систем/ Дж.Ортега. — М.: Мир, 1991. — 364 с.

5 Ярмушкин, С. В. Исследование параллельных алгоритмов решения трехдиагональных систем линейных уравнений [Текст]/С.В. Ярмушкин, Д.Л. Головашкин//Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки. — 2004. — вып. 26. — С. 78−82.

6 Hockney R.W., Jesshope C.R. Parallel Computers 2. Architecture, Programming and Algorithms. — Adam Hilger, Bristol and Philadelphia, 1988. — 625 p.

7 Головашкин, Д. Л. Параллельные алгоритмы решения сеточных уравнений трехдиагонального вида, основанные на методе встречных прогонок [Текст]/Д.Л. Головашкин // Математическое моделирование. — 2005. — Том 17, № 11. — С. 118−128.

8 Голуб, Дж. Матричные вычисления [Текст]/Дж. Голуб, Ч. Ван Лоун. — М.: Мир, 1999. — 548 с.

9 Бахвалов, Н. С. Численные методы [Текст]: учебное пособие для физмат. специальностей вузов/Н.С. Бахвалов, Н. П. Жидков, Г. М. Кобельков. — М.: Физматлит, 2002. — 630 с.

10 Головашкин, Д. Л. Параллельные алгоритмы метода циклической прогонки [Текст]/Д.Л. Головашкин, М.В. Филатов//Компьютерная Оптика. — ИСОИ РАН, Самараю — 2005. — вып. 27. — С. 123−130.

Приложение

Исходный код программы

#include

#include

#include

#include

#include «CTimer.h»

#include

#define M (500)

#define N (1024)

#define K (0)

#define D (0)

__global__ void reduction (float *a, float *b, float *c, float *f) {

int thread = (blockIdx.x/int (pow (2.0, D)))*N*int (pow (2.0, D))+threadIdx.x+1;

int block = blockIdx. x/int (pow (2.0, D));

for (int i=0; i

int min = thread-pow (2.0, i);

int max = thread+pow (2.0, i);

if (min

if (max>(block+1)*N*int (pow (2.0, D))) max=0;

float alfa=-a[thread]/b[min];

float betta=-c[thread]/b[max];

float zn0=alfa*a[min];

float zn1=b[thread]+alfa*c[min]+betta*a[max];

float zn2=betta*c[max];

float zn3=f[thread]+alfa*f[min]+betta*f[max];

__threadfence ();

a[thread]=zn0;

b[thread]=zn1;

c[thread]=zn2;

f[thread]=zn3;

}

__syncthreads ();

}

__global__ void sweep (float *a, float *b, float *c, float *f, float *test, int shareDim) {

float temp[N*M];

int block = blockIdx. x;

int first=(block/int (pow (2.0, K)))*N+(block%int (pow (2.0, K))+1);

int xDim = N/int (pow (2.0, K));

int redDim = int (pow (2.0, K));

for (int i=0; i

int next = first+i*redDim;

temp[i]=a[next];

temp[i+xDim ]=b[next];

temp[i+2*xDim ]=c[next];

temp[i+3*xDim ]=f[next];

}

temp[xDim*4]=-temp[2*xDim]/temp[xDim];

temp[xDim*4+xDim-1]=temp[3*xDim]/temp[xDim];

for (int i=1; i

temp[xDim*4+i]=-temp[2*xDim+i]/(temp[i]*temp[xDim*4+i-1]+temp[xDim+i]);

temp[xDim*4+xDim-1+i]=(temp[3*xDim+i]-temp[i]*temp[xDim*4+xDim-1+i-1])/(temp[i]*temp[xDim*4+i-1]+temp[xDim+i]);

}

int i =xDim-1;

temp[i]=(temp[3*xDim+i]-temp[i]*temp[xDim*4+xDim-1+i-1])/(temp[i]*temp[xDim*4+i-1]+temp[xDim+i]);

for (int i=xDim-2; i>-1; i—){

temp[i]=temp[i+1]*temp[xDim*4+i-1+1]+temp[xDim*4+xDim-1+i-1+1];

}

if (block==0){

for (int i=0; i< shareDim; i++) {

test[i]=temp[i];

}

}

syncthreads ();

}

int main (void){

CTimer Time;

CTimer Time1;

int shareDim = N/int (pow (2.0, K))*4+(N/int (pow (2.0, K))-1)*2;

float* a = (float*)malloc ((M*N+1)*sizeof (float));

float* b = (float*)malloc ((M*N+1)*sizeof (float));

float* c = (float*)malloc ((M*N+1)*sizeof (float));

float* f = (float*)malloc ((M*N+1)*sizeof (float));

float* test = (float*)malloc ((M*N)*sizeof (float));

a[0]=0; b[0]=1; c[0]=0; f[0]=0;

for (int i=1; i

a[i]=(rand ()%5 + 1);

}

for (int i=1; i

b[i]=(rand ()%9 + 10);

}

for (int i=1; i

c[i]=(rand ()%5 + 1);

}

for (int i=1; i

f[i]=(rand ()%9 + 1);

}

for (int i=0; i

a[1+i*N]=0; c[N+i*N]=0;

}

float* a1 = (float*)malloc ((M*N)*sizeof (float));

float* b1 = (float*)malloc ((M*N)*sizeof (float));

float* c1 = (float*)malloc ((M*N)*sizeof (float));

float* f1 = (float*)malloc ((M*N)*sizeof (float));

for (int i=1; i

a1[i-1]=a[i];

c1[i-1]=c[i];

b1[i-1]=b[i];

f1[i-1]=f[i];

}

//cudaEvent_t start3, stop3;

//cudaEventCreate (&start3);

//cudaEventCreate (&stop3);

float *dev_a, *dev_b, *dev_c, *dev_f, *dev_test;

cudaMalloc ((void**) &dev_test, (M*N)*sizeof (float));

cudaMemcpy (dev_test, test,(M*N)*sizeof (float), cudaMemcpyHostToDevice);

cudaMalloc ((void**) &dev_a, (N*M+1)*sizeof (float));

cudaMalloc ((void**) &dev_b, (N*M+1)*sizeof (float));

cudaMalloc ((void**) &dev_c, (N*M+1)*sizeof (float));

cudaMalloc ((void**) &dev_f, (N*M+1)*sizeof (float));

cudaMemcpy (dev_a, a,(N*M+1)*sizeof (float), cudaMemcpyHostToDevice);

cudaMemcpy (dev_b, b,(N*M+1)*sizeof (float), cudaMemcpyHostToDevice);

cudaMemcpy (dev_c, c,(N*M+1)*sizeof (float), cudaMemcpyHostToDevice);

cudaMemcpy (dev_f, f,(N*M+1)*sizeof (float), cudaMemcpyHostToDevice);

reduction<<>>(dev_a, dev_b, dev_c, dev_f);

Time.Start ();

sweep<<>>(dev_a, dev_b, dev_c, dev_f, dev_test, shareDim);

cudaThreadSynchronize ();

Time.End ();

cudaMemcpy (test, dev_test,(shareDim)*sizeof (float), cudaMemcpyDeviceToHost);

cudaMemcpy (a, dev_a,(N*M+1)*sizeof (float), cudaMemcpyDeviceToHost);

cudaMemcpy (b, dev_b,(N*M+1)*sizeof (float), cudaMemcpyDeviceToHost);

cudaMemcpy (c, dev_c,(N*M+1)*sizeof (float), cudaMemcpyDeviceToHost);

cudaMemcpy (f, dev_f,(N*M+1)*sizeof (float), cudaMemcpyDeviceToHost);

printf («%lf «, test[0]);

printf («n»);

printf («———————————n»);

float* x1 = (float*)malloc ((N*M)*sizeof (float));

float* alfa = (float*)malloc ((N-1)*sizeof (float));

float* betta = (float*)malloc ((N-1)*sizeof (float));

Time1.Start ();

for (int flag=0; flag

alfa[0]=-c1[flag*N]/b1[flag*N];

betta[0]=f1[flag*N]/b1[flag*N];

for (int j=1; j

alfa[j]=-c1[flag*N+j]/(alfa[+j-1]*a1[flag*N+j]+b1[flag*N+j]);

betta[j]=(f1[flag*N+j]-a1[flag*N+j]*betta[j-1])/(alfa[j-1]*a1[flag*N+j]+b1[flag*N+j]);

}

int j=N-1;

x1[flag*N+j]=(f1[flag*N+j]-a1[flag*N+j]*betta[j-1])/(alfa[j-1]*a1[flag*N+j]+b1[flag*N+j]);

for (int i=N-2; i>-1; i—)

x1[flag*N+i]=alfa[i]*x1[flag*N+i+1]+betta[i];

}

Time1.End ();

//———————————————————————————————————;

/*

printf («n»);

for (int i=0; i

printf («%lf—-%lf », x1[i*2], test[i]);

printf («n»);

}

*/

printf («————-n»);

printf («n»);

std:cout << Time. GetTimeInSeconds () << «seconds» << std: endl;

printf («n»);

std:cout << Time1. GetTimeInSeconds () << «seconds» << std: endl;

printf («n»);

printf («%lf », Time1. GetTimeInSeconds ()/Time.GetTimeInSeconds ());

//printf («### Total CPU time: t%.6fnn», 1.0*ev_time/CLOCKS_PER_SEC);

free (a);

free (b);

free (c);

free (f);

free (a1);

free (b1);

free (c1);

free (f1);

free (test);

free (x1);

free (alfa);

free (betta);

system («pause»);

}