Commit 0ee486a5 authored by MichelJuillard's avatar MichelJuillard
Browse files

more on perfectforesight solvers

parent 02f3a3de
Pipeline #5330 failed with stage
in 1 minute and 46 seconds
......@@ -100,7 +100,7 @@ function jacobian_time_vec!(y::AbstractVector{Float64},
k += 1
end
end
mul!(y, offset_y, ws[1].jacobian, (npred+n_current)*nvar + 1, nvar,
@inbounds mul!(y, offset_y, ws[1].jacobian, (npred+n_current)*nvar + 1, nvar,
nfwrd, ws[1].temp_vec, 1, 1.0, 1.0)
return 0
end
......@@ -170,7 +170,7 @@ function hh!(hh::AbstractMatrix{Float64}, h::AbstractMatrix{Float64},
n::Int64, work1::AbstractMatrix{Float64},
work2::AbstractMatrix{Float64})
computes hh = [h, h(1), h(2), ..., h(n-1)] and hh(i) = -h*f*h(-1)
computes hh = [h, h(1), h(2), ..., h(n-1)] and h(i) = -h*f*h(-1)
"""
function hh!(hh::AbstractMatrix{Float64}, h::AbstractMatrix{Float64},
f::AbstractMatrix{Float64}, hf::AbstractMatrix{Float64},
......@@ -233,6 +233,7 @@ function preconditioner!(rout::AbstractVector{Float64},
mul!(rout, ir, g, 1, m, m, rout, ir_m, 1, 1)
ir += m
end
return rout
end
struct LREprecond
......@@ -331,11 +332,11 @@ struct GmresWs
work.jacobian,
options,
LREWs)
g1_1 = LREresults.g1_1
@inbounds for i = 1:LREWs.backward_nbr
k = LREWs.backward_indices[i]
for j = 1:n
g[j, k] = LREresults.g1_1[j, i]
end
vg = view(g, :, LREWs.backward_indices[i])
vg1 = view(g1_1, :, i)
copy!(vg, vg1)
end
if any(isnan.(g))
throw(ArgumentError("NaN in g"))
......
......@@ -14,7 +14,7 @@ function make_one_period!(JA::Jacobian,
get_jacobian!(ws, endogenous, exogenous, JA.steadystate, md, t)
i, j, v = findnz(sparse(ws.jacobian))
r1 = (t - 1)*nnz_period + 1
for el = 1:length(i)
@inbounds for el = 1:length(i)
if j[el] <= maxcol
kjel = jacobian_columns[j[el]]
JA.I[r1] = i[el] + oc
......@@ -43,7 +43,8 @@ function makeJacobian!(
r = 1
jacobian_columns = [i for (i, x) in enumerate(transpose(md.lead_lag_incidence)) if x > 0]
i, j, v = findnz(sparse(ws[1].jacobian))
for el = 1:length(i)
fill!(JA.V, 0.0)
@inbounds for el = 1:length(i)
if j[el] <= maxcol
kjel = jacobian_columns[j[el]]
if kjel > nvar
......@@ -64,7 +65,7 @@ function makeJacobian!(
i, j, v = findnz(sparse(ws[1].jacobian))
oc = (periods - 1)*nvar
r = (periods - 1)*nnz_period + 1
for el = 1:length(i)
@inbounds for el = 1:length(i)
if j[el] <= maxcol
kjel = jacobian_columns[j[el]]
if kjel <= 2*nvar
......@@ -78,7 +79,7 @@ function makeJacobian!(
# terminal condition
ic = md.n_bkwrd + md.n_both + md.n_current .+ (1:md.n_fwrd+md.n_both)
g = context.results.model_results[1].linearrationalexpectations.g1_1
CmultG!(
@inbounds CmultG!(
JA.tmp_nvar_npred,
JA.tmp_nvar_nfwrd,
JA.tmp_nfwrd_npred,
......@@ -88,7 +89,7 @@ function makeJacobian!(
md.i_fwrd_b,
)
(i, j, v) = findnz(sparse(JA.tmp_nvar_npred))
needed_space = periods*nnz_period + length(i) - 1
needed_space = r + length(i) - 1
extra_space = needed_space - length(JA.I)
if extra_space > 0
resize!(JA.I, needed_space)
......@@ -98,7 +99,7 @@ function makeJacobian!(
resize!(JA.colval, needed_space)
resize!(JA.nzval, needed_space)
end
for el = 1:length(i)
@inbounds for el = 1:length(i)
JA.I[r] = i[el] + oc
JA.J[r] = jacobian_columns[j[el]] + oc
JA.V[r] = v[el]
......@@ -108,6 +109,7 @@ function makeJacobian!(
resize!(JA.J, r - 1)
resize!(JA.V, r - 1)
n = periods * nvar
#length(colptr) == n + 1 && colptr[end] - 1 == length(rowval) == length(nzval)
return SparseArrays.sparse!(
JA.I,
JA.J,
......
......@@ -60,11 +60,11 @@ struct Jacobian
nz = periods * md.NNZDerivatives[1]
nrow = periods * md.endogenous_nbr
I = Vector{Int64}(undef, nz)
fill!(I, 1)
@inbounds fill!(I, 1)
J = Vector{Int64}(undef, nz)
fill!(J, 1)
@inbounds fill!(J, 1)
V = Vector{Float64}(undef, nz)
fill!(V, 0.0)
@inbounds fill!(V, 0.0)
klasstouch = Vector{Int64}(undef, nz)
colptr = Vector{Int64}(undef, nrow + 1)
rowptr = Vector{Int64}(undef, nrow + 1)
......@@ -97,7 +97,8 @@ end
function get_dynamic_endogenous_variables!(y::AbstractVector{Float64},
data::AbstractVector{Float64},
lli::Matrix{Int64},
m::Model, period::Int64)
m::Model,
period::Int64)
m, n = size(lli)
p = (period - 2)*n
@inbounds for i = 1:m
......@@ -133,7 +134,7 @@ function compute_jacobian(ws::JacTimesVec,
period::Int64)
dynamic! = m.dynamic!.dynamic!
fill!(ws.jacobian, 0.0)
Base.invokelatest(dynamic!,
@inbounds Base.invokelatest(dynamic!,
ws.temp_vec,
ws.residuals,
ws.jacobian,
......
......@@ -3,7 +3,7 @@ using Dynare
using LinearAlgebra.BLAS
#include("../src/perfectforesight/gmres_solver.jl")
#include("../src/perfectforesight/perfectforesight_solvers.jl")
context = @dynare "test/models/irbc/irbc1.mod" "savemacro" "-DN=300"
context = @dynare "test/models/irbc/irbc1.mod" "savemacro" "-DN=200"
md = context.models[1]
......@@ -53,13 +53,15 @@ include("models/irbc/irbc1Dynamic.jl")
endo_steadystate,
period)
@btime Dynare.compute_jacobian(work, dynamic_variables, exogenous,
endo_steadystate, md, period)
wsJ = Dynare.JacTimesVec(context)
@btime Dynare.compute_jacobian(wsJ, exogenous,
endo_steadystate, md, period)
periods = 400;
periods = 200;
preconditioner_window = 3
res = zeros(periods*md.endogenous_nbr)
res[[3, 4]] .= 0.1
NN = 200*md.endogenous_nbr
res[1:NN] .= 0.05*randn(NN)
rout = zeros(periods*md.endogenous_nbr)
ws = Dynare.GmresWs(periods, preconditioner_window, context, "GS")
ws.endogenous .= repeat(endo_steadystate, periods + 2)
......@@ -69,56 +71,63 @@ n = md.endogenous_nbr
work = context.work
function f1(periods)
JA = Dynare.Jacobian(context, periods)
A = Dynare.makeJacobian!(JA, ws.endogenous, ws.exogenous, context, periods)
A = Dynare.makeJacobian!(JA, ws.endogenous, ws.exogenous, context, periods, ws_threaded)
y = A\res
return 0
end
function f2(periods, ws; verbose = false)
function f2(periods, preconditioner_window; verbose = false)
ws = Dynare.GmresWs(periods, preconditioner_window, context, "GS")
JA = Dynare.Jacobian(context, periods)
A = Dynare.makeJacobian!(JA, ws.endogenous, ws.exogenous, context, periods)
A = Dynare.makeJacobian!(JA, ws.endogenous, ws.exogenous, context, periods, ws_threaded)
Dynare.gmres!(rout, A, res, log=false,
verbose=verbose, Pr=ws.P, initially_zero=true)
Dynare.ldiv!(ws.P, rout)
return 0
end
function f3(periods, ws; verbose = false)
function f3(periods, preconditioner_window; verbose = false)
ws = Dynare.GmresWs(periods, preconditioner_window, context, "GS")
ws.endogenous .= repeat(endo_steadystate, periods + 2)
ws.exogenous .= zeros(periods + 2, md.exogenous_nbr)
Dynare.gmres!(rout, ws.LREMap, res,Pr=ws.P,
log=false, verbose=verbose, initially_zero=true)
return 0
end
BLAS.set_num_threads(3)
z = zeros(n*periods)
y = zeros(n*periods)
Dynare.jacobian_time_vec!(z, ws.dynamic_variables, ws.residuals,
ws.endogenous, ws.exogenous,
ws.steadystate, ws.presiduals, ws.g,
ws.temp_vec, context.work, md, periods)
ws_threaded = [Dynare.JacTimesVec(context) for i=1:Threads.nthreads()]
Dynare.jacobian_time_vec!(z, ws.residuals, ws.endogenous,
ws.exogenous, ws.steadystate, ws.g, md,
periods, ws_threaded)
@show "get_jacobian"
@btime Dynare.get_jacobian!(work, ws.endogenous, ws.exogenous, endo_steadystate, md, 2)
wsJ = Dynare.JacTimesVec(context)
@btime Dynare.get_jacobian!(wsJ, ws.endogenous, ws.exogenous, endo_steadystate, md, 2)
@show "get_abc"
@btime Dynare.get_abc!(ws.a, ws.b, ws.c, work.jacobian, md)
@show "makeJacobian"
JA = Dynare.Jacobian(context, periods)
@btime Dynare.makeJacobian!(JA, ws.endogenous, ws.exogenous, context, periods)
A = Dynare.makeJacobian!(JA, ws.endogenous, ws.exogenous, context, periods)
@btime Dynare.makeJacobian!(JA, ws.endogenous, ws.exogenous, context, periods, ws_threaded)
A = Dynare.makeJacobian!(JA, ws.endogenous, ws.exogenous, context, periods, ws_threaded)
@show "A\res"
@btime y = A\res
@show "f1"
@show "f1 \\"
@btime f1(periods)
@show "f2"
preconditioner_window = 50
@show "f2 gmres A"
fill!(rout, 0.0)
@btime f2(periods, ws)
f2(periods, ws, verbose=true)
@btime f2(periods, preconditioner_window)
f2(periods, preconditioner_window, verbose=true)
fill!(rout, 0.0)
@show "f3"
@show "f3 gmres LREMap"
fill!(rout, 0.0)
@btime f3(periods, ws)
@btime f3(periods, preconditioner_window)
f3(periods, preconditioner_window, verbose=true)
#=
for i = 1:4
......
Markdown is supported
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