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