Вычисление БПФ на GPU с использованием Java и CUDA

В последнее время у меня на работе встала задача создания спектрограмм для довольно большого объема данных. Главным инструментом для этого является алгоритм Быстрого Преобразования Фурье. На сегодняшний день для Java написано уже множество библиотек, реализующих вычисление БПФ и другие математические инструменты, например, мне нравится Commons Math от Apache project. Но для того, чтобы обрабатывать большие объемы данных с помощью этой библиотеки за более или менее разумное время, нужна довольно мощная и недешевая машина, которой у меня нет. Поэтому я и обратил свой взор на технологию NVIDIA CUDA, которая позволяет очень шустро производить расчеты с плавающей точкой на чипе видеокарты. А, как известно, БПФ это очень даже плавающая точка :) Вдвойне приятно то, что хорошие люди уже давно написали Java-интерфейс JCuda для нативной библиотеки CUDA. И ещё больше я обрадовался, когда обнаружил, что библиотека CUDA имеет в своем составе модуль CUFFT, который предоставляет средства для расчета БПФ, использующие по полной программе возможности видеокарты для выполнения параллельных вычислений с плавающей точкой.

CUDA поддерживается видеокартами NVIDIA GeForce начиная с чипа G80, т.е. с восьмой серии. Машина, на которой я тестирую в данный момент связку Java + CUDA довольно старенькая и имеет на борту процессор AMD Athlon X2 4200+, видеокарту NVIDIA GeForce 8800 GTS 640 Mb, всего 2 Gb оперативной памяти и управляется операционной системой openSUSE 11.4.

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

Процессоры Модель Гигафлопс AMD Phenom II X6 1100T Black Edition 3.3 ГГц 60.1 Intel Core i7-975 XE 3.33 ГГц 53.28 Intel Core 2 Quad Q8300 2.5 ГГц 40 AMD AMD ATHLON II X4 645 3.1 ГГц 38.44 Intel Core 2 Duo 2.4 ГГц 19.2 AMD Athlon 64 X2 4200+ 2.2 ГГц 13.2

Видеокарты Модель Гигафлопс GeForce GTX 590 2488.3 GeForce GTX 580 1581.056 GeForce GTX 480 1344.96 GeForce GTX 280 933.120 GeForce 9800 GTX 648 GeForce 8800 GTS 640mb 346

Как видно из таблиц, теоретическая производительность моей видеокарты в 26 раз превосходит производительность моего процессора и очень уж хочется это превосходство использовать при расчете БПФ. Посмотрим как это сделать.

На данном этапе я пока не разбирался глубоко в возможностях CUDA и использовал только готовый модуль CUFFT и интерфейс к нему из JCuda. Т.е. в моем случае мне не приходилось писать свои собственные ядра (kernels) для CUDA. Что такое ядро CUDA? Это файл с расширением PTX, исходный код которого пишется на языке, являющимся по сути ограниченной версией языка C, а затем компилируется с помощью специального компилятора NVCC. Файл ядра CUDA содержит функции, которые должны выполняться на GPU. С помощью JCuda можно, грубо говоря, загружать такие ядра в вашу java-программу и вызывать функции, описанные в них.

Подготовка

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

1. Скачать и установить "родные" проприетарные драйвера NVIDIA.
2. Скачать и установить CUDA Toolkit
3. Скачать и распаковать библиотеку JCuda
4. Добавить путь до jar-файлов JCuda в CLASSPATH или подключить её в вашей IDE как пользовательскую библиотеку
5. Установить путь до нативных библиотек JCuda (файлы *.so в linux или *.dll в windows) в переменной окружения LD_LIBRARY_PATH

Важные замечания

