From f455557376dad9f39eb75a5c159dd230a8db11d1 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Wed, 4 Apr 2012 16:38:22 +0200 Subject: [PATCH] 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