Fixed homotopy for perfect foresight models with block option, cosmetic changes.
Stéphane Adjemian (Charybdis)
parent
1faf755d56
commit
d19592f761
|
|
@ -110,10 +110,14 @@ if nargout>1
|
|||
[i_cols_J1,~,i_cols_1] = find(illi(:));
|
||||
i_cols_T = nonzeros(M_.lead_lag_incidence(1:2,:)');
|
||||
end
|
||||
residuals = perfect_foresight_problem(yy(:),str2func([M_.fname '_dynamic']), y0, yT, ...
|
||||
oo_.exo_simul,M_.params,oo_.steady_state, ...
|
||||
options_.periods,M_.endo_nbr,i_cols, ...
|
||||
i_cols_J1, i_cols_1, i_cols_T, i_cols_j, ...
|
||||
M_.NNZDerivatives(1));
|
||||
maxerror = max(max(abs(residuals)));
|
||||
if options_.block && ~options_.bytecode
|
||||
maxerror = oo_.deterministic_simulation.error;
|
||||
else
|
||||
residuals = perfect_foresight_problem(yy(:),str2func([M_.fname '_dynamic']), y0, yT, ...
|
||||
oo_.exo_simul,M_.params,oo_.steady_state, ...
|
||||
options_.periods,M_.endo_nbr,i_cols, ...
|
||||
i_cols_J1, i_cols_1, i_cols_T, i_cols_j, ...
|
||||
M_.NNZDerivatives(1));
|
||||
maxerror = max(max(abs(residuals)));
|
||||
end
|
||||
end
|
||||
|
|
|
|||
|
|
@ -80,7 +80,7 @@ correcting_factor=0.01;
|
|||
ilu_setup.droptol=1e-10;
|
||||
max_resa=1e100;
|
||||
reduced = 0;
|
||||
if(forward_backward)
|
||||
if forward_backward
|
||||
incr = 1;
|
||||
start = y_kmin+1;
|
||||
finish = periods+y_kmin;
|
||||
|
|
@ -89,140 +89,113 @@ else
|
|||
start = periods+y_kmin;
|
||||
finish = y_kmin+1;
|
||||
end
|
||||
%lambda=1;
|
||||
|
||||
for it_=start:incr:finish
|
||||
cvg=0;
|
||||
iter=0;
|
||||
g1=spalloc( Blck_size, Blck_size, nze);
|
||||
while ~(cvg==1 || iter>maxit_),
|
||||
if(is_dynamic)
|
||||
[r, y, g1, g2, g3] = feval(fname, y, x, params, steady_state, ...
|
||||
it_, 0);
|
||||
while ~(cvg==1 || iter>maxit_)
|
||||
if is_dynamic
|
||||
[r, y, g1, g2, g3] = feval(fname, y, x, params, steady_state, it_, 0);
|
||||
else
|
||||
[r, y, g1] = feval(fname, y, x, params);
|
||||
end;
|
||||
if(~isreal(r))
|
||||
end
|
||||
if ~isreal(r)
|
||||
max_res=(-(max(max(abs(r))))^2)^0.5;
|
||||
else
|
||||
max_res=max(max(abs(r)));
|
||||
end;
|
||||
%['max_res=' num2str(max_res) ' Block_Num=' int2str(Block_Num) ' it_=' int2str(it_)]
|
||||
%disp(['iteration : ' int2str(iter+1) ' => ' num2str(max_res) ' time = ' int2str(it_)]);
|
||||
|
||||
% fjac = zeros(Blck_size, Blck_size);
|
||||
% disp(['Blck_size=' int2str(Blck_size) ' it_=' int2str(it_)]);
|
||||
% dh = max(abs(y(it_, y_index_eq)),options_.gstep*ones(1, Blck_size))*eps^(1/3);
|
||||
% fvec = r;
|
||||
% for j = 1:Blck_size
|
||||
% ydh = y ;
|
||||
% ydh(it_, y_index_eq(j)) = ydh(it_, y_index_eq(j)) + dh(j) ;
|
||||
% [t, y1, g11, g21, g31]=feval(fname, ydh, x, params, it_, 0);
|
||||
% fjac(:,j) = (t - fvec)./dh(j) ;
|
||||
% end;
|
||||
% diff = g1 -fjac;
|
||||
% diff
|
||||
% disp('g1');
|
||||
% disp([num2str(g1,'%4.5f')]);
|
||||
% disp('fjac');
|
||||
% disp([num2str(fjac,'%4.5f')]);
|
||||
% [c_max, i_c_max] = max(abs(diff));
|
||||
% [l_c_max, i_r_max] = max(c_max);
|
||||
% disp(['maximum element row=' int2str(i_c_max(i_r_max)) ' and column=' int2str(i_r_max) ' value = ' num2str(l_c_max)]);
|
||||
% equation = i_c_max(i_r_max);
|
||||
% variable = i_r_max;
|
||||
% variable
|
||||
% mod(variable, Blck_size)
|
||||
% disp(['equation ' int2str(equation) ' and variable ' int2str(y_index_eq(variable)) ' ' M_.endo_names(y_index_eq(variable), :)]);
|
||||
% disp(['g1(' int2str(equation) ', ' int2str(variable) ')=' num2str(g1(equation, variable),'%3.10f') ' fjac(' int2str(equation) ', ' int2str(variable) ')=' num2str(fjac(equation, variable), '%3.10f') ' y(' int2str(it_) ', ' int2str(variable) ')=' num2str(y(it_, variable))]);
|
||||
% %return;
|
||||
% %g1 = fjac;
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
if(verbose==1)
|
||||
disp(['iteration : ' int2str(iter+1) ' => ' num2str(max_res) ' time = ' int2str(it_)]);
|
||||
if(is_dynamic)
|
||||
disp([M.endo_names(y_index_eq,:) num2str([y(it_,y_index_eq)' r g1])]);
|
||||
end
|
||||
if verbose==1
|
||||
disp(['iteration : ' int2str(iter+1) ' => ' num2str(max_res) ' time = ' int2str(it_)])
|
||||
if is_dynamic
|
||||
disp([M.endo_names(y_index_eq,:) num2str([y(it_,y_index_eq)' r g1])])
|
||||
else
|
||||
disp([M.endo_names(y_index_eq,:) num2str([y(y_index_eq) r g1])]);
|
||||
end;
|
||||
end;
|
||||
if(~isreal(max_res) || isnan(max_res))
|
||||
disp([M.endo_names(y_index_eq,:) num2str([y(y_index_eq) r g1])])
|
||||
end
|
||||
end
|
||||
if ~isreal(max_res) || isnan(max_res)
|
||||
cvg = 0;
|
||||
elseif(is_linear && iter>0)
|
||||
elseif is_linear && iter>0
|
||||
cvg = 1;
|
||||
else
|
||||
cvg=(max_res<solve_tolf);
|
||||
end;
|
||||
if(~cvg)
|
||||
if(iter>0)
|
||||
if(~isreal(max_res) || isnan(max_res) || (max_resa<max_res && iter>1))
|
||||
if(isnan(max_res) || (max_resa<max_res && iter>0))
|
||||
end
|
||||
if ~cvg
|
||||
if iter>0
|
||||
if ~isreal(max_res) || isnan(max_res) || (max_resa<max_res && iter>1)
|
||||
if isnan(max_res) || (max_resa<max_res && iter>0)
|
||||
detJ=det(g1a);
|
||||
if(abs(detJ)<1e-7)
|
||||
max_factor=max(max(abs(g1a)));
|
||||
ze_elem=sum(diag(g1a)<cutoff);
|
||||
disp([num2str(full(ze_elem),'%d') ' elements on the Jacobian diagonal are below the cutoff (' num2str(cutoff,'%f') ')']);
|
||||
if(correcting_factor<max_factor)
|
||||
if verbose
|
||||
disp([num2str(full(ze_elem),'%d') ' elements on the Jacobian diagonal are below the cutoff (' num2str(cutoff,'%f') ')'])
|
||||
end
|
||||
if correcting_factor<max_factor
|
||||
correcting_factor=correcting_factor*4;
|
||||
disp(['The Jacobain matrix is singular, det(Jacobian)=' num2str(detJ,'%f') '.']);
|
||||
disp([' trying to correct the Jacobian matrix:']);
|
||||
disp([' correcting_factor=' num2str(correcting_factor,'%f') ' max(Jacobian)=' num2str(full(max_factor),'%f')]);
|
||||
if verbose
|
||||
disp(['The Jacobain matrix is singular, det(Jacobian)=' num2str(detJ,'%f') '.'])
|
||||
disp([' trying to correct the Jacobian matrix:'])
|
||||
disp([' correcting_factor=' num2str(correcting_factor,'%f') ' max(Jacobian)=' num2str(full(max_factor),'%f')])
|
||||
end
|
||||
dx = - r/(g1+correcting_factor*speye(Blck_size));
|
||||
%dx = -b'/(g1+correcting_factor*speye(Blck_size))-ya_save;
|
||||
y(it_,y_index_eq)=ya_save+lambda*dx;
|
||||
continue;
|
||||
continue
|
||||
else
|
||||
disp('The singularity of the jacobian matrix could not be corrected');
|
||||
if verbose
|
||||
disp('The singularity of the jacobian matrix could not be corrected')
|
||||
end
|
||||
info = -Block_Num*10;
|
||||
return;
|
||||
end;
|
||||
end;
|
||||
elseif(lambda>1e-8)
|
||||
return
|
||||
end
|
||||
end
|
||||
elseif lambda>1e-8
|
||||
lambda=lambda/2;
|
||||
reduced = 1;
|
||||
disp(['reducing the path length: lambda=' num2str(lambda,'%f')]);
|
||||
if(is_dynamic)
|
||||
if verbose
|
||||
disp(['reducing the path length: lambda=' num2str(lambda,'%f')])
|
||||
end
|
||||
if is_dynamic
|
||||
y(it_,y_index_eq)=ya_save-lambda*dx;
|
||||
else
|
||||
y(y_index_eq)=ya_save-lambda*dx;
|
||||
end;
|
||||
continue;
|
||||
end
|
||||
continue
|
||||
else
|
||||
if(cutoff == 0)
|
||||
fprintf('Error in simul: Convergence not achieved in block %d, at time %d, after %d iterations.\n Increase "options_.simul.maxit".\n',Block_Num, it_, iter);
|
||||
else
|
||||
fprintf('Error in simul: Convergence not achieved in block %d, at time %d, after %d iterations.\n Increase "options_.simul.maxit" or set "cutoff=0" in model options.\n',Block_Num, it_, iter);
|
||||
end;
|
||||
if(is_dynamic)
|
||||
if verbose
|
||||
if cutoff==0
|
||||
fprintf('Error in simul: Convergence not achieved in block %d, at time %d, after %d iterations.\n Increase "options_.simul.maxit".\n',Block_Num, it_, iter);
|
||||
else
|
||||
fprintf('Error in simul: Convergence not achieved in block %d, at time %d, after %d iterations.\n Increase "options_.simul.maxit" or set "cutoff=0" in model options.\n',Block_Num, it_, iter);
|
||||
end
|
||||
end
|
||||
if is_dynamic
|
||||
oo_.deterministic_simulation.status = 0;
|
||||
oo_.deterministic_simulation.error = max_res;
|
||||
oo_.deterministic_simulation.iterations = iter;
|
||||
oo_.deterministic_simulation.block(Block_Num).status = 0;% Convergency failed.
|
||||
oo_.deterministic_simulation.block(Block_Num).error = max_res;
|
||||
oo_.deterministic_simulation.block(Block_Num).iterations = iter;
|
||||
end;
|
||||
end
|
||||
info = -Block_Num*10;
|
||||
return;
|
||||
end;
|
||||
return
|
||||
end
|
||||
else
|
||||
if(lambda<1)
|
||||
if lambda<1
|
||||
lambda=max(lambda*2, 1);
|
||||
end;
|
||||
end;
|
||||
end;
|
||||
if(is_dynamic)
|
||||
end
|
||||
end
|
||||
end
|
||||
if is_dynamic
|
||||
ya = y(it_,y_index_eq)';
|
||||
else
|
||||
ya = y(y_index_eq);
|
||||
end;
|
||||
end
|
||||
ya_save=ya;
|
||||
g1a=g1;
|
||||
if(~is_dynamic && options.solve_algo == 0)
|
||||
if (verbose == 1)
|
||||
disp('steady: fsolve');
|
||||
if ~is_dynamic && options.solve_algo==0
|
||||
if verbose
|
||||
disp('steady: fsolve')
|
||||
end
|
||||
if ~isoctave
|
||||
if ~user_has_matlab_license('optimization_toolbox')
|
||||
|
|
@ -248,31 +221,30 @@ for it_=start:incr:finish
|
|||
else
|
||||
yn = y(y_index_eq);
|
||||
exitval = 3;
|
||||
end;
|
||||
end
|
||||
end
|
||||
|
||||
y(y_index_eq) = yn;
|
||||
if exitval > 0
|
||||
info = 0;
|
||||
else
|
||||
info = -Block_Num*10;
|
||||
end
|
||||
elseif((~is_dynamic && options.solve_algo==2) || (is_dynamic && stack_solve_algo==4))
|
||||
if (verbose == 1 && ~is_dynamic)
|
||||
disp('steady: LU + lnsrch1');
|
||||
elseif (~is_dynamic && options.solve_algo==2) || (is_dynamic && stack_solve_algo==4)
|
||||
if verbose==1 && ~is_dynamic
|
||||
disp('steady: LU + lnsrch1')
|
||||
end
|
||||
lambda=1;
|
||||
stpmx = 100 ;
|
||||
if (is_dynamic)
|
||||
if is_dynamic
|
||||
stpmax = stpmx*max([sqrt(ya'*ya);size(y_index_eq,2)]);
|
||||
else
|
||||
stpmax = stpmx*max([sqrt(ya'*ya);size(y_index_eq,2)]);
|
||||
end;
|
||||
end
|
||||
nn=1:size(y_index_eq,2);
|
||||
g = (r'*g1)';
|
||||
f = 0.5*r'*r;
|
||||
p = -g1\r ;
|
||||
if (is_dynamic)
|
||||
if is_dynamic
|
||||
[ya,f,r,check]=lnsrch1(y(it_,:)',f,g,p,stpmax, ...
|
||||
'lnsrch1_wrapper_one_boundary',nn, ...
|
||||
y_index_eq, y_index_eq, fname, y, x, params, steady_state, it_);
|
||||
|
|
@ -281,126 +253,133 @@ for it_=start:incr:finish
|
|||
[ya,f,r,check]=lnsrch1(y,f,g,p,stpmax,fname,nn,y_index_eq,x, ...
|
||||
params, steady_state,0);
|
||||
dx = ya - y(y_index_eq);
|
||||
end;
|
||||
|
||||
if(is_dynamic)
|
||||
end
|
||||
if is_dynamic
|
||||
y(it_,:) = ya';
|
||||
else
|
||||
y = ya';
|
||||
end;
|
||||
elseif(~is_dynamic && options.solve_algo==3)
|
||||
if (verbose == 1)
|
||||
disp('steady: csolve');
|
||||
end
|
||||
elseif ~is_dynamic && options.solve_algo==3
|
||||
if verbose==1
|
||||
disp('steady: csolve')
|
||||
end
|
||||
[yn,info] = csolve(@local_fname, y(y_index_eq),@ ...
|
||||
local_fname,1e-6,500, x, params, steady_state, y, y_index_eq, fname, 1);
|
||||
dx = ya - yn;
|
||||
y(y_index_eq) = yn;
|
||||
elseif((stack_solve_algo==1 && is_dynamic) || (stack_solve_algo==0 && is_dynamic) || (~is_dynamic && (options.solve_algo==1 || options.solve_algo==6))),
|
||||
if (verbose == 1 && ~is_dynamic)
|
||||
disp('steady: Sparse LU ');
|
||||
elseif (stack_solve_algo==1 && is_dynamic) || (stack_solve_algo==0 && is_dynamic) || (~is_dynamic && (options.solve_algo==1 || options.solve_algo==6))
|
||||
if verbose==1 && ~is_dynamic
|
||||
disp('steady: Sparse LU ')
|
||||
end
|
||||
dx = g1\r;
|
||||
ya = ya - lambda*dx;
|
||||
if(is_dynamic)
|
||||
if is_dynamic
|
||||
y(it_,y_index_eq) = ya';
|
||||
else
|
||||
y(y_index_eq) = ya;
|
||||
end;
|
||||
elseif((stack_solve_algo==2 && is_dynamic) || (options.solve_algo==7 && ~is_dynamic)),
|
||||
end
|
||||
elseif (stack_solve_algo==2 && is_dynamic) || (options.solve_algo==7 && ~is_dynamic)
|
||||
flag1=1;
|
||||
if isoctave
|
||||
error('SOLVE_ONE_BOUNDARY: you can''t use solve_algo=7 since GMRES is not implemented in Octave')
|
||||
end
|
||||
if (verbose == 1 && ~is_dynamic)
|
||||
disp('steady: GMRES ');
|
||||
if verbose == 1 && ~is_dynamic
|
||||
disp('steady: GMRES ')
|
||||
end
|
||||
while(flag1>0)
|
||||
while flag1>0
|
||||
[L1, U1]=ilu(g1,ilu_setup);
|
||||
[dx,flag1] = gmres(g1,-r,Blck_size,1e-6,Blck_size,L1,U1);
|
||||
if (flag1>0 || reduced)
|
||||
if(flag1==1)
|
||||
disp(['Error in simul: No convergence inside GMRES after ' num2str(iter,'%6d') ' iterations, in block' num2str(Block_Num,'%3d')]);
|
||||
elseif(flag1==2)
|
||||
disp(['Error in simul: Preconditioner is ill-conditioned, in block' num2str(Block_Num,'%3d')]);
|
||||
elseif(flag1==3)
|
||||
disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block' num2str(Block_Num,'%3d')]);
|
||||
end;
|
||||
if flag1>0 || reduced
|
||||
if verbose
|
||||
if flag1==1
|
||||
disp(['Error in simul: No convergence inside GMRES after ' num2str(iter,'%6d') ' iterations, in block' num2str(Block_Num,'%3d')])
|
||||
elseif(flag1==2)
|
||||
disp(['Error in simul: Preconditioner is ill-conditioned, in block' num2str(Block_Num,'%3d')])
|
||||
elseif(flag1==3)
|
||||
disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block' num2str(Block_Num,'%3d')])
|
||||
end
|
||||
end
|
||||
ilu_setup.droptol = ilu_setup.droptol/10;
|
||||
reduced = 0;
|
||||
else
|
||||
ya = ya + lambda*dx;
|
||||
if(is_dynamic)
|
||||
if is_dynamic
|
||||
y(it_,y_index_eq) = ya';
|
||||
else
|
||||
y(y_index_eq) = ya';
|
||||
end;
|
||||
end;
|
||||
end;
|
||||
elseif((stack_solve_algo==3 && is_dynamic) || (options.solve_algo==8 && ~is_dynamic)),
|
||||
flag1=1;
|
||||
if (verbose == 1 && ~is_dynamic)
|
||||
disp('steady: BiCGStab');
|
||||
end
|
||||
end
|
||||
end
|
||||
while(flag1>0)
|
||||
elseif (stack_solve_algo==3 && is_dynamic) || (options.solve_algo==8 && ~is_dynamic)
|
||||
flag1=1;
|
||||
if verbose == 1 && ~is_dynamic
|
||||
disp('steady: BiCGStab')
|
||||
end
|
||||
while flag1>0
|
||||
[L1, U1]=ilu(g1,ilu_setup);
|
||||
phat = ya - U1 \ (L1 \ r);
|
||||
if(is_dynamic)
|
||||
if is_dynamic
|
||||
y(it_,y_index_eq) = phat;
|
||||
else
|
||||
y(y_index_eq) = phat;
|
||||
end;
|
||||
if(is_dynamic)
|
||||
end
|
||||
if is_dynamic
|
||||
[r, y, g1, g2, g3] = feval(fname, y, x, params, ...
|
||||
steady_state, it_, 0);
|
||||
else
|
||||
[r, y, g1] = feval(fname, y, x, params);
|
||||
end;
|
||||
if max(abs(r)) >= options.solve_tolf
|
||||
end
|
||||
if max(abs(r))>=options.solve_tolf
|
||||
[dx,flag1] = bicgstab(g1,-r,1e-7,Blck_size,L1,U1);
|
||||
else
|
||||
flag1 = 0;
|
||||
dx = phat - ya;
|
||||
end;
|
||||
if (flag1>0 || reduced)
|
||||
if(flag1==1)
|
||||
disp(['Error in simul: No convergence inside BICGSTAB after ' num2str(iter,'%6d') ' iterations, in block' num2str(Block_Num,'%3d')]);
|
||||
elseif(flag1==2)
|
||||
disp(['Error in simul: Preconditioner is ill-conditioned, in block' num2str(Block_Num,'%3d')]);
|
||||
elseif(flag1==3)
|
||||
disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block' num2str(Block_Num,'%3d')]);
|
||||
end;
|
||||
end
|
||||
if flag1>0 || reduced
|
||||
if verbose
|
||||
if(flag1==1)
|
||||
disp(['Error in simul: No convergence inside BICGSTAB after ' num2str(iter,'%6d') ' iterations, in block' num2str(Block_Num,'%3d')])
|
||||
elseif(flag1==2)
|
||||
disp(['Error in simul: Preconditioner is ill-conditioned, in block' num2str(Block_Num,'%3d')])
|
||||
elseif(flag1==3)
|
||||
disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block' num2str(Block_Num,'%3d')])
|
||||
end
|
||||
end
|
||||
ilu_setup.droptol = ilu_setup.droptol/10;
|
||||
reduced = 0;
|
||||
else
|
||||
ya = ya + lambda*dx;
|
||||
if(is_dynamic)
|
||||
if is_dynamic
|
||||
y(it_,y_index_eq) = ya';
|
||||
else
|
||||
y(y_index_eq) = ya';
|
||||
end;
|
||||
end;
|
||||
end;
|
||||
end
|
||||
end
|
||||
end
|
||||
else
|
||||
disp('unknown option : ');
|
||||
if(is_dynamic)
|
||||
disp(['options_.stack_solve_algo = ' num2str(stack_solve_algo) ' not implemented']);
|
||||
else
|
||||
disp(['options_.solve_algo = ' num2str(options.solve_algo) ' not implemented']);
|
||||
end;
|
||||
if verbose
|
||||
disp('unknown option : ')
|
||||
if is_dynamic
|
||||
disp(['options_.stack_solve_algo = ' num2str(stack_solve_algo) ' not implemented'])
|
||||
else
|
||||
disp(['options_.solve_algo = ' num2str(options.solve_algo) ' not implemented'])
|
||||
end
|
||||
end
|
||||
info = -Block_Num*10;
|
||||
return;
|
||||
end;
|
||||
return
|
||||
end
|
||||
iter=iter+1;
|
||||
max_resa = max_res;
|
||||
end
|
||||
end
|
||||
if cvg==0
|
||||
if(cutoff == 0)
|
||||
fprintf('Error in simul: Convergence not achieved in block %d, at time %d, after %d iterations.\n Increase "options_.simul.maxit\".\n',Block_Num, it_,iter);
|
||||
else
|
||||
fprintf('Error in simul: Convergence not achieved in block %d, at time %d, after %d iterations.\n Increase "options_.simul.maxit" or set "cutoff=0" in model options.\n',Block_Num, it_,iter);
|
||||
end;
|
||||
if verbose
|
||||
if cutoff == 0
|
||||
fprintf('Error in simul: Convergence not achieved in block %d, at time %d, after %d iterations.\n Increase "options_.simul.maxit\".\n',Block_Num, it_,iter);
|
||||
else
|
||||
fprintf('Error in simul: Convergence not achieved in block %d, at time %d, after %d iterations.\n Increase "options_.simul.maxit" or set "cutoff=0" in model options.\n',Block_Num, it_,iter);
|
||||
end
|
||||
end
|
||||
if(is_dynamic)
|
||||
oo_.deterministic_simulation.status = 0;
|
||||
oo_.deterministic_simulation.error = max_res;
|
||||
|
|
@ -410,10 +389,11 @@ for it_=start:incr:finish
|
|||
oo_.deterministic_simulation.block(Block_Num).iterations = iter;
|
||||
end;
|
||||
info = -Block_Num*10;
|
||||
return;
|
||||
return
|
||||
end
|
||||
end
|
||||
if(is_dynamic)
|
||||
|
||||
if is_dynamic
|
||||
info = 1;
|
||||
oo_.deterministic_simulation.status = 1;
|
||||
oo_.deterministic_simulation.error = max_res;
|
||||
|
|
@ -423,12 +403,11 @@ if(is_dynamic)
|
|||
oo_.deterministic_simulation.block(Block_Num).iterations = iter;
|
||||
else
|
||||
info = 0;
|
||||
end;
|
||||
return;
|
||||
end
|
||||
|
||||
function [err, G]=local_fname(yl, x, params, steady_state, y, y_index_eq, fname, is_csolve)
|
||||
y(y_index_eq) = yl;
|
||||
[err, y, G] = feval(fname, y, x, params, steady_state, 0);
|
||||
if(is_csolve)
|
||||
G = full(G);
|
||||
end;
|
||||
end
|
||||
|
|
@ -63,6 +63,8 @@ function [y, oo]= solve_two_boundaries(fname, y, x, params, steady_state, y_inde
|
|||
% You should have received a copy of the GNU General Public License
|
||||
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
verbose = options.verbosity;
|
||||
|
||||
cvg=0;
|
||||
iter=0;
|
||||
Per_u_=0;
|
||||
|
|
@ -80,7 +82,7 @@ max_resa=1e100;
|
|||
Jacobian_Size=Blck_size*(y_kmin+y_kmax_l +periods);
|
||||
g1=spalloc( Blck_size*periods, Jacobian_Size, nze*periods);
|
||||
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,options.periods);
|
||||
preconditioner = 2;
|
||||
g1a=g1(:, y_kmin*Blck_size+1:(periods+y_kmin)*Blck_size);
|
||||
|
|
@ -88,83 +90,86 @@ while ~(cvg==1 || iter>maxit_),
|
|||
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;
|
||||
[max_res, max_indx]=max(max(abs(r')));
|
||||
if(~isreal(r))
|
||||
if ~isreal(r)
|
||||
max_res = (-max_res^2)^0.5;
|
||||
end;
|
||||
% if(~isreal(r))
|
||||
% max_res=(-(max(max(abs(r))))^2)^0.5;
|
||||
% else
|
||||
% max_res=max(max(abs(r)));
|
||||
% end;
|
||||
if(~isreal(max_res) || isnan(max_res))
|
||||
if ~isreal(max_res) || isnan(max_res)
|
||||
cvg = 0;
|
||||
elseif(is_linear && iter>0)
|
||||
cvg = 1;
|
||||
else
|
||||
cvg=(max_res<solve_tolf);
|
||||
end;
|
||||
if(~cvg)
|
||||
if(iter>0)
|
||||
if(~isreal(max_res) || isnan(max_res) || (max_resa<max_res && iter>1))
|
||||
if(~isreal(max_res))
|
||||
end
|
||||
if ~cvg
|
||||
if iter>0
|
||||
if ~isreal(max_res) || isnan(max_res) || (max_resa<max_res && iter>1)
|
||||
if verbose && ~isreal(max_res)
|
||||
disp(['Variable ' M.endo_names(max_indx,:) ' (' int2str(max_indx) ') returns an undefined value']);
|
||||
end;
|
||||
if(isnan(max_res))
|
||||
end
|
||||
if isnan(max_res)
|
||||
detJ=det(g1aa);
|
||||
if(abs(detJ)<1e-7)
|
||||
if abs(detJ)<1e-7
|
||||
max_factor=max(max(abs(g1aa)));
|
||||
ze_elem=sum(diag(g1aa)<cutoff);
|
||||
disp([num2str(full(ze_elem),'%d') ' elements on the Jacobian diagonal are below the cutoff (' num2str(cutoff,'%f') ')']);
|
||||
if(correcting_factor<max_factor)
|
||||
if verbose
|
||||
disp([num2str(full(ze_elem),'%d') ' elements on the Jacobian diagonal are below the cutoff (' num2str(cutoff,'%f') ')']);
|
||||
end
|
||||
if correcting_factor<max_factor
|
||||
correcting_factor=correcting_factor*4;
|
||||
disp(['The Jacobain matrix is singular, det(Jacobian)=' num2str(detJ,'%f') '.']);
|
||||
disp([' trying to correct the Jacobian matrix:']);
|
||||
disp([' correcting_factor=' num2str(correcting_factor,'%f') ' max(Jacobian)=' num2str(full(max_factor),'%f')]);
|
||||
if verbose
|
||||
disp(['The Jacobain matrix is singular, det(Jacobian)=' num2str(detJ,'%f') '.']);
|
||||
disp([' trying to correct the Jacobian matrix:']);
|
||||
disp([' correcting_factor=' num2str(correcting_factor,'%f') ' max(Jacobian)=' num2str(full(max_factor),'%f')]);
|
||||
end
|
||||
dx = (g1aa+correcting_factor*speye(periods*Blck_size))\ba- ya;
|
||||
y(1+y_kmin:periods+y_kmin,y_index)=reshape((ya_save+lambda*dx)',length(y_index),periods)';
|
||||
continue;
|
||||
else
|
||||
disp('The singularity of the jacobian matrix could not be corrected');
|
||||
return;
|
||||
end;
|
||||
end;
|
||||
elseif(lambda>1e-8)
|
||||
return
|
||||
end
|
||||
end
|
||||
elseif lambda>1e-8
|
||||
lambda=lambda/2;
|
||||
reduced = 1;
|
||||
disp(['reducing the path length: lambda=' num2str(lambda,'%f')]);
|
||||
if verbose
|
||||
disp(['reducing the path length: lambda=' num2str(lambda,'%f')]);
|
||||
end
|
||||
y(1+y_kmin:periods+y_kmin,y_index)=reshape((ya_save+lambda*dx)',length(y_index),periods)';
|
||||
continue;
|
||||
continue
|
||||
else
|
||||
if(cutoff == 0)
|
||||
fprintf('Error in simul: Convergence not achieved in block %d, after %d iterations.\n Increase "options_.simul.maxit".\n',Block_Num, iter);
|
||||
else
|
||||
fprintf('Error in simul: Convergence not achieved in block %d, after %d iterations.\n Increase "options_.simul.maxit" or set "cutoff=0" in model options.\n',Block_Num, iter);
|
||||
end;
|
||||
if verbose
|
||||
if cutoff==0
|
||||
fprintf('Error in simul: Convergence not achieved in block %d, after %d iterations.\n Increase "options_.simul.maxit".\n',Block_Num, iter);
|
||||
else
|
||||
fprintf('Error in simul: Convergence not achieved in block %d, after %d iterations.\n Increase "options_.simul.maxit" or set "cutoff=0" in model options.\n',Block_Num, iter);
|
||||
end
|
||||
end
|
||||
oo.deterministic_simulation.status = 0;
|
||||
oo.deterministic_simulation.error = max_res;
|
||||
oo.deterministic_simulation.iterations = iter;
|
||||
oo.deterministic_simulation.block(Block_Num).status = 0;% Convergency failed.
|
||||
oo.deterministic_simulation.block(Block_Num).error = max_res;
|
||||
oo.deterministic_simulation.block(Block_Num).iterations = iter;
|
||||
return;
|
||||
end;
|
||||
return
|
||||
end
|
||||
else
|
||||
if(lambda<1)
|
||||
if lambda<1
|
||||
lambda=max(lambda*2, 1);
|
||||
end;
|
||||
end;
|
||||
end;
|
||||
end
|
||||
end
|
||||
end
|
||||
ya = reshape(y(y_kmin+1:y_kmin+periods,y_index)',1,periods*Blck_size)';
|
||||
ya_save=ya;
|
||||
g1aa=g1a;
|
||||
ba=b;
|
||||
max_resa=max_res;
|
||||
if(stack_solve_algo==0),
|
||||
if stack_solve_algo==0
|
||||
dx = g1a\b- ya;
|
||||
ya = ya + lambda*dx;
|
||||
y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)';
|
||||
elseif(stack_solve_algo==1),
|
||||
for t=1:periods;
|
||||
elseif stack_solve_algo==1
|
||||
for t=1:periods
|
||||
first_elem = (t-1)*Blck_size+1;
|
||||
last_elem = t*Blck_size;
|
||||
next_elem = (t+1)*Blck_size;
|
||||
|
|
@ -177,19 +182,19 @@ while ~(cvg==1 || iter>maxit_),
|
|||
g1a(Elem, Elem_1) = S1;
|
||||
b(Elem) = B1_inv * b(Elem);
|
||||
g1a(Elem, Elem) = ones(Blck_size, Blck_size);
|
||||
if (t < periods)
|
||||
if t<periods
|
||||
g1a(Elem_1, Elem_1) = g1a(Elem_1, Elem_1) - g1a(Elem_1, Elem) * S1;
|
||||
b(Elem_1) = b(Elem_1) - g1a(Elem_1, Elem) * b(Elem);
|
||||
g1a(Elem_1, Elem) = zeros(Blck_size, Blck_size);
|
||||
end;
|
||||
end;
|
||||
end
|
||||
end
|
||||
za = b(Elem);
|
||||
zaa = za;
|
||||
y_Elem = (periods - 1) * Blck_size + 1:(periods) * Blck_size;
|
||||
dx = ya;
|
||||
dx(y_Elem) = za - ya(y_Elem);
|
||||
ya(y_Elem) = ya(y_Elem) + lambda*dx(y_Elem);
|
||||
for t=periods-1:-1:1;
|
||||
for t=periods-1:-1:1
|
||||
first_elem = (t-1)*Blck_size+1;
|
||||
last_elem = t*Blck_size;
|
||||
next_elem = (t+1)*Blck_size;
|
||||
|
|
@ -201,29 +206,29 @@ while ~(cvg==1 || iter>maxit_),
|
|||
dx(y_Elem) = za - ya(y_Elem);
|
||||
ya(y_Elem) = ya(y_Elem) + lambda*dx(y_Elem);
|
||||
y(y_kmin + t, y_index) = ya(y_Elem);
|
||||
end;
|
||||
elseif(stack_solve_algo==2),
|
||||
end
|
||||
elseif stack_solve_algo==2
|
||||
flag1=1;
|
||||
while(flag1>0)
|
||||
if preconditioner == 2
|
||||
while flag1>0
|
||||
if preconditioner==2
|
||||
[L1, U1]=ilu(g1a,ilu_setup);
|
||||
elseif preconditioner == 3
|
||||
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
|
||||
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);
|
||||
|
|
@ -231,31 +236,33 @@ while ~(cvg==1 || iter>maxit_),
|
|||
U1 = kron(eye(ceil(periods/3)),U);
|
||||
L1 = L1(1:periods * Size, 1:periods * Size);
|
||||
U1 = U1(1:periods * Size, 1:periods * Size);
|
||||
end;
|
||||
end
|
||||
[za,flag1] = gmres(g1a,b,Blck_size,1e-6,Blck_size*periods,L1,U1);
|
||||
if (flag1>0 || reduced)
|
||||
if(flag1==1)
|
||||
disp(['Error in simul: No convergence inside GMRES after ' num2str(periods*10,'%6d') ' iterations, in block ' num2str(Blck_size,'%3d')]);
|
||||
elseif(flag1==2)
|
||||
disp(['Error in simul: Preconditioner is ill-conditioned, in block ' num2str(Blck_size,'%3d')]);
|
||||
elseif(flag1==3)
|
||||
disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block ' num2str(Blck_size,'%3d')]);
|
||||
end;
|
||||
if verbose
|
||||
if flag1==1
|
||||
disp(['Error in simul: No convergence inside GMRES after ' num2str(periods*10,'%6d') ' iterations, in block ' num2str(Blck_size,'%3d')]);
|
||||
elseif flag1==2
|
||||
disp(['Error in simul: Preconditioner is ill-conditioned, in block ' num2str(Blck_size,'%3d')]);
|
||||
elseif flag1==3
|
||||
disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block ' num2str(Blck_size,'%3d')]);
|
||||
end
|
||||
end
|
||||
ilu_setup.droptol = ilu_setup.droptol/10;
|
||||
reduced = 0;
|
||||
else
|
||||
dx = za - ya;
|
||||
ya = ya + lambda*dx;
|
||||
y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)';
|
||||
end;
|
||||
end;
|
||||
elseif(stack_solve_algo==3),
|
||||
end
|
||||
end
|
||||
elseif stack_solve_algo==3
|
||||
flag1=1;
|
||||
while(flag1>0)
|
||||
if preconditioner == 2
|
||||
while flag1>0
|
||||
if preconditioner==2
|
||||
[L1, U1]=ilu(g1a,ilu_setup);
|
||||
[za,flag1] = bicgstab(g1a,b,1e-7,Blck_size*periods,L1,U1);
|
||||
elseif (preconditioner == 3)
|
||||
elseif preconditioner==3
|
||||
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);
|
||||
|
|
@ -273,24 +280,26 @@ while ~(cvg==1 || iter>maxit_),
|
|||
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==1)
|
||||
disp(['Error in simul: No convergence inside BICGSTAB after ' num2str(periods*10,'%6d') ' iterations, in block ' num2str(Blck_size,'%3d')]);
|
||||
elseif(flag1==2)
|
||||
disp(['Error in simul: Preconditioner is ill-conditioned, in block ' num2str(Blck_size,'%3d')]);
|
||||
elseif(flag1==3)
|
||||
disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block ' num2str(Blck_size,'%3d')]);
|
||||
end;
|
||||
end
|
||||
if flag1>0 || reduced
|
||||
if verbose
|
||||
if flag1==1
|
||||
disp(['Error in simul: No convergence inside BICGSTAB after ' num2str(periods*10,'%6d') ' iterations, in block ' num2str(Blck_size,'%3d')]);
|
||||
elseif flag1==2
|
||||
disp(['Error in simul: Preconditioner is ill-conditioned, in block ' num2str(Blck_size,'%3d')]);
|
||||
elseif flag1==3
|
||||
disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block ' num2str(Blck_size,'%3d')]);
|
||||
end
|
||||
end
|
||||
ilu_setup.droptol = ilu_setup.droptol/10;
|
||||
reduced = 0;
|
||||
else
|
||||
dx = za - ya;
|
||||
ya = ya + lambda*dx;
|
||||
y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)';
|
||||
end;
|
||||
end;
|
||||
elseif(stack_solve_algo==4),
|
||||
end
|
||||
end
|
||||
elseif stack_solve_algo==4
|
||||
ra = reshape(r(:, y_kmin+1:periods+y_kmin),periods*Blck_size, 1);
|
||||
stpmx = 100 ;
|
||||
stpmax = stpmx*max([sqrt(ya'*ya);size(y_index,2)]);
|
||||
|
|
@ -304,22 +313,28 @@ while ~(cvg==1 || iter>maxit_),
|
|||
end
|
||||
end
|
||||
iter=iter+1;
|
||||
disp(['iteration: ' num2str(iter,'%d') ' error: ' num2str(max_res,'%e')]);
|
||||
end;
|
||||
if verbose
|
||||
disp(['iteration: ' num2str(iter,'%d') ' error: ' num2str(max_res,'%e')]);
|
||||
end
|
||||
end
|
||||
|
||||
if (iter>maxit_)
|
||||
disp(['No convergence after ' num2str(iter,'%4d') ' iterations in Block ' num2str(Block_Num,'%d')]);
|
||||
if verbose
|
||||
printline(41)
|
||||
%disp(['No convergence after ' num2str(iter,'%4d') ' iterations in Block ' num2str(Block_Num,'%d')])
|
||||
end
|
||||
oo.deterministic_simulation.status = 0;
|
||||
oo.deterministic_simulation.error = max_res;
|
||||
oo.deterministic_simulation.iterations = iter;
|
||||
oo.deterministic_simulation.block(Block_Num).status = 0;% Convergency failed.
|
||||
oo.deterministic_simulation.block(Block_Num).error = max_res;
|
||||
oo.deterministic_simulation.block(Block_Num).iterations = iter;
|
||||
return;
|
||||
return
|
||||
end
|
||||
|
||||
oo.deterministic_simulation.status = 1;
|
||||
oo.deterministic_simulation.error = max_res;
|
||||
oo.deterministic_simulation.iterations = iter;
|
||||
oo.deterministic_simulation.block(Block_Num).status = 1;% Convergency obtained.
|
||||
oo.deterministic_simulation.block(Block_Num).error = max_res;
|
||||
oo.deterministic_simulation.block(Block_Num).iterations = iter;
|
||||
return;
|
||||
oo.deterministic_simulation.block(Block_Num).iterations = iter;
|
||||
|
|
@ -1912,9 +1912,11 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
|
|||
<< " else" << endl
|
||||
<< " mthd='UNKNOWN';" << endl
|
||||
<< " end;" << endl
|
||||
<< " disp (['-----------------------------------------------------']) ;" << endl
|
||||
<< " disp (['MODEL SIMULATION: (method=' mthd ')']) ;" << endl
|
||||
<< " fprintf('\\n') ;" << endl
|
||||
<< " if options_.verbosity" << endl
|
||||
<< " printline(41)" << endl
|
||||
<< " disp(sprintf('MODEL SIMULATION (method=%s):',mthd))" << endl
|
||||
<< " skipline()" << endl
|
||||
<< " end" << endl
|
||||
<< " periods=options_.periods;" << endl
|
||||
<< " maxit_=options_.simul.maxit;" << endl
|
||||
<< " solve_tolf=options_.solve_tolf;" << endl
|
||||
|
|
@ -1951,7 +1953,7 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
|
|||
mDynamicModelFile << " g1=[];g2=[];g3=[];\n";
|
||||
mDynamicModelFile << " y=" << dynamic_basename << "_" << block + 1 << "(y, x, params, steady_state, 0, y_kmin, periods);\n";
|
||||
mDynamicModelFile << " tmp = y(:,M_.block_structure.block(" << block + 1 << ").variable);\n";
|
||||
mDynamicModelFile << " if any(isnan(tmp) | isinf(tmp))\n";
|
||||
mDynamicModelFile << " if any(isnan(tmp) | isinf(tmp))\n";
|
||||
mDynamicModelFile << " disp(['Inf or Nan value during the evaluation of block " << block <<"']);\n";
|
||||
mDynamicModelFile << " return;\n";
|
||||
mDynamicModelFile << " end;\n";
|
||||
|
|
|
|||
Loading…
Reference in New Issue