Нижегородский государственный университет им. Н.И. Лобачевского

Национальный исследовательский университет

Учебно-научный и инновационный комплекс

«Новые многофункциональные материалы и нанотехнологии»

Исследовательской школа «Наноматериалы и нанотехнологии»

Основная профессиональная образовательная программа аспирантуры

01.04.10 Физика полупроводников

Название дисциплины Параллельное программирование в физических исследованиях

Денисенко М.В., Сатанин А.М.

ПРИМЕНЕНИЕ ГЕТЕРОГЕННЫХ ВЫЧИСЛИТЕЛЬНЫХ СИСТЕМ И ТЕХНОЛОГИИ CUDA ДЛЯ МОДЕЛИРОВАНИЯ ФИЗИЧЕСКИХ ПРОЦЕССОВ

Электронное учебно-методическое пособие

Мероприятие 3.1: Развитие системы поддержки ведущих научно-педагогических коллективов, молодых ученых, преподавателей и специалистов

Нижний Новгород

2012

ПРИМЕНЕНИЕ ГЕТЕРОГЕННЫХ ВЫЧИСЛИТЕЛЬНЫХ СИСТЕМ И ТЕХНОЛОГИИ CUDA ДЛЯ МОДЕЛИРОВАНИЯ ФИЗИЧЕСКИХ ПРОЦЕССОВ

Денисенко М.В., Сатанин А.М. Электронное учебно-методическое пособие. – Нижний Новгород: Нижегородский госуниверситет, 2012. – 53 с.

Аннотация. В учебно-методическом пособии рассматриваются возможности применения технологии CUDA в научных исследованиях. Описываются принципы распараллеливания, и демонстрируется пример работающей программы – расчет диссипативной динамики квантового элемента памяти (кубита) методом квантовых траекторий (квантовым методом Монте-Карло). На базе примера показана эффективность работы графических ускорителей, масштабируемость и описывается возможность взаимодействия нескольких графических ускорителей с применением технологии MPI. Приведена схема вычислительного кластера НИФТИ ННГУ, дано описание и порядок выполнения лабораторных работ по изучению принципов гетерогенных вычислительных систем и реализации численных физических экспериментов. Пособие содержит задачи для самостоятельного решения и вопросы для анализировано результатов.

Электронное учебно-методическое пособие предназначено для аспирантов ННГУ, обучающихся по основной профессиональной образовательной программе аспирантуры

01.04.07 Физика конденсированного состояния, студентов старших курсов, изучающих курс параллельное программирование в физических исследованиях, а также для использования в УНИК «Новые многофункциональные материалы и нанотехнологии».

СОДЕРЖАНИЕ

ВВЕДЕНИЕ 4

1 Параллельное программирование на GPU 5

1.1 Общие принципы работы GPU-ускорителей 7

1.2 Типы памяти 10

1.3 Расширения языка C для работы с CUDA 12

1.4 Установка и настройка среды CUDA под Windows 15

1.5 Компиляция CUDA приложений 17

2 Вычислительная установка (схема кластера НИФТИ ННГУ) 19

3 Пример программы с применением технологии CUDA: расчет диссипативной динамики кубита квантовым методом Монте-Карло 21

3.1. Математическая модель 22

3.2 Схема распараллеливания 26

3.3 Результаты расчета и полученное ускорение 32

4 Указания по выполнению лабораторной работы 36

5 Индивидуальные дополнительные задания 37

ЗАКЛЮЧЕНИЕ 38

Литература 39

ПРИЛОЖЕНИЕ А. Пример программы на CUDA: расчет диссипативной динамики кубита квантовым методом Монте-Карло 41

ПРИЛОЖЕНИЕ Б. Пример программы CUDA+MPI: расчет интерференционных картин вероятности переходов кубита квантовым методом Монте-Карло 45

ВВЕДЕНИЕ

Для ускорения расчетов в последние годы активно применяются возможности современных графических карт, реализующие массово-параллельные вычисления общего назначения на мощных графических процессорах. Примером данной технологии является программно-аппаратная архитектура CUDA (от англ. Compute Unified Device Architecture), разработанная компанией Nvidia [1]. CUDA реализует аппаратный параллелизм, базируясь на принципах вычислений SIMD (от англ. Single Instruction Multiple Data), т.е. позволяет применять одни и те же команды параллельно к множеству данных.

Методическая работа базируется на оборудовании, закупленном в рамках программы развития ННГУ как национального исследовательского университета – вычислительного кластера НИФТИ ННГУ. В рамках данной работы описываются возможности применения технологии CUDA в научных исследованиях, описываются принципы распараллеливания и приводится пример работающей программы – расчет диссипативной динамики квантового элемента памяти (кубита) квантовым методом Монте-Карло. На базе примера демонстрируется эффективность работы графических ускорителей, масштабируемость и описывается возможность взаимодействия нескольких графических ускорителей с применением технологии MPI.

