Аннотация результатов, полученных в 2022 году
Проект посвящен разработке эффективных итерационных методов решения проблемы собственных значений для симметричных незнакоопределенных матриц на высокопроизводительных вычислительных системах. Нахождение собственных значений и собственных векторов матриц является одной из основных задач для многих разделов вычислительной физики и математики: при исследовании собственных колебаний различных механических систем, колебательных и электронных спектров молекул и кристаллов, в задачах квантовой механики. В настоящее время важной особенностью данных задач является высокий порядок их матриц – до миллиона и более неизвестных. Несмотря на разреженную структуру матриц, из-за высоких порядков матриц решение таких задач требует значительных вычислительных ресурсов, которые могут быть обеспечены современными высокопроизводительными вычислительными системами. В рамках проекта был проведен сравнительный анализ методов, основанных на подпространствах Крылова и предназначенных для решения как стандартной, так и обобщенной симметричной задач на собственные значения. Были исследованы следующие классы методов: обобщенный степенной метод, стандартный метод Ланцоша, метод Ланцоша со сдвигом и обращением, методы итерирования подпространства и метод обратных итераций. https://inftech.uginfo.sfedu.ru/sites/default/files/Abstract_SITO22_fin.pdf (стр.193-198) В конкретной вычислительной задаче требования достижения заранее заданной точности, надежности, скорости исполнения, разумных запросов к памяти, краткости программы и возможности эффективного распараллеливания не всегда бывает возможным удовлетворить в комплексе. Метод, наилучший для одной задачи, оказывается непригодным или малоэффективным для другой. В рамках проекта сначала была рассмотрена и решена стандартная задача на собственные значения с симметричной квазиопределенной (SQD) матрицей А, которая является невырожденной незнакоопределенной, а далее была исследована обобщенная задача на собственные значения, возникающая при конечно-элементном моделировании электроэластичных материалов. Для решения стандартной задачи на собственные значения наилучшие результаты по эффективности и точности вычислений показал метод Ланцоша. Было также апробировано применение спектральной трансформации «сдвиг-обращение», которая необходима для поиска внутренних точек спектра. Крайние (максимальные и минимальные) собственные числа, достаточно хорошо отделенные друг от друга, надежно находятся алгоритмом Ланцоша. Было проведено численное сравнение различных вариантов переортогонализации векторов Ланцоша. Расчеты показали, что как частичная, так и полная переортогонализация улучшают значительно точность вычислений, не требуя при этом существенных дополнительных временных затрат (с учетом распараллеливания ряда операций таких, как матрично-векторные произведения). Обобщенный степенной метод, метод обратных итераций и итерирования подпространства не зарекомендовали себя при решении стандартной задачи на собственные значения. Поэтому при переходе к более сложной, обобщенной задаче, выбор был сразу сделан в пользу методов типа Ланцоша. Проведенные исследования показали, метод Ланцоша неприменим непосредственно к обобщенной задаче. Точнее, его можно использовать, но скорость сходимости оказывается предельно низкой. Однако, выход есть – это использование спектральных трансформаций типа «сдвиг-обращение». Они значительно увеличивают скорость сходимости метода Ланцоша и позволяют находить как крайние (экстремальные), так и внутренние точки спектра. При использовании спектральной трансформации необходимо производить полную переортогонализацию векторов Ланцоша на каждом шаге алгоритма, поскольку, в противном случае, метод начинает создавать «копии» уже сошедшихся собственных значений. Как и для первой задачи, при эффективном распараллеливании матрично-векторных умножений, эта процедура не оказывается сильно затратной по времени вычислений. Для обеих рассмотренных задач критерием окончания итераций служит малая величина невязки, которая находится из полученной трехдиагональной матрицы небольшого размера. Размерность подпространства Крылова бралась, как правило, порядка 2n1/2, где n – размерность исходных матриц. Собственные значения трехдиагональной матрицы находились с использованием QR-алгоритма. Были проведены исследования эффективности обобщенного степенного метода в сравнении с другими методами: стандартный метод Ланцоша, метод Ланцоша со сдвигом и обращением, методы итерирования подпространства и метод обратных итераций. Для решения стандартной задачи на собственные значения наилучшие результаты по эффективности и точности вычислений показал метод Ланцоша. Было также апробировано применение спектральной трансформации «сдвиг-обращение», которая необходима для поиска внутренних точек спектра. Крайние (максимальные и минимальные) собственные числа, достаточно хорошо отделенные друг от друга, надежно находятся алгоритмом Ланцоша. Было проведено численное сравнение различных вариантов переортогонализации векторов Ланцоша. Расчеты показали, что как частичная, так и полная переортогонализация значительно улучшают точность вычислений, не требуя при этом существенных дополнительных временных затрат. При переходе к более сложной, обобщенной задаче, выбор следует делать в пользу методов типа Ланцоша. В рамках реализации проекта были разработаны программные продукты, реализующих созданные алгоритмы на параллельных компьютерах. Для эффективного использования высокопроизводительных рабочих станций с общей памятью и большим количеством ядер ЦПУ были реализованы параллельные версии расчетных модулей. Использовались два подхода: автоматическое распараллеливание циклов с помощью инструментов платформы .net и реализация параллельного умножения матрицы на вектор для формата CRS. В первом случае выигрыш во времени получить не удалось, так как большинство циклов, не задействованных в матричных операциях (поиск нормы вектора, умножение вектора на число) не играют существенной роли в общем времени вычислений и достаточно хорошо оптимизируются на этапе компиляции. Параллельное умножение матрицы было реализовано путем ручного разделения строк матрицы на группы. Для увеличения скорости доступа к данным, матрицы, которые хранились как верхне-треугольные с учетом симметричности, были дополнительно заполнены элементами из нижне-треугольной части. Это позволило независимо выделить строки матриц и разбить их на блоки по количеству потоков для обработки. Реализовав операцию извлечения ненулевых элементов строки для умножения с последующим умножением на вектор, удалось сократить время решения СЛАУ на 40–45% при использовании 8 потоков для матрицы порядка 2 тысяч строк и на 45-50% для матрицы порядка 20 тысяч строк. При этом использовалась единая структура данных для разреженной матрицы, чтение из которой осуществлялось из всех потоков. Аналогично, общие структуры данных использовались для участвующих в операциях векторов. Возможен дальнейший прирост производительности за счет использования блочной структуры матриц. Наиболее перспективным направлением дальнейших исследований является переход к векторным вычислениям с использованием технологии CUDА. В рамках текущих работ были подготовлены тестовые программы для представления разреженных матриц в формате для обработки на GPU и изучены существующие пакеты программ. На этапе численных экспериментов предстоит оценить переход от двойной точности, используемой в текущей реализации методов, к одинарной, позволяющей получить максимальный прирост скорости с использованием GPU.

