From 881d2819a3698fa581f3e34581c0c3e68d6f689d Mon Sep 17 00:00:00 2001 From: ferhat Date: Wed, 24 Dec 2008 16:29:51 +0000 Subject: [PATCH] - 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 --- matlab/dr1_sparse.m | 13 ++++++++++--- matlab/solve_one_boundary.m | 8 +++++--- matlab/solve_two_boundaries.m | 7 +++++-- 3 files changed, 20 insertions(+), 8 deletions(-) diff --git a/matlab/dr1_sparse.m b/matlab/dr1_sparse.m index 32dc13543..568a0cec9 100644 --- a/matlab/dr1_sparse.m +++ b/matlab/dr1_sparse.m @@ -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); diff --git a/matlab/solve_one_boundary.m b/matlab/solve_one_boundary.m index 0a72135ac..ed5a34ad1 100644 --- a/matlab/solve_one_boundary.m +++ b/matlab/solve_one_boundary.m @@ -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; diff --git a/matlab/solve_two_boundaries.m b/matlab/solve_two_boundaries.m index a652caa03..d9780c6c4 100644 --- a/matlab/solve_two_boundaries.m +++ b/matlab/solve_two_boundaries.m @@ -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