diff --git a/mex/sources/block_kalman_filter/block_kalman_filter.cc b/mex/sources/block_kalman_filter/block_kalman_filter.cc index 4608411c7..0f484b04f 100644 --- a/mex/sources/block_kalman_filter/block_kalman_filter.cc +++ b/mex/sources/block_kalman_filter/block_kalman_filter.cc @@ -187,9 +187,9 @@ BlockKalmanFilter::BlockKalmanFilter(int nlhs, mxArray *plhs[], int nrhs, const pY = mxDuplicateArray(prhs[5]); start = mxGetScalar(prhs[6]); /*Defining the initials values*/ - n = mxGetN(pT); // Number of state variables. - pp = mxGetM(pY); // Maximum number of observed variables. - smpl = mxGetN(pY); // Sample size. ; + n = mxGetN(pT); // Number of state variables. + pp = mxGetM(pY); // Maximum number of observed variables. + smpl = mxGetN(pY); // Sample size. ; mfd = mxGetPr(prhs[7]); kalman_tol = mxGetScalar(prhs[8]); riccati_tol = mxGetScalar(prhs[9]); @@ -204,9 +204,9 @@ BlockKalmanFilter::BlockKalmanFilter(int nlhs, mxArray *plhs[], int nrhs, const P = mxGetPr(pP); Y = mxGetPr(pY); - n = mxGetN(pT); // Number of state variables. - pp = mxGetM(pY); // Maximum number of observed variables. - smpl = mxGetN(pY); // Sample size. ; + n = mxGetN(pT); // Number of state variables. + pp = mxGetM(pY); // Maximum number of observed variables. + smpl = mxGetN(pY); // Sample size. ; n_state = n - pure_obs; /*mexPrintf("T\n"); @@ -224,27 +224,27 @@ BlockKalmanFilter::BlockKalmanFilter(int nlhs, mxArray *plhs[], int nrhs, const for (int i = 0; i < n; i++) i_nz_state_var[i] = nz_state_var[i]; - pa = mxCreateDoubleMatrix(n, 1, mxREAL); // State vector. + pa = mxCreateDoubleMatrix(n, 1, mxREAL); // State vector. a = mxGetPr(pa); tmp_a = std::make_unique(n); - dF = 0.0; // det(F). + dF = 0.0; // det(F). p_tmp1 = mxCreateDoubleMatrix(n, n_shocks, mxREAL); tmp1 = mxGetPr(p_tmp1); - t = 0; // Initialization of the time index. - plik = mxCreateDoubleMatrix(smpl, 1, mxREAL); + t = 0; // Initialization of the time index. + plik = mxCreateDoubleMatrix(smpl, 1, mxREAL); lik = mxGetPr(plik); - Inf = mxGetInf(); - LIK = 0.0; // Default value of the log likelihood. - notsteady = true; // Steady state flag. - F_singular = true; + Inf = mxGetInf(); + LIK = 0.0; // Default value of the log likelihood. + notsteady = true; // Steady state flag. + F_singular = true; v_pp = std::make_unique(pp); v_n = std::make_unique(n); mf = std::make_unique(pp); for (int i = 0; i < pp; i++) mf[i] = mfd[i] - 1; - /*compute QQ = R*Q*transpose(R)*/ // Variance of R times the vector of structural innovations.; + /*compute QQ = R*Q*transpose(R)*/ // Variance of R times the vector of structural innovations.; // tmp = R * Q; for (int i = 0; i < n; i++) for (int j = 0; j < n_shocks; j++) @@ -270,9 +270,9 @@ BlockKalmanFilter::BlockKalmanFilter(int nlhs, mxArray *plhs[], int nrhs, const pv = mxCreateDoubleMatrix(pp, 1, mxREAL); v = mxGetPr(pv); - pF = mxCreateDoubleMatrix(pp, pp, mxREAL); + pF = mxCreateDoubleMatrix(pp, pp, mxREAL); F = mxGetPr(pF); - piF = mxCreateDoubleMatrix(pp, pp, mxREAL); + piF = mxCreateDoubleMatrix(pp, pp, mxREAL); iF = mxGetPr(piF); lw = pp * 4; w = std::make_unique(lw); @@ -316,7 +316,7 @@ BlockKalmanFilter::block_kalman_filter_ss() { //v = Y(:,t)-a(mf); for (int i = 0; i < pp; i++) - v[i] = Y[i + t * pp] - a[mf[i]]; + v[i] = Y[i + t * pp] - a[mf[i]]; //a = T*(a+K*v); for (int i = pure_obs; i < n; i++) @@ -375,7 +375,7 @@ BlockKalmanFilter::block_kalman_filter(int nlhs, mxArray *plhs[]) for (auto i = d_index.begin(); i != d_index.end(); i++) { //mexPrintf("i_i=%d, omp_get_max_threads()=%d\n",i_i,omp_get_max_threads()); - v[i_i] = Y[*i + t * pp] - a[mf[*i]]; + v[i_i] = Y[*i + t * pp] - a[mf[*i]]; i_i++; } @@ -404,7 +404,7 @@ BlockKalmanFilter::block_kalman_filter(int nlhs, mxArray *plhs[]) //v = Y(:,t) - a(mf) for (int i = 0; i < pp; i++) - v[i] = Y[i + t * pp] - a[mf[i]]; + v[i] = Y[i + t * pp] - a[mf[i]]; //F = P(mf,mf) + H; if (H_size == 1) @@ -432,10 +432,10 @@ BlockKalmanFilter::block_kalman_filter(int nlhs, mxArray *plhs[]) mexPrintf("dgecon failure with error %d\n", static_cast(info)); if (rcond < kalman_tol) - if (not_all_abs_F_bellow_crit(F, size_d_index * size_d_index, kalman_tol)) //~all(abs(F(:)) > code_liste_type; +typedef vector> code_liste_type; typedef code_liste_type::const_iterator it_code_type; class GeneralExceptionHandling @@ -271,7 +271,7 @@ struct s_plan { string var, exo; int var_num, exo_num; - vector > per_value; + vector> per_value; vector value; }; @@ -310,7 +310,7 @@ public: vector jacobian_block, jacobian_other_endo_block, jacobian_exo_block, jacobian_det_exo_block; map TEF; map, double > TEFD; - map >, double > TEFDD; + map>, double > TEFDD; ExpressionType EQN_type; it_code_type it_code_expr; @@ -319,7 +319,7 @@ public: size_t /*unsigned int*/ endo_name_length, exo_name_length, param_name_length; unsigned int EQN_equation, EQN_block, EQN_block_number; unsigned int EQN_dvar1, EQN_dvar2, EQN_dvar3; - vector > > Variable_list; + vector>> Variable_list; inline ErrorMsg() @@ -955,7 +955,7 @@ public: Stack.pop(); if (compute) { - x[it_+lag+var*nb_row_x] = Stackf.top(); + x[it_+lag+var*nb_row_x] = Stackf.top(); Stackf.pop(); } break; @@ -1036,7 +1036,7 @@ public: Stack.pop(); if (compute) { - x[var] = Stackf.top(); + x[var] = Stackf.top(); Stackf.pop(); } break; @@ -2134,7 +2134,7 @@ public: #ifdef DEBUG double rr = Stackf.top(); mexPrintf("rr=%f\n", rr); - map >, double>::const_iterator it = TEFDD.find(make_pair(indx, make_pair(row-1, col-1))); + map>, double>::const_iterator it = TEFDD.find(make_pair(indx, make_pair(row-1, col-1))); mexPrintf("FSTP TEFDD[make_pair(indx, make_pair(row, col))]=%f done\n", it->second); mexEvalString("drawnow;"); #endif @@ -2157,13 +2157,13 @@ public: #ifdef DEBUG mexPrintf("FLDTEFD\n"); mexPrintf("indx=%d Stack.size()=%d\n", indx, Stack.size()); - map >, double>::const_iterator it = TEFDD.find(make_pair(indx, make_pair(row-1, col-1))); + map>, double>::const_iterator it = TEFDD.find(make_pair(indx, make_pair(row-1, col-1))); mexPrintf("FLD TEFD[make_pair(indx, make_pair(row, col))]=%f done\n", it->second); mexEvalString("drawnow;"); #endif if (compute) { - map >, double>::const_iterator it = TEFDD.find(make_pair(indx, make_pair(row-1, col-1))); + map>, double>::const_iterator it = TEFDD.find(make_pair(indx, make_pair(row-1, col-1))); Stackf.push(it->second); } tmp_out.str(""); diff --git a/mex/sources/bytecode/Evaluate.cc b/mex/sources/bytecode/Evaluate.cc index 2ac731eab..a488a0350 100644 --- a/mex/sources/bytecode/Evaluate.cc +++ b/mex/sources/bytecode/Evaluate.cc @@ -36,7 +36,7 @@ Evaluate::Evaluate() } 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) + print_it(print_it_arg), minimal_solving_periods(minimal_solving_periods_arg) { symbol_table_endo_nbr = 0; Block_List_Max_Lag = 0; @@ -45,7 +45,7 @@ Evaluate::Evaluate(const int y_size_arg, const int y_kmin_arg, const int y_kmax_ block = -1; y_size = y_size_arg; y_kmin = y_kmin_arg; - y_kmax = y_kmax_arg; + y_kmax = y_kmax_arg; periods = periods_arg; steady_state = steady_state_arg; slowc = slowc_arg; @@ -545,7 +545,7 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int case SymbolType::exogenous: var = static_cast(it_code->second)->get_pos(); lag = static_cast(it_code->second)->get_lead_lag(); - x[it_+lag+var*nb_row_x] = Stack.top(); + x[it_+lag+var*nb_row_x] = Stack.top(); #ifdef DEBUG tmp_out << "=>"; mexPrintf(" x[%d, %d](%f)=%s\n", it_+lag, var, x[it_+lag+var*nb_row_x], tmp_out.str().c_str()); @@ -591,7 +591,7 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int case SymbolType::exogenous: case SymbolType::exogenousDet: var = static_cast(it_code->second)->get_pos(); - x[var] = Stack.top(); + x[var] = Stack.top(); #ifdef DEBUG tmp_out << "=>"; mexPrintf(" x[%d, %d](%f)=%s\n", it_+lag, var, x[var], tmp_out.str().c_str()); diff --git a/mex/sources/bytecode/Evaluate.hh b/mex/sources/bytecode/Evaluate.hh index a6aa83da1..84e592411 100644 --- a/mex/sources/bytecode/Evaluate.hh +++ b/mex/sources/bytecode/Evaluate.hh @@ -73,7 +73,7 @@ protected: vector Block_Contain; int size; - int *index_vara; + int *index_vara; bool print_it, forward; int minimal_solving_periods; diff --git a/mex/sources/bytecode/Interpreter.cc b/mex/sources/bytecode/Interpreter.cc index cc11c15af..78e73d179 100644 --- a/mex/sources/bytecode/Interpreter.cc +++ b/mex/sources/bytecode/Interpreter.cc @@ -606,13 +606,13 @@ Interpreter::check_for_controlled_exo_validity(FBEGINBLOCK_ *fb, vector vector endogenous = fb->get_endogenous(); for (vector::iterator it = sconstrained_extended_path.begin(); it != sconstrained_extended_path.end(); it++) { - if ((find(endogenous.begin(), endogenous.end(), it->exo_num) != endogenous.end()) && (find(exogenous.begin(), exogenous.end(), it->var_num) == exogenous.end())) + if ((find(endogenous.begin(), endogenous.end(), it->exo_num) != endogenous.end()) && (find(exogenous.begin(), exogenous.end(), it->var_num) == exogenous.end())) { ostringstream tmp; tmp << "\n the conditional forecast involving as constrained variable " << get_variable(SymbolType::endogenous, it->exo_num) << " and as endogenized exogenous " << get_variable(SymbolType::exogenous, it->var_num) << " that do not appear in block=" << Block_Count+1 << ")\n You should not use block in model options\n"; throw FatalExceptionHandling(tmp.str()); } - else if ((find(endogenous.begin(), endogenous.end(), it->exo_num) != endogenous.end()) && (find(exogenous.begin(), exogenous.end(), it->var_num) != exogenous.end()) && ((fb->get_type() == EVALUATE_FORWARD) || (fb->get_type() != EVALUATE_BACKWARD))) + else if ((find(endogenous.begin(), endogenous.end(), it->exo_num) != endogenous.end()) && (find(exogenous.begin(), exogenous.end(), it->var_num) != exogenous.end()) && ((fb->get_type() == EVALUATE_FORWARD) || (fb->get_type() != EVALUATE_BACKWARD))) { ostringstream tmp; tmp << "\n the conditional forecast cannot be implemented for the block=" << Block_Count+1 << ") that has to be evaluated instead to be solved\n You should not use block in model options\n"; @@ -921,9 +921,9 @@ Interpreter::extended_path(string file_name, string bin_basename, bool evaluate, MainLoop(bin_basename, code, evaluate, block, true, true, sconstrained_extended_path, vector_table_conditional_local); for (int j = 0; j < y_size; j++) { - y_save[j + (t + y_kmin) * y_size] = y[ j + (y_kmin) * y_size]; + y_save[j + (t + y_kmin) * y_size] = y[j + y_kmin * y_size]; if (y_kmin > 0) - y[j ] = y[ j + (y_kmin) * y_size]; + y[j] = y[j + y_kmin * y_size]; } for (int j = 0; j < col_x; j++) { @@ -951,7 +951,7 @@ Interpreter::extended_path(string file_name, string bin_basename, bool evaluate, y_save[j + (k + y_kmin) * y_size] = y[ j + ( k - (nb_periods-1) + y_kmin) * y_size]; }*/ for (int i = 0; i < y_size * col_y; i++) - y[i] = y_save[i]; + y[i] = y_save[i]; for (int j = 0; j < col_x * nb_row_x; j++) x[j] = x_save[j]; if (Init_Code->second) diff --git a/mex/sources/bytecode/Mem_Mngr.cc b/mex/sources/bytecode/Mem_Mngr.cc index fae613bc4..9695b24bd 100644 --- a/mex/sources/bytecode/Mem_Mngr.cc +++ b/mex/sources/bytecode/Mem_Mngr.cc @@ -63,7 +63,7 @@ NonZeroElem * Mem_Mngr::mxMalloc_NZE() { long unsigned int i; - if (!Chunk_Stack.empty()) /*An unused block of memory available inside the heap*/ + if (!Chunk_Stack.empty()) /*An unused block of memory available inside the heap*/ { NonZeroElem *p1 = Chunk_Stack.back(); Chunk_Stack.pop_back(); @@ -74,23 +74,23 @@ Mem_Mngr::mxMalloc_NZE() i = CHUNK_heap_pos++; return (NZE_Mem_add[i]); } - else /*We have to allocate extra memory space*/ + else /*We have to allocate extra memory space*/ { CHUNK_SIZE += CHUNK_BLCK_SIZE; Nb_CHUNK++; - NZE_Mem = static_cast(mxMalloc(CHUNK_BLCK_SIZE*sizeof(NonZeroElem))); /*The block of memory allocated*/ + NZE_Mem = static_cast(mxMalloc(CHUNK_BLCK_SIZE*sizeof(NonZeroElem))); /*The block of memory allocated*/ error_msg.test_mxMalloc(NZE_Mem, __LINE__, __FILE__, __func__, CHUNK_BLCK_SIZE*sizeof(NonZeroElem)); NZE_Mem_Allocated.push_back(NZE_Mem); if (!NZE_Mem) mexPrintf("Not enough memory available\n"); if (NZE_Mem_add) { - NZE_Mem_add = static_cast(mxRealloc(NZE_Mem_add, CHUNK_SIZE*sizeof(NonZeroElem *))); /*We have to redefine the size of pointer on the memory*/ + NZE_Mem_add = static_cast(mxRealloc(NZE_Mem_add, CHUNK_SIZE*sizeof(NonZeroElem *))); /*We have to redefine the size of pointer on the memory*/ error_msg.test_mxMalloc(NZE_Mem_add, __LINE__, __FILE__, __func__, CHUNK_SIZE*sizeof(NonZeroElem *)); } else { - NZE_Mem_add = static_cast(mxMalloc(CHUNK_SIZE*sizeof(NonZeroElem *))); /*We have to define the size of pointer on the memory*/ + NZE_Mem_add = static_cast(mxMalloc(CHUNK_SIZE*sizeof(NonZeroElem *))); /*We have to define the size of pointer on the memory*/ error_msg.test_mxMalloc(NZE_Mem_add, __LINE__, __FILE__, __func__, CHUNK_SIZE*sizeof(NonZeroElem *)); } diff --git a/mex/sources/bytecode/Mem_Mngr.hh b/mex/sources/bytecode/Mem_Mngr.hh index 48c65b2d3..b754157fe 100644 --- a/mex/sources/bytecode/Mem_Mngr.hh +++ b/mex/sources/bytecode/Mem_Mngr.hh @@ -60,7 +60,7 @@ private: NonZeroElem *NZE_Mem; vector NZE_Mem_Allocated; int swp_f_b; - fstream SaveCode_swp; + fstream SaveCode_swp; string filename_mem; }; diff --git a/mex/sources/bytecode/SparseMatrix.cc b/mex/sources/bytecode/SparseMatrix.cc index 322b6943f..4af1ffa20 100644 --- a/mex/sources/bytecode/SparseMatrix.cc +++ b/mex/sources/bytecode/SparseMatrix.cc @@ -39,49 +39,49 @@ HINSTANCE hinstLib; # define UMFPACK_INFO 90 # define UMFPACK_CONTROL 20 /* used in all UMFPACK_report_* routines: */ -# define UMFPACK_PRL 0 /* print level */ +# define UMFPACK_PRL 0 /* print level */ /* returned by all routines that use Info: */ # define UMFPACK_OK (0) -# define UMFPACK_STATUS 0 /* UMFPACK_OK, or other result */ +# define UMFPACK_STATUS 0 /* UMFPACK_OK, or other result */ typedef void (*t_umfpack_dl_free_numeric)(void **Numeric); t_umfpack_dl_free_numeric umfpack_dl_free_numeric; typedef void (*t_umfpack_dl_free_symbolic)(void **Symbolic); t_umfpack_dl_free_symbolic umfpack_dl_free_symbolic; typedef int64_t (*t_umfpack_dl_solve)(int64_t sys, - const int64_t Ap [], - const int64_t Ai [], - const double Ax [], - double X [], - const double B [], + const int64_t Ap[], + const int64_t Ai[], + const double Ax[], + double X[], + const double B[], void *Numeric, - const double Control [UMFPACK_CONTROL], - double Info [UMFPACK_INFO]); + const double Control[UMFPACK_CONTROL], + double Info[UMFPACK_INFO]); t_umfpack_dl_solve umfpack_dl_solve; -typedef int64_t (*t_umfpack_dl_numeric)(const int64_t Ap [], - const int64_t Ai [], - const double Ax [], +typedef int64_t (*t_umfpack_dl_numeric)(const int64_t Ap[], + const int64_t Ai[], + const double Ax[], void *Symbolic, void **Numeric, - const double Control [UMFPACK_CONTROL], - double Info [UMFPACK_INFO]); + const double Control[UMFPACK_CONTROL], + double Info[UMFPACK_INFO]); t_umfpack_dl_numeric umfpack_dl_numeric; typedef int64_t (*t_umfpack_dl_symbolic)(int64_t n_row, int64_t n_col, - const int64_t Ap [], - const int64_t Ai [], - const double Ax [], + const int64_t Ap[], + const int64_t Ai[], + const double Ax[], void **Symbolic, - const double Control [UMFPACK_CONTROL], - double Info [UMFPACK_INFO]); + const double Control[UMFPACK_CONTROL], + double Info[UMFPACK_INFO]); t_umfpack_dl_symbolic umfpack_dl_symbolic; -typedef void (*t_umfpack_dl_report_info)(const double Control [UMFPACK_CONTROL], - const double Info [UMFPACK_INFO]); +typedef void (*t_umfpack_dl_report_info)(const double Control[UMFPACK_CONTROL], + const double Info[UMFPACK_INFO]); t_umfpack_dl_report_info umfpack_dl_report_info; -typedef void (*t_umfpack_dl_report_status)(const double Control [UMFPACK_CONTROL], +typedef void (*t_umfpack_dl_report_status)(const double Control[UMFPACK_CONTROL], int64_t status); t_umfpack_dl_report_status umfpack_dl_report_status; -typedef void (*t_umfpack_dl_defaults)(double Control [UMFPACK_CONTROL]); +typedef void (*t_umfpack_dl_defaults)(double Control[UMFPACK_CONTROL]); t_umfpack_dl_defaults umfpack_dl_defaults; #endif @@ -688,7 +688,7 @@ dynSparseMatrix::Simple_Init(int Size, map, int>, int> &IM, var = it4->first.first.second; eq = it4->first.first.first; lag = it4->first.second; - if (lag == 0) /*Build the index for sparse matrix containing the jacobian : u*/ + if (lag == 0) /*Build the index for sparse matrix containing the jacobian : u*/ { NbNZRow[eq]++; NbNZCol[var]++; @@ -796,7 +796,7 @@ dynSparseMatrix::Init_Matlab_Sparse_Simple(int Size, map, in var = it4->first.first.first; if (var != last_var) { - Aj[1+last_var ] = NZE; + Aj[1+last_var] = NZE; last_var = var; } eq = it4->first.second; @@ -894,7 +894,7 @@ dynSparseMatrix::Init_UMFPACK_Sparse_Simple(int Size, map, i ya[eq+it_*y_size] = y[eq+it_*y_size]; } #ifdef DEBUG - unsigned int max_nze = prior_nz;//mxGetNzmax(A_m); + unsigned int max_nze = prior_nz; //mxGetNzmax(A_m); #endif unsigned int NZE = 0; int last_var = 0; @@ -918,7 +918,7 @@ dynSparseMatrix::Init_UMFPACK_Sparse_Simple(int Size, map, i var = it4->first.first.first; if (var != last_var) { - (*Ap)[1+last_var ] = NZE; + (*Ap)[1+last_var] = NZE; last_var = var; } eq = it4->first.second; @@ -980,11 +980,11 @@ dynSparseMatrix::find_exo_num(vector sconstrained_extended_path, int val } int -dynSparseMatrix::find_int_date(vector > per_value, int value) +dynSparseMatrix::find_int_date(vector> per_value, int value) { int res = -1; int i = 0; - for (vector >::iterator it = per_value.begin(); it != per_value.end(); it++, i++) + for (vector>::iterator it = per_value.begin(); it != per_value.end(); it++, i++) if (it->first == value) { res = i; @@ -1158,7 +1158,7 @@ dynSparseMatrix::Init_UMFPACK_Sparse(int periods, int y_kmin, int y_kmax, int Si throw FatalExceptionHandling(tmp.str()); } #endif - u[k] -= jacob_exo[k + row_x*flip_exo] * x[t+y_kmin+flip_exo*nb_row_x]; + u[k] -= jacob_exo[k + row_x*flip_exo] * x[t+y_kmin+flip_exo*nb_row_x]; } } } @@ -1178,7 +1178,7 @@ dynSparseMatrix::Init_UMFPACK_Sparse(int periods, int y_kmin, int y_kmax, int Si ti_y_kmax = min(periods-(t +1), y_kmax); int ti_new_y_kmax = min(t, y_kmax); int ti_new_y_kmin = -min(periods-(t+1), y_kmin); - if (lag <= ti_new_y_kmax && lag >= ti_new_y_kmin) /*Build the index for sparse matrix containing the jacobian : u*/ + if (lag <= ti_new_y_kmax && lag >= ti_new_y_kmin) /*Build the index for sparse matrix containing the jacobian : u*/ { #ifdef DEBUG if (index < 0 || index >= u_count_alloc || index > Size + Size*Size) @@ -1248,10 +1248,10 @@ dynSparseMatrix::Init_UMFPACK_Sparse(int periods, int y_kmin, int y_kmax, int Si throw FatalExceptionHandling(tmp.str()); } #endif - (*b)[eq] += u[index+lag*u_count_init]*y[index_vara[var+Size*(y_kmin+t+lag)]]; + (*b)[eq] += u[index+lag*u_count_init]*y[index_vara[var+Size*(y_kmin+t+lag)]]; } } - else /* ...and store it in the u vector*/ + else /* ...and store it in the u vector*/ { #ifdef DEBUG if (index < 0 || index >= u_count_alloc) @@ -1267,7 +1267,7 @@ dynSparseMatrix::Init_UMFPACK_Sparse(int periods, int y_kmin, int y_kmax, int Si throw FatalExceptionHandling(tmp.str()); } #endif - (*b)[eq] += u[index]; + (*b)[eq] += u[index]; } it4++; } @@ -1369,7 +1369,7 @@ dynSparseMatrix::Init_CUDA_Sparse_Simple(int Size, map, int> var = it4->first.first.first; if (var != last_var) { - (*Ap)[1+last_var ] = NZE; + (*Ap)[1+last_var] = NZE; last_var = var; } eq = it4->first.second; @@ -1509,7 +1509,7 @@ dynSparseMatrix::Init_CUDA_Sparse(int periods, int y_kmin, int y_kmax, int Size, # endif Host_Ap[1+last_eq + t * Size] = NZE; if (preconditioner == 3 && t == 0) - Host_Ap_tild[1+last_eq ] = NZE_tild; + Host_Ap_tild[1+last_eq] = NZE_tild; last_eq = eq; } var = it4->first.second+Size*t; @@ -1519,7 +1519,7 @@ dynSparseMatrix::Init_CUDA_Sparse(int periods, int y_kmin, int y_kmax, int Size, { ti_y_kmin = -min(t, y_kmin); ti_y_kmax = min(periods-(t + 1), y_kmax); - if ((lag <= ti_y_kmax && lag >= ti_y_kmin) || preconditioner == 3) /*Build the index for sparse matrix containing the jacobian : u*/ + if ((lag <= ti_y_kmax && lag >= ti_y_kmin) || preconditioner == 3) /*Build the index for sparse matrix containing the jacobian : u*/ { # ifdef DEBUG if (index < 0 || index >= u_count_alloc || index > (periods-1)* IM.size() + Size * Size + periods * Size) @@ -1547,7 +1547,7 @@ dynSparseMatrix::Init_CUDA_Sparse(int periods, int y_kmin, int y_kmax, int Size, { if (lag > ti_y_kmax || lag < ti_y_kmin) { - Host_b[eq + t * Size] += u[index]*y[index_vara[var+Size*(y_kmin+lag)]]; + Host_b[eq + t * Size] += u[index]*y[index_vara[var+Size*(y_kmin+lag)]]; to_store = false; } if (t == 0) @@ -1593,10 +1593,10 @@ dynSparseMatrix::Init_CUDA_Sparse(int periods, int y_kmin, int y_kmax, int Size, throw FatalExceptionHandling(tmp.str()); } # endif - Host_b[eq + t * Size] += u[index]*y[index_vara[var+Size*(y_kmin+lag)]]; + Host_b[eq + t * Size] += u[index]*y[index_vara[var+Size*(y_kmin+lag)]]; } } - else // ...and store it in the u vector + else // ...and store it in the u vector { # ifdef DEBUG if (index < 0 || index >= u_count_alloc) @@ -1612,7 +1612,7 @@ dynSparseMatrix::Init_CUDA_Sparse(int periods, int y_kmin, int y_kmax, int Size, throw FatalExceptionHandling(tmp.str()); } # endif - Host_b[var] += u[index]; + Host_b[var] += u[index]; } it4++; } @@ -1680,15 +1680,15 @@ dynSparseMatrix::Init_CUDA_Sparse(int periods, int y_kmin, int y_kmax, int Size, } cudaChk(cudaMalloc((void **) A_tild, preconditioner_size * sizeof(double)), " in Init_Cuda_Sparse, can't allocate A_tild on the graphic card\n"); - cudaChk(cudaMemcpy(*x0, Host_x0, n * sizeof(double), cudaMemcpyHostToDevice), " in Init_CUDA_Sparse, cudaMemcpy x0 = Host_x0 failed"); - cudaChk(cudaMemcpy(*b, Host_b, n * sizeof(double), cudaMemcpyHostToDevice), " in Init_CUDA_Sparse, cudaMemcpy b = Host_b failed"); - cudaChk(cudaMemcpy(*Ap, Host_Ap, (n + 1) * sizeof(int), cudaMemcpyHostToDevice), " in Init_CUDA_Sparse, cudaMemcpy Ap = Host_Ap failed"); - cudaChk(cudaMemcpy(*Ai, Host_Ai, NZE * sizeof(int), cudaMemcpyHostToDevice), " in Init_CUDA_Sparse, cudaMemcpy Ai = Host_Ai failed"); - cudaChk(cudaMemcpy(*Ax, Host_Ax, NZE * sizeof(double), cudaMemcpyHostToDevice), " in Init_CUDA_Sparse, cudaMemcpy Ax = Host_Ax failed"); + cudaChk(cudaMemcpy(*x0, Host_x0, n * sizeof(double), cudaMemcpyHostToDevice), " in Init_CUDA_Sparse, cudaMemcpy x0 = Host_x0 failed"); + cudaChk(cudaMemcpy(*b, Host_b, n * sizeof(double), cudaMemcpyHostToDevice), " in Init_CUDA_Sparse, cudaMemcpy b = Host_b failed"); + cudaChk(cudaMemcpy(*Ap, Host_Ap, (n + 1) * sizeof(int), cudaMemcpyHostToDevice), " in Init_CUDA_Sparse, cudaMemcpy Ap = Host_Ap failed"); + cudaChk(cudaMemcpy(*Ai, Host_Ai, NZE * sizeof(int), cudaMemcpyHostToDevice), " in Init_CUDA_Sparse, cudaMemcpy Ai = Host_Ai failed"); + cudaChk(cudaMemcpy(*Ax, Host_Ax, NZE * sizeof(double), cudaMemcpyHostToDevice), " in Init_CUDA_Sparse, cudaMemcpy Ax = Host_Ax failed"); if (preconditioner == 3) { - cudaChk(cudaMemcpy(*Ap_tild, Host_Ap_tild, (n + 1) * sizeof(int), cudaMemcpyHostToDevice), " in Init_CUDA_Sparse, cudaMemcpy Ap_tild = Host_Ap_tild failed"); - cudaChk(cudaMemcpy(*Ai_tild, Host_Ai_tild, NZE_tild * sizeof(int), cudaMemcpyHostToDevice), " in Init_CUDA_Sparse, cudaMemcpy Ai_tild = Host_Ai_til failed"); + cudaChk(cudaMemcpy(*Ap_tild, Host_Ap_tild, (n + 1) * sizeof(int), cudaMemcpyHostToDevice), " in Init_CUDA_Sparse, cudaMemcpy Ap_tild = Host_Ap_tild failed"); + cudaChk(cudaMemcpy(*Ai_tild, Host_Ai_tild, NZE_tild * sizeof(int), cudaMemcpyHostToDevice), " in Init_CUDA_Sparse, cudaMemcpy Ai_tild = Host_Ai_til failed"); } cudaChk(cudaMemcpy(*A_tild, Host_A_tild, preconditioner_size * sizeof(double), cudaMemcpyHostToDevice), " in Init_CUDA_Sparse, cudaMemcpy A_tild = Host_A_tild failed"); } @@ -1800,7 +1800,7 @@ dynSparseMatrix::Init_Matlab_Sparse(int periods, int y_kmin, int y_kmax, int Siz ti_y_kmax = min(periods-(t +1), y_kmax); int ti_new_y_kmax = min(t, y_kmax); int ti_new_y_kmin = -min(periods-(t+1), y_kmin); - if (lag <= ti_new_y_kmax && lag >= ti_new_y_kmin) /*Build the index for sparse matrix containing the jacobian : u*/ + if (lag <= ti_new_y_kmax && lag >= ti_new_y_kmin) /*Build the index for sparse matrix containing the jacobian : u*/ { #ifdef DEBUG if (index < 0 || index >= u_count_alloc || index > Size + Size*Size) @@ -1842,10 +1842,10 @@ dynSparseMatrix::Init_Matlab_Sparse(int periods, int y_kmin, int y_kmax, int Siz throw FatalExceptionHandling(tmp.str()); } #endif - b[eq] += u[index+lag*u_count_init]*y[index_vara[var+Size*(y_kmin+t+lag)]]; + b[eq] += u[index+lag*u_count_init]*y[index_vara[var+Size*(y_kmin+t+lag)]]; } } - else /* ...and store it in the u vector*/ + else /* ...and store it in the u vector*/ { #ifdef DEBUG if (index < 0 || index >= u_count_alloc) @@ -1861,7 +1861,7 @@ dynSparseMatrix::Init_Matlab_Sparse(int periods, int y_kmin, int y_kmax, int Siz throw FatalExceptionHandling(tmp.str()); } #endif - b[eq] += u[index]; + b[eq] += u[index]; } it4++; } @@ -1941,7 +1941,7 @@ dynSparseMatrix::Init_GE(int periods, int y_kmin, int y_kmax, int Size, mapfirst.second; - if (lag <= ti_y_kmax && lag >= ti_y_kmin) /*Build the index for sparse matrix containing the jacobian : u*/ + if (lag <= ti_y_kmax && lag >= ti_y_kmin) /*Build the index for sparse matrix containing the jacobian : u*/ { nnz++; var += Size*t; @@ -1965,7 +1965,7 @@ dynSparseMatrix::Init_GE(int periods, int y_kmin, int y_kmax, int Size, mapsecond+u_count_init*t; u[b[eq]] += tmp_b; @@ -2062,7 +2062,7 @@ dynSparseMatrix::End_GE(int Size) } bool -dynSparseMatrix::compare(int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, long int nop4, int Size) +dynSparseMatrix::compare(int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, long int nop4, int Size) { long int i, j, nop = nop4/2; double r = 0.0; @@ -2981,7 +2981,7 @@ dynSparseMatrix::golden(double ax, double bx, double cx, double tol, double solv } void -dynSparseMatrix::Solve_Matlab_Relaxation(mxArray *A_m, mxArray *b_m, unsigned int Size, double slowc_l, bool is_two_boundaries, int it_) +dynSparseMatrix::Solve_Matlab_Relaxation(mxArray *A_m, mxArray *b_m, unsigned int Size, double slowc_l, bool is_two_boundaries, int it_) { mxArray *B1, *C1, *A2, *B2, *A3, *b1, *b2; double *b_m_d = mxGetPr(b_m); @@ -3133,7 +3133,7 @@ dynSparseMatrix::Solve_Matlab_Relaxation(mxArray *A_m, mxArray *b_m, unsigned in while (A3_var < Size) A3_j[++A3_var] = A3_nze; mxArray *d1 = NULL; - vector > triangular_form; + vector> triangular_form; double sumc = 0, C_sumc = 1000; mxArray *B1_inv = NULL; mxArray *B1_inv_t = NULL; @@ -3255,7 +3255,7 @@ dynSparseMatrix::Solve_Matlab_Relaxation(mxArray *A_m, mxArray *b_m, unsigned in } void -dynSparseMatrix::Solve_Matlab_LU_UMFPack(mxArray *A_m, mxArray *b_m, int Size, double slowc_l, bool is_two_boundaries, int it_) +dynSparseMatrix::Solve_Matlab_LU_UMFPack(mxArray *A_m, mxArray *b_m, int Size, double slowc_l, bool is_two_boundaries, int it_) { size_t n = mxGetM(A_m); mxArray *z; @@ -3305,12 +3305,12 @@ void dynSparseMatrix::Printfull_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, double *Ax, double *b, int n) { double A[n*n]; - for (int i = 0; i < n*n; i++) + for (int i = 0; i < n*n; i++) A[i] = 0; int k = 0; for (int i = 0; i < n; i++) for (int j = Ap[i]; j < Ap[i+1]; j++) - A[Ai[j] * n + i] = Ax[k++]; + A[Ai[j] * n + i] = Ax[k++]; for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) @@ -3329,11 +3329,11 @@ dynSparseMatrix::Print_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, doubl } void -dynSparseMatrix::Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, double *Ax, double *b, int n, int Size, double slowc_l, bool is_two_boundaries, int it_, vector_table_conditional_local_type vector_table_conditional_local) +dynSparseMatrix::Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, double *Ax, double *b, int n, int Size, double slowc_l, bool is_two_boundaries, int it_, vector_table_conditional_local_type vector_table_conditional_local) { SuiteSparse_long status, sys = 0; #ifndef _MSC_VER - double Control [UMFPACK_CONTROL], Info [UMFPACK_INFO], res [n]; + double Control[UMFPACK_CONTROL], Info[UMFPACK_INFO], res[n]; #else double *Control, *Info, *res; Control = (double *) mxMalloc(UMFPACK_CONTROL * sizeof(double)); @@ -3345,7 +3345,7 @@ dynSparseMatrix::Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, do #endif umfpack_dl_defaults(Control); - Control [UMFPACK_PRL] = 5; + Control[UMFPACK_PRL] = 5; status = 0; if (iter == 0) { @@ -3354,7 +3354,7 @@ dynSparseMatrix::Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, do { umfpack_dl_report_info(Control, Info); umfpack_dl_report_status(Control, status); - ostringstream Error; + ostringstream Error; Error << " umfpack_dl_symbolic failed\n"; throw FatalExceptionHandling(Error.str()); } @@ -3366,7 +3366,7 @@ dynSparseMatrix::Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, do { umfpack_dl_report_info(Control, Info); umfpack_dl_report_status(Control, status); - ostringstream Error; + ostringstream Error; Error << " umfpack_dl_numeric failed\n"; throw FatalExceptionHandling(Error.str()); } @@ -3375,7 +3375,7 @@ dynSparseMatrix::Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, do { umfpack_dl_report_info(Control, Info); umfpack_dl_report_status(Control, status); - ostringstream Error; + ostringstream Error; Error << " umfpack_dl_solve failed\n"; throw FatalExceptionHandling(Error.str()); } @@ -3393,14 +3393,14 @@ dynSparseMatrix::Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, do { int eq = index_vara[i+Size*(y_kmin)]; int flip_exo = vector_table_conditional_local[i].var_exo; - double yy = -(res[i] + x[y_kmin + flip_exo*nb_row_x]); + double yy = -(res[i] + x[y_kmin + flip_exo*nb_row_x]); direction[eq] = 0; x[flip_exo*nb_row_x + y_kmin] += slowc_l * yy; } else { int eq = index_vara[i+Size*(y_kmin)]; - double yy = -(res[i ] + y[eq]); + double yy = -(res[i] + y[eq]); direction[eq] = yy; y[eq] += slowc_l * yy; } @@ -3457,11 +3457,11 @@ dynSparseMatrix::Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, do } void -dynSparseMatrix::Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, double *Ax, double *b, int n, int Size, double slowc_l, bool is_two_boundaries, int it_) +dynSparseMatrix::Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, double *Ax, double *b, int n, int Size, double slowc_l, bool is_two_boundaries, int it_) { SuiteSparse_long status, sys = 0; #ifndef _MSC_VER - double Control [UMFPACK_CONTROL], Info [UMFPACK_INFO], res [n]; + double Control[UMFPACK_CONTROL], Info[UMFPACK_INFO], res[n]; #else double *Control, *Info, *res; Control = (double *) mxMalloc(UMFPACK_CONTROL * sizeof(double)); @@ -3473,7 +3473,7 @@ dynSparseMatrix::Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, do #endif umfpack_dl_defaults(Control); - Control [UMFPACK_PRL] = 5; + Control[UMFPACK_PRL] = 5; status = 0; if (iter == 0) { @@ -3482,7 +3482,7 @@ dynSparseMatrix::Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, do { umfpack_dl_report_info(Control, Info); umfpack_dl_report_status(Control, status); - ostringstream Error; + ostringstream Error; Error << " umfpack_dl_symbolic failed\n"; throw FatalExceptionHandling(Error.str()); } @@ -3494,7 +3494,7 @@ dynSparseMatrix::Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, do { umfpack_dl_report_info(Control, Info); umfpack_dl_report_status(Control, status); - ostringstream Error; + ostringstream Error; Error << " umfpack_dl_numeric failed\n"; throw FatalExceptionHandling(Error.str()); } @@ -3503,7 +3503,7 @@ dynSparseMatrix::Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, do { umfpack_dl_report_info(Control, Info); umfpack_dl_report_status(Control, status); - ostringstream Error; + ostringstream Error; Error << " umfpack_dl_solve failed\n"; throw FatalExceptionHandling(Error.str()); } @@ -3536,7 +3536,7 @@ dynSparseMatrix::Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, do } void -dynSparseMatrix::Solve_LU_UMFPack(mxArray *A_m, mxArray *b_m, int Size, double slowc_l, bool is_two_boundaries, int it_) +dynSparseMatrix::Solve_LU_UMFPack(mxArray *A_m, mxArray *b_m, int Size, double slowc_l, bool is_two_boundaries, int it_) { SuiteSparse_long n = mxGetM(A_m); @@ -3544,10 +3544,10 @@ dynSparseMatrix::Solve_LU_UMFPack(mxArray *A_m, mxArray *b_m, int Size, double s SuiteSparse_long *Ai = reinterpret_cast(mxGetIr(A_m)); double *Ax = mxGetPr(A_m); - double *B = mxGetPr(b_m); + double *B = mxGetPr(b_m); SuiteSparse_long status, sys = 0; #ifndef _MSC_VER - double Control [UMFPACK_CONTROL], Info [UMFPACK_INFO], res [n]; + double Control[UMFPACK_CONTROL], Info[UMFPACK_INFO], res[n]; #else double *Control, *Info, *res; Control = (double *) mxMalloc(UMFPACK_CONTROL * sizeof(double)); @@ -3600,7 +3600,7 @@ dynSparseMatrix::Solve_LU_UMFPack(mxArray *A_m, mxArray *b_m, int Size, double s #ifdef CUDA void -printM(int n, double *Ax, int *Ap, int *Ai, cusparseMatDescr_t descrA, cusparseHandle_t cusparse_handle) +printM(int n, double *Ax, int *Ap, int *Ai, cusparseMatDescr_t descrA, cusparseHandle_t cusparse_handle) { //cudaError_t cuda_error; //cusparseStatus_t cusparse_status; @@ -3627,8 +3627,8 @@ printM(int n, double *Ax, int *Ap, int *Ai, cusparseMatDescr_t descrA, cusparse void dynSparseMatrix::Solve_CUDA_BiCGStab_Free(double *tmp_vect_host, double *p, double *r, double *v, double *s, double *t, double *y_, double *z, double *tmp_, - int *Ai, double *Ax, int *Ap, double *x0, double *b, double *A_tild, int *A_tild_i, int *A_tild_p/*, double* Lx, int* Li, int* Lp, - double* Ux, int* Ui, int* Up, int* device_n*/, cusparseSolveAnalysisInfo_t infoL, cusparseSolveAnalysisInfo_t infoU, + int *Ai, double *Ax, int *Ap, double *x0, double *b, double *A_tild, int *A_tild_i, int *A_tild_p /*, double* Lx, int* Li, int* Lp, + double* Ux, int* Ui, int* Up, int* device_n*/, cusparseSolveAnalysisInfo_t infoL, cusparseSolveAnalysisInfo_t infoU, cusparseMatDescr_t descrL, cusparseMatDescr_t descrU, int preconditioner) { //cudaError_t cuda_error; @@ -3726,7 +3726,7 @@ Check(int n, double *Ax, int *Ap, int *Ai, double *b, double *x, bool Lower) if (k < i) sum += x[k] * Ax[j]; } - double err = b[i] - sum - x[i]; + double err = b[i] - sum - x[i]; if (abs(err) > 1e-10) mexPrintf("error at i=%d\n", i); } @@ -3742,7 +3742,7 @@ Check(int n, double *Ax, int *Ap, int *Ai, double *b, double *x, bool Lower) if (k >= i) sum += x[k] * Ax[j]; } - double err = b[i] - sum; + double err = b[i] - sum; if (abs(err) > 1e-10) mexPrintf("error at i=%d\n", i); } @@ -3752,11 +3752,11 @@ Check(int n, double *Ax, int *Ap, int *Ai, double *b, double *x, bool Lower) #ifdef CUDA int dynSparseMatrix::Solve_CUDA_BiCGStab(int *Ap, int *Ai, double *Ax, int *Ap_tild, int *Ai_tild, double *A_tild, double *b, double *x0, int n, int Size, double slowc_l, bool is_two_boundaries, - int it_, int nnz, int nnz_tild, int preconditioner, int max_iterations, int block) + int it_, int nnz, int nnz_tild, int preconditioner, int max_iterations, int block) { cusparseSolveAnalysisInfo_t info, infoL, infoU; cusparseMatDescr_t descrL, descrU; - const double tol = 1.0e-6;//1.0e-6; + const double tol = 1.0e-6; //1.0e-6; const double eps = 1.0e-16; double *p, *r, *r0, *v, *s, *t, *y_, *z, *tmp_; int *A_tild_i, *A_tild_p; @@ -4277,14 +4277,14 @@ dynSparseMatrix::Solve_CUDA_BiCGStab(int *Ap, int *Ai, double *Ax, int *Ap_tild, for (int t = 0; t < periods; t++) for (int i = 0; i < Size; i++) { - for (mwIndex l = Ujw_host[i]; l < Ujw_host[i+1]; l++) + for (mwIndex l = Ujw_host[i]; l < Ujw_host[i+1]; l++) { Wi[W_nnz] = Uiw_host[l] + t * Size; pW[W_nnz] = Ux_host[l]; //mexPrintf("Wj[%d] = %d, Wi[%d] = Uiw_host[%d] + t * Size = %d, pW[%d]=%f\n", i + t * Size, Wj[i + t * Size], W_nnz, l, Uiw_host[l] + t * Size, W_nnz, Ux_host[l]); W_nnz++; } - for (mwIndex l = Ljw_host[i]; l < Ljw_host[i+1]; l++) + for (mwIndex l = Ljw_host[i]; l < Ljw_host[i+1]; l++) { if (Liw_host[l] > i) { @@ -4405,7 +4405,7 @@ dynSparseMatrix::Solve_CUDA_BiCGStab(int *Ap, int *Ai, double *Ax, int *Ap_tild, " in Solve_Cuda_BiCGStab, cusparseSetMatDiagType has failed for descrU\n"); int host_nnz_tild; - if (preconditioner == 3) + if (preconditioner == 3) host_nnz_tild = W_nnz; else host_nnz_tild = nnz; @@ -4601,7 +4601,7 @@ dynSparseMatrix::Solve_CUDA_BiCGStab(int *Ap, int *Ai, double *Ax, int *Ap_tild, &tmp2), " in Solve_Cuda_BiCGStab, cublasDnrm2(y_) has failed\n"); mexPrintf("abs(alpha)*tmp1 = %f, alpha = %f, tmp1 = %f, tmp2 = %f, eps = %f\n", abs(alpha)*tmp1, alpha, tmp1, tmp2, eps); - if (abs(alpha)*tmp1 < eps * tmp2) + if (abs(alpha)*tmp1 < eps * tmp2) stagnation++; else stagnation = 0; @@ -4979,7 +4979,7 @@ dynSparseMatrix::Solve_Matlab_BiCGStab(mxArray *A_m, mxArray *b_m, int Size, dou double flags = 2; mxArray *z; z = NULL; - if (steady_state) /*Octave BicStab algorihtm involves a 0 division in case of a preconditionner equal to the LU decomposition of A matrix*/ + if (steady_state) /*Octave BicStab algorihtm involves a 0 division in case of a preconditionner equal to the LU decomposition of A matrix*/ { mxArray *res = mult_SAT_B(Sparse_transpose(A_m), x0_m); double *resid = mxGetPr(res); @@ -5129,7 +5129,7 @@ dynSparseMatrix::Singular_display(int block, int Size) { int k = first->u_index; int jj = first->r_index; - pind[ii * Size + jj ] = u[k]; + pind[ii * Size + jj] = u[k]; first = first->NZE_C_N; } } @@ -5314,8 +5314,8 @@ dynSparseMatrix::Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, i if (markovitz > markovitz_max) { piv = piv_v[j]; - pivj = pivj_v[j]; //Line number - pivk = pivk_v[j]; //positi + pivj = pivj_v[j]; //Line number + pivk = pivk_v[j]; //positi markovitz_max = markovitz; } } @@ -5341,8 +5341,8 @@ dynSparseMatrix::Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, i if (NR[j] == 1) { piv = piv_v[j]; - pivj = pivj_v[j]; //Line number - pivk = pivk_v[j]; //positi + pivj = pivj_v[j]; //Line number + pivk = pivk_v[j]; //positi markovitz_max = markovitz; } } @@ -5588,8 +5588,8 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo if (markovitz > markovitz_max) { piv = piv_v[j]; - pivj = pivj_v[j]; //Line number - pivk = pivk_v[j]; //positi + pivj = pivj_v[j]; //Line number + pivk = pivk_v[j]; //positi markovitz_max = markovitz; NR_max = NR[j]; } @@ -5616,8 +5616,8 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo if (NR[j] == 1) { piv = piv_v[j]; - pivj = pivj_v[j]; //Line number - pivk = pivk_v[j]; //positi + pivj = pivj_v[j]; //Line number + pivk = pivk_v[j]; //positi markovitz_max = markovitz; NR_max = NR[j]; } @@ -6458,7 +6458,6 @@ dynSparseMatrix::solve_linear(const int block_num, const int y_size, const int y void dynSparseMatrix::solve_non_linear(const int block_num, const int y_size, const int y_kmin, const int y_kmax, const int size) - { max_res_idx = 0; bool cvg = false; diff --git a/mex/sources/bytecode/SparseMatrix.hh b/mex/sources/bytecode/SparseMatrix.hh index 953d0a229..017497f6a 100644 --- a/mex/sources/bytecode/SparseMatrix.hh +++ b/mex/sources/bytecode/SparseMatrix.hh @@ -86,7 +86,7 @@ struct t_save_op_s int first, second; }; -const int IFLD = 0; +const int IFLD = 0; const int IFDIV = 1; const int IFLESS = 2; const int IFSUB = 3; @@ -122,7 +122,7 @@ public: double g0, gp0, glambda2; int try_at_iteration; int find_exo_num(vector sconstrained_extended_path, int value); - int find_int_date(vector > per_value, int value); + int find_int_date(vector> per_value, int value); private: void Init_GE(int periods, int y_kmin, int y_kmax, int Size, map, int>, int> &IM); @@ -140,14 +140,14 @@ private: bool golden(double ax, double bx, double cx, double tol, double solve_tolf, double *xmin); void Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bool symbolic, int Block_number); bool Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, 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_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 Print_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, double *Ax, int n); void Printfull_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, double *Ax, double *b, int n); void PrintM(int n, double *Ax, mwIndex *Ap, mwIndex *Ai); - void Solve_LU_UMFPack(mxArray *A_m, mxArray *b_m, int Size, double slowc_l, bool is_two_boundaries, int it_); - void Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, double *Ax, double *b, int n, int Size, double slowc_l, bool is_two_boundaries, int it_, vector_table_conditional_local_type vector_table_conditional_local); - void Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, double *Ax, double *b, int n, int Size, double slowc_l, bool is_two_boundaries, int it_); + void Solve_LU_UMFPack(mxArray *A_m, mxArray *b_m, int Size, double slowc_l, bool is_two_boundaries, int it_); + void Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, double *Ax, double *b, int n, int Size, double slowc_l, bool is_two_boundaries, int it_, vector_table_conditional_local_type vector_table_conditional_local); + void Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, double *Ax, double *b, int n, int Size, double slowc_l, bool is_two_boundaries, int it_); void End_Matlab_LU_UMFPack(); #ifdef CUDA @@ -155,7 +155,7 @@ private: int *Ai, double *Ax, int *Ap, double *x0, double *b, double *A_tild, int *A_tild_i, int *A_tild_p, cusparseSolveAnalysisInfo_t infoL, cusparseSolveAnalysisInfo_t infoU, cusparseMatDescr_t descrL, cusparseMatDescr_t descrU, int preconditioner); - int Solve_CUDA_BiCGStab(int *Ap, int *Ai, double *Ax, int *Ap_tild, int *Ai_tild, double *A_tild, double *b, double *x0, int n, int Size, double slowc_l, bool is_two_boundaries, int it_, int nnz, int nnz_tild, int preconditioner, int max_iterations, int block); + int Solve_CUDA_BiCGStab(int *Ap, int *Ai, double *Ax, int *Ap_tild, int *Ai_tild, double *A_tild, double *b, double *x0, int n, int Size, double slowc_l, bool is_two_boundaries, int it_, int nnz, int nnz_tild, int preconditioner, int max_iterations, int block); #endif void Solve_Matlab_GMRES(mxArray *A_m, mxArray *b_m, int Size, double slowc, int block, bool is_two_boundaries, int it_, mxArray *x0_m); void Solve_Matlab_BiCGStab(mxArray *A_m, mxArray *b_m, int Size, double slowc, int block, bool is_two_boundaries, int it_, mxArray *x0_m, int precond); @@ -164,7 +164,7 @@ private: bool solve_linear(const int block_num, const int y_size, const int y_kmin, const int y_kmax, const int size, const int iter); void solve_non_linear(const int block_num, const int y_size, const int y_kmin, const int y_kmax, const int size); string preconditioner_print_out(string s, int preconditioner, bool ss); - bool compare(int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, long int nop4, int Size + bool compare(int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, long int nop4, int Size #ifdef PROFILER , long int *ndiv, long int *nsub #endif diff --git a/mex/sources/bytecode/bytecode.cc b/mex/sources/bytecode/bytecode.cc index 971d2704f..ee81f0843 100644 --- a/mex/sources/bytecode/bytecode.cc +++ b/mex/sources/bytecode/bytecode.cc @@ -208,7 +208,7 @@ GPU_Test_and_Info(cublasHandle_t *cublas_handle, cusparseHandle_t *cusparse_hand mexPrintf("> Driver version:\n"); int cuda_version; cuda_error = cudaDriverGetVersion(&cuda_version); - if (cuda_error != cudaSuccess) + if (cuda_error != cudaSuccess) { ostringstream tmp; tmp << " cudaGetVersion has failed\n"; @@ -297,7 +297,7 @@ Get_Arguments_and_global_variables(int nrhs, col_y = mxGetN(prhs[i]); break; case 1: - *xd = mxGetPr(prhs[i]); + *xd = mxGetPr(prhs[i]); row_x = mxGetM(prhs[i]); col_x = mxGetN(prhs[i]); break; @@ -349,7 +349,7 @@ Get_Arguments_and_global_variables(int nrhs, pos = pos1 + 1; else pos += 5; - block = atoi(Get_Argument(prhs[i]).substr(pos, string::npos).c_str())-1; + block = atoi(Get_Argument(prhs[i]).substr(pos, string::npos).c_str())-1; } else if (Get_Argument(prhs[i]).substr(0, 13) == "extended_path") { @@ -369,7 +369,7 @@ Get_Arguments_and_global_variables(int nrhs, pos = pos1 + 1; else pos += 6; - *pfplan_struct_name = deblank(Get_Argument(prhs[i]).substr(pos, string::npos)); + *pfplan_struct_name = deblank(Get_Argument(prhs[i]).substr(pos, string::npos)); } else if (Get_Argument(prhs[i]).substr(0, 4) == "plan") { @@ -378,7 +378,7 @@ Get_Arguments_and_global_variables(int nrhs, pos = pos1 + 1; else pos += 4; - *plan_struct_name = deblank(Get_Argument(prhs[i]).substr(pos, string::npos)); + *plan_struct_name = deblank(Get_Argument(prhs[i]).substr(pos, string::npos)); } else { @@ -582,7 +582,7 @@ main(int nrhs, const char *prhs[]) int nb_controlled = 0; mxArray *options_cond_fcst_ = mxGetField(extended_path_struct, 0, "options_cond_fcst_"); mxArray *controlled_varexo = NULL; - if (options_cond_fcst_ != NULL) + if (options_cond_fcst_ != NULL) { controlled_varexo = mxGetField(options_cond_fcst_, 0, "controlled_varexo"); nb_controlled = mxGetM(controlled_varexo) * mxGetN(controlled_varexo); @@ -704,7 +704,7 @@ main(int nrhs, const char *prhs[]) string tmp = "Can not allocated memory to store the date_str in the extended path descriptor"; DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str()); } - dates.push_back(string(buf));//string(Dates[i]); + dates.push_back(string(buf)); //string(Dates[i]); mxFree(buf); } } @@ -726,7 +726,7 @@ main(int nrhs, const char *prhs[]) mxArray *tmp = mxGetField(plan_struct, i, "exo"); if (tmp) { - char name [100]; + char name[100]; mxGetString(tmp, name, 100); splan[i].var = name; SymbolType variable_type = SymbolType::endogenous; @@ -744,7 +744,7 @@ main(int nrhs, const char *prhs[]) tmp = mxGetField(plan_struct, i, "var"); if (tmp) { - char name [100]; + char name[100]; mxGetString(tmp, name, 100); splan[i].exo = name; SymbolType variable_type; @@ -778,7 +778,7 @@ main(int nrhs, const char *prhs[]) mexPrintf(" plan fliping var=%s (%d) exo=%s (%d) for the following periods and with the following values:\n", it->var.c_str(), it->var_num, it->exo.c_str(), it->exo_num); else mexPrintf(" plan shocks on var=%s for the following periods and with the following values:\n", it->var.c_str()); - for (vector >::iterator it1 = it->per_value.begin(); it1 != it->per_value.end(); it1++) + for (vector>::iterator it1 = it->per_value.begin(); it1 != it->per_value.end(); it1++) { mexPrintf(" %3d %10.5f\n", it1->first, it1->second); } @@ -804,7 +804,7 @@ main(int nrhs, const char *prhs[]) mxArray *tmp = mxGetField(pfplan_struct, i, "var"); if (tmp) { - char name [100]; + char name[100]; mxGetString(tmp, name, 100); spfplan[i].var = name; SymbolType variable_type = SymbolType::endogenous; @@ -822,7 +822,7 @@ main(int nrhs, const char *prhs[]) tmp = mxGetField(pfplan_struct, i, "exo"); if (tmp) { - char name [100]; + char name[100]; mxGetString(tmp, name, 100); spfplan[i].exo = name; SymbolType variable_type; @@ -856,7 +856,7 @@ main(int nrhs, const char *prhs[]) mexPrintf(" plan flipping var=%s (%d) exo=%s (%d) for the following periods and with the following values:\n", it->var.c_str(), it->var_num, it->exo.c_str(), it->exo_num); else mexPrintf(" plan shocks on var=%s (%d) for the following periods and with the following values:\n", it->var.c_str(), it->var_num); - for (vector >::iterator it1 = it->per_value.begin(); it1 != it->per_value.end(); it1++) + for (vector>::iterator it1 = it->per_value.begin(); it1 != it->per_value.end(); it1++) { mexPrintf(" %3d %10.5f\n", it1->first, it1->second); } @@ -1067,7 +1067,7 @@ main(int nrhs, const char *prhs[]) } for (i = 0; i < row_y*col_y; i++) { - y[i] = double (yd[i]); + y[i] = double (yd[i]); ya[i] = double (yd[i]); } size_t y_size = row_y; diff --git a/mex/sources/bytecode/testing/mex_interface.cc b/mex/sources/bytecode/testing/mex_interface.cc index 2eac721fd..5094185a9 100644 --- a/mex/sources/bytecode/testing/mex_interface.cc +++ b/mex/sources/bytecode/testing/mex_interface.cc @@ -347,7 +347,7 @@ read_struct(FILE *fid) fscanf(fid, "%d", &size_2); vector fieldnames; vector v_Array; - vector > vv_Array; + vector> vv_Array; for (unsigned int j = 0; j < size_1 * size_2; j++) { for (unsigned int i = 0; i < nfields; i++) diff --git a/mex/sources/bytecode/testing/mex_interface.hh b/mex/sources/bytecode/testing/mex_interface.hh index ede5e63d7..9e8989189 100644 --- a/mex/sources/bytecode/testing/mex_interface.hh +++ b/mex/sources/bytecode/testing/mex_interface.hh @@ -41,27 +41,27 @@ typedef unsigned int mwSize; enum mxData_type { - mxREAL, - mxCOMPLEX + mxREAL, + mxCOMPLEX }; enum mxArray_type { - mxDOUBLE_CLASS = 0, - mxSINGLE_CLASS = 1, - mxLOGICAL_CLASS = 2, - mxCHAR_CLASS = 3, - mxSPARSE_CLASS = 4, - mxSTRUCT_CLASS = 5, - mxCELL_CLASS = 6, - mxFUNCTION_CLASS = 7, - mxUINT8_CLASS = 8, - mxINT16_CLASS = 9, - mxUINT16_CLASS = 10, - mxINT32_CLASS = 11, - mxUINT32_CLASS = 12, - mxINT64_CLASS = 13, - mxUINT64_CLASS = 14 + mxDOUBLE_CLASS = 0, + mxSINGLE_CLASS = 1, + mxLOGICAL_CLASS = 2, + mxCHAR_CLASS = 3, + mxSPARSE_CLASS = 4, + mxSTRUCT_CLASS = 5, + mxCELL_CLASS = 6, + mxFUNCTION_CLASS = 7, + mxUINT8_CLASS = 8, + mxINT16_CLASS = 9, + mxUINT16_CLASS = 10, + mxINT32_CLASS = 11, + mxUINT32_CLASS = 12, + mxINT64_CLASS = 13, + mxUINT64_CLASS = 14 }; class mxArray_tag @@ -73,7 +73,7 @@ public: double *data; mxArray_type type; vector field_name; - vector > field_array; + vector> field_array; mxArray_tag(); }; diff --git a/mex/sources/k_order_perturbation/dynamic_abstract_class.hh b/mex/sources/k_order_perturbation/dynamic_abstract_class.hh index 901782670..e1f8c3709 100644 --- a/mex/sources/k_order_perturbation/dynamic_abstract_class.hh +++ b/mex/sources/k_order_perturbation/dynamic_abstract_class.hh @@ -29,7 +29,9 @@ class DynamicModelAC protected: int ntt; // Size of vector of temporary terms public: - DynamicModelAC(int ntt_arg) : ntt{ntt_arg} {}; + DynamicModelAC(int ntt_arg) : ntt{ntt_arg} + { + }; virtual ~DynamicModelAC() = default; virtual void eval(const Vector &y, const Vector &x, const Vector ¶ms, const Vector &ySteady, Vector &residual, std::vector &md) = 0; diff --git a/mex/sources/k_order_perturbation/dynamic_dll.hh b/mex/sources/k_order_perturbation/dynamic_dll.hh index 2e7d87195..693e8b2b9 100644 --- a/mex/sources/k_order_perturbation/dynamic_dll.hh +++ b/mex/sources/k_order_perturbation/dynamic_dll.hh @@ -35,7 +35,7 @@ #include "dynamic_abstract_class.hh" using dynamic_tt_fct = void (*)(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, double *T); -using dynamic_deriv_fct = void (*) (const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, const double *T, double *deriv); +using dynamic_deriv_fct = void (*)(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, const double *T, double *deriv); /** * creates pointer to Dynamic function inside _dynamic.dll @@ -47,7 +47,7 @@ private: std::vector dynamic_tt; std::vector dynamic_deriv; #if defined(_WIN32) || defined(__CYGWIN32__) - HINSTANCE dynamicHinstance; // DLL instance pointer in Windows + HINSTANCE dynamicHinstance; // DLL instance pointer in Windows #else void *dynamicHinstance; // and in Linux or Mac #endif diff --git a/mex/sources/k_order_perturbation/k_ord_dynare.cc b/mex/sources/k_order_perturbation/k_ord_dynare.cc index 80d05179c..5ccd442fb 100644 --- a/mex/sources/k_order_perturbation/k_ord_dynare.cc +++ b/mex/sources/k_order_perturbation/k_ord_dynare.cc @@ -236,4 +236,3 @@ DynareStateNameList::DynareStateNameList(const KordpDynare &dynare, const Dynare for (int i = 0; i < dynare.nexog(); i++) names.emplace_back(denl.getName(i)); } - diff --git a/mex/sources/k_order_perturbation/k_order_perturbation.cc b/mex/sources/k_order_perturbation/k_order_perturbation.cc index 03578a97a..0a44e8206 100644 --- a/mex/sources/k_order_perturbation/k_order_perturbation.cc +++ b/mex/sources/k_order_perturbation/k_order_perturbation.cc @@ -179,7 +179,7 @@ extern "C" { int ntt = std::accumulate(mxGetPr(dynamic_tmp_nbr_mx), mxGetPr(dynamic_tmp_nbr_mx)+kOrder+1, 0); // Extract various fields from dr - const mxArray *ys_mx = mxGetField(dr_mx, 0, "ys"); // and not in order of dr.order_var + const mxArray *ys_mx = mxGetField(dr_mx, 0, "ys"); // and not in order of dr.order_var if (!(ys_mx && mxIsDouble(ys_mx))) DYN_MEX_FUNC_ERR_MSG_TXT("dr.ys should be a double precision array"); Vector ySteady{ConstVector{ys_mx}}; @@ -312,4 +312,3 @@ extern "C" { plhs[0] = mxCreateDoubleScalar(0); } // end of mexFunction() } // end of extern C - diff --git a/mex/sources/k_order_perturbation/tests/k_order_test_main.cc b/mex/sources/k_order_perturbation/tests/k_order_test_main.cc index d924067ac..114e1e1e7 100644 --- a/mex/sources/k_order_perturbation/tests/k_order_test_main.cc +++ b/mex/sources/k_order_perturbation/tests/k_order_test_main.cc @@ -49,7 +49,7 @@ main(int argc, char *argv[]) 0.7000, 0.7870, 0.0200}; - Vector *modParams = new Vector(dparams, npar); + Vector *modParams = new Vector(dparams, npar); #ifdef DEBUG mexPrintf("k_ord_perturbation: nParams=%d .\n", npar); @@ -65,24 +65,24 @@ main(int argc, char *argv[]) #endif double d2Dparams[4] = { //(double *) mxGetData(mxFldp); - 0.1960e-3, 0.0, - 0.0, 0.0250e-3 + 0.1960e-3, 0.0, + 0.0, 0.0250e-3 }; npar = 2; //(int)mxGetN(mxFldp); - TwoDMatrix *vCov = new TwoDMatrix(npar, npar, (d2Dparams)); - double dYSparams [16] = { // 27 mxGetData(mxFldp); - // 1.0110, 2.2582, 5.8012, 0.5808, - 1.0110, 2.2582, 0.4477, 1.0000, - 4.5959, 1.0212, 5.8012, 0.8494, - 0.1872, 0.8604, 1.0030, 1.0080, - 0.5808, 1.0030, 2.2582, 0.4477 - //, 1.0110, 2.2582, 0.4477, 1.0000, 0.1872, 2.2582, 0.4477 + TwoDMatrix *vCov = new TwoDMatrix(npar, npar, (d2Dparams)); + double dYSparams[16] = { // 27 mxGetData(mxFldp); + // 1.0110, 2.2582, 5.8012, 0.5808, + 1.0110, 2.2582, 0.4477, 1.0000, + 4.5959, 1.0212, 5.8012, 0.8494, + 0.1872, 0.8604, 1.0030, 1.0080, + 0.5808, 1.0030, 2.2582, 0.4477 + //, 1.0110, 2.2582, 0.4477, 1.0000, 0.1872, 2.2582, 0.4477 }; const int nSteady = 16; //27 //31;//29, 16 (int)mxGetM(mxFldp); - Vector *ySteady = new Vector(dYSparams, nSteady); + Vector *ySteady = new Vector(dYSparams, nSteady); double nnzd[3] = { 77, 217, 0}; - const Vector *NNZD = new Vector(nnzd, 3); + const Vector *NNZD = new Vector(nnzd, 3); //mxFldp = mxGetField(dr, 0,"nstatic" ); const int nStat = 7; //(int)mxGetScalar(mxFldp); @@ -109,32 +109,32 @@ main(int argc, char *argv[]) int var_order[] //[18] = { - 5, 6, 8, 10, 11, 12, 14, 7, 13, 1, 2, 3, 4, 9, 15, 16 - // 5, 6, 8, 10, 11, 12, 16, 7, 13, 14, 15, 1, 2, 3, 4, 9, 17, 18 + 5, 6, 8, 10, 11, 12, 14, 7, 13, 1, 2, 3, 4, 9, 15, 16 + // 5, 6, 8, 10, 11, 12, 16, 7, 13, 14, 15, 1, 2, 3, 4, 9, 17, 18 }; //Vector * varOrder = new Vector(var_order, nEndo); vector *var_order_vp = new vector(nEndo); //nEndo)); for (int v = 0; v < nEndo; v++) (*var_order_vp)[v] = var_order[v]; - const double ll_incidence [] //[3][18] + const double ll_incidence[] //[3][18] = { - 1, 5, 21, - 2, 6, 22, - 0, 7, 23, - 0, 8, 24, - 0, 9, 0, - 0, 10, 0, - 3, 11, 0, - 0, 12, 0, - 0, 13, 25, - 0, 14, 0, - 0, 15, 0, - 0, 16, 0, - 4, 17, 0, - 0, 18, 0, - 0, 19, 26, - 0, 20, 27 + 1, 5, 21, + 2, 6, 22, + 0, 7, 23, + 0, 8, 24, + 0, 9, 0, + 0, 10, 0, + 3, 11, 0, + 0, 12, 0, + 0, 13, 25, + 0, 14, 0, + 0, 15, 0, + 0, 16, 0, + 4, 17, 0, + 0, 18, 0, + 0, 19, 26, + 0, 20, 27 }; TwoDMatrix *llincidence = new TwoDMatrix(3, nEndo, ll_incidence); @@ -166,7 +166,7 @@ main(int argc, char *argv[]) mexPrintf("k_ord_perturbation: ExoNameList[%d][0]= %s.\n", i, exoNamesMX[i]); } #endif - if ((nEndo != nendo) || (nExog != nexo)) //(nPar != npar) + if ((nEndo != nendo) || (nExog != nexo)) //(nPar != npar) { mexPrintf("Incorrect number of input parameters.\n"); return; @@ -220,7 +220,7 @@ main(int argc, char *argv[]) mexPrintf("k_order_perturbation: Calling dynare constructor.\n"); #endif // make KordpDynare object - KordpDynare dynare(endoNamesMX, nEndo, exoNamesMX, nExog, nPar, // paramNames, + KordpDynare dynare(endoNamesMX, nEndo, exoNamesMX, nExog, nPar, // paramNames, ySteady, vCov, modParams, nStat, nPred, nForw, nBoth, jcols, NNZD, nSteps, kOrder, journal, dynamicDLL, sstol, var_order_vp, //var_order llincidence, qz_criterium); @@ -314,7 +314,7 @@ main(int argc, char *argv[]) mexPrintf("k_order_perturbation: Filling MATLAB outputs.\n"); #endif - double *dgy, *dgu, *ysteady; + double *dgy, *dgu, *ysteady; int nb_row_x; ysteady = NULL; diff --git a/mex/sources/kalman_steady_state/kalman_steady_state.cc b/mex/sources/kalman_steady_state/kalman_steady_state.cc index 7e455e347..0a0db3181 100644 --- a/mex/sources/kalman_steady_state/kalman_steady_state.cc +++ b/mex/sources/kalman_steady_state/kalman_steady_state.cc @@ -95,7 +95,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) if (static_cast(q) != mxGetN(prhs[1])) DYN_MEX_FUNC_ERR_MSG_TXT("kalman_steady_state: The second input argument (QQ) must be a square matrix!"); if (mxIsNumeric(prhs[1]) == 0 || mxIsComplex(prhs[1]) == 1) - DYN_MEX_FUNC_ERR_MSG_TXT("kalman_steady_state: The second input argument (QQ) must be a real matrix!"); + DYN_MEX_FUNC_ERR_MSG_TXT("kalman_steady_state: The second input argument (QQ) must be a real matrix!"); if (q != n) DYN_MEX_FUNC_ERR_MSG_TXT("kalman_steady_state: The size of the second input argument (QQ) must match the size of the first argument (T)!"); lapack_int p = mxGetN(prhs[2]); diff --git a/mex/sources/local_state_space_iterations/local_state_space_iteration_2.cc b/mex/sources/local_state_space_iterations/local_state_space_iteration_2.cc index 15b422771..d0dc82e16 100644 --- a/mex/sources/local_state_space_iterations/local_state_space_iteration_2.cc +++ b/mex/sources/local_state_space_iterations/local_state_space_iteration_2.cc @@ -31,7 +31,7 @@ #include -#define FIRST_ORDER_LOOP 1// Comment out this line to use mkl-blas instead of loops when computing ghx*yhat and ghu*epsilon +#define FIRST_ORDER_LOOP 1 // Comment out this line to use mkl-blas instead of loops when computing ghx*yhat and ghu*epsilon std::tuple, std::vector, std::vector> set_vector_of_indices(int n, int r) @@ -216,7 +216,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) prhs[7] ghxu [double] m×nq array, second order reduced form. prhs[8] yhat_ [double] [OPTIONAL] n×s array, time t particles (pruning additional latent variables). prhs[9] ss [double] [OPTIONAL] m×1 array, steady state for the union of the states and the observed variables (needed for the pruning mode). - + prhs[9 or 11] [double] num of threads plhs[0] y [double] n×s array, time t+1 particles. @@ -231,14 +231,14 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) mexErrMsgTxt("Too many output arguments."); // Get dimensions. - size_t n = mxGetM(prhs[0]);// Number of states. - size_t s = mxGetN(prhs[0]);// Number of particles. - size_t q = mxGetM(prhs[1]);// Number of innovations. - size_t m = mxGetM(prhs[2]);// Number of elements in the union of states and observed variables. + size_t n = mxGetM(prhs[0]); // Number of states. + size_t s = mxGetN(prhs[0]); // Number of particles. + size_t q = mxGetM(prhs[1]); // Number of innovations. + size_t m = mxGetM(prhs[2]); // Number of elements in the union of states and observed variables. //mexPrintf("\n s (the number of column of yhat) is equal to %d.", s); //mexPrintf("\n The number of column of epsilon is %d.", mxGetN(prhs[1])); // Check the dimensions. - if (s != mxGetN(prhs[1]) // Number of columns for epsilon + if (s != mxGetN(prhs[1]) // Number of columns for epsilon || n != mxGetN(prhs[2]) // Number of columns for ghx || m != mxGetM(prhs[3]) // Number of rows for ghu || q != mxGetN(prhs[3]) // Number of columns for ghu @@ -248,12 +248,12 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) || m != mxGetM(prhs[6]) // Number of rows for ghuu || q*q != mxGetN(prhs[6]) // Number of columns for ghuu || m != mxGetM(prhs[7]) // Number of rows for ghxu - || n*q != mxGetN(prhs[7])) // Number of rows for ghxu + || n*q != mxGetN(prhs[7])) // Number of rows for ghxu mexErrMsgTxt("Input dimension mismatch!."); if (nrhs > 9) - if (n != mxGetM(prhs[8]) // Number of rows for yhat_ - || s != mxGetN(prhs[8]) // Number of columns for yhat_ - || m != mxGetM(prhs[9])) // Number of rows for ss + if (n != mxGetM(prhs[8]) // Number of rows for yhat_ + || s != mxGetN(prhs[8]) // Number of columns for yhat_ + || m != mxGetM(prhs[9])) // Number of rows for ss mexErrMsgTxt("Input dimension mismatch!."); // Get Input arrays. diff --git a/mex/sources/local_state_space_iterations/local_state_space_iteration_k.cc b/mex/sources/local_state_space_iterations/local_state_space_iteration_k.cc index dcc80f6ea..43e9a4779 100644 --- a/mex/sources/local_state_space_iterations/local_state_space_iteration_k.cc +++ b/mex/sources/local_state_space_iterations/local_state_space_iteration_k.cc @@ -37,13 +37,13 @@ struct ParticleWorker : public sthread::detach_thread { const int npred_both, exo_nbr; - const std::pair particle_range; + const std::pair particle_range; const ConstGeneralMatrix &yhat, ε const Vector &ys_reordered; const UnfoldDecisionRule &dr; GeneralMatrix &ynext; - ParticleWorker(int npred_both_arg, int exo_nbr_arg, std::pair particle_range_arg, + ParticleWorker(int npred_both_arg, int exo_nbr_arg, std::pair particle_range_arg, const ConstGeneralMatrix &yhat_arg, const ConstGeneralMatrix &epsilon_arg, const Vector &ys_reordered_arg, const UnfoldDecisionRule &dr_arg, GeneralMatrix &ynext_arg) @@ -52,7 +52,8 @@ struct ParticleWorker : public sthread::detach_thread ynext{ynext_arg} { } - void operator()(std::mutex &mut) override + void + operator()(std::mutex &mut) override { Vector dyu(npred_both+exo_nbr); Vector dy(dyu, 0, npred_both); diff --git a/mex/sources/ms-sbvar/mex_top_level.cc b/mex/sources/ms-sbvar/mex_top_level.cc index 7c2f6d56b..e4060c2be 100644 --- a/mex/sources/ms-sbvar/mex_top_level.cc +++ b/mex/sources/ms-sbvar/mex_top_level.cc @@ -48,16 +48,16 @@ mexFunction(int nlhs, mxArray *plhs[], /* * Allocate memory */ - maxnargs = static_cast (mxGetN(prhs[0])/2+1); - argument = static_cast (mxCalloc(mxGetN(prhs[0])+1, sizeof(char))); - args = static_cast (mxCalloc(maxnargs, sizeof(char *))); + maxnargs = static_cast(mxGetN(prhs[0])/2+1); + argument = static_cast(mxCalloc(mxGetN(prhs[0])+1, sizeof(char))); + args = static_cast(mxCalloc(maxnargs, sizeof(char *))); if (argument == NULL || args == NULL) DYN_MEX_FUNC_ERR_MSG_TXT("Error in MS-SBVAR MEX file: could not allocate memory. (1)"); /* * Create argument string from prhs and parse to create args / nargs */ - if (!(args[nargs] = static_cast (mxCalloc(strlen(mainarg)+1, sizeof(char))))) + if (!(args[nargs] = static_cast(mxCalloc(strlen(mainarg)+1, sizeof(char))))) DYN_MEX_FUNC_ERR_MSG_TXT("Error in MS-SBVAR MEX file: could not allocate memory. (2)"); strncpy(args[nargs++], mainarg, strlen(mainarg)); @@ -68,7 +68,7 @@ mexFunction(int nlhs, mxArray *plhs[], beginarg = &argument[0]; while ((n = strcspn(beginarg, " "))) { - if (!(args[nargs] = static_cast (mxCalloc(n+1, sizeof(char))))) + if (!(args[nargs] = static_cast(mxCalloc(n+1, sizeof(char))))) DYN_MEX_FUNC_ERR_MSG_TXT("Error in MS-SBVAR MEX file: could not allocate memory. (3)"); strncpy(args[nargs++], beginarg, n); diff --git a/mex/sources/ms-sbvar/modify_for_mex.cc b/mex/sources/ms-sbvar/modify_for_mex.cc index 42bbe64fa..bdcf2a85d 100644 --- a/mex/sources/ms-sbvar/modify_for_mex.cc +++ b/mex/sources/ms-sbvar/modify_for_mex.cc @@ -19,12 +19,12 @@ #if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE) -#include +# include -#ifdef __cplusplus +# ifdef __cplusplus extern "C" { -#endif +# endif int constant_seed; @@ -34,7 +34,7 @@ extern "C" throw "Error in MS-SBVAR MEX file.\n"; } -#ifdef __cplusplus +# ifdef __cplusplus } -#endif +# endif #endif diff --git a/mex/sources/perfect_foresight_problem/DynamicModelCaller.hh b/mex/sources/perfect_foresight_problem/DynamicModelCaller.hh index 4dd1d02c8..abfc019f1 100644 --- a/mex/sources/perfect_foresight_problem/DynamicModelCaller.hh +++ b/mex/sources/perfect_foresight_problem/DynamicModelCaller.hh @@ -39,7 +39,9 @@ public: // Used to store error messages (as exceptions cannot cross the OpenMP boundary) static std::string error_msg; - DynamicModelCaller(bool linear_arg, bool compute_jacobian_arg) : linear{linear_arg}, compute_jacobian{compute_jacobian_arg} {}; + DynamicModelCaller(bool linear_arg, bool compute_jacobian_arg) : linear{linear_arg}, compute_jacobian{compute_jacobian_arg} + { + }; virtual ~DynamicModelCaller() = default; virtual double &y(size_t i) const = 0; virtual double jacobian(size_t i) const = 0; @@ -55,7 +57,7 @@ private: static void *dynamic_mex; #endif using dynamic_tt_fct = void (*)(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, double *T); - using dynamic_deriv_fct = void (*) (const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, const double *T, double *deriv); + using dynamic_deriv_fct = void (*)(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, const double *T, double *deriv); static dynamic_tt_fct residual_tt_fct, g1_tt_fct; static dynamic_deriv_fct residual_fct, g1_fct; size_t nb_row_x; @@ -65,8 +67,16 @@ private: public: DynamicModelDllCaller(size_t ntt, mwIndex nx, mwIndex ny, size_t ndynvars, const double *x_arg, size_t nb_row_x_arg, const double *params_arg, const double *steady_state_arg, bool linear_arg, bool compute_jacobian_arg); virtual ~DynamicModelDllCaller() = default; - double &y(size_t i) const override { return y_p[i]; }; - double jacobian(size_t i) const override { return jacobian_p[i]; }; + double & + y(size_t i) const override + { + return y_p[i]; + }; + double + jacobian(size_t i) const override + { + return jacobian_p[i]; + }; void eval(int it, double *resid) override; static void load_dll(const std::string &basename); static void unload_dll(); @@ -85,13 +95,23 @@ private: public: DynamicModelMatlabCaller(std::string basename_arg, size_t ntt, size_t ndynvars, const mxArray *x_mx_arg, const mxArray *params_mx_arg, const mxArray *steady_state_mx_arg, bool linear_arg, bool compute_jacobian_arg); ~DynamicModelMatlabCaller() override; - double &y(size_t i) const override { return mxGetPr(y_mx)[i]; }; - double jacobian(size_t i) const override { return jacobian_mx ? mxGetPr(jacobian_mx)[i] : std::numeric_limits::quiet_NaN(); }; + double & + y(size_t i) const override + { + return mxGetPr(y_mx)[i]; + }; + double + jacobian(size_t i) const override + { + return jacobian_mx ? mxGetPr(jacobian_mx)[i] : std::numeric_limits::quiet_NaN(); + }; void eval(int it, double *resid) override; - class Exception { + class Exception + { public: const std::string msg; - Exception(std::string msg_arg) : msg{std::move(msg_arg)} {}; + Exception(std::string msg_arg) : msg{std::move(msg_arg)} + { + }; }; }; - diff --git a/mex/sources/sobol/gaussian.hh b/mex/sources/sobol/gaussian.hh index ea5377cb9..a259c66a9 100644 --- a/mex/sources/sobol/gaussian.hh +++ b/mex/sources/sobol/gaussian.hh @@ -45,36 +45,36 @@ icdf(const T uniform) { static T A[6] = { - -3.969683028665376e+01, - 2.209460984245205e+02, - -2.759285104469687e+02, - 1.383577518672690e+02, - -3.066479806614716e+01, - 2.506628277459239e+00 + -3.969683028665376e+01, + 2.209460984245205e+02, + -2.759285104469687e+02, + 1.383577518672690e+02, + -3.066479806614716e+01, + 2.506628277459239e+00 }; static T B[5] = { - -5.447609879822406e+01, - 1.615858368580409e+02, - -1.556989798598866e+02, - 6.680131188771972e+01, - -1.328068155288572e+01 + -5.447609879822406e+01, + 1.615858368580409e+02, + -1.556989798598866e+02, + 6.680131188771972e+01, + -1.328068155288572e+01 }; static T C[6] = { - -7.784894002430293e-03, - -3.223964580411365e-01, - -2.400758277161838e+00, - -2.549732539343734e+00, - 4.374664141464968e+00, - 2.938163982698783e+00 + -7.784894002430293e-03, + -3.223964580411365e-01, + -2.400758277161838e+00, + -2.549732539343734e+00, + 4.374664141464968e+00, + 2.938163982698783e+00 }; static T D[4] = { - 7.784695709041462e-03, - 3.224671290700398e-01, - 2.445134137142996e+00, - 3.754408661907416e+00 + 7.784695709041462e-03, + 3.224671290700398e-01, + 2.445134137142996e+00, + 3.754408661907416e+00 }; T gaussian = static_cast(0.0); if (0 < uniform && uniform < lb) @@ -148,15 +148,15 @@ usphere(int d, int n, T *U) { icdfm(n*d, U); #pragma omp parallel for - for (int j = 0; j < n; j++)// sequence index. + for (int j = 0; j < n; j++) // sequence index. { int k = j*d; double norm = 0.0; - for (int i = 0; i < d; i++)// dimension index. + for (int i = 0; i < d; i++) // dimension index. norm = norm + U[k+i]*U[k+i]; norm = sqrt(norm); - for (int i = 0; i < d; i++)// dimension index. + for (int i = 0; i < d; i++) // dimension index. U[k+i] = U[k+i]/norm; } } @@ -167,15 +167,15 @@ usphereRadius(int d, int n, double radius, T *U) { icdfm(n*d, U); #pragma omp parallel for - for (int j = 0; j < n; j++)// sequence index. + for (int j = 0; j < n; j++) // sequence index. { int k = j*d; double norm = 0.0; - for (int i = 0; i < d; i++)// dimension index. + for (int i = 0; i < d; i++) // dimension index. norm = norm + U[k+i]*U[k+i]; norm = sqrt(norm); - for (int i = 0; i < d; i++)// dimension index. + for (int i = 0; i < d; i++) // dimension index. U[k+i] = radius*U[k+i]/norm; } } diff --git a/mex/sources/sobol/sobol.hh b/mex/sources/sobol/sobol.hh index 1e80d6be1..a6f29d9f7 100644 --- a/mex/sources/sobol/sobol.hh +++ b/mex/sources/sobol/sobol.hh @@ -251,118 +251,118 @@ next_sobol(int dim_num, T1 *seed, T2 *quasi) T1 l = 0; static T1 poly[DIM_MAX] = { - 1, 3, 7, 11, 13, 19, 25, 37, 59, 47, - 61, 55, 41, 67, 97, 91, 109, 103, 115, 131, - 193, 137, 145, 143, 241, 157, 185, 167, 229, 171, - 213, 191, 253, 203, 211, 239, 247, 285, 369, 299, - 301, 333, 351, 355, 357, 361, 391, 397, 425, 451, - 463, 487, 501, 529, 539, 545, 557, 563, 601, 607, - 617, 623, 631, 637, 647, 661, 675, 677, 687, 695, - 701, 719, 721, 731, 757, 761, 787, 789, 799, 803, - 817, 827, 847, 859, 865, 875, 877, 883, 895, 901, - 911, 949, 953, 967, 971, 973, 981, 985, 995, 1001, - 1019, 1033, 1051, 1063, 1069, 1125, 1135, 1153, 1163, 1221, - 1239, 1255, 1267, 1279, 1293, 1305, 1315, 1329, 1341, 1347, - 1367, 1387, 1413, 1423, 1431, 1441, 1479, 1509, 1527, 1531, - 1555, 1557, 1573, 1591, 1603, 1615, 1627, 1657, 1663, 1673, - 1717, 1729, 1747, 1759, 1789, 1815, 1821, 1825, 1849, 1863, - 1869, 1877, 1881, 1891, 1917, 1933, 1939, 1969, 2011, 2035, - 2041, 2053, 2071, 2091, 2093, 2119, 2147, 2149, 2161, 2171, - 2189, 2197, 2207, 2217, 2225, 2255, 2257, 2273, 2279, 2283, - 2293, 2317, 2323, 2341, 2345, 2363, 2365, 2373, 2377, 2385, - 2395, 2419, 2421, 2431, 2435, 2447, 2475, 2477, 2489, 2503, - 2521, 2533, 2551, 2561, 2567, 2579, 2581, 2601, 2633, 2657, - 2669, 2681, 2687, 2693, 2705, 2717, 2727, 2731, 2739, 2741, - 2773, 2783, 2793, 2799, 2801, 2811, 2819, 2825, 2833, 2867, - 2879, 2881, 2891, 2905, 2911, 2917, 2927, 2941, 2951, 2955, - 2963, 2965, 2991, 2999, 3005, 3017, 3035, 3037, 3047, 3053, - 3083, 3085, 3097, 3103, 3159, 3169, 3179, 3187, 3205, 3209, - 3223, 3227, 3229, 3251, 3263, 3271, 3277, 3283, 3285, 3299, - 3305, 3319, 3331, 3343, 3357, 3367, 3373, 3393, 3399, 3413, - 3417, 3427, 3439, 3441, 3475, 3487, 3497, 3515, 3517, 3529, - 3543, 3547, 3553, 3559, 3573, 3589, 3613, 3617, 3623, 3627, - 3635, 3641, 3655, 3659, 3669, 3679, 3697, 3707, 3709, 3713, - 3731, 3743, 3747, 3771, 3791, 3805, 3827, 3833, 3851, 3865, - 3889, 3895, 3933, 3947, 3949, 3957, 3971, 3985, 3991, 3995, - 4007, 4013, 4021, 4045, 4051, 4069, 4073, 4179, 4201, 4219, - 4221, 4249, 4305, 4331, 4359, 4383, 4387, 4411, 4431, 4439, - 4449, 4459, 4485, 4531, 4569, 4575, 4621, 4663, 4669, 4711, - 4723, 4735, 4793, 4801, 4811, 4879, 4893, 4897, 4921, 4927, - 4941, 4977, 5017, 5027, 5033, 5127, 5169, 5175, 5199, 5213, - 5223, 5237, 5287, 5293, 5331, 5391, 5405, 5453, 5523, 5573, - 5591, 5597, 5611, 5641, 5703, 5717, 5721, 5797, 5821, 5909, - 5913, 5955, 5957, 6005, 6025, 6061, 6067, 6079, 6081, 6231, - 6237, 6289, 6295, 6329, 6383, 6427, 6453, 6465, 6501, 6523, - 6539, 6577, 6589, 6601, 6607, 6631, 6683, 6699, 6707, 6761, - 6795, 6865, 6881, 6901, 6923, 6931, 6943, 6999, 7057, 7079, - 7103, 7105, 7123, 7173, 7185, 7191, 7207, 7245, 7303, 7327, - 7333, 7355, 7365, 7369, 7375, 7411, 7431, 7459, 7491, 7505, - 7515, 7541, 7557, 7561, 7701, 7705, 7727, 7749, 7761, 7783, - 7795, 7823, 7907, 7953, 7963, 7975, 8049, 8089, 8123, 8125, - 8137, 8219, 8231, 8245, 8275, 8293, 8303, 8331, 8333, 8351, - 8357, 8367, 8379, 8381, 8387, 8393, 8417, 8435, 8461, 8469, - 8489, 8495, 8507, 8515, 8551, 8555, 8569, 8585, 8599, 8605, - 8639, 8641, 8647, 8653, 8671, 8675, 8689, 8699, 8729, 8741, - 8759, 8765, 8771, 8795, 8797, 8825, 8831, 8841, 8855, 8859, - 8883, 8895, 8909, 8943, 8951, 8955, 8965, 8999, 9003, 9031, - 9045, 9049, 9071, 9073, 9085, 9095, 9101, 9109, 9123, 9129, - 9137, 9143, 9147, 9185, 9197, 9209, 9227, 9235, 9247, 9253, - 9257, 9277, 9297, 9303, 9313, 9325, 9343, 9347, 9371, 9373, - 9397, 9407, 9409, 9415, 9419, 9443, 9481, 9495, 9501, 9505, - 9517, 9529, 9555, 9557, 9571, 9585, 9591, 9607, 9611, 9621, - 9625, 9631, 9647, 9661, 9669, 9679, 9687, 9707, 9731, 9733, - 9745, 9773, 9791, 9803, 9811, 9817, 9833, 9847, 9851, 9863, - 9875, 9881, 9905, 9911, 9917, 9923, 9963, 9973, 10003, 10025, - 10043, 10063, 10071, 10077, 10091, 10099, 10105, 10115, 10129, 10145, - 10169, 10183, 10187, 10207, 10223, 10225, 10247, 10265, 10271, 10275, - 10289, 10299, 10301, 10309, 10343, 10357, 10373, 10411, 10413, 10431, - 10445, 10453, 10463, 10467, 10473, 10491, 10505, 10511, 10513, 10523, - 10539, 10549, 10559, 10561, 10571, 10581, 10615, 10621, 10625, 10643, - 10655, 10671, 10679, 10685, 10691, 10711, 10739, 10741, 10755, 10767, - 10781, 10785, 10803, 10805, 10829, 10857, 10863, 10865, 10875, 10877, - 10917, 10921, 10929, 10949, 10967, 10971, 10987, 10995, 11009, 11029, - 11043, 11045, 11055, 11063, 11075, 11081, 11117, 11135, 11141, 11159, - 11163, 11181, 11187, 11225, 11237, 11261, 11279, 11297, 11307, 11309, - 11327, 11329, 11341, 11377, 11403, 11405, 11413, 11427, 11439, 11453, - 11461, 11473, 11479, 11489, 11495, 11499, 11533, 11545, 11561, 11567, - 11575, 11579, 11589, 11611, 11623, 11637, 11657, 11663, 11687, 11691, - 11701, 11747, 11761, 11773, 11783, 11795, 11797, 11817, 11849, 11855, - 11867, 11869, 11873, 11883, 11919, 11921, 11927, 11933, 11947, 11955, - 11961, 11999, 12027, 12029, 12037, 12041, 12049, 12055, 12095, 12097, - 12107, 12109, 12121, 12127, 12133, 12137, 12181, 12197, 12207, 12209, - 12239, 12253, 12263, 12269, 12277, 12287, 12295, 12309, 12313, 12335, - 12361, 12367, 12391, 12409, 12415, 12433, 12449, 12469, 12479, 12481, - 12499, 12505, 12517, 12527, 12549, 12559, 12597, 12615, 12621, 12639, - 12643, 12657, 12667, 12707, 12713, 12727, 12741, 12745, 12763, 12769, - 12779, 12781, 12787, 12799, 12809, 12815, 12829, 12839, 12857, 12875, - 12883, 12889, 12901, 12929, 12947, 12953, 12959, 12969, 12983, 12987, - 12995, 13015, 13019, 13031, 13063, 13077, 13103, 13137, 13149, 13173, - 13207, 13211, 13227, 13241, 13249, 13255, 13269, 13283, 13285, 13303, - 13307, 13321, 13339, 13351, 13377, 13389, 13407, 13417, 13431, 13435, - 13447, 13459, 13465, 13477, 13501, 13513, 13531, 13543, 13561, 13581, - 13599, 13605, 13617, 13623, 13637, 13647, 13661, 13677, 13683, 13695, - 13725, 13729, 13753, 13773, 13781, 13785, 13795, 13801, 13807, 13825, - 13835, 13855, 13861, 13871, 13883, 13897, 13905, 13915, 13939, 13941, - 13969, 13979, 13981, 13997, 14027, 14035, 14037, 14051, 14063, 14085, - 14095, 14107, 14113, 14125, 14137, 14145, 14151, 14163, 14193, 14199, - 14219, 14229, 14233, 14243, 14277, 14287, 14289, 14295, 14301, 14305, - 14323, 14339, 14341, 14359, 14365, 14375, 14387, 14411, 14425, 14441, - 14449, 14499, 14513, 14523, 14537, 14543, 14561, 14579, 14585, 14593, - 14599, 14603, 14611, 14641, 14671, 14695, 14701, 14723, 14725, 14743, - 14753, 14759, 14765, 14795, 14797, 14803, 14831, 14839, 14845, 14855, - 14889, 14895, 14909, 14929, 14941, 14945, 14951, 14963, 14965, 14985, - 15033, 15039, 15053, 15059, 15061, 15071, 15077, 15081, 15099, 15121, - 15147, 15149, 15157, 15167, 15187, 15193, 15203, 15205, 15215, 15217, - 15223, 15243, 15257, 15269, 15273, 15287, 15291, 15313, 15335, 15347, - 15359, 15373, 15379, 15381, 15391, 15395, 15397, 15419, 15439, 15453, - 15469, 15491, 15503, 15517, 15527, 15531, 15545, 15559, 15593, 15611, - 15613, 15619, 15639, 15643, 15649, 15661, 15667, 15669, 15681, 15693, - 15717, 15721, 15741, 15745, 15765, 15793, 15799, 15811, 15825, 15835, - 15847, 15851, 15865, 15877, 15881, 15887, 15899, 15915, 15935, 15937, - 15955, 15973, 15977, 16011, 16035, 16061, 16069, 16087, 16093, 16097, - 16121, 16141, 16153, 16159, 16165, 16183, 16189, 16195, 16197, 16201, - 16209, 16215, 16225, 16259, 16265, 16273, 16299, 16309, 16355, 16375, - 16381 + 1, 3, 7, 11, 13, 19, 25, 37, 59, 47, + 61, 55, 41, 67, 97, 91, 109, 103, 115, 131, + 193, 137, 145, 143, 241, 157, 185, 167, 229, 171, + 213, 191, 253, 203, 211, 239, 247, 285, 369, 299, + 301, 333, 351, 355, 357, 361, 391, 397, 425, 451, + 463, 487, 501, 529, 539, 545, 557, 563, 601, 607, + 617, 623, 631, 637, 647, 661, 675, 677, 687, 695, + 701, 719, 721, 731, 757, 761, 787, 789, 799, 803, + 817, 827, 847, 859, 865, 875, 877, 883, 895, 901, + 911, 949, 953, 967, 971, 973, 981, 985, 995, 1001, + 1019, 1033, 1051, 1063, 1069, 1125, 1135, 1153, 1163, 1221, + 1239, 1255, 1267, 1279, 1293, 1305, 1315, 1329, 1341, 1347, + 1367, 1387, 1413, 1423, 1431, 1441, 1479, 1509, 1527, 1531, + 1555, 1557, 1573, 1591, 1603, 1615, 1627, 1657, 1663, 1673, + 1717, 1729, 1747, 1759, 1789, 1815, 1821, 1825, 1849, 1863, + 1869, 1877, 1881, 1891, 1917, 1933, 1939, 1969, 2011, 2035, + 2041, 2053, 2071, 2091, 2093, 2119, 2147, 2149, 2161, 2171, + 2189, 2197, 2207, 2217, 2225, 2255, 2257, 2273, 2279, 2283, + 2293, 2317, 2323, 2341, 2345, 2363, 2365, 2373, 2377, 2385, + 2395, 2419, 2421, 2431, 2435, 2447, 2475, 2477, 2489, 2503, + 2521, 2533, 2551, 2561, 2567, 2579, 2581, 2601, 2633, 2657, + 2669, 2681, 2687, 2693, 2705, 2717, 2727, 2731, 2739, 2741, + 2773, 2783, 2793, 2799, 2801, 2811, 2819, 2825, 2833, 2867, + 2879, 2881, 2891, 2905, 2911, 2917, 2927, 2941, 2951, 2955, + 2963, 2965, 2991, 2999, 3005, 3017, 3035, 3037, 3047, 3053, + 3083, 3085, 3097, 3103, 3159, 3169, 3179, 3187, 3205, 3209, + 3223, 3227, 3229, 3251, 3263, 3271, 3277, 3283, 3285, 3299, + 3305, 3319, 3331, 3343, 3357, 3367, 3373, 3393, 3399, 3413, + 3417, 3427, 3439, 3441, 3475, 3487, 3497, 3515, 3517, 3529, + 3543, 3547, 3553, 3559, 3573, 3589, 3613, 3617, 3623, 3627, + 3635, 3641, 3655, 3659, 3669, 3679, 3697, 3707, 3709, 3713, + 3731, 3743, 3747, 3771, 3791, 3805, 3827, 3833, 3851, 3865, + 3889, 3895, 3933, 3947, 3949, 3957, 3971, 3985, 3991, 3995, + 4007, 4013, 4021, 4045, 4051, 4069, 4073, 4179, 4201, 4219, + 4221, 4249, 4305, 4331, 4359, 4383, 4387, 4411, 4431, 4439, + 4449, 4459, 4485, 4531, 4569, 4575, 4621, 4663, 4669, 4711, + 4723, 4735, 4793, 4801, 4811, 4879, 4893, 4897, 4921, 4927, + 4941, 4977, 5017, 5027, 5033, 5127, 5169, 5175, 5199, 5213, + 5223, 5237, 5287, 5293, 5331, 5391, 5405, 5453, 5523, 5573, + 5591, 5597, 5611, 5641, 5703, 5717, 5721, 5797, 5821, 5909, + 5913, 5955, 5957, 6005, 6025, 6061, 6067, 6079, 6081, 6231, + 6237, 6289, 6295, 6329, 6383, 6427, 6453, 6465, 6501, 6523, + 6539, 6577, 6589, 6601, 6607, 6631, 6683, 6699, 6707, 6761, + 6795, 6865, 6881, 6901, 6923, 6931, 6943, 6999, 7057, 7079, + 7103, 7105, 7123, 7173, 7185, 7191, 7207, 7245, 7303, 7327, + 7333, 7355, 7365, 7369, 7375, 7411, 7431, 7459, 7491, 7505, + 7515, 7541, 7557, 7561, 7701, 7705, 7727, 7749, 7761, 7783, + 7795, 7823, 7907, 7953, 7963, 7975, 8049, 8089, 8123, 8125, + 8137, 8219, 8231, 8245, 8275, 8293, 8303, 8331, 8333, 8351, + 8357, 8367, 8379, 8381, 8387, 8393, 8417, 8435, 8461, 8469, + 8489, 8495, 8507, 8515, 8551, 8555, 8569, 8585, 8599, 8605, + 8639, 8641, 8647, 8653, 8671, 8675, 8689, 8699, 8729, 8741, + 8759, 8765, 8771, 8795, 8797, 8825, 8831, 8841, 8855, 8859, + 8883, 8895, 8909, 8943, 8951, 8955, 8965, 8999, 9003, 9031, + 9045, 9049, 9071, 9073, 9085, 9095, 9101, 9109, 9123, 9129, + 9137, 9143, 9147, 9185, 9197, 9209, 9227, 9235, 9247, 9253, + 9257, 9277, 9297, 9303, 9313, 9325, 9343, 9347, 9371, 9373, + 9397, 9407, 9409, 9415, 9419, 9443, 9481, 9495, 9501, 9505, + 9517, 9529, 9555, 9557, 9571, 9585, 9591, 9607, 9611, 9621, + 9625, 9631, 9647, 9661, 9669, 9679, 9687, 9707, 9731, 9733, + 9745, 9773, 9791, 9803, 9811, 9817, 9833, 9847, 9851, 9863, + 9875, 9881, 9905, 9911, 9917, 9923, 9963, 9973, 10003, 10025, + 10043, 10063, 10071, 10077, 10091, 10099, 10105, 10115, 10129, 10145, + 10169, 10183, 10187, 10207, 10223, 10225, 10247, 10265, 10271, 10275, + 10289, 10299, 10301, 10309, 10343, 10357, 10373, 10411, 10413, 10431, + 10445, 10453, 10463, 10467, 10473, 10491, 10505, 10511, 10513, 10523, + 10539, 10549, 10559, 10561, 10571, 10581, 10615, 10621, 10625, 10643, + 10655, 10671, 10679, 10685, 10691, 10711, 10739, 10741, 10755, 10767, + 10781, 10785, 10803, 10805, 10829, 10857, 10863, 10865, 10875, 10877, + 10917, 10921, 10929, 10949, 10967, 10971, 10987, 10995, 11009, 11029, + 11043, 11045, 11055, 11063, 11075, 11081, 11117, 11135, 11141, 11159, + 11163, 11181, 11187, 11225, 11237, 11261, 11279, 11297, 11307, 11309, + 11327, 11329, 11341, 11377, 11403, 11405, 11413, 11427, 11439, 11453, + 11461, 11473, 11479, 11489, 11495, 11499, 11533, 11545, 11561, 11567, + 11575, 11579, 11589, 11611, 11623, 11637, 11657, 11663, 11687, 11691, + 11701, 11747, 11761, 11773, 11783, 11795, 11797, 11817, 11849, 11855, + 11867, 11869, 11873, 11883, 11919, 11921, 11927, 11933, 11947, 11955, + 11961, 11999, 12027, 12029, 12037, 12041, 12049, 12055, 12095, 12097, + 12107, 12109, 12121, 12127, 12133, 12137, 12181, 12197, 12207, 12209, + 12239, 12253, 12263, 12269, 12277, 12287, 12295, 12309, 12313, 12335, + 12361, 12367, 12391, 12409, 12415, 12433, 12449, 12469, 12479, 12481, + 12499, 12505, 12517, 12527, 12549, 12559, 12597, 12615, 12621, 12639, + 12643, 12657, 12667, 12707, 12713, 12727, 12741, 12745, 12763, 12769, + 12779, 12781, 12787, 12799, 12809, 12815, 12829, 12839, 12857, 12875, + 12883, 12889, 12901, 12929, 12947, 12953, 12959, 12969, 12983, 12987, + 12995, 13015, 13019, 13031, 13063, 13077, 13103, 13137, 13149, 13173, + 13207, 13211, 13227, 13241, 13249, 13255, 13269, 13283, 13285, 13303, + 13307, 13321, 13339, 13351, 13377, 13389, 13407, 13417, 13431, 13435, + 13447, 13459, 13465, 13477, 13501, 13513, 13531, 13543, 13561, 13581, + 13599, 13605, 13617, 13623, 13637, 13647, 13661, 13677, 13683, 13695, + 13725, 13729, 13753, 13773, 13781, 13785, 13795, 13801, 13807, 13825, + 13835, 13855, 13861, 13871, 13883, 13897, 13905, 13915, 13939, 13941, + 13969, 13979, 13981, 13997, 14027, 14035, 14037, 14051, 14063, 14085, + 14095, 14107, 14113, 14125, 14137, 14145, 14151, 14163, 14193, 14199, + 14219, 14229, 14233, 14243, 14277, 14287, 14289, 14295, 14301, 14305, + 14323, 14339, 14341, 14359, 14365, 14375, 14387, 14411, 14425, 14441, + 14449, 14499, 14513, 14523, 14537, 14543, 14561, 14579, 14585, 14593, + 14599, 14603, 14611, 14641, 14671, 14695, 14701, 14723, 14725, 14743, + 14753, 14759, 14765, 14795, 14797, 14803, 14831, 14839, 14845, 14855, + 14889, 14895, 14909, 14929, 14941, 14945, 14951, 14963, 14965, 14985, + 15033, 15039, 15053, 15059, 15061, 15071, 15077, 15081, 15099, 15121, + 15147, 15149, 15157, 15167, 15187, 15193, 15203, 15205, 15215, 15217, + 15223, 15243, 15257, 15269, 15273, 15287, 15291, 15313, 15335, 15347, + 15359, 15373, 15379, 15381, 15391, 15395, 15397, 15419, 15439, 15453, + 15469, 15491, 15503, 15517, 15527, 15531, 15545, 15559, 15593, 15611, + 15613, 15619, 15639, 15643, 15649, 15661, 15667, 15669, 15681, 15693, + 15717, 15721, 15741, 15745, 15765, 15793, 15799, 15811, 15825, 15835, + 15847, 15851, 15865, 15877, 15881, 15887, 15899, 15915, 15935, 15937, + 15955, 15973, 15977, 16011, 16035, 16061, 16069, 16087, 16093, 16097, + 16121, 16141, 16153, 16159, 16165, 16183, 16189, 16195, 16197, 16201, + 16209, 16215, 16225, 16259, 16265, 16273, 16299, 16309, 16355, 16375, + 16381 }; static T2 recipd; static T1 seed_save = -1;