Skip to content

Add Sylvester workspaces to the buffers #150

@jlperla

Description

@jlperla

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.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions