Computational Procedures for a Class of GI/D/ k Systems in Discrete Time

A class of discrete time GI/D/ k systems is considered for which the interarrival times have ﬁnite support and customers are served in ﬁrst-in ﬁrst-out (cid:2) FIFO (cid:3) order. The system is formulated as a single server queue with new general independent interarrival times and constant service duration by assuming cyclic assignment of customers to the identical servers. Then the queue length is set up as a quasi-birth-death (cid:2) QBD (cid:3) type Markov chain. It is shown that this transformed GI/D/1 system has special structures which make the computation of the matrix R simple and e ﬃ cient, thereby reducing the number of multiplications in each iteration signiﬁcantly. As a result we were able to keep the computation time very low. Moreover, use of the resulting structural properties makes the computation of the distribution of queue length of the transformed system e ﬃ cient. The computation of the distribution of waiting time is also shown to be simple by exploiting the special structures.


Introduction
In most communication systems studied in discrete time, the interarrival times of packets usually follow independent general distribution.Hence, this arrival of packets in an ATM switch can be represented as independent general distribution.On the other hand, multiple packets 53 byte cells are transmitted simultaneously through identical transmission lines in an ATM switch.Therefore, an ATM switch can be considered as a multiserver system with deterministic service time as the packets are of same size.Hence, the performance of an ATM switch can be analyzed by studying a GI/D/k system with finite support on interarrival times.
To the best of our knowledge GI/D/k systems in discrete time have never been analyzed with focus on computational efficiency in the literature.In this paper, matrixgeometric method is used to analyze a class of GI/D/k systems in discrete time with focus on computational efficiency.The class of GI/D/k systems considered in this paper has finite support for the interarrival times and the service is first-in first-out FIFO order.The idea used here to analyze the waiting time distribution of this multiserver system is due to Crommelin 1 .This idea is also used in analyzing the distribution of waiting time of multiserver systems with constant service time duration for both continuous and discrete times such as 2-20 .The idea is that the assignment of customers to servers in cyclic order does not change the distribution of waiting time in a multiserver system where customers are served in FIFO order in the identical servers with same constant service time.In this case, only an arbitrary server can be studied to compute the waiting time distribution.Crommelin 1 analyzed the waiting time distribution of an M/D/c system by converting the system to an E c /D/1 system.The work of Crommelin 1 is followed by some studies on deterministic service time multiserver systems in both continuous time and discrete time.Iversen 15 decomposed an M/D/rk queue with FIFO into E k /D/r queues with FIFO for analyzing waiting time distribution.Franx 8 developed expression for the waiting time distribution of the M/D/c queue by a full probabilistic analysis requiring neither generating functions nor Laplace transforms.Wittevrongel and Bruneel 19 used transformbased method to analyze multiserver system with constant service time and correlated arrival process.Nishimura 17 studied a MAP/D/N system in continuous time using spectral technique.Takine 18 analyzed a MAP/D/s system and computed the distribution of queue length by first characterizing the sojourn time distribution.Kim and Chaudhry 21 also computed the distribution of queue length of a discrete time GI/D/c queue after computing the distribution of waiting time by using distributional Little's law and transform-based method.Exact elapsed time from the time of last arrival is approximated with time average in 21 .Chaudhry et al. 7 analyzed the distribution of waiting time of a MAP/D/k system in discrete time.Later Alfa 3 gave a simple and more efficient computational scheme for the MAP/D/k system in discrete time.Alfa 2 also carried out algorithmic analysis for a discrete time BMAP/D/k system by exploiting its structural properties.Chaudhry et al. 7 and Alfa 2, 3 used matrix geometric method 22 for analysis.There is some work on the GI/D/c queues such as Wuyts and Bruneel 4 and Gao et al. 10-14 which used transformbased method.There are other algorithms for multiserver queues with constant service times such as GI/D/1 and GI/D/c queues by Chaudhry  In this paper, the cyclic assignment of customers to servers is used to model a class of GI/D/k systems in discrete time as a single server system with the assumption of first-in first-out FIFO service order.The modeling as single server system sets up the distribution of queue length as a quasi-birth-death QBD which has some structural properties.Analysis of the GI/D/k system is carried out efficiently by exploiting these structural properties in this paper.This paper has three contributions-reductions in the computational complexities of i the matrix R, ii the distribution of queue length of the transformed system, and iii the distribution of waiting time.The first contribution is that the time complexity of the computation of the matrix R is reduced by decreasing the number of multiplications in each iteration from O n 3 k 3 d 3 to O n 3 kd 2 while requiring the same number of iteration as the natural iteration.Here, n is the support of general interarrival times, k is the number of servers, and d is the duration of service.The second contribution is that the distribution of queue length of the transformed system is computed efficiently by exploiting the special structures of the system.The third and most important contribution is that the computation of the distribution of waiting time is simplified.Chaudhry et al. 7 and Alfa 2, 3 computed the distribution of waiting time using the technique for phase type service time distribution.However, it is shown in this paper that the computation of the distribution of waiting time for deterministic service time distribution in discrete time does not need the complicated steps of phase type service time distribution.The rest of the paper is organized as follows.Section 2 introduces the GI/D/k system and explains the modeling of this multiserver system into a single server system.Section 3 discusses the special structures of the matrices of our model and exploitation of those structures for efficient computation of the matrix R. The computation of the distributions of queue length and waiting time is explained in Sections 4 and 5, respectively.Some numerical examples are provided in Section 6 to show the suitability of our proposed method to compute the matrix R. Section 7 concludes the paper.

