Force precision in sylvester equation solution, particularly useful for...

Force precision in sylvester equation solution, particularly useful for different behaviour of the function schur under octave

matlab/gensylv/sylvester3amr.m

+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 <http://www.gnu.org/licenses/>.
+
+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

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,