Курсовая работа (т): Изучение возможностей массивно-параллельных вычислений в применении к задачам математической физики

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

4.2 Быстрое преобразование Фурье и библиотека CUFFT

В набор разработчика CUDA Toolkit, по обыкновению, входят две дополнительные библиотеки, связанные с наиболее часто встречающимися задачами. Первая из них, CUBLAS - реализация стандарта BLAS (Basic Linear Algebra Subprograms) - содержит в себе фундаментальные алгоритмы осуществления базовых операций линейной алгебры. О второй, CUFFT, и пойдёт речь в этом подразделе.

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

) операции с комплексными и вещественными данными;

) одномерные, двумерные и трёхмерные преобразования;

) выполнение преобразований над массивами размером до шестидесяти четырёх миллионов элементов с одинарной точностью и до ста двадцати восьми миллионов элементов с двойной точностью;

) поддержка in-place преобразований, при которых используется минимум места для хранения данных, а входной массив замещается результирующим, и out-of-place преобразований;

) двойная точность вычислений на некоторых совместимых видеокартах;

) поддержка потокового выполнения, позволяющего несинхронные вычисления и перемещение данных /17/.

Вообще быстрое преобразование Фурье (Fast Fourier Transform, FFT) - собирательное название группы алгоритмов вычисления дискретного преобразования Фурье. Главной отличительной чертой этих алгоритмов является значительное уменьшение числа шагов, необходимых для расчёта коэффициентов, по сравнению со стандартными способами преобразования. Впервые такой метод был предложен Джеймсом У.Кули из Исследовательского центра имени Томаса Уотсона корпорации IBM и Джоном У.Тьюки из Bell Telephone Laboratories /3/. Основная идея их работы заключалась в том, чтобы разделить сумму из  слагаемых на две суммы с числом слагаемых, равным , а затем вычислить обе суммы по отдельности. В свою очередь вычисление этих сумм производилось тем же самым способом, иначе говоря, рекурсивно. Постоянное деление оставшегося отрезка пополам в конце концов приводило к двухточечному преобразованию Фурье. Несложно заметить, что длина входного вектора  для реализации такого алгоритма должна быть степенью двойки, однако это не накладывает особого ограничения на входные данные, ибо всегда существует возможность дополнить вектор с  элементами до необходимой длины, добавив к нему  нулевых коэффициентов, наличие которых не повлияет на результат.

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

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

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

Очевидно, что коэффициенты  и  можно представить в виде:

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

то есть линейную комбинацию двух дискретных преобразований Фурье над «половинными» векторами. Продолжая рекурсивное разделение, приходим к двухточечному преобразованию, которое, с учётом выражения (4.2.4) и того факта, что значение  равно минус единице, будет вычисляться по формулам:

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

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

Интерфейс программирования CUFFT API (Application Programming Interface) базируется на свободно распространяемой библиотеке FFTW, которая считается одной из наиболее эффективных библиотек быстрого преобразования Фурье в среде CPU-программирования. FFTW предоставляет простой механизм конфигурации, так называемый plan, который полностью предопределяет оптимальную схему выполнения расчета для определённого размера и типа входных данных. При вызове функции вычисления преобразования производятся точно по предписанному плану. Преимущество такого подхода заключается в том, что единожды создав подобную схему, система сохраняет состояние, необходимое для множественного выполнения преобразований без пересчёта конфигурации.

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

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

где  - указатель на массив входных данных, созданный CPU. Далее нужно объявить объект типа :

в котором будет сохранён план, и определить его:

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

)  - преобразование вещественных данных в комплексные;

)  - преобразование комплексных данных в комплексные;

)  - комплексно-вещественное преобразование,

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

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

Здесь  и  - указатели на массивы комплексных данных в памяти GPU,  - исходный сигнал и  - полученный. В нашем случае роль этих аргументов будет играть преобразованная к типу указателя на  переменная  (выполняем in-place преобразование).

Параметром  типа  можно задать направление преобразования, его значение  даёт прямое преобразование Фурье,  - обратное.

Всё, что осталось сделать - это скопировать данные из памяти GPU в память, доступную CPU, «убить» план функцией  и освободить занимаемое место на видеокарте:

На этом знакомство с быстрым преобразованием Фурье и основным функциональным аппаратом библиотеки CUFFT будем считать законченным и перейдём непосредственно к их использованию в поставленной задаче.

4.3 Реализация FFT для задачи о колебаниях струны и мембраны

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

Начиная реализацию FFT для случая струны, прежде всего зададим функцию её начального положения:

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