Displays more details in case of singular system

mex/sources/bytecode/Interpreter.cc

@@ -71,7 +71,7 @@ double
 Interpreter::pow1(double a, double b)
 {
   double r = pow_(a, b);
-  if (isnan(r))
+  if (isnan(r)
     {
       res1 = NAN;
       r = 0.0000000000000000000000001;
@@ -1718,6 +1718,16 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
     }
 }
 
+
+void
+Interpreter::SingularDisplay(int Per_u_, bool evaluate, int Block_Count, int size, bool steady_state, it_code_type begining)
+{
+  it_code = begining;
+  compute_block_time(Per_u_, evaluate, Block_Count, size, steady_state);
+  Singular_display(Block_Count, size, steady_state, begining);
+}
+
+
 int
 Interpreter::simulate_a_block(const int size, const int type, string file_name, string bin_basename, bool Gaussian_Elimination, bool steady_state, bool print_it, int block_num,
                               const bool is_linear, const int symbol_table_endo_nbr, const int Block_List_Max_Lag, const int Block_List_Max_Lead, const int u_count_int)
@@ -1726,6 +1736,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
   int i;
   bool cvg;
   bool result = true;
+  bool singular_system;
   double *y_save;
   res1 = 0;
 #ifdef DEBUG
@@ -1979,7 +1990,10 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
                   if (cvg)
                     continue;
                   int prev_iter = iter;
-                  Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, print_it, cvg, iter, true, stack_solve_algo, solve_algo);
+                  singular_system = Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, print_it, cvg, iter, true, stack_solve_algo, solve_algo);
+                  if (singular_system)
+                    SingularDisplay(0, false, block_num, size, steady_state, begining);
+
                   iter++;
                   if (iter > prev_iter)
                     {
@@ -2024,7 +2038,9 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
                 }
               else
                 cvg = false;
-              Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, print_it, cvg, iter, true, stack_solve_algo, solve_algo);
+              singular_system = Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, print_it, cvg, iter, true, stack_solve_algo, solve_algo);
+              if (singular_system)
+                SingularDisplay(0, false, block_num, size, steady_state, begining);
               if (!result)
                 {
                   mexPrintf(" in Solve Forward complete, convergence not achieved in block %d\n", Block_Count+1);
@@ -2076,7 +2092,9 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
                       if (cvg)
                         continue;
                       int prev_iter = iter;
-                      Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, print_it, cvg, iter, false, stack_solve_algo, solve_algo);
+                      singular_system = Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, print_it, cvg, iter, false, stack_solve_algo, solve_algo);
+                      if (singular_system)
+                        SingularDisplay(0, false, block_num, size, steady_state, begining);
                       iter++;
                       if (iter > prev_iter)
                         {
@@ -2123,7 +2141,9 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
                     }
                   else
                     cvg = false;
-                  Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, print_it, cvg, iter, false, stack_solve_algo, solve_algo);
+                  singular_system = Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, print_it, cvg, iter, false, stack_solve_algo, solve_algo);
+                  if (singular_system)
+                    SingularDisplay(0, false, block_num, size, steady_state, begining);
                 }
             }
         }
@@ -2181,7 +2201,9 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
                   if (cvg)
                     continue;
                   int prev_iter = iter;
-                  Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, print_it, cvg, iter, true, stack_solve_algo, solve_algo);
+                  singular_system = Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, print_it, cvg, iter, true, stack_solve_algo, solve_algo);
+                  if (singular_system)
+                    SingularDisplay(0, false, block_num, size, steady_state, begining);
                   iter++;
                   if (iter > prev_iter)
                     {
@@ -2225,7 +2247,9 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
                 }
               else
                 cvg = false;