Аннотация результатов, полученных в 2023 году
Проект посвящен разработке эффективных итерационных методов решения проблемы собственных значений для симметричных незнакоопределенных матриц на высокопроизводительных вычислительных системах. В рамках реализации проекта за отчетный период проведено исследование методов Ланцоша со спектральной трансформацией «сдвиг-обращение» для решения обобщенной задачи на собственные значения, которая появляется в результате конечно-элементного моделирования электроэластичных материалов. Особенностью задачи является отсутствие положительной определённости у обеих матриц пучка. Матрица жёсткости, стоящая в левой части уравнения, является симметричной квази-определённой, а матрица масс, стоящая в правой части - положительно полуопределённой (вырожденной). Тем не менее, матричный пучок определён (существует положительно определённая линейная комбинация матриц для некоторых вещественных скаляров), все собственные значения задачи вещественны и положительны. Теория таких пучков близка к стандартной для эрмитовых проблем за исключением того, что некоторые собственные значения являются бесконечными. Метод Ланцоша не применим непосредственно к обобщенной задаче, поэтому используется спектральная трансформация. Существуют сложные вычислительные моменты решения данной задачи с матрицами больших размеров. А именно, эффективное решение незнакоопределённые системы линейных алгебраических уравнений (СЛАУ) во внутреннем цикле метода. Помимо использования стандартных методов для незнакоопределённых систем, предложены и исследованы альтернативные подходы. Один из них – это трансформация к линейной системе с несимметричной положительно определённой матрицей, спектр которой находится в правой полуплоскости и использование авторских ранее разработанных методов решения таких СЛАУ. https://inftech.uginfo.sfedu.ru/sites/default/files/AbstractSITO23_fin.pdf Другой подход – это понижение порядка СЛАУ с использованием дополнения Шура. В этом случае возникает необходимость решения обобщенной задачи на собственные значения с двумя положительно определёнными матрицами. Она решается предобусловленным методом сопряжённых градиентов для системы нормальных уравнений. Предобусловливателем является (1, 1) блок исходной матрицы жёсткости. Преимуществом данного подхода является как более низкий порядок решаемой задачи и определённость матриц пучка, так и отсутствие необходимости поиска эффективного предобусловливателя, т.к. он определяется матрицей исходной задачи. Третий подход к решению СЛАУ состоит в переупорядочивании неизвестных в задаче, что приводит к блочно-структурированным матрицам. Рассматривалось использование методов неполной факторизации (ILU, MILU, итерационного варианта ILU) как для решения структурированной задачи в целом, так и во внутреннем цикле метода Ланцоша со спектральной трансформацией. Проведено исследование критериев остановки внутренних итераций с заданным допуском и по заранее заданному фиксированному числу итераций. Другой особенностью обобщённой задачи на собственные значения является наличие (либо отсутствие) явного нулевого (2, 2) блока у вырожденной матрицы масс. Методы Ланцоша порождают базисные векторы с большими компонентами в ядре матрицы масс и, соответственно, с большой 2-нормой. Если матрица масс M имеет явный нулевой блок, то эти компоненты уничтожаются в M-полускалярном произведении, и v^T Mv≥0. Если вырожденная матрица M не имеет явного нулевого блока, то большие компоненты в ядре M могут доминировать и приводить к серьёзным ошибкам. Предложен способ устранения данной вычислительной сложности, связанный с выбором начального вектора для процесса Ланцоша. Исследованы варианты перезапусков метода Ланцоша, поскольку для нахождения собственных векторов задачи необходимо хранить весь набор полученных векторов. Это требует большого объема памяти при работе алгоритма, причем данный объем нельзя определить заранее, так как неизвестно число итераций, обеспечивающее вычисление собственных значений с требуемой точностью. В ходе исследований применимости предложенных методов для различных прикладных задач в области теории упругости и электро-упругости, решаемых методом конечных элементов, были проведены расчеты в комплексах ACELAN-COMPOS и ANSYS. Рассматривались задачи, требующие проведения модального анализа пьезоустройств, работающих в режиме колебаний. Для таких устройств модальный анализ позволяет получить формы колебаний и частоты резонанса и анти-резонанса, необходимые для последующего амплитудно-частотного анализа. Были проведены исследования для получаемых при решении упомянутых задач конечно-элементных сеток на основе тетраэдров и шестиугольных регулярных гексаэдров. Были проведены расчеты с линейными и квадратичными элементами, получены оценки производительности ключевых операций (матрично-векторное и матрично-матричное умножение) для различных форматов хранения разреженных матриц жесткости и масс, возникающих при различных способах перенумерации узлов конечно-элементных сеток и соответствующих им степеней свободы.