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

Beautified MEX source code 

git-svn-id: https://www.dynare.org/svn/dynare/trunk@3251 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
sebastien 2009-12-16 17:18:38 +00:00
parent 502e3e1df8
commit f3549b4c64
25 changed files with 2576 additions and 2575 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
mex/sources/build_matlab.m
View File

@ -36,9 +36,9 @@ if strcmpi('GLNX86', computer) || strcmpi('GLNXA64', computer) ...
elseif strcmpi('PCWIN', computer) || strcmpi('PCWIN64', computer) elseif strcmpi('PCWIN', computer) || strcmpi('PCWIN64', computer)
% Windows (x86-32 or x86-64) with Microsoft or gcc compiler % Windows (x86-32 or x86-64) with Microsoft or gcc compiler
if strcmpi('PCWIN', computer) if strcmpi('PCWIN', computer)
LIBRARY_PATH = [MATLAB_PATH '/extern/lib/win32/microsoft/']; LIBRARY_PATH = [MATLAB_PATH '/extern/lib/win32/microsoft/'];
else else
LIBRARY_PATH = [MATLAB_PATH '/extern/lib/win64/microsoft/']; LIBRARY_PATH = [MATLAB_PATH '/extern/lib/win64/microsoft/'];
end end
LAPACK_PATH = ['"' LIBRARY_PATH 'libmwlapack.lib"']; LAPACK_PATH = ['"' LIBRARY_PATH 'libmwlapack.lib"'];
if matlab_ver_less_than('7.5') if matlab_ver_less_than('7.5')
@ -59,7 +59,7 @@ COMPILE_OPTIONS = [ COMPILE_OPTIONS ' -DMATLAB_MEX_FILE -DMATLAB_VERSION=0x' spr
% Large array dims for 64 bits platforms appeared in Matlab 7.3 % Large array dims for 64 bits platforms appeared in Matlab 7.3
if (strcmpi('GLNXA64', computer) || strcmpi('PCWIN64', computer)) ... if (strcmpi('GLNXA64', computer) || strcmpi('PCWIN64', computer)) ...
&& ~matlab_ver_less_than('7.3') && ~matlab_ver_less_than('7.3')
COMPILE_OPTIONS = [ COMPILE_OPTIONS ' -largeArrayDims' ]; COMPILE_OPTIONS = [ COMPILE_OPTIONS ' -largeArrayDims' ];
end end

8
mex/sources/build_matlab_multithread.m
View File

@ -36,9 +36,9 @@ if strcmpi('GLNX86', computer) || strcmpi('GLNXA64', computer) ...
elseif strcmpi('PCWIN', computer) || strcmpi('PCWIN64', computer) elseif strcmpi('PCWIN', computer) || strcmpi('PCWIN64', computer)
% Windows (x86-32 or x86-64) with Microsoft or gcc compiler % Windows (x86-32 or x86-64) with Microsoft or gcc compiler
if strcmpi('PCWIN', computer) if strcmpi('PCWIN', computer)
LIBRARY_PATH = [MATLAB_PATH '/extern/lib/win32/microsoft/']; LIBRARY_PATH = [MATLAB_PATH '/extern/lib/win32/microsoft/'];
else else
LIBRARY_PATH = [MATLAB_PATH '/extern/lib/win64/microsoft/']; LIBRARY_PATH = [MATLAB_PATH '/extern/lib/win64/microsoft/'];
end end
LAPACK_PATH = ['"' LIBRARY_PATH 'libmwlapack.lib"']; LAPACK_PATH = ['"' LIBRARY_PATH 'libmwlapack.lib"'];
if matlab_ver_less_than('7.5') if matlab_ver_less_than('7.5')
@ -58,7 +58,7 @@ COMPILE_OPTIONS = [ COMPILE_OPTIONS ' -DMATLAB_MEX_FILE -DMATLAB_VERSION=0x' spr
% Large array dims for 64 bits platforms appeared in Matlab 7.3 % Large array dims for 64 bits platforms appeared in Matlab 7.3
if (strcmpi('GLNXA64', computer) || strcmpi('PCWIN64', computer)) ... if (strcmpi('GLNXA64', computer) || strcmpi('PCWIN64', computer)) ...
&& ~matlab_ver_less_than('7.3') && ~matlab_ver_less_than('7.3')
COMPILE_OPTIONS = [ COMPILE_OPTIONS ' -largeArrayDims' ]; COMPILE_OPTIONS = [ COMPILE_OPTIONS ' -largeArrayDims' ];
end end
@ -122,7 +122,7 @@ eval([ COMPILE_COMMAND ' -I. mjdgges/mjdgges.c ' LAPACK_PATH ]);
disp('Compiling sparse_hessian_times_B_kronecker_C...') disp('Compiling sparse_hessian_times_B_kronecker_C...')
eval([ COMPILE_COMMAND_OMP ' -I. kronecker/sparse_hessian_times_B_kronecker_C.cc' ]); eval([ COMPILE_COMMAND_OMP ' -I. kronecker/sparse_hessian_times_B_kronecker_C.cc' ]);
disp('Compiling A_times_B_kronecker_C...') disp('Compiling A_times_B_kronecker_C...')
eval([ COMPILE_COMMAND_OMP ' -I. kronecker/A_times_B_kronecker_C.cc ']); eval([ COMPILE_COMMAND_OMP ' -I. kronecker/A_times_B_kronecker_C.cc ']);

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

File diff suppressed because it is too large Load Diff

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

@ -31,54 +31,52 @@
# include "linbcg.hh" # include "linbcg.hh"
#endif #endif
#ifndef DEBUG_EX #ifndef DEBUG_EX
#include "mex.h" # include "mex.h"
#else #else
#include "mex_interface.hh" # include "mex_interface.hh"
#endif #endif
//#define DEBUGC //#define DEBUGC
using namespace std; using namespace std;
#define pow_ pow #define pow_ pow
typedef vector<pair<Tags, void* > >::const_iterator it_code_type; typedef vector<pair<Tags, void * > >::const_iterator it_code_type;
class Interpreter : SparseMatrix class Interpreter : SparseMatrix
{ {
protected : protected:
double pow1(double a, double b); double pow1(double a, double b);
double log1(double a); double log1(double a);
void compute_block_time(int Per_u_, bool evaluate, int block_num); void compute_block_time(int Per_u_, bool evaluate, int block_num);
void evaluate_a_block(const int size, const int type, string bin_basename, bool steady_state, int block_num, 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); 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);
bool simulate_a_block(const int size, const int type, string file_name, string bin_basename, bool Gaussian_Elimination, bool steady_state, int block_num, bool simulate_a_block(const int size, const int type, string file_name, string bin_basename, bool Gaussian_Elimination, 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); 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);
double *T; double *T;
vector<Block_contain_type> Block_Contain; vector<Block_contain_type> Block_Contain;
vector<pair<Tags, void* > > code_liste; vector<pair<Tags, void * > > code_liste;
it_code_type it_code; it_code_type it_code;
stack<double> Stack; stack<double> Stack;
int Block_Count, Per_u_, Per_y_; int Block_Count, Per_u_, Per_y_;
int it_, nb_row_x, nb_row_xd, maxit_, size_of_direction; int it_, nb_row_x, nb_row_xd, maxit_, size_of_direction;
double *g1, *r; double *g1, *r;
double solve_tolf; double solve_tolf;
bool GaussSeidel; bool GaussSeidel;
double *x, *params; double *x, *params;
double *steady_y, *steady_x; double *steady_y, *steady_x;
map<pair<pair<int, int> ,int>, int> IM_i; map<pair<pair<int, int>, int>, int> IM_i;
int equation, derivative_equation, derivative_variable; int equation, derivative_equation, derivative_variable;
string filename; string filename;
int minimal_solving_periods; int minimal_solving_periods;
public : public:
Interpreter(double *params_arg, double *y_arg, double *ya_arg, double *x_arg, double *steady_y_arg, double *steady_x_arg, 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, 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, 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); double slowc_arg, int y_decal_arg, double markowitz_c_arg, string &filename_arg, int minimal_solving_periods_arg);
bool compute_blocks(string file_name, string bin_basename, bool steady_state, bool evaluate); bool compute_blocks(string file_name, string bin_basename, bool steady_state, bool evaluate);
}; };
#endif #endif

79
mex/sources/bytecode/Mem_Mngr.cc
View File

@ -21,16 +21,16 @@
Mem_Mngr::Mem_Mngr() Mem_Mngr::Mem_Mngr()
{ {
swp_f=false; swp_f = false;
swp_f_b=0; swp_f_b = 0;
} }
void void
Mem_Mngr::Print_heap() Mem_Mngr::Print_heap()
{ {
int i; int i;
mexPrintf("i :"); mexPrintf("i :");
for (i=0;i<CHUNK_SIZE;i++) for (i = 0; i < CHUNK_SIZE; i++)
mexPrintf("%3d ",i); mexPrintf("%3d ", i);
mexPrintf("\n"); mexPrintf("\n");
} }
@ -38,90 +38,89 @@ void
Mem_Mngr::init_Mem() Mem_Mngr::init_Mem()
{ {
Chunk_Stack.clear(); Chunk_Stack.clear();
CHUNK_SIZE=0; CHUNK_SIZE = 0;
Nb_CHUNK=0; Nb_CHUNK = 0;
NZE_Mem=NULL; NZE_Mem = NULL;
NZE_Mem_add=NULL; NZE_Mem_add = NULL;
CHUNK_heap_pos=0; CHUNK_heap_pos = 0;
NZE_Mem_Allocated.clear(); NZE_Mem_Allocated.clear();
} }
void Mem_Mngr::fixe_file_name(string filename_arg) void
Mem_Mngr::fixe_file_name(string filename_arg)
{ {
filename=filename_arg; filename = filename_arg;
} }
void void
Mem_Mngr::init_CHUNK_BLCK_SIZE(int u_count) Mem_Mngr::init_CHUNK_BLCK_SIZE(int u_count)
{ {
CHUNK_BLCK_SIZE=u_count; CHUNK_BLCK_SIZE = u_count;
} }
NonZeroElem* NonZeroElem *
Mem_Mngr::mxMalloc_NZE() Mem_Mngr::mxMalloc_NZE()
{ {
long int i; long 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(); NonZeroElem *p1 = Chunk_Stack.back();
Chunk_Stack.pop_back(); Chunk_Stack.pop_back();
return(p1); return (p1);
} }
else if (CHUNK_heap_pos<CHUNK_SIZE) /*there is enough allocated memory space available we keep it at the top of the heap*/ else if (CHUNK_heap_pos < CHUNK_SIZE) /*there is enough allocated memory space available we keep it at the top of the heap*/
{ {
i=CHUNK_heap_pos++; i = CHUNK_heap_pos++;
return(NZE_Mem_add[i]); 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; CHUNK_SIZE += CHUNK_BLCK_SIZE;
Nb_CHUNK++; Nb_CHUNK++;
NZE_Mem=(NonZeroElem*)mxMalloc(CHUNK_BLCK_SIZE*sizeof(NonZeroElem)); /*The block of memory allocated*/ NZE_Mem = (NonZeroElem *) mxMalloc(CHUNK_BLCK_SIZE*sizeof(NonZeroElem)); /*The block of memory allocated*/
NZE_Mem_Allocated.push_back(NZE_Mem); NZE_Mem_Allocated.push_back(NZE_Mem);
if(!NZE_Mem) if (!NZE_Mem)
{ {
mexPrintf("Not enough memory available\n"); mexPrintf("Not enough memory available\n");
mexEvalString("drawnow;"); mexEvalString("drawnow;");
} }
NZE_Mem_add=(NonZeroElem**)mxRealloc(NZE_Mem_add, CHUNK_SIZE*sizeof(NonZeroElem*)); /*We have to redefine the size of pointer on the memory*/ NZE_Mem_add = (NonZeroElem **) mxRealloc(NZE_Mem_add, CHUNK_SIZE*sizeof(NonZeroElem *)); /*We have to redefine the size of pointer on the memory*/
if(!NZE_Mem_add) if (!NZE_Mem_add)
{ {
mexPrintf("Not enough memory available\n"); mexPrintf("Not enough memory available\n");
mexEvalString("drawnow;"); mexEvalString("drawnow;");
} }
for (i=CHUNK_heap_pos;i<CHUNK_SIZE;i++) for (i = CHUNK_heap_pos; i < CHUNK_SIZE; i++)
{ {
NZE_Mem_add[i]=(NonZeroElem*)(NZE_Mem+(i-CHUNK_heap_pos)); NZE_Mem_add[i] = (NonZeroElem *)(NZE_Mem+(i-CHUNK_heap_pos));
} }
i=CHUNK_heap_pos++; i = CHUNK_heap_pos++;
return(NZE_Mem_add[i]); return (NZE_Mem_add[i]);
} }
} }
void void
Mem_Mngr::mxFree_NZE(void* pos) Mem_Mngr::mxFree_NZE(void *pos)
{ {
int i; int i;
size_t gap; size_t gap;
for (i=0;i<Nb_CHUNK;i++) for (i = 0; i < Nb_CHUNK; i++)
{ {
gap=((size_t)(pos)-(size_t)(NZE_Mem_add[i*CHUNK_BLCK_SIZE]))/sizeof(NonZeroElem); gap = ((size_t)(pos)-(size_t)(NZE_Mem_add[i*CHUNK_BLCK_SIZE]))/sizeof(NonZeroElem);
if ((gap<CHUNK_BLCK_SIZE) && (gap>=0)) if ((gap < CHUNK_BLCK_SIZE) && (gap >= 0))
break; break;
} }
Chunk_Stack.push_back((NonZeroElem*)pos); Chunk_Stack.push_back((NonZeroElem *) pos);
} }
void void
Mem_Mngr::Free_All() Mem_Mngr::Free_All()
{ {
while(NZE_Mem_Allocated.size()) while (NZE_Mem_Allocated.size())
{ {
mxFree(NZE_Mem_Allocated.back()); mxFree(NZE_Mem_Allocated.back());
NZE_Mem_Allocated.pop_back(); NZE_Mem_Allocated.pop_back();
} }
mxFree(NZE_Mem_add); mxFree(NZE_Mem_add);
init_Mem(); init_Mem();
} }

53
mex/sources/bytecode/Mem_Mngr.hh
View File

@ -20,47 +20,46 @@
#ifndef MEM_MNGR_HH_INCLUDED #ifndef MEM_MNGR_HH_INCLUDED
#define MEM_MNGR_HH_INCLUDED #define MEM_MNGR_HH_INCLUDED
#include <vector> #include <vector>
#include <fstream> #include <fstream>
#ifndef DEBUG_EX #ifndef DEBUG_EX
#include "mex.h" # include "mex.h"
#else #else
#include "mex_interface.hh" # include "mex_interface.hh"
#endif #endif
using namespace std; using namespace std;
struct NonZeroElem struct NonZeroElem
{ {
int u_index; int u_index;
int r_index, c_index, lag_index; int r_index, c_index, lag_index;
NonZeroElem *NZE_R_N, *NZE_C_N; NonZeroElem *NZE_R_N, *NZE_C_N;
}; };
typedef vector<NonZeroElem*> v_NonZeroElem; typedef vector<NonZeroElem *> v_NonZeroElem;
class Mem_Mngr class Mem_Mngr
{ {
public: public:
void Print_heap(); void Print_heap();
void init_Mem(); void init_Mem();
void mxFree_NZE(void* pos); void mxFree_NZE(void *pos);
NonZeroElem* mxMalloc_NZE(); NonZeroElem *mxMalloc_NZE();
void init_CHUNK_BLCK_SIZE(int u_count); void init_CHUNK_BLCK_SIZE(int u_count);
void Free_All(); void Free_All();
Mem_Mngr(); Mem_Mngr();
void fixe_file_name(string filename_arg); void fixe_file_name(string filename_arg);
bool swp_f; bool swp_f;
private: private:
v_NonZeroElem Chunk_Stack; v_NonZeroElem Chunk_Stack;
int CHUNK_SIZE, CHUNK_BLCK_SIZE, Nb_CHUNK; int CHUNK_SIZE, CHUNK_BLCK_SIZE, Nb_CHUNK;
int CHUNK_heap_pos; int CHUNK_heap_pos;
NonZeroElem** NZE_Mem_add; NonZeroElem **NZE_Mem_add;
NonZeroElem* NZE_Mem; NonZeroElem *NZE_Mem;
vector<NonZeroElem*> NZE_Mem_Allocated; vector<NonZeroElem *> NZE_Mem_Allocated;
int swp_f_b; int swp_f_b;
fstream SaveCode_swp; fstream SaveCode_swp;
string filename; string filename;
}; };
#endif #endif

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

@ -141,7 +141,7 @@ SparseMatrix::At_Col(int c, int lag, NonZeroElem **first)
void void
SparseMatrix::Delete(const int r, const int c) SparseMatrix::Delete(const int r, const int c)
{ {
NonZeroElem *first = FNZE_R[r], *firsta = NULL; NonZeroElem *first = FNZE_R[r], *firsta = NULL;
while (first->c_index != c) while (first->c_index != c)
{ {
@ -165,7 +165,7 @@ SparseMatrix::Delete(const int r, const int c)
if (firsta != NULL) if (firsta != NULL)
firsta->NZE_C_N = first->NZE_C_N; firsta->NZE_C_N = first->NZE_C_N;
if (first == FNZE_C[c]) if (first == FNZE_C[c])
FNZE_C[c] = first->NZE_C_N; FNZE_C[c] = first->NZE_C_N;
u_liste.push_back(first->u_index); u_liste.push_back(first->u_index);
mem_mngr.mxFree_NZE(first); mem_mngr.mxFree_NZE(first);
@ -283,24 +283,24 @@ SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, i
mem_mngr.fixe_file_name(file_name); mem_mngr.fixe_file_name(file_name);
if (!SaveCode.is_open()) if (!SaveCode.is_open())
{ {
if(steady_state) if (steady_state)
SaveCode.open((file_name + "_static.bin").c_str(), ios::in | ios::binary); SaveCode.open((file_name + "_static.bin").c_str(), ios::in | ios::binary);
else else
SaveCode.open((file_name + "_dynamic.bin").c_str(), ios::in | ios::binary); SaveCode.open((file_name + "_dynamic.bin").c_str(), ios::in | ios::binary);
if (!SaveCode.is_open()) if (!SaveCode.is_open())
{ {
if(steady_state) if (steady_state)
mexPrintf("Error : Can't open file \"%s\" for reading\n", (file_name + "_static.bin").c_str()); mexPrintf("Error : Can't open file \"%s\" for reading\n", (file_name + "_static.bin").c_str());
else else
mexPrintf("Error : Can't open file \"%s\" for reading\n", (file_name + "_dynamic.bin").c_str()); mexPrintf("Error : Can't open file \"%s\" for reading\n", (file_name + "_dynamic.bin").c_str());
mexEvalString("st=fclose('all');clear all;"); mexEvalString("st=fclose('all');clear all;");
mexErrMsgTxt("Exit from Dynare"); mexErrMsgTxt("Exit from Dynare");
} }
} }
IM_i.clear(); IM_i.clear();
if(two_boundaries) if (two_boundaries)
{ {
for (i = 0; i < u_count_init-Size; i++) for (i = 0; i < u_count_init-Size; i++)
{ {
SaveCode.read(reinterpret_cast<char *>(&eq), sizeof(eq)); SaveCode.read(reinterpret_cast<char *>(&eq), sizeof(eq));
SaveCode.read(reinterpret_cast<char *>(&var), sizeof(var)); SaveCode.read(reinterpret_cast<char *>(&var), sizeof(var));
@ -308,12 +308,12 @@ SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, i
SaveCode.read(reinterpret_cast<char *>(&j), sizeof(j)); SaveCode.read(reinterpret_cast<char *>(&j), sizeof(j));
IM_i[make_pair(make_pair(eq, var), lag)] = j; IM_i[make_pair(make_pair(eq, var), lag)] = j;
} }
for (j=0;j<Size;j++) for (j = 0; j < Size; j++)
IM_i[make_pair(make_pair(j, Size*(periods+y_kmax)), 0)] = j; IM_i[make_pair(make_pair(j, Size*(periods+y_kmax)), 0)] = j;
} }
else else
{ {
for (i = 0; i < u_count_init; i++) for (i = 0; i < u_count_init; i++)
{ {
SaveCode.read(reinterpret_cast<char *>(&eq), sizeof(eq)); SaveCode.read(reinterpret_cast<char *>(&eq), sizeof(eq));
SaveCode.read(reinterpret_cast<char *>(&var), sizeof(var)); SaveCode.read(reinterpret_cast<char *>(&var), sizeof(var));
@ -321,7 +321,7 @@ SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, i
SaveCode.read(reinterpret_cast<char *>(&j), sizeof(j)); SaveCode.read(reinterpret_cast<char *>(&j), sizeof(j));
IM_i[make_pair(make_pair(eq, var), lag)] = j; IM_i[make_pair(make_pair(eq, var), lag)] = j;
} }
} }
index_vara = (int *) mxMalloc(Size*(periods+y_kmin+y_kmax)*sizeof(int)); index_vara = (int *) mxMalloc(Size*(periods+y_kmin+y_kmax)*sizeof(int));
for (j = 0; j < Size; j++) for (j = 0; j < Size; j++)
SaveCode.read(reinterpret_cast<char *>(&index_vara[j]), sizeof(*index_vara)); SaveCode.read(reinterpret_cast<char *>(&index_vara[j]), sizeof(*index_vara));
@ -415,7 +415,7 @@ SparseMatrix::Simple_Init(int it_, int y_kmin, int y_kmax, int Size, map<pair<pa
} }
//#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS"))) //#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
for (i = 0; i < Size; i++) for (i = 0; i < Size; i++)
b[i] = i; b[i] = i;
mxFree(temp_NZE_R); mxFree(temp_NZE_R);
mxFree(temp_NZE_C); mxFree(temp_NZE_C);
u_count = u_count1; u_count = u_count1;
@ -438,7 +438,7 @@ SparseMatrix::Init(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<
mem_mngr.init_CHUNK_BLCK_SIZE(u_count); mem_mngr.init_CHUNK_BLCK_SIZE(u_count);
g_save_op = NULL; g_save_op = NULL;
g_nop_all = 0; g_nop_all = 0;
i = (periods+y_kmax+1)*Size*sizeof(NonZeroElem*); i = (periods+y_kmax+1)*Size*sizeof(NonZeroElem *);
FNZE_R = (NonZeroElem **) mxMalloc(i); FNZE_R = (NonZeroElem **) mxMalloc(i);
FNZE_C = (NonZeroElem **) mxMalloc(i); FNZE_C = (NonZeroElem **) mxMalloc(i);
NonZeroElem **temp_NZE_R = (NonZeroElem **) mxMalloc(i); NonZeroElem **temp_NZE_R = (NonZeroElem **) mxMalloc(i);
@ -720,7 +720,7 @@ SparseMatrix::compare(int *save_op, int *save_opa, int *save_opaa, int beg_t, in
i = j = 0; i = j = 0;
while (i < nop4) while (i < nop4)
{ {
save_op_s = (t_save_op_s *) (&(save_op[i])); save_op_s = (t_save_op_s *)(&(save_op[i]));
up = &u[save_op_s->first+t*diff1[j]]; up = &u[save_op_s->first+t*diff1[j]];
switch (save_op_s->operat) switch (save_op_s->operat)
{ {
@ -751,7 +751,7 @@ SparseMatrix::compare(int *save_op, int *save_opa, int *save_opaa, int beg_t, in
i = j = 0; i = j = 0;
while (i < nop4) while (i < nop4)
{ {
save_op_s = (t_save_op_s *) (&(save_op[i])); save_op_s = (t_save_op_s *)(&(save_op[i]));
if (save_op_s->lag < (periods_beg_t-t)) if (save_op_s->lag < (periods_beg_t-t))
{ {
up = &u[save_op_s->first+t*diff1[j]]; up = &u[save_op_s->first+t*diff1[j]];
@ -1018,14 +1018,14 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
double piv_abs; double piv_abs;
NonZeroElem *first, *firsta, *first_suba; NonZeroElem *first, *firsta, *first_suba;
double *piv_v; double *piv_v;
int *pivj_v, *pivk_v, *NR; int *pivj_v, *pivk_v, *NR;
int l, N_max; int l, N_max;
bool one; bool one;
Clear_u(); Clear_u();
piv_v = (double*)mxMalloc(Size*sizeof(double)); piv_v = (double *) mxMalloc(Size*sizeof(double));
pivj_v = (int*)mxMalloc(Size*sizeof(int)); pivj_v = (int *) mxMalloc(Size*sizeof(int));
pivk_v = (int*)mxMalloc(Size*sizeof(int)); pivk_v = (int *) mxMalloc(Size*sizeof(int));
NR = (int*)mxMalloc(Size*sizeof(int)); NR = (int *) mxMalloc(Size*sizeof(int));
error_not_printed = true; error_not_printed = true;
u_count_alloc_save = u_count_alloc; u_count_alloc_save = u_count_alloc;
if (isnan(res1) || isinf(res1)) if (isnan(res1) || isinf(res1))
@ -1034,63 +1034,63 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
{ {
for (j = 0; j < y_size; j++) for (j = 0; j < y_size; j++)
{ {
bool select=false; bool select = false;
for(int i = 0; i<Size; i++) for (int i = 0; i < Size; i++)
if(j == index_vara[i]) if (j == index_vara[i])
{ {
select=true; select = true;
break; break;
} }
if(select) if (select)
mexPrintf("-> variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]); mexPrintf("-> variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]);
else else
mexPrintf(" variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]); mexPrintf(" variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]);
} }
mexPrintf("res1=%5.25\n",res1); mexPrintf("res1=%5.25\n", res1);
mexPrintf("The initial values of endogenous variables are too far from the solution.\n"); mexPrintf("The initial values of endogenous variables are too far from the solution.\n");
mexPrintf("Change them!\n"); mexPrintf("Change them!\n");
mexEvalString("drawnow;"); mexEvalString("drawnow;");
mxFree(piv_v); mxFree(piv_v);
mxFree(pivj_v); mxFree(pivj_v);
mxFree(pivk_v); mxFree(pivk_v);
mxFree(NR); mxFree(NR);
if(steady_state) if (steady_state)
return false; return false;
else else
{ {
mexEvalString("st=fclose('all');clear all;"); mexEvalString("st=fclose('all');clear all;");
filename += " stopped"; filename += " stopped";
mexErrMsgTxt(filename.c_str()); mexErrMsgTxt(filename.c_str());
} }
} }
if (slowc_save < 1e-8) if (slowc_save < 1e-8)
{ {
for (j = 0; j < y_size; j++) for (j = 0; j < y_size; j++)
{ {
bool select=false; bool select = false;
for(int i = 0; i<Size; i++) for (int i = 0; i < Size; i++)
if(j == index_vara[i]) if (j == index_vara[i])
{ {
select=true; select = true;
break; break;
} }
if(select) if (select)
mexPrintf("-> variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]); mexPrintf("-> variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]);
else else
mexPrintf(" variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]); mexPrintf(" variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]);
} }
mexPrintf("Dynare cannot improve the simulation in block %d at time %d (variable %d)\n", blck+1, it_+1, max_res_idx); mexPrintf("Dynare cannot improve the simulation in block %d at time %d (variable %d)\n", blck+1, it_+1, max_res_idx);
mexEvalString("drawnow;"); mexEvalString("drawnow;");
mxFree(piv_v); mxFree(piv_v);
mxFree(pivj_v); mxFree(pivj_v);
mxFree(pivk_v); mxFree(pivk_v);
mxFree(NR); mxFree(NR);
if(steady_state) if (steady_state)
return false; return false;
else else
{ {
mexEvalString("st=fclose('all');clear all;"); mexEvalString("st=fclose('all');clear all;");
filename += " stopped"; filename += " stopped";
mexErrMsgTxt(filename.c_str()); mexErrMsgTxt(filename.c_str());
} }
} }
@ -1098,24 +1098,24 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
mexPrintf("Error: Simulation diverging, trying to correct it using slowc=%f\n", slowc_save); mexPrintf("Error: Simulation diverging, trying to correct it using slowc=%f\n", slowc_save);
for (i = 0; i < y_size; i++) for (i = 0; i < y_size; i++)
y[i+it_*y_size] = ya[i+it_*y_size] + slowc_save*direction[i+it_*y_size]; y[i+it_*y_size] = ya[i+it_*y_size] + slowc_save*direction[i+it_*y_size];
mxFree(piv_v); mxFree(piv_v);
mxFree(pivj_v); mxFree(pivj_v);
mxFree(pivk_v); mxFree(pivk_v);
mxFree(NR); mxFree(NR);
iter--; iter--;
return true; return true;
} }
if (cvg) if (cvg)
{ {
mxFree(piv_v); mxFree(piv_v);
mxFree(pivj_v); mxFree(pivj_v);
mxFree(pivk_v); mxFree(pivk_v);
mxFree(NR); mxFree(NR);
return (true); return (true);
} }
Simple_Init(it_, y_kmin, y_kmax, Size, IM_i); Simple_Init(it_, y_kmin, y_kmax, Size, IM_i);
NonZeroElem **bc; NonZeroElem **bc;
bc = (NonZeroElem**)mxMalloc(Size*sizeof(*bc)); bc = (NonZeroElem **) mxMalloc(Size*sizeof(*bc));
for (i = 0; i < Size; i++) for (i = 0; i < Size; i++)
{ {
/*finding the max-pivot*/ /*finding the max-pivot*/
@ -1183,7 +1183,7 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
{ {
for (j = 0; j < l; j++) for (j = 0; j < l; j++)
{ {
markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double (NR[j])/double(N_max))); markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double (NR[j])/double (N_max)));
if (markovitz > markovitz_max && NR[j] == 1) if (markovitz > markovitz_max && NR[j] == 1)
{ {
piv = piv_v[j]; piv = piv_v[j];
@ -1199,14 +1199,14 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
line_done[pivj] = true; line_done[pivj] = true;
if (piv_abs < eps) if (piv_abs < eps)
{ {
mexPrintf("Error: singular system in Simulate_NG in block %d\n",blck+1); mexPrintf("Error: singular system in Simulate_NG in block %d\n", blck+1);
mexEvalString("drawnow;"); mexEvalString("drawnow;");
mxFree(piv_v); mxFree(piv_v);
mxFree(pivj_v); mxFree(pivj_v);
mxFree(pivk_v); mxFree(pivk_v);
mxFree(NR); mxFree(NR);
mxFree(bc); mxFree(bc);
if(steady_state) if (steady_state)
return false; return false;
else else
{ {
@ -1236,7 +1236,7 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
} }
//#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS"))) //#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
for (j = 0; j < nb_eq_todo; j++) for (j = 0; j < nb_eq_todo; j++)
{ {
first = bc[j]; first = bc[j];
int row = first->r_index; int row = first->r_index;
double first_elem = u[first->u_index]; double first_elem = u[first->u_index];
@ -1310,18 +1310,18 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
} }
} }
u[b[row]] -= u[b[pivj]]*first_elem; u[b[row]] -= u[b[pivj]]*first_elem;
} }
} }
double slowc_lbx = slowc, res1bx; double slowc_lbx = slowc, res1bx;
for (i = 0; i < y_size; i++) for (i = 0; i < y_size; i++)
ya[i+it_*y_size] = y[i+it_*y_size]; ya[i+it_*y_size] = y[i+it_*y_size];
slowc_save = slowc; slowc_save = slowc;
res1bx = simple_bksub(it_, Size, slowc_lbx); res1bx = simple_bksub(it_, Size, slowc_lbx);
End(Size); End(Size);
mxFree(piv_v); mxFree(piv_v);
mxFree(pivj_v); mxFree(pivj_v);
mxFree(pivk_v); mxFree(pivk_v);
mxFree(NR); mxFree(NR);
mxFree(bc); mxFree(bc);
return true; return true;
} }
@ -1382,8 +1382,8 @@ SparseMatrix::CheckIt(int y_size, int y_kmin, int y_kmax, int Size, int periods,
mexPrintf("G1a red done\n"); mexPrintf("G1a red done\n");
SaveResult >> row; SaveResult >> row;
mexPrintf("row(2)=%d\n", row); mexPrintf("row(2)=%d\n", row);
double *B; double *B;
B = (double*)mxMalloc(row*sizeof(double)); B = (double *) mxMalloc(row*sizeof(double));
for (int i = 0; i < row; i++) for (int i = 0; i < row; i++)
SaveResult >> B[i]; SaveResult >> B[i];
SaveResult.close(); SaveResult.close();
@ -1393,7 +1393,7 @@ SparseMatrix::CheckIt(int y_size, int y_kmin, int y_kmax, int Size, int periods,
{ {
if (abs(u[b[i]]+B[i]) > epsilon) if (abs(u[b[i]]+B[i]) > epsilon)
mexPrintf("Problem at i=%d u[b[i]]=%f B[i]=%f\n", i, u[b[i]], B[i]); mexPrintf("Problem at i=%d u[b[i]]=%f B[i]=%f\n", i, u[b[i]], B[i]);
} }
mxFree(B); mxFree(B);
} }
@ -1470,9 +1470,9 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
{ {
for (j = 0; j < y_size; j++) for (j = 0; j < y_size; j++)
mexPrintf("variable %d at time %d = %f\n", j+1, it_, y[j+it_*y_size]); mexPrintf("variable %d at time %d = %f\n", j+1, it_, y[j+it_*y_size]);
for(j = 0; j < Size; j++) for (j = 0; j < Size; j++)
mexPrintf("residual(%d)=%5.25f\n",j, u[j]); mexPrintf("residual(%d)=%5.25f\n", j, u[j]);
mexPrintf("res1=%5.25f\n",res1); mexPrintf("res1=%5.25f\n", res1);
mexPrintf("The initial values of endogenous variables are too far from the solution.\n"); mexPrintf("The initial values of endogenous variables are too far from the solution.\n");
mexPrintf("Change them!\n"); mexPrintf("Change them!\n");
mexEvalString("drawnow;"); mexEvalString("drawnow;");
@ -1526,7 +1526,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
} }
else else
{ {
start_compare = max( y_kmin, minimal_solving_periods); start_compare = max(y_kmin, minimal_solving_periods);
restart = 0; restart = 0;
} }
res1a = res1; res1a = res1;
@ -1549,10 +1549,10 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
Init(periods, y_kmin, y_kmax, Size, IM_i); Init(periods, y_kmin, y_kmax, Size, IM_i);
double *piv_v; double *piv_v;
int *pivj_v, *pivk_v, *NR; int *pivj_v, *pivk_v, *NR;
piv_v = (double*)mxMalloc(Size*sizeof(double)); piv_v = (double *) mxMalloc(Size*sizeof(double));
pivj_v = (int*)mxMalloc(Size*sizeof(int)); pivj_v = (int *) mxMalloc(Size*sizeof(int));
pivk_v = (int*)mxMalloc(Size*sizeof(int)); pivk_v = (int *) mxMalloc(Size*sizeof(int));
NR = (int*)mxMalloc(Size*sizeof(int)); NR = (int *) mxMalloc(Size*sizeof(int));
for (int t = 0; t < periods; t++) for (int t = 0; t < periods; t++)
{ {
if (record && symbolic) if (record && symbolic)
@ -1670,7 +1670,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
nopa = int (1.5*nopa); nopa = int (1.5*nopa);
save_op = (int *) mxRealloc(save_op, nopa*sizeof(int)); save_op = (int *) mxRealloc(save_op, nopa*sizeof(int));
} }
save_op_s = (t_save_op_s *) (&(save_op[nop])); save_op_s = (t_save_op_s *)(&(save_op[nop]));
save_op_s->operat = IFLD; save_op_s->operat = IFLD;
save_op_s->first = pivk; save_op_s->first = pivk;
save_op_s->lag = 0; save_op_s->lag = 0;
@ -1687,7 +1687,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
/*divide all the non zeros elements of the line pivj by the max_pivot*/ /*divide all the non zeros elements of the line pivj by the max_pivot*/
int nb_var = At_Row(pivj, &first); int nb_var = At_Row(pivj, &first);
NonZeroElem **bb; NonZeroElem **bb;
bb = (NonZeroElem**)mxMalloc(nb_var*sizeof(first)); bb = (NonZeroElem **) mxMalloc(nb_var*sizeof(first));
for (j = 0; j < nb_var; j++) for (j = 0; j < nb_var; j++)
{ {
bb[j] = first; bb[j] = first;
@ -1707,13 +1707,13 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
nopa = int (1.5*nopa); nopa = int (1.5*nopa);
save_op = (int *) mxRealloc(save_op, nopa*sizeof(int)); save_op = (int *) mxRealloc(save_op, nopa*sizeof(int));
} }
save_op_s = (t_save_op_s *) (&(save_op[nop+j*2])); save_op_s = (t_save_op_s *)(&(save_op[nop+j*2]));
save_op_s->operat = IFDIV; save_op_s->operat = IFDIV;
save_op_s->first = first->u_index; save_op_s->first = first->u_index;
save_op_s->lag = first->lag_index; save_op_s->lag = first->lag_index;
} }
} }
} }
mxFree(bb); mxFree(bb);
nop += nb_var*2; nop += nb_var*2;
u[b[pivj]] /= piv; u[b[pivj]] /= piv;
@ -1726,7 +1726,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
nopa = int (1.5*nopa); nopa = int (1.5*nopa);
save_op = (int *) mxRealloc(save_op, nopa*sizeof(int)); save_op = (int *) mxRealloc(save_op, nopa*sizeof(int));
} }
save_op_s = (t_save_op_s *) (&(save_op[nop])); save_op_s = (t_save_op_s *)(&(save_op[nop]));
save_op_s->operat = IFDIV; save_op_s->operat = IFDIV;
save_op_s->first = b[pivj]; save_op_s->first = b[pivj];
save_op_s->lag = 0; save_op_s->lag = 0;
@ -1739,7 +1739,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
int nb_var_piva = At_Row(pivj, &first_piva); int nb_var_piva = At_Row(pivj, &first_piva);
NonZeroElem **bc; NonZeroElem **bc;
bc = (NonZeroElem**)mxMalloc(nb_eq*sizeof(first)); bc = (NonZeroElem **) mxMalloc(nb_eq*sizeof(first));
int nb_eq_todo = 0; int nb_eq_todo = 0;
for (j = 0; j < nb_eq && first; j++) for (j = 0; j < nb_eq && first; j++)
{ {
@ -1766,7 +1766,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
save_op = (int *) mxRealloc(save_op, nopa*sizeof(int)); save_op = (int *) mxRealloc(save_op, nopa*sizeof(int));
} }
} }
save_op_s_l = (t_save_op_s *) (&(save_op[nop])); save_op_s_l = (t_save_op_s *)(&(save_op[nop]));
save_op_s_l->operat = IFLD; save_op_s_l->operat = IFLD;
save_op_s_l->first = first->u_index; save_op_s_l->first = first->u_index;
save_op_s_l->lag = abs(first->lag_index); save_op_s_l->lag = abs(first->lag_index);
@ -1817,7 +1817,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
save_op = (int *) mxRealloc(save_op, nopa*sizeof(int)); save_op = (int *) mxRealloc(save_op, nopa*sizeof(int));
} }
} }
save_op_s_l = (t_save_op_s *) (&(save_op[nop])); save_op_s_l = (t_save_op_s *)(&(save_op[nop]));
save_op_s_l->operat = IFLESS; save_op_s_l->operat = IFLESS;
save_op_s_l->first = tmp_u_count; save_op_s_l->first = tmp_u_count;
save_op_s_l->second = first_piv->u_index; save_op_s_l->second = first_piv->u_index;
@ -1875,7 +1875,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
save_op = (int *) mxRealloc(save_op, nopa*sizeof(int)); save_op = (int *) mxRealloc(save_op, nopa*sizeof(int));
} }
} }
save_op_s_l = (t_save_op_s *) (&(save_op[nop])); save_op_s_l = (t_save_op_s *)(&(save_op[nop]));
save_op_s_l->operat = IFSUB; save_op_s_l->operat = IFSUB;
save_op_s_l->first = first_sub->u_index; save_op_s_l->first = first_sub->u_index;
save_op_s_l->second = first_piv->u_index; save_op_s_l->second = first_piv->u_index;
@ -1912,7 +1912,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
save_op = (int *) mxRealloc(save_op, nopa*sizeof(int)); save_op = (int *) mxRealloc(save_op, nopa*sizeof(int));
} }
} }
save_op_s_l = (t_save_op_s *) (&(save_op[nop])); save_op_s_l = (t_save_op_s *)(&(save_op[nop]));
save_op_s_l->operat = IFSUB; save_op_s_l->operat = IFSUB;
save_op_s_l->first = b[row]; save_op_s_l->first = b[row];
save_op_s_l->second = b[pivj]; save_op_s_l->second = b[pivj];
@ -1920,7 +1920,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
} }
nop += 3; nop += 3;
} }
} }
mxFree(bc); mxFree(bc);
} }
if (symbolic) if (symbolic)
@ -1977,10 +1977,10 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
nop1 = nop; nop1 = nop;
} }
} }
mxFree(piv_v); mxFree(piv_v);
mxFree(pivj_v); mxFree(pivj_v);
mxFree(pivk_v); mxFree(pivk_v);
mxFree(NR); mxFree(NR);
} }
nop_all += nop; nop_all += nop;
if (symbolic) if (symbolic)

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

@ -37,35 +37,40 @@ using namespace std;
extern unsigned long _nan[2]; extern unsigned long _nan[2];
extern double NAN; extern double NAN;
inline bool isnan(double value) inline bool
isnan(double value)
{ {
return _isnan(value); return _isnan(value);
} }
inline bool isinf(double value) inline bool
isinf(double value)
{ {
return (std::numeric_limits<double>::has_infinity && return (std::numeric_limits<double>::has_infinity
value == std::numeric_limits<double>::infinity()); && value == std::numeric_limits<double>::infinity());
} }
template<typename T> template<typename T>
inline T asinh(T x) inline T
asinh(T x)
{ {
return log(x+sqrt(x*x+1)); return log(x+sqrt(x*x+1));
} }
template<typename T> template<typename T>
inline T acosh(T x) inline T
acosh(T x)
{ {
if (!(x>=1.0)) if (!(x >= 1.0))
return sqrt(-1.0); return sqrt(-1.0);
return log(x+sqrt(x*x-1.0)); return log(x+sqrt(x*x-1.0));
} }
template<typename T> template<typename T>
inline T atanh(T x) inline T
atanh(T x)
{ {
if(!(x>-1.0 && x<1.0)) if (!(x > -1.0 && x < 1.0))
return sqrt(-1.0); return sqrt(-1.0);
return log((1.0+x)/(1.0-x))/2.0; return log((1.0+x)/(1.0-x))/2.0;
} }
@ -78,108 +83,105 @@ struct t_save_op_s
int first, second; int first, second;
}; };
const int IFLD =0; const int IFLD = 0;
const int IFDIV =1; const int IFDIV = 1;
const int IFLESS=2; const int IFLESS = 2;
const int IFSUB =3; const int IFSUB = 3;
const int IFLDZ =4; const int IFLDZ = 4;
const int IFMUL =5; const int IFMUL = 5;
const int IFSTP =6; const int IFSTP = 6;
const int IFADD =7; const int IFADD = 7;
const double eps=1e-10; const double eps = 1e-10;
const double very_big=1e24; const double very_big = 1e24;
const int alt_symbolic_count_max=1; const int alt_symbolic_count_max = 1;
class SparseMatrix class SparseMatrix
{ {
public: public:
SparseMatrix(); SparseMatrix();
int simulate_NG1(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 simulate_NG1(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);
bool simulate_NG(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); bool simulate_NG(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);
void Direct_Simulate(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it, int iter); void Direct_Simulate(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it, int iter);
void fixe_u(double **u, int u_count_int, int max_lag_plus_max_lead_plus_1); void 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); void Read_SparseMatrix(string file_name, const int Size, int periods, int y_kmin, int y_kmax, bool steady_state, bool two_boundaries);
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 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);
private: private:
void Init(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int> ,int>, int> &IM); void Init(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM);
void ShortInit(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int> ,int>, int> &IM); void ShortInit(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM);
void Simple_Init(int it_, int y_kmin, int y_kmax, int Size, std::map<std::pair<std::pair<int, int> ,int>, int> &IM); void Simple_Init(int it_, int y_kmin, int y_kmax, int Size, std::map<std::pair<std::pair<int, int>, int>, int> &IM);
void End(int Size); void End(int Size);
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 #ifdef PROFILER
, long int *ndiv, long int *nsub , long int *ndiv, long int *nsub
#endif #endif
); );
void Insert(const int r, const int c, const int u_index, const int lag_index); void Insert(const int r, const int c, const int u_index, const int lag_index);
void Delete(const int r,const int c); void Delete(const int r, const int c);
int At_Row(int r, NonZeroElem **first); int At_Row(int r, NonZeroElem **first);
int At_Pos(int r, int c, NonZeroElem **first); int At_Pos(int r, int c, NonZeroElem **first);
int At_Col(int c, NonZeroElem **first); int At_Col(int c, NonZeroElem **first);
int At_Col(int c, int lag, NonZeroElem **first); int At_Col(int c, int lag, NonZeroElem **first);
int NRow(int r); int NRow(int r);
int NCol(int c); int NCol(int c);
int Union_Row(int row1, int row2); int Union_Row(int row1, int row2);
void Print(int Size,int *b); void Print(int Size, int *b);
int Get_u(); int Get_u();
void Delete_u(int pos); void Delete_u(int pos);
void Clear_u(); void Clear_u();
void Print_u(); void Print_u();
void CheckIt(int y_size, int y_kmin, int y_kmax, int Size, int periods, int iter); void CheckIt(int y_size, int y_kmin, int y_kmax, int Size, int periods, int iter);
void Check_the_Solution(int periods, int y_kmin, int y_kmax, int Size, double *u, int* pivot, int* b); 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); int complete(int beg_t, int Size, int periods, int *b);
double bksub( int tbreak, int last_period, int Size, double slowc_l double bksub(int tbreak, int last_period, int Size, double slowc_l
#ifdef PROFILER #ifdef PROFILER
, long int *nmul , long int *nmul
#endif #endif
); );
double simple_bksub(int it_, int Size, double slowc_l); double simple_bksub(int it_, int Size, double slowc_l);
stack<double> Stack; stack<double> Stack;
int nb_prologue_table_u, nb_first_table_u, nb_middle_table_u, nb_last_table_u; 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 nb_prologue_table_y, nb_first_table_y, nb_middle_table_y, nb_last_table_y;
int middle_count_loop; int middle_count_loop;
char type; char type;
fstream SaveCode; fstream SaveCode;
string filename; string filename;
int max_u, min_u; int max_u, min_u;
clock_t time00; clock_t time00;
Mem_Mngr mem_mngr; Mem_Mngr mem_mngr;
vector<int> u_liste; vector<int> u_liste;
map<pair<int, int>,NonZeroElem*> Mapped_Array; map<pair<int, int>, NonZeroElem *> Mapped_Array;
int *NbNZRow, *NbNZCol; int *NbNZRow, *NbNZCol;
NonZeroElem **FNZE_R, **FNZE_C; NonZeroElem **FNZE_R, **FNZE_C;
int nb_endo, u_count_init; int nb_endo, u_count_init;
int *pivot, *pivotk, *pivot_save; int *pivot, *pivotk, *pivot_save;
double *pivotv, *pivotva; double *pivotv, *pivotva;
int *b; int *b;
bool *line_done; bool *line_done;
bool symbolic, alt_symbolic; bool symbolic, alt_symbolic;
int alt_symbolic_count; int alt_symbolic_count;
int *g_save_op; int *g_save_op;
int first_count_loop; int first_count_loop;
int g_nop_all; int g_nop_all;
int u_count_alloc, u_count_alloc_save; int u_count_alloc, u_count_alloc_save;
double markowitz_c_s; double markowitz_c_s;
double res1a; double res1a;
long int nop_all, nop1, nop2; long int nop_all, nop1, nop2;
map<pair<pair<int, int> ,int>, int> IM_i; map<pair<pair<int, int>, int>, int> IM_i;
protected: protected:
double *u, *y, *ya; double *u, *y, *ya;
double res1, res2, max_res, max_res_idx; double res1, res2, max_res, max_res_idx;
double slowc, slowc_save, markowitz_c; double slowc, slowc_save, markowitz_c;
int y_kmin, y_kmax, y_size, periods, y_decal; int y_kmin, y_kmax, y_size, periods, y_decal;
int *index_vara, *index_equa; int *index_vara, *index_equa;
int u_count, tbreak_g; int u_count, tbreak_g;
int iter; int iter;
double *direction; double *direction;
int start_compare; int start_compare;
int restart; int restart;
bool error_not_printed; bool error_not_printed;
}; };
#endif #endif

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

@ -26,18 +26,17 @@
#include "Interpreter.hh" #include "Interpreter.hh"
#ifndef DEBUG_EX #ifndef DEBUG_EX
#include "mex.h" # include "mex.h"
#else #else
#include "mex_interface.hh" # include "mex_interface.hh"
#endif #endif
#include "Mem_Mngr.hh" #include "Mem_Mngr.hh"
#ifdef DEBUG_EX #ifdef DEBUG_EX
using namespace std; using namespace std;
#include <sstream> # include <sstream>
string string
Get_Argument(const char *argv) Get_Argument(const char *argv)
{ {
@ -45,23 +44,22 @@ Get_Argument(const char *argv)
return f; return f;
} }
int int
main( int argc, const char* argv[] ) main(int argc, const char *argv[])
{ {
FILE *fid; FILE *fid;
bool steady_state = false; bool steady_state = false;
bool evaluate = false; bool evaluate = false;
printf("argc=%d\n",argc); printf("argc=%d\n", argc);
if(argc<2) if (argc < 2)
{ {
mexPrintf("model filename expected\n"); mexPrintf("model filename expected\n");
mexEvalString("st=fclose('all');clear all;"); mexEvalString("st=fclose('all');clear all;");
mexErrMsgTxt("Exit from Dynare"); mexErrMsgTxt("Exit from Dynare");
} }
float f_tmp; float f_tmp;
ostringstream tmp_out(""); ostringstream tmp_out("");
tmp_out << argv[1] << "_options.txt"; tmp_out << argv[1] << "_options.txt";
cout << tmp_out.str().c_str() << "\n"; cout << tmp_out.str().c_str() << "\n";
int nb_params; int nb_params;
int i, row_y, col_y, row_x, col_x; int i, row_y, col_y, row_x, col_x;
@ -71,13 +69,13 @@ main( int argc, const char* argv[] )
string file_name(argv[1]); string file_name(argv[1]);
for(i=2;i<argc; i++) for (i = 2; i < argc; i++)
{ {
if(Get_Argument(argv[i])=="static") if (Get_Argument(argv[i]) == "static")
steady_state = true; steady_state = true;
else if(Get_Argument(argv[i])=="dynamic") else if (Get_Argument(argv[i]) == "dynamic")
steady_state = false; steady_state = false;
else if(Get_Argument(argv[i])=="evaluate") else if (Get_Argument(argv[i]) == "evaluate")
evaluate = true; evaluate = true;
else else
{ {
@ -89,110 +87,110 @@ main( int argc, const char* argv[] )
mexErrMsgTxt(f.c_str()); mexErrMsgTxt(f.c_str());
} }
} }
fid = fopen(tmp_out.str().c_str(),"r"); fid = fopen(tmp_out.str().c_str(), "r");
int periods; int periods;
fscanf(fid,"%d",&periods); fscanf(fid, "%d", &periods);
int maxit_; int maxit_;
fscanf(fid,"%d",&maxit_); fscanf(fid, "%d", &maxit_);
fscanf(fid,"%f",&f_tmp); fscanf(fid, "%f", &f_tmp);
double slowc = f_tmp; double slowc = f_tmp;
fscanf(fid,"%f",&f_tmp); fscanf(fid, "%f", &f_tmp);
double markowitz_c = f_tmp; double markowitz_c = f_tmp;
fscanf(fid,"%f",&f_tmp); fscanf(fid, "%f", &f_tmp);
double solve_tolf = f_tmp; double solve_tolf = f_tmp;
fscanf(fid,"%d",&minimal_solving_periods); fscanf(fid, "%d", &minimal_solving_periods);
fclose(fid); fclose(fid);
tmp_out.str(""); tmp_out.str("");
tmp_out << argv[1] << "_M.txt"; tmp_out << argv[1] << "_M.txt";
fid = fopen(tmp_out.str().c_str(),"r"); fid = fopen(tmp_out.str().c_str(), "r");
int y_kmin; int y_kmin;
fscanf(fid,"%d",&y_kmin); fscanf(fid, "%d", &y_kmin);
int y_kmax; int y_kmax;
fscanf(fid,"%d",&y_kmax); fscanf(fid, "%d", &y_kmax);
int y_decal; int y_decal;
fscanf(fid,"%d",&y_decal); fscanf(fid, "%d", &y_decal);
fscanf(fid,"%d",&nb_params); fscanf(fid, "%d", &nb_params);
fscanf(fid,"%d",&row_x); fscanf(fid, "%d", &row_x);
fscanf(fid,"%d",&col_x); fscanf(fid, "%d", &col_x);
fscanf(fid,"%d",&row_y); fscanf(fid, "%d", &row_y);
fscanf(fid,"%d",&col_y); fscanf(fid, "%d", &col_y);
int steady_row_y, steady_col_y; int steady_row_y, steady_col_y;
int steady_row_x, steady_col_x; int steady_row_x, steady_col_x;
fscanf(fid,"%d",&steady_row_y); fscanf(fid, "%d", &steady_row_y);
fscanf(fid,"%d",&steady_col_y); fscanf(fid, "%d", &steady_col_y);
fscanf(fid,"%d",&steady_row_x); fscanf(fid, "%d", &steady_row_x);
fscanf(fid,"%d",&steady_col_x); fscanf(fid, "%d", &steady_col_x);
int nb_row_xd; int nb_row_xd;
fscanf(fid,"%d",&nb_row_xd); fscanf(fid, "%d", &nb_row_xd);
double * params = (double*)malloc(nb_params*sizeof(params[0])); double *params = (double *) malloc(nb_params*sizeof(params[0]));
for(i=0; i < nb_params; i++) for (i = 0; i < nb_params; i++)
{ {
fscanf(fid,"%f",&f_tmp); fscanf(fid, "%f", &f_tmp);
params[i] = f_tmp; params[i] = f_tmp;
} }
fclose(fid); fclose(fid);
yd = (double*)malloc(row_y*col_y*sizeof(yd[0])); yd = (double *) malloc(row_y*col_y*sizeof(yd[0]));
xd = (double*)malloc(row_x*col_x*sizeof(xd[0])); xd = (double *) malloc(row_x*col_x*sizeof(xd[0]));
tmp_out.str(""); tmp_out.str("");
tmp_out << argv[1] << "_oo.txt"; tmp_out << argv[1] << "_oo.txt";
fid = fopen(tmp_out.str().c_str(),"r"); fid = fopen(tmp_out.str().c_str(), "r");
for(i=0; i < col_y*row_y; i++) for (i = 0; i < col_y*row_y; i++)
{ {
fscanf(fid,"%f",&f_tmp); fscanf(fid, "%f", &f_tmp);
yd[i] = f_tmp; yd[i] = f_tmp;
} }
for(i=0; i < col_x*row_x; i++) for (i = 0; i < col_x*row_x; i++)
{ {
fscanf(fid,"%f",&f_tmp); fscanf(fid, "%f", &f_tmp);
xd[i] = f_tmp; xd[i] = f_tmp;
} }
double *steady_yd, *steady_xd; double *steady_yd, *steady_xd;
steady_yd = (double*)malloc(steady_row_y*steady_col_y*sizeof(steady_yd[0])); steady_yd = (double *) malloc(steady_row_y*steady_col_y*sizeof(steady_yd[0]));
steady_xd = (double*)malloc(steady_row_x*steady_col_x*sizeof(steady_xd[0])); steady_xd = (double *) malloc(steady_row_x*steady_col_x*sizeof(steady_xd[0]));
for(i=0; i < steady_row_y*steady_col_y; i++) for (i = 0; i < steady_row_y*steady_col_y; i++)
{ {
fscanf(fid,"%f",&f_tmp); fscanf(fid, "%f", &f_tmp);
steady_yd[i] = f_tmp; steady_yd[i] = f_tmp;
} }
for(i=0; i < steady_row_x*steady_col_x; i++) for (i = 0; i < steady_row_x*steady_col_x; i++)
{ {
fscanf(fid,"%f",&f_tmp); fscanf(fid, "%f", &f_tmp);
steady_xd[i] = f_tmp; steady_xd[i] = f_tmp;
} }
fclose(fid); fclose(fid);
int size_of_direction=col_y*row_y*sizeof(double); int size_of_direction = col_y*row_y*sizeof(double);
double * y=(double*)mxMalloc(size_of_direction); double *y = (double *) mxMalloc(size_of_direction);
double * ya=(double*)mxMalloc(size_of_direction); double *ya = (double *) mxMalloc(size_of_direction);
direction=(double*)mxMalloc(size_of_direction); direction = (double *) mxMalloc(size_of_direction);
memset(direction,0,size_of_direction); memset(direction, 0, size_of_direction);
double * x=(double*)mxMalloc(col_x*row_x*sizeof(double)); double *x = (double *) mxMalloc(col_x*row_x*sizeof(double));
for (i=0;i<row_x*col_x;i++) for (i = 0; i < row_x*col_x; i++)
x[i]=double(xd[i]); x[i] = double (xd[i]);
for (i=0;i<row_y*col_y;i++) for (i = 0; i < row_y*col_y; i++)
y[i]=double(yd[i]); y[i] = double (yd[i]);
free(yd); free(yd);
free(xd); free(xd);
int y_size=row_y; int y_size = row_y;
int nb_row_x=row_x; int nb_row_x = row_x;
clock_t t0= clock(); 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); 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);
string f(file_name); string f(file_name);
interprete.compute_blocks(f, f, steady_state, evaluate); interprete.compute_blocks(f, f, steady_state, evaluate);
clock_t t1= clock(); clock_t t1 = clock();
if(!evaluate) if (!evaluate)
mexPrintf("Simulation Time=%f milliseconds\n",1000.0*(double(t1)-double(t0))/double(CLOCKS_PER_SEC)); mexPrintf("Simulation Time=%f milliseconds\n", 1000.0*(double (t1)-double (t0))/double (CLOCKS_PER_SEC));
if(x) if (x)
mxFree(x); mxFree(x);
if(y) if (y)
mxFree(y); mxFree(y);
if(ya) if (ya)
mxFree(ya); mxFree(ya);
if(direction) if (direction)
mxFree(direction); mxFree(direction);
free(params); free(params);
} }
@ -203,9 +201,9 @@ string
Get_Argument(const mxArray *prhs) Get_Argument(const mxArray *prhs)
{ {
const mxArray *mxa = prhs; const mxArray *mxa = prhs;
int buflen=mxGetM(mxa) * mxGetN(mxa) + 1; int buflen = mxGetM(mxa) * mxGetN(mxa) + 1;
char *first_argument; char *first_argument;
first_argument=(char*)mxCalloc(buflen, sizeof(char)); first_argument = (char *) mxCalloc(buflen, sizeof(char));
int status = mxGetString(mxa, first_argument, buflen); int status = mxGetString(mxa, first_argument, buflen);
if (status != 0) if (status != 0)
mexWarnMsgTxt("Not enough space. The first argument is truncated."); mexWarnMsgTxt("Not enough space. The first argument is truncated.");
@ -214,25 +212,25 @@ Get_Argument(const mxArray *prhs)
return f; return f;
} }
/* The gateway routine */ /* The gateway routine */
void void
mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{ {
mxArray *M_, *oo_, *options_; mxArray *M_, *oo_, *options_;
int i, row_y, col_y, row_x, col_x, nb_row_xd; int i, row_y, col_y, row_x, col_x, nb_row_xd;
int steady_row_y, steady_col_y, steady_row_x, steady_col_x, steady_nb_row_xd; int steady_row_y, steady_col_y, steady_row_x, steady_col_x, steady_nb_row_xd;
int y_kmin=0, y_kmax=0, y_decal=0, periods=1; int y_kmin = 0, y_kmax = 0, y_decal = 0, periods = 1;
double * pind ; double *pind;
double *direction; double *direction;
bool steady_state = false; bool steady_state = false;
bool evaluate = false; bool evaluate = false;
for(i=0;i<nrhs; i++) for (i = 0; i < nrhs; i++)
{ {
if(Get_Argument(prhs[i])=="static") if (Get_Argument(prhs[i]) == "static")
steady_state = true; steady_state = true;
else if(Get_Argument(prhs[i])=="dynamic") else if (Get_Argument(prhs[i]) == "dynamic")
steady_state = false; steady_state = false;
else if(Get_Argument(prhs[i])=="evaluate") else if (Get_Argument(prhs[i]) == "evaluate")
evaluate = true; evaluate = true;
else else
{ {
@ -243,130 +241,128 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
mexErrMsgTxt(f.c_str()); mexErrMsgTxt(f.c_str());
} }
} }
M_ = mexGetVariable("global","M_"); M_ = mexGetVariable("global", "M_");
if (M_ == NULL ) if (M_ == NULL)
{ {
mexPrintf("Global variable not found : "); mexPrintf("Global variable not found : ");
mexErrMsgTxt("M_ \n"); mexErrMsgTxt("M_ \n");
} }
/* Gets variables and parameters from global workspace of Matlab */ /* Gets variables and parameters from global workspace of Matlab */
oo_ = mexGetVariable("global","oo_"); oo_ = mexGetVariable("global", "oo_");
if (oo_ == NULL ) if (oo_ == NULL)
{ {
mexPrintf("Global variable not found : "); mexPrintf("Global variable not found : ");
mexErrMsgTxt("oo_ \n"); mexErrMsgTxt("oo_ \n");
} }
options_ = mexGetVariable("global","options_"); options_ = mexGetVariable("global", "options_");
if (options_ == NULL ) if (options_ == NULL)
{ {
mexPrintf("Global variable not found : "); mexPrintf("Global variable not found : ");
mexErrMsgTxt("options_ \n"); mexErrMsgTxt("options_ \n");
} }
double * params = mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,"params"))); double *params = mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "params")));
double *yd, *xd , *steady_yd = NULL, *steady_xd = NULL; double *yd, *xd, *steady_yd = NULL, *steady_xd = NULL;
if(!steady_state) if (!steady_state)
{ {
yd= mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,"endo_simul"))); yd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "endo_simul")));
row_y=mxGetM(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")));; col_y = mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "endo_simul")));;
xd= mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,"exo_simul"))); xd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_simul")));
row_x=mxGetM(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"))); col_x = mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_simul")));
nb_row_xd=int(floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,"exo_det_nbr")))))); nb_row_xd = int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "exo_det_nbr"))))));
y_kmin=int(floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,"maximum_lag")))))); 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_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"))))))); y_decal = max(0, y_kmin-int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "maximum_endo_lag")))))));
periods=int(floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,"periods")))))); periods = int (floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "periods"))))));
steady_yd= mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,"steady_state"))); steady_yd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "steady_state")));
steady_row_y=mxGetM(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")));; steady_col_y = mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "steady_state")));;
steady_xd= mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,"exo_steady_state"))); steady_xd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_steady_state")));
steady_row_x=mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,"exo_steady_state"))); steady_row_x = mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_steady_state")));
steady_col_x=mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,"exo_steady_state"))); steady_col_x = mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_steady_state")));
steady_nb_row_xd=int(floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,"exo_det_nbr")))))); steady_nb_row_xd = int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "exo_det_nbr"))))));
} }
else else
{ {
yd= mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,"steady_state"))); yd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "steady_state")));
row_y=mxGetM(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")));; col_y = mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "steady_state")));;
xd= mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,"exo_steady_state"))); xd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_steady_state")));
row_x=mxGetM(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"))); col_x = mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_steady_state")));
nb_row_xd=int(floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,"exo_det_nbr")))))); nb_row_xd = int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "exo_det_nbr"))))));
} }
int maxit_=int(floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,"maxit_")))))); int maxit_ = int (floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "maxit_"))))));
double slowc=double(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,"slowc"))))); double slowc = double (*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "slowc")))));
double markowitz_c=double(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,"markowitz"))))); 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 minimal_solving_periods = int (*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "minimal_solving_periods")))));
double solve_tolf; double solve_tolf;
if(steady_state) if (steady_state)
solve_tolf=*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,"solve_tolf")))); solve_tolf = *(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "solve_tolf"))));
else else
solve_tolf=*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,"dynatol")))); solve_tolf = *(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "dynatol"))));
mxArray *mxa=mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,"fname")); mxArray *mxa = mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "fname"));
int buflen=mxGetM(mxa) * mxGetN(mxa) + 1; int buflen = mxGetM(mxa) * mxGetN(mxa) + 1;
char *fname; char *fname;
fname=(char*)mxCalloc(buflen, sizeof(char)); fname = (char *) mxCalloc(buflen, sizeof(char));
string file_name=fname; string file_name = fname;
int status = mxGetString(mxa, fname, buflen); int status = mxGetString(mxa, fname, buflen);
if (status != 0) if (status != 0)
mexWarnMsgTxt("Not enough space. Filename is truncated."); mexWarnMsgTxt("Not enough space. Filename is truncated.");
int size_of_direction = col_y*row_y*sizeof(double);
int size_of_direction=col_y*row_y*sizeof(double); double *y = (double *) mxMalloc(size_of_direction);
double * y=(double*)mxMalloc(size_of_direction); double *ya = (double *) mxMalloc(size_of_direction);
double * ya=(double*)mxMalloc(size_of_direction); direction = (double *) mxMalloc(size_of_direction);
direction=(double*)mxMalloc(size_of_direction); memset(direction, 0, size_of_direction);
memset(direction,0,size_of_direction); double *x = (double *) mxMalloc(col_x*row_x*sizeof(double));
double * x=(double*)mxMalloc(col_x*row_x*sizeof(double)); for (i = 0; i < row_x*col_x; i++)
for (i=0;i<row_x*col_x;i++) x[i] = double (xd[i]);
x[i]=double(xd[i]); for (i = 0; i < row_y*col_y; i++)
for (i=0;i<row_y*col_y;i++)
{ {
y[i] = double(yd[i]); y[i] = double (yd[i]);
ya[i] = double(yd[i]); ya[i] = double (yd[i]);
} }
int y_size=row_y; int y_size = row_y;
int nb_row_x=row_x; int nb_row_x = row_x;
clock_t t0 = clock();
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); 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);
string f(fname); string f(fname);
bool result = interprete.compute_blocks(f, f, steady_state, evaluate); bool result = interprete.compute_blocks(f, f, steady_state, evaluate);
clock_t t1= clock(); clock_t t1 = clock();
if(!steady_state && !evaluate) if (!steady_state && !evaluate)
mexPrintf("Simulation Time=%f milliseconds\n",1000.0*(double(t1)-double(t0))/double(CLOCKS_PER_SEC)); mexPrintf("Simulation Time=%f milliseconds\n", 1000.0*(double (t1)-double (t0))/double (CLOCKS_PER_SEC));
if (nlhs>0) if (nlhs > 0)
{ {
plhs[0] = mxCreateDoubleMatrix(row_y, col_y, mxREAL); plhs[0] = mxCreateDoubleMatrix(row_y, col_y, mxREAL);
pind = mxGetPr(plhs[0]); pind = mxGetPr(plhs[0]);
if(evaluate) if (evaluate)
for (i=0;i<row_y*col_y;i++) for (i = 0; i < row_y*col_y; i++)
pind[i]=y[i]-ya[i]; pind[i] = y[i]-ya[i];
else else
for (i=0;i<row_y*col_y;i++) for (i = 0; i < row_y*col_y; i++)
pind[i]=y[i]; pind[i] = y[i];
if(nlhs>1) if (nlhs > 1)
{ {
plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL); plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
pind = mxGetPr(plhs[1]); pind = mxGetPr(plhs[1]);
if(result) if (result)
pind[0] = 0; pind[0] = 0;
else else
pind[0] = 1; pind[0] = 1;
} }
} }
if(x) if (x)
mxFree(x); mxFree(x);
if(y) if (y)
mxFree(y); mxFree(y);
if(ya) if (ya)
mxFree(ya); mxFree(ya);
if(direction) if (direction)
mxFree(direction); mxFree(direction);
} }
#endif #endif

15
mex/sources/bytecode/testing/mex_interface.cc
View File

@ -10,9 +10,9 @@ mexPrintf(const char *str, ...)
va_list args; va_list args;
int retval; int retval;
va_start (args, str); va_start(args, str);
retval = vprintf (str, args); retval = vprintf(str, args);
va_end (args); va_end(args);
return retval; return retval;
} }
@ -25,24 +25,23 @@ mexErrMsgTxt(const string str)
} }
void void
mxFree(void* to_release) mxFree(void *to_release)
{ {
free(to_release); free(to_release);
} }
void* void *
mxMalloc(int amount) mxMalloc(int amount)
{ {
return malloc(amount); return malloc(amount);
} }
void* void *
mxRealloc(void* to_extend, int amount) mxRealloc(void *to_extend, int amount)
{ {
return realloc(to_extend, amount); return realloc(to_extend, amount);
} }
void void
mexEvalString(const string str) mexEvalString(const string str)
{ {

8
mex/sources/bytecode/testing/mex_interface.hh
View File

@ -5,9 +5,9 @@
#include <stdarg.h> #include <stdarg.h>
using namespace std; using namespace std;
int mexPrintf(/*const string*/const char* str, ...); int mexPrintf(/*const string*/ const char *str, ...);
void mexErrMsgTxt(const string str); void mexErrMsgTxt(const string str);
void* mxMalloc(int amount); void *mxMalloc(int amount);
void* mxRealloc(void* to_extend, int amount); void *mxRealloc(void *to_extend, int amount);
void mxFree(void* to_release); void mxFree(void *to_release);
void mexEvalString(const string str); void mexEvalString(const string str);

73
mex/sources/bytecode/testing/simulate_debug.m
View File

@ -1,37 +1,36 @@
function simulate_debug() function simulate_debug()
global M_ oo_ options_; global M_ oo_ options_;
fid = fopen([M_.fname '_options.txt'],'wt'); fid = fopen([M_.fname '_options.txt'],'wt');
fprintf(fid,'%d\n',options_.periods); fprintf(fid,'%d\n',options_.periods);
fprintf(fid,'%d\n',options_.maxit_); fprintf(fid,'%d\n',options_.maxit_);
fprintf(fid,'%6.20f\n',options_.slowc); fprintf(fid,'%6.20f\n',options_.slowc);
fprintf(fid,'%6.20f\n',options_.markowitz); fprintf(fid,'%6.20f\n',options_.markowitz);
fprintf(fid,'%6.20f\n',options_.dynatol); fprintf(fid,'%6.20f\n',options_.dynatol);
fprintf(fid,'%d\n',options_.minimal_solving_periods); fprintf(fid,'%d\n',options_.minimal_solving_periods);
fclose(fid); fclose(fid);
fid = fopen([M_.fname '_M.txt'],'wt'); fid = fopen([M_.fname '_M.txt'],'wt');
fprintf(fid,'%d\n',M_.maximum_lag); fprintf(fid,'%d\n',M_.maximum_lag);
fprintf(fid,'%d\n',M_.maximum_lead); fprintf(fid,'%d\n',M_.maximum_lead);
fprintf(fid,'%d\n',M_.maximum_endo_lag); fprintf(fid,'%d\n',M_.maximum_endo_lag);
fprintf(fid,'%d\n',M_.param_nbr); fprintf(fid,'%d\n',M_.param_nbr);
fprintf(fid,'%d\n',size(oo_.exo_simul, 1)); fprintf(fid,'%d\n',size(oo_.exo_simul, 1));
fprintf(fid,'%d\n',size(oo_.exo_simul, 2)); fprintf(fid,'%d\n',size(oo_.exo_simul, 2));
fprintf(fid,'%d\n',M_.endo_nbr); fprintf(fid,'%d\n',M_.endo_nbr);
fprintf(fid,'%d\n',size(oo_.endo_simul, 2)); fprintf(fid,'%d\n',size(oo_.endo_simul, 2));
fprintf(fid,'%d\n',M_.exo_det_nbr); fprintf(fid,'%d\n',M_.exo_det_nbr);
fprintf(fid,'%d\n',size(oo_.steady_state,1)); fprintf(fid,'%d\n',size(oo_.steady_state,1));
fprintf(fid,'%d\n',size(oo_.steady_state,2)); fprintf(fid,'%d\n',size(oo_.steady_state,2));
fprintf(fid,'%d\n',size(oo_.exo_steady_state,1)); fprintf(fid,'%d\n',size(oo_.exo_steady_state,1));
fprintf(fid,'%d\n',size(oo_.exo_steady_state,2)); fprintf(fid,'%d\n',size(oo_.exo_steady_state,2));
fprintf(fid,'%6.20f\n',M_.params); fprintf(fid,'%6.20f\n',M_.params);
fclose(fid); fclose(fid);
fid = fopen([M_.fname '_oo.txt'],'wt'); fid = fopen([M_.fname '_oo.txt'],'wt');
fprintf(fid,'%6.20f\n',oo_.endo_simul); fprintf(fid,'%6.20f\n',oo_.endo_simul);
fprintf(fid,'%6.20f\n',oo_.exo_simul); fprintf(fid,'%6.20f\n',oo_.exo_simul);
fprintf(fid,'%6.20f\n',oo_.steady_state); fprintf(fid,'%6.20f\n',oo_.steady_state);
fprintf(fid,'%6.20f\n',oo_.exo_steady_state); fprintf(fid,'%6.20f\n',oo_.exo_steady_state);
fclose(fid); fclose(fid);

60
mex/sources/dynblas.h
View File

@ -30,77 +30,77 @@
#define _DYNBLAS_H #define _DYNBLAS_H
/* Starting from version 7.8, MATLAB BLAS expects ptrdiff_t arguments for integers */ /* Starting from version 7.8, MATLAB BLAS expects ptrdiff_t arguments for integers */
# if defined(MATLAB_MEX_FILE) && MATLAB_VERSION >= 0x0708 #if defined(MATLAB_MEX_FILE) && MATLAB_VERSION >= 0x0708
# ifdef __cplusplus
# include <cstdlib>
# else
# include <stdlib.h>
# endif
typedef ptrdiff_t blas_int;
# else
typedef int blas_int;
# endif
# if defined(MATLAB_MEX_FILE) && defined(_WIN32)
# define FORTRAN_WRAPPER(x) x
# else
# define FORTRAN_WRAPPER(x) x ## _
# endif
# ifdef __cplusplus # ifdef __cplusplus
extern "C" { # include <cstdlib>
# else
# include <stdlib.h>
# endif # endif
typedef ptrdiff_t blas_int;
#else
typedef int blas_int;
#endif
#if defined(MATLAB_MEX_FILE) && defined(_WIN32)
# define FORTRAN_WRAPPER(x) x
#else
# define FORTRAN_WRAPPER(x) x ## _
#endif
#ifdef __cplusplus
extern "C" {
#endif
typedef const char *BLCHAR; typedef const char *BLCHAR;
typedef const blas_int *CONST_BLINT; typedef const blas_int *CONST_BLINT;
typedef const double *CONST_BLDOU; typedef const double *CONST_BLDOU;
typedef double *BLDOU; typedef double *BLDOU;
# define dgemm FORTRAN_WRAPPER(dgemm) #define dgemm FORTRAN_WRAPPER(dgemm)
void dgemm(BLCHAR transa, BLCHAR transb, CONST_BLINT m, CONST_BLINT n, void dgemm(BLCHAR transa, BLCHAR transb, CONST_BLINT m, CONST_BLINT n,
CONST_BLINT k, CONST_BLDOU alpha, CONST_BLDOU a, CONST_BLINT lda, CONST_BLINT k, CONST_BLDOU alpha, CONST_BLDOU a, CONST_BLINT lda,
CONST_BLDOU b, CONST_BLINT ldb, CONST_BLDOU beta, CONST_BLDOU b, CONST_BLINT ldb, CONST_BLDOU beta,
BLDOU c, CONST_BLINT ldc); BLDOU c, CONST_BLINT ldc);
# define dgemv FORTRAN_WRAPPER(dgemv) #define dgemv FORTRAN_WRAPPER(dgemv)
void dgemv(BLCHAR trans, CONST_BLINT m, CONST_BLINT n, CONST_BLDOU alpha, void dgemv(BLCHAR trans, CONST_BLINT m, CONST_BLINT n, CONST_BLDOU alpha,
CONST_BLDOU a, CONST_BLINT lda, CONST_BLDOU x, CONST_BLINT incx, CONST_BLDOU a, CONST_BLINT lda, CONST_BLDOU x, CONST_BLINT incx,
CONST_BLDOU beta, BLDOU y, CONST_BLINT incy); CONST_BLDOU beta, BLDOU y, CONST_BLINT incy);
# define dtrsv FORTRAN_WRAPPER(dtrsv) #define dtrsv FORTRAN_WRAPPER(dtrsv)
void dtrsv(BLCHAR uplo, BLCHAR trans, BLCHAR diag, CONST_BLINT n, void dtrsv(BLCHAR uplo, BLCHAR trans, BLCHAR diag, CONST_BLINT n,
CONST_BLDOU a, CONST_BLINT lda, BLDOU x, CONST_BLINT incx); CONST_BLDOU a, CONST_BLINT lda, BLDOU x, CONST_BLINT incx);
# define dtrmv FORTRAN_WRAPPER(dtrmv) #define dtrmv FORTRAN_WRAPPER(dtrmv)
void dtrmv(BLCHAR uplo, BLCHAR trans, BLCHAR diag, CONST_BLINT n, void dtrmv(BLCHAR uplo, BLCHAR trans, BLCHAR diag, CONST_BLINT n,
CONST_BLDOU a, CONST_BLINT lda, BLDOU x, CONST_BLINT incx); CONST_BLDOU a, CONST_BLINT lda, BLDOU x, CONST_BLINT incx);
# define daxpy FORTRAN_WRAPPER(daxpy) #define daxpy FORTRAN_WRAPPER(daxpy)
void daxpy(CONST_BLINT n, CONST_BLDOU a, CONST_BLDOU x, CONST_BLINT incx, void daxpy(CONST_BLINT n, CONST_BLDOU a, CONST_BLDOU x, CONST_BLINT incx,
BLDOU y, CONST_BLINT incy); BLDOU y, CONST_BLINT incy);
# define dcopy FORTRAN_WRAPPER(dcopy) #define dcopy FORTRAN_WRAPPER(dcopy)
void dcopy(CONST_BLINT n, CONST_BLDOU x, CONST_BLINT incx, void dcopy(CONST_BLINT n, CONST_BLDOU x, CONST_BLINT incx,
BLDOU y, CONST_BLINT incy); BLDOU y, CONST_BLINT incy);
# define zaxpy FORTRAN_WRAPPER(zaxpy) #define zaxpy FORTRAN_WRAPPER(zaxpy)
void zaxpy(CONST_BLINT n, CONST_BLDOU a, CONST_BLDOU x, CONST_BLINT incx, void zaxpy(CONST_BLINT n, CONST_BLDOU a, CONST_BLDOU x, CONST_BLINT incx,
BLDOU y, CONST_BLINT incy); BLDOU y, CONST_BLINT incy);
# define dscal FORTRAN_WRAPPER(dscal) #define dscal FORTRAN_WRAPPER(dscal)
void dscal(CONST_BLINT n, CONST_BLDOU a, BLDOU x, CONST_BLINT incx); void dscal(CONST_BLINT n, CONST_BLDOU a, BLDOU x, CONST_BLINT incx);
# define dtrsm FORTRAN_WRAPPER(dtrsm) #define dtrsm FORTRAN_WRAPPER(dtrsm)
void dtrsm(BLCHAR side, BLCHAR uplo, BLCHAR transa, BLCHAR diag, CONST_BLINT m, void dtrsm(BLCHAR side, BLCHAR uplo, BLCHAR transa, BLCHAR diag, CONST_BLINT m,
CONST_BLINT n, CONST_BLDOU alpha, CONST_BLDOU a, CONST_BLINT lda, CONST_BLINT n, CONST_BLDOU alpha, CONST_BLDOU a, CONST_BLINT lda,
BLDOU b, CONST_BLINT ldb); BLDOU b, CONST_BLINT ldb);
# define ddot FORTRAN_WRAPPER(ddot) #define ddot FORTRAN_WRAPPER(ddot)
double ddot(CONST_BLINT n, CONST_BLDOU x, CONST_BLINT incx, CONST_BLDOU y, double ddot(CONST_BLINT n, CONST_BLDOU x, CONST_BLINT incx, CONST_BLDOU y,
CONST_BLINT incy); CONST_BLINT incy);
# ifdef __cplusplus #ifdef __cplusplus
} /* extern "C" */ } /* extern "C" */
# endif #endif
#endif /* _DYNBLAS_H */ #endif /* _DYNBLAS_H */

70
mex/sources/dynlapack.h
View File

@ -6,7 +6,7 @@
* *
* When used in the context of a MATLAB MEX file, you must define MATLAB_MEX_FILE * When used in the context of a MATLAB MEX file, you must define MATLAB_MEX_FILE
* and MATLAB_VERSION (for version 7.4, define it to 0x0704). * and MATLAB_VERSION (for version 7.4, define it to 0x0704).
* *
* *
* Copyright (C) 2009 Dynare Team * Copyright (C) 2009 Dynare Team
* *
@ -30,88 +30,88 @@
#define _DYNLAPACK_H #define _DYNLAPACK_H
/* Starting from version 7.8, MATLAB LAPACK expects ptrdiff_t arguments for integers */ /* Starting from version 7.8, MATLAB LAPACK expects ptrdiff_t arguments for integers */
# if defined(MATLAB_MEX_FILE) && MATLAB_VERSION >= 0x0708 #if defined(MATLAB_MEX_FILE) && MATLAB_VERSION >= 0x0708
# ifdef __cplusplus
# include <cstdlib>
# else
# include <stdlib.h>
# endif
typedef ptrdiff_t lapack_int;
# else
typedef int lapack_int;
# endif
# if defined(MATLAB_MEX_FILE) && defined(_WIN32)
# define FORTRAN_WRAPPER(x) x
# else
# define FORTRAN_WRAPPER(x) x ## _
# endif
# ifdef __cplusplus # ifdef __cplusplus
extern "C" { # include <cstdlib>
# else
# include <stdlib.h>
# endif # endif
typedef ptrdiff_t lapack_int;
#else
typedef int lapack_int;
#endif
#if defined(MATLAB_MEX_FILE) && defined(_WIN32)
# define FORTRAN_WRAPPER(x) x
#else
# define FORTRAN_WRAPPER(x) x ## _
#endif
#ifdef __cplusplus
extern "C" {
#endif
typedef const char *LACHAR; typedef const char *LACHAR;
typedef const lapack_int *CONST_LAINT; typedef const lapack_int *CONST_LAINT;
typedef lapack_int *LAINT; typedef lapack_int *LAINT;
typedef const double *CONST_LADOU; typedef const double *CONST_LADOU;
typedef double *LADOU; typedef double *LADOU;
typedef lapack_int (*DGGESCRIT)(const double*, const double*, const double*); typedef lapack_int (*DGGESCRIT)(const double *, const double *, const double *);
# define dgetrs FORTRAN_WRAPPER(dgetrs) #define dgetrs FORTRAN_WRAPPER(dgetrs)
void dgetrs(LACHAR trans, CONST_LAINT n, CONST_LAINT nrhs, CONST_LADOU a, CONST_LAINT lda, CONST_LAINT ipiv, void dgetrs(LACHAR trans, CONST_LAINT n, CONST_LAINT nrhs, CONST_LADOU a, CONST_LAINT lda, CONST_LAINT ipiv,
LADOU b, CONST_LAINT ldb, LAINT info); LADOU b, CONST_LAINT ldb, LAINT info);
# define dgetrf FORTRAN_WRAPPER(dgetrf) #define dgetrf FORTRAN_WRAPPER(dgetrf)
void dgetrf(CONST_LAINT m, CONST_LAINT n, LADOU a, void dgetrf(CONST_LAINT m, CONST_LAINT n, LADOU a,
CONST_LAINT lda, LAINT ipiv, LAINT info); CONST_LAINT lda, LAINT ipiv, LAINT info);
# define dgees FORTRAN_WRAPPER(dgees) #define dgees FORTRAN_WRAPPER(dgees)
void dgees(LACHAR jobvs, LACHAR sort, const void* select, void dgees(LACHAR jobvs, LACHAR sort, const void *select,
CONST_LAINT n, LADOU a, CONST_LAINT lda, LAINT sdim, CONST_LAINT n, LADOU a, CONST_LAINT lda, LAINT sdim,
LADOU wr, LADOU wi, LADOU vs, CONST_LAINT ldvs, LADOU wr, LADOU wi, LADOU vs, CONST_LAINT ldvs,
LADOU work, CONST_LAINT lwork, const void* bwork, LAINT info); LADOU work, CONST_LAINT lwork, const void *bwork, LAINT info);
# define dgecon FORTRAN_WRAPPER(dgecon) #define dgecon FORTRAN_WRAPPER(dgecon)
void dgecon(LACHAR norm, CONST_LAINT n, CONST_LADOU a, CONST_LAINT lda, void dgecon(LACHAR norm, CONST_LAINT n, CONST_LADOU a, CONST_LAINT lda,
CONST_LADOU anorm, LADOU rnorm, LADOU work, LAINT iwork, CONST_LADOU anorm, LADOU rnorm, LADOU work, LAINT iwork,
LAINT info); LAINT info);
# define dtrexc FORTRAN_WRAPPER(dtrexc) #define dtrexc FORTRAN_WRAPPER(dtrexc)
void dtrexc(LACHAR compq, CONST_LAINT n, LADOU t, CONST_LAINT ldt, void dtrexc(LACHAR compq, CONST_LAINT n, LADOU t, CONST_LAINT ldt,
LADOU q, CONST_LAINT ldq, LAINT ifst, LAINT ilst, LADOU work, LADOU q, CONST_LAINT ldq, LAINT ifst, LAINT ilst, LADOU work,
LAINT info); LAINT info);
# define dtrsyl FORTRAN_WRAPPER(dtrsyl) #define dtrsyl FORTRAN_WRAPPER(dtrsyl)
void dtrsyl(LACHAR trana, LACHAR tranb, CONST_LAINT isgn, CONST_LAINT m, void dtrsyl(LACHAR trana, LACHAR tranb, CONST_LAINT isgn, CONST_LAINT m,
CONST_LAINT n, CONST_LADOU a, CONST_LAINT lda, CONST_LADOU b, CONST_LAINT n, CONST_LADOU a, CONST_LAINT lda, CONST_LADOU b,
CONST_LAINT ldb, LADOU c, CONST_LAINT ldc, LADOU scale, CONST_LAINT ldb, LADOU c, CONST_LAINT ldc, LADOU scale,
LAINT info); LAINT info);
# define dpotrf FORTRAN_WRAPPER(dpotrf) #define dpotrf FORTRAN_WRAPPER(dpotrf)
void dpotrf(LACHAR uplo, CONST_LAINT n, LADOU a, CONST_LAINT lda, void dpotrf(LACHAR uplo, CONST_LAINT n, LADOU a, CONST_LAINT lda,
LAINT info); LAINT info);
# define dgges FORTRAN_WRAPPER(dgges) #define dgges FORTRAN_WRAPPER(dgges)
void dgges(LACHAR jobvsl, LACHAR jobvsr, LACHAR sort, DGGESCRIT delztg, void dgges(LACHAR jobvsl, LACHAR jobvsr, LACHAR sort, DGGESCRIT delztg,
CONST_LAINT n, LADOU a, CONST_LAINT lda, LADOU b, CONST_LAINT ldb, CONST_LAINT n, LADOU a, CONST_LAINT lda, LADOU b, CONST_LAINT ldb,
LAINT sdim, LADOU alphar, LADOU alphai, LADOU beta, LAINT sdim, LADOU alphar, LADOU alphai, LADOU beta,
LADOU vsl, CONST_LAINT ldvsl, LADOU vsr, CONST_LAINT ldvsr, LADOU vsl, CONST_LAINT ldvsl, LADOU vsr, CONST_LAINT ldvsr,
LADOU work, CONST_LAINT lwork, LAINT bwork, LAINT info); LADOU work, CONST_LAINT lwork, LAINT bwork, LAINT info);
# define dsyev FORTRAN_WRAPPER(dsyev) #define dsyev FORTRAN_WRAPPER(dsyev)
void dsyev(LACHAR jobz, LACHAR uplo, CONST_LAINT n, LADOU a, CONST_LAINT lda, void dsyev(LACHAR jobz, LACHAR uplo, CONST_LAINT n, LADOU a, CONST_LAINT lda,
LADOU w, LADOU work, CONST_LAINT lwork, LAINT info); LADOU w, LADOU work, CONST_LAINT lwork, LAINT info);
# define dsyevr FORTRAN_WRAPPER(dsyevr) #define dsyevr FORTRAN_WRAPPER(dsyevr)
void dsyevr(LACHAR jobz, LACHAR range, LACHAR uplo, CONST_LAINT n, LADOU a, void dsyevr(LACHAR jobz, LACHAR range, LACHAR uplo, CONST_LAINT n, LADOU a,
CONST_LAINT lda, LADOU lv, LADOU vu, CONST_LAINT il, CONST_LAINT iu, CONST_LAINT lda, LADOU lv, LADOU vu, CONST_LAINT il, CONST_LAINT iu,
CONST_LADOU abstol, LAINT m, LADOU w, LADOU z, CONST_LAINT ldz, CONST_LADOU abstol, LAINT m, LADOU w, LADOU z, CONST_LAINT ldz,
LAINT isuppz, LADOU work, CONST_LAINT lwork, LAINT iwork, CONST_LAINT liwork, LAINT isuppz, LADOU work, CONST_LAINT lwork, LAINT iwork, CONST_LAINT liwork,
LAINT info); LAINT info);
# ifdef __cplusplus #ifdef __cplusplus
} /* extern "C" */ } /* extern "C" */
# endif #endif
#endif /* _DYNLAPACK_H */ #endif /* _DYNLAPACK_H */

12
mex/sources/dynmex.h
View File

@ -20,16 +20,16 @@
#ifndef _DYNMEX_H #ifndef _DYNMEX_H
#define _DYNMEX_H #define _DYNMEX_H
# if !defined(MATLAB_MEX_FILE) && !defined(OCTAVE_MEX_FILE) #if !defined(MATLAB_MEX_FILE) && !defined(OCTAVE_MEX_FILE)
# error You must define either MATLAB_MEX_FILE or OCTAVE_MEX_FILE # error You must define either MATLAB_MEX_FILE or OCTAVE_MEX_FILE
# endif #endif
# include <mex.h> #include <mex.h>
/* mwSize, mwIndex and mwSignedIndex appeared in MATLAB 7.3 */ /* mwSize, mwIndex and mwSignedIndex appeared in MATLAB 7.3 */
# if defined(MATLAB_MEX_FILE) && MATLAB_VERSION < 0x0703 #if defined(MATLAB_MEX_FILE) && MATLAB_VERSION < 0x0703
typedef int mwIndex; typedef int mwIndex;
typedef int mwSize; typedef int mwSize;
# endif #endif
#endif #endif

22
mex/sources/k_order_perturbation/dynamic_dll.cpp
View File

@ -27,8 +27,8 @@
* <model>_dynamic () function * <model>_dynamic () function
**************************************/ **************************************/
DynamicModelDLL::DynamicModelDLL(const string &modName, const int y_length, const int j_cols, DynamicModelDLL::DynamicModelDLL(const string &modName, const int y_length, const int j_cols,
const int n_max_lag, const int n_exog, const string &sExt) throw (DynareException) const int n_max_lag, const int n_exog, const string &sExt) throw (DynareException) :
: length(y_length), jcols(j_cols), nMax_lag(n_max_lag), nExog(n_exog) length(y_length), jcols(j_cols), nMax_lag(n_max_lag), nExog(n_exog)
{ {
string fName; string fName;
#if !defined(__CYGWIN32__) && !defined(_WIN32) #if !defined(__CYGWIN32__) && !defined(_WIN32)
@ -42,12 +42,12 @@ DynamicModelDLL::DynamicModelDLL(const string &modName, const int y_length, cons
dynamicHinstance = LoadLibrary(fName.c_str()); dynamicHinstance = LoadLibrary(fName.c_str());
if (dynamicHinstance == NULL) if (dynamicHinstance == NULL)
throw 1; throw 1;
Dynamic = (DynamicFn)GetProcAddress(dynamicHinstance, "Dynamic"); Dynamic = (DynamicFn) GetProcAddress(dynamicHinstance, "Dynamic");
if (Dynamic == NULL) if (Dynamic == NULL)
{ {
FreeLibrary(dynamicHinstance); // Free the library FreeLibrary(dynamicHinstance); // Free the library
throw 2; throw 2;
} }
#else // Linux or Mac #else // Linux or Mac
dynamicHinstance = dlopen(fName.c_str(), RTLD_NOW); dynamicHinstance = dlopen(fName.c_str(), RTLD_NOW);
if ((dynamicHinstance == NULL) || dlerror()) if ((dynamicHinstance == NULL) || dlerror())
@ -58,7 +58,7 @@ DynamicModelDLL::DynamicModelDLL(const string &modName, const int y_length, cons
Dynamic = (DynamicFn) dlsym(dynamicHinstance, "Dynamic"); Dynamic = (DynamicFn) dlsym(dynamicHinstance, "Dynamic");
if ((Dynamic == NULL) || dlerror()) if ((Dynamic == NULL) || dlerror())
{ {
dlclose(dynamicHinstance); // Free the library dlclose(dynamicHinstance); // Free the library
cerr << dlerror() << endl; cerr << dlerror() << endl;
throw 2; throw 2;
} }
@ -120,16 +120,16 @@ DynamicModelDLL::eval(const Vector &y, const TwoDMatrix &x, const Vector *modPa
dg1 = const_cast<double *>(g1->base()); dg1 = const_cast<double *>(g1->base());
} }
if (g2 != NULL) if (g2 != NULL)
dg2 = const_cast<double *>(g2->base()); dg2 = const_cast<double *>(g2->base());
dresidual = const_cast<double *>(residual.base()); dresidual = const_cast<double *>(residual.base());
if (g3 != NULL) if (g3 != NULL)
dg3 = const_cast<double *>(g3->base()); dg3 = const_cast<double *>(g3->base());
dresidual = const_cast<double *>(residual.base()); dresidual = const_cast<double *>(residual.base());
double *dy = const_cast<double *>(y.base()); double *dy = const_cast<double *>(y.base());
double *dx = const_cast<double *>(x.base()); double *dx = const_cast<double *>(x.base());
double *dbParams = const_cast<double *>(modParams->base()); double *dbParams = const_cast<double *>(modParams->base());
Dynamic(dy, dx, nExog, dbParams, it_, dresidual, dg1, dg2, dg3); Dynamic(dy, dx, nExog, dbParams, it_, dresidual, dg1, dg2, dg3);
} }
void void

326
mex/sources/k_order_perturbation/k_ord_dynare.cpp
View File

@ -36,35 +36,35 @@
KordpDynare::KordpDynare(const char **endo, int num_endo, KordpDynare::KordpDynare(const char **endo, int num_endo,
const char **exo, int nexog, int npar, //const char** par, const char **exo, int nexog, int npar, //const char** par,
Vector *ysteady, TwoDMatrix *vcov, Vector *inParams, int nstat, Vector *ysteady, TwoDMatrix *vcov, Vector *inParams, int nstat,
int npred, int nforw, int nboth, const int jcols, const Vector *nnzd, int npred, int nforw, int nboth, const int jcols, const Vector *nnzd,
const int nsteps, int norder, //const char* modName, const int nsteps, int norder, //const char* modName,
Journal &jr, DynamicModelDLL &dynamicDLL, double sstol, Journal &jr, DynamicModelDLL &dynamicDLL, double sstol,
const vector<int> *var_order, const TwoDMatrix *llincidence, double criterium) throw (TLException) const vector<int> *var_order, const TwoDMatrix *llincidence, double criterium) throw (TLException) :
: nStat(nstat), nBoth(nboth), nPred(npred), nForw(nforw), nExog(nexog), nPar(npar), nStat(nstat), nBoth(nboth), nPred(npred), nForw(nforw), nExog(nexog), nPar(npar),
nYs(npred + nboth), nYss(nboth + nforw), nY(num_endo), nJcols(jcols), NNZD(nnzd), nSteps(nsteps), nYs(npred + nboth), nYss(nboth + nforw), nY(num_endo), nJcols(jcols), NNZD(nnzd), nSteps(nsteps),
nOrder(norder), journal(jr), dynamicDLL(dynamicDLL), ySteady(ysteady), vCov(vcov), params(inParams), nOrder(norder), journal(jr), dynamicDLL(dynamicDLL), ySteady(ysteady), vCov(vcov), params(inParams),
md(1), dnl(NULL), denl(NULL), dsnl(NULL), ss_tol(sstol), varOrder(var_order), md(1), dnl(NULL), denl(NULL), dsnl(NULL), ss_tol(sstol), varOrder(var_order),
ll_Incidence(llincidence), qz_criterium(criterium) ll_Incidence(llincidence), qz_criterium(criterium)
{ {
dnl = new DynareNameList(*this, endo); dnl = new DynareNameList(*this, endo);
denl = new DynareExogNameList(*this, exo); denl = new DynareExogNameList(*this, exo);
dsnl = new DynareStateNameList(*this, *dnl, *denl); dsnl = new DynareStateNameList(*this, *dnl, *denl);
JacobianIndices = ReorderDynareJacobianIndices(varOrder); JacobianIndices = ReorderDynareJacobianIndices(varOrder);
// Initialise ModelDerivativeContainer(*this, this->md, nOrder); // Initialise ModelDerivativeContainer(*this, this->md, nOrder);
for (int iord = 1; iord <= nOrder; iord++) for (int iord = 1; iord <= nOrder; iord++)
{ {
FSSparseTensor *t = new FSSparseTensor(iord, nY+nYs+nYss+nExog, nY); FSSparseTensor *t = new FSSparseTensor(iord, nY+nYs+nYss+nExog, nY);
md.insert(t); md.insert(t);
} }
} }
KordpDynare::KordpDynare(const KordpDynare &dynare) KordpDynare::KordpDynare(const KordpDynare &dynare) :
: nStat(dynare.nStat), nBoth(dynare.nBoth), nPred(dynare.nPred), nStat(dynare.nStat), nBoth(dynare.nBoth), nPred(dynare.nPred),
nForw(dynare.nForw), nExog(dynare.nExog), nPar(dynare.nPar), nForw(dynare.nForw), nExog(dynare.nExog), nPar(dynare.nPar),
nYs(dynare.nYs), nYss(dynare.nYss), nY(dynare.nY), nJcols(dynare.nJcols), nYs(dynare.nYs), nYss(dynare.nYss), nY(dynare.nY), nJcols(dynare.nJcols),
NNZD(dynare.NNZD), nSteps(dynare.nSteps), nOrder(dynare.nOrder), NNZD(dynare.NNZD), nSteps(dynare.nSteps), nOrder(dynare.nOrder),
journal(dynare.journal), dynamicDLL(dynare.dynamicDLL), journal(dynare.journal), dynamicDLL(dynare.dynamicDLL),
ySteady(NULL), params(NULL), vCov(NULL), md(dynare.md), ySteady(NULL), params(NULL), vCov(NULL), md(dynare.md),
dnl(NULL), denl(NULL), dsnl(NULL), ss_tol(dynare.ss_tol), dnl(NULL), denl(NULL), dsnl(NULL), ss_tol(dynare.ss_tol),
@ -100,7 +100,7 @@ KordpDynare::~KordpDynare()
if (ll_Incidence) if (ll_Incidence)
delete ll_Incidence; delete ll_Incidence;
if (NNZD) if (NNZD)
delete NNZD; delete NNZD;
} }
/** This clears the container of model derivatives and initializes it /** This clears the container of model derivatives and initializes it
@ -143,8 +143,8 @@ KordpDynare::evaluateSystem(Vector &out, const Vector &yym, const Vector &yy,
void void
KordpDynare::calcDerivatives(const Vector &yy, const Vector &xx) throw (DynareException) KordpDynare::calcDerivatives(const Vector &yy, const Vector &xx) throw (DynareException)
{ {
TwoDMatrix *g2 = 0;// NULL; TwoDMatrix *g2 = 0; // NULL;
TwoDMatrix *g3 = 0;// NULL; TwoDMatrix *g3 = 0; // NULL;
TwoDMatrix *g1 = new TwoDMatrix(nY, nJcols); // generate g1 for jacobian TwoDMatrix *g1 = new TwoDMatrix(nY, nJcols); // generate g1 for jacobian
g1->zeros(); g1->zeros();
@ -154,21 +154,21 @@ KordpDynare::calcDerivatives(const Vector &yy, const Vector &xx) throw (DynareEx
if (nOrder > 1) if (nOrder > 1)
{ {
// generate g2 space for sparse Hessian 3x NNZH = 3x NNZD[1] // generate g2 space for sparse Hessian 3x NNZH = 3x NNZD[1]
g2 = new TwoDMatrix((int) (*NNZD)[1],3); g2 = new TwoDMatrix((int)(*NNZD)[1], 3);
g2->zeros(); g2->zeros();
} }
if (nOrder > 2) if (nOrder > 2)
{ {
g3 = new TwoDMatrix((int) (*NNZD)[2],3); g3 = new TwoDMatrix((int)(*NNZD)[2], 3);
g3->zeros(); g3->zeros();
} }
Vector out(nY); Vector out(nY);
out.zeros(); out.zeros();
const Vector *llxYYp; // getting around the constantness const Vector *llxYYp; // getting around the constantness
if ((nJcols - nExog) > yy.length()) if ((nJcols - nExog) > yy.length())
llxYYp = (LLxSteady(yy)); llxYYp = (LLxSteady(yy));
else else
llxYYp = &yy; llxYYp = &yy;
const Vector &llxYY = *(llxYYp); const Vector &llxYY = *(llxYYp);
dynamicDLL.eval(llxYY, xx, params, out, g1, g2, g3); dynamicDLL.eval(llxYY, xx, params, out, g1, g2, g3);
@ -180,13 +180,13 @@ KordpDynare::calcDerivatives(const Vector &yy, const Vector &xx) throw (DynareEx
delete g1; delete g1;
if (nOrder > 1) if (nOrder > 1)
{ {
populateDerivativesContainer(g2, 2, JacobianIndices); populateDerivativesContainer(g2, 2, JacobianIndices);
delete g2; delete g2;
} }
if (nOrder > 2) if (nOrder > 2)
{ {
populateDerivativesContainer(g3, 3, JacobianIndices); populateDerivativesContainer(g3, 3, JacobianIndices);
delete g3; delete g3;
} }
delete llxYYp; delete llxYYp;
} }
@ -200,8 +200,8 @@ KordpDynare::calcDerivativesAtSteady() throw (DynareException)
} }
/******************************************************************************* /*******************************************************************************
* populateDerivatives to sparse Tensor and fit it in the Derivatives Container * populateDerivatives to sparse Tensor and fit it in the Derivatives Container
*******************************************************************************/ *******************************************************************************/
void void
KordpDynare::populateDerivativesContainer(TwoDMatrix *g, int ord, const vector<int> *vOrder) KordpDynare::populateDerivativesContainer(TwoDMatrix *g, int ord, const vector<int> *vOrder)
{ {
@ -213,80 +213,80 @@ KordpDynare::populateDerivativesContainer(TwoDMatrix *g, int ord, const vector<i
if (ord == 1) if (ord == 1)
{ {
for (int i = 0; i < g->ncols(); i++) for (int i = 0; i < g->ncols(); i++)
{ {
for (int j = 0; j < g->nrows(); j++) for (int j = 0; j < g->nrows(); j++)
{ {
double x; double x;
if (s[0] < nJcols-nExog) if (s[0] < nJcols-nExog)
x = g->get(j, (*vOrder)[s[0]]); x = g->get(j, (*vOrder)[s[0]]);
else else
x = g->get(j, s[0]); x = g->get(j, s[0]);
if (x != 0.0) if (x != 0.0)
mdTi->insert(s, j, x); mdTi->insert(s, j, x);
} }
s[0]++; s[0]++;
} }
} }
else if (ord == 2) else if (ord == 2)
{ {
int nJcols1 = nJcols-nExog; int nJcols1 = nJcols-nExog;
vector<int> revOrder(nJcols1); vector<int> revOrder(nJcols1);
for (int i = 0; i < nJcols1; i++) for (int i = 0; i < nJcols1; i++)
revOrder[(*vOrder)[i]] = i; revOrder[(*vOrder)[i]] = i;
for (int i = 0; i < g->nrows(); i++) for (int i = 0; i < g->nrows(); i++)
{ {
int j = (int)g->get(i,0)-1; // hessian indices start with 1 int j = (int) g->get(i, 0)-1; // hessian indices start with 1
int i1 = (int)g->get(i,1) -1; int i1 = (int) g->get(i, 1) -1;
int s0 = (int)floor(((double) i1)/((double) nJcols)); int s0 = (int) floor(((double) i1)/((double) nJcols));
int s1 = i1- (nJcols*s0); int s1 = i1- (nJcols*s0);
if (s0 < nJcols1) if (s0 < nJcols1)
s[0] = revOrder[s0]; s[0] = revOrder[s0];
else else
s[0] = s0; s[0] = s0;
if (s1 < nJcols1) if (s1 < nJcols1)
s[1] = revOrder[s1]; s[1] = revOrder[s1];
else else
s[1] = s1; s[1] = s1;
if (s[1] >= s[0]) if (s[1] >= s[0])
{ {
double x = g->get(i,2); double x = g->get(i, 2);
mdTi->insert(s, j, x); mdTi->insert(s, j, x);
} }
} }
} }
else if (ord == 3) else if (ord == 3)
{ {
int nJcols1 = nJcols-nExog; int nJcols1 = nJcols-nExog;
int nJcols2 = nJcols*nJcols; int nJcols2 = nJcols*nJcols;
vector<int> revOrder(nJcols1); vector<int> revOrder(nJcols1);
for (int i = 0; i < nJcols1; i++) for (int i = 0; i < nJcols1; i++)
revOrder[(*vOrder)[i]] = i; revOrder[(*vOrder)[i]] = i;
for (int i = 0; i < g->nrows(); i++) for (int i = 0; i < g->nrows(); i++)
{ {
int j = (int)g->get(i,0)-1; int j = (int) g->get(i, 0)-1;
int i1 = (int)g->get(i,1) -1; int i1 = (int) g->get(i, 1) -1;
int s0 = (int)floor(((double) i1)/((double) nJcols2)); int s0 = (int) floor(((double) i1)/((double) nJcols2));
int i2 = i1 - nJcols2*s0; int i2 = i1 - nJcols2*s0;
int s1 = (int)floor(((double) i2)/((double) nJcols)); int s1 = (int) floor(((double) i2)/((double) nJcols));
int s2 = i2 - nJcols*s1; int s2 = i2 - nJcols*s1;
if (s0 < nJcols1) if (s0 < nJcols1)
s[0] = revOrder[s0]; s[0] = revOrder[s0];
else else
s[0] = s0; s[0] = s0;
if (s1 < nJcols1) if (s1 < nJcols1)
s[1] = revOrder[s1]; s[1] = revOrder[s1];
else else
s[1] = s1; s[1] = s1;
if (s2 < nJcols1) if (s2 < nJcols1)
s[2] = revOrder[s2]; s[2] = revOrder[s2];
else else
s[2] = s2; s[2] = s2;
if ((s[2] >= s[1]) && (s[1] >= s[0])) if ((s[2] >= s[1]) && (s[1] >= s[0]))
{ {
double x = g->get(i,2); double x = g->get(i, 2);
mdTi->insert(s, j, x); mdTi->insert(s, j, x);
} }
} }
} }
// md container // md container
@ -298,20 +298,20 @@ void
KordpDynare::writeModelInfo(Journal &jr) const KordpDynare::writeModelInfo(Journal &jr) const
{ {
// write info on variables // write info on variables
JournalRecordPair rp(journal); JournalRecordPair rp(journal);
rp << "Information on variables" << endrec; rp << "Information on variables" << endrec;
JournalRecord rec1(journal); JournalRecord rec1(journal);
rec1 << "Number of endogenous: " << ny() << endrec; rec1 << "Number of endogenous: " << ny() << endrec;
JournalRecord rec2(journal); JournalRecord rec2(journal);
rec2 << "Number of exogenous: " << nexog() << endrec; rec2 << "Number of exogenous: " << nexog() << endrec;
JournalRecord rec3(journal); JournalRecord rec3(journal);
rec3 << "Number of static: " << nstat() << endrec; rec3 << "Number of static: " << nstat() << endrec;
JournalRecord rec4(journal); JournalRecord rec4(journal);
rec4 << "Number of predetermined: " << npred()+nboth() << endrec; rec4 << "Number of predetermined: " << npred()+nboth() << endrec;
JournalRecord rec5(journal); JournalRecord rec5(journal);
rec5 << "Number of forward looking: " << nforw()+nboth() << endrec; rec5 << "Number of forward looking: " << nforw()+nboth() << endrec;
JournalRecord rec6(journal); JournalRecord rec6(journal);
rec6 << "Number of both: " << nboth() << endrec; rec6 << "Number of both: " << nboth() << endrec;
} }
/********************************************************* /*********************************************************
@ -328,35 +328,35 @@ KordpDynare::LLxSteady(const Vector &yS) throw (DynareException, TLException)
// create temporary square 2D matrix size nEndo x nEndo (sparse) // create temporary square 2D matrix size nEndo x nEndo (sparse)
// for the lag, current and lead blocks of the jacobian // for the lag, current and lead blocks of the jacobian
Vector *llxSteady = new Vector(nJcols-nExog); Vector *llxSteady = new Vector(nJcols-nExog);
for (int ll_row = 0; ll_row < ll_Incidence->nrows(); ll_row++) for (int ll_row = 0; ll_row < ll_Incidence->nrows(); ll_row++)
{
// populate (non-sparse) vector with ysteady values
for (int i = 0; i < nY; i++)
{ {
// populate (non-sparse) vector with ysteady values if (ll_Incidence->get(ll_row, i))
for (int i = 0; i < nY; i++) (*llxSteady)[((int) ll_Incidence->get(ll_row, i))-1] = yS[i];
{
if (ll_Incidence->get(ll_row, i))
(*llxSteady)[((int) ll_Incidence->get(ll_row, i))-1] = yS[i];
}
} }
}
return llxSteady; return llxSteady;
} }
/************************************ /************************************
* Reorder DynareJacobianIndices of variables in a vector according to * Reorder DynareJacobianIndices of variables in a vector according to
* given int * varOrder together with lead & lag incidence matrix and * given int * varOrder together with lead & lag incidence matrix and
* any the extra columns for exogenous vars, and then, * any the extra columns for exogenous vars, and then,
* reorders its blocks given by the varOrder and the Dynare++ expectations: * reorders its blocks given by the varOrder and the Dynare++ expectations:
* extra nboth+ npred (t-1) lags * extra nboth+ npred (t-1) lags
* varOrder * varOrder
static: static:
pred pred
both both
forward forward
* extra both + nforw (t+1) leads, and * extra both + nforw (t+1) leads, and
* extra exogen * extra exogen
* so to match the jacobian organisation expected by the Appoximation class * so to match the jacobian organisation expected by the Appoximation class
both + nforw (t+1) leads both + nforw (t+1) leads
static static
pred pred
@ -375,37 +375,37 @@ KordpDynare::ReorderDynareJacobianIndices(const vector<int> *varOrder) throw (TL
vector <int> tmp(nY); vector <int> tmp(nY);
int i, j, rjoff = nJcols-nExog-1; int i, j, rjoff = nJcols-nExog-1;
for (int ll_row = 0; ll_row < ll_Incidence->nrows(); ll_row++) for (int ll_row = 0; ll_row < ll_Incidence->nrows(); ll_row++)
{ {
// reorder in orde-var order & populate temporary nEndo (sparse) vector with // reorder in orde-var order & populate temporary nEndo (sparse) vector with
// the lag, current and lead blocks of the jacobian respectively // the lag, current and lead blocks of the jacobian respectively
for (i = 0; i < nY; i++) for (i = 0; i < nY; i++)
tmp[i] = ((int) ll_Incidence->get(ll_row, (*varOrder)[i]-1)); tmp[i] = ((int) ll_Incidence->get(ll_row, (*varOrder)[i]-1));
// write the reordered blocks back to the jacobian // write the reordered blocks back to the jacobian
// in reverse order // in reverse order
for (j = nY-1; j >= 0; j--) for (j = nY-1; j >= 0; j--)
if (tmp[j]) if (tmp[j])
{ {
(*JacobianIndices)[rjoff] = tmp[j] -1; (*JacobianIndices)[rjoff] = tmp[j] -1;
rjoff--; rjoff--;
if (rjoff < 0) if (rjoff < 0)
break; break;
} }
} }
//add the indices for the nExog exogenous jacobians //add the indices for the nExog exogenous jacobians
for (j = nJcols-nExog; j < nJcols; j++) for (j = nJcols-nExog; j < nJcols; j++)
(*JacobianIndices)[j] = j; (*JacobianIndices)[j] = j;
return JacobianIndices; return JacobianIndices;
} }
/************************************ /************************************
* Reorder first set of columns of variables in a (jacobian) matrix * Reorder first set of columns of variables in a (jacobian) matrix
* according to order given in varsOrder together with the extras * according to order given in varsOrder together with the extras
* assuming tdx ncols() - nExog is eaqual or less than length of varOrder and * assuming tdx ncols() - nExog is eaqual or less than length of varOrder and
* of any of its elements too. * of any of its elements too.
************************************/ ************************************/
void void
KordpDynare::ReorderCols(TwoDMatrix *tdx, const vector<int> *vOrder) throw (DynareException, TLException) KordpDynare::ReorderCols(TwoDMatrix *tdx, const vector<int> *vOrder) throw (DynareException, TLException)
@ -419,8 +419,8 @@ KordpDynare::ReorderCols(TwoDMatrix *tdx, const vector<int> *vOrder) throw (Dyna
tdx->zeros(); // empty original matrix tdx->zeros(); // empty original matrix
// reorder the columns // reorder the columns
for (int i = 0; i < tdx->ncols(); i++) for (int i = 0; i < tdx->ncols(); i++)
tdx->copyColumn(tmpR, (*vOrder)[i], i); tdx->copyColumn(tmpR, (*vOrder)[i], i);
} }
void void
@ -431,15 +431,15 @@ KordpDynare::ReorderCols(TwoDMatrix *tdx, const int *vOrder) throw (TLException)
TwoDMatrix &tmpR = tmp; TwoDMatrix &tmpR = tmp;
tdx->zeros(); // empty original matrix tdx->zeros(); // empty original matrix
// reorder the columns // reorder the columns
for (int i = 0; i < tdx->ncols(); i++) for (int i = 0; i < tdx->ncols(); i++)
tdx->copyColumn(tmpR, vOrder[i], i); tdx->copyColumn(tmpR, vOrder[i], i);
} }
/*********************************************************************** /***********************************************************************
* Recursive hierarchical block reordering of the higher order, input model * Recursive hierarchical block reordering of the higher order, input model
* derivatives inc. Hessian * derivatives inc. Hessian
* This is now obsolete but kept in in case it is needed * This is now obsolete but kept in in case it is needed
***********************************************************************/ ***********************************************************************/
void void
KordpDynare::ReorderBlocks(TwoDMatrix *tdx, const vector<int> *vOrder) throw (DynareException, TLException) KordpDynare::ReorderBlocks(TwoDMatrix *tdx, const vector<int> *vOrder) throw (DynareException, TLException)
@ -476,8 +476,8 @@ KordpDynare::ReorderBlocks(TwoDMatrix *tdx, const vector<int> *vOrder) throw (Dy
throw DynareException(__FILE__, __LINE__, "Size of order var is too small"); throw DynareException(__FILE__, __LINE__, "Size of order var is too small");
// reorder the columns // reorder the columns
for (int i = 0; i < tdx->ncols(); i++) for (int i = 0; i < tdx->ncols(); i++)
tdx->copyColumn(tmpR, (*vOrder)[i], i); tdx->copyColumn(tmpR, (*vOrder)[i], i);
} }
} }
@ -514,7 +514,7 @@ DynareNameList::selectIndices(const vector<const char *> &ns) const throw (Dynar
DynareNameList::DynareNameList(const KordpDynare &dynare) DynareNameList::DynareNameList(const KordpDynare &dynare)
{ {
for (int i = 0; i < dynare.ny(); i++) for (int i = 0; i < dynare.ny(); i++)
names.push_back(dynare.dnl->getName(i)); names.push_back(dynare.dnl->getName(i));
} }
DynareNameList::DynareNameList(const KordpDynare &dynare, const char **namesp) DynareNameList::DynareNameList(const KordpDynare &dynare, const char **namesp)
{ {
@ -538,7 +538,7 @@ DynareStateNameList::DynareStateNameList(const KordpDynare &dynare, const Dynare
const DynareExogNameList &denl) const DynareExogNameList &denl)
{ {
for (int i = 0; i < dynare.nys(); i++) for (int i = 0; i < dynare.nys(); i++)
names.push_back(dnl.getName(i+dynare.nstat())); names.push_back(dnl.getName(i+dynare.nstat()));
for (int i = 0; i < dynare.nexog(); i++) for (int i = 0; i < dynare.nexog(); i++)
names.push_back(denl.getName(i)); names.push_back(denl.getName(i));
} }

8
mex/sources/k_order_perturbation/k_ord_dynare.h
View File

@ -144,10 +144,10 @@ public:
KordpDynare(const char **endo, int num_endo, KordpDynare(const char **endo, int num_endo,
const char **exo, int num_exo, int num_par, const char **exo, int num_exo, int num_par,
Vector *ySteady, TwoDMatrix *vCov, Vector *params, int nstat, int nPred, Vector *ySteady, TwoDMatrix *vCov, Vector *params, int nstat, int nPred,
int nforw, int nboth, const int nJcols, const Vector *NNZD, int nforw, int nboth, const int nJcols, const Vector *NNZD,
const int nSteps, const int ord, const int nSteps, const int ord,
Journal &jr, DynamicModelDLL &dynamicDLL, double sstol, Journal &jr, DynamicModelDLL &dynamicDLL, double sstol,
const vector<int> *varOrder, const TwoDMatrix *ll_Incidence, const vector<int> *varOrder, const TwoDMatrix *ll_Incidence,
double qz_criterium) throw (TLException); double qz_criterium) throw (TLException);
/** Makes a deep copy of the object. */ /** Makes a deep copy of the object. */
@ -291,8 +291,8 @@ class KordpVectorFunction : public ogu::VectorFunction
protected: protected:
KordpDynare &d; KordpDynare &d;
public: public:
KordpVectorFunction(KordpDynare &dyn) KordpVectorFunction(KordpDynare &dyn) :
: d(dyn) d(dyn)
{ {
} }
virtual ~KordpVectorFunction() virtual ~KordpVectorFunction()

75
mex/sources/k_order_perturbation/k_order_perturbation.cpp
View File

@ -18,24 +18,24 @@
*/ */
/****************************************************** /******************************************************
// k_order_perturbation.cpp : Defines the entry point for the k-order perturbation application DLL. // k_order_perturbation.cpp : Defines the entry point for the k-order perturbation application DLL.
// //
// called from Dynare dr1_k_order.m, (itself called form resol.m instead of regular dr1.m) // called from Dynare dr1_k_order.m, (itself called form resol.m instead of regular dr1.m)
// if options_.order < 2 % 1st order only // if options_.order < 2 % 1st order only
// [ysteady, ghx_u]=k_order_perturbation(dr,task,M_,options_, oo_ , ['.' mexext]); // [ysteady, ghx_u]=k_order_perturbation(dr,task,M_,options_, oo_ , ['.' mexext]);
// else % 2nd order // else % 2nd order
// [ysteady, ghx_u, g_2]=k_order_perturbation(dr,task,M_,options_, oo_ , ['.' mexext]); // [ysteady, ghx_u, g_2]=k_order_perturbation(dr,task,M_,options_, oo_ , ['.' mexext]);
// inputs: // inputs:
// dr, - Dynare structure // dr, - Dynare structure
// task, - check or not, not used // task, - check or not, not used
// M_ - Dynare structure // M_ - Dynare structure
// options_ - Dynare structure // options_ - Dynare structure
// oo_ - Dynare structure // oo_ - Dynare structure
// ['.' mexext] Matlab dll extension // ['.' mexext] Matlab dll extension
// returns: // returns:
// ysteady steady state // ysteady steady state
// ghx_u - first order rules packed in one matrix // ghx_u - first order rules packed in one matrix
// g_2 - 2nd order rules packed in one matrix // g_2 - 2nd order rules packed in one matrix
**********************************************************/ **********************************************************/
#include "k_ord_dynare.h" #include "k_ord_dynare.h"
@ -84,7 +84,7 @@ DynareMxArrayToString(const char *cNamesCharStr, const int len, const int width)
{ {
char **cNamesMX; char **cNamesMX;
cNamesMX = (char **) calloc(len, sizeof(char *)); cNamesMX = (char **) calloc(len, sizeof(char *));
for(int i = 0; i < len; i++) for (int i = 0; i < len; i++)
cNamesMX[i] = (char *) calloc(width+1, sizeof(char)); cNamesMX[i] = (char *) calloc(width+1, sizeof(char));
for (int i = 0; i < width; i++) for (int i = 0; i < width; i++)
@ -99,7 +99,7 @@ DynareMxArrayToString(const char *cNamesCharStr, const int len, const int width)
else cNamesMX[j][i] = '\0'; else cNamesMX[j][i] = '\0';
} }
} }
return (const char **)cNamesMX; return (const char **) cNamesMX;
} }
////////////////////////////////////////////////////// //////////////////////////////////////////////////////
@ -112,7 +112,7 @@ DynareMxArrayToString(const mxArray *mxFldp, const int len, const int width)
{ {
char *cNamesCharStr = mxArrayToString(mxFldp); char *cNamesCharStr = mxArrayToString(mxFldp);
const char **ret = DynareMxArrayToString(cNamesCharStr, len, width); const char **ret = DynareMxArrayToString(cNamesCharStr, len, width);
return ret; return ret;
} }
@ -127,7 +127,7 @@ extern "C" {
mexErrMsgTxt("Must have at least 5 input parameters."); mexErrMsgTxt("Must have at least 5 input parameters.");
if (nlhs == 0) if (nlhs == 0)
mexErrMsgTxt("Must have at least 1 output parameter."); mexErrMsgTxt("Must have at least 1 output parameter.");
const mxArray *dr = prhs[0]; const mxArray *dr = prhs[0];
const int check_flag = (int) mxGetScalar(prhs[1]); const int check_flag = (int) mxGetScalar(prhs[1]);
const mxArray *M_ = prhs[2]; const mxArray *M_ = prhs[2];
@ -136,7 +136,7 @@ extern "C" {
mxArray *mFname = mxGetField(M_, 0, "fname"); mxArray *mFname = mxGetField(M_, 0, "fname");
if (!mxIsChar(mFname)) if (!mxIsChar(mFname))
mexErrMsgTxt("Input must be of type char."); mexErrMsgTxt("Input must be of type char.");
string fName = mxArrayToString(mFname); string fName = mxArrayToString(mFname);
const mxArray *mexExt = prhs[5]; const mxArray *mexExt = prhs[5];
string dfExt = mxArrayToString(mexExt); //Dynamic file extension, e.g.".dll" or .mexw32; string dfExt = mxArrayToString(mexExt); //Dynamic file extension, e.g.".dll" or .mexw32;
@ -147,12 +147,12 @@ extern "C" {
kOrder = (int) mxGetScalar(mxFldp); kOrder = (int) mxGetScalar(mxFldp);
else else
kOrder = 1; kOrder = 1;
if (kOrder == 1 && nlhs != 1) if (kOrder == 1 && nlhs != 1)
mexErrMsgTxt("k_order_perturbation at order 1 requires exactly 1 argument in output"); mexErrMsgTxt("k_order_perturbation at order 1 requires exactly 1 argument in output");
else if (kOrder > 1 && nlhs != kOrder+1) else if (kOrder > 1 && nlhs != kOrder+1)
mexErrMsgTxt("k_order_perturbation at order > 1 requires exactly order + 1 argument in output"); mexErrMsgTxt("k_order_perturbation at order > 1 requires exactly order + 1 argument in output");
double qz_criterium = 1+1e-6; double qz_criterium = 1+1e-6;
mxFldp = mxGetField(options_, 0, "qz_criterium"); mxFldp = mxGetField(options_, 0, "qz_criterium");
if (mxIsNumeric(mxFldp)) if (mxIsNumeric(mxFldp))
@ -201,11 +201,11 @@ extern "C" {
mxFldp = mxGetField(dr, 0, "order_var"); mxFldp = mxGetField(dr, 0, "order_var");
dparams = (double *) mxGetData(mxFldp); dparams = (double *) mxGetData(mxFldp);
npar = (int) mxGetM(mxFldp); npar = (int) mxGetM(mxFldp);
if (npar != nEndo) //(nPar != npar) if (npar != nEndo) //(nPar != npar)
mexErrMsgTxt("Incorrect number of input var_order vars."); mexErrMsgTxt("Incorrect number of input var_order vars.");
vector<int> *var_order_vp = (new vector<int>(nEndo)); vector<int> *var_order_vp = (new vector<int>(nEndo));
for (int v = 0; v < nEndo; v++) for (int v = 0; v < nEndo; v++)
(*var_order_vp)[v] = (int)(*(dparams++)); (*var_order_vp)[v] = (int)(*(dparams++));
// the lag, current and lead blocks of the jacobian respectively // the lag, current and lead blocks of the jacobian respectively
mxFldp = mxGetField(M_, 0, "lead_lag_incidence"); mxFldp = mxGetField(M_, 0, "lead_lag_incidence");
@ -213,15 +213,14 @@ extern "C" {
npar = (int) mxGetN(mxFldp); npar = (int) mxGetN(mxFldp);
int nrows = (int) mxGetM(mxFldp); int nrows = (int) mxGetM(mxFldp);
TwoDMatrix *llincidence = new TwoDMatrix(nrows, npar, dparams); TwoDMatrix *llincidence = new TwoDMatrix(nrows, npar, dparams);
if (npar != nEndo) if (npar != nEndo)
mexErrMsgIdAndTxt("dynare:k_order_perturbation", "Incorrect length of lead lag incidences: ncol=%d != nEndo=%d.", npar, nEndo); mexErrMsgIdAndTxt("dynare:k_order_perturbation", "Incorrect length of lead lag incidences: ncol=%d != nEndo=%d.", npar, nEndo);
//get NNZH =NNZD(2) = the total number of non-zero Hessian elements //get NNZH =NNZD(2) = the total number of non-zero Hessian elements
mxFldp = mxGetField(M_, 0, "NNZDerivatives"); mxFldp = mxGetField(M_, 0, "NNZDerivatives");
dparams = (double *) mxGetData(mxFldp); dparams = (double *) mxGetData(mxFldp);
Vector *NNZD = new Vector (dparams, (int) mxGetM(mxFldp)); Vector *NNZD = new Vector(dparams, (int) mxGetM(mxFldp));
const int jcols = nExog+nEndo+nsPred+nsForw; // Num of Jacobian columns const int jcols = nExog+nEndo+nsPred+nsForw; // Num of Jacobian columns
@ -235,7 +234,7 @@ extern "C" {
const int widthExog = (int) mxGetN(mxFldp); const int widthExog = (int) mxGetN(mxFldp);
const char **exoNamesMX = DynareMxArrayToString(mxFldp, nexo, widthExog); const char **exoNamesMX = DynareMxArrayToString(mxFldp, nexo, widthExog);
if ((nEndo != nendo) || (nExog != nexo)) //(nPar != npar) if ((nEndo != nendo) || (nExog != nexo)) //(nPar != npar)
mexErrMsgTxt("Incorrect number of input parameters."); mexErrMsgTxt("Incorrect number of input parameters.");
/* Fetch time index */ /* Fetch time index */
@ -259,11 +258,11 @@ extern "C" {
// make KordpDynare object // 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, ySteady, vCov, modParams, nStat, nPred, nForw, nBoth,
jcols, NNZD, nSteps, kOrder, journal, dynamicDLL, jcols, NNZD, nSteps, kOrder, journal, dynamicDLL,
sstol, var_order_vp, llincidence, qz_criterium); sstol, var_order_vp, llincidence, qz_criterium);
// construct main K-order approximation class // construct main K-order approximation class
Approximation app(dynare, journal, nSteps, false, qz_criterium); Approximation app(dynare, journal, nSteps, false, qz_criterium);
// run stochastic steady // run stochastic steady
app.walkStochSteady(); app.walkStochSteady();
@ -301,11 +300,11 @@ extern "C" {
if (kOrder == 1) if (kOrder == 1)
{ {
/* Set the output pointer to the output matrix ysteady. */ /* Set the output pointer to the output matrix ysteady. */
map<string, ConstTwoDMatrix>::const_iterator cit = mm.begin(); map<string, ConstTwoDMatrix>::const_iterator cit = mm.begin();
++cit; ++cit;
plhs[0] = mxCreateDoubleMatrix((*cit).second.numRows(), (*cit).second.numCols(), mxREAL); plhs[0] = mxCreateDoubleMatrix((*cit).second.numRows(), (*cit).second.numCols(), mxREAL);
TwoDMatrix g((*cit).second.numRows(), (*cit).second.numCols(), mxGetPr(plhs[0])); TwoDMatrix g((*cit).second.numRows(), (*cit).second.numCols(), mxGetPr(plhs[0]));
g = (const TwoDMatrix &)(*cit).second; g = (const TwoDMatrix &)(*cit).second;
} }
if (kOrder >= 2) if (kOrder >= 2)
{ {

258
mex/sources/k_order_perturbation/tests/first_order.m
View File

@ -1,129 +1,129 @@
function [gy]=first_order(M_, dr, jacobia) function [gy]=first_order(M_, dr, jacobia)
% Emulation of Dynare++ c++ first_order.cpp for testing pruposes % Emulation of Dynare++ c++ first_order.cpp for testing pruposes
% Copyright (C) 2009 Dynare Team % Copyright (C) 2009 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
% Dynare is free software: you can redistribute it and/or modify % Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by % it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or % the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version. % (at your option) any later version.
% %
% Dynare is distributed in the hope that it will be useful, % Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of % but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details. % GNU General Public License for more details.
% %
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>. % along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% fd = jacobia_ % fd = jacobia_
% reorder jacobia_ % reorder jacobia_
[fd]=k_reOrderedJacobia(M_, jacobia) [fd]=k_reOrderedJacobia(M_, jacobia)
%ypart=dr; %ypart=dr;
ypart.ny=M_.endo_nbr; ypart.ny=M_.endo_nbr;
ypart.nyss=dr.nfwrd+dr.nboth; ypart.nyss=dr.nfwrd+dr.nboth;
ypart.nys=dr.npred; ypart.nys=dr.npred;
ypart.npred=dr.npred-dr.nboth; ypart.npred=dr.npred-dr.nboth;
ypart.nboth=dr.nboth; ypart.nboth=dr.nboth;
ypart.nforw=dr.nfwrd; ypart.nforw=dr.nfwrd;
ypart.nstat =dr.nstatic ypart.nstat =dr.nstatic
nu=M_.exo_nbr nu=M_.exo_nbr
off= 1; % = 0 in C off= 1; % = 0 in C
fyplus = fd(:,off:off+ypart.nyss-1); fyplus = fd(:,off:off+ypart.nyss-1);
off= off+ypart.nyss; off= off+ypart.nyss;
fyszero= fd(:,off:off+ypart.nstat-1); fyszero= fd(:,off:off+ypart.nstat-1);
off= off+ypart.nstat; off= off+ypart.nstat;
fypzero= fd(:,off:off+ypart.npred-1); fypzero= fd(:,off:off+ypart.npred-1);
off= off+ypart.npred; off= off+ypart.npred;
fybzero= fd(:,off:off+ypart.nboth-1); fybzero= fd(:,off:off+ypart.nboth-1);
off= off+ypart.nboth; off= off+ypart.nboth;
fyfzero= fd(:,off:off+ypart.nforw-1); fyfzero= fd(:,off:off+ypart.nforw-1);
off= off+ypart.nforw; off= off+ypart.nforw;
fymins= fd(:,off:off+ypart.nys-1); fymins= fd(:,off:off+ypart.nys-1);
off= off+ypart.nys; off= off+ypart.nys;
fuzero= fd(:,off:off+nu-1); fuzero= fd(:,off:off+nu-1);
off=off+ nu; off=off+ nu;
n= ypart.ny+ypart.nboth; n= ypart.ny+ypart.nboth;
%TwoDMatrix %TwoDMatrix
matD=zeros(n,n); matD=zeros(n,n);
% matD.place(fypzero,0,0); % matD.place(fypzero,0,0);
matD(1:n-ypart.nboth,1:ypart.npred)= fypzero; matD(1:n-ypart.nboth,1:ypart.npred)= fypzero;
% matD.place(fybzero,0,ypart.npred); % matD.place(fybzero,0,ypart.npred);
matD(1:n-ypart.nboth,ypart.npred+1:ypart.npred+ypart.nboth)=fybzero; matD(1:n-ypart.nboth,ypart.npred+1:ypart.npred+ypart.nboth)=fybzero;
% matD.place(fyplus,0,ypart.nys()+ypart.nstat); % matD.place(fyplus,0,ypart.nys()+ypart.nstat);
matD(1:n-ypart.nboth,ypart.nys+ypart.nstat+1:ypart.nys+ypart.nstat+ypart.nyss)=fyplus; matD(1:n-ypart.nboth,ypart.nys+ypart.nstat+1:ypart.nys+ypart.nstat+ypart.nyss)=fyplus;
for i=1:ypart.nboth for i=1:ypart.nboth
matD(ypart.ny()+i,ypart.npred+i)= 1.0; matD(ypart.ny()+i,ypart.npred+i)= 1.0;
end end
matE=[fymins, fyszero, zeros(n-ypart.nboth,ypart.nboth), fyfzero; zeros(ypart.nboth,n)]; matE=[fymins, fyszero, zeros(n-ypart.nboth,ypart.nboth), fyfzero; zeros(ypart.nboth,n)];
% matE.place(fymins; % matE.place(fymins;
% matE.place(fyszero,0,ypart.nys()); % matE.place(fyszero,0,ypart.nys());
% matE.place(fyfzero,0,ypart.nys()+ypart.nstat+ypart.nboth); % matE.place(fyfzero,0,ypart.nys()+ypart.nstat+ypart.nboth);
for i= 1:ypart.nboth for i= 1:ypart.nboth
matE(ypart.ny()+i,ypart.nys()+ypart.nstat+i)= -1.0; matE(ypart.ny()+i,ypart.nys()+ypart.nstat+i)= -1.0;
end end
matE=-matE; %matE.mult(-1.0); matE=-matE; %matE.mult(-1.0);
% vsl=zeros(n,n); % vsl=zeros(n,n);
% vsr=zeros(n,n); % vsr=zeros(n,n);
% lwork= 100*n+16; % lwork= 100*n+16;
% work=zeros(1,lwork); % work=zeros(1,lwork);
% bwork=zeros(1,n); % bwork=zeros(1,n);
%int info; %int info;
% LAPACK_dgges("N","V","S",order_eigs,&n,matE.getData().base(),&n, % LAPACK_dgges("N","V","S",order_eigs,&n,matE.getData().base(),&n,
% matD.getData().base(),&n,&sdim,alphar.base(),alphai.base(), % matD.getData().base(),&n,&sdim,alphar.base(),alphai.base(),
% beta.base(),vsl.getData().base(),&n,vsr.getData().base(),&n, % beta.base(),vsl.getData().base(),&n,vsr.getData().base(),&n,
% work.base(),&lwork,&(bwork[0]),&info); % work.base(),&lwork,&(bwork[0]),&info);
[matE1,matD1,vsr,sdim,dr.eigval,info] = mjdgges(matE,matD,1); [matE1,matD1,vsr,sdim,dr.eigval,info] = mjdgges(matE,matD,1);
bk_cond= (sdim==ypart.nys); bk_cond= (sdim==ypart.nys);
% ConstGeneralMatrix z11(vsr,0,0,ypart.nys(),ypart.nys()); % ConstGeneralMatrix z11(vsr,0,0,ypart.nys(),ypart.nys());
z11=vsr(1:ypart.nys,1:ypart.nys); z11=vsr(1:ypart.nys,1:ypart.nys);
% ConstGeneralMatrix z12(vsr,0,ypart.nys(),ypart.nys(),n-ypart.nys()); % ConstGeneralMatrix z12(vsr,0,ypart.nys(),ypart.nys(),n-ypart.nys());
z12=vsr(1:ypart.nys(),ypart.nys+1:end);%, n-ypart.nys); z12=vsr(1:ypart.nys(),ypart.nys+1:end);%, n-ypart.nys);
% ConstGeneralMatrix z21(vsr,ypart.nys(),0,n-ypart.nys(),ypart.nys()); % ConstGeneralMatrix z21(vsr,ypart.nys(),0,n-ypart.nys(),ypart.nys());
z21=vsr(ypart.nys+1:end,1:ypart.nys); z21=vsr(ypart.nys+1:end,1:ypart.nys);
% ConstGeneralMatrix z22(vsr,ypart.nys(),ypart.nys(),n-ypart.nys(),n-ypart.nys()); % ConstGeneralMatrix z22(vsr,ypart.nys(),ypart.nys(),n-ypart.nys(),n-ypart.nys());
z22=vsr(ypart.nys+1:end,ypart.nys+1:end); z22=vsr(ypart.nys+1:end,ypart.nys+1:end);
% GeneralMatrix sfder(z12,"transpose"); % GeneralMatrix sfder(z12,"transpose");
sfder=z12';%,"transpose"); sfder=z12';%,"transpose");
% z22.multInvLeftTrans(sfder); % z22.multInvLeftTrans(sfder);
sfder=z22'\sfder; sfder=z22'\sfder;
sfder=-sfder;% .mult(-1); sfder=-sfder;% .mult(-1);
%s11(matE,0,0,ypart.nys(),ypart.nys()); %s11(matE,0,0,ypart.nys(),ypart.nys());
s11=matE1(1:ypart.nys,1:ypart.nys); s11=matE1(1:ypart.nys,1:ypart.nys);
% t11=(matD1,0,0,ypart.nys(),ypart.nys()); % t11=(matD1,0,0,ypart.nys(),ypart.nys());
t11=matD1(1:ypart.nys,1:ypart.nys); t11=matD1(1:ypart.nys,1:ypart.nys);
dumm=(s11');%,"transpose"); dumm=(s11');%,"transpose");
%z11.multInvLeftTrans(dumm); %z11.multInvLeftTrans(dumm);
dumm=z11'\dumm; dumm=z11'\dumm;
preder=(dumm');%,"transpose"); preder=(dumm');%,"transpose");
%t11.multInvLeft(preder); %t11.multInvLeft(preder);
preder=t11\preder; preder=t11\preder;
%preder.multLeft(z11); %preder.multLeft(z11);
preder= z11*preder; preder= z11*preder;
% gy.place(preder,ypart.nstat,0); % gy.place(preder,ypart.nstat,0);
% gy=(zeros(ypart.nstat,size(preder,2)) ;preder); % gy=(zeros(ypart.nstat,size(preder,2)) ;preder);
% sder(sfder,0,0,ypart.nstat,ypart.nys()); % sder(sfder,0,0,ypart.nstat,ypart.nys());
sder=sfder(1:ypart.nstat,1:ypart.nys); sder=sfder(1:ypart.nstat,1:ypart.nys);
% gy.place(sder,0,0); % gy.place(sder,0,0);
% gy(1:ypart.nstat, 1:ypart.nys)=sder; % gy(1:ypart.nstat, 1:ypart.nys)=sder;
% gy=[sder;preder]; % gy=[sder;preder];
% fder(sfder,ypart.nstat+ypart.nboth,0,ypart.nforw,ypart.nys()); % fder(sfder,ypart.nstat+ypart.nboth,0,ypart.nforw,ypart.nys());
fder=sfder(ypart.nstat+ypart.nboth+1:ypart.nstat+ypart.nboth+ypart.nforw,1:ypart.nys); fder=sfder(ypart.nstat+ypart.nboth+1:ypart.nstat+ypart.nboth+ypart.nforw,1:ypart.nys);
% gy.place(fder,ypart.nstat+ypart.nys(),0); % gy.place(fder,ypart.nstat+ypart.nys(),0);
% gy(ypart.nstat+ypart.nys,1)=fder; % gy(ypart.nstat+ypart.nys,1)=fder;
gy=[sder;preder;fder]; gy=[sder;preder;fder];

6
mex/sources/k_order_perturbation/tests/k_order_test_main.cpp
View File

@ -81,7 +81,7 @@ main(int argc, char *argv[])
const int nSteady = 16; //27 //31;//29, 16 (int)mxGetM(mxFldp); 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}; 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" ); //mxFldp = mxGetField(dr, 0,"nstatic" );
@ -111,7 +111,7 @@ main(int argc, char *argv[])
= { = {
5, 6, 8, 10, 11, 12, 14, 7, 13, 1, 2, 3, 4, 9, 15, 16 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, 16, 7, 13, 14, 15, 1, 2, 3, 4, 9, 17, 18
}; };
//Vector * varOrder = new Vector(var_order, nEndo); //Vector * varOrder = new Vector(var_order, nEndo);
vector<int> *var_order_vp = new vector<int>(nEndo); //nEndo)); vector<int> *var_order_vp = new vector<int>(nEndo); //nEndo));
for (int v = 0; v < nEndo; v++) for (int v = 0; v < nEndo; v++)
@ -135,7 +135,7 @@ main(int argc, char *argv[])
0, 18, 0, 0, 18, 0,
0, 19, 26, 0, 19, 26,
0, 20, 27 0, 20, 27
}; };
TwoDMatrix *llincidence = new TwoDMatrix(3, nEndo, ll_incidence); TwoDMatrix *llincidence = new TwoDMatrix(3, nEndo, ll_incidence);
const int jcols = nExog+nEndo+nsPred+nsForw; // Num of Jacobian columns const int jcols = nExog+nEndo+nsPred+nsForw; // Num of Jacobian columns

198
mex/sources/kronecker/A_times_B_kronecker_C.cc
View File

@ -18,8 +18,8 @@
*/ */
/* /*
* This mex file computes A*kron(B,C) or A*kron(B,B) without explicitely building kron(B,C) or kron(B,B), so that * This mex file computes A*kron(B,C) or A*kron(B,B) without explicitely building kron(B,C) or kron(B,B), so that
* one can consider large matrices B and/or C. * one can consider large matrices B and/or C.
*/ */
#include <string.h> #include <string.h>
@ -28,105 +28,107 @@
#include <dynblas.h> #include <dynblas.h>
#ifdef USE_OMP #ifdef USE_OMP
#include <omp.h> # include <omp.h>
#endif #endif
void full_A_times_kronecker_B_C(double *A, double *B, double *C, double *D, void
blas_int mA, blas_int nA, blas_int mB, blas_int nB, blas_int mC, blas_int nC) full_A_times_kronecker_B_C(double *A, double *B, double *C, double *D,
blas_int mA, blas_int nA, blas_int mB, blas_int nB, blas_int mC, blas_int nC)
{ {
#if USE_OMP #if USE_OMP
#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS"))) # pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
for(long int colD = 0; colD<nB*nC; colD++) for (long int colD = 0; colD < nB*nC; colD++)
{ {
#if DEBUG_OMP # if DEBUG_OMP
mexPrintf("%d thread number is %d (%d).\n",colD,omp_get_thread_num(),omp_get_num_threads()); mexPrintf("%d thread number is %d (%d).\n", colD, omp_get_thread_num(), omp_get_num_threads());
#endif # endif
unsigned int colB = colD/nC; unsigned int colB = colD/nC;
unsigned int colC = colD%nC; unsigned int colC = colD%nC;
for(int colA = 0; colA<nA; colA++) for (int colA = 0; colA < nA; colA++)
{ {
unsigned int rowB = colA/mC; unsigned int rowB = colA/mC;
unsigned int rowC = colA%mC; unsigned int rowC = colA%mC;
unsigned int idxA = colA*mA; unsigned int idxA = colA*mA;
unsigned int idxD = colD*mA; unsigned int idxD = colD*mA;
double BC = B[colB*mB+rowB]*C[colC*mC+rowC]; double BC = B[colB*mB+rowB]*C[colC*mC+rowC];
for(int rowD = 0; rowD<mA; rowD++) for (int rowD = 0; rowD < mA; rowD++)
{ {
D[idxD+rowD] += A[idxA+rowD]*BC; D[idxD+rowD] += A[idxA+rowD]*BC;
} }
} }
} }
#else #else
const unsigned long shiftA = mA*mC ; const unsigned long shiftA = mA*mC;
const unsigned long shiftD = mA*nC ; const unsigned long shiftD = mA*nC;
unsigned long int kd = 0, ka = 0 ; unsigned long int kd = 0, ka = 0;
char transpose[2] = "N"; char transpose[2] = "N";
double one = 1.0 ; double one = 1.0;
for(unsigned long int col=0; col<nB; col++) for (unsigned long int col = 0; col < nB; col++)
{ {
ka = 0 ; ka = 0;
for(unsigned long int row=0; row<mB; row++) for (unsigned long int row = 0; row < mB; row++)
{ {
dgemm(transpose, transpose, &mA, &nC, &mC, &B[mB*col+row], &A[ka], &mA, &C[0], &mC, &one, &D[kd], &mA); dgemm(transpose, transpose, &mA, &nC, &mC, &B[mB*col+row], &A[ka], &mA, &C[0], &mC, &one, &D[kd], &mA);
ka += shiftA; ka += shiftA;
} }
kd += shiftD; kd += shiftD;
} }
#endif #endif
} }
void
void full_A_times_kronecker_B_B(double *A, double *B, double *D, blas_int mA, blas_int nA, blas_int mB, blas_int nB) full_A_times_kronecker_B_B(double *A, double *B, double *D, blas_int mA, blas_int nA, blas_int mB, blas_int nB)
{ {
#if USE_OMP #if USE_OMP
#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS"))) # pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
for(long int colD = 0; colD<nB*nB; colD++) for (long int colD = 0; colD < nB*nB; colD++)
{ {
#if DEBUG_OMP # if DEBUG_OMP
mexPrintf("%d thread number is %d (%d).\n",colD,omp_get_thread_num(),omp_get_num_threads()); mexPrintf("%d thread number is %d (%d).\n", colD, omp_get_thread_num(), omp_get_num_threads());
#endif # endif
unsigned int j1B = colD/nB; unsigned int j1B = colD/nB;
unsigned int j2B = colD%nB; unsigned int j2B = colD%nB;
for(int colA = 0; colA<nA; colA++) for (int colA = 0; colA < nA; colA++)
{ {
unsigned int i1B = colA/mB; unsigned int i1B = colA/mB;
unsigned int i2B = colA%mB; unsigned int i2B = colA%mB;
unsigned int idxA = colA*mA; unsigned int idxA = colA*mA;
unsigned int idxD = colD*mA; unsigned int idxD = colD*mA;
double BB = B[j1B*mB+i1B]*B[j2B*mB+i2B]; double BB = B[j1B*mB+i1B]*B[j2B*mB+i2B];
for(int rowD = 0; rowD<mA; rowD++) for (int rowD = 0; rowD < mA; rowD++)
{ {
D[idxD+rowD] += A[idxA+rowD]*BB; D[idxD+rowD] += A[idxA+rowD]*BB;
} }
} }
} }
#else #else
const unsigned long int shiftA = mA*mB ; const unsigned long int shiftA = mA*mB;
const unsigned long int shiftD = mA*nB ; const unsigned long int shiftD = mA*nB;
unsigned long int kd = 0, ka = 0 ; unsigned long int kd = 0, ka = 0;
char transpose[2] = "N"; char transpose[2] = "N";
double one = 1.0; double one = 1.0;
for(unsigned long int col=0; col<nB; col++) for (unsigned long int col = 0; col < nB; col++)
{ {
ka = 0; ka = 0;
for(unsigned long int row=0; row<mB; row++) for (unsigned long int row = 0; row < mB; row++)
{ {
dgemm(transpose, transpose, &mA, &nB, &mB, &B[mB*col+row], &A[ka], &mA, &B[0], &mB, &one, &D[kd], &mA); dgemm(transpose, transpose, &mA, &nB, &mB, &B[mB*col+row], &A[ka], &mA, &B[0], &mB, &one, &D[kd], &mA);
ka += shiftA; ka += shiftA;
} }
kd += shiftD; kd += shiftD;
} }
#endif #endif
} }
void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ) void
mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{ {
// Check input and output: // Check input and output:
if ( (nrhs > 3) || (nrhs <2) ) if ((nrhs > 3) || (nrhs < 2))
{ {
mexErrMsgTxt("Two or Three input arguments required."); mexErrMsgTxt("Two or Three input arguments required.");
} }
if (nlhs>1) if (nlhs > 1)
{ {
mexErrMsgTxt("Too many output arguments."); mexErrMsgTxt("Too many output arguments.");
} }
@ -136,21 +138,21 @@ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] )
nA = mxGetN(prhs[0]); nA = mxGetN(prhs[0]);
mB = mxGetM(prhs[1]); mB = mxGetM(prhs[1]);
nB = mxGetN(prhs[1]); nB = mxGetN(prhs[1]);
if (nrhs == 3)// A*kron(B,C) is to be computed. if (nrhs == 3) // A*kron(B,C) is to be computed.
{ {
mC = mxGetM(prhs[2]); mC = mxGetM(prhs[2]);
nC = mxGetN(prhs[2]); nC = mxGetN(prhs[2]);
if (mB*mC != nA) if (mB*mC != nA)
{ {
mexErrMsgTxt("Input dimension error!"); mexErrMsgTxt("Input dimension error!");
} }
} }
else// A*kron(B,B) is to be computed. else // A*kron(B,B) is to be computed.
{ {
if (mB*mB != nA) if (mB*mB != nA)
{ {
mexErrMsgTxt("Input dimension error!"); mexErrMsgTxt("Input dimension error!");
} }
} }
// Get input matrices: // Get input matrices:
double *B, *C, *A; double *B, *C, *A;
@ -164,11 +166,11 @@ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] )
double *D; double *D;
if (nrhs == 3) if (nrhs == 3)
{ {
plhs[0] = mxCreateDoubleMatrix(mA,nB*nC,mxREAL); plhs[0] = mxCreateDoubleMatrix(mA, nB*nC, mxREAL);
} }
else else
{ {
plhs[0] = mxCreateDoubleMatrix(mA,nB*nB,mxREAL); plhs[0] = mxCreateDoubleMatrix(mA, nB*nB, mxREAL);
} }
D = mxGetPr(plhs[0]); D = mxGetPr(plhs[0]);
// Computational part: // Computational part:

191
mex/sources/kronecker/sparse_hessian_times_B_kronecker_C.cc
View File

@ -19,7 +19,7 @@
/* /*
* This mex file computes A*kron(B,C) or A*kron(B,B) without explicitly building kron(B,C) or kron(B,B), so that * This mex file computes A*kron(B,C) or A*kron(B,B) without explicitly building kron(B,C) or kron(B,B), so that
* one can consider large matrices A, B and/or C, and assuming that A is a the hessian of a dsge model * one can consider large matrices A, B and/or C, and assuming that A is a the hessian of a dsge model
* (dynare format). This mex file should not be used outside dr1.m. * (dynare format). This mex file should not be used outside dr1.m.
*/ */
@ -31,77 +31,79 @@
# include <omp.h> # include <omp.h>
#endif #endif
void sparse_hessian_times_B_kronecker_B(mwIndex *isparseA, mwIndex *jsparseA, double *vsparseA, void
double *B, double *D, mwSize mA, mwSize nA, mwSize mB, mwSize nB) sparse_hessian_times_B_kronecker_B(mwIndex *isparseA, mwIndex *jsparseA, double *vsparseA,
double *B, double *D, mwSize mA, mwSize nA, mwSize mB, mwSize nB)
{ {
/* /*
** Loop over the columns of kron(B,B) (or of the result matrix D). ** Loop over the columns of kron(B,B) (or of the result matrix D).
** This loop is splitted into two nested loops because we use the ** This loop is splitted into two nested loops because we use the
** symmetric pattern of the hessian matrix. ** symmetric pattern of the hessian matrix.
*/ */
#if USE_OMP #if USE_OMP
#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS"))) # pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
#endif #endif
for(int j1B=0; j1B<nB; j1B++) for (int j1B = 0; j1B < nB; j1B++)
{ {
#if DEBUG_OMP #if DEBUG_OMP
mexPrintf("%d thread number is %d (%d).\n",j1B,omp_get_thread_num(),omp_get_num_threads()); mexPrintf("%d thread number is %d (%d).\n", j1B, omp_get_thread_num(), omp_get_num_threads());
#endif #endif
for(unsigned int j2B=j1B; j2B<nB; j2B++) for (unsigned int j2B = j1B; j2B < nB; j2B++)
{ {
unsigned long int jj = j1B*nB+j2B;// column of kron(B,B) index. unsigned long int jj = j1B*nB+j2B; // column of kron(B,B) index.
unsigned long int iv =0; unsigned long int iv = 0;
unsigned int nz_in_column_ii_of_A = 0; unsigned int nz_in_column_ii_of_A = 0;
unsigned int k1 = 0; unsigned int k1 = 0;
unsigned int k2 = 0; unsigned int k2 = 0;
/* /*
** Loop over the rows of kron(B,B) (column jj). ** Loop over the rows of kron(B,B) (column jj).
*/ */
for(unsigned long int ii=0; ii<nA; ii++) for (unsigned long int ii = 0; ii < nA; ii++)
{ {
k1 = jsparseA[ii]; k1 = jsparseA[ii];
k2 = jsparseA[ii+1]; k2 = jsparseA[ii+1];
if (k1 < k2)// otherwise column ii of A does not have non zero elements (and there is nothing to compute). if (k1 < k2) // otherwise column ii of A does not have non zero elements (and there is nothing to compute).
{ {
++nz_in_column_ii_of_A; ++nz_in_column_ii_of_A;
unsigned int i1B = (ii/mB); unsigned int i1B = (ii/mB);
unsigned int i2B = (ii%mB); unsigned int i2B = (ii%mB);
double bb = B[j1B*mB+i1B]*B[j2B*mB+i2B]; double bb = B[j1B*mB+i1B]*B[j2B*mB+i2B];
/* /*
** Loop over the non zero entries of A(:,ii). ** Loop over the non zero entries of A(:,ii).
*/ */
for(unsigned int k=k1; k<k2; k++) for (unsigned int k = k1; k < k2; k++)
{ {
unsigned int kk = isparseA[k]; unsigned int kk = isparseA[k];
D[jj*mA+kk] = D[jj*mA+kk] + bb*vsparseA[iv]; D[jj*mA+kk] = D[jj*mA+kk] + bb*vsparseA[iv];
iv++; iv++;
} }
} }
} }
if (nz_in_column_ii_of_A>0) if (nz_in_column_ii_of_A > 0)
{ {
memcpy(&D[(j2B*nB+j1B)*mA],&D[jj*mA],mA*sizeof(double)); memcpy(&D[(j2B*nB+j1B)*mA], &D[jj*mA], mA*sizeof(double));
} }
} }
} }
} }
void sparse_hessian_times_B_kronecker_C(mwIndex *isparseA, mwIndex *jsparseA, double *vsparseA, void
double *B, double *C, double *D, sparse_hessian_times_B_kronecker_C(mwIndex *isparseA, mwIndex *jsparseA, double *vsparseA,
mwSize mA, mwSize nA, mwSize mB, mwSize nB, mwSize mC, mwSize nC) double *B, double *C, double *D,
mwSize mA, mwSize nA, mwSize mB, mwSize nB, mwSize mC, mwSize nC)
{ {
/* /*
** Loop over the columns of kron(B,B) (or of the result matrix D). ** Loop over the columns of kron(B,B) (or of the result matrix D).
*/ */
#if USE_OMP #if USE_OMP
#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS"))) # pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
#endif #endif
for(long int jj=0; jj<nB*nC; jj++)// column of kron(B,C) index. for (long int jj = 0; jj < nB*nC; jj++) // column of kron(B,C) index.
{ {
// Uncomment the following line to check if all processors are used. // Uncomment the following line to check if all processors are used.
#if DEBUG_OMP #if DEBUG_OMP
mexPrintf("%d thread number is %d (%d).\n",jj,omp_get_thread_num(),omp_get_num_threads()); mexPrintf("%d thread number is %d (%d).\n", jj, omp_get_thread_num(), omp_get_num_threads());
#endif #endif
unsigned int jB = jj/nC; unsigned int jB = jj/nC;
unsigned int jC = jj%nC; unsigned int jC = jj%nC;
unsigned int k1 = 0; unsigned int k1 = 0;
@ -111,38 +113,39 @@ void sparse_hessian_times_B_kronecker_C(mwIndex *isparseA, mwIndex *jsparseA, do
/* /*
** Loop over the rows of kron(B,C) (column jj). ** Loop over the rows of kron(B,C) (column jj).
*/ */
for(unsigned long int ii=0; ii<nA; ii++) for (unsigned long int ii = 0; ii < nA; ii++)
{ {
k1 = jsparseA[ii]; k1 = jsparseA[ii];
k2 = jsparseA[ii+1]; k2 = jsparseA[ii+1];
if (k1 < k2)// otherwise column ii of A does not have non zero elements (and there is nothing to compute). if (k1 < k2) // otherwise column ii of A does not have non zero elements (and there is nothing to compute).
{ {
++nz_in_column_ii_of_A; ++nz_in_column_ii_of_A;
unsigned int iC = (ii%mB); unsigned int iC = (ii%mB);
unsigned int iB = (ii/mB); unsigned int iB = (ii/mB);
double cb = C[jC*mC+iC]*B[jB*mB+iB]; double cb = C[jC*mC+iC]*B[jB*mB+iB];
/* /*
** Loop over the non zero entries of A(:,ii). ** Loop over the non zero entries of A(:,ii).
*/ */
for(unsigned int k=k1; k<k2; k++) for (unsigned int k = k1; k < k2; k++)
{ {
unsigned int kk = isparseA[k]; unsigned int kk = isparseA[k];
D[jj*mA+kk] = D[jj*mA+kk] + cb*vsparseA[iv]; D[jj*mA+kk] = D[jj*mA+kk] + cb*vsparseA[iv];
iv++; iv++;
} }
} }
} }
} }
} }
void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ) void
mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{ {
// Check input and output: // Check input and output:
if ( (nrhs > 3) || (nrhs <2) ) if ((nrhs > 3) || (nrhs < 2))
{ {
mexErrMsgTxt("Two or Three input arguments required."); mexErrMsgTxt("Two or Three input arguments required.");
} }
if (nlhs>1) if (nlhs > 1)
{ {
mexErrMsgTxt("Too many output arguments."); mexErrMsgTxt("Too many output arguments.");
} }
@ -156,21 +159,21 @@ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] )
nA = mxGetN(prhs[0]); nA = mxGetN(prhs[0]);
mB = mxGetM(prhs[1]); mB = mxGetM(prhs[1]);
nB = mxGetN(prhs[1]); nB = mxGetN(prhs[1]);
if (nrhs == 3)// A*kron(B,C) is to be computed. if (nrhs == 3) // A*kron(B,C) is to be computed.
{ {
mC = mxGetM(prhs[2]); mC = mxGetM(prhs[2]);
nC = mxGetN(prhs[2]); nC = mxGetN(prhs[2]);
if (mB*mC != nA) if (mB*mC != nA)
{ {
mexErrMsgTxt("Input dimension error!"); mexErrMsgTxt("Input dimension error!");
} }
} }
else// A*kron(B,B) is to be computed. else // A*kron(B,B) is to be computed.
{ {
if (mB*mB != nA) if (mB*mB != nA)
{ {
mexErrMsgTxt("Input dimension error!"); mexErrMsgTxt("Input dimension error!");
} }
} }
// Get input matrices: // Get input matrices:
double *B, *C; double *B, *C;
@ -180,23 +183,23 @@ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] )
C = mxGetPr(prhs[2]); C = mxGetPr(prhs[2]);
} }
// Sparse (dynare) hessian matrix. // Sparse (dynare) hessian matrix.
mwIndex *isparseA = (mwIndex*)mxGetIr(prhs[0]); mwIndex *isparseA = (mwIndex *) mxGetIr(prhs[0]);
mwIndex *jsparseA = (mwIndex*)mxGetJc(prhs[0]); mwIndex *jsparseA = (mwIndex *) mxGetJc(prhs[0]);
double *vsparseA = mxGetPr(prhs[0]); double *vsparseA = mxGetPr(prhs[0]);
// Initialization of the ouput: // Initialization of the ouput:
double *D; double *D;
if (nrhs == 3) if (nrhs == 3)
{ {
plhs[0] = mxCreateDoubleMatrix(mA,nB*nC,mxREAL); plhs[0] = mxCreateDoubleMatrix(mA, nB*nC, mxREAL);
} }
else else
{ {
plhs[0] = mxCreateDoubleMatrix(mA,nB*nB,mxREAL); plhs[0] = mxCreateDoubleMatrix(mA, nB*nB, mxREAL);
} }
D = mxGetPr(plhs[0]); D = mxGetPr(plhs[0]);
// Computational part: // Computational part:
if (nrhs == 2) if (nrhs == 2)
{ {
sparse_hessian_times_B_kronecker_B(isparseA, jsparseA, vsparseA, B, D, mA, nA, mB, nB); sparse_hessian_times_B_kronecker_B(isparseA, jsparseA, vsparseA, B, D, mA, nA, mB, nB);
} }
else else

142
mex/sources/kronecker/tests/test_kron.m
View File

@ -16,82 +16,82 @@ function test_kron(test)
% %
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>. % along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if ~nargin
test = 3;
end
if test == 1
if ~nargin percentage_of_non_zero_elements = 10e-4;
test = 3; NumberOfVariables = 549;%100;
NumberOfEquations = 256;%50
NumberOfColsInB = 50 ;
A = zeros(NumberOfEquations,NumberOfVariables^2);
for i = 1:NumberOfEquations
for j = 1:NumberOfVariables
for k = j:NumberOfVariables
if rand<percentage_of_non_zero_elements
A(i,(j-1)*NumberOfVariables+k) = randn;
end
end
for h = j+1:NumberOfVariables
A(i,(h-1)*NumberOfVariables+j) = A(i,(j-1)*NumberOfVariables+h);
end
end
end end
A = sparse(A);
B = randn(NumberOfVariables,NumberOfColsInB);
disp('')
disp('Computation of A*kron(B,B) with the mex file (v1):')
tic
D1 = sparse_hessian_times_B_kronecker_C(A,B);
toc
disp('')
disp('Computation of A*kron(B,B) with the mex file (v2):')
tic
D2 = sparse_hessian_times_B_kronecker_C(A,B,B);
toc
if test == 1 size(D1)
percentage_of_non_zero_elements = 10e-4; disp('');
NumberOfVariables = 549;%100; disp(['Difference between D1 and D2 = ' num2str(max(max(abs(D1-D2))))]);
NumberOfEquations = 256;%50
NumberOfColsInB = 50 ; return
A = zeros(NumberOfEquations,NumberOfVariables^2); disp(' ')
for i = 1:NumberOfEquations disp('Computation of A*kron(B,B) with two nested loops:')
for j = 1:NumberOfVariables tic
for k = j:NumberOfVariables D3 = zeros(NumberOfEquations,NumberOfColsInB*NumberOfColsInB);
if rand<percentage_of_non_zero_elements k = 1;
A(i,(j-1)*NumberOfVariables+k) = randn; for i1 = 1:NumberOfColsInB
end for i2 = 1:NumberOfColsInB
end D3(:,k) = A*kron(B(:,i1),B(:,i2));
for h = j+1:NumberOfVariables k = k+1;
A(i,(h-1)*NumberOfVariables+j) = A(i,(j-1)*NumberOfVariables+h);
end
end
end end
A = sparse(A); end
B = randn(NumberOfVariables,NumberOfColsInB); toc
disp('') disp('');
disp('Computation of A*kron(B,B) with the mex file (v1):') disp(['Difference between D1 and D3 = ' num2str(max(max(abs(D1-D3))))]);
tic
D1 = sparse_hessian_times_B_kronecker_C(A,B);
toc
disp('')
disp('Computation of A*kron(B,B) with the mex file (v2):')
tic
D2 = sparse_hessian_times_B_kronecker_C(A,B,B);
toc
size(D1)
disp('');
disp(['Difference between D1 and D2 = ' num2str(max(max(abs(D1-D2))))]);
return
disp(' ')
disp('Computation of A*kron(B,B) with two nested loops:')
tic
D3 = zeros(NumberOfEquations,NumberOfColsInB*NumberOfColsInB);
k = 1;
for i1 = 1:NumberOfColsInB
for i2 = 1:NumberOfColsInB
D3(:,k) = A*kron(B(:,i1),B(:,i2));
k = k+1;
end
end
toc
disp('');
disp(['Difference between D1 and D3 = ' num2str(max(max(abs(D1-D3))))]);
disp(' ') disp(' ')
disp('Direct computation of A*kron(B,B):') disp('Direct computation of A*kron(B,B):')
tic tic
try try
D4 = A*kron(B,B); D4 = A*kron(B,B);
notest = 0; notest = 0;
catch catch
notest = 1; notest = 1;
disp('Out of memory') disp('Out of memory')
end end
toc toc
if ~notest if ~notest
disp(''); disp('');
disp(['Difference between D1 and D4 = ' num2str(max(max(abs(D1-D4))))]); disp(['Difference between D1 and D4 = ' num2str(max(max(abs(D1-D4))))]);
end end
end end
@ -113,13 +113,13 @@ if test > 1
disp('') disp('')
disp('Computation of A*kron(B,B) with the mex file (v1):') disp('Computation of A*kron(B,B) with the mex file (v1):')
tic tic
D1 = sparse_hessian_times_B_kronecker_C(hessian,zx); D1 = sparse_hessian_times_B_kronecker_C(hessian,zx);
toc toc
disp('') disp('')
disp('Computation of A*kron(B,B) with the mex file (v2):') disp('Computation of A*kron(B,B) with the mex file (v2):')
tic tic
D2 = sparse_hessian_times_B_kronecker_C(hessian,zx,zx); D2 = sparse_hessian_times_B_kronecker_C(hessian,zx,zx);
toc toc
disp(''); disp('');

