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