diff --git a/mex/sources/bytecode/Interpreter.cc b/mex/sources/bytecode/Interpreter.cc index d8dc69987..276a2b72d 100644 --- a/mex/sources/bytecode/Interpreter.cc +++ b/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); } } } diff --git a/mex/sources/bytecode/Interpreter.hh b/mex/sources/bytecode/Interpreter.hh index 8b971e227..22778b14d 100644 --- a/mex/sources/bytecode/Interpreter.hh +++ b/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; code_liste_type code_liste; it_code_type it_code; diff --git a/mex/sources/bytecode/SparseMatrix.cc b/mex/sources/bytecode/SparseMatrix.cc index a7cf37698..db577ac11 100644 --- a/mex/sources/bytecode/SparseMatrix.cc +++ b/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, int>, int> &IM, bool &zero_solution) +SparseMatrix::Simple_Init(int Size, map, int>, int> &IM, bool &zero_solution) { int i, eq, var, lag; map, 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 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 diff --git a/mex/sources/bytecode/SparseMatrix.hh b/mex/sources/bytecode/SparseMatrix.hh index dc7bb3861..e6df80e94 100644 --- a/mex/sources/bytecode/SparseMatrix.hh +++ b/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, int>, int> &IM); void Init_Matlab_Sparse(int periods, int y_kmin, int y_kmax, int Size, map, int>, int> &IM, mxArray *A_m, mxArray *b_m, mxArray *x0_m); void Init_Matlab_Sparse_Simple(int Size, map, 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, int>, int> &IM, bool &zero_solution); + void Simple_Init(int Size, std::map, 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);