-              Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, print_it, cvg, iter, true, stack_solve_algo, solve_algo);
+              singular_system = Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, print_it, cvg, iter, true, stack_solve_algo, solve_algo);
+              if (singular_system)
+                SingularDisplay(0, false, block_num, size, steady_state, begining);
               if (!result)
                 {
                   mexPrintf(" in Solve Backward complete, convergence not achieved in block %d\n", Block_Count+1);
@@ -2277,7 +2301,9 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
                       if (cvg)
                         continue;
                       int prev_iter = iter;
-                      Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, print_it, cvg, iter, false, stack_solve_algo, solve_algo);
+                      singular_system = Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, print_it, cvg, iter, false, stack_solve_algo, solve_algo);
+                      if (singular_system)
+                        SingularDisplay(0, false, block_num, size, steady_state, begining);
                       iter++;
                       if (iter > prev_iter)
                         {
@@ -2320,7 +2346,9 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
                     }
                   else
                     cvg = false;
-                  Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, print_it, cvg, iter, false, stack_solve_algo, solve_algo);
+                  singular_system = Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, print_it, cvg, iter, false, stack_solve_algo, solve_algo);
+                  if (singular_system)
+                    SingularDisplay(0, false, block_num, size, steady_state, begining);
                 }
             }
         }

mex/sources/bytecode/Interpreter.hh

@@ -61,6 +61,7 @@ protected:
   void print_a_block(const int size, const int type, string bin_basename, bool steady_state, int block_num,
                      const bool is_linear, const int symbol_table_endo_nbr, const int Block_List_Max_Lag,
                      const int Block_List_Max_Lead, const int u_count_int, int block);
+  void SingularDisplay(int Per_u_, bool evaluate, int Block_Count, int size, bool steady_state, it_code_type begining);
   vector<Block_contain_type> Block_Contain;
   code_liste_type code_liste;
   it_code_type it_code;

mex/sources/bytecode/SparseMatrix.cc

@@ -366,7 +366,7 @@ SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, i
 }
 
 void
