"slowc" option is now available to solve a simple equation in bytecode

time-shift
Ferhat Mihoubi 2013-06-10 16:45:57 +02:00
parent b257c1393a
commit b1da2f2b7d
4 changed files with 27 additions and 4 deletions

View File

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

View File

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

View File

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

View File

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