Приволжский научно-образовательный центр суперкомпьютерных технологий
Исследования
  • Основные направления
  • Гранты, проекты, договора

    • Методы распараллеливания вычислений для решения задач квантовой механики

    Проект «Методы распараллеливания вычислений для решения задач квантовой механики» по гранту Российского Фонда Фундаментальных Исследований №18-37-00277.

    1. Цель и задачи фундаментального исследования

    Основной целью проекта является разработка новых параллельных алгоритмов и их высокопроизводительных программных реализаций для численного моделирования динамики квантовых систем на современных суперкомпьютерах, узлы которых содержат мультиядерные процессоры, включая процессоры Intel Xeon Phi второго поколения. Рассмотрены следующие методы численного моделирования: метод Магнуса, метод квантовых траекторий и метод перехода к решению системы ОДУ в базисе обобщенных матриц Гелл-Манна.

    Достижение этой цели требует решения следующих задач:

    1. Аналитический обзор и анализ методов эффективного использования мультиядерных процессоров для численного решения задач квантовой механики.
    2. Апробация разработанных ранее коллективом авторов программных реализаций метода Магнуса, метода квантовых траекторий и метода перехода к решению ОДУ в базисе обобщенных матриц Гелл-Манна при решении вычислительно-трудоемких задач на суперкомпьютерах, построенных на основе мультиядерных процессоров.
    3. Разработка методов балансировки нагрузки для эффективного использования параллелизма на трех уровнях вычислительной системы (распределенная память, общая память, инструкции SIMD).
    4. Программная реализация разработанных методов.
    5. Анализ корректности подготовленных реализаций путем анализа результатов моделирования на наборе тестовых задач.
    6. Анализ производительности и масштабируемости разработанных программных реализаций.

    2. Важнейшие результаты, полученные при реализации Проекта

    В работе рассматривались задачи исследования динамики, как закрытых квантовых систем, описываемых уравнением Шредингера, так и открытых квантовых систем с использованием марковского кинетического уравнения в форме «Горини-Коссаковского-Сударшана-Линдблада» (далее просто «уравнение Линдблада») [7].

    Для исследования динамики закрытых систем необходимо было провести апробацию разработанной ранее программной реализация метода Магнуса [8]. Для решения уравнения Линдблада необходимо было использовать разработанную ранее программную реализацию метода квантовых траекторий [9], а также разработать новую программную реализацию метода перехода к решению ОДУ в базисе обобщенных матриц Гелл-Манна [1].

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

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

    В результате выполнения проекта был подготовлен обзор методов численного решения задач квантовой механики. Обзор показал наличие высокого интереса к разработке и применению новых методов, способных интегрировать системы общего вида, а также к созданию методов, учитывающих специфику конкретных приложений. Увеличение производительности и объема памяти современных суперкомпьютеров открывает новые перспективы численного моделирования для исследования динамики сложных многокомпонентных квантовых систем. Для этого требуется разработка методов эффективного использования ресурсов суперкомпьютеров. С учетом выполненного обзора литературы была проиведена разработка трех программных реализаций методов численного моделирования: метода Магнуса, метода квантовых траекторий и метода перехода к решению системы ОДУ в базисе обобщенных матриц Гелл-Манна.

    2.2. Апробация разработанных ранее коллективом авторов программных реализаций при решении вычислительно-трудоемких задач на суперкомпьютерах, построенных на основе мультиядерных процессоров

    Метод Магнуса и метод квантовых траекторий

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

    Для метода Магнуса выполнялся расчет Флоке-состояний для димера Бозе-Хаббарда с N = 1000 бозонами, измерялось полное время работы метода (меньше лучше).

    В методе квантовых траекторий выполнялось интегрирование 1000 траекторий для системы с N = 1000 состояниями, измерялось число траекторий, проинтегрированных на 1 период за 1 секунду (больше лучше).

    Тесты запускались на двух системах:

    1. Суперкомпьютере Intel Endeavor. Вычислительные узлы обладали следующими характеристиками: 2x Intel Xeon Gold 6148 (40 ядер, 2.40 GHz), более чем 128 GB оперативной памяти.
    2. Суперкомпьютере МСЦ РАН. Вычислительные узлы обладали следующими характеристиками: Intel Xeon Phi 7290 (72 ядра, 1.50 GHz).

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

    Для анализа производительности строилась Roofline-модель. Сбор данных для Roofline-модели производился с помощью инструмента Intel Advisor из пакета Parallel Studio XE 2018.

    Анализ корректности разработанных программных реализаций производился как задел, до начала выполнения проекта. Некоторые результаты по анализу корректности программных реализаций приведены в разделе 2.5.

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

    Метод перехода к решению системы ОДУ в базисе обобщенных матриц Гелл-Манна

    Для апробация программной реализации метода перехода к решению ОДУ в базисе обобщенных матриц Гелл-Манна в начале проверялась ее корректность. Для этого были построены бифуркационные диаграммы для квантовой модулированной системы размера N = 1000  при достижении 10 периодов по времени методом квантовых траекторий и рассматриваемым методом. Бифуркационные диаграммы использовались для качественного сравнения получаемых результатов. Численное сравнение получаемых результатов было выполнено при решении систем меньшего размера (N = 200 и N = 400). Результаты проверки корректности полученной программной реализации представлен в разделе 2.5.

    Кроме анализа корректности для метода решения ОДУ в базисе обобщенных матриц Гелл-Манна анализировались:

    1. Затраты памяти в зависимости от размера задачи N.
    2. Масштабируемость полученных параллельных программных реализаций.

    Численные эксперименты проводились на суперкомпьютерах:

    1. «Лобачевский», установленном в ННГУ им. Н.И. Лобачевского. Вычислительные узлы обладали следующими характеристиками: 2×Intel Xeon E5 2660 (8 ядер, 2.2 GHz), 64 GB оперативной памяти, операционная система CentOS 6.4.
    2. Суперкомпьютере Intel Endeavor. Вычислительные узлы обладали следующими характеристиками: 2x Intel Xeon Gold 6148 (40 ядер, 2.40 GHz), более чем 128 GB оперативной памяти.
    3. Суперкомпьютере МСЦ РАН. Вычислительные узлы обладали следующими характеристиками: Intel Xeon Phi 7290 (72 ядра, 1.50 GHz).

    Использовались программные пакеты Intel Math Kernel Library, Intel C/C++ Compiler и Intel MPI, а также Intel Advisor в составе Intel Parallel Studio.

    Результаты вычислительных экспериментов приведены в разделе 2.6.

    2.3. Разработка методов балансировки нагрузки для эффективного использования параллелизма на трех уровнях вычислительной системы

    Метод Магнуса

    Реализация метода подробна изложена в [8]. Приведем основные этапы метода Магнуса:

    1. Вычисление базиса собственных векторов матрицы, вспомогательных матриц, коэффициентов функции Бесселя и коммутаторов для разложения Магнуса.
    2. Интегрирование единичной матрицы на время периода T посредством применения последовательности M одношаговых разложений Магнуса на временных интервалах размера T/M.  Используемый в разложении Магнуса экспоненциальный оператор, в свою очередь, вычислялся через полиномиальное разложение Чебышева.
    3. Численная диагонализация матрицы, полученной на этапе 2.

    Схематично реализация метода представлена на рис. 1.

    Рис. 1. Вычислительная схема метода Магнуса

    Наиболее вычислительно-трудоемким является этап приближения матричной экспоненты полиномиальным представлением Чебышева, ядром которого является операция умножения матрицы на вектор.

    В работе [8] представлен подход, позволяющий перейти от матрично-векторных операций к матричным операциям. В указанном алгоритме распараллеливание было выполнено для расчетов на кластерных системах. Задействованы три уровня параллелизма: распределенная память, общая память, векторные вычисления (команды SIMD). Для разделения вычислительной нагрузки в системах с распределенной памятью используются стандартные механизмы технологии MPI. В рамках шага реализовано равномерное распределение блоков матриц, состоящих из базисных векторов, на несколько узлов. На общей памяти и на уровне инструкций SIMD использовались оптимизированные реализации функций BLAS из библиотеки Intel Math Kernel Library (MKL). В качестве основных операций использовались операции матричного умножения и сложения векторов. Указанный способ распределения вычислительной нагрузки позволяет равномерно загрузить все вычислительные устройства вычислительной системы.

    В работе [4] представлена разработанная программная реализация метода Магнуса для графических ускорителей.  При реализации метода Магнуса на графических ускорителях интегрируемая единичная матрица базисных векторов разбивается на блоки столбцов равного размера согласно общему числу используемых GPU. На каждый GPU загружается его блок начальных векторов. Метод и алгоритм адаптированы таким образом, что все шаги являются операциями над комплексными матрицами, которые выполняются непосредственно на GPU с помощью функций cuBLAS. По окончании интегрирования на время периода результат копируется из памяти GPU в RAM и собирается на одном процессе с использованием механизмов из технологии MPI.

    Метод квантовых траекторий

    Идея метода заключается в усреднении по многим реализациям результатов интегрирования во времени волновой функции с помощью неэрмитового гамильтониана специального вида [9]. В случае кусочно-линейной правой части уравнения Линдблада использовались две ускоряющие модификации метода:

    1. Интегрирование во времени волновой функции с помощью предварительно вычисленного набора экспоненциальных операторов.
    2. Одновременное интегрирование во времени групп траекторий, позволяющее использовать вместо операций умножения матрицы на вектор матричные умножения. В этом случае 95% времени вычислений выполняются операции умножения комплексных матриц.

    Общая вычислительная схема программной реализации представлена на рис.2.

    Рис. 2. Вычислительная схема метода квантовых траекторий

    Распараллеливание выполнено для расчетов на кластерных системах, задействованы три уровня параллелизма: распределенная память, общая память, векторные вычисления (команды SIMD). Благодаря естественному параллелизму алгоритмов Монте-Карло модифицированный метод распараллеливается на несколько узлов с линейной эффективностью. На уровне общей памяти и инструкций SIMD использовались оптимизированные реализации функций BLAS из библиотеки Intel Math Kernel Library (MKL).

    Метод перехода к решению системы ОДУ в базисе обобщенных матриц Гелл-Манна

    Реализация метода перехода к решению системы ОДУ в базисе обобщенных матриц Гелл-Манна содержит следующие основные этапы:

    1. Инициализация данных. В рамках первого этапа определяется решаемая система ОДУ. Этап не оказывает сильного влияния на общее время вычислений.
    2. Подготовка данных. В рамках этапа исходная система ОДУ переводиться в базис обобщенных матриц Гелл-Манна.
    3. Интегрирование системы ОДУ. Этап занимает основное время вычислений.
    4. Завершение вычислений. На этапе завершения вычислений полученное на третьем этапе решение переводится в исходный базис. Этап не оказывает сильного влияния на общее время вычислений.

    Отметим следующие характерные особенности алгоритма:

    Подготовка данных (этап 2) занимает незначительную долю времени относительно времени работы всего алгоритма. При этом данный этап очень требователен по памяти, что и ограничивает размер возможных решаемых систем. В рамках проекта предложен способ распараллеливания данного этапа, основной целью которого является распределение данных по узлам суперкомпьютера при решении задач большой размерности [1]. Достигаемое при этом ускорение вычислений также является важным, но не первостепенным фактором, поскольку подготовка данных выполняется только один раз, а основные расчеты производятся на этапе интегрирования системы (этап 3).

    Алгоритм интегрирования системы ОДУ методом Рунге-Кутта (этап 3) вносит основной вклад во время решения задачи. В рамках проекта предполагалось, что матрицы, описывающие получаемую систему ОДУ, являются существенно разреженными. Вместе с тем наличие зависимостей по данным в ходе численного интегрирования на каждом шаге метода потребовало разработки параллельного алгоритма, учитывающего эти зависимости.

    В реализации метода Рунге-Кутта основные вычисления сосредоточены в матрично-векторных умножениях и в суммировании векторов. Распараллеливание выполнено для расчетов на кластерных системах, задействованы три уровня параллелизма: распределенная память, общая память, векторные вычисления (команды SIMD).

    Разделение вычислительной нагрузки в системах с распределенной памятью использует стандартные механизмы технологии MPI, на общей памяти и на уровне инструкций SIMD использовались оптимизированные реализации функций BLAS из библиотеки Intel Math Kernel Library (MKL).

    Для организации распределенных вычислений использовался следующий подход. В начале матрицы решаемых систем ОДУ равномерно распределялись между вычислительными процессами. Разделение матриц осуществлялось по строкам. Далее использовалась библиотека ParMetis [11] для поиска такого разделения строк между процессами, которое минимизирует количество пересылок. Суть данной оптимизации заключается в следующем. На каждом шаге метода Рунге-Кутта процессы получают свою часть интегрируемого вектора. Чтобы продолжить вычисления, процессы обмениваются вычисленными значениями. Поскольку матрицы являются разреженными, то каждому процессу нужны не все вычисленные значения. Перераспределение строк между процессами призвано уменьшить объем передаваемых данных.

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

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

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

    2.4. Программная реализация разработанных методов

     Метод Магнуса

    Метод Магнуса реализован для вычислительных систем как с общей памятью, так и с распределенной памятью. Программная реализация допускает использования Intel Xeon Phi и GPU. Детали реализации изложены в публикациях [2,4,8].

    Метод квантовых траекторий

    Метод квантовых траекторий реализован для вычислительных систем как с общей памятью, так и с распределенной памятью. Программная реализация допускает использования ускорителей Intel Xeon Phi. Детали реализации изложены в публикациях [2,9]. В 2019 году по программной реализации получено свидетельство о регистрации программы для ЭВМ [6].

    Метод перехода к решению системы ОДУ в базисе обобщенных матриц Гелл-Манна

    В рамках проекта разработан новый метод перехода в базис обобщенных матриц Гелл-Манна (см. раздел 2.3). На основе метода разработана программная реализация, допускающая вычисления на системах как с общей, так и с распределенной памятью. Детали реализации представлены в публикации индексируемой Web Of Science [1]. В 2019 году по программной реализации получено свидетельство о регистрации программы для ЭВМ [5].

    2.5. Анализ корректности подготовленных реализаций путем анализа результатов моделирования на наборе тестовых задач

    Метод Магнуса

    Корректность реализации метода была подтверждена при подготовке публикации [8]. Реализация метода Магнуса для вычислений на графических ускорителях была протестирована путем сопоставления результатов полученной реализации и существующей корректной реализации для центральных процессоров.

     Метод квантовых траекторий

    Корректность реализации метода была подтверждена при подготовке публикации [10].

    Реализация метода квантовых траекторий способна работать с квантовыми системами с кусочно-постоянным драйвингом. У таких систем при небольшом размере можно получить решение построением супероператора. Для проверки корректности разработанной программной реализации проверялось, сходится ли решение, полученное методом квантовых траекторий к решению, полученному методом супероператора для квантового димера размера 100. Так же проверялся порядок сходимости схемы. Метод квантовых траекторий имеет теоретическую сходимость порядка 1/ \sqrt{M_r}, где M_r используемое число траекторий. Результаты проверки корректности представлены на рис. 3.


    Рис. 3. Матричная ошибка метода квантовых траекторий относительно метода супероператора для задачи квантового димера размера 100 в зависимости от изменения количества траекторий и времени интегрирования системы

    По осям отложено количество периодов интегрирования (t/T), количество квантовых траекторий, используемых в расчете (M_r), и матричная ошибка полученного результата от эталонного решения (ε). Матричная ошибка вычисляется как максимальный модуль собственных чисел разницы матриц.

    Исходя из рис 3. видно, что метод квантовых траекторий сходится со скоростью 1/ \sqrt{M_r} к решению, полученному методом супероператора, а величина ошибки сохраняется при интегрировании системы во времени.

    Метод перехода к решению системы ОДУ в базисе обобщенных матриц Гелл-Манна

    Проверка корректности разработанной программной реализации выполнялась в несколько этапов:

    1. Качественное сравнение решений, получаемых разными методами.
    2. Численное сравнение решений, получаемых разными методами.

    Для качественного сравнения решений была построена бифуркационная диаграмма для квантовой модулированной системы размера N = 1000 при достижении 10 периодов по времени методом квантовых скачков и рассматриваемым методом (см. рис. 4).

    (а) (б)

    Рис.4. Однопараметрическая бифрукационая диаграмма для квантовой модулированной системы полученной (а, слева) методом квантовых траекторий (8000 реализаций) и (б, справа) рассматриваемый метод

    Качественный анализ позволяет сделать следующие выводы. С вычислительной точки зрения недостаток метода перехода к решению ОДУ в базисе обобщенных матриц Гелл-Манна состоит в том, что расчеты до достижения аттрактора занимают продолжительное время и в основном определяются временем интегрирования системы ОДУ. Напротив, в методе квантовых траекторий стробоскопическое накопление позволяет производить вычисления достаточно далеко во времени существенно быстрее. Достоинство метода перехода к решению ОДУ в базисе обобщенных матриц Гелл-Манна в сравнении с методом квантовых траекторий заключается в его точности. Данный факт можно пронаблюдать на рис. 4. Как видно из картинок, они имеют общую структуру. На обоих графиках совпадает количество и время возникновения новых максимумов при переходе к квантовому хаотическому аттрактору. При этом диаграмма, посчитанная рассматриваемым методом перехода к решению ОДУ в базисе обобщенных матриц Гелл-Манна, очевидно, является более непрерывной.

    Численное сравнение решений, для подтверждения корректности реализаций было выполнено при подготовке публикации [1].

    Для проверки корректности сравнивалась сходимость двух методов: метода квантовых траекторий (QJ) и метода перехода к решению ОДУ в базисе обобщенных матриц Гелл-Манна (SU). В качестве тестовой задачи использовался квантовый димер для размеров 200 и 400. Результаты проверки корректности представлены на рис. 5.


    Рис.5. Сравнение сходимости методов по матричной норме для задачи квантового димера размеров 200 и 400

    Стоит отметить, что количество траекторий в методе квантовых траекторий и шаг интегрирования при решении системы в базисе обобщенных матриц Гелл-Манна линейно связаны с временем работы программ. При сравнении методов сопоставлялись решения, которые можно получить за одно и тоже время счета. По оси Х отложено время решения задачи в минутах. По оси Y отложена матричная ошибка относительно эталонного решения. Матричная ошибка вычисляется как максимальный модуль собственных чисел разницы матриц. Эталонное решение было получено благодаря длительному интегрированию системы по времени.

    Исходя из рис 5. видно, что решение системы, полученное в базисе обобщенных матриц Гелл-Манна, имеет точность на 8-10 порядков выше, чем решения, полученные методом квантовых траекторий за тоже время вычислений. Так же видно, что решение в базисе обобщенных матриц Гелл-Манна сходится к эталонному решению как \delta t^4, что соответствует методу Рунге-Кутты 4-го порядка.

    2.6. Анализ производительности и масштабируемости разработанных программных реализаций

    Метод Магнуса и метод квантовых траекторий

    Напомним, что для метода Магнуса выполнялся расчет Флоке-состояний для димера Бозе-Хаббарда с N = 1000 бозонами, измерялось полное время работы метода (меньше лучше). В методе квантовых траекторий выполнялось интегрирование 1000 траекторий для системы с N = 1000 состояниями, измерялось число траекторий, проинтегрированных на 1 период за 1 секунду (больше лучше).

    Результаты экспериментов частично представлены в Таблицах 1-2 (лучшие и показательные значения). Показатели Roofline приводятся для запуска с максимальной производительностью.

    Таблица 1. Время выполнения расчета состояний Флоке методом Магнуса (сек.)

    # процессов # MKL потоков Intel Xeon Gold Intel Xeon Phi
    1 40/72 498,9 614,2
    2 20/36 529,8 697,6
    4 10/18 467,2 673,5
    8 5/9 470,5 729,7

    Таблица 2. Число траекторий, проинтегрированных на 1 период за 1 секунду
    метода квантовых траекторий

    # процессов # OpenMP потоков # MKL потоков Intel Xeon Gold Intel Xeon Phi
    1 1 40/72 8,09 6,19
    1 40/72 1 15,31 6,22
    1 20/36 2 16,80 9,14
    2 1 20/36 15,04 7,50
    2 20/36 1 17,52 11,12
    4 10/18 1 13,75 10,63
    20/36 1 2 10,96 12,79
    20/36 2 1 11,84 13,44
    40/72 1 1 16,20 13,14

    В начале отметим максимальную производительность вычислительных систем используемых при проведении экспериментов – 3046,28 GFlops (Xeon Gold), 3230,5 GFlops (Xeon Phi).

    При использовании метода Магнуса Roofline-модель показывает арифметическую интенсивность (AI) реализации 1,32 FLOP/Byte и производительность до 1674,52 GFlops для системы Xeon Gold, 1,46 FLOP/Byte и 1027,52 GFlops для Xeon Phi.  Для сравнения, используемая реализация умножения комплексных матриц, функция zgemm() библиотеки MKL, для N=1000 показывает характеристики 1,34 FLOP/Byte и 1867,00 GFlops (Xeon Gold), 1,50 FLOP/Byte и 1239,37 GFlops (Xeon Phi). 

    При реализации метода квантовых траекторий отметим следующие характеристики реализации: 0,96 FLOP/Byte и 1406,22 GFlops на системе Intel Xeon Gold, 0,65 FLOP/Byte и 1284,78 GFlops на системе Intel Xeon Phi.

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

    Эксперименты с реализацией метода Магнуса для вычислительных систем с графическим ускорителями выполнялись на 16 узлах 2x Intel E5-2660, 2.2 GHz, 8 ядер, 64 GB RAM, 3x NVIDIA Kepler K20Х (2688 потоковых процессоров, 6 GB GDDR5), QDR Infiniband, Intel Parallel Studio XE 2017, CUDA Toolkit 8.0 [4].

    Выполнялся расчет Флоке-состояний для димера Бозе-Хаббарда с N бозонами.  Проведена серия экспериментов для систем с различной размерностью N, шагом по времени h = T/5000, количеством полиномов в разложении Чебышева L = 50.  Измерялось общее время выполнения расчета (столбец «Итого»), а также суммарные времена вычисления разложений Магнуса и Чебышева.  Ввиду того, что целью проведения экспериментов являлась оценка производительности, выполнялось 100 шагов интегрирования по времени вместо 5000.  Время работы для произвольного числа шагов можно оценить, линейно масштабируя время этапа интегрирования. Результаты представлены в Таблице 3.

    Таблица 3. Время выполнения расчета состояний Флоке (сек.)

    N 1 MPI-процесс, 1 GPU 1 MPI-процесс, 3 GPU 16 MPI-процессов, 3 GPU
    Магнус Чебышев Итого Магнус Чебышев Итого Магнус Чебышев Итого
    256 0.02 1.05 3.21 0.02 0.69 3.82 0.02 0.63 3.79
    512 0.08 6.10 9.42 0.08 2.79 7.11 0.08 1.27 5.49
    1024 0.29 41.32 48.30 0.29 15.31 22.99 0.29 2.57 10.46
    2048 1.14 308.37 333.28 1.14 107.02 133.43 1.14 12.08 38.47
    4096 2.54 1024.62 1094.46 2.53 344.49 419.19 2.54 25.12 102.67
    5120 Недостаточно памяти GPU 7.03 1584.63 1851.1 7.03 123.03 388.90

    Результаты экспериментов показывают, что в реальных расчетах этап вычисления Флоке-оператора будет занимать около 95% времени, эффективность его масштабирования при максимальных N составляет 75-80%.  Производительность каждого GPU на данном этапе достигает 1,135 TFlops (86,6% от теоретически возможных 1,31 TFlops). 

    Метод перехода к решению системы ОДУ в базисе обобщенных матриц Гелл-Манна

    Для разработанной программной реализации метода решения системы ОДУ в базисе обобщенных матриц Гелл-Манна [1] были выполнены следующие эксперименты:

    1. Проведен анализ потребляемой памяти в процессе вычислений.
    2. Проведена оценка масштабируемости программной реализации.

    Сначала рассмотрим результаты вычислительных экспериментов, показывающие затраты по памяти, необходимых для решения систем ОДУ в базисе обобщенных матриц Гелл-Манна. Решалась задача интегрирования по времени системы ОДУ описывающих димер.  Пиковое потребление памяти при решении задач размером от N = 100 до N = 2000 показано на рис. 6. Результаты экспериментов позволяют сделать вывод, что при решении задачи размерности N = 1600 требования к памяти на узел снижаются с 54,5 ГБ при использовании 5 узлов до 16,5 ГБ при использовании 25 узлов (эффективность масштабирования составляет 66,5%). Кроме того, мы смогли выполнить вычисления для N = 2000, для чего потребовалось 31 ГБ памяти на каждом из 25 узлов кластера. Время вычислений является значительным, но не критичным для этапа подготовки данных для разреженной задачи. Тем не менее, отметим, что при использовании 5 узлов кластера время вычислений уменьшается примерно вдвое по сравнению с последовательным алгоритмом, а затем время уменьшается незначительно.


    Рис.6. Количество потребляемой памяти на этапе подготовки данных в зависимости от количества используемых вычислительных узлов (исходная система ОДУ разреженная)

    Кроме ситуации, когда исходная система ОДУ описывалась разреженными матрицами, также решались системы ОДУ с плотными матрицами. Для плотной модели нам удалось осуществить переход в новый базис при N = 200. На рис. 7 слева показано, как уменьшаются затраты памяти на узел за счет увеличения числа узлов с 1 до 25. Обратите внимание, что в отличие от случая разреженной модели время, затрачиваемое на подготовку данных, уменьшается почти линейно, что является дополнительным преимуществом разработанной программной реализации.


    Рис.7. Количество потребляемой памяти на этапе подготовки данных в зависимости от количества используемых вычислительных узлов (исходная система ОДУ плотная)

    На рис. 8. показано время интегрирования системы в зависимости от количества MPI процессов. Отметим, что при проведении вычислительных экспериментов в рамках одного узла были задействованы все доступные вычислительные ядра. В случае, когда исходная система ОДУ была разреженной, использовалась задача с размерностью N = 1600, а в случае плотных систем N = 150. Исходя из результатов экспериментов можно сделать вывод о том, что в случае разреженных систем ОДУ использование больше 4 узлов является не целесообразным, а в случае плотных систем не целесообразно использование больше 2 узлов. На производительность алгоритмов при большом числе узлов существенную роль оказывают временные затраты на передачу данных.


    Рис.8. Время вычислений в зависимости от количества используемых вычислительных узлов

    3. Сопоставление результатов, полученных при реализации с мировым уровнем

    Результаты проведенных исследований сопоставимы с мировым уровнем, что подтверждается экспертизой результатов, проведенной при публикации в научных журналах и при представлении докладов на российских и международных научных конференциях [1-6].

    Эффективность и достоверность разработанных методов, математического и программного обеспечения показана благодаря выполненной репрезентативной серии вычислительных экспериментов.

    4. Методы и подходы, использованные при реализации Проекта

    Метод Магнуса

    Базовая реализация метода Магнуса предполагает формирование множества базовых векторов, а затем интегрирование полученного множества векторов по времени. Интегрирование по времени каждого вектора может быть сведено к матрично-векторным операциям. Для метода Магнуса была разработана модификация, позволяющая интегрировать по времени сразу множество базовых векторов. Благодаря модификации удалось перейти от операций матрично-векторного умножения к операциям матричного умножения. Для операций матричного умножения использовались оптимизированные реализации из Intel MKL и cuBLAS. Распределение базовых векторов по вычислительным узлам и сбор решений осуществлялся с применением технологии MPI.

    Метод квантовых траекторий

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

    Метод перехода к решению системы ОДУ в базисе обобщенных матриц Гелл-Манна

    В методе решения систем ОДУ в базисе обобщенных матриц Гелл-Манна были выделены два основных этапа:

    1. Этап подготовки данных.
    2. Этап численного интегрирования.

    Первый этап требует большого объема оперативной памяти. Второй этап занимает основное время вычислений.

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

    На этапе численного интегрирования был реализован распределенный метод Рунге-Кутта 4-го порядка. В случае, когда исходная система ОДУ является разреженной, произведена оптимизация затрат на передачу данных. Для оптимизации затрат исходная матрица системы уравнений ОДУ была переупорядочена с использованием библиотеки ParMetis. Также для снижения затрат на передачу данных при реализации активно использовались коллективные операции.

    Заключение

    В качестве заключения приведем основные результаты, полученные в рамках проекта.

    Во-первых, подготовлен обзор методов численного решения задач квантовой механики.

    Во-вторых, проведена апробация ранее разработанных программных реализаций метода Магнуса, метода квантовых траекторий. Проведенные эксперименты на суперкомпьютерах, построенных на основе мультиядерных процессоров, включая Intel Xeon Phi архитектуры Knights Landing показали эффективность данных реализаций.

    В-третьих, для метода перехода к решению ОДУ в базисе обобщенных матриц Гелл-Манна разработана как последовательная, так и параллельная реализации.

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

    В-пятых, написана и опубликована научная статья [1] в журнале Physical Review E индексируемым WoS. Написаны и опубликованы статья и тезисы [2-4]. Зарегистрированы две программы для ЭВМ [5,6].

    В-шестых, результаты выполненных исследований представлены на конференциях:

    • Конференция «Суперкомпьютерные дни в России», 24-25 сентября 2018 г., Москва.
    • Конференция «Суперкомпьютерные дни в России», 23-24 сентября 2019 г., Москва.
    • Конференция «Параллельные вычислительные технологии (ПаВТ) 2019», 2-4 апреля 2019 г., Калининград.

    Публикации

    1. Liniov A.V., Meyerov I.B., Kozinov E.A., Volokitin V.D., Yusipov I.I., Ivanchenko M.V., Denisov S. Unfolding a quantum master equation into a system of real-valued equations: Computationally effective expansion over the basis of SU (N) generators // Physical Review E. № 100. 2019. P. 053305.
    2. Волокитин В.Д., Козинов Е.А., Лаптева Т.В., Линев А.В. Применение Roofline-модели для анализа производительности двух алгоритмов численного моделирования квантовых систем// Суперкомпьютерные дни в России: Труды международной конференции (24-25 сентября 2018 г., г. Москва). – М.: Изд-во МГУ, 2018. – С. 1009-1011
    3. Козинов Е.А., Линев А.В., Волокитин В.Д., Иванченко М.В., Мееров И.Б., Денисов С. Параллельный алгоритм исследования динамики открытых квантовых систем: решение основного кинетического уравнения в базисе обобщенных матриц гелл-манна // Параллельные вычислительные технологии (ПАВТ’2019). Издательство: Издательский центр ЮУрГУ (Челябинск). 2019. С. 263-274.
    4. Красикова Е.А., Волокитин В.Д., Иванченко М.В., Лаптева Т.В., Линев А.В., Мееров И.Б. Параллельный алгоритм вычисления Флоке-состояний квантовой системы на гибридных кластерах с GPU // Суперкомпьютерные дни в России: Труды международной конференции. 23-24 сентября 2019 г., г. Москва / Под. ред. Вл.В. Воеводина – Москва : МАКС Пресс, 2019. C. 217-219
    5. Волокитин В.Д., Козинов Е.А. Численное решение основного кинетического уравнения в базисе обобщенных матриц Гелл-Манна на кластерных системах // Свидетельство о регистрации программы для ЭВМ №2019667407 от 23.12.2019.
    6. Волокитин В.Д., Козинов Е.А. Параллельная реализация метода квантовых траекторий для систем с распределенной памятью // Свидетельство о регистрации программы для ЭВМ №2019666749 от 13.12.2019.

    Литература

    1. Lindblad G. On the generators of quantum dynamical semigroups //Communications in Mathematical Physics. – 1976. – Т. 48. – №. 2. – С. 119-130.
    2. Laptyeva T. V. et al. Calculating Floquet states of large quantum systems: A parallelization strategy and its cluster implementation //Computer Physics Communications. – 2016. – Т. 201. – С. 85-94.
    3. Liniov A. et al. Increasing Performance of the Quantum Trajectory Method by Grouping Trajectories //Russian Supercomputing Days. – Springer, Cham, 2017. – С. 136-150.
    4. Volokitin V. et al. Computation of the asymptotic states of modulated open quantum systems with a numerically exact realization of the quantum trajectory method //Physical Review E. – 2017. – Т. 96. – №. 5. – С. 053313.
    5. Karypis G., Kumar V. ParMetis: Parallel graph partitioning and sparse matrix ordering library // Tech. Rep. TR 97-060, University Of Minnesota, Department Of Computer Science. 1997.