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

- Jacobian matrix stored only for blocks -> jacobian matrix only as output argument 

git-svn-id: https://www.dynare.org/svn/dynare/trunk@2342 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
ferhat 2008-12-24 16:29:51 +00:00
parent 487afce68d
commit 881d2819a3
3 changed files with 20 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
matlab/dr1_sparse.m
View File

@ -209,7 +209,7 @@ function [dr,info,M_,options_,oo_] = dr1_sparse(dr,task,M_,options_,oo_)
elseif options_.model_mode==1
if options_.order == 1
[junk,jacobia_] = feval([M_.fname '_dynamic'],ones(M_.maximum_lag+M_.maximum_lead+1,1)*dr.ys',[oo_.exo_simul ...
[junk,derivate] = feval([M_.fname '_dynamic'],ones(M_.maximum_lag+M_.maximum_lead+1,1)*dr.ys',[oo_.exo_simul ...
oo_.exo_det_simul], M_.params, it_);
%full(jacobia_)
dr.eigval = [];
@ -222,14 +222,21 @@ function [dr,info,M_,options_,oo_] = dr1_sparse(dr,task,M_,options_,oo_)
M_.block_structure.block(i).dr=set_state_space(M_.block_structure.block(i).dr,M_.block_structure.block(i));
col_selector=repmat(M_.block_structure.block(i).variable,1,M_.block_structure.block(i).maximum_endo_lag+M_.block_structure.block(i).maximum_endo_lead+1)+kron([M_.maximum_endo_lag-M_.block_structure.block(i).maximum_endo_lag:M_.maximum_endo_lag+M_.block_structure.block(i).maximum_endo_lead],M_.endo_nbr*ones(1,M_.block_structure.block(i).endo_nbr));
row_selector = M_.block_structure.block(i).equation;
jcb_=jacobia_(row_selector,col_selector);
%jcb_=jacobia_(row_selector,col_selector);
jcb_=derivate(i).g1;
%disp('jcb_');
%full(jcb_)
%M_.block_structure.block(i).lead_lag_incidence'
jcb_ = jcb_(:,find(M_.block_structure.block(i).lead_lag_incidence')) ;
if M_.block_structure.block(i).exo_nbr>0
col_selector = [ first_col_exo + ...
repmat(M_.block_structure.block(i).exogenous,1,M_.block_structure.block(i).maximum_exo_lag+M_.block_structure.block(i).maximum_exo_lead+1)+kron([M_.maximum_exo_lag-M_.block_structure.block(i).maximum_exo_lag:M_.maximum_exo_lag+M_.block_structure.block(i).maximum_exo_lead],M_.exo_nbr*ones(1,M_.block_structure.block(i).exo_nbr))];
end
%derivate(i).g1
%derivate(i).g1_x
%col_selector
jcb_ = [ jcb_ jacobia_(row_selector,col_selector)];
%jcb_ = [ jcb_ jacobia_(row_selector,col_selector)];
jcb_ = [ jcb_ derivate(i).g1_x];
%full(jcb_)
hss_=0; %hessian(M_.block_structure.block(i).equation,M_.block_structure.block(i).variable);

8
matlab/solve_one_boundary.m
View File

@ -42,7 +42,7 @@ function y = solve_one_boundary(fname, y, x, params, y_index_eq, nze, periods, i
% y [matrix] All endogenous variables of the model
%
% ALGORITHM
% Newton with LU or GMRES or BicGstab
% Newton with LU or GMRES or BicGstab for dynamic block
%
% SPECIAL REQUIREMENTS
% none.
@ -90,9 +90,9 @@ function y = solve_one_boundary(fname, y, x, params, y_index_eq, nze, periods, i
g1=spalloc( Blck_size, Blck_size, nze);
while ~(cvg==1 | iter>maxit_),
if(is_dynamic)
[r, g1, g2, g3, b] = feval(fname, y, x, params, it_, 0, g1, g2, g3);
[r, g1, g2, g3] = feval(fname, y, x, params, it_, 0);
else
[r, g1, g2, g3, b] = feval(fname, y, x, params, 0);
[r, g1, g2, g3] = feval(fname, y, x, params, 0);
end;
if(~isreal(r))
max_res=(-(max(max(abs(r))))^2)^0.5;
@ -204,8 +204,10 @@ function y = solve_one_boundary(fname, y, x, params, y_index_eq, nze, periods, i
f = 0.5*r'*r;
p = -g1\r ;
[y,f,r,check]=lnsrch1(y,f,g,p,stpmax,fname,nn,y_index_eq,x, params, 0);
dx = ya - y(y_index_eq);
elseif(~is_dynamic & options_.solve_algo==3)
[yn,info] = csolve(@local_fname, y(y_index_eq),@local_fname,1e-6,500, x, params, y, y_index_eq, fname, 1);
dx = ya - yn;
y(y_index_eq) = yn;
elseif((simulation_method==0 & is_dynamic) | (~is_dynamic & options_.solve_algo==1)),
dx = g1\r;

7
matlab/solve_two_boundaries.m
View File

@ -72,9 +72,12 @@ function y = solve_two_boundaries(fname, y, x, params, y_index, nze, periods, y_
g1=spalloc( Blck_size*periods, Jacobian_Size, nze*periods);
reduced = 0;
while ~(cvg==1 | iter>maxit_),
[r, g1, g2, g3, b]=feval(fname, y, x, params, periods, 0, g1, g2, g3, y_kmin, Blck_size);
[r, g1, g2, g3, b]=feval(fname, y, x, params, periods, 0, y_kmin, Blck_size);
g1a=g1(:, y_kmin*Blck_size+1:(periods+y_kmin)*Blck_size);
b = b' -g1(:, 1+(y_kmin-y_kmin_l)*Blck_size:y_kmin*Blck_size)*reshape(y(1+y_kmin-y_kmin_l:y_kmin,y_index)',1,y_kmin_l*Blck_size)'-g1(:, (periods+y_kmin)*Blck_size+1:(periods+y_kmin+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)';
%disp(['size(g1)=' int2str(size(g1))]);
%disp(['g1(:,' int2str(1:y_kmin_l*Blck_size) ')']);
%disp(['g1(:,' int2str((periods+y_kmin_l)*Blck_size+1:(periods+y_kmin_l+y_kmax_l)*Blck_size) ')']);
b = b' -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)'-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)';
if(~isreal(r))
max_res=(-(max(max(abs(r))))^2)^0.5;
else