stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Major update of bytecode: 

- Iterative linear solvers using CUDA
 - interpreter.cc decomposed
time-shift
Ferhat Mihoubi 2013-03-22 15:44:34 +01:00
parent 7a8b407380
commit 03e487a092
8 changed files with 5104 additions and 2921 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
mex/build/bytecode.am
View File

@ -1,6 +1,6 @@
noinst_PROGRAMS = bytecode
bytecode_CPPFLAGS = $(AM_CPPFLAGS) -I$(top_srcdir)/../../sources/bytecode -I$(top_srcdir)/../../../preprocessor
bytecode_CPPFLAGS = $(AM_CPPFLAGS) -I$(top_srcdir)/../../sources -I$(top_srcdir)/../../sources/bytecode -I$(top_srcdir)/../../../preprocessor
TOPDIR = $(top_srcdir)/../../sources/bytecode
@ -9,8 +9,10 @@ nodist_bytecode_SOURCES = \
$(TOPDIR)/Interpreter.cc \
$(TOPDIR)/Mem_Mngr.cc \
$(TOPDIR)/SparseMatrix.cc \
$(TOPDIR)/Evaluate.cc \
$(TOPDIR)/Interpreter.hh \
$(TOPDIR)/Mem_Mngr.hh \
$(TOPDIR)/SparseMatrix.hh \
$(TOPDIR)/Evaluate.hh \
$(TOPDIR)/ErrorHandling.hh

253
mex/sources/bytecode/ErrorHandling.hh
View File

