CC BY
88
Информационные технологии Вестник Нижегородского университета и м. Н.И. Лобаче вского, 2012, № 4 (1), с. 238-246
УДК 519.612
РЕШЕНИЕ ПЛОХООБУСЛОВЛЕННЫХ СИММЕТРИЧНЫХ СЛАУ ДЛЯ ЗАДАЧ СТРОИТЕЛЬНОЙ МЕХАНИКИ ПАРАЛЛЕЛЬНЫМИ ИТЕРАЦИОННЫМИ МЕТОДАМИ
© 2012 г. В.Л. Якушев 1, В.Н. Симбиркин 2, А.В. Филимонов 1,2, П.А. Новиков 1,
И.Н. Коньшин 3’4, Г.Б. Сушко 4, С.А. Харченко 3,4
1 Институт автоматизации проектирования РАН, Москва 2 ООО «Еврософт», Москва 3 Вычислительный центр им. А.А. Дородницына РАН, Москва
4 ООО «ТЕСИС», Москва
П8ступела ч редакцею 30.12.2011
Приведены данные по решению итерационными методами плохообусловленных симметричных систем линейных алгебраических уравнений (СЛАУ) для набора задач строительной механики. При построении предобусловливания использовалась приближенная треугольная факторизация второго порядка точности. Для организации распараллеливания вычислений по процессорам и ядрам использовались библиотеки MPI и TBB соответственно. Численные эксперименты проведены на суперкомпьютере «Ломоносов» в различных параллельных режимах.
Ключечые сл8ча: задачи строительной механики, система линейных алгебраических уравнений, итерационные методы, неполная треугольная факторизация, параллельные вычисления, MPI, распараллеливание по нитям, суперкомпьютер «Ломоносов».
Введение
В России развернута работа по реализации крупномасштабных строительных проектов. К их числу относится, например, строительство комплекса «Москва-Сити», высотных зданий в Москве, Санкт-Петербурге, Сочи, Новороссийске и других крупных городах, а также уникальных спортивных сооружений (Олимпийские объекты в Сочи, новый футбольный стадион на Крестовском острове в Санкт-Петербурге и т.д.). При этом отметим, что за последние годы произошел ряд аварий и катастроф с разрушением строительных объектов, например, Аквапарка и Басманного рынка в Москве, машинного зала Саяно-Шушенской ГЭС и др. Это способствовало ужесточению требований к расчетным моделям сооружений и точности их реализации численными методами.
В ходе научно-исследовательских и проектных работ возникла необходимость существенного развития вычислительных средств и информационных технологий для обеспечения надежности и экономичности строительства и эксплуатации указанных объектов. Сейчас предъявляются всё большие требования к подробности конечно-элементных моделей сооружений, и, следовательно, растет размерность моделей. Актуальной проблемой также является достижение высокой скорости решения задач.
Следует заметить, что проектирование многих крупных строительных объектов в России уже сейчас выполняется проектно-строительными фирмами Европы, США и других стран. Эти фирмы имеют соответствующий опыт, в том числе в применении программных продуктов и вычислительных систем большой мощности. К сожалению, в России разработки САЕ- и САD-систем для автоматизации инженерных расчетов и конструирования в строительстве никем не финансируются, разработка программного обеспечения ведется небольшими фирмами за счет средств, получаемых, в основном, от другой коммерческой деятельности. В дальнейшем отсутствие отечественных исследований и разработок в данном направлении приведет к тотальному господству в инженерной практике западных программных продуктов и технологий, как это уже произошло в ряде отраслей.
Разработка отечественного САЕ-продукта, в том числе для мощных ЭВМ с параллельной архитектурой, соответствует современным тенденциям развития автоматизации исследований и проектирования в строительстве. При этом могут быть более полно реализованы перспективные отечественные подходы к расчету строительных конструкций с использованием единых трехмерных моделей всех несущих конструкций сооружения и основания, учитывающих пространственный характер воздействий,
физическую и геометрическую нелинейность, технологию и последовательность возведения конструкций и другие важные факторы.
Данная статья является результатом совместной работы нескольких коллективов по созданию эффективных итерационных методов решения плохообусловленных симметричных систем линейных алгебраических уравнений (СЛАУ). Для генерации входных данных был использован программный комплекс STARK ES, разработанный в ООО «Еврософт» (Россия, Москва) и предназначенный для расчета пространственных строительных конструкций на прочность, устойчивость и колебания [1,2]. При проведении численных экспериментов и анализе эффективности использовались алгоритмы и программы, разработанные ООО «ТЕСИС» (Россия, Москва).
Моделирование задач строительной механики
Для моделирования плосконапряженно-
го/плоскодеформированного и изгибного напряженного состояния используются плоские четырехугольные и треугольные конечные элементы.
При моделировании плосконапряженно-го/плоскодеформированного состояния в каждом узле используются три степени свободы в плоскости элемента: две поступательные и одна вращательная. Используются три типа конечных элементов: две разновидности гибридных конечных элементов и один вид элементов метода перемещений. При построении гибридных элементов использован функционал Рейсснера, а элементов метода перемещений - функционал Лагранжа. Первый тип гибридного элемента получен с использованием квадратичной аппроксимации напряжений по площади элемента (для четырехугольных элементов) или линейной аппроксимации напряжений по площади элемента (для треугольных элементов) и квадратичной аппроксимации перемещений по границам элемента (схема Олмана для перемещений). Второй тип гибридного элемента получен с использованием квадратичной аппроксимации напряжений по площади элемента (для четырехугольных элементов) или линейной аппроксимации напряжений по площади элемента (для треугольных элементов) и кубической аппроксимации перемещений по границам элемента (схема Пиана). Элементы метода перемещений получены по схеме Олмана.
Для моделирования изгибного напряженного состояния также используются три типа конечных элементов. При построении гибридных
элементов использован функционал Рейсснера, а для элементов метода перемещений - функционал Лагранжа. Г ибридные элементы первого типа строятся на основе теории толстых плит Миндлина-Рейсснера и имеют в узле три степени свободы: поперечное перемещение и два угла поворота нормали к срединной плоскости. В изгибном элементе поперечное перемещение на границе элемента аппроксимируется полиномом второго порядка, а углы поворота нормали
- полиномом первого порядка. Использована квадратичная аппроксимация напряжений по площади элемента (для четырехугольных элементов) и линейная аппроксимация напряжений по площади элемента (для треугольных элементов). Эти элементы свободны от эффекта «сдвигового запирания» и могут использоваться для расчета толстых и тонких плит. Гибридные элементы второго типа получены с использованием квадратичной аппроксимации напряжений по площади элемента (для четырехугольных элементов) или линейной аппроксимации напряжений по площади элемента (для треугольных элементов) и кубической аппроксимации поперечных перемещений по границам элемента (схема Пиана). Элементы метода перемещений получены по схеме Батоша на основе теории толстых плит Миндлина-Рейсснера [3].
Для моделирования напряженного состояния массивных тел используются современные объемные конечные элементы. Используются три типа объемных конечных элементов: два типа гибридных конечных элементов и один тип элементов метода перемещений. При построении гибридных элементов используется функционал Рейсснера, а элементов метода перемещений - функционал Лагранжа. Гибридные элементы первого типа являются равновесными, а второго типа - неравновесными. Элементы метода перемещений являются изопарамет-рическими. Все типы элементов имеют в узле три степени свободы. Равновесные гибридные элементы свободны от эффекта «запирания по толщине», обладают высокой точностью по напряжениям и перемещениям и могут использоваться для расчета как толстых, так и тонких тел [3].
Для дальнейшего анализа выбраны гибридные элементы первого типа.
Опесанее расчетным М8делей. Для анализа и проведения экспериментов отобран ряд конечно-элементных моделей проектируемых строительных объектов из практики института ЦНИИСК им. В.А. Кучеренко. Основные характеристики расчетных моделей сооружений представлены в таблице 1, а на рисунке 1 приведен общий вид некоторых из них.
Таблеца 1
Описание моделей строительных конструкций
№ п/п Наименование объекта Кол-во узлов Кол-во конечных элементов Кол-во степеней свободы Кол-во ненулевых элементов в матрице жесткости
1 Оболочка покрытия 2-й сцены Мариинского театра, С.-Петербург 6 670 22 773 40 020 871 093
2 Спортивно-досуговый центр 14 797 15 801 88 782 1 189 954
3 Нефтедобывающая платформа 17 312 28 623 98 650 13 842 114
4 Горнолыжный спуск «Воробьевы Г оры», Москва 18 760 34 822 112 560 2 826 091
5 Высотное здание комплекса «Новосити», Новороссийск 32 626 33 612 195 024 2 736 424
6 22-эт. здание по ул. Параллельная, Сочи 41 257 48 267 247 542 3 494 364
7 Жилой комплекс «Королевский парк», Сочи 42 762 62 408 256 572 5 470 229
8 Модельная задача Куб 50x50x50 132 651 125 000 397 950 13 842 114
9 Высотное здание по ул. Мосфильмовской (жилой комплекс Воробьевы горы) 71 169 88 377 427 014 4 327 831
10 Высотное здание «Москва-Сити», Н=380 м 125 932 178 059 755 592 11 228 013
11 Стадион «Газпром-Арена», С.-Петербург, модель I 159 686 230 386 798 230 15 566 477
12 Большая ледовая арена, Сочи 148 315 172 285 889 890 16 035 641
13 Высотное здание 1 270 695 284 000 1 623 600 22 203 726
14 Стадион «Газпром-Арена», С.-Петербург, модель II 426 071 577 500 2 534 446 49 500 433
15 Высотное здание 2 541 295 568 000 3 247 200 44 404 241
16 Высотное здание 3 811 895 852 000 4 870 800 66 764 339
17 Высотное здание 4 1 339 015 1 403 570 8 033 520 110 021 813
модель! модель 4 модель 10
Рис. 1. Общий вид конечно-элементных моделей сооружений
Особенности СЛАУ задач строительной механики
Эффективное решение систем линейных уравнений, возникающих в задачах строительной механики, предполагает использование самой современной вычислительной техники, поскольку число неизвестных в задачах быстро растет и уже сейчас доходит до нескольких десятков миллионов. Среди других различных типов систем линейных уравнений эти задачи выделяются гигантской обусловленностью (109-1013) даже при относительно небольшой размерности, что приводит к тому, что при их решении почти всегда используются прямые методы. При попытках использования итерационных методов, как правило, возникают существенные проблемы со сходимостью алгоритмов. С другой стороны, основная общемировая тенденция развития алгоритмов решения СЛАУ состоит в использовании итерационных алгоритмов, в том числе по причине их высокой эффективности на самой современной вычислительной технике для широкого круга приложений. В данной работе сделана попытка применить существующие в ООО «ТЕСИС» разработки в области итерационных алгоритмов при решении СЛАУ задач строительной механики.
Параллельный итерационный алгоритм решения СЛАУ
Комбинированный режим параллельных вычислений. Современные суперкомпьютерные вычислительные системы, как правило, имеют неоднородную архитектуру. С одной стороны, имеется набор вычислительных узлов с распределенной памятью, обмен данными между которыми может быть осуществлен по быстрой обменной сети. С другой стороны, каждый узел представляет собой многопроцессорный/многоядерный компьютер с общим доступом к оперативной памяти. Программная реализация вычислительных алгоритмов на компьютерах подобной архитектуры предполагает [4] использование стандарта MPI при распараллеливании по распределенной памяти (по узлам вычислений), а по общей памяти узла распараллеливание по процессорам/ядрам естественно осуществлять на основе стандартов работы с потоками с встроенными механизмами динамической балансировки нагрузки, таких, как OpenMP или Intel® Threading Building Blocks (TBB). При этом предполагается (рис. 2), что на каждом узле имеется только один MPI процесс, который порождает на этом узле нужное количество одновременно работающих потоков вычислений.
Рис. 2. Комбинированная MPI+threads организация параллельных вычислений
Шаблонная реализация, типы данных, мелкоблочные вычисления. Многообразие приложений алгоритмов решения СЛАУ диктует различные форматы представления данных. В том числе, необходимо решать задачи с симметричными и несимметричными матрицами, вещественными и комплексными, с постоянным размером блока мелкоблочного биения (1, 3, 5 и т.д.), решать плохообусловленные задачи со сверхвысокой точностью и т.д. Единая реализация алгоритмов решения СЛАУ для достаточно широкого многообразия типов данных обеспечивается шаблонной реализацией всех алгоритмов в терминах шаблонов C++. В частности:
а) поддерживаются три типа данных (short, int, long long) для хранения целочисленных массивов;
б) поддерживаются два основных типа данных для хранения вещественных чисел (float, double) при хранении данных матрицы и предо-бусловливателя;
в) вещественные данные матрицы и пред-обусловливателя могут храниться в блочном формате k х k, где k - постоянный размер блока;
г) для поддержки несимметричных неполных треугольных разложений используется сдвоенный формат хранения блоков мелкоблочного биения;
д) поддерживаются 4 типа вещественных векторных данных (double, double long, quad и octal), два последних типа данных - эмуляция четверной и восьмерной точности с помощью библиотеки QD [5].
Предобусловливание на основе неполных разложений второго порядка точности. В работе [6] предложен алгоритм построения предобу-словливания второго порядка точности, для которого имеют место соотношения:
A + E = LU + L(I)U(II) + L(II)U(I), (1)
где A - матрица коэффициентов системы уравнений, LI'I) и UJ) - соответственно нижняя и верхняя треугольная части предобусловливате-
Рис. 3. Геометрическая объемная декомпозиция области Рис. 4. Геометрическая декомпозиция квазиповерхности
ля (элементы первого порядка точности), Ll'11) и и®, соответственно, нижняя и верхняя треугольная части элементов второго порядка точности предобусловливания, E - матрица ошибок. Предобусловленные итерации производятся с матрицей А^^и®)-1. В процессе неполной факторизации в памяти одновременно хранятся только все элементы матриц A, Ll'1) и U'1) и некоторая небольшая часть элементов матриц Ll'11) и Ц'11). Элементы матрицы ошибок Е отбрасываются с корректировкой соответствующих диагональных значений. В работе [6] приводятся некоторые теоретические оценки качества предобусловливания в симметричном случае. Этот алгоритм показал высокую эффективность и надёжность при последовательном решении широкого класса задач с плохо обусловленными матрицами коэффициентов. В случае мелкоблочных вычислений операции с элементами заменяются на соответствующие операции с блоками.
Параллельный алгоритм построения неполных треугольных разложений на основе геометрической редукции. Основой одного из возможных способов распараллеливания алгоритма построения предобусловливания типа (1) является так называемая геометрическая редукция задачи [7]. Её смысл состоит в следующем:
1. Строится геометрическая объёмная декомпозиция области нар частей (рис. 3).
2. В объемах выделяются строки, в столбцах которых нет элементов из других объемов, т.е. исходная матрица приводится к окаймленному блочнодиагональному виду
0 С,
0 Ар Ср
В1 • Вр D
(2)
где Ai - независимые диагональные блоки, Bi и Ci - соответствующие им блоки окаймления, г = = 1,.. .,р, а D - последний диагональный блок.
3. Для матрицы дополнения по Шуру
р
S = D - ^ В1А-1С1 строится ее структурное при-
г=1
ближение вида $р^=$р ^ +^Г^р(В1 )*Sp(Zi),
г=1
где Sp(Zi) = Sp[(Aг. )ЫСТСХЕ С1 ], КСУСЬЕ > 1. Через Sp(Z) здесь обозначена структура разреженности матрицы ^
4. В строках матрицы D выделим наборы строк Му, 1 - г - Р, 1 - 7 - Р, г < 7, такие, что в каждой строке из набора Му в столбце есть элементы как из г-й, так и из 7-й подобласти исходной геометрической декомпозиции области, а элементы из других исходных объемов отсутствуют. Множества ячеек из наборов Му имеют в 3D случае квазиразмерность 2.
5. По подматрице Sp(Dij) матрицы Sp(D)
для строк и столбцов из множества Мг7 проводится геометрическая декомпозиция на две части с выделением границы между ними (рис. 4). Граница раздела имеет квазиразмерность 1 - это некоторая квазикривая в 3D случае.
6. Аналогично в строках матрицы D выделяются тройные наборы строк М^, квазикривые в 3D случае, для которых проводится декомпозиция на 3 части.
7. Каждое множество из больше чем 3 соседних объемных блоков в строке рассматривается как целое - квазиточка в 3D случае.
Таким образом, строится некоторое блочное биение задачи. Упорядочивание блоков осуществляется следующим образом: сначала все объемы, затем все квазиповерхности, затем все квазикривые, затем все квазиточки. Распределение по процессорам осуществляется с учетом
геометрической локальности: объемы - каждый своему процессору, квазиповерхности - каждая из частей одному из процессоров в примыкающих объемах, квазикривые - каждая из частей одному из трех примыкающих процессоров, и т.д. Параллельный алгоритм построения предо-бусловливания строится на основе геометрической редукции как «блочный по позициям» с неполным разложением «по значениям» второго порядка точности внутри каждого блока.
Итерационный алгоритм SOFGMRES(m). Рассмотрим систему линейных уравнений вида Ах = Ь, (3)
где А е RИхИ - заданная невырожденная, возможно незнакоопpеделенная несимметричная матрица, Ь е RИ - заданный вектор правой части, х е RИ - неизвестный вектор. В дальнейшем будем предполагать, что ^едобусловливатель и ненулевое начальное приближение, если таковые имеются, уже учтены соответствующим образом в А и Ь.
Основой алгоритма SOFGMRES(m) являются матричные соотношения [8] вида:
' АУк=^кЯк,
‘УТ¥к=1к, (4)
У^к=1к,
где Ук е Rтк , Wk е Rтк , Як е Rкхк, матрица Rk - верхняя треугольная. Соотношения (4) будем называть матричными соотношениями в QR форме. В процессе итераций матричные соотношения (4) достраиваются, строится новое приближение к решению, затем построенные расширенные матричные соотношения фильтруются, новые направления пополняют начальное подпространство, и т.д. Оказывается, что если начальное подпространство и фильтрации отсутствуют, то при некоторых необремени-тельных предположениях алгоритм SOFGMRES(m) математически эквивалентен алгоритму GMRES без рестартов. Оказывается также, что сходимость алгоритма SOFGMRES(m) определяется монотонной характеристикой подпространства -обобщенным обратным числом обусловленности предобусловленной матрицы на подпространстве. В симметричном положительно определенном случае если подпространство пусто, эта характеристика совпадает с обратным числом обусловленности предобусловленной матрицы в обычном смысле. Кроме того, эта характеристика в общем случае принимает значения в интервале [0,1], монотонно не убывает при расширении подпространства, и чем она ближе к 1, тем
быстрее сходимость алгоритма SOFGMRES(rn) к решению системы уравнений (3).
Существует множество способов построения начального подпространства: а) итерационный, алгоритм SOFGMRES(m); б) на основе априорной информации; в) последовательным уточнением подпространства при решении систем уравнений с последовательно возникающими правыми частями; г) через предысторию вычислений - финальное подпространство для предыдущей системы уравнений можно свести к набору полей моделируемых физических величин, а затем на этой основе построить начальное подпространство для следующей системы линейных уравнений, в том числе и при небольших изменениях расчетной области.
Высокоточное решение плохообусловленных СЛАУ. Идея высокоточного решения плохообусловленных СЛАУ состоит в использовании повышенной точности для небольшой части вычислений. Высокоточные вычисления при решении плохообусловленных СЛАУ проводятся на основе программной эмуляции в пакете QD [5]. Основная идея этого пакета состоит в том, чтобы хранить и проводить вычисления с «четверной» double-double и «восьмерной» quad-double точностью как с не вычисленной явным образом суммой двух или четырех вещественных чисел двойной точности. К сожалению, явное использование программной эмуляции повышенной точности во всех вычислениях является слишком дорогой операцией (~ 10-кратное замедление для double-double и —100-кратное замедление для quad-double), поэтому для эффективности высокоточных алгоритмов необходимо минимизировать вычисления с кратной точностью.
Подобная методика подробно рассматривается в работе [9], где описывается двухуровневый алгоритм, в котором повышенная точность используется только на верхнем уровне итерационной схемы, на котором вычисляется текущая невязка в высокой точности. Далее невязка преобразуется в двойную точность, и с ней, как с обычной правой частью, с сохранением подпространства на входе/выходе, проводятся вычисления в double для нахождения в double нового пробного вектора решения. Этот пробный вектор преобразуется назад в высокую точность, строится новое приближение к решению в высокой точности и новый вектор невязки, и т.д. При такой организации вычислений основные затраты - решение системы уравнений в double для последовательности правых частей -текущих невязок системы уравнений.
Таблица 2
Основные параметры решения тестовых задач
№ модели Кол-во неизвестных Размер блока DShift T1 T2 Время, c
1 40 020 3 1є-6 1e-3 1e-6 0.42
2 88 782 3 1e-5 3e-3 1e-5 3.65
3 98 б50 2 1e-6 1e-3 1e-6 8.59
4 112 5б0 3 1e-6 1e-3 1e-6 4.10
5 195 024 3 1e-6 1e-3 1e-6 10.4
б 247 542 3 1e-5 3e-3 1e-5 17.0
7 25б 572 3 1e-5 3e-3 1e-5 16.8
8 397 950 3 1e-4 1e-2 1e-4 53.6
9 427 014 3 1e-6 1e-3 1e-6 41.6
10 755 592 3 1e-7 3e-4 1e-7 101.
11 798 230 2 3e-4 3e-3 1e-5 242.
12 889 890 3 1e-5 3e-3 1e-5 72.
13 1 б23 б00 3 1e-10 1e-5 1e-10 113.
14 2 534 44б 2 1e-6 1e-3 1e-6 431.
15 3 247 200 3 1e-10 1e-5 1e-10 351.
1б 4 870 800 3 1e-12 1e-6 1e-12 823.
17 8 033 520 3 1e-8 1e-6 1e-10 766. (*)
Таблица 3
Влияние выбора параметра размера блока на сходимость итераций для задачи № 13
Размер блока, b nzA(b)/nzA(1) nzU(b)/nzAU(b) Кол-во итераций Время(Ргес), c Время(Вет), c Время, c
1 1.00 7.40 58 393.3 146.0 539.3
2 1.60 7.47 55 83.8 81.9 165.7
3 2.05 7.28 60 46.0 77.0 123.0
4 2.62 7.24 61 51.3 84.7 136.0
5 3.20 7.05 60 58.6 79.6 138.2
6 2.16 7.46 65 60.1 96.1 156.2
8 3.87 6.37 54 119.3 96.2 215.5
10 4.42 6.25 55 129.7 107.0 236.7
Результаты численных экспериментов
В таблице 2 для набора тестовых задач представлены оптимальные параметры предобу-словливания при решении СЛАУ с первой правой частью в режиме 1 MPI х 4 threads на одном узле суперкомпьютера «Ломоносов». Параметр «Размер блока» в таблице означает размер блока мелкоблочного биения матрицы и предо-бусловливателя при решении СЛАУ, параметр «DShift» означает диагональный сдвиг после мелкоблочного масштабирования, параметры т и т2 обозначают мелкоблочные пороги фильтрации первого и второго порядков точности,
соответственно. Критерий остановки итерационной схемы - уменьшение относительной нормы невязки в 109 раз. Результат для наибольшей по размерности задачи № 17 приведен для расчета на 4 узлах, т.к. памяти одного узла не хватает для построения предобусловливателя.
В таблице 3 для задачи № 13 приведены результаты решения СЛАУ с оптимальным набором параметров предобусловливания для различных вариантов выбора размера блока мелкоблочного биения матрицы. В таблице 3 в том числе приведены значения относительного заполнения матрицы (п2А(Ь)/п2А(1)) по отношению к хранению данных в точечном формате и
Таблица 4
Время решения задачи № 13 для одной и двух правых частей различными итерационными методами
Итерационный метод № правой части Кол-во итераций Время, с
CG 1 400 199.3
2 402 199.4
S0FGMRES(10) 1 337 348.6
2 (с сохранением подпространства от 1-й) 45 48.5
Таблица 5
Время расчета в многопроцессорном режиме
N° модели 1 (1*1) 2 (1*2) 4 (1*4) 8 (2*4) 16 (4*4) 32 (8*4)
1 0.45+0.73(38) 0.23+0.83(66) 0.13+0.51(66) 0.08+0.97(187) - -
2 2.5+4.0(62) 1.2+5.1(126) 0.7+3.3(126) 0.3+2.9(197) - -
3 7.9+11.5(107) 3.8+15.3(248) 2.0+10.0(248) - - -
4 3.1+7.5(90) 1.6+6.3(123) 0.8+4.3(124) 0.5+4.0(209) - -
5 11+8(48) 6+15(150) 3+10(150) 2+12(243) 1+9(363) 0.4+7.8(473)
б 13+17(88) 8+13(117) 5+9(117) 2+9(200) - -
7 20+16(82) 10+11(96) 6+7(96) 3+6(139) 1+7(221) -
8 35+74(198) 18+55(198) 10+34(198) 5+18(204) - -
9 30+55(160) 19+53(218) 12+31(218) - - -
10 135+168(180) 75+140(240) 49+100(240) 23+92(400) 11+91(630) 6+78(840)
11 174+191(173) 88+135(177) 116+126(197) - - -
12 39+75(98) 19+52(99) 12+42(101) 8+20(99) - -
13 165+60(30) 86+63(46) 44+56(46) 21+37(46) 12+23(45) -
14 443+476(141) 256+536(181) 132+349(206) - - -
15 332+754(160) 169+580(160) 90+508(164) 48+287(176) 36+157(208) 23+151(367)
16 - - 137+686(139) - - -
17 - - - - 96+670(350) -
предобусловливания по отношению к хранению матрицы в блочном формате (nzU(b)/nzAU(b)). Из таблицы следует, что для задачи № 13 оптимальным размером блока мелкоблочного биения является значение Ь = 3.
В таблице 4 для задачи № 13 приведено время итерационной схемы при решении задачи с несколькими правыми частями. Из таблицы видно, что алгоритм SOFGMRES при меньшем числе итераций, чем в CG, для первой правой части требует большего времени счета из-за дополнительных затрат на ортогонализации. С другой стороны, для второй правой части, при условии сохранения подпространства, алгоритм SOFGMRES находит решение за значительно меньшее время. При большом числе правых частей, что часто имеет место в задачах строительной механики, использование алгоритма SOFGMRES может быть более предпочтительным.
Для некоторых задач были оценены числа обусловленности отмасштабированных к единичной диагонали матриц жесткости: для задачи № 9 получена оценка — 3-109, для задачи № 13 — 2-1010, для задачи № 15 — 8-1011, для задачи № 16 — 4-1012, для задачи № 17 — 3-1013. Алгоритм вычисления оценок числа обусловленности для симметричных матриц см. в [9]. При таких больших значениях чисел обусловленности матриц, достаточно близких к обратной точности представления числа в формате double, трудно ожидать высокой точности вектора решения СЛАУ даже при вычислении этого решения прямыми методами в формате double. По видимому, для получения высокоточного решения подобных СЛАУ необходимо использовать повышенную точность представления чисел, в том числе при построении исходной матрицы коэффициентов.
В таблице 5 для всего набора тестовых задач представлены результаты некоторых парал-
лельных расчетов на 1, 2, 4 и 8 узлах суперкомпьютера «Ломоносов» с распараллеливанием по нитям вычислений (от 1 до 4 нитей). Приводятся времена построения предобуслов-ливателя плюс время итерационной схемы (в скобках указано количество итераций).
Выводы
В работе приведены данные по решению итерационными методами плохообусловленных симметричных СЛАУ для набора задач строительной механики в различных параллельных режимах на суперкомпьютере «Ломоносов». Из полученных оценок числа обусловленности для некоторых матриц следует, что эти задачи являются чрезвычайно трудными для решения итерационными методами. Результаты экспериментов показывают, что тем не менее все задачи были успешно решены. Кроме того, для некоторых задач рассмотренные итерационные методы демонстрируют неплохую вычислительную эффективность и параллельную масштабируемость. С другой стороны, для ряда самых трудных задач не удалось получить желаемую эффективность параллельных итерационных алгоритмов. По мнению авторов, это связано со значительным ухудшением качества предобусловливания при априорном блочном ограничении заполнения предобусловливателя при распараллеливании по процессорам/нитям. Кроме того, не использованы в полной мере возможности улучшения эффективности алгоритмов за счет применения блочных (многовекторных) итерационных алгоритмов. Использование мелкоблочного формата хранения матрицы и предобусловливателя, а также блочных итерационных алгоритмов позволяет рассчитывать, что предложенные итерационные алгоритмы могут быть эффективно реализованы на гибридных вычислительных системах с ускорителями, например, на основе GPU. В последующих работах авторы планируют развивать итерационные методы в этих направлениях.
Список литературы
1. Симбиркин В.Н., Фросин А.В. Развитие и взаимодействие компонентов автоматизированной системы инженерного анализа СТАРКОН для массового применения в строительстве // Строительная механика и расчет сооружений. 2010. № 5. С. 67-71.
2. Жук Ю.Н., Симбиркин В.Н. Программный комплекс STARK ES // В кн.: Современное высотное строительство. М.: ИТЦ Москомархитектуры, 2007. 464 с.
3. Программный комплекс для расчета строительных конструкций на прочность, устойчивость и колебания STARK ES. Версия 4х4 (2007). Руководство пользователя. М.: Еврософт, 2008. 399 с.
4. Сушко Г.Б., Харченко С.А. Экспериментальное исследование на СКИФ МГУ «Чебышев» комбинированной MPI+threads реализации алгоритма решения систем линейных уравнений, возникающих во FlowVision при моделировании задач вычислительной гидродинамики // Материалы Международной научной конференции «Параллельные вычислительные технологии» (ПаВТ'2009), Н. Новгород, 30 марта
- 3 апреля 2009 г. Челябинск: Изд-во ЮУрГУ, 2009. C. 316-324.
5. http://crd.lbl.gov/—dhbailey/mpdist/ (дата обращения: 03.11.2011)
6. Kaporin I.E. High quality preconditioning of a general symmetric positive definite matrix based on its UTU + UTR + R U decomposition // Numer. Linear Algebra Appl. 1998. V. 5. P. 483-509.
7. Харченко С.А. Влияние распараллеливания вычислений с поверхностными межпроцессорными границами на масштабируемость параллельного итерационного алгоритма решения систем линейных уравнений на примере уравнений вычислительной гидродинамики // Материалы Международной научной конференции «Параллельные вычислительные технологии» (ПаВТ'2008), СПб., 28 января-1 февраля 2008 г. Челябинск: Изд-во ЮУрГУ, 2008. C. 494-499.
8. Харченко С.А. Параллельная реализация итерационного алгоритма решения несимметричных систем линейных уравнений с частичным сохранением спектральной/сингулярной информации при явных рестартах // Вычислительные методы и программирование. 2010. T. 11. С. 373-381.
9. Харченко С.А. Высокоточный двухуровневый параллельный итерационный алгоритм решения плохо-обусловленных СЛАУ // Труды Международной супер-компьютерной конференции «Научный сервис в сети Интернет: экзафлопсное будущее». 2011. С. 50-55.
SOLUTION OF ILL-CONDITIONED SYMMETRIC SLAE FOR STRUCTURAL MECHANICS PROBLEMS BY PARALLEL ITERATIVE METHODS
V.L. Yakushev, V.N. Simbirkin, A V. Filimonov, P.A Novikov,
I.N. Konshin, G.B. Sushko, S.A Kharchenko
The results of numerical experiments using the iterative methods are presented for the solution of ill-conditioned symmetric SLAE for a number of structural mechanics problems. The second order incomplete triangular factorization has been used to construct preconditionings. MPI and TBB technologies have been exploited to parallelize computations. Numerical experiments in different parallel regimes have been performed on the Lomonosov supercomputer.
Keywords: structural mechanics problems, system of linear algebraic equations, iterative methods, incomplete triangular factorization, parallel computing, message passing interface (MPI), multithreading, the Lomonosov supercomputer.
