SKopecz / PositiveIntegrators.jl

A Julia library of positivity-preserving time integration methods
https://skopecz.github.io/PositiveIntegrators.jl/
MIT License
13 stars 4 forks source link

Specializations for `SparseMatrix` #98

Closed SKopecz closed 2 months ago

SKopecz commented 2 months ago

The benchmark results in the "Linear Advection Tutorial" show that using a SparseMatrix is not much faster compared to a dense matrix. Moreover, the number of allocations is increased tremendously. See https://skopecz.github.io/PositiveIntegrators.jl/previews/PR96/linear_advection/ or the local results for N=100 below.

julia> @benchmark solve(prob, MPRK43I(1.0, 0.5); save_everystep = false)
BenchmarkTools.Trial: 30 samples with 1 evaluation.
 Range (min … max):  160.886 ms … 189.764 ms  ┊ GC (min … max): 0.91% … 1.00%
 Time  (median):     169.169 ms               ┊ GC (median):    0.96%
 Time  (mean ± σ):   171.651 ms ±   7.304 ms  ┊ GC (mean ± σ):  1.01% ± 0.33%

        ▁   ▄  ▁▁  ▁                  █     ▁                    
  ▆▁▁▆▆▁█▆▁▁█▁▁██▁▆█▆▁▁▁▁▆▆▁▁▁▁▁▁▁▁▁▁▆█▆▁▆▁▁█▆▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▆ ▁
  161 ms           Histogram: frequency by time          190 ms <

 Memory estimate: 83.61 MiB, allocs estimate: 3282.

julia> @benchmark solve(prob_sparse, MPRK43I(1.0, 0.5); save_everystep = false)
BenchmarkTools.Trial: 24 samples with 1 evaluation.
 Range (min … max):  209.740 ms … 226.565 ms  ┊ GC (min … max): 0.53% … 0.48%
 Time  (median):     214.642 ms               ┊ GC (median):    0.54%
 Time  (mean ± σ):   215.483 ms ±   3.314 ms  ┊ GC (mean ± σ):  0.75% ± 0.39%

             ▃ █  █▃             █                               
  ▇▁▁▁▁▁▇▁▁▁▁█▁█▇▁██▇▁▇▇▇▁▁▁▇▁▁▁▇█▁▁▇▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▇ ▁
  210 ms           Histogram: frequency by time          227 ms <

 Memory estimate: 91.61 MiB, allocs estimate: 43393.

To see the efficiency of sparse matrices we need to add specializations for general sparse matrices.

Another big problem related to the linear advection tutorial is that p_prototype in this case does not include the diagonal of the matrix in its sparsity pattern. Using this sparsity pattern to build the MPRK matrix requires adding values to the diagonal, which is inefficient. To avoid this problem we need to make sure, that the diagonal is included in the sparsity pattern of p_prototype.

SKopecz commented 2 months ago

Is there an efficient way to add the diagonal to the sparsity pattern of a given matrix @ranocha? Or shall we just add a sparse diagonal matrix to p_protoype?

ranocha commented 2 months ago

Is there an efficient way to add the diagonal to the sparsity pattern of a given matrix @ranocha? Or shall we just add a sparse diagonal matrix to p_protoype?

This is only required once when setting up the prototype, isn't it? I would try to add I from LinearAlgebra or an appropriate spdiagm