How it works

All pairs, from one sparse factorization.

An algebraic path problem asks for the star of a graph's weighted adjacency matrix, A* = I ⊕ A ⊕ A² ⊕ ⋯, over a semiring. With ⊕ = min and ⊗ = +, A*[i, j] is the shortest distance from i to j. Other semirings ask other questions with the same formula: widest paths (max, min), bottleneck paths (min, max), reachability (or, and), sums over all walks (+, ×).

Floyd–Warshall computes A* with the loops of Gaussian elimination: n³ updates on a dense matrix. PATHPACK.jl uses a sparse direct method instead. It orders the vertices to limit fill, factorizes A into sparse triangular factors L and U over the semiring, and then solves against the identity, row by row, to get A*.

The steps below run on a small street grid. Every order, factor and distance in them is computed in your browser for the graph shown.

Step 1 of 5

The graph is a matrix

144 street corners and 232 streets. Row i of the matrix A holds the streets out of corner i: A[i, j] is the length of the street from i to j. An empty entry is ∞: no street.

PATHPACK.jl computes the star A* = I ⊕ A ⊕ A² ⊕ ⋯. With ⊕ = min and ⊗ = +, its entry A*[i, j] is the distance from i to j. Floyd–Warshall finds it with n³ updates on a dense matrix. PATHPACK.jl keeps A sparse for as long as it can.

The corners are numbered in the order a data file might list them, so the nonzeros look scattered. Hover a corner or an entry.

Graph144 corners, 232 streets

MatrixA, input order

Step 1 of 5

In code

One call does all of it. To answer many single-source questions, keep the factorization: each query is two triangular solves.

using PATHPACK, SparseArrays, LinearAlgebra
using PATHPACK.CPU: MinPlus, mstar, mlu
D = mstar(MinPlus(), A)                   # all of A*: order, factorize, solve
F = mlu(MinPlus(), A)                     # or keep the factors
B = fill(Inf, 1, size(A, 1)); B[s] = 0.0  # the unit row vector at s
rmul!(B, F)                               # row s of A*: the distances from s

The keyword alg of mlu and mstar picks the ordering, for example METIS() (nested dissection) from CliqueTrees.jl. On NVIDIA GPUs, apsp_gpu(A) runs the same method: CPU threads factorize the many small parts, the GPU the large separator blocks, and the rows of A* are split across one GPU or several.