Commit d32dd992 authored by ferhat's avatar ferhat
Browse files

- extension of normalization of equations to nonlinear equations

- mfs: new option for 'steady' and 'model' commands. Determines the equation belonging to the set of feedback variables.
  mfs = 0 => all variables are considered as feedback variables (default value)
  mfs = 1 => using only naturally normalized equation as potential recursive equations (all variables assigned to unnormalized equations are considered as feedback variable)
  mfs = 2 => adding to the set of potential recursive equation with mfs = 1 the linear equation in endogenous variable normalized (all variables assigned to nonlinear unnormalized equations are considered as feedback variable)
  mfs = 3 => adding to the set of potential recursive equation with mfs = 2 the non linear equation in endogenous variable normalized
- correction of few buggs in simulate.dll
- block_mfs_dll: new option for 'steady' command. Use simulate.dll to solve the steady state model (speedup the computation of the steady-state and the homotopy)

git-svn-id: https://www.dynare.org/svn/dynare/trunk@2866 ac1d8469-bf42-47a9-8791-bf33cf982152
parent 5d1fdab1
......@@ -225,3 +225,6 @@ function global_initialization()
% block decomposition + minimum feedback set for steady state computation
options_.block_mfs = 0;
% block decomposition + minimum feedback set for steady state computation
% using simulate.dll
options_.block_mfs_dll = 0;
......@@ -24,23 +24,23 @@ function model_info;
if(isfield(M_,'block_structure'))
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 | %30s |\n','Block n','Size','Block Type','Equation','Dependent variable');
fprintf('|============|============|================================|================|================================|\n');
fprintf('===============================================================================================================\n');
fprintf('| %10s | %10s | %30s | %14s | %31s |\n','Block n','Size','Block Type','E quation','Dependent variable');
fprintf('|============|============|================================|================|=================================|\n');
for i=1:nb_blocks
size_block=length(M_.block_structure.block(i).equation);
if(i>1)
fprintf('|------------|------------|--------------------------------|----------------|--------------------------------|\n');
fprintf('|------------|------------|--------------------------------|----------------|---------------------------------|\n');
end;
for j=1:size_block
if(j==1)
fprintf('| %3d (%4d) | %10d | %30s | %14d | %30s |\n',i,M_.block_structure.block(i).num,size_block,Sym_type(M_.block_structure.block(i).Simulation_Type),M_.block_structure.block(i).equation(j),M_.endo_names(M_.block_structure.block(i).variable(j),:));
fprintf('| %3d (%4d) | %10d | %30s | %14d | %-6d %24s |\n',i,M_.block_structure.block(i).num,size_block,Sym_type(M_.block_structure.block(i).Simulation_Type),M_.block_structure.block(i).equation(j),M_.block_structure.block(i).variable(j),M_.endo_names(M_.block_structure.block(i).variable(j),:));
else
fprintf('| %10s | %10s | %30s | %14d | %30s |\n','','','',M_.block_structure.block(i).equation(j),M_.endo_names(M_.block_structure.block(i).variable(j),:));
fprintf('| %10s | %10s | %30s | %14d | %-6d %24s |\n','','','',M_.block_structure.block(i).equation(j),M_.block_structure.block(i).variable(j),M_.endo_names(M_.block_structure.block(i).variable(j),:));
end;
end;
end;
fprintf('==============================================================================================================\n');
fprintf('===============================================================================================================\n');
fprintf('\n');
for k=1:M_.maximum_endo_lag+M_.maximum_endo_lead+1
if(k==M_.maximum_endo_lag+1)
......
......@@ -53,12 +53,12 @@ function resid(period)
addpath(M_.fname);
end;
for it_=M_.maximum_lag+1:period+M_.maximum_lag
if(options_.model_mode == 1 || options_.model_mode == 3)
z(:,it_-M_.maximum_lag) = feval(fh,oo_.endo_simul',oo_.exo_simul, M_.params, it_);
else
%if(options_.model_mode == 1 || options_.model_mode == 3)
% z(:,it_-M_.maximum_lag) = feval(fh,oo_.endo_simul',oo_.exo_simul, M_.params, it_);
%else
z(:,it_-M_.maximum_lag) = feval(fh,y(iyr0),oo_.exo_simul, M_.params, it_);
iyr0 = iyr0 + n;
end;
%end;
end
if(options_.model_mode == 1 || options_.model_mode == 3)
rmpath(M_.fname);
......
......@@ -100,8 +100,42 @@ function [y, info] = solve_one_boundary(fname, y, x, params, y_index_eq, nze, pe
max_res=max(max(abs(r)));
end;
%['max_res=' num2str(max_res) ' Block_Num=' int2str(Block_Num) ' it_=' int2str(it_)]
disp(['iteration : ' int2str(iter+1) ' => ' num2str(max_res) ' time = ' int2str(it_)]);
% fjac = zeros(Blck_size, Blck_size);
% disp(['Blck_size=' int2str(Blck_size) ' it_=' int2str(it_)]);
% dh = max(abs(y(it_, y_index_eq)),options_.gstep*ones(1, Blck_size))*eps^(1/3);
% fvec = r;
% for j = 1:Blck_size
% ydh = y ;
% ydh(it_, y_index_eq(j)) = ydh(it_, y_index_eq(j)) + dh(j) ;
% [t, y1, g11, g21, g31]=feval(fname, ydh, x, params, it_, 0);
% fjac(:,j) = (t - fvec)./dh(j) ;
% end;
% diff = g1 -fjac;
% diff
% disp('g1');
% disp([num2str(g1,'%4.5f')]);
% disp('fjac');
% disp([num2str(fjac,'%4.5f')]);
% [c_max, i_c_max] = max(abs(diff));
% [l_c_max, i_r_max] = max(c_max);
% disp(['maximum element row=' int2str(i_c_max(i_r_max)) ' and column=' int2str(i_r_max) ' value = ' num2str(l_c_max)]);
% equation = i_c_max(i_r_max);
% variable = i_r_max;
% variable
% mod(variable, Blck_size)
% disp(['equation ' int2str(equation) ' and variable ' int2str(y_index_eq(variable)) ' ' M_.endo_names(y_index_eq(variable), :)]);
% disp(['g1(' int2str(equation) ', ' int2str(variable) ')=' num2str(g1(equation, variable),'%3.10f') ' fjac(' int2str(equation) ', ' int2str(variable) ')=' num2str(fjac(equation, variable), '%3.10f') ' y(' int2str(it_) ', ' int2str(variable) ')=' num2str(y(it_, variable))]);
% %return;
% %g1 = fjac;
if(verbose==1)
disp(['iteration : ' int2str(iter) ' => ' num2str(max_res) ' time = ' int2str(it_)]);
disp(['iteration : ' int2str(iter+1) ' => ' num2str(max_res) ' time = ' int2str(it_)]);
if(is_dynamic)
disp([M_.endo_names(y_index_eq,:) num2str([y(it_,y_index_eq)' r g1])]);
else
......@@ -118,7 +152,7 @@ function [y, info] = solve_one_boundary(fname, y, x, params, y_index_eq, nze, pe
if(~cvg)
if(iter>0)
if(~isreal(max_res) | isnan(max_res) | (max_resa<max_res && iter>1))
if(isnan(max_res))
if(isnan(max_res)| (max_resa<max_res && iter>0))
detJ=det(g1a);
if(abs(detJ)<1e-7)
max_factor=max(max(abs(g1a)));
......@@ -198,15 +232,23 @@ function [y, info] = solve_one_boundary(fname, y, x, params, y_index_eq, nze, pe
else
info = -Block_Num*10;
end
elseif(~is_dynamic & options_.solve_algo==2)
elseif((~is_dynamic & options_.solve_algo==2) || (is_dynamic & options_.solve_algo==4))
lambda=1;
stpmx = 100 ;
if (is_dynamic)
stpmax = stpmx*max([sqrt(y*y');size(y_index_eq,2)]);
else
stpmax = stpmx*max([sqrt(y'*y);size(y_index_eq,2)]);
end;
nn=1:size(y_index_eq,2);
g = (r'*g1)';
f = 0.5*r'*r;
p = -g1\r ;
if (is_dynamic)
[y,f,r,check]=lnsrch1(y,f,g,p,stpmax,fname,nn,y_index_eq,x, params, it_, 0);
else
[y,f,r,check]=lnsrch1(y,f,g,p,stpmax,fname,nn,y_index_eq,x, params, 0);
end;
dx = ya - y(y_index_eq);
elseif(~is_dynamic & options_.solve_algo==3)
[yn,info] = csolve(@local_fname, y(y_index_eq),@local_fname,1e-6,500, x, params, y, y_index_eq, fname, 1);
......@@ -283,6 +325,7 @@ function [y, info] = solve_one_boundary(fname, y, x, params, y_index_eq, nze, pe
return;
end;
iter=iter+1;
max_resa = max_res;
end
end
if cvg==0
......
......@@ -58,7 +58,7 @@ function y = solve_two_boundaries(fname, y, x, params, y_index, nze, periods, y_
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global oo_ M_ T9025 T1149 T11905;
global options_ oo_ M_ T9025 T1149 T11905;
cvg=0;
iter=0;
Per_u_=0;
......@@ -73,11 +73,48 @@ function y = solve_two_boundaries(fname, y, x, params, y_index, nze, periods, y_
reduced = 0;
while ~(cvg==1 | iter>maxit_),
[r, y, g1, g2, g3, b]=feval(fname, y, x, params, periods, 0, y_kmin, Blck_size);
% fjac = zeros(Blck_size, Blck_size*(y_kmin_l+1+y_kmax_l));
% disp(['Blck_size=' int2str(Blck_size) ' size(y_index)=' int2str(size(y_index,2))]);
% dh = max(abs(y(y_kmin+1-y_kmin_l:y_kmin+1+y_kmax_l, y_index)),options_.gstep*ones(y_kmin_l+1+y_kmax_l, Blck_size))*eps^(1/3);
% fvec = r;
% %for i = y_kmin+1-y_kmin_l:y_kmin+1+y_kmax_l
% i = y_kmin+1;
% i
% for j = 1:Blck_size
% ydh = y ;
% ydh(i, y_index(j)) = ydh(i, y_index(j)) + dh(i, j) ;
% if(j==11 && i==2)
% disp(['y(i,y_index(11)=' int2str(y_index(11)) ')= ' num2str(y(i,y_index(11))) ' ydh(i, y_index(j))=' num2str(ydh(i, y_index(j))) ' dh(i,j)= ' num2str(dh(i,j))]);
% end;
% [t, y1, g11, g21, g31, b1]=feval(fname, ydh, x, params, periods, 0, y_kmin, Blck_size);
% fjac(:,j+(i-(y_kmin+1-y_kmin_l))*Blck_size) = (t(:, 1+y_kmin) - fvec(:, 1+y_kmin))./dh(i, j) ;
% if(j==11 && i==2)
% disp(['fjac(:,' int2str(j+(i-(y_kmin+1-y_kmin_l))*Blck_size) ')=']);
% disp([num2str(fjac(:,j+(i-(y_kmin+1-y_kmin_l))*Blck_size))]);
% end;
% end;
% % end
% %diff = g1(1:Blck_size, 1:Blck_size*(y_kmin_l+1+y_kmax_l)) -fjac;
% diff = g1(1:Blck_size, y_kmin_l*Blck_size+1:(y_kmin_l+1)*Blck_size) -fjac(1:Blck_size, y_kmin_l*Blck_size+1:(y_kmin_l+1)*Blck_size);
% disp(diff);
% [c_max, i_c_max] = max(abs(diff));
% [l_c_max, i_r_max] = max(c_max);
% disp(['maximum element row=' int2str(i_c_max(i_r_max)) ' and column=' int2str(i_r_max) ' value = ' num2str(l_c_max)]);
% equation = i_c_max(i_r_max);
% variable = i_r_max;
% variable
% disp(['equation ' int2str(equation) ' and variable ' int2str(y_index(mod(variable, Blck_size))) ' ' M_.endo_names(y_index(mod(variable, Blck_size)), :)]);
% disp(['g1(' int2str(equation) ', ' int2str(variable) ')=' num2str(g1(equation, y_kmin_l*Blck_size+variable),'%3.10f') ' fjac(' int2str(equation) ', ' int2str(variable) ')=' num2str(fjac(equation, y_kmin_l*Blck_size+variable), '%3.10f')]);
% return;
% for i=1:periods;
% disp([sprintf('%5.14f ',[T9025(i) T1149(i) T11905(i)])]);
% end;
% return;
residual = r(:,y_kmin+1:y_kmin+1+y_kmax_l);
%residual = r(:,y_kmin+1:y_kmin+1+y_kmax_l);
%num2str(residual,' %1.6f')
%jac_ = g1(1:(y_kmin)*Blck_size, 1:(y_kmin+1+y_kmax_l)*Blck_size);
%jac_
......
......@@ -79,6 +79,8 @@ function steady_()
[oo_.exo_steady_state; ...
oo_.exo_det_steady_state], M_.params);
end
elseif options_.block_mfs_dll
[oo_.steady_state,check] = simulate('steady_state');
else
[oo_.steady_state,check] = dynare_solve([M_.fname '_static'],...
oo_.steady_state,...
......
This diff is collapsed.
......@@ -67,8 +67,9 @@ class Interpreter : SparseMatrix
{
protected :
double pow1(double a, double b);
double log1(double a);
void compute_block_time(int Per_u_);
void simulate_a_block(int size,int type, string file_name, string bin_basename, bool Gaussian_Elimination);
bool simulate_a_block(int size,int type, string file_name, string bin_basename, bool Gaussian_Elimination, bool steady_state, int block_num);
double *T;
vector<Block_contain_type> Block_Contain;
vector<Block_type> Block;
......@@ -83,6 +84,7 @@ class Interpreter : SparseMatrix
double *x, *params;
//double *y, *ya, *x, *direction;
map<pair<pair<int, int> ,int>, int> IM_i;
int equation, derivative_equation, derivative_variable;
string filename;
public :
//ReadBinFile read_bin_file;
......@@ -90,7 +92,7 @@ class Interpreter : SparseMatrix
Interpreter(double *params_arg, double *y_arg, double *ya_arg, double *x_arg, double *direction_arg, int y_size_arg, int nb_row_x_arg,
int nb_row_xd_arg, int periods_arg, int y_kmin_arg, int y_kmax_arg, int maxit_arg_, double solve_tolf_arg, int size_o_direction_arg,
double slowc_arg, int y_decal_arg, double markowitz_c_arg, string &filename_arg);
void compute_blocks(string file_name, string bin_basename);
bool compute_blocks(string file_name, string bin_basename, bool steady_state);
};
......
......@@ -34,7 +34,7 @@ struct NonZeroElem
{
int u_index;
int r_index, c_index, lag_index;
NonZeroElem *NZE_R_N, *NZE_C_N;
NonZeroElem *NZE_R_N, *NZE_C_N/*, *NZE_C_P*/;
};
typedef vector<NonZeroElem*> v_NonZeroElem;
......
This diff is collapsed.
......@@ -33,15 +33,13 @@
#endif
#define NEW_ALLOC
#define MARKOVITZ
//#define PROFILER
//#define MEMORY_LEAKS
using namespace std;
struct t_save_op_s
{
short int lag, operat;
long int first, second;
int first, second;
};
const int IFLD =0;
......@@ -62,11 +60,10 @@ class SparseMatrix
public:
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 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 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 initialize(int periods_arg, int nb_endo_arg, int y_kmin_arg, int y_kmax_arg, int y_size_arg, int u_count_arg, int u_count_init_arg, double *u_arg, double *y_arg, double *ya_arg, double slowc_arg, int y_decal_arg, double markowitz_c_arg, double res1_arg, double res2_arg, double max_res_arg);
void Read_SparseMatrix(string file_name, int Size, int periods, int y_kmin, int y_kmax);
void Read_SparseMatrix(string file_name, int Size, int periods, int y_kmin, int y_kmax, bool steady_state);
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:
......@@ -80,7 +77,7 @@ class SparseMatrix
#endif
);
void Insert(const int r, const int c, const int u_index, const int lag_index);
void Delete(const int r,const int c, const int Size);
void Delete(const int r,const int c);
int At_Row(int r, NonZeroElem **first);
int At_Pos(int r, int c, NonZeroElem **first);
int At_Col(int c, NonZeroElem **first);
......@@ -102,9 +99,6 @@ class SparseMatrix
#endif
);
double simple_bksub(int it_, int Size, double slowc_l);
void run_triangular(int nop_all,int *op_all);
void run_it(int nop_all,int *op_all);
void run_u_period1(int periods);
void close_swp_file();
stack<double> Stack;
int nb_prologue_table_u, nb_first_table_u, nb_middle_table_u, nb_last_table_u;
......@@ -118,6 +112,7 @@ class SparseMatrix
Mem_Mngr mem_mngr;
vector<int> u_liste;
map<pair<int, int>,NonZeroElem*> Mapped_Array;
int *NbNZRow, *NbNZCol;
NonZeroElem **FNZE_R, **FNZE_C;
int nb_endo, u_count_init;
......@@ -148,6 +143,7 @@ protected:
double *direction;
int start_compare;
int restart;
bool error_not_printed;
};
......
......@@ -43,7 +43,6 @@ max(int a, int b)
#ifdef DEBUG_EX
/*The Matlab c++ interface*/
using namespace std;
#include <sstream>
......@@ -53,7 +52,14 @@ int
main( int argc, const char* argv[] )
{
FILE *fid;
bool steady_state = false;
printf("argc=%d\n",argc);
if(argc<2)
{
mexPrintf("model filename expected\n");
mexEvalString("st=fclose('all');clear all;");
mexErrMsgTxt("Exit from Dynare");
}
float f_tmp;
ostringstream tmp_out("");
tmp_out << argv[1] << "_options.txt";
......@@ -64,51 +70,48 @@ main( int argc, const char* argv[] )
double *direction;
string file_name(argv[1]);
//mexPrintf("file_name=%s\n",file_name.c_str());
if(argc>2)
{
string f(argv[1]);
if(f == "steady_state")
steady_state = true;
}
fid = fopen(tmp_out.str().c_str(),"r");
int periods;
fscanf(fid,"%d",&periods);
int maxit_;
fscanf(fid,"%d",&maxit_);
fscanf(fid,"%f",&f_tmp);
slowc = f_tmp;
//mexPrintf("slowc_save=%f\n",slowc_save);
double slowc = f_tmp;
fscanf(fid,"%f",&f_tmp);
markowitz_c = f_tmp;
double markowitz_c = f_tmp;
fscanf(fid,"%f",&f_tmp);
solve_tolf = f_tmp;
double solve_tolf = f_tmp;
fclose(fid);
tmp_out.str("");
tmp_out << argv[1] << "_M.txt";
//printf("%s\n",tmp_out.str().c_str());
fid = fopen(tmp_out.str().c_str(),"r");
int y_kmin;
fscanf(fid,"%d",&y_kmin);
//printf("y_kmin=%d\n",y_kmin);
int y_kmax;
fscanf(fid,"%d",&y_kmax);
//printf("y_kmax=%d\n",y_kmax);
int y_decal;
fscanf(fid,"%d",&y_decal);
//printf("y_decal=%d\n",y_decal);
fscanf(fid,"%d",&nb_params);
//printf("nb_params=%d\n",nb_params);
fscanf(fid,"%d",&row_x);
//printf("row_x=%d\n",row_x);
fscanf(fid,"%d",&col_x);
//printf("col_x=%d\n",col_x);
fscanf(fid,"%d",&row_y);
//printf("row_y=%d\n",row_y);
fscanf(fid,"%d",&col_y);
//printf("col_y=%d\n",col_y);
int nb_row_xd;
fscanf(fid,"%d",&nb_row_xd);
//printf("nb_row_xd=%d\n",nb_row_xd);
params = (double*)malloc(nb_params*sizeof(params[0]));
//printf("OK1\n");
double * params = (double*)malloc(nb_params*sizeof(params[0]));
for(i=0; i < nb_params; i++)
{
fscanf(fid,"%f",&f_tmp);
params[i] = f_tmp;
//printf("param[%d]=%f\n",i,params[i]);
}
//printf("OK2\n");
fclose(fid);
yd = (double*)malloc(row_y*col_y*sizeof(yd[0]));
xd = (double*)malloc(row_x*col_x*sizeof(xd[0]));
......@@ -127,12 +130,12 @@ main( int argc, const char* argv[] )
}
fclose(fid);
size_of_direction=col_y*row_y*sizeof(double);
y=(double*)mxMalloc(size_of_direction);
ya=(double*)mxMalloc(size_of_direction);
int size_of_direction=col_y*row_y*sizeof(double);
double * y=(double*)mxMalloc(size_of_direction);
double * ya=(double*)mxMalloc(size_of_direction);
direction=(double*)mxMalloc(size_of_direction);
memset(direction,0,size_of_direction);
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++)
x[i]=double(xd[i]);
for (i=0;i<row_y*col_y;i++)
......@@ -140,36 +143,14 @@ main( int argc, const char* argv[] )
free(yd);
free(xd);
y_size=row_y;
x_size=col_x/*row_x*/;
nb_row_x=row_x;
/*for(int i=0; i<y_size; i++)
{
for(int it_=0; it_<8;it_++)
mexPrintf("y[t=%d, var=%d]=%f ",it_+1, i+1, y[(it_)*y_size+i]);
mexPrintf("\n");
}
for(int i=0; i<col_x; i++)
{
for(int it_=0; it_<8;it_++)
mexPrintf("x[t=%d, var=%d]=%f ",it_, i+1, x[it_+i*nb_row_x]);
mexPrintf("\n");
}*/
t0= clock();
int y_size=row_y;
int nb_row_x=row_x;
clock_t t0= clock();
Interpreter interprete(params, y, ya, x, 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);
string f(file_name);
interprete.compute_blocks(f+"_dynamic", f);
t1= clock();
interprete.compute_blocks(f, f, steady_state);
clock_t t1= clock();
mexPrintf("Simulation Time=%f milliseconds\n",1000.0*(double(t1)-double(t0))/double(CLOCKS_PER_SEC));
/*if (nlhs>0)
{
plhs[0] = mxCreateDoubleMatrix(row_y, col_y, mxREAL);
pind = mxGetPr(plhs[0]);
for (i=0;i<row_y*col_y;i++)
pind[i]=y[i];
}*/
if(x)
mxFree(x);
if(y)
......@@ -187,9 +168,24 @@ void
mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
mxArray *M_, *oo_, *options_;
int i, row_y, col_y, row_x, col_x;
int i, row_y, col_y, row_x, col_x, nb_row_xd;
int y_kmin=0, y_kmax=0, y_decal=0, periods=1;
double * pind ;
double *direction;
bool steady_state = false;
if(nrhs>0)
{
const mxArray *mxa = prhs[0];
int buflen=mxGetM(mxa) * mxGetN(mxa) + 1;
char *first_argument;
first_argument=(char*)mxCalloc(buflen, sizeof(char));
int status = mxGetString(mxa, first_argument, buflen);
if (status != 0)
mexWarnMsgTxt("Not enough space. The first argument is truncated.");
string f(first_argument);
if(f == "steady_state")
steady_state = true;
}
M_ = mexGetVariable("global","M_");
if (M_ == NULL )
{
......@@ -213,22 +209,41 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
mexErrMsgTxt("options_ \n");
}
//mexPrintf("ok0\n");
params = mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,"params")));
double * params = mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,"params")));
double *yd, *xd;
if(!steady_state)
{
yd= mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,"endo_simul")));
row_y=mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,"endo_simul")));
col_y=mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,"endo_simul")));;
xd= mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,"exo_simul")));
row_x=mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,"exo_simul")));
col_x=mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,"exo_simul")));
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_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")))))));
periods=int(floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,"periods"))))));
maxit_=int(floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,"maxit_"))))));
slowc=double(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,"slowc")))));
//slowc_save=slowc;
markowitz_c=double(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,"markowitz")))));
}
else
{
yd= mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,"steady_state")));
row_y=mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,"steady_state")));
col_y=mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,"steady_state")));;
xd= mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,"exo_steady_state")));
row_x=mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,"exo_steady_state")));
col_x=mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,"exo_steady_state")));
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_"))))));
double slowc=double(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,"slowc")))));
double markowitz_c=double(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,"markowitz")))));
double solve_tolf;
if(steady_state)
solve_tolf=*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,"solve_tolf"))));
else
solve_tolf=*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,"dynatol"))));
mxArray *mxa=mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,"fname"));
int buflen=mxGetM(mxa) * mxGetN(mxa) + 1;
char *fname;
......@@ -237,41 +252,31 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
int status = mxGetString(mxa, fname, buflen);
if (status != 0)
mexWarnMsgTxt("Not enough space. Filename is truncated.");
col_y=mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,"endo_simul")));;
solve_tolf=*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,"dynatol"))));
size_of_direction=col_y*row_y*sizeof(double);
y=(double*)mxMalloc(size_of_direction);
ya=(double*)mxMalloc(size_of_direction);
int size_of_direction=col_y*row_y*sizeof(double);
double * y=(double*)mxMalloc(size_of_direction);
double * ya=(double*)mxMalloc(size_of_direction);
direction=(double*)mxMalloc(size_of_direction);
memset(direction,0,size_of_direction);
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++)
x[i]=double(xd[i]);
for (i=0;i<row_y*col_y;i++)
y[i]=double(yd[i]);
y_size=row_y;
x_size=col_x/*row_x*/;
nb_row_x=row_x;
/*for(int i=0; i<y_size; i++)
{
for(int it_=0; it_<8;it_++)
mexPrintf("y[t=%d, var=%d]=%f ",it_+1, i+1, y[(it_)*y_size+i]);
mexPrintf("\n");
}
int y_size=row_y;
int nb_row_x=row_x;
for(int i=0; i<col_x; i++)
{
for(int it_=0; it_<8;it_++)
mexPrintf("x[t=%d, var=%d]=%f ",it_, i+1, x[it_+i*nb_row_x]);
mexPrintf("\n");
}*/
/*int it_ = y_kmin;
for (int j = 0; j < y_size; j++)
mexPrintf(" variable %d at time %d and %d = %f\n", j+1, it_, it_+1, y[j+it_*y_size]);*/
t0= clock();
clock_t t0= clock();
Interpreter interprete(params, y, ya