1. Сначала я установил Developer Drivers for Linux (260.19.26), ссылка на которые находится на той же странице, где и ссылка на CUDA Toolkit 3.2. Уж не знаю является ли установка именно этой версии драйверов обязательной, но после того, как я убедился, что CUDA работает, я установил последние на тот момент драйвера NVIDIA версии 280.13, которые при установке сообщили о том, что сносят драйвера 260.19.26, потому что сами всё умеют. CUDA при этом осталась работоспособной.
2. У меня 32-разрядная версия openSUSE, поэтому мне были нужны соответствующие версии CUDA и JCuda. Не особо запариваясь, я взял последнюю доступную собранную для Linux 32-bit версию JCuda 0.3.2a и соответствующую ей версию CUDA Toolkit 3.2. CUDA Toolkit 3.2 была доступна для версии openSUSE 11.2, но вполне себе благополучно встала на 11.4. Входе установки CUDA Toolkit, предложит выбрать путь установки и по завершению скажет о том, что нужно будет добавить в переменные окружения PATH и LD_LIBRARY_PATH. Сделайте это.
3. В Eclipse JCuda можно подключить как пользовательскую библиотеку: Window -> Preferences -> Java -> Build Path -> User Libraries, New, AddJARs... и добавить все файлы *.jar из каталога jcuda. После этого библиотеку можно добавить в проект.
4. Необходимо указать путь к нативным библиотекам JCuda. В Eclipse можно поступить двумя способами. Первый способ заключается в том, чтобы при добавлении библиотеки JCuda в User Libraries, для каждого jar-файла указать путь до каталога, в котором лежат нативные библиотеки JCuda (Native Library Location). Второй способ это создать каталог в проекте eclipse и скопировать в него все файлы нативных библиотек JCuda. В Run Configurations проекта добавить переменную окружения LD_LIBRARY_PATH и указать в ней путь до этого каталога. На скриншоте ниже второй способ показан наглядно.

От слов к делу

CUFFT базируется на очень популярной и эффективной библиотеке FFTW. Довольно интересная концепция этой библиотеки заключается в том, что прежде чем осуществить преобразование какого-либо массива данных, нужно сформировать так называемый план (plan). Такой подход позволяет не тратить время на пересчет конфигурации алгоритма преобразования при множественных однотипных преобразованиях. Чтобы создать план используя JCuda, сначала нужно создать объект cufftHandle, который представляет указатель на план, а потом, используя этот указатель, создать план с помощью какой-либо из функций: cufftPlan1d, cufftPlan2d или cufftPlan3d. Эти методы создают план для выполнения одномерного, двумерного и трехмерного БПФ соответственно. В моем случае необходимо осуществить одномерное (1D) преобразование вещественных чисел и на выходе получить также вещественные числа. В CUFFT есть несколько вариантов проведения БПФ:

• CUFFT_C2C - преобразование комплексных чисел в комплексные (одинарная точность)
• CUFFT_C2R - преобразование комплексных чисел в вещественные (одинарная точность)
• CUFFT_D2Z - преобразование вещественных чисел в комплексные (двойная точность)
• CUFFT_R2C - преобразование вещественных чисел в комплексные (одинарная точность)
• CUFFT_Z2D - преобразование комплексных чисел в вещественные (двойная точность)
• CUFFT_Z2Z - преобразование комплексных чисел в комплексные (двойная точность)

В Java тип с одинарной точностью это float, а с двойной - double. Массив комплексных чисел в CUFFT представляет чередующиеся реальные и мнимые части числа. Т.е. для представления одного комплексного числа требуется два элемента массива. Например, массив размером 10, может содержать 5 комплексных чисел.

Как видно выше, в CUFFT нет возможности подать на вход вещественные числа и на выходе получить тоже вещественные, поэтому можно использовать преобразования CUFFT_R2C или CUFFT_D2Z, а затем самостоятельно найти модули полученных комплексных чисел, это уже не такое дорогое действие как выполнение БПФ в плане времени.

Еще оговорюсь, что, к сожалению, моя GeForce 8800 GTS не поддерживает операции с двойной точностью, поэтому мне приходится использовать CUFFT_R2C, хотя хотелось бы CUFFT_D2Z, но это уже мелочи которые лечатся покупкой новой видеокарты.

Пример создания плана для одномерного БПФ вещественных чисел в комплексные:

cufftHandle plan = new cufftHandle();
JCufft.cufftPlan1d(plan, inputDataSize, cufftType.CUFFT_R2C, 1);

• plan - указатель на план
• inputDataSize - количество входных значений; если входные данные являются вещественными, то этот параметр равняется длине входного массива, а если входные данные это массив комплексных чисел, то значение параметра равно половине длины входного массива, т.к. одно комплексное число представлено двумя элементами массива
• cufftType.CUFFT_R2C - тип преобразования
• 1 - количество подобных преобразований

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

