improve performance with sparse matrices (#101)

* implement sparse optimization for out-of-place build_mprk_matrix!

* fix broadcasting

* format

* format

* add package versions to tutorial

* add warning to docs

* fix nonconservative systems

* fix comment

* mention KLUFactorization

* add LinearSolve

* skip broken tests

* Update docs/src/faq.md

Co-authored-by: Stefan Kopecz <[email protected]>

---------

Co-authored-by: Stefan Kopecz <[email protected]>

Author: ranocha

docs/Project.toml

@@ -1,13 +1,19 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
 BenchmarkTools = "1"
 Documenter = "1"
+InteractiveUtils = "1"
+LinearSolve = "2.21"
 OrdinaryDiffEq = "6.59"
+Pkg = "1"
 Plots = "1"
 SparseArrays = "1.7"

docs/make.jl

@@ -76,6 +76,7 @@ makedocs(modules = [PositiveIntegrators],
              "Tutorials" => [
                  "Linear Advection" => "linear_advection.md",
              ],
+             "Troubleshooting, FAQ" => "faq.md",
              "API reference" => "api_reference.md",
              "Contributing" => "contributing.md",
              "Code of conduct" => "code_of_conduct.md",

docs/src/faq.md

@@ -0,0 +1,30 @@
+# Troubleshooting and frequently asked questions
+
+## Sparse matrices
+
+You can use sparse matrices for the linear systems arising in
+[PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl),
+as described, e.g., in the [tutorial on linear advection](@ref tutorial-linear-advection).
+However, you need to make sure that you do not change the sparsity pattern
+of the production term matrix since we assume that the structural nonzeros
+are kept fixed. This is a [known issue](https://github.com/JuliaSparse/SparseArrays.jl/issues/190).
+For example, you should avoid something like
+
+```@repl
+using SparseArrays
+p = spdiagm(0 => ones(4), 1 => zeros(3))
+p .= 2 * p
+```
+
+Instead, you should be able to use a pattern like the following, where the function `nonzeros` is used to modify the values of a sparse matrix.
+
+```@repl
+using SparseArrays
+p = spdiagm(0 => ones(4), 1 => zeros(3))
+for j in axes(p, 2)
+    for idx in nzrange(p, j)
+        i = rowvals(p)[idx]
+        nonzeros(p)[idx] = 10 * i + j # value p[i, j]
+    end
+end; p
+```

docs/src/linear_advection.md

@@ -1,7 +1,7 @@
-# Tutorial: Solution of the linear advection equation
+# [Tutorial: Solution of the linear advection equation](@id tutorial-linear-advection)
 
-This tutorial is about the efficient solution of production-destruction systems (PDS) with a large number of differential equations. 
-We will explore several ways to represent such large systems and assess their efficiency. 
+This tutorial is about the efficient solution of production-destruction systems (PDS) with a large number of differential equations.
+We will explore several ways to represent such large systems and assess their efficiency.
 
 ## Definition of the production-destruction system
 
@@ -11,7 +11,7 @@ One example of the occurrence of a PDS with a large number of equations is the s
 \partial_t u(t,x)=-a\partial_x u(t,x),\quad u(0,x)=u_0(x)
 ```
 
-with ``a>0``, ``t≥ 0``, ``x\in[0,1]`` and periodic boundary conditions. To keep things as simple as possible, we 
+with ``a>0``, ``t≥ 0``, ``x\in[0,1]`` and periodic boundary conditions. To keep things as simple as possible, we
 discretize the space domain as ``0=x_0<x_1\dots <x_{N-1}<x_N=1`` with ``x_i = i Δ x`` for ``i=0,\dots,N`` and ``Δx=1/N``. An upwind discretization of the spatial derivative yields the ODE system
 
 ```math
@@ -22,13 +22,13 @@ discretize the space domain as ``0=x_0<x_1\dots <x_{N-1}<x_N=1`` with ``x_i = i
 ```
 
 where ``u_i(t)`` is an approximation of ``u(t,x_i)`` for ``i=1,\dots, N``.
-This system can also be written as ``\partial_t \mathbf u(t)=\mathbf A\mathbf u(t)`` with ``\mathbf u(t)=(u_1(t),\dots,u_N(t))`` and 
+This system can also be written as ``\partial_t \mathbf u(t)=\mathbf A\mathbf u(t)`` with ``\mathbf u(t)=(u_1(t),\dots,u_N(t))`` and
 
 ```math
 \mathbf A= \frac{a}{Δ x}\begin{bmatrix}-1&0&\dots&0&1\\1&-1&\ddots&&0\\0&\ddots&\ddots&\ddots&\vdots\\ \vdots&\ddots&\ddots&\ddots&0\\0&\dots&0&1&-1\end{bmatrix}.
 ```
 
-In particular the matrix ``\mathbf A`` shows that there is a single production term and a single destruction term per equation. 
+In particular the matrix ``\mathbf A`` shows that there is a single production term and a single destruction term per equation.
 Furthermore, the system is conservative as ``\mathbf A`` has column sum zero.
 To be precise, the production matrix ``\mathbf P = (p_{i,j})`` of this conservative PDS is given by
 
@@ -67,11 +67,15 @@ As mentioned above, we will try different approaches to solve this PDS and compa
 
 ### Standard in-place implementation
 
+By default, we will use dense matrices to store the production terms
+and to setup/solve the linear systems arising in MPRK methods. Of course,
+this is not efficient for large and sparse systems like in this case.
+
 ```@example LinearAdvection
 using PositiveIntegrators # load ConservativePDSProblem
 
 function lin_adv_P!(P, u, p, t)
-    P .= 0.0
+    fill!(P, 0.0)
     N = length(u)
     dx = 1 / N
     P[1, N] = u[N] / dx
@@ -97,7 +101,9 @@ plot!(x, last(sol.u); label = "u")
 
 ### Using sparse matrices
 
-TODO: Some text
+To use different matrix types for the production terms and linear systems,
+you can use the keyword argument `p_prototype` of
+[`ConservativePDSProblem`](@ref) and [`PDSProblem`](@ref).
 
 ```@example LinearAdvection
 using SparseArrays
@@ -127,4 +133,30 @@ using BenchmarkTools
 
 ```@example LinearAdvection
 @benchmark solve(prob_sparse, MPRK43I(1.0, 0.5); save_everystep = false)
-```
+```
+
+By default, we use an LU factorization for the linear systems. At the time of
+writing, Julia uses [SparseArrays.jl](https://github.com/JuliaSparse/SparseArrays.jl)
+defaulting to UMFPACK from SuiteSparse in this case. However, the linear systems
+do not necessarily have the structure for which UMFPACK is optimized for. Thus,
+it is often possible to gain performance by switching to KLU instead.
+
+```@example LinearAdvection
+using LinearSolve
+@benchmark solve(prob_sparse, MPRK43I(1.0, 0.5; linsolve = KLUFactorization()); save_everystep = false)
+```
+
+
+## Package versions
+
+These results were obtained using the following versions.
+```@example LinearAdvection
+using InteractiveUtils
+versioninfo()
+println()
+
+using Pkg
+Pkg.status(["PositiveIntegrators", "SparseArrays", "KLU", "LinearSolve", "OrdinaryDiffEq"],
+           mode=PKGMODE_MANIFEST)
+nothing # hide
+```

src/PositiveIntegrators.jl

@@ -2,7 +2,8 @@ module PositiveIntegrators
 
 # 1. Load dependencies
 using LinearAlgebra: LinearAlgebra, Tridiagonal, I, diag, diagind, mul!
-using SparseArrays: SparseArrays, AbstractSparseMatrix
+using SparseArrays: SparseArrays, AbstractSparseMatrix,
+                    issparse, nonzeros, nzrange, rowvals, spdiagm
 using StaticArrays: SVector, MVector, SMatrix, StaticArray, @SVector, @SMatrix
 
 using FastBroadcast: @..