The GI/D/k System
The class of GI/D/k system in discrete time considered here has finite support for interarrival time which is of general distribution.Moreover, the customers are served in FIFO service order in the identical servers which have constant service time.
The interarrival times, A, of this multiserver system are general, independent, and identically distributed iid with distribution a v Pr {A v}, v 1, 2, . . . .
We let λ 1/E A be the mean arrival rate of customers to the system with 0 < λ ≤ 1.The interarrival times can be represented as a Markov chain where each state represents the elapsed time using the approach followed in the analysis of discrete time GI/G/1 and GI X /G/1 systems by Alfa 24 .General independent interarrival times can be represented as an absorbing Markov chain with a transition matrix P A T t 0 1,n 1 where T is a square matrix of order n, t is a column vector of order n, t 1 n − T1 n , 0 i,j is an i × j matrix of zeros, and 1 i is a column vector of ones of order i.The case of infinite interarrival times can be captured for n ∞.However, the interarrival times are finite if data is collected from practical systems.Therefore, general independent interarrival times with finite support n < ∞ are considered in this paper.The probabilities of interarrival times can be represented by a vector a a 1 , a 2 , . . ., a n .Here, a1 n 1 and a 0 0. This general arrival process can be defined by phase type distribution α, T of dimension n in elapsed time representation.The parameters α, T, and t are given as α 1, 0, 0 . . ., 0 , T i,j 1 − i j 0 a j / 1 − i−1 v 0 a v a i , for 1 ≤ i < n, j i 1, T i,j 0, for all j / i 1 and t a 1 a 2 • • • a n−1 1 where a i 1 − a i for 1 ≤ i ≤ n − 1 and W is the transpose of matrix W.
The general independent arrival process can be represented by two square matrices D 0 and D 1 of order n where D l ij l 0, 1 and 1 ≤ i, j ≤ n represents a transition from phase i to phase j with l arrivals.Now, D 0 and D 1 can be represented in terms of phase type distribution α, T .Here, D 0 T and D 1 tα t 0 n,n−1 .The first column of D 1 is nonzero and this column is the vector t.The matrices D 0 and D 1 are similar to the matrices of zero or one arrival for MAP distributed interarrival times.However, they are special cases because the matrices D 0 and D 1 only capture general and independently distributed interarrival times.
A discrete time Markov chain { N t , J t ; t 0, 1, 2, . . .; N t ∈ {0, 1, 2, . ..};J t ∈ {1, 2, 3, . . ., n}} can be considered to represent this general arrival process where N t is the number of arrivals up to and including time t and J t is the phase of the next arrival at time t.
The matrix D D 0 D 1 is stochastic.If π π 1 , π 2 , . . ., π n is the invariant probability vector for arrival in n different phases then π πD and π1 n 1.The arrival rate is given by λ πD 1 1 n πt.
There are k identical servers in a GI/D/k system.The service times are of constant duration d units d ≥ 1 .Service time can be represented as phase type distribution in elapsed time form with representation β, S of dimension d where β 1, 0, 0, . . ., 0 and , I j is an identity matrix of order j.Another column vector s 1 d − S1 d can be defined for this phase type distribution β, S .
The condition that ensures the stability of GI/D/k system is ρ < 1 where ρ λd/k.We assume that the system is stable i.e., ρ < 1 .

