From 06f665e231c1a66863612eb75802ff4f0b698c57 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Wed, 15 Jun 2022 14:46:44 +0200 Subject: [PATCH] Perfect foresight: LBJ now available under stack_solve_algo=1 (with/without block/bytecode) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Previously, LBJ was available: – under stack_solve_algo=6 when neither block nor bytecode were present – under stack_solve_algo=1 with either block or bytecode (but the documentation was not making it clear that it was LBJ) This commit merges the two values for the option, and makes them interchangeable. LBJ should now be invoked with stack_solve_algo=1 (but stack_solve_algo=6 is kept for compatibility, and is a synonymous). --- doc/manual/source/the-model-file.rst | 23 +++++++++++-------- .../perfect_foresight_solver_core.m | 2 +- .../private/check_input_arguments.m | 12 ++++------ .../solve_block_decomposed_problem.m | 4 ++-- matlab/solve_one_boundary.m | 9 ++------ matlab/solve_two_boundaries.m | 2 +- mex/sources/bytecode/SparseMatrix.cc | 20 +++++++++------- tests/deterministic_simulations/lbj/rbc.mod | 21 +++++++++++++++++ tests/run_block_byte_tests_matlab.m | 6 ++--- tests/run_block_byte_tests_octave.m | 6 ++--- 10 files changed, 62 insertions(+), 43 deletions(-) diff --git a/doc/manual/source/the-model-file.rst b/doc/manual/source/the-model-file.rst index 06f898335..a54250b23 100644 --- a/doc/manual/source/the-model-file.rst +++ b/doc/manual/source/the-model-file.rst @@ -3533,10 +3533,9 @@ method to solve the simultaneous equation system. Because the resulting Jacobian is in the order of ``n`` by ``T`` and hence will be very large for long simulations with many variables, Dynare makes use of the sparse matrix capacities of MATLAB/Octave. A slower but -potentially less memory consuming alternative (``stack_solve_algo=6``) +potentially less memory consuming alternative (``stack_solve_algo=1``) is based on a Newton-type algorithm first proposed by *Laffargue -(1990)* and *Boucekkine (1995)*, which uses relaxation -techniques. Thereby, the algorithm avoids ever storing the full +(1990)* and *Boucekkine (1995)*, which avoids ever storing the full Jacobian. The details of the algorithm can be found in *Juillard (1996)*. The third type of algorithms makes use of block decomposition techniques (divide-and-conquer methods) that exploit the structure of @@ -3654,9 +3653,15 @@ speed-up on large models. ``1`` - Use a Newton algorithm with a sparse LU solver at each - iteration (requires ``bytecode`` and/or ``block`` - option, see :ref:`model-decl`). + Use the Laffargue-Boucekkine-Juillard (LBJ) algorithm proposed + in *Juillard (1996)*. It is slower than ``stack_solve_algo=0``, + but may be less memory consuming on big models. Note that if the + ``block`` option is used (see :ref:`model-decl`), a simple + Newton algorithm with sparse matrices is used for blocks which + are purely backward or forward (of type ``SOLVE BACKWARD`` or + ``SOLVE FORWARD``, see :comm:`model_info`), since LBJ only makes + sense on blocks with both leads and lags (of type ``SOLVE TWO + BOUNDARIES``). ``2`` @@ -3686,10 +3691,8 @@ speed-up on large models. ``6`` - Use the historical algorithm proposed in *Juillard - (1996)*: it is slower than ``stack_solve_algo=0``, but - may be less memory consuming on big models (not - available with ``bytecode`` and/or ``block`` options). + Synonymous for ``stack_solve_algo=1``. Kept for historical + reasons. ``7`` diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m index d37ccf524..f0efa5fc5 100644 --- a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m +++ b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m @@ -94,7 +94,7 @@ else [oo_.endo_simul, oo_.deterministic_simulation] = ... sim1(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_); end - case 6 + case {1 6} if options_.linear_approximation error('Invalid value of stack_solve_algo option!') end diff --git a/matlab/perfect-foresight-models/private/check_input_arguments.m b/matlab/perfect-foresight-models/private/check_input_arguments.m index 166163b46..94586ebcf 100644 --- a/matlab/perfect-foresight-models/private/check_input_arguments.m +++ b/matlab/perfect-foresight-models/private/check_input_arguments.m @@ -3,7 +3,7 @@ function check_input_arguments(DynareOptions, DynareModel, DynareResults) %Conducts checks for inconsistent/missing inputs to deterministic %simulations -% Copyright © 2015-2017 Dynare Team +% Copyright © 2015-2022 Dynare Team % % This file is part of Dynare. % @@ -25,19 +25,15 @@ if DynareOptions.stack_solve_algo < 0 || DynareOptions.stack_solve_algo > 7 end if ~DynareOptions.block && ~DynareOptions.bytecode && DynareOptions.stack_solve_algo ~= 0 ... - && DynareOptions.stack_solve_algo ~= 6 && DynareOptions.stack_solve_algo ~= 7 - error('perfect_foresight_solver:ArgCheck','PERFECT_FORESIGHT_SOLVER: you must use stack_solve_algo=0 or stack_solve_algo=6 when not using block nor bytecode option') + && DynareOptions.stack_solve_algo ~= 1 && DynareOptions.stack_solve_algo ~= 6 ... + && DynareOptions.stack_solve_algo ~= 7 + error('perfect_foresight_solver:ArgCheck','PERFECT_FORESIGHT_SOLVER: you must use stack_solve_algo={0,1,6,7} when not using block nor bytecode option') end if DynareOptions.block && ~DynareOptions.bytecode && DynareOptions.stack_solve_algo == 5 error('perfect_foresight_solver:ArgCheck','PERFECT_FORESIGHT_SOLVER: you can''t use stack_solve_algo = 5 without bytecode option') end -if (DynareOptions.block || DynareOptions.bytecode) && DynareOptions.stack_solve_algo == 6 - error('perfect_foresight_solver:ArgCheck','PERFECT_FORESIGHT_SOLVER: you can''t use stack_solve_algo = 6 with block or bytecode option') -end - - if isempty(DynareResults.endo_simul) || any(size(DynareResults.endo_simul) ~= [ DynareModel.endo_nbr, DynareModel.maximum_lag+DynareOptions.periods+DynareModel.maximum_lead ]) if DynareOptions.initval_file diff --git a/matlab/perfect-foresight-models/solve_block_decomposed_problem.m b/matlab/perfect-foresight-models/solve_block_decomposed_problem.m index 7b495d83e..0bf59162c 100644 --- a/matlab/perfect-foresight-models/solve_block_decomposed_problem.m +++ b/matlab/perfect-foresight-models/solve_block_decomposed_problem.m @@ -22,8 +22,8 @@ cutoff = 1e-15; if options_.stack_solve_algo==0 mthd='Sparse LU'; -elseif options_.stack_solve_algo==1 - mthd='Relaxation'; +elseif options_.stack_solve_algo==1 || options_.stack_solve_algo==6 + mthd='LBJ'; elseif options_.stack_solve_algo==2 mthd='GMRES'; elseif options_.stack_solve_algo==3 diff --git a/matlab/solve_one_boundary.m b/matlab/solve_one_boundary.m index 5b7e218df..9c1a9ac1b 100644 --- a/matlab/solve_one_boundary.m +++ b/matlab/solve_one_boundary.m @@ -23,12 +23,7 @@ function [y, T, oo_, info] = solve_one_boundary(fname, y, x, params, steady_stat % solve_tolf [double] convergence criteria % cutoff [double] cutoff to correct the direction in Newton in case % of singular jacobian matrix -% stack_solve_algo [integer] linear solver method used in the -% Newton algorithm : -% - 1 sparse LU -% - 2 GMRES -% - 3 BicGStab -% - 4 Optimal path length +% stack_solve_algo [integer] linear solver method used in the Newton algorithm % is_forward [logical] Whether the block has to be solved forward % If false, the block is solved backward % is_dynamic [logical] If true, the block belongs to the dynamic file @@ -212,7 +207,7 @@ for it_=start:incr:finish y(y_index_eq, it_) = ya; %% Recompute temporary terms, since they are not given as output of lnsrch1 [~, ~, T(:, it_)] = feval(fname, Block_Num, dynvars_from_endo_simul(y, it_, M), x, params, steady_state, T(:, it_), it_, false); - elseif (is_dynamic && (stack_solve_algo==1 || stack_solve_algo==0)) || (~is_dynamic && options.solve_algo==6) + elseif (is_dynamic && (stack_solve_algo==1 || stack_solve_algo==0 || stack_solve_algo==6)) || (~is_dynamic && options.solve_algo==6) if verbose && ~is_dynamic disp('steady: Sparse LU ') end diff --git a/matlab/solve_two_boundaries.m b/matlab/solve_two_boundaries.m index 3da38a702..46b793ca7 100644 --- a/matlab/solve_two_boundaries.m +++ b/matlab/solve_two_boundaries.m @@ -182,7 +182,7 @@ while ~(cvg || iter>maxit_) dx = g1a\b- ya; ya = ya + lambda*dx; y(y_index, y_kmin+(1:periods))=reshape(ya',length(y_index),periods); - elseif stack_solve_algo==1 + elseif stack_solve_algo==1 || stack_solve_algo==6 for t=1:periods first_elem = (t-1)*Blck_size+1; last_elem = t*Blck_size; diff --git a/mex/sources/bytecode/SparseMatrix.cc b/mex/sources/bytecode/SparseMatrix.cc index 14488dd8e..8e0cb2b3d 100644 --- a/mex/sources/bytecode/SparseMatrix.cc +++ b/mex/sources/bytecode/SparseMatrix.cc @@ -319,7 +319,8 @@ dynSparseMatrix::Read_SparseMatrix(const string &file_name, int Size, int period for (int j = 0; j < Size; j++) IM_i[{ j, Size*(periods+y_kmax), 0 }] = j; } - else if (stack_solve_algo >= 0 && stack_solve_algo <= 4) + else if ((stack_solve_algo >= 0 && stack_solve_algo <= 4) + || stack_solve_algo == 6) { for (int i = 0; i < u_count_init-Size; i++) { @@ -350,7 +351,8 @@ dynSparseMatrix::Read_SparseMatrix(const string &file_name, int Size, int period IM_i[{ eq, var, lag }] = val; } } - else if (((stack_solve_algo >= 0 && stack_solve_algo <= 4) && !steady_state) + else if ((((stack_solve_algo >= 0 && stack_solve_algo <= 4) + || stack_solve_algo == 6) && !steady_state) || ((solve_algo >= 6 || solve_algo <= 8) && steady_state)) { for (int i = 0; i < u_count_init; i++) @@ -3799,7 +3801,8 @@ dynSparseMatrix::Simulate_One_Boundary(int block_num, int y_size, int y_kmin, in if (!x0_m) throw FatalExceptionHandling(" in Simulate_One_Boundary, can't allocate x0_m vector\n"); if (!((solve_algo == 6 && steady_state) - || ((stack_solve_algo == 0 || stack_solve_algo == 1 || stack_solve_algo == 4) && !steady_state))) + || ((stack_solve_algo == 0 || stack_solve_algo == 1 || stack_solve_algo == 4 + || stack_solve_algo == 6) && !steady_state))) { Init_Matlab_Sparse_Simple(size, IM_i, A_m, b_m, zero_solution, x0_m); A_m_save = mxDuplicateArray(A_m); @@ -3839,7 +3842,7 @@ dynSparseMatrix::Simulate_One_Boundary(int block_num, int y_size, int y_kmin, in Solve_Matlab_GMRES(A_m, b_m, size, slowc, block_num, false, it_, x0_m); else if ((solve_algo == 8 && steady_state) || (stack_solve_algo == 3 && !steady_state)) Solve_Matlab_BiCGStab(A_m, b_m, size, slowc, block_num, false, it_, x0_m, preconditioner); - else if ((solve_algo == 6 && steady_state) || ((stack_solve_algo == 0 || stack_solve_algo == 1 || stack_solve_algo == 4) && !steady_state)) + else if ((solve_algo == 6 && steady_state) || ((stack_solve_algo == 0 || stack_solve_algo == 1 || stack_solve_algo == 4 || stack_solve_algo == 6) && !steady_state)) Solve_LU_UMFPack(Ap, Ai, Ax, b, size, size, slowc, false, it_); } return singular_system; @@ -3898,7 +3901,7 @@ dynSparseMatrix::Simulate_Newton_One_Boundary(bool forward) test_mxMalloc(r, __LINE__, __FILE__, __func__, size*sizeof(double)); iter = 0; if ((solve_algo == 6 && steady_state) - || ((stack_solve_algo == 0 || stack_solve_algo == 1 || stack_solve_algo == 4) && !steady_state)) + || ((stack_solve_algo == 0 || stack_solve_algo == 1 || stack_solve_algo == 4 || stack_solve_algo == 6) && !steady_state)) { Ap_save = static_cast(mxMalloc((size + 1) * sizeof(SuiteSparse_long))); test_mxMalloc(Ap_save, __LINE__, __FILE__, __func__, (size + 1) * sizeof(SuiteSparse_long)); @@ -3937,7 +3940,7 @@ dynSparseMatrix::Simulate_Newton_One_Boundary(bool forward) solve_linear(block_num, y_size, y_kmin, y_kmax, size, 0); } if ((solve_algo == 6 && steady_state) - || ((stack_solve_algo == 0 || stack_solve_algo == 1 || stack_solve_algo == 4) && !steady_state)) + || ((stack_solve_algo == 0 || stack_solve_algo == 1 || stack_solve_algo == 4 || stack_solve_algo == 6) && !steady_state)) { mxFree(Ap_save); mxFree(Ai_save); @@ -4126,7 +4129,8 @@ dynSparseMatrix::Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin mexPrintf("MODEL SIMULATION: (method=Sparse LU)\n"); break; case 1: - mexPrintf("MODEL SIMULATION: (method=Relaxation)\n"); + case 6: + mexPrintf("MODEL SIMULATION: (method=LBJ)\n"); break; case 2: mexPrintf(preconditioner_print_out("MODEL SIMULATION: (method=GMRES)\n", preconditioner, false).c_str()); @@ -4178,7 +4182,7 @@ dynSparseMatrix::Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin } if (stack_solve_algo == 0 || stack_solve_algo == 4) Solve_LU_UMFPack(Ap, Ai, Ax, b, Size * periods, Size, slowc, true, 0, vector_table_conditional_local); - else if (stack_solve_algo == 1) + else if (stack_solve_algo == 1 || stack_solve_algo == 6) Solve_Matlab_Relaxation(A_m, b_m, Size, slowc, 0); else if (stack_solve_algo == 2) Solve_Matlab_GMRES(A_m, b_m, Size, slowc, blck, true, 0, x0_m); diff --git a/tests/deterministic_simulations/lbj/rbc.mod b/tests/deterministic_simulations/lbj/rbc.mod index 8d31197e7..488870f16 100644 --- a/tests/deterministic_simulations/lbj/rbc.mod +++ b/tests/deterministic_simulations/lbj/rbc.mod @@ -79,6 +79,27 @@ end oo0 = oo_; +perfect_foresight_setup(periods=400); +perfect_foresight_solver(stack_solve_algo=1); + +if ~oo_.deterministic_simulation.status + error('Perfect foresight simulation failed') +end + +oo1 = oo_; + +maxabsdiff = max(max(abs(oo0.endo_simul-oo1.endo_simul))); + +if max(max(abs(oo0.endo_simul-oo1.endo_simul)))>options_.dynatol.x + error('stack_solve_algo={0,1} return different paths for the endogenous variables!') +else + skipline() + fprintf('Maximum (absolute) differrence between paths is %s', num2str(maxabsdiff)) + skipline() +end + +% Also test stack_solve_algo=6, which is a synonymous for stack_solve_algo=1 + perfect_foresight_setup(periods=400); perfect_foresight_solver(stack_solve_algo=6); diff --git a/tests/run_block_byte_tests_matlab.m b/tests/run_block_byte_tests_matlab.m index b0119c7d2..c352784c1 100644 --- a/tests/run_block_byte_tests_matlab.m +++ b/tests/run_block_byte_tests_matlab.m @@ -49,13 +49,13 @@ for blockFlag = 0:1 default_stack_solve_algo = 0; if ~blockFlag && storageFlag ~= 2 solve_algos = [1:4 9]; - stack_solve_algos = [0 6]; + stack_solve_algos = [0 1 6]; elseif blockFlag && storageFlag ~= 2 solve_algos = [1:4 6:9]; - stack_solve_algos = 0:4; + stack_solve_algos = [0:4 6]; else solve_algos = 1:9; - stack_solve_algos = 0:5; + stack_solve_algos = 0:6; end if has_optimization_toolbox solve_algos = [ solve_algos 0 ]; diff --git a/tests/run_block_byte_tests_octave.m b/tests/run_block_byte_tests_octave.m index f0bd25a90..5299732d9 100644 --- a/tests/run_block_byte_tests_octave.m +++ b/tests/run_block_byte_tests_octave.m @@ -45,13 +45,13 @@ for blockFlag = 0:1 default_stack_solve_algo = 0; if !blockFlag && storageFlag != 2 solve_algos = [0:4 9]; - stack_solve_algos = [0 6]; + stack_solve_algos = [0 1 6]; elseif blockFlag && storageFlag != 2 solve_algos = [0:4 6:9]; - stack_solve_algos = 0:4; + stack_solve_algos = [0:4 6]; else solve_algos = 0:9; - stack_solve_algos = 0:5; + stack_solve_algos = 0:6; endif # Workaround for strange race condition related to the static/dynamic