@ -23,24 +23,138 @@
#include <cstring>
#include <iostream>
#include <sstream>
#include <map>
#define BYTE_CODE
#include "CodeInterpreter.hh"
#ifdef DEBUG_EX
# include <math>
# include <math.h>
# include "mex_interface.hh"
#endif
#ifdef OCTAVE_MEX_FILE
# define CHAR_LENGTH 1
#else
# define CHAR_LENGTH 2
#endif
#ifdef _MSC_VER
#include <limits>
#define M_E 2.71828182845904523536
#define M_LOG2E 1.44269504088896340736
#define M_LOG10E 0.434294481903251827651
#define M_LN2 0.693147180559945309417
#define M_LN10 2.30258509299404568402
#define M_PI 3.14159265358979323846
#define M_PI_2 1.57079632679489661923
#define M_PI_4 0.785398163397448309616
#define M_1_PI 0.318309886183790671538
#define M_2_PI 0.636619772367581343076
#define M_1_SQRTPI 0.564189583547756286948
#define M_2_SQRTPI 1.12837916709551257390
#define M_SQRT2 1.41421356237309504880
#define M_SQRT_2 0.707106781186547524401
#define NAN numeric_limits<double>::quiet_NaN()
#define isnan(x) _isnan(x)
#define isinf(x) (!_finite(x))
#define fpu_error(x) (isinf(x) || isnan(x))
class MSVCpp_missings
{
public:
inline double
asinh(double x) const
{
if(x==0.0)
return 0.0;
double ax = abs(x);
return log(x+ax*sqrt(1.+1./(ax*ax)));
}
inline double
acosh(double x) const
{
if(x==0.0)
return 0.0;
double ax = abs(x);
return log(x+ax*sqrt(1.-1./(ax*ax)));
}
inline double
atanh(double x) const
{
return log((1+x)/(1-x))/2;
}
inline double
erf(double x) const
{
const double a1 = -1.26551223, a2 = 1.00002368,
a3 = 0.37409196, a4 = 0.09678418,
a5 = -0.18628806, a6 = 0.27886807,
a7 = -1.13520398, a8 = 1.48851587,
a9 = -0.82215223, a10 = 0.17087277;
double v = 1;
double z = abs(x);
if (z <= 0)
return v;
double t = 1 / (1 + 0.5 * z);
v = t*exp((-z*z) +a1+t*(a2+t*(a3+t*(a4+t*(a5+t*(a6+t*(a7+t*(a8+t*(a9+t*a10)))))))));
if (x < 0)
v = 2 - v;
return 1 - v;
}
inline double
nearbyint(double x) const
{
return floor(x + 0.5);
}
inline double
fmax(double x, double y) const
{
if (x > y)
return x;
else
return y;
}
inline double
fmin(double x, double y) const
{
if (x < y)
return x;
else
return y;
}
};
#endif
//#define DEBUG
using namespace std;
const int NO_ERROR_ON_EXIT = 0;
const int ERROR_ON_EXIT = 1;
typedef vector<pair<Tags, void * > > code_liste_type;
typedef code_liste_type::const_iterator it_code_type;
class GeneralExceptionHandling
{
string ErrorMsg;
public:
#ifdef _MSC_VER_
~GeneralExceptionHandling()
{
FreeLibrary(hinstLib);
};
#endif
GeneralExceptionHandling(string ErrorMsg_arg) : ErrorMsg(ErrorMsg_arg)
{
};
@ -121,6 +235,16 @@ public:
};
};
class UserExceptionHandling : public GeneralExceptionHandling
{
double value;
public:
UserExceptionHandling() : GeneralExceptionHandling("Fatal error in bytecode:")
{
completeErrorMsg(" User break\n");
};
};
class FatalExceptionHandling : public GeneralExceptionHandling
{
public:
@ -133,16 +257,40 @@ public:
};
};
class ErrorMsg
struct s_plan
{
string var, exo;
int var_num, exo_num;
vector<pair<int, double> > per_value;
};
#ifdef MATLAB_MEX_FILE
extern "C" bool utIsInterruptPending();
#else
#include <octave/oct.h>
#include <octave/unwind-prot.h>
#endif
#ifdef _MSC_VER
class ErrorMsg : public MSVCpp_missings
#else
class ErrorMsg
#endif
{
private:
bool is_load_variable_list;
public:
double *y, *ya;
int y_size;
double *T;
int nb_row_xd, nb_row_x, y_size;
int nb_row_xd, nb_row_x;
int y_kmin, y_kmax, periods;
double *x, *params;
double *u, *y, *ya;
double *u;
double *steady_y, *steady_x;
double *g2, *g1, *r;
double *g2, *g1, *r, *res;
vector<s_plan> splan, spfplan;
vector<mxArray *> jacobian_block, jacobian_other_endo_block, jacobian_exo_block, jacobian_det_exo_block;
map<unsigned int, double> TEF;
map<pair<unsigned int, unsigned int>, double > TEFD;
@ -150,11 +298,12 @@ public:
ExpressionType EQN_type;
it_code_type it_code_expr;
unsigned int nb_endo, nb_exo, nb_param;
/*unsigned int*/size_t nb_endo, nb_exo, nb_param;
char *P_endo_names, *P_exo_names, *P_param_names;
unsigned int endo_name_length, exo_name_length, param_name_length;
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<pair<string, pair<SymbolType, unsigned int> > > Variable_list;
inline
ErrorMsg()
@ -169,6 +318,7 @@ public:
nb_param = mxGetM(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "param_names")));
param_name_length = mxGetN(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "param_names")));
P_param_names = (char *) mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "param_names")));
is_load_variable_list = false;
}
inline string
@ -184,9 +334,9 @@ public:
else
{
if (str[i] == '$')
pos1 = temp.length();
pos1 = int(temp.length());
else
pos2 = temp.length();
pos2 = int(temp.length());
if (pos1 >= 0 && pos2 >= 0)
{
tmp_n.erase(pos1, pos2-pos1+1);
@ -199,6 +349,50 @@ public:
return temp;
}
inline void
load_variable_list()
{
ostringstream res;
for (unsigned int variable_num = 0; variable_num < (unsigned int)nb_endo; variable_num++)
{
for (unsigned int i = 0; i < endo_name_length; i++)
if (P_endo_names[CHAR_LENGTH*(variable_num+i*nb_endo)] != ' ')
res << P_endo_names[CHAR_LENGTH*(variable_num+i*nb_endo)];
Variable_list.push_back(make_pair(res.str(), make_pair(eEndogenous, variable_num)));
}
for (unsigned int variable_num = 0; variable_num < (unsigned int)nb_exo; variable_num++)
{
for (unsigned int i = 0; i < exo_name_length; i++)
if (P_exo_names[CHAR_LENGTH*(variable_num+i*nb_exo)] != ' ')
res << P_exo_names[CHAR_LENGTH*(variable_num+i*nb_exo)];
Variable_list.push_back(make_pair(res.str(), make_pair(eExogenous, variable_num)));
}
}
inline int
get_ID(const string variable_name, SymbolType *variable_type)
{
if (!is_load_variable_list)
{
load_variable_list();
is_load_variable_list = true;
}
size_t n = Variable_list.size();
int i = 0;
bool notfound = true;
while (notfound && i < n)
{
if (variable_name == Variable_list[i].first)
{
notfound = false;
*variable_type = Variable_list[i].second.first;
return Variable_list[i].second.second;
}
i++;
}
return(-1);
}
inline string
get_variable(const SymbolType variable_type, const unsigned int variable_num) const
{
@ -293,7 +487,6 @@ public:
break;
default:
return ("???");
break;
}
else
switch (EQN_type)
@ -342,7 +535,6 @@ public:
break;
default:
return ("???");
break;
}
it_code_type it_code_ret;
Error_loc << endl << add_underscore_to_fpe(" " + print_expression(it_code_expr, evaluate, size, block_num, steady_state, Per_u_, it_, it_code_ret, true));
@ -378,6 +570,12 @@ public:
while (go_on)
{
#ifdef OCTAVE_MEX_FILE
OCTAVE_QUIT;
#else
if ( utIsInterruptPending() )
throw UserExceptionHandling();
#endif
switch (it_code->first)
{
case FNUMEXPR:
@ -441,7 +639,9 @@ public:
case eParameter:
var = ((FLDV_ *) it_code->second)->get_pos();
#ifdef DEBUG
mexPrintf("FLDV_ Param var=%d", var);
mexPrintf("FLDV_ Param var=%d\n", var);
mexPrintf("get_variable(eParameter, var)=%s\n",get_variable(eParameter, var).c_str());
mexEvalString("drawnow;");
#endif
Stack.push(get_variable(eParameter, var));
if (compute)
@ -451,7 +651,10 @@ public:
var = ((FLDV_ *) it_code->second)->get_pos();
lag = ((FLDV_ *) it_code->second)->get_lead_lag();
#ifdef DEBUG
mexPrintf("FLDV_ endo var=%d, lag=%d", var, lag);
mexPrintf("FLDV_ endo var=%d, lag=%d\n", var, lag);
mexPrintf("get_variable(eEndogenous, var)=%s, compute=%d\n",get_variable(eEndogenous, var).c_str(), compute);
mexPrintf("it_=%d, lag=%d, y_size=%d, var=%d, y=%x\n", it_, lag, y_size, var, y);
mexEvalString("drawnow;");
#endif
tmp_out.str("");
if (lag > 0)
@ -1250,7 +1453,7 @@ public:
Stack.pop();
if (compute)
{
int derivOrder = nearbyint(Stackf.top());
int derivOrder = int(nearbyint(Stackf.top()));
Stackf.pop();
if (fabs(v1f) < NEAR_ZERO && v2f > 0
&& derivOrder > v2f
@ -1570,7 +1773,11 @@ public:
}
tmp_out.str("");
tmp_out << function_name << "(";
#ifndef _MSC_VER
string ss[nb_input_arguments];
#else
vector<string> ss(nb_input_arguments);
#endif
for (unsigned int i = 0; i < nb_input_arguments; i++)
{
ss[nb_input_arguments-i-1] = Stack.top();
@ -1624,7 +1831,11 @@ public:
tmp_out.str("");
tmp_out << function_name << "(";
tmp_out << arg_func_name.c_str() << ", " << fc->get_row() << ", {";
#ifndef _MSC_VER
string ss[nb_add_input_arguments];
#else
vector<string> ss(nb_input_arguments);
#endif
for (unsigned int i = 0; i < nb_add_input_arguments; i++)
{
ss[nb_add_input_arguments-i-1] = Stack.top();
@ -1655,7 +1866,11 @@ public:
}
tmp_out.str("");
tmp_out << function_name << "(";
#ifndef _MSC_VER
string ss[nb_input_arguments];
#else
vector<string> ss(nb_input_arguments);
#endif
for (unsigned int i = 0; i < nb_input_arguments; i++)
{
ss[nb_input_arguments-i-1] = Stack.top();
@ -1708,7 +1923,11 @@ public:
tmp_out.str("");
tmp_out << function_name << "(";
tmp_out << arg_func_name.c_str() << ", " << fc->get_row() << ", " << fc->get_col() << ", {";
#ifndef _MSC_VER
string ss[nb_add_input_arguments];
#else
vector<string> ss(nb_input_arguments);
#endif
for (unsigned int i = 0; i < nb_add_input_arguments; i++)
{
ss[nb_add_input_arguments-i-1] = Stack.top();
@ -1739,7 +1958,11 @@ public:
}
tmp_out.str("");
tmp_out << function_name << "(";
#ifndef _MSC_VER
string ss[nb_input_arguments];
#else
vector<string> ss(nb_input_arguments);
#endif
for (unsigned int i = 0; i < nb_input_arguments; i++)
{
ss[nb_input_arguments-i-1] = Stack.top();
@ -1965,7 +2188,7 @@ public:
it_code++;
}
#ifdef DEBUG
mexPrintf("print_expression end\n"); mexEvalString("drawnow;");
mexPrintf("print_expression end tmp_out.str().c_str()=%s\n", tmp_out.str().c_str()); mexEvalString("drawnow;");
#endif
it_code_ret = it_code;
return (tmp_out.str());

2436
mex/sources/bytecode/Interpreter.cc
View File

File diff suppressed because it is too large Load Diff

54
mex/sources/bytecode/Interpreter.hh
View File

@ -27,6 +27,7 @@
#define BYTE_CODE
#include "CodeInterpreter.hh"
#include "SparseMatrix.hh"
#include "Evaluate.hh"
#ifdef LINBCG
# include "linbcg.hh"
#endif
@ -40,50 +41,29 @@
using namespace std;
#define pow_ pow
class Interpreter : public SparseMatrix
class Interpreter : public dynSparseMatrix
{
private:
unsigned int EQN_dvar1, EQN_dvar2, EQN_dvar3;
int EQN_lag1, EQN_lag2, EQN_lag3;
mxArray *GlobalTemporaryTerms;
protected:
double pow1(double a, double b);
double divide(double a, double b);
double log1(double a);
double log10_1(double a);
void compute_block_time(int Per_u_, bool evaluate, int block_num, int size, bool steady_state);
void evaluate_a_block(const int size, const int type, string bin_basename, bool steady_state, int block_num,
const bool is_linear = false, const int symbol_table_endo_nbr = 0, const int Block_List_Max_Lag = 0, const int Block_List_Max_Lead = 0, const int u_count_int = 0, int block = -1);
int 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 = false, const int symbol_table_endo_nbr = 0, const int Block_List_Max_Lag = 0, const int Block_List_Max_Lead = 0, const int u_count_int = 0);
void print_a_block(const int size, const int type, string bin_basename, bool steady_state, int block_num,
const bool is_linear, const int symbol_table_endo_nbr, const int Block_List_Max_Lag,
const int Block_List_Max_Lead, const int u_count_int, int block);
void SingularDisplay(int Per_u_, bool evaluate, int Block_Count, int size, bool steady_state, it_code_type begining);
vector<Block_contain_type> Block_Contain;
code_liste_type code_liste;
it_code_type it_code;
int Block_Count, Per_u_, Per_y_;
int it_, maxit_, size_of_direction;
double solve_tolf;
bool GaussSeidel;
map<pair<pair<int, int>, int>, int> IM_i;
int equation, derivative_equation, derivative_variable;
string filename;
int minimal_solving_periods;
int stack_solve_algo, solve_algo;
bool global_temporary_terms;
bool print, print_error;
void evaluate_a_block();
int simulate_a_block();
void print_a_block();
public:
~Interpreter();
Interpreter(double *params_arg, double *y_arg, double *ya_arg, double *x_arg, double *steady_y_arg, double *steady_x_arg,
double *direction_arg, int y_size_arg, int nb_row_x_arg,
int nb_row_xd_arg, int periods_arg, int y_kmin_arg, int y_kmax_arg, int maxit_arg_, double solve_tolf_arg, int size_o_direction_arg,
double slowc_arg, int y_decal_arg, double markowitz_c_arg, string &filename_arg, int minimal_solving_periods_arg, int stack_solve_algo_arg, int solve_algo_arg,
bool global_temporary_terms_arg, bool print_arg, bool print_error_arg, mxArray *GlobalTemporaryTerms_arg);
bool compute_blocks(string file_name, string bin_basename, bool steady_state, bool evaluate, int block, int &nb_blocks, bool print_it);
double *direction_arg, size_t y_size_arg,
size_t nb_row_x_arg, size_t nb_row_xd_arg, int periods_arg, int y_kmin_arg, int y_kmax_arg,
int maxit_arg_, double solve_tolf_arg, size_t size_of_direction_arg, double slowc_arg, int y_decal_arg, double markowitz_c_arg,
string &filename_arg, int minimal_solving_periods_arg, int stack_solve_algo_arg, int solve_algo_arg,
bool global_temporary_terms_arg, bool print_arg, bool print_error_arg, mxArray *GlobalTemporaryTerms_arg,
bool steady_state_arg, bool print_it_arg
#ifdef CUDA
, const int CUDA_device, cublasHandle_t cublas_handle_arg, cusparseHandle_t cusparse_handle_arg, cusparseMatDescr_t descr_arg
#endif
);
bool compute_blocks(string file_name, string bin_basename, bool evaluate, int block, int &nb_blocks);
inline mxArray *
get_jacob(int block_num)
{

4474
mex/sources/bytecode/SparseMatrix.cc
View File

File diff suppressed because it is too large Load Diff

133
mex/sources/bytecode/SparseMatrix.hh
View File

@ -19,21 +19,62 @@
#ifndef SPARSEMATRIX_HH_INCLUDED
#define SPARSEMATRIX_HH_INCLUDED
#define PRINTF_ printf
#include <fstream>
#include <stack>
#include <cmath>
#include <map>
#include <ctime>
#include "dynblas.h"
#if !(defined _MSC_VER)
#include "dynumfpack.h"
#endif
#ifdef OCTAVE_MEX_FILE
# define CHAR_LENGTH 1
#else
# define CHAR_LENGTH 2
#ifdef CUDA
#include "cuda.h"
#include "cuda_runtime_api.h"
#include "cublas_v2.h"
#include "cusparse_v2.h"
#endif
#include "Mem_Mngr.hh"
#include "ErrorHandling.hh"
//#include "Interpreter.hh"
#include "Evaluate.hh"
#define cudaChk(x, y) \
{ \
cudaError_t cuda_error = x; \
if (cuda_error != cudaSuccess) \
{ \
ostringstream tmp; \
tmp << y; \
throw FatalExceptionHandling(tmp.str()); \
} \
};
#define cusparseChk(x, y) \
{ \
cusparseStatus_t cusparse_status = x; \
if (cusparse_status != CUSPARSE_STATUS_SUCCESS) \
{ \
ostringstream tmp; \
tmp << y; \
throw FatalExceptionHandling(tmp.str()); \
} \
};
#define cublasChk(x, y) \
{ \
cublasStatus_t cublas_status = x; \
if (cublas_status != CUBLAS_STATUS_SUCCESS) \
{ \
ostringstream tmp; \
tmp << y; \
throw FatalExceptionHandling(tmp.str()); \
} \
};
#define NEW_ALLOC
#define MARKOVITZ
@ -53,41 +94,76 @@ const int IFLDZ = 4;
const int IFMUL = 5;
const int IFSTP = 6;
const int IFADD = 7;
const double eps = 1e-10;
const double eps = 1e-15;
const double very_big = 1e24;
const int alt_symbolic_count_max = 1;
const double mem_increasing_factor = 1.1;
class SparseMatrix : public ErrorMsg
class dynSparseMatrix : public Evaluate
{
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)*/;
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);
#if (defined _MSC_VER)
typedef int64_t SuiteSparse_long;
#endif
dynSparseMatrix();
dynSparseMatrix(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
#ifdef CUDA
,const int CUDA_device_arg, cublasHandle_t cublas_handle_arg, cusparseHandle_t cusparse_handle_arg, cusparseMatDescr_t descr_arg
#endif
);
void Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin, int y_kmax, int Size, int periods, bool cvg, int minimal_solving_periods, int stack_solve_algo, unsigned int endo_name_length, char *P_endo_names);
void Simulate_Newton_One_Boundary(bool forward);
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_SparseMatrix(string file_name, const int Size, int periods, int y_kmin, int y_kmax, 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;
void Singular_display(int block, int Size);
void End_Solver();
double g0, gp0, glambda2;
int try_at_iteration;
private:
void Init_GE(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM);
void Init_Matlab_Sparse(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM, mxArray *A_m, mxArray *b_m, mxArray *x0_m);
void Init_UMFPACK_Sparse(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM, SuiteSparse_long **Ap, SuiteSparse_long **Ai, double **Ax, double **b, mxArray *x0_m);
#ifdef CUDA
void Init_CUDA_Sparse(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM, int **Ap, int **Ai, double **Ax, int **Ap_tild, int **Ai_tild, double **A_tild, double **b, double **x0, mxArray *x0_m, int *nnz, int *nnz_tild, int preconditioner);
#endif
void Init_Matlab_Sparse_Simple(int Size, map<pair<pair<int, int>, int>, int> &IM, mxArray *A_m, mxArray *b_m, bool &zero_solution, mxArray *x0_m);
void Init_UMFPACK_Sparse_Simple(int Size, map<pair<pair<int, int>, int>, int> &IM, SuiteSparse_long **Ap, SuiteSparse_long **Ai, double **Ax, double **b, bool &zero_solution, mxArray *x0_m);
void Init_CUDA_Sparse_Simple(int Size, map<pair<pair<int, int>, int>, int> &IM, SuiteSparse_long **Ap, SuiteSparse_long **Ai, double **Ax, double **b, double **x0, bool &zero_solution, mxArray *x0_m);
void Simple_Init(int Size, std::map<std::pair<std::pair<int, int>, int>, int> &IM, bool &zero_solution);
void End_GE(int Size);
bool mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb, double *fc);
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, bool steady_state, int it_);
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_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);
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, bool steady_state);
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_);
void End_Matlab_LU_UMFPack();
#ifdef CUDA
void 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,
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);
#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);
void Check_and_Correct_Previous_Iteration(int block_num, int y_size, int size, double crit_opt_old);
bool Simulate_One_Boundary(int blck, int y_size, int y_kmin, int y_kmax, int Size, bool cvg);
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 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
);
void Grad_f_product(int n, mxArray *b_m, double* vectr, mxArray *A_m, SuiteSparse_long *Ap, SuiteSparse_long *Ai, double* Ax, double *b);
void Insert(const int r, const int c, const int u_index, const int lag_index);
void Delete(const int r, const int c);
int At_Row(int r, NonZeroElem **first);
@ -102,7 +178,8 @@ private:
void Delete_u(int pos);
void Clear_u();
void Print_u();
void CheckIt(int y_size, int y_kmin, int y_kmax, int Size, int periods, int iter);
void *Symbolic, *Numeric ;
void CheckIt(int y_size, int y_kmin, int y_kmax, int Size, int periods);
void Check_the_Solution(int periods, int y_kmin, int y_kmax, int Size, double *u, int *pivot, int *b);
int complete(int beg_t, int Size, int periods, int *b);
void bksub(int tbreak, int last_period, int Size, double slowc_l
@ -118,12 +195,18 @@ private:
mxArray *Sparse_substract_SA_SB(mxArray *A_m, mxArray *B_m);
mxArray *Sparse_substract_A_SB(mxArray *A_m, mxArray *B_m);
mxArray *substract_A_B(mxArray *A_m, mxArray *B_m);
#ifdef CUDA
int CUDA_device;
cublasHandle_t cublas_handle;
cusparseHandle_t cusparse_handle;
cusparseMatDescr_t CUDA_descr;
#endif
protected:
stack<double> Stack;
int nb_prologue_table_u, nb_first_table_u, nb_middle_table_u, nb_last_table_u;
int nb_prologue_table_y, nb_first_table_y, nb_middle_table_y, nb_last_table_y;
int middle_count_loop;
char type;
//char type;
fstream SaveCode;
string filename;
int max_u, min_u;
@ -154,19 +237,19 @@ protected:
int u_count_alloc, u_count_alloc_save;
vector<double *> jac;
double *jcb;
double res1, res2, max_res;
int max_res_idx;
double slowc, slowc_save, prev_slowc_save, markowitz_c;
int y_decal;
int *index_vara, *index_equa;
int *index_equa;
int u_count, tbreak_g;
int iter;
double *direction;
int start_compare;
int restart;
bool error_not_printed;
double g_lambda1, g_lambda2, gp_0;
double lu_inc_tol;
//private:
SuiteSparse_long *Ap_save, *Ai_save;
double *Ax_save, *b_save;
mxArray *A_m_save, *b_m_save;
};
#endif

