From b1da2f2b7d94482ac8aed9f314cf7d538ae72a78 Mon Sep 17 00:00:00 2001
From: Ferhat Mihoubi <ferhat.mihoubi@univ-evry.fr>
Date: Mon, 10 Jun 2013 16:45:57 +0200
Subject: [PATCH] "slowc" option is now available to solve a simple equation in
 bytecode

---
 mex/sources/bytecode/ErrorHandling.hh |  1 +
 mex/sources/bytecode/Evaluate.cc      | 25 +++++++++++++++++++++++--
 mex/sources/bytecode/Evaluate.hh      |  3 ++-
 mex/sources/bytecode/Interpreter.cc   |  2 +-
 4 files changed, 27 insertions(+), 4 deletions(-)

diff --git a/mex/sources/bytecode/ErrorHandling.hh b/mex/sources/bytecode/ErrorHandling.hh
index b09c62deb..b89fc88f2 100644
--- a/mex/sources/bytecode/ErrorHandling.hh
+++ b/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
 {
diff --git a/mex/sources/bytecode/Evaluate.cc b/mex/sources/bytecode/Evaluate.cc
index 20d392843..68d3172c9 100644
--- a/mex/sources/bytecode/Evaluate.cc
+++ b/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;
diff --git a/mex/sources/bytecode/Evaluate.hh b/mex/sources/bytecode/Evaluate.hh
index 491e81be6..4d95b0428 100644
--- a/mex/sources/bytecode/Evaluate.hh
+++ b/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);
diff --git a/mex/sources/bytecode/Interpreter.cc b/mex/sources/bytecode/Interpreter.cc
index c4b677c82..2e3f8690b 100644
--- a/mex/sources/bytecode/Interpreter.cc
+++ b/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