# copyright Barry L. Nelson and Michael R. Taaffe 2000

with(linalg):

PhPhInfNetSymbolic := proc(A, lambda, B, mu, K, Pr)

#
# this procedure returns a symbolic representation
# of the first and second moment differential
# equations for a K-node network of Ph_t/Ph_t/Infinity queues 
# with a single arrival process 
# 
# initial conditions for the queue are empty and idle
#
# this version is designed to work with Maple V release 5
#              
# last update 7/28/00
#
# Inputs:
# A = Markov chain transition matrix for the arrival process
# lambda = arrival rate vector
# B = Markov chain transition matrix for the service process
#     [vector of matrices of dimension K]
# mu = service rate vector
#     [vector of vectors of dimension K]
# K = number of  nodes in the network
# Pr = Markov routing matrix
#

	local BB, Bmu, i, j, g, h, v, r, offi, offj, mB;

# first check to make sure it is not a single-node problem
	if K = 1 then
	   RETURN(PhPhInfSymbolic(A, lambda, B, mu));

# network problem; first extract key dimensions
	else
	   v := vector(K);
	   v[1] := rowdim(B[1]) - 1;
	   mB := v[1];
	   for r from 2 to K do
	      v[r] := rowdim(B[r]) - 1;
	      mB := mB + v[r];
	   od;

	   BB := matrix(mB+1, mB+1, 0);

# build the composite service-rate vector
	   Bmu := vector(mB);
           offi := 0;
           for i from 1 to K do
	      for g from 1 to v[i] do
	         Bmu[offi + g] := mu[i][g];
 	      od;
	      offi := offi + v[i]
	   od; 

# build the composite phase transition matrix
	   offi := 0; offj := 0;

#      build rows and columns 1 through mB
	   for i from 1 to K do
	      for j from 1 to K do
	         for g from 1 to v[i] do
	            for h from 1 to v[j] do
	               BB[offi + g, offj + h] := Pr[i,j]*B[i][g, v[i]+1]*B[j][v[j]+1,h];
	            od;
	         od;
	         offj := offj + v[j];
              od;
	      offi := offi + v[i];
	      offj := 0;
           od;

#      build the last row
	   offj := 0;
	   for j from 1 to K do
	      for h from 1 to v[j] do
	         BB[mB+1, offj + h] := Pr[K+1,j]*B[j][v[j]+1,h];
              od;
              offj := offj + v[j];
           od;

#      build the last column
	   offi := 0;
	   for i from 1 to K do
	      for g from 1 to v[i] do
	         BB[offi + g, mB+1] := Pr[i, K+1]*B[i][g, v[i]+1];
              od;
	      offi := offi + v[i];
           od;

	   BB[mB+1, mB+1] := 0;

#      add the terms representing immediate feedback
	   BB := matadd(BB, MatDiag(B, K));

	   RETURN(PhPhInfSymbolic(A, lambda, BB, Bmu));
	fi;
	end;




MatDiag := proc(x, k)

# auxiliary proc to build a matrix with other matrices on the diagonal

	local m, lower, i, j, l, d, M;

	m := rowdim(x[1]) - 1;
  	for l from 2 to k do
	  m := m + rowdim(x[l]) - 1;
	od;
	M := matrix(m+1, m+1, 0);

	lower := 0;

	for l from 1 to k do
	   d := rowdim(x[l]) - 1;
	   for i from 1 to d do
	      for j from 1 to d do
	         M[lower + i, lower + j] := x[l][i,j];
              od;
	   od;
	   lower := lower + d;
	od;
	RETURN(M);
	end;