Развитая в работе техника расчета диссипативной динамики кубитов может быть естественным образом распространена на более сложные системы. Например, используемый квантовый метод Монте-Карло применим для моделирования работы многих других квантовых объектов (например, квантовых ям, квантовых точек, ионов в ловушках и т.д.) в реальных условиях взаимодействия с внешним окружением, а также для изучения работы генератора одиночных фотонов, полупроводниковых транзисторов и генераторов на основе квантовых точек и др.

1. Параллельное программирование на GPU

В настоящее время высокопроизводительные вычисления являются необходимым составляющим звеном во многих видах человеческой деятельности: научных исследованиях, промышленном производстве, образовании [2]. Уровень развития вычислительной техники и методов математического моделирования позволяет выполнять симуляцию и анализ природных и технологических процессов, что обеспечивает эффективное решение научных и производственных задач, предоставляя промышленному производству конкурентное преимущество, а исследователям – лидирующие позиции в современной науке. Спектр применения суперкомпьютерных технологий чрезвычайно широк: авиа-, судо- и машиностроение, строительство, нефтегазодобывающая промышленность, моделирование физических, химических и климатических процессов, медицина, биотехнологии и многое другое.

Современные суперкомпьютеры обеспечивают выполнение до 1016 операций с плавающей точкой в секунду (10 Петафлоп) при запусках тестов Linpack (но существуют другие подходы к оценке их производительности) [3]. Достаточно ли подобной вычислительной мощности для решения современных задач? Конечно, любая имеющаяся производительность будет полностью потреблена пользователями, но даже грубая оценка необходимых объемов вычислений показывает, что для решения многих задач требования к производительности отличаются на порядки от возможностей существующих ресурсов, например:

- моделирование номинальных и переходных режимов работы ядерного реактора: для расчета процессов в реакторной установке требуется порядка 100 часов на машине производительностью 1 Эфлопс [4];

- проектирование и разработка изделий в авиастроении: при использовании метода прямого численного моделирования требуется использовать сетку размером 1016 элементов и выполнять не менее 106 шагов моделирования по времени [4];

- квантово-химические расчеты малых систем для одного состояния методом связанных кластеров (Coupled Cluster Singles and Doubles, CCSD) для системы из нескольких десятков (~50) атомов требуется порядка 1014 операций; при решении оптимизационных задач требуется выполнение расчетов для десятков тысяч состояний (при большом количестве измерений – миллионов).

Однако перенос расчетов на гетерогенные вычислительные устройства, содержащие GPU узлы, позволяет существенно ускорить расчет. Графический процессор изначально создавался как многоядерная структура, в которой количество ядер может достигать сотен, а современные CPU содержат несколько ядер (на большинстве современных систем от 2 до 6, по состоянию на 2012 г.). Чтобы понять, какие преимущества приносит перенос расчётов на GPU, приведём усреднённые цифры, полученные исследователями по всему миру. В среднем, при переносе вычислений на GPU, во многих задачах достигается ускорение в 5-30 раз по сравнению с быстрыми универсальными процессорами [5]. Самые большие цифры (порядка 100-кратного ускорения и даже более!) достигаются на коде, который не очень хорошо подходит для расчётов при помощи блоков SSE (Streaming SIMD Extensions, потоковое SIMD-расширение процессора), но вполне удобен для GPU. Рассмотрим здесь некоторые примеры ускорений при использовании синтетического кода на GPU против SSE-векторизованного кода на CPU (по данным NVIDIA):

- флуоресцентная микроскопия: в 12 раз по сравнению с расчетом на CPU (12x);

- молекулярная динамика (non-bonded force calc): 8-16x;

- электростатика (прямое и многоуровневое суммирование Кулона): 40-120x и 7x.

Рис. 1 Примеры использования GPU-ускорителей в научных исследованиях и полученные ускорения по сравнению с CPU. Данные взяты с официального сайта NVidia 2008 г.

1. Общие принципы работы GPU-ускорителей

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

Графический процессор представляется в виде набора независимых мультипроцессоров (multiprocessors). Каждый мультипроцессор состоит из нескольких CUDA-ядер (CUDA cores), нескольких модулей для вычисления математических функций (SFU), конвейера, а также разделяемой памяти (shared memory) и, кэша (для определенных видов памяти).

Графический процессор основан на так называемой архитектуре SIMT (Single Instruction, Multiple Thread). Технология CUDA позволяет определять специальные функции – ядра (kernels), которые выполняются параллельно на CPU в виде множества различных потоков (threads). Таким образом, ядро является аналогом потоковой функции. Каждый поток исполняется на одном CUDA-ядре, используя собственный стек инструкций и локальную память.

Отдельные потоки группируются в блоки потоков (thread block) одинакового размера, при этом каждый блок потоков выполняется на отдельном мультипроцессоре. Количество потоков в блоке ограничено (максимальные значения для конкретных устройств могут быть найдены в [7] или получены во время выполнения при помощи функций CUDA API). Потоки внутри блока потоков могут эффективно взаимодействовать между собой с помощью общих данных в разделяемой памяти и синхронизации. Кроме того, потоки могут взаимодействовать при помощи глобальной памяти и атомарных операций.