The Model of Arrivals to an Arbitrary Server
Let us assume that customers are served in first-in first-out FIFO order in this GI/D/k system.It can also be assumed that the servers are assigned with customers in a cyclic manner particularly during the idle times.This assignment of customers to servers does not affect the waiting times experienced by the customers as the servers are identical.Using the same approach of Alfa 2, 3 we can consider the case of the jth 1 ≤ j ≤ k server.In this approach if the first customer is assigned to the first server then the jth server will be assigned with j, j k, j 2k, . ..th customers.In this case we only need to study one server, an arbitrary server j 1 ≤ j ≤ k .This arbitrary server sees the arrival of customers to be according to another independent general process which we call k GI.This k GI arrival process is a k-fold convolution of original arrival process to the multiserver system.k GI can be defined by two square matrices C 0 and C 1 of order nk.The element C v ij v 0, 1 and 1 ≤ i, j ≤ nk represents transition from phase i to phase j with v arrivals.Both C 0 and C 1 can be considered to comprise of k 2 square block matrices of order n.In this representation C 0 i,j D 0 for i j, 1 ≤ i ≤ k, C 0 i,j D 1 for j i 1, 1 ≤ i < k and C 0 i,j 0 n,n for 1 ≤ i ≤ k and i / j or j / i 1. Similarly C 1 i,j D 1 for i k, j 1 and C 1 i,j 0 n,n in all other cases for 1 ≤ i, j ≤ k.The first column and the last row of C 0 are zero as the first column and last row of D 0 are zero.Only the last n rows of C 1 are nonzero and only the first element of each of these n rows is nonzero.
The matrices C 0 and C 1 exhibit behavior similar to the matrices D 0 and D 1 except that C v v 0, 1 is a new matrix representing the arrival of v supercustomers where a supercustomer is the last customer to form a group of k customers or the arrival of v regular customers to an arbitrary server.On the other hand, D v v 0, 1 represents arrival of v regular customers to the GI/D/k system.Now, this new arrival process k GI to an arbitrary server can be represented by a discrete time Markov chain { N t , J t ; t 0, 1, 2, . . .; N t ∈ {0, 1, 2, . ..};J t ∈ {1, 2, . . ., nk}} where N t is the number of supercustomers that have arrived by time t and J t is the phase of the next arrival at time t.

The New Model k GI/D/1 System
A model can be developed for a single server queue with arrival k GI and deterministic service time of d units following the techniques of Alfa 3 .This single server queue is called k GI/D/1 system and the resulting model is of QBD process type for which standard methods can be used 22 .Let us assume that L t , J t , and S t represent the number of supercustomers, the arrival phase of the next supercustomer, and the elapsed time of the supercustomer in service, respectively, at time t in this QBD model.The state space can be represented by {L t , J t , S t } for L t ≥ 1 and the state space is {0, J t } for L t 0. The total state space is { 0, j ∪ i, j, l : i ≥ 1, 1 ≤ j ≤ nk, 1 ≤ l ≤ d}.The transition matrix for this k GI/D/1 system is where . .be the stationary distribution of P where x xP, x 0 1 nk ∞ i 1 x i 1 nkd 1 and x j j ≥ 0 is the probability vector representing j supercustomers in the k GI/D/1 system.During the computation of the stationary distribution of queue length x, a matrix R is computed first which is followed by the computations of x 0 and x 1 .Then x i i ≥ 2 is computed using the matrix-geometric result represents the expected number of visits into m 1, v , starting from m, u , before the first return to level m m ≥ 1 .Here, m is the number of supercustomers in the system and u and v are any of the nkd phases of the system.

