CC BY
69
Математическое моделирование. Оптимальное управление Вестник Нижегородского университета им. Н.И. Лобачевского, 2012, № 5 (2), с. 42-48
УДК 004.942, 004.272.26, 530.145
АЛГОРИТМ РАСЧЕТА ДИССИПАТИВНОЙ ДИНАМИКИ КУБИТОВ НА СУПЕРКОМПЬЮТЕРНОМ КОМПЛЕКСЕ С ИСПОЛЬЗОВАНИЕМ GPU-УСКОРИТЕЛЕЙ
© 2012 г. А.В. Волков, М.В. Денисенко, А.М. Сатанин
Нижегородский госуниверситет им. Н.И. Лобачевского
Поступила в редакцию 10.09.2012
Квантовым методом Монте-Карло изучена диссипативная динамика кубита в сильном электромагнитном поле. Алгоритм моделирования реализован на гетерогенной суперкомпьютерной системе с применением технологий СИЭА и МРІ. Предложен метод распараллеливания численной схемы, и продемонстрировано ускорение ~170 раз по сравнению с последовательной версией на СРИ. Разработанный программный комплекс позволяет рассчитывать и предсказывать результаты дорогостоящих реальных экспериментов над кубитами.
Ключевые слова: квантовый метод Монте-Карло, GPU, технология СИЭА, параллельные вычисления, кубит.
Введение
В последнее время благодаря успехам нанотехнологий наблюдается тенденция к минимизации размеров затворов транзисторов, что, несомненно, ведет к предельной точке, когда будут достигнуты атомные размеры, и тем самым к нарушению известного закона Мура [1, 2]. В связи с этим возрастает интерес к изучению поведения во внешних полях «искусственных атомов» - систем, представляющих собой новый тип квантовой памяти. Управление и контроль состояний таких систем представляется перспективным для дальнейшего развития информационно-телекоммуникационных и суперком-пьютерных систем, поскольку активно обсуждается возможность использования «искусственных атомов» в качестве логических элементов (кубитов) для квантового компьютера [3]; также развивается новое направление создания атомных чипов, оптоэлектронных устройств минимальных размеров, позволяющее в устройствах из нескольких атомов организовывать хранение и передачу информации, оперируя системой на временах порядка миллисекунд.
Определяющую роль в работе подобного рода вычислительных устройств играют квантовые релаксационные процессы, поэтому актуальной задачей на сегодняшний день является экспериментальное и теоретическое изучение работы квантовых приборов в реальных условиях взаимодействия с внешней средой. Так, например, для кубитов слабая связь с внешним
окружением, с одной стороны, предохраняет их от быстрой потери когерентности, а с другой -существенно затрудняет связь кубитов между собой. Без такой связи нельзя приготовить запутанные состояния (entanglement states) нескольких кубитов и в полной мере использовать потенциал квантовых вычислений. Современные методы спектроскопии и динамического контроля позволяют также осуществлять «однократные» измерения (single-shot measurement) над кубитами [4-6], то есть мониторинг состояний открытой квантовой системы в реальном времени в условиях связи с внешней средой. Благодаря этому стало возможно исследовать динамику и характеристики отдельных квантовых систем как для каждой реализации, так и в среднем по ансамблю измерений.
В данной работе основное внимание уделяется разработке гетерогенного параллельного алгоритма для моделирования работы потокового кубита [7], находящегося в сильных электромагнитных полях при реальных уровнях шума. Потоковый кубит представляет собой сверхпроводящую петлю с тремя джозефсонов-скими переходами, имеющую макроскопический размер и возможность управления за счет электромагнитных полей.
Для моделирования диссипативной динамики в работе используется квантовый метод Монте-Карло (МК) [8], который идеально подходит для решения задачи на многопроцессорных системах, особенно для систем на базе GPU-ускорителей. Результатом моделирования
является получение единичных квантовых траекторий (аналог «однократных измерений») и интерференционных картин вероятностей переходов в сильных электромагнитных полях, которые позволяют извлечь дополнительную информацию о параметрах кубита и шума. Главной целью является разработка и оптимизация алгоритма с использованием ОРИ-ускорителей и современных технологий параллельной обработки данных (МР1, СИБЛ) для расчета диссипативной динамики квантовых приборов, а также получение физических результатов с помощью данной программы для конкретного прибора - потокового кубита.
Квантово-механическая модель кубита и алгоритм моделирования диссипативной динамики
С точки зрения квантовой механики кубит представляет собой обычную двухуровневую систему, гамильтониан которой имеет вид [7]:
Н = 2(вст * + Дст х X (1)
где е - управляющий параметр кубита, имеющий смысл магнитного поля, Д - туннельное расщепление уровней, ст2 и ох - матрицы Паули. Для изучения когерентной динамики в сильных полях в последнее время применяется новая экспериментальная методика - амплитудная спектроскопия [9-11], в основе которой лежит метод получения информации с помощью функции отклика по амплитудам постоянного и переменного поля сигнала при фиксированной частоте. В рамках данной методики считается, что на кубит подается внешнее управляющее поле е(?) = е0 + Асо8шt , где е0 - амплитуда постоянного магнитного поля, А- амплитуда переменного высокочастотного поля частоты га = = 2тс/Т, Т - период внешнего магнитного поля.
При рассмотрении диссипативной динамики кубита существенное влияние оказывает взаимодействие с резервуаром. Заметим, что эксперименты над кубитами проводятся при низких температурах (~мК), поэтому вклад шума ква-зичастичных состояний (фермионных возбуждений) мал, следовательно, процессы релаксации могут описываться только фононным резервуаром с большим числом степеней свободы. Управление состояниями потокового кубита осуществляется внешними полями, которые и являются главным источником шума. Экспериментально показано [9], что время поперечной релаксации велико по сравнению с временем продольной, поэтому уравнение для матрицы
плотности имеет вид:
р = - 7[Н*, Р] +1 г(^грст+ -рХ (2)
п 2
где Г - скорость поперечной релаксации.
Обычно для численного расчета динамики квантовых систем с учетом диссипации используют уравнение для оператора матрицы плотно-
2
сти (2), содержащее N переменных, где N -размерность гамильтониана. Другие методы (Гейзенберга - Ланжевена, функции распределения) используются, в основном, для аналитических расчетов [12].
Использование квантового метода МК позволяет свести задачу к решению уравнения для волновой функции системы |Т(у)(?)) (/ - номер
траектории), содержащего N переменных и описывающего динамику системы с эффективным гамильтонианом Вигнера - Вайскопфа Не^ =
= Н* - ('ЙГ/4)I (I - единичная матрица)
¥| ч»)(')) = - '-пН «НЧшМ). (3)
В процессе эволюции динамика системы может испытывать квантовые скачки, что соответствует диссипативной динамике [8]. Таким образом, в момент времени t' = t + Дt вектор состояний определяется оператором матрицы Паули <зг, который отвечает за дефазировку системы
|Т( + Дф = Ха г|Т( ^)), (4)
где X - нормировочная константа. Вероятность квантовых скачков определяется следующим выражением:
Г
Р( у)(0 = - Д/|Т( у ^ )|2. (5)
Считая, что время дискретно и меняется с интервалом Дt, численный алгоритм для вычисления одной у-й квантовой траектории, представляющей одну реализацию реального эксперимента - диссипативную динамику одного кубита, можно представить на схеме рис. 1.
Поскольку процесс релаксации является случайным, каждая траектория уникальна. Полученное решение необходимо усреднить по М реализациям, обычно М порядка 1000-10000, для достижения точности квантового метода
МК, определяемой 1/4М. Реализации статистически независимы, поэтому есть возможность запускать каждую реализацию независимо в отдельном потоке (на отдельном процессоре), собирая затем данные и усредняя. Отсутствует необходимость обмена данными между потоками в процессе работы программы, что позволяет использовать многопроцессорные
системы, особенно системы на базе ОРИ-ускорителей, практически с максимальной эффективностью .
Схема распараллеливания
Как сказано выше, на основе метода МК рассчитывается динамика единичной реалии-зации на отдельных ОРИ с применением технологии СИБЛ [13]. Поскольку реализации статистически независимы, отсутствует необходимость обмена данными между блоками и потоками внутри отдельных блоков, создаваемых на ОРИ, на каждом процессоре выполняется одна и та же подпрограмма, что позволяет обеспечить параллелизм на уровне данных (БІМТ).
Известно, что в реальных экспериментах получают данные о поведении квантового прибора в зависимости от начальных условий и параметров задачи [9-11]. Например, для потокового кубита это интерференционные зависимости вероятностей населенности Р^(А, е0) возбужденного уровня кубита от амплитуд постоянного е0 и переменного А магнитных полей, которые в качестве примера рассчитываются в рамках данной работы. Так как значения вероятностей при каждых значениях параметров е0 и А независимы, то возможно повышение производительности работы программы за счет использования кластерной системы, узлы которой содержат ОРИ, на каждом из которых значение вероятности населенности рассчитывается независимо. Для обеспечения работы нескольких ОРИ была использована технология программирования для систем с распределенной памятью МРІ.
Для использования всех имеющихся графических ускорителей работа с данными организована следующим образом (рис. 2):
- перед началом выполнения алгоритма создаются структуры, которые хранят данные для работы генератора случайных чисел (начальное смещение в генерируемой последовательности для конкретного устройства и общие данные генератора), и делается несколько копий структуры с информацией о гамильтониане, которые используются на отдельных ОРИ;
- населённости уровней кубита рассчитываются в рамках вычислительного потока на графическом адаптере и, следовательно, не требуют выполнения дополнительных действий при переходе к множеству ОРИ;
- промежуточные результаты накапливаются в памяти каждого из используемых графиче-
ских адаптеров, а после окончания основного алгоритма суммируются на центральном процессоре.
Оптимизация работы программы
Одним из важных условий эффективности работы программы на ОРИ является оптимальное использование имеющихся типов памяти: глобальной, разделяемой и памяти констант. В частности, для решения данной задачи необходимы следующие наборы данных: информация о гамильтониане (хранится в регистрах разделяемой памяти), константы (хранятся в кешируемой памяти констант), текущие населённости уровней (хранятся в регистрах разделяемой памяти) и результирующие данные (имеют большой объем, хранятся в глобальной памяти). Регистры являются быстрым типом памяти, но имеют небольшой объем, поэтому использование их для работы с результирующими даными не представляется возможным. Глобальная память, напротив, позволяет хранить большие объемы данных, но является относительно медленной. Однако в данной задаче обращений к результирующим данным массивов, находящимся в глобальной памяти, немного, поэтому хранение в глобальной памяти не оказывает существенного влияния на производительность программы. Также отметим, что если при усовершенствовании архитектуры ОРИ объем разделяемой памяти будет увеличен, то существует возможность хранить результирующие данные в данной памяти ОРИ, что позволит ещё более эффективно использовать обращение к памяти и ускорить расчет.
Заметим, что точность квантового метода
МК зависит от числа реализаций ~ 1/4М. В силу наличия стохастических процессов двойная точность для расчётов не требуется, а использование операций с одинарной точностью на ОРИ существенно ускоряет работу приложения, т.к. вычисления с двойной точностью на несколько порядков медленнее [14].
Для обеспечения сходимости результатов моделирования диссипативной динамики необходимо наличие большого числа реализаций (М > 103) и генерирование последовательностей случайных чисел непосредственно для каждой точки временной траектории. Генерация необходимого массива случайных чисел производилась на ОРИ с помощью возможностей СИБЛ-библиотеки - СИЯЛКО. Использование данной библиотеки позволило существенно увеличить скорость выполнения расчётов, так как скорость
Рис. 1. Алгоритм моделирования диссипативной динамики потокового кубита, основанного на квантовом методе Монте-Карло
Управляющий CFU/<y~^
Запись в Global рассчитанных данных
Выделение памяти Global и Constant
ори ори ори ори
ори
Выделение Shared памяти Выделение Shared памяти
задание основного ... задание основного
состояния в t = 0 состояния в t = 0
Инициализация потоков на каждом ЄРІ/ Число реализаций - 1 ... М
Динамика Динамика Динамика Динамика Динамика
системы (МК) системы (МК) системы (МК) системы (МК) системы (МК)
Сбор данных по реализациям
Рис. 2. Алгоритм моделирования вероятности населенности возбужденного уровня кубита, где блоки «Динамика системы (МК)» реализуются согласно алгоритму, приведенному на рис. 1. Пунктирной линией выделена схема алгоритма для вычисления на ОРИ-ускорителе одной усредненной квантовой траектории
генерации случайных чисел на ОРИ в 50 раз больше, чем на СРИ [15].
Результаты расчетов и ускорение при использовании архитектуры СиБЛ
Для расчёта единичных усредненных траекторий была выбрана длительность импульса т = = 5Т, где Т - период внешнего поля, а шаг по времени М = 0.003Т выбирался, чтобы как удовлетворить условиям метода МК [8], так и обеспечить точность метода Рунге - Кутта четвертого порядка, применявшегося при численном моделировании динамики между скачками. При данных параметрах наличие статистических ошибок не должно превышать 5%, что приемлемо для извлечения необходимой физической информации из такого рода систем.
На рис. 3 приведен пример квантовой траектории, описывающей динамику сверхпроводящего кубита в сильном поле с учётом релаксационных эффектов при фиксированных значениях управляющего параметра е0, ответственного за постоянное магнитное поле в системе, и амплитуде А переменного поля, измеряемых в ГГц.
Р
0.8 0.6 0.4 0.2
0 1 2 3 4 t/T
Рис. 3. Населенность верхнего уровня кубита (усредненная квантовая траектория)
В табл. 1 представлены характерные времена, требуемые для вычисления с одинарной точностью населенности верхнего уровня кубита на CPU (Intel Core i7 960, использовалось 1 ядро) и на GPU (Nvidia Tesla С1060). Версия для GPU создана с использованием технологии CUDA. Видно, что графические ускорители позволяют существенно сократить время расчета, и поэтому их использование является целесообразным при решении данной задачи. Как видно из табл. 1, ускорение относительно CPU растет в зависимости от числа реализаций и достигает примерно 170 раз. Рост ускорения от реализаций обусловлен тем, что время расчета одной траектории мало, а при «перегрузке» графического ускорителя наступает эффективная загруженность и достигается оптимальная работа.
Таблица 1
Время расчета единичной усредненной по реализациям квантовой траектории на СРи и ОРИ
Число реализаций Время расчета на 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
Для расчета интерференционной картины (рис. 4), наблюдавшейся в недавних экспериментах [9-11], была выбрана 2Б-сетка (А, е0) 100^100 точек; таким образом, было необходимо рассчитывать 10000 усредненных квантовых траекторий (каждая из которых была усреднена по 15360 реализациям) и, соответственно, вероятность населенности возбужденного уровня в момент окончания длительности импульса.
Во, ГГц
5 0
-5
-10
0 2 4 6 8 10 12 А, ГГц
0 0.5 1.0
Рис. 4. Результаты расчёта интерференционной картины для возбужденного уровня системы
На рис. 5 и в табл. 2 представлены данные о приросте производительности вычислений в зависимости от количества используемых графических ускорителей.
Количество графических ускорителей
Рис. 5. Ускорение работы алгоритма в зависимости от числа использованных ОРИ по сравнению с одним ОРИ
Зависимость аппроксимируется линейной, что говорит о хорошем масштабировании задачи при переходе к множеству ОРИ, эффективность при использовании 20 ОРИ составляет 97%.
Таблица 2
Время расчета и ускорение алгоритма при одновременном расчёте 10000 единичных усредненных траекторий на ОРИ с применением технологии МР1
Тип устройства Число устройств Время, с Ускорение
GPU (Tesla C1060) 1 6020 ~ в 66 раз отн. СРИ
GPU (Tesla C1060) + MPI 5 1228 ~ в 323 раза отн. СРИ, ~ в 5 раз отн. ОРИ
GPU (Tesla C1060) + MPI 20 311 ~ в 1273 раз отн. СРИ, ~ в 19.4 раза отн. ОРИ
Выводы
В работе представлены математическая модель и алгоритмы компьютерного моделирования, позволяющие описать процессы возбуждения кубитов сильными электромагнитными полями. Создан программный комплекс, который реализует параллельные вычисления, имитирующие процессы измерения вероятности нахождения кубита в возбужденном состоянии, на многопроцессорном кластере (с использованием протокола МРІ) и графических процессорных устройствах (с использованием технологии СИБЛ). Используя разработанный программный комплекс, выполнено численное моделирование диссипативной динамики одиночных квантовых систем на основе метода квантовых траекторий (квантового метода Монте-Карло) на примере кубитов во внешних полях. Комплекс поддерживает суперкомпьютер-ные технологии на базе ОРИ-ускорителей и кластерное распараллеливание (МРІ). Разработанные методики позволяют из первых принципов промоделировать результаты недавних экспериментов (например, «однократные» измерения с помощью нелинейного бифуркационного осциллятора [4-6]), а также рассмотреть влияние квантового и классического шума на дина-
мику электронных состояний в квантовых точках и предложить методы управления и томографии состояния кубита. На основе результатов теоретических расчетов и компьютерного моделирования может быть разработана методика оптимального подбора параметров сигналов (частоты, амплитуды, длительности импульса) для необходимых процессов записи и считывания информации из кубита [16, 17].
Работа выполнена при поддержке грантов РФФИ № 12-07-00546-а и №12-07-31144-мол_а; госконтракта Минобрнауки № 07.514.11.4147 и соглашения ФЦП «Кадры» №14.132.21.1399.
Список литературы
1. Moore G.E. Cramming more components onto integrated circuits // Electronics. Volume 38. Number 8. April 19, 1965. URL: ftp://download.intel.com/muse-um/Moores_Law/Articles-Press_Releases/Gordon_Moo-re_ 1965_rticle.pdf
2. Wikipedia - Moore's_law. URL: http://en.wiki-pedia.org/wiki/Moore's_law
3. Нильсен М., Чанг И. Квантовые вычисления и квантовая информация. М.: Мир, 2006. 824 с.
4. Neumann P., Beck J., Steiner M. et al. // Science. 2010. V. 329. P. 542.
5. Vamivakas A.N., Lu C.-Y., Matthiesen C. // Nature. 2010. V. 467. P. 297.
6. Vijay R., Slichter D.H., Siddiqi I. // Phys. Rev. Lett. 2011. V. 106. P. 110502.
7. You J.Q., Nori F. // Nature. 2011. V. 474. P. 589.
8. Plenio M.B., Knight P.L. // Rev. Mod. Phys. 1998. V. 70. P. 101.
9. Oliver W.D. et al. // Science. 2005. V. 310. P. 1653.
10. Berns D.M. et al. // Phys. Rev. Lett. 2006. V. 97. P. 150502.
11. Oliver W. et al. // Quantum Inf. Process. 2009. V. 8. P. 261.
12. Скалли M.O., Зубайри M.C. Квантовая оптика. М.: Физматлит, 2003.
13. Боресков А.В., Харламов А. А. Основы работы с технологией CUDA. Изд-во «ДМК Пресс», 2010.
14. http://www.ixbt.com/video3/gt200-part1.shtml (дата обращения: 14.10.2012).
15. http://www. nag.co.uk/numeric/GPU s/benchmarks. asp (дата обращения: 14.10.2012).
16. Denisenko M.V., Satanin A.M., Ashhab S. et al. // Phys. Rev. B. 2012. V. 85. P. 184524.
17. Гельман А.И., Сатанин А.М. // Письма в ЖЭТФ. 2010. Т. 91. С. 584-589.
AN ALGORITHM FOR CALCULATING QUBIT DISSIPATIVE DYNAMICS ON A SUPERCOMPUTER COMPLEX USING GPU ACCELERATORS
A. V. Volkov, M.V. Denisenko, A.M. Satanin
The dissipative qubit dynamics in a strong electromagnetic field has been studied by the quantum Monte-Carlo method. The simulation algorithm has been realized on a heterogeneous supercomputing system using CUDA and MPI technologies. A method proposed to parallelize and accelerate the calculation scheme has demonstrated ~ 170 times greater acceleration as compared to the CPU sequential version. The developed software complex allows calculating and predicting the results of expensive experiments with qubits.
Keywords: quantum Monte-Carlo method, GPU, CUDA technology, parallel computing, qubit.
