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-Mackey 算法

Cuthill-McKee 算法对图形的顶点重新编号，使其邻接矩阵（我们正在研究的矩阵）的带宽最小。反向 Cuthill-McKee 算法（RCM，reverse Cuthill-McKee）只是以相反的方向重新编号顶点索引。

RCM 算法很早以前就出现了，你可以使用 CurthillMcKee.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 因式分解或 Cholecki 分解的预排序。它既适用于对称矩阵，也适用于非对称矩阵。

原始邻接矩阵的整个区域可能都 "散布 "着非零点。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] 的 Choletsky 乘数比原来的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

同样，我们可以找到一种排列方式，在这种排列方式中，后续分解算法引入的新非零较少，同时又能保持对称性（这对力学和物理学中的问题很重要）。

结论

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

其他算法可在Metis.jl 和IterativeSolvers.jl 软件包中找到。前者包含图形重排序算法，后者包含处理稀疏矩阵的附加方法。我们还建议研究SuiteSparse - 它是对经过良好测试的工业 C 代码的封装，但缺少我们需要的一些功能。