Structures of the Matrices
This section explains the structures of the block-matrices B 0 , B 1 , E, A 2 , A 1 , and A 0 and how their structures are exploited for efficient computation of the matrix R. The structures of these matrices are also fully exploited in computing the stationary distribution of the queue length in Section 4.
Here, B 0 and B 1 are matrices of order nk × nk and nk × nkd, respectively.The first column and the last row of the matrix B 0 are zero.On the other hand, only the last n rows of the matrix B 1 are nonzero and among these last n rows only the first column is nonzero which is t.The matrix B 1 can be represented as B 1 The matrix E is of order nkd × nk which can be represented as where Moreover, L has n − 1 nonzero columns and each column contains only one nonzero element.L has n−1 nonzero elements in jd, j 1 th 1 ≤ j ≤ n − 1 positions.M D 1 ⊗ s is also a matrix of order nd × n matrix.Furthermore, only the first column of M is nonzero which contains n nonzero elements.M has nonzero elements in jd, 1 th 1 ≤ j ≤ n positions.
The matrices A 2 , A 1 , and A 0 are of size nkd × nkd.The matrix A 2 can be represented as where L 1 and M 1 are both matrices of order nd × nd, There are n − 1 nonzero columns in L 1 and each column contains only one nonzero element.L 1 has these nonzero elements in jd, jd 1 th 1 ≤ j ≤ n − 1 positions.Only the first column of M 1 is nonzero with n nonzero elements in jd, 1 th 1 ≤ j ≤ n positions.
The matrix A 1 can be represented as where . . .
Here, L 2 and M 2 are both matrices of order nd × nd.

