1function [x,u] = lyapunov_symm(a,b,lyapunov_fixed_point_tol,qz_criterium,lyapunov_complex_threshold,method,debug)  % --*-- Unitary tests --*--
2% Solves the Lyapunov equation x-a*x*a' = b, for b and x symmetric matrices.
3% If a has some unit roots, the function computes only the solution of the stable subsystem.
4%
5% INPUTS:
6%   a                           [double]    n*n matrix.
7%   b                           [double]    n*n matrix.
8%   qz_criterium                [double]    unit root threshold for eigenvalues
9%   lyapunov_fixed_point_tol    [double]    convergence criteria for fixed_point algorithm.
10%   lyapunov_complex_threshold  [double]    scalar, complex block threshold for the upper triangular matrix T.
11%   method                      [integer]   Scalar, if method=0 [default] then U, T, n and k are not persistent.
12%                                                      method=1 then U, T, n and k are declared as persistent
13%                                                               variables and the Schur decomposition is triggered.
14%                                                      method=2 then U, T, n and k are declared as persistent
15%                                                               variables and the Schur decomposition is not performed.
16%                                                      method=3 fixed point method
17% OUTPUTS
18%   x:      [double]    m*m solution matrix of the lyapunov equation, where m is the dimension of the stable subsystem.
19%   u:      [double]    Schur vectors associated with unit roots
20%
21% ALGORITHM
22%   Uses reordered Schur decomposition (Bartels-Stewart algorithm)
23%   [method<3] or a fixed point algorithm (method==3)
24%
25% SPECIAL REQUIREMENTS
26%   None
27
28% Copyright (C) 2006-2017 Dynare Team
29%
30% This file is part of Dynare.
31%
32% Dynare is free software: you can redistribute it and/or modify
33% it under the terms of the GNU General Public License as published by
34% the Free Software Foundation, either version 3 of the License, or
35% (at your option) any later version.
36%
37% Dynare is distributed in the hope that it will be useful,
38% but WITHOUT ANY WARRANTY; without even the implied warranty of
39% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
40% GNU General Public License for more details.
41%
42% You should have received a copy of the GNU General Public License
43% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
44
45if nargin<6 || isempty(method)
46    method = 0;
47end
48
49if nargin<7
50    debug = 0;
51end
52
53if method == 3
54    persistent X method1;
55    if ~isempty(method1)
56        method = method1;
57    end
58    if debug
59        fprintf('lyapunov_symm:: [method=%d] \n',method);
60    end
61    if method == 3
62        %tol = 1e-10;
63        it_fp = 0;
64        evol = 100;
65        if isempty(X) || length(X)~=length(b)
66            X = b;
67            max_it_fp = 2000;
68        else
69            max_it_fp = 300;
70        end
71        at = a';
72        %fixed point iterations
73        while evol >  lyapunov_fixed_point_tol && it_fp < max_it_fp
74            X_old = X;
75            X = a * X * at + b;
76            evol = max(sum(abs(X - X_old))); %norm_1
77                                             %evol = max(sum(abs(X - X_old)')); %norm_inf
78            it_fp = it_fp + 1;
79        end
80        if debug
81            fprintf('lyapunov_symm:: lyapunov fixed_point iterations=%d norm=%g\n',it_fp,evol);
82        end
83        if it_fp >= max_it_fp
84            disp(['lyapunov_symm:: convergence not achieved in solution of Lyapunov equation after ' int2str(it_fp) ' iterations, switching method from 3 to 0']);
85            method1 = 0;
86            method = 0;
87        else
88            method1 = 3;
89            x = X;
90            return
91        end
92    end
93end
94
95if method
96    persistent U T k n
97else
98    %    if exist('U','var')
99    %        clear('U','T','k','n')
100    %    end
101end
102
103u = [];
104
105if size(a,1) == 1
106    x=b/(1-a*a);
107    return
108end
109
110if method<2
111    [U,T] = schur(a);
112    e1 = abs(ordeig(T)) > 2-qz_criterium;
113    k = sum(e1);       % Number of unit roots.
114    n = length(e1)-k;  % Number of stationary variables.
115    if k > 0
116        % Selects stable roots
117        [U,T] = ordschur(U,T,e1);
118        T = T(k+1:end,k+1:end);
119    end
120end
121
122B = U(:,k+1:end)'*b*U(:,k+1:end);
123
124x = zeros(n,n);
125i = n;
126
127while i >= 2
128    if abs(T(i,i-1))<lyapunov_complex_threshold
129        if i == n
130            c = zeros(n,1);
131        else
132            c = T(1:i,:)*(x(:,i+1:end)*T(i,i+1:end)') + ...
133                T(i,i)*T(1:i,i+1:end)*x(i+1:end,i);
134        end
135        q = eye(i)-T(1:i,1:i)*T(i,i);
136        x(1:i,i) = q\(B(1:i,i)+c);
137        x(i,1:i-1) = x(1:i-1,i)';
138        i = i - 1;
139    else
140        if i == n
141            c = zeros(n,1);
142            c1 = zeros(n,1);
143        else
144            c = T(1:i,:)*(x(:,i+1:end)*T(i,i+1:end)') + ...
145                T(i,i)*T(1:i,i+1:end)*x(i+1:end,i) + ...
146                T(i,i-1)*T(1:i,i+1:end)*x(i+1:end,i-1);
147            c1 = T(1:i,:)*(x(:,i+1:end)*T(i-1,i+1:end)') + ...
148                 T(i-1,i-1)*T(1:i,i+1:end)*x(i+1:end,i-1) + ...
149                 T(i-1,i)*T(1:i,i+1:end)*x(i+1:end,i);
150        end
151        q = [  eye(i)-T(1:i,1:i)*T(i,i) ,  -T(1:i,1:i)*T(i,i-1) ; ...
152               -T(1:i,1:i)*T(i-1,i)     ,   eye(i)-T(1:i,1:i)*T(i-1,i-1) ];
153        z =  q\[ B(1:i,i)+c ; B(1:i,i-1) + c1 ];
154        x(1:i,i) = z(1:i);
155        x(1:i,i-1) = z(i+1:end);
156        x(i,1:i-1) = x(1:i-1,i)';
157        x(i-1,1:i-2) = x(1:i-2,i-1)';
158        i = i - 2;
159    end
160end
161if i == 1
162    c = T(1,:)*(x(:,2:end)*T(1,2:end)') + T(1,1)*T(1,2:end)*x(2:end,1);
163    x(1,1) = (B(1,1)+c)/(1-T(1,1)*T(1,1));
164end
165x = U(:,k+1:end)*x*U(:,k+1:end)';
166u = U(:,1:k);