-SparseMatrix::Simple_Init(int it_, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM, bool &zero_solution)
+SparseMatrix::Simple_Init(int Size, map<pair<pair<int, int>, int>, int> &IM, bool &zero_solution)
 {
   int i, eq, var, lag;
   map<pair<pair<int, int>, int>, int>::iterator it4;
@@ -2153,6 +2153,94 @@ SparseMatrix::Solve_Matlab_BiCGStab(mxArray *A_m, mxArray *b_m, int Size, double
 }
 
 void
+SparseMatrix::Singular_display(int block, int Size, bool steady_state, it_code_type it_code)
+{
+  bool zero_solution;
+  Simple_Init(Size, IM_i, zero_solution);
+  NonZeroElem *first;
+  mxArray *rhs[1];
+  rhs[0] = mxCreateDoubleMatrix(Size, Size, mxREAL);
+  double *pind;
+  pind = mxGetPr(rhs[0]);
+  for (int j = 0; j < Size * Size; j++)
+    pind[j] = 0.0;
+  for (int ii = 0; ii < Size; ii++)
+    {
+      int nb_eq = At_Col(ii, &first);
+      for (int j = 0; j < nb_eq; j++)
+        {
+          int k = first->u_index;
+          int jj = first->r_index;
+          pind[ii * Size + jj ] = u[k];
+          first = first->NZE_C_N;
+        }
+    }
+  mxArray *lhs[3];
+  mexCallMATLAB(3, lhs, 1, rhs, "svd");
+  mxArray* SVD_u = lhs[0];
+  mxArray* SVD_s = lhs[1];
+  mxArray* SVD_v = lhs[2];
+  double *SVD_ps = mxGetPr(SVD_s);
+  double *SVD_pu = mxGetPr(SVD_u);
+  for (int i = 0; i < Size; i++)
+    {
+      if (abs(SVD_ps[i * (1 + Size)]) < 1e-12)
+        {
+            mexPrintf(" The following equations form a linear combination:\n    ");
+            double max_u = 0;
+            for (int j = 0; j < Size; j++)
+              if (abs(SVD_pu[j + i * Size]) > abs(max_u))
+                max_u = SVD_pu[j + i * Size];
+            vector<int> equ_list;
+            for (int j = 0; j < Size; j++)
+              {
+                double rr = SVD_pu[j + i * Size] / max_u;
+                if ( rr < -1e-10)
+                  {
+                    equ_list.push_back(j);
+                    if (rr != -1)
+                      mexPrintf(" - %3.2f*Dequ_%d_dy",abs(rr),j+1);
+                    else
+                      mexPrintf(" - Dequ_%d_dy",j+1);
+                  }
+                else if (rr > 1e-10)
+                  {
+                    equ_list.push_back(j);
+                    if (j > 0)
+                      if (rr != 1)
+                        mexPrintf(" + %3.2f*Dequ_%d_dy",rr,j+1);
+                      else
+                        mexPrintf(" + Dequ_%d_dy",j+1);
+                    else
+                      if (rr != 1)
+                        mexPrintf(" %3.2f*Dequ_%d_dy",rr,j+1);
+                      else
+                        mexPrintf(" Dequ_%d_dy",j+1);
+                  }
+              }
+            mexPrintf(" = 0\n");
+            /*mexPrintf(" with:\n");
+            it_code = get_begin_block(block);
+            for (int j=0; j < Size; j++)
+              {
+                if (find(equ_list.begin(), equ_list.end(), j) != equ_list.end())
+                  mexPrintf("  equ_%d: %s\n",j, print_expression(it_code_expr, false, Size, block, steady_state, 0, 0, it_code, true).c_str());
+              }*/
+        }
+    }
+  mxDestroyArray(lhs[0]);
+  mxDestroyArray(lhs[1]);
+  mxDestroyArray(lhs[2]);
+  ostringstream tmp;
+  if (block > 1)
+    tmp << " in Solve_ByteCode_Sparse_GaussianElimination, singular system in block " << block+1 << "\n";
+  else
+    tmp << " in Solve_ByteCode_Sparse_GaussianElimination, singular system\n";
+  throw FatalExceptionHandling(tmp.str());
+}
+
+
+bool
 SparseMatrix::Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, bool steady_state, int it_)
 {
   bool one;
@@ -2219,7 +2307,6 @@ SparseMatrix::Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, bool
         }
       if (piv_abs < eps)
         {
-          mexEvalString("drawnow;");
           mxFree(piv_v);
           mxFree(pivj_v);
           mxFree(pivk_v);
@@ -2231,7 +2318,7 @@ SparseMatrix::Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, bool
                 mexPrintf("Error: singular system in Simulate_NG in block %d\n", blck+1);
               else
                 mexPrintf("Error: singular system in Simulate_NG\n");
-              return;
+              return true;
             }
           else
             {
@@ -2411,6 +2498,7 @@ SparseMatrix::Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, bool
   mxFree(pivk_v);
   mxFree(NR);
   mxFree(bc);
+  return false;
 }
 
 void
@@ -2898,13 +2986,14 @@ SparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bool
   End_GE(Size);
 }
 