Computation of the Matrix R
The matrix R can be computed using any of the three linearly convergent formulae from Neuts 22 -R i 1 : which are known as natural, traditional, and U-based iterations, respectively.Here, R i is the ith iteration value of R, R 0 : 0 nkd,nkd , U is a square matrix of the same order as the matrix R and represents the taboo probability of starting from m, u for m ≥ 1 the eventual visit of the Markov chain to level m by visiting m, v under the taboo of level m − 1 .Any of these three methods terminate in the ith i ≥ 1 iteration for R i − R i − 1 ∞ < ε where ε is a very small positive constant.In each step, natural itearation requires two matrix multiplications whereas traditional iteration requires three matrix multiplications and one matrix inversion.On the other hand, U-based iteration requires two matrix multiplications and one matrix inversion in each step.Traditional iteration requires less number of iterations than natural iteration does while the number of iterations required for U-based iteration is less than the numbers of iterations for other two methods.Alfa and Xue 25 developed efficient inversion technique for U-based iteration for the matrix R in analyzing the queue length of a GI/G/1 system in discrete time.
There are several efficient methods for computing the matrix R for a QBD such as logarithmic reduction algorithm of Latouche and Ramaswami 26 , cyclic reduction technique of Bini and Meini 27 , and invariant subspace approach of Akar and Sohraby 28 .These quadratically convergent algorithms compute a matrix G which is the minimal nonnegative solution to G Here, G is an m × m matrix for m × m matrices A 2 , A 1 , and A 0 where G i,j 1 ≤ i, j ≤ m represents the probability that starting from state N 1, i the Markov chain eventually visits the level N N ≥ 1 and does so by visiting the state N, j .The matrices R, U, and G are related and these relations can be found in 26 .Similar to the computation of the matrix R, the matrix G can be computed using three different iterations-natural, traditional, and U-based.Generally a zero matrix is used for the initial value for iterations of the matrix G i.e., G 0 m,m .Latouche and Ramaswami 26 showed that if S G ε iterative steps are needed to converge for a very small positive constant , then logarithmic reduction requires not more than log 2 S G ε iterations and overall complexity of their algorithm is O m 3 log 2 S G ε .They also developed logarithmic reduction technique for the matrix R.However, the quadratic convergent algorithms 26-28 are not considered here for computing the matrix R as these algorithms involve inversion of matrices which cannot exploit the special structures of the matrices A 2 , A 1 , and A 0 .
The complexity of computing the matrix R can be reduced by using the structure of A 0 .The first nd k −1 rows are zero in A 0 .Therefore, the first nd k −1 rows are zero in R.Among the last nd rows of A 0 , n rows are zero which are rows nd k − 1 d, nd k − 1 2d, nd k − 1 3d, . . ., nkd.The corresponding rows are zero for R. The matrix R can be represented as follows: Here, R j 1 ≤ j ≤ k is a matrix of order nd × nd which contains d, 2d, 3d, . . ., ndth rows as zero rows.
The computations of R i A 2 and R i A 1 require 2k − 1 and 2k block-matrix multiplications, respectively, where each block matrix is of order nd × nd.Similarly, the computation of R 2 i A 2 requires multiplication of R i and R i A 2 which involves k blockmatrix multiplications.Inversions of both the matrices I nkd − A 1 and I nkd − A 1 − R i A 2 generate k 2 block matrices each of order nd × nd.Inversion of the matrix I nkd − A 1 needs to be computed only once for traditional iteration whereas inversion of the matrix I nkd − A 1 − R i A 2 needs to be computed in each iteration for U-based iteration.On the other hand, the multiplication of A 0 R 2 i A 2 and I nkd − A 1 −1 in traditional iteration and the multiplication of A 0 and I nkd − A 1 − R i A 2 −1 in U-based iteration involve k 2 and k blockmatrix multiplications, respectively, where each block matrix is of order nd×nd.Each iteration of natural method requires substantially less multiplications than each step of traditional and U-based methods do.Therefore, we use natural iteration for developing our algorithm for computing the matrix R in this paper.
The matrix R j 1 ≤ j ≤ k of 3.7 can be computed using the natural iteration with the following equation: where R j i i ≥ 0, 1 ≤ j ≤ k is the ith iteration value of R j with the initial value R j 0 : 0 nd,nd 1 ≤ j ≤ k .The terminating condition for this iteration is same as the terminating condition for natural, traditional, and U-based methods but needs to consider only The number of multiplications required to compute Similarly, the number of multiplications required to compute both R j i L 1 and In the same way, the number of multiplications required to compute Hence, the numbers of multiplications required for computing R 1 i 1 from R 1 i and respectively.Therefore, the total number of multiplications required to compute all k block matrices R j i 1 's 1 ≤ j ≤ k of the matrix R i 1 in one iteration is By exploiting the special structures of the matrices the computational complexity is decreased in each iteration from O n 3 k 3 d 3 to O n 3 kd 2 which is substantial reduction in computations.Moreover, the memory requirement for our method is also less than that of natural, traditional, U-based methods and quadratic convergent algorithms.

Stationary Distribution of P for k GI/D/1 System
Two methods for computing stationary distribution x for the Markov chain represented by the matrix P are briefly presented here which are followed by Section 4.1.Section 4.1 explains how the second method can be made efficient by utilizing the special structures of the matrices.
In the first method an invariant probability vector z z 0 , z 1 needs to be computed for a stochastic matrix The row vector z is computed from z P z and z 0 1 nk z 1 1 nkd 1.Then x 0 and x 1 can be expressed as x 0 v −1 z 0 and In the second method, we need to compute h, the invariant probability vector of a stochastic matrix H defined as On the other hand, h is computed using the relations h hH and h1 nkd 1. x 0 and x 1 can be computed as

Exploiting Structures for Stationary Distribution
The stationary distribution of x can be computed by first computing the probability vector x 1 which can be made efficient by exploiting the special structures of the matrices H, E I nk − B 0 −1 1 nk and I nkd −R −1 1 nkd .The special structures of the matrices , and H are due to the special structures of the matrices B 0 , B 1 , E, and R. Here, how the computations of these matrices can be carried out efficiently is explained. Computing requires O n 3 kd 3 multiplications by using block-matrix inversion 29 where Υ −1 can be represented as Here, V i and The matrix H of can be represented in the following form: where Here,