// создание случайных входных данных размером 1024
float[] inputData = createRandomData(1024);
// массив для сохранения результатов преобразования
float[] fftResults = new float[inputData.length + 2];
// выполнение преобразования
JCufft.cufftExecR2C(plan, inputData, fftResults);

В этом коде можно заметить некоторую странность. Почему размер массива для сохранения результатов преобразования на 2 элемента больше, чем размер входных данных? Важно в этом разобраться. На входе мы имеем массив из N = 1024 вещественных чисел. А результатом у нас будут комплексные числа, т.е. по 2 элемента массива на одно число. Но в тоже время в документации CUFFT сказано, что результатом преобразования является массив без избыточных значений и его размер равен N/2 + 1. Вот тут для меня небольшая загадка. Почему +1? Зачем он нужен? На сколько я знаю, полезными результатами БПФ является массив размером равным ровно половине размера входного массива, т.е. N/2. Поэтому этот +1 элемент может является особенностью алгоритма, а может это пробел в моих знаниях :) Но это и не важно. Важно помнить то, что результатами преобразований CUFFT является массив который содержит N/2 + 1 значений. Но почему тогда в примере выше размер выходного массива равен не N/2 + 1, а N + 2? Это потому, что на выходе у нас комплексные числа, и для их представления, как было уже много раз упомянуто выше, необходимо два элемента массива. Если писать полностью, то правильно было бы вычислять размер выходного массива так: (inputData.length/2 + 1) * 2. Если раскрыть скобки, то мы и получим inputData.length + 2.

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

// создание случайных входных данных размером 1024
float[] inputData = createRandomData(1024);
// массив для сохранения результатов преобразования
float[] fftResults = new float[inputData.length + 2];
// указатель на устройство
Pointer deviceDataIn = new Pointer();
// выделение памяти на устройстве для входных данных
JCuda.cudaMalloc(deviceDataIn, inputData.length * 4);
// копирование данных из оперативной памяти в память видеокарты
JCuda.cudaMemcpy(deviceDataIn, Pointer.to(inputData), inputData.length * 4, 
          cudaMemcpyKind.cudaMemcpyHostToDevice);
     
Pointer deviceDataOut = new Pointer();
// выделение памяти на устройстве для результатов
JCuda.cudaMalloc(deviceDataOut, fftResults.length * 4);

// создание плана     
cufftHandle plan = new cufftHandle();
JCufft.cufftPlan1d(plan, inputData.length, cufftType.CUFFT_R2C, 1);

// преобразование
JCufft.cufftExecR2C(plan, deviceDataIn, deviceDataOut);

// копирование данных из памяти видеокарты в оперативную память     
JCuda.cudaMemcpy(Pointer.to(fftResults), deviceDataOut, fftResults.length * 4, 
             cudaMemcpyKind.cudaMemcpyDeviceToHost);

В данном примере метод JCuda.cudaMalloc(Pointer devPtr, long size) осуществляет выделение памяти размером size байт на видеокарте и возвращает на неё указатель в параметре devPtr. JCuda.cudaMemcpy(Pointer dst, Pointer src, long count, int cudaMemcpyKind_kind) осуществляет копирование данных между хостом и видеокартой, где dst - указатель на место назначения, src - указатель на источник, count - количество байт для копирования, cudaMemcpyKind_kind - направление копирования (например, cudaMemcpyKind.cudaMemcpyDeviceToHost для копирования данных с памяти видеокарты в оперативную память). С помощью Pointer.to(float[] values) создается указатель на массив значений.

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

JCufft.cufftDestroy(plan);
JCuda.cudaFree(deviceDataIn);
JCuda.cudaFree(deviceDataOut);

Пример программы

Теперь приведу пример рабочей программы, которая будет генерировать массив случайных значений размером 2 в 23 степени (кратность степени двойки необходима для библиотеки от Apache, для CUFFT это условие не обязательно) и производить БПФ этого массива с помощью Apache Commons Math и CUDA, замеряя при этом время выполнения каждого преобразования. В конце программа выведет первое, среднее и последнее значение для каждого из результатов преобразования, чтобы можно было убедиться в правильности расчетов. Для того, чтобы программа успешно выполнилась, необходимо увеличить размер кучи у JVM с помощью параметров командной строки -Xms512m -Xmx768m или же уменьшить размер генерируемого входного массива. В Eclipse указать параметры командной строки можно в Run Configurations во вкладке Arguments, раздел VM arguments.

