%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    Copyright (c) 2004.  Takayuki Osogami.  All Rights Reserved.   %
%                                                                   %
% Permission to use, copy, modify, and distribute this source code  %
% and its documentation for any purpose, without fee, and without   %
% a written agreement, is hereby granted, provided that the above   %
% copyright notice, this paragraph and the following two paragraphs %
% appear in all copies, modifications, and distributions.           %
% Created by Takayuki Osogami, Department of Computer Science,      %
% Carnegie Mellon University.                                       %
%                                                                   %
% IN NO EVENT SHALL THE AUTHOR BE LIABLE TO ANY PARTY FOR ANY       %
% DAMAGES, INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS   %
% SOURCE CODE AND ITS DOCUMENTATION, EVEN IF THE AUTHOR HAS BEEN    %
% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.                        %
%                                                                   %
% THE AUTHOR SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING,      %
% BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY     %
% AND FITNESS FOR A PARTICULAR PURPOSE.  THE SOURCE CODE AND        %
% ACCOMPANYING DOCUMENTATION, IF ANY, PROVIDED HEREUNDER IS         %
% PROVIDED "AS IS."  THE AUTHOR HAS NO OBLIGATION TO PROVIDE        %
% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Usage: [G] = Gmatrix(A0,A1,A2,epsilon)
% 
% This program iteratively computes the matrix G
% given A2, A1, and A0 (of the repeating part).
%
% INPUT: A0,A1,A2: A matrices for the repeating part
%
% OUTPUT: G: G matrix for the repeating part
%
%
% author: Takayuki Osogami
%         Department of Computer Science
%         Carnegie Mellon University
%         osogami@cs.cmu.edu
% date: April 28, 2004

function [G] = Gmatrix(A0,A1,A2,epsilon)

S = size(A0);
dim = S(1);
G = zeros(dim,dim);

%% exit after this many iteration: did not converge %%
maxiterate = 100000;

%% use quadratic convergent algorithm %%
QUADRATIC = 1;

if( QUADRATIC == 0 )

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % linear congergent algorithm % 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 i = 1;
 D = 1000;

 while( D > epsilon )

    i = i + 1;
    if( i > maxiterate )
     maxiterate_Gmatrix = 1
     break;
    end

    old = G;
    G = - (A1)^(-1) * A0 * G * G  - (A1)^(-1) * A2;
    diff = 0;
    for k = 1:dim
      for l = 1:dim
        diff = max(diff, abs(old(k,l) - G(k,l)));
      end
    end
    D = diff;

 end

else

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % quadratic convergent algorithm %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 i = 1;
 D = 1000;
 H = (-A1)^(-1) * A0;
 L = (-A1)^(-1) * A2;
 G = L;
 T = H;

 while( D > epsilon )

    if( i > maxiterate )
     maxiterate_Gmatrix = 1
     break;
    end

    U = H*L + L*H;
    M = H^2;
    S = size(U);
    E = eye(S(1));
    H = (E-U)^(-1) * M;
    M = L^2;
    L = (E-U)^(-1) * M;
    G = G + T*L;
    T = T*H;

    S = size(G);
    I = ones(S(1),1);
    error = I - G*I;

    diff = 0;
    for k = 1:S
      diff = max(diff,abs(error(k)));
    end
    D = diff;

 end

end