Merge pull request #424 from FerhatMihoubi/master

1 commit

mex/sources/bytecode/ErrorHandling.hh

@@ -59,6 +59,7 @@
 #define isinf(x) (!_finite(x))
 #define fpu_error(x) (isinf(x) || isnan(x))
 
+#define finite(x) _finite(x)
 
 class MSVCpp_missings
 {

mex/sources/bytecode/Evaluate.cc

@@ -19,7 +19,7 @@ Evaluate::Evaluate()
   block = -1;
 }
 
-Evaluate::Evaluate(const int y_size_arg, const int y_kmin_arg, const int y_kmax_arg, const bool print_it_arg, const bool steady_state_arg, const int periods_arg, const int minimal_solving_periods_arg):
+Evaluate::Evaluate(const int y_size_arg, const int y_kmin_arg, const int y_kmax_arg, const bool print_it_arg, const bool steady_state_arg, const int periods_arg, const int minimal_solving_periods_arg, const double slowc_arg):
 print_it(print_it_arg),  minimal_solving_periods(minimal_solving_periods_arg)
 {
   symbol_table_endo_nbr = 0;
@@ -32,6 +32,7 @@ print_it(print_it_arg),  minimal_solving_periods(minimal_solving_periods_arg)
   y_kmax  = y_kmax_arg;
   periods = periods_arg;
   steady_state = steady_state_arg;
+  slowc = slowc_arg;
 }
 
 double
@@ -1504,11 +1505,30 @@ Evaluate::solve_simple_one_periods()
 {
   bool cvg = false;
   int iter = 0;
+  double ya ;
+  double slowc_save = slowc;
+  res1 = 0;
   while (!(cvg || (iter > maxit_)))
     {
       it_code = start_code;
       Per_y_ = it_*y_size;
+      ya = y[Block_Contain[0].Variable + Per_y_];
       compute_block_time(0, false, false);
+      if (!finite(res1))
+        {
+          res1 = NAN;
+          while ((isinf(res1) || isnan(res1)) && (slowc > 1e-9) )
+            {
+              it_code = start_code;
+              compute_block_time(0, false, false);
+              if (!finite(res1))
+                {
+                  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]);
+                }
+            }
+        }
       double rr;
       rr = r[0];
       cvg = (fabs(rr) < solve_tolf);
@@ -1517,7 +1537,7 @@ Evaluate::solve_simple_one_periods()
         continue;
       try
         {
-          y[Block_Contain[0].Variable + Per_y_] += -divide(rr, g1[0]);
+          y[Block_Contain[0].Variable + Per_y_] += - slowc * divide(rr, g1[0]);
         }
       catch (FloatingPointExceptionHandling &fpeh)
         {
@@ -1526,6 +1546,7 @@ Evaluate::solve_simple_one_periods()
         }
       iter++;
     }
+  slowc = slowc_save;
   if (!cvg)
     {
       ostringstream tmp;

mex/sources/bytecode/Evaluate.hh

@@ -82,8 +82,9 @@ protected:
   bool Gaussian_Elimination, is_linear;
 public:
   bool steady_state;
+  double slowc;
   Evaluate();
-  Evaluate(const int y_size_arg, const int y_kmin_arg, const int y_kmax_arg, const bool print_it_arg, const bool steady_state_arg, const int periods_arg, const int minimal_solving_periods_arg);
+  Evaluate(const int y_size_arg, const int y_kmin_arg, const int y_kmax_arg, const bool print_it_arg, const bool steady_state_arg, const int periods_arg, const int minimal_solving_periods_arg, const double slowc);
   //typedef  void (Interpreter::*InterfpreterMemFn)(const int block_num, const int size, const bool steady_state, int it);
   void set_block(const int size_arg, const int type_arg, string file_name_arg, string bin_base_name_arg, const int block_num_arg,
           const bool is_linear_arg, const int symbol_table_endo_nbr_arg, const int Block_List_Max_Lag_arg, const int Block_List_Max_Lead_arg, const int u_count_int_arg, const int block_arg);

mex/sources/bytecode/Interpreter.cc

@@ -39,7 +39,7 @@ Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, doub
                          , const int CUDA_device_arg, cublasHandle_t cublas_handle_arg, cusparseHandle_t cusparse_handle_arg, cusparseMatDescr_t descr_arg
 #endif
                          )
-                         : dynSparseMatrix(y_size_arg, y_kmin_arg, y_kmax_arg, print_it_arg, steady_state_arg, periods_arg, minimal_solving_periods_arg
+                         : dynSparseMatrix(y_size_arg, y_kmin_arg, y_kmax_arg, print_it_arg, steady_state_arg, periods_arg, minimal_solving_periods_arg, slowc_arg
 #ifdef CUDA
                                         , CUDA_device_arg, cublas_handle_arg, cusparse_handle_arg, descr_arg
 #endif