commit
a5e9c7bef4
|
@ -107,12 +107,20 @@ fclose(fid);
|
||||||
% Set newline code (ok for *nix, check for mac and windows)
|
% Set newline code (ok for *nix, check for mac and windows)
|
||||||
if isunix
|
if isunix
|
||||||
newline_code = 10;
|
newline_code = 10;
|
||||||
|
elseif ispc
|
||||||
|
newline_code = 13;
|
||||||
|
elseif ismac
|
||||||
|
newline_code = 10;
|
||||||
else
|
else
|
||||||
error('readcsv:: Not implemented for your OS!')
|
error('readcsv:: Not implemented for your OS!')
|
||||||
end
|
end
|
||||||
|
|
||||||
% Get the positions of the end-of-line code;
|
% Get the positions of the end-of-line code;
|
||||||
end_of_line_locations = find(bfile==newline_code);
|
end_of_line_locations = find(bfile==newline_code);
|
||||||
|
if ispc && isempty(end_of_line_locations)
|
||||||
|
newline_code=10;
|
||||||
|
end_of_line_locations = find(bfile==newline_code);
|
||||||
|
end;
|
||||||
tmp = find(bfile==newline_code);
|
tmp = find(bfile==newline_code);
|
||||||
|
|
||||||
% Get the number of lines in the file.
|
% Get the number of lines in the file.
|
||||||
|
|
|
@ -102,6 +102,26 @@ switch (extension)
|
||||||
end
|
end
|
||||||
dyn_data_01(:,dyn_i_01) = dyn_tmp_01;
|
dyn_data_01(:,dyn_i_01) = dyn_tmp_01;
|
||||||
end
|
end
|
||||||
|
case '.csv'
|
||||||
|
[freq,init,data,varlist] = load_csv_file_data(fullname);
|
||||||
|
disp('size(data)');
|
||||||
|
size(data)
|
||||||
|
% for i=1:length(varlist)
|
||||||
|
% if isnan(varlist)
|
||||||
|
% varlist(1,i) = ' ';
|
||||||
|
% end
|
||||||
|
% end
|
||||||
|
%var_names_01 = deblank(var_names_01);
|
||||||
|
|
||||||
|
for dyn_i_01=1:var_size_01
|
||||||
|
iv = strmatch(deblank(var_names_01(dyn_i_01,:)),varlist,'exact') + 1;
|
||||||
|
dyn_tmp_01 = [data(2:end,iv)]';
|
||||||
|
if length(dyn_tmp_01) > dyn_size_01 && dyn_size_01 > 0
|
||||||
|
cd(old_pwd)
|
||||||
|
error('data size is too large')
|
||||||
|
end
|
||||||
|
dyn_data_01(:,dyn_i_01) = dyn_tmp_01;
|
||||||
|
end
|
||||||
otherwise
|
otherwise
|
||||||
cd(old_pwd)
|
cd(old_pwd)
|
||||||
error(['Unsupported extension for datafile: ' extension])
|
error(['Unsupported extension for datafile: ' extension])
|
||||||
|
|
|
@ -74,69 +74,11 @@ g1=spalloc( Blck_size*periods, Jacobian_Size, nze*periods);
|
||||||
reduced = 0;
|
reduced = 0;
|
||||||
while ~(cvg==1 || iter>maxit_),
|
while ~(cvg==1 || iter>maxit_),
|
||||||
[r, y, g1, g2, g3, b]=feval(fname, y, x, params, steady_state, periods, 0, y_kmin, Blck_size);
|
[r, y, g1, g2, g3, b]=feval(fname, y, x, params, steady_state, periods, 0, y_kmin, Blck_size);
|
||||||
% fjac = zeros(Blck_size, Blck_size*(y_kmin_l+1+y_kmax_l));
|
preconditioner = 2;
|
||||||
% 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);
|
|
||||||
%num2str(residual,' %1.6f')
|
|
||||||
%jac_ = g1(1:(y_kmin)*Blck_size, 1:(y_kmin+1+y_kmax_l)*Blck_size);
|
|
||||||
%jac_
|
|
||||||
|
|
||||||
g1a=g1(:, y_kmin*Blck_size+1:(periods+y_kmin)*Blck_size);
|
g1a=g1(:, y_kmin*Blck_size+1:(periods+y_kmin)*Blck_size);
|
||||||
term1 = g1(:, 1:y_kmin_l*Blck_size)*reshape(y(1+y_kmin-y_kmin_l:y_kmin,y_index)',1,y_kmin_l*Blck_size)';
|
term1 = g1(:, 1:y_kmin_l*Blck_size)*reshape(y(1+y_kmin-y_kmin_l:y_kmin,y_index)',1,y_kmin_l*Blck_size)';
|
||||||
term2 = g1(:, (periods+y_kmin_l)*Blck_size+1:(periods+y_kmin_l+y_kmax_l)*Blck_size)*reshape(y(periods+y_kmin+1:periods+y_kmin+y_kmax_l,y_index)',1,y_kmax_l*Blck_size)';
|
term2 = g1(:, (periods+y_kmin_l)*Blck_size+1:(periods+y_kmin_l+y_kmax_l)*Blck_size)*reshape(y(periods+y_kmin+1:periods+y_kmin+y_kmax_l,y_index)',1,y_kmax_l*Blck_size)';
|
||||||
b = b - term1 - term2;
|
b = b - term1 - term2;
|
||||||
|
|
||||||
% fid = fopen(['result' num2str(iter)],'w');
|
|
||||||
% fg1a = full(g1a);
|
|
||||||
% fprintf(fid,'%d\n',size(fg1a,1));
|
|
||||||
% fprintf(fid,'%d\n',size(fg1a,2));
|
|
||||||
% fprintf(fid,'%5.14f\n',fg1a);
|
|
||||||
% fprintf(fid,'%d\n',size(b,1));
|
|
||||||
% fprintf(fid,'%5.14f\n',b);
|
|
||||||
% fclose(fid);
|
|
||||||
% return;
|
|
||||||
%ipconfigb_ = b(1:(1+y_kmin)*Blck_size);
|
|
||||||
%b_
|
|
||||||
|
|
||||||
|
|
||||||
[max_res, max_indx]=max(max(abs(r')));
|
[max_res, max_indx]=max(max(abs(r')));
|
||||||
if(~isreal(r))
|
if(~isreal(r))
|
||||||
max_res = (-max_res^2)^0.5;
|
max_res = (-max_res^2)^0.5;
|
||||||
|
@ -255,7 +197,33 @@ while ~(cvg==1 || iter>maxit_),
|
||||||
elseif(stack_solve_algo==2),
|
elseif(stack_solve_algo==2),
|
||||||
flag1=1;
|
flag1=1;
|
||||||
while(flag1>0)
|
while(flag1>0)
|
||||||
[L1, U1]=luinc(g1a,luinc_tol);
|
if preconditioner == 2
|
||||||
|
[L1, U1]=luinc(g1a,luinc_tol);
|
||||||
|
elseif preconditioner == 3
|
||||||
|
Size = Blck_size;
|
||||||
|
gss1 = g1a(Size + 1: 2*Size,Size + 1: 2*Size) + g1a(Size + 1: 2*Size,2*Size+1: 3*Size);
|
||||||
|
[L1, U1]=lu(gss1);
|
||||||
|
L(1:Size,1:Size) = L1;
|
||||||
|
U(1:Size,1:Size) = U1;
|
||||||
|
gss2 = g1a(Size + 1: 2*Size,1: Size) + g1a(Size + 1: 2*Size,Size+1: 2*Size) + g1a(Size + 1: 2*Size,2*Size+1: 3*Size);
|
||||||
|
[L2, U2]=lu(gss2);
|
||||||
|
L(Size+1:(periods-1)*Size,Size+1:(periods-1)*Size) = kron(eye(periods-2), L2);
|
||||||
|
U(Size+1:(periods-1)*Size,Size+1:(periods-1)*Size) = kron(eye(periods-2), U2);
|
||||||
|
gss2 = g1a(Size + 1: 2*Size,1: Size) + g1a(Size + 1: 2*Size,Size+1: 2*Size);
|
||||||
|
[L3, U3]=lu(gss2);
|
||||||
|
L((periods-1)*Size+1:periods*Size,(periods-1)*Size+1:periods*Size) = L3;
|
||||||
|
U((periods-1)*Size+1:periods*Size,(periods-1)*Size+1:periods*Size) = U3;
|
||||||
|
L1 = L;
|
||||||
|
U1 = U;
|
||||||
|
elseif preconditioner == 4
|
||||||
|
Size = Blck_size;
|
||||||
|
gss1 = g1a(1: 3*Size, 1: 3*Size);
|
||||||
|
[L, U] = lu(gss1);
|
||||||
|
L1 = kron(eye(ceil(periods/3)),L);
|
||||||
|
U1 = kron(eye(ceil(periods/3)),U);
|
||||||
|
L1 = L1(1:periods * Size, 1:periods * Size);
|
||||||
|
U1 = U1(1:periods * Size, 1:periods * Size);
|
||||||
|
end;
|
||||||
[za,flag1] = gmres(g1a,b,Blck_size,1e-6,Blck_size*periods,L1,U1);
|
[za,flag1] = gmres(g1a,b,Blck_size,1e-6,Blck_size*periods,L1,U1);
|
||||||
if (flag1>0 || reduced)
|
if (flag1>0 || reduced)
|
||||||
if(flag1==1)
|
if(flag1==1)
|
||||||
|
@ -276,8 +244,25 @@ while ~(cvg==1 || iter>maxit_),
|
||||||
elseif(stack_solve_algo==3),
|
elseif(stack_solve_algo==3),
|
||||||
flag1=1;
|
flag1=1;
|
||||||
while(flag1>0)
|
while(flag1>0)
|
||||||
[L1, U1]=luinc(g1a,luinc_tol);
|
if (preconditioner == 3)
|
||||||
[za,flag1] = bicgstab(g1a,b,1e-7,Blck_size*periods,L1,U1);
|
Size = Blck_size;
|
||||||
|
gss0 = g1a(Size + 1: 2*Size,1: Size) + g1a(Size + 1: 2*Size,Size+1: 2*Size) + g1a(Size + 1: 2*Size,2*Size+1: 3*Size);
|
||||||
|
[L1, U1]=lu(gss0);
|
||||||
|
P1 = eye(size(gss0));
|
||||||
|
Q1 = eye(size(gss0));
|
||||||
|
L = kron(eye(periods),L1);
|
||||||
|
U = kron(eye(periods),U1);
|
||||||
|
P = kron(eye(periods),P1);
|
||||||
|
Q = kron(eye(periods),Q1);
|
||||||
|
[za,flag1] = bicgstab1(g1a,b,1e-7,Blck_size*periods,L,U, P, Q);
|
||||||
|
else
|
||||||
|
Size = Blck_size;
|
||||||
|
gss0 = g1a(Size + 1: 2*Size,1: Size) + g1a(Size + 1: 2*Size,Size+1: 2*Size) + g1a(Size + 1: 2*Size,2*Size+1: 3*Size);
|
||||||
|
[L1, U1]=lu(gss0);
|
||||||
|
L1 = kron(eye(periods),L1);
|
||||||
|
U1 = kron(eye(periods),U1);
|
||||||
|
[za,flag1] = bicgstab(g1a,b,1e-7,Blck_size*periods,L1,U1);
|
||||||
|
end;
|
||||||
if (flag1>0 || reduced)
|
if (flag1>0 || reduced)
|
||||||
if(flag1==1)
|
if(flag1==1)
|
||||||
disp(['Error in simul: No convergence inside BICGSTAB after ' num2str(periods*10,'%6d') ' iterations, in block' num2str(Block_Size,'%3d')]);
|
disp(['Error in simul: No convergence inside BICGSTAB after ' num2str(periods*10,'%6d') ' iterations, in block' num2str(Block_Size,'%3d')]);
|
||||||
|
|
|
@ -1,6 +1,6 @@
|
||||||
noinst_PROGRAMS = bytecode
|
noinst_PROGRAMS = bytecode
|
||||||
|
|
||||||
bytecode_CPPFLAGS = $(AM_CPPFLAGS) -I$(top_srcdir)/../../sources/bytecode -I$(top_srcdir)/../../../preprocessor
|
bytecode_CPPFLAGS = $(AM_CPPFLAGS) -I$(top_srcdir)/../../sources -I$(top_srcdir)/../../sources/bytecode -I$(top_srcdir)/../../../preprocessor
|
||||||
|
|
||||||
TOPDIR = $(top_srcdir)/../../sources/bytecode
|
TOPDIR = $(top_srcdir)/../../sources/bytecode
|
||||||
|
|
||||||
|
@ -9,8 +9,10 @@ nodist_bytecode_SOURCES = \
|
||||||
$(TOPDIR)/Interpreter.cc \
|
$(TOPDIR)/Interpreter.cc \
|
||||||
$(TOPDIR)/Mem_Mngr.cc \
|
$(TOPDIR)/Mem_Mngr.cc \
|
||||||
$(TOPDIR)/SparseMatrix.cc \
|
$(TOPDIR)/SparseMatrix.cc \
|
||||||
|
$(TOPDIR)/Evaluate.cc \
|
||||||
$(TOPDIR)/Interpreter.hh \
|
$(TOPDIR)/Interpreter.hh \
|
||||||
$(TOPDIR)/Mem_Mngr.hh \
|
$(TOPDIR)/Mem_Mngr.hh \
|
||||||
$(TOPDIR)/SparseMatrix.hh \
|
$(TOPDIR)/SparseMatrix.hh \
|
||||||
|
$(TOPDIR)/Evaluate.hh \
|
||||||
$(TOPDIR)/ErrorHandling.hh
|
$(TOPDIR)/ErrorHandling.hh
|
||||||
|
|
||||||
|
|
|
@ -26,10 +26,13 @@
|
||||||
#endif
|
#endif
|
||||||
#include "block_kalman_filter.h"
|
#include "block_kalman_filter.h"
|
||||||
using namespace std;
|
using namespace std;
|
||||||
//#define BLAS
|
#define BLAS
|
||||||
#define DIRECT
|
//#define CUBLAS
|
||||||
|
|
||||||
|
|
||||||
|
#ifdef CUBLAS
|
||||||
|
#include <cuda_runtime.h>
|
||||||
|
#include <cublas_v2.h>
|
||||||
|
#endif
|
||||||
void
|
void
|
||||||
mexDisp(mxArray* P)
|
mexDisp(mxArray* P)
|
||||||
{
|
{
|
||||||
|
@ -157,7 +160,7 @@ BlockKalmanFilter::BlockKalmanFilter(int nlhs, mxArray *plhs[], int nrhs, const
|
||||||
if (missing_observations)
|
if (missing_observations)
|
||||||
{
|
{
|
||||||
if (! mxIsCell (prhs[0]))
|
if (! mxIsCell (prhs[0]))
|
||||||
DYN_MEX_FUNC_ERR_MSG_TXT("the first input argument of block_missing_observations_kalman_filter must be a Call Array.");
|
DYN_MEX_FUNC_ERR_MSG_TXT("the first input argument of block_missing_observations_kalman_filter must be a Cell Array.");
|
||||||
pdata_index = prhs[0];
|
pdata_index = prhs[0];
|
||||||
if (! mxIsDouble (prhs[1]))
|
if (! mxIsDouble (prhs[1]))
|
||||||
DYN_MEX_FUNC_ERR_MSG_TXT("the second input argument of block_missing_observations_kalman_filter must be a scalar.");
|
DYN_MEX_FUNC_ERR_MSG_TXT("the second input argument of block_missing_observations_kalman_filter must be a scalar.");
|
||||||
|
@ -234,14 +237,13 @@ BlockKalmanFilter::BlockKalmanFilter(int nlhs, mxArray *plhs[], int nrhs, const
|
||||||
*a = mxGetPr(pa);
|
*a = mxGetPr(pa);
|
||||||
tmp_a = (double*)mxMalloc(n * sizeof(double));
|
tmp_a = (double*)mxMalloc(n * sizeof(double));
|
||||||
dF = 0.0; // det(F).
|
dF = 0.0; // det(F).
|
||||||
p_tmp = mxCreateDoubleMatrix(n, n_state, mxREAL);
|
|
||||||
*tmp = mxGetPr(p_tmp);
|
|
||||||
p_tmp1 = mxCreateDoubleMatrix(n, n_shocks, mxREAL);
|
p_tmp1 = mxCreateDoubleMatrix(n, n_shocks, mxREAL);
|
||||||
tmp1 = mxGetPr(p_tmp1);
|
tmp1 = mxGetPr(p_tmp1);
|
||||||
t = 0; // Initialization of the time index.
|
t = 0; // Initialization of the time index.
|
||||||
plik = mxCreateDoubleMatrix(smpl, 1, mxREAL);
|
plik = mxCreateDoubleMatrix(smpl, 1, mxREAL);
|
||||||
lik = mxGetPr(plik);
|
lik = mxGetPr(plik);
|
||||||
Inf = mxGetInf();
|
Inf = mxGetInf();
|
||||||
LIK = 0.0; // Default value of the log likelihood.
|
LIK = 0.0; // Default value of the log likelihood.
|
||||||
notsteady = true; // Steady state flag.
|
notsteady = true; // Steady state flag.
|
||||||
F_singular = true;
|
F_singular = true;
|
||||||
|
@ -287,6 +289,22 @@ BlockKalmanFilter::BlockKalmanFilter(int nlhs, mxArray *plhs[], int nrhs, const
|
||||||
iw = (lapack_int*)mxMalloc(pp * sizeof(lapack_int));
|
iw = (lapack_int*)mxMalloc(pp * sizeof(lapack_int));
|
||||||
ipiv = (lapack_int*)mxMalloc(pp * sizeof(lapack_int));
|
ipiv = (lapack_int*)mxMalloc(pp * sizeof(lapack_int));
|
||||||
info = 0;
|
info = 0;
|
||||||
|
#ifdef BLAS || CUBLAS
|
||||||
|
p_tmp = mxCreateDoubleMatrix(n, n, mxREAL);
|
||||||
|
*tmp = mxGetPr(p_tmp);
|
||||||
|
p_P_t_t1 = mxCreateDoubleMatrix(n, n, mxREAL);
|
||||||
|
*P_t_t1 = mxGetPr(p_P_t_t1);
|
||||||
|
pK = mxCreateDoubleMatrix(n, n, mxREAL);
|
||||||
|
*K = mxGetPr(pK);
|
||||||
|
p_K_P = mxCreateDoubleMatrix(n, n, mxREAL);
|
||||||
|
*K_P = mxGetPr(p_K_P);
|
||||||
|
oldK = (double*)mxMalloc(n * n * sizeof(double));
|
||||||
|
*P_mf = (double*)mxMalloc(n * n * sizeof(double));
|
||||||
|
for (int i = 0; i < n * n; i++)
|
||||||
|
oldK[i] = Inf;
|
||||||
|
#else
|
||||||
|
p_tmp = mxCreateDoubleMatrix(n, n_state, mxREAL);
|
||||||
|
*tmp = mxGetPr(p_tmp);
|
||||||
p_P_t_t1 = mxCreateDoubleMatrix(n_state, n_state, mxREAL);
|
p_P_t_t1 = mxCreateDoubleMatrix(n_state, n_state, mxREAL);
|
||||||
*P_t_t1 = mxGetPr(p_P_t_t1);
|
*P_t_t1 = mxGetPr(p_P_t_t1);
|
||||||
pK = mxCreateDoubleMatrix(n, pp, mxREAL);
|
pK = mxCreateDoubleMatrix(n, pp, mxREAL);
|
||||||
|
@ -297,6 +315,7 @@ BlockKalmanFilter::BlockKalmanFilter(int nlhs, mxArray *plhs[], int nrhs, const
|
||||||
*P_mf = (double*)mxMalloc(n * pp * sizeof(double));
|
*P_mf = (double*)mxMalloc(n * pp * sizeof(double));
|
||||||
for (int i = 0; i < n * pp; i++)
|
for (int i = 0; i < n * pp; i++)
|
||||||
oldK[i] = Inf;
|
oldK[i] = Inf;
|
||||||
|
#endif
|
||||||
}
|
}
|
||||||
|
|
||||||
void
|
void
|
||||||
|
@ -424,17 +443,17 @@ BlockKalmanFilter::block_kalman_filter(int nlhs, mxArray *plhs[], double *P_mf,
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
/* Computes the norm of iF */
|
/* Computes the norm of iF */
|
||||||
double anorm = dlange("1", &size_d_index, &size_d_index, iF, &size_d_index, w);
|
double anorm = dlange("1", &size_d_index, &size_d_index, iF, &size_d_index, w);
|
||||||
//mexPrintf("anorm = %f\n",anorm);
|
//mexPrintf("anorm = %f\n",anorm);
|
||||||
|
|
||||||
/* Modifies F in place with a LU decomposition */
|
/* Modifies F in place with a LU decomposition */
|
||||||
dgetrf(&size_d_index, &size_d_index, iF, &size_d_index, ipiv, &info);
|
dgetrf(&size_d_index, &size_d_index, iF, &size_d_index, ipiv, &info);
|
||||||
if (info != 0) fprintf(stderr, "dgetrf failure with error %d\n", (int) info);
|
if (info != 0) mexPrintf("dgetrf failure with error %d\n", (int) info);
|
||||||
|
|
||||||
/* Computes the reciprocal norm */
|
/* Computes the reciprocal norm */
|
||||||
dgecon("1", &size_d_index, iF, &size_d_index, &anorm, &rcond, w, iw, &info);
|
dgecon("1", &size_d_index, iF, &size_d_index, &anorm, &rcond, w, iw, &info);
|
||||||
if (info != 0) fprintf(stderr, "dgecon failure with error %d\n", (int) info);
|
if (info != 0) mexPrintf("dgecon failure with error %d\n", (int) info);
|
||||||
|
|
||||||
if (rcond < kalman_tol)
|
if (rcond < kalman_tol)
|
||||||
if (not_all_abs_F_bellow_crit(F, size_d_index * size_d_index, kalman_tol)) //~all(abs(F(:))<kalman_tol)
|
if (not_all_abs_F_bellow_crit(F, size_d_index * size_d_index, kalman_tol)) //~all(abs(F(:))<kalman_tol)
|
||||||
|
@ -506,7 +525,7 @@ BlockKalmanFilter::block_kalman_filter(int nlhs, mxArray *plhs[], double *P_mf,
|
||||||
//iF = inv(F);
|
//iF = inv(F);
|
||||||
//int lwork = 4/*2*/* pp;
|
//int lwork = 4/*2*/* pp;
|
||||||
dgetri(&size_d_index, iF, &size_d_index, ipiv, w, &lw, &info);
|
dgetri(&size_d_index, iF, &size_d_index, ipiv, w, &lw, &info);
|
||||||
if (info != 0) fprintf(stderr, "dgetri failure with error %d\n", (int) info);
|
if (info != 0) mexPrintf("dgetri failure with error %d\n", (int) info);
|
||||||
|
|
||||||
//lik(t) = log(dF)+transpose(v)*iF*v;
|
//lik(t) = log(dF)+transpose(v)*iF*v;
|
||||||
#ifdef USE_OMP
|
#ifdef USE_OMP
|
||||||
|
@ -628,14 +647,126 @@ BlockKalmanFilter::block_kalman_filter(int nlhs, mxArray *plhs[], double *P_mf,
|
||||||
double one = 1.0;
|
double one = 1.0;
|
||||||
double zero = 0.0;
|
double zero = 0.0;
|
||||||
memcpy(P, QQ, n * n *sizeof(double));
|
memcpy(P, QQ, n * n *sizeof(double));
|
||||||
dsymm("R", "U", &n, &n,
|
blas_int n_b = n;
|
||||||
&one, P_t_t1, &n,
|
/*mexPrintf("sizeof(n_b)=%d, n_b=%d, sizeof(n)=%d, n=%d\n",sizeof(n_b),n_b,sizeof(n),n);
|
||||||
T, &n, &zero,
|
mexEvalString("drawnow;");*/
|
||||||
tmp, &n);
|
dsymm("R", "U", &n_b, &n_b,
|
||||||
dgemm("N", "T", &n, &n,
|
&one, P_t_t1, &n_b,
|
||||||
&n, &one, tmp, &n,
|
T, &n_b, &zero,
|
||||||
T, &n, &one,
|
tmp, &n_b);
|
||||||
P, &n);
|
dgemm("N", "T", &n_b, &n_b,
|
||||||
|
&n_b, &one, tmp, &n_b,
|
||||||
|
T, &n_b, &one,
|
||||||
|
P, &n_b);
|
||||||
|
#else
|
||||||
|
#ifdef CUBLAS
|
||||||
|
for (int i = 0; i < n; i++)
|
||||||
|
for (int j = i; j < n; j++)
|
||||||
|
{
|
||||||
|
double res = 0.0;
|
||||||
|
//int j_pp = j * pp;
|
||||||
|
for (int k = 0; k < size_d_index; k++)
|
||||||
|
res += K[i + k * n] * P_mf[k + j * size_d_index];
|
||||||
|
K_P[i * n + j] = K_P[j * n + i] = res;
|
||||||
|
}
|
||||||
|
//#pragma omp parallel for shared(P, K_P, P_t_t1)
|
||||||
|
for (int i = size_d_index; i < n; i++)
|
||||||
|
for (int j = i; j < n; j++)
|
||||||
|
{
|
||||||
|
unsigned int k = i * n + j;
|
||||||
|
P_t_t1[j * n + i] = P_t_t1[k] = P[k] - K_P[k];
|
||||||
|
}
|
||||||
|
mexPrintf("CudaBLAS\n");
|
||||||
|
mexEvalString("drawnow;");
|
||||||
|
double one = 1.0;
|
||||||
|
double zero = 0.0;
|
||||||
|
cublasStatus_t status;
|
||||||
|
cublasHandle_t handle;
|
||||||
|
status = cublasCreate(&handle);
|
||||||
|
if (status != CUBLAS_STATUS_SUCCESS)
|
||||||
|
{
|
||||||
|
mexPrintf("!!!! CUBLAS initialization error\n");
|
||||||
|
return false;
|
||||||
|
}
|
||||||
|
/*int device;
|
||||||
|
cudaGetDevice(&device);*/
|
||||||
|
int n2 = n * n;
|
||||||
|
double* d_A = 0;
|
||||||
|
double* d_B = 0;
|
||||||
|
double* d_C = 0;
|
||||||
|
double* d_D = 0;
|
||||||
|
// Allocate device memory for the matrices
|
||||||
|
if (cudaMalloc((void**)&d_A, n2 * sizeof(double)) != cudaSuccess)
|
||||||
|
{
|
||||||
|
mexPrintf("!!!! device memory allocation error (allocate A)\n");
|
||||||
|
return false;
|
||||||
|
}
|
||||||
|
if (cudaMalloc((void**)&d_B, n2 * sizeof(d_B[0])) != cudaSuccess)
|
||||||
|
{
|
||||||
|
mexPrintf("!!!! device memory allocation error (allocate B)\n");
|
||||||
|
return false;
|
||||||
|
}
|
||||||
|
if (cudaMalloc((void**)&d_C, n2 * sizeof(d_C[0])) != cudaSuccess)
|
||||||
|
{
|
||||||
|
mexPrintf("!!!! device memory allocation error (allocate C)\n");
|
||||||
|
return false;
|
||||||
|
}
|
||||||
|
if (cudaMalloc((void**)&d_D, n2 * sizeof(d_D[0])) != cudaSuccess)
|
||||||
|
{
|
||||||
|
mexPrintf("!!!! device memory allocation error (allocate D)\n");
|
||||||
|
return false;
|
||||||
|
}
|
||||||
|
// Initialize the device matrices with the host matrices
|
||||||
|
status = cublasSetVector(n2, sizeof(P_t_t1[0]), P_t_t1, 1, d_A, 1);
|
||||||
|
if (status != CUBLAS_STATUS_SUCCESS)
|
||||||
|
{
|
||||||
|
mexPrintf("!!!! device access error (write A)\n");
|
||||||
|
return false;
|
||||||
|
}
|
||||||
|
status = cublasSetVector(n2, sizeof(T[0]), T, 1, d_B, 1);
|
||||||
|
if (status != CUBLAS_STATUS_SUCCESS)
|
||||||
|
{
|
||||||
|
mexPrintf("!!!! device access error (write B)\n");
|
||||||
|
return false;
|
||||||
|
}
|
||||||
|
status = cublasSetVector(n2, sizeof(tmp[0]), tmp, 1, d_C, 1);
|
||||||
|
if (status != CUBLAS_STATUS_SUCCESS)
|
||||||
|
{
|
||||||
|
mexPrintf("!!!! device access error (write C)\n");
|
||||||
|
return false;
|
||||||
|
}
|
||||||
|
mexPrintf("just before calling\n");
|
||||||
|
mexEvalString("drawnow;");
|
||||||
|
status = cublasSetVector(n2, sizeof(QQ[0]), QQ, 1, d_D, 1);
|
||||||
|
if (status != CUBLAS_STATUS_SUCCESS)
|
||||||
|
{
|
||||||
|
mexPrintf("!!!! device access error (write D)\n");
|
||||||
|
return false;
|
||||||
|
}
|
||||||
|
|
||||||
|
// Performs operation using plain C code
|
||||||
|
|
||||||
|
cublasDsymm(handle, CUBLAS_SIDE_RIGHT, CUBLAS_FILL_MODE_UPPER, n, n,
|
||||||
|
&one, d_A, n,
|
||||||
|
d_B, n, &zero,
|
||||||
|
d_C, n);
|
||||||
|
/*dgemm("N", "T", &n_b, &n_b,
|
||||||
|
&n_b, &one, tmp, &n_b,
|
||||||
|
T, &n_b, &one,
|
||||||
|
P, &n_b);*/
|
||||||
|
cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, n, n,
|
||||||
|
n, &one, d_C, n,
|
||||||
|
d_B, n, &one,
|
||||||
|
d_D, n);
|
||||||
|
//double_symm(n, &one, h_A, h_B, &zero, h_C);
|
||||||
|
|
||||||
|
status = cublasGetVector(n2, sizeof(P[0]), d_D, 1, P, 1);
|
||||||
|
if (status != CUBLAS_STATUS_SUCCESS)
|
||||||
|
{
|
||||||
|
mexPrintf("!!!! device access error (read P)\n");
|
||||||
|
return false;
|
||||||
|
}
|
||||||
|
|
||||||
#else
|
#else
|
||||||
#ifdef USE_OMP
|
#ifdef USE_OMP
|
||||||
#pragma omp parallel for shared(K_P) num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
|
#pragma omp parallel for shared(K_P) num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
|
||||||
|
@ -717,7 +848,7 @@ BlockKalmanFilter::block_kalman_filter(int nlhs, mxArray *plhs[], double *P_mf,
|
||||||
P[i + j * n] = P[j + i * n];
|
P[i + j * n] = P[j + i * n];
|
||||||
}
|
}
|
||||||
#endif
|
#endif
|
||||||
|
#endif
|
||||||
if (t >= no_more_missing_observations)
|
if (t >= no_more_missing_observations)
|
||||||
{
|
{
|
||||||
double max_abs = 0.0;
|
double max_abs = 0.0;
|
||||||
|
|
|
@ -23,24 +23,138 @@
|
||||||
#include <cstring>
|
#include <cstring>
|
||||||
#include <iostream>
|
#include <iostream>
|
||||||
#include <sstream>
|
#include <sstream>
|
||||||
|
#include <map>
|
||||||
|
#define BYTE_CODE
|
||||||
#include "CodeInterpreter.hh"
|
#include "CodeInterpreter.hh"
|
||||||
#ifdef DEBUG_EX
|
#ifdef DEBUG_EX
|
||||||
# include <math>
|
# include <math.h>
|
||||||
# include "mex_interface.hh"
|
# include "mex_interface.hh"
|
||||||
#endif
|
#endif
|
||||||
|
|
||||||
|
#ifdef OCTAVE_MEX_FILE
|
||||||
|
# define CHAR_LENGTH 1
|
||||||
|
#else
|
||||||
|
# define CHAR_LENGTH 2
|
||||||
|
#endif
|
||||||
|
|
||||||
|
#ifdef _MSC_VER
|
||||||
|
#include <limits>
|
||||||
|
#define M_E 2.71828182845904523536
|
||||||
|
#define M_LOG2E 1.44269504088896340736
|
||||||
|
#define M_LOG10E 0.434294481903251827651
|
||||||
|
#define M_LN2 0.693147180559945309417
|
||||||
|
#define M_LN10 2.30258509299404568402
|
||||||
|
#define M_PI 3.14159265358979323846
|
||||||
|
#define M_PI_2 1.57079632679489661923
|
||||||
|
#define M_PI_4 0.785398163397448309616
|
||||||
|
#define M_1_PI 0.318309886183790671538
|
||||||
|
#define M_2_PI 0.636619772367581343076
|
||||||
|
#define M_1_SQRTPI 0.564189583547756286948
|
||||||
|
#define M_2_SQRTPI 1.12837916709551257390
|
||||||
|
#define M_SQRT2 1.41421356237309504880
|
||||||
|
#define M_SQRT_2 0.707106781186547524401
|
||||||
|
#define NAN numeric_limits<double>::quiet_NaN()
|
||||||
|
|
||||||
|
#define isnan(x) _isnan(x)
|
||||||
|
#define isinf(x) (!_finite(x))
|
||||||
|
#define fpu_error(x) (isinf(x) || isnan(x))
|
||||||
|
|
||||||
|
|
||||||
|
class MSVCpp_missings
|
||||||
|
{
|
||||||
|
public:
|
||||||
|
inline double
|
||||||
|
asinh(double x) const
|
||||||
|
{
|
||||||
|
if(x==0.0)
|
||||||
|
return 0.0;
|
||||||
|
double ax = abs(x);
|
||||||
|
return log(x+ax*sqrt(1.+1./(ax*ax)));
|
||||||
|
}
|
||||||
|
|
||||||
|
inline double
|
||||||
|
acosh(double x) const
|
||||||
|
{
|
||||||
|
if(x==0.0)
|
||||||
|
return 0.0;
|
||||||
|
double ax = abs(x);
|
||||||
|
return log(x+ax*sqrt(1.-1./(ax*ax)));
|
||||||
|
}
|
||||||
|
|
||||||
|
inline double
|
||||||
|
atanh(double x) const
|
||||||
|
{
|
||||||
|
return log((1+x)/(1-x))/2;
|
||||||
|
}
|
||||||
|
|
||||||
|
inline double
|
||||||
|
erf(double x) const
|
||||||
|
{
|
||||||
|
const double a1 = -1.26551223, a2 = 1.00002368,
|
||||||
|
a3 = 0.37409196, a4 = 0.09678418,
|
||||||
|
a5 = -0.18628806, a6 = 0.27886807,
|
||||||
|
a7 = -1.13520398, a8 = 1.48851587,
|
||||||
|
a9 = -0.82215223, a10 = 0.17087277;
|
||||||
|
double v = 1;
|
||||||
|
double z = abs(x);
|
||||||
|
if (z <= 0)
|
||||||
|
return v;
|
||||||
|
double t = 1 / (1 + 0.5 * z);
|
||||||
|
v = t*exp((-z*z) +a1+t*(a2+t*(a3+t*(a4+t*(a5+t*(a6+t*(a7+t*(a8+t*(a9+t*a10)))))))));
|
||||||
|
if (x < 0)
|
||||||
|
v = 2 - v;
|
||||||
|
return 1 - v;
|
||||||
|
}
|
||||||
|
|
||||||
|
inline double
|
||||||
|
nearbyint(double x) const
|
||||||
|
{
|
||||||
|
return floor(x + 0.5);
|
||||||
|
}
|
||||||
|
|
||||||
|
inline double
|
||||||
|
fmax(double x, double y) const
|
||||||
|
{
|
||||||
|
if (x > y)
|
||||||
|
return x;
|
||||||
|
else
|
||||||
|
return y;
|
||||||
|
}
|
||||||
|
|
||||||
|
inline double
|
||||||
|
fmin(double x, double y) const
|
||||||
|
{
|
||||||
|
if (x < y)
|
||||||
|
return x;
|
||||||
|
else
|
||||||
|
return y;
|
||||||
|
}
|
||||||
|
|
||||||
|
};
|
||||||
|
#endif
|
||||||
|
|
||||||
|
|
||||||
//#define DEBUG
|
//#define DEBUG
|
||||||
using namespace std;
|
using namespace std;
|
||||||
|
|
||||||
const int NO_ERROR_ON_EXIT = 0;
|
const int NO_ERROR_ON_EXIT = 0;
|
||||||
const int ERROR_ON_EXIT = 1;
|
const int ERROR_ON_EXIT = 1;
|
||||||
|
|
||||||
|
|
||||||
typedef vector<pair<Tags, void * > > code_liste_type;
|
typedef vector<pair<Tags, void * > > code_liste_type;
|
||||||
typedef code_liste_type::const_iterator it_code_type;
|
typedef code_liste_type::const_iterator it_code_type;
|
||||||
|
|
||||||
|
|
||||||
class GeneralExceptionHandling
|
class GeneralExceptionHandling
|
||||||
{
|
{
|
||||||
string ErrorMsg;
|
string ErrorMsg;
|
||||||
public:
|
public:
|
||||||
|
#ifdef _MSC_VER_
|
||||||
|
~GeneralExceptionHandling()
|
||||||
|
{
|
||||||
|
FreeLibrary(hinstLib);
|
||||||
|
};
|
||||||
|
#endif
|
||||||
GeneralExceptionHandling(string ErrorMsg_arg) : ErrorMsg(ErrorMsg_arg)
|
GeneralExceptionHandling(string ErrorMsg_arg) : ErrorMsg(ErrorMsg_arg)
|
||||||
{
|
{
|
||||||
};
|
};
|
||||||
|
@ -121,6 +235,16 @@ public:
|
||||||
};
|
};
|
||||||
};
|
};
|
||||||
|
|
||||||
|
class UserExceptionHandling : public GeneralExceptionHandling
|
||||||
|
{
|
||||||
|
double value;
|
||||||
|
public:
|
||||||
|
UserExceptionHandling() : GeneralExceptionHandling("Fatal error in bytecode:")
|
||||||
|
{
|
||||||
|
completeErrorMsg(" User break\n");
|
||||||
|
};
|
||||||
|
};
|
||||||
|
|
||||||
class FatalExceptionHandling : public GeneralExceptionHandling
|
class FatalExceptionHandling : public GeneralExceptionHandling
|
||||||
{
|
{
|
||||||
public:
|
public:
|
||||||
|
@ -133,16 +257,40 @@ public:
|
||||||
};
|
};
|
||||||
};
|
};
|
||||||
|
|
||||||
class ErrorMsg
|
struct s_plan
|
||||||
{
|
{
|
||||||
|
string var, exo;
|
||||||
|
int var_num, exo_num;
|
||||||
|
vector<pair<int, double> > per_value;
|
||||||
|
};
|
||||||
|
|
||||||
|
#ifdef MATLAB_MEX_FILE
|
||||||
|
extern "C" bool utIsInterruptPending();
|
||||||
|
#else
|
||||||
|
#include <octave/oct.h>
|
||||||
|
#include <octave/unwind-prot.h>
|
||||||
|
#endif
|
||||||
|
|
||||||
|
#ifdef _MSC_VER
|
||||||
|
class ErrorMsg : public MSVCpp_missings
|
||||||
|
#else
|
||||||
|
class ErrorMsg
|
||||||
|
#endif
|
||||||
|
{
|
||||||
|
private:
|
||||||
|
bool is_load_variable_list;
|
||||||
|
|
||||||
public:
|
public:
|
||||||
|
double *y, *ya;
|
||||||
|
int y_size;
|
||||||
double *T;
|
double *T;
|
||||||
int nb_row_xd, nb_row_x, y_size;
|
int nb_row_xd, nb_row_x;
|
||||||
int y_kmin, y_kmax, periods;
|
int y_kmin, y_kmax, periods;
|
||||||
double *x, *params;
|
double *x, *params;
|
||||||
double *u, *y, *ya;
|
double *u;
|
||||||
double *steady_y, *steady_x;
|
double *steady_y, *steady_x;
|
||||||
double *g2, *g1, *r;
|
double *g2, *g1, *r, *res;
|
||||||
|
vector<s_plan> splan, spfplan;
|
||||||
vector<mxArray *> jacobian_block, jacobian_other_endo_block, jacobian_exo_block, jacobian_det_exo_block;
|
vector<mxArray *> jacobian_block, jacobian_other_endo_block, jacobian_exo_block, jacobian_det_exo_block;
|
||||||
map<unsigned int, double> TEF;
|
map<unsigned int, double> TEF;
|
||||||
map<pair<unsigned int, unsigned int>, double > TEFD;
|
map<pair<unsigned int, unsigned int>, double > TEFD;
|
||||||
|
@ -150,11 +298,12 @@ public:
|
||||||
|
|
||||||
ExpressionType EQN_type;
|
ExpressionType EQN_type;
|
||||||
it_code_type it_code_expr;
|
it_code_type it_code_expr;
|
||||||
unsigned int nb_endo, nb_exo, nb_param;
|
/*unsigned int*/size_t nb_endo, nb_exo, nb_param;
|
||||||
char *P_endo_names, *P_exo_names, *P_param_names;
|
char *P_endo_names, *P_exo_names, *P_param_names;
|
||||||
unsigned int endo_name_length, exo_name_length, param_name_length;
|
size_t/*unsigned int*/ endo_name_length, exo_name_length, param_name_length;
|
||||||
unsigned int EQN_equation, EQN_block, EQN_block_number;
|
unsigned int EQN_equation, EQN_block, EQN_block_number;
|
||||||
unsigned int EQN_dvar1, EQN_dvar2, EQN_dvar3;
|
unsigned int EQN_dvar1, EQN_dvar2, EQN_dvar3;
|
||||||
|
vector<pair<string, pair<SymbolType, unsigned int> > > Variable_list;
|
||||||
|
|
||||||
inline
|
inline
|
||||||
ErrorMsg()
|
ErrorMsg()
|
||||||
|
@ -169,6 +318,7 @@ public:
|
||||||
nb_param = mxGetM(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "param_names")));
|
nb_param = mxGetM(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "param_names")));
|
||||||
param_name_length = mxGetN(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "param_names")));
|
param_name_length = mxGetN(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "param_names")));
|
||||||
P_param_names = (char *) mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "param_names")));
|
P_param_names = (char *) mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "param_names")));
|
||||||
|
is_load_variable_list = false;
|
||||||
}
|
}
|
||||||
|
|
||||||
inline string
|
inline string
|
||||||
|
@ -184,9 +334,9 @@ public:
|
||||||
else
|
else
|
||||||
{
|
{
|
||||||
if (str[i] == '$')
|
if (str[i] == '$')
|
||||||
pos1 = temp.length();
|
pos1 = int(temp.length());
|
||||||
else
|
else
|
||||||
pos2 = temp.length();
|
pos2 = int(temp.length());
|
||||||
if (pos1 >= 0 && pos2 >= 0)
|
if (pos1 >= 0 && pos2 >= 0)
|
||||||
{
|
{
|
||||||
tmp_n.erase(pos1, pos2-pos1+1);
|
tmp_n.erase(pos1, pos2-pos1+1);
|
||||||
|
@ -199,6 +349,50 @@ public:
|
||||||
return temp;
|
return temp;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
inline void
|
||||||
|
load_variable_list()
|
||||||
|
{
|
||||||
|
ostringstream res;
|
||||||
|
for (unsigned int variable_num = 0; variable_num < (unsigned int)nb_endo; variable_num++)
|
||||||
|
{
|
||||||
|
for (unsigned int i = 0; i < endo_name_length; i++)
|
||||||
|
if (P_endo_names[CHAR_LENGTH*(variable_num+i*nb_endo)] != ' ')
|
||||||
|
res << P_endo_names[CHAR_LENGTH*(variable_num+i*nb_endo)];
|
||||||
|
Variable_list.push_back(make_pair(res.str(), make_pair(eEndogenous, variable_num)));
|
||||||
|
}
|
||||||
|
for (unsigned int variable_num = 0; variable_num < (unsigned int)nb_exo; variable_num++)
|
||||||
|
{
|
||||||
|
for (unsigned int i = 0; i < exo_name_length; i++)
|
||||||
|
if (P_exo_names[CHAR_LENGTH*(variable_num+i*nb_exo)] != ' ')
|
||||||
|
res << P_exo_names[CHAR_LENGTH*(variable_num+i*nb_exo)];
|
||||||
|
Variable_list.push_back(make_pair(res.str(), make_pair(eExogenous, variable_num)));
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
inline int
|
||||||
|
get_ID(const string variable_name, SymbolType *variable_type)
|
||||||
|
{
|
||||||
|
if (!is_load_variable_list)
|
||||||
|
{
|
||||||
|
load_variable_list();
|
||||||
|
is_load_variable_list = true;
|
||||||
|
}
|
||||||
|
size_t n = Variable_list.size();
|
||||||
|
int i = 0;
|
||||||
|
bool notfound = true;
|
||||||
|
while (notfound && i < n)
|
||||||
|
{
|
||||||
|
if (variable_name == Variable_list[i].first)
|
||||||
|
{
|
||||||
|
notfound = false;
|
||||||
|
*variable_type = Variable_list[i].second.first;
|
||||||
|
return Variable_list[i].second.second;
|
||||||
|
}
|
||||||
|
i++;
|
||||||
|
}
|
||||||
|
return(-1);
|
||||||
|
}
|
||||||
|
|
||||||
inline string
|
inline string
|
||||||
get_variable(const SymbolType variable_type, const unsigned int variable_num) const
|
get_variable(const SymbolType variable_type, const unsigned int variable_num) const
|
||||||
{
|
{
|
||||||
|
@ -293,7 +487,6 @@ public:
|
||||||
break;
|
break;
|
||||||
default:
|
default:
|
||||||
return ("???");
|
return ("???");
|
||||||
break;
|
|
||||||
}
|
}
|
||||||
else
|
else
|
||||||
switch (EQN_type)
|
switch (EQN_type)
|
||||||
|
@ -342,7 +535,6 @@ public:
|
||||||
break;
|
break;
|
||||||
default:
|
default:
|
||||||
return ("???");
|
return ("???");
|
||||||
break;
|
|
||||||
}
|
}
|
||||||
it_code_type it_code_ret;
|
it_code_type it_code_ret;
|
||||||
Error_loc << endl << add_underscore_to_fpe(" " + print_expression(it_code_expr, evaluate, size, block_num, steady_state, Per_u_, it_, it_code_ret, true));
|
Error_loc << endl << add_underscore_to_fpe(" " + print_expression(it_code_expr, evaluate, size, block_num, steady_state, Per_u_, it_, it_code_ret, true));
|
||||||
|
@ -378,6 +570,12 @@ public:
|
||||||
|
|
||||||
while (go_on)
|
while (go_on)
|
||||||
{
|
{
|
||||||
|
#ifdef OCTAVE_MEX_FILE
|
||||||
|
OCTAVE_QUIT;
|
||||||
|
#else
|
||||||
|
if ( utIsInterruptPending() )
|
||||||
|
throw UserExceptionHandling();
|
||||||
|
#endif
|
||||||
switch (it_code->first)
|
switch (it_code->first)
|
||||||
{
|
{
|
||||||
case FNUMEXPR:
|
case FNUMEXPR:
|
||||||
|
@ -441,7 +639,9 @@ public:
|
||||||
case eParameter:
|
case eParameter:
|
||||||
var = ((FLDV_ *) it_code->second)->get_pos();
|
var = ((FLDV_ *) it_code->second)->get_pos();
|
||||||
#ifdef DEBUG
|
#ifdef DEBUG
|
||||||
mexPrintf("FLDV_ Param var=%d", var);
|
mexPrintf("FLDV_ Param var=%d\n", var);
|
||||||
|
mexPrintf("get_variable(eParameter, var)=%s\n",get_variable(eParameter, var).c_str());
|
||||||
|
mexEvalString("drawnow;");
|
||||||
#endif
|
#endif
|
||||||
Stack.push(get_variable(eParameter, var));
|
Stack.push(get_variable(eParameter, var));
|
||||||
if (compute)
|
if (compute)
|
||||||
|
@ -451,7 +651,10 @@ public:
|
||||||
var = ((FLDV_ *) it_code->second)->get_pos();
|
var = ((FLDV_ *) it_code->second)->get_pos();
|
||||||
lag = ((FLDV_ *) it_code->second)->get_lead_lag();
|
lag = ((FLDV_ *) it_code->second)->get_lead_lag();
|
||||||
#ifdef DEBUG
|
#ifdef DEBUG
|
||||||
mexPrintf("FLDV_ endo var=%d, lag=%d", var, lag);
|
mexPrintf("FLDV_ endo var=%d, lag=%d\n", var, lag);
|
||||||
|
mexPrintf("get_variable(eEndogenous, var)=%s, compute=%d\n",get_variable(eEndogenous, var).c_str(), compute);
|
||||||
|
mexPrintf("it_=%d, lag=%d, y_size=%d, var=%d, y=%x\n", it_, lag, y_size, var, y);
|
||||||
|
mexEvalString("drawnow;");
|
||||||
#endif
|
#endif
|
||||||
tmp_out.str("");
|
tmp_out.str("");
|
||||||
if (lag > 0)
|
if (lag > 0)
|
||||||
|
@ -1250,7 +1453,7 @@ public:
|
||||||
Stack.pop();
|
Stack.pop();
|
||||||
if (compute)
|
if (compute)
|
||||||
{
|
{
|
||||||
int derivOrder = nearbyint(Stackf.top());
|
int derivOrder = int(nearbyint(Stackf.top()));
|
||||||
Stackf.pop();
|
Stackf.pop();
|
||||||
if (fabs(v1f) < NEAR_ZERO && v2f > 0
|
if (fabs(v1f) < NEAR_ZERO && v2f > 0
|
||||||
&& derivOrder > v2f
|
&& derivOrder > v2f
|
||||||
|
@ -1570,7 +1773,11 @@ public:
|
||||||
}
|
}
|
||||||
tmp_out.str("");
|
tmp_out.str("");
|
||||||
tmp_out << function_name << "(";
|
tmp_out << function_name << "(";
|
||||||
|
#ifndef _MSC_VER
|
||||||
string ss[nb_input_arguments];
|
string ss[nb_input_arguments];
|
||||||
|
#else
|
||||||
|
vector<string> ss(nb_input_arguments);
|
||||||
|
#endif
|
||||||
for (unsigned int i = 0; i < nb_input_arguments; i++)
|
for (unsigned int i = 0; i < nb_input_arguments; i++)
|
||||||
{
|
{
|
||||||
ss[nb_input_arguments-i-1] = Stack.top();
|
ss[nb_input_arguments-i-1] = Stack.top();
|
||||||
|
@ -1624,7 +1831,11 @@ public:
|
||||||
tmp_out.str("");
|
tmp_out.str("");
|
||||||
tmp_out << function_name << "(";
|
tmp_out << function_name << "(";
|
||||||
tmp_out << arg_func_name.c_str() << ", " << fc->get_row() << ", {";
|
tmp_out << arg_func_name.c_str() << ", " << fc->get_row() << ", {";
|
||||||
|
#ifndef _MSC_VER
|
||||||
string ss[nb_add_input_arguments];
|
string ss[nb_add_input_arguments];
|
||||||
|
#else
|
||||||
|
vector<string> ss(nb_input_arguments);
|
||||||
|
#endif
|
||||||
for (unsigned int i = 0; i < nb_add_input_arguments; i++)
|
for (unsigned int i = 0; i < nb_add_input_arguments; i++)
|
||||||
{
|
{
|
||||||
ss[nb_add_input_arguments-i-1] = Stack.top();
|
ss[nb_add_input_arguments-i-1] = Stack.top();
|
||||||
|
@ -1655,7 +1866,11 @@ public:
|
||||||
}
|
}
|
||||||
tmp_out.str("");
|
tmp_out.str("");
|
||||||
tmp_out << function_name << "(";
|
tmp_out << function_name << "(";
|
||||||
|
#ifndef _MSC_VER
|
||||||
string ss[nb_input_arguments];
|
string ss[nb_input_arguments];
|
||||||
|
#else
|
||||||
|
vector<string> ss(nb_input_arguments);
|
||||||
|
#endif
|
||||||
for (unsigned int i = 0; i < nb_input_arguments; i++)
|
for (unsigned int i = 0; i < nb_input_arguments; i++)
|
||||||
{
|
{
|
||||||
ss[nb_input_arguments-i-1] = Stack.top();
|
ss[nb_input_arguments-i-1] = Stack.top();
|
||||||
|
@ -1708,7 +1923,11 @@ public:
|
||||||
tmp_out.str("");
|
tmp_out.str("");
|
||||||
tmp_out << function_name << "(";
|
tmp_out << function_name << "(";
|
||||||
tmp_out << arg_func_name.c_str() << ", " << fc->get_row() << ", " << fc->get_col() << ", {";
|
tmp_out << arg_func_name.c_str() << ", " << fc->get_row() << ", " << fc->get_col() << ", {";
|
||||||
|
#ifndef _MSC_VER
|
||||||
string ss[nb_add_input_arguments];
|
string ss[nb_add_input_arguments];
|
||||||
|
#else
|
||||||
|
vector<string> ss(nb_input_arguments);
|
||||||
|
#endif
|
||||||
for (unsigned int i = 0; i < nb_add_input_arguments; i++)
|
for (unsigned int i = 0; i < nb_add_input_arguments; i++)
|
||||||
{
|
{
|
||||||
ss[nb_add_input_arguments-i-1] = Stack.top();
|
ss[nb_add_input_arguments-i-1] = Stack.top();
|
||||||
|
@ -1739,7 +1958,11 @@ public:
|
||||||
}
|
}
|
||||||
tmp_out.str("");
|
tmp_out.str("");
|
||||||
tmp_out << function_name << "(";
|
tmp_out << function_name << "(";
|
||||||
|
#ifndef _MSC_VER
|
||||||
string ss[nb_input_arguments];
|
string ss[nb_input_arguments];
|
||||||
|
#else
|
||||||
|
vector<string> ss(nb_input_arguments);
|
||||||
|
#endif
|
||||||
for (unsigned int i = 0; i < nb_input_arguments; i++)
|
for (unsigned int i = 0; i < nb_input_arguments; i++)
|
||||||
{
|
{
|
||||||
ss[nb_input_arguments-i-1] = Stack.top();
|
ss[nb_input_arguments-i-1] = Stack.top();
|
||||||
|
@ -1965,7 +2188,7 @@ public:
|
||||||
it_code++;
|
it_code++;
|
||||||
}
|
}
|
||||||
#ifdef DEBUG
|
#ifdef DEBUG
|
||||||
mexPrintf("print_expression end\n"); mexEvalString("drawnow;");
|
mexPrintf("print_expression end tmp_out.str().c_str()=%s\n", tmp_out.str().c_str()); mexEvalString("drawnow;");
|
||||||
#endif
|
#endif
|
||||||
it_code_ret = it_code;
|
it_code_ret = it_code;
|
||||||
return (tmp_out.str());
|
return (tmp_out.str());
|
||||||
|
|
File diff suppressed because it is too large
Load Diff
|
@ -0,0 +1,97 @@
|
||||||
|
/*
|
||||||
|
* Copyright (C) 2007-2012 Dynare Team
|
||||||
|
*
|
||||||
|
* This file is part of Dynare.
|
||||||
|
*
|
||||||
|
* Dynare is free software: you can redistribute it and/or modify
|
||||||
|
* it under the terms of the GNU General Public License as published by
|
||||||
|
* the Free Software Foundation, either version 3 of the License, or
|
||||||
|
* (at your option) any later version.
|
||||||
|
*
|
||||||
|
* Dynare is distributed in the hope that it will be useful,
|
||||||
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||||
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||||
|
* GNU General Public License for more details.
|
||||||
|
*
|
||||||
|
* You should have received a copy of the GNU General Public License
|
||||||
|
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
||||||
|
*/
|
||||||
|
|
||||||
|
#ifndef EVALUATE_HH_INCLUDED
|
||||||
|
#define EVALUATE_HH_INCLUDED
|
||||||
|
|
||||||
|
#include <stack>
|
||||||
|
#include <vector>
|
||||||
|
#include <string>
|
||||||
|
#include <cmath>
|
||||||
|
#define BYTE_CODE
|
||||||
|
#include "CodeInterpreter.hh"
|
||||||
|
#ifdef LINBCG
|
||||||
|
# include "linbcg.hh"
|
||||||
|
#endif
|
||||||
|
#ifndef DEBUG_EX
|
||||||
|
# include <dynmex.h>
|
||||||
|
#else
|
||||||
|
# include "mex_interface.hh"
|
||||||
|
#endif
|
||||||
|
#include "ErrorHandling.hh"
|
||||||
|
|
||||||
|
#define pow_ pow
|
||||||
|
|
||||||
|
class Evaluate : public ErrorMsg
|
||||||
|
{
|
||||||
|
private:
|
||||||
|
unsigned int EQN_dvar1, EQN_dvar2, EQN_dvar3;
|
||||||
|
int EQN_lag1, EQN_lag2, EQN_lag3;
|
||||||
|
protected:
|
||||||
|
mxArray *GlobalTemporaryTerms;
|
||||||
|
it_code_type start_code, end_code;
|
||||||
|
double pow1(double a, double b);
|
||||||
|
double divide(double a, double b);
|
||||||
|
double log1(double a);
|
||||||
|
double log10_1(double a);
|
||||||
|
void evaluate_over_periods(const bool forward);
|
||||||
|
void solve_simple_one_periods();
|
||||||
|
void solve_simple_over_periods(const bool forward);
|
||||||
|
void compute_block_time(const int Per_u_, const bool evaluate, const bool no_derivatives);
|
||||||
|
code_liste_type code_liste;
|
||||||
|
it_code_type it_code;
|
||||||
|
int Block_Count, Per_u_, Per_y_;
|
||||||
|
int it_;
|
||||||
|
int maxit_, size_of_direction;
|
||||||
|
double *direction;
|
||||||
|
double solve_tolf;
|
||||||
|
bool GaussSeidel;
|
||||||
|
map<pair<pair<int, int>, int>, int> IM_i;
|
||||||
|
int equation, derivative_equation, derivative_variable;
|
||||||
|
string filename;
|
||||||
|
int stack_solve_algo, solve_algo;
|
||||||
|
bool global_temporary_terms;
|
||||||
|
bool print, print_error;
|
||||||
|
double res1, res2, max_res;
|
||||||
|
int max_res_idx;
|
||||||
|
vector<Block_contain_type> Block_Contain;
|
||||||
|
|
||||||
|
int size;
|
||||||
|
int *index_vara;
|
||||||
|
|
||||||
|
bool print_it, forward;
|
||||||
|
int minimal_solving_periods;
|
||||||
|
int type, block_num, symbol_table_endo_nbr, Block_List_Max_Lag, Block_List_Max_Lead, u_count_int, block;
|
||||||
|
string file_name, bin_base_name;
|
||||||
|
bool Gaussian_Elimination, is_linear;
|
||||||
|
public:
|
||||||
|
bool steady_state;
|
||||||
|
Evaluate();
|
||||||
|
Evaluate(const int y_size_arg, const int y_kmin_arg, const int y_kmax_arg, const bool print_it_arg, const bool steady_state_arg, const int periods_arg, const int minimal_solving_periods_arg);
|
||||||
|
//typedef void (Interpreter::*InterfpreterMemFn)(const int block_num, const int size, const bool steady_state, int it);
|
||||||
|
void set_block(const int size_arg, const int type_arg, string file_name_arg, string bin_base_name_arg, const int block_num_arg,
|
||||||
|
const bool is_linear_arg, const int symbol_table_endo_nbr_arg, const int Block_List_Max_Lag_arg, const int Block_List_Max_Lead_arg, const int u_count_int_arg, const int block_arg);
|
||||||
|
void evaluate_complete(const bool no_derivatives);
|
||||||
|
bool compute_complete(const bool no_derivatives, double &res1, double &res2, double &max_res, int &max_res_idx);
|
||||||
|
void compute_complete_2b(const bool no_derivatives, double *_res1, double *_res2, double *_max_res, int *_max_res_idx);
|
||||||
|
|
||||||
|
bool compute_complete(double lambda, double *crit);
|
||||||
|
};
|
||||||
|
|
||||||
|
#endif
|
File diff suppressed because it is too large
Load Diff
|
@ -27,6 +27,7 @@
|
||||||
#define BYTE_CODE
|
#define BYTE_CODE
|
||||||
#include "CodeInterpreter.hh"
|
#include "CodeInterpreter.hh"
|
||||||
#include "SparseMatrix.hh"
|
#include "SparseMatrix.hh"
|
||||||
|
#include "Evaluate.hh"
|
||||||
#ifdef LINBCG
|
#ifdef LINBCG
|
||||||
# include "linbcg.hh"
|
# include "linbcg.hh"
|
||||||
#endif
|
#endif
|
||||||
|
@ -40,50 +41,29 @@
|
||||||
|
|
||||||
using namespace std;
|
using namespace std;
|
||||||
|
|
||||||
#define pow_ pow
|
|
||||||
|
|
||||||
class Interpreter : public SparseMatrix
|
class Interpreter : public dynSparseMatrix
|
||||||
{
|
{
|
||||||
private:
|
private:
|
||||||
unsigned int EQN_dvar1, EQN_dvar2, EQN_dvar3;
|
|
||||||
int EQN_lag1, EQN_lag2, EQN_lag3;
|
|
||||||
mxArray *GlobalTemporaryTerms;
|
|
||||||
protected:
|
protected:
|
||||||
double pow1(double a, double b);
|
void evaluate_a_block();
|
||||||
double divide(double a, double b);
|
int simulate_a_block();
|
||||||
double log1(double a);
|
void print_a_block();
|
||||||
double log10_1(double a);
|
|
||||||
void compute_block_time(int Per_u_, bool evaluate, int block_num, int size, bool steady_state);
|
|
||||||
void evaluate_a_block(const int size, const int type, string bin_basename, bool steady_state, int block_num,
|
|
||||||
const bool is_linear = false, const int symbol_table_endo_nbr = 0, const int Block_List_Max_Lag = 0, const int Block_List_Max_Lead = 0, const int u_count_int = 0, int block = -1);
|
|
||||||
int simulate_a_block(const int size, const int type, string file_name, string bin_basename, bool Gaussian_Elimination, bool steady_state, bool print_it, int block_num,
|
|
||||||
const bool is_linear = false, const int symbol_table_endo_nbr = 0, const int Block_List_Max_Lag = 0, const int Block_List_Max_Lead = 0, const int u_count_int = 0);
|
|
||||||
void print_a_block(const int size, const int type, string bin_basename, bool steady_state, int block_num,
|
|
||||||
const bool is_linear, const int symbol_table_endo_nbr, const int Block_List_Max_Lag,
|
|
||||||
const int Block_List_Max_Lead, const int u_count_int, int block);
|
|
||||||
void SingularDisplay(int Per_u_, bool evaluate, int Block_Count, int size, bool steady_state, it_code_type begining);
|
|
||||||
vector<Block_contain_type> Block_Contain;
|
|
||||||
code_liste_type code_liste;
|
|
||||||
it_code_type it_code;
|
|
||||||
int Block_Count, Per_u_, Per_y_;
|
|
||||||
int it_, maxit_, size_of_direction;
|
|
||||||
double solve_tolf;
|
|
||||||
bool GaussSeidel;
|
|
||||||
map<pair<pair<int, int>, int>, int> IM_i;
|
|
||||||
int equation, derivative_equation, derivative_variable;
|
|
||||||
string filename;
|
|
||||||
int minimal_solving_periods;
|
|
||||||
int stack_solve_algo, solve_algo;
|
|
||||||
bool global_temporary_terms;
|
|
||||||
bool print, print_error;
|
|
||||||
public:
|
public:
|
||||||
~Interpreter();
|
~Interpreter();
|
||||||
Interpreter(double *params_arg, double *y_arg, double *ya_arg, double *x_arg, double *steady_y_arg, double *steady_x_arg,
|
Interpreter(double *params_arg, double *y_arg, double *ya_arg, double *x_arg, double *steady_y_arg, double *steady_x_arg,
|
||||||
double *direction_arg, int y_size_arg, int nb_row_x_arg,
|
double *direction_arg, size_t y_size_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,
|
size_t nb_row_x_arg, size_t nb_row_xd_arg, int periods_arg, int y_kmin_arg, int y_kmax_arg,
|
||||||
double slowc_arg, int y_decal_arg, double markowitz_c_arg, string &filename_arg, int minimal_solving_periods_arg, int stack_solve_algo_arg, int solve_algo_arg,
|
int maxit_arg_, double solve_tolf_arg, size_t size_of_direction_arg, double slowc_arg, int y_decal_arg, double markowitz_c_arg,
|
||||||
bool global_temporary_terms_arg, bool print_arg, bool print_error_arg, mxArray *GlobalTemporaryTerms_arg);
|
string &filename_arg, int minimal_solving_periods_arg, int stack_solve_algo_arg, int solve_algo_arg,
|
||||||
bool compute_blocks(string file_name, string bin_basename, bool steady_state, bool evaluate, int block, int &nb_blocks, bool print_it);
|
bool global_temporary_terms_arg, bool print_arg, bool print_error_arg, mxArray *GlobalTemporaryTerms_arg,
|
||||||
|
bool steady_state_arg, bool print_it_arg
|
||||||
|
#ifdef CUDA
|
||||||
|
, const int CUDA_device, cublasHandle_t cublas_handle_arg, cusparseHandle_t cusparse_handle_arg, cusparseMatDescr_t descr_arg
|
||||||
|
#endif
|
||||||
|
);
|
||||||
|
bool compute_blocks(string file_name, string bin_basename, bool evaluate, int block, int &nb_blocks);
|
||||||
|
|
||||||
inline mxArray *
|
inline mxArray *
|
||||||
get_jacob(int block_num)
|
get_jacob(int block_num)
|
||||||
{
|
{
|
||||||
|
|
File diff suppressed because it is too large
Load Diff
|
@ -19,21 +19,62 @@
|
||||||
|
|
||||||
#ifndef SPARSEMATRIX_HH_INCLUDED
|
#ifndef SPARSEMATRIX_HH_INCLUDED
|
||||||
#define SPARSEMATRIX_HH_INCLUDED
|
#define SPARSEMATRIX_HH_INCLUDED
|
||||||
|
#define PRINTF_ printf
|
||||||
|
|
||||||
#include <fstream>
|
|
||||||
#include <stack>
|
#include <stack>
|
||||||
#include <cmath>
|
#include <cmath>
|
||||||
#include <map>
|
#include <map>
|
||||||
#include <ctime>
|
#include <ctime>
|
||||||
|
#include "dynblas.h"
|
||||||
|
#if !(defined _MSC_VER)
|
||||||
|
#include "dynumfpack.h"
|
||||||
|
#endif
|
||||||
|
|
||||||
#ifdef OCTAVE_MEX_FILE
|
#ifdef CUDA
|
||||||
# define CHAR_LENGTH 1
|
#include "cuda.h"
|
||||||
#else
|
#include "cuda_runtime_api.h"
|
||||||
# define CHAR_LENGTH 2
|
#include "cublas_v2.h"
|
||||||
|
#include "cusparse_v2.h"
|
||||||
#endif
|
#endif
|
||||||
|
|
||||||
#include "Mem_Mngr.hh"
|
#include "Mem_Mngr.hh"
|
||||||
#include "ErrorHandling.hh"
|
#include "ErrorHandling.hh"
|
||||||
|
//#include "Interpreter.hh"
|
||||||
|
#include "Evaluate.hh"
|
||||||
|
|
||||||
|
#define cudaChk(x, y) \
|
||||||
|
{ \
|
||||||
|
cudaError_t cuda_error = x; \
|
||||||
|
if (cuda_error != cudaSuccess) \
|
||||||
|
{ \
|
||||||
|
ostringstream tmp; \
|
||||||
|
tmp << y; \
|
||||||
|
throw FatalExceptionHandling(tmp.str()); \
|
||||||
|
} \
|
||||||
|
};
|
||||||
|
|
||||||
|
#define cusparseChk(x, y) \
|
||||||
|
{ \
|
||||||
|
cusparseStatus_t cusparse_status = x; \
|
||||||
|
if (cusparse_status != CUSPARSE_STATUS_SUCCESS) \
|
||||||
|
{ \
|
||||||
|
ostringstream tmp; \
|
||||||
|
tmp << y; \
|
||||||
|
throw FatalExceptionHandling(tmp.str()); \
|
||||||
|
} \
|
||||||
|
};
|
||||||
|
|
||||||
|
#define cublasChk(x, y) \
|
||||||
|
{ \
|
||||||
|
cublasStatus_t cublas_status = x; \
|
||||||
|
if (cublas_status != CUBLAS_STATUS_SUCCESS) \
|
||||||
|
{ \
|
||||||
|
ostringstream tmp; \
|
||||||
|
tmp << y; \
|
||||||
|
throw FatalExceptionHandling(tmp.str()); \
|
||||||
|
} \
|
||||||
|
};
|
||||||
|
|
||||||
#define NEW_ALLOC
|
#define NEW_ALLOC
|
||||||
#define MARKOVITZ
|
#define MARKOVITZ
|
||||||
|
|
||||||
|
@ -53,41 +94,76 @@ const int IFLDZ = 4;
|
||||||
const int IFMUL = 5;
|
const int IFMUL = 5;
|
||||||
const int IFSTP = 6;
|
const int IFSTP = 6;
|
||||||
const int IFADD = 7;
|
const int IFADD = 7;
|
||||||
const double eps = 1e-10;
|
const double eps = 1e-15;
|
||||||
const double very_big = 1e24;
|
const double very_big = 1e24;
|
||||||
const int alt_symbolic_count_max = 1;
|
const int alt_symbolic_count_max = 1;
|
||||||
const double mem_increasing_factor = 1.1;
|
const double mem_increasing_factor = 1.1;
|
||||||
|
|
||||||
class SparseMatrix : public ErrorMsg
|
|
||||||
|
|
||||||
|
class dynSparseMatrix : public Evaluate
|
||||||
{
|
{
|
||||||
public:
|
public:
|
||||||
SparseMatrix();
|
#if (defined _MSC_VER)
|
||||||
void Simulate_Newton_Two_Boundaries(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it, bool cvg, int &iter, int minimal_solving_periods, int stack_solve_algo, unsigned int endo_name_length, char *P_endo_names) /*throw(ErrorHandlingException)*/;
|
typedef int64_t SuiteSparse_long;
|
||||||
bool Simulate_Newton_One_Boundary(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, int stack_solve_algo, int solve_algo);
|
#endif
|
||||||
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);
|
dynSparseMatrix();
|
||||||
|
dynSparseMatrix(const int y_size_arg, const int y_kmin_arg, const int y_kmax_arg, const bool print_it_arg, const bool steady_state_arg, const int periods_arg, const int minimal_solving_periods_arg
|
||||||
|
#ifdef CUDA
|
||||||
|
,const int CUDA_device_arg, cublasHandle_t cublas_handle_arg, cusparseHandle_t cusparse_handle_arg, cusparseMatDescr_t descr_arg
|
||||||
|
#endif
|
||||||
|
);
|
||||||
|
void Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin, int y_kmax, int Size, int periods, bool cvg, int minimal_solving_periods, int stack_solve_algo, unsigned int endo_name_length, char *P_endo_names);
|
||||||
|
void Simulate_Newton_One_Boundary(bool forward);
|
||||||
void fixe_u(double **u, int u_count_int, int max_lag_plus_max_lead_plus_1);
|
void fixe_u(double **u, int u_count_int, int max_lag_plus_max_lead_plus_1);
|
||||||
void Read_SparseMatrix(string file_name, const int Size, int periods, int y_kmin, int y_kmax, bool steady_state, bool two_boundaries, int stack_solve_algo, int solve_algo);
|
void Read_SparseMatrix(string file_name, const int Size, int periods, int y_kmin, int y_kmax, bool two_boundaries, int stack_solve_algo, int solve_algo);
|
||||||
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);
|
||||||
void Singular_display(int block, int Size, bool steady_state, it_code_type it_code);
|
void Singular_display(int block, int Size);
|
||||||
double g0, gp0, glambda2, try_at_iteration;
|
void End_Solver();
|
||||||
|
double g0, gp0, glambda2;
|
||||||
|
int try_at_iteration;
|
||||||
|
|
||||||
private:
|
private:
|
||||||
void Init_GE(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM);
|
void Init_GE(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM);
|
||||||
void Init_Matlab_Sparse(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM, mxArray *A_m, mxArray *b_m, mxArray *x0_m);
|
void Init_Matlab_Sparse(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM, mxArray *A_m, mxArray *b_m, mxArray *x0_m);
|
||||||
|
void Init_UMFPACK_Sparse(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM, SuiteSparse_long **Ap, SuiteSparse_long **Ai, double **Ax, double **b, mxArray *x0_m);
|
||||||
|
#ifdef CUDA
|
||||||
|
void Init_CUDA_Sparse(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM, int **Ap, int **Ai, double **Ax, int **Ap_tild, int **Ai_tild, double **A_tild, double **b, double **x0, mxArray *x0_m, int *nnz, int *nnz_tild, int preconditioner);
|
||||||
|
#endif
|
||||||
void Init_Matlab_Sparse_Simple(int Size, map<pair<pair<int, int>, int>, int> &IM, mxArray *A_m, mxArray *b_m, bool &zero_solution, mxArray *x0_m);
|
void Init_Matlab_Sparse_Simple(int Size, map<pair<pair<int, int>, int>, int> &IM, mxArray *A_m, mxArray *b_m, bool &zero_solution, mxArray *x0_m);
|
||||||
|
void Init_UMFPACK_Sparse_Simple(int Size, map<pair<pair<int, int>, int>, int> &IM, SuiteSparse_long **Ap, SuiteSparse_long **Ai, double **Ax, double **b, bool &zero_solution, mxArray *x0_m);
|
||||||
|
void Init_CUDA_Sparse_Simple(int Size, map<pair<pair<int, int>, int>, int> &IM, SuiteSparse_long **Ap, SuiteSparse_long **Ai, double **Ax, double **b, double **x0, bool &zero_solution, mxArray *x0_m);
|
||||||
void Simple_Init(int Size, std::map<std::pair<std::pair<int, int>, int>, int> &IM, bool &zero_solution);
|
void Simple_Init(int Size, std::map<std::pair<std::pair<int, int>, int>, int> &IM, bool &zero_solution);
|
||||||
void End_GE(int Size);
|
void End_GE(int Size);
|
||||||
|
bool mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb, double *fc);
|
||||||
|
bool golden(double ax, double bx, double cx, double tol, double solve_tolf, double *xmin);
|
||||||
void Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bool symbolic, int Block_number);
|
void Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bool symbolic, int Block_number);
|
||||||
bool Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, bool steady_state, int it_);
|
bool Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, int it_);
|
||||||
void Solve_Matlab_Relaxation(mxArray *A_m, mxArray *b_m, unsigned int Size, double slowc_l, bool is_two_boundaries, int it_);
|
void Solve_Matlab_Relaxation(mxArray *A_m, mxArray *b_m, unsigned int Size, double slowc_l, bool is_two_boundaries, int it_);
|
||||||
void Solve_Matlab_LU_UMFPack(mxArray *A_m, mxArray *b_m, int Size, double slowc_l, bool is_two_boundaries, int it_);
|
void Solve_Matlab_LU_UMFPack(mxArray *A_m, mxArray *b_m, int Size, double slowc_l, bool is_two_boundaries, int it_);
|
||||||
void Solve_Matlab_GMRES(mxArray *A_m, mxArray *b_m, int Size, double slowc, int block, bool is_two_boundaries, int it_, bool steady_state, mxArray *x0_m);
|
void Solve_LU_UMFPack(mxArray *A_m, mxArray *b_m, int Size, double slowc_l, bool is_two_boundaries, int it_);
|
||||||
void Solve_Matlab_BiCGStab(mxArray *A_m, mxArray *b_m, int Size, double slowc, int block, bool is_two_boundaries, int it_, mxArray *x0_m, bool steady_state);
|
void Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, double *Ax, double *b, int n, int Size, double slowc_l, bool is_two_boundaries, int it_);
|
||||||
|
void End_Matlab_LU_UMFPack();
|
||||||
|
#ifdef CUDA
|
||||||
|
void Solve_CUDA_BiCGStab_Free(double* tmp_vect_host, double* p, double* r, double* v, double* s, double* t, double* y_, double* z, double* tmp_,
|
||||||
|
int* Ai, double* Ax, int* Ap, double* x0, double* b, double* A_tild, int* A_tild_i, int* A_tild_p,
|
||||||
|
cusparseSolveAnalysisInfo_t infoL, cusparseSolveAnalysisInfo_t infoU,
|
||||||
|
cusparseMatDescr_t descrL, cusparseMatDescr_t descrU, int preconditioner);
|
||||||
|
int Solve_CUDA_BiCGStab(int *Ap, int *Ai, double *Ax, int *Ap_tild, int *Ai_tild, double *A_tild, double *b, double *x0, int n, int Size, double slowc_l, bool is_two_boundaries, int it_, int nnz, int nnz_tild, int preconditioner, int max_iterations, int block);
|
||||||
|
#endif
|
||||||
|
void Solve_Matlab_GMRES(mxArray *A_m, mxArray *b_m, int Size, double slowc, int block, bool is_two_boundaries, int it_, mxArray *x0_m);
|
||||||
|
void Solve_Matlab_BiCGStab(mxArray *A_m, mxArray *b_m, int Size, double slowc, int block, bool is_two_boundaries, int it_, mxArray *x0_m, int precond);
|
||||||
|
void Check_and_Correct_Previous_Iteration(int block_num, int y_size, int size, double crit_opt_old);
|
||||||
|
bool Simulate_One_Boundary(int blck, int y_size, int y_kmin, int y_kmax, int Size, bool cvg);
|
||||||
|
bool solve_linear(const int block_num, const int y_size, const int y_kmin, const int y_kmax, const int size, const int iter);
|
||||||
|
void solve_non_linear(const int block_num, const int y_size, const int y_kmin, const int y_kmax, const int size);
|
||||||
|
string preconditioner_print_out(string s, int preconditioner);
|
||||||
bool compare(int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, long int nop4, int Size
|
bool compare(int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, long int nop4, int Size
|
||||||
#ifdef PROFILER
|
#ifdef PROFILER
|
||||||
, long int *ndiv, long int *nsub
|
, long int *ndiv, long int *nsub
|
||||||
#endif
|
#endif
|
||||||
);
|
);
|
||||||
|
void Grad_f_product(int n, mxArray *b_m, double* vectr, mxArray *A_m, SuiteSparse_long *Ap, SuiteSparse_long *Ai, double* Ax, double *b);
|
||||||
void Insert(const int r, const int c, const int u_index, const int lag_index);
|
void Insert(const int r, const int c, const int u_index, const int lag_index);
|
||||||
void Delete(const int r, const int c);
|
void Delete(const int r, const int c);
|
||||||
int At_Row(int r, NonZeroElem **first);
|
int At_Row(int r, NonZeroElem **first);
|
||||||
|
@ -102,7 +178,8 @@ private:
|
||||||
void Delete_u(int pos);
|
void Delete_u(int pos);
|
||||||
void Clear_u();
|
void Clear_u();
|
||||||
void Print_u();
|
void Print_u();
|
||||||
void CheckIt(int y_size, int y_kmin, int y_kmax, int Size, int periods, int iter);
|
void *Symbolic, *Numeric ;
|
||||||
|
void CheckIt(int y_size, int y_kmin, int y_kmax, int Size, int periods);
|
||||||
void Check_the_Solution(int periods, int y_kmin, int y_kmax, int Size, double *u, int *pivot, int *b);
|
void Check_the_Solution(int periods, int y_kmin, int y_kmax, int Size, double *u, int *pivot, int *b);
|
||||||
int complete(int beg_t, int Size, int periods, int *b);
|
int complete(int beg_t, int Size, int periods, int *b);
|
||||||
void bksub(int tbreak, int last_period, int Size, double slowc_l
|
void bksub(int tbreak, int last_period, int Size, double slowc_l
|
||||||
|
@ -118,12 +195,18 @@ private:
|
||||||
mxArray *Sparse_substract_SA_SB(mxArray *A_m, mxArray *B_m);
|
mxArray *Sparse_substract_SA_SB(mxArray *A_m, mxArray *B_m);
|
||||||
mxArray *Sparse_substract_A_SB(mxArray *A_m, mxArray *B_m);
|
mxArray *Sparse_substract_A_SB(mxArray *A_m, mxArray *B_m);
|
||||||
mxArray *substract_A_B(mxArray *A_m, mxArray *B_m);
|
mxArray *substract_A_B(mxArray *A_m, mxArray *B_m);
|
||||||
|
#ifdef CUDA
|
||||||
|
int CUDA_device;
|
||||||
|
cublasHandle_t cublas_handle;
|
||||||
|
cusparseHandle_t cusparse_handle;
|
||||||
|
cusparseMatDescr_t CUDA_descr;
|
||||||
|
#endif
|
||||||
|
protected:
|
||||||
stack<double> Stack;
|
stack<double> Stack;
|
||||||
int nb_prologue_table_u, nb_first_table_u, nb_middle_table_u, nb_last_table_u;
|
int nb_prologue_table_u, nb_first_table_u, nb_middle_table_u, nb_last_table_u;
|
||||||
int nb_prologue_table_y, nb_first_table_y, nb_middle_table_y, nb_last_table_y;
|
int nb_prologue_table_y, nb_first_table_y, nb_middle_table_y, nb_last_table_y;
|
||||||
int middle_count_loop;
|
int middle_count_loop;
|
||||||
char type;
|
//char type;
|
||||||
fstream SaveCode;
|
fstream SaveCode;
|
||||||
string filename;
|
string filename;
|
||||||
int max_u, min_u;
|
int max_u, min_u;
|
||||||
|
@ -154,19 +237,19 @@ protected:
|
||||||
int u_count_alloc, u_count_alloc_save;
|
int u_count_alloc, u_count_alloc_save;
|
||||||
vector<double *> jac;
|
vector<double *> jac;
|
||||||
double *jcb;
|
double *jcb;
|
||||||
double res1, res2, max_res;
|
|
||||||
int max_res_idx;
|
|
||||||
double slowc, slowc_save, prev_slowc_save, markowitz_c;
|
double slowc, slowc_save, prev_slowc_save, markowitz_c;
|
||||||
int y_decal;
|
int y_decal;
|
||||||
int *index_vara, *index_equa;
|
int *index_equa;
|
||||||
int u_count, tbreak_g;
|
int u_count, tbreak_g;
|
||||||
int iter;
|
int iter;
|
||||||
double *direction;
|
|
||||||
int start_compare;
|
int start_compare;
|
||||||
int restart;
|
int restart;
|
||||||
bool error_not_printed;
|
|
||||||
double g_lambda1, g_lambda2, gp_0;
|
double g_lambda1, g_lambda2, gp_0;
|
||||||
double lu_inc_tol;
|
double lu_inc_tol;
|
||||||
|
//private:
|
||||||
|
SuiteSparse_long *Ap_save, *Ai_save;
|
||||||
|
double *Ax_save, *b_save;
|
||||||
|
mxArray *A_m_save, *b_m_save;
|
||||||
};
|
};
|
||||||
|
|
||||||
#endif
|
#endif
|
||||||
|
|
|
@ -0,0 +1,121 @@
|
||||||
|
/*
|
||||||
|
* Copyright (C) 2007-2012 Dynare Team
|
||||||
|
*
|
||||||
|
* This file is part of Dynare.
|
||||||
|
*
|
||||||
|
* Dynare is free software: you can redistribute it and/or modify
|
||||||
|
* it under the terms of the GNU General Public License as published by
|
||||||
|
* the Free Software Foundation, either version 3 of the License, or
|
||||||
|
* (at your option) any later version.
|
||||||
|
*
|
||||||
|
* Dynare is distributed in the hope that it will be useful,
|
||||||
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||||
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||||
|
* GNU General Public License for more details.
|
||||||
|
*
|
||||||
|
* You should have received a copy of the GNU General Public License
|
||||||
|
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
||||||
|
*/
|
||||||
|
|
||||||
|
#ifndef SPARMATRIX_KERNEL
|
||||||
|
#define SPARMATRIX_KERNEL
|
||||||
|
|
||||||
|
// Kernel definition of vector division
|
||||||
|
__global__ void
|
||||||
|
VecDiv(double* A, double* B, double* C, int n)
|
||||||
|
{
|
||||||
|
int i = blockIdx.x * 1024 + threadIdx.x;
|
||||||
|
if (i < n)
|
||||||
|
C[i] = (B[i] != 0.0 ? A[i] / B[i] : A[i]);
|
||||||
|
}
|
||||||
|
|
||||||
|
__global__ void
|
||||||
|
VecAdd(double* res, double* r, double alpha, double* x, int n)
|
||||||
|
{
|
||||||
|
int i = blockIdx.x * 1024 + threadIdx.x;
|
||||||
|
if (i < n)
|
||||||
|
res[i] = r[i] + alpha * x[i];
|
||||||
|
}
|
||||||
|
|
||||||
|
__global__ void
|
||||||
|
VecInc(double* res, double alpha, double* x, int n)
|
||||||
|
{
|
||||||
|
int i = blockIdx.x * 1024 + threadIdx.x;
|
||||||
|
if (i < n)
|
||||||
|
res[i] += alpha * x[i];
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
__global__ void
|
||||||
|
update_x(double* x, double alpha, double* y, double omega, double *z)
|
||||||
|
{
|
||||||
|
int i = threadIdx.x;
|
||||||
|
x[i] += alpha * y[i] + omega * z[i];
|
||||||
|
}
|
||||||
|
|
||||||
|
__global__ void
|
||||||
|
Get_LU_dim(int *n, int* A_tild_i, int *A_tild_p, int *nnz_l, int *nnz_u)
|
||||||
|
{
|
||||||
|
nnz_u[0] = 0;
|
||||||
|
nnz_l[0] = 0;
|
||||||
|
for (int i = 0; i < n[0]; i++)
|
||||||
|
{
|
||||||
|
for (int j = A_tild_p[i]; j < A_tild_p[i+1]; j++)
|
||||||
|
{
|
||||||
|
if (A_tild_i[j] < i)
|
||||||
|
nnz_l[0]++;
|
||||||
|
else if (A_tild_i[j] == i)
|
||||||
|
{
|
||||||
|
nnz_u[0]++;
|
||||||
|
//nnz_l[0]++;
|
||||||
|
}
|
||||||
|
else
|
||||||
|
nnz_u[0]++;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
__global__ void
|
||||||
|
Get_LU1_dim(int* n, int *nnz_l, int *nnz_u)
|
||||||
|
{
|
||||||
|
nnz_u[0] = 3+n[0];
|
||||||
|
nnz_l[0] = 1+n[0];
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
__global__ void
|
||||||
|
Get_L_and_U(int *n, double* A_tild_x, int* A_tild_i, int *A_tild_p, double* Lx, int* Li, int *Lp, double* Ux, int* Ui, int* Up)
|
||||||
|
{
|
||||||
|
int nnz_u = 0, nnz_l = 0;
|
||||||
|
Lp[0] = 0;
|
||||||
|
Up[0] = 0;
|
||||||
|
for (int i = 0; i < n[0]; i++)
|
||||||
|
{
|
||||||
|
for (int j = A_tild_p[i]; j < A_tild_p[i+1]; j++)
|
||||||
|
{
|
||||||
|
if (A_tild_i[j] < i)
|
||||||
|
{
|
||||||
|
Lx[nnz_l] = A_tild_x[j];
|
||||||
|
Li[nnz_l] = A_tild_i[j];
|
||||||
|
nnz_l++;
|
||||||
|
}
|
||||||
|
else if (A_tild_i[j] == i)
|
||||||
|
{
|
||||||
|
Ux[nnz_u] = A_tild_x[j];
|
||||||
|
Lx[nnz_l] = 1.0;
|
||||||
|
Li[nnz_l] = Ui[nnz_u] = A_tild_i[j];
|
||||||
|
nnz_u++;
|
||||||
|
//nnz_l++;
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
Ux[nnz_u] = A_tild_x[j];
|
||||||
|
Ui[nnz_u] = A_tild_i[j];
|
||||||
|
nnz_u++;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
Lp[i+1] = nnz_l;
|
||||||
|
Up[i+1] = nnz_u;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
#endif
|
|
@ -18,12 +18,18 @@
|
||||||
*/
|
*/
|
||||||
#include <cstring>
|
#include <cstring>
|
||||||
#include "Interpreter.hh"
|
#include "Interpreter.hh"
|
||||||
|
#include "ErrorHandling.hh"
|
||||||
|
#include <ctime>
|
||||||
|
#include <math.h>
|
||||||
|
|
||||||
#ifdef DEBUG_EX
|
#ifdef DEBUG_EX
|
||||||
|
|
||||||
using namespace std;
|
using namespace std;
|
||||||
# include <sstream>
|
# include <sstream>
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
string
|
string
|
||||||
Get_Argument(const char *argv)
|
Get_Argument(const char *argv)
|
||||||
{
|
{
|
||||||
|
@ -33,14 +39,17 @@ Get_Argument(const char *argv)
|
||||||
|
|
||||||
#else
|
#else
|
||||||
|
|
||||||
|
void (*prev_fn)(int);
|
||||||
|
|
||||||
|
|
||||||
string
|
string
|
||||||
Get_Argument(const mxArray *prhs)
|
Get_Argument(const mxArray *prhs)
|
||||||
{
|
{
|
||||||
const mxArray *mxa = prhs;
|
const mxArray *mxa = prhs;
|
||||||
int buflen = mxGetM(mxa) * mxGetN(mxa) + 1;
|
mwSize buflen = mwSize(mxGetM(mxa) * mxGetN(mxa) + 1);
|
||||||
char *first_argument;
|
char *first_argument;
|
||||||
first_argument = (char *) mxCalloc(buflen, sizeof(char));
|
first_argument = (char *) mxCalloc(buflen, sizeof(char));
|
||||||
int status = mxGetString(mxa, first_argument, buflen);
|
size_t status = mxGetString(mxa, first_argument, buflen);
|
||||||
if (status != 0)
|
if (status != 0)
|
||||||
mexWarnMsgTxt("Not enough space. The first argument is truncated.");
|
mexWarnMsgTxt("Not enough space. The first argument is truncated.");
|
||||||
string f(first_argument);
|
string f(first_argument);
|
||||||
|
@ -49,6 +58,178 @@ Get_Argument(const mxArray *prhs)
|
||||||
}
|
}
|
||||||
#endif
|
#endif
|
||||||
|
|
||||||
|
|
||||||
|
//#include <windows.h>
|
||||||
|
#include <stdio.h>
|
||||||
|
|
||||||
|
|
||||||
|
#ifdef CUDA
|
||||||
|
int
|
||||||
|
GPU_Test_and_Info(cublasHandle_t *cublas_handle, cusparseHandle_t *cusparse_handle, cusparseMatDescr_t *descr)
|
||||||
|
{
|
||||||
|
cudaDeviceProp deviceProp;
|
||||||
|
int device_count, device, version, version_max = 0;
|
||||||
|
cublasStatus_t cublas_status;
|
||||||
|
cudaError_t cuda_error;
|
||||||
|
*descr=0;
|
||||||
|
|
||||||
|
/* ask cuda how many devices it can find */
|
||||||
|
cudaGetDeviceCount(&device_count);
|
||||||
|
if (device_count < 1)
|
||||||
|
{
|
||||||
|
/* if it couldn't find any fail out */
|
||||||
|
ostringstream tmp;
|
||||||
|
tmp << " Unable to find a CUDA device. Unable to implement CUDA solvers\n";
|
||||||
|
throw FatalExceptionHandling(tmp.str());
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
mexPrintf("-----------------------------------------\n");
|
||||||
|
for (int i = 0; i < device_count; i++)
|
||||||
|
{
|
||||||
|
cudaSetDevice(i);
|
||||||
|
// Statistics about the GPU device
|
||||||
|
cuda_error = cudaGetDeviceProperties(&deviceProp, i);
|
||||||
|
if (cuda_error != cudaSuccess)
|
||||||
|
{
|
||||||
|
ostringstream tmp;
|
||||||
|
tmp << " bytecode cudaGetDeviceProperties failed\n";
|
||||||
|
throw FatalExceptionHandling(tmp.str());
|
||||||
|
}
|
||||||
|
mexPrintf("> GPU device %d: \"%s\" has:\n - %d Multi-Processors,\n - %d threads per multiprocessor,\n", i, deviceProp.name, deviceProp.multiProcessorCount, deviceProp.maxThreadsPerMultiProcessor);
|
||||||
|
mexEvalString("drawnow;");
|
||||||
|
version = (deviceProp.major * 0x10 + deviceProp.minor);
|
||||||
|
if (version >= version_max)
|
||||||
|
{
|
||||||
|
device = i;
|
||||||
|
version_max = version;
|
||||||
|
}
|
||||||
|
mexPrintf(" - %4.2fMhz clock rate,\n - %2.0fMb of memory,\n - %d.%d compute capabilities.\n", double(deviceProp.clockRate) / (1024 * 1024), double(deviceProp.totalGlobalMem) / (1024 * 1024), deviceProp.major, deviceProp.minor);
|
||||||
|
mexEvalString("drawnow;");
|
||||||
|
}
|
||||||
|
}
|
||||||
|
mexPrintf("> Device %d selected\n", device);
|
||||||
|
mexEvalString("drawnow;");
|
||||||
|
|
||||||
|
cuda_error = cudaSetDevice(device);
|
||||||
|
if (cuda_error != cudaSuccess)
|
||||||
|
{
|
||||||
|
ostringstream tmp;
|
||||||
|
tmp << " bytecode cudaSetDevice failed\n";
|
||||||
|
throw FatalExceptionHandling(tmp.str());
|
||||||
|
}
|
||||||
|
|
||||||
|
if(version_max < 0x11)
|
||||||
|
{
|
||||||
|
ostringstream tmp;
|
||||||
|
tmp << " bytecode requires a minimum CUDA compute 1.1 capability\n";
|
||||||
|
cudaDeviceReset();
|
||||||
|
throw FatalExceptionHandling(tmp.str());
|
||||||
|
}
|
||||||
|
|
||||||
|
// Initialize CuBlas library
|
||||||
|
cublas_status = cublasCreate(cublas_handle);
|
||||||
|
if (cublas_status != CUBLAS_STATUS_SUCCESS)
|
||||||
|
{
|
||||||
|
ostringstream tmp;
|
||||||
|
switch(cublas_status)
|
||||||
|
{
|
||||||
|
case CUBLAS_STATUS_NOT_INITIALIZED:
|
||||||
|
tmp << " the CUBLAS initialization failed.\n";
|
||||||
|
break;
|
||||||
|
case CUBLAS_STATUS_ALLOC_FAILED:
|
||||||
|
tmp << " the resources could not be allocated.\n";
|
||||||
|
break;
|
||||||
|
default:
|
||||||
|
tmp << " unknown error during the initialization of cusparse library.\n";
|
||||||
|
}
|
||||||
|
throw FatalExceptionHandling(tmp.str());
|
||||||
|
}
|
||||||
|
|
||||||
|
// Initialize the CuSparse library
|
||||||
|
cusparseStatus_t cusparse_status;
|
||||||
|
cusparse_status = cusparseCreate(cusparse_handle);
|
||||||
|
if (cusparse_status != CUSPARSE_STATUS_SUCCESS)
|
||||||
|
{
|
||||||
|
ostringstream tmp;
|
||||||
|
switch(cusparse_status)
|
||||||
|
{
|
||||||
|
case CUSPARSE_STATUS_NOT_INITIALIZED:
|
||||||
|
tmp << " the CUDA Runtime initialization failed.\n";
|
||||||
|
break;
|
||||||
|
case CUSPARSE_STATUS_ALLOC_FAILED:
|
||||||
|
tmp << " the resources could not be allocated.\n";
|
||||||
|
break;
|
||||||
|
case CUSPARSE_STATUS_ARCH_MISMATCH:
|
||||||
|
tmp << " the device compute capability (CC) is less than 1.1. The CC of at least 1.1 is required.\n";
|
||||||
|
break;
|
||||||
|
default:
|
||||||
|
tmp << " unknown error during the initialization of cusparse library.\n";
|
||||||
|
}
|
||||||
|
throw FatalExceptionHandling(tmp.str());
|
||||||
|
}
|
||||||
|
|
||||||
|
// Create and setup matrix descriptor
|
||||||
|
cusparse_status = cusparseCreateMatDescr(descr);
|
||||||
|
if (cusparse_status != CUSPARSE_STATUS_SUCCESS)
|
||||||
|
{
|
||||||
|
ostringstream tmp;
|
||||||
|
tmp << " Matrix descriptor initialization failed\n";
|
||||||
|
throw FatalExceptionHandling(tmp.str());
|
||||||
|
}
|
||||||
|
cusparseSetMatType(*descr, CUSPARSE_MATRIX_TYPE_GENERAL);
|
||||||
|
cusparseSetMatIndexBase(*descr, CUSPARSE_INDEX_BASE_ZERO);
|
||||||
|
|
||||||
|
mexPrintf("> Driver version:\n");
|
||||||
|
int cuda_version;
|
||||||
|
cuda_error = cudaDriverGetVersion(&cuda_version);
|
||||||
|
if (cuda_error != cudaSuccess)
|
||||||
|
{
|
||||||
|
ostringstream tmp;
|
||||||
|
tmp << " cudaGetVersion has failed\n";
|
||||||
|
throw FatalExceptionHandling(tmp.str());
|
||||||
|
}
|
||||||
|
mexPrintf(" - CUDA version %5.3f\n", double(cuda_version) / 1000);
|
||||||
|
int cublas_version;
|
||||||
|
cublas_status = cublasGetVersion(*cublas_handle, &cublas_version);
|
||||||
|
if (cublas_status != CUBLAS_STATUS_SUCCESS)
|
||||||
|
{
|
||||||
|
ostringstream tmp;
|
||||||
|
tmp << " cublasGetVersion has failed\n";
|
||||||
|
throw FatalExceptionHandling(tmp.str());
|
||||||
|
}
|
||||||
|
mexPrintf(" - CUBLAS version %5.3f\n", double(cublas_version) / 1000);
|
||||||
|
int cusparse_version;
|
||||||
|
cusparse_status = cusparseGetVersion(*cusparse_handle, &cusparse_version);
|
||||||
|
if (cusparse_status != CUSPARSE_STATUS_SUCCESS)
|
||||||
|
{
|
||||||
|
ostringstream tmp;
|
||||||
|
tmp << " cusparseGetVersion has failed\n";
|
||||||
|
throw FatalExceptionHandling(tmp.str());
|
||||||
|
}
|
||||||
|
mexPrintf(" - CUSPARSE version %5.3f\n", double(cusparse_version) / 1000);
|
||||||
|
mexPrintf("-----------------------------------------\n");
|
||||||
|
return device;
|
||||||
|
}
|
||||||
|
|
||||||
|
void
|
||||||
|
GPU_close(cublasHandle_t cublas_handle, cusparseHandle_t cusparse_handle, cusparseMatDescr_t descr)
|
||||||
|
{
|
||||||
|
cublasChk(cublasDestroy(cublas_handle),"in bytecode cublasDestroy failed\n");
|
||||||
|
cusparseChk(cusparseDestroyMatDescr(descr), "in bytecode cusparseDestroyMatDescr failed\n");
|
||||||
|
cusparseChk(cusparseDestroy(cusparse_handle),"in bytecode cusparseDestroy failed\n");
|
||||||
|
}
|
||||||
|
|
||||||
|
#endif
|
||||||
|
string
|
||||||
|
deblank(string x)
|
||||||
|
{
|
||||||
|
for(int i = 0; i < x.length(); i++)
|
||||||
|
if (x[i] == ' ')
|
||||||
|
x.erase(i--, 1);
|
||||||
|
return x;
|
||||||
|
}
|
||||||
|
|
||||||
void
|
void
|
||||||
Get_Arguments_and_global_variables(int nrhs,
|
Get_Arguments_and_global_variables(int nrhs,
|
||||||
#ifndef DEBUG_EX
|
#ifndef DEBUG_EX
|
||||||
|
@ -57,10 +238,10 @@ Get_Arguments_and_global_variables(int nrhs,
|
||||||
const char *prhs[],
|
const char *prhs[],
|
||||||
#endif
|
#endif
|
||||||
int &count_array_argument,
|
int &count_array_argument,
|
||||||
double *yd[], unsigned int &row_y, unsigned int &col_y,
|
double *yd[], size_t &row_y, size_t &col_y,
|
||||||
double *xd[], unsigned int &row_x, unsigned int &col_x,
|
double *xd[], size_t &row_x, size_t &col_x,
|
||||||
double *params[],
|
double *params[],
|
||||||
double *steady_yd[], unsigned int &steady_row_y, unsigned int &steady_col_y,
|
double *steady_yd[], size_t &steady_row_y, size_t &steady_col_y,
|
||||||
unsigned int &periods,
|
unsigned int &periods,
|
||||||
#ifndef DEBUG_EX
|
#ifndef DEBUG_EX
|
||||||
mxArray *block_structur[],
|
mxArray *block_structur[],
|
||||||
|
@ -69,8 +250,10 @@ Get_Arguments_and_global_variables(int nrhs,
|
||||||
mxArray *M_[], mxArray *oo_[], mxArray *options_[], bool &global_temporary_terms,
|
mxArray *M_[], mxArray *oo_[], mxArray *options_[], bool &global_temporary_terms,
|
||||||
bool &print,
|
bool &print,
|
||||||
bool &print_error,
|
bool &print_error,
|
||||||
mxArray *GlobalTemporaryTerms[])
|
mxArray *GlobalTemporaryTerms[],
|
||||||
|
string *plan_struct_name, string *pfplan_struct_name)
|
||||||
{
|
{
|
||||||
|
size_t pos;
|
||||||
#ifdef DEBUG_EX
|
#ifdef DEBUG_EX
|
||||||
for (int i = 2; i < nrhs; i++)
|
for (int i = 2; i < nrhs; i++)
|
||||||
#else
|
#else
|
||||||
|
@ -101,7 +284,7 @@ Get_Arguments_and_global_variables(int nrhs,
|
||||||
steady_col_y = mxGetN(prhs[i]);
|
steady_col_y = mxGetN(prhs[i]);
|
||||||
break;
|
break;
|
||||||
case 4:
|
case 4:
|
||||||
periods = mxGetScalar(prhs[i]);
|
periods = int(mxGetScalar(prhs[i]));
|
||||||
break;
|
break;
|
||||||
case 5:
|
case 5:
|
||||||
*block_structur = mxDuplicateArray(prhs[i]);
|
*block_structur = mxDuplicateArray(prhs[i]);
|
||||||
|
@ -111,7 +294,7 @@ Get_Arguments_and_global_variables(int nrhs,
|
||||||
*GlobalTemporaryTerms = mxDuplicateArray(prhs[i]);
|
*GlobalTemporaryTerms = mxDuplicateArray(prhs[i]);
|
||||||
break;
|
break;
|
||||||
default:
|
default:
|
||||||
//mexPrintf("Unknown argument count_array_argument=%d\n",count_array_argument);
|
mexPrintf("Unknown argument count_array_argument=%d\n",count_array_argument);
|
||||||
break;
|
break;
|
||||||
}
|
}
|
||||||
count_array_argument++;
|
count_array_argument++;
|
||||||
|
@ -132,16 +315,34 @@ Get_Arguments_and_global_variables(int nrhs,
|
||||||
print_error = false;
|
print_error = false;
|
||||||
else
|
else
|
||||||
{
|
{
|
||||||
int pos = Get_Argument(prhs[i]).find("block");
|
;
|
||||||
if (pos != (int) string::npos)
|
if ((pos = Get_Argument(prhs[i]).find("block")) != (int) string::npos)
|
||||||
{
|
{
|
||||||
int pos1 = Get_Argument(prhs[i]).find("=", pos+5);
|
size_t pos1 = Get_Argument(prhs[i]).find("=", pos+5);
|
||||||
if (pos1 != (int) string::npos)
|
if (pos1 != (int) string::npos)
|
||||||
pos = pos1 + 1;
|
pos = pos1 + 1;
|
||||||
else
|
else
|
||||||
pos += 5;
|
pos += 5;
|
||||||
block = atoi(Get_Argument(prhs[i]).substr(pos, string::npos).c_str())-1;
|
block = atoi(Get_Argument(prhs[i]).substr(pos, string::npos).c_str())-1;
|
||||||
}
|
}
|
||||||
|
else if ((pos = Get_Argument(prhs[i]).find("pfplan")) != (int) string::npos)
|
||||||
|
{
|
||||||
|
size_t pos1 = Get_Argument(prhs[i]).find("=", pos+6);
|
||||||
|
if (pos1 != (int) string::npos)
|
||||||
|
pos = pos1 + 1;
|
||||||
|
else
|
||||||
|
pos += 6;
|
||||||
|
*pfplan_struct_name = deblank(Get_Argument(prhs[i]).substr(pos, string::npos));
|
||||||
|
}
|
||||||
|
else if ((pos = Get_Argument(prhs[i]).find("plan")) != (int) string::npos)
|
||||||
|
{
|
||||||
|
size_t pos1 = Get_Argument(prhs[i]).find("=", pos+4);
|
||||||
|
if (pos1 != (int) string::npos)
|
||||||
|
pos = pos1 + 1;
|
||||||
|
else
|
||||||
|
pos += 4;
|
||||||
|
*plan_struct_name = deblank(Get_Argument(prhs[i]).substr(pos, string::npos));
|
||||||
|
}
|
||||||
else
|
else
|
||||||
{
|
{
|
||||||
ostringstream tmp;
|
ostringstream tmp;
|
||||||
|
@ -185,6 +386,7 @@ Get_Arguments_and_global_variables(int nrhs,
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
#ifdef DEBUG_EX
|
#ifdef DEBUG_EX
|
||||||
int
|
int
|
||||||
main(int nrhs, const char *prhs[])
|
main(int nrhs, const char *prhs[])
|
||||||
|
@ -203,9 +405,9 @@ main(int nrhs, const char *prhs[])
|
||||||
char *plhs[1];
|
char *plhs[1];
|
||||||
load_global((char *) prhs[1]);
|
load_global((char *) prhs[1]);
|
||||||
#endif
|
#endif
|
||||||
//ErrorHandlingException error_handling;
|
mxArray *plan_struct = NULL, *pfplan_struct = NULL;
|
||||||
unsigned int i, row_y = 0, col_y = 0, row_x = 0, col_x = 0, nb_row_xd = 0;
|
size_t i, row_y = 0, col_y = 0, row_x = 0, col_x = 0, nb_row_xd = 0;
|
||||||
unsigned int steady_row_y, steady_col_y;
|
size_t steady_row_y, steady_col_y;
|
||||||
int y_kmin = 0, y_kmax = 0, y_decal = 0;
|
int y_kmin = 0, y_kmax = 0, y_decal = 0;
|
||||||
unsigned int periods = 1;
|
unsigned int periods = 1;
|
||||||
double *direction;
|
double *direction;
|
||||||
|
@ -218,13 +420,22 @@ main(int nrhs, const char *prhs[])
|
||||||
bool global_temporary_terms = false;
|
bool global_temporary_terms = false;
|
||||||
bool print = false, print_error = true, print_it = false;
|
bool print = false, print_error = true, print_it = false;
|
||||||
double *steady_yd = NULL, *steady_xd = NULL;
|
double *steady_yd = NULL, *steady_xd = NULL;
|
||||||
|
string plan, pfplan;
|
||||||
|
|
||||||
|
vector<s_plan> splan, spfplan;
|
||||||
|
|
||||||
|
#ifdef CUDA
|
||||||
|
int CUDA_device = -1;
|
||||||
|
cublasHandle_t cublas_handle;
|
||||||
|
cusparseHandle_t cusparse_handle;
|
||||||
|
cusparseMatDescr_t descr;
|
||||||
|
#endif
|
||||||
try
|
try
|
||||||
{
|
{
|
||||||
Get_Arguments_and_global_variables(nrhs, prhs, count_array_argument,
|
Get_Arguments_and_global_variables(nrhs, prhs, count_array_argument,
|
||||||
&yd, row_y, col_y,
|
&yd, row_y, col_y,
|
||||||
&xd, row_x, col_x,
|
&xd, row_x, col_x,
|
||||||
¶ms,
|
¶ms,
|
||||||
&steady_yd, steady_row_y, steady_col_y,
|
&steady_yd, steady_row_y, steady_col_y,
|
||||||
periods,
|
periods,
|
||||||
#ifndef DEBUG_EX
|
#ifndef DEBUG_EX
|
||||||
|
@ -232,123 +443,401 @@ main(int nrhs, const char *prhs[])
|
||||||
#endif
|
#endif
|
||||||
steady_state, evaluate, block,
|
steady_state, evaluate, block,
|
||||||
&M_, &oo_, &options_, global_temporary_terms,
|
&M_, &oo_, &options_, global_temporary_terms,
|
||||||
print, print_error, &GlobalTemporaryTerms);
|
print, print_error, &GlobalTemporaryTerms,
|
||||||
|
&plan, &pfplan);
|
||||||
}
|
}
|
||||||
catch (GeneralExceptionHandling &feh)
|
catch (GeneralExceptionHandling &feh)
|
||||||
{
|
{
|
||||||
DYN_MEX_FUNC_ERR_MSG_TXT(feh.GetErrorMsg().c_str());
|
DYN_MEX_FUNC_ERR_MSG_TXT(feh.GetErrorMsg().c_str());
|
||||||
}
|
}
|
||||||
|
|
||||||
if (!count_array_argument)
|
if (!count_array_argument)
|
||||||
params = mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "params")));
|
{
|
||||||
|
int field = mxGetFieldNumber(M_, "params");
|
||||||
|
if (field < 0)
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT("params is not a field of M_");
|
||||||
|
params = mxGetPr(mxGetFieldByNumber(M_, 0, field));
|
||||||
|
}
|
||||||
|
|
||||||
|
ErrorMsg emsg;
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
if (plan.length()>0)
|
||||||
|
{
|
||||||
|
mxArray* plan_struct = mexGetVariable("base", plan.c_str());
|
||||||
|
if (plan_struct == NULL)
|
||||||
|
{
|
||||||
|
string tmp = plan;
|
||||||
|
tmp.insert(0,"Can't find the plan: ");
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
|
||||||
|
}
|
||||||
|
size_t n_plan = mxGetN(plan_struct);
|
||||||
|
splan.resize(n_plan);
|
||||||
|
for (int i = 0; i < n_plan; i++)
|
||||||
|
{
|
||||||
|
splan[i].var = "";
|
||||||
|
splan[i].exo = "";
|
||||||
|
mxArray* tmp = mxGetField(plan_struct, i, "exo");
|
||||||
|
if (tmp)
|
||||||
|
{
|
||||||
|
char name [100];
|
||||||
|
mxGetString(tmp, name, 100);
|
||||||
|
splan[i].var = name;
|
||||||
|
SymbolType variable_type;
|
||||||
|
int exo_num = emsg.get_ID(name, &variable_type);
|
||||||
|
if (variable_type == eExogenous || variable_type == eExogenousDet)
|
||||||
|
splan[i].var_num = exo_num;
|
||||||
|
else
|
||||||
|
{
|
||||||
|
string tmp = name;
|
||||||
|
tmp.insert(0,"the variable '");
|
||||||
|
tmp.append("' defined as var in plan is not an exogenous or a deterministic exogenous\n");
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
|
||||||
|
}
|
||||||
|
}
|
||||||
|
tmp = mxGetField(plan_struct, i, "var");
|
||||||
|
if (tmp)
|
||||||
|
{
|
||||||
|
char name [100];
|
||||||
|
mxGetString(tmp, name, 100);
|
||||||
|
splan[i].exo = name;
|
||||||
|
SymbolType variable_type;
|
||||||
|
int exo_num = emsg.get_ID(name, &variable_type);
|
||||||
|
if (variable_type == eEndogenous)
|
||||||
|
splan[i].exo_num = exo_num;
|
||||||
|
else
|
||||||
|
{
|
||||||
|
string tmp = name;
|
||||||
|
tmp.insert(0,"the variable '");
|
||||||
|
tmp.append("' defined as exo in plan is not an endogenous variable\n");
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
|
||||||
|
}
|
||||||
|
}
|
||||||
|
tmp = mxGetField(plan_struct, i, "per_value");
|
||||||
|
if (tmp)
|
||||||
|
{
|
||||||
|
size_t num_shocks = mxGetM(tmp);
|
||||||
|
(splan[i]).per_value.resize(num_shocks);
|
||||||
|
double * per_value = mxGetPr(tmp);
|
||||||
|
for (int j = 0; j < num_shocks; j++)
|
||||||
|
(splan[i]).per_value[j] = make_pair(ceil(per_value[j]), per_value[j + num_shocks]);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
int i;
|
||||||
|
for (vector<s_plan>::iterator it = splan.begin(); it != splan.end(); it++)
|
||||||
|
{
|
||||||
|
mexPrintf("----------------------------------------------------------------------------------------------------\n");
|
||||||
|
mexPrintf("suprise n°%d\n", i+1);
|
||||||
|
if (it->exo.length())
|
||||||
|
mexPrintf(" plan fliping var=%s (%d) exo=%s (%d) for the following periods and with the following values:\n", it->var.c_str(), it->var_num, it->exo.c_str(), it->exo_num);
|
||||||
|
else
|
||||||
|
mexPrintf(" plan shocks on var=%s for the following periods and with the following values:\n", it->var.c_str());
|
||||||
|
for (vector<pair<int, double> >::iterator it1 = it->per_value.begin(); it1 != it->per_value.end(); it1++)
|
||||||
|
{
|
||||||
|
mexPrintf(" %3d %10.5f\n",it1->first, it1->second);
|
||||||
|
}
|
||||||
|
i++;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
if (pfplan.length()>0)
|
||||||
|
{
|
||||||
|
pfplan_struct = mexGetVariable("base", pfplan.c_str());
|
||||||
|
if (!pfplan_struct)
|
||||||
|
{
|
||||||
|
string tmp = pfplan;
|
||||||
|
tmp.insert(0,"Can't find the pfplan: ");
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
|
||||||
|
}
|
||||||
|
size_t n_plan = mxGetN(pfplan_struct);
|
||||||
|
spfplan.resize(n_plan);
|
||||||
|
for (int i = 0; i < n_plan; i++)
|
||||||
|
{
|
||||||
|
spfplan[i].var = "";
|
||||||
|
spfplan[i].exo = "";
|
||||||
|
mxArray* tmp = mxGetField(pfplan_struct, i, "var");
|
||||||
|
if (tmp)
|
||||||
|
{
|
||||||
|
char name [100];
|
||||||
|
mxGetString(tmp, name, 100);
|
||||||
|
spfplan[i].var = name;
|
||||||
|
SymbolType variable_type;
|
||||||
|
int exo_num = emsg.get_ID(name, &variable_type);
|
||||||
|
if (variable_type == eExogenous || variable_type == eExogenousDet)
|
||||||
|
splan[i].var_num = exo_num;
|
||||||
|
else
|
||||||
|
{
|
||||||
|
string tmp = name;
|
||||||
|
tmp.insert(0,"the variable '");
|
||||||
|
tmp.append("' defined as var in pfplan is not an exogenous or a deterministic exogenous\n");
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
|
||||||
|
}
|
||||||
|
}
|
||||||
|
tmp = mxGetField(pfplan_struct, i, "exo");
|
||||||
|
if (tmp)
|
||||||
|
{
|
||||||
|
char name [100];
|
||||||
|
mxGetString(tmp, name, 100);
|
||||||
|
spfplan[i].exo = name;
|
||||||
|
SymbolType variable_type;
|
||||||
|
int exo_num = emsg.get_ID(name, &variable_type);
|
||||||
|
if (variable_type == eEndogenous)
|
||||||
|
spfplan[i].exo_num = exo_num;
|
||||||
|
else
|
||||||
|
{
|
||||||
|
string tmp = name;
|
||||||
|
tmp.insert(0,"the variable '");
|
||||||
|
tmp.append("' defined as exo in pfplan is not an endogenous variable\n");
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
|
||||||
|
}
|
||||||
|
}
|
||||||
|
tmp = mxGetField(pfplan_struct, i, "per_value");
|
||||||
|
if (tmp)
|
||||||
|
{
|
||||||
|
size_t num_shocks = mxGetM(tmp);
|
||||||
|
double * per_value = mxGetPr(tmp);
|
||||||
|
(spfplan[i]).per_value.resize(num_shocks);
|
||||||
|
for (int j = 0; j < num_shocks; j++)
|
||||||
|
spfplan[i].per_value[j] = make_pair(ceil(per_value[j]), per_value[j+ num_shocks]);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
int i;
|
||||||
|
for (vector<s_plan>::iterator it = spfplan.begin(); it != spfplan.end(); it++)
|
||||||
|
{
|
||||||
|
mexPrintf("----------------------------------------------------------------------------------------------------\n");
|
||||||
|
mexPrintf("perfect foresight n°%d\n", i+1);
|
||||||
|
if (it->exo.length())
|
||||||
|
mexPrintf(" plan flipping var=%s (%d) exo=%s (%d) for the following periods and with the following values:\n", it->var.c_str(), it->var_num, it->exo.c_str(), it->exo_num);
|
||||||
|
else
|
||||||
|
mexPrintf(" plan shocks on var=%s (%d) for the following periods and with the following values:\n", it->var.c_str(), it->var_num);
|
||||||
|
for (vector<pair<int, double> >::iterator it1 = it->per_value.begin(); it1 != it->per_value.end(); it1++)
|
||||||
|
{
|
||||||
|
mexPrintf(" %3d %10.5f\n",it1->first, it1->second);
|
||||||
|
}
|
||||||
|
i++;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
int field_steady_state = mxGetFieldNumber(oo_, "steady_state");
|
||||||
|
if (field_steady_state < 0)
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT("steady_state is not a field of oo_");
|
||||||
|
int field_exo_steady_state = mxGetFieldNumber(oo_, "exo_steady_state");
|
||||||
|
if (field_exo_steady_state < 0)
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT("exo_steady_state is not a field of oo_");
|
||||||
|
|
||||||
if (!steady_state)
|
if (!steady_state)
|
||||||
{
|
{
|
||||||
|
int field_endo_simul = mxGetFieldNumber(oo_, "endo_simul");
|
||||||
|
if (field_endo_simul < 0)
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT("endo_simul is not a field of oo_");
|
||||||
|
|
||||||
|
int field_exo_simul = mxGetFieldNumber(oo_, "exo_simul");
|
||||||
|
if (field_exo_simul < 0)
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT("exo_simul is not a field of oo_");
|
||||||
|
|
||||||
if (!count_array_argument)
|
if (!count_array_argument)
|
||||||
{
|
{
|
||||||
yd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "endo_simul")));
|
mxArray* endo_sim_arr = mxGetFieldByNumber(oo_, 0, field_endo_simul);
|
||||||
row_y = mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "endo_simul")));
|
yd = mxGetPr(endo_sim_arr);
|
||||||
col_y = mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "endo_simul")));
|
row_y = mxGetM(endo_sim_arr);
|
||||||
xd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_simul")));
|
col_y = mxGetN(endo_sim_arr);
|
||||||
row_x = mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_simul")));
|
mxArray* exo_sim_arr = mxGetFieldByNumber(oo_, 0, field_exo_simul);
|
||||||
col_x = mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_simul")));
|
xd = mxGetPr(exo_sim_arr);
|
||||||
|
row_x = mxGetM(exo_sim_arr);
|
||||||
|
col_x = mxGetN(exo_sim_arr);
|
||||||
nb_row_xd = row_x;
|
nb_row_xd = row_x;
|
||||||
}
|
}
|
||||||
|
int field = mxGetFieldNumber(M_, "maximum_lag");
|
||||||
|
if (field >= 0)
|
||||||
|
y_kmin = int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, field)))));
|
||||||
|
else
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT("maximum_lag is not a field of M_");
|
||||||
|
field = mxGetFieldNumber(M_, "maximum_lead");
|
||||||
|
if (field >= 0)
|
||||||
|
y_kmax = int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, field)))));
|
||||||
|
else
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT("maximum_lead is not a field of M_");
|
||||||
|
field = mxGetFieldNumber(M_, "maximum_endo_lag");
|
||||||
|
if (field >= 0)
|
||||||
|
y_decal = max(0, y_kmin-int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, field))))));
|
||||||
|
else
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT("maximum_endo_lag is not a field of M_");
|
||||||
|
|
||||||
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")))))));
|
|
||||||
if (!count_array_argument)
|
if (!count_array_argument)
|
||||||
periods = int (floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "periods"))))));
|
{
|
||||||
|
int field = mxGetFieldNumber(options_, "periods");
|
||||||
|
if (field >= 0)
|
||||||
|
periods = int (floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, field)))));
|
||||||
|
else
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT("options_ is not a field of options_");
|
||||||
|
}
|
||||||
|
|
||||||
if (!steady_yd )
|
if (!steady_yd )
|
||||||
{
|
{
|
||||||
steady_yd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "steady_state")));
|
mxArray* steady_state_arr = mxGetFieldByNumber(oo_, 0, field_steady_state);
|
||||||
steady_row_y = mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "steady_state")));
|
steady_yd = mxGetPr(steady_state_arr);
|
||||||
steady_col_y = mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "steady_state")));;
|
steady_row_y = mxGetM(steady_state_arr);
|
||||||
|
steady_col_y = mxGetN(steady_state_arr);
|
||||||
}
|
}
|
||||||
steady_xd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_steady_state")));
|
steady_xd = mxGetPr(mxGetFieldByNumber(oo_, 0, field_exo_steady_state));
|
||||||
}
|
}
|
||||||
else
|
else
|
||||||
{
|
{
|
||||||
if (!count_array_argument)
|
if (!count_array_argument)
|
||||||
{
|
{
|
||||||
yd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "steady_state")));
|
mxArray* steady_state_arr = mxGetFieldByNumber(oo_, 0, field_steady_state);
|
||||||
row_y = mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "steady_state")));
|
yd = mxGetPr(steady_state_arr);
|
||||||
col_y = mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "steady_state")));;
|
row_y = mxGetM(steady_state_arr);
|
||||||
|
col_y = mxGetN(steady_state_arr);
|
||||||
|
|
||||||
xd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_steady_state")));
|
mxArray* exo_steady_state_arr = mxGetFieldByNumber(oo_, 0, field_exo_steady_state);
|
||||||
row_x = mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_steady_state")));
|
xd = mxGetPr(exo_steady_state_arr);
|
||||||
col_x = mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_steady_state")));
|
row_x = mxGetM(exo_steady_state_arr);
|
||||||
|
col_x = mxGetN(exo_steady_state_arr);
|
||||||
nb_row_xd = row_x;
|
nb_row_xd = row_x;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
int verbose= int(*mxGetPr((mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "verbosity")))));
|
int field = mxGetFieldNumber(options_, "verbosity");
|
||||||
|
int verbose = 0;
|
||||||
|
if (field >= 0)
|
||||||
|
verbose = int(*mxGetPr((mxGetFieldByNumber(options_, 0, field))));
|
||||||
|
else
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT("verbosity is not a field of options_");
|
||||||
if (verbose)
|
if (verbose)
|
||||||
print_it = true;
|
print_it = true;
|
||||||
int maxit_ = int (floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "maxit_"))))));
|
field = mxGetFieldNumber(options_, "maxit_");
|
||||||
double slowc = double (*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "slowc")))));
|
if (field < 0)
|
||||||
double markowitz_c = double (*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "markowitz")))));
|
DYN_MEX_FUNC_ERR_MSG_TXT("maxit_ is not a field of options_");
|
||||||
int minimal_solving_periods = int (*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "minimal_solving_periods")))));
|
int maxit_ = int (floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, field)))));
|
||||||
int stack_solve_algo = int (*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "stack_solve_algo")))));
|
field = mxGetFieldNumber(options_, "slowc");
|
||||||
|
if (field < 0)
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT("slows is not a field of options_");
|
||||||
|
double slowc = double (*(mxGetPr(mxGetFieldByNumber(options_, 0, field))));
|
||||||
|
field = mxGetFieldNumber(options_, "markowitz");
|
||||||
|
if (field < 0)
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT("markowitz is not a field of options_");
|
||||||
|
double markowitz_c = double (*(mxGetPr(mxGetFieldByNumber(options_, 0, field))));
|
||||||
|
field = mxGetFieldNumber(options_, "minimal_solving_periods");
|
||||||
|
if (field < 0)
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT("minimal_solving_periods is not a field of options_");
|
||||||
|
int minimal_solving_periods = int (*(mxGetPr(mxGetFieldByNumber(options_, 0, field))));
|
||||||
|
field = mxGetFieldNumber(options_, "stack_solve_algo");
|
||||||
|
if (field < 0)
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT("stack_solve_algo is not a field of options_");
|
||||||
|
int stack_solve_algo = int (*(mxGetPr(mxGetFieldByNumber(options_, 0, field))));
|
||||||
int solve_algo;
|
int solve_algo;
|
||||||
double solve_tolf;
|
double solve_tolf;
|
||||||
|
|
||||||
if (steady_state)
|
if (steady_state)
|
||||||
{
|
{
|
||||||
solve_algo = int (*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "solve_algo")))));
|
int field = mxGetFieldNumber(options_, "solve_algo");
|
||||||
solve_tolf = *(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "solve_tolf"))));
|
if (field >= 0)
|
||||||
|
solve_algo = int (*(mxGetPr(mxGetFieldByNumber(options_, 0, field))));
|
||||||
|
else
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT("solve_algo is not a field of options_");
|
||||||
|
field = mxGetFieldNumber(options_, "solve_tolf");
|
||||||
|
if (field >= 0)
|
||||||
|
solve_tolf = *(mxGetPr(mxGetFieldByNumber(options_, 0, field)));
|
||||||
|
else
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT("solve_tolf is not a field of options_");
|
||||||
}
|
}
|
||||||
else
|
else
|
||||||
{
|
{
|
||||||
solve_algo = stack_solve_algo;
|
solve_algo = stack_solve_algo;
|
||||||
mxArray *dynatol = mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "dynatol"));
|
int field = mxGetFieldNumber(options_, "dynatol");
|
||||||
solve_tolf= *mxGetPr((mxGetFieldByNumber(dynatol, 0, mxGetFieldNumber(dynatol, "f"))));
|
mxArray *dynatol;
|
||||||
|
if (field >= 0)
|
||||||
|
dynatol = mxGetFieldByNumber(options_, 0, field);
|
||||||
|
else
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT("dynatol is not a field of options_");
|
||||||
|
field = mxGetFieldNumber(dynatol, "f");
|
||||||
|
if (field >= 0)
|
||||||
|
solve_tolf= *mxGetPr((mxGetFieldByNumber(dynatol, 0, field)));
|
||||||
|
else
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT("f is not a field of options_.dynatol");
|
||||||
}
|
}
|
||||||
|
field = mxGetFieldNumber(M_, "fname");
|
||||||
mxArray *mxa = mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "fname"));
|
mxArray *mxa;
|
||||||
int buflen = mxGetM(mxa) * mxGetN(mxa) + 1;
|
if (field >= 0)
|
||||||
|
mxa = mxGetFieldByNumber(M_, 0, field);
|
||||||
|
else
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT("fname is not a field of M_");
|
||||||
|
size_t buflen = mxGetM(mxa) * mxGetN(mxa) + 1;
|
||||||
char *fname;
|
char *fname;
|
||||||
fname = (char *) mxCalloc(buflen+1, sizeof(char));
|
fname = (char *) mxCalloc(buflen+1, sizeof(char));
|
||||||
int status = mxGetString(mxa, fname, buflen);
|
size_t status = mxGetString(mxa, fname, int(buflen));
|
||||||
fname[buflen] = ' ';
|
fname[buflen] = ' ';
|
||||||
if (status != 0)
|
if (status != 0)
|
||||||
mexWarnMsgTxt("Not enough space. Filename is truncated.");
|
mexWarnMsgTxt("Not enough space. Filename is truncated.");
|
||||||
string file_name = fname;
|
string file_name = fname;
|
||||||
|
|
||||||
int size_of_direction = col_y*row_y*sizeof(double);
|
#ifdef CUDA
|
||||||
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);
|
|
||||||
|
|
||||||
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]);
|
|
||||||
ya[i] = double (yd[i]);
|
|
||||||
}
|
|
||||||
int y_size = row_y;
|
|
||||||
int nb_row_x = row_x;
|
|
||||||
clock_t t0 = clock();
|
|
||||||
|
|
||||||
Interpreter interprete(params, y, ya, x, steady_yd, steady_xd, direction, y_size, nb_row_x, nb_row_xd, periods, y_kmin, y_kmax, maxit_, solve_tolf, size_of_direction, slowc, y_decal, markowitz_c, file_name, minimal_solving_periods, stack_solve_algo, solve_algo, global_temporary_terms, print, print_error, GlobalTemporaryTerms);
|
|
||||||
|
|
||||||
string f(fname);
|
|
||||||
mxFree(fname);
|
|
||||||
int nb_blocks = 0;
|
|
||||||
double *pind;
|
|
||||||
bool no_error = true;
|
|
||||||
|
|
||||||
try
|
try
|
||||||
{
|
{
|
||||||
interprete.compute_blocks(f, f, steady_state, evaluate, block, nb_blocks,print_it);
|
if (stack_solve_algo == 7 && !steady_state)
|
||||||
|
CUDA_device = GPU_Test_and_Info(&cublas_handle, &cusparse_handle, &descr);
|
||||||
}
|
}
|
||||||
catch (GeneralExceptionHandling &feh)
|
catch (GeneralExceptionHandling &feh)
|
||||||
{
|
{
|
||||||
DYN_MEX_FUNC_ERR_MSG_TXT(feh.GetErrorMsg().c_str());
|
DYN_MEX_FUNC_ERR_MSG_TXT(feh.GetErrorMsg().c_str());
|
||||||
}
|
}
|
||||||
|
#else
|
||||||
|
if (stack_solve_algo == 7 && !steady_state)
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT("bytecode has not been compiled with CUDA option. Bytecode Can't use options_.stack_solve_algo=7\n");
|
||||||
|
#endif
|
||||||
|
|
||||||
|
size_t 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);
|
||||||
|
double *x = (double *) mxMalloc(col_x*row_x*sizeof(double));
|
||||||
|
#ifdef USE_OMP
|
||||||
|
#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
|
||||||
|
#endif
|
||||||
|
for (i = 0; i < row_x*col_x; i++)
|
||||||
|
{
|
||||||
|
x[i] = double (xd[i]);
|
||||||
|
}
|
||||||
|
|
||||||
|
#ifdef USE_OMP
|
||||||
|
#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
|
||||||
|
#endif
|
||||||
|
for (i = 0; i < row_y*col_y; i++)
|
||||||
|
{
|
||||||
|
y[i] = double (yd[i]);
|
||||||
|
ya[i] = double (yd[i]);
|
||||||
|
}
|
||||||
|
size_t y_size = row_y;
|
||||||
|
size_t nb_row_x = row_x;
|
||||||
|
clock_t t0 = clock();
|
||||||
|
Interpreter interprete(params, y, ya, x, steady_yd, steady_xd, direction, y_size, nb_row_x, nb_row_xd, periods, y_kmin, y_kmax, maxit_, solve_tolf, size_of_direction, slowc, y_decal,
|
||||||
|
markowitz_c, file_name, minimal_solving_periods, stack_solve_algo, solve_algo, global_temporary_terms, print, print_error, GlobalTemporaryTerms, steady_state,
|
||||||
|
print_it
|
||||||
|
#ifdef CUDA
|
||||||
|
, CUDA_device, cublas_handle, cusparse_handle, descr
|
||||||
|
#endif
|
||||||
|
);
|
||||||
|
string f(fname);
|
||||||
|
mxFree(fname);
|
||||||
|
int nb_blocks = 0;
|
||||||
|
double *pind;
|
||||||
|
bool no_error = true;
|
||||||
|
try
|
||||||
|
{
|
||||||
|
interprete.compute_blocks(f, f, evaluate, block, nb_blocks);
|
||||||
|
}
|
||||||
|
catch (GeneralExceptionHandling &feh)
|
||||||
|
{
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT(feh.GetErrorMsg().c_str());
|
||||||
|
}
|
||||||
|
|
||||||
|
#ifdef CUDA
|
||||||
|
if (stack_solve_algo == 7 && !steady_state)
|
||||||
|
GPU_close(cublas_handle, cusparse_handle, descr);
|
||||||
|
#endif
|
||||||
|
|
||||||
clock_t t1 = clock();
|
clock_t t1 = clock();
|
||||||
if (!steady_state && !evaluate && no_error && print)
|
if (!steady_state && !evaluate && no_error && print)
|
||||||
|
@ -370,14 +859,14 @@ main(int nrhs, const char *prhs[])
|
||||||
if (evaluate)
|
if (evaluate)
|
||||||
{
|
{
|
||||||
vector<double> residual = interprete.get_residual();
|
vector<double> residual = interprete.get_residual();
|
||||||
plhs[1] = mxCreateDoubleMatrix(residual.size()/col_y, col_y, mxREAL);
|
plhs[1] = mxCreateDoubleMatrix(int(residual.size()/double(col_y)), int(col_y), mxREAL);
|
||||||
pind = mxGetPr(plhs[1]);
|
pind = mxGetPr(plhs[1]);
|
||||||
for (i = 0; i < residual.size(); i++)
|
for (i = 0; i < residual.size(); i++)
|
||||||
pind[i] = residual[i];
|
pind[i] = residual[i];
|
||||||
}
|
}
|
||||||
else
|
else
|
||||||
{
|
{
|
||||||
plhs[1] = mxCreateDoubleMatrix(row_y, col_y, mxREAL);
|
plhs[1] = mxCreateDoubleMatrix(int(row_y), int(col_y), mxREAL);
|
||||||
pind = mxGetPr(plhs[1]);
|
pind = mxGetPr(plhs[1]);
|
||||||
for (i = 0; i < row_y*col_y; i++)
|
for (i = 0; i < row_y*col_y; i++)
|
||||||
pind[i] = y[i];
|
pind[i] = y[i];
|
||||||
|
@ -385,7 +874,7 @@ main(int nrhs, const char *prhs[])
|
||||||
}
|
}
|
||||||
else
|
else
|
||||||
{
|
{
|
||||||
plhs[1] = mxCreateDoubleMatrix(row_y, col_y, mxREAL);
|
plhs[1] = mxCreateDoubleMatrix(int(row_y), int(col_y), mxREAL);
|
||||||
pind = mxGetPr(plhs[1]);
|
pind = mxGetPr(plhs[1]);
|
||||||
if (evaluate)
|
if (evaluate)
|
||||||
{
|
{
|
||||||
|
@ -409,7 +898,7 @@ main(int nrhs, const char *prhs[])
|
||||||
jacob_exo_field_number = 1;
|
jacob_exo_field_number = 1;
|
||||||
jacob_exo_det_field_number = 2;
|
jacob_exo_det_field_number = 2;
|
||||||
jacob_other_endo_field_number = 3;
|
jacob_other_endo_field_number = 3;
|
||||||
mwSize dims[1] = {nb_blocks };
|
mwSize dims[1] = {(mwSize)nb_blocks };
|
||||||
plhs[2] = mxCreateStructArray(1, dims, 4, field_names);
|
plhs[2] = mxCreateStructArray(1, dims, 4, field_names);
|
||||||
}
|
}
|
||||||
else if (!mxIsStruct(block_structur))
|
else if (!mxIsStruct(block_structur))
|
||||||
|
@ -456,15 +945,15 @@ main(int nrhs, const char *prhs[])
|
||||||
}
|
}
|
||||||
if (nlhs > 3)
|
if (nlhs > 3)
|
||||||
{
|
{
|
||||||
plhs[3] = mxCreateDoubleMatrix(row_y, col_y, mxREAL);
|
plhs[3] = mxCreateDoubleMatrix(int(row_y), int(col_y), mxREAL);
|
||||||
pind = mxGetPr(plhs[3]);
|
pind = mxGetPr(plhs[3]);
|
||||||
for (i = 0; i < row_y*col_y; i++)
|
for (i = 0; i < row_y*col_y; i++)
|
||||||
pind[i] = y[i];
|
pind[i] = y[i];
|
||||||
if (nlhs > 4)
|
if (nlhs > 4)
|
||||||
{
|
{
|
||||||
mxArray *GlobalTemporaryTerms = interprete.get_Temporary_Terms();
|
mxArray *GlobalTemporaryTerms = interprete.get_Temporary_Terms();
|
||||||
unsigned int nb_temp_terms = mxGetM(GlobalTemporaryTerms);
|
size_t nb_temp_terms = mxGetM(GlobalTemporaryTerms);
|
||||||
plhs[4] = mxCreateDoubleMatrix(nb_temp_terms, 1, mxREAL);
|
plhs[4] = mxCreateDoubleMatrix(int(nb_temp_terms), 1, mxREAL);
|
||||||
pind = mxGetPr(plhs[4]);
|
pind = mxGetPr(plhs[4]);
|
||||||
double *tt = mxGetPr(GlobalTemporaryTerms);
|
double *tt = mxGetPr(GlobalTemporaryTerms);
|
||||||
for (i = 0; i < nb_temp_terms; i++)
|
for (i = 0; i < nb_temp_terms; i++)
|
||||||
|
@ -486,4 +975,8 @@ main(int nrhs, const char *prhs[])
|
||||||
mxFree(ya);
|
mxFree(ya);
|
||||||
if (direction)
|
if (direction)
|
||||||
mxFree(direction);
|
mxFree(direction);
|
||||||
|
#ifdef _MSC_VER_
|
||||||
|
/*fFreeResult =*/ FreeLibrary(hinstLib);
|
||||||
|
#endif
|
||||||
|
return;
|
||||||
}
|
}
|
||||||
|
|
|
@ -41,7 +41,7 @@ typedef ptrdiff_t blas_int;
|
||||||
typedef int blas_int;
|
typedef int blas_int;
|
||||||
#endif
|
#endif
|
||||||
|
|
||||||
#if defined(MATLAB_MEX_FILE) && defined(_WIN32)
|
#if defined(MATLAB_MEX_FILE) && defined(_WIN32) && !defined(_MSC_VER)
|
||||||
# define FORTRAN_WRAPPER(x) x
|
# define FORTRAN_WRAPPER(x) x
|
||||||
#else
|
#else
|
||||||
# define FORTRAN_WRAPPER(x) x ## _
|
# define FORTRAN_WRAPPER(x) x ## _
|
||||||
|
|
|
@ -0,0 +1,122 @@
|
||||||
|
/*
|
||||||
|
* Defines the prototypes for BLAS Fortran functions.
|
||||||
|
*
|
||||||
|
* Also defines a typedef blas_int to be used for all integers passed to BLAS
|
||||||
|
* functions.
|
||||||
|
*
|
||||||
|
* When used in the context of a MATLAB MEX file, you must define MATLAB_MEX_FILE
|
||||||
|
* and MATLAB_VERSION (for version 7.4, define it to 0x0704).
|
||||||
|
*
|
||||||
|
*
|
||||||
|
* Copyright (C) 2009-2011 Dynare Team
|
||||||
|
*
|
||||||
|
* This file is part of Dynare.
|
||||||
|
*
|
||||||
|
* Dynare is free software: you can redistribute it and/or modify
|
||||||
|
* it under the terms of the GNU General Public License as published by
|
||||||
|
* the Free Software Foundation, either version 3 of the License, or
|
||||||
|
* (at your option) any later version.
|
||||||
|
*
|
||||||
|
* Dynare is distributed in the hope that it will be useful,
|
||||||
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||||
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||||
|
* GNU General Public License for more details.
|
||||||
|
*
|
||||||
|
* You should have received a copy of the GNU General Public License
|
||||||
|
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
||||||
|
*/
|
||||||
|
|
||||||
|
#ifndef _DYNUMFPACK_H
|
||||||
|
#define _DYNUMFPACK_H
|
||||||
|
|
||||||
|
/* Starting from version 7.8, MATLAB BLAS expects ptrdiff_t arguments for integers */
|
||||||
|
#if defined(MATLAB_MEX_FILE) && MATLAB_VERSION >= 0x0708
|
||||||
|
# ifdef __cplusplus
|
||||||
|
# include <cstddef>
|
||||||
|
# else
|
||||||
|
# include <stddef.h>
|
||||||
|
# endif
|
||||||
|
#endif
|
||||||
|
/*
|
||||||
|
#if defined(MATLAB_MEX_FILE) && defined(_WIN32) && !defined(_MSC_VER)
|
||||||
|
# define FORTRAN_WRAPPER(x) x
|
||||||
|
#else
|
||||||
|
# define FORTRAN_WRAPPER(x) x
|
||||||
|
#endif
|
||||||
|
*/
|
||||||
|
#ifdef __cplusplus
|
||||||
|
extern "C" {
|
||||||
|
#endif
|
||||||
|
|
||||||
|
|
||||||
|
/* -------------------------------------------------------------------------- */
|
||||||
|
/* size of Info and Control arrays */
|
||||||
|
/* -------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
/* These might be larger in future versions, since there are only 3 unused
|
||||||
|
* entries in Info, and no unused entries in Control. */
|
||||||
|
|
||||||
|
#define UMFPACK_INFO 90
|
||||||
|
#define UMFPACK_CONTROL 20
|
||||||
|
/* used in all UMFPACK_report_* routines: */
|
||||||
|
#define UMFPACK_PRL 0 /* print level */
|
||||||
|
/* returned by all routines that use Info: */
|
||||||
|
#define UMFPACK_OK (0)
|
||||||
|
#define UMFPACK_STATUS 0 /* UMFPACK_OK, or other result */
|
||||||
|
|
||||||
|
|
||||||
|
typedef long long int SuiteSparse_long;
|
||||||
|
|
||||||
|
#define umfpack_dl_defaults FORTRAN_WRAPPER(umfpack_dl_defaults)
|
||||||
|
void umfpack_dl_defaults(double Control[UMFPACK_CONTROL]);
|
||||||
|
|
||||||
|
#define umfpack_dl_symbolic FORTRAN_WRAPPER(umfpack_dl_symbolic)
|
||||||
|
SuiteSparse_long umfpack_dl_symbolic(SuiteSparse_long n_row, SuiteSparse_long n_col,
|
||||||
|
const SuiteSparse_long Ap [ ], const SuiteSparse_long Ai [ ],
|
||||||
|
const double Ax [ ], void **Symbolic,
|
||||||
|
const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO]);
|
||||||
|
|
||||||
|
#define umfpack_dl_numeric FORTRAN_WRAPPER(umfpack_dl_numeric)
|
||||||
|
SuiteSparse_long umfpack_dl_numeric(const SuiteSparse_long Ap [ ], const SuiteSparse_long Ai [ ],
|
||||||
|
const double Ax [ ], void *Symbolic, void **Numeric,
|
||||||
|
const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO]);
|
||||||
|
|
||||||
|
#define umfpack_dl_solve FORTRAN_WRAPPER(umfpack_dl_solve)
|
||||||
|
SuiteSparse_long umfpack_dl_solve(SuiteSparse_long sys, const SuiteSparse_long Ap [ ],
|
||||||
|
const SuiteSparse_long Ai [ ], const double Ax [ ],
|
||||||
|
double X [ ], const double B [ ], void *Numeric,
|
||||||
|
const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO]);
|
||||||
|
|
||||||
|
#define umfpack_dl_report_info FORTRAN_WRAPPER(umfpack_dl_report_info)
|
||||||
|
void umfpack_dl_report_info(const double Control [UMFPACK_CONTROL],
|
||||||
|
const double Info [UMFPACK_INFO]);
|
||||||
|
|
||||||
|
#define umfpack_dl_report_status FORTRAN_WRAPPER(umfpack_dl_report_status)
|
||||||
|
void umfpack_dl_report_status(const double Control [UMFPACK_CONTROL],
|
||||||
|
SuiteSparse_long status);
|
||||||
|
|
||||||
|
#define umfpack_dl_free_symbolic FORTRAN_WRAPPER(umfpack_dl_free_symbolic)
|
||||||
|
void umfpack_dl_free_symbolic(void **Symbolic);
|
||||||
|
|
||||||
|
#define umfpack_dl_free_numeric FORTRAN_WRAPPER(umfpack_dl_free_numeric)
|
||||||
|
void umfpack_dl_free_numeric(void **Numeric);
|
||||||
|
|
||||||
|
|
||||||
|
#define umfpack_dl_load_symbolic FORTRAN_WRAPPER(umfpack_dl_load_symbolic )
|
||||||
|
SuiteSparse_long umfpack_dl_load_symbolic (void **Symbolic, char *filename) ;
|
||||||
|
|
||||||
|
#define umfpack_dl_load_numeric FORTRAN_WRAPPER(umfpack_dl_load_numeric )
|
||||||
|
SuiteSparse_long umfpack_dl_load_numeric (void **Numeric, char *filename) ;
|
||||||
|
|
||||||
|
#define umfpack_dl_save_symbolic FORTRAN_WRAPPER(umfpack_dl_save_symbolic )
|
||||||
|
SuiteSparse_long umfpack_dl_save_symbolic (void *Symbolic, char *filename) ;
|
||||||
|
|
||||||
|
#define umfpack_dl_save_numeric FORTRAN_WRAPPER(umfpack_dl_save_numeric )
|
||||||
|
SuiteSparse_long umfpack_dl_save_numeric (void *Numeric, char *filename) ;
|
||||||
|
|
||||||
|
|
||||||
|
#ifdef __cplusplus
|
||||||
|
} /* extern "C" */
|
||||||
|
#endif
|
||||||
|
|
||||||
|
#endif /* DYNUMFPACK */
|
|
@ -1232,7 +1232,7 @@ public:
|
||||||
CompileCode.write(reinterpret_cast<char *>(&row), sizeof(row));
|
CompileCode.write(reinterpret_cast<char *>(&row), sizeof(row));
|
||||||
CompileCode.write(reinterpret_cast<char *>(&col), sizeof(col));
|
CompileCode.write(reinterpret_cast<char *>(&col), sizeof(col));
|
||||||
CompileCode.write(reinterpret_cast<char *>(&function_type), sizeof(function_type));
|
CompileCode.write(reinterpret_cast<char *>(&function_type), sizeof(function_type));
|
||||||
int size = func_name.size();
|
size_t size = func_name.size();
|
||||||
CompileCode.write(reinterpret_cast<char *>(&size), sizeof(int));
|
CompileCode.write(reinterpret_cast<char *>(&size), sizeof(int));
|
||||||
const char *name = func_name.c_str();
|
const char *name = func_name.c_str();
|
||||||
CompileCode.write(reinterpret_cast<const char *>(name), func_name.size());
|
CompileCode.write(reinterpret_cast<const char *>(name), func_name.size());
|
||||||
|
@ -1607,7 +1607,7 @@ class CodeLoad
|
||||||
private:
|
private:
|
||||||
uint8_t *code;
|
uint8_t *code;
|
||||||
unsigned int nb_blocks;
|
unsigned int nb_blocks;
|
||||||
vector<unsigned int> begin_block;
|
vector<size_t> begin_block;
|
||||||
public:
|
public:
|
||||||
|
|
||||||
inline unsigned int
|
inline unsigned int
|
||||||
|
@ -1616,7 +1616,7 @@ public:
|
||||||
return nb_blocks;
|
return nb_blocks;
|
||||||
};
|
};
|
||||||
|
|
||||||
unsigned int inline
|
size_t inline
|
||||||
get_begin_block(int block)
|
get_begin_block(int block)
|
||||||
{
|
{
|
||||||
return begin_block[block];
|
return begin_block[block];
|
||||||
|
|
|
@ -3715,6 +3715,19 @@ DynamicModel::testTrendDerivativesEqualToZero(const eval_context_t &eval_context
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
void
|
||||||
|
DynamicModel::print_trend_vars()
|
||||||
|
{
|
||||||
|
for (trend_symbols_map_t::const_iterator it = nonstationary_symbols_map.begin();
|
||||||
|
it != nonstationary_symbols_map.end(); it++)
|
||||||
|
{
|
||||||
|
cout << "it->first:" << symbol_table.getName(it->first) << " ";
|
||||||
|
it->second->print_deflator();
|
||||||
|
cout << endl;
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
void
|
void
|
||||||
DynamicModel::writeParamsDerivativesFile(const string &basename) const
|
DynamicModel::writeParamsDerivativesFile(const string &basename) const
|
||||||
{
|
{
|
||||||
|
|
|
@ -230,6 +230,7 @@ public:
|
||||||
virtual int getDerivID(int symb_id, int lag) const throw (UnknownDerivIDException);
|
virtual int getDerivID(int symb_id, int lag) const throw (UnknownDerivIDException);
|
||||||
virtual int getDynJacobianCol(int deriv_id) const throw (UnknownDerivIDException);
|
virtual int getDynJacobianCol(int deriv_id) const throw (UnknownDerivIDException);
|
||||||
virtual void addAllParamDerivId(set<int> &deriv_id_set);
|
virtual void addAllParamDerivId(set<int> &deriv_id_set);
|
||||||
|
void print_trend_vars();
|
||||||
|
|
||||||
//! Returns true indicating that this is a dynamic model
|
//! Returns true indicating that this is a dynamic model
|
||||||
virtual bool
|
virtual bool
|
||||||
|
|
|
@ -175,6 +175,64 @@ ExprNode::writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output
|
||||||
// Nothing to do
|
// Nothing to do
|
||||||
}
|
}
|
||||||
|
|
||||||
|
void
|
||||||
|
ExprNode::print_deflator()
|
||||||
|
{
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
void
|
||||||
|
VariableNode::print_deflator()
|
||||||
|
{
|
||||||
|
cout << datatree.symbol_table.getName(symb_id);
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
void
|
||||||
|
UnaryOpNode::print_deflator()
|
||||||
|
{
|
||||||
|
arg->print_deflator();
|
||||||
|
}
|
||||||
|
|
||||||
|
void
|
||||||
|
BinaryOpNode::print_deflator()
|
||||||
|
{
|
||||||
|
arg1->print_deflator();
|
||||||
|
arg2->print_deflator();
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
void
|
||||||
|
TrinaryOpNode::print_deflator()
|
||||||
|
{
|
||||||
|
arg1->print_deflator();
|
||||||
|
arg2->print_deflator();
|
||||||
|
arg3->print_deflator();
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
void
|
||||||
|
ExternalFunctionNode::print_deflator()
|
||||||
|
{
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
void
|
||||||
|
FirstDerivExternalFunctionNode::print_deflator()
|
||||||
|
{
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
void
|
||||||
|
SecondDerivExternalFunctionNode::print_deflator()
|
||||||
|
{
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
void
|
void
|
||||||
ExprNode::compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
|
ExprNode::compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
|
||||||
bool lhs_rhs, const temporary_terms_t &temporary_terms,
|
bool lhs_rhs, const temporary_terms_t &temporary_terms,
|
||||||
|
@ -297,6 +355,12 @@ NumConstNode::collectTemporary_terms(const temporary_terms_t &temporary_terms, t
|
||||||
temporary_terms_inuse.insert(idx);
|
temporary_terms_inuse.insert(idx);
|
||||||
}
|
}
|
||||||
|
|
||||||
|
void
|
||||||
|
NumConstNode::print_deflator()
|
||||||
|
{
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
void
|
void
|
||||||
NumConstNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
|
NumConstNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
|
||||||
const temporary_terms_t &temporary_terms,
|
const temporary_terms_t &temporary_terms,
|
||||||
|
|
|
@ -284,6 +284,7 @@ public:
|
||||||
//! Returns the relative period of the most forward term in this expression
|
//! Returns the relative period of the most forward term in this expression
|
||||||
/*! A negative value means that the expression contains only lagged variables */
|
/*! A negative value means that the expression contains only lagged variables */
|
||||||
virtual int maxLead() const = 0;
|
virtual int maxLead() const = 0;
|
||||||
|
virtual void print_deflator();
|
||||||
|
|
||||||
//! Returns a new expression where all the leads/lags have been shifted backwards by the same amount
|
//! Returns a new expression where all the leads/lags have been shifted backwards by the same amount
|
||||||
/*!
|
/*!
|
||||||
|
@ -455,6 +456,7 @@ public:
|
||||||
virtual expr_t detrend(int symb_id, expr_t trend) const;
|
virtual expr_t detrend(int symb_id, expr_t trend) const;
|
||||||
virtual expr_t cloneDynamic(DataTree &dynamic_datatree) const;
|
virtual expr_t cloneDynamic(DataTree &dynamic_datatree) const;
|
||||||
virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
|
virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
|
||||||
|
virtual void print_deflator();
|
||||||
};
|
};
|
||||||
|
|
||||||
//! Symbol or variable node
|
//! Symbol or variable node
|
||||||
|
@ -483,6 +485,7 @@ public:
|
||||||
virtual double eval(const eval_context_t &eval_context) const throw (EvalException, EvalExternalFunctionException);
|
virtual double eval(const eval_context_t &eval_context) const throw (EvalException, EvalExternalFunctionException);
|
||||||
virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic, deriv_node_temp_terms_t &tef_terms) const;
|
virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic, deriv_node_temp_terms_t &tef_terms) const;
|
||||||
virtual expr_t toStatic(DataTree &static_datatree) const;
|
virtual expr_t toStatic(DataTree &static_datatree) const;
|
||||||
|
virtual void print_deflator();
|
||||||
SymbolType
|
SymbolType
|
||||||
get_type() const
|
get_type() const
|
||||||
{
|
{
|
||||||
|
@ -554,6 +557,7 @@ public:
|
||||||
virtual void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const;
|
virtual void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const;
|
||||||
static double eval_opcode(UnaryOpcode op_code, double v) throw (EvalException, EvalExternalFunctionException);
|
static double eval_opcode(UnaryOpcode op_code, double v) throw (EvalException, EvalExternalFunctionException);
|
||||||
virtual double eval(const eval_context_t &eval_context) const throw (EvalException, EvalExternalFunctionException);
|
virtual double eval(const eval_context_t &eval_context) const throw (EvalException, EvalExternalFunctionException);
|
||||||
|
virtual void print_deflator();
|
||||||
virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic, deriv_node_temp_terms_t &tef_terms) const;
|
virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic, deriv_node_temp_terms_t &tef_terms) const;
|
||||||
//! Returns operand
|
//! Returns operand
|
||||||
expr_t
|
expr_t
|
||||||
|
@ -633,6 +637,7 @@ public:
|
||||||
virtual double eval(const eval_context_t &eval_context) const throw (EvalException, EvalExternalFunctionException);
|
virtual double eval(const eval_context_t &eval_context) const throw (EvalException, EvalExternalFunctionException);
|
||||||
virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic, deriv_node_temp_terms_t &tef_terms) const;
|
virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic, deriv_node_temp_terms_t &tef_terms) const;
|
||||||
virtual expr_t Compute_RHS(expr_t arg1, expr_t arg2, int op, int op_type) const;
|
virtual expr_t Compute_RHS(expr_t arg1, expr_t arg2, int op, int op_type) const;
|
||||||
|
virtual void print_deflator();
|
||||||
//! Returns first operand
|
//! Returns first operand
|
||||||
expr_t
|
expr_t
|
||||||
get_arg1() const
|
get_arg1() const
|
||||||
|
@ -733,6 +738,7 @@ public:
|
||||||
virtual int maxEndoLag() const;
|
virtual int maxEndoLag() const;
|
||||||
virtual int maxExoLag() const;
|
virtual int maxExoLag() const;
|
||||||
virtual int maxLead() const;
|
virtual int maxLead() const;
|
||||||
|
virtual void print_deflator();
|
||||||
virtual expr_t decreaseLeadsLags(int n) const;
|
virtual expr_t decreaseLeadsLags(int n) const;
|
||||||
virtual expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const;
|
virtual expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const;
|
||||||
//! Creates another TrinaryOpNode with the same opcode, but with a possibly different datatree and arguments
|
//! Creates another TrinaryOpNode with the same opcode, but with a possibly different datatree and arguments
|
||||||
|
@ -797,6 +803,7 @@ public:
|
||||||
bool lhs_rhs, const temporary_terms_t &temporary_terms,
|
bool lhs_rhs, const temporary_terms_t &temporary_terms,
|
||||||
const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
|
const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
|
||||||
deriv_node_temp_terms_t &tef_terms) const;
|
deriv_node_temp_terms_t &tef_terms) const;
|
||||||
|
virtual void print_deflator();
|
||||||
|
|
||||||
virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic, deriv_node_temp_terms_t &tef_terms) const;
|
virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic, deriv_node_temp_terms_t &tef_terms) const;
|
||||||
virtual expr_t toStatic(DataTree &static_datatree) const;
|
virtual expr_t toStatic(DataTree &static_datatree) const;
|
||||||
|
@ -855,6 +862,7 @@ public:
|
||||||
bool lhs_rhs, const temporary_terms_t &temporary_terms,
|
bool lhs_rhs, const temporary_terms_t &temporary_terms,
|
||||||
const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
|
const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
|
||||||
deriv_node_temp_terms_t &tef_terms) const;
|
deriv_node_temp_terms_t &tef_terms) const;
|
||||||
|
virtual void print_deflator();
|
||||||
};
|
};
|
||||||
|
|
||||||
class SecondDerivExternalFunctionNode : public ExternalFunctionNode
|
class SecondDerivExternalFunctionNode : public ExternalFunctionNode
|
||||||
|
@ -880,6 +888,7 @@ public:
|
||||||
virtual void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
|
virtual void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
|
||||||
const temporary_terms_t &temporary_terms,
|
const temporary_terms_t &temporary_terms,
|
||||||
deriv_node_temp_terms_t &tef_terms) const;
|
deriv_node_temp_terms_t &tef_terms) const;
|
||||||
|
virtual void print_deflator();
|
||||||
};
|
};
|
||||||
|
|
||||||
#endif
|
#endif
|
||||||
|
|
|
@ -107,6 +107,7 @@ ModFile::addStatementAtFront(Statement *st)
|
||||||
void
|
void
|
||||||
ModFile::checkPass()
|
ModFile::checkPass()
|
||||||
{
|
{
|
||||||
|
dynamic_model.print_trend_vars();
|
||||||
for (vector<Statement *>::iterator it = statements.begin();
|
for (vector<Statement *>::iterator it = statements.begin();
|
||||||
it != statements.end(); it++)
|
it != statements.end(); it++)
|
||||||
(*it)->checkPass(mod_file_struct, warnings);
|
(*it)->checkPass(mod_file_struct, warnings);
|
||||||
|
@ -374,8 +375,8 @@ ModFile::computingPass(bool no_tmp_terms)
|
||||||
// Mod file may have no equation (for example in a standalone BVAR estimation)
|
// Mod file may have no equation (for example in a standalone BVAR estimation)
|
||||||
if (dynamic_model.equation_number() > 0)
|
if (dynamic_model.equation_number() > 0)
|
||||||
{
|
{
|
||||||
if (nonstationary_variables)
|
/*if (nonstationary_variables)
|
||||||
trend_dynamic_model.runTrendTest(global_eval_context);
|
trend_dynamic_model.runTrendTest(global_eval_context);*/
|
||||||
|
|
||||||
// Compute static model and its derivatives
|
// Compute static model and its derivatives
|
||||||
dynamic_model.toStatic(static_model);
|
dynamic_model.toStatic(static_model);
|
||||||
|
|
|
@ -1002,14 +1002,21 @@ ModelTree::computeJacobian(const set<int> &vars)
|
||||||
{
|
{
|
||||||
for (set<int>::const_iterator it = vars.begin();
|
for (set<int>::const_iterator it = vars.begin();
|
||||||
it != vars.end(); it++)
|
it != vars.end(); it++)
|
||||||
for (int eq = 0; eq < (int) equations.size(); eq++)
|
{
|
||||||
{
|
int prev_deriv = NNZDerivatives[0];
|
||||||
expr_t d1 = equations[eq]->getDerivative(*it);
|
for (int eq = 0; eq < (int) equations.size(); eq++)
|
||||||
if (d1 == Zero)
|
{
|
||||||
continue;
|
expr_t d1 = equations[eq]->getDerivative(*it);
|
||||||
first_derivatives[make_pair(eq, *it)] = d1;
|
if (d1 == Zero)
|
||||||
++NNZDerivatives[0];
|
continue;
|
||||||
}
|
first_derivatives[make_pair(eq, *it)] = d1;
|
||||||
|
++NNZDerivatives[0];
|
||||||
|
}
|
||||||
|
if (NNZDerivatives[0] == prev_deriv)
|
||||||
|
{
|
||||||
|
cout << "the derivatives w.r. to " << symbol_table.getName(*it) << " is always equal to 0\n";
|
||||||
|
}
|
||||||
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
void
|
void
|
||||||
|
|
|
@ -13,6 +13,7 @@ rho = 0.7;
|
||||||
psi = 0.787;
|
psi = 0.787;
|
||||||
del = 0.02;
|
del = 0.02;
|
||||||
|
|
||||||
|
//model(block, bytecode);
|
||||||
model;
|
model;
|
||||||
dA = exp(gam+e_a);
|
dA = exp(gam+e_a);
|
||||||
log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
|
log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
|
||||||
|
@ -56,19 +57,39 @@ steady;
|
||||||
|
|
||||||
check;
|
check;
|
||||||
|
|
||||||
stoch_simul(irf=0);
|
//stoch_simul(irf=0);
|
||||||
|
|
||||||
conditional_forecast_paths;
|
conditional_forecast_paths;
|
||||||
|
var gp_obs;
|
||||||
|
periods 1 2:5;
|
||||||
|
//values 0.05;
|
||||||
|
//values 0.98 1.00797;
|
||||||
|
values 0.98 0.99;
|
||||||
|
//expectation perfect_foresight;
|
||||||
var gy_obs;
|
var gy_obs;
|
||||||
periods 1 2 3:5;
|
periods 1 2 3:5;
|
||||||
values 0.01 -0.02 0;
|
//values 0.01 -0.02 0;
|
||||||
var gp_obs;
|
//values 0.85 0.85 0.95;
|
||||||
periods 1:5;
|
values 0.95 0.95 0.99;
|
||||||
values 0.05;
|
//expectation perfect_foresight;
|
||||||
end;
|
end;
|
||||||
|
|
||||||
conditional_forecast(parameter_set=calibration, controlled_varexo=(e_a,e_m));
|
options_.stack_solve_algo = 0;
|
||||||
|
options_.maxit_ = 50;
|
||||||
|
|
||||||
if ~(exist('OCTAVE_VERSION') && octave_ver_less_than('3.4.0'))
|
conditional_forecast(parameter_set=calibration, controlled_varexo=(e_m,e_a), simulation_type = deterministic);
|
||||||
plot_conditional_forecast(periods=10) gy_obs gp_obs;
|
|
||||||
end
|
/*shocks;
|
||||||
|
var e_a;
|
||||||
|
periods 1 2 3 4 5;
|
||||||
|
values -0.0109 -0.0122 -0.0137 -0.0154 -0.0173;
|
||||||
|
var e_m;
|
||||||
|
periods 1 2 3 4 5;
|
||||||
|
values -0.1242 -0.0386 -0.0392 -0.0398 -0.0405;
|
||||||
|
end;
|
||||||
|
simul(periods=40);*/
|
||||||
|
rplot gy_obs;
|
||||||
|
rplot gp_obs;
|
||||||
|
//if ~(exist('OCTAVE_VERSION') && octave_ver_less_than('3.4.0'))
|
||||||
|
//plot_conditional_forecast(periods=10) gy_obs gp_obs;
|
||||||
|
//end
|
||||||
|
|
Loading…
Reference in New Issue