From 602ebab4a5b7d0e582aed545c42480ef6d9f8fea Mon Sep 17 00:00:00 2001 From: Thore Kockerols Date: Thu, 5 Jun 2025 10:53:17 +0100 Subject: [PATCH 1/2] Integrate cached block decomposition --- benchmark/benchmarks.jl | 14 +++-- src/MacroModelling.jl | 59 +++++++++++++------- src/filter/inversion.jl | 9 +-- src/filter/kalman.jl | 3 +- src/macros.jl | 1 + src/moments.jl | 18 +++--- src/perturbation.jl | 95 ++++++++++++++++++++++++-------- src/structures.jl | 1 + test/test_standalone_function.jl | 2 +- 9 files changed, 142 insertions(+), 60 deletions(-) diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl index 78091152b..4760159c4 100644 --- a/benchmark/benchmarks.jl +++ b/benchmark/benchmarks.jl @@ -43,13 +43,19 @@ function run_benchmarks!(𝓂::β„³, SUITE::BenchmarkGroup) SUITE[𝓂.model_name]["qme"] = BenchmarkGroup() - sol, qme_sol, solved = calculate_first_order_solution(βˆ‡β‚; T = 𝓂.timings, opts = merge_calculation_options(quadratic_matrix_equation_algorithm = :schur)) + sol, qme_sol, solved = calculate_first_order_solution(βˆ‡β‚; T = 𝓂.timings, + opts = merge_calculation_options(quadratic_matrix_equation_algorithm = :schur), + decomposition = 𝓂.solution.perturbation.first_order_block_decomposition) clear_solution_caches!(𝓂, :first_order) - SUITE[𝓂.model_name]["qme"]["schur"] = @benchmarkable calculate_first_order_solution($βˆ‡β‚; T = $𝓂.timings, opts = merge_calculation_options(quadratic_matrix_equation_algorithm = :schur)) setup = clear_solution_caches!($𝓂, :first_order) - - SUITE[𝓂.model_name]["qme"]["doubling"] = @benchmarkable calculate_first_order_solution($βˆ‡β‚; T = $𝓂.timings, opts = merge_calculation_options(quadratic_matrix_equation_algorithm = :doubling)) setup = clear_solution_caches!($𝓂, :first_order) + SUITE[𝓂.model_name]["qme"]["schur"] = @benchmarkable calculate_first_order_solution($βˆ‡β‚; T = $𝓂.timings, + opts = merge_calculation_options(quadratic_matrix_equation_algorithm = :schur), + decomposition = $𝓂.solution.perturbation.first_order_block_decomposition) setup = clear_solution_caches!($𝓂, :first_order) + + SUITE[𝓂.model_name]["qme"]["doubling"] = @benchmarkable calculate_first_order_solution($βˆ‡β‚; T = $𝓂.timings, + opts = merge_calculation_options(quadratic_matrix_equation_algorithm = :doubling), + decomposition = $𝓂.solution.perturbation.first_order_block_decomposition) setup = clear_solution_caches!($𝓂, :first_order) A = @views sol[:, 1:𝓂.timings.nPast_not_future_and_mixed] * β„’.diagm(ones(𝓂.timings.nVars))[𝓂.timings.past_not_future_and_mixed_idx,:] diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index 7d542fef1..f02966dbc 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -5042,10 +5042,11 @@ function calculate_second_order_stochastic_steady_state(parameters::Vector{M}, # @timeit_debug timer "Calculate first order solution" begin - 𝐒₁, qme_sol, solved = calculate_first_order_solution(βˆ‡β‚; - T = 𝓂.timings, + 𝐒₁, qme_sol, solved = calculate_first_order_solution(βˆ‡β‚; + T = 𝓂.timings, opts = opts, - initial_guess = 𝓂.solution.perturbation.qme_solution) + initial_guess = 𝓂.solution.perturbation.qme_solution, + decomposition = 𝓂.solution.perturbation.first_order_block_decomposition) if solved 𝓂.solution.perturbation.qme_solution = qme_sol end @@ -5369,10 +5370,11 @@ function calculate_third_order_stochastic_steady_state( parameters::Vector{M}, βˆ‡β‚ = calculate_jacobian(parameters, SS_and_pars, 𝓂)# |> Matrix - 𝐒₁, qme_sol, solved = calculate_first_order_solution(βˆ‡β‚; - T = 𝓂.timings, + 𝐒₁, qme_sol, solved = calculate_first_order_solution(βˆ‡β‚; + T = 𝓂.timings, opts = opts, - initial_guess = 𝓂.solution.perturbation.qme_solution) + initial_guess = 𝓂.solution.perturbation.qme_solution, + decomposition = 𝓂.solution.perturbation.first_order_block_decomposition) if solved 𝓂.solution.perturbation.qme_solution = qme_sol end @@ -5754,11 +5756,11 @@ function solve!(𝓂::β„³; # @timeit_debug timer "Calculate first order solution" begin - S₁, qme_sol, solved = calculate_first_order_solution(βˆ‡β‚; - T = 𝓂.timings, + S₁, qme_sol, solved = calculate_first_order_solution(βˆ‡β‚; + T = 𝓂.timings, opts = opts, - initial_guess = 𝓂.solution.perturbation.qme_solution) - + initial_guess = 𝓂.solution.perturbation.qme_solution, + decomposition = 𝓂.solution.perturbation.first_order_block_decomposition) if solved 𝓂.solution.perturbation.qme_solution = qme_sol end # end # timeit_debug @@ -5776,11 +5778,12 @@ function solve!(𝓂::β„³; βˆ‡Μ‚β‚ = calculate_jacobian(𝓂.parameter_values, SS_and_pars, 𝓂)# |> Matrix - Ŝ₁, qme_sol, solved = calculate_first_order_solution(βˆ‡Μ‚β‚; - T = 𝓂.timings, + Ŝ₁, qme_sol, solved = calculate_first_order_solution(βˆ‡Μ‚β‚; + T = 𝓂.timings, opts = opts, - initial_guess = 𝓂.solution.perturbation.qme_solution) - + initial_guess = 𝓂.solution.perturbation.qme_solution, + decomposition = 𝓂.solution.perturbation.first_order_block_decomposition) + if solved 𝓂.solution.perturbation.qme_solution = qme_sol end write_parameters_input!(𝓂, :activeα΅’α΅‡αΆœshocks => 0, verbose = false) @@ -6622,6 +6625,23 @@ function write_functions_mapping!(𝓂::β„³, max_perturbation_order::Int; 𝓂.jacobian_SS_and_pars = buffer_SS_and_pars, func_βˆ‡β‚_SS_and_pars + # precompute block decomposition for first order solution + rows, cols, _ = findnz(βˆ‡β‚_dyn) + pattern_full = sparse(rows, cols, ones(Int, length(rows)), size(βˆ‡β‚_dyn,1), size(βˆ‡β‚_dyn,2)) + + T = 𝓂.timings + dynIndex = T.nPresent_only+1:T.nVars + comb = union(T.future_not_past_and_mixed_idx, T.past_not_future_idx) + sort!(comb) + f_in = indexin(T.future_not_past_and_mixed_idx, comb) + Ir = β„’.I(length(comb)) + + Aplus_pattern = pattern_full[dynIndex,1:T.nFuture_not_past_and_mixed] * Ir[f_in,:] + A0_pattern = pattern_full[dynIndex,T.nFuture_not_past_and_mixed .+ range(1,T.nVars)][:,comb] + block_pattern = (Aplus_pattern .+ A0_pattern) .> 0 + Q, P, R, _, _ = BlockTriangularForm.order(block_pattern) + 𝓂.solution.perturbation.first_order_block_decomposition = (Vector{Int}(Q), Vector{Int}(P), Vector{Int}(R)) + @@ -8358,11 +8378,12 @@ function get_relevant_steady_state_and_state_update(::Val{:first_order}, βˆ‡β‚ = calculate_jacobian(parameter_values, SS_and_pars, 𝓂) # , timer = timer)# |> Matrix - 𝐒₁, qme_sol, solved = calculate_first_order_solution(βˆ‡β‚; - T = TT, - # timer = timer, - initial_guess = 𝓂.solution.perturbation.qme_solution, - opts = opts) + 𝐒₁, qme_sol, solved = calculate_first_order_solution(βˆ‡β‚; + T = TT, + # timer = timer, + initial_guess = 𝓂.solution.perturbation.qme_solution, + opts = opts, + decomposition = 𝓂.solution.perturbation.first_order_block_decomposition) if solved 𝓂.solution.perturbation.qme_solution = qme_sol end diff --git a/src/filter/inversion.jl b/src/filter/inversion.jl index 7427b8d4b..1ed26014d 100644 --- a/src/filter/inversion.jl +++ b/src/filter/inversion.jl @@ -3457,10 +3457,11 @@ function filter_data_with_model(𝓂::β„³, βˆ‡β‚ = calculate_jacobian(𝓂.parameter_values, SS_and_pars, 𝓂)# |> Matrix - 𝐒₁, qme_sol, solved = calculate_first_order_solution(βˆ‡β‚; - T = T, - initial_guess = 𝓂.solution.perturbation.qme_solution, - opts = opts) + 𝐒₁, qme_sol, solved = calculate_first_order_solution(βˆ‡β‚; + T = T, + initial_guess = 𝓂.solution.perturbation.qme_solution, + opts = opts, + decomposition = 𝓂.solution.perturbation.first_order_block_decomposition) if solved 𝓂.solution.perturbation.qme_solution = qme_sol end diff --git a/src/filter/kalman.jl b/src/filter/kalman.jl index 767810b00..1bf749fef 100644 --- a/src/filter/kalman.jl +++ b/src/filter/kalman.jl @@ -601,7 +601,8 @@ function filter_and_smooth(𝓂::β„³, βˆ‡β‚ = calculate_jacobian(parameters, SS_and_pars, 𝓂)# |> Matrix - sol, qme_sol, solved = calculate_first_order_solution(βˆ‡β‚; T = 𝓂.timings, opts = opts) + sol, qme_sol, solved = calculate_first_order_solution(βˆ‡β‚; T = 𝓂.timings, opts = opts, + decomposition = 𝓂.solution.perturbation.first_order_block_decomposition) if solved 𝓂.solution.perturbation.qme_solution = qme_sol end diff --git a/src/macros.jl b/src/macros.jl index 01d2ed004..a47b7a6ba 100644 --- a/src/macros.jl +++ b/src/macros.jl @@ -949,6 +949,7 @@ macro model(𝓂,ex...) third_order_perturbation_solution([], (x,y)->nothing, (x,y)->nothing), third_order_perturbation_solution([], (x,y)->nothing, (x,y)->nothing), zeros(0,0), # 1st order sol + (Int[], Int[], Int[]), SparseMatrixCSC{Float64, Int64}(β„’.I,0,0), # 2nd order sol SparseMatrixCSC{Float64, Int64}(β„’.I,0,0), # 3rd order sol auxilliary_indices(Int[],Int[],Int[],Int[],Int[]), diff --git a/src/moments.jl b/src/moments.jl index 4d5982988..8f8c38b9b 100644 --- a/src/moments.jl +++ b/src/moments.jl @@ -11,10 +11,11 @@ function calculate_covariance(parameters::Vector{R}, βˆ‡β‚ = calculate_jacobian(parameters, SS_and_pars, 𝓂) - sol, qme_sol, solved = calculate_first_order_solution(βˆ‡β‚; - T = 𝓂.timings, - initial_guess = 𝓂.solution.perturbation.qme_solution, - opts = opts) + sol, qme_sol, solved = calculate_first_order_solution(βˆ‡β‚; + T = 𝓂.timings, + initial_guess = 𝓂.solution.perturbation.qme_solution, + opts = opts, + decomposition = 𝓂.solution.perturbation.first_order_block_decomposition) if solved 𝓂.solution.perturbation.qme_solution = qme_sol end @@ -56,10 +57,11 @@ function calculate_mean(parameters::Vector{T}, else βˆ‡β‚ = calculate_jacobian(parameters, SS_and_pars, 𝓂)# |> Matrix - 𝐒₁, qme_sol, solved = calculate_first_order_solution(βˆ‡β‚; - T = 𝓂.timings, - initial_guess = 𝓂.solution.perturbation.qme_solution, - opts = opts) + 𝐒₁, qme_sol, solved = calculate_first_order_solution(βˆ‡β‚; + T = 𝓂.timings, + initial_guess = 𝓂.solution.perturbation.qme_solution, + opts = opts, + decomposition = 𝓂.solution.perturbation.first_order_block_decomposition) if !solved mean_of_variables = SS_and_pars[1:𝓂.timings.nVars] diff --git a/src/perturbation.jl b/src/perturbation.jl index 5ed5c3592..9289722ed 100644 --- a/src/perturbation.jl +++ b/src/perturbation.jl @@ -1,9 +1,62 @@ @stable default_mode = "disable" begin -function calculate_first_order_solution(βˆ‡β‚::Matrix{R}; - T::timings, +function _solve_first_order_blocks(Ap, A0, Am, Q, P, R, T, initial_guess, opts) + n_blocks = length(R) - 1 + + Ap_perm = @views Ap[Q, P] + A0_perm = @views A0[Q, P] + Am_perm = @views Am[Q, P] + + ig_perm = if length(initial_guess) > 0 + @views initial_guess[Q, P] + else + zeros(eltype(Ap), 0, 0) + end + + sol = zeros(eltype(Ap), size(Ap)) + + for i in 1:n_blocks + rng = R[i]:(R[i+1]-1) + Apb = @views Ap_perm[rng, rng] + A0b = @views A0_perm[rng, rng] + Amb = @views Am_perm[rng, rng] + igb = if length(ig_perm) > 0 + @views ig_perm[rng, rng] + else + zeros(eltype(Ap), 0, 0) + end + + if size(Apb,1) == 1 && abs(Apb[1]) < opts.tol.qme_tol + sol[Q[rng], P[rng]] .= -Amb[1] / A0b[1] + elseif size(Apb,1) == 1 + a = Apb[1]; b = A0b[1]; c = Amb[1] + disc = b^2 - 4a*c + root1 = (-b + sqrt(disc)) / (2a) + root2 = (-b - sqrt(disc)) / (2a) + sol[Q[rng], P[rng]] .= abs(root1) < 1 ? root1 : root2 + else + Xblock, solved = solve_quadratic_matrix_equation(Apb, A0b, Amb, T, + initial_guess = igb, + quadratic_matrix_equation_algorithm = opts.quadratic_matrix_equation_algorithm, + tol = opts.tol.qme_tol, + acceptance_tol = opts.tol.qme_acceptance_tol, + verbose = opts.verbose) + if !solved + return sol, false + end + sol[Q[rng], P[rng]] .= Xblock + end + end + + return sol, true +end + + +function calculate_first_order_solution(βˆ‡β‚::Matrix{R}; + T::timings, opts::CalculationOptions = merge_calculation_options(), - initial_guess::AbstractMatrix{R} = zeros(0,0))::Tuple{Matrix{R}, Matrix{R}, Bool} where R <: AbstractFloat + initial_guess::AbstractMatrix{R} = zeros(0,0), + decomposition::Tuple{Vector{Int},Vector{Int},Vector{Int}})::Tuple{Matrix{R}, Matrix{R}, Bool} where R <: AbstractFloat # @timeit_debug timer "Calculate 1st order solution" begin # @timeit_debug timer "Preprocessing" begin @@ -43,12 +96,9 @@ function calculate_first_order_solution(βˆ‡β‚::Matrix{R}; # end # timeit_debug # @timeit_debug timer "Quadratic matrix equation solve" begin - sol, solved = solve_quadratic_matrix_equation(AΜƒβ‚Š, AΜƒβ‚€, AΜƒβ‚‹, T, - initial_guess = initial_guess, - quadratic_matrix_equation_algorithm = opts.quadratic_matrix_equation_algorithm, - tol = opts.tol.qme_tol, - acceptance_tol = opts.tol.qme_acceptance_tol, - verbose = opts.verbose) + Q, P, R = decomposition + + sol, solved = _solve_first_order_blocks(AΜƒβ‚Š, AΜƒβ‚€, AΜƒβ‚‹, Q, P, R, T, initial_guess, opts) if !solved if opts.verbose println("Quadratic matrix equation solution failed.") end @@ -117,11 +167,12 @@ end end # dispatch_doctor -function rrule(::typeof(calculate_first_order_solution), +function rrule(::typeof(calculate_first_order_solution), βˆ‡β‚::Matrix{R}; - T::timings, + T::timings, opts::CalculationOptions = merge_calculation_options(), - initial_guess::AbstractMatrix{R} = zeros(0,0)) where R <: AbstractFloat + initial_guess::AbstractMatrix{R} = zeros(0,0), + decomposition::Tuple{Vector{Int},Vector{Int},Vector{Int}}) where R <: AbstractFloat # Forward pass to compute the output and intermediate values needed for the backward pass # @timeit_debug timer "Calculate 1st order solution" begin # @timeit_debug timer "Preprocessing" begin @@ -162,12 +213,9 @@ function rrule(::typeof(calculate_first_order_solution), # end # timeit_debug # @timeit_debug timer "Quadratic matrix equation solve" begin - sol, solved = solve_quadratic_matrix_equation(AΜƒβ‚Š, AΜƒβ‚€, AΜƒβ‚‹, T, - initial_guess = initial_guess, - quadratic_matrix_equation_algorithm = opts.quadratic_matrix_equation_algorithm, - tol = opts.tol.qme_tol, - acceptance_tol = opts.tol.qme_acceptance_tol, - verbose = opts.verbose) + Q, P, R = decomposition + + sol, solved = _solve_first_order_blocks(AΜƒβ‚Š, AΜƒβ‚€, AΜƒβ‚‹, Q, P, R, T, initial_guess, opts) if !solved return (zeros(T.nVars,T.nPast_not_future_and_mixed + T.nExo), sol, false), x -> NoTangent(), NoTangent(), NoTangent() @@ -276,10 +324,11 @@ end @stable default_mode = "disable" begin -function calculate_first_order_solution(βˆ‡β‚::Matrix{β„±.Dual{Z,S,N}}; - T::timings, +function calculate_first_order_solution(βˆ‡β‚::Matrix{β„±.Dual{Z,S,N}}; + T::timings, opts::CalculationOptions = merge_calculation_options(), - initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0))::Tuple{Matrix{β„±.Dual{Z,S,N}}, Matrix{Float64}, Bool} where {Z,S,N} + initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), + decomposition::Tuple{Vector{Int},Vector{Int},Vector{Int}})::Tuple{Matrix{β„±.Dual{Z,S,N}}, Matrix{Float64}, Bool} where {Z,S,N} βˆ‡Μ‚β‚ = β„±.value.(βˆ‡β‚) expand = [β„’.I(T.nVars)[T.future_not_past_and_mixed_idx,:], β„’.I(T.nVars)[T.past_not_future_and_mixed_idx,:]] @@ -287,7 +336,7 @@ function calculate_first_order_solution(βˆ‡β‚::Matrix{β„±.Dual{Z,S,N}}; A = βˆ‡Μ‚β‚[:,1:T.nFuture_not_past_and_mixed] * expand[1] B = βˆ‡Μ‚β‚[:,T.nFuture_not_past_and_mixed .+ range(1,T.nVars)] - 𝐒₁, qme_sol, solved = calculate_first_order_solution(βˆ‡Μ‚β‚; T = T, opts = opts, initial_guess = initial_guess) + 𝐒₁, qme_sol, solved = calculate_first_order_solution(βˆ‡Μ‚β‚; T = T, opts = opts, initial_guess = initial_guess, decomposition = decomposition) if !solved return βˆ‡β‚, qme_sol, false @@ -363,7 +412,7 @@ function calculate_first_order_solution(βˆ‡β‚::Matrix{β„±.Dual{Z,S,N}}; B = -((βˆ‡β‚Š * x * Jm + βˆ‡β‚€) \ βˆ‡β‚‘) return hcat(x, B), qme_sol, solved -end +end diff --git a/src/structures.jl b/src/structures.jl index 65eacf772..b23303aac 100644 --- a/src/structures.jl +++ b/src/structures.jl @@ -217,6 +217,7 @@ mutable struct perturbation third_order::third_order_perturbation_solution pruned_third_order::third_order_perturbation_solution qme_solution::Matrix{Float64} + first_order_block_decomposition::Tuple{Vector{Int},Vector{Int},Vector{Int}} second_order_solution::AbstractMatrix{Float64} third_order_solution::AbstractMatrix{Float64} auxilliary_indices::auxilliary_indices diff --git a/test/test_standalone_function.jl b/test/test_standalone_function.jl index 146c2de2c..89e05c6c0 100644 --- a/test/test_standalone_function.jl +++ b/test/test_standalone_function.jl @@ -72,7 +72,7 @@ get_irf(RBC_CME, algorithm = :pruned_second_order) T = timings([:R, :y], [:Pi, :c], [:k, :z_delta], [:A], [:A, :Pi, :c], [:A, :k, :z_delta], [:A, :Pi, :c, :k, :z_delta], [:A], [:k, :z_delta], [:A], [:delta_eps, :eps_z], [:A, :Pi, :R, :c, :k, :y, :z_delta], Symbol[], Symbol[], 2, 1, 3, 3, 5, 7, 2, [3, 6], [1, 2, 4, 5, 7], [1, 2, 4], [2, 3], [1, 5, 7], [1], [1], [5, 7], [5, 6, 1, 7, 3, 2, 4], [3, 4, 5, 1, 2]) -first_order_solution, qme_sol, solved = calculate_first_order_solution(βˆ‡β‚; T = T)# |> Matrix{Float32} +first_order_solution, qme_sol, solved = calculate_first_order_solution(βˆ‡β‚; T = T, decomposition = (Int[], Int[], Int[])) second_order_solution, solved2 = calculate_second_order_solution(βˆ‡β‚, βˆ‡β‚‚, first_order_solution, RBC_CME.solution.perturbation.second_order_auxilliary_matrices, From 941adec0134edd0709da2beb8eb3cb3e042a0e24 Mon Sep 17 00:00:00 2001 From: Thore Kockerols Date: Thu, 5 Jun 2025 10:55:05 +0000 Subject: [PATCH 2/2] fix naming in order to avoid conflicts with other vars --- src/perturbation.jl | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/src/perturbation.jl b/src/perturbation.jl index 9289722ed..17154f004 100644 --- a/src/perturbation.jl +++ b/src/perturbation.jl @@ -1,14 +1,14 @@ @stable default_mode = "disable" begin -function _solve_first_order_blocks(Ap, A0, Am, Q, P, R, T, initial_guess, opts) - n_blocks = length(R) - 1 +function _solve_first_order_blocks(Ap, A0, Am, 𝒬, 𝒫, β„›, T, initial_guess, opts) + n_blocks = length(β„›) - 1 - Ap_perm = @views Ap[Q, P] - A0_perm = @views A0[Q, P] - Am_perm = @views Am[Q, P] + Ap_perm = @views Ap[𝒬, 𝒫] + A0_perm = @views A0[𝒬, 𝒫] + Am_perm = @views Am[𝒬, 𝒫] ig_perm = if length(initial_guess) > 0 - @views initial_guess[Q, P] + @views initial_guess[𝒬, 𝒫] else zeros(eltype(Ap), 0, 0) end @@ -16,7 +16,7 @@ function _solve_first_order_blocks(Ap, A0, Am, Q, P, R, T, initial_guess, opts) sol = zeros(eltype(Ap), size(Ap)) for i in 1:n_blocks - rng = R[i]:(R[i+1]-1) + rng = β„›[i]:(β„›[i+1]-1) Apb = @views Ap_perm[rng, rng] A0b = @views A0_perm[rng, rng] Amb = @views Am_perm[rng, rng] @@ -27,13 +27,13 @@ function _solve_first_order_blocks(Ap, A0, Am, Q, P, R, T, initial_guess, opts) end if size(Apb,1) == 1 && abs(Apb[1]) < opts.tol.qme_tol - sol[Q[rng], P[rng]] .= -Amb[1] / A0b[1] + sol[𝒬[rng], 𝒫[rng]] .= -Amb[1] / A0b[1] elseif size(Apb,1) == 1 a = Apb[1]; b = A0b[1]; c = Amb[1] disc = b^2 - 4a*c root1 = (-b + sqrt(disc)) / (2a) root2 = (-b - sqrt(disc)) / (2a) - sol[Q[rng], P[rng]] .= abs(root1) < 1 ? root1 : root2 + sol[𝒬[rng], 𝒫[rng]] .= abs(root1) < 1 ? root1 : root2 else Xblock, solved = solve_quadratic_matrix_equation(Apb, A0b, Amb, T, initial_guess = igb, @@ -44,7 +44,7 @@ function _solve_first_order_blocks(Ap, A0, Am, Q, P, R, T, initial_guess, opts) if !solved return sol, false end - sol[Q[rng], P[rng]] .= Xblock + sol[𝒬[rng], 𝒫[rng]] .= Xblock end end @@ -96,9 +96,9 @@ function calculate_first_order_solution(βˆ‡β‚::Matrix{R}; # end # timeit_debug # @timeit_debug timer "Quadratic matrix equation solve" begin - Q, P, R = decomposition + 𝒬, 𝒫, β„› = decomposition - sol, solved = _solve_first_order_blocks(AΜƒβ‚Š, AΜƒβ‚€, AΜƒβ‚‹, Q, P, R, T, initial_guess, opts) + sol, solved = _solve_first_order_blocks(AΜƒβ‚Š, AΜƒβ‚€, AΜƒβ‚‹, 𝒬, 𝒫, β„›, T, initial_guess, opts) if !solved if opts.verbose println("Quadratic matrix equation solution failed.") end @@ -213,9 +213,9 @@ function rrule(::typeof(calculate_first_order_solution), # end # timeit_debug # @timeit_debug timer "Quadratic matrix equation solve" begin - Q, P, R = decomposition + 𝒬, 𝒫, β„› = decomposition - sol, solved = _solve_first_order_blocks(AΜƒβ‚Š, AΜƒβ‚€, AΜƒβ‚‹, Q, P, R, T, initial_guess, opts) + sol, solved = _solve_first_order_blocks(AΜƒβ‚Š, AΜƒβ‚€, AΜƒβ‚‹, 𝒬, 𝒫, β„›, T, initial_guess, opts) if !solved return (zeros(T.nVars,T.nPast_not_future_and_mixed + T.nExo), sol, false), x -> NoTangent(), NoTangent(), NoTangent() @@ -405,7 +405,7 @@ function calculate_first_order_solution(βˆ‡β‚::Matrix{β„±.Dual{Z,S,N}}; Jm = @view(β„’.diagm(ones(S,T.nVars))[T.past_not_future_and_mixed_idx,:]) - βˆ‡β‚Š = βˆ‡β‚[:,1:T.nFuture_not_past_and_mixed] * β„’.diagm(ones(S,T.nVars))[T.future_not_past_and_mixed_idx,:] + βˆ‡β‚Š = βˆ‡β‚[:,1:T.nFuture_not_past_and_mixed]βˆ‡β‚Š = βˆ‡β‚ βˆ‡β‚€ = βˆ‡β‚[:,T.nFuture_not_past_and_mixed .+ range(1,T.nVars)] βˆ‡β‚‘ = βˆ‡β‚[:,(T.nFuture_not_past_and_mixed + T.nVars + T.nPast_not_future_and_mixed + 1):end]