Графическое представление разреженных матриц¶
В этом примере мы создадим сетку для расчета в конечных элементах аэродинамического профиля NASA - крыло и два закрылка. Полученный граф мы пропустим через разные алгоритмы, чтобы показать, как можно оптимизировать представление матрицы смежности графа.
Загрузка данных¶
Для дальнейшей работы нам потребуются библиотеки для визуализации графов, специальные функции из библиотеки SparseArrays (написанной на Си и оптимизированной), а также две библиотеки функций, написанные на Julia.
Pkg.add( ["Graphs", "GraphPlot", "AMD", "SymRCM"] )
Данные хранятся в файле airfoil.jld2 и содержат:
- 4253 пары координат (x,y) точек сетки
- Массив из 12 289 пар индексов (i,j), задающих соединения между точками сетки
Загрузим файл данных в рабочую область.
using FileIO
i,j,x,y = load( "$(@__DIR__)/airfoil.jld2", "i", "j", "x", "y" );
Визуализация конечно-элементной сетки¶
Создадим разреженную матрицу смежности из вектора соединений (i,j) и сделаем её положительно определенной. Затем построим график матрицы смежности, используя (x,y) как координаты вершин (точек сетки).
Для дополнительного ускорения проинструктируем систему сохранять графики в формате PNG и не использовать интерактивную визуализацию.
gr( fmt = "png", size=(600,400) );
using SparseArrays, Graphs, GraphPlot
n = convert(Int32, max( maximum(i), maximum(j) ));
A = sparse( vec(i), vec(j), -1.0, n, n );
A = A .+ A';
gplot( Graph( A ), vec(x), .-vec(y), EDGELINEWIDTH=0.3,
edgestrokec = "deepskyblue", nodefillc = "deepskyblue",
plot_size = (16cm, 12cm), title="Профиль крыла и закрылков" )
Визуализация матрицы смежности¶
Функция spy позволяет визуализировать ненулевые элементы матрицы, что особенно полезно для анализа структуры разреженных матриц. spy(A) отображает структуру ненулевых элементов матрицы A.
spy( A, c="deepskyblue",
title="Матрица смежности графа профиля крыла",
titlefont=font(10) )
Симметричное переупорядочивание - обратный алгоритм Катхилла-Макки¶
Алгоритм Катхилл—Макки перенумеровывает вершины графа с тем, чтобы его матрица смежности (изучаемая нами) имела наименьшую ширину ленты. Обратный алгоритм (RCM, reverse Cuthill—McKee) просто нумерует индексы вершин в обратную сторону.
Алгоритмы для RCM появились давно, можно использовать и более старую имплементацию в пакете CurthillMcKee.jl.
symrcm реализует обратное преобразование Катхилла-Макки для переупорядочивания матрицы смежности. Команда r = symrcm(A) возвращает вектор перестановки r, при котором матрица A[r,r] имеет диагональные элементы, расположенные ближе к главной диагонали, чем в исходной матрице A.
using SymRCM
r = symrcm(A)
spy( A[r,r], c="deepskyblue",
title="Упорядочивание Катхилла-Макки (обратное)",
titlefont=font(10) )
Такой формат полезен как предварительное упорядочивание для LU-факторизации или разложения Холецкого. Работает как для симметричных, так и для несимметричных матриц.
Исходная матрица смежности может иметь ненули "размазанные" по всей ее площади. RCM-перестановка упорядочивает точки так, что ненули концентрируются у диагонали, после чего вычисления могут вестись более эффективно.
Симметричное переупорядочивание - перестановка столбцов¶
В Julia нет функции по умолчанию для получения вектора перестановки, упорядочивающего столбцы разреженной матрицы A в порядке неубывания количества ненулевых элементов. Такой подход тоже иногда полезен как предварительное упорядочивание для LU-разложения, например lu(A[:,j]).
Поскольку суть операции colperm крайне проста, легче написать собственную реализацию:
function colperm(S::SparseMatrixCSC)
nnz_per_col = sum(S .!= 0, dims=1) # Количество ненулевых элементов в каждом столбце
return sortperm(vec(nnz_per_col)) # Возвращает перестановку индексов
end
j = colperm( A )
spy( A[j,j], c="deepskyblue",
title="Упорядочивание по количеству элементов в столбцах",
titlefont=font(10) )
Эта операция позволяет добиться более эффективного сжатия матрицы, а также может уменьшить количество ненулевых элементов, образуемых командами LU разложения.
Симметричный приближенный минимальный порядок¶
Команда amd из соответствующего пакета выполняет симметричную перестановку по методу приближенного минимального порядка. Для симметричной положительно определенной матрицы A команда p = amd(A) возвращает вектор перестановки p, при котором матрица A[p,p] имеет более разреженный множитель Холецкого, чем исходная A. Иногда amd хорошо работает и для симметричных неопределенных матриц.
using AMD
m = amd(A);
spy( A[m,m], c="deepskyblue",
title="Приближенный минимальный порядок",
titlefont=font(10) )
Опять же, мы находим перестановку, при которой последующие алгоритмы разложения вводят меньше новых ненулей, при этом сохраняется симметрия (важно для задач механики и физики).
Заключение¶
Все эти алгоритмы помогают выполнить подготовку данных для быстрого и более стабильного решения СЛАУ, которые возникают в расчетах прочности крыла, аэродинамическом моделировании, тепловом анализе и т.д. Применив оптимизацию мы, скорее всего, сможем сократить время решения задачи на 2-3 порядка.
Дополнительные алгоритмы есть, например, в пакетах Metis.jl и IterativeSolvers.jl. В первом собраны алгоритмы переупорядочивания графов, во втором дополнительные методы для работы с разреженными матрицами. Рекомендуем также изучить SuiteSparse - обертку над хорошо проверенным промышленным кодом на Си, тем не менее лишена некоторых функций, которые нам потребовались.