CUDA: Работа с памятью. Часть I.

В процессе работы с CUDA я практически не касался вопросов об использовании памяти видеокарты. Настало время убрать этот пробел.

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

Видеокарта и типы памяти

При использовании GPU разработчику доступно несколько видов памяти: регистры, локальная, глобальная, разделяемая, константная и текстурная память. Каждая из этих типов памяти имеет определенное назначение, которое обуславливается её техническими параметрами (скорость работы, уровень доступа на чтение и запись). Иерархия типов памяти представлена на рис. 1.


Рис. 1. Типы памяти видеокарты

1. Регистровая память (register) является самой быстрой из всех видов. Определить количество регистров доступных GPU можно с помощью уже хорошо известной функции cudaGetDeviceProperties. Рассчитать количество регистров, доступных одной нити GPU, так же не составляет труда, для этого необходимо разделить общее число регистров на произведение количества нитей в блоке и количества блоков в гриде. Все регистры GPU 32 разрядные. В CUDA нет явных способов использования регистровой памяти, всю работу по размещению данных в регистрах берет на себя компилятор.
2. Локальная память (local memory) может быть использована компилятором при большом количестве локальных переменных в какой-либо функции. По скоростным характеристикам локальная память значительно медленнее, чем регистровая. В документации от nVidia рекомендуется использовать локальную память только в самых необходимых случаях. Явных средств, позволяющих блокировать использование локальной памяти, не предусмотрено, поэтому при падении производительности стоит тщательно проанализировать код и исключить лишние локальные переменные.
3. Глобальная память (global memory) – самый медленный тип памяти, из доступных GPU. Глобальные переменные можно выделить с помощью спецификатора __global__, а так же динамически, с помощью функций из семейства cudMallocXXX. Глобальная память в основном служит для хранения больших объемов данных, поступивших на device с host’а, данное перемещение осуществляется с использованием функций cudaMemcpyXXX. В алгоритмах, требующих высокой производительности, количество операций с глобальной памятью необходимо свести к минимуму.
4. Разделяемая память (shared memory) относиться к быстрому типу памяти. Разделяемую память рекомендуется использовать для минимизации обращение к глобальной памяти, а так же для хранения локальных переменных функций. Адресация разделяемой памяти между нитями потока одинакова в пределах одного блока, что может быть использовано для обмена данными между потоками в пределах одного блока. Для размещения данных в разделяемой памяти используется спецификатор __shared__.
5. Константная память (constant memory) является достаточно быстрой из доступных GPU. Отличительной особенностью константной памяти является возможность записи данных с хоста, но при этом в пределах GPU возможно лишь чтение из этой памяти, что и обуславливает её название. Для размещения данных в константной памяти предусмотрен спецификатор __constant__. Если необходимо использовать массив в константной памяти, то его размер необходимо указать заранее, так как динамическое выделение в отличие от глобальной памяти в константной не поддерживается. Для записи с хоста в константную память используется функция cudaMemcpyToSymbol, и для копирования с device’а на хост cudaMemcpyFromSymbol, как видно этот подход несколько отличается от подхода при работе с глобальной памятью.
6. Текстурная память (texture memory), как и следует из названия, предназначена главным образом для работы с текстурами. Текстурная память имеет специфические особенности в адресации, чтении и записи данных. Более подробно о текстурной памяти я расскажу при рассмотрении вопросов обработки изображений на GPU.

Пример использования разделяемой памяти

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

Перед тем, как приступить к написанию основного кода, приведу небольшой способ отладки. Как известно, функции из CUDA runtime API могут возвращать различные коды ошибок, но в предыдущий раз я ни как это не учитывал. Чтобы упростить себе жизнь можно использовать следующий макрос для отлова ошибок:

#define CUDA_DEBUG

#ifdef CUDA_DEBUG

#define CUDA_CHECK_ERROR(err)           \
if (err != cudaSuccess) {          \
printf("Cuda error: %s\n", cudaGetErrorString(err));    \
printf("Error in file: %s, line: %i\n", __FILE__, __LINE__);  \
}                 \

#else

#define CUDA_CHECK_ERROR(err)

#endif

* This source code was highlighted with Source Code Highlighter.

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

Приступаем к основной задаче.

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

// Функция транспонирования матрицы без использования разделяемой памяти
//
// inputMatrix - указатель на исходную матрицу
// outputMatrix - указатель на матрицу результат
// width - ширина исходной матрицы (она же высота матрицы-результата)
// height - высота исходной матрицы (она же ширина матрицы-результата)
//
__global__ void transposeMatrixSlow(float* inputMatrix, float* outputMatrix, int width, int height)
{
  int xIndex = blockDim.x * blockIdx.x + threadIdx.x;
  int yIndex = blockDim.y * blockIdx.y + threadIdx.y;

  if ((xIndex < width) && (yIndex < height))
  {
    //Линейный индекс элемента строки исходной матрицы  
    int inputIdx = xIndex + width * yIndex;  

    //Линейный индекс элемента столбца матрицы-результата
    int outputIdx = yIndex + height * xIndex;

    outputMatrix[outputIdx] = inputMatrix[inputIdx];
  }
}