4.4
We can recursively define Let us define another invariant probability vector u u 1 , u 2 , u 3 , . . ., u k of the stochastic matrix H where u i 1 ≤ i ≤ k is a row vector having nd scalar elements.Here, u is computed from x 1 can be computed as x 1 u −1 u using O n 3 kd 3 multiplications instead of O n 3 k 3 d 3 multiplications using the special structure of H. Let us partition x 0 and x 1 as x 0 κ 1 , κ 2 , κ 3 , . . ., κ k and x 1 ν 1 , ν 2 , ν 3 , . . ., ν k , respectively, for efficient computation of x 0 .Here, κ i and ν i for 1 ≤ i ≤ k are of size 1 × n and 1 × nd, respectively.κ i 1 ≤ i ≤ k is computed as κ i i j 1 ν i−j 1 V j and these vector-matrix multiplications can be made efficient by considering only nonzero rows of V j 1 ≤ j ≤ k .
x i i ≥ 2 can be computed as

Waiting-Time Distribution of k GI/D/1 System
The waiting time distribution of an arbitrary customer in k GI/D/1 system is the same as the waiting time distribution of a customer in GI/D/k system.Let w r be the probability that a customer i.e., the kth customer of a supercustomer group has to wait r units of time.w r is computed after obtaining the probability vectors y i 's where y i i ≥ 0 is the stationary probability vector that a supercustomer sees i supercustomers ahead of him.y 0 can be expressed as y 0 y 0,1 , y 0,2 , . . ., y 0,nk .Similarly, y i i ≥ 1 can be expressed as y i y i,1 , y i,2 , . . ., y i,nk and y i,v can be expressed as Two approaches for computing the distribution of waiting time are presented in Sections 5.1 and 5.2, respectively.The first method uses the approach for phase type service time distribution and this approach is used in the analysis of Alfa 3 and Chaudhry et al. 7 .
In the second method we develop a computationally simpler approach.This is a simplified version of Method 1.It simplifies the computations of the probability vectors y i 's for i ≥ 0 and w r , the probability mass function of the distribution of waiting time.

Method 1
The expressions for the stationary distribution of an entering supercustomer to find different number of supercustomers ahead of him in the system can be written as

5.1
In this method, the probability that the waiting time is r units of time is computed as where Ω i r is a square matrix of order d that represents the probability distribution of r unit of work in system at the arrival of a supercustomer who finds that i supercustomers are ahead of him.For example, the element Ω i r u,v represents the probability that the service of the arbitrary supercustomer, who arrives to find i supercustomers in the system with remaining workload of r units begins in phase v 1 ≤ v ≤ d , given that service of the supercustomer who was in service at the arrival of the arbitrary supercustomer was in phase u 1 ≤ u ≤ d .The matrix Ω i r is computed in the following way: 5.3

Method 2
This method simplifies the computation of y 0 and y i i ≥ 1 .Moreover, this method does not need to compute the matrix Ω i r i, r ≥ 0 which simplifies the computation of the probability mass function of the distribution of waiting time.
The special structures of C 1 , s, S, and β can be exploited to reduce the number of operations to compute y 0 and y i i ≥ 1 .Only the first columns of the matrices C 1 and C 1 ⊗ s are nonzero.Therefore, the first element of y 0 is nonzero and the remaining nk − 1 elements of y 0 are equal to zero.This phenomenon is quite intuitive as the arriving supercustomer does not find anybody ahead of him, he can start his service in phase 1 instantly and the arrival process for next supercustomer starts at phase 1 of the nk phases.The first element y 0,1 of y 0 is computed as x 0, k−1 n j t j x 1, k−1 n j,d t j .

