1function X = fastgensylv(A, B, C, D, tol,maxit,X0)
2
3%@info:
4%! @deftypefn {Function File} {[@var{X1}, @var{info}] =} fastgensylv (@var{A},@var{B},@var{C},@var{tol},@var{maxit},@var{X0})
5%! @anchor{fastgensylv}
6%! @sp 1
7%! Solves the Sylvester equation A * X + B * X * C + D = 0 for X.
8%! @sp 2
9%! @strong{Inputs}
10%! @sp 1
11%! @table @ @var
12%! @item A
13%! Square matrix of doubles, n*n.
14%! @item B
15%! Square matrix of doubles, n*n.
16%! @item C
17%! Square matrix of doubles, n*n.
18%! @item tol
19%! Scalar double, tolerance parameter.
20%! @item maxit
21%! Integer scalar, maximum number of iterations.
22%! @item X0
23%! Square matrix of doubles, n*n, initial condition.
24%! @end table
25%! @sp 1
26%! @strong{Outputs}
27%! @sp 1
28%! @table @ @var
29%! @item X
30%! Square matrix of doubles, n*n, solution of the matrix equation.
31%! @item info
32%! Scalar integer, if nonzero the algorithm failed in finding the solution of the matrix equation.
33%! @end table
34%! @sp 2
35%! @strong{This function is called by:}
36%! @sp 2
37%! @strong{This function calls:}
38%! @sp 2
39%! @end deftypefn
40%@eod:
41
42% Copyright (C) 2012-2017 Dynare Team
43%
44% This file is part of Dynare.
45%
46% Dynare is free software: you can redistribute it and/or modify
47% it under the terms of the GNU General Public License as published by
48% the Free Software Foundation, either version 3 of the License, or
49% (at your option) any later version.
50%
51% Dynare is distributed in the hope that it will be useful,
52% but WITHOUT ANY WARRANTY; without even the implied warranty of
53% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
54% GNU General Public License for more details.
55%
56% You should have received a copy of the GNU General Public License
57% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
58
59if size(A,1)~=size(D,1) || size(A,1)~=size(B,1) || size(C,2)~=size(D,2)
60    error('fastgensylv:: Dimension error!')
61end
62
63if nargin<7 || isempty(X0)
64    X = zeros(size(A,2),size(C,1));
65elseif nargin==7 && ~isempty(X0)
66    X = X0;
67end
68
69kk = 0;
70cc = 1+tol;
71
72iA = inv(A);
73Z = - (B * X * C + D);
74
75while kk<=maxit && cc>tol
76    X = iA * Z;
77    Z_old = Z;
78    Z = - (B * X * C + D);
79    cc = max(sum(abs(Z-Z_old)));
80    kk = kk + 1;
81end
82
83if kk==maxit && cc>tol
84    error(['fastgensylv:: Convergence not achieved in fixed point solution of Sylvester equation after ' int2str(maxit) ' iterations']);
85end
86
87
88
89
90% function X = fastgensylv(A, B, C, D)
91% Solve the Sylvester equation:
92% A * X + B * X * C + D = 0
93% INPUTS
94%   A
95%   B
96%   C
97%   D
98%   block : block number (for storage purpose)
99%   tol : convergence criteria
100% OUTPUTS
101%   X solution
102%
103% ALGORITHM
104%   fixed point method
105%   MARLLINY MONSALVE (2008): "Block linear method for large scale
106%   Sylvester equations", Computational & Applied Mathematics, Vol 27, n°1,
107%   p47-59
108%   ||A^-1||.||B||.||C|| < 1 is a suffisant condition:
109%    - to get a unique solution for the Sylvester equation
110%    - to get a convergent fixed-point algorithm
111%
112% SPECIAL REQUIREMENTS
113%   none.
114