From 3ca9485eddb8026fdffac2c81d2324ca5831894b Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Wed, 4 Apr 2012 10:35:18 +0200 Subject: [PATCH 1/7] bug fix when calling gsa set default order=1 set analytic derivation =1 --- matlab/dynare_identification.m | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/matlab/dynare_identification.m b/matlab/dynare_identification.m index b60ffec94..b9384bda9 100644 --- a/matlab/dynare_identification.m +++ b/matlab/dynare_identification.m @@ -61,7 +61,7 @@ options_ident = set_default_option(options_ident,'advanced',0); options_ident = set_default_option(options_ident,'normalize_jacobians',1); options_ident = set_default_option(options_ident,'lik_init',1); if options_ident.gsa_sample_file, - GSAFolder = checkpath('gsa'); + GSAFolder = checkpath('gsa',M_.dname); if options_ident.gsa_sample_file==1, load([GSAFolder,filesep,fname_,'_prior'],'lpmat','lpmat0','istable'); elseif options_ident.gsa_sample_file==2, @@ -105,11 +105,13 @@ nlags = options_ident.ar; periods = options_ident.periods; replic = options_ident.replic; useautocorr = options_ident.useautocorr; +options_.order=1; options_.ar=nlags; options_.prior_mc = options_ident.prior_mc; options_.options_ident = options_ident; options_.Schur_vec_tol = 1.e-8; options_.nomoments=0; +options_.analytic_derivation=1; options_ = set_default_option(options_,'datafile',[]); options_.mode_compute = 0; From fe6aa88c2745f2f1374e791489766295bc20c7d7 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Wed, 4 Apr 2012 10:36:32 +0200 Subject: [PATCH 2/7] 1) check if AHess is positive semi-definite; 2) fix memory issues related to chh 3) fix issue for the rare case of empty indok (octave) --- matlab/identification_analysis.m | 23 +++++++++++++++++------ 1 file changed, 17 insertions(+), 6 deletions(-) diff --git a/matlab/identification_analysis.m b/matlab/identification_analysis.m index 2347418f8..695b2d44b 100644 --- a/matlab/identification_analysis.m +++ b/matlab/identification_analysis.m @@ -138,6 +138,10 @@ if info(1)==0, [fval,cost_flag,ys,trend_coeff,info,M_,options_,bayestopt_,oo_,DLIK,AHess] = dsge_likelihood(params',data_info,options_,M_,estim_params_,bayestopt_,oo_,derivatives_info); % fval = DsgeLikelihood(xparam1,data_info,options_,M_,estim_params_,bayestopt_,oo_); AHess=-AHess; + if min(eig(AHess))<0, + error('Analytic Hessian is not positive semi-definite!') + end +% chol(AHess); ide_hess.AHess= AHess; deltaM = sqrt(diag(AHess)); iflag=any((deltaM.*deltaM)==0); @@ -153,7 +157,9 @@ if info(1)==0, cmm = NaN(size(siJ,1),size(siJ,1)); ind1=find(ide_hess.ind0); cmm = siJ(:,ind1)*((AHess(ind1,ind1))\siJ(:,ind1)'); - chh = siH(:,ind1)*((AHess(ind1,ind1))\siH(:,ind1)'); + temp1=((AHess(ind1,ind1))\siH(:,ind1)'); + diag_chh=sum(siH(:,ind1)'.*temp1)'; +% chh = siH(:,ind1)*((AHess(ind1,ind1))\siH(:,ind1)'); ind1=ind1(ind1>offset); clre = siLRE(:,ind1-offset)*((AHess(ind1,ind1))\siLRE(:,ind1-offset)'); rhoM=sqrt(1./diag(inv(tildaM(indok,indok)))); @@ -190,12 +196,16 @@ if info(1)==0, % rhoM=sqrt(1-1./diag(inv(tildaM))); % rhoM=(1-1./diag(inv(tildaM))); ind1=find(ide_hess.ind0); - chh = siH(:,ind1)*((MIM(ind1,ind1))\siH(:,ind1)'); + temp1=((MIM(ind1,ind1))\siH(:,ind1)'); + diag_chh=sum(siH(:,ind1)'.*temp1)'; +% chh = siH(:,ind1)*((MIM(ind1,ind1))\siH(:,ind1)'); ind1=ind1(ind1>offset); clre = siLRE(:,ind1-offset)*((MIM(ind1,ind1))\siLRE(:,ind1-offset)'); - rhoM(indok)=sqrt(1./diag(inv(tildaM(indok,indok)))); - normaliz(indok) = (sqrt(diag(inv(tildaM(indok,indok))))./deltaM(indok))'; %sqrt(diag(inv(MIM(indok,indok))))'; -% deltaM = deltaM.*abs(params'); + if ~isempty(indok), + rhoM(indok)=sqrt(1./diag(inv(tildaM(indok,indok)))); + normaliz(indok) = (sqrt(diag(inv(tildaM(indok,indok))))./deltaM(indok))'; %sqrt(diag(inv(MIM(indok,indok))))'; + end + % deltaM = deltaM.*abs(params') flag_score=0; end ide_strength_J(indok) = (1./(normaliz(indok)'./abs(params(indok)'))); @@ -212,7 +222,8 @@ if info(1)==0, % isok = find((abs(TAU)>=1.e-8)); % quant(isok,:) = siH(isok,:)./repmat(TAU(isok,1),1,nparam); % quant(inok,:) = siH(inok,:)./repmat(mean(abs(TAU)),length(inok),nparam); - quant = siH./repmat(sqrt(diag(chh)),1,nparam); +% quant = siH./repmat(sqrt(diag(chh)),1,nparam); + quant = siH./repmat(sqrt(diag_chh),1,nparam); siHnorm = vnorm(quant).*normaliz1; % siHnorm = vnorm(siH./repmat(TAU,1,nparam)).*normaliz; quant=[]; From a49e5504480a689da7ddf4203ec72d9b68ae2a7e Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Wed, 4 Apr 2012 10:36:54 +0200 Subject: [PATCH 3/7] typo fixed --- matlab/simulated_moment_uncertainty.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/matlab/simulated_moment_uncertainty.m b/matlab/simulated_moment_uncertainty.m index f65e80f70..f40be9497 100644 --- a/matlab/simulated_moment_uncertainty.m +++ b/matlab/simulated_moment_uncertainty.m @@ -22,7 +22,7 @@ global options_ oo_ mm=zeros(length(indx),replic); -disp('Evaluting simulated moment uncertainty ... please wait') +disp('Evaluating simulated moment uncertainty ... please wait') disp(['Doing ',int2str(replic),' replicas of length ',int2str(periods),' periods.']) noprint0 = options_.noprint; for j=1:replic; From 76093094a65e45c01153368301954bad5d9ab124 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Wed, 4 Apr 2012 10:39:19 +0200 Subject: [PATCH 4/7] Force precision in sylvester equation solution, particularly useful for different behaviour of the function schur under octave --- matlab/gensylv/sylvester3amr.m | 38 ++++++++++++++++++++++++++++++++++ matlab/getH.m | 8 +++++++ 2 files changed, 46 insertions(+) create mode 100644 matlab/gensylv/sylvester3amr.m diff --git a/matlab/gensylv/sylvester3amr.m b/matlab/gensylv/sylvester3amr.m new file mode 100644 index 000000000..4eb696296 --- /dev/null +++ b/matlab/gensylv/sylvester3amr.m @@ -0,0 +1,38 @@ +function [x0, flag]=sylvester3amr(x0,a,b,c,dd) +% solves iteratively ax+bxc=d + +% Copyright (C) 2012 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 . + +a_1 = inv(a); +b = a_1*b; +flag=0; +for j=1:size(dd,3), + d = a_1*dd(:,:,j); + e = 1; + iter = 1; + while e > 1e-8 && iter < 500 + x = d-b*x0(:,:,j)*c; + e = max(max(abs(x-x0(:,:,j)))); + x0(:,:,j) = x; + iter = iter + 1; + end + if iter == 500 + sprintf('sylvester3amr : Only accuracy of %10.8f is achieved after 500 iterations',e); + flag=1; + end +end \ No newline at end of file diff --git a/matlab/getH.m b/matlab/getH.m index a457969bb..5d9316ed4 100644 --- a/matlab/getH.m +++ b/matlab/getH.m @@ -356,6 +356,10 @@ else % generalized sylvester equation d(:,:,j) = Dg2(:,:,j)-elem(:,:,j)*A; end xx=sylvester3mr(a,b,c,d); + flag=1; + while flag, + [xx, flag]=sylvester3amr(xx,a,b,c,d); + end H=zeros(m*m+m*(m+1)/2,param_nbr+length(indexo)); if nargout>1, dOm = zeros(m,m,param_nbr+length(indexo)); @@ -436,6 +440,10 @@ if nargout > 5, end end xx2=sylvester3mr(a,b,c,d); + flag=1; + while flag, + [xx2, flag]=sylvester3amr(xx2,a,b,c,d); + end jcount = 0; for j=1:param_nbr, for i=j:param_nbr, From f455557376dad9f39eb75a5c159dd230a8db11d1 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Wed, 4 Apr 2012 16:38:22 +0200 Subject: [PATCH 5/7] compact all sylvester routines, with generic dimension of d argument --- matlab/gensylv/sylvester3.m | 61 +++++++++++++------- matlab/gensylv/sylvester3a.m | 30 +++++----- matlab/gensylv/sylvester3amr.m | 38 ------------- matlab/gensylv/sylvester3mr.m | 100 --------------------------------- 4 files changed, 57 insertions(+), 172 deletions(-) delete mode 100644 matlab/gensylv/sylvester3amr.m delete mode 100644 matlab/gensylv/sylvester3mr.m diff --git a/matlab/gensylv/sylvester3.m b/matlab/gensylv/sylvester3.m index 389ff7fb7..b2b18f594 100644 --- a/matlab/gensylv/sylvester3.m +++ b/matlab/gensylv/sylvester3.m @@ -1,5 +1,5 @@ -function x=sylvester3(a,b,c,d) -% solves a*x+b*x*c=d +function x=sylvester3mr(a,b,c,d) +% solves a*x+b*x*c=d where d is [n x m x p] % Copyright (C) 2005-2009 Dynare Team % @@ -20,53 +20,72 @@ function x=sylvester3(a,b,c,d) n = size(a,1); m = size(c,1); +p = size(d,3); if n == 1 - x=d./(a*ones(1,m)+b*c); + for j=1:p, + x(:,:,j)=d(:,:,j)./(a*ones(1,m)+b*c); + end return end if m == 1 - x = (a+c*b)\d; + for j=1:p, + x(:,:,j) = (a+c*b)\d(:,:,j); + end return; end -x=zeros(n,m); +x=zeros(n,m,p); [u,t]=schur(c); if exist('OCTAVE_VERSION') [aa,bb,qq,zz]=qz(full(a),full(b)); - d=qq'*d*u; + for j=1:p, + d(:,:,j)=qq'*d(:,:,j)*u; + end else [aa,bb,qq,zz]=qz(full(a),full(b),'real'); % available in Matlab version 6.0 - d=qq*d*u; + for j=1:p, + d(:,:,j)=qq*d(:,:,j)*u; + end end i = 1; +c = zeros(n,1,p); while i < m if t(i+1,i) == 0 if i == 1 - c = zeros(n,1); + c = zeros(n,1,p); else - c = bb*(x(:,1:i-1)*t(1:i-1,i)); + for j=1:p, + c(:,:,j) = bb*(x(:,1:i-1,j)*t(1:i-1,i)); + end end - x(:,i)=(aa+bb*t(i,i))\(d(:,i)-c); + x(:,i,:)=(aa+bb*t(i,i))\squeeze(d(:,i,:)-c); i = i+1; else if i == n - c = zeros(n,1); - c1 = zeros(n,1); + c = zeros(n,1,p); + c1 = zeros(n,1,p); else - c = bb*(x(:,1:i-1)*t(1:i-1,i)); - c1 = bb*(x(:,1:i-1)*t(1:i-1,i+1)); + for j=1:p, + c(:,:,j) = bb*(x(:,1:i-1,j)*t(1:i-1,i)); + c1(:,:,j) = bb*(x(:,1:i-1,j)*t(1:i-1,i+1)); + end end - z = [aa+bb*t(i,i) bb*t(i+1,i); bb*t(i,i+1) aa+bb*t(i+1,i+1)]... - \[d(:,i)-c;d(:,i+1)-c1]; - x(:,i) = z(1:n); - x(:,i+1) = z(n+1:end); + bigmat = ([aa+bb*t(i,i) bb*t(i+1,i); bb*t(i,i+1) aa+bb*t(i+1,i+1)]); + z = bigmat\squeeze([d(:,i,:)-c;d(:,i+1,:)-c1]); + x(:,i,:) = z(1:n,:); + x(:,i+1,:) = z(n+1:end,:); i = i + 2; end end if i == m - c = bb*(x(:,1:m-1)*t(1:m-1,m)); - x(:,m)=(aa+bb*t(m,m))\(d(:,m)-c); + for j=1:p, + c(:,:,j) = bb*(x(:,1:m-1,j)*t(1:m-1,m)); + end + aabbt = (aa+bb*t(m,m)); + x(:,m,:)=aabbt\squeeze(d(:,m,:)-c); +end +for j=1:p, + x(:,:,j)=zz*x(:,:,j)*u'; end -x=zz*x*u'; % 01/25/03 MJ corrected bug for i==m (sign of c in x determination) % 01/31/03 MJ added 'real' to qz call diff --git a/matlab/gensylv/sylvester3a.m b/matlab/gensylv/sylvester3a.m index f4aafe4b5..2200ef4d1 100644 --- a/matlab/gensylv/sylvester3a.m +++ b/matlab/gensylv/sylvester3a.m @@ -1,7 +1,7 @@ -function x=sylvester3a(x0,a,b,c,d) +function [x0, flag]=sylvester3a(x0,a,b,c,dd) % solves iteratively ax+bxc=d -% Copyright (C) 2005-2009 Dynare Team +% Copyright (C) 2005-2012 Dynare Team % % This file is part of Dynare. % @@ -20,15 +20,19 @@ function x=sylvester3a(x0,a,b,c,d) a_1 = inv(a); b = a_1*b; -d = a_1*d; -e = 1; -iter = 1; -while e > 1e-8 & iter < 500 - x = d-b*x0*c; - e = max(max(abs(x-x0))); - x0 = x; - iter = iter + 1; -end -if iter == 500 - warning('sylvester3a : Only accuracy of %10.8f is achieved after 500 iterations') +flag=0; +for j=1:size(dd,3), + d = a_1*dd(:,:,j); + e = 1; + iter = 1; + while e > 1e-8 && iter < 500 + x = d-b*x0(:,:,j)*c; + e = max(max(abs(x-x0(:,:,j)))); + x0(:,:,j) = x; + iter = iter + 1; + end + if iter == 500 + sprintf('sylvester3amr : Only accuracy of %10.8f is achieved after 500 iterations',e); + flag=1; + end end \ No newline at end of file diff --git a/matlab/gensylv/sylvester3amr.m b/matlab/gensylv/sylvester3amr.m deleted file mode 100644 index 4eb696296..000000000 --- a/matlab/gensylv/sylvester3amr.m +++ /dev/null @@ -1,38 +0,0 @@ -function [x0, flag]=sylvester3amr(x0,a,b,c,dd) -% solves iteratively ax+bxc=d - -% Copyright (C) 2012 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 . - -a_1 = inv(a); -b = a_1*b; -flag=0; -for j=1:size(dd,3), - d = a_1*dd(:,:,j); - e = 1; - iter = 1; - while e > 1e-8 && iter < 500 - x = d-b*x0(:,:,j)*c; - e = max(max(abs(x-x0(:,:,j)))); - x0(:,:,j) = x; - iter = iter + 1; - end - if iter == 500 - sprintf('sylvester3amr : Only accuracy of %10.8f is achieved after 500 iterations',e); - flag=1; - end -end \ No newline at end of file diff --git a/matlab/gensylv/sylvester3mr.m b/matlab/gensylv/sylvester3mr.m deleted file mode 100644 index 671fb4b21..000000000 --- a/matlab/gensylv/sylvester3mr.m +++ /dev/null @@ -1,100 +0,0 @@ -function x=sylvester3mr(a,b,c,d) -% solves a*x+b*x*c=d where d is [n x m x p] - -% Copyright (C) 2005-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 . - -n = size(a,1); -m = size(c,1); -if length(size(d))==2, - x=sylvester3(a,b,c,d); - return -end -p = size(d,3); -if n == 1 - for j=1:p, - x(:,:,j)=d(:,:,j)./(a*ones(1,m)+b*c); - end - return -end -if m == 1 - invacb = inv(a+c*b); - x = invacb*d; - return; -end -x=zeros(n,m,p); -[u,t]=schur(c); -if exist('OCTAVE_VERSION') - [aa,bb,qq,zz]=qz(full(a),full(b)); - for j=1:p, - d(:,:,j)=qq'*d(:,:,j)*u; - end -else - [aa,bb,qq,zz]=qz(full(a),full(b),'real'); % available in Matlab version 6.0 - for j=1:p, - d(:,:,j)=qq*d(:,:,j)*u; - end -end -i = 1; -c = zeros(n,1,p); -while i < m - if t(i+1,i) == 0 - if i == 1 - c = zeros(n,1,p); - else - for j=1:p, - c(:,:,j) = bb*(x(:,1:i-1,j)*t(1:i-1,i)); - end - end - % aabbtinv = inv(aa+bb*t(i,i)); - % x(:,i,:)=aabbtinv*squeeze(d(:,i,:)-c); - x(:,i,:)=(aa+bb*t(i,i))\squeeze(d(:,i,:)-c); - i = i+1; - else - if i == n - c = zeros(n,1,p); - c1 = zeros(n,1,p); - else - for j=1:p, - c(:,:,j) = bb*(x(:,1:i-1,j)*t(1:i-1,i)); - c1(:,:,j) = bb*(x(:,1:i-1,j)*t(1:i-1,i+1)); - end - end - % bigmatinv = inv([aa+bb*t(i,i) bb*t(i+1,i); bb*t(i,i+1) aa+bb*t(i+1,i+1)]); - % z = bigmatinv * squeeze([d(:,i,:)-c;d(:,i+1,:)-c1]); - bigmat = ([aa+bb*t(i,i) bb*t(i+1,i); bb*t(i,i+1) aa+bb*t(i+1,i+1)]); - z = bigmat\squeeze([d(:,i,:)-c;d(:,i+1,:)-c1]); - x(:,i,:) = z(1:n,:); - x(:,i+1,:) = z(n+1:end,:); - i = i + 2; - end -end -if i == m - for j=1:p, - c(:,:,j) = bb*(x(:,1:m-1,j)*t(1:m-1,m)); - end - % aabbtinv = inv(aa+bb*t(m,m)); - % x(:,m,:)=aabbtinv*squeeze(d(:,m,:)-c); - aabbt = (aa+bb*t(m,m)); - x(:,m,:)=aabbt\squeeze(d(:,m,:)-c); -end -for j=1:p, - x(:,:,j)=zz*x(:,:,j)*u'; -end - -% 01/25/03 MJ corrected bug for i==m (sign of c in x determination) -% 01/31/03 MJ added 'real' to qz call From 32de0fd675d6ee620d26237a015a5499a2e26f83 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Wed, 4 Apr 2012 16:41:49 +0200 Subject: [PATCH 6/7] fixed typos --- matlab/gensylv/sylvester3.m | 4 ++-- matlab/gensylv/sylvester3a.m | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/matlab/gensylv/sylvester3.m b/matlab/gensylv/sylvester3.m index b2b18f594..fcc5d8e78 100644 --- a/matlab/gensylv/sylvester3.m +++ b/matlab/gensylv/sylvester3.m @@ -1,4 +1,4 @@ -function x=sylvester3mr(a,b,c,d) +function x=sylvester3(a,b,c,d) % solves a*x+b*x*c=d where d is [n x m x p] % Copyright (C) 2005-2009 Dynare Team @@ -21,6 +21,7 @@ function x=sylvester3mr(a,b,c,d) n = size(a,1); m = size(c,1); p = size(d,3); +x=zeros(n,m,p); if n == 1 for j=1:p, x(:,:,j)=d(:,:,j)./(a*ones(1,m)+b*c); @@ -33,7 +34,6 @@ if m == 1 end return; end -x=zeros(n,m,p); [u,t]=schur(c); if exist('OCTAVE_VERSION') [aa,bb,qq,zz]=qz(full(a),full(b)); diff --git a/matlab/gensylv/sylvester3a.m b/matlab/gensylv/sylvester3a.m index 2200ef4d1..00724dc38 100644 --- a/matlab/gensylv/sylvester3a.m +++ b/matlab/gensylv/sylvester3a.m @@ -32,7 +32,7 @@ for j=1:size(dd,3), iter = iter + 1; end if iter == 500 - sprintf('sylvester3amr : Only accuracy of %10.8f is achieved after 500 iterations',e); + sprintf('sylvester3a : Only accuracy of %10.8f is achieved after 500 iterations',e); flag=1; end end \ No newline at end of file From d4723038becc8fc8849fe93188bfe2c7d1c93c9b Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Wed, 4 Apr 2012 16:43:44 +0200 Subject: [PATCH 7/7] 1) fixed infinite loop 2) fixed names of sylvester functions --- matlab/dynare_identification.m | 2 +- matlab/getH.m | 16 ++++++++++------ 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/matlab/dynare_identification.m b/matlab/dynare_identification.m index b9384bda9..aaf6ac279 100644 --- a/matlab/dynare_identification.m +++ b/matlab/dynare_identification.m @@ -154,7 +154,7 @@ end SampleSize = options_ident.prior_mc; -if ~(exist('sylvester3mr','file')==2), +if ~(exist('sylvester3','file')==2), dynareroot = strrep(which('dynare'),'dynare.m',''); addpath([dynareroot 'gensylv']) diff --git a/matlab/getH.m b/matlab/getH.m index 5d9316ed4..71a012437 100644 --- a/matlab/getH.m +++ b/matlab/getH.m @@ -355,10 +355,12 @@ else % generalized sylvester equation elem(:,:,j) = (Dg0(:,:,j)-Dg1(:,:,j)*A); d(:,:,j) = Dg2(:,:,j)-elem(:,:,j)*A; end - xx=sylvester3mr(a,b,c,d); + xx=sylvester3(a,b,c,d); flag=1; - while flag, - [xx, flag]=sylvester3amr(xx,a,b,c,d); + icount=0; + while flag && icount<4, + [xx, flag]=sylvester3a(xx,a,b,c,d); + icount=icount+1; end H=zeros(m*m+m*(m+1)/2,param_nbr+length(indexo)); if nargout>1, @@ -439,10 +441,12 @@ if nargout > 5, d(:,:,jcount) = elem1+elem2; end end - xx2=sylvester3mr(a,b,c,d); + xx2=sylvester3(a,b,c,d); flag=1; - while flag, - [xx2, flag]=sylvester3amr(xx2,a,b,c,d); + icount=0; + while flag && icount<4, + [xx2, flag]=sylvester3a(xx2,a,b,c,d); + icount = icount + 1; end jcount = 0; for j=1:param_nbr,