На аппаратном уровне потоки блока группируются в так называемые варпы (warps) по 32 элемента (на всех текущих устройствах), внутри которых все потоки параллельно выполняют одинаковые инструкции (по принципу SIMD – Single Instruction, Multiple Data). Важным моментом является то, что потоки фактически выполняют одну и ту же команды, но каждая со своими данными. Поэтому если внутри варпа происходит ветвление (например в результате выполнения оператора if), то все нити варпа выполняют все возникающие при этом ветви. По этой причине операции ветвления могут негативно сказываться на производительности – различные пути не могут выполняться параллельно (в тоже время потоки одного варпа, выполняющие один путь, работают параллельно).

В свою очередь, блоки потоков объединяются в решетки блоков потоков (grid of thread blocks). Следует отметить, что взаимодействие потоков из разных блоков во время работы ядра затруднено: отсутствуют явные инструкции синхронизации, взаимодействие возможно через глобальную память и использованием атомарных функций (другим вариантом является разбиение ядра на несколько ядер без внутреннего взаимодействия между потоками разных блоков).

Рис. 2 Иерархия потоков CUDA

Каждый поток внутри блока потоков имеет свои координаты (одно-, двух- или трехмерные), которые доступны через встроенную переменную threadIdx. В свою очередь, координаты блока потоков (одно-, двух- или трехмерные) внутри решетки определяются встроенной переменной blockIdx. Пример иерархии потоков приведен на рис. 2. Данные встроенные переменные являются структурами с полями.x,.y,.z.

2. Типы памяти

Кроме иерархии потоков существует также несколько различных типов памяти. Быстродействие приложения очень сильно зависит от скорости работы с памятью. Именно поэтому в традиционных CPU большую часть кристалла занимают различные кэши, предназначенные для ускорения работы с памятью (в то время как для GPU основную часть кристалла занимают ALU).

В CUDA для GPU существует несколько различных типов памяти, доступных нитям, сильно различающихся между собой (см. табл. 1) [1].

Таблица 1. Типы памяти в CUDA

Тип памяти	Доступ	Уровень выделения	Скорость работы
регистры (registers)	R/W	per-thread	высокая (on chip)
local	R/W	per-thread	низкая (DRAM)
shared	R/W	per-block	высокая (on-chip)
global	R/W	per-grid	низкая(DRAM)
constant	R/O	per-grid	высокая(on chip L1 cache)
texture	R/O	per-grid	высокая(on chip L1 cache)

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

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

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

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

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

Естественно, что глобальная, локальная, текстурная и память констант – это физически одна и та же память, известная как локальная видеопамять видеокарты. Их отличия в различных алгоритмах кэширования и моделях доступа. При этом центральный процессор (CPU) имеет R/W доступ только к глобальной, константной и текстурной памяти (находящейся в DRAM GPU) и только через функции копирования памяти между CPU и GPU (предоставляемые CUDA API).

3. Расширения языка C для работы с CUDA

Программы для CUDA (соответствующие файлы обычно имеют расширение .cu) пишутся на «расширенном» С и компилируются при помощи команды nvcc.

Вводимые в CUDA расширения языка С состоят из

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

Таблица 2. Спецификаторы функций в CUDA

Спецификатор	Выполняется на	Может вызываться из
__device__	device	device
__global__	device	host
__host__	host	host

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

Спецификатор __global__ обозначает ядро и соответствующая функция должна возвращать значение типа void.

__global__ void myKernel ( float * a, float * b, float * c )

{

int index = threadIdx.x;

c [i] = a [i] * b [i];

}

На функции, выполняемые на GPU (__device__ и __global__) накладываются следующие ограничения:

• нельзя брать их адрес (за исключением __global__ функций);
• не поддерживается рекурсия;
• не поддерживаются static-переменные внутри функции;
• не поддерживается переменное число входных аргументов.

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

• эти спецификаторы не могут быть применены к полям структуры (struct или union);
• соответствующие переменные могут использоваться только в пределах одного файла, их нельзя объявлять как extern;
• запись в переменные типа __constant__ может осуществляться только CPU при помощи специальных функций;
• __shared__ переменные не могут инициализироваться при объявлении.

Добавленные переменные

В язык добавлены следующие специальные переменные

• gridDim - размер grid'а (имеет тип dim3);
• blockDim - размер блока (имеет тип dim3);
• blockIdx - индекс текущего блока в grid'е (имеет тип uint3);
• threadIdx - индекс текущей нити в блоке (имеет тип uint3);
• warpSize - размер warp'а (имеет тип int).

В язык добавляются 1/2/3/4-мерные вектора из базовых типов - char1, char2, char3, char4, uchar1, uchar2, uchar3, uchar4, short1, short2, short3, short4, ushort1, ushort2, ushort3, ushort4, int1,int2, int3, int4, uint1, uint2, uint3, uint4, long1, long2, long3, long4, ulong1, ulong2, ulong3, ulong4, float1, float2, float3, float2, и double2.

