And then there is SparsityDetection.jl which can detect the sparsity pattern automatically:
using SparsityDetection
input = rand(n)
output = similar(input)
sparsity_pattern = jacobian_sparsity(ffd,output,input)
jac = Float64.(sparse(sparsity_pattern))
colors = matrix_colors(jac)
forwarddiff_color_jacobian!(jac, ffd, X, colorvec = colors)
resulting in
According to the test, the sparsity detection algorithm implemented exhibits O(N) behavior. In this example it is beaten by the atomic assembly. However, e.g. in an implementation of Newton’s algorithm and in many other cases one would perform the sparsity detection and matrix coloring only once and re-use the results during the computation of series of jacobians, so the relative overhead from the detection step would decrease.
I find the O(N) behavior of the sparsity detection step quite amazing. The algorithm behind this is described here and uses more deep Julia stuff…