import java.util.Date;
import java.util.Random;
import org.apache.commons.math.complex.Complex;
import org.apache.commons.math.transform.FastFourierTransformer;
import jcuda.Pointer;
import jcuda.jcufft.JCufft;
import jcuda.jcufft.cufftHandle;
import jcuda.jcufft.cufftType;
import jcuda.runtime.JCuda;
import jcuda.runtime.cudaMemcpyKind;

public class JCufftDemo {

 public static void main(String[] args) {

  double[] fftResults;  
  int dataSize = 1<<23;

  System.out.println("Генерация входных данных размером "+dataSize+" значений...\n");
  float[] inputData = createRandomData(dataSize);
  
  System.out.println("1D БПФ с использованием apache commons math...");
  fftResults = commonsTransform(floatDataToDoubleData(inputData.clone()));
  printSomeValues(fftResults);
  
  System.out.println();
  System.out.println("1D БПФ JCufft (данные в оперативной памяти)...");
  fftResults = jcudaTransformHostMemory(inputData.clone());
  printSomeValues(fftResults);
   
  System.out.println();
  System.out.println("1D БПФ JCufft (данные в памяти видеокарты)...");
  fftResults = jcudaTransformDeviceMemory(inputData.clone());
  printSomeValues(fftResults);
 }
 
 /**
  * Генерирует массив случайных чисел
  * 
  * @param dataSize - размер генерируемого массива
  * @return массив случайных чисел
  */
 public static float[] createRandomData(int dataSize){
        Random random = new Random();
        float data[] = new float[dataSize];
        
        for (int i = 0; i < dataSize; i++)
         data[i] = random.nextFloat();

        return data;
 }
 
 /**
  * Конвертирует массив значений типа float в массив значений double
  * 
  * @param data - массив который нужно конвертировать
  * @return - сконвертированный массив
  */
 public static double[] floatDataToDoubleData(float[] data){ 
  
  double[] doubleData = new double[data.length];  
  for(int i=0; i < data.length; i++) doubleData[i] = data[i];
  
  return doubleData;
 }
 
 /**
  * Выполняет БПФ массива значений с помощью CUDA, осуществляя
  * операции с данными в оперативной памяти
  * 
  * @param inputData - массив входных значений
  * @return массив с результатами БПФ
  */
 public static double[] jcudaTransformHostMemory(float[] inputData){
  float[] fftResults = new float[inputData.length + 2];
     // создание плана
     cufftHandle plan = new cufftHandle();
     JCufft.cufftPlan1d(plan, inputData.length, cufftType.CUFFT_R2C, 1);
     // выполнение БПФ
     long timeStart = new Date().getTime();
     JCufft.cufftExecR2C(plan, inputData, fftResults);
     System.out.println("Время преобразования: " + (new Date().getTime() - timeStart)/1000.0+" сек");
     // уничтожение плана
     JCufft.cufftDestroy(plan);

  return cudaComplexToDouble(fftResults);
 }
 
 /**
  * Выполняет БПФ массива значений с помощью CUDA, осуществляя
  * операции с данными в памяти видеокарты
  * 
  * @param inputData - массив входных значений
  * @return массив с результатами БПФ
  */
 public static double[] jcudaTransformDeviceMemory(float[] inputData){

     float[] fftResults = new float[inputData.length + 2];

     // указатель на устройство
     Pointer deviceDataIn = new Pointer();
     // выделение памяти на видеокарте для входных данных
     JCuda.cudaMalloc(deviceDataIn, inputData.length * 4);
     // копирование входных данных в память видеокарты
     JCuda.cudaMemcpy(deviceDataIn, Pointer.to(inputData), inputData.length * 4, 
          cudaMemcpyKind.cudaMemcpyHostToDevice);
     
     Pointer deviceDataOut = new Pointer();
     // выделение памяти на видеокарте для результатов преобразования
     JCuda.cudaMalloc(deviceDataOut, fftResults.length * 4);

     // создание плана
     cufftHandle plan = new cufftHandle();
     JCufft.cufftPlan1d(plan, inputData.length, cufftType.CUFFT_R2C, 1);

     // выполнение БПФ
     long timeStart = new Date().getTime();
     JCufft.cufftExecR2C(plan, deviceDataIn, deviceDataOut);
     System.out.println("Время преобразования: " + (new Date().getTime() - timeStart)/1000.+" сек");
     
     // копирование результатов из памяти видеокарты в оперативную память
     JCuda.cudaMemcpy(Pointer.to(fftResults), deviceDataOut, fftResults.length * 4, 
             cudaMemcpyKind.cudaMemcpyDeviceToHost);

     // освобождение ресурсов
     JCufft.cufftDestroy(plan);
     JCuda.cudaFree(deviceDataIn);
     JCuda.cudaFree(deviceDataOut);

  return cudaComplexToDouble(fftResults);
 }
 