Обращение к компонентам вектора идет по именам - x, y, z и w. Для создания значений-векторов заданного типа служит конструкция вида make_<typeName>.

int2 a = make_int2 ( 1, 7 );

float3 u = make_float3 ( 1, 2, 3.4f );

Обратите внимание, что для этих типов (в отличии от шейдерных языков GLSL/Cg/HLSL) не поддерживаются векторные покомпонентные операции, т.е. нельзя просто сложить два вектора при помощи оператора "+" - это необходимо явно делать для каждой компоненты.

Также для задания размерности служит тип dim3, основанный на типе uint3, но обладающий нормальным конструктором, инициализирующим все не заданные компоненты единицами.

Директива вызова ядра

Для запуска ядра на GPU используется следующая конструкция:

kernelName <<<Dg,Db,Ns,S>>> ( args )

Здесь kernelName это имя (адрес) соответствующей __global__ функции, Dg - переменная (или значение) типа dim3, задающая размерность и размер grid'a (в блоках), Db - переменная (или значение) типа dim3, задающая размерность и размер блока (в нитях), Ns - переменная (или значение) типа size_t, задающая дополнительный объем shared-памяти, которая должна быть динамически выделена (к уже статически выделенной shared-памяти), S - переменная (или значение) типа cudaStream_t задает поток (CUDA stream), в котором должен произойти вызов, по умолчанию используется поток 0. Через args обозначены аргументы вызова функции kernelName.

Также в язык С добавлена функция __syncthreads, осуществляющая синхронизацию всех нитей блока. Управление из нее будет возвращено только тогда, когда все нити данного блока вызовут эту функцию. Т.е. когда весь код, идущий перед этим вызовом, уже выполнен (и, значит, на его результаты можно смело рассчитывать). Эта функция очень удобная для организации бесконфликтной работы сshared-памятью.

Также CUDA поддерживает все математические функции из стандартной библиотеки С, однако с точки зрения быстродействия лучше использовать их float-аналоги (а не double) - например sinf. Кроме этого CUDA предоставляет дополнительный набор математических функций (__sinf, __powf и т.д.) обеспечивающие более низкую точность, но заметно более высокое быстродействие чем sinf, powf и т.п.

4. Установка и настройка среды CUDA под Windows

CUDA – это программно-аппаратная вычислительная технология от компании NVidia, поэтому перед началом установки необходимо проверить, поддерживает ли ваша видеокарта технологию CUDA (это можно сделать на официальном сайте: https://developer.nvidia.com/cuda-gpus). Для разработки и запуска приложений необходимо:

• графический процессор, поддерживающий CUDA;
• драйвер для устройства и CUDA Toolkit;
• программное обеспечение NVidia, которое можно бесплатно загрузить с сайта https://developer.nvidia.com/cuda-downloads
• среда разработки программ (Microsoft Visual Studio, Netbeans, etc.).

В данном методическом пособии будет рассмотрена интеграция и установка CUDA (версия 4.0) с помощью Microsoft Visual Studio. Отметим, что у Вас должна иметься уже установленная программа Microsoft Visual Studio 2010.

Этапы установки:

1. Необходимо загрузить и установить драйвер для работы CUDA, который интегрирован в последние версии NVidia ForceWare (см. http://www.nvidia.com/drivers). Нужно обратить внимание, что версия драйвера должна соответствовать требованиям, указанным в примечаниях к релизу (Release Notes) CUDA Toolkit.
2. Скачать драйверы для разработчиков (CUDA Developer Drivers) находятся в разделе CUDA Toolkit на сайте NVidia: https://developer.nvidia.com/cuda-toolkit. Этот программный пакет содержит инструменты, библиотеки, заголовочные файлы для компиляции программ с помощью Microsoft Visual Studio. GPU Computing SDK содержит в себе демонстрационные примеры, которые уже предварительно сконфигурированы для удобной работы в среде Microsoft Visual Studio, и является исключительно полезным при изучении технологии, но не требуется для создания и запуска приложений.
3. Запустите установочный пакет (CUDA Toolkit) и поэтапно следуйте процессу установки. Инструментарий CUDA по умолчанию устанавливается в каталог C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v4.0. В процессе установки автоматически создаются соответствующие переменные окружения (в том числе CUDA_LIB_PATH, CUDA_INC_PATH), для настройки можно использовать их.
4. Для компиляции и сборки проектов в настройках среды (Microsoft Visual Studio) необходимо выбрать вкладку

Tools | Options | Projects and Solutions | VC++ Directories.

Во вкладку Include files добавить следующие каталоги:

C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v4.0\include

Во вкладку Library files добавить следующие каталоги:

C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v4.0\lib\Win32 (или C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v4.0\lib\x64 для 64-битной сборки)

5. Запустите установочный пакет GPU Computing SDK и поэтапно следуйте процессу установки. Данный программный комплект содержит не только исходные коды программ для множества примеров и задач, но и шаблоны для Microsoft Visual Studio.
6. Для компиляции с использованием CUDA необходимо использовать соответствующие правила сборки (они автоматически добавляются в общий список при установке CUDA Toolkit): пункт меню

Project | Custom Build Rules…, далее выбрать CUDA Runtime API Build Rule. При этом в настройках проекта будет добавлена вкладка CUDA Runtime API с настройками, специфичными для CUDA-части. Компилятор nvcc может быть вызван и отдельно, через командную строку.

5. Компиляция CUDA приложений

Для проверки правильной работоспособности CUDA следует произвести компиляцию и запуск некоторых из включённых в SDK примеров. GPU Computing. Рекомендуется запуск программы bandwidthTest, расположенной в C:\ Documents and Settings\All Users\Application Data\NVIDIA Corporation\NVIDIA GPU Computing SDK #.#\C\bin\win32\ Release в Windows XP и в %ProgramData%\NVIDIA Corporation\ NVIDIA GPU Computing SDK #.#\C\bin\win32 в случае с Windows Vista и позже. (На 64-х разрядных версиях Windows путь к директории заканчивается на win64\Release). В случае если установка была произведена корректно, вы увидите на экране терминала примерно следующую картину см. рис. 3.

Рис. 3 Допустимые результаты тестовой программы bandwidthTest 

Название устройства (вторая строка) и пропускная способность отличаются от системы к системе. Здесь важно содержимое пары строк: второй, которая отображает имя модели устройства, и предпоследней, которая сообщает, что все требуемые тесты пройдены.

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

Замечания.

1) Убедитесь, что при установке CUDA Toolkit вы не забыли поставить галочку напротив CUDA Toolkit Visual Studio Integration (должна быть установлена по умолчанию).

2) CUDA Toolkit версии 3.1 и ниже устанавливаются в C:\CUDA по умолчанию, предварительно требуя деинсталяции более ранних версий CUDA Toolkit. С версии 3.2 появилась возможность одновременной установки сразу нескольких версий пакета.

