Troubleshooting and frequently asked questions

Sparse matrices

You can use sparse matrices for the linear systems arising in PositiveIntegrators.jl, as described, e.g., in the tutorial on 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. For example, you should avoid something like

julia> using SparseArrays
julia> p = spdiagm(0 => ones(4), 1 => zeros(3))4×4 SparseArrays.SparseMatrixCSC{Float64, Int64} with 7 stored entries:
 1.0  0.0   ⋅    ⋅
  ⋅   1.0  0.0   ⋅
  ⋅    ⋅   1.0  0.0
  ⋅    ⋅    ⋅   1.0
julia> p .= 2 * p4×4 SparseArrays.SparseMatrixCSC{Float64, Int64} with 4 stored entries:
 2.0   ⋅    ⋅    ⋅
  ⋅   2.0   ⋅    ⋅
  ⋅    ⋅   2.0   ⋅
  ⋅    ⋅    ⋅   2.0

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.

julia> using SparseArrays
julia> p = spdiagm(0 => ones(4), 1 => zeros(3))4×4 SparseArrays.SparseMatrixCSC{Float64, Int64} with 7 stored entries:
 1.0  0.0   ⋅    ⋅
  ⋅   1.0  0.0   ⋅
  ⋅    ⋅   1.0  0.0
  ⋅    ⋅    ⋅   1.0
julia> 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; p4×4 SparseArrays.SparseMatrixCSC{Float64, Int64} with 7 stored entries:
 11.0  12.0    ⋅     ⋅
   ⋅   22.0  23.0    ⋅
   ⋅     ⋅   33.0  34.0
   ⋅     ⋅     ⋅   44.0