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

Uses the initial method to manage the floating exceptions

time-shift
Ferhat MIHOUBI 2010-02-05 18:50:57 +01:00
parent f7ac31b58a
commit cea26af06e
3 changed files with 23 additions and 54 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -16,8 +16,8 @@
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#define _GLIBCXX_USE_C99_FENV_TR1 1
#include <cfenv>
//#define _GLIBCXX_USE_C99_FENV_TR1 1
//#include <cfenv>
#include <cstring>
#include <sstream>
@ -112,7 +112,7 @@ double
Interpreter::pow1(double a, double b)
{
double r = pow_(a, b);
if (fetestexcept(FE_INVALID))
if (isnan(r) || isinf(r))
{
if (error_not_printed)
{
@ -120,7 +120,6 @@ Interpreter::pow1(double a, double b)
error_not_printed = false;
r = 0.0000000000000000000000001;
}
feclearexcept (FE_ALL_EXCEPT);
res1 = NAN;
return r;
}
@ -131,7 +130,7 @@ double
Interpreter::log1(double a)
{
double r = log(a);
if (fetestexcept(FE_INVALID | FE_DIVBYZERO))
if (isnan(r) || isinf(r))
{
if (error_not_printed)
{
@ -139,7 +138,6 @@ Interpreter::log1(double a)
error_not_printed = false;
r = 0.0000000000000000000000001;
}
feclearexcept (FE_ALL_EXCEPT);
res1 = NAN;
return r;
}
@ -150,7 +148,7 @@ double
Interpreter::log10_1(double a)
{
double r = log(a);
if (fetestexcept(FE_INVALID | FE_DIVBYZERO))
if (isnan(r) || isinf(r))
{
if (error_not_printed)
{
@ -158,7 +156,6 @@ Interpreter::log10_1(double a)
error_not_printed = false;
r = 0.0000000000000000000000001;
}
feclearexcept (FE_ALL_EXCEPT);
res1 = NAN;
return r;
}
@ -655,7 +652,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
case oDivide:
double r;
r = v1 / v2;
if (fetestexcept(FE_DIVBYZERO))
if (isinf(r))
{
if (error_not_printed)
{
@ -663,7 +660,6 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
error_not_printed = false;
r = 1e70;
}
feclearexcept (FE_ALL_EXCEPT);
res1 = NAN;
}
Stack.push(r);
@ -1086,14 +1082,13 @@ 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 (isinf(y[Block_Contain[0].Variable]))
{
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++;
@ -1126,14 +1121,13 @@ 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 (isinf(y[Per_y_+Block_Contain[0].Variable]))
{
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++;
@ -1167,14 +1161,13 @@ 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 (isinf(y[Block_Contain[0].Variable]))
{
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++;
@ -1206,14 +1199,13 @@ 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 (isinf(y[Per_y_+Block_Contain[0].Variable]))
{
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++;
@ -1273,7 +1265,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, true/*false*/, cvg, iter, true, EQN_block_number);
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)
{

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

@ -17,8 +17,8 @@
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#define _GLIBCXX_USE_C99_FENV_TR1 1
#include <cfenv>
//#define _GLIBCXX_USE_C99_FENV_TR1 1
//#include <cfenv>
#include <cstring>
#include <ctime>
@ -1718,7 +1718,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
}
else
markovitz = fabs(piv_v[j])/piv_abs;
if (fetestexcept(FE_DIVBYZERO))
/*if (fetestexcept(FE_DIVBYZERO))
{
if (error_not_printed)
{
@ -1728,7 +1728,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
}
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)
{
@ -1758,7 +1758,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
}
else
markovitz = fabs(piv_v[j])/piv_abs;
if (fetestexcept(FE_DIVBYZERO))
/*if (fetestexcept(FE_DIVBYZERO))
{
if (error_not_printed)
{
@ -1768,7 +1768,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
}
feclearexcept (FE_ALL_EXCEPT);
res1 = NAN;
}
}*/
if (/*markovitz > markovitz_max &&*/ NR[j] == 1)
{
piv = piv_v[j];

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

@ -21,8 +21,8 @@
// simulate.cc //
// simulate file designed for GNU GCC C++ compiler //
////////////////////////////////////////////////////////////////////////
#define _GLIBCXX_USE_C99_FENV_TR1 1
#include <cfenv>
//#define _GLIBCXX_USE_C99_FENV_TR1 1
//#include <cfenv>
#include <cstring>
@ -55,14 +55,14 @@ main(int argc, const char *argv[])
bool steady_state = false;
bool evaluate = false;
printf("argc=%d\n", argc);
fexcept_t *flagp;
/*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);
feclearexcept (FE_ALL_EXCEPT);*/
if (argc < 2)
{
mexPrintf("model filename expected\n");
@ -229,29 +229,6 @@ 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
@ -265,14 +242,14 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
double *direction;
bool steady_state = false;
bool evaluate = false;
fexcept_t *flagp;
/*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);
feclearexcept (FE_ALL_EXCEPT);*/
for (i = 0; i < nrhs; i++)
{
if (Get_Argument(prhs[i]) == "static")