diff --git a/doc/manual.xml b/doc/manual.xml index 06f426a3a..d19391db8 100644 --- a/doc/manual.xml +++ b/doc/manual.xml @@ -1836,7 +1836,7 @@ steady(homotopy_mode = 1, homotopy_steps = 50); 0: Newton method to solve simultaneously all the equations for every period, see . (Default) 1: use a Newton algorithm with a sparse LU solver at each iteration. - 2: use a Newton algorithm with a Generalized Minimal Residual (GMRES) solver at each iteration. + 2: use a Newton algorithm with a Generalized Minimal Residual (GMRES) solver at each iteration. This option is not available under Octave. 3: use a Newton algorithm with a Stabilized Bi-Conjugate Gradient (BICGSTAB) solver at each iteration. 4: use a Newton algorithm with a optimal path length at each iteration. 5: use a Newton algorithm with a sparse Gaussian elimination (SPE) solver at each iteration. diff --git a/matlab/lnsrch1_wrapper_one_boundary.m b/matlab/lnsrch1_wrapper_one_boundary.m new file mode 100644 index 000000000..fcc5c465a --- /dev/null +++ b/matlab/lnsrch1_wrapper_one_boundary.m @@ -0,0 +1,43 @@ +function r = lnsrch1_wrapper_one_boundary(ya, y_index, fname, y, x, params, it_) +% wrapper for solve_one_boundary m-file when it is used with a dynamic +% model +% +% INPUTS +% ya [vector] The endogenous of the current block +% y_index [vector of int] The index of the endogenous variables of +% the block +% fname [string] name of the file containing the block +% to simulate +% y [matrix] All the endogenous variables of the model +% x [matrix] All the exogenous variables of the model +% params [vector] All the parameters of the model +% OUTPUTS +% r [vector] The residuals of the current block +% +% ALGORITHM +% none. +% +% SPECIAL REQUIREMENTS +% none. +% + +% Copyright (C) 2009 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 0) [L1, U1]=luinc(g1,luinc_tol); - [za,flag1] = gmres(g1,-r,Blck_size,1e-6,Blck_size,L1,U1); - %[za,flag1] = gmres(-g1,b',Blck_size,1e-6,Blck_size,L1,U1); + [dx,flag1] = gmres(g1,-r,Blck_size,1e-6,Blck_size,L1,U1); if (flag1>0 | reduced) if(flag1==1) disp(['Error in simul: No convergence inside GMRES after ' num2str(iter,'%6d') ' iterations, in block' num2str(Block_Num,'%3d')]); @@ -280,7 +284,6 @@ for it_=start:incr:finish luinc_tol = luinc_tol/10; reduced = 0; else - dx = za - ya; ya = ya + lambda*dx; if(is_dynamic) y(it_,y_index_eq) = ya'; @@ -293,8 +296,7 @@ for it_=start:incr:finish flag1=1; while(flag1>0) [L1, U1]=luinc(g1,luinc_tol); - [za,flag1] = bicgstab(g1,-r,1e-7,Blck_size,L1,U1); - %[za,flag1] = bicgstab(-g1,b',1e-7,Blck_size,L1,U1); + [dx,flag1] = bicgstab(g1,-r,1e-7,Blck_size,L1,U1); if (flag1>0 | reduced) if(flag1==1) disp(['Error in simul: No convergence inside BICGSTAB after ' num2str(iter,'%6d') ' iterations, in block' num2str(Block_Num,'%3d')]); @@ -306,7 +308,6 @@ for it_=start:incr:finish luinc_tol = luinc_tol/10; reduced = 0; else - dx = za - ya; ya = ya + lambda*dx; if(is_dynamic) y(it_,y_index_eq) = ya'; diff --git a/matlab/solve_two_boundaries.m b/matlab/solve_two_boundaries.m index 9b54d530c..fa7dc79a2 100644 --- a/matlab/solve_two_boundaries.m +++ b/matlab/solve_two_boundaries.m @@ -213,14 +213,6 @@ while ~(cvg==1 | iter>maxit_), dx = g1a\b- ya; ya = ya + lambda*dx; y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)'; - % v = ''; - % for i=1:(size(y_index,2)) - % v = [v ' %1.6f']; - % end; - % v = [v '\n']; - % v - % sprintf(v,y(:,y_index)') - % return; elseif(stack_solve_algo==2), flag1=1; while(flag1>0) @@ -264,11 +256,20 @@ while ~(cvg==1 | iter>maxit_), end; end; elseif(stack_solve_algo==4), - error('SOLVE_TWO_BOUNDARIES: stack_solve_algo=4 not implemented') - end; + ra = reshape(r(:, y_kmin+1:periods+y_kmin),periods*Blck_size, 1); + stpmx = 100 ; + stpmax = stpmx*max([sqrt(ya'*ya);size(y_index,2)]); + nn=1:size(ra,1); + g = (ra'*g1a)'; + f = 0.5*ra'*ra; + p = -g1a\ra; + [yn,f,ra,check]=lnsrch1(ya,f,g,p,stpmax,'lnsrch1_wrapper_two_boundaries',nn,nn, fname, y, y_index, x, params, periods, y_kmin, Blck_size); + dx = ya - yn; + y(1+y_kmin:periods+y_kmin,y_index)=reshape(yn',length(y_index),periods)'; + end end iter=iter+1; - disp(['iteration: ' num2str(iter,'%d') ' error: ' num2str(max_res,'%e')]); + disp(['iteration: ' num2str(iter,'%d') ' error: ' num2str(max_res,'%e') ' stack_solve_algo=' num2str(stack_solve_algo)]); end; if (iter>maxit_) disp(['No convergence after ' num2str(iter,'%4d') ' iterations in Block ' num2str(Block_Num,'%d')]); diff --git a/preprocessor/DynamicModel.cc b/preprocessor/DynamicModel.cc index 26b7c2c91..1bfb09d63 100644 --- a/preprocessor/DynamicModel.cc +++ b/preprocessor/DynamicModel.cc @@ -1995,7 +1995,7 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative evaluateAndReduceJacobian(eval_context, contemporaneous_jacobian, static_jacobian, dynamic_jacobian, cutoff, false); - computePossiblySingularNormalization(contemporaneous_jacobian, cutoff == 0); + computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian, dynamic_jacobian); computePrologueAndEpilogue(static_jacobian, equation_reordered, variable_reordered, prologue, epilogue); diff --git a/preprocessor/ModelTree.cc b/preprocessor/ModelTree.cc index 0b2853159..d7d6f1e28 100644 --- a/preprocessor/ModelTree.cc +++ b/preprocessor/ModelTree.cc @@ -32,8 +32,8 @@ using namespace boost; using namespace MFS; -void -ModelTree::computeNormalization(const set > &endo_eqs_incidence) throw (NormalizationException) +bool +ModelTree::computeNormalization(const jacob_map &contemporaneous_jacobian, bool verbose) { const int n = equation_number(); @@ -50,8 +50,8 @@ ModelTree::computeNormalization(const set > &endo_eqs_incidence) // Fill in the graph set > endo; - for (set >::const_iterator it = endo_eqs_incidence.begin(); it != endo_eqs_incidence.end(); it++) - add_edge(it->first + n, it->second, g); + for (jacob_map::const_iterator it = contemporaneous_jacobian.begin(); it != contemporaneous_jacobian.end(); it++) + add_edge(it->first.first + n, it->first.second, g); // Compute maximum cardinality matching vector mate_map(2*n); @@ -125,62 +125,95 @@ ModelTree::computeNormalization(const set > &endo_eqs_incidence) // Check if all variables are normalized vector::const_iterator it = find(mate_map.begin(), mate_map.begin() + n, graph_traits::null_vertex()); if (it != mate_map.begin() + n) - throw NormalizationException(symbol_table.getID(eEndogenous, it - mate_map.begin())); + { + if (verbose) + cerr << "ERROR: Could not normalize the model. Variable " + << symbol_table.getName(symbol_table.getID(eEndogenous, it - mate_map.begin())) + << " is not in the maximum cardinality matching." << endl; + check = false; + } + return check; } void -ModelTree::computePossiblySingularNormalization(const jacob_map &contemporaneous_jacobian, bool try_symbolic) +ModelTree::computeNonSingularNormalization(jacob_map &contemporaneous_jacobian, double cutoff, jacob_map &static_jacobian, dynamic_jacob_map &dynamic_jacobian) { + bool check = false; + cout << "Normalizing the model..." << endl; - set > endo_eqs_incidence; + int n = equation_number(); - for (jacob_map::const_iterator it = contemporaneous_jacobian.begin(); - it != contemporaneous_jacobian.end(); it++) - endo_eqs_incidence.insert(make_pair(it->first.first, it->first.second)); + // compute the maximum value of each row of the contemporaneous Jacobian matrix + //jacob_map normalized_contemporaneous_jacobian; + jacob_map normalized_contemporaneous_jacobian(contemporaneous_jacobian); + vector max_val(n, 0.0); + for (jacob_map::const_iterator iter = contemporaneous_jacobian.begin(); iter != contemporaneous_jacobian.end(); iter++) + if (fabs(iter->second) > max_val[iter->first.first]) + max_val[iter->first.first] = fabs(iter->second); - try + for (jacob_map::iterator iter = normalized_contemporaneous_jacobian.begin(); iter != normalized_contemporaneous_jacobian.end(); iter++) + iter->second /= max_val[iter->first.first]; + + //We start with the highest value of the cutoff and try to normalize the model + double current_cutoff = 0.99999999; + + int suppressed = 0; + while (!check && current_cutoff > 1e-19) { - computeNormalization(endo_eqs_incidence); - return; - } - catch (NormalizationException &e) - { - if (try_symbolic) - cout << "Normalization failed with cutoff, trying symbolic normalization..." << endl; - else + jacob_map tmp_normalized_contemporaneous_jacobian; + int suppress = 0; + for (jacob_map::iterator iter = normalized_contemporaneous_jacobian.begin(); iter != normalized_contemporaneous_jacobian.end(); iter++) + if (fabs(iter->second) > max(current_cutoff, cutoff)) + tmp_normalized_contemporaneous_jacobian[make_pair(iter->first.first, iter->first.second)] = iter->second; + else + suppress++; + + if (suppress != suppressed) + check = computeNormalization(tmp_normalized_contemporaneous_jacobian, false); + suppressed = suppress; + if (!check) { - cerr << "ERROR: Could not normalize the model. Variable " - << symbol_table.getName(e.symb_id) - << " is not in the maximum cardinality matching. Try to decrease the cutoff." << endl; - exit(EXIT_FAILURE); + current_cutoff /= 2; + // In this last case try to normalize with the complete jacobian + if (current_cutoff <= 1e-19) + check = computeNormalization(normalized_contemporaneous_jacobian, false); } } - // If no non-singular normalization can be found, try to find a normalization even with a potential singularity - if (try_symbolic) + if (!check) { - endo_eqs_incidence.clear(); + cout << "Normalization failed with cutoff, trying symbolic normalization..." << endl; + //if no non-singular normalization can be found, try to find a normalization even with a potential singularity + jacob_map tmp_normalized_contemporaneous_jacobian; set > endo; - for (int i = 0; i < equation_number(); i++) + for (int i = 0; i < n; i++) { endo.clear(); equations[i]->collectEndogenous(endo); for (set >::const_iterator it = endo.begin(); it != endo.end(); it++) - endo_eqs_incidence.insert(make_pair(i, it->first)); + tmp_normalized_contemporaneous_jacobian[make_pair(i, it->first)] = 1; } + check = computeNormalization(tmp_normalized_contemporaneous_jacobian, true); + if (check) + { + // Update the jacobian matrix + for (jacob_map::const_iterator it = tmp_normalized_contemporaneous_jacobian.begin(); it != tmp_normalized_contemporaneous_jacobian.end(); it++) + { + if (static_jacobian.find(make_pair(it->first.first, it->first.second)) == static_jacobian.end()) + static_jacobian[make_pair(it->first.first, it->first.second)] = 0; + if (dynamic_jacobian.find(make_pair(0, make_pair(it->first.first, it->first.second))) == dynamic_jacobian.end()) + dynamic_jacobian[make_pair(0, make_pair(it->first.first, it->first.second))] = 0; + if (contemporaneous_jacobian.find(make_pair(it->first.first, it->first.second)) == contemporaneous_jacobian.end()) + contemporaneous_jacobian[make_pair(it->first.first, it->first.second)] = 0; + } + } + } - try - { - computeNormalization(endo_eqs_incidence); - } - catch (NormalizationException &e) - { - cerr << "ERROR: Could not normalize the model even with zero cutoff. Variable " - << symbol_table.getName(e.symb_id) - << " is not in the maximum cardinality matching." << endl; - exit(EXIT_FAILURE); - } + if (!check) + { + cerr << "No normalization could be computed. Aborting." << endl; + exit(EXIT_FAILURE); } } diff --git a/preprocessor/ModelTree.hh b/preprocessor/ModelTree.hh index 7b1d7036a..4614a12c6 100644 --- a/preprocessor/ModelTree.hh +++ b/preprocessor/ModelTree.hh @@ -136,24 +136,21 @@ protected: //! for each block contains pair< max_lag, max_lead> t_lag_lead_vector block_lag_lead; - //! Exception thrown when normalization fails - class NormalizationException - { - public: - //! A variable missing from the maximum cardinal matching - int symb_id; - NormalizationException(int symb_id_arg) : symb_id(symb_id_arg) - { - } - }; - //! Compute the matching between endogenous and variable using the jacobian contemporaneous_jacobian - /*! \param endo_eqs_incidence A set indicating which endogenous appear in which equation. First element of pairs is equation number, second is type specific endo ID */ - void computeNormalization(const set > &endo_eqs_incidence) throw (NormalizationException); + /*! + \param contemporaneous_jacobian Jacobian used as an incidence matrix: all elements declared in the map (even if they are zero), are used as vertices of the incidence matrix + \return True if a complete normalization has been achieved + */ + bool computeNormalization(const jacob_map &contemporaneous_jacobian, bool verbose); + //! Try to compute the matching between endogenous and variable using a decreasing cutoff - /*! applied to the jacobian contemporaneous_jacobian and stop when a matching is found.*/ - /*! if no matching is found with a cutoff close to zero an error message is printout */ - void computePossiblySingularNormalization(const jacob_map &contemporaneous_jacobian, bool try_symbolic); + /*! + Applied to the jacobian contemporaneous_jacobian and stop when a matching is found. + If no matching is found using a strictly positive cutoff, then a zero cutoff is applied (i.e. use a symbolic normalization); in that case, the method adds zeros in the jacobian matrices to reflect all the edges in the symbolic incidence matrix. + If no matching is found with a zero cutoff close to zero an error message is printout. + */ + void computeNonSingularNormalization(jacob_map &contemporaneous_jacobian, double cutoff, jacob_map &static_jacobian, dynamic_jacob_map &dynamic_jacobian); + //! Try to normalized each unnormalized equation (matched endogenous variable only on the LHS) void computeNormalizedEquations(multimap &endo2eqs) const; //! Evaluate the jacobian and suppress all the elements below the cutoff diff --git a/preprocessor/StaticModel.cc b/preprocessor/StaticModel.cc index 784403e72..6bce4e2b1 100644 --- a/preprocessor/StaticModel.cc +++ b/preprocessor/StaticModel.cc @@ -45,7 +45,7 @@ StaticModel::StaticModel(SymbolTable &symbol_table_arg, } void -StaticModel::compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, map_idx_type &map_idx) const +StaticModel::compileDerivative(ofstream &code_file, int eq, int symb_id, map_idx_type &map_idx) const { first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, symbol_table.getID(eEndogenous, symb_id))); if (it != first_derivatives.end()) @@ -539,7 +539,7 @@ StaticModel::writeModelEquationsCodeOrdered(const string file_name, const string { case SOLVE_BACKWARD_SIMPLE: case SOLVE_FORWARD_SIMPLE: - compileDerivative(code_file, getBlockEquationID(block, 0), getBlockVariableID(block, 0), 0, map_idx); + compileDerivative(code_file, getBlockEquationID(block, 0), getBlockVariableID(block, 0), map_idx); { FSTPG_ fstpg(0); fstpg.write(code_file); @@ -727,7 +727,7 @@ StaticModel::computingPass(const eval_context_type &eval_context, bool no_tmp_te evaluateAndReduceJacobian(eval_context, contemporaneous_jacobian, static_jacobian, dynamic_jacobian, cutoff, false); - computePossiblySingularNormalization(contemporaneous_jacobian, cutoff == 0); + computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian, dynamic_jacobian); computePrologueAndEpilogue(static_jacobian, equation_reordered, variable_reordered, prologue, epilogue); diff --git a/preprocessor/StaticModel.hh b/preprocessor/StaticModel.hh index 5b82ab541..1dbfb07d1 100644 --- a/preprocessor/StaticModel.hh +++ b/preprocessor/StaticModel.hh @@ -74,7 +74,7 @@ private: void computeTemporaryTermsOrdered(); //! Write derivative code of an equation w.r. to a variable - void compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, map_idx_type &map_idx) const; + void compileDerivative(ofstream &code_file, int eq, int symb_id, map_idx_type &map_idx) const; //! Write chain rule derivative code of an equation w.r. to a variable void compileChainRuleDerivative(ofstream &code_file, int eq, int var, int lag, map_idx_type &map_idx) const; diff --git a/tests/Makefile.am b/tests/Makefile.am index 3ead51951..41322bb31 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -13,6 +13,8 @@ OCTAVE_MODS = \ predetermined_variables.mod \ block_bytecode/fs2000_simk.mod \ block_bytecode/fs2000_lu.mod \ + block_bytecode/fs2000_bicgstab.mod \ + block_bytecode/fs2000_optpath.mod \ block_bytecode/fs2000_bytecode.mod \ block_bytecode/ramst.mod \ block_bytecode/ireland.mod \ @@ -48,7 +50,9 @@ MODS = $(OCTAVE_MODS) \ AIM/ls2003_2L0L.mod \ AIM/ls2003_2L0L_AIM.mod \ AIM/ls2003_2L2L.mod \ - AIM/ls2003_2L2L_AIM.mod + AIM/ls2003_2L2L_AIM.mod \ + block_bytecode/fs2000_gmres.mod \ + block_bytecode/ramst_a.mod EXTRA_DIST = $(MODS) \ run_test_octave.m \