1function [X, info] = cycle_reduction(A0, A1, A2, cvg_tol, ch) % --*-- Unitary tests --*-- 2 3%@info: 4%! @deftypefn {Function File} {[@var{X}, @var{info}] =} cycle_reduction (@var{A0},@var{A1},@var{A2},@var{cvg_tol},@var{ch}) 5%! @anchor{cycle_reduction} 6%! @sp 1 7%! Solves the quadratic matrix equation A2*X^2 + A1*X + A0 = 0. 8%! @sp 2 9%! @strong{Inputs} 10%! @sp 1 11%! @table @ @var 12%! @item A0 13%! Square matrix of doubles, n*n. 14%! @item A1 15%! Square matrix of doubles, n*n. 16%! @item A2 17%! Square matrix of doubles, n*n. 18%! @item cvg_tol 19%! Scalar double, tolerance parameter. 20%! @item ch 21%! Any matlab object, if not empty the solution is checked. 22%! @end table 23%! @sp 1 24%! @strong{Outputs} 25%! @sp 1 26%! @table @ @var 27%! @item X 28%! Square matrix of doubles, n*n, solution of the matrix equation. 29%! @item info 30%! Scalar integer, if nonzero the algorithm failed in finding the solution of the matrix equation. 31%! @end table 32%! @sp 2 33%! @strong{This function is called by:} 34%! @sp 2 35%! @strong{This function calls:} 36%! @sp 2 37%! @strong{References:} 38%! @sp 1 39%! D.A. Bini, G. Latouche, B. Meini (2002), "Solving matrix polynomial equations arising in queueing problems", Linear Algebra and its Applications 340, pp. 222-244 40%! @sp 1 41%! D.A. Bini, B. Meini (1996), "On the solution of a nonlinear matrix equation arising in queueing problems", SIAM J. Matrix Anal. Appl. 17, pp. 906-926. 42%! @sp 2 43%! @end deftypefn 44%@eod: 45 46% Copyright (C) 2013-2017 Dynare Team 47% 48% This file is part of Dynare. 49% 50% Dynare is free software: you can redistribute it and/or modify 51% it under the terms of the GNU General Public License as published by 52% the Free Software Foundation, either version 3 of the License, or 53% (at your option) any later version. 54% 55% Dynare is distributed in the hope that it will be useful, 56% but WITHOUT ANY WARRANTY; without even the implied warranty of 57% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 58% GNU General Public License for more details. 59% 60% You should have received a copy of the GNU General Public License 61% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 62 63max_it = 300; 64it = 0; 65info = 0; 66X = []; 67crit = Inf; 68A0_0 = A0; 69Ahat1 = A1; 70if (nargin == 5 && ~isempty(ch) ) 71 A1_0 = A1; 72 A2_0 = A2; 73end 74n = length(A0); 75id0 = 1:n; 76id2 = id0+n; 77 78cont = 1; 79while cont 80 tmp = ([A0; A2]/A1)*[A0 A2]; 81 A1 = A1 - tmp(id0,id2) - tmp(id2,id0); 82 A0 = -tmp(id0,id0); 83 A2 = -tmp(id2,id2); 84 Ahat1 = Ahat1 -tmp(id2,id0); 85 crit = norm(A0,1); 86 if crit < cvg_tol 87 % keep iterating until condition on A2 is met 88 if norm(A2,1) < cvg_tol 89 cont = 0; 90 end 91 elseif isnan(crit) || it == max_it 92 if crit < cvg_tol 93 info(1) = 4; 94 info(2) = log(norm(A2,1)); 95 else 96 info(1) = 3; 97 info(2) = log(norm(A1,1)); 98 end 99 return 100 end 101 it = it + 1; 102end 103 104X = -Ahat1\A0_0; 105 106if (nargin == 5 && ~isempty(ch) ) 107 %check the solution 108 res = A0_0 + A1_0 * X + A2_0 * X * X; 109 if (sum(sum(abs(res))) > cvg_tol) 110 disp(['the norm residual of the residu ' num2str(res) ' compare to the tolerance criterion ' num2str(cvg_tol)]); 111 end 112end 113 114%@test:1 115%$ 116%$ t = zeros(3,1); 117%$ 118%$ % Set the dimension of the problem to be solved. 119%$ n = 500; 120%$ 121%$ % Set the equation to be solved 122%$ A = eye(n); 123%$ B = diag(30*ones(n,1)); B(1,1) = 20; B(end,end) = 20; B = B - diag(10*ones(n-1,1),-1); B = B - diag(10*ones(n-1,1),1); 124%$ C = diag(15*ones(n,1)); C = C - diag(5*ones(n-1,1),-1); C = C - diag(5*ones(n-1,1),1); 125%$ 126%$ % Solve the equation with the cycle reduction algorithm 127%$ try 128%$ tic; X1 = cycle_reduction(C,B,A,1e-7); elapsedtime = toc; 129%$ disp(['Elapsed time for cycle reduction algorithm is: ' num2str(elapsedtime) ' (n=' int2str(n) ').']) 130%$ t(1) = 1; 131%$ catch 132%$ % nothing to do here. 133%$ end 134%$ 135%$ % Solve the equation with the logarithmic reduction algorithm 136%$ try 137%$ tic; X2 = logarithmic_reduction(A,B,C,1e-16,100); elapsedtime = toc; 138%$ disp(['Elapsed time for logarithmic reduction algorithm is: ' num2str(elapsedtime) ' (n=' int2str(n) ').']) 139%$ t(2) = 1; 140%$ catch 141%$ % nothing to do here. 142%$ end 143%$ 144%$ % Check the results. 145%$ if t(1) && t(2) 146%$ t(3) = dassert(X1,X2,1e-12); 147%$ end 148%$ 149%$ T = all(t); 150%@eof:1