1function X = gensylv_fp(A, B, C, D, block, tol)
2% function X = gensylv_fp(A, B, C, D)
3% Solve the Sylvester equation:
4% A * X + B * X * C + D = 0
5% INPUTS
6%   A
7%   B
8%   C
9%   D
10%   block : block number (for storage purpose)
11%   tol : convergence criteria
12% OUTPUTS
13%   X solution
14%
15% ALGORITHM
16%   fixed point method
17%   MARLLINY MONSALVE (2008): "Block linear method for large scale
18%   Sylvester equations", Computational & Applied Mathematics, Vol 27, n°1,
19%   p47-59
20%   ||A^-1||.||B||.||C|| < 1 is a suffisant condition:
21%    - to get a unique solution for the Sylvester equation
22%    - to get a convergent fixed-point algorithm
23%
24% SPECIAL REQUIREMENTS
25%   none.
26% Copyright (C) 1996-2017 Dynare Team
27%
28% This file is part of Dynare.
29%
30% Dynare is free software: you can redistribute it and/or modify
31% it under the terms of the GNU General Public License as published by
32% the Free Software Foundation, either version 3 of the License, or
33% (at your option) any later version.
34%
35% Dynare is distributed in the hope that it will be useful,
36% but WITHOUT ANY WARRANTY; without even the implied warranty of
37% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
38% GNU General Public License for more details.
39%
40% You should have received a copy of the GNU General Public License
41% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
42
43
44%tol = 1e-12;
45evol = 100;
46A1 = inv(A);
47eval(['persistent hxo_' int2str(block) ';']);
48hxo = eval(['hxo_' int2str(block) ';']);
49if isempty(hxo)
50    X = zeros(size(B, 2), size(C, 1));
51else
52    X = hxo;
53end
54it_fp = 0;
55maxit_fp = 1000;
56Z = - (B * X * C + D);
57while it_fp < maxit_fp && evol > tol
58    %X_old = X;
59    %X = - A1 * ( B * X * C + D);
60    %evol = max(max(abs(X - X_old)));
61    X = A1 * Z;
62    Z_old = Z;
63    Z = - (B * X * C + D);
64    evol = max(sum(abs(Z - Z_old))); %norm_1
65                                     %evol = max(sum(abs(Z - Z_old)')); %norm_inf
66    it_fp = it_fp + 1;
67end
68%fprintf('sylvester it_fp=%d evol=%g | ',it_fp,evol);
69if evol < tol
70    eval(['hxo_' int2str(block) ' = X;']);
71else
72    error(['convergence not achieved in fixed point solution of Sylvester equation after ' int2str(it_fp) ' iterations']);
73end