Scalar method requires element-wise evaluation of array-valued functions

[elisno]

I have a function definition that solves a system of equations by calling gesv!.
If the system is large, gesv! is called multiple times for each solution element.
This is because the current FFT call in '_wcoeff' only handles scalar valued functions.

I can make a seperate method that deals with array valued functions by storing the Laguerre-coefficients in higher dimensional arrays and applying a FFT along the first dimension of those arrays.

For simplicity, I'd start with vector valued functions.

[jlapeyre]

Can you post some example code, so I can understand the problem better ?

[elisno]

In the following system of differential equations, U(s) is determined with gesv!.
For this minimal example, we assume that the operator A is time-independent.

using InverseLaplace
using LinearAlgebra.LAPACK: gesv!
using BenchmarkTools
using Profile

const N = 64

function resolvent(s)
    A = rand(Complex{Float64},N,N)
    u₀ = rand(Complex{Float64},N)
    return gesv!(s .* one(A) - A,u₀)[1]
end


function Weeksresolvent()
    Weeksvector = Vector{Weeks{Complex{Float64}}}(undef,N)
    for i = 1:N
        Weeksvector[i] = Weeks(s -> resolvent(s)[i],32,1.0,1.0,datatype=Complex)
    end
    return Weeksvector
end


@btime Weeksresolvent()
@code_warntype Weeksresolvent()

### -  Juno Profiler
#Profile.clear()
#@profiler Weeksresolvent()

This is the output of the code above.

julia> include("gesvcalls.jl")
  1.276 s (47041 allocations: 1.01 GiB)
Body::Array{Weeks{Complex{Float64}},1}
16 1 ── %1  = $(Expr(:foreigncall, :(:jl_alloc_array_1d), Array{Weeks{Complex{Float64}},1}, svec(Any, Int64), :(:ccall), 2, Array{Weeks{Complex{Float64}},1}, 64, 64))::Array{Weeks{Complex{Float64}},1}        │╻    Type
17 │    %2  = Main.N::Core.Compiler.Const(64, false)                                                                                                                                                            │    
   │          (Base.ifelse)(true, %2, 0)                                                                                                                                                                        │╻╷   Colon
   │    %4  = (Base.slt_int)(64, 1)::Bool                                                                                                                                                                       ││╻╷╷  isempty
   └───       goto #3 if not %4                                                                                                                                                                                 ││   
   2 ──       goto #4                                                                                                                                                                                           ││   
   3 ──       goto #4                                                                                                                                                                                           ││   
   4 ┄─ %8  = φ (#2 => true, #3 => false)::Bool                                                                                                                                                                 │    
   │    %9  = φ (#3 => 1)::Int64                                                                                                                                                                                │    
   │    %10 = φ (#3 => 1)::Int64                                                                                                                                                                                │    
   │    %11 = (Base.not_int)(%8)::Bool                                                                                                                                                                          │    
   └───       goto #10 if not %11                                                                                                                                                                               │    
   5 ┄─ %13 = φ (#4 => %9, #9 => %26)::Int64                                                                                                                                                                    │    
   │    %14 = φ (#4 => %10, #9 => %27)::Int64                                                                                                                                                                   │    
18 │    %15 = %new(getfield(Main, Symbol("##3#4")){Int64}, %13)::getfield(Main, Symbol("##3#4")){Int64}                                                                                                         │    
   │    %16 = Main.Complex::Core.Compiler.Const(Complex, false)                                                                                                                                                 │    
   │    %17 = %new(NamedTuple{(:datatype,),Tuple{UnionAll}}, %16)::NamedTuple{(:datatype,),Tuple{UnionAll}}                                                                                                     ││╻╷   Type
   │    %18 = invoke Core.kwfunc(Main.Weeks::Any)::Core.Compiler.Const(getfield(Core, Symbol("#kw#Type"))(), false)                                                                                             │    
   │    %19 = invoke %18(%17::NamedTuple{(:datatype,),Tuple{UnionAll}}, Main.Weeks::Type{Weeks}, %15::Function, 32::Int64, 1.0::Float64, 1.0::Float64)::Weeks{_1} where _1                                      │    
   │          (Base.setindex!)(%1, %19, %13)                                                                                                                                                                    │    
   │    %21 = (%14 === 64)::Bool                                                                                                                                                                                ││╻    ==
   └───       goto #7 if not %21                                                                                                                                                                                ││   
   6 ──       goto #8                                                                                                                                                                                           ││   
   7 ── %24 = (Base.add_int)(%14, 1)::Int64                                                                                                                                                                     ││╻    +
   └───       goto #8                                                                                                                                                                                           │╻    iterate
   8 ┄─ %26 = φ (#7 => %24)::Int64                                                                                                                                                                              │    
   │    %27 = φ (#7 => %24)::Int64                                                                                                                                                                              │    
   │    %28 = φ (#6 => true, #7 => false)::Bool                                                                                                                                                                 │    
   │    %29 = (Base.not_int)(%28)::Bool                                                                                                                                                                         │    
   └───       goto #10 if not %29                                                                                                                                                                               │    
   9 ──       goto #5                                                                                                                                                                                           │    
20 10 ─       return %1    

Using the profiler, I see that FF0 = map(F, s) in _wcoeff is the most expensive part of the call to Weeks (the majority of the calls are to gesv!).

Why is Weeks{_1} where _1 a non-concrete type when working with complex numbers? (in %19 = invoke %18(%17::NamedTuple{(:datatype,),Tuple{UnionAll}}, Main.Weeks::Type{Weeks}, %15::Function, 32::Int64, 1.0::Float64, 1.0::Float64)::Weeks{_1} where _1)
This isn't the case when I call the Weeks' method for floating points.

I'll try to post the vector-valued method tomorrow for comparison.

[elisno]

The previous code was run using InverseLaplace from your most recent master branch.

This following code uses InverseLaplace based on #3.

using InverseLaplace
using LinearAlgebra.LAPACK: gesv!
using BenchmarkTools
using Profile

const N = 64

function resolvent(s)
    A = rand(Complex{Float64},N,N)
    u₀ = rand(Complex{Float64},N)
    return gesv!(s .* one(A) - A,u₀)[1]
end


function Weeksresolvent()
    Weeksvector = Vector{Weeks{Complex{Float64}}}(undef,N)
    for i = 1:N
        Weeksvector[i] = Weeks(s -> resolvent(s)[i],32,1.0,1.0,datatype=Complex)
    end
    return Weeksvector
end


function Weeksvectorresolvent()
    return Weeks(resolvent,32,1.0,1.0,datatype=Complex)
end

function comparescalarandvector()
    println("")
    println("Iterative")
    @btime Weeksresolvent()

    println("")
    println("Non-iterative")
    @btime Weeksvectorresolvent()
end

### -  Juno Profiler
#comparescalarandvector()
#Profile.clear()
#@profiler comparescalarandvector()

The output:

julia> comparescalarandvector()

Iterative
  1.151 s (51713 allocations: 1.02 GiB)

Non-iterative
  16.612 ms (878 allocations: 16.86 MiB)
Weeks{Complex{Float64}}(Nterms=32,sigma=1.0,b=1.0)

The relative speed difference is around 70x, It only increases for increasing N.
gesv! is still a bottleneck, but only due to multiple calls for large Nterms.

Comparing the element-wise versions, I've just now noticed an increase in the number of allocations.