5.4
d − 1 columns of C 1 ⊗ S are nonzero which are column 2 to column d.On the other hand, only the first column of C 1 ⊗ sβ is nonzero.Therefore, the first d elements of y i i ≥ 1 are nonzero.It is quite intuitive that the first d elements of y i i ≥ 1 are nonzero i.e., y i,1 is nonzero and y i,v 0 1,d 2 ≤ v ≤ nk as the arriving supercustomer finds i supercustomers ahead of him, the arrival process of next supercustomer starts at phase 1 of nk phases and the ongoing service is in the start of phase j 1 ≤ j ≤ d .The d nonzero elements y i,1,j 's i ≥ 1, 1 ≤ j ≤ d of the vector y i can be computed using

5.5
The special structures of 1 nk ⊗ I d and y i for i ≥ 1 result in

5.6
If a supercustomer sees i supercustomers ahead of him then he has to wait or id 1, id 2, . . ., id d amount of time for d > 1 or i units of time for d 1 before getting service.
First we consider the case of d 1.In this case, S 0 , s 1 , and β 1 .For d 1, 5.3 can be written as

5.7
For i < r we can assume i i 1 r and i 1 > 0. In this case, Ω i r For i r, we have Ω i r 1 from 5.7 .For i > r we can assume i r i 1 and i 1 > 0. In this case, Ω i r 0. Therefore, we can express for d 1, Now, we consider the case d ≥ 2. Let us define e j 1 ≤ j ≤ d as the jth column of an identity matrix of order d.It is found that for i ≥ 1 and d ≥ 2

5.10
Alfa 2, 3 pointed out that Ω i r u,v 0 for i ≥ 1, 2 ≤ v ≤ d.However, he did not provide any explicit expression for Ω i r 1 d and his method needs to compute the first column of the matrix Ω i r i ≥ 1 recursively as our method does not need to compute.Use of 5.10 to compute Ω i r 1 d can significantly reduce the complexity of computing  5.10 reduces the usage of memory as there is no need to compute and store the matrix Ω i r for different values of i and r.Now, if an arriving customer has to wait r units of time r ≥ 0 then there are r/d supercustomers ahead of him and the customer, who is in service in the server, is in phase d − r mod d of service time.Therefore, the probability of a customer waiting r r ≥ 0 units of time can be expressed by using value of Ω i r 1 d of 5.8 and 5.10 in 5.2 5.11  1 and 2 show the number of iterations and time required to compute the matrices R, G, and U for the interarrival times having probability vector 0.1, 0.2, 0.2, 0.5 .The means of interarrival time and arrival rate for this probability vector are 3.1 time slots and 0.3226/time slot, respectively, and ρ 0.8602.Tables 3 and 4 show the number of iterations and time required for the interarrival times having probability vector 0.1, 0.2, 0.3, 0.4 .The means of interarrival time and arrival rate for the probability vector 0.1, 0.2, 0.3, 0.4 are 3.0 time slots and 0.3333/time slot, respectively, and ρ 0.8889.Our methods of computing the matrix R are found to require least time for these two examples among the methods used here though logarithmic reduction technique and U-based methods require less number of iterations.Tables 5 and 6 show the number of iterations and time required for the interarrival times having probability vector 0.1, 0.3, 0.4, 0.2 .The means of interarrival time and arrival rate for this probability vector 0.1, 0.3, 0.4, 0.2 of interarrival times are 2.7 time slots and 0.3704/time slot, respectively, and ρ 0.9877.Table 6 shows that for the interarrival times probability vector 0.1, 0.3, 0.4, 0.2 , Logarithmic Reduction techniques for the matrices G and R require less time to converge than our method of exploiting the structure of the matrices A 2 , A 1 , and A 0 .The relative advantage of our method of exploiting the structures to Logarithmic Reduction techniques decreases with respect to time with increasing value of ρ.