2. Вычислительная установка (схема кластера НИФТИ ННГУ)

Методическая работа базируется на оборудовании, закупленном в рамках программы развития ННГУ как национального исследовательского университета – вычислительного кластера НИФТИ ННГУ. Кластер состоит из девяти рабочих станций (см. схему на рис. 4), восемь из которых содержат графические ускорители nVidia® Tesla.

Рис. 4. Схема кластера лаборатории «Теория наноструктур» НИФТИ ННГУ им. Лобачевского.

Краткое техническое описание оборудования и сети:

1. Сервер Flagman QD820 (8 процессоров AMD® Opteron™ SixCore, 16 х DIMM 4096Mb DDR-II, 4 х HDD 300Gb SerialATA 10000rpm);
2. процессоров AMD® Opteron™ SixCore, 16 х DIMM 4096Mb DDR-II, 4 х HDD 300Gb SerialATA 10000rpm);
3. 2 станции Flagman WX240T.2 (2 процессора Intel® Xeon® X5550, 12 х DIMM 2048Mb DDR-III, 4 вычислителя nVidia® Tesla® C1060 4096Mb DDR-III);
4. 6 станций Flagman WP120N.2 (процессор Intel® Core™ i7 I7-950, 6 х 2048Mb DDR-III, 2 вычислителя nVidia® Tesla® C1060 4096Mb DDR-III).

Сетевое оборудование: Вычислительная сеть: 20 Gbit/s InfiniBand (коммутратор Mellanox MTS3600); Служебная сеть: 1 Gbit/s (коммутатор 3Com, 24 ports).

Для удобства доступа к кластеру организован удаленный доступ с помощью стандартной программы ОС Windows – Удаленный рабочий стол, из аудиторий физического факультета ННГУ им. Н.И. Лобачевского №419 (лаборатория «Теория наноструктур») и на №537 (терминал класс кафедры теоретической физики, физического факультета).

Для выполнения лабораторной работы в оборудованном терминал классе (аудитория №537), необходимо произвести следующие действия:

Пуск Все программы Стандартные Подключение к удаленному рабочему столу

или

C:\Windows\system32\mstsc.exe

Для подключения к выбранному компьютеру укажите IP-адрес и порт:

Для подключения к Lab25-6 – 85.143.6.98:3389

Для подключения к Lab25-7 – 85.143.6.98:3390

Для подключения к Lab25-8 – 85.143.6.98:3391

Введите имя учетной записи и пароль:

Имя учетной записи: cluster\Theorlab

Пароль: Theorlab

Для завершения работы просто закройте окно Удаленного рабочего стола.

3. Пример программы с применением технологии CUDA: расчет диссипативной динамики кубита квантовым методом Монте-Карло

В рамках данной работы описываются возможности применения технологии CUDA в научных исследованиях, описываются принципы распараллеливания и демонстрируется пример работающей программы – расчет диссипативной динамики квантового элемента памяти (кубита) квантовым методом Монте-Карло (МК).