* This source code was highlighted with Source Code Highlighter.

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

Пишем функцию транспонирования, которая использует разделяемую память:

#define BLOCK_DIM 16

// Функция транспонирования матрицы c использования разделяемой памяти
//
// inputMatrix - указатель на исходную матрицу
// outputMatrix - указатель на матрицу результат
// width - ширина исходной матрицы (она же высота матрицы-результата)
// height - высота исходной матрицы (она же ширина матрицы-результата)
//
__global__ void transposeMatrixFast(float* inputMatrix, float* outputMatrix, int width, int height)
{
  __shared__ float temp[BLOCK_DIM][BLOCK_DIM];

  int xIndex = blockIdx.x * blockDim.x + threadIdx.x;
  int yIndex = blockIdx.y * blockDim.y + threadIdx.y;

  if ((xIndex < width) && (yIndex < height))
  {
    // Линейный индекс элемента строки исходной матрицы  
    int idx = yIndex * width + xIndex;

    //Копируем элементы исходной матрицы
    temp[threadIdx.y][threadIdx.x] = inputMatrix[idx];
  }

  //Синхронизируем все нити в блоке
  __syncthreads();

  xIndex = blockIdx.y * blockDim.y + threadIdx.x;
  yIndex = blockIdx.x * blockDim.x + threadIdx.y;

  if ((xIndex < height) && (yIndex < width))
  {
    // Линейный индекс элемента строки исходной матрицы  
    int idx = yIndex * height + xIndex;

    //Копируем элементы исходной матрицы
     outputMatrix[idx] = temp[threadIdx.x][threadIdx.y];
  }
}

* This source code was highlighted with Source Code Highlighter.

В этой функции я использую разделяемую память в виде двумерного массива.
Как уже было сказано, адресация разделяемой памяти в пределах одного блока одинакова для всех потоков, поэтому, чтобы избежать коллизий при доступе и записи, каждому элементу в массиве соответствует одна нить в блоке.
После копирования элементов исходной матрицы в буфер temp, вызывается функция __syncthreads. Эта функция синхронизирует потоки в пределах блока. Её отличие от других способов синхронизации заключаеться в том, что она выполняеться только на GPU.
В конце происходит копирование сохраненных элементов исходной матрицы в матрицу-результат, в соответствии с правилом транспонирования.
Может показаться, что эта функция должна выполняться медленне, чем её версия без разделяемой памяти, где нет никаких посредников. Но на самом деле копирование из глобальной памяти в глобальную работает значительно медленее, чем связка глобальная память – разделяемая память – глобальная память.
Хочу заметить, что проверять границы массивов матриц стоит вручную, в GPU нет аппаратных средств для слежения за границами массивов.

Ну и напоследок напишем функцию транспонирования, которая исполняется только на CPU:

// Функция транспонирования матрицы, выполняемая на CPU
__host__ void transposeMatrixCPU(float* inputMatrix, float* outputMatrix, int width, int height)
{
  for (int y = 0; y < height; y++)
  {
    for (int x = 0; x < width; x++)
    {
      outputMatrix[x * height + y] = inputMatrix[y * width + x];
    }
  }
}

* This source code was highlighted with Source Code Highlighter.

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

#define GPU_SLOW 1
#define GPU_FAST 2
#define CPU 3

#define ITERATIONS 20    //Количество нагрузочных циклов