-void
+bool
 SparseMatrix::Simulate_Newton_One_Boundary(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 stack_solve_algo, int solve_algo)
 {
   int i, j;
   mxArray *b_m = NULL, *A_m = NULL, *x0_m = NULL;
   Clear_u();
   error_not_printed = true;
+  bool singular_system = false;
   u_count_alloc_save = u_count_alloc;
   if (isnan(res1) || isinf(res1) || (res2 > 12*g0 && iter > 0))
     {
@@ -2931,7 +3020,7 @@ SparseMatrix::Simulate_Newton_One_Boundary(int blck, int y_size, int it_, int y_
               else
                 mexPrintf(" dynare cannot improve the simulation in block %d at time %d (variable %d)\n", blck+1, it_+1, index_vara[max_res_idx]+1);
               mexEvalString("drawnow;");
-              return;
+              return singular_system;
             }
           else
             {
@@ -2972,11 +3061,11 @@ SparseMatrix::Simulate_Newton_One_Boundary(int blck, int y_size, int it_, int y_
       for (i = 0; i < y_size; i++)
         y[i+it_*y_size] = ya[i+it_*y_size] + slowc_save*direction[i+it_*y_size];
       iter--;
-      return;
+      return singular_system;
     }
   if (cvg)
     {
-      return;
+      return singular_system;
     }
   if (print_it)
     {
@@ -3024,7 +3113,7 @@ SparseMatrix::Simulate_Newton_One_Boundary(int blck, int y_size, int it_, int y_
     }
   bool zero_solution;
   if ((solve_algo == 5 && steady_state) || (stack_solve_algo == 5 && !steady_state))
-    Simple_Init(it_, y_kmin, y_kmax, Size, IM_i, zero_solution);
+    Simple_Init(Size, IM_i, zero_solution);
   else
     {
       b_m = mxCreateDoubleMatrix(Size, 1, mxREAL);
@@ -3063,7 +3152,7 @@ SparseMatrix::Simulate_Newton_One_Boundary(int blck, int y_size, int it_, int y_
   else
     {
       if ((solve_algo == 5 && steady_state) || (stack_solve_algo == 5 && !steady_state))
-        Solve_ByteCode_Sparse_GaussianElimination(Size, blck, steady_state, it_);
+        singular_system = Solve_ByteCode_Sparse_GaussianElimination(Size, blck, steady_state, it_);
       else if ((solve_algo == 7 && steady_state) || (stack_solve_algo == 2 && !steady_state))
         Solve_Matlab_GMRES(A_m, b_m, Size, slowc, blck, false, it_, steady_state, x0_m);
       else if ((solve_algo == 8 && steady_state) || (stack_solve_algo == 3 && !steady_state))
@@ -3071,7 +3160,7 @@ SparseMatrix::Simulate_Newton_One_Boundary(int blck, int y_size, int it_, int y_
       else if ((solve_algo == 6 && steady_state) || ((stack_solve_algo == 0 || stack_solve_algo == 1) && !steady_state))
         Solve_Matlab_LU_UMFPack(A_m, b_m, Size, slowc, false, it_);
     }
-  return;
+  return singular_system;
 }
 
 void

mex/sources/bytecode/SparseMatrix.hh

@@ -63,21 +63,22 @@ class SparseMatrix : public ErrorMsg
 public:
   SparseMatrix();
   void Simulate_Newton_Two_Boundaries(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 stack_solve_algo, unsigned int endo_name_length, char *P_endo_names) /*throw(ErrorHandlingException)*/;
-  void Simulate_Newton_One_Boundary(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 stack_solve_algo, int solve_algo);
+  bool Simulate_Newton_One_Boundary(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 stack_solve_algo, int solve_algo);
   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, int stack_solve_algo, int solve_algo);
   void Read_file(string file_name, int periods, int u_size1, int y_size, int y_kmin, int y_kmax, int &nb_endo, int &u_count, int &u_count_init, double *u);
+  void Singular_display(int block, int Size, bool steady_state, it_code_type it_code);
   double g0, gp0, glambda2, try_at_iteration;
 
 private:
   void Init_GE(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM);
   void Init_Matlab_Sparse(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM, mxArray *A_m, mxArray *b_m, mxArray *x0_m);
   void Init_Matlab_Sparse_Simple(int Size, map<pair<pair<int, int>, int>, int> &IM, mxArray *A_m, mxArray *b_m, bool &zero_solution, mxArray *x0_m);
-  void Simple_Init(int it_, int y_kmin, int y_kmax, int Size, std::map<std::pair<std::pair<int, int>, int>, int> &IM, bool &zero_solution);
+  void Simple_Init(int Size, std::map<std::pair<std::pair<int, int>, int>, int> &IM, bool &zero_solution);
   void End_GE(int Size);
   void Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bool symbolic, int Block_number);
-  void Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, bool steady_state, int it_);
+  bool Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, bool steady_state, int it_);
   void Solve_Matlab_Relaxation(mxArray *A_m, mxArray *b_m, unsigned int Size, double slowc_l, bool is_two_boundaries, int  it_);
   void Solve_Matlab_LU_UMFPack(mxArray *A_m, mxArray *b_m, int Size, double slowc_l, bool is_two_boundaries, int it_);
   void Solve_Matlab_GMRES(mxArray *A_m, mxArray *b_m, int Size, double slowc, int block, bool is_two_boundaries, int it_, bool steady_state, mxArray *x0_m);