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