diff --git a/mex/sources/bytecode/ErrorHandling.hh b/mex/sources/bytecode/ErrorHandling.hh
index 0af66aee3453e2b0c3b5426802bd96fb4e6810c5..00e6c8cdb8569b0194cde357020c7bc80211663b 100644
--- a/mex/sources/bytecode/ErrorHandling.hh
+++ b/mex/sources/bytecode/ErrorHandling.hh
@@ -33,57 +33,59 @@ struct GeneralException
struct FloatingPointException : public GeneralException
{
- FloatingPointException(const string &details) :
- GeneralException {"Floating point error: " + details}
+ const string location;
+ FloatingPointException(const string &details, string location_arg) :
+ GeneralException {"Floating point error: " + details},
+ location {move(location_arg)}
{
}
};
struct LogException : public FloatingPointException
{
- LogException(double value) :
+ LogException(double value, string location_arg) :
FloatingPointException { [=]
{
// We don’t use std::to_string(), because it uses fixed formatting
ostringstream s;
s << "log(X) with X=" << defaultfloat << value;
return s.str();
- }() }
+ }(), move(location_arg) }
{
}
};
struct Log10Exception : public FloatingPointException
{
- Log10Exception(double value) :
+ Log10Exception(double value, string location_arg) :
FloatingPointException { [=]
{
// We don’t use std::to_string(), because it uses fixed formatting
ostringstream s;
s << "log10(X) with X=" << defaultfloat << value;
return s.str();
- }() }
+ }(), move(location_arg) }
{
}
};
struct DivideException : public FloatingPointException
{
- DivideException(double divisor) :
+ DivideException(double v1, double v2, string location_arg) :
FloatingPointException { [=]
{
// We don’t use std::to_string(), because it uses fixed formatting
ostringstream s;
- s << "a/X with X=" << defaultfloat << divisor;
+ s << "a/X with a=" << defaultfloat << v1 << " and X= " << v2;
return s.str();
- }() }
+ }(), move(location_arg) }
{
}
};
struct PowException : public FloatingPointException
{
- PowException(double base, double exponent) :
+ PowException(double base, double exponent, string location_arg) :
FloatingPointException { [=]
{
// We don’t use std::to_string(), because it uses fixed formatting
@@ -92,7 +94,7 @@ struct PowException : public FloatingPointException
if (fabs(base) <= 1e-10)
s << " and a=" << exponent;
return s.str();
- }() }
+ }(), move(location_arg) }
{
}
};
diff --git a/mex/sources/bytecode/Evaluate.cc b/mex/sources/bytecode/Evaluate.cc
index 0772fe2c999dc020db025429456c247ed1472d49..05ba368a36d9bf2983a3c732334361231a715d08 100644
--- a/mex/sources/bytecode/Evaluate.cc
+++ b/mex/sources/bytecode/Evaluate.cc
@@ -21,6 +21,7 @@
#include <cmath>
#include <limits>
#include <stack>
+#include <cfenv>
#include <dynmex.h>
@@ -987,64 +988,8 @@ Evaluate::print_expression(const Evaluate::it_code_type &expr_begin, const optio
return { get<0>(Stack.top()), it_code };
}
-double
-Evaluate::pow1(double a, double b)
-{
- double r = pow(a, b);
- if (isnan(r) || isinf(r))
- {
- res1 = std::numeric_limits<double>::quiet_NaN();
- r = 0.0000000000000000000000001;
- if (print_error)
- throw PowException{a, b};
- }
- return r;
-}
-
-double
-Evaluate::divide(double a, double b)
-{
- double r = a / b;
- if (isnan(r) || isinf(r))
- {
- res1 = std::numeric_limits<double>::quiet_NaN();
- r = 1e70;
- if (print_error)
- throw DivideException{b};
- }
- return r;
-}
-
-double
-Evaluate::log1(double a)
-{
- double r = log(a);
- if (isnan(r) || isinf(r))
- {
- res1 = std::numeric_limits<double>::quiet_NaN();
- r = -1e70;
- if (print_error)
- throw LogException{a};
- }
- return r;
-}
-
-double
-Evaluate::log10_1(double a)
-{
- double r = log(a);
- if (isnan(r) || isinf(r))
- {
- res1 = std::numeric_limits<double>::quiet_NaN();
- r = -1e70;
- if (print_error)
- throw Log10Exception{a};
- }
- return r;
-}
-
void
-Evaluate::compute_block_time(int Per_u_, bool evaluate, bool no_derivative)
+Evaluate::evaluateBlock(int Per_u_, bool evaluate, bool no_derivative)
{
auto it_code { currentBlockBeginning() };
int var{0}, lag{0};
@@ -1619,20 +1564,13 @@ Evaluate::compute_block_time(int Per_u_, bool evaluate, bool no_derivative)
#endif
break;
case BinaryOpcode::divide:
- double tmp;
-#ifdef DEBUG
- mexPrintf("v1=%f / v2=%f\n", v1, v2);
-#endif
- try
- {
- tmp = divide(v1, v2);
- }
- catch (FloatingPointException &fpeh)
- {
- mexPrintf("%s\n %s\n", fpeh.message.c_str(), error_location(it_code_expr, it_code, steady_state, it_).c_str());
- go_on = false;
- }
- Stack.push(tmp);
+ {
+ feclearexcept(FE_ALL_EXCEPT);
+ double tmp {v1 / v2};
+ if (fetestexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW))
+ throw DivideException{v1, v2, error_location(it_code_expr, it_code, steady_state, it_)};
+ Stack.push(tmp);
+ }
#ifdef DEBUG
tmp_out << " |" << v1 << "/" << v2 << "|";
#endif
@@ -1674,20 +1612,13 @@ Evaluate::compute_block_time(int Per_u_, bool evaluate, bool no_derivative)
#endif
break;
case BinaryOpcode::power:
-#ifdef DEBUG
- mexPrintf("pow\n");
-#endif
- try
- {
- tmp = pow1(v1, v2);
- }
- catch (FloatingPointException &fpeh)
- {
- mexPrintf("%s\n %s\n", fpeh.message.c_str(), error_location(it_code_expr, it_code, steady_state, it_).c_str());
- go_on = false;
- }
- Stack.push(tmp);
-
+ {
+ feclearexcept(FE_ALL_EXCEPT);
+ double tmp {pow(v1, v2)};
+ if (fetestexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW))
+ throw PowException{v1, v2, error_location(it_code_expr, it_code, steady_state, it_)};
+ Stack.push(tmp);
+ }
#ifdef DEBUG
tmp_out << " |" << v1 << "^" << v2 << "|";
#endif
@@ -1704,7 +1635,10 @@ Evaluate::compute_block_time(int Per_u_, bool evaluate, bool no_derivative)
Stack.push(0.0);
else
{
- double dxp = pow1(v1, v2-derivOrder);
+ feclearexcept(FE_ALL_EXCEPT);
+ double dxp {pow(v1, v2-derivOrder)};
+ if (fetestexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW))
+ throw PowException{v1, v2-derivOrder, error_location(it_code_expr, it_code, steady_state, it_)};
for (int i = 0; i < derivOrder; i++)
dxp *= v2--;
Stack.push(dxp);
@@ -1767,33 +1701,25 @@ Evaluate::compute_block_time(int Per_u_, bool evaluate, bool no_derivative)
#endif
break;
case UnaryOpcode::log:
- double tmp;
- try
- {
- tmp = log1(v1);
- }
- catch (FloatingPointException &fpeh)
- {
- mexPrintf("%s\n %s\n", fpeh.message.c_str(), error_location(it_code_expr, it_code, steady_state, it_).c_str());
- go_on = false;
- }
- Stack.push(tmp);
-
+ {
+ feclearexcept(FE_ALL_EXCEPT);
+ double tmp {log(v1)};
+ if (fetestexcept(FE_DIVBYZERO | FE_INVALID))
+ throw LogException{v1, error_location(it_code_expr, it_code, steady_state, it_)};
+ Stack.push(tmp);
+ }
#ifdef DEBUG
tmp_out << " |log(" << v1 << ")|";
#endif
break;
case UnaryOpcode::log10:
- try
- {
- tmp = log10_1(v1);
- }
- catch (FloatingPointException &fpeh)
- {
- mexPrintf("%s\n %s\n", fpeh.message.c_str(), error_location(it_code_expr, it_code, steady_state, it_).c_str());
- go_on = false;
- }
- Stack.push(tmp);
+ {
+ feclearexcept(FE_ALL_EXCEPT);
+ double tmp {log10(v1)};
+ if (fetestexcept(FE_DIVBYZERO | FE_INVALID))
+ throw Log10Exception{v1, error_location(it_code_expr, it_code, steady_state, it_)};
+ Stack.push(tmp);
+ }
#ifdef DEBUG
tmp_out << " |log10(" << v1 << ")|";
#endif
diff --git a/mex/sources/bytecode/Evaluate.hh b/mex/sources/bytecode/Evaluate.hh
index 16b22aa326218daa8c6543ee9dd44b411f915052..d6a37e673e08cacce9cdba022dea8dd91e1ff370 100644
--- a/mex/sources/bytecode/Evaluate.hh
+++ b/mex/sources/bytecode/Evaluate.hh
@@ -89,14 +89,8 @@ protected:
double *g1, *r, *res;
vector<mxArray *> jacobian_block, jacobian_exo_block, jacobian_det_exo_block;
mxArray *GlobalTemporaryTerms;
- double pow1(double a, double b);
- double divide(double a, double b);
- double log1(double a);
- double log10_1(double a);
- void compute_block_time(int Per_u_, bool evaluate, bool no_derivatives);
+ void evaluateBlock(int Per_u_, bool evaluate, bool no_derivatives);
int it_;
- bool print_error;
- double res1;
int block_num; // Index of the current block
int size; // Size of the current block
diff --git a/mex/sources/bytecode/Interpreter.cc b/mex/sources/bytecode/Interpreter.cc
index 8f4e069710dcdeda96f90fefbe381db7d2bbe5f1..cba6fecd04d7915ee57b4a524dba2dfb90d42caf 100644
--- a/mex/sources/bytecode/Interpreter.cc
+++ b/mex/sources/bytecode/Interpreter.cc
@@ -22,6 +22,7 @@
#include <cstring>
#include <filesystem>
#include <numeric>
+#include <cfenv>
#include "Interpreter.hh"
@@ -34,7 +35,7 @@ Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, doub
string &filename_arg, int minimal_solving_periods_arg, int stack_solve_algo_arg, int solve_algo_arg,
bool global_temporary_terms_arg, bool print_arg, bool print_error_arg, mxArray *GlobalTemporaryTerms_arg,
bool steady_state_arg, bool block_decomposed_arg, bool print_it_arg, int col_x_arg, int col_y_arg, BasicSymbolTable &symbol_table_arg)
-: dynSparseMatrix {y_size_arg, y_kmin_arg, y_kmax_arg, print_it_arg, steady_state_arg, block_decomposed_arg, periods_arg, minimal_solving_periods_arg, symbol_table_arg}
+: dynSparseMatrix {y_size_arg, y_kmin_arg, y_kmax_arg, print_it_arg, steady_state_arg, block_decomposed_arg, periods_arg, minimal_solving_periods_arg, symbol_table_arg, print_error_arg}
{
params = params_arg;
y = y_arg;
@@ -61,7 +62,6 @@ Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, doub
col_x = col_x_arg;
col_y = col_y_arg;
GlobalTemporaryTerms = GlobalTemporaryTerms_arg;
- print_error = print_error_arg;
print_it = print_it_arg;
}
@@ -110,7 +110,14 @@ Interpreter::solve_simple_one_periods()
{
slowc /= 1.5;
mexPrintf("Reducing the path length in Newton step slowc=%f\n", slowc);
- y[Block_Contain[0].Variable + Per_y_] = ya - slowc * divide(r[0], g1[0]);
+ feclearexcept(FE_ALL_EXCEPT);
+ y[Block_Contain[0].Variable + Per_y_] = ya - slowc * (r[0] / g1[0]);
+ if (fetestexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW))
+ {
+ res1 = numeric_limits<double>::quiet_NaN();
+ if (print_error)
+ mexPrintf(" Singularity in block %d", block_num+1);
+ }
}
}
}
@@ -119,14 +126,13 @@ Interpreter::solve_simple_one_periods()
cvg = (fabs(rr) < solve_tolf);
if (cvg)
continue;
- try
+ feclearexcept(FE_ALL_EXCEPT);
+ y[Block_Contain[0].Variable + Per_y_] += -slowc * (rr / g1[0]);
+ if (fetestexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW))
{
- y[Block_Contain[0].Variable + Per_y_] += -slowc *divide(rr, g1[0]);
- }
- catch (FloatingPointException &fpeh)
- {
- mexPrintf("%s\n \n", fpeh.message.c_str());
- mexPrintf(" Singularity in block %d", block_num+1);
+ res1 = numeric_limits<double>::quiet_NaN();
+ if (print_error)
+ mexPrintf(" Singularity in block %d", block_num+1);
}
iter++;
}
diff --git a/mex/sources/bytecode/SparseMatrix.cc b/mex/sources/bytecode/SparseMatrix.cc
index 303124ccc9c6b4885c1c23cdd7226a816140742b..1b10ac0766966ea6469ed2dc11f7c8dd0b71b945 100644
--- a/mex/sources/bytecode/SparseMatrix.cc
+++ b/mex/sources/bytecode/SparseMatrix.cc
@@ -25,11 +25,13 @@
#include "SparseMatrix.hh"
dynSparseMatrix::dynSparseMatrix(int y_size_arg, int y_kmin_arg, int y_kmax_arg, bool print_it_arg, bool steady_state_arg, bool block_decomposed_arg, int periods_arg,
- int minimal_solving_periods_arg, BasicSymbolTable &symbol_table_arg) :
+ int minimal_solving_periods_arg, BasicSymbolTable &symbol_table_arg,
+ bool print_error_arg) :
Evaluate {y_size_arg, y_kmin_arg, y_kmax_arg, steady_state_arg, periods_arg, symbol_table_arg},
block_decomposed {block_decomposed_arg},
minimal_solving_periods {minimal_solving_periods_arg},
- print_it {print_it_arg}
+ print_it {print_it_arg},
+ print_error {print_error_arg}
{
pivotva = nullptr;
g_save_op = nullptr;
@@ -1818,6 +1820,21 @@ dynSparseMatrix::Sparse_transpose(const mxArray *A_m)
return C_m;
}
+void
+dynSparseMatrix::compute_block_time(int Per_u_, bool evaluate, bool no_derivatives)
+{
+ try
+ {
+ evaluateBlock(Per_u_, evaluate, no_derivatives);
+ }
+ catch (FloatingPointException &e)
+ {
+ res1 = numeric_limits<double>::quiet_NaN();
+ if (print_error)
+ mexPrintf("%s\n %s\n", e.message.c_str(), e.location.c_str());
+ }
+}
+
bool
dynSparseMatrix::compute_complete(bool no_derivatives, double &_res1, double &_res2, double &_max_res, int &_max_res_idx)
{
diff --git a/mex/sources/bytecode/SparseMatrix.hh b/mex/sources/bytecode/SparseMatrix.hh
index 12c5c76d141f33768808dcb5aac112eef3c19108..0f018e702b5227458df5d8251d571e194364af85 100644
--- a/mex/sources/bytecode/SparseMatrix.hh
+++ b/mex/sources/bytecode/SparseMatrix.hh
@@ -68,7 +68,7 @@ constexpr double mem_increasing_factor = 1.1;
class dynSparseMatrix : public Evaluate
{
public:
- dynSparseMatrix(int y_size_arg, int y_kmin_arg, int y_kmax_arg, bool print_it_arg, bool steady_state_arg, bool block_decomposed_arg, int periods_arg, int minimal_solving_periods_arg, BasicSymbolTable &symbol_table_arg);
+ dynSparseMatrix(int y_size_arg, int y_kmin_arg, int y_kmax_arg, bool print_it_arg, bool steady_state_arg, bool block_decomposed_arg, int periods_arg, int minimal_solving_periods_arg, BasicSymbolTable &symbol_table_arg, bool print_error_arg);
void Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin, int y_kmax, int Size, int periods, bool cvg, int minimal_solving_periods, int stack_solve_algo, const vector_table_conditional_local_type &vector_table_conditional_local);
void Simulate_Newton_One_Boundary(bool forward);
void fixe_u(double **u, int u_count_int, int max_lag_plus_max_lead_plus_1);
@@ -198,13 +198,16 @@ protected:
int maxit_;
double *direction;
double solve_tolf;
- double res2, max_res;
+ double res1, res2, max_res;
int max_res_idx;
int *index_vara;
+ void compute_block_time(int Per_u_, bool evaluate, bool no_derivatives);
bool compute_complete(bool no_derivatives, double &res1, double &res2, double &max_res, int &max_res_idx);
bool compute_complete(double lambda, double *crit);
+
+ bool print_error; // Whether to stop processing on floating point exceptions
};
#endif // _SPARSEMATRIX_HH