Метод Монте-Карло на графических процессорах

Министерство образования и науки Российской Федерации
Уральский федеральный университет
имени первого Президента России Б. Н. Ельцина

К. А. Некрасов, С. И. Поташников, 
А. С. Боярченков, А. Я. Купряжкин

Метод Монте-Карло 
на графических процессорах

Учебное пособие

Рекомендовано методическим советом УрФУ
для студентов, обучающихся по направлениям подготовки
14.04.02 — Ядерная физика и технологии;
09.04.02 — Информационные системы и технологии;
14.04.01 — Ядерные реакторы и материалы

Екатеринбург
Издательство Уральского университета
2016

УДК 519.245:004(075.8)
ББК 22я73+32.97я73
          М54
Авторы:
Некрасов К. А., Поташников С. И., Боярченков А. С., Купряжкин А. Я.

Рецензенты:
Институт теплофизики УрО РАН (д‑р физ.‑мат. наук, проф. В. Г. Байдаков); гл. науч. сотр. лаборатории математического моделирования Института промышленной экологии УрО РАН д‑р физ.‑мат. наук, проф. 
А. Н. Вадаксин

М54
    Метод Монте-Карло на графических процессорах : учебное пособие / К. А. Некрасов, С. И. Поташников, А. С. Боярченков, А. Я. Купряжкин. — Екатеринбург : Изд‑во Урал. ун‑та, 2016. — 60 с.

ISBN 978‑5‑7996‑1723‑3

В учебном пособии изложены основные идеи метода Монте‑Карло, принципы 
распараллеливания расчетов, организации высокоскоростных параллельных вычислений методом Монте‑Карло на графических процессорах NVIDIA архитектуры CUDA. Подробно проанализирован пример моделирования методом Монте‑Карло диффузии нейтронов через пластину.
Пособие предназначено для проведения практических занятий по программированию графических процессоров для магистрантов по направлениям подготовки 09.04.02, 14.04.02, специалистов по направлению подготовки 141401.

Библиогр.: 15 назв. Рис. 15. Прил. 1.

УДК 519.245:004(075.8)
ББК 22я73+32.97я73

ISBN 978‑5‑7996‑1723‑3 
© Уральский федеральный
 
     университет, 2016

Введение

Одним из широко используемых подходов к численному моделированию физических систем является метод Монте‑Карло [1]–[2], для 
применения которого необходим анализ большого количества случайных состояний и вариантов поведения исследуемой системы, что связано с большими объемами вычислений.
Варианты случайного поведения системы, рассматриваемые методом Монте‑Карло, обычно генерируются и обрабатываются по одному и тому же общему алгоритму, будучи при этом независимыми друг 
от друга. Это позволяет эффективно реализовывать метод Монте‑Карло на поточно‑параллельных вычислительных системах.
В последние годы появилась практическая возможность проведения поточно‑параллельного физико‑математического моделирования 
на графических процессорах (GPU) — вычислительных устройствах 
параллельной архитектуры, изначально предназначенных для отображения графики на персональных компьютерах и рабочих станциях [3].
В настоящее время ведущие производители графических процессоров сами позиционируют новые GPU как универсальные системы, 
предназначенные не только для графических вычислений, но и для расчетов общего назначения. В частности, компанией NVIDIA выпускаются графические процессоры архитектуры CUDA [4] и соответствующее программное обеспечение [5], которые позволяют сравнительно 
легко программировать GPU на языке, являющемся расширением 
языка C, сохраняющем синтаксис последнего и включающем все его 
основные возможности.
Современные графические процессоры имеют в своем составе 
до нескольких сотен параллельных процессоров, каждый из которых 
по вычислительной производительности сравним с центральным процессором персонального компьютера [6]– [7].

ВВедение

Опыт показывает, что при решении хорошо распараллеливаемых 
задач реальная производительность GPU оказывается не менее чем 
на два порядка выше, чем у современных центральных процессоров 
ПК. Отметим, что затраты на приобретение и эксплуатацию одного 
вычислительного комплекса с GPU примерно совпадают с затратами на обычный персональный компьютер, тогда как кластер эквивалентной производительности занимал бы намного больше места, будучи при этом намного сложнее в эксплуатации и на порядки дороже.
В работах [8]–[10] представлены результаты применения графических процессоров для моделирования кристаллов диоксида урана 
методом молекулярной динамики. Использование графических процессоров позволило ускорить решение данной задачи (по сравнению 
с центральным процессором) более чем в 600 раз. Настоящее пособие 
посвящено оценке возможностей графических процессоров архитектуры CUDA (Compute Unified Device Architecture [3]–[5]) на примере метода Монте‑Карло, в частности при моделировании одномерной диффузии нейтронов через пластину.
Архитектуры графических процессоров, выпускаемых ведущими 
производителями — компаниями NVIDIA и AMD/ATI, несколько 
различны, что особенно проявляется при неграфических вычислениях. Наиболее удобной для таких вычислений мы считаем архитектуру 
CUDA [3]–[5], предлагаемую компанией NVIDIA.
В настоящем пособии будут рассмотрены особенности параллельной реализации метода Монте‑Карло для высокопроизводительных 
вычислений на графических процессорах. При этом мы ограничимся 
программированием для CUDA, так как в настоящее время эта платформа представляется наиболее эффективной, а реализованный в ней 
подход — наиболее перспективным.

