HighDimensionalEconLab / DifferentiableStateSpaceModels.jl

MIT License
46 stars 1 forks source link

Add Sylvester workspaces to the buffers #150

Open jlperla opened 2 years ago

jlperla commented 2 years ago

If #127 is successful, then consider caching the workspaces.

In particular, the workspace calls in https://github.com/HighDimensionalEconLab/DifferentiableStateSpaceModels.jl/blob/v0.4.19/src/generate_perturbation_derivatives.jl#L280 and https://github.com/HighDimensionalEconLab/DifferentiableStateSpaceModels.jl/blob/v0.4.19/src/generate_perturbation_derivatives.jl#L143 and https://github.com/HighDimensionalEconLab/DifferentiableStateSpaceModels.jl/blob/v0.4.19/src/generate_perturbation.jl#L338 could themselves be cached.

There are only 2 variations of parameters, I believe.

I would check with Michiel to see if they need to be zero'd out each time. If not, then move them to https://github.com/HighDimensionalEconLab/DifferentiableStateSpaceModels.jl/blob/v0.4.19/src/types.jl#L59-L157

Basically, you would want to add in a generic type parameter for the workspace buffer, then call that constructor for GeneralizedSylvester, then just use the workspace each time.

For example, maybe the firstordersolverbuffers becomes

struct FirstOrderSolverBuffers{WType}
    A::Matrix{Complex{Float64}}
    B::Matrix{Complex{Float64}}
    Z::Matrix{Float64}  # real version, transposed relative to the schur
    Z_ll::Matrix{Float64}
    S_bb::UpperTriangular{Float64,Matrix{Float64}}
    T_bb::UpperTriangular{Float64,Matrix{Float64}}
    sylvester_ws_2::WType
end
function FirstOrderSolverBuffers(n_y, n_x, n_p_d, n_ϵ, n_z)
    return FirstOrderSolverBuffers(zeros(Complex{Float64}, n_x + n_y, n_x + n_y),
                                   zeros(Complex{Float64}, n_x + n_y, n_x + n_y),
                                   zeros(n_x + n_y, n_x + n_y),
                                   zeros(n_y, n_y),
                                   UpperTriangular(zeros(n_x, n_x)),
                                   UpperTriangular(zeros(n_x, n_x)),
                                   GeneralizedSylvesterWs(n, n, n_x, 2))
end

Then replace https://github.com/HighDimensionalEconLab/DifferentiableStateSpaceModels.jl/blob/v0.4.19/src/generate_perturbation.jl#L336-L339 with just

            generalized_sylvester_solver!(buff.A, buff.C, c.h_x, buff.E, 2, buff.sylvester_ws_2)

Or something along the lines. If it turns out that we need to zero out the workspace buffers between calls, then add in some calls to do that somewhere.

The same thing can be done for the other workspace uses.