Commit e39adc5e authored by MichelJuillard's avatar MichelJuillard
Browse files

separate extended_lyapd_core!

parent 6836fcb5
......@@ -4,5 +4,6 @@ include("linear_rational_expectations.jl")
export LinearRationalExpectationsWs, LinearRationalExpectationsResults, first_order_solver!
include("extended_lyapunov.jl")
export LyapdWs, extended_lyapd!
export LyapdWs, extended_lyapd!, extended_lyapd_core!
end
......@@ -133,42 +133,8 @@ function extended_lyapd!(Σ, A, B, ws)
end
end
fill!(Σ, 0.0)
row = n
while row >= rowu
if row == 1 || ws.AA[row, row - 1] == 0
solve_one_row!(Σ, ws.AA, ws.BB, n, row, ws)
vB = view(ws.BB, 1:row - 1, 1:row - 1)
vtemp = view(ws.temp1, 1:row - 1, 1)
= view(Σ, 1:row -1 , row)
vA = view(ws.AA, 1:row - 1, 1:row - 1)
mul!(vtemp, vA, )
vA = view(ws.AA, 1:row - 1, row)
mul!(vB, vtemp, vA', 1.0, 1.0)
= view(Σ, 1:row, row)
vA = view(ws.AA, 1:row - 1, 1:row)
mul!(vtemp, vA, )
vA = view(ws.AA, 1:row - 1, row)
vB .= vB .+ vA.*vtemp'
row -= 1
else
solve_two_rows!(Σ, ws.AA, ws.BB, n, row, ws)
vB = view(ws.BB, 1:row - 2, 1:row - 2)
vtemp = view(ws.temp1, 1:row - 2, 1:2)
= view(Σ, 1:row - 2, row - 1:row)
vA = view(ws.AA, 1:row - 2, 1:row - 2)
mul!(vtemp, vA, )
vA = view(ws.AA, 1:row - 2, row - 1:row)
mul!(vB, vtemp, vA', 1.0, 1.0)
= view(Σ, 1:row, row - 1:row)
vA = view(ws.AA, 1:row - 2, 1:row)
mul!(vtemp, vA, )
vA = view(ws.AA, 1:row - 2, row - 1:row)
mul!(vB, vA, vtemp', 1.0, 1.0)
row -= 2
end
end
extended_lyapd_core!(Σ, ws.AA, ws.BB, ws)
mul!(ws.temp1, ws.dgees_ws.vs, Σ)
mul!(Σ, ws.temp1, ws.dgees_ws.vs')
......@@ -191,3 +157,42 @@ function extended_lyapd!(Σ, A, B, ws)
end
end
function extended_lyapd_core!(Σ, A, B, ws)
fill!(Σ, 0.0)
n = size(A, 1)
row = n
while row >= 1
if row == 1 || A[row, row - 1] == 0
solve_one_row!(Σ, A, B, n, row, ws)
vB = view(B, 1:row - 1, 1:row - 1)
vtemp = view(ws.temp1, 1:row - 1, 1)
= view(Σ, 1:row -1 , row)
vA = view(A, 1:row - 1, 1:row - 1)
mul!(vtemp, vA, )
vA = view(A, 1:row - 1, row)
mul!(vB, vtemp, vA', 1.0, 1.0)
= view(Σ, 1:row, row)
vA = view(A, 1:row - 1, 1:row)
mul!(vtemp, vA, )
vA = view(A, 1:row - 1, row)
vB .= vB .+ vA.*vtemp'
row -= 1
else
solve_two_rows!(Σ, A, B, n, row, ws)
vB = view(B, 1:row - 2, 1:row - 2)
vtemp = view(ws.temp1, 1:row - 2, 1:2)
= view(Σ, 1:row - 2, row - 1:row)
vA = view(A, 1:row - 2, 1:row - 2)
mul!(vtemp, vA, )
vA = view(A, 1:row - 2, row - 1:row)
mul!(vB, vtemp, vA', 1.0, 1.0)
= view(Σ, 1:row, row - 1:row)
vA = view(A, 1:row - 2, 1:row)
mul!(vtemp, vA, )
vA = view(A, 1:row - 2, row - 1:row)
mul!(vB, vA, vtemp', 1.0, 1.0)
row -= 2
end
end
end
......@@ -36,9 +36,8 @@ end
for i = 1:100
s = Int64(floor(1000*rand()))
@show s
Random.seed!(s)
n = 5
local n = 5
X1 = randn(n, n)
X = X1'*X1
A = randn(n, n)
......@@ -56,7 +55,7 @@ for i = 1:100
s = Int64(floor(1000*rand()))
@show s
Random.seed!(s)
n = 2
local n = 2
X1 = randn(n, n)
X = X1'*X1
A = randn(n, n)
......@@ -68,6 +67,9 @@ for i = 1:100
AA = vcat(hcat([ 1.0 -0.5 0; 0 1.0 0], randn(m, n)),
[1 -1 0 0 0],
hcat(zeros(n, m+1), A))
if abs(det(AA)) < 1e-12
continue
end
B = X - A*X*A'
BB = vcat(hcat(I(m), zeros(m, n+1)),
[0 0 0.2 0 0],
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment