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

Adaptation of simulate to matlab 64 bits 

git-svn-id: https://www.dynare.org/svn/dynare/trunk@2874 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
ferhat 2009-08-29 16:04:06 +00:00
parent ed6cda0ca1
commit 7a6df1c342
7 changed files with 208 additions and 196 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
mex/sources/simulate/Interpreter.cc
View File

@ -18,7 +18,7 @@
*/ */
#include <cstring> #include <cstring>
#include <sstream> #include <sstream>
#include "Interpreter.hh" #include "Interpreter.hh"
#define BIG 1.0e+8; #define BIG 1.0e+8;
#define SMALL 1.0e-5; #define SMALL 1.0e-5;
//#define DEBUG //#define DEBUG
@ -811,7 +811,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
Block_List_Max_Lead=get_code_int; Block_List_Max_Lead=get_code_int;
u_count_int=get_code_int; u_count_int=get_code_int;
fixe_u(&u, u_count_int, u_count_int); fixe_u(&u, u_count_int, u_count_int);
Read_SparseMatrix(bin_basename, size, 1, 0, 0, steady_state); Read_SparseMatrix(bin_basename, size, 1, 0, 0, steady_state, false);
g1=(double*)mxMalloc(size*size*sizeof(double)); g1=(double*)mxMalloc(size*size*sizeof(double));
r=(double*)mxMalloc(size*sizeof(double)); r=(double*)mxMalloc(size*sizeof(double));
begining=get_code_pointer; begining=get_code_pointer;
@ -966,7 +966,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
Block_List_Max_Lead=get_code_int; Block_List_Max_Lead=get_code_int;
u_count_int=get_code_int; u_count_int=get_code_int;
fixe_u(&u, u_count_int, u_count_int); fixe_u(&u, u_count_int, u_count_int);
Read_SparseMatrix(bin_basename, size, 1, 0, 0, steady_state); Read_SparseMatrix(bin_basename, size, 1, 0, 0, steady_state, false);
g1=(double*)mxMalloc(size*size*sizeof(double)); g1=(double*)mxMalloc(size*size*sizeof(double));
r=(double*)mxMalloc(size*sizeof(double)); r=(double*)mxMalloc(size*sizeof(double));
begining=get_code_pointer; begining=get_code_pointer;
@ -1121,7 +1121,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
Block_List_Max_Lead=get_code_int; Block_List_Max_Lead=get_code_int;
u_count_int=get_code_int; u_count_int=get_code_int;
fixe_u(&u, u_count_int, u_count_int); fixe_u(&u, u_count_int, u_count_int);
Read_SparseMatrix(bin_basename, size, periods, y_kmin, y_kmax, steady_state); Read_SparseMatrix(bin_basename, size, periods, y_kmin, y_kmax, steady_state, true);
u_count=u_count_int*(periods+y_kmax+y_kmin); u_count=u_count_int*(periods+y_kmax+y_kmin);
r=(double*)mxMalloc(size*sizeof(double)); r=(double*)mxMalloc(size*sizeof(double));
y_save=(double*)mxMalloc(y_size*sizeof(double)*(periods+y_kmax+y_kmin)); y_save=(double*)mxMalloc(y_size*sizeof(double)*(periods+y_kmax+y_kmin));
@ -1230,7 +1230,8 @@ Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_s
{ {
ifstream CompiledCode; ifstream CompiledCode;
bool result = true; bool result = true;
int Code_Size, var; streamoff Code_Size;
int var;
if(steady_state) if(steady_state)
file_name += "_static"; file_name += "_static";
else else
@ -1271,7 +1272,7 @@ Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_s
Block.clear(); Block.clear();
Block_Contain.clear(); Block_Contain.clear();
Block_contain_type lBlock_Contain; Block_contain_type lBlock_Contain;
lBlock.begin=get_code_pos-(uint64_t)Init_Code; lBlock.begin=get_code_pos-(long int*)Init_Code;
lBlock.size=get_code_int; lBlock.size=get_code_int;
lBlock.type=get_code_int; lBlock.type=get_code_int;
Block.push_back(lBlock); Block.push_back(lBlock);
@ -1302,7 +1303,7 @@ Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_s
T=(double*)mxMalloc(var*sizeof(double)); T=(double*)mxMalloc(var*sizeof(double));
break; break;
default : default :
mexPrintf("Unknow command : %d at pos %d !!\n",(long int)(code),(uint64_t*)(get_code_pos)-(uint64_t*)(Init_Code)); mexPrintf("Unknow command : %d at pos %d !!\n",(long int)(code),(long int*)(get_code_pos)-(long int*)(Init_Code));
mexEvalString("st=fclose('all');clear all;"); mexEvalString("st=fclose('all');clear all;");
mexEvalString("drawnow;"); mexEvalString("drawnow;");
mexErrMsgTxt("End of simulate"); mexErrMsgTxt("End of simulate");

2
mex/sources/simulate/Interpreter.hh
View File

@ -58,7 +58,7 @@ struct Block_type
#define get_code_pdouble (double*)Code; Code+=sizeof(double); #define get_code_pdouble (double*)Code; Code+=sizeof(double);
#define get_code_bool *((bool*)(Code++)) #define get_code_bool *((bool*)(Code++))
#define get_code_char *((char*)(Code++)) #define get_code_char *((char*)(Code++))
#define get_code_pos (long int)Code #define get_code_pos (long int*)Code
#define get_code_pointer Code #define get_code_pointer Code
#define set_code_pointer(pos) Code=pos #define set_code_pointer(pos) Code=pos

127
mex/sources/simulate/Mem_Mngr.cc
View File

@ -16,7 +16,7 @@
* 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 <stdint.h> //#include "stdint.h"
#include "Mem_Mngr.hh" #include "Mem_Mngr.hh"
Mem_Mngr::Mem_Mngr() Mem_Mngr::Mem_Mngr()
@ -103,10 +103,11 @@ Mem_Mngr::mxMalloc_NZE()
void void
Mem_Mngr::mxFree_NZE(void* pos) Mem_Mngr::mxFree_NZE(void* pos)
{ {
int i, gap; int i;
size_t gap;
for (i=0;i<Nb_CHUNK;i++) for (i=0;i<Nb_CHUNK;i++)
{ {
gap=((uint64_t)(pos)-(uint64_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;
} }
@ -114,126 +115,6 @@ Mem_Mngr::mxFree_NZE(void* pos)
} }
void
Mem_Mngr::write_swp_f(int *save_op_all,long int *nop_all)
{
swp_f=true;
swp_f_b++;
mexPrintf("writing the block %d with size=%d\n",swp_f_b,*nop_all);
if (!SaveCode_swp.is_open())
{
mexPrintf("open the swp file for writing\n");
SaveCode_swp.open((filename + ".swp").c_str(), std::ios::out | std::ios::binary);
if (!SaveCode_swp.is_open())
{
mexPrintf("Error : Can't open file \"%s\" for writing\n", (filename + ".swp").c_str());
mexEvalString("st=fclose('all');clear all;");
mexErrMsgTxt("Exit from Dynare");
}
}
SaveCode_swp.write(reinterpret_cast<char *>(nop_all), sizeof(*nop_all));
SaveCode_swp.write(reinterpret_cast<char *>(save_op_all), (*nop_all)*sizeof(int));
(*nop_all)=0;
}
bool
Mem_Mngr::read_swp_f(int **save_op_all,long int *nop_all)
{
int j;
swp_f=true;
if (!SaveCode_swp.is_open())
{
mexPrintf("open the file %s\n",(filename + ".swp").c_str());
SaveCode_swp.open((filename + ".swp").c_str(), std::ios::in | std::ios::binary);
j=SaveCode_swp.is_open();
mexPrintf("is_open()=%d\n",j);
if (!SaveCode_swp.is_open())
{
mexPrintf("Error : Can't open file \"%s\" for reading\n", (filename + ".swp").c_str());
mexEvalString("st=fclose('all');clear all;");
mexErrMsgTxt("Exit from Dynare");
}
SaveCode_swp.seekg(0);
}
j=SaveCode_swp.tellg();
SaveCode_swp.read(reinterpret_cast<char *>(nop_all), sizeof(*nop_all));
(*save_op_all)=(int*)mxMalloc((*nop_all)*sizeof(int));
SaveCode_swp.read(reinterpret_cast<char *>(*save_op_all), (*nop_all)*sizeof(int));
return(SaveCode_swp.good());
}
void
Mem_Mngr::close_swp_f()
{
if (SaveCode_swp.is_open())
{
SaveCode_swp.close();
mexPrintf("close the swp file\n");
}
}
int*
Mem_Mngr::malloc_std(long int nop)
{
return((int*)malloc(nop*sizeof(int)));
}
int*
Mem_Mngr::realloc_std(int* save_op_o, long int &nopa)
{
int *save_op=(int*)realloc(save_op_o,nopa*sizeof(int));
if (!save_op)
{
int nopag=int(nopa/3);
nopa=nopa-nopag;
while (!save_op && nopag>0)
{
nopag=int(nopag*0.66);
save_op=(int*)realloc(save_op_o,nopa*sizeof(int));
}
if (!save_op)
{
mexPrintf("Memory exhausted\n");
mexEvalString("drawnow;");
mexEvalString("st=fclose('all');clear all;");
filename+=" stopped";
mexErrMsgTxt(filename.c_str());
}
}
return(save_op);
}
void
Mem_Mngr::chk_avail_mem(int **save_op_all,long int *nop_all,long int *nopa_all,int add, int t)
{
mexPrintf("Error: out of save_op_all[%d] nopa_all=%d t=%d\n",(*nop_all)+add,(*nopa_all),t);
int tmp_nopa_all=int(1.5*(*nopa_all));
int *tmp_i;
if (tmp_nopa_all*sizeof(int)<1024*1024)
{
mexPrintf("allocate %d bites save_op_all=%x\n",tmp_nopa_all*sizeof(int),*save_op_all);
tmp_i=(int*)mxRealloc(*save_op_all,tmp_nopa_all*sizeof(int));
mexPrintf("tmp_i=");
mexPrintf("%x\n",tmp_i);
}
else
tmp_i=NULL;
if (!tmp_i)
{
write_swp_f((*save_op_all),nop_all);
}
else
{
mexPrintf("allocated\n");
(*save_op_all)=tmp_i;
(*nopa_all)=tmp_nopa_all;
}
mexPrintf("end of chk\n");
}
void void
Mem_Mngr::Free_All() Mem_Mngr::Free_All()
{ {

8
mex/sources/simulate/Mem_Mngr.hh
View File

@ -42,9 +42,9 @@ typedef vector<NonZeroElem*> v_NonZeroElem;
class Mem_Mngr class Mem_Mngr
{ {
public: public:
void write_swp_f(int *save_op_all,long int *nop_all); /*void write_swp_f(int *save_op_all,long int *nop_all);
bool read_swp_f(int **save_op_all,long int *nop_all); bool read_swp_f(int **save_op_all,long int *nop_all);
void close_swp_f(); void close_swp_f();*/
void Print_heap(); void Print_heap();
void init_Mem(); void init_Mem();
void mxFree_NZE(void* pos); void mxFree_NZE(void* pos);
@ -53,9 +53,9 @@ public:
void Free_All(); void Free_All();
Mem_Mngr(); Mem_Mngr();
void fixe_file_name(string filename_arg); void fixe_file_name(string filename_arg);
int* malloc_std(long int nop); /*int* malloc_std(long int nop);
int* realloc_std(int* save_op_o, long int &nopa); int* realloc_std(int* save_op_o, long int &nopa);
void chk_avail_mem(int **save_op_all,long int *nop_all,long int *nopa_all,int add, int t); void chk_avail_mem(int **save_op_all,long int *nop_all,long int *nopa_all,int add, int t);*/
bool swp_f; bool swp_f;
//bool verbose; //bool verbose;
private: private:

202
mex/sources/simulate/SparseMatrix.cc
View File

@ -39,7 +39,12 @@ SparseMatrix::SparseMatrix()
tbreak_g = 0; tbreak_g = 0;
start_compare = 0; start_compare = 0;
restart = 0; restart = 0;
IM_i.clear(); IM_i.clear();
#ifdef _MSC_VER
nan__[0] = 0xffffffff;
nan__[1] = 0x7fffffff;
NAN = *( double* )nan__;
#endif
} }
int int
@ -339,7 +344,7 @@ SparseMatrix::Insert(const int r, const int c, const int u_index, const int lag_
} }
void void
SparseMatrix::Read_SparseMatrix(string file_name, int Size, int periods, int y_kmin, int y_kmax, bool steady_state) SparseMatrix::Read_SparseMatrix(string file_name, int Size, int periods, int y_kmin, int y_kmax, bool steady_state, bool two_boundaries)
{ {
int i, j, eq, var, lag; int i, j, eq, var, lag;
filename = file_name; filename = file_name;
@ -361,14 +366,47 @@ SparseMatrix::Read_SparseMatrix(string file_name, int Size, int periods, int y_k
} }
} }
IM_i.clear(); IM_i.clear();
for (i = 0; i < u_count_init; i++) if(two_boundaries)
{ {
SaveCode.read(reinterpret_cast<char *>(&eq), sizeof(eq)); for (i = 0; i < u_count_init-Size; i++)
SaveCode.read(reinterpret_cast<char *>(&var), sizeof(var)); {
SaveCode.read(reinterpret_cast<char *>(&lag), sizeof(lag)); SaveCode.read(reinterpret_cast<char *>(&eq), sizeof(eq));
SaveCode.read(reinterpret_cast<char *>(&j), sizeof(j)); SaveCode.read(reinterpret_cast<char *>(&var), sizeof(var));
IM_i[make_pair(make_pair(eq, var), lag)] = j; //mexPrintf("var=%d\n",var);
} SaveCode.read(reinterpret_cast<char *>(&lag), sizeof(lag));
SaveCode.read(reinterpret_cast<char *>(&j), sizeof(j));
IM_i[make_pair(make_pair(eq, var), lag)] = j;
}
/*
int eqr1=j;
int varr=Size*(block_triangular.periods
+block_triangular.incidencematrix.Model_Max_Lead_Endo);
int k1=0;
SaveCode.write(reinterpret_cast<char *>(&eqr1), sizeof(eqr1));
SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr));
SaveCode.write(reinterpret_cast<char *>(&k1), sizeof(k1));
SaveCode.write(reinterpret_cast<char *>(&eqr1), sizeof(eqr1));
u_count_int++;
*/
for (j=0;j<Size;j++)
{
//mexPrintf("var = %d\n",Size*(periods+y_kmax));
IM_i[make_pair(make_pair(j, Size*(periods+y_kmax)), 0)] = j;
//u_count_init++;
}
}
else
{
for (i = 0; i < u_count_init; i++)
{
SaveCode.read(reinterpret_cast<char *>(&eq), sizeof(eq));
SaveCode.read(reinterpret_cast<char *>(&var), sizeof(var));
//mexPrintf("var=%d\n",var);
SaveCode.read(reinterpret_cast<char *>(&lag), sizeof(lag));
SaveCode.read(reinterpret_cast<char *>(&j), sizeof(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));
@ -985,12 +1023,12 @@ SparseMatrix::complete(int beg_t, int Size, int periods, int *b)
return (beg_t); return (beg_t);
} }
void /*void
SparseMatrix::close_swp_file() SparseMatrix::close_swp_file()
{ {
mem_mngr.close_swp_f(); mem_mngr.close_swp_f();
} }
*/
double double
SparseMatrix::bksub(int tbreak, int last_period, int Size, double slowc_l) SparseMatrix::bksub(int tbreak, int last_period, int Size, double slowc_l)
{ {
@ -1074,10 +1112,15 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
int pivj = 0, pivk = 0; int pivj = 0, pivk = 0;
double piv_abs/*, first_elem*/; double piv_abs/*, first_elem*/;
NonZeroElem *first, *firsta, *first_suba; NonZeroElem *first, *firsta, *first_suba;
double piv_v[Size]; double *piv_v;
int pivj_v[Size], pivk_v[Size], NR[Size], l, N_max; int *pivj_v, *pivk_v, *NR;
int l, N_max;
bool one; bool one;
Clear_u(); Clear_u();
piv_v = (double*)mxMalloc(Size*sizeof(double));
pivj_v = (int*)mxMalloc(Size*sizeof(int));
pivk_v = (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))
@ -1101,7 +1144,11 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
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(pivj_v);
mxFree(pivk_v);
mxFree(NR);
if(steady_state) if(steady_state)
return false; return false;
else else
@ -1129,21 +1176,29 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
} }
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;");
if(steady_state) mxFree(piv_v);
return false; mxFree(pivj_v);
else mxFree(pivk_v);
{ mxFree(NR);
if(steady_state)
return false;
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());
} }
} }
slowc_save /= 2; slowc_save /= 2;
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];
/*for (i = 0; i < y_size*(periods+y_kmin); i++) /*for (i = 0; i < y_size*(periods+y_kmin); i++)
y[i] = ya[i]+slowc_save*direction[i];*/ y[i] = ya[i]+slowc_save*direction[i];*/
mxFree(piv_v);
mxFree(pivj_v);
mxFree(pivk_v);
mxFree(NR);
iter--; iter--;
return true; return true;
} }
@ -1156,9 +1211,17 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
mexPrintf(" abs. error=%.10e \n", double (res1)); mexPrintf(" abs. error=%.10e \n", double (res1));
mexPrintf("-----------------------------------\n"); mexPrintf("-----------------------------------\n");
}*/ }*/
if (cvg) if (cvg)
return (true); {
Simple_Init(it_, y_kmin, y_kmax, Size, IM_i); mxFree(piv_v);
mxFree(pivj_v);
mxFree(pivk_v);
mxFree(NR);
return (true);
}
Simple_Init(it_, y_kmin, y_kmax, Size, IM_i);
NonZeroElem **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*/
@ -1243,15 +1306,20 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
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(pivj_v);
mxFree(pivk_v);
mxFree(NR);
mxFree(bc);
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());
} }
} }
/*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);
@ -1265,14 +1333,13 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
nb_eq = At_Col(i, &first); nb_eq = At_Col(i, &first);
NonZeroElem *first_piva; NonZeroElem *first_piva;
int nb_var_piva = At_Row(pivj, &first_piva); int nb_var_piva = At_Row(pivj, &first_piva);
NonZeroElem *bc[nb_eq]; 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++)
{ {
if (!line_done[first->r_index]) if (!line_done[first->r_index])
bc[nb_eq_todo++] = first; bc[nb_eq_todo++] = first;
first = first->NZE_C_N; first = first->NZE_C_N;
} }
//#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++)
{ {
@ -1352,14 +1419,19 @@ 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;
first = first->NZE_C_N; first = first->NZE_C_N;
} }
} }
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(pivj_v);
mxFree(pivk_v);
mxFree(NR);
mxFree(bc);
return true; return true;
} }
@ -1419,7 +1491,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[row]; double *B;
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();
@ -1429,7 +1502,8 @@ 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);
} }
void void
@ -1581,16 +1655,22 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
} }
else else
{ {
Init(periods, y_kmin, y_kmax, Size, IM_i); Init(periods, y_kmin, y_kmax, Size, IM_i);
double *piv_v;
int *pivj_v, *pivk_v, *NR;
piv_v = (double*)mxMalloc(Size*sizeof(double));
pivj_v = (int*)mxMalloc(Size*sizeof(int));
pivk_v = (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)
{ {
if (save_op) ; if (save_op)
{ {
mxFree(save_op); mxFree(save_op);
save_op = NULL; save_op = NULL;
} }
save_op = (int *) mxMalloc(nop*sizeof(int)); save_op = (int *) mxMalloc(nop*sizeof(int));
nopa = nop; nopa = nop;
} }
@ -1604,8 +1684,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
int nb_eq = At_Col(i, 0, &first); int nb_eq = At_Col(i, 0, &first);
if ((symbolic && t <= start_compare) || !symbolic) if ((symbolic && t <= start_compare) || !symbolic)
{ {
double piv_v[Size]; int l = 0, N_max = 0;
int pivj_v[Size], pivk_v[Size], NR[Size], l = 0, N_max = 0;
bool one = false; bool one = false;
piv_abs = 0; piv_abs = 0;
for (j = 0; j < nb_eq; j++) for (j = 0; j < nb_eq; j++)
@ -1666,7 +1745,7 @@ SparseMatrix::simulate_NG1(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(NR[j]/N_max)); markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log((double)(NR[j]/N_max)));
if (markovitz > markovitz_max && NR[j] == 1) if (markovitz > markovitz_max && NR[j] == 1)
{ {
piv = piv_v[j]; piv = piv_v[j];
@ -1717,10 +1796,14 @@ 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[nb_var]; NonZeroElem **bb;
bb = (NonZeroElem**)mxMalloc(nb_var*sizeof(first));
//mexPrintf("nb_var=%d\n",nb_var);
for (j = 0; j < nb_var; j++) for (j = 0; j < nb_var; j++)
{ {
bb[j] = first; bb[j] = first;
/*if(nb_var==122)
mexPrintf("j=%d first=%x\n",j,first);*/
first = first->NZE_R_N; first = first->NZE_R_N;
} }
@ -1743,7 +1826,8 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
save_op_s->lag = first->lag_index; save_op_s->lag = first->lag_index;
} }
} }
} }
mxFree(bb);
nop += nb_var*2; nop += nb_var*2;
u[b[pivj]] /= piv; u[b[pivj]] /= piv;
if (symbolic) if (symbolic)
@ -1767,7 +1851,9 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
NonZeroElem *first_piva; NonZeroElem *first_piva;
int nb_var_piva = At_Row(pivj, &first_piva); int nb_var_piva = At_Row(pivj, &first_piva);
NonZeroElem *bc[nb_eq]; NonZeroElem **bc;
bc = (NonZeroElem**)mxMalloc(nb_eq*sizeof(first));
//NonZeroElem *bc[nb_eq];
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++)
{ {
@ -1949,7 +2035,8 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
} }
nop += 3; nop += 3;
} }
} }
mxFree(bc);
} }
if (symbolic) if (symbolic)
{ {
@ -2006,6 +2093,10 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
} }
//record = true; //record = true;
} }
mxFree(piv_v);
mxFree(pivj_v);
mxFree(pivk_v);
mxFree(NR);
} }
nop_all += nop; nop_all += nop;
if (symbolic) if (symbolic)
@ -2017,7 +2108,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
if (save_opaa) if (save_opaa)
mxFree(save_opaa); mxFree(save_opaa);
} }
close_swp_file(); //close_swp_file();
/*The backward substitution*/ /*The backward substitution*/
double slowc_lbx = slowc, res1bx; double slowc_lbx = slowc, res1bx;
@ -2034,7 +2125,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
mexEvalString("drawnow;"); mexEvalString("drawnow;");
} }
close_swp_file(); //close_swp_file();
time00 = clock(); time00 = clock();
if (tbreak_g == 0) if (tbreak_g == 0)
tbreak_g = periods; tbreak_g = periods;
@ -2043,7 +2134,6 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
/*Check_the_Solution(periods, y_kmin, y_kmax, Size, ua, pivot, b); /*Check_the_Solution(periods, y_kmin, y_kmax, Size, ua, pivot, b);
mxFree(ua);*/ mxFree(ua);*/
return (0); return (0);
} }

