diff --git a/matlab/perfect-foresight-models/private/simulation_core.m b/matlab/perfect-foresight-models/private/simulation_core.m
index 79f91604864097cd58002a4bea6ebdc2986af643..d09b51b585093c75e8e17489153772b08bdff50c 100644
--- a/matlab/perfect-foresight-models/private/simulation_core.m
+++ b/matlab/perfect-foresight-models/private/simulation_core.m
@@ -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
diff --git a/matlab/solve_one_boundary.m b/matlab/solve_one_boundary.m
index 304b81892851dbbbcbbd1fc3d0733b2e076762ce..8c3d8cca4947d5162be39308c85b7f18bbfd07f9 100644
--- a/matlab/solve_one_boundary.m
+++ b/matlab/solve_one_boundary.m
@@ -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)),
+ 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');
+ if verbose == 1 && ~is_dynamic
+ disp('steady: BiCGStab')
end
- while(flag1>0)
+ 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
\ No newline at end of file
diff --git a/matlab/solve_two_boundaries.m b/matlab/solve_two_boundaries.m
index cddc8f2da14f76021559165b1143c65e3b00925a..cf00fdacb359b74d33a02e08b0edd27249e249bf 100644
--- a/matlab/solve_two_boundaries.m
+++ b/matlab/solve_two_boundaries.m
@@ -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;
\ No newline at end of file
diff --git a/preprocessor/DynamicModel.cc b/preprocessor/DynamicModel.cc
index f4efef49b893894bb0db0dd287a4d5029a293a62..645fa6469322bf1d413d898da7bc49100408999f 100644
--- a/preprocessor/DynamicModel.cc
+++ b/preprocessor/DynamicModel.cc
@@ -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";