Conclusion
In this paper, a class of GI/D/k systems in discrete time is analyzed by converting it to a single server queue problem with a convoluted arrival process.Then the stationary distribution of the length of the queue of the transformed system is analyzed using QBD approach.The special structures of the matrices make the computation of the matrix R and the distribution of queue length of the transformed system efficient.Numerical examples show that our proposed method for computing the matrix R can be as efficient as quadratic convergent algorithms such as logarithmic reduction technique.It is also shown that for deterministic service time the waiting time distribution does not need any cumbersome computation like phase type service time distribution.
23 , GI X /D/c queue by Chaudhry and Kim 5 , M X /D/c queue by Chaudhry et al. 6 and Franx 9 , Ph/D/c queues by van Hoorn 16 , and an E k /D/r queue by Xerocostas and Demertzes 20 .

L 2 D 2 .
0 ⊗ S and M 2 D 1 ⊗ S.There are n − 1 d − 1 nonzero columns, n − 1 d − 1 nonzero rows, and n − 1 d − 1 nonzero elements in L 2 .Moreover, each nonzero row or each nonzero column of L 2 contains only one nonzero element.On the other hand, there are d − 1 nonzero columns in M 2 from column 2 to d each having n nonzero elements.The matrix A 0 can be represented as A 0 M A 0 has only n d − 1 nonzero rows which are in its last nd rows and d − 1 nonzero columns from column 2 to column d.
and D 0 I n − D 0 .The computation of I nk − B 0 −1 requires O n 2 k nk 2 instead of O n 3 k 3 multiplications by utilizing the special structure of I nk − B 0 .The matrices E I nk − B 0 −1 and E I nk − B 0 −1 B 1 can be represented as are matrices of order nd × n and nd × nd, respectively.Here, multiplying E by I nk − B 0 −1 and E I nk − B 0 −1 by B 1 both require O n 2 k multiplications due to utilizing special structures of matrices.If there were no special structures of E, I nk − B 0 −1 and B 1 then the multiplying E by I nk − B 0 −1 and E I nk − B 0 −1 by B 1 would require O n 3 k 3 d and O n 3 k 2 d 3 multiplications, respectively.
positive real number u is computed as u u η ζ where η E I nk − B 0 −1 1 nk and ζ I nkd − R −1 1 nkd .η and ζ are both of size nkd×1.The column vectors η and ζ can be expressed as η d for j 0, 1, 2, . . ., i − 1 d,Ω i i − 1 d j S j−1 sβ e d−j 1 0 d,d−1 for j 1,2, . . ., d, Ω i id j 0 d,d for j ≥ 1, 5.9 which are proved by induction in 30 .If the matrix Ω i r is nonzero then it contains nonzero elements in the first column.The nonzero first column in nonzero Ω i r is due to elapsed time representation of service time where service starts at phase 1.Now, e d−j 1 0 d,d−1 1 d e d−j 1 .Therefore, we can express for d ≥ 2,

Table 1 :
Number of iterations required for different methods for interarrival times having probability The program was executed in SunOS 5.10 31 on an i386 processor which operates at 2660 MHz and has an i387 compatible floating point processor.The number of iterations and time in seconds is recorded in Tables1-6for each method from the execution of the program."LR" is used to denote logarithmic reduction technique in those tables.Tables This section compares our proposed method of computing the matrix R with natural, traditional, and U-based methods and Logarithmic Reduction technique 26 for the matrices R and G. Three different matrices R, G, and U are computed in each of these methods.The values of the convergence constant ε used are 10 −6 , 10 −7 , 10 −8 , and 10 −9 .Three identical servers i.e., k 3 are considered with duration of service as eight time units i.e., d 8 and the support for interarrival times is assumed to be four i.e., n 4 .Three

Table 2 :
Time in seconds required for different methods to execute for interarrival times having probability

Table 3 :
Number of iterations required for different methods for interarrival times having probability vector a 0.1, 0.2, 0.3, 0.4 .

Table 4 :
Time in seconds required for different methods to execute for interarrival times having probability vector a 0.1, 0.2, 0.3, 0.4 .

Table 5 :
Number of iterations required for different methods for interarrival times having probability vector a 0.1, 0.3, 0.4, 0.2 .

Table 6 :
Time in seconds required for different methods to execute for interarrival times having probability vector a 0.1, 0.3, 0.4, 0.2 .