__host__ int main()
{  
  int width = 2048;    //Ширина матрицы
  int height = 1536;    //Высота матрицы

  int matrixSize = width * height;
  int byteSize = matrixSize * sizeof(float);

  //Выделяем память под матрицы на хосте
  float* inputMatrix = new float[matrixSize];      
  float* outputMatrix = new float[matrixSize];

  //Заполняем исходную матрицу данными
  for (int i = 0; i < matrixSize; i++)
  {
    inputMatrix[i] = i;
  }

  //Выбираем способ расчета транспонированной матрицы
  printf("Select compute mode: 1 - Slow GPU, 2 - Fast GPU, 3 - CPU\n");
  int mode;
  scanf("%i", &mode);

  //Записываем исходную матрицу в файл
  printMatrixToFile("before.txt", inputMatrix, width, height);
  
  if (mode == CPU)    //Если используеться только CPU
  {    
    int start = GetTickCount();        
    for (int i = 0; i < ITERATIONS; i++)
    {
      transposeMatrixCPU(inputMatrix, outputMatrix, width, height);
    }
    //Выводим время выполнения функции на CPU (в миллиекундах)
    printf ("CPU compute time: %i\n", GetTickCount() - start);
  }  
  else  //В случае расчета на GPU
  {
    float* devInputMatrix;
    float* devOutputMatrix;

    //Выделяем глобальную память для храния данных на девайсе
    CUDA_CHECK_ERROR(cudaMalloc((void**)&devInputMatrix, byteSize));
    CUDA_CHECK_ERROR(cudaMalloc((void**)&devOutputMatrix, byteSize));

    //Копируем исходную матрицу с хоста на девайс
    CUDA_CHECK_ERROR(cudaMemcpy(devInputMatrix, inputMatrix, byteSize, cudaMemcpyHostToDevice));

    //Конфигурация запуска ядра
    dim3 gridSize = dim3(width / BLOCK_DIM, height / BLOCK_DIM, 1);
    dim3 blockSize = dim3(BLOCK_DIM, BLOCK_DIM, 1);

    cudaEvent_t start;
    cudaEvent_t stop;
    
    //Создаем event'ы для синхронизации и замера времени работы GPU
    CUDA_CHECK_ERROR(cudaEventCreate(&start));
    CUDA_CHECK_ERROR(cudaEventCreate(&stop));

    //Отмечаем старт расчетов на GPU
    cudaEventRecord(start, 0);

    if (mode == GPU_SLOW)    //Используеться функция без разделяемой памяти
    {
      for (int i = 0; i < ITERATIONS; i++)
      {
      
        transposeMatrixSlow<<<gridSize, blockSize>>>(devInputMatrix, devOutputMatrix, width, height);      
      }
    }
    else if (mode == GPU_FAST)  //Используеться функция с разделяемой памятью
    {
      for (int i = 0; i < ITERATIONS; i++)
      {
      
        transposeMatrixFast<<<gridSize, blockSize>>>(devInputMatrix, devOutputMatrix, width, height);      
      }
    }

    //Отмечаем окончание расчета
    cudaEventRecord(stop, 0);

    float time = 0;
    //Синхронизируемя с моментом окончания расчетов
    cudaEventSynchronize(stop);
    //Рассчитываем время работы GPU
    cudaEventElapsedTime(&time, start, stop);

    //Выводим время расчета в консоль
    printf("GPU compute time: %.0f\n", time);

    //Копируем результат с девайса на хост
    CUDA_CHECK_ERROR(cudaMemcpy(outputMatrix, devOutputMatrix, byteSize, cudaMemcpyDeviceToHost));

    //
    //Чистим ресурсы на видеокарте
    //

    CUDA_CHECK_ERROR(cudaFree(devInputMatrix));
    CUDA_CHECK_ERROR(cudaFree(devOutputMatrix));

    CUDA_CHECK_ERROR(cudaEventDestroy(start));
    CUDA_CHECK_ERROR(cudaEventDestroy(stop));
  }

  //Записываем матрицу-результат в файл
  printMatrixToFile("after.txt", outputMatrix, height, width);

  //Чистим память на хосте
  delete[] inputMatrix;
  delete[] outputMatrix;

  return 0;
}

* This source code was highlighted with Source Code Highlighter.

В случае если расчеты выполняются только на CPU, то для замера времени расчетов используется функция GetTickCount(), которая подключается из windows.h. Для замера времени расчетов на GPU используеться функция cudaEventElapsedTime, прототип которой имеет следующий вид:

cudaError_t cudaEventElapsedTime( float* time, cudaEvent_t start, cudaEvent_t end ), где

1. time – указатель на float, для записи времени между event’ами start и end (в миллисекундах),
2. start – хендл первого event’а,
3. end – хендл второго event’а.

Возвращает:

1. cudaSuccess – в случае успеха
2. cudaErrorInvalidValue – неверное значение
3. cudaErrorInitializationError – ошибка инициализации
4. cudaErrorPriorLaunchFailure – ошибка при предыдущем асинхронном запуске функции
5. cudaErrorInvalidResourceHandle – неверный хендл event’а

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

__host__ void printMatrixToFile(char* fileName, float* matrix, int width, int height)
{
  FILE* file = fopen(fileName, "wt");
  for (int y = 0; y < height; y++)
  {
    for (int x = 0; x < width; x++)
    {
      fprintf(file, "%.0f\t", matrix[y * width + x]);
    }
    fprintf(file, "\n");
  }
  fclose(file);
}

* This source code was highlighted with Source Code Highlighter.

Если матрицы очень большие, то вывод данных в файлы может сильно замедлить выполнение программы.

Заключение

В процессе тестирования я использовал матрицы размерностью 2048 * 1536= 3145728 элементов и 20 итераций в нагрузочных циклах. После результатов замеров у меня получились следующие результаты (рис. 2).


Рис. 2. Время расчетов. (меньше –лучше).

Как видно, GPU версия с разделяемой памятью выполняется почти в 20 раз быстрее, чем версия на CPU. Так же стоит отметить, что при использовании разделяемой памяти расчет выполняется примерно в 4 раза быстрее, чем без неё.
В своем примере я не учитываю время копирования данных с хоста на девайс и обратно, но в реальных приложениях их так же необходимо брать в расчет. Количество перемещений данных между CPU и GPU по-возможности необходимо свести к минимуму.

P.S. Надеюсь, вам понравился прирост производительности, который можно получить с помощью GPU.