Commit aa23ed73 authored by sebastien's avatar sebastien
Browse files

preprocessor + bytecode DLL: various enhancements to block and bytecode options (changes by Ferhat)


git-svn-id: https://www.dynare.org/svn/dynare/trunk@3244 ac1d8469-bf42-47a9-8791-bf33cf982152
parent 9f57f73d
......@@ -25,7 +25,7 @@ function model_info;
nb_blocks=length(M_.block_structure.block);
fprintf('The model has %d equations and is decomposed in %d blocks as follow:\n',M_.endo_nbr,nb_blocks);
fprintf('===============================================================================================================\n');
fprintf('| %10s | %10s | %30s | %14s | %31s |\n','Block n','Size','Block Type','E quation','Dependent variable');
fprintf('| %10s | %10s | %30s | %14s | %31s |\n','Block no','Size','Block Type',' Equation','Dependent variable');
fprintf('|============|============|================================|================|=================================|\n');
for i=1:nb_blocks
size_block=length(M_.block_structure.block(i).equation);
......@@ -44,7 +44,7 @@ function model_info;
fprintf('\n');
for k=1:M_.maximum_endo_lag+M_.maximum_endo_lead+1
if(k==M_.maximum_endo_lag+1)
fprintf('%-30s %s','the variable','is used in equations contemporously');
fprintf('%-30s %s','the variable','is used in equations ontemporaneously');
elseif(k<M_.maximum_endo_lag+1)
fprintf('%-30s %s %d','the variable','is used in equations with lag ',M_.maximum_endo_lag+1-k);
else
......
......@@ -110,7 +110,7 @@ function steady_()
ss(M_.blocksMFS{b}), ...
options_.jacobian_flag, b);
if check ~= 0
error('STEADY: convergence problems')
error(['STEADY: convergence problems in block ' int2str(b)])
end
ss(M_.blocksMFS{b}) = y;
end
......
......@@ -11,7 +11,6 @@ nodist_bytecode_SOURCES = \
Interpreter.cc \
Mem_Mngr.cc \
SparseMatrix.cc \
bytecode.hh \
Interpreter.hh \
Mem_Mngr.hh \
SparseMatrix.hh
......@@ -21,8 +21,6 @@
#include "Interpreter.hh"
#define BIG 1.0e+8;
#define SMALL 1.0e-5;
//#define DEBUG
Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, double *x_arg, double *steady_y_arg, double *steady_x_arg,
double *direction_arg, int y_size_arg,
......@@ -59,30 +57,12 @@ Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, doub
double
Interpreter::pow1(double a, double b)
{
/*double r;
if(a>=0)
r=pow_(a,b);
else
{
//r=0;
//max_res=res1=res2=BIG;
if(error_not_printed)
{
mexPrintf("Error: X^a with X<0\n");
error_not_printed = false;
}
//r = BIG;
//r = -pow_(-a, b);
//r = 0;
//r = SMALL;
//r = pow_(-a, b);
}*/
double r = pow_(a, b);
if (isnan(r) || isinf(r))
{
if(a<0 && error_not_printed)
{
mexPrintf("Error: X^a with X=%5.25f\n",a);
mexPrintf("Error: X^a with X=%5.15f\n",a);
error_not_printed = false;
r = 0.0000000000000000000000001;
}
......@@ -96,30 +76,12 @@ Interpreter::pow1(double a, double b)
double
Interpreter::log1(double a)
{
/*double r;
if(a>=0)
r=pow_(a,b);
else
{
//r=0;
//max_res=res1=res2=BIG;
if(error_not_printed)
{
mexPrintf("Error: X^a with X<0\n");
error_not_printed = false;
}
//r = BIG;
//r = -pow_(-a, b);
//r = 0;
//r = SMALL;
//r = pow_(-a, b);
}*/
double r = log(a);
if (isnan(r) || isinf(r))
{
if(a<=0 && error_not_printed)
{
mexPrintf("Error: log(X) with X<=0\n");
mexPrintf("Error: log(X) with X=%5.15f\n",a);
error_not_printed = false;
}
res1=NAN;
......@@ -129,10 +91,8 @@ Interpreter::log1(double a)
return r;
}
void
Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num) /*throw(EvalException)*/
Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
{
int var, lag = 0, op;
ostringstream tmp_out;
......@@ -161,9 +121,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num) /*thro
Stack.push(ya[(it_+lag)*y_size+var]);
else
{
/*mexPrintf(" y[%d, %d]=",(it_+lag)*y_size, var );
mexPrintf("%f\n",y[(it_+lag)*y_size+var]);*/
Stack.push(y[(it_+lag)*y_size+var]);
Stack.push(y[(it_+lag)*y_size+var]);
}
#ifdef DEBUG
tmp_out << " y[" << it_+lag << ", " << var << "](" << y[(it_+lag)*y_size+var] << ")";
......@@ -175,7 +133,6 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num) /*thro
Stack.push(x[it_+lag+var*nb_row_x]);
#ifdef DEBUG
tmp_out << " x[" << it_+lag << ", " << var << "](" << x[it_+lag+var*nb_row_x] << ")";
//mexPrintf(" x[%d, %d](%f)\n", it_+lag, var, x[it_+lag+var*nb_row_x]);
#endif
break;
case eExogenousDet :
......@@ -681,7 +638,6 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num) /*thro
#endif
break;
default:
/*throw EvalException();*/
;
}
break;
......@@ -703,14 +659,12 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num) /*thro
if (Stack.size()>0)
{
mexPrintf("error: Stack not empty!\n");
mexEvalString("st=fclose('all');clear all;");
mexErrMsgTxt("End of simulate");
mexErrMsgTxt("End of bytecode");
}
break;
default :
mexPrintf("Unknown opcode %d!! FENDEQU=%d\n",it_code->first,FENDEQU);
mexEvalString("st=fclose('all');clear all;");
mexErrMsgTxt("End of simulate");
mexErrMsgTxt("End of bytecode");
break;
}
it_code++;
......@@ -964,8 +918,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
if (!cvg)
{
mexPrintf("Convergence not achieved in block %d, at time %d after %d iterations\n",Block_Count,it_,iter);
mexEvalString("st=fclose('all');clear all;");
mexErrMsgTxt("End of simulate");
mexErrMsgTxt("End of bytecode");
}
}
}
......@@ -1025,8 +978,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
if (!cvg)
{
mexPrintf("Convergence not achieved in block %d, at time %d after %d iterations\n",Block_Count,it_,iter);
mexEvalString("st=fclose('all');clear all;");
mexErrMsgTxt("End of simulate");
mexErrMsgTxt("End of bytecode");
}
}
}
......@@ -1080,7 +1032,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
}
if (!cvg)
{
mexPrintf("Convergence not achieved in block %d, at time %d after %d iterations\n", Block_Count, it_, iter);
mexPrintf("Convergence not achieved in block %d, after %d iterations\n", Block_Count, iter);
return false;
}
}
......@@ -1143,8 +1095,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
if (!cvg)
{
mexPrintf("Convergence not achieved in block %d, at time %d after %d iterations\n", Block_Count, it_, iter);
mexEvalString("st=fclose('all');clear all;");
mexErrMsgTxt("End of simulate");
mexErrMsgTxt("End of bytecode");
}
}
}
......@@ -1216,7 +1167,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
}
if (!cvg)
{
mexPrintf("Convergence not achieved in block %d, at time %d after %d iterations\n", Block_Count, it_, iter);
mexPrintf("Convergence not achieved in block %d, after %d iterations\n", Block_Count, iter);
return false;
}
}
......@@ -1279,8 +1230,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
if (!cvg)
{
mexPrintf("Convergence not achieved in block %d, at time %d after %d iterations\n", Block_Count, it_, iter);
mexEvalString("st=fclose('all');clear all;");
mexErrMsgTxt("End of simulate");
mexErrMsgTxt("End of bytecode");
}
}
}
......@@ -1364,8 +1314,6 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
cvg = false;
else
cvg=(max_res<solve_tolf);
/*if(cvg)
continue;*/
u_count=u_count_saved;
simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter, minimal_solving_periods);
iter++;
......@@ -1373,8 +1321,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
if (!cvg)
{
mexPrintf("Convergence not achieved in block %d, after %d iterations\n",Block_Count, iter);
mexEvalString("st=fclose('all');clear all;");
mexErrMsgTxt("End of simulate");
mexErrMsgTxt("End of bytecode");
}
}
else
......@@ -1411,9 +1358,8 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
break;
default:
mexPrintf("Unknown type =%d\n",type);
mexEvalString("st=fclose('all');clear all;");
mexEvalString("drawnow;");
mexErrMsgTxt("End of simulate");
mexErrMsgTxt("End of bytecode");
}
return true;
}
......@@ -1436,7 +1382,6 @@ Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_s
{
mexPrintf("%s.cod Cannot be opened\n",file_name.c_str());
mexEvalString("drawnow;");
mexEvalString("st=fclose('all');clear all;");
filename+=" stopped";
mexEvalString("drawnow;");
mexErrMsgTxt(filename.c_str());
......@@ -1500,9 +1445,8 @@ Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_s
break;
default :
mexPrintf("Unknown command \n");
mexEvalString("st=fclose('all');clear all;");
mexEvalString("drawnow;");
mexErrMsgTxt("End of simulate");
mexErrMsgTxt("End of bytecode");
break;
}
}
......
......@@ -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 "Mem_Mngr.hh"
Mem_Mngr::Mem_Mngr()
......@@ -74,7 +74,6 @@ Mem_Mngr::mxMalloc_NZE()
}
else /*We have to allocate extra memory space*/
{
//mexPrintf("CHUNK_SIZE=%d CHUNK_BLCK_SIZE=%d Nb_CHUNK=%d\n",CHUNK_SIZE,CHUNK_BLCK_SIZE,Nb_CHUNK);
CHUNK_SIZE+=CHUNK_BLCK_SIZE;
Nb_CHUNK++;
NZE_Mem=(NonZeroElem*)mxMalloc(CHUNK_BLCK_SIZE*sizeof(NonZeroElem)); /*The block of memory allocated*/
......
......@@ -20,7 +20,7 @@
#ifndef MEM_MNGR_HH_INCLUDED
#define MEM_MNGR_HH_INCLUDED
//#include <stack>
#include <vector>
#include <fstream>
#ifndef DEBUG_EX
......@@ -34,7 +34,7 @@ struct NonZeroElem
{
int u_index;
int r_index, c_index, lag_index;
NonZeroElem *NZE_R_N, *NZE_C_N/*, *NZE_C_P*/;
NonZeroElem *NZE_R_N, *NZE_C_N;
};
typedef vector<NonZeroElem*> v_NonZeroElem;
......@@ -42,9 +42,6 @@ typedef vector<NonZeroElem*> v_NonZeroElem;
class Mem_Mngr
{
public:
/*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 Print_heap();
void init_Mem();
void mxFree_NZE(void* pos);
......@@ -53,15 +50,11 @@ public:
void Free_All();
Mem_Mngr();
void fixe_file_name(string filename_arg);
/*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);*/
bool swp_f;
//bool verbose;
private:
v_NonZeroElem Chunk_Stack;
int CHUNK_SIZE, CHUNK_BLCK_SIZE, Nb_CHUNK;
int CHUNK_heap_pos/*, CHUNK_heap_max_size*/;
int CHUNK_heap_pos;
NonZeroElem** NZE_Mem_add;
NonZeroElem* NZE_Mem;
vector<NonZeroElem*> NZE_Mem_Allocated;
......
......@@ -143,24 +143,6 @@ SparseMatrix::Delete(const int r, const int c)
{
NonZeroElem *first = FNZE_R[r], *firsta = NULL;
//mexPrintf("Delete r=%d, c=%d\n", r, c);
/*map<pair<int, int>,NonZeroElem*>::const_iterator it;
it = Mapped_Array.find(make_pair(r, c));
if(it==Mapped_Array.end())
mexPrintf("Not Found\n");
first = it->second;
if(it != Mapped_Array.begin())
{
it--;
if(it->first.first == r)
firsta = it->second;
else
firsta = NULL;
}
else
firsta = NULL;
Mapped_Array.erase(make_pair(r, c));*/
while (first->c_index != c)
{
firsta = first;
......@@ -179,14 +161,6 @@ SparseMatrix::Delete(const int r, const int c)
firsta = first;
first = first->NZE_C_N;
}
//firsta = first->NZE_C_P;
/*if(first->NZE_C_N != NULL)
{
//mexPrintf("ultime\n");
(first->NZE_C_N)->NZE_C_P = firsta;
}*/
//mexPrintf("ok1\n");
if (firsta != NULL)
firsta->NZE_C_N = first->NZE_C_N;
......@@ -196,24 +170,6 @@ SparseMatrix::Delete(const int r, const int c)
u_liste.push_back(first->u_index);
mem_mngr.mxFree_NZE(first);
NbNZCol[c]--;
/*Check the deletition*/
/*int nb_var=NbNZRow[r];
first=FNZE_R[r];
for(int j=0;j<nb_var;j++)
{
if(!first)
mexPrintf("Error in Delete (Row) r=%d and c=%d \n",r,c);
first=first->NZE_R_N;
}
nb_var=NbNZCol[c];
first=FNZE_C[c];
for(int j=0;j<nb_var;j++)
{
if(!first)
mexPrintf("Error in Delete (Col) r=%d and c=%d \n",r,c);
first=first->NZE_C_N;
}*/
//mexPrintf("done\n");
}
void
......@@ -288,7 +244,7 @@ SparseMatrix::Insert(const int r, const int c, const int u_index, const int lag_
firsta->NZE_R_N = firstn;
firstn->NZE_R_N = first;
}
else /*first.c_index<c*/
else
{
first->NZE_R_N = firstn;
firstn->NZE_R_N = NULL;
......@@ -309,42 +265,20 @@ SparseMatrix::Insert(const int r, const int c, const int u_index, const int lag_
firsta->NZE_C_N = firstn;
firstn->NZE_C_N = first;
}
else /*first.r_index<r*/
else
{
first->NZE_C_N = firstn;
firstn->NZE_C_N = NULL;
}
//if (firsta != NULL)
/*firstn->NZE_C_P = firsta;
if(first != NULL)
first->NZE_C_P = firstn;*/
NbNZCol[c]++;
//Mapped_Array[make_pair(r, c)] = firstn;
/*Check the insertion*/
/*int nb_var=NbNZRow[r];
first=FNZE_R[r];
for(int j=0;j<nb_var;j++)
{
if(!first)
mexPrintf("Error in insert (Row) r=%d and c=%d \n",r,c);
first=first->NZE_R_N;
}
nb_var=NbNZCol[c];
first=FNZE_C[c];
for(int j=0;j<nb_var;j++)
{
if(!first)
mexPrintf("Error in insert (Col) r=%d and c=%d \n",r,c);
first=first->NZE_C_N;
}*/
}
void
SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, int y_kmin, int y_kmax, bool steady_state, bool two_boundaries)
{
int i, j, eq, var, lag;
unsigned int eq, var;
int i, j, lag;
filename = file_name;
mem_mngr.fixe_file_name(file_name);
if (!SaveCode.is_open())
......@@ -370,28 +304,12 @@ SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, 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++;
}
IM_i[make_pair(make_pair(j, Size*(periods+y_kmax)), 0)] = j;
}
else
{
......@@ -399,7 +317,6 @@ SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, 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;
......@@ -478,7 +395,7 @@ SparseMatrix::Simple_Init(int it_, int y_kmin, int y_kmax, int Size, map<pair<pa
first = mem_mngr.mxMalloc_NZE();
first->NZE_C_N = NULL;
first->NZE_R_N = NULL;
first->u_index = u_count1 /*it4->second+u_count_init*it_*/;
first->u_index = u_count1;
first->r_index = eq;
first->c_index = var;
first->lag_index = lag;
......@@ -492,7 +409,6 @@ SparseMatrix::Simple_Init(int it_, int y_kmin, int y_kmax, int Size, map<pair<pa
temp_NZE_C[var]->NZE_C_N = first;
temp_NZE_R[eq] = first;
temp_NZE_C[var] = first;
//Mapped_Array[make_pair(eq,var)]=first;
u_count1++;
}
it4++;
......@@ -502,7 +418,7 @@ SparseMatrix::Simple_Init(int it_, int y_kmin, int y_kmax, int Size, map<pair<pa
b[i] = i;
mxFree(temp_NZE_R);
mxFree(temp_NZE_C);
u_count = u_count1/*+Size*/;
u_count = u_count1;
}
void
......@@ -572,7 +488,6 @@ SparseMatrix::Init(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<
NbNZRow[eq]++;
NbNZCol[var]++;
first = mem_mngr.mxMalloc_NZE();
//first->NZE_C_P = temp_NZE_C[var];
first->NZE_C_N = NULL;
first->NZE_R_N = NULL;
first->u_index = it4->second+u_count_init*t;
......@@ -589,7 +504,6 @@ SparseMatrix::Init(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<
temp_NZE_C[var]->NZE_C_N = first;
temp_NZE_R[eq] = first;
temp_NZE_C[var] = first;
//Mapped_Array[make_pair(eq,var)]=first;
}
else /*Build the additive terms ooutside the simulation periods related to the first lags and the last leads...*/
{
......@@ -777,7 +691,6 @@ SparseMatrix::compare(int *save_op, int *save_opa, int *save_opaa, int beg_t, in
}
j++;
}
//mexPrintf("beg_t=%d, periods=%d, j=%d OK=%d y_kmin=%d\n",beg_t, periods, j, OK, y_kmin);
// the same pivot for all remaining periods
if (OK)
{
......@@ -1021,12 +934,6 @@ SparseMatrix::complete(int beg_t, int Size, int periods, int *b)
return (beg_t);
}
/*void
SparseMatrix::close_swp_file()
{
mem_mngr.close_swp_f();
}
*/
double
SparseMatrix::bksub(int tbreak, int last_period, int Size, double slowc_l)
{
......@@ -1108,7 +1015,7 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
{
int i, j, k;
int pivj = 0, pivk = 0;
double piv_abs/*, first_elem*/;
double piv_abs;
NonZeroElem *first, *firsta, *first_suba;
double *piv_v;
int *pivj_v, *pivk_v, *NR;
......@@ -1191,8 +1098,6 @@ 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);
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];*/
mxFree(piv_v);
mxFree(pivj_v);
mxFree(pivk_v);
......@@ -1200,15 +1105,6 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
iter--;
return true;
}
/*if (Size>1)
{
mexPrintf("-----------------------------------\n");
mexPrintf(" Simulate iteration %d \n", iter+1);
mexPrintf(" max. error=%.10e \n", double (max_res));
mexPrintf(" sqr. error=%.10e \n", double (res2));
mexPrintf(" abs. error=%.10e \n", double (res1));
mexPrintf("-----------------------------------\n");
}*/
if (cvg)
{
mxFree(piv_v);
......@@ -1341,8 +1237,6 @@ 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")))
for (j = 0; j < nb_eq_todo; j++)
{
//t_save_op_s *save_op_s_l;
//NonZeroElem *
first = bc[j];
int row = first->r_index;
double first_elem = u[first->u_index];
......@@ -1376,7 +1270,7 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
piv_c_index = Size;
l_piv++;
}
else //first_sub->c_index==first_piv->c_index
else
{
if (i == sub_c_index)
{
......@@ -1640,7 +1534,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
if (print_it)
{
mexPrintf("-----------------------------------\n");
mexPrintf(" Simulate iteration %d \n", iter+1);
mexPrintf(" Simulate iteration no %d \n", iter+1);
mexPrintf(" max. error=%.10e \n", double (max_res));
mexPrintf(" sqr. error=%.10e \n", double (res2));
mexPrintf(" abs. error=%.10e \n", double (res1));
......@@ -1652,12 +1546,12 @@ 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);
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));
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++)
{
......@@ -1787,20 +1681,16 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
{
mexPrintf("Error: singular system in Simulate_NG1\n");
mexEvalString("d