sylvester3

 solves a*x+b*x*c=d where d is [n x m x p]

0001 function x=sylvester3(a,b,c,d)
0002 % solves a*x+b*x*c=d where d is [n x m x p]
0003 
0004 % Copyright (C) 2005-2009 Dynare Team
0005 %
0006 % This file is part of Dynare.
0007 %
0008 % Dynare is free software: you can redistribute it and/or modify
0009 % it under the terms of the GNU General Public License as published by
0010 % the Free Software Foundation, either version 3 of the License, or
0011 % (at your option) any later version.
0012 %
0013 % Dynare is distributed in the hope that it will be useful,
0014 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0015 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0016 % GNU General Public License for more details.
0017 %
0018 % You should have received a copy of the GNU General Public License
0019 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0020 
0021 n = size(a,1);
0022 m = size(c,1);
0023 p = size(d,3);
0024 x=zeros(n,m,p);
0025 if n == 1
0026     for j=1:p,
0027         x(:,:,j)=d(:,:,j)./(a*ones(1,m)+b*c);
0028     end
0029     return
0030 end
0031 if m == 1
0032     for j=1:p,
0033         x(:,:,j) = (a+c*b)\d(:,:,j);
0034     end
0035     return;
0036 end
0037 [u,t]=schur(c);
0038 if exist('OCTAVE_VERSION')
0039     [aa,bb,qq,zz]=qz(full(a),full(b));
0040     for j=1:p,
0041         if octave_ver_less_than('3.4.0')
0042             d(:,:,j)=qq'*d(:,:,j)*u;
0043         else
0044             d(:,:,j)=qq*d(:,:,j)*u;
0045         end
0046     end
0047 else
0048     [aa,bb,qq,zz]=qz(full(a),full(b),'real'); % available in Matlab version 6.0
0049     for j=1:p,
0050         d(:,:,j)=qq*d(:,:,j)*u;
0051     end
0052 end
0053 i = 1;
0054 c = zeros(n,1,p);
0055 while i < m
0056     if t(i+1,i) == 0
0057         if i == 1
0058             c = zeros(n,1,p);
0059         else
0060             for j=1:p,
0061                 c(:,:,j) = bb*(x(:,1:i-1,j)*t(1:i-1,i));
0062             end
0063         end
0064         x(:,i,:)=(aa+bb*t(i,i))\squeeze(d(:,i,:)-c);
0065         i = i+1;
0066     else
0067         if i == n
0068             c = zeros(n,1,p);
0069             c1 = zeros(n,1,p);
0070         else
0071             for j=1:p,
0072                 c(:,:,j) = bb*(x(:,1:i-1,j)*t(1:i-1,i));
0073                 c1(:,:,j) = bb*(x(:,1:i-1,j)*t(1:i-1,i+1));
0074             end
0075         end
0076         bigmat = ([aa+bb*t(i,i) bb*t(i+1,i); bb*t(i,i+1) aa+bb*t(i+1,i+1)]);
0077         z = bigmat\squeeze([d(:,i,:)-c;d(:,i+1,:)-c1]);
0078         x(:,i,:) = z(1:n,:);
0079         x(:,i+1,:) = z(n+1:end,:);
0080         i = i + 2;
0081     end
0082 end
0083 if i == m
0084     for j=1:p,
0085         c(:,:,j) = bb*(x(:,1:m-1,j)*t(1:m-1,m));
0086     end
0087     aabbt = (aa+bb*t(m,m));
0088     x(:,m,:)=aabbt\squeeze(d(:,m,:)-c);
0089 end
0090 for j=1:p,
0091     x(:,:,j)=zz*x(:,:,j)*u';
0092 end
0093 
0094 % 01/25/03 MJ corrected bug for i==m (sign of c in x determination)
0095 % 01/31/03 MJ added 'real' to qz call