stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

- In bytecode, the MATLAB function "mexErrMsgTxt" has been replaced by an exception handling see ticket #137

time-shift
Ferhat Mihoubi 2010-09-24 12:52:58 +02:00
parent bc089f732e
commit df1b1e4ed0
9 changed files with 1365 additions and 1256 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
matlab/resid.m
View File

@ -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,...

6
matlab/simul.m
View File

@ -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')

6
matlab/steady_.m
View File

@ -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,...

122
mex/sources/bytecode/ErrorHandling.hh Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef ERROR_HANDLING
#define ERROR_HANDLING
#include <cstring>
#include <iostream>
#include <sstream>
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

325
mex/sources/bytecode/Interpreter.cc
View File

@ -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);

9
mex/sources/bytecode/Interpreter.hh
View File

@ -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_type> Block_Contain;

2052
mex/sources/bytecode/SparseMatrix.cc
View File

File diff suppressed because it is too large Load Diff

7
mex/sources/bytecode/SparseMatrix.hh
View File

@ -26,6 +26,7 @@
#include <map>
#include <ctime>
#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<pair<pair<int, int>, int>, int> &IM, mxArray *A_m, mxArray *b_m);
void Simple_Init(int it_, int y_kmin, int y_kmax, int Size, std::map<std::pair<std::pair<int, int>, int>, int> &IM);
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_);

91
mex/sources/bytecode/bytecode.cc
View File

@ -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++)
{