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