Обычно для численного расчета динамики квантовых систем с учетом диссипации используют уравнение для оператора плотности (Master equation), содержащее n2 переменных, где n – размерность гамильтониана. Другие методы (Гейзенберга-Ланжевена, функции распределения) используются, в основном, для аналитических расчетов. За счет использования квантового метода МК, удается свести задачу к решению уравнения для волновой функции системы, содержащее n переменных, однако полученное решение необходимо усреднить по m реализациям, обычно m порядка 1000-10000 для получения нужной точности. Реализации статистически независимы, поэтому есть возможность запускать каждую реализацию в отдельном потоке (на отдельном процессоре), собирая затем данные и усредняя. Отсутствует необходимость в обмене данными между потоками в процессе работы программы, что позволяет использовать многопроцессорные системы, особенно системы на базе GPU-ускорителей практически со 100% эффективностью. Таким образом, ускорение по сравнению со стандартными методами расчета может достигать n раз.

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

3.1. Математическая модель

В последнее десятилетие в связи с развитием нанотехнологии и квантовой информации возник интерес к изучению поведения во внешних полях «искусственных атомов»: полупроводниковых квантовых точек, квантовых ям, ионов в магнитных ловушках, макроскопических сверхпроводящих контуров и т.п. Управление и контроль квантовых состояний таких систем представляется важным, поскольку активно обсуждается возможность использования «искусственных атомов» в качестве логических элементов (кубитов) для квантового компьютера, а также для устройств кодирования и передачи информации [8, 9].

Роль кубита может играть любая двухуровневая система. При этом связь с внешней средой и экспериментальной установкой приводит к обратному шумовому воздействию на систему. Неконтролируемое воздействие шума нарушает когерентность и затрудняет реализацию квантовых вычислений [8], таким образом, необходимым условием является устойчивость таких систем к шумам и возможность производить «неразрушающие измерения» [10].

Одна из перспективных реализаций кубита представляет собой сверхпроводящую петлю, прерванную тремя джозефсоновскими переходами [25] (рис. 5). Они обладают малой диссипацией, хорошей устойчивостью к шумам и относительно простым способом управления состояниями, что важно в квантовой информатике при передаче, хранении и обработке информации. Логическими состояниями кубита являются направления незатухающего тока в сверхпроводящей петле. Для детектирования состояний в настоящее время предложено несколько экспериментальных методик. Первоначально группой Дж. Моои (J.E. Mooij) был предложен метод с использованием СКВИДа (см. рис. 5 (а)), обладающего двумя дополнительными джозефсоновскими переходами [11-13]. В дальнейшем для изучения квантовых когерентных эффектов и мониторинга состояний был предложен Е. Ильичёвым новый радиочастотный (РЧ) метод [14,15], в котором кубит связывается с высокочастотным резонансным колебательным контуром, т.е. помещается в центр сверхпроводящей ниобиевой планарной катушки (см. рис. 5 (б)).

Рис. 5. Примеры реализации потокового кубита (а) индуктивно связанного с СКВИДом [11] и (б) связанного со сверхпроводящей управляемой катушкой [15]. На вставке справа представлен кубит с тремя слабыми связями (желтые кривые), расположенный в центре катушки.

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

, (1)

где - управляющий параметр кубита, имеющий смысл магнитного поля, - туннельное расщепление уровней, а и - матрицы Паули. Для изучения когерентной динамики в сильных полях в последнее время применяется новая экспериментальная методика – амплитудная спектроскопия [16-18], в основе которой лежит метод получения информации с помощью функции отклика по амплитудам постоянного и переменного поля сигнала при фиксированной частоте. В рамках данной методики считается, что на кубит подается внешнее управляющее поле , где – амплитуда постоянного магнитного поля, – амплитуда переменного высокочастотного поля частоты , где T – период внешнего магнитного поля.

При рассмотрении диссипативной динамики кубита существенное влияние оказывает взаимодействие с резервуаром. Заметим, что эксперименты над кубитами проводятся при низких температурах (~мК), поэтому вклад шума квазичастичных состояний (фермионных возбуждений) мал, следовательно, процессы релаксации могут описываться только фононным резервуаром с большим числом степеней свободы. Управление состояниями потокового кубита осуществляется внешними полями, которые и являются главным источником шума. Экспериментально показано [19], что время поперечной релаксации велико по сравнению с временем продольной, поэтому уравнение для матрицы плотности имеет вид:

, (2)

где - скорость поперечной релаксации.

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

. (3)

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

, (4)

где – нормировочная константа. Вероятность квантовых скачков определяется следующим выражением:

. (4)

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

Учитывая, что процесс релаксации является случайным, каждая траектория уникальна. Полученное решение необходимо усреднить по M реализациям, обычно M порядка 1000-10000, для достижения точности квантового МК, определяемой . Реализации статистически независимы, поэтому есть возможность запускать каждую реализацию независимо в отдельном потоке (на отдельном процессоре), собирая затем данные и усредняя. Отсутствует необходимость в обмене данными между потоками в процессе работы программы, что позволяет использовать многопроцессорные системы, особенно системы на базе GPU-ускорителей практически с максимальной эффективностью, см. Приложение А.

