diff --git a/doc/manual/source/the-model-file.rst b/doc/manual/source/the-model-file.rst
index 06f8983352046134c4dcfe11108565a458269ef6..a54250b231cf83400da82cc2fbef5e2cc885fb8a 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 d37ccf524a7a8a32bbd1986018795f5d22520cfd..f0efa5fc5e89cc0c656e9e8d1f1c20e1aecf564e 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 166163b46b1a902507feb41650cf446314978650..94586ebcf3d45efea86834f4be5e2d93c2fb58b9 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 7b495d83efd6e534b3bbc41477cf286863083b63..0bf59162c58e0590a355b9f37075d0d1adafc7db 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 5b7e218dfdad5fa0661461230b0d28d1af0e09fc..9c1a9ac1b31a527455817e0ec87b0189021654e8 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 3da38a702a0a71d32a98b3ae71a430f1abd32960..46b793ca7408ba62eeb1dedb60ca3633f6bc28fd 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 14488dd8e03d8c56bc7e42738ef906989ce29892..8e0cb2b3d4424334bf1b94538e0d57d06216d137 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<SuiteSparse_long *>(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 8d31197e714a36948784d386e4274530450abedf..488870f16baf005447c9f2b2b4ed04217f122bb7 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 b0119c7d28ef224a6f33ca2f9718b0b1f298a565..c352784c12a7b0d7a1f0f42829adbc663f9007d5 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 f0bd25a90f6f2fd370778ed2a4d4d5bb48f27ff9..5299732d9b50f17b44c8b96977faf869c727821c 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