1. Применение метода Монте-Карло  
для моделирования физических систем

1.1. Содержание метода Монте-Карло

1.1.1. Метод Монте-Карло на примере вычисления площадей

Простейшей иллюстрацией метода Монте‑Карло [1] является вычисление площади плоской фигуры (S) по принципу, показанному 
на рис. 1.1 [2].
Фигура заключена в квадрат, по которому 
случайным образом разбрасывается большое 
количество точек. Если число точек достаточно велико, то отношение площади S к площади квадрата L 2 будет с заданной точностью равняться отношению количества точек 
ў
N , 
попавших внутрь фигуры, к полному количеству точек N

S
L N

N
»
ў
2
.

Точность оценки площади возрастает вместе с количеством точек N. Погрешность вычислений обычно пропорциональна N ‑1/2 .
Площадь фигуры можно интерпретировать как ее объем в двухмерном пространстве. Очевидно, что метод, показанный на рис. 1.1, 
может быть применен и для вычисления объема трехмерной фигуры, 
заключенной внутрь куба. Аналогично можно рассчитывать объемы 
многомерных фигур в пространствах бόльшей размерности. Этот подход может быть без модификации применен и для многомерного численного интегрирования, поскольку определенные интегралы могут 
быть представлены как объемы многомерных тел.

L

S

Рис. 1.1. Вычисление 
площади плоской фигуры с криволинейной 
границей методом Монте‑Карло [2]

1. Применение метода монте-Карло для моделироВания физичесКих систем 

1.1.2. Распределения случайных величин

Одним из основных этапов любой вариации метода Монте‑Карло является генерирование случайных величин, распределенных 
по необходимому закону. При вычислении площади фигуры в параграфе 1.1.1 этими случайными величинами были координаты точек, 
разбрасываемых по площади квадрата. Каждой точке соответствовали две координаты (x, y), каждая из которых должна была быть случайной величиной, равномерно распределенной на интервале (0; L).
Пусть некоторая случайная величина x определена на интервале 
(a; b). Плотностью распределения x называют функцию w (x), такую, что

W x
w x dx
P
x

a

x
( ) =
( )
=
<
(
)
т
x
 —

вероятность того, что x окажется меньше, чем x. Сам интеграл W(x) называется интегральной (кумулятивной) функцией распределения величины x, значения W(x) всегда принадлежат интервалу (0; 1). 
Поскольку P
b
g Ј
(
)  = 1, для любой плотности распределения выполняется условие нормировки

 
w x dx

a

b
( )
=
т
1. 
 (1.1)

Несколько нестрого равномерно распределенными можно назвать 
те случайные величины, все допустимые значения которых равновероятны. Соответственно равномерно распределенные величины характеризуются плотностью распределения, не зависящей от x:

w x
b
a
x
a b

x
a b
( ) =
О(
)

П(
)

м

нп

оп

1
, при
;

0,
при
;

;

.

Такое распределение имеет прямоугольную форму, как это показано для интервала (0; 1) на рис. 1.2.

Кроме равномерно распределенных, 
широко используются и часто рассматриваются нормально распределенные (распределенные по Гауссу) случайные веw(x)

1

x
0
1
0

Рис. 1.2. Плотность вероятности для равномерного распределения на интервале (0; 1)

1.1. содержание метода монте-Карло

личины. Нормальной (гауссовой) случайной величиной называется 
случайная величина x, определенная на всей числовой оси ‑Ґ Ґ
(
)
;
 
с плотностью распределения

 
w x
x
( ) =
(
)
й

л

к
к

щ

ы

ъ
ъ
1

2
2

2

2
s
p

x

s
exp
,

где s — дисперсия x , характеризующая разброс значений относительно среднего; x  — математическое ожидание (среднее значение) x.
Среднее значение произвольной непрерывной случайной величины x с плотностью распределения w (x) может быть вычислено по формуле

 
x =
( )
т xw x dx

0

1

,

а дисперсия определяется по выражению

 
s
x
x
2
2

0

1
2
[ ] =
( )
‑( )
т x w x dx
.

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

Стандартные генераторы случайных чисел (см. подглаву 1.2) обычно создают числовые последовательности с равномерным распределением. В то же время для решения конкретных задач методом Монте‑Карло могут потребоваться случайные величины, распределенные 
по самым разным законам. Существуют достаточно простые способы 
перехода от равномерно распределенной случайной последовательности к последовательностям с любыми другими распределениями. 
Остановимся на одном из них.
Пусть необходимо генерировать значения случайной величины x, 
распределенной в интервале (a; b) с плотностью вероятности w(x) 
и функцией распределения W(x). Значения функции распределения 
W(x) обязательно приходятся на интервал (0; 1), что позволяет сопоставить ей случайную величину g, равномерно распределенную на этом же 
интервале.