CC BY
25
УДК 004.021, 004.04
Н. Ю. Зятьков \ А. А. Айзенберг 2
1 Новосибирский государственный университет ул. Пирогова, 2, Новосибирск, 630090, Россия
2 Бергенский университет Аллегатен, 41, 7803 №-5020, Берген, Норвегия
nikolay.zyatkov@gmail. com; [email protected]
ВЫСОКООПТИМИЗИРОВАННАЯ РЕАЛИЗАЦИЯ ВЫЧИСЛЕНИЯ МАТРИЦЫ ТЕНИ ДЛЯ МОДЕЛИРОВАНИЯ КАСКАДНОЙ ДИФРАКЦИИ В ГЕОЛОГИЧЕСКИХ СЛОЯХ *
Представлены алгоритмы высокооптимизированного вычисления и хранения матрицы виртуальной тени. Процедура, численно реализующая эти алгоритмы, является частью разрабатываемого программного комплекса МНКВ (метод наложения концевых волн) для дифракционного моделирования волновых полей в слоистых и блоковых средах со сложными границами. Матрица виртуальной тени - важная компонента МНКВ, позволяющая вычислять физически реализуемое фундаментальное решение для данной области (блок или слой) однородной среды с учетом каскадной дифракции на границе сложной криволинейной формы. Описан алгоритм построения семейства зон виртуальной тени, его оптимизации и адаптации для GPU и GPU-кластера. Также описан алгоритм хранения матрицы тени в сжатом состоянии. Приведены результаты тестирования разработанных подходов.
Ключевые слова: оптимизация программ, высокопроизводительные вычисления, GPU, дифракционное моделирование.
Введение
Сейсмический метод исследования недр Земли использует отраженные волны для послойного восстановления геометрических и материальных параметров реальной геологической среды по наблюдаемым данным. Многие методы восстановления искомых параметров целевых глубоких геологических структур базируются на математическом моделировании волновых полей в покрывающей среде известного строения, которое имитирует наблюдаемые волновые поля [1-6]. В случае сложного геологического строения покрывающей среды на первый план выходят методы моделирования, использующие математический аппарат фундаментальных решений систем уравнений, описывающих механические колебания в плавно-неоднородных областях покрывающей среды. Из-за неединственности фундаментальных ре-
* Работа проводилась при частичной поддержке Шведского фонда по международному сотрудничеству в науке и высшем образовании и выполнена с использованием ресурсов информационно--вычислительного центра НГУ, а также ресурсов суперкомпьютерного комплекса МГУ им. М. В. Ломоносова.
Зятьков Н. Ю., Айзенберг А. А. Высокооптимизированная реализация вычисления матрицы тени для моделирования каскадной дифракции в геологических слоях // Вестн. Новосиб. гос. ун-та. Серия: Информационные технологии. 2016. Т. 14, № 2. С. 17-37.
ISSN 1818-7900. Вестник НГУ. Серия: Информационные технологии. 2016. Том 14, № 2 © Н. Ю. Зятьков, А. А. Айзенберг, 2016
шений проблема поиска решения, которое дает разумный баланс между вычислительной скоростью и аналитической точностью, остается актуальной проблемой математической волновой теории [1; 7; 8].
Недавно было показано [9], что физически реализуемое фундаментальное решение в плавно-неоднородной области может быть представлено в виде суммы расходящегося фундаментального решения для безграничной среды и каскадной дифракции, которая представляется в виде бесконечного ряда Неймана из дифракционных поправок нарастающего порядка. Дифракционная поправка каждого порядка представляет собой действие композитного оператора распространения-поглощения на дифракционную поправку предыдущего порядка. В ядре оператора поглощения встроена целочисленная функция виртуальной тени, принимающая либо значение 0, либо 1, которая для заданной криволинейной границы слоя определяет семейство освещенных и затененных точек границы относительно каждой точки этой границы. Оператор поглощения, за счет встроенной в него функции тени, распространяет только нефизическое волновое поле между зонами геометрических теней границы. Это позволяет моделировать каскадную дифракцию, возникающую на выпукло-вогнутых частях границ слоя.
В данной работе речь пойдет о высокооптимизированной реализации вычисления дискретного аналога функции виртуальной тени - матрицы тени. Эта процедура является важной составляющей разрабатываемого программного комплекса МНКВ (метод наложения концевых волн) для вычисления дифракционных полей внутри трехмерных геологических слоев со сложными границами [7; 8]. Работа состоит из введения, пяти содержательных разделов и заключения. В разделе 1 дано строгое определение функции виртуальной тени и предложен алгоритм вычисления ее дискретного аналога. В разделе 2 приведен алгоритм хранения матрицы виртуальной тени в памяти ЭВМ. В разделе 3 описана оптимизация вычисления матрицы виртуальной тени. В разделе 4 показаны алгоритмы адаптации вычисления матрицы тени для GPU и GPU-кластера. В разделе 5 приведены результаты тестирования разработанных процедур. В заключение приведены оценки разработанных алгоритмов и дальнейшие перспективы их развития.
1. Определение и алгоритм вычисления матрицы виртуальной тени
Рассмотрим рис. 1, на котором изображена однородная область D в форме безграничного полупространства, ограниченного криволинейной границей S и бесконечно удаленной сферической поверхностью S+м . Для границы S определим двумерную функцию виртуальной тени h(s, s') следующей формулой:
Здесь Ь(- отрезок, который соединяет две внутренние точки 5 и s', предельные к границе 5". Формуле (1) можно придать следующую «оптическую» интерпретацию, если полагать отрезок Ь( 5,5') оптическим лучом. Если из точки 5 «оптически видна» точка я', то отрезок
Ь(5, лежит полностью в области В. Тогда отрезок Ь(5, 5') будем считать физически реализуемым лучом и к(= 0 . Если из точки 5 «оптически не видна» точка я', то отрезок Ь(5, лежит не полностью в области В. Тогда отрезок Ь(5, 5') будем считать физически нереализуемым лучом и к (= 1. Описанная интерпретация соответствует рис. 1.
Используя описанную функцию, можно определять так называемые освещенные и затененные зоны на границе для моделирования каскадной дифракции внутри заданного полупространства, блока или слоя.
S (s,s')n S = 0, L (s, s'jnS
(1)
Рис. 1. Определение функции виртуальной тени для криволинейной границы: случай «оптической видимости» точки s' точкой ^ и случай «оптической невидимости» точки 5' точкой 5
Если границу дискретизировать на N треугольников, то функцию я') можно представить в виде квадратной матрицы [Нк, J размерности N х N . Отметим основные свойства матрицы виртуальной тени [Нк, J . Во-первых, матрица Нк, J - симметричная относительно
главной диагонали, что следует из определения функции виртуальной тени (1). Так как транспонирование симметричной матрицы не меняет ее, то выполняется равенство
Нц = Н9к (2)
для любого номера столбца к и строки , . Симметрия матрицы виртуальной тени Нк, J позволяет вычислять лишь ее верхнетреугольную часть, а ее нижнетреугольную часть получать из вычисленной верхнетреугольной по формуле (2). Отметим, что матрица [ Нк, J заполняется
только нулями и единицами.
Рассмотрим возможные ситуации при вычислении отдельного элемента матрицы виртуальной тени [. Выделим пару произвольных малых треугольных элементов АБк и АБ,. Эта
пара соответствует элементу матрицы Нк, с двойным индексом кц . Пусть на элементе границы АБк задана внутренняя нормаль пк и на элементе границы АБ, задана внутренняя нормаль п. Соединим центр элемента границы А5к и центр элемента границы отрезком Ь (5к, ) . Для пары элементов границы АБк и в зависимости от расположения их
нормалей относительно луча, соединяющего центральные точки элементов, существуют четыре различные ситуации:
1) если углы между обеими нормалями пк и п, и отрезком Ь (5к, ) являются тупыми,
то элементы А5к и «не видят» друг друга и скалярный элемент матрицы виртуальной тени Нк, = 1 (рис. 2, а);
2) если угол между нормалью пк и отрезком Ь (5к, ) является тупым, а между другой нормалью п, и отрезком Ь (5к, ) острым, то элементы АБк и «не видят» друг друга и скалярный элемент матрицы виртуальной тени Нк, = 1 (рис. 2, б);
3) если угол между нормалью пк и отрезком Ь (5к, ) является острым, а между другой нормалью п, и отрезком Ь (5к, ) тупым, то элементы АБк и «не видят» друг друга и скалярный элемент матрицы виртуальной тени Нк, = 1 (рис. 2, в);
б
а
в
Рис. 2. Элементы А5к и А5д «не видят» друг друга
г
д
4) если углы между обеими нормалями пк и пд и отрезком Ь (5к, ^ ) являются острыми, то элементы А5к и А5д могут как «видеть», так и «не видеть» друг друга, и значение скалярного элемента матрицы виртуальной тени Икд зависит от того, есть пересечение отрезка Ь (5к, 5д) с хотя бы одним третьим элементом границы А5п или нет такого пересечения.
Поэтому в ситуации 4) возможны два варианта:
4а) если отрезок Ь (8к , 8д) между элементами А5к и А5д пересекает какой-то промежуточный элемент АБп, то элементы границы «не видят» друг друга, тогда скалярный элемент матрицы виртуальной тени Икд = 1 (рис. 2, г);
4б) если отрезок Ь (5к, 5д) между элементами А5к и А5д не пересекает ни один промежуточный элемент АБп, то элементы границы «видят» друг друга, тогда скалярный элемент матрицы виртуальной тени Икд = 0 (рис. 2, д).
Используя приведенный выше логический анализ, опишем возможный алгоритм вычисления матрицы тени:
// Изначально массив [Икд J должен быть заполнен нулями, т. е.
И кд = 0, к, д е {1, 2,..., К].
// Двойной цикл по элементам границы
// (заполняем только верхнетреугольную часть матрицы;
// нижнетреугольную часть получаем из верхнетреугольной; // из условия симметрии матрицы).
Цикл к = 1, 2,..., N
{
Цикл д = 1, 2,..., N
{
1. Вычислить вектор Ь (чк, ) .
2. Вычислить косинус угла С08(пк, Ь) между нормалью элемента
А^к и вект°р°м Ь (к, Яд) .
3. Вычислить косинус угла С08(пд, - Ь) между нормалью элемента
Мд и вект°р°м - Ь (чк, ) .
Если С08(пк, Ь) < 0 или С08(Пд, — Ь) < 0 {
1. Нк9 = 1;
2. Перейти к выполнению новой итерации цикла по д .
}
Иначе
{
// Новый цикл по элементам границы, с целью определить // пересекает ли вектор Ь (чк, ) границу хотя бы в одной
// точке. Если пересекает, то элемент с номером к «не // видит» элемент с номером д и Нкд = 1, иначе к «видит»
// д и Нк9 = о.
Цикл п = 1, 2,..., N
{
1. Находим точку пересечения О вектора Ь (чк, ) с
плоскостью, которую образует элемент границы с номером п .
2. Определяем, принадлежит ли точка О элементу с номером п. Если да, то найдена точка пересечения вектора Ь (чк, ) с границей, элемент с номером к «не
видит» элемент с номером д и Нкд = 1. Переход к выполнению новой итерации цикла по д .
}
}
Отметим основные проблемы реализации приведенного алгоритма.
• Матрица виртуальной тени имеет размерность N х N, где N - количество треугольных элементов границы. К примеру, объем необходимой памяти для хранения одной такой матрицы для одной границы слоистой среды при N = 100 000 составляет 100 000 х 100 000 = 9,3 Гб, а при N = 150 000 - уже 21 Гб, что неприемлемо много.
}
}
• Из алгоритма, приведенного выше, видно, что его трудоемкость составляет о(N3) , что при больших размерностях матрицы виртуальной тени также неприемлемо.
Для нормальной работы алгоритма необходима оптимизация вычисления матрицы виртуальной тени как по производительности, так и по требуемой памяти.
2. Оптимизация хранения матрицы виртуальной тени
Перед тем как описывать алгоритм оптимизации хранения матрицы виртуальной тени, необходимо описать принцип реализации МНКВ, так как процедура вычисления матрицы тени является неотъемлемой компонентой этого метода. МНКВ для заданной слоистой модели среды, источника и приемников вычисляет компоненты волнового поля, в зависимости от заданного волнового кода [7; 8]. Эти вычисления численно реализуются в виде последовательного перемножения слоевых матриц на векторы волнового поля, в соответствии с указанным волновым кодом. Каждая граница слоистой модели разбивается на N треугольников с помощью процедуры триангуляции. Все векторы волнового поля на границах имеют размерность N . Все слоевые матрицы распространения-поглощения имеют размерность N х N . Матрица виртуальной тени также имеет размерность N х N. При увеличении параметра N время расчета алгоритма и объем требуемый памяти для его нормальной работы будут увеличиваться экспоненциально. Следует отметить, что все матрицы и векторы являются плотными. Уже при N ~ 105 объем необходимой памяти для хранения одной такой матрицы составит ~ 102 Гб , что практически не позволяет ее целиком хранить в оперативной памяти ЭВМ. Эта проблема была решена посредством разбиения исходной задачи на подзадачи и их последовательного решения. Было решено разрезать все матрицы размерности N х N на горизонтальные полосы, при этом выделять память только для одной полосы матрицы и затем последовательно обрабатывать каждую полосу.
При реализации алгоритма МНКВ оптимальным подходом является заполнение всех матриц виртуальных теней для каждой границы среды до запуска основного алгоритма и хранение их в сжатом состоянии, поскольку одна такая матрица зависит только от формы границы -контакта двух сред. Также было бы идеальным хранить матрицу виртуальной тени для данной границы в виде набора сжатых по определенному алгоритму полос, на которые бы разрезалась матрица.
Так как матрица виртуальной тени симметричная, то запускать алгоритм для заполнения ее нижнетреугольной части не оптимально. Выгодно нижнетреугольную часть заполнять из вычисленной ранее верхнетреугольной части. Но так как память будет выделяться только для одной полосы матрицы и после ее заполнения эта полоса будет сжиматься и храниться в сжатом состоянии, то возникает проблема выбора способа заполнения части полосы, принадлежащей нижнему треугольнику матрицы (рис. 3). Для оптимального выбора необходимо знать информацию о столбцах всех полос, заполненных на предыдущих итерациях, которые хранятся в сжатом состоянии.
НСотр[ 0]
НСотр[ 1] ->
НСотр[гагЬ — 1]
N
N
Рис. 3. Проблема, возникающая при заполнении матрицы тени, в случае ее хранения в виде набора сжатых полос
Поэтому, чтобы заполнить часть полосы из нижнего треугольника матрицы, необходимо разжимать все полосы, заполненные на предыдущих шагах и извлекать из них соответствующие столбцы. Из рассуждений нетрудно понять, что этот подход не оптимален.
Было решено разрезать матрицу виртуальной тени на полосы, а каждую полосу на квадраты (рис. 4). Все квадраты, принадлежащие верхнему треугольнику матрицы, заполняются по алгоритму расчета матрицы виртуальной тени, который приведен выше, а все квадраты, принадлежащие нижнему треугольнику матрицы, заполняются путем транспонирования соответствующих квадратов из верхнего треугольника.
НСотр[ 0]
НСотр[ 1] ->
НСотр[гагЬ — 1] ->
N
Рис. 4. Алгоритм хранения матрицы тени
Далее приведен модифицированный алгоритм расчета и хранения матрицы тени. // Цикл по квадратным подматрицам из верхнетреугольной части // матрицы тени
Цикл (1 = 1, 2,..., гахЬ) {
Цикл (у = 1,1 +1,..., гахЬ ) {
// Если квадратная подматрица диагональная.
Если(1 == у)
{
1. Заполнить верхнетреугольную часть квадратной подматрицы по алгоритму расчета матрицы тени.
2. Получить нижнетреугольную часть квадратной подматрицы из верхнетреугольной, используя свойство симметричности матрицы тени.
3. Сжать заполненную квадратную подматрицу и добавить в полосу 1 (рис. 5, а).
}
иначе
{
1. Заполнить квадратную подматрицу (1, у) по алгоритму расчета матрицы тени.
2. Заполнить квадратную подматрицу (у, 1) путем транспонирования подматрицы (1, у).
3. Сжать заполненные подматрицы и добавить в полосы 1 и у соответственно (рис. 5, б).
}
}
}
ш
НСотр[0] - —>
НСотр[ 1] -s \
\
HComp[razb — 1] - \
N
НСотр[0] НСотр[1]
HComp[razb — 1]
□
X
N
N N
Рис. 5. Алгоритм заполнения матрицы виртуальной тени: а - случай диагональной подматрицы; б - случай не диагональной подматрицы
Сжатие подматриц выполнялось в следующем виде. Считалось количество нулей матрицы, идущих подряд до первой единицы, и это количество сохранялось в виде целого числа. Затем считалось количество подряд идущих единиц до первого нуля и сохранялось в виде целого числа, и т. д. (рис. 6).
8 10 3
Рис. 6. Алгоритм сжатия матрицы виртуальной тени
б
а
Этот алгоритм сжатия является вариацией известного алгоритма Run-length encoding (RLE) [10]. Таким способом удалось достигнуть сжатия матрицы виртуальной тени более чем в 100 раз.
3. Оптимизация вычисления матрицы виртуальной тени
Обратимся к алгоритму вычисления матрицы виртуальной тени, описанному в разделе 1. Трудоемкость данного алгоритма, как уже было сказано, составляет O(N3). Реализация алгоритма представляет собой двойной цикл по парам элементов границы и третьего внутреннего цикла опять же по элементам границы. Причем количество запусков третьего цикла зависит от формы границы, для которой вычисляется матрица виртуальной тени. Чем больше пар элементов границы окажется такими, что оба угла, образующиеся линией соединяющей центры элементов и нормалями элементов, будут являться острыми, тем чаще будет запускаться внутренний цикл и тем более трудоемким станет алгоритм вычисления. К примеру, для любых вогнутых в среду границ внутренний цикл не запустится ни одного раза, так как нормали каждой пары элементов не будут направлены друг на друга и производительность алгоритма для таких типов границ окажется максимальной (рис. 7, а), а трудоемкость составит O(N2).
Матрицей тени в данном случае окажется единичная матрица, диагональные элементы которой равны нулю (каждый элемент сам себя «видит», так как Hkk = 0, Vk = 1 ... N по определению).
Другой предельный случай алгоритма вычисления матрицы виртуальной тени - если нормали всех пар элементов границы направлены друг на друга. В качестве подобной границы могут выступать любые выпуклые границы (рис. 7, б).
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
N
б
Рис. 7. Пример предельного случая типа границы и ее матрицы виртуальной тени, для которой производительность вычисления матрицы: а - максимальная; б - минимальная
Производительность алгоритма для таких типов границ окажется минимальной, а трудоёмкость составит O(N3) . Матрица виртуальной тени в данном случае окажется нулевой (зон геометрической тени для таких границ нет).
Видно, что при такой высокой трудоёмкости 0(Ы2) и тем более O(N3) для приемлемого по времени вычисления матрицы виртуальной тени для заданной границы требуются большие вычислительные ресурсы. Попробуем оптимизировать алгоритм, описанный в разделе 1. Больше всего нас будет интересовать третий внутренний цикл построения теневой матрицы, так как основные трудоёмкие вычисления производятся именно в нём. В этом цикле для заданной пары элементов проверяется, не пересекает ли отрезок, соединяющий центры этих элементов, границу хотя бы в одной точке. Если такая точка найдётся, то заданная пара элементов друг друга «не видит» и матрица виртуальной тени для этой пары заполняется 1. Заметим, что реализация третьего внутреннего цикла частично пересекается с известным алгоритмом трассировки лучей [11]. В алгоритме трассировки лучей, как и в предлагаемом алгоритме построения участков тени, самой трудоёмкой частью является поиск пересечения луча с примитивами 3Б сцены. В нашем случае в качестве примитивов сцены выступает граница, состоящая из N треугольников. Было принято решение оптимизировать работу внутреннего цикла вычисления матрицы виртуальной тени по аналогии С уже исследованной задачей.
Проверка пересечения луча с треугольником. В качестве алгоритма проверки пересечения луча с треугольником рассматривалось два варианта. Первый вариант - это так называемый «барицентрический тест» [12]. Суть его состоит в том, что, имея три точки на плоскости, можно выразить любую другую точку через ее три барицентрические координаты. Если каждая из этих координат будет больше или равна нулю, то искомая точка принадлежит треугольнику. В противном случае - не принадлежит. По определению, каждая барицентрическая координата
представляет собой отношение площади малых треугольников, на которые делит искомый треугольник рассматриваемая точка, к площади этого треугольника (рис. 8, a).
Второй вариант проверки, лежит ли точка внутри треугольника, - это так называемый «единичный тест Вупа» (англ. Woop's unit test) [13]. Основная идея данного алгоритма заключается в том, чтобы вычислить такую матрицу аффинных преобразований рассматриваемого треугольника, что после применения этого преобразования (поворот и параллельный перенос в пространстве) он преобразуется в некоторый единичный треугольник с вершинами (1,0,0); (0,1,0); (0,0,0) и нормалью (0,0,1) (рис. 8, б).
В
А
Р
а
б
Рис. 8. Алгоритм проверки пересечения луча с треугольником: а - «барицентрический тест»; б - «единичный тест Вупа»
После вычисления матрицы преобразования рассматриваемый луч преобразовывается по этим правилам, и далее делается проверка: принадлежит ли точка пересечения луча и единичного треугольника этому треугольнику. Для единичного треугольника сделать подобную проверку проще, чем для произвольного треугольника.
Проанализировав 2 описанных алгоритма на предмет производительности для задачи определения освещенных и затененных участков триангулированной границы, было решено остановиться на первом способе. В случае реализации алгоритма для GPU этот способ дал лучшие результаты, по сравнению со вторым.
Исследование ускоряющих структур для метода трассировки лучей. В алгоритме трассировки лучей, с целью ускорения поиска точек пересечения луча с примитивами сцены, как правило, используются ускоряющие структуры. Это нужно для того, чтобы не перебирать все примитивы, а делать проверку только в некоторой окрестности луча. В методе трассировки лучей выделяют в основном два типа ускоряющих структур:
1) разбиение сцены на регулярную сетку;
2) иерархические представления, разбиение сцены на нерегулярную сетку.
Смысл обеих структур в том, что 3D сцена разбивается на равномерную или неравномерную сетку, и затем рассматриваются только те кубы сетки (воксели), через которые прошел луч. Если все примитивы сцены имеют примерно одинаковый размер, то логично разбивать
трехмерное пространство на равномерную сетку. Это можно сделать трехмерным алгоритмом Брезенхема [14] или алгоритмом Фуджимото [15] (рис. 9, а).
а
г - \ —
ч ч - ^ >
S ч ч >
б
Рис. 9. Процесс «траверса» луча в иерархической сетке. Красным отмечены примитивы (треугольники), которые проверяются на предмет пересечения их с лучом
В случае когда примитивы сцены друг с другом имеют большую разницу в размерах (например, сцена «чайник на стадионе»), равномерное разбиение сцены работает неэффективно, так как один примитив (треугольник) может попасть стразу в несколько вокселей и алгоритм будет вынужден проверять пересечение его с лучом несколько раз. С целью уменьшить подобные повторные вычисления, используются иерархические представления или неравномерные сетки. Одним из таких алгоритмов разбиения является структура бинарного пространственного разбиения, называемая kd-дерево [16]. Структура представляет собой бинарное дерево ограничивающих примитивы сцены параллелепипедов или боксов, вложенных друг в друга. Каждый такой бокс может быть разбит плоскостью на два дочерних бокса. Вся сцена целиком содержится внутри корневого бокса, но, продолжая рекурсивное разбиение сцены на боксы, можно прийти к тому, что листовые узлы такого бинарного дерева будут содержать лишь небольшое количество примитивов. Этот подход позволяет использовать алгоритм бинарного поиска для нахождения примитивов, пресекаемых лучом, и решает проблему разномасштабных примитивов сцены (рис. 9, б).
Несмотря на то, что 3D сцену (в нашем случае триангулированную границу) можно разбить на равномерную или неравномерную сетку и хранить ее в таком виде до запуска основного алгоритма, существует следующая проблема описанных ускоряющих структур. В ходе выполнения поиска пересечения луча с треугольником приходится выполнять «траверс» или по-другому - определять воксели, через которые прошел луч. И на это также требуются время и ресурсы вычислительной машины.
Изучив ускоряющие структуры, используемые при реализации метода трассировки лучей, было решено применить подход, описанный в следующем разделе.
4. Адаптация вычисления матрицы виртуальной тени на GPU и GPU-кластер
Адаптация вычисления матрицы виртуальной тени на GPU. Авторы исходили из того, что вычисление матрицы виртуальной тени будет происходить с использованием GPU-ускорителей, матрица будет считаться один раз до запуска основного алгоритма МНКВ и результат будет сохраняться на жесткий диск, который затем можно переиспользовать. Как известно, на GPU часто выгодны избыточные операции, если они позволяют выполнить меньше запросов в память и ядро (функция, исполняемая GPU) имеет достаточно простую структу-
ру. Было принято решение реализовать третий внутренний цикл алгоритма построения теневой матрицы методом так называемой «грубой силы» на GPU, тем самым избегая «траверса», ветвлений и возможных сложностей, возникающих при программировании подобной задачи с использованием графических ускорителей. При этом авторы стремились делать трудоемкую проверку пересечения луча с треугольником как можно реже. Для этого каждый отрезок, соединяющий два элемента, которые «смотрят» друг на друга, обрамлялся некоторым замкнутым боксом (рис. 10).
Каждой GPU-нити назначался свой элемент триангулированной границы. Если этот элемент не лежал хотя бы частично внутри бокса, то дальнейшая его обработка GPU-нитью прекращалась. Если элемент полностью или частично лежал внутри бокса, то вычислялся тест, принадлежит ли точка пересечения луча с плоскостью треугольника этому треугольнику. Таким образом, количество операций поиска пересечения луча с треугольниками было значительно уменьшено.
Рис. 10. Схема ускоряющей структуры для алгоритма построения матрицы виртуальной тени.
Красным отмечены примитивы (треугольники), которые проверяются на предмет пересечения их с лучом
В листинге 1 показано ядро, которое реализует на языке CUDA C третий внутренний цикл алгоритма поиска освещенных и затененных зон на GPU. В этом ядре происходит проверка, лежит ли треугольник внутри бокса. Если результат отрицателен, то запускается алгоритм «барицентрический тест» проверки, пересекает ли луч треугольник. Нужно отметить, что после того, как найден хотя бы один элемент, с которым пересекается луч, дальнейшие вычисления можно не продолжать.
Листинг 1. Ядро на языке CUDA C, реализующее барицентрический тест для определения затененных и освещенных зон.
_global_ void ShadowMatrixLoop_kernel (Triangle *tr_g, int N, Coord
a, Coord b,
Coord r_dir, Coord r_o, char *h_g, int idx)
{
int n = threadIdx.x + blockIdx.x * blockDim.x;
Coord p, t, q; double bary_x, bary_y; double det; Coord e1, e2;
if (n < N) {
// если треугольник не попал в обрамляющий бокс, то прекращаем
// дальнейшую проверку
if (tr_g[n].center.x < a.x || tr_g[n].center.x > b.x) return;
if (tr_g[n].center.y < a.y || tr_g[n].center.y > b.y)
и e2
Равна
та. Равна
нию e1 и p
угольник
return;
if (tr_g[n].center.z < a.z || tr_g[n].center.z > b.z) return;
// el = v1 - v0 - вектор-сторона треугольника el.x = tr_g[n].p2.x - tr_g[n].p1.x; el.y = tr_g[n].p2.y - tr_g[n].pl.y; el.z = tr_g[n].p2.z - tr_g[n].p1.z;
// e2 = v2 - v0 - вектор-сторона треугольника e2.x = tr_g[n].p3.x - tr_g[n].pl.x; e2.y = tr_g[n].p3.y - tr_g[n].pl.y; e2.z = tr_g[n].p3.z - tr_g[n].pl.z;
// t = r_o - vo - ray origine минус вершина треугольника v0 t.x = r_o.x - tr_g[n].pl.x; t.y = r_o.y - tr_g[n].pl.y; t.z = r_o.z - tr_g[n].pl.z;
// p = cross(r_dir, e2) - векторное произведение ray direction
p.x = r_dir.y * e2.z - r_dir.z * e2.y; p.y = r_dir.z * e2.x - r_dir.x * e2.z; p.z = r_dir.x * e2.y - r_dir.y * e2.x;
// q = cross(t, el) - векторное произведение t и el q.x = t.y * el.z - t.z * el.y; q.y = t.z * el.x - t.x * el.z; q.z = t.x * el.y - t.y * el.x;
// bary_x = dot(t, p) - первая барицентрическая координата.
// скалярному произведению t и p bary_x = t.x*p.x + t.y*p.y + t.z*p.z;
// bary_y = dot(r_dir, q) - вторая барицентрическая координа-
// скалярному произведению ray direction и q bary_y = r_dir.x*q.x + r_dir.y*q.y + r_dir.z*q.z;
// det = dot(el, p) - детерминант. Равен скалярному произведе-
det = el.x*p.x + el.y*p.y + el.z*p.z;
// барицентрический тест - проверка, пересекает ли луч тре-
if ( bary_x < 0 || bary_y < 0 || (bary_x + bary_y > det) ) return;
// если пересекает, то элемент матрицы виртуальной тени равен h_g[idx] = l;
l
}
}
Адаптация вычисления матрицы виртуальной тени на GPU-кластер. В разделе 2 была показана оптимизация хранения теневой матрицы. Напомним, что в ходе оптимизации эта матрица разбивается на квадратные подматрицы. Каждая подматрица может вычисляться независимо от других подматриц. Поэтому, располагая GPU-кластером (набором GPUs), было решено воспользоваться именно этим свойством: каждому GPU назначался свой набор подматриц, которые он заполнял, используя описанный выше параллельный алгоритм построения матрицы виртуальной тени. На рис. 11 схематично показано разбиение матрицы виртуальной тени на квадратные подматрицы и алгоритм назначения подматриц каждому GPU, на примере использования трех графических ускорителей. Видно, что вычисления проводятся только для квадратов из верхнетреугольной части матрицы.
Рис. 11. Схема распределения вычислений матрицы виртуальной тени на примере трех GPU-ускорителей
Схема реализации заполнения матрицы виртуальной тени для GPU-кластера выглядит следующим образом.
1. Если заданному GPU назначена диагональная подматрица, то:
а) ее верхнетреугольная часть заполняется алгоритмом построения матрицы виртуальной тени, а нижнетреугольная получится из верхнетреугольной из свойства симметрии матрицы виртуальной тени;
б) подматрица сжимается по алгоритму сжатия матрицы виртуальной тени;
в) сжатая подматрица отправляется процессу с рангом 0.
2. Если заданному GPU назначена не диагональная подматрица, то:
а) подматрица вычисляется по алгоритму построения матрицы виртуальной тени;
б) подматрица транспонируется, результат сохраняется в виде новой матрицы;
в) обе подматрицы сжимаются по алгоритму сжатия матрицы виртуальной тени;
г) обе сжатые подматрицы отправляются процессу с рангом 0.
3. Для процесса, имеющего ранг 0:
а) занимается аналогичными вычислениями, как и остальные процессы, но при этом является приемником сжатых подматриц, получаемых от остальных процессов;
б) формирует сжатую результирующую матрицу виртуальной тени заданной границы в виде трехмерного массива HComp[k][i][ j], где к - номер полосы разбиения матрицы, i -номер квадратной подматрицы внутри заданной полосы, j - номер элемента подматрицы HComp[k][i] ;
в) сжатая матрица HComp сохраняется на жесткий диск.
Алгоритм построения матрицы виртуальной тени был реализован на GPU-кластере с использованием связки MPI и CUDA.
5. Тестирование алгоритма вычисления матрицы виртуальной тени
Первый тест вычисления матрицы виртуальной тени был запущен для границы с небольшим числом элементов для анализа и демонстрации качества вычислений. В качестве тестовой
_ 2 _ 2
границы была взята аналитическая функция /(х, у) = е_х е~у , х е[_2; 2], у е [_2; 2], которая показана на рис. 12. Данная граница была дискретизирована на 2 048 треугольников.
Рис. 12. Вид поверхности, выбранной для тестирования алгоритма построения матрицы виртуальной тени
0000000000000000000000000000000000000000000000000000000000000000 0000000000000000000000000000000000000000000000000000000000000000 0000000000000000000000000000000000000000000000000000000000000000 1111011101000000000000000000000000000000000000000000000000000000 1111111111111111110101000000000000000000000000000000000000000000 1111111111111111111111110100000000000000000001010111111111111101 1111111111111111111111111111000000000001111111111111111111111111 1111111111111111111111111111111100011111111111111111111111111111 111111111111111111111111111111110111111111111111111111111111111 1111111111111111111111111111111100111111111111111111111111111111 1111111111111111111111111111111100111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111110 1111111111111111111111111111111111111111111111111111111111111100 1111111111111111111111111111111111111111111111111111111111111100 1111111111111111111111111111111111111111111111111111111111111000 1111111111111111111111111111111111111111111111111111111111111000 1111111111111111111111111111111111111111111111111111111111110000 1111111111111111111111111111111111111111111111111111111111100000 1111111111111111111111111111111111111111111111111111111111000000 1111111111111111111111111111111111111111110000111111111110000000 1111111111111111111111111111111111111111000000001111111100000000 1111111111111111111111111111111111111111000000000011110000000000 11111111111111111111111111111111111111110000000000#0 00000 000000 1111111111111111111111111111111111111111110000000111000000000000 1111111111111111111111111111111111111111111111111111100000000000 1111111111111111111111111111111111111111111111111111100000000000 1111111111111111111111111111111111111111111111111111110000000000 1111111111111111111111111111111111111111111111111111110000000000 1111111111111111111111111111111111111111111111111111110000000000 1111111111111111111111111111111111111111111111111111110000000000 1111111111111111111111111111111111111111111111111111100000000000 1111111111111111111111111111111111111111111111111111100000000000 1111111111111111111111111111111111111111111111111111100000000000 1111111111111111111111111111111111111111111111111111000000000000 1111111111111111111111111111111111111111111111111111000000000000 1111111111111111111111111111111111111111111111111111000000000000
1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111100000111111111111111111111 1111111111111111111111111111111111000000000000111111111111111111 1111111111111111111111111111111100000000000000001111111111111111 1111111111111111111111111111111100000000000000000111111111111111 1111111111111111111111111111111100000000000000000011111111111111 1111111111111111111111111111111100000000000000000001111111111111 1111111111111111111111111111111100000000000000000001111111111110 1111111111111111111111111111111101000000000000000001111111100000 1111111111111111111111111111111111010000000000000001111001100000 1111111111111111111111111111111111111100000000000000000000000000 111111111111111111111111111111111111111111010110300000000000000 1111111111111111111111111111111111111111111111000000000000000000 1111111111111111111111111111111111111111111111000000000000000000 1111111111111111111111111111111111111111111000000000000000000000 1111111111111111111111111111111111111111101000000000000000000000 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111ф1111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111
Рис. 13. Результат работы заполнения матрицы виртуальной тени
на примере ее четырех векторов
На рис. 13 показаны 4 вектора матрицы виртуальной тени, как результат тестирования реализованного алгоритма.
Каждый из векторов, изображенных на рис. 13, можно интерпретировать как проекцию функции /(х, у) на плоскость ху . Красными точками отмечены точки на границе, с которых ведется наблюдение на всю поверхность. Нулями отмечены освещенные, относительно этих
Рис. 14. Границы типа «диск» (а), «катушка» (в), «вулкан» (Э). Три вектора матрицы виртуальной тени для этих границ (б, г, е соответсвенно)
точек, зоны, а единицами - затененные. Последний вектор изображает ситуацию, когда ведется наблюдение с точки максимума функции / (х, у) на всю поверхность. В этом случае, при данной границе, указанная точка не может «видеть» ни одну другую точку на поверхности, кроме себя. Визуальная оценка результатов вычисления матрицы виртуальной тени для приведенной границы указывает на достаточную точность реализации процедуры построения матрицы виртуальной тени.Следующие три теста были запущены для искривленных границ с реалистичным числом покрывающих элементов и были нацелены на анализ производительности вычислений и степени сжатия матрицы виртуальной тени. Тестовые границы представлены на рис. 14, а, в, д).
Каждая из этих границ образуется в результате различных способов поворота двупарабо-лической кривой, показанной на рис. 15, в трехмерном пространстве. Эта кривая представляется следующей аналитической формулой:
х(г + 0,5) = 4-4 ^/(7+0,5)^-0,5
(3)
Граница на рис. 14, а представляется в виде формулы
:(у, г + 0,5) = 1 + 13 - 4 г + 0,5)2 - 0,5
- у2 , -1,7 <г < + 0,7,-1,0 < у <+1,0
и получена путем поворота кривой (3) вокруг вертикальной оси х = 1, у = 0 . Граница на рис. 14, в представляется в виде формулы
:(у, г + 0,5) = 6 + I2 + 4 Гд/(г + 0,5)2 - 0,5
- у2 , -1,7 <г <+0,7,-1,0 < у <+1,0
и получена путем поворота кривой (3) вокруг вертикальной оси х = 6, у = 0 . Граница на рис. 14, д представляется в виде формулы
х(у, г) = 4 - 4\д/(г + 0,5)2 + у2 - 0,5
-1,7 <г <+ 0,7,-1,0 <у <+1,0
и получена путем поворота кривой (3) вокруг горизонтальной оси г = -0,5, у = 0,0 .
На рис. 14, б, г, е показано по 3 вектора матриц виртуальных теней для соответствующих границ. Каждый из векторов, изображенных на рисунке, можно интерпретировать как проекцию функции х (у, г) на плоскость уг . Красными точками отмечены точки на границе, с которых ведется наблюдение на всю поверхность. Данные точки были выбраны на кривой, которая образуется в результате пересечения тестируемых границ плоскостью у = 0 . Положение точек схематично указано на рис. 15. Точка 1 имеет координаты (х = 3,0; у = 0,0; г = 0,5), точка 2 имеет координаты (х = 3,84; у = 0,0; г = -0,2), точка 3 имеет координаты (х = 3,75; у = 0,0; г = -0,75). Более темные пятна соответствуют освещенным, относительно этих точек, зонам, а более светлым - затененные.
2
-3
012345678 X (Кт)
Рис. 15. Разрез границ, изображенных на рис. 14, а, в, д, плоскостью у = 0
В таблице указаны параметры границ и результаты тестирования высокооптимизирован-ного алгоритма построения матриц виртуальных теней, реализованного для ОРИ-кластера. Показана зависимость времени вычисления матриц от количества задействованных ОРШ. Также показана степень сжатия данных матриц. По таблице видно, что матрицы виртуальных теней удается сжимать более чем в 100 раз, а в некоторых случаях и более чем в 500 раз. Объемы сжатых матриц составляют порядка сотни Мб, что является приемлемым результатом для использования данной процедуры в реалистических задачах дифракционного моделирования.
Результаты тестирования алгоритма построения матрицы виртуальной тени
для разных типов границ
Тип границы Количество элементов Размер матрицы до сжатия, Гб Размер матрицы после сжатия, Мб Коэффициент сжатия Время счета (мин)
2 GPUs 4 GPUs 6 GPUs 12 GPUs
Диск 240 000 53,6 Гб 132 Мб 416 30 15 10 5
Катушка 240 000 53,6 Гб 419 Мб 131 28 15 10 5
Вулкан 252 000 59,1 Гб 90 Мб 673 31,3 15,8 10,5 5,5
Заключение
В данной работе были представлены алгоритмы вычисления, оптимизации хранения и вычисления, а также адаптации на GPU и GPU-кластер матрицы виртуальной тени. Учитывая, что алгоритм построения семейства зон виртуальной тени частично пересекается с известным алгоритмом трассировки лучей, его оптимизация вычислений была произведена по аналогии с уже исследованной задачей. Также процедура была адаптирована под GPU-кластер, что в совокупности дало существенный прирост ее производительности. За счет построенного алгоритма хранения матрицы тени авторам удалось добиться ее сжатия в 100 и более раз. Полученные результаты дают перспективу оптимально строить, хранить и использовать матрицы виртуальной тени для трехмерных границ реалистичных сред.
* * *
Авторы выражают благодарность А. М. Айзенбергу, А. А. Дучкову (Институт нефтегазовой геологии и геофизики им. А. А. Трофимука СО РАН, Новосибирск, Россия), А. А. Рома-ненко (Новосибирский государственный университет, Новосибирск, Россия), Ф. Андерссону (Университет Лунда, г. Лунд, Швеция), А. Стовасу (Норвежский университет естественных и технических наук, г. Тронхейм, Норвегия) и М. А. Айзенберг (Исследовательский Центр фирмы Статойл, г. Берген, Норвегия) за полезные обсуждения и замечания в ходе работы.
Список литературы
1. Carcione J. M., Herman G. C., ten Kroode A.P.E. Seismic modeling. Geophysics. 2002, 67, 4, P.1304-1325.
2. Chandler-Wilde S. N., Graham I. G., Langdon S., Spence E. A. Numerical-asymptotic boundary integral methods in high-frequency acoustic scattering. Acta Numerica, Cambridge Univ. Press, 2012. P.89-305.
3. Goldin S. V. Estimation of reflection coefficient under migration of converted and monotype waves // Russian Geology and Geophysics. 1992. 33, 4. P. 76-90.
4. Gray S. H. Seismic imaging // Geophysics. 2001. 66. P. 15-17.
5. Virieux J., Operto S. An overview of full-waveform inversion in exploration geophysics // Geophysics. 2009. 74, 6, WCC127-WCC152.
6. Virieux J., Calandra H., Plessix R-E. A review of the spectral, pseudo-spectral, finite -difference and finite-element modelling techniques for geophysical imaging // Geophysical Prospecting. 2011. 59. P. 794-813.
7. Aizenberg A. M., AyzenbergM. A., Klem-Musatov K. D. Seismic diffraction modeling with the tip-wave superposition method // Extended Abstracts of the 73-th EAGE Conference & Exhibition, Austria, Vienna, 23-26 May 2011, B018.
8. Aizenberg A. M., Klem-Musatov K. D. Progress in seismic diffraction theory - From edge and tip waves to multiple reflections-transmissions with diffractions // Extended Abstracts of the 72th EAGE Conference & Exhibition, Spain, Barcelona, 14-17 June 2010, G034.
9. Aizenberg A. M., Ayzenberg A. A. Feasible fundamental solution of the multiphysics wave equation in inhomogeneous domain of complex shape // Wave Motion. 2015, 53. P. 66-79.
10. James A. Storer. Data Compression: methods and theory. Computer Science Press, 1988.
11. Whitted Turner. An improved illumination model for shaded display. Communications of the ACM. 1980.Vol. 23, No. 6. P. 343-349,
12. Möller T., Trumbore B. Fast, minimum storage ray/triangle intersection // Journal of Graphics Tools (jgt) 1997. 2, 1, P. 21-28.
13. Woop S., Benthin C., Wald I. Watertight Ray/Triangle Intersection // Journal of Computer Graphics Techniques. 2013.Vol. 2, No. 1.
14. Роджерс Д. Алгоритмические основы машинной графики. М.: Мир, 1989.
15. Fujimoto A., Tanaka Takayuki, Iwata K. Tutorial: computer graphics; image synthesis / Eds. I. Joy Kenneth, W. Grant Charles, L. Max Nelson, Lansing Hatfield. New York, NY, USA: Computer Science Press, Inc., 1988. P. 148-159.
16. Bentley J. L. Multidimensional binary search trees used for associative searching. Communications of the ACM. 1975. Vol. 18, No. 9. P. 509-517.
Материал поступил в редколлегию 15.04.2016
36
H. 3mKOB, A. A. А^íзен6ерг
1 2 N. Yu. Zyatkov , A. A. Ayzenberg
1 Novosibirsk State University 2 Pirogov Str., Novosibirsk, 630090, Russian Federation
2 University of Bergen 41 Allegaten, Bergen, 78 N-5020, Norway
nikolay.zyatkov@gmail. com, Alena.Ayzenberg@uib. no
HIGHLY-OPTIMIZED REALIZATION OF SHADOW MATRIX COMPUTATION FOR CASCADE DIFFRACTION MODELING IN GEOLOGICAL LAYERS
We provide algorithms for the virtual shadow matrix highly-optimized computation and its storing. Procedure numerically implement these algorithms are part of TWSM software package (tip-wave superposition method) for diffraction modeling of wavefields in layered and block media with complex-shaped boundaries. The shadow matrix is an important procedure of TWSM and allows calculating the feasible fundamental solution for the area (block or layer) of homogeneous medium considering cascade diffraction at the boundary of the complex convex-concave form. We describe algorithms for the family of virtual shadow zones construction, their optimization and adaptation for GPU and GPU-cluster. We also provide an algorithm of the shadow matrix storage in the compressed state. Finally, we give the test results of the developed approaches.
Keywords: software optimization, high performance computing, GPU, diffraction modeling.
References
1. Carcione J. M., Herman G. C., ten Kroode A. P. E. Seismic modeling. Geophysics, 2002, 67, 4,1304-1325.
2. Chandler-Wilde S. N., Graham I. G., Langdon S., Spence E. A. Numerical-asymptotic boundary integral methods in high-frequency acoustic scattering. Acta Numerica, Cambridge University Press, 2012, 89-305.
3. Goldin S. V. Estimation of reflection coefficient under migration of converted and monotype waves. Russian Geology and Geophysics, 1992, 33, 4, 76-90.
4. Gray S. H. Seismic imaging. Geophysics, 2001, 66, 15-17.
5. Virieux J., Operto S. An overview of full-waveform inversion in exploration geophysics. Geophysics, 2009, 74, 6, WCC127-WCC152.
6. Virieux J., Calandra H., Plessix R-E. A review of the spectral, pseudo-spectral, finite-difference and finite-element modelling techniques for geophysical imaging. Geophysical Prospecting, 2011, 59, 794-813.
7. Aizenberg A. M., Ayzenberg M. A., Klem-Musatov K. D. Seismic diffraction modeling with the tip-wave superposition method. Extended Abstracts of the 73-th EAGE Conference & Exhibition, Austria, Vienna, 23-26 May 2011, B018.
8. Aizenberg A. M., Klem-Musatov K. D. Progress in seismic diffraction theory - From edge and tip waves to multiple reflections-transmissions with diffractions. Extended Abstracts of the 72-th EAGE Conference & Exhibition, Spain, Barcelona, 14-17 June 2010, G034.
9. Aizenberg A. M., Ayzenberg A. A. Feasible fundamental solution of the multiphysics wave equation in inhomogeneous domain of complex shape. Wave Motion, 2015, 53, 66-79.
10. Storer James A. Data Compression: methods and theory. Computer Science Press, 1988.
11. Whitted Turner. An improved illumination model for shaded display. Communications of the ACM. 1980. Vol. 23, No. 6, 343-349.
12. Möller T., Trumbore B. Fast, minimum storage ray/triangle intersection. Journal of graphics tools (jgt). 1997. 2, 1, 21-28.
13. Woop S., Benthin C., Wald I. Watertight Ray/Triangle Intersection. Journal of Computer Graphics Techniques. 2013. Vol. 2, No. 1.
14. Rogers David F. Procedural Elements for Computer Graphics. New York, McGraw-Hill, 1985. (In Russ.)
15. Fujimoto A., Tanaka Takayuki, Iwata K. Tutorial: computer graphics; image synthesis. (Eds.) I. Joy Kenneth, W. Grant Charles, L. Max Nelson, Lansing Hatfield. New York, NY, USA: Computer Science Press, Inc., 1988. P. 148-159.
16. Bentley J. L. Multidimensional binary search trees used for associative searching. Communications of the ACM. 1975. Vol. 18, No. 9, 509-517.
