Engee documentation
Notebook

Graphical representation of sparse matrices

In this example, we will create a mesh for a finite element computation of a NASA airfoil - a wing and two flaps. We will run the resulting graph through different algorithms to show how the representation of the graph adjacency matrix can be optimised.

Data loading

For further work we need libraries for graph visualisation, special functions from the SparseArrays library (written in C and optimised), and two function libraries written in Julia.

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

The data is stored in the file airfoil.jld2 and contains:

  • 4253 pairs of coordinates (x,y) of grid points
  • An array of 12,289 pairs of indices (i,j) defining connections between grid points

Let's load the data file into the workspace.

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

Visualisation of the finite element mesh

Create a sparse adjacency matrix from the connection vector (i,j) and make it positively defined. Then plot the adjacency matrix using (x,y) as coordinates of the vertices (grid points).

For added speed, instruct the system to save the graphs in PNG format and not to use interactive visualisation.

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]:
Профиль крыла и закрылков

Visualisation of the adjacency matrix

The function spy allows you to visualise the non-zero elements of a matrix, which is particularly useful for analysing the structure of sparse matrices. spy(A) displays the structure of the non-zero elements of matrix A.

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

Symmetric reordering - inverse Cuthill-Mackey algorithm

The Cuthill-McKee algorithm renumbers the vertices of a graph so that its adjacency matrix (the one we are studying) has the smallest tape width. The reverse Cuthill-McKee algorithm (RCM, reverse Cuthill-McKee) simply renumbers the vertex indices in the opposite direction.

Algorithms for RCM appeared long ago, you can use an older implementation in the package CurthillMcKee.jl.

symrcm implements the inverse Cuthill-McKee transform to reorder the adjacency matrix. The command r = symrcm(A) returns the permutation vector r, where the matrix A[r,r] has diagonal elements closer to the main diagonal than in the original matrix 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

This format is useful as a preorder for LU factorisation or Cholecki decomposition. It works for both symmetric and asymmetric matrices.

The original adjacency matrix may have non-zeros "smeared" throughout its area. RCM rearrangement arranges the points so that the non-zeros are concentrated at the diagonal, after which computations can be performed more efficiently.

Symmetric reordering - column rearrangement

There is no default function in Julia to obtain a permutation vector that orders the columns of a sparse matrix A in non-decreasing order of the number of non-zero elements. This approach is also sometimes useful as a preorder for LU decomposition, e.g. lu(A[:,j]).

Since the essence of the colperm operation is extremely simple, it is easier to write your own implementation:

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

This operation can achieve more efficient matrix compression, and can also reduce the number of non-zero elements formed by LU decomposition commands.

Symmetric approximate minimum order

The command amd from the corresponding package performs symmetric permutation using the approximate minimum order method. For a symmetric positive definite matrix A, the command p = amd(A) returns a permutation vector p, where the matrix A[p,p] has a more sparse Choletsky multiplier than the original A. Sometimes amd works well for symmetric indefinite matrices as well.

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

Again, we find a permutation in which subsequent decomposition algorithms introduce fewer new nonzeros while preserving symmetry (important for problems in mechanics and physics).

Conclusion

All these algorithms help to prepare data for faster and more stable solution of SLAEs that arise in wing strength calculations, aerodynamic modelling, thermal analysis, etc. By applying optimisation we are likely to be able to reduce the time to solve the problem by 2-3 orders of magnitude.

Additional algorithms are available, for example, in the packages Metis.jl and IterativeSolvers.jl. The former contains algorithms for graph reordering, the latter contains additional methods for working with sparse matrices. We also recommend to study SuiteSparse - it is a wrapper over well-tested industrial C code, nevertheless it lacks some features we needed.