Рис. 6 Алгоритм моделирования диссипативной динамики потокового кубита, основанного на квантовом методе Монте-Карло.

3.2 Схема распараллеливания

Как сказано выше, на основе метода МК рассчитывается динамика единичной реализации на отдельных GPU с применением технологии CUDA [1,5]. Поскольку реализации статистически независимы, отсутствует необходимость в обмене данными между блоками и потоками внутри отдельных блоков, создаваемых на GPU, на каждом процессоре выполняется одна и та же подпрограмма, что позволяет обеспечить параллелизм на уровне данных (SIMT).

Известно, что в реальных экспериментах получают данные о поведении квантового прибора в зависимости от начальных условий и параметров задачи [17-19]. Например, для потокового кубита – это интерференционные зависимости вероятностей населенности возбужденного уровня кубита от амплитуд постоянного и переменного А магнитных полей, которые в качестве примера рассчитываются в рамках данной работы. Так как значения вероятностей при каждых значениях параметров и А независимы, то возможно повышение производительности работы программы за счет использования Кластерной системы, узлы которой содержат GPU, на каждом из которых значение вероятности населенности рассчитываться независимо. Для обеспечения работы нескольких GPU была использована технология программирования для систем с распределенной памятью MPI [21, 22], см. Приложение Б.

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

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

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

Рис. 7 Алгоритм моделирования вероятности населенности возбужденного уровня кубита, где блоки «Динамика системы (МК)» реализуются согласно алгоритму, приведенному на рис.6. Пунктирной линей представлена схема алгоритма для вычисления усредненной единичной квантовой траектории, выполняемой на каждом из GPU-ускорителях.

Приведем более подробное описание реализованного алгоритма взаимодействия нескольких GPU-ускорителей. С программным кодом более подробно можно ознакомиться в Приложении Б. Рассылка данных осуществляется от главного потока Master, запускаемом на CPU, всем остальным доступным видеокартам, на которых уже запускается копия программа, представленная в Приложении А.

Таблица 3. Схема работы MPI приложения для взаимодействия нескольких GPU устройств

