function [M,status] = DoublingAlgorithmLyapunov(S,W,MaxIter,ConvValue)
% DoublingAlgorithmLyapunov: Approximates the solution to the simplified version of the Sylvester equations.
%
% USAGE:
%
%       [M,status] = DoublingAlgorithmLyapunov(S,W,maxIter,ConvValue)
%
% The simplified version of the Sylvester equations (the Lyapunov equations) is
%
%    M = W + S*M*S'
%
% The matrix M that solves these equation is here approximated through the Doubling Algorithm. The algorithm
% is simply:
%
%    gamma(k+1) = gamma(k) + alpha(k)*gamma(k)*alpha(k)'
%      alpha(k) = alpha(k-1)*alpha(k-1)
%
% where the initial values are:
%
%      gamma(0) = W
%      alpha(0) = S
%
% For k large enough, gamma(k+1) is approximately equal to M.
%
% The convergence criterion used by the DoublingAlgorithmLyapunov function is simply
%
%    norm( gamma(k+1) - gamma(k) ) <= ConvValue
%
% REQUIRED INPUTS:  S (nxn) matrix
%
%                   W (nxn) matrix
%
%                   MaxIter, a positive integer measuring the maximum number of iterations
%
%                   ConvValue, a positive real number measuring convergence
%
% REQUIRED OUTPUTS: M (nxn) matrix that approximates the solution to the Lyapunov equations
%                     status, a boolean variable taking on 0 if the iterations have converged, and 1
%                     otherwise.
%
%                   status (Boolean) takes the value 0 if the algorithm has converged and 1 otherwise.
%
%
%                       Written by: Anders Warne
%                                   New Area Wide Model Project
%                                   DG-R/EMO
%                                   European Central Bank (ECB)
%                                   Email: anders.warne@ecb.europa.eu
%                                   Copyright © 2006-2008 European Central Bank.
%
%                       First version: July 17, 2006.
%                        This version: December 5, 2008.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% LICENSE INFORMATION:
%
%      YADA is free software: you can redistribute it and/or modify
%      it under the terms of the GNU General Public License as published by
%      the Free Software Foundation, either version 3 of the License, or
%      (at your option) any later version.
%
%      This program is distributed in the hope that it will be useful,
%      but WITHOUT ANY WARRANTY; without even the implied warranty of
%      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%      GNU General Public License for more details.
%
%      You should have received a copy of the GNU General Public License
%      along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%      YADA is released under the GNU General Public License, Version 3,
%      29 June 2007 <http://www.gnu.org/licenses/>. The current release of
%      the program was last modified by the ECB on the "This version" date
%      above.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% CHANGELOG:
%
% * 19-7-2006:  Updated function documentation.
%
% * 24-10-2006: Changes the convergence criterion from max,max,abs to norm.
%
% * 25-10-2006: Updated the documentation.
%
% * 2-11-2007:  Updated the documentation.
%
% * 23-5-2008:  Updated the documentation.
%
% * 5-12-2008:  Updated the documentation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

status = 1;
%
% Initialize gamma and alpha
%
gamma_old = W;
alpha = S;
%
% start up the iterations
%
for i=1:MaxIter;
   gamma_new = gamma_old+(alpha*gamma_old*alpha');
   alpha = alpha*alpha;
   M = gamma_new;
   %
   % check for convergence
   %
   if norm(M-gamma_old)<=ConvValue;
%   if (max(max(abs(M-gamma_old))))<=ConvValue;
      status = 0;
      break;
   end;
   gamma_old = gamma_new;
end;

%
% end of DoublingAlgorithmLyapunov.m
%
