diff --git a/mex/sources/bytecode/Interpreter.cc b/mex/sources/bytecode/Interpreter.cc
index df982f9a5..399f87619 100644
--- a/mex/sources/bytecode/Interpreter.cc
+++ b/mex/sources/bytecode/Interpreter.cc
@@ -16,6 +16,9 @@
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see .
*/
+#define _GLIBCXX_USE_C99_FENV_TR1 1
+#include
+
#include
#include
#include "Interpreter.hh"
@@ -23,6 +26,10 @@
#define SMALL 1.0e-5;
//#define DEBUG
+Interpreter::~Interpreter()
+{
+}
+
Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, double *x_arg, double *steady_y_arg, double *steady_x_arg,
double *direction_arg, int y_size_arg,
int nb_row_x_arg, int nb_row_xd_arg, int periods_arg, int y_kmin_arg, int y_kmax_arg,
@@ -60,37 +67,35 @@ Interpreter::error_location()
{
string tmp;
mxArray *M_;
- double * P;
- int R, C;
- M_ = mexGetVariable("global", "M_");
+ char * P;
+ unsigned int R, C;
stringstream Error_loc("in ");
switch(EQN_type)
{
case TemporaryTerm:
- if (EQN_block_number>1)
+ if (EQN_block_number > 1)
Error_loc << "temporary term " << EQN_equation+1 << " in block " << EQN_block+1;
else
Error_loc << "temporary term " << EQN_equation+1;
break;
case ModelEquation:
- if (EQN_block_number>1)
+ if (EQN_block_number > 1)
Error_loc << "equation " << EQN_equation+1 << " in block " << EQN_block+1;
else
Error_loc << "equation " << EQN_equation+1;
break;
case FirstEndoDerivative:
- if (EQN_block_number>1)
- Error_loc << "first order derivative " << EQN_equation+1 << " in block " << EQN_block+1 << " with respect to endogenous variable ";
+ M_ = mexGetVariable("global", "M_");
+ if (EQN_block_number > 1)
+ Error_loc << "first order derivative of equation " << EQN_equation+1 << " in block " << EQN_block+1 << " with respect to endogenous variable ";
else
- Error_loc << "first order derivative " << EQN_equation+1 << " with respect to endogenous variable ";
- R = mxGetN(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "endo_names")));
- C = mxGetM(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "endo_names")));
- P = (double*)mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "endo_names")));
- if(EQN_dvar1first)
@@ -622,7 +648,20 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
#endif
break;
case oDivide:
- Stack.push(v1 / v2);
+ double r;
+ r = v1 / v2;
+ if (fetestexcept(FE_DIVBYZERO))
+ {
+ if (error_not_printed)
+ {
+ mexPrintf("--------------------------------------\n Error: Divide by zero with %5.15f/%5.15f\n in %s \n--------------------------------------\n", v1, v2,error_location().c_str());
+ error_not_printed = false;
+ r = 1e70;
+ }
+ feclearexcept (FE_ALL_EXCEPT);
+ res1 = NAN;
+ }
+ Stack.push(r);
#ifdef DEBUG
tmp_out << " |" << v1 << "/" << v2 << "|";
#endif
@@ -713,7 +752,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
#endif
break;
case oLog10:
- Stack.push(log10(v1));
+ Stack.push(log10_1(v1));
#ifdef DEBUG
tmp_out << " |log10(" << v1 << ")|";
#endif
@@ -1042,6 +1081,16 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
if (cvg)
continue;
y[Block_Contain[0].Variable] += -r[0]/g1[0];
+ if (fetestexcept(FE_DIVBYZERO))
+ {
+ 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;
+ }
+ feclearexcept (FE_ALL_EXCEPT);
+ res1 = NAN;
+ }
iter++;
}
if (!cvg)
@@ -1072,6 +1121,16 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
if (cvg)
continue;
y[Per_y_+Block_Contain[0].Variable] += -r[0]/g1[0];
+ if (fetestexcept(FE_DIVBYZERO))
+ {
+ 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;
+ }
+ feclearexcept (FE_ALL_EXCEPT);
+ res1 = NAN;
+ }
iter++;
}
if (!cvg)
@@ -1103,6 +1162,16 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
if (cvg)
continue;
y[Block_Contain[0].Variable] += -r[0]/g1[0];
+ if (fetestexcept(FE_DIVBYZERO))
+ {
+ 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;
+ }
+ feclearexcept (FE_ALL_EXCEPT);
+ res1 = NAN;
+ }
iter++;
}
if (!cvg)
@@ -1132,6 +1201,16 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
if (cvg)
continue;
y[Per_y_+Block_Contain[0].Variable] += -r[0]/g1[0];
+ if (fetestexcept(FE_DIVBYZERO))
+ {
+ 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;
+ }
+ feclearexcept (FE_ALL_EXCEPT);
+ res1 = NAN;
+ }
iter++;
}
if (!cvg)
@@ -1186,8 +1265,13 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
cvg = false;
if (cvg)
continue;
- result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true);
+ int prev_iter = iter;
+ result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, EQN_block_number);
iter++;
+ if (iter > prev_iter)
+ {
+
+ }
}
if (!cvg or !result)
{
@@ -1204,7 +1288,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);
cvg = false;
- result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true);
+ result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, EQN_block_number);
if (!result)
{
mexPrintf("Convergence not achieved in block %d\n", Block_Count);
@@ -1253,7 +1337,8 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
cvg = false;
if (cvg)
continue;
- result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false);
+ 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);
iter++;
}
if (!cvg)
@@ -1274,7 +1359,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);
cvg = false;
- result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false);
+ result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, EQN_block_number);
}
}
}
@@ -1326,7 +1411,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
cvg = false;
if (cvg)
continue;
- result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true);
+ result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, EQN_block_number);
iter++;
}
if (!cvg or !result)
@@ -1344,7 +1429,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);
cvg = false;
- result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true);
+ result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, EQN_block_number);
if (!result)
{
mexPrintf("Convergence not achieved in block %d\n", Block_Count);
@@ -1393,7 +1478,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
cvg = false;
if (cvg)
continue;
- result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false);
+ result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, EQN_block_number);
iter++;
}
if (!cvg)
@@ -1412,7 +1497,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);
cvg = false;
- result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false);
+ result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, EQN_block_number);
}
}
}
@@ -1485,7 +1570,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
else
cvg = (max_res < solve_tolf);
u_count = u_count_saved;
- simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter, minimal_solving_periods);
+ simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter, minimal_solving_periods, EQN_block_number);
iter++;
}
if (!cvg)
@@ -1517,7 +1602,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);
+ simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter, minimal_solving_periods, EQN_block_number);
}
mxFree(r);
mxFree(y_save);
diff --git a/mex/sources/bytecode/Interpreter.hh b/mex/sources/bytecode/Interpreter.hh
index c2eb071f5..90dc1a013 100644
--- a/mex/sources/bytecode/Interpreter.hh
+++ b/mex/sources/bytecode/Interpreter.hh
@@ -54,6 +54,7 @@ private:
protected:
double pow1(double a, double b);
double log1(double a);
+ double log10_1(double a);
string error_location();
void compute_block_time(int Per_u_, bool evaluate, int block_num);
void evaluate_a_block(const int size, const int type, string bin_basename, bool steady_state, int block_num,
@@ -77,7 +78,7 @@ private:
string filename;
int minimal_solving_periods;
public:
-
+ ~Interpreter();
Interpreter(double *params_arg, double *y_arg, double *ya_arg, double *x_arg, double *steady_y_arg, double *steady_x_arg,
double *direction_arg, int y_size_arg, int nb_row_x_arg,
int nb_row_xd_arg, int periods_arg, int y_kmin_arg, int y_kmax_arg, int maxit_arg_, double solve_tolf_arg, int size_o_direction_arg,
diff --git a/mex/sources/bytecode/SparseMatrix.cc b/mex/sources/bytecode/SparseMatrix.cc
index 7f67fad06..13a53dea8 100644
--- a/mex/sources/bytecode/SparseMatrix.cc
+++ b/mex/sources/bytecode/SparseMatrix.cc
@@ -17,6 +17,9 @@
* along with Dynare. If not, see .
*/
+#define _GLIBCXX_USE_C99_FENV_TR1 1
+#include
+
#include
#include
#include
@@ -994,7 +997,7 @@ SparseMatrix::simple_bksub(int it_, int Size, double slowc_l)
}
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)
+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 i, j, k;
int pivj = 0, pivk = 0;
@@ -1082,7 +1085,12 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
}
}
}
+#ifdef GLOBAL_CONVERGENCE
+ if (iter == 1)
+
+#else
slowc_save /= 2;
+#endif
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];
@@ -1152,42 +1160,12 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
}
first = first->NZE_C_N;
}
- double markovitz = 0, markovitz_max = -9e70;
- if (!one)
- {
- for (j = 0; j < l; j++)
- {
- markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double (NR[j])/double (N_max)));
- if (markovitz > markovitz_max)
- {
- piv = piv_v[j];
- pivj = pivj_v[j]; //Line number
- pivk = pivk_v[j]; //positi
- markovitz_max = markovitz;
- }
- }
- }
- else
- {
- for (j = 0; j < l; j++)
- {
- markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double (NR[j])/double (N_max)));
- if (markovitz > markovitz_max && NR[j] == 1)
- {
- piv = piv_v[j];
- pivj = pivj_v[j]; //Line number
- pivk = pivk_v[j]; //positi
- markovitz_max = markovitz;
- }
- }
- }
- pivot[i] = pivj;
- pivotk[i] = pivk;
- pivotv[i] = piv;
- line_done[pivj] = true;
if (piv_abs < eps)
{
- mexPrintf("Error: singular system in Simulate_NG in block %d\n", blck+1);
+ 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);
@@ -1203,6 +1181,98 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
mexErrMsgTxt(filename.c_str());
}
}
+ double markovitz = 0, markovitz_max = -9e70;
+ int NR_max = 0;
+ //mexPrintf("l=%d\n",l);
+ 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 (fetestexcept(FE_DIVBYZERO))
+ {
+ if (error_not_printed)
+ {
+ mexPrintf("--------------------------------------\n Error: Divide by zero in simul_NG(1) piv_abs=%f and N_max=%d \n--------------------------------------\n", piv_abs, N_max);
+ error_not_printed = false;
+ markovitz = 1e70;
+ }
+ feclearexcept (FE_ALL_EXCEPT);
+ res1 = NAN;
+ }
+ //mexPrintf("piv_v[j]=%f NR[j]=%d markovitz=%f markovitz_max=%f\n", piv_v[j], NR[j], markovitz, markovitz_max);
+ 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 (fetestexcept(FE_DIVBYZERO))
+ {
+ if (error_not_printed)
+ {
+ mexPrintf("--------------------------------------\n Error: Divide by zero in simul_NG(2) piv_abs=%f and N_max=%d \n--------------------------------------\n", piv_abs, N_max);
+ error_not_printed = false;
+ markovitz = 1e70;
+ }
+ feclearexcept (FE_ALL_EXCEPT);
+ res1 = NAN;
+ }
+ if (/*markovitz > markovitz_max &&*/ 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];
+ //mexPrintf("forced piv = %f NR_max =%d\n",piv, NR_max);
+ }
+ }
+ }
+ 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;
+ 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++)
@@ -1426,7 +1496,7 @@ SparseMatrix::Check_the_Solution(int periods, int y_kmin, int y_kmax, int Size,
}
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)
+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)
{
/*Triangularisation at each period of a block using a simple gaussian Elimination*/
t_save_op_s *save_op_s;
@@ -1606,17 +1676,45 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
first = first->NZE_C_N;
}
double markovitz = 0, markovitz_max = -9e70;
+ int NR_max = 0;
+ //mexPrintf("l=%d\n",l);
if (!one)
{
for (j = 0; j < l; j++)
{
- markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double (NR[j])/double (N_max)));
+ 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 (fetestexcept(FE_DIVBYZERO))
+ {
+ if (error_not_printed)
+ {
+ mexPrintf("--------------------------------------\n Error: Divide by zero in simul_NG(1) piv_abs=%f and N_max=%d \n--------------------------------------\n", piv_abs, N_max);
+ error_not_printed = false;
+ markovitz = 1e70;
+ }
+ feclearexcept (FE_ALL_EXCEPT);
+ res1 = NAN;
+ }
+ //mexPrintf("piv_v[j]=%f NR[j]=%d markovitz=%f markovitz_max=%f\n", piv_v[j], NR[j], markovitz, markovitz_max);
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];
}
}
}
@@ -1624,16 +1722,46 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
{
for (j = 0; j < l; j++)
{
- markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log((double)(NR[j]/N_max)));
- if (markovitz > markovitz_max && NR[j] == 1)
+ 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 (fetestexcept(FE_DIVBYZERO))
+ {
+ if (error_not_printed)
+ {
+ mexPrintf("--------------------------------------\n Error: Divide by zero in simul_NG(2) piv_abs=%f and N_max=%d \n--------------------------------------\n", piv_abs, N_max);
+ error_not_printed = false;
+ markovitz = 1e70;
+ }
+ feclearexcept (FE_ALL_EXCEPT);
+ res1 = NAN;
+ }
+ if (/*markovitz > markovitz_max &&*/ NR[j] == 1)
{
piv = piv_v[j];
pivj = pivj_v[j]; //Line number
- pivk = pivk_v[j]; //position
+ pivk = pivk_v[j]; //positi
markovitz_max = markovitz;
+ NR_max = NR[j];
+ //mexPrintf("forced piv = %f NR_max =%d\n",piv, NR_max);
}
}
}
+ 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;
@@ -1667,7 +1795,10 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
}
if (piv_abs < eps)
{
- mexPrintf("Error: singular system in Simulate_NG1\n");
+ 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());
diff --git a/mex/sources/bytecode/SparseMatrix.hh b/mex/sources/bytecode/SparseMatrix.hh
index fd24b22f3..f2cf06dc5 100644
--- a/mex/sources/bytecode/SparseMatrix.hh
+++ b/mex/sources/bytecode/SparseMatrix.hh
@@ -99,8 +99,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);
- 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 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);
+ 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);
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);
@@ -182,6 +182,7 @@ protected:
int start_compare;
int restart;
bool error_not_printed;
+ double g_lambda1, g_lambda2, gp_0;
};
#endif
diff --git a/mex/sources/bytecode/bytecode.cc b/mex/sources/bytecode/bytecode.cc
index cef457ec4..59292040b 100644
--- a/mex/sources/bytecode/bytecode.cc
+++ b/mex/sources/bytecode/bytecode.cc
@@ -21,9 +21,11 @@
// simulate.cc //
// simulate file designed for GNU GCC C++ compiler //
////////////////////////////////////////////////////////////////////////
+#define _GLIBCXX_USE_C99_FENV_TR1 1
+#include
+
#include
-
#include "Interpreter.hh"
#ifndef DEBUG_EX
# include "mex.h"
@@ -33,10 +35,12 @@
#include "Mem_Mngr.hh"
+
#ifdef DEBUG_EX
using namespace std;
# include
+
string
Get_Argument(const char *argv)
{
@@ -213,6 +217,30 @@ Get_Argument(const mxArray *prhs)
return f;
}
+void fpe_handler(int) {
+ int e;
+ mexPrintf("caught FPE, exiting.\n");
+ e = fetestexcept(FE_ALL_EXCEPT);
+ if (!e) {
+ mexPrintf("no exception information set\n");
+ }
+ if (e & FE_DIVBYZERO) {
+ mexPrintf("divide by zero\n");
+ }
+ if (e & FE_INVALID) {
+ mexPrintf("invalide operand\n");
+ }
+ if (e & FE_UNDERFLOW) {
+ mexPrintf("underflow\n");
+ }
+ if (e & FE_OVERFLOW) {
+ mexPrintf("overflow\n");
+ }
+ exit(1);
+}
+
+
+
/* The gateway routine */
void
mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
@@ -225,6 +253,14 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
double *direction;
bool steady_state = false;
bool evaluate = false;
+ fexcept_t *flagp;
+ flagp = (fexcept_t*) mxMalloc(sizeof(fexcept_t));
+ if(fegetexceptflag(flagp, FE_ALL_EXCEPT))
+ mexPrintf("fegetexceptflag failed\n");
+ if(fesetexceptflag(flagp,FE_INVALID | FE_DIVBYZERO))
+ mexPrintf("fesetexceptflag failed\n");
+ mxFree(flagp);
+ feclearexcept (FE_ALL_EXCEPT);
for (i = 0; i < nrhs; i++)
{
if (Get_Argument(prhs[i]) == "static")