Commit f4555573 by Marco Ratto

### compact all sylvester routines, with generic dimension of d argument

parent 76093094
 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
 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
 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
 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
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!