1function x=sylvester3(a,b,c,d) 2% solves a*x+b*x*c=d where d is [n x m x p] 3 4% Copyright (C) 2005-2017 Dynare Team 5% 6% This file is part of Dynare. 7% 8% Dynare is free software: you can redistribute it and/or modify 9% it under the terms of the GNU General Public License as published by 10% the Free Software Foundation, either version 3 of the License, or 11% (at your option) any later version. 12% 13% Dynare is distributed in the hope that it will be useful, 14% but WITHOUT ANY WARRANTY; without even the implied warranty of 15% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 16% GNU General Public License for more details. 17% 18% You should have received a copy of the GNU General Public License 19% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 20 21n = size(a,1); 22m = size(c,1); 23p = size(d,3); 24x=zeros(n,m,p); 25if n == 1 26 for j=1:p 27 x(:,:,j)=d(:,:,j)./(a*ones(1,m)+b*c); 28 end 29 return 30end 31if m == 1 32 for j=1:p 33 x(:,:,j) = (a+c*b)\d(:,:,j); 34 end 35 return 36end 37[u,t]=schur(c); 38if isoctave 39 [aa,bb,qq,zz]=qz(full(a),full(b)); 40 for j=1:p 41 d(:,:,j)=qq*d(:,:,j)*u; 42 end 43else 44 [aa,bb,qq,zz]=qz(full(a),full(b),'real'); % available in Matlab version 6.0 45 for j=1:p 46 d(:,:,j)=qq*d(:,:,j)*u; 47 end 48end 49i = 1; 50c = zeros(n,1,p); 51c1 = zeros(n,1,p); 52while i < m 53 if t(i+1,i) == 0 54 if i == 1 55 c = zeros(n,1,p); 56 else 57 for j=1:p 58 c(:,:,j) = bb*(x(:,1:i-1,j)*t(1:i-1,i)); 59 end 60 end 61 x(:,i,:)=(aa+bb*t(i,i))\squeeze(d(:,i,:)-c); 62 i = i+1; 63 else 64 if i == n 65 c = zeros(n,1,p); 66 c1 = zeros(n,1,p); 67 else 68 for j=1:p 69 c(:,:,j) = bb*(x(:,1:i-1,j)*t(1:i-1,i)); 70 c1(:,:,j) = bb*(x(:,1:i-1,j)*t(1:i-1,i+1)); 71 end 72 end 73 bigmat = ([aa+bb*t(i,i) bb*t(i+1,i); bb*t(i,i+1) aa+bb*t(i+1,i+1)]); 74 z = bigmat\squeeze([d(:,i,:)-c;d(:,i+1,:)-c1]); 75 x(:,i,:) = z(1:n,:); 76 x(:,i+1,:) = z(n+1:end,:); 77 i = i + 2; 78 end 79end 80if i == m 81 for j=1:p 82 c(:,:,j) = bb*(x(:,1:m-1,j)*t(1:m-1,m)); 83 end 84 aabbt = (aa+bb*t(m,m)); 85 x(:,m,:)=aabbt\squeeze(d(:,m,:)-c); 86end 87for j=1:p 88 x(:,:,j)=zz*x(:,:,j)*u'; 89end 90 91% 01/25/03 MJ corrected bug for i==m (sign of c in x determination) 92% 01/31/03 MJ added 'real' to qz call 93