299
mex/sources/mjdgges/mjdgges.c
View File

@ -1,147 +1,152 @@
/* /*
* Copyright (C) 2006-2008 Dynare Team * Copyright (C) 2006-2008 Dynare Team
* *
* This file is part of Dynare. * This file is part of Dynare.
* *
* Dynare is free software: you can redistribute it and/or modify * Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by * it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or * the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version. * (at your option) any later version.
* *
* Dynare is distributed in the hope that it will be useful, * Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of * but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details. * GNU General Public License for more details.
* *
* You should have received a copy of the GNU General Public License * You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>. * along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/ */
#include <string.h> #include <string.h>
#include <dynmex.h> #include <dynmex.h>
#include <dynlapack.h> #include <dynlapack.h>
double criterium; double criterium;
lapack_int lapack_int
my_criteria(const double *alphar, const double *alphai, const double *beta) my_criteria(const double *alphar, const double *alphai, const double *beta)
{ {
return( (*alphar * *alphar + *alphai * *alphai) < criterium * *beta * *beta); return ((*alphar * *alphar + *alphai * *alphai) < criterium * *beta * *beta);
} }
void mjdgges(double *a, double *b, double *z, double *n, double *sdim, double *eval_r, double *eval_i, double *info) void
{ mjdgges(double *a, double *b, double *z, double *n, double *sdim, double *eval_r, double *eval_i, double *info)
lapack_int i_n, i_info, i_sdim, one, lwork; {
double *alphar, *alphai, *beta, *work, *par, *pai, *pb, *per, *pei; lapack_int i_n, i_info, i_sdim, one, lwork;
double *junk; double *alphar, *alphai, *beta, *work, *par, *pai, *pb, *per, *pei;
lapack_int *bwork; double *junk;
lapack_int *bwork;
one = 1;
i_n = (lapack_int)*n; one = 1;
alphar = mxCalloc(i_n,sizeof(double)); i_n = (lapack_int)*n;
alphai = mxCalloc(i_n,sizeof(double)); alphar = mxCalloc(i_n, sizeof(double));
beta = mxCalloc(i_n,sizeof(double)); alphai = mxCalloc(i_n, sizeof(double));
lwork = 16*i_n+16; beta = mxCalloc(i_n, sizeof(double));
work = mxCalloc(lwork,sizeof(double)); lwork = 16*i_n+16;
bwork = mxCalloc(i_n,sizeof(lapack_int)); work = mxCalloc(lwork, sizeof(double));
/* made necessary by bug in Lapack */ bwork = mxCalloc(i_n, sizeof(lapack_int));
junk = mxCalloc(i_n*i_n,sizeof(double)); /* made necessary by bug in Lapack */
junk = mxCalloc(i_n*i_n, sizeof(double));
dgges("N", "V", "S", my_criteria, &i_n, a, &i_n, b, &i_n, &i_sdim, alphar, alphai, beta, junk, &i_n, z, &i_n, work, &lwork, bwork, &i_info);
dgges("N", "V", "S", my_criteria, &i_n, a, &i_n, b, &i_n, &i_sdim, alphar, alphai, beta, junk, &i_n, z, &i_n, work, &lwork, bwork, &i_info);
*sdim = i_sdim;
*info = i_info; *sdim = i_sdim;
*info = i_info;
par = alphar;
pai = alphai; par = alphar;
pb = beta; pai = alphai;
pei = eval_i; pb = beta;
for(per = eval_r; per <= &eval_r[i_n-1]; ++per) pei = eval_i;
{ for (per = eval_r; per <= &eval_r[i_n-1]; ++per)
*per = *par / *pb; {
*pei = *pai / *pb; *per = *par / *pb;
++par; *pei = *pai / *pb;
++pai; ++par;
++pb; ++pai;
++pei; ++pb;
} ++pei;
} }
}
/* MATLAB interface */
void mexFunction( int nlhs, mxArray *plhs[], /* MATLAB interface */
int nrhs, const mxArray *prhs[] ) void
mexFunction(int nlhs, mxArray *plhs[],
{ int nrhs, const mxArray *prhs[])
unsigned int m1,n1,m2,n2;
double *s, *t, *z, *sdim, *eval_r, *eval_i, *info, *a, *b; {
double n; unsigned int m1, n1, m2, n2;
double *s, *t, *z, *sdim, *eval_r, *eval_i, *info, *a, *b;
/* Check for proper number of arguments */ double n;
if (nrhs < 2 || nrhs > 3) { /* Check for proper number of arguments */
mexErrMsgTxt("MJDGGES: two or three input arguments are required.");
} else if (nlhs > 6) { if (nrhs < 2 || nrhs > 3)
mexErrMsgTxt("MJDGGES: too many output arguments."); {
} mexErrMsgTxt("MJDGGES: two or three input arguments are required.");
}
/* Check that A and B are real matrices of the same dimension.*/ else if (nlhs > 6)
{
m1 = mxGetM(prhs[0]); mexErrMsgTxt("MJDGGES: too many output arguments.");
n1 = mxGetN(prhs[0]); }
m2 = mxGetM(prhs[1]);
n2 = mxGetN(prhs[1]); /* Check that A and B are real matrices of the same dimension.*/
if (!mxIsDouble(prhs[0]) || mxIsComplex(prhs[0]) ||
!mxIsDouble(prhs[1]) || mxIsComplex(prhs[1]) || m1 = mxGetM(prhs[0]);
(m1 != n1) || (m2!= n1) || (m2 != n2)) { n1 = mxGetN(prhs[0]);
mexErrMsgTxt("MJDGGES requires two square real matrices of the same dimension."); m2 = mxGetM(prhs[1]);
} n2 = mxGetN(prhs[1]);
if (!mxIsDouble(prhs[0]) || mxIsComplex(prhs[0])
/* Create a matrix for the return argument */ || !mxIsDouble(prhs[1]) || mxIsComplex(prhs[1])
plhs[0] = mxCreateDoubleMatrix(n1, n1, mxREAL); || (m1 != n1) || (m2 != n1) || (m2 != n2))
plhs[1] = mxCreateDoubleMatrix(n1, n1, mxREAL); {
plhs[2] = mxCreateDoubleMatrix(n1, n1, mxREAL); mexErrMsgTxt("MJDGGES requires two square real matrices of the same dimension.");
plhs[3] = mxCreateDoubleMatrix(1, 1, mxREAL); }
plhs[4] = mxCreateDoubleMatrix(n1, 1, mxCOMPLEX);
plhs[5] = mxCreateDoubleMatrix(1, 1, mxREAL); /* Create a matrix for the return argument */
plhs[0] = mxCreateDoubleMatrix(n1, n1, mxREAL);
/* Assign pointers to the various parameters */ plhs[1] = mxCreateDoubleMatrix(n1, n1, mxREAL);
s = mxGetPr(plhs[0]); plhs[2] = mxCreateDoubleMatrix(n1, n1, mxREAL);
t = mxGetPr(plhs[1]); plhs[3] = mxCreateDoubleMatrix(1, 1, mxREAL);
z = mxGetPr(plhs[2]); plhs[4] = mxCreateDoubleMatrix(n1, 1, mxCOMPLEX);
sdim = mxGetPr(plhs[3]); plhs[5] = mxCreateDoubleMatrix(1, 1, mxREAL);
eval_r = mxGetPr(plhs[4]);
eval_i = mxGetPi(plhs[4]); /* Assign pointers to the various parameters */
info = mxGetPr(plhs[5]); s = mxGetPr(plhs[0]);
t = mxGetPr(plhs[1]);
a = mxGetPr(prhs[0]); z = mxGetPr(plhs[2]);
b = mxGetPr(prhs[1]); sdim = mxGetPr(plhs[3]);
eval_r = mxGetPr(plhs[4]);
/* set criterium for stable eigenvalues */ eval_i = mxGetPi(plhs[4]);
if ( nrhs == 3) info = mxGetPr(plhs[5]);
{
criterium = *mxGetPr(prhs[2]); a = mxGetPr(prhs[0]);
} b = mxGetPr(prhs[1]);
else
{ /* set criterium for stable eigenvalues */
criterium = 1+1e-6; if (nrhs == 3)
} {
criterium = *mxGetPr(prhs[2]);
/* keep a and b intact */ }
memcpy(s,a,sizeof(double)*n1*n1); else
memcpy(t,b,sizeof(double)*n1*n1); {
criterium = 1+1e-6;
n = n1; }
/* Do the actual computations in a subroutine */ /* keep a and b intact */
mjdgges(s, t, z, &n, sdim, eval_r, eval_i, info); memcpy(s, a, sizeof(double)*n1*n1);
memcpy(t, b, sizeof(double)*n1*n1);
return; n = n1;
} /* Do the actual computations in a subroutine */
mjdgges(s, t, z, &n, sdim, eval_r, eval_i, info);
/*
07/30/03 MJ added user set criterium for stable eigenvalues return;
corrected error messages in mexfunction()
*/ }
/*
07/30/03 MJ added user set criterium for stable eigenvalues
corrected error messages in mexfunction()
*/