Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

improve performance with sparse matrices #101

Merged
merged 12 commits into from
Jul 18, 2024
6 changes: 6 additions & 0 deletions docs/Project.toml
Original file line number Diff line number Diff line change
@@ -1,14 +1,20 @@
[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"
PositiveIntegrators = "d1b20bf0-b083-4985-a874-dc5121669aa5"
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"
1 change: 1 addition & 0 deletions docs/make.jl
Original file line number Diff line number Diff line change
@@ -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",
30 changes: 30 additions & 0 deletions docs/src/faq.md
Original file line number Diff line number Diff line change
@@ -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.
ranocha marked this conversation as resolved.
Show resolved Hide resolved

```@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
```
50 changes: 41 additions & 9 deletions docs/src/linear_advection.md
Original file line number Diff line number Diff line change
@@ -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.

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are you saying that this is the reason for the large number of allocations we see in the linear advection tutorial?

There's also something else we should look at. When we initialize the linear solvers we use something like

linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
                        assumptions = LinearSolve.OperatorAssumptions(true))

The statement alias_A=true allows the linear solver to reuse the memory of the system matrix and modify it. The idea behind this was that in an MPRK scheme we no longer need the system matrix once the system is solved, so the linear solver can change it at will.

For a sparse matrix this modification of the system matrix could become problematic, since the sparsity pattern of L and U does not match that of A, especially in the linear advection tutorial, where A is a banded matrix with an entry in the upper right corner, for which in genaral L or U could become a full triangular matrix.

All I want to say is that we should try alias_A=false and check if we see any difference.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Shall we use alias_A = !issparse(the_matrix_we_put_in) for the general case?

```@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
```
3 changes: 2 additions & 1 deletion src/PositiveIntegrators.jl
Original file line number Diff line number Diff line change
@@ -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: @..
Loading