diff --git a/matlab/resid.m b/matlab/resid.m index e88b20a82..90162f0c3 100644 --- a/matlab/resid.m +++ b/matlab/resid.m @@ -70,7 +70,8 @@ if options_.block && ~options_.bytecode z(idx) = r; end elseif options_.block && options_.bytecode - [z,check] = bytecode('evaluate','static'); + [check, z] = bytecode('evaluate','static'); + mexErrCheck('bytecode', check); else z = feval([M_.fname '_static'],... oo_.steady_state,... diff --git a/matlab/simul.m b/matlab/simul.m index 26c9bad3a..067a389a9 100644 --- a/matlab/simul.m +++ b/matlab/simul.m @@ -73,13 +73,15 @@ end if(options_.block) if(options_.bytecode) - oo_.endo_simul=bytecode('dynamic'); + [info, oo_.endo_simul] = bytecode('dynamic'); + mexErrCheck('bytecode', info); else eval([M_.fname '_dynamic']); end; else if(options_.bytecode) - oo_.endo_simul=bytecode('dynamic'); + [info, oo_.endo_simul]=bytecode('dynamic'); + mexErrCheck('bytecode', info); else if M_.maximum_endo_lead == 0 error('SIMUL: purely backward models are not supported') diff --git a/matlab/steady_.m b/matlab/steady_.m index 74195c931..501bbe88d 100644 --- a/matlab/steady_.m +++ b/matlab/steady_.m @@ -76,9 +76,10 @@ if options_.steadystate_flag check1 = 1; end elseif options_.block && options_.bytecode - [residuals, check1] = bytecode('evaluate','static',oo_.steady_state,... + [check1, residuals] = bytecode('evaluate','static',oo_.steady_state,... [oo_.exo_steady_state; ... oo_.exo_det_steady_state], M_.params, 1); + mexErrCheck('bytecode', check1); else check1 = 0; check1 = max(abs(feval([M_.fname '_static'],... @@ -122,7 +123,8 @@ elseif options_.block && ~options_.bytecode oo_.exo_det_steady_state], M_.params); end elseif options_.bytecode - [oo_.steady_state,check] = bytecode('static'); + [check, oo_.steady_state] = bytecode('static'); + mexErrCheck('bytecode', check); else [oo_.steady_state,check] = dynare_solve([M_.fname '_static'],... oo_.steady_state,... diff --git a/mex/sources/bytecode/ErrorHandling.hh b/mex/sources/bytecode/ErrorHandling.hh new file mode 100644 index 000000000..95142c6ea --- /dev/null +++ b/mex/sources/bytecode/ErrorHandling.hh @@ -0,0 +1,122 @@ +/* + * Copyright (C) 2007-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 . + */ + +#ifndef ERROR_HANDLING +#define ERROR_HANDLING + +#include +#include +#include + +using namespace std; + +const int NO_ERROR_ON_EXIT = 0; +const int ERROR_ON_EXIT = 1; +class GeneralExceptionHandling +{ + string ErrorMsg; +public: + GeneralExceptionHandling(string ErrorMsg_arg) : ErrorMsg(ErrorMsg_arg){}; + inline string + GetErrorMsg() + { + return ErrorMsg; + } + inline void + completeErrorMsg(string ErrorMsg_arg) + { + ErrorMsg += ErrorMsg_arg; + } +}; + +class FloatingPointExceptionHandling : public GeneralExceptionHandling +{ +public: + FloatingPointExceptionHandling(string value) : GeneralExceptionHandling(string("Floating point error in bytrecode: " + value)) + {}; +}; + + +class LogExceptionHandling : public FloatingPointExceptionHandling +{ + double value; +public: + LogExceptionHandling(double value_arg) : FloatingPointExceptionHandling("log(X)"), + value(value_arg) + { + ostringstream tmp; + tmp << " with X=" << value << "\n"; + completeErrorMsg(tmp.str()); + }; +}; + +class Log10ExceptionHandling : public FloatingPointExceptionHandling +{ + double value; +public: + Log10ExceptionHandling(double value_arg) : FloatingPointExceptionHandling("log10(X)"), + value(value_arg) + { + ostringstream tmp; + tmp << " with X=" << value << "\n"; + completeErrorMsg(tmp.str()); + }; +}; + +class DivideExceptionHandling : public FloatingPointExceptionHandling +{ + double value1, value2; + DivideExceptionHandling(double value1_arg, double value2_arg) : FloatingPointExceptionHandling("a/X"), + value1(value1_arg), + value2(value2_arg) + { + ostringstream tmp; + tmp << " with X=" << value2 << "\n"; + completeErrorMsg(tmp.str()); + }; +}; + + +class PowExceptionHandling : public FloatingPointExceptionHandling +{ + double value1, value2; +public: + PowExceptionHandling(double value1_arg, double value2_arg) : FloatingPointExceptionHandling("X^a"), + value1(value1_arg), + value2(value2_arg) + { + ostringstream tmp; + tmp << " with X=" << value1 << "\n"; + completeErrorMsg(tmp.str()); + }; +}; + + +class FatalExceptionHandling : public GeneralExceptionHandling +{ +public: + FatalExceptionHandling(string ErrorMsg_arg) : GeneralExceptionHandling("Fatal error in bytrecode:") + { + completeErrorMsg(ErrorMsg_arg); + }; + FatalExceptionHandling() : GeneralExceptionHandling("") + { + }; +}; +#endif diff --git a/mex/sources/bytecode/Interpreter.cc b/mex/sources/bytecode/Interpreter.cc index 362f51bf1..93aa194cd 100644 --- a/mex/sources/bytecode/Interpreter.cc +++ b/mex/sources/bytecode/Interpreter.cc @@ -199,60 +199,58 @@ Interpreter::error_location(bool evaluate, bool steady_state) return("???"); break; } - Error_loc << endl << add_underscore_to_fpe(print_expression(it_code_expr, evaluate)); + Error_loc << endl << add_underscore_to_fpe(" " + print_expression(it_code_expr, evaluate)); return(Error_loc.str()); } double -Interpreter::pow1(double a, double b, bool evaluate, bool steady_state) +Interpreter::pow1(double a, double b) { double r = pow_(a, b); if (isnan(r)) { - if (error_not_printed) - { - mexPrintf("--------------------------------------\nError: X^a with X=%5.15f\n in %s \n--------------------------------------\n", a, error_location(evaluate, steady_state).c_str()); - error_not_printed = false; - r = 0.0000000000000000000000001; - } res1 = NAN; - return r; + r = 0.0000000000000000000000001; + throw PowExceptionHandling(a, b); } return r; } double -Interpreter::log1(double a, bool evaluate, bool steady_state) +Interpreter::divide(double a, double b) { - double r = log(a); - if (isnan(r)) + double r = a/ b; + if (isinf(r)) { - if (error_not_printed) - { - mexPrintf("--------------------------------------\nError: log(X) with X=%5.15f\n in %s \n--------------------------------------\n", a, error_location(evaluate, steady_state).c_str()); - error_not_printed = false; - r = 0.0000000000000000000000001; - } res1 = NAN; - return r; + r = 1e70; + throw PowExceptionHandling(a, b); } return r; } double -Interpreter::log10_1(double a, bool evaluate, bool steady_state) +Interpreter::log1(double a) { double r = log(a); if (isnan(r)) { - if (error_not_printed) - { - mexPrintf("--------------------------------------\nError: log10(X) with X=%5.15f\n in %s \n--------------------------------------\n", a, error_location(evaluate, steady_state).c_str()); - error_not_printed = false; - r = 0.0000000000000000000000001; - } res1 = NAN; - return r; + r = -1e70; + throw LogExceptionHandling(a); + } + return r; +} + +double +Interpreter::log10_1(double a) +{ + double r = log(a); + if (isnan(r)) + { + res1 = NAN; + r = -1e70; + throw Log10ExceptionHandling(a); } return r; } @@ -381,8 +379,9 @@ Interpreter::print_expression(it_code_type it_code, bool evaluate) dvar3 = ((FNUMEXPR_ *) it_code->second)->get_dvariable2(); break; default: - mexPrintf("Dérivatives %d not implemented yet\n", it_code->first); - mexErrMsgTxt("end of bytecode\n"); + ostringstream tmp; + tmp << " in print_expression, derivatives " << it_code->first << " not implemented yet\n"; + throw FatalExceptionHandling(tmp.str()); } break; case FLDV: @@ -1087,9 +1086,9 @@ Interpreter::print_expression(it_code_type it_code, bool evaluate) go_on = false; break; default: - mexPrintf("Unknown opcode %d!! FENDEQU=%d\n", it_code->first, FENDEQU); - mexErrMsgTxt("End of bytecode"); - break; + ostringstream tmp; + tmp << " in print_expression, unknown opcode " << it_code->first << "!! FENDEQU=" << FENDEQU << "\n"; + throw FatalExceptionHandling(tmp.str()); } it_code++; } @@ -1354,11 +1353,6 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num, int si #ifdef DEBUG mexPrintf("FLDSV y[var=%d]\n",var); tmp_out << " y[" << var << "](" << y[var] << ")"; - if(var<0 || var>= y_size) - { - mexPrintf("y[%d]=",var); - mexErrMsgTxt("End of bytecode"); - } #endif if (evaluate) Stack.push(ya[var]); @@ -1699,8 +1693,9 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num, int si jacob_exo_det[eq + size*pos_col] = rr; break; default: - mexPrintf("Variable %d not used yet\n", EQN_type); - mexErrMsgTxt("end of bytecode\n"); + ostringstream tmp; + tmp << " in compute_block_time, variable " << EQN_type << " not used yet\n"; + throw FatalExceptionHandling(tmp.str()); } #ifdef DEBUG tmp_out << "=>"; @@ -1737,20 +1732,17 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num, int si #endif break; case oDivide: - double r; - r = v1 / v2; - if (isinf(r)) + double tmp; + try { - if (error_not_printed) - { - mexPrintf("--------------------------------------\n Error: Divide by zero with %5.15f/%5.15f\n in %s \n--------------------------------------\n", v1, v2,error_location(evaluate, steady_state).c_str()); - error_not_printed = false; - r = 1e70; - } - res1 = NAN; + tmp = divide(v1 , v2); + } + catch(FloatingPointExceptionHandling &fpeh) + { + mexPrintf("%s %s\n",fpeh.GetErrorMsg().c_str(),error_location(evaluate, steady_state).c_str()); go_on = false; } - Stack.push(r); + Stack.push(tmp); #ifdef DEBUG tmp_out << " |" << v1 << "/" << v2 << "|"; #endif @@ -1792,10 +1784,16 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num, int si #endif break; case oPower: - r = pow1(v1, v2, evaluate, steady_state); - if (isnan(res1)) - go_on = false; - Stack.push(r); + try + { + tmp = pow1(v1, v2); + } + catch(FloatingPointExceptionHandling &fpeh) + { + mexPrintf("%s %s\n",fpeh.GetErrorMsg().c_str(),error_location(evaluate, steady_state).c_str()); + go_on = false; + } + Stack.push(tmp); #ifdef DEBUG tmp_out << " |" << v1 << "^" << v2 << "|"; @@ -1839,17 +1837,34 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num, int si #endif break; case oLog: - Stack.push(log1(v1, evaluate, steady_state)); - if (isnan(res1)) - go_on = false; + double tmp; + try + { + tmp = log1(v1); + } + catch(FloatingPointExceptionHandling &fpeh) + { + mexPrintf("%s %s\n",fpeh.GetErrorMsg().c_str(),error_location(evaluate, steady_state).c_str()); + go_on = false; + } + Stack.push(tmp); + //if (isnan(res1)) + #ifdef DEBUG tmp_out << " |log(" << v1 << ")|"; #endif break; case oLog10: - Stack.push(log10_1(v1, evaluate, steady_state)); - if (isnan(res1)) - go_on = false; + try + { + tmp = log10_1(v1); + } + catch(FloatingPointExceptionHandling &fpeh) + { + mexPrintf("%s %s\n",fpeh.GetErrorMsg().c_str(),error_location(evaluate, steady_state).c_str()); + go_on = false; + } + Stack.push(tmp); #ifdef DEBUG tmp_out << " |log10(" << v1 << ")|"; #endif @@ -1936,7 +1951,8 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num, int si Stack.push(erf(v1)); #ifdef DEBUG tmp_out << " |erf(" << v1 << ")|"; -#endif + +# endif break; default: ; @@ -2004,14 +2020,15 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num, int si op = ((FOK_ *) it_code->second)->get_arg(); if (Stack.size() > 0) { - mexPrintf("error: Stack not empty!\n"); - mexErrMsgTxt("End of bytecode"); + ostringstream tmp; + tmp << " in compute_block_time, stack not empty\n"; + throw FatalExceptionHandling(tmp.str()); } break; default: - mexPrintf("Unknown opcode %d!! FENDEQU=%d\n", it_code->first, FENDEQU); - mexErrMsgTxt("End of bytecode"); - break; + ostringstream tmp; + tmp << " in compute_block_time, unknown opcode " << it_code->first << "\n"; + throw FatalExceptionHandling(tmp.str()); } it_code++; } @@ -2174,7 +2191,7 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam } } -bool +int Interpreter::simulate_a_block(const int size, const int type, string file_name, string bin_basename, bool Gaussian_Elimination, 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) { @@ -2245,23 +2262,23 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, cvg = (fabs(rr) < solve_tolf); if (cvg) continue; - y[Block_Contain[0].Variable] += -r[0]/g1[0]; - if (isinf(y[Block_Contain[0].Variable])) + + try { - if (error_not_printed) - { - mexPrintf("--------------------------------------\n Error: Divide by zero with %5.15f/%5.15f\nSingularity in block %d\n--------------------------------------\n", r[0], g1[0], block_num); - error_not_printed = false; - } - res1 = NAN; + y[Block_Contain[0].Variable] += - divide(r[0],g1[0]); + } + catch(FloatingPointExceptionHandling &fpeh) + { + mexPrintf("%s \n",fpeh.GetErrorMsg().c_str()); + mexPrintf(" Singularity in block %d", block_num+1); } iter++; } if (!cvg) { - mexPrintf("In Solve Forward simple, convergence not achieved in block %d, after %d iterations\n", Block_Count+1, iter); - mexPrintf("r[0]= %f\n", r[0]); - return false; + ostringstream tmp; + tmp << " in Solve Forward simple, convergence not achieved in block " << Block_Count+1 << ", after " << iter << " iterations\n"; + throw FatalExceptionHandling(tmp.str()); } } else @@ -2284,22 +2301,22 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, cvg = (fabs(rr) < solve_tolf); if (cvg) continue; - y[Per_y_+Block_Contain[0].Variable] += -r[0]/g1[0]; - if (isinf(y[Per_y_+Block_Contain[0].Variable])) + try { - if (error_not_printed) - { - mexPrintf("--------------------------------------\n Error: Divide by zero with %5.15f/%5.15f\nSingularity in block %d\n--------------------------------------\n", r[0], g1[0], block_num); - error_not_printed = false; - } - res1 = NAN; + y[Per_y_+Block_Contain[0].Variable] += - divide(r[0], g1[0]); + } + catch(FloatingPointExceptionHandling &fpeh) + { + mexPrintf("%s \n",fpeh.GetErrorMsg().c_str()); + mexPrintf(" Singularity in block %d", block_num+1); } iter++; } if (!cvg) { - mexPrintf("In Solve Forward simple, convergence not achieved in block %d, at time %d after %d iterations\n", Block_Count, it_, iter); - mexErrMsgTxt("End of bytecode"); + ostringstream tmp; + tmp << " in Solve Forward simple, convergence not achieved in block " << Block_Count+1 << ", at time " << it_ << ", after " << iter << " iterations\n"; + throw FatalExceptionHandling(tmp.str()); } } } @@ -2327,8 +2344,16 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, cvg = (fabs(rr) < solve_tolf); if (cvg) continue; - y[Block_Contain[0].Variable] += -r[0]/g1[0]; - if (isinf(y[Block_Contain[0].Variable])) + try + { + y[Block_Contain[0].Variable] += - divide(r[0], g1[0]); + } + catch(FloatingPointExceptionHandling &fpeh) + { + mexPrintf("%s \n",fpeh.GetErrorMsg().c_str()); + mexPrintf(" Singularity in block %d", block_num+1); + } + /*if (isinf(y[Block_Contain[0].Variable])) { if (error_not_printed) { @@ -2336,13 +2361,14 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, error_not_printed = false; } res1 = NAN; - } + }*/ iter++; } if (!cvg) { - mexPrintf("In Solve Backward simple, convergence not achieved in block %d, after %d iterations\n", Block_Count+1, iter); - return false; + ostringstream tmp; + tmp << " in Solve Backward simple, convergence not achieved in block " << Block_Count+1 << ", after " << iter << " iterations\n"; + throw FatalExceptionHandling(tmp.str()); } } else @@ -2365,8 +2391,17 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, cvg = (fabs(rr) < solve_tolf); if (cvg) continue; - y[Per_y_+Block_Contain[0].Variable] += -r[0]/g1[0]; - if (isinf(y[Per_y_+Block_Contain[0].Variable])) + try + { + y[Per_y_+Block_Contain[0].Variable] += - divide(r[0], g1[0]); + } + catch(FloatingPointExceptionHandling &fpeh) + { + mexPrintf("%s \n",fpeh.GetErrorMsg().c_str()); + mexPrintf(" Singularity in block %d", block_num+1); + } + + /*if (isinf(y[Per_y_+Block_Contain[0].Variable])) { if (error_not_printed) { @@ -2374,13 +2409,14 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, error_not_printed = false; } res1 = NAN; - } + }*/ iter++; } if (!cvg) { - mexPrintf("In Solve Backward simple, convergence not achieved in block %d, at time %d after %d iterations\n", Block_Count+1, it_, iter); - mexErrMsgTxt("End of bytecode"); + ostringstream tmp; + tmp << " in Solve Backward simple, convergence not achieved in block " << Block_Count+1 << ", at time " << it_ << ", after " << iter << " iterations\n"; + throw FatalExceptionHandling(tmp.str()); } } } @@ -2435,7 +2471,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, if (cvg) continue; int prev_iter = iter; - result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, EQN_block_number, solve_algo); + Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, solve_algo); iter++; if (iter > prev_iter) { @@ -2446,8 +2482,9 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, } if (!cvg || !result) { - mexPrintf("In Solve Forward complete, convergence not achieved in block %d, after %d iterations\n", Block_Count+1, iter); - return false; + ostringstream tmp; + tmp << " in Solve Forward complete, convergence not achieved in block " << Block_Count+1 << ", after " << iter << " iterations\n"; + throw FatalExceptionHandling(tmp.str()); } } else @@ -2462,11 +2499,11 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, error_not_printed = true; compute_block_time(0, false, block_num, size, steady_state); cvg = false; - result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, EQN_block_number, solve_algo); + Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, solve_algo); if (!result) { - mexPrintf("In Solve Forward complete, convergence not achieved in block %d\n", Block_Count+1); - return false; + mexPrintf(" in Solve Forward complete, convergence not achieved in block %d\n", Block_Count+1); + return ERROR_ON_EXIT; } } } @@ -2514,7 +2551,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, if (cvg) continue; int prev_iter = iter; - result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, EQN_block_number, solve_algo); + Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, solve_algo); iter++; if (iter > prev_iter) { @@ -2525,8 +2562,9 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, } if (!cvg) { - mexPrintf("In Solve Forward complete, convergence not achieved in block %d, at time %d after %d iterations\n", Block_Count+1, it_, iter); - mexErrMsgTxt("End of bytecode"); + ostringstream tmp; + tmp << " in Solve Forward complete, convergence not achieved in block " << Block_Count+1 << ", at time " << it_ << ", after " << iter << " iterations\n"; + throw FatalExceptionHandling(tmp.str()); } } } @@ -2543,7 +2581,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, error_not_printed = true; compute_block_time(0, false, block_num, size, steady_state); cvg = false; - result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, EQN_block_number, solve_algo); + Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, solve_algo); } } } @@ -2601,7 +2639,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, if (cvg) continue; int prev_iter = iter; - result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, EQN_block_number, solve_algo); + Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, solve_algo); iter++; if (iter > prev_iter) { @@ -2612,8 +2650,9 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, } if (!cvg || !result) { - mexPrintf("In Solve Backward complete, convergence not achieved in block %d, after %d iterations\n", Block_Count+1, iter); - return false; + ostringstream tmp; + tmp << " in Solve Backward complete, convergence not achieved in block " << Block_Count+1 << ", after " << iter << " iterations\n"; + throw FatalExceptionHandling(tmp.str()); } } else @@ -2627,11 +2666,11 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, error_not_printed = true; compute_block_time(0, false, block_num, size, steady_state); cvg = false; - result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, EQN_block_number, solve_algo); + Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, solve_algo); if (!result) { - mexPrintf("In Solve Backward complete, convergence not achieved in block %d\n", Block_Count+1); - return false; + mexPrintf(" in Solve Backward complete, convergence not achieved in block %d\n", Block_Count+1); + return ERROR_ON_EXIT; } } } @@ -2679,7 +2718,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, if (cvg) continue; int prev_iter = iter; - result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, EQN_block_number, solve_algo); + Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, solve_algo); iter++; if (iter > prev_iter) { @@ -2690,8 +2729,9 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, } if (!cvg) { - mexPrintf("In Solve Backward complete, convergence not achieved in block %d, at time %d after %d iterations\n", Block_Count+1, it_, iter); - mexErrMsgTxt("End of bytecode"); + ostringstream tmp; + tmp << " in Solve Backward complete, convergence not achieved in block " << Block_Count+1 << ", at time " << it_ << ", after " << iter << " iterations\n"; + throw FatalExceptionHandling(tmp.str()); } } } @@ -2704,7 +2744,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, error_not_printed = true; compute_block_time(0, false, block_num, size, steady_state); cvg = false; - result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, EQN_block_number, solve_algo); + Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, solve_algo); } } } @@ -2723,7 +2763,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, if (steady_state) { mexPrintf("SOLVE_TWO_BOUNDARIES in a steady state model: impossible case\n"); - return false; + return ERROR_ON_EXIT; } fixe_u(&u, u_count_int, u_count_int); Read_SparseMatrix(bin_basename, size, periods, y_kmin, y_kmax, steady_state, true, stack_solve_algo, solve_algo); @@ -2731,9 +2771,6 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, r = (double *) mxMalloc(size*sizeof(double)); y_save = (double *) mxMalloc(y_size*sizeof(double)*(periods+y_kmax+y_kmin)); begining = it_code; - if (!Gaussian_Elimination) - { - } giter = 0; iter = 0; if (!is_linear) @@ -2784,7 +2821,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, cvg = (max_res < solve_tolf); u_count = u_count_saved; int prev_iter = iter; - simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter, minimal_solving_periods, EQN_block_number, stack_solve_algo, endo_name_length, P_endo_names); + Simulate_Newton_Two_Boundaries(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter, minimal_solving_periods, stack_solve_algo, endo_name_length, P_endo_names); iter++; if (iter > prev_iter) { @@ -2796,8 +2833,9 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, } if (!cvg) { - mexPrintf("In Solve two boundaries, convergence not achieved in block %d, after %d iterations\n", Block_Count+1, iter); - mexErrMsgTxt("End of bytecode"); + ostringstream tmp; + tmp << " in Solve two boundaries, convergence not achieved in block " << Block_Count+1 << ", after " << iter << " iterations\n"; + throw FatalExceptionHandling(tmp.str()); } } else @@ -2825,7 +2863,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, } } cvg = false; - simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter, minimal_solving_periods, EQN_block_number, stack_solve_algo, endo_name_length, P_endo_names); + Simulate_Newton_Two_Boundaries(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter, minimal_solving_periods, stack_solve_algo, endo_name_length, P_endo_names); } mxFree(r); mxFree(y_save); @@ -2835,17 +2873,17 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, memset(direction, 0, size_of_direction); break; default: - mexPrintf("Unknown type =%d\n", type); - mexEvalString("drawnow;"); - mexErrMsgTxt("End of bytecode"); + ostringstream tmp; + tmp << " in simulate_a_block, Unknown type = " << type << "\n"; + throw FatalExceptionHandling(tmp.str()); + return ERROR_ON_EXIT; } - return true; + return NO_ERROR_ON_EXIT; } bool Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_state, bool evaluate, bool block, int &nb_blocks) { - bool result = true; int var; @@ -2859,11 +2897,9 @@ Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_s EQN_block_number = code.get_block_number(); if (!code_liste.size()) { - mexPrintf("%s.cod Cannot be opened\n", file_name.c_str()); - mexEvalString("drawnow;"); - filename += " stopped"; - mexEvalString("drawnow;"); - mexErrMsgTxt(filename.c_str()); + ostringstream tmp; + tmp << " in compute_blocks, " << file_name.c_str() << " cannot be opened\n"; + throw FatalExceptionHandling(tmp.str()); } //The big loop on intructions Block_Count = -1; @@ -2901,11 +2937,15 @@ Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_s fb->get_is_linear(), fb->get_endo_nbr(), fb->get_Max_Lag(), fb->get_Max_Lead(), fb->get_u_count_int()); } else - result = simulate_a_block(fb->get_size(), fb->get_type(), file_name, bin_basename, true, steady_state, Block_Count, - fb->get_is_linear(), fb->get_endo_nbr(), fb->get_Max_Lag(), fb->get_Max_Lead(), fb->get_u_count_int()); + { + result = simulate_a_block(fb->get_size(), fb->get_type(), file_name, bin_basename, true, steady_state, Block_Count, + fb->get_is_linear(), fb->get_endo_nbr(), fb->get_Max_Lag(), fb->get_Max_Lead(), fb->get_u_count_int()); + if (result == ERROR_ON_EXIT) + return ERROR_ON_EXIT; + } delete fb; } - if (!result) + if (result == ERROR_ON_EXIT) go_on = false; break; case FEND: @@ -2936,10 +2976,9 @@ Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_s it_code++; break; default: - mexPrintf("Unknown command %d\n",it_code->first); - mexEvalString("drawnow;"); - mexErrMsgTxt("End of bytecode"); - break; + ostringstream tmp; + tmp << " in compute_blocks, unknown command " << it_code->first << "\n"; + throw FatalExceptionHandling(tmp.str()); } } mxFree(Init_Code->second); diff --git a/mex/sources/bytecode/Interpreter.hh b/mex/sources/bytecode/Interpreter.hh index b7ae9367a..3b88194d6 100644 --- a/mex/sources/bytecode/Interpreter.hh +++ b/mex/sources/bytecode/Interpreter.hh @@ -57,9 +57,10 @@ private: int EQN_lag1, EQN_lag2, EQN_lag3; it_code_type it_code_expr; protected: - double pow1(double a, double b, bool evaluate, bool steady_state); - double log1(double a, bool evaluate, bool steady_state); - double log10_1(double a, bool evaluate, bool steady_state); + double pow1(double a, double b); + double divide(double a, double b); + double log1(double a); + double log10_1(double a); /*string remove_white(string str);*/ string add_underscore_to_fpe(const string &str); string get_variable(const SymbolType variable_type, const unsigned int variable_num); @@ -68,7 +69,7 @@ private: string print_expression(it_code_type it_code, bool evaluate); void evaluate_a_block(const int size, const int type, string bin_basename, bool steady_state, int block_num, const bool is_linear = false, const int symbol_table_endo_nbr = 0, const int Block_List_Max_Lag = 0, const int Block_List_Max_Lead = 0, const int u_count_int = 0); - bool simulate_a_block(const int size, const int type, string file_name, string bin_basename, bool Gaussian_Elimination, bool steady_state, int block_num, + int simulate_a_block(const int size, const int type, string file_name, string bin_basename, bool Gaussian_Elimination, bool steady_state, int block_num, const bool is_linear = false, const int symbol_table_endo_nbr = 0, const int Block_List_Max_Lag = 0, const int Block_List_Max_Lead = 0, const int u_count_int = 0); double *T; vector Block_Contain; diff --git a/mex/sources/bytecode/SparseMatrix.cc b/mex/sources/bytecode/SparseMatrix.cc index 3bf5154e8..d2088a354 100644 --- a/mex/sources/bytecode/SparseMatrix.cc +++ b/mex/sources/bytecode/SparseMatrix.cc @@ -289,12 +289,12 @@ SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, i SaveCode.open((file_name + "_dynamic.bin").c_str(), ios::in | ios::binary); if (!SaveCode.is_open()) { + ostringstream tmp; if (steady_state) - mexPrintf("Error : Can't open file \"%s\" for reading\n", (file_name + "_static.bin").c_str()); + tmp << " in Read_SparseMatrix, " << file_name << "_static.bin cannot be opened\n"; else - mexPrintf("Error : Can't open file \"%s\" for reading\n", (file_name + "_dynamic.bin").c_str()); - mexEvalString("st=fclose('all');clear all;"); - mexErrMsgTxt("Exit from Dynare"); + tmp << " in Read_SparseMatrix, " << file_name << "_dynamic.bin cannot be opened\n"; + throw FatalExceptionHandling(tmp.str()); } } IM_i.clear(); @@ -449,26 +449,30 @@ SparseMatrix::Init_Matlab_Sparse_Simple(int Size, map, int>, double *b = mxGetPr(b_m); if (!b) { - mexPrintf("Can't retrieve b vector in LU solver\n"); - mexErrMsgTxt("end of bytecode\n"); + ostringstream tmp; + tmp << " in Init_Matlab_Sparse_Simple, can't retrieve b vector\n"; + throw FatalExceptionHandling(tmp.str()); } mwIndex *Ai = mxGetIr(A_m); if (!Ai) { - mexPrintf("Can't allocate Ai index vector in LU solver\n"); - mexErrMsgTxt("end of bytecode\n"); + ostringstream tmp; + tmp << " in Init_Matlab_Sparse_Simple, can't allocate Ai index vector\n"; + throw FatalExceptionHandling(tmp.str()); } mwIndex *Aj = mxGetJc(A_m); if (!Aj) { - mexPrintf("Can't allocate Aj index vector in LU solver\n"); - mexErrMsgTxt("end of bytecode\n"); + ostringstream tmp; + tmp << " in Init_Matlab_Sparse_Simple, can't allocate Aj index vector\n"; + throw FatalExceptionHandling(tmp.str()); } double *A = mxGetPr(A_m); if (!A) { - mexPrintf("Can't retrieve A matrix in LU solver\n"); - mexErrMsgTxt("end of bytecode\n"); + ostringstream tmp; + tmp << " in Init_Matlab_Sparse_Simple, can't retrieve A matrix\n"; + throw FatalExceptionHandling(tmp.str()); } map, int>, int>::iterator it4; for (i = 0; i < y_size*(periods+y_kmin); i++) @@ -496,13 +500,15 @@ SparseMatrix::Init_Matlab_Sparse_Simple(int Size, map, int>, #if DEBUG if (index<0 || index >= u_count_alloc || index > Size + Size*Size) { - mexPrintf("index (%d) out of range for u vector (0) max %d allocated %d\n",index, Size+Size*Size, u_count_alloc); - mexErrMsgTxt("end of bytecode\n"); + ostringstream tmp; + tmp << " in Init_Matlab_Sparse_Simple, index (" << index << ") out of range for u vector max = " << Size+Size*Size << " allocated = " << u_count_alloc << "\n"; + throw FatalExceptionHandling(tmp.str()); } if (NZE >= max_nze) { - mexPrintf("Exceeds the capacity of A_m sparse matrix in LU (0)\n"); - mexErrMsgTxt("end of bytecode\n"); + ostringstream tmp; + tmp << " in Init_Matlab_Sparse_Simple, exceeds the capacity of A_m sparse matrix\n"; + throw FatalExceptionHandling(tmp.str()); } #endif A[NZE] = u[index]; @@ -511,18 +517,21 @@ SparseMatrix::Init_Matlab_Sparse_Simple(int Size, map, int>, #if DEBUG if (eq < 0 || eq >= Size) { - mexPrintf("index (%d) out of range for b vector (0)\n",eq); - mexErrMsgTxt("end of bytecode\n"); + ostringstream tmp; + tmp << " in Init_Matlab_Sparse_Simple, index (" << eq << ") out of range for b vector\n"; + throw FatalExceptionHandling(tmp.str()); } if (var < 0 || var >= Size) { - mexPrintf("index (%d) out of range for index_vara vector (0)\n",var); - mexErrMsgTxt("end of bytecode\n"); + ostringstream tmp; + tmp << " in Init_Matlab_Sparse_Simple, index (" << var << ") out of range for index_vara vector\n"; + throw FatalExceptionHandling(tmp.str()); } if (index_vara[var] < 0 || index_vara[var] >= y_size) { - mexPrintf("index (%d) out of range for y vector max=%d (0)\n",index_vara[var], y_size); - mexErrMsgTxt("end of bytecode\n"); + ostringstream tmp; + tmp << " in Init_Matlab_Sparse_Simple, index (" << index_vara[var] << ") out of range for y vector max=" << y_size << " (0)\n"; + throw FatalExceptionHandling(tmp.str()); } #endif it4++; @@ -539,26 +548,30 @@ SparseMatrix::Init_Matlab_Sparse(int periods, int y_kmin, int y_kmax, int Size, double *b = mxGetPr(b_m); if (!b) { - mexPrintf("Can't retrieve b vector in LU solver\n"); - mexErrMsgTxt("end of bytecode\n"); + ostringstream tmp; + tmp << " in Init_Matlab_Sparse, can't retrieve b vector\n"; + throw FatalExceptionHandling(tmp.str()); } mwIndex *Ai = mxGetIr(A_m); if (!Ai) { - mexPrintf("Can't allocate Ai index vector in LU solver\n"); - mexErrMsgTxt("end of bytecode\n"); + ostringstream tmp; + tmp << " in Init_Matlab_Sparse, can't allocate Ai index vector\n"; + throw FatalExceptionHandling(tmp.str()); } mwIndex *Aj = mxGetJc(A_m); if (!Aj) { - mexPrintf("Can't allocate Aj index vector in LU solver\n"); - mexErrMsgTxt("end of bytecode\n"); + ostringstream tmp; + tmp << " in Init_Matlab_Sparse, can't allocate Aj index vector\n"; + throw FatalExceptionHandling(tmp.str()); } double *A = mxGetPr(A_m); if (!A) { - mexPrintf("Can't retrieve A matrix in LU solver\n"); - mexErrMsgTxt("end of bytecode\n"); + ostringstream tmp; + tmp << " in Init_Matlab_Sparse, can't retrieve A matrix\n"; + throw FatalExceptionHandling(tmp.str()); } map, int>, int>::iterator it4; for (i = 0; i < y_size*(periods+y_kmin); i++) @@ -595,16 +608,18 @@ SparseMatrix::Init_Matlab_Sparse(int periods, int y_kmin, int y_kmax, int Size, if (lag <= ti_new_y_kmax && lag >= ti_new_y_kmin) /*Build the index for sparse matrix containing the jacobian : u*/ { #if DEBUG - if (index<0 || index >= u_count_alloc) - { - mexPrintf("index (%d) out of range for u vector (0)\n",index); - mexErrMsgTxt("end of bytecode\n"); - } - if (NZE >= max_nze) - { - mexPrintf("Exceeds the capacity of A_m sparse matrix in LU (0)\n"); - mexErrMsgTxt("end of bytecode\n"); - } + if (index<0 || index >= u_count_alloc || index > Size + Size*Size) + { + ostringstream tmp; + tmp << " in Init_Matlab_Sparse, index (" << index << ") out of range for u vector max = " << Size+Size*Size << " allocated = " << u_count_alloc << "\n"; + throw FatalExceptionHandling(tmp.str()); + } + if (NZE >= max_nze) + { + ostringstream tmp; + tmp << " in Init_Matlab_Sparse, exceeds the capacity of A_m sparse matrix\n"; + throw FatalExceptionHandling(tmp.str()); + } #endif A[NZE] = u[index]; Ai[NZE] = eq - lag * Size; @@ -613,25 +628,23 @@ SparseMatrix::Init_Matlab_Sparse(int periods, int y_kmin, int y_kmax, int Size, if (lag > ti_y_kmax || lag < ti_y_kmin) { #if DEBUG - if ((index+lag*u_count_init) < 0 || (index+lag*u_count_init) >= u_count_alloc) + if (eq < 0 || eq >= Size * periods) { - mexPrintf("index (%d) out of range for u vector (1)\n",index+lag*u_count_init); - mexErrMsgTxt("end of bytecode\n"); - } - if (eq < 0 || eq >= (Size*periods)) - { - mexPrintf("index (%d) out of range for b vector (0)\n",eq); - mexErrMsgTxt("end of bytecode\n"); - } + ostringstream tmp; + tmp << " in Init_Matlab_Sparse, index (" << eq << ") out of range for b vector\n"; + throw FatalExceptionHandling(tmp.str()); + } if (var+Size*(y_kmin+t+lag) < 0 || var+Size*(y_kmin+t+lag) >= Size*(periods+y_kmin+y_kmax)) { - mexPrintf("index (%d) out of range for index_vara vector (0)\n",var+Size*(y_kmin+t+lag)); - mexErrMsgTxt("end of bytecode\n"); + ostringstream tmp; + tmp << " in Init_Matlab_Sparse, index (" << var+Size*(y_kmin+t+lag) << ") out of range for index_vara vector\n"; + throw FatalExceptionHandling(tmp.str()); } if (index_vara[var+Size*(y_kmin+t+lag)] < 0 || index_vara[var+Size*(y_kmin+t+lag)] >= y_size*(periods+y_kmin+y_kmax)) { - mexPrintf("index (%d) out of range for y vector max=%d (0)\n",index_vara[var+Size*(y_kmin+t+lag)], y_size*(periods+y_kmin+y_kmax)); - mexErrMsgTxt("end of bytecode\n"); + ostringstream tmp; + tmp << " in Init_Matlab_Sparse, index (" << index_vara[var+Size*(y_kmin+t+lag)] << ") out of range for y vector max=" << y_size*(periods+y_kmin+y_kmax) << "\n"; + throw FatalExceptionHandling(tmp.str()); } #endif b[eq] += u[index+lag*u_count_init]*y[index_vara[var+Size*(y_kmin+t+lag)]]; @@ -642,13 +655,15 @@ SparseMatrix::Init_Matlab_Sparse(int periods, int y_kmin, int y_kmax, int Size, #if DEBUG if (index < 0 || index >= u_count_alloc) { - mexPrintf("index (%d) out of range for u vector (2)\n",index); - mexErrMsgTxt("end of bytecode\n"); + ostringstream tmp; + tmp << " in Init_Matlab_Sparse, index (" << index << ") out of range for u vector\n"; + throw FatalExceptionHandling(tmp.str()); } if (eq < 0 || eq >= (Size*periods)) { - mexPrintf("index (%d) out of range for b vector (1)\n",eq); - mexErrMsgTxt("end of bytecode\n"); + ostringstream tmp; + tmp << " in Init_Matlab_Sparse, index (" << eq << ") out of range for b vector\n"; + throw FatalExceptionHandling(tmp.str()); } #endif b[eq] += u[index]; @@ -793,9 +808,9 @@ SparseMatrix::Get_u() u = (double *) mxRealloc(u, u_count_alloc*sizeof(double)); if (!u) { - mexPrintf("Error in Get_u: memory exhausted (realloc(%d))\n", u_count_alloc*sizeof(double)); - mexEvalString("st=fclose('all');clear all;"); - mexErrMsgTxt("Exit from Dynare"); + ostringstream tmp; + tmp << " in Get_u, memory exhausted (realloc(" << u_count_alloc*sizeof(double) << "))\n"; + throw FatalExceptionHandling(tmp.str()); } int i = u_count; u_count++; @@ -881,11 +896,9 @@ SparseMatrix::compare(int *save_op, int *save_opa, int *save_opaa, int beg_t, in i += 3; break; default: - mexPrintf("unknown operator = %d ", save_op_s->operat); - mexEvalString("st=fclose('all');clear all;"); - filename += " stopped"; - mexErrMsgTxt(filename.c_str()); - break; + ostringstream tmp; + tmp << " in compare, unknown operator = " << save_op_s->operat << "\n"; + throw FatalExceptionHandling(tmp.str()); } j++; } @@ -907,9 +920,9 @@ SparseMatrix::compare(int *save_op, int *save_opa, int *save_opaa, int beg_t, in u = (double *) mxRealloc(u, u_count_alloc*sizeof(double)); if (!u) { - mexPrintf("Error in Get_u: memory exhausted (realloc(%d))\n", u_count_alloc*sizeof(double)); - mexEvalString("st=fclose('all');clear all;"); - mexErrMsgTxt("Exit from Dynare"); + ostringstream tmp; + tmp << " in compare, memory exhausted (realloc(" << u_count_alloc*sizeof(double) << "))\n"; + throw FatalExceptionHandling(tmp.str()); } } double *up; @@ -1199,417 +1212,6 @@ SparseMatrix::simple_bksub(int it_, int Size, double slowc_l) return res1; } -bool -SparseMatrix::simulate_NG(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 Block_number, int solve_algo) -{ - int i, j, k; - int pivj = 0, pivk = 0; - double piv_abs; - NonZeroElem *first, *firsta, *first_suba; - double *piv_v; - int *pivj_v, *pivk_v, *NR; - int l, N_max; - bool one; - mxArray *b_m, *A_m; - Clear_u(); - piv_v = (double *) mxMalloc(Size*sizeof(double)); - pivj_v = (int *) mxMalloc(Size*sizeof(int)); - pivk_v = (int *) mxMalloc(Size*sizeof(int)); - NR = (int *) mxMalloc(Size*sizeof(int)); - error_not_printed = true; - u_count_alloc_save = u_count_alloc; - if (isnan(res1) || isinf(res1) || (res2 > 12*g0 && iter>0)) - { - if (iter == 0) - { - for (j = 0; j < y_size; j++) - { - bool select = false; - for (int i = 0; i < Size; i++) - if (j == index_vara[i]) - { - select = true; - break; - } - if (select) - mexPrintf("-> variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]); - else - mexPrintf(" variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]); - } - mexPrintf("res1=%5.10f\n", res1); - mexPrintf("The initial values of endogenous variables are too far from the solution.\n"); - mexPrintf("Change them!\n"); - mexEvalString("drawnow;"); - mxFree(piv_v); - mxFree(pivj_v); - mxFree(pivk_v); - mxFree(NR); - if (steady_state) - return false; - else - { - mexEvalString("st=fclose('all');clear all;"); - filename += " stopped"; - mexErrMsgTxt(filename.c_str()); - } - } - if (fabs(slowc_save) < 1e-8) - { - for (j = 0; j < y_size; j++) - { - bool select = false; - for (int i = 0; i < Size; i++) - if (j == index_vara[i]) - { - select = true; - break; - } - if (select) - mexPrintf("-> variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]); - else - mexPrintf(" variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]); - } - mexPrintf("Dynare cannot improve the simulation in block %d at time %d (variable %d)\n", blck+1, it_+1, max_res_idx); - mexEvalString("drawnow;"); - mxFree(piv_v); - mxFree(pivj_v); - mxFree(pivk_v); - mxFree(NR); - if (steady_state) - return false; - else - { - mexEvalString("st=fclose('all');clear all;"); - filename += " stopped"; - mexErrMsgTxt(filename.c_str()); - } - } - if(!(isnan(res1) || isinf(res1)) && !(isnan(g0) || isinf(g0))) - { - if (try_at_iteration == 0) - { - prev_slowc_save = slowc_save; - slowc_save = max( - gp0 / (2 * (res2 - g0 - gp0)) , 0.1); - } - else - { - double t1 = res2 - gp0 * slowc_save - g0; - double t2 = glambda2 - gp0 * prev_slowc_save - g0; - double a = (1/(slowc_save * slowc_save) * t1 - 1/(prev_slowc_save * prev_slowc_save) * t2) / (slowc_save - prev_slowc_save); - double b = (-prev_slowc_save/(slowc_save * slowc_save) * t1 + slowc_save/(prev_slowc_save * prev_slowc_save) * t2) / (slowc_save - prev_slowc_save); - prev_slowc_save = slowc_save; - slowc_save = max(min( -b + sqrt(b*b - 3 * a * gp0) / (3 * a), 0.5 * slowc_save), 0.1 * slowc_save); - } - glambda2 = res2; - try_at_iteration ++; - } - else - { - prev_slowc_save = slowc_save; - slowc_save /= 1.1; - } - if (print_it) - mexPrintf("Error: Simulation diverging, trying to correct it using slowc=%f\n", slowc_save); - for (i = 0; i < y_size; i++) - y[i+it_*y_size] = ya[i+it_*y_size] + slowc_save*direction[i+it_*y_size]; - mxFree(piv_v); - mxFree(pivj_v); - mxFree(pivk_v); - mxFree(NR); - iter--; - return true; - } - if (cvg) - { - mxFree(piv_v); - mxFree(pivj_v); - mxFree(pivk_v); - mxFree(NR); - return (true); - } - if (print_it ) - { - mexPrintf("solwc=%f g0=%f res2=%f glambda2=%f\n",slowc_save,g0, res2, glambda2); - mexPrintf("-----------------------------------\n"); - mexPrintf(" Simulate iteration no %d \n", iter+1); - mexPrintf(" max. error=%.10e \n", double (max_res)); - mexPrintf(" sqr. error=%.10e \n", double (res2)); - mexPrintf(" abs. error=%.10e \n", double (res1)); - mexPrintf("-----------------------------------\n"); - } - - if (solve_algo == 5) - Simple_Init(it_, y_kmin, y_kmax, Size, IM_i); - else - { - b_m = mxCreateDoubleMatrix(Size,1,mxREAL); - if (!b_m) - { - mexPrintf("Can't allocate b_m matrix in LU solver\n"); - mexErrMsgTxt("end of bytecode\n"); - } - A_m = mxCreateSparse(Size, Size, min(int(IM_i.size()*2), Size*Size), mxREAL); - if (!A_m) - { - mexPrintf("Can't allocate A_m matrix in LU solver\n"); - mexErrMsgTxt("end of bytecode\n"); - } - Init_Matlab_Sparse_Simple(Size, IM_i, A_m, b_m); - } - if (solve_algo == 5) - { - NonZeroElem **bc; - bc = (NonZeroElem **) mxMalloc(Size*sizeof(*bc)); - for (i = 0; i < Size; i++) - { - /*finding the max-pivot*/ - double piv = piv_abs = 0; - int nb_eq = At_Col(i, &first); - l = 0; N_max = 0; - one = false; - piv_abs = 0; - for (j = 0; j < nb_eq; j++) - { - if (!line_done[first->r_index]) - { - k = first->u_index; - int jj = first->r_index; - int NRow_jj = NRow(jj); - - piv_v[l] = u[k]; - double piv_fabs = fabs(u[k]); - pivj_v[l] = jj; - pivk_v[l] = k; - NR[l] = NRow_jj; - if (NRow_jj == 1 && !one) - { - one = true; - piv_abs = piv_fabs; - N_max = NRow_jj; - } - if (!one) - { - if (piv_fabs > piv_abs) - piv_abs = piv_fabs; - if (NRow_jj > N_max) - N_max = NRow_jj; - } - else - { - if (NRow_jj == 1) - { - if (piv_fabs > piv_abs) - piv_abs = piv_fabs; - if (NRow_jj > N_max) - N_max = NRow_jj; - } - } - l++; - } - first = first->NZE_C_N; - } - if (piv_abs < eps) - { - if (Block_number > 1) - mexPrintf("Error: singular system in Simulate_NG in block %d\n", blck+1); - else - mexPrintf("Error: singular system in Simulate_NG\n"); - mexEvalString("drawnow;"); - mxFree(piv_v); - mxFree(pivj_v); - mxFree(pivk_v); - mxFree(NR); - mxFree(bc); - if (steady_state) - return false; - else - { - mexEvalString("st=fclose('all');clear all;"); - filename += " stopped"; - mexErrMsgTxt(filename.c_str()); - } - } - double markovitz = 0, markovitz_max = -9e70; - int NR_max = 0; - if (!one) - { - for (j = 0; j < l; j++) - { - if (N_max > 0 && NR[j] > 0) - { - if (fabs(piv_v[j]) > 0) - { - if (markowitz_c > 0) - markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double (NR[j])/double (N_max))); - else - markovitz = fabs(piv_v[j])/piv_abs; - } - else - markovitz = 0; - } - else - markovitz = fabs(piv_v[j])/piv_abs; - if (markovitz > markovitz_max) - { - piv = piv_v[j]; - pivj = pivj_v[j]; //Line number - pivk = pivk_v[j]; //positi - markovitz_max = markovitz; - NR_max = NR[j]; - } - } - } - else - { - for (j = 0; j < l; j++) - { - if (N_max > 0 && NR[j] > 0) - { - if (fabs(piv_v[j]) > 0) - { - if (markowitz_c > 0) - markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double (NR[j])/double (N_max))); - else - markovitz = fabs(piv_v[j])/piv_abs; - } - else - markovitz = 0; - } - else - markovitz = fabs(piv_v[j])/piv_abs; - if (NR[j] == 1) - { - piv = piv_v[j]; - pivj = pivj_v[j]; //Line number - pivk = pivk_v[j]; //positi - markovitz_max = markovitz; - NR_max = NR[j]; - } - } - } - pivot[i] = pivj; - pivotk[i] = pivk; - pivotv[i] = piv; - line_done[pivj] = true; - - /*divide all the non zeros elements of the line pivj by the max_pivot*/ - int nb_var = At_Row(pivj, &first); - for (j = 0; j < nb_var; j++) - { - u[first->u_index] /= piv; - first = first->NZE_R_N; - } - u[b[pivj]] /= piv; - /*substract the elements on the non treated lines*/ - nb_eq = At_Col(i, &first); - NonZeroElem *first_piva; - int nb_var_piva = At_Row(pivj, &first_piva); - int nb_eq_todo = 0; - for (j = 0; j < nb_eq && first; j++) - { - if (!line_done[first->r_index]) - bc[nb_eq_todo++] = first; - first = first->NZE_C_N; - } - //#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS"))) - for (j = 0; j < nb_eq_todo; j++) - { - first = bc[j]; - int row = first->r_index; - double first_elem = u[first->u_index]; - - int nb_var_piv = nb_var_piva; - NonZeroElem *first_piv = first_piva; - NonZeroElem *first_sub; - int nb_var_sub = At_Row(row, &first_sub); - int l_sub = 0, l_piv = 0; - int sub_c_index = first_sub->c_index, piv_c_index = first_piv->c_index; - while (l_sub < nb_var_sub || l_piv < nb_var_piv) - { - if (l_sub < nb_var_sub && (sub_c_index < piv_c_index || l_piv >= nb_var_piv)) - { - first_sub = first_sub->NZE_R_N; - if (first_sub) - sub_c_index = first_sub->c_index; - else - sub_c_index = Size; - l_sub++; - } - else if (sub_c_index > piv_c_index || l_sub >= nb_var_sub) - { - int tmp_u_count = Get_u(); - Insert(row, first_piv->c_index, tmp_u_count, 0); - u[tmp_u_count] = -u[first_piv->u_index]*first_elem; - first_piv = first_piv->NZE_R_N; - if (first_piv) - piv_c_index = first_piv->c_index; - else - piv_c_index = Size; - l_piv++; - } - else - { - if (i == sub_c_index) - { - firsta = first; - first_suba = first_sub->NZE_R_N; - Delete(first_sub->r_index, first_sub->c_index); - first = firsta->NZE_C_N; - first_sub = first_suba; - if (first_sub) - sub_c_index = first_sub->c_index; - else - sub_c_index = Size; - l_sub++; - first_piv = first_piv->NZE_R_N; - if (first_piv) - piv_c_index = first_piv->c_index; - else - piv_c_index = Size; - l_piv++; - } - else - { - u[first_sub->u_index] -= u[first_piv->u_index]*first_elem; - first_sub = first_sub->NZE_R_N; - if (first_sub) - sub_c_index = first_sub->c_index; - else - sub_c_index = Size; - l_sub++; - first_piv = first_piv->NZE_R_N; - if (first_piv) - piv_c_index = first_piv->c_index; - else - piv_c_index = Size; - l_piv++; - } - } - } - u[b[row]] -= u[b[pivj]]*first_elem; - } - } - double slowc_lbx = slowc, res1bx; - for (i = 0; i < y_size; i++) - ya[i+it_*y_size] = y[i+it_*y_size]; - slowc_save = slowc; - res1bx = simple_bksub(it_, Size, slowc_lbx); - End_GE(Size); - mxFree(piv_v); - mxFree(pivj_v); - mxFree(pivk_v); - mxFree(NR); - mxFree(bc); - } - else if (solve_algo == 1 || solve_algo == 4 || solve_algo == 0) - Solve_Matlab_LU_UMFPack(A_m, b_m, Size, slowc, false, it_); - else if (solve_algo == 2) - Solve_Matlab_GMRES(A_m, b_m, Size, slowc, blck, false, it_); - else if (solve_algo == 3) - Solve_Matlab_BiCGStab(A_m, b_m, Size, slowc, blck, false, it_); - return true; -} - void SparseMatrix::CheckIt(int y_size, int y_kmin, int y_kmax, int Size, int periods, int iter) { @@ -1620,9 +1222,9 @@ SparseMatrix::CheckIt(int y_size, int y_kmin, int y_kmax, int Size, int periods, SaveResult.open(out.str().c_str(), ios::in); if (!SaveResult.is_open()) { - mexPrintf("Error : Can't open file \"%s\" for reading\n", "Result"); - mexEvalString("st=fclose('all');clear all;"); - mexErrMsgTxt("Exit from Dynare"); + ostringstream tmp; + tmp << " in CheckIt, Result file cannot be opened\n"; + throw FatalExceptionHandling(tmp.str()); } mexPrintf("Reading Result..."); int row, col; @@ -1717,8 +1319,6 @@ SparseMatrix::Check_the_Solution(int periods, int y_kmin, int y_kmax, int Size, if (abs(res) > epsilon) mexPrintf("Error for equation %d => res=%f y[%d]=%f u[b[%d]]=%f somme(y*u)=%f\n", pos, res, pos, y[index_vara[pos]], pos, u[b[pos]], tmp_); } - filename += " stopped"; - mexErrMsgTxt(filename.c_str()); } void @@ -1729,26 +1329,30 @@ SparseMatrix::Solve_Matlab_Relaxation(mxArray* A_m, mxArray* b_m, unsigned int S double* b_m_d = mxGetPr(b_m); if (!b_m_d) { - mexPrintf("Can't retrieve b_m vector in LU solver\n"); - mexErrMsgTxt("end of bytecode\n"); + ostringstream tmp; + tmp << " in Solve_Matlab_Relaxation, can't retrieve b_m vector\n"; + throw FatalExceptionHandling(tmp.str()); } mwIndex *A_m_i = mxGetIr(A_m); if (!A_m_i) { - mexPrintf("Can't allocate A_m_i index vector in LU solver\n"); - mexErrMsgTxt("end of bytecode\n"); + ostringstream tmp; + tmp << " in Solve_Matlab_Relaxation, can't allocate A_m_i index vector\n"; + throw FatalExceptionHandling(tmp.str()); } mwIndex *A_m_j = mxGetJc(A_m); if (!A_m_j) { - mexPrintf("Can't allocate A_m_j index vector in LU solver\n"); - mexErrMsgTxt("end of bytecode\n"); + ostringstream tmp; + tmp << " in Solve_Matlab_Relaxation, can't allocate A_m_j index vector\n"; + throw FatalExceptionHandling(tmp.str()); } double *A_m_d = mxGetPr(A_m); if (!A_m_d) { - mexPrintf("Can't retrieve A matrix in LU solver\n"); - mexErrMsgTxt("end of bytecode\n"); + ostringstream tmp; + tmp << " in Solve_Matlab_Relaxation, can't retrieve A matrix\n"; + throw FatalExceptionHandling(tmp.str()); } unsigned int max_nze = A_m_j[Size*periods]; unsigned int nze = 0; @@ -2198,29 +1802,891 @@ SparseMatrix::Solve_Matlab_BiCGStab(mxArray* A_m, mxArray* b_m, int Size, double mxDestroyArray(flag); } +void +SparseMatrix::Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, bool steady_state, int it_) +{ + bool one; + int pivj = 0, pivk = 0; + double *piv_v; + int *pivj_v, *pivk_v, *NR; + int l, N_max; + NonZeroElem *first, *firsta, *first_suba; + double piv_abs; + NonZeroElem **bc; + bc = (NonZeroElem **) mxMalloc(Size*sizeof(*bc)); + piv_v = (double *) mxMalloc(Size*sizeof(double)); + pivj_v = (int *) mxMalloc(Size*sizeof(int)); + pivk_v = (int *) mxMalloc(Size*sizeof(int)); + NR = (int *) mxMalloc(Size*sizeof(int)); -int -SparseMatrix::simulate_NG1(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 Block_number, int stack_solve_algo, unsigned int endo_name_length, char *P_endo_names) + for (int i = 0; i < Size; i++) + { + /*finding the max-pivot*/ + double piv = piv_abs = 0; + int nb_eq = At_Col(i, &first); + l = 0; N_max = 0; + one = false; + piv_abs = 0; + for (int j = 0; j < nb_eq; j++) + { + if (!line_done[first->r_index]) + { + int k = first->u_index; + int jj = first->r_index; + int NRow_jj = NRow(jj); + + piv_v[l] = u[k]; + double piv_fabs = fabs(u[k]); + pivj_v[l] = jj; + pivk_v[l] = k; + NR[l] = NRow_jj; + if (NRow_jj == 1 && !one) + { + one = true; + piv_abs = piv_fabs; + N_max = NRow_jj; + } + if (!one) + { + if (piv_fabs > piv_abs) + piv_abs = piv_fabs; + if (NRow_jj > N_max) + N_max = NRow_jj; + } + else + { + if (NRow_jj == 1) + { + if (piv_fabs > piv_abs) + piv_abs = piv_fabs; + if (NRow_jj > N_max) + N_max = NRow_jj; + } + } + l++; + } + first = first->NZE_C_N; + } + if (piv_abs < eps) + { + mexEvalString("drawnow;"); + mxFree(piv_v); + mxFree(pivj_v); + mxFree(pivk_v); + mxFree(NR); + mxFree(bc); + if (steady_state) + { + if (blck > 1) + mexPrintf("Error: singular system in Simulate_NG in block %d\n", blck+1); + else + mexPrintf("Error: singular system in Simulate_NG\n"); + return; + } + else + { + ostringstream tmp; + if (blck > 1) + tmp << " in Solve_ByteCode_Sparse_GaussianElimination, singular system in block " << blck+1 << "\n"; + else + tmp << " in Solve_ByteCode_Sparse_GaussianElimination, singular system\n"; + throw FatalExceptionHandling(tmp.str()); + } + } + double markovitz = 0, markovitz_max = -9e70; + int NR_max = 0; + if (!one) + { + for (int j = 0; j < l; j++) + { + if (N_max > 0 && NR[j] > 0) + { + if (fabs(piv_v[j]) > 0) + { + if (markowitz_c > 0) + markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double (NR[j])/double (N_max))); + else + markovitz = fabs(piv_v[j])/piv_abs; + } + else + markovitz = 0; + } + else + markovitz = fabs(piv_v[j])/piv_abs; + if (markovitz > markovitz_max) + { + piv = piv_v[j]; + pivj = pivj_v[j]; //Line number + pivk = pivk_v[j]; //positi + markovitz_max = markovitz; + NR_max = NR[j]; + } + } + } + else + { + for (int j = 0; j < l; j++) + { + if (N_max > 0 && NR[j] > 0) + { + if (fabs(piv_v[j]) > 0) + { + if (markowitz_c > 0) + markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double (NR[j])/double (N_max))); + else + markovitz = fabs(piv_v[j])/piv_abs; + } + else + markovitz = 0; + } + else + markovitz = fabs(piv_v[j])/piv_abs; + if (NR[j] == 1) + { + piv = piv_v[j]; + pivj = pivj_v[j]; //Line number + pivk = pivk_v[j]; //positi + markovitz_max = markovitz; + NR_max = NR[j]; + } + } + } + pivot[i] = pivj; + pivotk[i] = pivk; + pivotv[i] = piv; + line_done[pivj] = true; + + /*divide all the non zeros elements of the line pivj by the max_pivot*/ + int nb_var = At_Row(pivj, &first); + for (int j = 0; j < nb_var; j++) + { + u[first->u_index] /= piv; + first = first->NZE_R_N; + } + u[b[pivj]] /= piv; + /*substract the elements on the non treated lines*/ + nb_eq = At_Col(i, &first); + NonZeroElem *first_piva; + int nb_var_piva = At_Row(pivj, &first_piva); + int nb_eq_todo = 0; + for (int j = 0; j < nb_eq && first; j++) + { + if (!line_done[first->r_index]) + bc[nb_eq_todo++] = first; + first = first->NZE_C_N; + } + //#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS"))) + for (int j = 0; j < nb_eq_todo; j++) + { + first = bc[j]; + int row = first->r_index; + double first_elem = u[first->u_index]; + + int nb_var_piv = nb_var_piva; + NonZeroElem *first_piv = first_piva; + NonZeroElem *first_sub; + int nb_var_sub = At_Row(row, &first_sub); + int l_sub = 0, l_piv = 0; + int sub_c_index = first_sub->c_index, piv_c_index = first_piv->c_index; + while (l_sub < nb_var_sub || l_piv < nb_var_piv) + { + if (l_sub < nb_var_sub && (sub_c_index < piv_c_index || l_piv >= nb_var_piv)) + { + first_sub = first_sub->NZE_R_N; + if (first_sub) + sub_c_index = first_sub->c_index; + else + sub_c_index = Size; + l_sub++; + } + else if (sub_c_index > piv_c_index || l_sub >= nb_var_sub) + { + int tmp_u_count = Get_u(); + Insert(row, first_piv->c_index, tmp_u_count, 0); + u[tmp_u_count] = -u[first_piv->u_index]*first_elem; + first_piv = first_piv->NZE_R_N; + if (first_piv) + piv_c_index = first_piv->c_index; + else + piv_c_index = Size; + l_piv++; + } + else + { + if (i == sub_c_index) + { + firsta = first; + first_suba = first_sub->NZE_R_N; + Delete(first_sub->r_index, first_sub->c_index); + first = firsta->NZE_C_N; + first_sub = first_suba; + if (first_sub) + sub_c_index = first_sub->c_index; + else + sub_c_index = Size; + l_sub++; + first_piv = first_piv->NZE_R_N; + if (first_piv) + piv_c_index = first_piv->c_index; + else + piv_c_index = Size; + l_piv++; + } + else + { + u[first_sub->u_index] -= u[first_piv->u_index]*first_elem; + first_sub = first_sub->NZE_R_N; + if (first_sub) + sub_c_index = first_sub->c_index; + else + sub_c_index = Size; + l_sub++; + first_piv = first_piv->NZE_R_N; + if (first_piv) + piv_c_index = first_piv->c_index; + else + piv_c_index = Size; + l_piv++; + } + } + } + u[b[row]] -= u[b[pivj]]*first_elem; + } + } + double slowc_lbx = slowc, res1bx; + for (int i = 0; i < y_size; i++) + ya[i+it_*y_size] = y[i+it_*y_size]; + slowc_save = slowc; + res1bx = simple_bksub(it_, Size, slowc_lbx); + End_GE(Size); + mxFree(piv_v); + mxFree(pivj_v); + mxFree(pivk_v); + mxFree(NR); + mxFree(bc); +} + +void +SparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bool symbolic, int Block_number) { /*Triangularisation at each period of a block using a simple gaussian Elimination*/ t_save_op_s *save_op_s; - bool record = false; int *save_op = NULL, *save_opa = NULL, *save_opaa = NULL; long int nop = 0, nopa = 0; - int tbreak = 0, last_period = periods; - int i, j, k; - int pivj = 0, pivk = 0; - int tmp_u_count, lag; - NonZeroElem *first; + bool record = false; + double *piv_v; double piv_abs; + int *pivj_v, *pivk_v, *NR; + int pivj = 0, pivk = 0; + NonZeroElem *first; + int tmp_u_count, lag; + int tbreak = 0, last_period = periods; + + piv_v = (double *) mxMalloc(Size*sizeof(double)); + pivj_v = (int *) mxMalloc(Size*sizeof(int)); + pivk_v = (int *) mxMalloc(Size*sizeof(int)); + NR = (int *) mxMalloc(Size*sizeof(int)); + + for (int t = 0; t < periods; t++) + { + if (record && symbolic) + { + if (save_op) + { + mxFree(save_op); + save_op = NULL; + } + save_op = (int *) mxMalloc(nop*sizeof(int)); + nopa = nop; + } + nop = 0; + Clear_u(); + int ti = t*Size; + for (int i = ti; i < Size+ti; i++) + { + /*finding the max-pivot*/ + double piv = piv_abs = 0; + int nb_eq = At_Col(i, 0, &first); + if ((symbolic && t <= start_compare) || !symbolic) + { + int l = 0, N_max = 0; + bool one = false; + piv_abs = 0; + for (int j = 0; j < nb_eq; j++) + { + if (!line_done[first->r_index]) + { + int k = first->u_index; + int jj = first->r_index; + int NRow_jj = NRow(jj); + piv_v[l] = u[k]; + double piv_fabs = fabs(u[k]); + pivj_v[l] = jj; + pivk_v[l] = k; + NR[l] = NRow_jj; + if (NRow_jj == 1 && !one) + { + one = true; + piv_abs = piv_fabs; + N_max = NRow_jj; + } + if (!one) + { + if (piv_fabs > piv_abs) + piv_abs = piv_fabs; + if (NRow_jj > N_max) + N_max = NRow_jj; + } + else + { + if (NRow_jj == 1) + { + if (piv_fabs > piv_abs) + piv_abs = piv_fabs; + if (NRow_jj > N_max) + N_max = NRow_jj; + } + } + l++; + } + first = first->NZE_C_N; + } + double markovitz = 0, markovitz_max = -9e70; + int NR_max = 0; + if (!one) + { + for (int j = 0; j < l; j++) + { + if (N_max > 0 && NR[j] > 0) + { + if (fabs(piv_v[j]) > 0) + { + if (markowitz_c > 0) + markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double (NR[j])/double (N_max))); + else + markovitz = fabs(piv_v[j])/piv_abs; + } + else + markovitz = 0; + } + else + markovitz = fabs(piv_v[j])/piv_abs; + if (markovitz > markovitz_max) + { + piv = piv_v[j]; + pivj = pivj_v[j]; //Line number + pivk = pivk_v[j]; //positi + markovitz_max = markovitz; + NR_max = NR[j]; + } + } + } + else + { + for (int j = 0; j < l; j++) + { + if (N_max > 0 && NR[j] > 0) + { + if (fabs(piv_v[j]) > 0) + { + if (markowitz_c > 0) + markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double (NR[j])/double (N_max))); + else + markovitz = fabs(piv_v[j])/piv_abs; + } + else + markovitz = 0; + } + else + markovitz = fabs(piv_v[j])/piv_abs; + if (NR[j] == 1) + { + piv = piv_v[j]; + pivj = pivj_v[j]; //Line number + pivk = pivk_v[j]; //positi + markovitz_max = markovitz; + NR_max = NR[j]; + } + } + } + if (fabs(piv) < eps) + mexPrintf("==> Error NR_max=%d, N_max=%d and piv=%f, piv_abs=%f, markovitz_max=%f\n",NR_max, N_max, piv, piv_abs, markovitz_max); + if (NR_max == 0) + mexPrintf("==> Error NR_max=0 and piv=%f, markovitz_max=%f\n",piv, markovitz_max); + pivot[i] = pivj; + pivot_save[i] = pivj; + pivotk[i] = pivk; + pivotv[i] = piv; + } + else + { + pivj = pivot[i-Size]+Size; + pivot[i] = pivj; + At_Pos(pivj, i, &first); + pivk = first->u_index; + piv = u[pivk]; + piv_abs = fabs(piv); + } + line_done[pivj] = true; + if (symbolic) + { + if (record) + { + if (nop+1 >= nopa) + { + nopa = long (mem_increasing_factor*(double)nopa); + save_op = (int *) mxRealloc(save_op, nopa*sizeof(int)); + } + save_op_s = (t_save_op_s *)(&(save_op[nop])); + save_op_s->operat = IFLD; + save_op_s->first = pivk; + save_op_s->lag = 0; + } + nop += 2; + } + if (piv_abs < eps) + { + ostringstream tmp; + if (Block_number>1) + tmp << " in Solve_ByteCode_Symbolic_Sparse_GaussianElimination, singular system in block " << Block_number+1 << "\n"; + else + tmp << " in Solve_ByteCode_Symbolic_Sparse_GaussianElimination, singular system\n"; + throw FatalExceptionHandling(tmp.str()); + } + /*divide all the non zeros elements of the line pivj by the max_pivot*/ + int nb_var = At_Row(pivj, &first); + NonZeroElem **bb; + bb = (NonZeroElem **) mxMalloc(nb_var*sizeof(first)); + for (int j = 0; j < nb_var; j++) + { + bb[j] = first; + first = first->NZE_R_N; + } + + for (int j = 0; j < nb_var; j++) + { + first = bb[j]; + u[first->u_index] /= piv; + if (symbolic) + { + if (record) + { + if (nop+j*2+1 >= nopa) + { + nopa = long (mem_increasing_factor*(double)nopa); + save_op = (int *) mxRealloc(save_op, nopa*sizeof(int)); + } + save_op_s = (t_save_op_s *)(&(save_op[nop+j*2])); + save_op_s->operat = IFDIV; + save_op_s->first = first->u_index; + save_op_s->lag = first->lag_index; + } + } + } + mxFree(bb); + nop += nb_var*2; + u[b[pivj]] /= piv; + if (symbolic) + { + if (record) + { + if (nop+1 >= nopa) + { + nopa = long (mem_increasing_factor*(double)nopa); + save_op = (int *) mxRealloc(save_op, nopa*sizeof(int)); + } + save_op_s = (t_save_op_s *)(&(save_op[nop])); + save_op_s->operat = IFDIV; + save_op_s->first = b[pivj]; + save_op_s->lag = 0; + } + nop += 2; + } + /*substract the elements on the non treated lines*/ + nb_eq = At_Col(i, &first); + NonZeroElem *first_piva; + int nb_var_piva = At_Row(pivj, &first_piva); + + NonZeroElem **bc; + bc = (NonZeroElem **) mxMalloc(nb_eq*sizeof(first)); + int nb_eq_todo = 0; + for (int j = 0; j < nb_eq && first; j++) + { + if (!line_done[first->r_index]) + bc[nb_eq_todo++] = first; + first = first->NZE_C_N; + } + //#pragma omp parallel for num_threads(2) shared(nb_var_piva, first_piva, nopa, nop, save_op, record) + for (int j = 0; j < nb_eq_todo; j++) + { + t_save_op_s *save_op_s_l; + first = bc[j]; + int row = first->r_index; + double first_elem = u[first->u_index]; + if (symbolic) + { + if (record) + { + if (nop+1 >= nopa) + { + nopa = long (mem_increasing_factor*(double)nopa); + save_op = (int *) mxRealloc(save_op, nopa*sizeof(int)); + } + save_op_s_l = (t_save_op_s *)(&(save_op[nop])); + save_op_s_l->operat = IFLD; + save_op_s_l->first = first->u_index; + save_op_s_l->lag = abs(first->lag_index); + } + nop += 2; + } + + int nb_var_piv = nb_var_piva; + NonZeroElem *first_piv = first_piva; + NonZeroElem *first_sub; + int nb_var_sub = At_Row(row, &first_sub); + int l_sub = 0; + int l_piv = 0; + int sub_c_index = first_sub->c_index; + int piv_c_index = first_piv->c_index; + int tmp_lag = first_sub->lag_index; + while (l_sub < nb_var_sub || l_piv < nb_var_piv) + { + if (l_sub < nb_var_sub && (sub_c_index < piv_c_index || l_piv >= nb_var_piv)) + { + //There is no nonzero element at row pivot for this column=> Nothing to do for the current element got to next column + first_sub = first_sub->NZE_R_N; + if (first_sub) + sub_c_index = first_sub->c_index; + else + sub_c_index = Size*periods; + l_sub++; + } + else if (sub_c_index > piv_c_index || l_sub >= nb_var_sub) + { + // There is an nonzero element at row pivot but not at the current row=> insert a negative element in the current row + tmp_u_count = Get_u(); + lag = first_piv->c_index/Size-row/Size; + //#pragma omp critical + { + Insert(row, first_piv->c_index, tmp_u_count, lag); + } + u[tmp_u_count] = -u[first_piv->u_index]*first_elem; + if (symbolic) + { + if (record) + { + if (nop+2 >= nopa) + { + nopa = long (mem_increasing_factor*(double)nopa); + save_op = (int *) mxRealloc(save_op, nopa*sizeof(int)); + } + save_op_s_l = (t_save_op_s *)(&(save_op[nop])); + save_op_s_l->operat = IFLESS; + save_op_s_l->first = tmp_u_count; + save_op_s_l->second = first_piv->u_index; + save_op_s_l->lag = max(first_piv->lag_index, abs(tmp_lag)); + } + nop += 3; + } + first_piv = first_piv->NZE_R_N; + if (first_piv) + piv_c_index = first_piv->c_index; + else + piv_c_index = Size*periods; + l_piv++; + } + else /*first_sub->c_index==first_piv->c_index*/ + { + if (i == sub_c_index) + { + NonZeroElem *firsta = first; + NonZeroElem *first_suba = first_sub->NZE_R_N; + Delete(first_sub->r_index, first_sub->c_index); + first = firsta->NZE_C_N; + first_sub = first_suba; + if (first_sub) + sub_c_index = first_sub->c_index; + else + sub_c_index = Size*periods; + l_sub++; + first_piv = first_piv->NZE_R_N; + if (first_piv) + piv_c_index = first_piv->c_index; + else + piv_c_index = Size*periods; + l_piv++; + } + else + { + u[first_sub->u_index] -= u[first_piv->u_index]*first_elem; + if (symbolic) + { + if (record) + { + if (nop+3 >= nopa) + { + nopa = long (mem_increasing_factor*(double)nopa); + save_op = (int *) mxRealloc(save_op, nopa*sizeof(int)); + } + save_op_s_l = (t_save_op_s *)(&(save_op[nop])); + save_op_s_l->operat = IFSUB; + save_op_s_l->first = first_sub->u_index; + save_op_s_l->second = first_piv->u_index; + save_op_s_l->lag = max(abs(tmp_lag), first_piv->lag_index); + } + nop += 3; + } + first_sub = first_sub->NZE_R_N; + if (first_sub) + sub_c_index = first_sub->c_index; + else + sub_c_index = Size*periods; + l_sub++; + first_piv = first_piv->NZE_R_N; + if (first_piv) + piv_c_index = first_piv->c_index; + else + piv_c_index = Size*periods; + l_piv++; + } + } + } + u[b[row]] -= u[b[pivj]]*first_elem; + + if (symbolic) + { + if (record) + { + if (nop+3 >= nopa) + { + nopa = long (mem_increasing_factor*(double)nopa); + save_op = (int *) mxRealloc(save_op, nopa*sizeof(int)); + } + save_op_s_l = (t_save_op_s *)(&(save_op[nop])); + save_op_s_l->operat = IFSUB; + save_op_s_l->first = b[row]; + save_op_s_l->second = b[pivj]; + save_op_s_l->lag = abs(tmp_lag); + } + nop += 3; + } + } + mxFree(bc); + } + if (symbolic) + { + if (record && (nop == nop1)) + { + if (save_opa && save_opaa) + { + if (compare(save_op, save_opa, save_opaa, t, periods, nop, Size)) + { + tbreak = t; + tbreak_g = tbreak; + break; + } + } + if (save_opa) + { + if (save_opaa) + { + mxFree(save_opaa); + save_opaa = NULL; + } + save_opaa = (int *) mxMalloc(nop1*sizeof(int)); + memcpy(save_opaa, save_opa, nop1*sizeof(int)); + } + if (save_opa) + { + mxFree(save_opa); + save_opa = NULL; + } + save_opa = (int *) mxMalloc(nop*sizeof(int)); + memcpy(save_opa, save_op, nop*sizeof(int)); + } + else + { + if (nop == nop1) + record = true; + else + { + record = false; + if (save_opa) + { + mxFree(save_opa); + save_opa = NULL; + } + if (save_opaa) + { + mxFree(save_opaa); + save_opaa = NULL; + } + } + } + nop2 = nop1; + nop1 = nop; + } + } + mxFree(piv_v); + mxFree(pivj_v); + mxFree(pivk_v); + mxFree(NR); + nop_all += nop; + if (symbolic) + { + if (save_op) + mxFree(save_op); + if (save_opa) + mxFree(save_opa); + if (save_opaa) + mxFree(save_opaa); + } + + /*The backward substitution*/ + double slowc_lbx = slowc, res1bx; + for (int i = 0; i < y_size*(periods+y_kmin); i++) + ya[i] = y[i]; + slowc_save = slowc; + res1bx = bksub(tbreak, last_period, Size, slowc_lbx); + End_GE(Size); +} + + +void +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 solve_algo) +{ + int i, j; + mxArray *b_m = NULL, *A_m = NULL; + Clear_u(); + error_not_printed = true; + u_count_alloc_save = u_count_alloc; + if (isnan(res1) || isinf(res1) || (res2 > 12*g0 && iter>0)) + { + if (iter == 0 || fabs(slowc_save) < 1e-8) + { + for (j = 0; j < y_size; j++) + { + bool select = false; + for (int i = 0; i < Size; i++) + if (j == index_vara[i]) + { + select = true; + break; + } + if (select) + mexPrintf("-> variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]); + else + mexPrintf(" variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]); + } + if (steady_state) + { + if (iter == 0) + mexPrintf(" the initial values of endogenous variables are too far from the solution.\nChange them!\n"); + 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; + } + else + { + ostringstream tmp; + if (iter == 0) + tmp << " in Simulate_Newton_One_Boundary, The initial values of endogenous variables are too far from the solution.\nChange them!\n"; + else + tmp << " in Simulate_Newton_One_Boundary, Dynare cannot improve the simulation in block " << blck+1 << " at time " << it_+1 << " (variable " << index_vara[max_res_idx]+1 << "%d)\n"; + throw FatalExceptionHandling(tmp.str()); + } + } + if(!(isnan(res1) || isinf(res1)) && !(isnan(g0) || isinf(g0))) + { + if (try_at_iteration == 0) + { + prev_slowc_save = slowc_save; + slowc_save = max( - gp0 / (2 * (res2 - g0 - gp0)) , 0.1); + } + else + { + double t1 = res2 - gp0 * slowc_save - g0; + double t2 = glambda2 - gp0 * prev_slowc_save - g0; + double a = (1/(slowc_save * slowc_save) * t1 - 1/(prev_slowc_save * prev_slowc_save) * t2) / (slowc_save - prev_slowc_save); + double b = (-prev_slowc_save/(slowc_save * slowc_save) * t1 + slowc_save/(prev_slowc_save * prev_slowc_save) * t2) / (slowc_save - prev_slowc_save); + prev_slowc_save = slowc_save; + slowc_save = max(min( -b + sqrt(b*b - 3 * a * gp0) / (3 * a), 0.5 * slowc_save), 0.1 * slowc_save); + } + glambda2 = res2; + try_at_iteration ++; + } + else + { + prev_slowc_save = slowc_save; + slowc_save /= 1.1; + } + if (print_it) + mexPrintf("Error: Simulation diverging, trying to correct it using slowc=%f\n", slowc_save); + 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; + } + if (cvg) + { + return; + } + if (print_it ) + { + mexPrintf("solwc=%f g0=%f res2=%f glambda2=%f\n",slowc_save,g0, res2, glambda2); + mexPrintf("-----------------------------------\n"); + mexPrintf(" Simulate iteration no %d \n", iter+1); + mexPrintf(" max. error=%.10e \n", double (max_res)); + mexPrintf(" sqr. error=%.10e \n", double (res2)); + mexPrintf(" abs. error=%.10e \n", double (res1)); + mexPrintf("-----------------------------------\n"); + } + + if (solve_algo == 5) + Simple_Init(it_, y_kmin, y_kmax, Size, IM_i); + else + { + b_m = mxCreateDoubleMatrix(Size,1,mxREAL); + if (!b_m) + { + ostringstream tmp; + tmp << " in Simulate_Newton_One_Boundary, can't allocate b_m vector\n"; + throw FatalExceptionHandling(tmp.str()); + } + A_m = mxCreateSparse(Size, Size, min(int(IM_i.size()*2), Size*Size), mxREAL); + if (!A_m) + { + ostringstream tmp; + tmp << " in Simulate_Newton_One_Boundary, can't allocate A_m matrix\n"; + throw FatalExceptionHandling(tmp.str()); + } + Init_Matlab_Sparse_Simple(Size, IM_i, A_m, b_m); + } + + if (solve_algo == 1 || solve_algo == 4 || solve_algo == 0) + Solve_Matlab_LU_UMFPack(A_m, b_m, Size, slowc, false, it_); + else if (solve_algo == 2) + Solve_Matlab_GMRES(A_m, b_m, Size, slowc, blck, false, it_); + else if (solve_algo == 3) + Solve_Matlab_BiCGStab(A_m, b_m, Size, slowc, blck, false, it_); + else if (solve_algo == 5) + Solve_ByteCode_Sparse_GaussianElimination(Size, blck, steady_state, it_); + return; +} + + + +void +SparseMatrix::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) +{ if (start_compare == 0) start_compare = y_kmin; u_count_alloc_save = u_count_alloc; - clock_t t01; clock_t t1 = clock(); nop1 = 0; error_not_printed = true; - mxArray *b_m, *A_m; + mxArray *b_m = NULL, *A_m = NULL; if (iter > 0) { mexPrintf("Sim : %f ms\n", (1000.0*(double (clock())-double (time00)))/double (CLOCKS_PER_SEC)); @@ -2229,9 +2695,9 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax } if (isnan(res1) || isinf(res1) || (res2 > 12*g0 && iter>0)) { - if (iter == 0) + if (iter == 0 || fabs(slowc_save) < 1e-8) { - for (j = 0; j < y_size; j++) + for (int j = 0; j < y_size; j++) { ostringstream res; for (unsigned int i = 0; i < endo_name_length; i++) @@ -2249,37 +2715,13 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax else mexPrintf(" variable %s (%d) at time %d = %f direction = %f\n", res.str().c_str(), j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]); } - mexPrintf("res1=%5.10f\n", res1); - mexPrintf("The initial values of endogenous variables are too far from the solution.\n"); - mexPrintf("Change them!\n"); - mexEvalString("drawnow;"); - filename += " stopped"; - mexErrMsgTxt(filename.c_str()); - } - if (fabs(slowc_save) < 1e-8) - { - for (j = 0; j < y_size; j++) - { - ostringstream res(""); - for (unsigned int i = 0; i < endo_name_length; i++) - if (P_endo_names[2*(j+i*y_size)] != ' ') - res << P_endo_names[2*(j+i*y_size)]; - bool select = false; - for (int i = 0; i < Size; i++) - if (j == index_vara[i]) - { - select = true; - break; - } - if (select) - mexPrintf("-> variable %s (%d) at time %d = %f direction = %f\n", res.str().c_str(), j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]); - else - mexPrintf(" variable %s (%d) at time %d = %f direction = %f\n", res.str().c_str(), j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]); - } - mexPrintf("Dynare cannot improve the simulation in block %d at time %d (variable %d)\n", blck+1, it_+1, max_res_idx); - mexEvalString("drawnow;"); - filename += " stopped"; - mexErrMsgTxt(filename.c_str()); + ostringstream Error; + if (iter == 0 ) + Error << " in Simulate_Newton_Two_Boundaries, the initial values of endogenous variables are too far from the solution.\nChange them!\n"; + else + Error << " in Simulate_Newton_Two_Boundaries, dynare cannot improve the simulation in block " << blck+1 << " at time " << it_+1 << " (variable " << index_vara[max_res_idx]+1 << ")\n"; + //Error << filename << " stopped"; + throw FatalExceptionHandling(Error.str()); } if(!(isnan(res1) || isinf(res1)) && !(isnan(g0) || isinf(g0)) && (stack_solve_algo == 4 || stack_solve_algo == 5)) { @@ -2301,13 +2743,13 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax try_at_iteration ++; if (slowc_save<=0.1) { - for (i = 0; i < y_size*(periods+y_kmin); i++) + for (int i = 0; i < y_size*(periods+y_kmin); i++) y[i] = ya[i]+direction[i]; g0 = res2; gp0 = -res2; try_at_iteration = 0; iter--; - return 0; + return; } } else @@ -2318,50 +2760,17 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax if (print_it) { if (isnan(res1) || isinf(res1)) - mexPrintf("Error: The model cannot be evaluated, trying to correct it using slowc=%f\n", slowc_save); + mexPrintf("The model cannot be evaluated, trying to correct it using slowc=%f\n", slowc_save); else - mexPrintf("Error: Simulation diverging, trying to correct it using slowc=%f\n", slowc_save); + mexPrintf("Simulation diverging, trying to correct it using slowc=%f\n", slowc_save); } - for (i = 0; i < y_size*(periods+y_kmin); i++) + for (int i = 0; i < y_size*(periods+y_kmin); i++) y[i] = ya[i]+slowc_save*direction[i]; iter--; - return 0; - } - if (isnan(res1) || isinf(res1)) - { - if (iter == 0) - { - for (j = 0; j < y_size; j++) - mexPrintf("variable %d at time %d = %f\n", j+1, it_, y[j+it_*y_size]); - for (j = 0; j < Size; j++) - mexPrintf("residual(%d)=%5.25f\n", j, u[j]); - mexPrintf("res1=%5.25f\n", res1); - mexPrintf("The initial values of endogenous variables are too far from the solution.\n"); - mexPrintf("Change them!\n"); - mexEvalString("drawnow;"); - mexEvalString("st=fclose('all');clear all;"); - filename += " stopped"; - mexErrMsgTxt(filename.c_str()); - } - if (slowc_save < 1e-8) - { - mexPrintf("slowc_save=%g\n", slowc_save); - for (j = 0; j < y_size; j++) - mexPrintf("variable %d at time %d = %f direction = %f variable at last step = %f b = %f\n", j+1, it_+1, y[j+it_*y_size], direction[j+it_*y_size], ya[j+it_*y_size], u[pivot[j+it_*y_size]]); - mexPrintf("Dynare cannot improve the simulation in block %d at time %d (variable %d max_res = %f, res1 = %f)\n", blck+1, it_+1, max_res_idx, max_res, res1); - mexEvalString("drawnow;"); - mexEvalString("st=fclose('all');clear all;"); - filename += " stopped"; - mexErrMsgTxt(filename.c_str()); - } - slowc_save /= 2; - mexPrintf("Error: Simulation diverging, trying to correct it using slowc=%f\n", slowc_save); - for (i = 0; i < y_size*(periods+y_kmin); i++) - y[i] = ya[i]+slowc_save*direction[i]; - iter--; - return (0); + return; } + u_count += u_count_init; if (stack_solve_algo == 5) { @@ -2436,7 +2845,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax } if (cvg) { - return (0); + return; } else { @@ -2447,498 +2856,30 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax b_m = mxCreateDoubleMatrix(periods*Size,1,mxREAL); if (!b_m) { - mexPrintf("Can't allocate b_m matrix in LU solver\n"); - mexErrMsgTxt("end of bytecode\n"); + ostringstream tmp; + tmp << " in Simulate_Newton_Two_Boundaries, can't allocate b_m vector\n"; + throw FatalExceptionHandling(tmp.str()); } A_m = mxCreateSparse(periods*Size, periods*Size, IM_i.size()* periods*2, mxREAL); if (!A_m) { - mexPrintf("Can't allocate A_m matrix in LU solver\n"); - mexErrMsgTxt("end of bytecode\n"); + ostringstream tmp; + tmp << " in Simulate_Newton_Two_Boundaries, can't allocate A_m matrix\n"; + throw FatalExceptionHandling(tmp.str()); } Init_Matlab_Sparse(periods, y_kmin, y_kmax, Size, IM_i, A_m, b_m); } - if (stack_solve_algo == 5) - { - double *piv_v; - int *pivj_v, *pivk_v, *NR; - piv_v = (double *) mxMalloc(Size*sizeof(double)); - pivj_v = (int *) mxMalloc(Size*sizeof(int)); - pivk_v = (int *) mxMalloc(Size*sizeof(int)); - NR = (int *) mxMalloc(Size*sizeof(int)); - for (int t = 0; t < periods; t++) - { - if (record && symbolic) - { - if (save_op) - { - mxFree(save_op); - save_op = NULL; - } - save_op = (int *) mxMalloc(nop*sizeof(int)); - nopa = nop; - } - nop = 0; - Clear_u(); - int ti = t*Size; - for (i = ti; i < Size+ti; i++) - { - /*finding the max-pivot*/ - double piv = piv_abs = 0; - int nb_eq = At_Col(i, 0, &first); - if ((symbolic && t <= start_compare) || !symbolic) - { - int l = 0, N_max = 0; - bool one = false; - piv_abs = 0; - for (j = 0; j < nb_eq; j++) - { - if (!line_done[first->r_index]) - { - k = first->u_index; - int jj = first->r_index; - int NRow_jj = NRow(jj); - piv_v[l] = u[k]; - double piv_fabs = fabs(u[k]); - pivj_v[l] = jj; - pivk_v[l] = k; - NR[l] = NRow_jj; - if (NRow_jj == 1 && !one) - { - one = true; - piv_abs = piv_fabs; - N_max = NRow_jj; - } - if (!one) - { - if (piv_fabs > piv_abs) - piv_abs = piv_fabs; - if (NRow_jj > N_max) - N_max = NRow_jj; - } - else - { - if (NRow_jj == 1) - { - if (piv_fabs > piv_abs) - piv_abs = piv_fabs; - if (NRow_jj > N_max) - N_max = NRow_jj; - } - } - l++; - } - first = first->NZE_C_N; - } - double markovitz = 0, markovitz_max = -9e70; - int NR_max = 0; - if (!one) - { - for (j = 0; j < l; j++) - { - if (N_max > 0 && NR[j] > 0) - { - if (fabs(piv_v[j]) > 0) - { - if (markowitz_c > 0) - markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double (NR[j])/double (N_max))); - else - markovitz = fabs(piv_v[j])/piv_abs; - } - else - markovitz = 0; - } - else - markovitz = fabs(piv_v[j])/piv_abs; - if (markovitz > markovitz_max) - { - piv = piv_v[j]; - pivj = pivj_v[j]; //Line number - pivk = pivk_v[j]; //positi - markovitz_max = markovitz; - NR_max = NR[j]; - } - } - } - else - { - for (j = 0; j < l; j++) - { - if (N_max > 0 && NR[j] > 0) - { - if (fabs(piv_v[j]) > 0) - { - if (markowitz_c > 0) - markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double (NR[j])/double (N_max))); - else - markovitz = fabs(piv_v[j])/piv_abs; - } - else - markovitz = 0; - } - else - markovitz = fabs(piv_v[j])/piv_abs; - if (NR[j] == 1) - { - piv = piv_v[j]; - pivj = pivj_v[j]; //Line number - pivk = pivk_v[j]; //positi - markovitz_max = markovitz; - NR_max = NR[j]; - } - } - } - if (fabs(piv) < eps) - mexPrintf("==> Error NR_max=%d, N_max=%d and piv=%f, piv_abs=%f, markovitz_max=%f\n",NR_max, N_max, piv, piv_abs, markovitz_max); - if (NR_max == 0) - mexPrintf("==> Error NR_max=0 and piv=%f, markovitz_max=%f\n",piv, markovitz_max); - pivot[i] = pivj; - pivot_save[i] = pivj; - pivotk[i] = pivk; - pivotv[i] = piv; - } - else - { - pivj = pivot[i-Size]+Size; - pivot[i] = pivj; - At_Pos(pivj, i, &first); - pivk = first->u_index; - piv = u[pivk]; - piv_abs = fabs(piv); - } - line_done[pivj] = true; - if (symbolic) - { - if (record) - { - if (nop+1 >= nopa) - { - nopa = long (mem_increasing_factor*(double)nopa); - save_op = (int *) mxRealloc(save_op, nopa*sizeof(int)); - } - save_op_s = (t_save_op_s *)(&(save_op[nop])); - save_op_s->operat = IFLD; - save_op_s->first = pivk; - save_op_s->lag = 0; - } - nop += 2; - } - if (piv_abs < eps) - { - if (Block_number>1) - mexPrintf("Error: singular system in Simulate_NG1 in block %d\n", blck); - else - mexPrintf("Error: singular system in Simulate_NG1\n"); - mexEvalString("drawnow;"); - filename += " stopped"; - mexErrMsgTxt(filename.c_str()); - } - /*divide all the non zeros elements of the line pivj by the max_pivot*/ - int nb_var = At_Row(pivj, &first); - NonZeroElem **bb; - bb = (NonZeroElem **) mxMalloc(nb_var*sizeof(first)); - for (j = 0; j < nb_var; j++) - { - bb[j] = first; - first = first->NZE_R_N; - } - for (j = 0; j < nb_var; j++) - { - first = bb[j]; - u[first->u_index] /= piv; - if (symbolic) - { - if (record) - { - if (nop+j*2+1 >= nopa) - { - nopa = long (mem_increasing_factor*(double)nopa); - save_op = (int *) mxRealloc(save_op, nopa*sizeof(int)); - } - save_op_s = (t_save_op_s *)(&(save_op[nop+j*2])); - save_op_s->operat = IFDIV; - save_op_s->first = first->u_index; - save_op_s->lag = first->lag_index; - } - } - } - mxFree(bb); - nop += nb_var*2; - u[b[pivj]] /= piv; - if (symbolic) - { - if (record) - { - if (nop+1 >= nopa) - { - nopa = long (mem_increasing_factor*(double)nopa); - save_op = (int *) mxRealloc(save_op, nopa*sizeof(int)); - } - save_op_s = (t_save_op_s *)(&(save_op[nop])); - save_op_s->operat = IFDIV; - save_op_s->first = b[pivj]; - save_op_s->lag = 0; - } - nop += 2; - } - /*substract the elements on the non treated lines*/ - nb_eq = At_Col(i, &first); - NonZeroElem *first_piva; - int nb_var_piva = At_Row(pivj, &first_piva); - - NonZeroElem **bc; - bc = (NonZeroElem **) mxMalloc(nb_eq*sizeof(first)); - int nb_eq_todo = 0; - for (j = 0; j < nb_eq && first; j++) - { - if (!line_done[first->r_index]) - bc[nb_eq_todo++] = first; - first = first->NZE_C_N; - } - //#pragma omp parallel for num_threads(2) shared(nb_var_piva, first_piva, nopa, nop, save_op, record) - for (j = 0; j < nb_eq_todo; j++) - { - t_save_op_s *save_op_s_l; - first = bc[j]; - int row = first->r_index; - double first_elem = u[first->u_index]; - if (symbolic) - { - if (record) - { - if (nop+1 >= nopa) - { - nopa = long (mem_increasing_factor*(double)nopa); - save_op = (int *) mxRealloc(save_op, nopa*sizeof(int)); - } - save_op_s_l = (t_save_op_s *)(&(save_op[nop])); - save_op_s_l->operat = IFLD; - save_op_s_l->first = first->u_index; - save_op_s_l->lag = abs(first->lag_index); - } - nop += 2; - } - - int nb_var_piv = nb_var_piva; - NonZeroElem *first_piv = first_piva; - NonZeroElem *first_sub; - int nb_var_sub = At_Row(row, &first_sub); - int l_sub = 0; - int l_piv = 0; - int sub_c_index = first_sub->c_index; - int piv_c_index = first_piv->c_index; - int tmp_lag = first_sub->lag_index; - while (l_sub < nb_var_sub || l_piv < nb_var_piv) - { - if (l_sub < nb_var_sub && (sub_c_index < piv_c_index || l_piv >= nb_var_piv)) - { - //There is no nonzero element at row pivot for this column=> Nothing to do for the current element got to next column - first_sub = first_sub->NZE_R_N; - if (first_sub) - sub_c_index = first_sub->c_index; - else - sub_c_index = Size*periods; - l_sub++; - } - else if (sub_c_index > piv_c_index || l_sub >= nb_var_sub) - { - // There is an nonzero element at row pivot but not at the current row=> insert a negative element in the current row - tmp_u_count = Get_u(); - lag = first_piv->c_index/Size-row/Size; - //#pragma omp critical - { - Insert(row, first_piv->c_index, tmp_u_count, lag); - } - u[tmp_u_count] = -u[first_piv->u_index]*first_elem; - if (symbolic) - { - if (record) - { - if (nop+2 >= nopa) - { - nopa = long (mem_increasing_factor*(double)nopa); - save_op = (int *) mxRealloc(save_op, nopa*sizeof(int)); - } - save_op_s_l = (t_save_op_s *)(&(save_op[nop])); - save_op_s_l->operat = IFLESS; - save_op_s_l->first = tmp_u_count; - save_op_s_l->second = first_piv->u_index; - save_op_s_l->lag = max(first_piv->lag_index, abs(tmp_lag)); - } - nop += 3; - } - first_piv = first_piv->NZE_R_N; - if (first_piv) - piv_c_index = first_piv->c_index; - else - piv_c_index = Size*periods; - l_piv++; - } - else /*first_sub->c_index==first_piv->c_index*/ - { - if (i == sub_c_index) - { - NonZeroElem *firsta = first; - NonZeroElem *first_suba = first_sub->NZE_R_N; - Delete(first_sub->r_index, first_sub->c_index); - first = firsta->NZE_C_N; - first_sub = first_suba; - if (first_sub) - sub_c_index = first_sub->c_index; - else - sub_c_index = Size*periods; - l_sub++; - first_piv = first_piv->NZE_R_N; - if (first_piv) - piv_c_index = first_piv->c_index; - else - piv_c_index = Size*periods; - l_piv++; - } - else - { - u[first_sub->u_index] -= u[first_piv->u_index]*first_elem; - if (symbolic) - { - if (record) - { - if (nop+3 >= nopa) - { - nopa = long (mem_increasing_factor*(double)nopa); - save_op = (int *) mxRealloc(save_op, nopa*sizeof(int)); - } - save_op_s_l = (t_save_op_s *)(&(save_op[nop])); - save_op_s_l->operat = IFSUB; - save_op_s_l->first = first_sub->u_index; - save_op_s_l->second = first_piv->u_index; - save_op_s_l->lag = max(abs(tmp_lag), first_piv->lag_index); - } - nop += 3; - } - first_sub = first_sub->NZE_R_N; - if (first_sub) - sub_c_index = first_sub->c_index; - else - sub_c_index = Size*periods; - l_sub++; - first_piv = first_piv->NZE_R_N; - if (first_piv) - piv_c_index = first_piv->c_index; - else - piv_c_index = Size*periods; - l_piv++; - } - } - } - u[b[row]] -= u[b[pivj]]*first_elem; - - if (symbolic) - { - if (record) - { - if (nop+3 >= nopa) - { - nopa = long (mem_increasing_factor*(double)nopa); - save_op = (int *) mxRealloc(save_op, nopa*sizeof(int)); - } - save_op_s_l = (t_save_op_s *)(&(save_op[nop])); - save_op_s_l->operat = IFSUB; - save_op_s_l->first = b[row]; - save_op_s_l->second = b[pivj]; - save_op_s_l->lag = abs(tmp_lag); - } - nop += 3; - } - } - mxFree(bc); - } - if (symbolic) - { - if (record && (nop == nop1)) - { - if (save_opa && save_opaa) - { - if (compare(save_op, save_opa, save_opaa, t, periods, nop, Size)) - { - tbreak = t; - tbreak_g = tbreak; - break; - } - } - if (save_opa) - { - if (save_opaa) - { - mxFree(save_opaa); - save_opaa = NULL; - } - save_opaa = (int *) mxMalloc(nop1*sizeof(int)); - memcpy(save_opaa, save_opa, nop1*sizeof(int)); - } - if (save_opa) - { - mxFree(save_opa); - save_opa = NULL; - } - save_opa = (int *) mxMalloc(nop*sizeof(int)); - memcpy(save_opa, save_op, nop*sizeof(int)); - } - else - { - if (nop == nop1) - record = true; - else - { - record = false; - if (save_opa) - { - mxFree(save_opa); - save_opa = NULL; - } - if (save_opaa) - { - mxFree(save_opaa); - save_opaa = NULL; - } - } - } - nop2 = nop1; - nop1 = nop; - } - } - mxFree(piv_v); - mxFree(pivj_v); - mxFree(pivk_v); - mxFree(NR); - nop_all += nop; - if (symbolic) - { - if (save_op) - mxFree(save_op); - if (save_opa) - mxFree(save_opa); - if (save_opaa) - mxFree(save_opaa); - } - - /*The backward substitution*/ - double slowc_lbx = slowc, res1bx; - for (i = 0; i < y_size*(periods+y_kmin); i++) - ya[i] = y[i]; - slowc_save = slowc; - res1bx = bksub(tbreak, last_period, Size, slowc_lbx); - t01 = clock(); - End_GE(Size); - } + if (stack_solve_algo == 0 || stack_solve_algo == 4) + Solve_Matlab_LU_UMFPack(A_m, b_m, Size, slowc, true, 0); else if (stack_solve_algo == 1) Solve_Matlab_Relaxation(A_m, b_m, Size, slowc, true, 0); - else if (stack_solve_algo == 0 || stack_solve_algo == 4) - Solve_Matlab_LU_UMFPack(A_m, b_m, Size, slowc, true, 0); else if (stack_solve_algo == 2) Solve_Matlab_GMRES(A_m, b_m, Size, slowc, blck, true, 0); else if (stack_solve_algo == 3) Solve_Matlab_BiCGStab(A_m, b_m, Size, slowc, blck, true, 0); + else if (stack_solve_algo == 5) + Solve_ByteCode_Symbolic_Sparse_GaussianElimination(Size, symbolic, blck); } if (print_it) { @@ -2950,8 +2891,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax time00 = clock(); if (tbreak_g == 0) tbreak_g = periods; - - return (0); + return; } void diff --git a/mex/sources/bytecode/SparseMatrix.hh b/mex/sources/bytecode/SparseMatrix.hh index b26e4df1f..287c35126 100644 --- a/mex/sources/bytecode/SparseMatrix.hh +++ b/mex/sources/bytecode/SparseMatrix.hh @@ -26,6 +26,7 @@ #include #include #include "Mem_Mngr.hh" +#include "ErrorHandling.hh" #define NEW_ALLOC #define MARKOVITZ @@ -54,8 +55,8 @@ class SparseMatrix { public: SparseMatrix(); - int simulate_NG1(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 Block_number, int stack_solve_algo, unsigned int endo_name_length, char *P_endo_names); - bool simulate_NG(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 Block_number, int solve_algo); + 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 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); @@ -68,6 +69,8 @@ private: void Init_Matlab_Sparse_Simple(int Size, map, int>, int> &IM, mxArray *A_m, mxArray *b_m); void Simple_Init(int it_, int y_kmin, int y_kmax, int Size, std::map, int>, int> &IM); 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_); 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_); diff --git a/mex/sources/bytecode/bytecode.cc b/mex/sources/bytecode/bytecode.cc index 92054f04b..5415a8920 100644 --- a/mex/sources/bytecode/bytecode.cc +++ b/mex/sources/bytecode/bytecode.cc @@ -64,6 +64,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) #else load_global((char*)prhs[1]); #endif + //ErrorHandlingException error_handling; int i, row_y = 0, col_y = 0, row_x = 0, col_x = 0, nb_row_xd = 0; int steady_row_y, steady_col_y, steady_row_x, steady_col_x, steady_nb_row_xd; int y_kmin = 0, y_kmax = 0, y_decal = 0, periods = 1; @@ -121,40 +122,43 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) block = true; else { - mexPrintf("Unknown argument : "); - string f; - f = Get_Argument(prhs[i]); - f.append("\n"); - mexErrMsgTxt(f.c_str()); + ostringstream tmp; + tmp << " in main, unknown argument : " << Get_Argument(prhs[i]) << "\n"; + throw FatalExceptionHandling(tmp.str()); } } if (count_array_argument > 0 && count_array_argument < 4) { - mexPrintf("Missing arguments. All the following arguments have to be indicated y, x, params, it_\n"); - mexErrMsgTxt("end of bytecode\n"); + ostringstream tmp; + tmp << " in main, missing arguments. All the following arguments have to be indicated y, x, params, it_\n"; + throw FatalExceptionHandling(tmp.str()); } M_ = mexGetVariable("global", "M_"); if (M_ == NULL) { - mexPrintf("Global variable not found : "); - mexErrMsgTxt("M_ \n"); + ostringstream tmp; + tmp << " in main, global variable not found: M_\n"; + throw FatalExceptionHandling(tmp.str()); } /* Gets variables and parameters from global workspace of Matlab */ oo_ = mexGetVariable("global", "oo_"); if (oo_ == NULL) { - mexPrintf("Global variable not found : "); - mexErrMsgTxt("oo_ \n"); + ostringstream tmp; + tmp << " in main, global variable not found: oo_\n"; + throw FatalExceptionHandling(tmp.str()); } options_ = mexGetVariable("global", "options_"); if (options_ == NULL) { - mexPrintf("Global variable not found : "); - mexErrMsgTxt("options_ \n"); + ostringstream tmp; + tmp << " in main, global variable not found: options_\n"; + throw FatalExceptionHandling(tmp.str()); } if (!count_array_argument) params = mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "params"))); + double *steady_yd = NULL, *steady_xd = NULL; if (!steady_state) { @@ -246,30 +250,40 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) string f(fname); mxFree(fname); int nb_blocks = 0; - bool result = interprete.compute_blocks(f, f, steady_state, evaluate, block, nb_blocks); + double *pind; + bool result = true, no_error = true; + + try + { + result = interprete.compute_blocks(f, f, steady_state, evaluate, block, nb_blocks); + } + catch (GeneralExceptionHandling &feh) + { + DYN_MEX_FUNC_ERR_MSG_TXT(feh.GetErrorMsg().c_str()); + } + clock_t t1 = clock(); - if (!steady_state && !evaluate) + if (!steady_state && !evaluate && no_error) mexPrintf("Simulation Time=%f milliseconds\n", 1000.0*(double (t1)-double (t0))/double (CLOCKS_PER_SEC)); #ifndef DEBUG_EX - double *pind; if (nlhs > 0) { - plhs[0] = mxCreateDoubleMatrix(row_y, col_y, mxREAL); + plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL); pind = mxGetPr(plhs[0]); - if (evaluate) - for (i = 0; i < row_y*col_y; i++) - pind[i] = y[i]-ya[i]; + if (no_error) + pind[0] = 0; else - for (i = 0; i < row_y*col_y; i++) - pind[i] = y[i]; + pind[0] = 1; if (nlhs > 1) { - plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL); + plhs[1] = mxCreateDoubleMatrix(row_y, col_y, mxREAL); pind = mxGetPr(plhs[1]); - if (result) - pind[0] = 0; + if (evaluate) + for (i = 0; i < row_y*col_y; i++) + pind[i] = y[i]-ya[i]; else - pind[0] = 1; + for (i = 0; i < row_y*col_y; i++) + pind[i] = y[i]; if (nlhs > 2) { int jacob_field_number = 0, jacob_exo_field_number = 0, jacob_exo_det_field_number = 0, jacob_other_endo_field_number = 0; @@ -285,37 +299,22 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) mexPrintf("the structure has been created\n"); } else if (!mxIsStruct(block_structur)) - { - mexPrintf("The third output argument must be a structure\n"); - mexErrMsgTxt("end of bytecode\n"); - } + DYN_MEX_FUNC_ERR_MSG_TXT("Fatal error in bytecode: in main, the third output argument must be a structure\n"); else { mexPrintf("Adding Fields\n"); jacob_field_number = mxAddField(block_structur, "jacob"); if (jacob_field_number == -1) - { - mexPrintf("Cannot add extra field to the structArray\n"); - mexErrMsgTxt("end of bytecode\n"); - } + DYN_MEX_FUNC_ERR_MSG_TXT("Fatal error in bytecode: in main, cannot add extra field jacob to the structArray\n"); jacob_exo_field_number = mxAddField(block_structur, "jacob_exo"); if (jacob_exo_field_number == -1) - { - mexPrintf("Cannot add extra field to the structArray\n"); - mexErrMsgTxt("end of bytecode\n"); - } + DYN_MEX_FUNC_ERR_MSG_TXT("Fatal error in bytecode: in main, cannot add extra field jacob_exo to the structArray\n"); jacob_exo_det_field_number = mxAddField(block_structur, "jacob_exo_det"); if (jacob_exo_det_field_number == -1) - { - mexPrintf("Cannot add extra field to the structArray\n"); - mexErrMsgTxt("end of bytecode\n"); - } + DYN_MEX_FUNC_ERR_MSG_TXT("Fatal error in bytecode: in main, cannot add extra field jacob_exo_det to the structArray\n"); jacob_other_endo_field_number = mxAddField(block_structur, "jacob_other_endo"); if (jacob_other_endo_field_number == -1) - { - mexPrintf("Cannot add extra field to the structArray\n"); - mexErrMsgTxt("end of bytecode\n"); - } + DYN_MEX_FUNC_ERR_MSG_TXT("Fatal error in bytecode: in main, cannot add extra field jacob_other_endo to the structArray\n"); } for (int i = 0; i < nb_blocks; i++) {