matrixfunctions / GraphMatFun.jl

Computation graphs for matrix functions
MIT License
12 stars 0 forks source link

Sparse matrices support #83

Open jarlebring opened 5 months ago

jarlebring commented 5 months ago

Consider

julia> using LinearAlgebra, SparseArrays;
julia> (g,_)=graph_sastre_exp(3);
julia> n=500; A=spdiagm(-1 => 1:n, 1 => 1:n);

Problem 1: eval_graph does not work:

julia> eval_graph(g,A)
ERROR: UndefVarError: `SparseMatrixCSC` not defined
Stacktrace:
 [1] init_vals_eval_graph!(graph::Compgraph{Float64}, x::SparseMatrixCSC{Int64, Int64}, vals::Nothing, input::Symbol)
   @ GraphMatFun ~/.julia/packages/GraphMatFun/3g8TS/src/eval.jl:114
 [2] eval_graph(graph::Compgraph{Float64}, x::SparseMatrixCSC{Int64, Int64}; vals::Nothing, input::Symbol, output::Int64, comporder::Nothing)
   @ GraphMatFun ~/.julia/packages/GraphMatFun/3g8TS/src/eval.jl:49
 [3] eval_graph(graph::Compgraph{Float64}, x::SparseMatrixCSC{Int64, Int64})
   @ GraphMatFun ~/.julia/packages/GraphMatFun/3g8TS/src/eval.jl:41
 [4] top-level scope
   @ REPL[26]:1

Problem 2: Generated code does not work efficiently.

julia> gen_code("myfun.jl",g);
julia> include("myfun.jl");
julia> @btime dummy(A)
  328.665 ms (70 allocations: 1.40 MiB)
julia> @btime A*A*A*A;
  50.880 μs (19 allocations: 182.45 KiB)

Since g has three multiplications they should be essentially the same CPU-time. It's because of memory allocation. Not sure we can think of recycling memory slots as the current generated code is doing since the sparsity pattern will change.

jarlebring commented 5 months ago

Problem 2 might be solvable together with #76 .

jarlebring commented 5 months ago

Problem 2:

Option 1: Allocate new memory for every operation. The code will be very similar to LangMatlab()-code

Option 2: Do preallocation of memslots where with the memslots have the sparsity pattern of the final matrix. This leads to extra computation (multiplications with zero) and the efficiency will be highly pattern dependent. Not sure how to compute the sparsity pattern without computing the products.