Master (myid = 0)	Slove (myid = 1 … N, где N – число доступных GPU)
Описание функций
void master () { float data[2];//массив нач. данных float res[3]; //массив результ. данных, первые два элемента это data[2] init_io(); init_net(); while (load_data ((float *) data)) { if (! is_free ()) { receive_result (res); save_result (res); } send_data ((float *) data); } while (! are_all_free ()) { receive_result (res); save_result (res); } close_io(); send_end (); }	void slave () { float data[2]; float res[3]; int j; j = 0; int myid; MPI_Comm_rank (MPI_COMM_WORLD, &myid); int device = myid % 2; cudaSetDevice(device); cudaDeviceProp deviceProp; cudaGetDeviceProperties(&deviceProp, device); while (receive_data ((float *) data)) { calc (data, res,deviceProp.maxThreadsDim[0], deviceProp.multiProcessorCount); send_result (res); j = j + 1; } }
Описание рассылки
init_io() – инициализация начальной точки, с которой начинается проход по сетке. init_net() – устанавливает число выполняющих расчёт процессов n, в начале программы n = 0, так как ещё не один из slove-процессов не считает. load_data - отвечает за проход по сетке (перемещение от узла к узлу) if (! is_free ()) - проверка на то, есть ли посчитавшие процессы, готовые отправить уже рассчитанные данные для записи их в файл. Но так как ещё все slove – процессы свободны (не получили начальные данные), то условие не выполняется и переходим к процедуре send_data – рассылка начальных данных. Принимаем запрос с тегом TAG_WMSG от всех slove-процессов (MPI_ANY_SOURCE) Далее master определяет ранг каждого процесса, который отправил запрос, и отправляет ему сообщение с текстом msg = W_OK, что означает «Готов отправить данные» Далее начинается отправка данных (массива data). Отметим, что сообщения с начальными данными имеют тег TAG_DATA. После каждой порции отправки число «занятых» (выполняющих расчёт) процессов: n = n + 1.	Функция receive_data выполняет запрос к master: «Готов ли он дать начальные данные?». Для этого slove-процесс отправляет сообщение (msg) с тегом TAG_WMSG, именно по тегу сообщения и различается данный запрос главным процессом. Принимает сообщение msg = W_OK и «встает» на приём данных.
На все процессы (slove) – устройства GPU – разосланы начальные задания. Число «занятых» процессов n = N.
Master находится в ожидании, функция is_free() позволила данному процессу перейти к receive_result, то есть «встать» на прием данных от всех slove-процессов (MPI_ANY_SOURCE), а также после приема число «занятых» процессов: n = n - 1. Затем полученные данные сохраняются в файл, для этого вызывается специальная процедура save_result.	Выполняется передача данных в модуль calc.cpp и производится расчёт на GPU c номером (myid % 2) на локальной машине при вызове фнкции void montecarlo( singlecomplex *a_dev, singlecomplex *b_dev, float *rand_dev,int N, int T, float Am, float eps,int maxThreadsPerBlock, int NumMultiprocessors) из файла gpu.cu Затем результат передаётся master. Для этого вызывается процедура send_result. Заметим, что сообщения с результатом имеют тег TAG_RESULT.
Когда последняя партия начальных данных будет разослана, то master выйдет из своего цикла while, а расчёт данных может ещё продолжатьсяю и необходимо принять последнюю порцию посчитанных результатов. Рассмотрим эту ситуацию подробнее.
are_all_free () – позволяет master «стоять» на приёме сообщений до тех пор пока n = 0 (n – число считающих в данный момент процессов) Выполняются аналогичные процедуры: receive_result и save_result. close_io() – закрытие файла для записи, когда n = 0. send_end () – отправка всем slove-процессам сообщение: «Не буду отправлять данные». То есть msg = W_NO c тегом TAG_WMSG.	Получив сообщение с msg = W_NO и тегом TAG_WMSG, функция receive_data останавливает работу программы.

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

Заметим, что точность квантового метода МК зависит от числа реализаций ~ , где M – число реализаций. В силу наличия стохастических процессов двойная точность для расчётов не требуется, а использование операций с одинарной точностью на GPU, существенно ускоряет работу приложения, т.к. вычисления с двойной точностью на несколько порядков медленнее [23]. Для обеспечения сходимости результатов моделирования диссипативной динамики необходимо наличие большого числа реализаций (M > 103) и генерирования последовательностей случайных чисел непосредственно для каждой точки временной траектории. Генерация необходимого массива случайных чисел производилась на GPU с помощью возможностей CUDA библиотеки – CURAND. Использование данной библиотеки существенно позволило увеличить скорость выполнения расчётов, так как скорость генерации случайных чисел на GPU в 50 раз больше, чем на СPU [24].

3.3 Результаты расчета и полученное ускорение

Для расчёта единичных усредненных траекторий было выбрана длительность импульса = 5 T, где Т – период внешнего поля, а шаг по времени выбирался, чтобы удовлетворить как условиям метода МК [20], так и обеспечить точность метода Рунге–Кутта четвертого порядка применявшегося при численном моделировании динамики между скачками. При данных параметрах наличие статистических ошибок не должно превышать 5%, что приемлемо для извлечения необходимой физической информации из такого рода систем.

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

Рис. 8 Населенность верхнего уровня кубита (усредненная квантовая траектория).

В Табл. 3. представлены характерные времена, требуемые для вычисления c одинарной точностью населенности верхнего уровня кубита на CPU (Intel Core i7 960 использовалось 1 ядро) и на GPU (Nvidia Tesla С1060). Версия для GPU создана с использованием технологии CUDA. Видно, что графические ускорители позволяют существенно сократить время расчета, и поэтому их использование является целесообразным при решении данной задачи. Как видно из табл. 4 ускорение относительно CPU растет в зависимости от числа реализаций и достигает примерно 170 раз. Рост ускорения от реализаций обусловлен тем, что время расчета одной траектории мало, а при «перегрузке» графического ускорителя наступает эффективная загруженность и достигается оптимальная работа.

Таблица 4. Время расчета единичной усредненной по реализациям квантовой траектории на CPU и GPU

Число реализаций	Время расчета на CPU, с.	Время расчета на GPU, с.
980	17,082	0,609
1980	35,303	0,656
3840	68,166	0,686
7680	136,454	0,874
15360	272,705	1,654
30720	545,458	3,260

Для расчета интерференционной картины (см. рис.9), наблюдаемой в недавних экспериментах [17-19], была выбрана 2D-сетка (А, 0) 100х100 точек, таким образом, было необходимо рассчитывать 10000 усредненных квантовых траекторий (каждая из которых была усреднена по 15360 реализациям) и, соответственно, вероятность населенности возбужденного уровня в момент окончания длительности импульса.

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

На рис. 10 и в табл. 5 представлены данные о приросте производительности вычислений в зависимости от количества используемых графических ускорителей, данная зависимость аппроксимируется линейной, что говорит о хорошем масштабировании задачи при переходе к множеству GPU, эффективность при использовании 20 GPU составляет 97%.

Рис. 10 Ускорение работы алгоритма в зависимости от числа использованных GPU по сравнению с одним GPU.

Таблица 5. Время расчета и ускорение алгоритма при одновременном расчёте 10000 единичных усредненных траекторий на GPU c применением технологии MPI

Расчёт интерференционных картин (10000 траекторий) на GPU
Тип устройства	Число устройств	Время, с	Ускорение
GPU (Tesla C1060)	1	6020	~ в 66 раза отн. CPU
GPU (Tesla C1060) + MPI	5	1228	~ в 323 раза отн. CPU ~ в 5 раза отн. GPU
GPU (Tesla C1060) + MPI	20	311	~ в 1273 раза отн. CPU ~ в 19,4 раза отн. GPU