Engee 文档
Notebook

稀疏矩阵的图形表示

在这个例子中，我们将创建一个网格，用于计算NASA翼型的最终元素-一个机翼和两个襟翼。 我们将通过不同的算法运行生成的图，以显示如何优化图的邻接矩阵的表示。

上传数据

为了进一步的工作，我们将需要用于图形可视化的库，SparseArrays库中的特殊函数（用C编写并优化），以及两个用Julia编写的函数库。

In [ ]: 
Pkg.add( ["Graphs", "GraphPlot", "AMD", "SymRCM"] )

数据存储在文件中 airfoil.jld2 并包含:

*4253对网格点的坐标（x，y）
*定义网格点之间连接的12,289对索引(i,j)的数组

将数据文件上传到工作区。

In [ ]: 
using FileIO
i,j,x,y = load( "$(@__DIR__)/airfoil.jld2", "i", "j", "x", "y" );

有限元网格的可视化

让我们从连接向量（i，j）创建一个稀疏邻接矩阵，并使其正定。 然后我们使用（x，y）作为顶点（网格点）的坐标绘制邻接矩阵。

为了获得额外的加速，我们指示系统以PNG格式保存图形，而不是使用交互式可视化。

In [ ]: 
gr( fmt = "png", size=(600,400) );
In [ ]: 
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="Профиль крыла и закрылков" )
Out[0]:
Профиль крыла и закрылков

邻接矩阵的可视化

功能 spy 它允许您可视化非零矩阵元素，这对于分析稀疏矩阵的结构特别有用。 spy(A) 显示矩阵a的非零元素的结构。

In [ ]: 
spy( A, c="deepskyblue",
     title="Матрица смежности графа профиля крыла",
     titlefont=font(10) )
Out[0]:
No description has been provided for this image

对称重排序是Cuthill-McKee算法的逆

Cuthill-McKee算法重新计算图的顶点，以便其邻接矩阵（我们正在研究）具有最小的带状宽度。 反向算法（RCM，reverse Cuthill-McKee）简单地对相反方向的顶点索引进行编号。

RCM的算法很久以前就出现了，你可以在[CurthillMcKee中使用较旧的实现。jl]包（https://github.com/rleegates/CuthillMcKee.jl )。

symrcm 实现逆Cuthill-McKee变换以重新排列邻接矩阵。 专责小组 r = symrcm(A) 返回置换向量 r，其中矩阵 A[r,r] 它具有位于比原始矩阵A中更靠近主对角线的对角元素。

In [ ]: 
using SymRCM

r = symrcm(A)
spy( A[r,r], c="deepskyblue",
     title="Упорядочивание Катхилла-Макки (обратное)",
     titlefont=font(10) )
Out[0]:
No description has been provided for this image

这种格式作为LU因子分解或Cholesky分解的预排序是有用的。 它适用于对称和非对称矩阵。

原始邻接矩阵可能在其整个区域上具有非零"涂抹"。 RCM置换排列点，使非零集中在对角线上，之后可以更有效地执行计算。

对称重新排序-重新排列列

Julia没有用于获取对稀疏矩阵的列进行排序的置换向量的默认函数。 A 非零元素的数量不递减的顺序。 这种方法有时也可用作LU分解的预排序，例如 lu(A[:,j]).

因为手术的本质是 colperm 非常简单，编写自己的实现更容易。:

In [ ]: 
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) )
Out[0]:
No description has been provided for this image

此操作允许更有效的矩阵压缩，并且还可以减少由LU分解命令生成的非零元素的数量。

对称近似最小阶

专责小组 amd 使用来自相应包的近似最小阶方法执行对称置换。 对于对称正定矩阵 A 专责小组 p = amd(A) 返回置换向量 p，其中矩阵 A[p,p] 有一个比原来的更稀疏的Cholesky乘数 A. 有时 amd 它也适用于对称不定矩阵。

In [ ]: 
using AMD

m = amd(A);
spy( A[m,m], c="deepskyblue",
     title="Приближенный минимальный порядок",
     titlefont=font(10) )
Out[0]:
No description has been provided for this image

再次，我们发现了一种置换，其中后续分解算法引入更少的新nonzeros，同时保持对称性（对于力学和物理学中的问题很重要）。

结论

所有这些算法都有助于准备数据，以便更快，更稳定地解决机翼强度计算，空气动力学建模，热分析等方面出现的弱点。 通过应用优化，我们很可能能够将解决问题所需的时间减少2-3个数量级。

其他算法可用，例如，在包中 Metis.jl 和 IterativeSolvers.jl. 第一个包含用于图重新排序的算法，第二个包含用于处理稀疏矩阵的其他方法。 我们也建议你学习 SuiteSparse -对经过良好测试的工业C代码进行包装，但它缺少我们所需要的一些功能。