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

- Correction of several bugs with Octave 

- Add wrapper needed to compute the steady-state using bytecode and solve_algo = 0, ..., 4
time-shift
Ferhat Mihoubi 2010-10-18 17:28:21 +02:00
parent cab8941c29
commit 1a09426706
18 changed files with 654 additions and 163 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

27
matlab/block_bytecode_mfs_steadystate.m Normal file
View File

@ -0,0 +1,27 @@
function [r, g1] = block_bytecode_mfs_steadystate(y, b)
% Wrapper around the *_static.m file, for use with dynare_solve,
% when block_mfs option is given to steady.
% Copyright (C) 2009 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/>.
global M_ oo_
ss = oo_.steady_state;
indx = M_.blocksMFS{b};
ss(indx) = y;
x = [oo_.exo_steady_state; oo_.exo_det_steady_state];
eval(['[chk, r, g1] = bytecode( ss, x, M_.params, 1, ss, ''evaluate'', ''static'', ''block = ' int2str(b) ''');']);

24
matlab/bytecode_steadystate.m Normal file
View File

@ -0,0 +1,24 @@
function [r, g1] = bytecode_steadystate(y)
% Wrapper around the *_static.m file, for use with dynare_solve,
% when block_mfs option is given to steady.
% Copyright (C) 2009 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/>.
global M_ oo_
x = [oo_.exo_steady_state; oo_.exo_det_steady_state];
eval('[chk, r, g1] = bytecode( y, x, M_.params, 1, x, ''evaluate'', ''static'', ''block = 1'');');

14
matlab/solve_one_boundary.m
View File

@ -1,4 +1,4 @@
function [y, info] = solve_one_boundary(fname, y, x, params, y_index_eq, nze, periods, is_linear, Block_Num, y_kmin, maxit_, solve_tolf, lambda, cutoff, stack_solve_algo, forward_backward, is_dynamic, verbose)
function [y, info] = solve_one_boundary(fname, y, x, params, y_index_eq, nze, periods, is_linear, Block_Num, y_kmin, maxit_, solve_tolf, lambda, cutoff, stack_solve_algo, forward_backward, is_dynamic, verbose, indirect_call)
% Computes the deterministic simulation of a block of equation containing
% lead or lag variables
%
@ -38,9 +38,13 @@ function [y, info] = solve_one_boundary(fname, y, x, params, y_index_eq, nze, pe
% field remains unchanged
% verbose [integer] (0) iterations are not printed
% (1) iterations are printed
%
% OUTPUTS
% indirect_call [integer] (0) direct call to the fname
% (1) indirect call via the
% local_fname wrapper
% OUTPUTS
% y [matrix] All endogenous variables of the model
% info [integer] >=0 no error
% <0 error
%
% ALGORITHM
% Newton with LU or GMRES or BicGstab for dynamic block
@ -84,7 +88,7 @@ else
start = periods+y_kmin;
finish = y_kmin+1;
end
lambda=1;
%lambda=1;
for it_=start:incr:finish
cvg=0;
iter=0;
@ -260,7 +264,7 @@ for it_=start:incr:finish
[yn,info] = csolve(@local_fname, y(y_index_eq),@local_fname,1e-6,500, x, params, y, y_index_eq, fname, 1);
dx = ya - yn;
y(y_index_eq) = yn;
elseif((stack_solve_algo==1 & is_dynamic) | (~is_dynamic & options_.solve_algo==1)),
elseif((stack_solve_algo==1 & is_dynamic) | (stack_solve_algo==0 & is_dynamic) | (~is_dynamic & options_.solve_algo==1)),
dx = g1\r;
ya = ya - lambda*dx;
if(is_dynamic)

3
matlab/solve_two_boundaries.m
View File

@ -247,7 +247,6 @@ while ~(cvg==1 | iter>maxit_),
Elem = first_elem:last_elem;
za = b(Elem) - g1a(Elem, Elem_1) * zaa;
zaa = za;
%y_Elem_1 = Blck_size * (t)+1:Blck_size * (t+1);
y_Elem = Blck_size * (t-1)+1:Blck_size * (t);
dx(y_Elem) = za - ya(y_Elem);
ya(y_Elem) = ya(y_Elem) + lambda*dx(y_Elem);
@ -309,7 +308,7 @@ while ~(cvg==1 | iter>maxit_),
end
end
iter=iter+1;
disp(['iteration: ' num2str(iter,'%d') ' error: ' num2str(max_res,'%e') ' stack_solve_algo=' num2str(stack_solve_algo)]);
disp(['iteration: ' num2str(iter,'%d') ' error: ' num2str(max_res,'%e')]);
end;
if (iter>maxit_)
disp(['No convergence after ' num2str(iter,'%4d') ' iterations in Block ' num2str(Block_Num,'%d')]);

14
matlab/steady_.m
View File

@ -33,8 +33,8 @@ if options_.bytecode && ...
(options_.solve_algo < 0 || options_.solve_algo > 8)
error('STEADY: for the moment, you must use solve_algo=1, 2, 3, 4, 5, 6, 7, 8 with bytecode option')
end
if ~options_.bytecode && options_.solve_algo == 5
error('STEADY: you can''t yet use solve_algo=5 without bytecode option')
if ~options_.bytecode && options_.solve_algo == 8
error('STEADY: you can''t yet use solve_algo=8 without bytecode option')
end
if options_.steadystate_flag
@ -126,7 +126,7 @@ elseif options_.bytecode
if options_.solve_algo > 4
[check, oo_.steady_state] = bytecode('static');
mexErrCheck('bytecode', check);
else
elseif options_.block
for b = 1:size(M_.blocksMFS,1)
n = size(M_.blocksMFS{b}, 1);
ss = oo_.steady_state;
@ -144,6 +144,14 @@ elseif options_.bytecode
oo_.exo_det_steady_state], M_.params, 1, options_.solve_algo, 'static', 'evaluate', ['block = ' int2str(b)]);
end;
end
else
[y, check] = dynare_solve('bytecode_steadystate', ...
oo_.steady_state, ...
options_.jacobian_flag);
if check ~= 0
error('STEADY: convergence problems')
end
oo_.steady_state = y;
end
else
[oo_.steady_state,check] = dynare_solve([M_.fname '_static'],...

28
mex/sources/bytecode/Interpreter.cc
View File

@ -111,8 +111,8 @@ Interpreter::get_variable(const SymbolType variable_type, const unsigned int var
if (variable_num < nb_endo)
{
for (unsigned int i = 0; i < endo_name_length; i++)
if (P_endo_names[2*(variable_num+i*nb_endo)] != ' ')
res << P_endo_names[2*(variable_num+i*nb_endo)];
if (P_endo_names[CHAR_LENGTH*(variable_num+i*nb_endo)] != ' ')
res << P_endo_names[CHAR_LENGTH*(variable_num+i*nb_endo)];
}
else
mexPrintf("=> Unknown endogenous variable n° %d",variable_num);
@ -122,8 +122,8 @@ Interpreter::get_variable(const SymbolType variable_type, const unsigned int var
if (variable_num < nb_exo)
{
for (unsigned int i = 0; i < exo_name_length; i++)
if (P_exo_names[2*(variable_num+i*nb_exo)] != ' ')
res << P_exo_names[2*(variable_num+i*nb_exo)];
if (P_exo_names[CHAR_LENGTH*(variable_num+i*nb_exo)] != ' ')
res << P_exo_names[CHAR_LENGTH*(variable_num+i*nb_exo)];
}
else
mexPrintf("=> Unknown exogenous variable n° %d",variable_num);
@ -132,8 +132,8 @@ Interpreter::get_variable(const SymbolType variable_type, const unsigned int var
if (variable_num < nb_param)
{
for (unsigned int i = 0; i < param_name_length; i++)
if (P_param_names[2*(variable_num+i*nb_param)] != ' ')
res << P_param_names[2*(variable_num+i*nb_param)];
if (P_param_names[CHAR_LENGTH*(variable_num+i*nb_param)] != ' ')
res << P_param_names[CHAR_LENGTH*(variable_num+i*nb_param)];
}
else
mexPrintf("=> Unknown parameter n° %d",variable_num);
@ -2597,7 +2597,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
if (cvg)
continue;
int prev_iter = iter;
Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, solve_algo);
Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, stack_solve_algo, solve_algo);
iter++;
if (iter > prev_iter)
{
@ -2625,7 +2625,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
error_not_printed = true;
compute_block_time(0, false, block_num, size, steady_state);
cvg = false;
Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, solve_algo);
Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, stack_solve_algo, solve_algo);
if (!result)
{
mexPrintf(" in Solve Forward complete, convergence not achieved in block %d\n", Block_Count+1);
@ -2677,7 +2677,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
if (cvg)
continue;
int prev_iter = iter;
Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, solve_algo);
Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, stack_solve_algo, solve_algo);
iter++;
if (iter > prev_iter)
{
@ -2707,7 +2707,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
error_not_printed = true;
compute_block_time(0, false, block_num, size, steady_state);
cvg = false;
Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, solve_algo);
Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, stack_solve_algo, solve_algo);
}
}
}
@ -2765,7 +2765,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
if (cvg)
continue;
int prev_iter = iter;
Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, solve_algo);
Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, stack_solve_algo, solve_algo);
iter++;
if (iter > prev_iter)
{
@ -2792,7 +2792,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
error_not_printed = true;
compute_block_time(0, false, block_num, size, steady_state);
cvg = false;
Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, solve_algo);
Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, stack_solve_algo, solve_algo);
if (!result)
{
mexPrintf(" in Solve Backward complete, convergence not achieved in block %d\n", Block_Count+1);
@ -2844,7 +2844,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
if (cvg)
continue;
int prev_iter = iter;
Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, solve_algo);
Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, stack_solve_algo, solve_algo);
iter++;
if (iter > prev_iter)
{
@ -2870,7 +2870,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
error_not_printed = true;
compute_block_time(0, false, block_num, size, steady_state);
cvg = false;
Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, solve_algo);
Simulate_Newton_One_Boundary(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, stack_solve_algo, solve_algo);
}
}
}

511
mex/sources/bytecode/SparseMatrix.cc
View File

@ -42,7 +42,6 @@ SparseMatrix::SparseMatrix()
restart = 0;
IM_i.clear();
lu_inc_tol = 1e-10;
reduced = true;
}
int
@ -366,7 +365,7 @@ SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, i
}
void
SparseMatrix::Simple_Init(int it_, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM)
SparseMatrix::Simple_Init(int it_, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM, bool &zero_solution)
{
int i, eq, var, lag;
map<pair<pair<int, int>, int>, int>::iterator it4;
@ -435,15 +434,24 @@ SparseMatrix::Simple_Init(int it_, int y_kmin, int y_kmax, int Size, map<pair<pa
it4++;
}
//#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
double cum_abs_sum = 0;
for (i = 0; i < Size; i++)
b[i] = i;
{
b[i] = i;
cum_abs_sum += fabs(u[i]);
}
if (cum_abs_sum < 1e-20)
zero_solution = true;
else
zero_solution = false;
mxFree(temp_NZE_R);
mxFree(temp_NZE_C);
u_count = u_count1;
}
void
SparseMatrix::Init_Matlab_Sparse_Simple(int Size, map<pair<pair<int, int>, int>, int> &IM, mxArray *A_m, mxArray *b_m)
SparseMatrix::Init_Matlab_Sparse_Simple(int Size, map<pair<pair<int, int>, int>, int> &IM, mxArray *A_m, mxArray *b_m, bool &zero_solution)
{
int i, eq, var;
double *b = mxGetPr(b_m);
@ -482,10 +490,19 @@ SparseMatrix::Init_Matlab_Sparse_Simple(int Size, map<pair<pair<int, int>, int>,
#endif
unsigned int NZE = 0;
int last_var = 0;
double cum_abs_sum = 0;
for (i = 0; i < Size; i++)
b[i] = u[i];
{
b[i] = u[i];
cum_abs_sum += fabs(b[i]);
}
if (cum_abs_sum < 1e-20)
zero_solution = true;
else
zero_solution = false;
Aj[0] = 0;
last_var = 0;
last_var = -1;
it4 = IM.begin();
while (it4 != IM.end())
{
@ -1321,10 +1338,307 @@ SparseMatrix::Check_the_Solution(int periods, int y_kmin, int y_kmax, int Size,
}
}
mxArray*
SparseMatrix::substract_A_B(mxArray* A_m, mxArray* B_m)
{
unsigned int n_A = mxGetN(A_m);
unsigned int m_A = mxGetM(A_m);
double* A_d = mxGetPr(A_m);
unsigned int n_B = mxGetN(B_m);
double* B_d = mxGetPr(B_m);
mxArray* C_m = mxCreateDoubleMatrix(m_A, n_B, mxREAL);
double* C_d = mxGetPr(C_m);
for (unsigned int j = 0; j < n_A; j++)
for (unsigned int i = 0; i < m_A; i++)
{
unsigned int index = j*m_A+i;
C_d[index] = A_d[index] - B_d[index];
}
return C_m;
}
mxArray*
SparseMatrix::Sparse_substract_A_SB(mxArray* A_m, mxArray* B_m)
{
unsigned int n_B = mxGetN(B_m);
unsigned int m_B = mxGetM(B_m);
mwIndex* B_i = mxGetIr(B_m);
mwIndex* B_j = mxGetJc(B_m);
unsigned int total_nze_B = B_j[n_B];
double* B_d = mxGetPr(B_m);
mxArray* C_m = mxDuplicateArray(A_m);
double* C_d = mxGetPr(C_m);
unsigned int nze_B = 0;
unsigned int B_col = 0;
while (nze_B < total_nze_B)
{
while (nze_B >= (unsigned int)B_j[B_col+1] && (nze_B < total_nze_B))
B_col++;
C_d[B_col*m_B+B_i[nze_B]] -= B_d[nze_B];
nze_B++;
}
return C_m;
}
mxArray*
SparseMatrix::Sparse_substract_SA_SB(mxArray* A_m, mxArray* B_m)
{
unsigned int n_A = mxGetN(A_m);
unsigned int m_A = mxGetM(A_m);
mwIndex* A_i = mxGetIr(A_m);
mwIndex* A_j = mxGetJc(A_m);
unsigned int total_nze_A = A_j[n_A];
double* A_d = mxGetPr(A_m);
unsigned int n_B = mxGetN(B_m);
mwIndex* B_i = mxGetIr(B_m);
mwIndex* B_j = mxGetJc(B_m);
unsigned int total_nze_B = B_j[n_B];
double* B_d = mxGetPr(B_m);
mxArray* C_m = mxCreateSparse(m_A, n_B, m_A*n_B, mxREAL);
mwIndex* C_i = mxGetIr(C_m);
mwIndex* C_j = mxGetJc(C_m);
double* C_d = mxGetPr(C_m);
unsigned int nze_B = 0, nze_C = 0, nze_A = 0;
unsigned int A_col = 0, B_col = 0, C_col = 0;
C_j[C_col] = 0;
while (nze_A < total_nze_A || nze_B < total_nze_B)
{
while (nze_A >= (unsigned int)A_j[A_col+1] && (nze_A < total_nze_A))
A_col++;
int A_row = A_i[nze_A];
while (nze_B >= (unsigned int)B_j[B_col+1] && (nze_B < total_nze_B))
B_col++;
int B_row = B_i[nze_B];
if (A_col == B_col)
{
if (A_row == B_row && (nze_B < total_nze_B && nze_A < total_nze_A))
{
C_d[nze_C] = A_d[nze_A++] - B_d[nze_B++];
C_i[nze_C] = A_row;
while (C_col < A_col)
C_j[++C_col] = nze_C;
C_j[A_col+1] = nze_C++;
C_col = A_col;
}
else if (A_row < B_row || (nze_B >= total_nze_B && nze_A < total_nze_A))
{
C_d[nze_C] = A_d[nze_A++];
C_i[nze_C] = A_row;
while (C_col < A_col)
C_j[++C_col] = nze_C;
C_j[A_col+1] = nze_C++;
C_col = A_col;
}
else
{
C_d[nze_C] = - B_d[nze_B++];
C_i[nze_C] = B_row;
while (C_col < B_col)
C_j[++C_col] = nze_C;
C_j[B_col+1] = nze_C++;
C_col = B_col;
}
}
else if (A_col < B_col || (nze_B >= total_nze_B && nze_A < total_nze_A))
{
C_d[nze_C] = A_d[nze_A++];
C_i[nze_C] = A_row;
while (C_col < A_col)
C_j[++C_col] = nze_C;
C_j[A_col+1] = nze_C++;
C_col = A_col;
}
else
{
C_d[nze_C] = - B_d[nze_B++];
C_i[nze_C] = B_row;
while (C_col < B_col)
C_j[++C_col] = nze_C;
C_j[B_col+1] = nze_C++;
C_col = B_col;
}
}
while (C_col < n_B)
C_j[++C_col] = nze_C;
mxSetNzmax(C_m, nze_C);
return C_m;
}
mxArray*
SparseMatrix::mult_SAT_B(mxArray* A_m, mxArray* B_m)
{
unsigned int n_A = mxGetN(A_m);
unsigned int m_A = mxGetM(A_m);
mwIndex* A_i = mxGetIr(A_m);
mwIndex* A_j = mxGetJc(A_m);
double* A_d = mxGetPr(A_m);
unsigned int n_B = mxGetN(B_m);
double* B_d = mxGetPr(B_m);
mxArray* C_m = mxCreateDoubleMatrix(m_A, n_B, mxREAL);
double* C_d = mxGetPr(C_m);
unsigned int nze_A = 0;
for (unsigned int j = 0; j < n_B; j++)
{
for (unsigned int i = 0; i < n_A; i++)
{
double sum = 0;
nze_A = A_j[i];
while ( nze_A < (unsigned int)A_j[i+1] )
{
unsigned int i_A = A_i[nze_A];
sum += A_d[nze_A++] * B_d[i_A];
}
C_d[j*n_A+i] = sum;
}
}
return C_m;
}
mxArray*
SparseMatrix::Sparse_mult_SAT_B(mxArray* A_m, mxArray* B_m)
{
unsigned int n_A = mxGetN(A_m);
unsigned int m_A = mxGetM(A_m);
mwIndex* A_i = mxGetIr(A_m);
mwIndex* A_j = mxGetJc(A_m);
double* A_d = mxGetPr(A_m);
unsigned int n_B = mxGetN(B_m);
unsigned int m_B = mxGetM(B_m);
double* B_d = mxGetPr(B_m);
mxArray* C_m = mxCreateSparse(m_A, n_B, m_A*n_B, mxREAL);
mwIndex* C_i = mxGetIr(C_m);
mwIndex* C_j = mxGetJc(C_m);
double* C_d = mxGetPr(C_m);
unsigned int nze_C = 0, nze_A = 0;
unsigned int C_col = 0;
C_j[C_col] = 0;
for (unsigned int j = 0; j < n_B; j++)
{
for (unsigned int i = 0; i < n_A; i++)
{
double sum = 0;
nze_A = A_j[i];
while ( nze_A < (unsigned int)A_j[i+1] )
{
unsigned int i_A = A_i[nze_A];
sum += A_d[nze_A++] * B_d[i_A];
}
if (fabs(sum) > 1e-10)
{
C_d[nze_C] = sum;
C_i[nze_C] = i;
while (C_col < j)
C_j[++C_col] = nze_C;
nze_C++;
}
}
}
while (C_col < m_B)
C_j[++C_col] = nze_C;
mxSetNzmax(C_m, nze_C);
return C_m;
}
mxArray*
SparseMatrix::Sparse_mult_SAT_SB(mxArray* A_m, mxArray* B_m)
{
unsigned int n_A = mxGetN(A_m);
unsigned int m_A = mxGetM(A_m);
mwIndex* A_i = mxGetIr(A_m);
mwIndex* A_j = mxGetJc(A_m);
double* A_d = mxGetPr(A_m);
unsigned int n_B = mxGetN(B_m);
mwIndex* B_i = mxGetIr(B_m);
mwIndex* B_j = mxGetJc(B_m);
double* B_d = mxGetPr(B_m);
mxArray* C_m = mxCreateSparse(m_A, n_B, m_A*n_B, mxREAL);
mwIndex* C_i = mxGetIr(C_m);
mwIndex* C_j = mxGetJc(C_m);
double* C_d = mxGetPr(C_m);
unsigned int nze_B = 0, nze_C = 0, nze_A = 0;
unsigned int C_col = 0;
C_j[C_col] = 0;
for (unsigned int j = 0; j < n_B; j++)
{
for (unsigned int i = 0; i < n_A; i++)
{
double sum = 0;
nze_B = B_j[j];
nze_A = A_j[i];
while ( nze_A < (unsigned int)A_j[i+1] && nze_B < (unsigned int)B_j[j+1])
{
unsigned int i_A = A_i[nze_A];
unsigned int i_B = B_i[nze_B];
if ( i_A == i_B)
sum += A_d[nze_A++] * B_d[nze_B++];
else if (i_A < i_B)
nze_A++;
else
nze_B++;
}
if (fabs(sum) > 1e-10)
{
C_d[nze_C] = sum;
C_i[nze_C] = i;
while (C_col < j)
C_j[++C_col] = nze_C;
nze_C++;
}
}
}
while (C_col < n_B)
C_j[++C_col] = nze_C;
mxSetNzmax(C_m, nze_C);
return C_m;
}
mxArray*
SparseMatrix::Sparse_transpose(mxArray* A_m)
{
unsigned int n_A = mxGetN(A_m);
unsigned int m_A = mxGetM(A_m);
mwIndex* A_i = mxGetIr(A_m);
mwIndex* A_j = mxGetJc(A_m);
unsigned int total_nze_A = A_j[n_A];
double* A_d = mxGetPr(A_m);
mxArray* C_m = mxCreateSparse(n_A, m_A, total_nze_A, mxREAL);
mwIndex* C_i = mxGetIr(C_m);
mwIndex* C_j = mxGetJc(C_m);
double* C_d = mxGetPr(C_m);
unsigned int nze_C = 0, nze_A = 0;
memset(C_j, 0, m_A);
map<pair<unsigned int, unsigned int>, double> B2;
for (unsigned int i = 0; i < n_A; i++)
{
while ( nze_A < (unsigned int)A_j[i+1])
{
C_j[A_i[nze_A]+1]++;
B2[make_pair(A_i[nze_A], i)] = A_d[nze_A];
nze_A++;
}
}
for (unsigned int i = 0; i < m_A; i++)
C_j[i+1] += C_j[i];
for (map<pair<unsigned int, unsigned int>, double>::const_iterator it = B2.begin(); it != B2.end(); it++)
{
C_d[nze_C] = it->second;
C_i[nze_C++] = it->first.second;
}
return C_m;
}
void
SparseMatrix::Solve_Matlab_Relaxation(mxArray* A_m, mxArray* b_m, unsigned int Size, double slowc_l, bool is_two_boundaries, int it_)
{
mexPrintf("Solve_Matlab_Relaxation Size=%d\n", Size);
mxArray *B1, *C1, *A2, *B2, *A3, *b1, *b2;
double* b_m_d = mxGetPr(b_m);
if (!b_m_d)
@ -1407,7 +1721,7 @@ SparseMatrix::Solve_Matlab_Relaxation(mxArray* A_m, mxArray* b_m, unsigned int S
A3*/
while (var < 2*Size && nze < max_nze)
{
if (A_m_j[var+1] <= nze)
if ((unsigned int)A_m_j[var+1] <= nze)
{
if (var < Size)
b1_d[var] = b_m_d[var];
@ -1474,20 +1788,19 @@ SparseMatrix::Solve_Matlab_Relaxation(mxArray* A_m, mxArray* b_m, unsigned int S
B2_j[++B2_var] = B2_nze;
while(A3_var < Size)
A3_j[++A3_var] = A3_nze;
mxArray *d1;
mxArray *val[2];
vector<pair<mxArray*, mxArray*> > triangular_form;
//mxArray* C_B1_inv;
int last_t = 0;
double sumc=0, C_sumc = 1000;
mxArray *B1_inv;
mxArray *B1_inv = NULL;
mxArray *B1_inv_t = NULL;
for (int t = 1; t <= periods; t++)
{
if (abs(sumc / C_sumc -1) > 1e-10*res1)
{
C_sumc = sumc;
mxDestroyArray(B1_inv);
if (B1_inv)
mxDestroyArray(B1_inv);
mexCallMATLAB(1, &B1_inv, 1, &B1, "inv");
mwIndex *B_inv_j = mxGetJc(B1_inv);
unsigned int B_inv_nze = B_inv_j[Size];
@ -1497,42 +1810,31 @@ SparseMatrix::Solve_Matlab_Relaxation(mxArray* A_m, mxArray* b_m, unsigned int S
sumc += fabs(B_inv_d[i]);
last_t = t;
}
B1_inv_t = Sparse_transpose(B1_inv);
mxArray *S1 = Sparse_mult_SAT_SB(B1_inv_t, C1);
mxArray *S1;
val[0] = B1_inv;
val[1] = C1;
mexCallMATLAB(1, &S1, 2, val, "*");
val[1] = b1;
mexCallMATLAB(1, &d1, 2, val, "*");
d1 = mult_SAT_B(B1_inv_t, b1);
if (t < periods)
//Computation for the next lines
{
val[0] = A2;
val[1] = S1;
mxDestroyArray(B1_inv_t);
mxArray *A2_t = Sparse_transpose(A2);
mxDestroyArray(A2);
mxArray *tmp = Sparse_mult_SAT_SB(A2_t, S1);
mxDestroyArray(B1);
mxArray *tmp;
mexCallMATLAB(1, &tmp, 2, val, "*");
val[0] = B2;
val[1] = tmp;
mexCallMATLAB(1, &B1, 2, val, "-");
B1 = Sparse_substract_SA_SB(B2, tmp);
mxDestroyArray(tmp);
val[0] = A2;
val[1] = d1;
mxDestroyArray(b1);
mexCallMATLAB(1, &tmp, 2, val, "*");
val[0] = b2;
val[1] = tmp;
mexCallMATLAB(1, &b1, 2, val, "-");
tmp = mult_SAT_B(A2_t, d1);
b1 = substract_A_B(b2, tmp);
mxDestroyArray(tmp);
triangular_form.push_back(make_pair(S1, d1));
mxDestroyArray(A2_t);
}
mxDestroyArray(A2);
A2 = mxDuplicateArray(A3);
//I S1
//0 B1 C1 =>B1 =
// A2 B2 => A2 = A3
@ -1544,7 +1846,7 @@ SparseMatrix::Solve_Matlab_Relaxation(mxArray* A_m, mxArray* b_m, unsigned int S
nze--;
while (var < (t+2)*Size && nze < max_nze)
{
if (A_m_j[var+1] <= nze)
if ((unsigned int)A_m_j[var+1] <= nze)
{
b2_d[var - (t+1) * Size] = b_m_d[var];
var++;
@ -1568,7 +1870,6 @@ SparseMatrix::Solve_Matlab_Relaxation(mxArray* A_m, mxArray* b_m, unsigned int S
nze++;
}
}
//mexPrintf("last_t=%d\n", last_t);
double *d1_d = mxGetPr(d1);
for(unsigned i=0; i<Size; i++)
{
@ -1579,18 +1880,15 @@ SparseMatrix::Solve_Matlab_Relaxation(mxArray* A_m, mxArray* b_m, unsigned int S
}
pair<mxArray*, mxArray*> tf;
for (int t = periods-2; t >= 0; t--)
{
mxArray* tmp;
tf = triangular_form.back();
triangular_form.pop_back();
val[0] = tf.first;
val[1] = d1;
mexCallMATLAB(1, &tmp, 2, val, "*");
val[0] = tf.second;
val[1] = tmp;
mexCallMATLAB(1, &d1, 2, val, "-");
mxArray* tf_first_t = Sparse_transpose(tf.first);
mxDestroyArray(tf.first);
tmp = mult_SAT_B(tf_first_t, d1);
d1 = substract_A_B(tf.second, tmp);
d1_d = mxGetPr(d1);
mxDestroyArray(tmp);
for(unsigned i=0; i<Size; i++)
@ -1600,7 +1898,7 @@ SparseMatrix::Solve_Matlab_Relaxation(mxArray* A_m, mxArray* b_m, unsigned int S
direction[eq] = yy;
y[eq] += slowc_l * yy;
}
mxDestroyArray(tf.first);
mxDestroyArray(tf_first_t);
mxDestroyArray(tf.second);
}
mxDestroyArray(B1);
@ -1646,10 +1944,17 @@ SparseMatrix::Solve_Matlab_LU_UMFPack(mxArray* A_m, mxArray* b_m, int Size, doub
}
void
SparseMatrix::Solve_Matlab_GMRES(mxArray* A_m, mxArray* b_m, int Size, double slowc, int block, bool is_two_boundaries, int it_)
SparseMatrix::Solve_Matlab_GMRES(mxArray* A_m, mxArray* b_m, int Size, double slowc, int block, bool is_two_boundaries, int it_, bool steady_state)
{
#ifdef OCTAVE_MEX_FILE
ostringstream tmp;
if (steady_state)
tmp << " GMRES method is not implemented in Octave. You cannot use solve_algo=6, change solve_algo.\n";
else
tmp << " GMRES method is not implemented in Octave. You cannot use stack_solve_algo=2, change stack_solve_algo.\n";
throw FatalExceptionHandling(tmp.str());
#endif
int n = mxGetM(A_m);
/*[L1, U1]=luinc(g1a,luinc_tol);*/
mxArray *lhs0[2];
mxArray *rhs0[2];
rhs0[0] = A_m;
@ -1677,26 +1982,25 @@ SparseMatrix::Solve_Matlab_GMRES(mxArray* A_m, mxArray* b_m, int Size, double sl
mxDestroyArray(rhs[4]);
mxDestroyArray(rhs[5]);
mxDestroyArray(rhs[6]);
if (*flag1 > 0 || reduced)
if (*flag1 > 0)
{
ostringstream tmp;
if (*flag1 == 1)
{
tmp << "Error in simul: No convergence inside GMRES, in block " << block;
tmp << "Error in bytecode: No convergence inside GMRES, in block " << block+1;
mexWarnMsgTxt(tmp.str().c_str());
}
else if (*flag1 == 2)
{
tmp << "Error in simul: Preconditioner is ill-conditioned, in block " << block;
tmp << "Error in bytecode: Preconditioner is ill-conditioned, in block " << block+1;
mexWarnMsgTxt(tmp.str().c_str());
}
else if (*flag1 == 3)
{
tmp << "Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block " << block;
tmp << "Error in bytecode: GMRES stagnated (Two consecutive iterates were the same.), in block " << block+1;
mexWarnMsgTxt(tmp.str().c_str());
}
lu_inc_tol /= 10;
reduced = false;
}
else
{
@ -1714,8 +2018,8 @@ SparseMatrix::Solve_Matlab_GMRES(mxArray* A_m, mxArray* b_m, int Size, double sl
{
int eq = index_vara[i];
double yy = - (res[i] + y[eq+it_*y_size]);
direction[eq] = yy;
y[eq] += slowc * yy;
direction[eq] = yy;
y[eq+it_*y_size] += slowc * yy;
}
}
mxDestroyArray(A_m);
@ -1755,26 +2059,25 @@ SparseMatrix::Solve_Matlab_BiCGStab(mxArray* A_m, mxArray* b_m, int Size, double
mxDestroyArray(rhs[3]);
mxDestroyArray(rhs[4]);
mxDestroyArray(rhs[5]);
if (*flag1 > 0 || reduced)
if (*flag1 > 0)
{
ostringstream tmp;
if (*flag1 == 1)
{
tmp << "Error in simul: No convergence inside BiCGStab, in block " << block;
tmp << "Error in bytecode: No convergence inside BiCGStab, in block " << block+1;
mexWarnMsgTxt(tmp.str().c_str());
}
else if (*flag1 == 2)
{
tmp << "Error in simul: Preconditioner is ill-conditioned, in block " << block;
tmp << "Error in bytecode: Preconditioner is ill-conditioned, in block " << block+1;
mexWarnMsgTxt(tmp.str().c_str());
}
else if (*flag1 == 3)
{
tmp << "Error in simul: BiCGStab stagnated (Two consecutive iterates were the same.), in block " << block;
tmp << "Error in bytecode: BiCGStab stagnated (Two consecutive iterates were the same.), in block " << block+1;
mexWarnMsgTxt(tmp.str().c_str());
}
lu_inc_tol /= 10;
reduced = false;
}
else
{
@ -1793,7 +2096,7 @@ SparseMatrix::Solve_Matlab_BiCGStab(mxArray* A_m, mxArray* b_m, int Size, double
int eq = index_vara[i];
double yy = - (res[i] + y[eq+it_*y_size]);
direction[eq] = yy;
y[eq] += slowc * yy;
y[eq+it_*y_size] += slowc * yy;
}
}
mxDestroyArray(A_m);
@ -2553,7 +2856,7 @@ SparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bool
void
SparseMatrix::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 solve_algo)
SparseMatrix::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)
{
int i, j;
mxArray *b_m = NULL, *A_m = NULL;
@ -2634,7 +2937,41 @@ SparseMatrix::Simulate_Newton_One_Boundary(int blck, int y_size, int it_, int y_
}
if (print_it )
{
mexPrintf("solwc=%f g0=%f res2=%f glambda2=%f\n",slowc_save,g0, res2, glambda2);
//mexPrintf("solwc=%f g0=%f res2=%f glambda2=%f\n",slowc_save,g0, res2, glambda2);
if (steady_state)
{
switch (solve_algo)
{
case 0:
mexPrintf("MODEL STEADY STATE: MATLAB fsolve\n");
break;
case 1:
mexPrintf("MODEL STEADY STATE: MATLAB solve1\n");
break;
case 2:
case 4:
mexPrintf("MODEL STEADY STATE: block decomposition + MATLAB solve1\n");
break;
case 3:
mexPrintf("MODEL STEADY STATE: MATLAB csolve\n");
break;
case 5:
mexPrintf("MODEL STEADY STATE: Sparse LU\n");
break;
case 6:
mexPrintf("MODEL SIMULATION: (method=GMRES)\n");
break;
case 7:
mexPrintf("MODEL SIMULATION: (method=BiCGStab)\n");
break;
case 8:
mexPrintf("MODEL SIMULATION: (method=ByteCode own solver)\n");
break;
default:
mexPrintf("MODEL SIMULATION: (method=Unknown - %d - )\n", stack_solve_algo);
}
}
mexPrintf("-----------------------------------\n");
mexPrintf(" Simulate iteration no %d \n", iter+1);
mexPrintf(" max. error=%.10e \n", double (max_res));
@ -2642,9 +2979,9 @@ SparseMatrix::Simulate_Newton_One_Boundary(int blck, int y_size, int it_, int y_
mexPrintf(" abs. error=%.10e \n", double (res1));
mexPrintf("-----------------------------------\n");
}
if (solve_algo == 8)
Simple_Init(it_, y_kmin, y_kmax, Size, IM_i);
bool zero_solution;
if ((solve_algo == 8 && steady_state) || (stack_solve_algo == 5 && !steady_state))
Simple_Init(it_, y_kmin, y_kmax, Size, IM_i, zero_solution);
else
{
b_m = mxCreateDoubleMatrix(Size,1,mxREAL);
@ -2661,17 +2998,29 @@ SparseMatrix::Simulate_Newton_One_Boundary(int blck, int y_size, int it_, int y_
tmp << " in Simulate_Newton_One_Boundary, can't allocate A_m matrix\n";
throw FatalExceptionHandling(tmp.str());
}
Init_Matlab_Sparse_Simple(Size, IM_i, A_m, b_m);
Init_Matlab_Sparse_Simple(Size, IM_i, A_m, b_m, zero_solution);
}
if (zero_solution)
{
for (int i = 0; i < Size; i++)
{
int eq = index_vara[i];
double yy = - (y[eq+it_*y_size]);
direction[eq] = yy;
y[eq+it_*y_size] += slowc * yy;
}
}
else
{
if ((solve_algo == 8 && steady_state) || (stack_solve_algo == 5 && !steady_state))
Solve_ByteCode_Sparse_GaussianElimination(Size, blck, steady_state, it_);
else if ((solve_algo == 6 && steady_state) || (stack_solve_algo == 2 && !steady_state))
Solve_Matlab_GMRES(A_m, b_m, Size, slowc, blck, false, it_, steady_state);
else if ((solve_algo == 7 && steady_state) || (stack_solve_algo == 3 && !steady_state))
Solve_Matlab_BiCGStab(A_m, b_m, Size, slowc, blck, false, it_);
else if ((solve_algo == 5 && steady_state) || ((stack_solve_algo == 0 || stack_solve_algo == 1) && !steady_state))
Solve_Matlab_LU_UMFPack(A_m, b_m, Size, slowc, false, it_);
}
if (solve_algo == 5)
Solve_Matlab_LU_UMFPack(A_m, b_m, Size, slowc, false, it_);
else if (solve_algo == 6)
Solve_Matlab_GMRES(A_m, b_m, Size, slowc, blck, false, it_);
else if (solve_algo == 7)
Solve_Matlab_BiCGStab(A_m, b_m, Size, slowc, blck, false, it_);
else if (solve_algo == 8)
Solve_ByteCode_Sparse_GaussianElimination(Size, blck, steady_state, it_);
return;
}
@ -2701,8 +3050,8 @@ SparseMatrix::Simulate_Newton_Two_Boundaries(int blck, int y_size, int it_, int
{
ostringstream res;
for (unsigned int i = 0; i < endo_name_length; i++)
if (P_endo_names[2*(j+i*y_size)] != ' ')
res << P_endo_names[2*(j+i*y_size)];
if (P_endo_names[CHAR_LENGTH*(j+i*y_size)] != ' ')
res << P_endo_names[CHAR_LENGTH*(j+i*y_size)];
bool select = false;
for (int i = 0; i < Size; i++)
if (j == index_vara[i])
@ -2875,7 +3224,7 @@ SparseMatrix::Simulate_Newton_Two_Boundaries(int blck, int y_size, int it_, int
else if (stack_solve_algo == 1)
Solve_Matlab_Relaxation(A_m, b_m, Size, slowc, true, 0);
else if (stack_solve_algo == 2)
Solve_Matlab_GMRES(A_m, b_m, Size, slowc, blck, true, 0);
Solve_Matlab_GMRES(A_m, b_m, Size, slowc, blck, true, 0, false);
else if (stack_solve_algo == 3)
Solve_Matlab_BiCGStab(A_m, b_m, Size, slowc, blck, true, 0);
else if (stack_solve_algo == 5)

26
mex/sources/bytecode/SparseMatrix.hh
View File

@ -25,6 +25,14 @@
#include <cmath>
#include <map>
#include <ctime>
#ifdef OCTAVE_MEX_FILE
#define CHAR_LENGTH 1
#else
#define CHAR_LENGTH 2
#endif
#include "Mem_Mngr.hh"
#include "ErrorHandling.hh"
#define NEW_ALLOC
@ -56,7 +64,7 @@ class SparseMatrix
public:
SparseMatrix();
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)*/;
void 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 solve_algo);
void 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);
void Direct_Simulate(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it, int iter);
void fixe_u(double **u, int u_count_int, int max_lag_plus_max_lead_plus_1);
void Read_SparseMatrix(string file_name, const int Size, int periods, int y_kmin, int y_kmax, bool steady_state, bool two_boundaries, int stack_solve_algo, int solve_algo);
@ -66,14 +74,14 @@ public:
private:
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);
void Init_Matlab_Sparse_Simple(int Size, map<pair<pair<int, int>, int>, int> &IM, mxArray *A_m, mxArray *b_m);
void Simple_Init(int it_, int y_kmin, int y_kmax, int Size, std::map<std::pair<std::pair<int, int>, int>, int> &IM);
void Init_Matlab_Sparse_Simple(int Size, map<pair<pair<int, int>, int>, int> &IM, mxArray *A_m, mxArray *b_m, bool &zero_solution);
void Simple_Init(int it_, int y_kmin, int y_kmax, int Size, std::map<std::pair<std::pair<int, int>, int>, int> &IM, bool &zero_solution);
void End_GE(int Size);
void Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bool symbolic, int Block_number);
void Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, bool steady_state, 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_GMRES(mxArray* A_m, mxArray* b_m, int Size, double slowc, int block, 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);
void Solve_Matlab_BiCGStab(mxArray* A_m, mxArray* b_m, int Size, double slowc, int block, bool is_two_boundaries, int it_);
bool compare(int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, long int nop4, int Size
#ifdef PROFILER
@ -103,6 +111,15 @@ private:
#endif
);
double simple_bksub(int it_, int Size, double slowc_l);
mxArray* Sparse_transpose(mxArray* A_m);
mxArray* Sparse_mult_SAT_SB(mxArray* A_m, mxArray* B_m);
mxArray* Sparse_mult_SAT_B(mxArray* A_m, mxArray* B_m);
mxArray* mult_SAT_B(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* substract_A_B(mxArray* A_m, mxArray* B_m);
stack<double> Stack;
int nb_prologue_table_u, nb_first_table_u, nb_middle_table_u, nb_last_table_u;
int nb_prologue_table_y, nb_first_table_y, nb_middle_table_y, nb_last_table_y;
@ -152,7 +169,6 @@ protected:
bool error_not_printed;
double g_lambda1, g_lambda2, gp_0;
double lu_inc_tol;
bool reduced;
};
#endif

11
mex/sources/bytecode/testing/mex_interface.cc
View File

@ -131,7 +131,7 @@ mxCreateSparse(unsigned int rows, unsigned int cols, unsigned int nz_max, mxData
Array->Nzmax = nz_max;
Array->data = (double*)mxMalloc(nz_max*sizeof(double));
Array->Ir = (mwIndex*)mxMalloc(nz_max*sizeof(mwIndex));
Array->Jc = (mwIndex*)mxMalloc(cols*sizeof(mwIndex));
Array->Jc = (mwIndex*)mxMalloc((cols+1)*sizeof(mwIndex));
return(Array);
}
@ -176,6 +176,15 @@ mxDestroyArray(mxArray* A_m)
mxFree(A_m);
}
mxArray*
mxSetNzmax(mxArray* A_m, mwSize nz_max)
{
A_m->Nzmax = nz_max;
A_m->data = (double*)mxRealloc(A_m->data, nz_max*sizeof(double));
A_m->Ir = (mwIndex*)mxRealloc(A_m->Ir, nz_max*sizeof(mwIndex));
return(A_m);
}
void
mexCallMATLAB(unsigned int n_lhs, mxArray* matrix_lhs[], unsigned int n_rhs, mxArray* matrix_rhs[], const char* function)
{

1
mex/sources/bytecode/testing/mex_interface.hh
View File

@ -93,6 +93,7 @@ void mexEvalString(const string str);
double* mxGetPr(const mxArray *b_m);
inline mwIndex* mxGetIr(const mxArray *A_m) {return (mwIndex*)A_m->Ir;};
inline mwIndex* mxGetJc(const mxArray *A_m) {return (mwIndex*)A_m->Jc;};
mxArray* mxSetNzmax(mxArray* A_m, mwSize nz_max);
inline int mxGetNzmax(const mxArray *A_m) {return A_m->Nzmax;};
inline int mxGetM(const mxArray *A_m) {return A_m->size_1;};
inline int mxGetN(const mxArray *A_m) {return A_m->size_2;};

51
preprocessor/DynamicModel.cc
View File

@ -2231,21 +2231,21 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
for (lag_var_t::const_iterator it = other_endo_block[block].begin(); it != other_endo_block[block].end(); it++)
for(var_t::const_iterator it1 = it->second.begin(); it1 != it->second.end(); it1++)
other_endogenous.insert(*it1);
output << "M_.block_structure.block(" << block+1 << ").Simulation_Type = " << simulation_type << ";\n";
output << "M_.block_structure.block(" << block+1 << ").maximum_lag = " << max_lag << ";\n";
output << "M_.block_structure.block(" << block+1 << ").maximum_lead = " << max_lead << ";\n";
output << "M_.block_structure.block(" << block+1 << ").maximum_endo_lag = " << max_lag_endo << ";\n";
output << "M_.block_structure.block(" << block+1 << ").maximum_endo_lead = " << max_lead_endo << ";\n";
output << "M_.block_structure.block(" << block+1 << ").maximum_exo_lag = " << max_lag_exo << ";\n";
output << "M_.block_structure.block(" << block+1 << ").maximum_exo_lead = " << max_lead_exo << ";\n";
output << "M_.block_structure.block(" << block+1 << ").maximum_exo_det_lag = " << max_lag_exo_det << ";\n";
output << "M_.block_structure.block(" << block+1 << ").maximum_exo_det_lead = " << max_lead_exo_det << ";\n";
output << "M_.block_structure.block(" << block+1 << ").endo_nbr = " << block_size << ";\n";
output << "M_.block_structure.block(" << block+1 << ").mfs = " << getBlockMfs(block) << ";\n";
output << "M_.block_structure.block(" << block+1 << ").equation = [" << tmp_s_eq.str() << "];\n";
output << "M_.block_structure.block(" << block+1 << ").variable = [" << tmp_s.str() << "];\n";
output << "M_.block_structure.block(" << block+1 << ").exo_nbr = " << getBlockExoSize(block) << ";\n";
output << "M_.block_structure.block(" << block+1 << ").exogenous = [";
output << "block_structure.block(" << block+1 << ").Simulation_Type = " << simulation_type << ";\n";
output << "block_structure.block(" << block+1 << ").maximum_lag = " << max_lag << ";\n";
output << "block_structure.block(" << block+1 << ").maximum_lead = " << max_lead << ";\n";
output << "block_structure.block(" << block+1 << ").maximum_endo_lag = " << max_lag_endo << ";\n";
output << "block_structure.block(" << block+1 << ").maximum_endo_lead = " << max_lead_endo << ";\n";
output << "block_structure.block(" << block+1 << ").maximum_exo_lag = " << max_lag_exo << ";\n";
output << "block_structure.block(" << block+1 << ").maximum_exo_lead = " << max_lead_exo << ";\n";
output << "block_structure.block(" << block+1 << ").maximum_exo_det_lag = " << max_lag_exo_det << ";\n";
output << "block_structure.block(" << block+1 << ").maximum_exo_det_lead = " << max_lead_exo_det << ";\n";
output << "block_structure.block(" << block+1 << ").endo_nbr = " << block_size << ";\n";
output << "block_structure.block(" << block+1 << ").mfs = " << getBlockMfs(block) << ";\n";
output << "block_structure.block(" << block+1 << ").equation = [" << tmp_s_eq.str() << "];\n";
output << "block_structure.block(" << block+1 << ").variable = [" << tmp_s.str() << "];\n";
output << "block_structure.block(" << block+1 << ").exo_nbr = " << getBlockExoSize(block) << ";\n";
output << "block_structure.block(" << block+1 << ").exogenous = [";
int i = 0;
for (set<int>::iterator it_exogenous = exogenous.begin(); it_exogenous != exogenous.end(); it_exogenous++)
if (*it_exogenous >= 0)
@ -2254,8 +2254,8 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
i++;
}
output << "];\n";
output << "M_.block_structure.block(" << block+1 << ").exo_det_nbr = " << i << ";\n";
output << "M_.block_structure.block(" << block+1 << ").exogenous_det = [";
output << "block_structure.block(" << block+1 << ").exo_det_nbr = " << i << ";\n";
output << "block_structure.block(" << block+1 << ").exogenous_det = [";
i = 0;
for (set<int>::iterator it_exogenous_det = exogenous_det.begin(); it_exogenous_det != exogenous_det.end(); it_exogenous_det++)
if (*it_exogenous_det >= 0)
@ -2265,8 +2265,8 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
}
output << "];\n";
output << "M_.block_structure.block(" << block+1 << ").other_endo_nbr = " << i << ";\n";
output << "M_.block_structure.block(" << block+1 << ").other_endogenous = [";
output << "block_structure.block(" << block+1 << ").other_endo_nbr = " << i << ";\n";
output << "block_structure.block(" << block+1 << ").other_endogenous = [";
i = 0;
for (set<int>::iterator it_other_endogenous = other_endogenous.begin(); it_other_endogenous != other_endogenous.end(); it_other_endogenous++)
if (*it_other_endogenous >= 0)
@ -2281,7 +2281,7 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
dynamic_jacob_map_t reordered_dynamic_jacobian;
for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != blocks_derivatives[block].end(); it++)
reordered_dynamic_jacobian[make_pair(it->second.first, make_pair(it->first.second, it->first.first))] = it->second.second;
output << "M_.block_structure.block(" << block+1 << ").lead_lag_incidence = [];\n";
output << "block_structure.block(" << block+1 << ").lead_lag_incidence = [];\n";
int last_var = -1;
for (int lag = -max_lag_endo; lag < max_lead_endo+1; lag++)
{
@ -2301,14 +2301,15 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
}
for (int i = last_var + 1; i < block_size; i++)
tmp_s << " 0";
output << "M_.block_structure.block(" << block+1 << ").lead_lag_incidence = [ M_.block_structure.block(" << block+1 << ").lead_lag_incidence; " << tmp_s.str() << "]; %lag = " << lag << "\n";
output << "block_structure.block(" << block+1 << ").lead_lag_incidence = [ block_structure.block(" << block+1 << ").lead_lag_incidence; " << tmp_s.str() << "]; %lag = " << lag << "\n";
tmp_s.str("");
}
output << "M_.block_structure.block(" << block+1 << ").n_static = " << block_col_type[block].first.first << ";\n";
output << "M_.block_structure.block(" << block+1 << ").n_forward = " << block_col_type[block].first.second << ";\n";
output << "M_.block_structure.block(" << block+1 << ").n_backward = " << block_col_type[block].second.first << ";\n";
output << "M_.block_structure.block(" << block+1 << ").n_mixed = " << block_col_type[block].second.second << ";\n";
output << "block_structure.block(" << block+1 << ").n_static = " << block_col_type[block].first.first << ";\n";
output << "block_structure.block(" << block+1 << ").n_forward = " << block_col_type[block].first.second << ";\n";
output << "block_structure.block(" << block+1 << ").n_backward = " << block_col_type[block].second.first << ";\n";
output << "block_structure.block(" << block+1 << ").n_mixed = " << block_col_type[block].second.second << ";\n";
}
output << "M_.block_structure.block = block_structure.block;\n";
string cst_s;
int nb_endo = symbol_table.endo_nbr();
output << "M_.block_structure.variable_reordered = [";

53
preprocessor/StaticModel.cc
View File

@ -465,6 +465,12 @@ StaticModel::writeModelEquationsCode(const string file_name, const string bin_ba
FENDEQU_ fendequ;
fendequ.write(code_file, instruction_number);
// Get the current code_file position and jump if eval = true
streampos pos1 = code_file.tellp();
FJMPIFEVAL_ fjmp_if_eval(0);
fjmp_if_eval.write(code_file, instruction_number);
int prev_instruction_number = instruction_number;
vector<vector<pair<int, int> > > derivatives;
derivatives.resize(symbol_table.endo_nbr());
count_u = symbol_table.endo_nbr();
@ -515,6 +521,53 @@ StaticModel::writeModelEquationsCode(const string file_name, const string bin_ba
FSTPSU_ fstpsu(i);
fstpsu.write(code_file, instruction_number);
}
// Get the current code_file position and jump = true
streampos pos2 = code_file.tellp();
FJMP_ fjmp(0);
fjmp.write(code_file, instruction_number);
// Set code_file position to previous JMPIFEVAL_ and set the number of instructions to jump
streampos pos3 = code_file.tellp();
code_file.seekp(pos1);
FJMPIFEVAL_ fjmp_if_eval1(instruction_number - prev_instruction_number);
fjmp_if_eval1.write(code_file, instruction_number);
code_file.seekp(pos3);
prev_instruction_number = instruction_number ;
temporary_terms_t tt2;
tt2.clear();
temporary_terms_t tt3;
tt3.clear();
// The Jacobian if we have to solve the block determinsitic bloc
for (first_derivatives_t::const_iterator it = first_derivatives.begin();
it != first_derivatives.end(); it++)
{
int deriv_id = it->first.second;
if (getTypeByDerivID(deriv_id) == eEndogenous)
{
expr_t d1 = it->second;
unsigned int eq = it->first.first;
int symb = getSymbIDByDerivID(deriv_id);
unsigned int var = symbol_table.getTypeSpecificID(symb);
FNUMEXPR_ fnumexpr(FirstEndoDerivative, eq, var);
fnumexpr.write(code_file, instruction_number);
if (!derivatives[eq].size())
derivatives[eq].clear();
derivatives[eq].push_back(make_pair(var, count_u));
d1->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
FSTPG2_ fstpg2(eq,var);
fstpg2.write(code_file, instruction_number);
}
}
// Set codefile position to previous JMP_ and set the number of instructions to jump
pos1 = code_file.tellp();
code_file.seekp(pos2);
FJMP_ fjmp1(instruction_number - prev_instruction_number);
fjmp1.write(code_file, instruction_number);
code_file.seekp(pos1);
FENDBLOCK_ fendblock;
fendblock.write(code_file, instruction_number);
FEND_ fend;

2
tests/block_bytecode/fs2000_block_bytecode.mod
View File

@ -3,6 +3,6 @@
@#define bytecode = 1
@#define block = 1
@#define solve_algo = 5
@#define solve_algo = 8
@#define stack_solve_algo = 5
@#include "fs2000_common.mod"

16
tests/block_bytecode/fs2000_bytecode.mod
View File

@ -1,8 +1,8 @@
// Tests option bytecode without block
// Must be launched after fs2000_simk.mod
@#define bytecode = 1
@#define block = 0
@#define solve_algo = 5
@#define stack_solve_algo = 5
@#include "fs2000_common.mod"
// Tests option bytecode without block
// Must be launched after fs2000_simk.mod
@#define bytecode = 1
@#define block = 0
@#define solve_algo = 8
@#define stack_solve_algo = 5
@#include "fs2000_common.mod"

16
tests/block_bytecode/fs2000_gmres.mod
View File

@ -1,8 +1,8 @@
// Tests option block + stack_solve_algo = 2 + solve_algo = 3
// Must be launched after fs2000_simk.mod
@#define block = 1
@#define bytecode = 0
@#define solve_algo = 3
@#define stack_solve_algo = 2
@#include "fs2000_common.mod"
// Tests option block + stack_solve_algo = 2 + solve_algo = 3
// Must be launched after fs2000_simk.mod
@#define block = 1
@#define bytecode = 0
@#define solve_algo = 2
@#define stack_solve_algo = 2
@#include "fs2000_common.mod"

16
tests/block_bytecode/fs2000_lu.mod
View File

@ -1,8 +1,8 @@
// Tests option block + stack_solve_algo = 1 + solve_algo = 2
// Must be launched after fs2000_simk.mod
@#define block = 1
@#define bytecode = 0
@#define solve_algo = 2
@#define stack_solve_algo = 1
@#include "fs2000_common.mod"
// Tests option block + stack_solve_algo = 1 + solve_algo = 2
// Must be launched after fs2000_simk.mod
@#define block = 1
@#define bytecode = 0
@#define solve_algo = 2
@#define stack_solve_algo = 0
@#include "fs2000_common.mod"

2
tests/block_bytecode/ireland.mod
View File

@ -23,7 +23,7 @@ scy = 0.0040;
shy = 0.0015;
shc = 0.0010;
model(block,cutoff=0);
model(bytecode, block,cutoff=0);
exp(y) = exp(a)*exp(k(-1))^theta*exp(h)^(1-theta);
a = (1-rho)*aa+rho*a(-1)+e;
exp(y) = exp(c) + exp(i);

2
tests/block_bytecode/mod1.mod
View File

@ -18,7 +18,7 @@ var e_x; stderr 0.01;
var e_y; stderr 0.01;
end;
steady;
check;
//check;
model_info;
simul(periods=50);
/*stoch_simul(order=1,periods=1000,irf=0,nomoments);