 /**
  * Выполняет БПФ массива значений с помощью Apache Commons Math
  * 
  * @param inputData - массив входных значений
  * @return массив с результатами БПФ
  */
 public static double[] commonsTransform(double[] inputData){
  
  FastFourierTransformer fft = new FastFourierTransformer();
  long timeStart = new Date().getTime();
  Complex[] cmx = fft.transform(inputData);
  System.out.println("Время преобразования: " + (new Date().getTime() - timeStart)/1000.+" сек");

  double[] fftReults = new double[inputData.length/2 + 1];
  for(int i = 0; i < fftReults.length; i++){
   fftReults[i] = cmx[i].abs();
  }
  
  return fftReults;
 }
 
    /**
     * Метод осуществляет преобразование массива комплексных чисел в
     * массив, содержащий их модули
     * 
     * @param complexData - массив комплексных чисел
     * @return массив модулей комплексных чисел
     */
    public static double[] cudaComplexToDouble(float[] complexData){

     double[] result = new double[complexData.length/2];
     int j=0;
        for(int i=0; i < complexData.length-1; i++) {         
         result[j++] = Math.sqrt(complexData[i]*complexData[i] + complexData[i+1]*complexData[i+1]);
         i++;
        }

     return result;  
    }
    
    /**
     * Выводит на стандартный вывод первое, среднее 
     * и последнее значения массива
     * 
     * @param data - массив, значения которого необходимо вывести
     */
    public static void printSomeValues(double[] data){    
     System.out.println("[0]: "+data[0]);
     System.out.println("["+(data.length/2)+"]: "+data[data.length/2]);
     System.out.println("["+(data.length-1)+"]: "+data[data.length-1]);
    }
}

Результат выполнения программы:

Генерация входных данных размером 8388608 значений...

1D БПФ с использованием apache commons math...

Время преобразования: 92.851 сек

[0]: 4195107.645827353

[2097152]: 1411.8392665442454

[4194304]: 114.60465306043625

1D БПФ JCufft (данные в оперативной памяти)...

Время преобразования: 0.085 сек

[0]: 4195107.922955976

[2097152]: 1411.791370918522

[4194304]: 114.75

1D БПФ JCufft (данные в памяти видеокарты)...

Время преобразования: 0.0 сек

[0]: 4195107.922955976

[2097152]: 1411.791370918522

[4194304]: 114.75

Как видно из результатов выполнения программы, выполнение БПФ с помощью CUDA на несколько порядков быстрее, чем при использовании Apache Commons Math, а при использовании памяти видеокарты оно вообще происходит меньше миллисекунды, что,  согласитесь, просто "сносит башню"! :) Но тут необходимо принять во внимание тот факт, что в программе измеряется только выполнение самого преобразования и не учитываются временные издержки на нахождение модуля комплексных данных, выделение памяти под массивы, копирование данных между видеокартой и хостом и т.п., а эти издержки имеют место быть. Поэтому если у вас есть задача в которой вам пригодилось бы использование CUFFT, то лучше сначала провести исследование стоит ли вам использовать память видеокарты или нет. Однако, если нужна двойная точность, то вариантов кроме как копировать данные на видеокарту и обратно нет. Кстати, по поводу двойной точности. В результате выполнения программы можно заметить, что результаты преобразование с помощью Apache Commons Math несколько отличаются от результатов CUFFT. Это происходит как раз из-за того, что в Commons используется двойная точность, а в CUFFT мне пришлось использовать одинарную по причине староватости видеокарты.