CC BY
69
Математическое моделирование. Оптимальное управление Вестник Нижегородского университета им. Н.И. Лобачевского, 23010, № 5 (1), с. 193-200
УДК 517.968:519.612:004.021
РАСПАРАЛЛЕЛИВАНИЕ АЛГОРИТМОВ РЕШЕНИЯ ЛИНЕЙНОЙ ОБРАТНОЙ ЗАДАЧИ ГРАВИМЕТРИИ НА МВС-1000 И ГРАФИЧЕСКИХ ПРОЦЕССОРАХ *
© 2010 г. Е.Н. Акимова, Д.В. Белоусов
Институт математики и механики УрО РАН, Екатеринбург aen15@yandex.ru
Поступила в редакцию 05.05.2010
Для решения линейной обратной задачи гравиметрии о восстановлении переменной плотности в слое численно реализованы на многопроцессорном вычислительном комплексе МВС-1000 и графических процессорах (видеокартах NVIDIA) регулярные параллельные итерационные алгоритмы. Для решения модельной задачи гравиметрии и задачи с реальными данными проведено сравнение времени счета параллельных алгоритмов на видеоускорителях GeForce NVIDIA и МВС-1000 с анализом эффективности и ускорения.
Ключевые слова: параллельные алгоритмы, обратная задача гравиметрии, эффективность и ускорение, многопроцессорный комплекс, графические процессоры.
Введение
Рассматривается линейная обратная задача гравиметрии о восстановлении переменной плотности в горизонтальном или криволинейном слое. После предварительной обработки гравитационных данных по методике, предложенной П.С. Мартышко и И.Л. Пруткиным в работе [1], задача нахождения неизвестной плотности сводится к решению линейного двумерного интегрального уравнения Фредгольма первого рода [2].
Задача гравиметрии является существенно некорректной задачей, решение которой обладает сильной чувствительностью к погрешности правой части, полученной в результате измерений и предварительной обработки геофизических данных. При разработке методов решения задач используются идеи итеративной регуляризации [3].
Для решения линейной обратной задачи гравиметрии о восстановлении переменной плотности в слое численно реализованы на многопроцессорном вычислительном комплексе МВС-1000 и графических процессорах регулярные параллельные итерационные алгоритмы.
Для решения модельной задачи гравиметрии и задачи с реальными данными проведено срав-
Статья рекомендована к печати программным комитетом Международной научной конференции «Параллельные вычислительные технологии 2010» (http://agora.guru.ru/pavt).
нение времени счета параллельных алгоритмов на видеоускорителях GeForce GTX 285, GeForce GTX 260 и МВС-1000/64 с анализом эффективности и ускорения.
Параллельные алгоритмы решения линейной обратной задачи гравиметрии
о восстановлении плотности в слое
Одной из важнейших моделей строения земной коры является модель горизонтальной слоистой среды.
Рассматривается задача о нахождении переменной плотности а = а(х, у) в горизонтальном слое
П = {(x,y,z) е R3 : (x,y) е D, Hx < z < H2}, где Hj, H2 - константы, либо криволинейном слое
ni=i(x’y z Xе R 3: (x’ y)e D
(Я! (x, y)< z < H2 {x, y) по гравитационным данным, измеренным на площади D = {(x,y) е R2: a < x < b, c < y < d}
земной поверхности. Используется априорная информация об отсутствии аномалий плотности вне слоя с криволинейными границами Hj = H1(x, y) и H2 = H2 (x, y), такими, что
Hj < H2 V(x, y), и выполняется условие
H (x, y) ^ h = const. При этом предполагает-
X——±^ y—±x
ся, что распределение плотности a(x, y) внутри слоя не зависит от z (ось z направлена вниз) (рис. 1).
G X
а(х,У) " z=HJx.y)^
b d
1
[(je - [ )2 + (y - y ' )2 + H ( ', y' )]
j(x ', y ' ) ' dy ' =
(1)
ентного типа [4]: метод минимальных невязок, метод наискорейшего спуска, метод минимальной ошибки и метод простой итерации (МПИ) в виде
(4)
где Ятах - максимальное собственное значение матрицы А + аЕ (симметричный случай).
Условие останова итерационных процессов следующее:
Рис. 1. Задача о нахождении плотности в слое
Задача нахождения неизвестной плотности сводится к решению линейного двумерного интегрального уравнения Фредгольма первого рода:
[( - [2 + (У - У' )2 + я2(х'> У' )]2,
= Д^(х, у ),
где / - гравитационая постоянная, А^(х, у) -гравитационный эффект, порождаемый источниками в горизонтальном или криволинейном слое.
После дискретизации уравнения на сетке, где задана (х, у), и аппроксимации инте-
грального оператора по квадратурным формулам задача (1) сводится к решению системы линейных алгебраических уравнений (СЛАУ) либо с симметричной положительно определенной матрицей (горизонтальный слой), либо с несимметричной матрицей (криволинейный слой). Так как уравнение (1) относится к классу некорректно поставленных задач, то СЛАУ, возникающая в результате дискретизации уравнения, является плохо обусловленной и преобразуется к виду (схема Лаврентьева)
(А + а Е )г = Ь, (2)
где а - параметр регуляризации.
В случае криволинейного слоя исходная матрица СЛАУ несимметрична, поэтому эта система предварительно преобразуется к виду (схема Тихонова)
(ЛтА + а'Е^ = АтЬ , (3)
т
где А - транспонированная матрица, а' - параметр регуляризации.
Для решения уравнений (2) и (3) используются регулярные итерационные методы гради-
< є.
В работах [5, 6] численно реализованы регулярные параллельные итерационные алгоритмы решения линейной обратной задачи гравиметрии (1) с помощью библиотеки MPI [7] на языке Фортран на многопроцессорном вычислительном комплексе МВС-1000 — российском массивно-параллельном суперкомпьютере кластерного типа с распределенной памятью, установленном в Институте математики и механики УрО РАН.
Алгоритмы были реализованы на следующих вычислителях:
1. МВС-1000/17ЕК (UM32), состоящем из 16 2-процессорных модулей Xeon 2.4 ГГц, интерфейса GigabitEthernet и 68 Гбайт оперативной памяти;
2. МВС-1000/64 (UM64), состоящем из 14 2-процессорных 2-ядерных модулей AMD Opte-ron 64 bit (2.6 ГГц), интерфейса GbitEthernet и 112 Гб оперативной памяти.
Распараллеливание итерационных методов градиентного типа основано на разбиении матрицы A горизонтальными полосами на m блоков, а вектора решения z и вектора правой части b СЛАУ на m частей так, что n = m*L, где n - размерность системы уравнений, m -число процессоров, L - число строк матрицы в блоке. На каждой итерации каждый из m процессоров вычисляет свою часть вектора решения. В случае умножения матрицы A на вектор z каждый из m процессоров умножает свою часть строк матрицы A на вектор z. В случае матричного умножения A A каждый из m процессоров умножает свою часть строк транспонированной матрицы A на всю матрицу A. Host-процессор отвечает за пересылки данных и также вычисляет свою часть вектора решения (рис. 2). Для метода простой итерации (4) максимальное собственное значение Xmax матрицы A+aE находится с помощью степенного метода с использованием параллельного алгоритма умножения матрицы на вектор.
а с
Рис. 2. Схема распределения данных по процессорам
В предыдущих работах для решения линейной задачи гравиметрии с различными гравитационными данными проведен анализ ускорения и эффективности параллельных алгоритмов
^ =Т 1 Тш > Ет = ^ 1 т ,
где Тт - время выполнения параллельного алгоритма на МВС-1000 с числом процессоров т (т>1), Т1 - время выполнения последовательного алгоритма на одном процессоре. Время Тт представляет собой совокупность чистого времени счета и накладных расходов на межпроцессорные обмены Тт = Тс + Т0. Число процессоров т соответствует упомянутому разбиению векторов на т частей и разбиению исходной области на т подобластей.
В общем случае эффективность распараллеливания меняется в пределах 0<Ет<1. В идеальном случае при равномерной и сбалансированной загрузке процессоров и минимальном времени обменов между ними Ет близка к единице, но при решении практических задач она уменьшается за счет накладных расходов.
При решении задачи о восстановлении плотности в слое на МВС-1000 с помощью параллельных алгоритмов матрица СЛАУ большой размерности формируется и хранится в памяти каждого процессора по частям, что (в некоторых случаях) дает эффективность распараллеливания Ет >1.
Результаты численных экспериментов на МВС-1000
Задача 1. На многопроцессорном вычислительном комплексе МВС-1000/64 решена модельная задача о восстановлении плотности в горизонтальном слое между глубинами И1=1 км и И2 =1.5 км для области £ь имеющей размеры 10^10 км . Шаги сетки:
Ах «Ду « 0.091 км. Гравитационная постоянная f = 6.67 -10-8 см3/г-с2.
После дискретизации исходного уравнения на сетке, где задана функция Дg(х, у), и аппроксимации интегрального оператора по квадратурным формулам задача (1) сводится к СЛАУ с симметричной плохо обусловленной заполненной матрицей 12100*12100.
Для решения задачи о восстановлении плотности в слое использовался параллельный итеративно регуляризованный МПИ (число итераций - 280) с параметром регуляризации а = 0.0001. При этом относительная норма невяз-
ки уменьшилась от 1.0 до
Ак - Ъ
= 3.1-10
4
На рис. 3 изображено исходное аномальное гравитационное поле для области £1, предоставленное Д.Е. Кокшаровым (ИГФ УрО РАН).
На рис. 4 изображены линии уровня и распределение плотности в слое, восстановленной по выделенному аномальному полю для области £1.
Задача 2. На многопроцессорном вычислительном комплексе МВС-1000/64 решена задача с реальными данными, предоставленными
В.А. Пьянковым (ИГФ УрО РАН), о восстановлении плотности в горизонтальном слое между глубинами Их = 10 км и Н2 = 20 км для области £2, имеющей размеры 120*220 км2. Шаги сетки: Дх да 0.6 км, Ау » 1.1 км. Гравитационная постоянная/= 6.6710"8 см3/г-с2.
После дискретизации исходного уравнения на сетке, где задана функция А§(х, у), и аппроксимации интегрального оператора по квадратурным формулам задача (1) сводится к СЛАУ с симметричной плохо обусловленной заполненной матрицей 40000*40000.
Для решения задачи использовался параллельный итеративно регуляризованный
10000+т
Рис. 4. Линии уровня и распределение восстановленной плотности в слое для области
МПИ (число итераций — 430) с параметром регуляризации а = 0.001. При этом относительная норма невязки уменьшилась от 1.0
\\Лгк - ¿11 до Л------!!« 4.10-3
И 410'
На рис. 5 изображено исходное аномальное гравитационное поле для области £2.
На рис. 6 изображены линии уровня и распределение плотности в слое, восстановленной по выделенному аномальному полю для области £2.
Результаты решения задачи 2 переданы специалистам в Институт геофизики УрО РАН для геофизической интерпретации.
В табл. 1 приведены времена счета и коэффициенты ускорения и эффективности решения задачи гравиметрии с реальными данными с использованием параллельного алгоритма МПИ (число итераций — 430) на МВ С-1000/64 для 200*200 точек сетки.
Отметим, что матрица СЛАУ размерности 40000*40000 формируется и хранится в памяти
Распараллеливание на видеоускорителях с помощью технологии CUDA
и результаты численных экспериментов
Для организации параллельных вычислений актуальным в настоящее время является использование видеоускорителей (GPU) компании NVIDIA (рис. 7) [8]. Основой распараллеливания служит архитектура графических процессоров. В видеочипах NVIDIA базовым блоком является мультипроцессор, содержащий несколько десятков ядер, несколько сотен арифметико-логических устройств (ALU), несколько тысяч регистров и небольшое количество разделяемой общей памяти. Видеоускорители NVIDIA содержат быструю глобальную память с возможностью доступа к ней всех мультипроцессоров, локальную память в каждом мультипроцессоре, а также специальную память для констант.
Работа нескольких ядер мультипроцессора основана на архитектуре типа SIMD, т.е. каждый из процессоров выполняет одну и ту же команду над разными элементами данных.
Видеочипы GPU состоят из массивов исполнительных блоков, управляющих потоками блоков, разделяемой памяти небольшого объема и контроллеров памяти на несколько каналов. Такая архитектура позволяет чипу обрабатывать нескольких тысяч потоков данных, требующих высокой пропускной способности памяти.
Для поддержки параллельных вычислений компания NVIDIA разработала технологию CUDA [9] — среду разработки программ на языке Си, позволяющую создавать программное
каждого процессора по частям.
Рис. 5. Аномальное гравитационное поле Ag (x, y) для области S2
40000 80000 120000
Рис. 6. Линии уровня и распределение восстановленной плотности в слое для области $2
Решение задачи гравиметрии 2 о восстановлении плотности в слое
Таблица 1
Число процессоров m Время мин Ускорение Sm Эффективность Em
1 55.82 — —
2 32.96 1.69 0.85
3 20.83 2.80 0.93
4 15.83 3.53 0.88
5 12.40 4.50 0.90
10 6.26 8.92 0.89
15 4.16 13.4 0.89
20 3.28 17.0 0.85
30 2.12 26.3 0.88
40 1.80 30.0 0.76
60 1.19 46.9 0.78
80 0.88 63.4 0.79
обеспечение для решения сложных вычислительных задач.
Модель программирования в CUDA основывается на группировании потоков. Потоки объединяются в блоки потоков (thread block) — одномерные или двумерные сетки потоков, взаимодействующих между собой при помощи разделяемой памяти и точек синхронизации. Программа (kernel) выполняется над сеткой (grid) блоков потоков (thread blocks) (рис. 8). Каждый блок может быть одно-, дву- или трехмерным по форме и состоять из 512 потоков на текущем аппаратном обеспечении.
Основной процесс приложения CUDA работает на универсальном компьютере (Host).
CPU-код инициализирует GPU, распределяет память на видеокарте, копирует константы в память видеокарты, запускает несколько копий процессов kernel на видеокарте, копирует полученный результат из видеопамяти, освобождает память и завершает работу.
Параллельная часть кода программы, написанной на языке CUDA, выполняется как множество
нитей (потоков). Для оптимизации работы с памятью нити группируют в блоки фиксированного размера. Блоки объединяются в группы блоков. Параллельная процедура выполняется над группой блоков. Типовая процедура (kernel) параллельной обработки массива Data выглядит следующим образом: kernel (Data) <<<blocks,
threads>>>, где blocks — количество блоков в сетке, threads — количество потоков в блоке.
Описанные выше линейные задачи гравиметрии 1 и 2 были решены на видеоускорителях GeForce GTX 285 (GPU-1) и GeForce GTX 260 (GPU-2).
Результаты решения задач представлены на рис. 4 и 6.
В табл. 2 приводятся технические характеристики используемых систем.
Отметим, что для эффективных расчетов объем памяти на Host-компьютере (CPU) должен превышать суммарный объем используемых в задаче массивов данных. Размер требуемой оперативной памяти для матрицы размерности M х N вещественных чисел одинарной
точности вычисляется по формуле: size_ t s = M X x N X sizeof (float).
Рассмотрим распараллеливание алгоритмов решения обратной задачи гравиметрии (1) о восстановлении переменной плотности в слое методом простой итерации с помощью CUDA. Базовыми операциями для распараллеливания итерационных процессов являются реализации параллельных функций над матрицами и векторами: параллельное умножение матрицы на вектор, параллельное умножение матриц и т.д.
На рис. 9 показано разбиение данных в процессе умножения матрицы на вектор.
Для оптимизации работы с памятью при вычислениях было использовано два приема.
1. Для сеток не очень большой размерности (размер СЛАУ 12100*12100), когда данные вхо-
Рис. 7. Видеоускоритель GeForce GTX 285
дят в память видеокарты, матрица A порядка N и вектор Z размерности N расширяются до размерности M и дополняются нулями таким образом, чтобы M было кратно числу блоков. Размер блока BLOCK_SIZE (threads) выбирается кратным 16, поскольку в одном блоке группируются до 512
Host Device
Grid 1
Kömel Block Slock Block
1 <0.0) (1,0) (2.0)
BlocK/' Bloch і Block
(Ш.М (1,1) 1 (2.1 і
/ / і \
/ ■ ' 1 Grid 2 і
/ І і
Kernel -■ ■ I ї і
Slock ПИ)
IhrL^U J Cfl. 01 Hir* jd IM» 1ЪГ£ ltd a.oi Ihr L'iidl .Jvftl Ihfod U.1I)
ItirL'Ull №11 llirrad (1. li Ilirv Jtl <2. 1 1 1 hrt1 Hil (І. I ) IlirTHd 14.1)
ІЬпміІ (№, 21 Itirf-ed tl.2) Ihrr-dri 42.21 Thrf-ad (j.ij llirr^d (4,2)
Рис. 8. Модель параллельной обработки данных
Таблица 2
Технические характеристики систем (GPU-1,2 и CPU)
Характеристики GPU-1 GPU-2
Количество процессорных ядер 240 192
Частота ядра (МГц) 648 576
Частота процессора (МГц) 1476 1242
Количество видеопамяти (Мб) 1024 896
Характеристики СРи Intel Core I5-750, 4-ядер.
Частота процессора (ГГц) 2.66
Оперативная память (Гб) 8
Разрядность ОС (Бит) 64
Вкчк ТІ1Г Bloch гі£р Fiki.: і. нп
4 —- №
. і 04*СИ SiZf!
і 1 . СНс-сЬ -iize г
1 |ь Block ііій f '
М
Рис. 9. Умножение матрицы на вектор
Таблица 3
Времена решения задач гравиметрии 1 и 2
Модельная задача гравиметрии 1 Метод — МПИ (280 итер.) Матрица СЛАУ 12100x12100 Относительная норма невязки: II Лг1 — ¿11 Host Intel Core I5-750 (2.66 ГГц) 84.0 сек
Host + GeForce GTX 285 (240 ядер) 14.0 сек
Host + GeForce GTX 260 (192 ядра) 19.5 сек
МВС-1000/64 (1 проц., 2.6 ГГц) 157.8 сек
МВС-1000/64 (2 проц.) 91.8 сек
И "3110' МВС-1000/64 (8 проц.) 20.4 сек
МВС-1000/64 (10 проц.) 16.8 сек
Задача гравиметрии 2 Метод — МПИ (430 итер.) Матрица СЛАУ 40000x40000 Относительная норма невязки: II Azк - ¿1 II II ~ 4 1 О-3 Host Intel Core I5-750 (2.66 ГГц) 31.46 мин
Host + GeForce GTX 285 (240 ядер) 4.08 мин
Host + GeForce GTX 260 (192 ядра) 5.28 мин
МВС-1000/64 (1 проц., 2.6 ГГц) 55.82 мин
МВС-1000/64 (2 проц.) 32.96 мин
н 410' МВС-1000/64 (4 проц.) 15.83 мин
МВС-1000/64 (5 проц.) 12.40 мин
МВС-1000/64 (10 проц.) 6.26 мин
МВС-1000/64 (15 проц.) 4.16 мин
потоков. Тогда количество блоков вычисляется по формуле: blocks = M/BLOCK_SIZE. Вычисления производятся без выгрузки данных в память Host-процессора. Данные находятся только в памяти видеокарты.
2. Для сеток довольно большой размерности (размер СЛАУ 40000x40000), когда данные не входят в память видеокарты, наилучшим по быстродействию оказывается метод вычисления элементов матрицы A «на лету», т.е. вычисление значения элемента матрицы происходит в момент обращения к этому элементу без сохранения его в память видеокарты. Это позволяет существенно снизить количество обращений к памяти видеокарты и заметно ускорить процесс вычислений по сравнению с хранением матрицы A в памяти Host-процессора и порционной загрузкой в видеоускоритель для вычислений.
В табл. 3 приводятся времена решения модельной задачи гравиметрии 1 и задачи гравиметрии с реальными данными 2 на Host-
процессоре Intel Core I5-750 без использования и с использованием видеоускорителей GeForce на разных сетках: 110x110 (матрица СЛАУ 12100x12100) и 200x200 (матрица СЛАУ 40000x40000). Приводится сравнение времени счета решения задач гравиметрии на видеокартах с временем счета задач на МВС-1000/64.
Таким образом, время решения линейной задачи гравиметрии довольно большой размерности методом МПИ на видеоускорителях GeForce сравнимо с временем решения линейной задачи гравиметрии на 15 процессорах МВС-1000/64.
Заключение
Для решения линейной обратной задачи гравиметрии о восстановлении переменной плотности в слое численно реализованы на многопроцессорном вычислительном комплексе
МВС-1000 и графических процессорах (видео-
картах NVIDIA) регулярные параллельные итерационные алгоритмы.
Для решения модельной задачи гравиметрии и задачи с реальными данными проведено сравнение времени счета параллельных алгоритмов на видеоускорителях GeForce GTX 285, GeForce GTX 260 NVIDIA и МВС-1000/64 с анализом эффективности и ускорения.
Авторы выражают признательность за поддержку и внимание к работе члену-коррес-понденту РАН В.В. Васину.
Работа выполнена при поддержке гранта РФФИ (проект № 09-01-00053), Междисциплинарного проекта УрО РАН, Программы фундаментальных исследований Президиума РАН № 14 (проект № 09-П-1-1003).
Список литературы
1. Мартышко П.С., Пруткин И. Л. Технология разделения источников гравитационного поля по глубине // Геофизический журнал. 2003. Т. 25. № 3. С. 159-168.
2. Мартышко П.С., Кокшаров Д.Е. Об определении плотности в слоистой среде по гравитационным данным // Геофизический журнал. 2005. Т. 27. № 4. С. 678-684.
3. Васин В.В., Агеев А.Л. Некорректные задачи с априорной информацией. Екатеринбург: Наука, 1993. 262 с.
4. Васин В.В., Еремин И.И. Операторы и итерационные процессы фейеровского типа. Теория и приложения. Екатеринбург: УрО РАН, 2005. 210 с.
5. Акимова Е.Н., Гемайдинов Д.В. Параллельные алгоритмы решения задачи гравиметрии о восстановлении плотности в слое // Труды Института математики и механики УрО РАН. 2007. Т. 13. № 3. С. 3-21.
6. Акимова Е.Н. Параллельные алгоритмы решения обратных задач гравиметрии и магнитометрии на МВС-1000 // Вестник ННГУ. 2009. № 4.
С. 181-189.
7. Baranov A.V., Latsis A.O., Sazhin C.V., Khramt-sov M.Yu. The MVS-1000 System User's Guide. URL: http://parallel.ru/mvs/user.html.
8. URL: http://www.nvidia.ru/
9. URL: http://www.ixbt.com/video3/cuda-1.shtml
PARALLELIZATION OF ALGORITHMS FOR SOLVING LINEAR INVERSE GRAVIMETRIC PROBLEM ON MVS-1000 AND GRAPHIC PROCESSORS
E.N. Akimova, D.V. Belousov
For solving the linear inverse gravimetric problem on restoration of a variable density in a layer, regular parallel iterative algorithms have been numerically realized on the multiprocessor system MVS-1000 and graphic processors (NVIDIA video cards). For solving the model gravimetric problem and the gravimetric problem with real gravitational data, the comparison has been carried out of the execution time for the parallel algorithms on NVIDIA video accelerators GeForce and MVS-1000 with the analysis of efficiency and acceleration.
Keywords: parallel algorithms, inverse gravimetric problem, efficiency and acceleration, multiprocessor complex, graphic processors.
