diff --git a/matlab/dr_block.m b/matlab/dr_block.m
index 7f7df97cc83b1bd3bf2f92050d6cf3fab2c7c920..e7d7d0398f8817467b7c0818e5dbc634a3dfa57d 100644
--- a/matlab/dr_block.m
+++ b/matlab/dr_block.m
@@ -567,7 +567,9 @@ for i = 1:Size;
if block_type == 5
vghx_other = - inv(kron(eye(size(D_,2)), A_) + kron(C_', B_)) * vec(D_);
ghx_other = reshape(vghx_other, size(D_,1), size(D_,2));
- else
+ elseif options_.sylvester_fp == 1
+ ghx_other = gensylv_fp(A_, B_, C_, D_, i);
+ else
[err, ghx_other] = gensylv(1, A_, B_, C_, -D_);
end;
if options_.aim_solver ~= 1 && options_.use_qzdiv
diff --git a/matlab/dsge_likelihood.m b/matlab/dsge_likelihood.m
index fbb2d1eab32148bd45b1813bcbe388a62591ca74..282a8abddf715eedda08f6f28e378379efd08c2e 100644
--- a/matlab/dsge_likelihood.m
+++ b/matlab/dsge_likelihood.m
@@ -366,7 +366,11 @@ switch DynareOptions.lik_init
% Use standard kalman filter except if the univariate filter is explicitely choosen.
kalman_algo = 1;
end
- Pstar = lyapunov_symm(T,R*Q*R',DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold);
+ if DynareOptions.lyapunov_fp == 1
+ Pstar = lyapunov_symm(T,Q,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, 4, R);
+ else
+ Pstar = lyapunov_symm(T,R*Q*R',DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold);
+ end;
Pinf = [];
a = zeros(mm,1);
Zflag = 0;
diff --git a/matlab/gensylv_fp.m b/matlab/gensylv_fp.m
new file mode 100644
index 0000000000000000000000000000000000000000..f56a65a9183b7939b5ef45d0ef67a529e3dc71bb
--- /dev/null
+++ b/matlab/gensylv_fp.m
@@ -0,0 +1,73 @@
+function X = gensylv_fp(A, B, C, D, block)
+% function X = gensylv_fp(A, B, C, D)
+% Solve the Sylvester equation:
+% A * X + B * X * C + D = 0
+% INPUTS
+% A
+% B
+% C
+% D
+% block : block number (for storage purpose)
+% OUTPUTS
+% X solution
+%
+% ALGORITHM
+% fixed point method
+% MARLLINY MONSALVE (2008): "Block linear method for large scale
+% Sylvester equations", Computational & Applied Mathematics, Vol 27, n�1,
+% p47-59
+% ||A^-1||.||B||.||C|| < 1 is a suffisant condition:
+% - to get a unique solution for the Sylvester equation
+% - to get a convergent fixed-point algorithm
+%
+% SPECIAL REQUIREMENTS
+% none.
+% Copyright (C) 1996-2010 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
+
+%tol = 1e-07;
+%tol = 1e-13;
+tol = 1e-12;
+evol = 100;
+A1 = inv(A);
+eval(['persistent hxo_' int2str(block) ';']);
+hxo = eval(['hxo_' int2str(block) ';']);
+if isempty(hxo)
+ X = zeros(size(B, 2), size(C, 1));
+else
+ X = hxo;
+end;
+it_fp = 0;
+maxit_fp = 1000;
+Z = - (B * X * C + D);
+while it_fp < maxit_fp && evol > tol;
+ %X_old = X;
+ %X = - A1 * ( B * X * C + D);
+ %evol = max(max(abs(X - X_old)));
+ X = A1 * Z;
+ Z_old = Z;
+ Z = - (B * X * C + D);
+ evol = max(sum(abs(Z - Z_old))); %norm_1
+ %evol = max(sum(abs(Z - Z_old)')); %norm_inf
+ it_fp = it_fp + 1;
+end;
+%fprintf('sylvester it_fp=%d evol=%g | ',it_fp,evol);
+if evol < tol
+ eval(['hxo_' int2str(block) ' = X;']);
+else
+ error(['convergence not achieved in fixed point solution of Sylvester equation after ' int2str(it_fp) ' iterations']);
+end;
\ No newline at end of file
diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index 2ee53d90f3d7a75d351085da5a8509ce5c1d7be0..0e5dd3eefc9c00d3585ca5ca60331254b2208de6 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -382,6 +382,14 @@ options_.use_dll = 0;
% model evaluated using bytecode.dll
options_.bytecode = 0;
+% use a fixed point method to solve Sylvester equation (for large scale
+% models)
+options_.sylvester_fp = 0;
+
+% use a fixed point method to solve Lyapunov equation (for large scale
+% models)
+options_.lyapunov_fp = 0;
+
% dates for historical time series
options_.initial_date.freq = 1;
options_.initial_date.period = 1;
diff --git a/matlab/lyapunov_symm.m b/matlab/lyapunov_symm.m
index d32f644394bb8f537abe97ad2b64768e4669c169..afa763c82f43036ff0332c643a79609bf1d7438d 100644
--- a/matlab/lyapunov_symm.m
+++ b/matlab/lyapunov_symm.m
@@ -1,4 +1,4 @@
-function [x,u] = lyapunov_symm(a,b,qz_criterium,lyapunov_complex_threshold,method)
+function [x,u] = lyapunov_symm(a,b,qz_criterium,lyapunov_complex_threshold,method, R)
% Solves the Lyapunov equation x-a*x*a' = b, for b and x symmetric matrices.
% If a has some unit roots, the function computes only the solution of the stable subsystem.
%
@@ -12,6 +12,7 @@ function [x,u] = lyapunov_symm(a,b,qz_criterium,lyapunov_complex_threshold,metho
% variables and the schur decomposition is triggered.
% method=2 then U, T, n and k are declared as persistent
% variables and the schur decomposition is not performed.
+% method=3 fixed point method
% OUTPUTS
% x: [double] m*m solution matrix of the lyapunov equation, where m is the dimension of the stable subsystem.
% u: [double] Schur vectors associated with unit roots
@@ -38,10 +39,55 @@ function [x,u] = lyapunov_symm(a,b,qz_criterium,lyapunov_complex_threshold,metho
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
-
if nargin<5
method = 0;
end
+
+if method == 3
+ persistent X method1;
+ if ~isempty(method1)
+ method = method1;
+ end;
+ fprintf(' [methode=%d] ',method);
+ if method == 3
+ %tol = 1e-8;
+ tol = 1e-10;
+ it_fp = 0;
+ evol = 100;
+ if isempty(X)
+ X = b;
+ max_it_fp = 2000;
+ else
+ max_it_fp = 300;
+ end;
+ at = a';
+ %fixed point iterations
+ while evol > tol && it_fp < max_it_fp;
+ X_old = X;
+ X = a * X * at + b;
+ evol = max(sum(abs(X - X_old))); %norm_1
+ %evol = max(sum(abs(X - X_old)')); %norm_inf
+ it_fp = it_fp + 1;
+ end;
+ fprintf('lyapunov it_fp=%d evol=%g\n',it_fp,evol);
+ if it_fp >= max_it_fp
+ disp(['convergence not achieved in solution of Lyapunov equation after ' int2str(it_fp) ' iterations, switching method from 3 to 0']);
+ method1 = 0;
+ method = 0;
+ else
+ method1 = 3;
+ x = X;
+ return;
+ end;
+ end;
+elseif method == 4
+ % works only with Matlab System Control toolbox
+ chol_b = R*chol(b,'lower');
+ Rx = dlyapchol(a,chol_b);
+ x = Rx' * Rx;
+ return;
+end;
+
if method
persistent U T k n
else
@@ -113,4 +159,4 @@ if i == 1
x(1,1) = (B(1,1)+c)/(1-T(1,1)*T(1,1));
end
x = U(:,k+1:end)*x*U(:,k+1:end)';
-u = U(:,1:k);
\ No newline at end of file
+u = U(:,1:k);
diff --git a/mex/sources/bytecode/ErrorHandling.hh b/mex/sources/bytecode/ErrorHandling.hh
index 9856b4338c719f3039411c8462397873d4ed4202..a27e5fe44b6af308fb2496d5bb829788a3a57dd0 100644
--- a/mex/sources/bytecode/ErrorHandling.hh
+++ b/mex/sources/bytecode/ErrorHandling.hh
@@ -113,7 +113,10 @@ public:
value2(value2_arg)
{
ostringstream tmp;
- tmp << " with X=" << value1 << "\n";
+ if (abs(value1) > 1e-10 )
+ tmp << " with X=" << value1 << "\n";
+ else
+ tmp << " with X=" << value1 << " and a=" << value2 << "\n";
completeErrorMsg(tmp.str());
};
};
diff --git a/mex/sources/bytecode/Interpreter.cc b/mex/sources/bytecode/Interpreter.cc
index d8dc69987276d48f8be780e71d25a4a75d02dc03..a75943b5acd05eb3c5ee6d036d3d8ca9f6c78275 100644
--- a/mex/sources/bytecode/Interpreter.cc
+++ b/mex/sources/bytecode/Interpreter.cc
@@ -71,7 +71,7 @@ double
Interpreter::pow1(double a, double b)
{
double r = pow_(a, b);
- if (isnan(r))
+ if (isnan(r) || isinf(r))
{
res1 = NAN;
r = 0.0000000000000000000000001;
@@ -85,7 +85,7 @@ double
Interpreter::divide(double a, double b)
{
double r = a / b;
- if (isinf(r))
+ if (isnan(r) || isinf(r))
{
res1 = NAN;
r = 1e70;
@@ -99,7 +99,7 @@ double
Interpreter::log1(double a)
{
double r = log(a);
- if (isnan(r))
+ if (isnan(r) || isinf(r))
{
res1 = NAN;
r = -1e70;
@@ -113,7 +113,7 @@ double
Interpreter::log10_1(double a)
{
double r = log(a);
- if (isnan(r))
+ if (isnan(r) || isinf(r))
{
res1 = NAN;
r = -1e70;
@@ -1718,6 +1718,16 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
}
}
+
+void
+Interpreter::SingularDisplay(int Per_u_, bool evaluate, int Block_Count, int size, bool steady_state, it_code_type begining)
+{
+ it_code = begining;
+ compute_block_time(Per_u_, evaluate, Block_Count, size, steady_state);
+ Singular_display(Block_Count, size, steady_state, begining);
+}
+
+
int
Interpreter::simulate_a_block(const int size, const int type, string file_name, string bin_basename, bool Gaussian_Elimination, bool steady_state, bool print_it, int block_num,
const bool is_linear, const int symbol_table_endo_nbr, const int Block_List_Max_Lag, const int Block_List_Max_Lead, const int u_count_int)
@@ -1726,6 +1736,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
int i;
bool cvg;
bool result = true;
+ bool singular_system;
double *y_save;
res1 = 0;
#ifdef DEBUG
@@ -1979,7 +1990,10 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
if (cvg)
continue;
int prev_iter = iter;
- Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, print_it, cvg, iter, true, stack_solve_algo, solve_algo);
+ singular_system = Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, print_it, cvg, iter, true, stack_solve_algo, solve_algo);
+ if (singular_system)
+ SingularDisplay(0, false, block_num, size, steady_state, begining);
+
iter++;
if (iter > prev_iter)
{
@@ -2024,7 +2038,9 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
}
else
cvg = false;
- Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, print_it, cvg, iter, true, stack_solve_algo, solve_algo);
+ singular_system = Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, print_it, cvg, iter, true, stack_solve_algo, solve_algo);
+ if (singular_system)
+ SingularDisplay(0, false, block_num, size, steady_state, begining);
if (!result)
{
mexPrintf(" in Solve Forward complete, convergence not achieved in block %d\n", Block_Count+1);
@@ -2076,7 +2092,9 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
if (cvg)
continue;
int prev_iter = iter;
- Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, print_it, cvg, iter, false, stack_solve_algo, solve_algo);
+ singular_system = Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, print_it, cvg, iter, false, stack_solve_algo, solve_algo);
+ if (singular_system)
+ SingularDisplay(0, false, block_num, size, steady_state, begining);
iter++;
if (iter > prev_iter)
{
@@ -2123,7 +2141,9 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
}
else
cvg = false;
- Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, print_it, cvg, iter, false, stack_solve_algo, solve_algo);
+ singular_system = Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, print_it, cvg, iter, false, stack_solve_algo, solve_algo);
+ if (singular_system)
+ SingularDisplay(0, false, block_num, size, steady_state, begining);
}
}
}
@@ -2181,7 +2201,9 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
if (cvg)
continue;
int prev_iter = iter;
- Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, print_it, cvg, iter, true, stack_solve_algo, solve_algo);
+ singular_system = Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, print_it, cvg, iter, true, stack_solve_algo, solve_algo);
+ if (singular_system)
+ SingularDisplay(0, false, block_num, size, steady_state, begining);
iter++;
if (iter > prev_iter)
{
@@ -2225,7 +2247,9 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
}
else
cvg = false;
- Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, print_it, cvg, iter, true, stack_solve_algo, solve_algo);
+ singular_system = Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, print_it, cvg, iter, true, stack_solve_algo, solve_algo);
+ if (singular_system)
+ SingularDisplay(0, false, block_num, size, steady_state, begining);
if (!result)
{
mexPrintf(" in Solve Backward complete, convergence not achieved in block %d\n", Block_Count+1);
@@ -2277,7 +2301,9 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
if (cvg)
continue;
int prev_iter = iter;
- Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, print_it, cvg, iter, false, stack_solve_algo, solve_algo);
+ singular_system = Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, print_it, cvg, iter, false, stack_solve_algo, solve_algo);
+ if (singular_system)
+ SingularDisplay(0, false, block_num, size, steady_state, begining);
iter++;
if (iter > prev_iter)
{
@@ -2320,7 +2346,9 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
}
else
cvg = false;
- Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, print_it, cvg, iter, false, stack_solve_algo, solve_algo);
+ singular_system = Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, print_it, cvg, iter, false, stack_solve_algo, solve_algo);
+ if (singular_system)
+ SingularDisplay(0, false, block_num, size, steady_state, begining);
}
}
}
diff --git a/mex/sources/bytecode/Interpreter.hh b/mex/sources/bytecode/Interpreter.hh
index 8b971e22710697e16e2b9e5e76273bfee7c1538e..22778b14dbaf97193635c1db1e513e20e1111353 100644
--- a/mex/sources/bytecode/Interpreter.hh
+++ b/mex/sources/bytecode/Interpreter.hh
@@ -61,6 +61,7 @@ protected:
void print_a_block(const int size, const int type, string bin_basename, bool steady_state, int block_num,
const bool is_linear, const int symbol_table_endo_nbr, const int Block_List_Max_Lag,
const int Block_List_Max_Lead, const int u_count_int, int block);
+ void SingularDisplay(int Per_u_, bool evaluate, int Block_Count, int size, bool steady_state, it_code_type begining);
vector<Block_contain_type> Block_Contain;
code_liste_type code_liste;
it_code_type it_code;
diff --git a/mex/sources/bytecode/SparseMatrix.cc b/mex/sources/bytecode/SparseMatrix.cc
index a7cf3769882b27ab9ee1014ac0ef4dc6208a9b39..db577ac11b97f2c87d7e885f022bdc805d6b0c67 100644
--- a/mex/sources/bytecode/SparseMatrix.cc
+++ b/mex/sources/bytecode/SparseMatrix.cc
@@ -366,7 +366,7 @@ SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, i
}
void
-SparseMatrix::Simple_Init(int it_, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM, bool &zero_solution)
+SparseMatrix::Simple_Init(int Size, map<pair<pair<int, int>, int>, int> &IM, bool &zero_solution)
{
int i, eq, var, lag;
map<pair<pair<int, int>, int>, int>::iterator it4;
@@ -2153,6 +2153,94 @@ SparseMatrix::Solve_Matlab_BiCGStab(mxArray *A_m, mxArray *b_m, int Size, double
}
void
+SparseMatrix::Singular_display(int block, int Size, bool steady_state, it_code_type it_code)
+{
+ bool zero_solution;
+ Simple_Init(Size, IM_i, zero_solution);
+ NonZeroElem *first;
+ mxArray *rhs[1];
+ rhs[0] = mxCreateDoubleMatrix(Size, Size, mxREAL);
+ double *pind;
+ pind = mxGetPr(rhs[0]);
+ for (int j = 0; j < Size * Size; j++)
+ pind[j] = 0.0;
+ for (int ii = 0; ii < Size; ii++)
+ {
+ int nb_eq = At_Col(ii, &first);
+ for (int j = 0; j < nb_eq; j++)
+ {
+ int k = first->u_index;
+ int jj = first->r_index;
+ pind[ii * Size + jj ] = u[k];
+ first = first->NZE_C_N;
+ }
+ }
+ mxArray *lhs[3];
+ mexCallMATLAB(3, lhs, 1, rhs, "svd");
+ mxArray* SVD_u = lhs[0];
+ mxArray* SVD_s = lhs[1];
+ mxArray* SVD_v = lhs[2];
+ double *SVD_ps = mxGetPr(SVD_s);
+ double *SVD_pu = mxGetPr(SVD_u);
+ for (int i = 0; i < Size; i++)
+ {
+ if (abs(SVD_ps[i * (1 + Size)]) < 1e-12)
+ {
+ mexPrintf(" The following equations form a linear combination:\n ");
+ double max_u = 0;
+ for (int j = 0; j < Size; j++)
+ if (abs(SVD_pu[j + i * Size]) > abs(max_u))
+ max_u = SVD_pu[j + i * Size];
+ vector<int> equ_list;
+ for (int j = 0; j < Size; j++)
+ {
+ double rr = SVD_pu[j + i * Size] / max_u;
+ if ( rr < -1e-10)
+ {
+ equ_list.push_back(j);
+ if (rr != -1)
+ mexPrintf(" - %3.2f*Dequ_%d_dy",abs(rr),j+1);
+ else
+ mexPrintf(" - Dequ_%d_dy",j+1);
+ }
+ else if (rr > 1e-10)
+ {
+ equ_list.push_back(j);
+ if (j > 0)
+ if (rr != 1)
+ mexPrintf(" + %3.2f*Dequ_%d_dy",rr,j+1);
+ else
+ mexPrintf(" + Dequ_%d_dy",j+1);
+ else
+ if (rr != 1)
+ mexPrintf(" %3.2f*Dequ_%d_dy",rr,j+1);
+ else
+ mexPrintf(" Dequ_%d_dy",j+1);
+ }
+ }
+ mexPrintf(" = 0\n");
+ /*mexPrintf(" with:\n");
+ it_code = get_begin_block(block);
+ for (int j=0; j < Size; j++)
+ {
+ if (find(equ_list.begin(), equ_list.end(), j) != equ_list.end())
+ mexPrintf(" equ_%d: %s\n",j, print_expression(it_code_expr, false, Size, block, steady_state, 0, 0, it_code, true).c_str());
+ }*/
+ }
+ }
+ mxDestroyArray(lhs[0]);
+ mxDestroyArray(lhs[1]);
+ mxDestroyArray(lhs[2]);
+ ostringstream tmp;
+ if (block > 1)
+ tmp << " in Solve_ByteCode_Sparse_GaussianElimination, singular system in block " << block+1 << "\n";
+ else
+ tmp << " in Solve_ByteCode_Sparse_GaussianElimination, singular system\n";
+ throw FatalExceptionHandling(tmp.str());
+}
+
+
+bool
SparseMatrix::Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, bool steady_state, int it_)
{
bool one;
@@ -2219,7 +2307,6 @@ SparseMatrix::Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, bool
}
if (piv_abs < eps)
{
- mexEvalString("drawnow;");
mxFree(piv_v);
mxFree(pivj_v);
mxFree(pivk_v);
@@ -2231,7 +2318,7 @@ SparseMatrix::Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, bool
mexPrintf("Error: singular system in Simulate_NG in block %d\n", blck+1);
else
mexPrintf("Error: singular system in Simulate_NG\n");
- return;
+ return true;
}
else
{
@@ -2411,6 +2498,7 @@ SparseMatrix::Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, bool
mxFree(pivk_v);
mxFree(NR);
mxFree(bc);
+ return false;
}
void
@@ -2898,13 +2986,14 @@ SparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bool
End_GE(Size);
}
-void
+bool
SparseMatrix::Simulate_Newton_One_Boundary(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, bool print_it, bool cvg, int &iter, bool steady_state, int stack_solve_algo, int solve_algo)
{
int i, j;
mxArray *b_m = NULL, *A_m = NULL, *x0_m = NULL;
Clear_u();
error_not_printed = true;
+ bool singular_system = false;
u_count_alloc_save = u_count_alloc;
if (isnan(res1) || isinf(res1) || (res2 > 12*g0 && iter > 0))
{
@@ -2931,7 +3020,7 @@ SparseMatrix::Simulate_Newton_One_Boundary(int blck, int y_size, int it_, int y_
else
mexPrintf(" dynare cannot improve the simulation in block %d at time %d (variable %d)\n", blck+1, it_+1, index_vara[max_res_idx]+1);
mexEvalString("drawnow;");
- return;
+ return singular_system;
}
else
{
@@ -2972,11 +3061,11 @@ SparseMatrix::Simulate_Newton_One_Boundary(int blck, int y_size, int it_, int y_
for (i = 0; i < y_size; i++)
y[i+it_*y_size] = ya[i+it_*y_size] + slowc_save*direction[i+it_*y_size];
iter--;
- return;
+ return singular_system;
}
if (cvg)
{
- return;
+ return singular_system;
}
if (print_it)
{
@@ -3024,7 +3113,7 @@ SparseMatrix::Simulate_Newton_One_Boundary(int blck, int y_size, int it_, int y_
}
bool zero_solution;
if ((solve_algo == 5 && steady_state) || (stack_solve_algo == 5 && !steady_state))
- Simple_Init(it_, y_kmin, y_kmax, Size, IM_i, zero_solution);
+ Simple_Init(Size, IM_i, zero_solution);
else
{
b_m = mxCreateDoubleMatrix(Size, 1, mxREAL);
@@ -3063,7 +3152,7 @@ SparseMatrix::Simulate_Newton_One_Boundary(int blck, int y_size, int it_, int y_
else
{
if ((solve_algo == 5 && steady_state) || (stack_solve_algo == 5 && !steady_state))
- Solve_ByteCode_Sparse_GaussianElimination(Size, blck, steady_state, it_);
+ singular_system = Solve_ByteCode_Sparse_GaussianElimination(Size, blck, steady_state, it_);
else if ((solve_algo == 7 && steady_state) || (stack_solve_algo == 2 && !steady_state))
Solve_Matlab_GMRES(A_m, b_m, Size, slowc, blck, false, it_, steady_state, x0_m);
else if ((solve_algo == 8 && steady_state) || (stack_solve_algo == 3 && !steady_state))
@@ -3071,7 +3160,7 @@ SparseMatrix::Simulate_Newton_One_Boundary(int blck, int y_size, int it_, int y_
else if ((solve_algo == 6 && steady_state) || ((stack_solve_algo == 0 || stack_solve_algo == 1) && !steady_state))
Solve_Matlab_LU_UMFPack(A_m, b_m, Size, slowc, false, it_);
}
- return;
+ return singular_system;
}
void
diff --git a/mex/sources/bytecode/SparseMatrix.hh b/mex/sources/bytecode/SparseMatrix.hh
index dc7bb3861d28396df2b020cc3f994707f4b1b3f7..e6df80e9492ae7c159b6b6f81009ba104d648d88 100644
--- a/mex/sources/bytecode/SparseMatrix.hh
+++ b/mex/sources/bytecode/SparseMatrix.hh
@@ -63,21 +63,22 @@ class SparseMatrix : public ErrorMsg
public:
SparseMatrix();
void Simulate_Newton_Two_Boundaries(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it, bool cvg, int &iter, int minimal_solving_periods, int stack_solve_algo, unsigned int endo_name_length, char *P_endo_names) /*throw(ErrorHandlingException)*/;
- void Simulate_Newton_One_Boundary(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, bool print_it, bool cvg, int &iter, bool steady_state, int stack_solve_algo, int solve_algo);
+ bool Simulate_Newton_One_Boundary(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, bool print_it, bool cvg, int &iter, bool steady_state, int stack_solve_algo, int solve_algo);
void Direct_Simulate(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it, int iter);
void fixe_u(double **u, int u_count_int, int max_lag_plus_max_lead_plus_1);
void Read_SparseMatrix(string file_name, const int Size, int periods, int y_kmin, int y_kmax, bool steady_state, bool two_boundaries, int stack_solve_algo, int solve_algo);
void Read_file(string file_name, int periods, int u_size1, int y_size, int y_kmin, int y_kmax, int &nb_endo, int &u_count, int &u_count_init, double *u);
+ void Singular_display(int block, int Size, bool steady_state, it_code_type it_code);
double g0, gp0, glambda2, try_at_iteration;
private:
void Init_GE(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM);
void Init_Matlab_Sparse(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM, mxArray *A_m, mxArray *b_m, mxArray *x0_m);
void Init_Matlab_Sparse_Simple(int Size, map<pair<pair<int, int>, int>, int> &IM, mxArray *A_m, mxArray *b_m, bool &zero_solution, mxArray *x0_m);
- void Simple_Init(int it_, int y_kmin, int y_kmax, int Size, std::map<std::pair<std::pair<int, int>, int>, int> &IM, bool &zero_solution);
+ void Simple_Init(int Size, std::map<std::pair<std::pair<int, int>, int>, int> &IM, bool &zero_solution);
void End_GE(int Size);
void Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bool symbolic, int Block_number);
- void Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, bool steady_state, int it_);
+ bool Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, bool steady_state, int it_);
void Solve_Matlab_Relaxation(mxArray *A_m, mxArray *b_m, unsigned int Size, double slowc_l, bool is_two_boundaries, int it_);
void Solve_Matlab_LU_UMFPack(mxArray *A_m, mxArray *b_m, int Size, double slowc_l, bool is_two_boundaries, int it_);
void Solve_Matlab_GMRES(mxArray *A_m, mxArray *b_m, int Size, double slowc, int block, bool is_two_boundaries, int it_, bool steady_state, mxArray *x0_m);
diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy
index 6a5a2055b873cebe307af96a71ee7a4932a564c0..40857352d314936f48c81608ed7f7e783bc19946 100644
--- a/preprocessor/DynareBison.yy
+++ b/preprocessor/DynareBison.yy
@@ -100,6 +100,7 @@ class ParsingDriver;
%token END ENDVAL EQUAL ESTIMATION ESTIMATED_PARAMS ESTIMATED_PARAMS_BOUNDS ESTIMATED_PARAMS_INIT
%token FILENAME FILTER_STEP_AHEAD FILTERED_VARS FIRST_OBS LAST_OBS SET_TIME
%token <string_val> FLOAT_NUMBER
+%token DEFAULT FIXED_POINT
%token FORECAST K_ORDER_SOLVER INSTRUMENTS PRIOR SHIFT MEAN STDEV VARIANCE MODE INTERVAL SHAPE DOMAINN
%token GAMMA_PDF GRAPH CONDITIONAL_VARIANCE_DECOMPOSITION NOCHECK STD
%token HISTVAL HOMOTOPY_SETUP HOMOTOPY_MODE HOMOTOPY_STEPS HP_FILTER HP_NGRID
@@ -108,7 +109,7 @@ class ParsingDriver;
%token <string_val> DATE_NUMBER
%token INV_GAMMA_PDF INV_GAMMA1_PDF INV_GAMMA2_PDF IRF IRF_SHOCKS
%token KALMAN_ALGO KALMAN_TOL SUBSAMPLES OPTIONS
-%token LABELS LAPLACE LIK_ALGO LIK_INIT LINEAR LOAD_IDENT_FILES LOAD_MH_FILE LOAD_PARAMS_AND_STEADY_STATE LOGLINEAR
+%token LABELS LAPLACE LIK_ALGO LIK_INIT LINEAR LOAD_IDENT_FILES LOAD_MH_FILE LOAD_PARAMS_AND_STEADY_STATE LOGLINEAR LYAPUNOV
%token MARKOWITZ MARGINAL_DENSITY MAX MAXIT
%token MFS MH_DROP MH_INIT_SCALE MH_JSCALE MH_MODE MH_NBLOCKS MH_REPLIC MH_RECOVER MIN MINIMAL_SOLVING_PERIODS
%token MODE_CHECK MODE_COMPUTE MODE_FILE MODEL MODEL_COMPARISON MODEL_INFO MSHOCKS ABS SIGN
@@ -123,7 +124,7 @@ class ParsingDriver;
%token QZ_CRITERIUM FULL DSGE_VAR DSGE_VARLAG DSGE_PRIOR_WEIGHT
%token RELATIVE_IRF REPLIC RPLOT SAVE_PARAMS_AND_STEADY_STATE
%token SHOCKS SHOCK_DECOMPOSITION SIGMA_E SIMUL SIMUL_ALGO SIMUL_SEED SMOOTHER STACK_SOLVE_ALGO STEADY_STATE_MODEL SOLVE_ALGO
-%token STDERR STEADY STOCH_SIMUL
+%token STDERR STEADY STOCH_SIMUL SYLVESTER
%token TEX RAMSEY_POLICY PLANNER_DISCOUNT DISCRETIONARY_POLICY
%token <string_val> TEX_NAME
%token UNIFORM_PDF UNIT_ROOT_VARS USE_DLL USEAUTOCORR GSA_SAMPLE_FILE
@@ -694,8 +695,8 @@ svar_identification_elem : EXCLUSION LAG INT_NUMBER ';' svar_equation_list
{ driver.add_constants_exclusion(); }
| RESTRICTION EQUATION INT_NUMBER COMMA
{ driver.add_restriction_equation_nbr($3);}
- restriction_expression EQUAL
- {driver.add_restriction_equal();}
+ restriction_expression EQUAL
+ {driver.add_restriction_equal();}
restriction_expression ';'
| UPPER_CHOLESKY ';'
{ driver.add_upper_cholesky(); }
@@ -925,6 +926,7 @@ stoch_simul_options : o_dr_algo
| o_conditional_variance_decomposition
| o_k_order_solver
| o_pruning
+ | o_sylvester
;
symbol_list : symbol_list symbol
@@ -1339,6 +1341,8 @@ estimation_options : o_datafile
| o_cova_compute
| o_irf_shocks
| o_sub_draws
+ | o_sylvester
+ | o_lyapunov
;
list_optim_option : QUOTED_STRING COMMA QUOTED_STRING
@@ -2072,6 +2076,10 @@ o_aim_solver: AIM_SOLVER {driver.option_num("aim_solver", "1"); };
o_partial_information : PARTIAL_INFORMATION {driver.option_num("partial_information", "1"); };
o_sub_draws: SUB_DRAWS EQUAL INT_NUMBER {driver.option_num("sub_draws",$3);};
o_planner_discount : PLANNER_DISCOUNT EQUAL expression { driver.declare_optimal_policy_discount_factor_parameter($3); };
+o_sylvester : SYLVESTER EQUAL FIXED_POINT {driver.option_num("sylvester_fp", "1"); }
+ | SYLVESTER EQUAL DEFAULT {driver.option_num("sylvester_fp", "0"); };
+o_lyapunov : LYAPUNOV EQUAL FIXED_POINT {driver.option_num("lyapunov_fp", "1"); }
+ | LYAPUNOV EQUAL DEFAULT {driver.option_num("lyapunov_fp", "0"); };
o_bvar_prior_tau : BVAR_PRIOR_TAU EQUAL signed_number { driver.option_num("bvar_prior_tau", $3); };
o_bvar_prior_decay : BVAR_PRIOR_DECAY EQUAL non_negative_number { driver.option_num("bvar_prior_decay", $3); };
diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll
index 2d91d096e1330661937395ceac979b689a4cf6e0..a5010b5bdc198b1b411b8b09f1c68e324504268a 100644
--- a/preprocessor/DynareFlex.ll
+++ b/preprocessor/DynareFlex.ll
@@ -290,6 +290,8 @@ string eofbuff;
<DYNARE_STATEMENT>dummy_obs {return token::DUMMY_OBS;}
<DYNARE_STATEMENT>nstates {return token::NSTATES;}
<DYNARE_STATEMENT>indxscalesstates {return token::INDXSCALESSTATES;}
+<DYNARE_STATEMENT>fixed_point {return token::FIXED_POINT;}
+<DYNARE_STATEMENT>default {return token::DEFAULT;}
<DYNARE_STATEMENT>alpha {
yylval->string_val = new string(yytext);
return token::ALPHA;
@@ -481,6 +483,8 @@ string eofbuff;
<DYNARE_STATEMENT>stack_solve_algo {return token::STACK_SOLVE_ALGO;}
<DYNARE_STATEMENT>drop {return token::DROP;}
<DYNARE_STATEMENT>order {return token::ORDER;}
+<DYNARE_STATEMENT>sylvester {return token::SYLVESTER;}
+<DYNARE_STATEMENT>lyapunov {return token::LYAPUNOV;}
<DYNARE_STATEMENT>replic {return token::REPLIC;}
<DYNARE_STATEMENT>ar {return token::AR;}
<DYNARE_STATEMENT>nofunctions {return token::NOFUNCTIONS;}
diff --git a/preprocessor/ParsingDriver.cc b/preprocessor/ParsingDriver.cc
index 80e01abd517a18ad095e0fc8e1e2ce377e557e27..db60330e4a0a308531d3c59158a5777e0700b8cc 100644
--- a/preprocessor/ParsingDriver.cc
+++ b/preprocessor/ParsingDriver.cc
@@ -2264,7 +2264,7 @@ ParsingDriver::add_model_var_or_external_function(string *function_name, bool in
else
{ //First time encountering this external function i.e., not previously declared or encountered
if (in_model_block)
- error("To use an external function within the model block, you must first declare it via the external_function() statement.");
+ error("To use an external function (" + *function_name + ") within the model block, you must first declare it via the external_function() statement.");
declare_symbol(function_name, eExternalFunction, NULL);
current_external_function_options.nargs = stack_external_function_args.top().size();