49
mex/sources/simulate/SparseMatrix.hh
View File

@ -30,11 +30,17 @@
// Nothing to single thread version. // Nothing to single thread version.
#else #else
#include <omp.h> #include <omp.h>
#endif
#ifdef _MSC_VER
#include <limits>
#endif #endif
#define NEW_ALLOC #define NEW_ALLOC
#define MARKOVITZ #define MARKOVITZ
using namespace std; using namespace std;
struct t_save_op_s struct t_save_op_s
{ {
@ -63,8 +69,45 @@ class SparseMatrix
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, int Size, int periods, int y_kmin, int y_kmax, bool steady_state); void Read_SparseMatrix(string file_name, 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);
#ifdef _MSC_VER
unsigned long nan__[2];
double NAN;
inline bool isnan(double value)
{
return value != value;
}
inline bool isinf(double value)
{
return (std::numeric_limits<double>::has_infinity &&
value == std::numeric_limits<double>::infinity());
}
inline double asinh(double x)
{
return log(x+sqrt(x*x+1));
}
template<typename T>
inline T acosh(T x)
{
if(!(x>=1.0)) return sqrt(-1.0);
return log(x+sqrt(x*x-1.0));
}
template<typename T>
inline T atanh(T x)
{
if(!(x>-1.0 && x<1.0)) return sqrt(-1.0);
return log((1.0+x)/(1.0-x))/2.0;
}
#endif
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);
@ -99,7 +142,7 @@ class SparseMatrix
#endif #endif
); );
double simple_bksub(int it_, int Size, double slowc_l); double simple_bksub(int it_, int Size, double slowc_l);
void close_swp_file(); /*void close_swp_file();*/
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;
@ -148,6 +191,4 @@ protected:
#endif #endif

1
mex/sources/simulate/simulate.cc
View File

@ -24,7 +24,6 @@
//////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////
#include <cstring> #include <cstring>
#include "simulate.hh" #include "simulate.hh"
#include "Interpreter.hh" #include "Interpreter.hh"
#include "Mem_Mngr.hh" #include "Mem_Mngr.hh"