From 7a6df1c3427ed3ae563fa5edb65ea54285ad5ee9 Mon Sep 17 00:00:00 2001
From: ferhat <ferhat@ac1d8469-bf42-47a9-8791-bf33cf982152>
Date: Sat, 29 Aug 2009 16:04:06 +0000
Subject: [PATCH] Adaptation of simulate to matlab 64 bits
git-svn-id: https://www.dynare.org/svn/dynare/trunk@2874 ac1d8469-bf42-47a9-8791-bf33cf982152
---
mex/sources/simulate/Interpreter.cc | 15 +-
mex/sources/simulate/Interpreter.hh | 2 +-
mex/sources/simulate/Mem_Mngr.cc | 127 +----------------
mex/sources/simulate/Mem_Mngr.hh | 8 +-
mex/sources/simulate/SparseMatrix.cc | 202 +++++++++++++++++++--------
mex/sources/simulate/SparseMatrix.hh | 49 ++++++-
mex/sources/simulate/simulate.cc | 1 -
7 files changed, 208 insertions(+), 196 deletions(-)
diff --git a/mex/sources/simulate/Interpreter.cc b/mex/sources/simulate/Interpreter.cc
index e455d14ef..6441a339d 100644
--- a/mex/sources/simulate/Interpreter.cc
+++ b/mex/sources/simulate/Interpreter.cc
@@ -18,7 +18,7 @@
*/
#include <cstring>
#include <sstream>
-#include "Interpreter.hh"
+#include "Interpreter.hh"
#define BIG 1.0e+8;
#define SMALL 1.0e-5;
//#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;
u_count_int=get_code_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));
r=(double*)mxMalloc(size*sizeof(double));
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;
u_count_int=get_code_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));
r=(double*)mxMalloc(size*sizeof(double));
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;
u_count_int=get_code_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);
r=(double*)mxMalloc(size*sizeof(double));
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;
bool result = true;
- int Code_Size, var;
+ streamoff Code_Size;
+ int var;
if(steady_state)
file_name += "_static";
else
@@ -1271,7 +1272,7 @@ Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_s
Block.clear();
Block_Contain.clear();
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.type=get_code_int;
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));
break;
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("drawnow;");
mexErrMsgTxt("End of simulate");
diff --git a/mex/sources/simulate/Interpreter.hh b/mex/sources/simulate/Interpreter.hh
index 66d7592cd..10f9256d7 100644
--- a/mex/sources/simulate/Interpreter.hh
+++ b/mex/sources/simulate/Interpreter.hh
@@ -58,7 +58,7 @@ struct Block_type
#define get_code_pdouble (double*)Code; Code+=sizeof(double);
#define get_code_bool *((bool*)(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 set_code_pointer(pos) Code=pos
diff --git a/mex/sources/simulate/Mem_Mngr.cc b/mex/sources/simulate/Mem_Mngr.cc
index ef2d92826..2b077d560 100644
--- a/mex/sources/simulate/Mem_Mngr.cc
+++ b/mex/sources/simulate/Mem_Mngr.cc
@@ -16,7 +16,7 @@
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
-#include <stdint.h>
+//#include "stdint.h"
#include "Mem_Mngr.hh"
Mem_Mngr::Mem_Mngr()
@@ -103,10 +103,11 @@ Mem_Mngr::mxMalloc_NZE()
void
Mem_Mngr::mxFree_NZE(void* pos)
{
- int i, gap;
+ int i;
+ size_t gap;
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))
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
Mem_Mngr::Free_All()
{
diff --git a/mex/sources/simulate/Mem_Mngr.hh b/mex/sources/simulate/Mem_Mngr.hh
index 5d1139608..7ade97786 100644
--- a/mex/sources/simulate/Mem_Mngr.hh
+++ b/mex/sources/simulate/Mem_Mngr.hh
@@ -42,9 +42,9 @@ typedef vector<NonZeroElem*> v_NonZeroElem;
class Mem_Mngr
{
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);
- void close_swp_f();
+ void close_swp_f();*/
void Print_heap();
void init_Mem();
void mxFree_NZE(void* pos);
@@ -53,9 +53,9 @@ public:
void Free_All();
Mem_Mngr();
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);
- 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 verbose;
private:
diff --git a/mex/sources/simulate/SparseMatrix.cc b/mex/sources/simulate/SparseMatrix.cc
index 20fa39b98..7617576a8 100644
--- a/mex/sources/simulate/SparseMatrix.cc
+++ b/mex/sources/simulate/SparseMatrix.cc
@@ -39,7 +39,12 @@ SparseMatrix::SparseMatrix()
tbreak_g = 0;
start_compare = 0;
restart = 0;
- IM_i.clear();
+ IM_i.clear();
+#ifdef _MSC_VER
+ nan__[0] = 0xffffffff;
+ nan__[1] = 0x7fffffff;
+ NAN = *( double* )nan__;
+#endif
}
int
@@ -339,7 +344,7 @@ SparseMatrix::Insert(const int r, const int c, const int u_index, const int lag_
}
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;
filename = file_name;
@@ -361,14 +366,47 @@ SparseMatrix::Read_SparseMatrix(string file_name, int Size, int periods, int y_k
}
}
IM_i.clear();
- for (i = 0; i < u_count_init; i++)
- {
- SaveCode.read(reinterpret_cast<char *>(&eq), sizeof(eq));
- SaveCode.read(reinterpret_cast<char *>(&var), sizeof(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;
- }
+ if(two_boundaries)
+ {
+ for (i = 0; i < u_count_init-Size; 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;
+ }
+ /*
+ 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));
for (j = 0; j < Size; j++)
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);
}
-void
+/*void
SparseMatrix::close_swp_file()
{
mem_mngr.close_swp_f();
}
-
+*/
double
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;
double piv_abs/*, first_elem*/;
NonZeroElem *first, *firsta, *first_suba;
- double piv_v[Size];
- int pivj_v[Size], pivk_v[Size], NR[Size], l, N_max;
+ double *piv_v;
+ int *pivj_v, *pivk_v, *NR;
+ int l, N_max;
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;
u_count_alloc_save = u_count_alloc;
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("The initial values of endogenous variables are too far from the solution.\n");
mexPrintf("Change them!\n");
- mexEvalString("drawnow;");
+ mexEvalString("drawnow;");
+ mxFree(piv_v);
+ mxFree(pivj_v);
+ mxFree(pivk_v);
+ mxFree(NR);
if(steady_state)
return false;
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);
mexEvalString("drawnow;");
- if(steady_state)
- return false;
- else
- {
+ mxFree(piv_v);
+ mxFree(pivj_v);
+ mxFree(pivk_v);
+ mxFree(NR);
+ if(steady_state)
+ return false;
+ else
+ {
mexEvalString("st=fclose('all');clear all;");
- filename += " stopped";
- mexErrMsgTxt(filename.c_str());
- }
+ filename += " stopped";
+ mexErrMsgTxt(filename.c_str());
+ }
}
slowc_save /= 2;
mexPrintf("Error: Simulation diverging, trying to correct it using slowc=%f\n", slowc_save);
for (i = 0; i < y_size; i++)
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++)
- 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--;
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("-----------------------------------\n");
}*/
- if (cvg)
- return (true);
- Simple_Init(it_, y_kmin, y_kmax, Size, IM_i);
+ if (cvg)
+ {
+ 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++)
{
/*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)
{
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)
return false;
- else
- {
+ else
+ {
mexEvalString("st=fclose('all');clear all;");
filename += " stopped";
mexErrMsgTxt(filename.c_str());
- }
+ }
}
/*divide all the non zeros elements of the line pivj by the max_pivot*/
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);
NonZeroElem *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++)
- {
+ {
if (!line_done[first->r_index])
bc[nb_eq_todo++] = first;
first = first->NZE_C_N;
- }
+ }
//#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
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;
first = first->NZE_C_N;
- }
- }
+ }
+ }
double slowc_lbx = slowc, res1bx;
for (i = 0; i < y_size; i++)
ya[i+it_*y_size] = y[i+it_*y_size];
slowc_save = slowc;
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;
}
@@ -1419,7 +1491,8 @@ SparseMatrix::CheckIt(int y_size, int y_kmin, int y_kmax, int Size, int periods,
mexPrintf("G1a red done\n");
SaveResult >> 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++)
SaveResult >> B[i];
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)
mexPrintf("Problem at i=%d u[b[i]]=%f B[i]=%f\n", i, u[b[i]], B[i]);
- }
+ }
+ mxFree(B);
}
void
@@ -1581,16 +1655,22 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
}
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++)
{
if (record && symbolic)
{
- if (save_op) ;
- {
- mxFree(save_op);
- save_op = NULL;
- }
+ if (save_op)
+ {
+ mxFree(save_op);
+ save_op = NULL;
+ }
save_op = (int *) mxMalloc(nop*sizeof(int));
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);
if ((symbolic && t <= start_compare) || !symbolic)
{
- double piv_v[Size];
- int pivj_v[Size], pivk_v[Size], NR[Size], l = 0, N_max = 0;
+ int l = 0, N_max = 0;
bool one = false;
piv_abs = 0;
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++)
{
- 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)
{
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*/
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++)
{
bb[j] = first;
+ /*if(nb_var==122)
+ mexPrintf("j=%d first=%x\n",j,first);*/
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;
}
}
- }
+ }
+ mxFree(bb);
nop += nb_var*2;
u[b[pivj]] /= piv;
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;
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;
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;
}
- }
+ }
+ mxFree(bc);
}
if (symbolic)
{
@@ -2006,6 +2093,10 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
}
//record = true;
}
+ mxFree(piv_v);
+ mxFree(pivj_v);
+ mxFree(pivk_v);
+ mxFree(NR);
}
nop_all += nop;
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)
mxFree(save_opaa);
}
- close_swp_file();
+ //close_swp_file();
/*The backward substitution*/
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;");
}
- close_swp_file();
+ //close_swp_file();
time00 = clock();
if (tbreak_g == 0)
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);
mxFree(ua);*/
-
return (0);
}
diff --git a/mex/sources/simulate/SparseMatrix.hh b/mex/sources/simulate/SparseMatrix.hh
index 3954fae7a..28eb9da07 100644
--- a/mex/sources/simulate/SparseMatrix.hh
+++ b/mex/sources/simulate/SparseMatrix.hh
@@ -30,11 +30,17 @@
// Nothing to single thread version.
#else
#include <omp.h>
+#endif
+#ifdef _MSC_VER
+ #include <limits>
#endif
#define NEW_ALLOC
#define MARKOVITZ
using namespace std;
+
+
+
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);
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 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);
+
+#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:
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
);
double simple_bksub(int it_, int Size, double slowc_l);
- void close_swp_file();
+ /*void close_swp_file();*/
stack<double> Stack;
int nb_prologue_table_u, nb_first_table_u, nb_middle_table_u, nb_last_table_u;
int nb_prologue_table_y, nb_first_table_y, nb_middle_table_y, nb_last_table_y;
@@ -148,6 +191,4 @@ protected:
-
-
#endif
diff --git a/mex/sources/simulate/simulate.cc b/mex/sources/simulate/simulate.cc
index f311df18a..e40a49c8d 100644
--- a/mex/sources/simulate/simulate.cc
+++ b/mex/sources/simulate/simulate.cc
@@ -24,7 +24,6 @@
////////////////////////////////////////////////////////////////////////
#include <cstring>
-
#include "simulate.hh"
#include "Interpreter.hh"
#include "Mem_Mngr.hh"
--
GitLab