%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    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: [F,P0,P1,P2] = FPmatrix(A0,A1,A2)
%
% Given A matrices in a level of a QBD process,
% FP matrix returns F,P0,P1,P2, defined below:
%
%        [ P11*F11(x) P12*F12(x) ... ]
%        [ P21*F21(x)                ]
% A(x) = [    :                      ]
%        [                           ]
%        [                           ]
%
% is the definition of A(x) in a general semi-Markov chain.
% When it is a continuous time Markov chain,
%
%        [ F1(x)  0     ... 0 ]   [ P11 P12 ... ]
%        [  0    F2(x)        ]   [ P21 P22 ... ]
% A(x) = [  :                 ] * [  :          ]
%        [  :                 ]   [  :          ]
%        [  0    ...          ]   [             ]
%
% where Fi(x) is the exponential distribution function.
%
% We obtain P and F from the infinitesimal generators.
% F just records the rates instead of distribution functions;
% the rates are sufficient to figure out distribution functions.

function [F,P0,P1,P2] = FPmatrix(A0,A1,A2)

S = size(A1);
dim = S(1);

for i = 1:dim
 totalrate = - A1(i,i);

 F(i,i) = totalrate;

 P0(i,:) = A0(i,:) / totalrate;
 P2(i,:) = A2(i,:) / totalrate;

 P1(i,:) = A1(i,:) / totalrate;
 P1(i,i) = 0;

end