669
mex/sources/bytecode/bytecode.cc
View File

@ -18,12 +18,18 @@
*/
#include <cstring>
#include "Interpreter.hh"
#include "ErrorHandling.hh"
#include <ctime>
#include <math.h>
#ifdef DEBUG_EX
using namespace std;
# include <sstream>
string
Get_Argument(const char *argv)
{
@ -33,14 +39,17 @@ Get_Argument(const char *argv)
#else
void (*prev_fn)(int);
string
Get_Argument(const mxArray *prhs)
{
const mxArray *mxa = prhs;
int buflen = mxGetM(mxa) * mxGetN(mxa) + 1;
mwSize buflen = mwSize(mxGetM(mxa) * mxGetN(mxa) + 1);
char *first_argument;
first_argument = (char *) mxCalloc(buflen, sizeof(char));
int status = mxGetString(mxa, first_argument, buflen);
size_t status = mxGetString(mxa, first_argument, buflen);
if (status != 0)
mexWarnMsgTxt("Not enough space. The first argument is truncated.");
string f(first_argument);
@ -49,6 +58,178 @@ Get_Argument(const mxArray *prhs)
}
#endif
//#include <windows.h>
#include <stdio.h>
#ifdef CUDA
int
GPU_Test_and_Info(cublasHandle_t *cublas_handle, cusparseHandle_t *cusparse_handle, cusparseMatDescr_t *descr)
{
cudaDeviceProp deviceProp;
int device_count, device, version, version_max = 0;
cublasStatus_t cublas_status;
cudaError_t cuda_error;
*descr=0;
/* ask cuda how many devices it can find */
cudaGetDeviceCount(&device_count);
if (device_count < 1)
{
/* if it couldn't find any fail out */
ostringstream tmp;
tmp << " Unable to find a CUDA device. Unable to implement CUDA solvers\n";
throw FatalExceptionHandling(tmp.str());
}
else
{
mexPrintf("-----------------------------------------\n");
for (int i = 0; i < device_count; i++)
{
cudaSetDevice(i);
// Statistics about the GPU device
cuda_error = cudaGetDeviceProperties(&deviceProp, i);
if (cuda_error != cudaSuccess)
{
ostringstream tmp;
tmp << " bytecode cudaGetDeviceProperties failed\n";
throw FatalExceptionHandling(tmp.str());
}
mexPrintf("> GPU device %d: \"%s\" has:\n - %d Multi-Processors,\n - %d threads per multiprocessor,\n", i, deviceProp.name, deviceProp.multiProcessorCount, deviceProp.maxThreadsPerMultiProcessor);
mexEvalString("drawnow;");
version = (deviceProp.major * 0x10 + deviceProp.minor);
if (version >= version_max)
{
device = i;
version_max = version;
}
mexPrintf(" - %4.2fMhz clock rate,\n - %2.0fMb of memory,\n - %d.%d compute capabilities.\n", double(deviceProp.clockRate) / (1024 * 1024), double(deviceProp.totalGlobalMem) / (1024 * 1024), deviceProp.major, deviceProp.minor);
mexEvalString("drawnow;");
}
}
mexPrintf("> Device %d selected\n", device);
mexEvalString("drawnow;");
cuda_error = cudaSetDevice(device);
if (cuda_error != cudaSuccess)
{
ostringstream tmp;
tmp << " bytecode cudaSetDevice failed\n";
throw FatalExceptionHandling(tmp.str());
}
if(version_max < 0x11)
{
ostringstream tmp;
tmp << " bytecode requires a minimum CUDA compute 1.1 capability\n";
cudaDeviceReset();
throw FatalExceptionHandling(tmp.str());
}
// Initialize CuBlas library
cublas_status = cublasCreate(cublas_handle);
if (cublas_status != CUBLAS_STATUS_SUCCESS)
{
ostringstream tmp;
switch(cublas_status)
{
case CUBLAS_STATUS_NOT_INITIALIZED:
tmp << " the CUBLAS initialization failed.\n";
break;
case CUBLAS_STATUS_ALLOC_FAILED:
tmp << " the resources could not be allocated.\n";
break;
default:
tmp << " unknown error during the initialization of cusparse library.\n";
}
throw FatalExceptionHandling(tmp.str());
}
// Initialize the CuSparse library
cusparseStatus_t cusparse_status;
cusparse_status = cusparseCreate(cusparse_handle);
if (cusparse_status != CUSPARSE_STATUS_SUCCESS)
{
ostringstream tmp;
switch(cusparse_status)
{
case CUSPARSE_STATUS_NOT_INITIALIZED:
tmp << " the CUDA Runtime initialization failed.\n";
break;
case CUSPARSE_STATUS_ALLOC_FAILED:
tmp << " the resources could not be allocated.\n";
break;
case CUSPARSE_STATUS_ARCH_MISMATCH:
tmp << " the device compute capability (CC) is less than 1.1. The CC of at least 1.1 is required.\n";
break;
default:
tmp << " unknown error during the initialization of cusparse library.\n";
}
throw FatalExceptionHandling(tmp.str());
}
// Create and setup matrix descriptor
cusparse_status = cusparseCreateMatDescr(descr);
if (cusparse_status != CUSPARSE_STATUS_SUCCESS)
{
ostringstream tmp;
tmp << " Matrix descriptor initialization failed\n";
throw FatalExceptionHandling(tmp.str());
}
cusparseSetMatType(*descr, CUSPARSE_MATRIX_TYPE_GENERAL);
cusparseSetMatIndexBase(*descr, CUSPARSE_INDEX_BASE_ZERO);
mexPrintf("> Driver version:\n");
int cuda_version;
cuda_error = cudaDriverGetVersion(&cuda_version);
if (cuda_error != cudaSuccess)
{
ostringstream tmp;
tmp << " cudaGetVersion has failed\n";
throw FatalExceptionHandling(tmp.str());
}
mexPrintf(" - CUDA version %5.3f\n", double(cuda_version) / 1000);
int cublas_version;
cublas_status = cublasGetVersion(*cublas_handle, &cublas_version);
if (cublas_status != CUBLAS_STATUS_SUCCESS)
{
ostringstream tmp;
tmp << " cublasGetVersion has failed\n";
throw FatalExceptionHandling(tmp.str());
}
mexPrintf(" - CUBLAS version %5.3f\n", double(cublas_version) / 1000);
int cusparse_version;
cusparse_status = cusparseGetVersion(*cusparse_handle, &cusparse_version);
if (cusparse_status != CUSPARSE_STATUS_SUCCESS)
{
ostringstream tmp;
tmp << " cusparseGetVersion has failed\n";
throw FatalExceptionHandling(tmp.str());
}
mexPrintf(" - CUSPARSE version %5.3f\n", double(cusparse_version) / 1000);
mexPrintf("-----------------------------------------\n");
return device;
}
void
GPU_close(cublasHandle_t cublas_handle, cusparseHandle_t cusparse_handle, cusparseMatDescr_t descr)
{
cublasChk(cublasDestroy(cublas_handle),"in bytecode cublasDestroy failed\n");
cusparseChk(cusparseDestroyMatDescr(descr), "in bytecode cusparseDestroyMatDescr failed\n");
cusparseChk(cusparseDestroy(cusparse_handle),"in bytecode cusparseDestroy failed\n");
}
#endif
string
deblank(string x)
{
for(int i = 0; i < x.length(); i++)
if (x[i] == ' ')
x.erase(i--, 1);
return x;
}
void
Get_Arguments_and_global_variables(int nrhs,
#ifndef DEBUG_EX
@ -57,10 +238,10 @@ Get_Arguments_and_global_variables(int nrhs,
const char *prhs[],
#endif
int &count_array_argument,
double *yd[], unsigned int &row_y, unsigned int &col_y,
double *xd[], unsigned int &row_x, unsigned int &col_x,
double *params[],
double *steady_yd[], unsigned int &steady_row_y, unsigned int &steady_col_y,
double *yd[], size_t &row_y, size_t &col_y,
double *xd[], size_t &row_x, size_t &col_x,
double *params[],
double *steady_yd[], size_t &steady_row_y, size_t &steady_col_y,
unsigned int &periods,
#ifndef DEBUG_EX
mxArray *block_structur[],
@ -69,8 +250,10 @@ Get_Arguments_and_global_variables(int nrhs,
mxArray *M_[], mxArray *oo_[], mxArray *options_[], bool &global_temporary_terms,
bool &print,
bool &print_error,
mxArray *GlobalTemporaryTerms[])
mxArray *GlobalTemporaryTerms[],
string *plan_struct_name, string *pfplan_struct_name)
{
size_t pos;
#ifdef DEBUG_EX
for (int i = 2; i < nrhs; i++)
#else
@ -101,7 +284,7 @@ Get_Arguments_and_global_variables(int nrhs,
steady_col_y = mxGetN(prhs[i]);
break;
case 4:
periods = mxGetScalar(prhs[i]);
periods = int(mxGetScalar(prhs[i]));
break;
case 5:
*block_structur = mxDuplicateArray(prhs[i]);
@ -111,7 +294,7 @@ Get_Arguments_and_global_variables(int nrhs,
*GlobalTemporaryTerms = mxDuplicateArray(prhs[i]);
break;
default:
//mexPrintf("Unknown argument count_array_argument=%d\n",count_array_argument);
mexPrintf("Unknown argument count_array_argument=%d\n",count_array_argument);
break;
}
count_array_argument++;
@ -132,16 +315,34 @@ Get_Arguments_and_global_variables(int nrhs,
print_error = false;
else
{
int pos = Get_Argument(prhs[i]).find("block");
if (pos != (int) string::npos)
;
if ((pos = Get_Argument(prhs[i]).find("block")) != (int) string::npos)
{
int pos1 = Get_Argument(prhs[i]).find("=", pos+5);
size_t pos1 = Get_Argument(prhs[i]).find("=", pos+5);
if (pos1 != (int) string::npos)
pos = pos1 + 1;
else
pos += 5;
block = atoi(Get_Argument(prhs[i]).substr(pos, string::npos).c_str())-1;
}
else if ((pos = Get_Argument(prhs[i]).find("pfplan")) != (int) string::npos)
{
size_t pos1 = Get_Argument(prhs[i]).find("=", pos+6);
if (pos1 != (int) string::npos)
pos = pos1 + 1;
else
pos += 6;
*pfplan_struct_name = deblank(Get_Argument(prhs[i]).substr(pos, string::npos));
}
else if ((pos = Get_Argument(prhs[i]).find("plan")) != (int) string::npos)
{
size_t pos1 = Get_Argument(prhs[i]).find("=", pos+4);
if (pos1 != (int) string::npos)
pos = pos1 + 1;
else
pos += 4;
*plan_struct_name = deblank(Get_Argument(prhs[i]).substr(pos, string::npos));
}
else
{
ostringstream tmp;
@ -185,6 +386,7 @@ Get_Arguments_and_global_variables(int nrhs,
}
}
#ifdef DEBUG_EX
int
main(int nrhs, const char *prhs[])
@ -203,9 +405,9 @@ main(int nrhs, const char *prhs[])
char *plhs[1];
load_global((char *) prhs[1]);
#endif
//ErrorHandlingException error_handling;
unsigned int i, row_y = 0, col_y = 0, row_x = 0, col_x = 0, nb_row_xd = 0;
unsigned int steady_row_y, steady_col_y;
mxArray *plan_struct = NULL, *pfplan_struct = NULL;
size_t i, row_y = 0, col_y = 0, row_x = 0, col_x = 0, nb_row_xd = 0;
size_t steady_row_y, steady_col_y;
int y_kmin = 0, y_kmax = 0, y_decal = 0;
unsigned int periods = 1;
double *direction;
@ -218,13 +420,22 @@ main(int nrhs, const char *prhs[])
bool global_temporary_terms = false;
bool print = false, print_error = true, print_it = false;
double *steady_yd = NULL, *steady_xd = NULL;
string plan, pfplan;
vector<s_plan> splan, spfplan;
#ifdef CUDA
int CUDA_device = -1;
cublasHandle_t cublas_handle;
cusparseHandle_t cusparse_handle;
cusparseMatDescr_t descr;
#endif
try
{
Get_Arguments_and_global_variables(nrhs, prhs, count_array_argument,
&yd, row_y, col_y,
&xd, row_x, col_x,
&params,
&params,
&steady_yd, steady_row_y, steady_col_y,
periods,
#ifndef DEBUG_EX
@ -232,123 +443,401 @@ main(int nrhs, const char *prhs[])
#endif
steady_state, evaluate, block,
&M_, &oo_, &options_, global_temporary_terms,
print, print_error, &GlobalTemporaryTerms);
print, print_error, &GlobalTemporaryTerms,
&plan, &pfplan);
}
catch (GeneralExceptionHandling &feh)
{
DYN_MEX_FUNC_ERR_MSG_TXT(feh.GetErrorMsg().c_str());
}
if (!count_array_argument)
params = mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "params")));
{
int field = mxGetFieldNumber(M_, "params");
if (field < 0)
DYN_MEX_FUNC_ERR_MSG_TXT("params is not a field of M_");
params = mxGetPr(mxGetFieldByNumber(M_, 0, field));
}
ErrorMsg emsg;
if (plan.length()>0)
{
mxArray* plan_struct = mexGetVariable("base", plan.c_str());
if (plan_struct == NULL)
{
string tmp = plan;
tmp.insert(0,"Can't find the plan: ");
DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
}
size_t n_plan = mxGetN(plan_struct);
splan.resize(n_plan);
for (int i = 0; i < n_plan; i++)
{
splan[i].var = "";
splan[i].exo = "";
mxArray* tmp = mxGetField(plan_struct, i, "exo");
if (tmp)
{
char name [100];
mxGetString(tmp, name, 100);
splan[i].var = name;
SymbolType variable_type;
int exo_num = emsg.get_ID(name, &variable_type);
if (variable_type == eExogenous || variable_type == eExogenousDet)
splan[i].var_num = exo_num;
else
{
string tmp = name;
tmp.insert(0,"the variable '");
tmp.append("' defined as var in plan is not an exogenous or a deterministic exogenous\n");
DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
}
}
tmp = mxGetField(plan_struct, i, "var");
if (tmp)
{
char name [100];
mxGetString(tmp, name, 100);
splan[i].exo = name;
SymbolType variable_type;
int exo_num = emsg.get_ID(name, &variable_type);
if (variable_type == eEndogenous)
splan[i].exo_num = exo_num;
else
{
string tmp = name;
tmp.insert(0,"the variable '");
tmp.append("' defined as exo in plan is not an endogenous variable\n");
DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
}
}
tmp = mxGetField(plan_struct, i, "per_value");
if (tmp)
{
size_t num_shocks = mxGetM(tmp);
(splan[i]).per_value.resize(num_shocks);
double * per_value = mxGetPr(tmp);
for (int j = 0; j < num_shocks; j++)
(splan[i]).per_value[j] = make_pair(ceil(per_value[j]), per_value[j + num_shocks]);
}
}
int i;
for (vector<s_plan>::iterator it = splan.begin(); it != splan.end(); it++)
{
mexPrintf("----------------------------------------------------------------------------------------------------\n");
mexPrintf("suprise n°%d\n", i+1);
if (it->exo.length())
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<pair<int, double> >::iterator it1 = it->per_value.begin(); it1 != it->per_value.end(); it1++)
{
mexPrintf(" %3d %10.5f\n",it1->first, it1->second);
}
i++;
}
}
if (pfplan.length()>0)
{
pfplan_struct = mexGetVariable("base", pfplan.c_str());
if (!pfplan_struct)
{
string tmp = pfplan;
tmp.insert(0,"Can't find the pfplan: ");
DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
}
size_t n_plan = mxGetN(pfplan_struct);
spfplan.resize(n_plan);
for (int i = 0; i < n_plan; i++)
{
spfplan[i].var = "";
spfplan[i].exo = "";
mxArray* tmp = mxGetField(pfplan_struct, i, "var");
if (tmp)
{
char name [100];
mxGetString(tmp, name, 100);
spfplan[i].var = name;
SymbolType variable_type;
int exo_num = emsg.get_ID(name, &variable_type);
if (variable_type == eExogenous || variable_type == eExogenousDet)
splan[i].var_num = exo_num;
else
{
string tmp = name;
tmp.insert(0,"the variable '");
tmp.append("' defined as var in pfplan is not an exogenous or a deterministic exogenous\n");
DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
}
}
tmp = mxGetField(pfplan_struct, i, "exo");
if (tmp)
{
char name [100];
mxGetString(tmp, name, 100);
spfplan[i].exo = name;
SymbolType variable_type;
int exo_num = emsg.get_ID(name, &variable_type);
if (variable_type == eEndogenous)
spfplan[i].exo_num = exo_num;
else
{
string tmp = name;
tmp.insert(0,"the variable '");
tmp.append("' defined as exo in pfplan is not an endogenous variable\n");
DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
}
}
tmp = mxGetField(pfplan_struct, i, "per_value");
if (tmp)
{
size_t num_shocks = mxGetM(tmp);
double * per_value = mxGetPr(tmp);
(spfplan[i]).per_value.resize(num_shocks);
for (int j = 0; j < num_shocks; j++)
spfplan[i].per_value[j] = make_pair(ceil(per_value[j]), per_value[j+ num_shocks]);
}
}
int i;
for (vector<s_plan>::iterator it = spfplan.begin(); it != spfplan.end(); it++)
{
mexPrintf("----------------------------------------------------------------------------------------------------\n");
mexPrintf("perfect foresight n°%d\n", i+1);
if (it->exo.length())
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<pair<int, double> >::iterator it1 = it->per_value.begin(); it1 != it->per_value.end(); it1++)
{
mexPrintf(" %3d %10.5f\n",it1->first, it1->second);
}
i++;
}
}
int field_steady_state = mxGetFieldNumber(oo_, "steady_state");
if (field_steady_state < 0)
DYN_MEX_FUNC_ERR_MSG_TXT("steady_state is not a field of oo_");
int field_exo_steady_state = mxGetFieldNumber(oo_, "exo_steady_state");
if (field_exo_steady_state < 0)
DYN_MEX_FUNC_ERR_MSG_TXT("exo_steady_state is not a field of oo_");
if (!steady_state)
{
int field_endo_simul = mxGetFieldNumber(oo_, "endo_simul");
if (field_endo_simul < 0)
DYN_MEX_FUNC_ERR_MSG_TXT("endo_simul is not a field of oo_");
int field_exo_simul = mxGetFieldNumber(oo_, "exo_simul");
if (field_exo_simul < 0)
DYN_MEX_FUNC_ERR_MSG_TXT("exo_simul is not a field of oo_");
if (!count_array_argument)
{
yd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "endo_simul")));
row_y = mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "endo_simul")));
col_y = mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "endo_simul")));
xd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_simul")));
row_x = mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_simul")));
col_x = mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_simul")));
mxArray* endo_sim_arr = mxGetFieldByNumber(oo_, 0, field_endo_simul);
yd = mxGetPr(endo_sim_arr);
row_y = mxGetM(endo_sim_arr);
col_y = mxGetN(endo_sim_arr);
mxArray* exo_sim_arr = mxGetFieldByNumber(oo_, 0, field_exo_simul);
xd = mxGetPr(exo_sim_arr);
row_x = mxGetM(exo_sim_arr);
col_x = mxGetN(exo_sim_arr);
nb_row_xd = row_x;
}
int field = mxGetFieldNumber(M_, "maximum_lag");
if (field >= 0)
y_kmin = int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, field)))));
else
DYN_MEX_FUNC_ERR_MSG_TXT("maximum_lag is not a field of M_");
field = mxGetFieldNumber(M_, "maximum_lead");
if (field >= 0)
y_kmax = int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, field)))));
else
DYN_MEX_FUNC_ERR_MSG_TXT("maximum_lead is not a field of M_");
field = mxGetFieldNumber(M_, "maximum_endo_lag");
if (field >= 0)
y_decal = max(0, y_kmin-int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, field))))));
else
DYN_MEX_FUNC_ERR_MSG_TXT("maximum_endo_lag is not a field of M_");
y_kmin = int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "maximum_lag"))))));
y_kmax = int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "maximum_lead"))))));
y_decal = max(0, y_kmin-int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "maximum_endo_lag")))))));
if (!count_array_argument)
periods = int (floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "periods"))))));
{
int field = mxGetFieldNumber(options_, "periods");
if (field >= 0)
periods = int (floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, field)))));
else
DYN_MEX_FUNC_ERR_MSG_TXT("options_ is not a field of options_");
}
if (!steady_yd )
{
steady_yd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "steady_state")));
steady_row_y = mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "steady_state")));
steady_col_y = mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "steady_state")));;
mxArray* steady_state_arr = mxGetFieldByNumber(oo_, 0, field_steady_state);
steady_yd = mxGetPr(steady_state_arr);
steady_row_y = mxGetM(steady_state_arr);
steady_col_y = mxGetN(steady_state_arr);
}
steady_xd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_steady_state")));
steady_xd = mxGetPr(mxGetFieldByNumber(oo_, 0, field_exo_steady_state));
}
else
{
if (!count_array_argument)
{
yd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "steady_state")));
row_y = mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "steady_state")));
col_y = mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "steady_state")));;
mxArray* steady_state_arr = mxGetFieldByNumber(oo_, 0, field_steady_state);
yd = mxGetPr(steady_state_arr);
row_y = mxGetM(steady_state_arr);
col_y = mxGetN(steady_state_arr);
xd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_steady_state")));
row_x = mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_steady_state")));
col_x = mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_steady_state")));
mxArray* exo_steady_state_arr = mxGetFieldByNumber(oo_, 0, field_exo_steady_state);
xd = mxGetPr(exo_steady_state_arr);
row_x = mxGetM(exo_steady_state_arr);
col_x = mxGetN(exo_steady_state_arr);
nb_row_xd = row_x;
}
}
int verbose= int(*mxGetPr((mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "verbosity")))));
int field = mxGetFieldNumber(options_, "verbosity");
int verbose = 0;
if (field >= 0)
verbose = int(*mxGetPr((mxGetFieldByNumber(options_, 0, field))));
else
DYN_MEX_FUNC_ERR_MSG_TXT("verbosity is not a field of options_");
if (verbose)
print_it = true;
int maxit_ = int (floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "maxit_"))))));
double slowc = double (*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "slowc")))));
double markowitz_c = double (*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "markowitz")))));
int minimal_solving_periods = int (*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "minimal_solving_periods")))));
int stack_solve_algo = int (*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "stack_solve_algo")))));
field = mxGetFieldNumber(options_, "maxit_");
if (field < 0)
DYN_MEX_FUNC_ERR_MSG_TXT("maxit_ is not a field of options_");
int maxit_ = int (floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, field)))));
field = mxGetFieldNumber(options_, "slowc");
if (field < 0)
DYN_MEX_FUNC_ERR_MSG_TXT("slows is not a field of options_");
double slowc = double (*(mxGetPr(mxGetFieldByNumber(options_, 0, field))));
field = mxGetFieldNumber(options_, "markowitz");
if (field < 0)
DYN_MEX_FUNC_ERR_MSG_TXT("markowitz is not a field of options_");
double markowitz_c = double (*(mxGetPr(mxGetFieldByNumber(options_, 0, field))));
field = mxGetFieldNumber(options_, "minimal_solving_periods");
if (field < 0)
DYN_MEX_FUNC_ERR_MSG_TXT("minimal_solving_periods is not a field of options_");
int minimal_solving_periods = int (*(mxGetPr(mxGetFieldByNumber(options_, 0, field))));
field = mxGetFieldNumber(options_, "stack_solve_algo");
if (field < 0)
DYN_MEX_FUNC_ERR_MSG_TXT("stack_solve_algo is not a field of options_");
int stack_solve_algo = int (*(mxGetPr(mxGetFieldByNumber(options_, 0, field))));
int solve_algo;
double solve_tolf;
if (steady_state)
{
solve_algo = int (*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "solve_algo")))));
solve_tolf = *(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "solve_tolf"))));
int field = mxGetFieldNumber(options_, "solve_algo");
if (field >= 0)
solve_algo = int (*(mxGetPr(mxGetFieldByNumber(options_, 0, field))));
else
DYN_MEX_FUNC_ERR_MSG_TXT("solve_algo is not a field of options_");
field = mxGetFieldNumber(options_, "solve_tolf");
if (field >= 0)
solve_tolf = *(mxGetPr(mxGetFieldByNumber(options_, 0, field)));
else
DYN_MEX_FUNC_ERR_MSG_TXT("solve_tolf is not a field of options_");
}
else
{
solve_algo = stack_solve_algo;
mxArray *dynatol = mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "dynatol"));
solve_tolf= *mxGetPr((mxGetFieldByNumber(dynatol, 0, mxGetFieldNumber(dynatol, "f"))));
int field = mxGetFieldNumber(options_, "dynatol");
mxArray *dynatol;
if (field >= 0)
dynatol = mxGetFieldByNumber(options_, 0, field);
else
DYN_MEX_FUNC_ERR_MSG_TXT("dynatol is not a field of options_");
field = mxGetFieldNumber(dynatol, "f");
if (field >= 0)
solve_tolf= *mxGetPr((mxGetFieldByNumber(dynatol, 0, field)));
else
DYN_MEX_FUNC_ERR_MSG_TXT("f is not a field of options_.dynatol");
}
mxArray *mxa = mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "fname"));
int buflen = mxGetM(mxa) * mxGetN(mxa) + 1;
field = mxGetFieldNumber(M_, "fname");
mxArray *mxa;
if (field >= 0)
mxa = mxGetFieldByNumber(M_, 0, field);
else
DYN_MEX_FUNC_ERR_MSG_TXT("fname is not a field of M_");
size_t buflen = mxGetM(mxa) * mxGetN(mxa) + 1;
char *fname;
fname = (char *) mxCalloc(buflen+1, sizeof(char));
int status = mxGetString(mxa, fname, buflen);
size_t status = mxGetString(mxa, fname, int(buflen));
fname[buflen] = ' ';
if (status != 0)
mexWarnMsgTxt("Not enough space. Filename is truncated.");
string file_name = fname;
int size_of_direction = col_y*row_y*sizeof(double);
double *y = (double *) mxMalloc(size_of_direction);
double *ya = (double *) mxMalloc(size_of_direction);
direction = (double *) mxMalloc(size_of_direction);
memset(direction, 0, size_of_direction);
double *x = (double *) mxMalloc(col_x*row_x*sizeof(double));
for (i = 0; i < row_x*col_x; i++)
x[i] = double (xd[i]);
for (i = 0; i < row_y*col_y; i++)
{
y[i] = double (yd[i]);
ya[i] = double (yd[i]);
}
int y_size = row_y;
int nb_row_x = row_x;
clock_t t0 = clock();
Interpreter interprete(params, y, ya, x, steady_yd, steady_xd, direction, y_size, nb_row_x, nb_row_xd, periods, y_kmin, y_kmax, maxit_, solve_tolf, size_of_direction, slowc, y_decal, markowitz_c, file_name, minimal_solving_periods, stack_solve_algo, solve_algo, global_temporary_terms, print, print_error, GlobalTemporaryTerms);
string f(fname);
mxFree(fname);
int nb_blocks = 0;
double *pind;
bool no_error = true;
#ifdef CUDA
try
{
interprete.compute_blocks(f, f, steady_state, evaluate, block, nb_blocks,print_it);
if (stack_solve_algo == 7 && !steady_state)
CUDA_device = GPU_Test_and_Info(&cublas_handle, &cusparse_handle, &descr);
}
catch (GeneralExceptionHandling &feh)
{
DYN_MEX_FUNC_ERR_MSG_TXT(feh.GetErrorMsg().c_str());
}
#else
if (stack_solve_algo == 7 && !steady_state)
DYN_MEX_FUNC_ERR_MSG_TXT("bytecode has not been compiled with CUDA option. Bytecode Can't use options_.stack_solve_algo=7\n");
#endif
size_t size_of_direction = col_y*row_y*sizeof(double);
double *y = (double *) mxMalloc(size_of_direction);
double *ya = (double *) mxMalloc(size_of_direction);
direction = (double *) mxMalloc(size_of_direction);
memset(direction, 0, size_of_direction);
double *x = (double *) mxMalloc(col_x*row_x*sizeof(double));
#ifdef USE_OMP
#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
#endif
for (i = 0; i < row_x*col_x; i++)
{
x[i] = double (xd[i]);
}
#ifdef USE_OMP
#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
#endif
for (i = 0; i < row_y*col_y; i++)
{
y[i] = double (yd[i]);
ya[i] = double (yd[i]);
}
size_t y_size = row_y;
size_t nb_row_x = row_x;
clock_t t0 = clock();
Interpreter interprete(params, y, ya, x, steady_yd, steady_xd, direction, y_size, nb_row_x, nb_row_xd, periods, y_kmin, y_kmax, maxit_, solve_tolf, size_of_direction, slowc, y_decal,
markowitz_c, file_name, minimal_solving_periods, stack_solve_algo, solve_algo, global_temporary_terms, print, print_error, GlobalTemporaryTerms, steady_state,
print_it
#ifdef CUDA
, CUDA_device, cublas_handle, cusparse_handle, descr
#endif
);
string f(fname);
mxFree(fname);
int nb_blocks = 0;
double *pind;
bool no_error = true;
try
{
interprete.compute_blocks(f, f, evaluate, block, nb_blocks);
}
catch (GeneralExceptionHandling &feh)
{
DYN_MEX_FUNC_ERR_MSG_TXT(feh.GetErrorMsg().c_str());
}
#ifdef CUDA
if (stack_solve_algo == 7 && !steady_state)
GPU_close(cublas_handle, cusparse_handle, descr);
#endif
clock_t t1 = clock();
if (!steady_state && !evaluate && no_error && print)
@ -370,14 +859,14 @@ main(int nrhs, const char *prhs[])
if (evaluate)
{
vector<double> residual = interprete.get_residual();
plhs[1] = mxCreateDoubleMatrix(residual.size()/col_y, col_y, mxREAL);
plhs[1] = mxCreateDoubleMatrix(int(residual.size()/double(col_y)), int(col_y), mxREAL);
pind = mxGetPr(plhs[1]);
for (i = 0; i < residual.size(); i++)
pind[i] = residual[i];
}
else
{
plhs[1] = mxCreateDoubleMatrix(row_y, col_y, mxREAL);
plhs[1] = mxCreateDoubleMatrix(int(row_y), int(col_y), mxREAL);
pind = mxGetPr(plhs[1]);
for (i = 0; i < row_y*col_y; i++)
pind[i] = y[i];
@ -385,7 +874,7 @@ main(int nrhs, const char *prhs[])
}
else
{
plhs[1] = mxCreateDoubleMatrix(row_y, col_y, mxREAL);
plhs[1] = mxCreateDoubleMatrix(int(row_y), int(col_y), mxREAL);
pind = mxGetPr(plhs[1]);
if (evaluate)
{
@ -409,7 +898,7 @@ main(int nrhs, const char *prhs[])
jacob_exo_field_number = 1;
jacob_exo_det_field_number = 2;
jacob_other_endo_field_number = 3;
mwSize dims[1] = {nb_blocks };
mwSize dims[1] = {(mwSize)nb_blocks };
plhs[2] = mxCreateStructArray(1, dims, 4, field_names);
}
else if (!mxIsStruct(block_structur))
@ -456,15 +945,15 @@ main(int nrhs, const char *prhs[])
}
if (nlhs > 3)
{
plhs[3] = mxCreateDoubleMatrix(row_y, col_y, mxREAL);
plhs[3] = mxCreateDoubleMatrix(int(row_y), int(col_y), mxREAL);
pind = mxGetPr(plhs[3]);
for (i = 0; i < row_y*col_y; i++)
pind[i] = y[i];
if (nlhs > 4)
{
mxArray *GlobalTemporaryTerms = interprete.get_Temporary_Terms();
unsigned int nb_temp_terms = mxGetM(GlobalTemporaryTerms);
plhs[4] = mxCreateDoubleMatrix(nb_temp_terms, 1, mxREAL);
size_t nb_temp_terms = mxGetM(GlobalTemporaryTerms);
plhs[4] = mxCreateDoubleMatrix(int(nb_temp_terms), 1, mxREAL);
pind = mxGetPr(plhs[4]);
double *tt = mxGetPr(GlobalTemporaryTerms);
for (i = 0; i < nb_temp_terms; i++)
@ -486,4 +975,8 @@ main(int nrhs, const char *prhs[])
mxFree(ya);
if (direction)
mxFree(direction);
#ifdef _MSC_VER_
/*fFreeResult =*/ FreeLibrary(hinstLib);
#endif
return;
}

2
mex/sources/dynblas.h
View File

@ -41,7 +41,7 @@ typedef ptrdiff_t blas_int;
typedef int blas_int;
#endif
#if defined(MATLAB_MEX_FILE) && defined(_WIN32)
#if defined(MATLAB_MEX_FILE) && defined(_WIN32) && !defined(_MSC_VER)
# define FORTRAN_WRAPPER(x) x
#else
# define FORTRAN_WRAPPER(x) x ## _