Stablizing solution of large sprase system

[andrewwinters5000]

As the sparse matrix for, in particular, the BicubicBSpline grows in size the built in \ operator may fail. From the docstring of the \ operator

When A is sparse, a similar polyalgorithm is used. For indefinite matrices, the LDLt factorization does not use pivoting during the numerical factorization and therefore the procedure can fail even for invertible matrices.

As the matrix for the Bicubic B-spline is indefinite, the lack of pivoting is an issue. The current strategy from #48 is to compute the qr factorization prior to the \ call. This serves to stabilize the computation of the B-splines coefficients, but may not be optimal.

One could add LinearSolve.jl as a dependency which opens the possibility of many sparse solvers. Then we would replace a line like

    x = qr(A) \ b

with

    Ax_b = LinearProblem(A, b)
    x = solve(Ax_b, selected_solver()).u

Some preliminary testing for an approximately 100_000 x 100_000 sparse system for particular choices of selected_solver() are

• QRFactorization() perform equivalent to qr(A) \ b, which is not surprising as LinearSolve.jl calls the underlying built in qr factorization from LinearAlgebra
• UMFPACKFactorization() fails
• KLUFactorization() works but is about 20x time slower than the QR strategy. The QR takes ≈13.5s and KLU takes ≈242s
• Many iterative solver options are available like KrylovJL_GMRES or KrylovJL_BICGSTAB but testing would need to be done to find a "good" preconditioner and initial guess.

Including the LinearSolve.jl as a dependency gives flexibility with respect to solver selection but also some pre-compilation overhead. Whether or not this is "worth it" remains to be seen and requires more testing. In particular, with respect to Krylov-type methods for very large datasets. However, for modest datasets (up to ≈ 500_000 x 500_000) the qr strategy has proven sufficient.

Finally, below is a spy plot of the sparsity pattern for the approximately 100_000 x 100_000 sparse matrix for the Bicubic B-spline with "not-a-knot" boundary conditions. There is a predominant block-tridiagonal structure but the boundary coupling at the far right may introduce some complexity. Perhaps seeing the structure could help others more experienced with numerical linear algebra guide the choice of a sparse solver.