BeginPackage["queue4`"]

queue4::usage="This is a collection of queueing theory programs used 
to model computer systems." 
Exact::usage="Exact[ Nc_?VectorQ, Tc_?VectorQ, Dck_?MatrixQ ] 
computes the performance statistics for a closed, multiclass BCMP 
queueing system where Nc is the population vector, Tc is the 
vector of think times and Dck is the matrix of service demands 
by workload class.  The output consists of four performance vectors: 
Wc the vector of response times by class, Lambdac the vector of 
throughputs by class, Lk the vector of number of customers at each 
service center, and Rhok the vector of the total utilization of each 
service center." 
Approx::usage=" Approx[ Nc_?VectorQ, Tc_?VectorQ, Dck_?MatrixQ, 
epsilon_Real ] computes the approximate performance statistics for 
a closed, multiclass BCMP queueing system where Nc is the population 
vector, Tc is the vector of think times and Dck is the matrix of 
service demands by workload class.  The output consists of four 
performance vectors: Wc the vector of response times by class, Lambdac
the vector of throughputs by class, Lk the vector of number of 
customers at each service center, and Rhok the vector of the total 
utilization of each service center."

Begin["queue4`private`"]

FixPerm[ numC_, PVc_, Nc_ ] := 
Block[ {i, m = PVc }, 
       For[i=numC, i>1, i--, 
         If[m[[i]] > Nc[[i]], 
            m[[i-1]]=m[[i-1]]+m[[i]]-Nc[[i]] ; 
            m[[i]]=Nc[[i]] ]];    
       If[ m[[1]] > Nc[[1]], {}, m] 
]

FirstPerm[ numC_, Nc_, n_ ] := 
Block[ {m}, 
    m = Table[ 0, { numC } ] ; 
    m[[numC]] = n ;  
    FixPerm[numC, m, Nc ]  
]

NextPerm[ numC_, Nc_, PVc_ ] := 
Block[ {m=PVc, i=numC, j}, 

      While[ m[[i]] == 0, i-- ] ; 
      If[i==1, Return[{}] ] ; 

      m[[i]]-- ; 
      i-- ; 
      While[ (i >= 1) && (m[[i]] == Nc[[i]]), i-- ] ; 
      If[i < 1, Return[{}] ] ; 
      m[[i]]++ ; 

      For[j=i+1, j<numC, j++, 
          m[[numC]] = m[[numC]] + m[[j]] ;
          m[[j]] = 0 ] ;
      FixPerm[numC, m, Nc ]  
]

Exact[ Nc_?VectorQ, Tc_?VectorQ, Dck_?MatrixQ ] := 

Block[ { n, PVc, PVcm1, Wck, Wc, Lambdac, Lck1, Lck2, Ltmp, Rhok, Lk,     
         numC = Length[Nc], numK = Dimensions[Dck][[2]],
         Ntot, Wtot, zVectorK }, 

   zVectorK = N[ Table[0, {numK}] ] ;
   Ntot = Sum[ Nc[[c]], {c, 1, numC} ] ; 
   Lck1[ Table[0, {numC}] ] = zVectorK ; 

   For[n=1, n <= Ntot, n++,  
       PVc = FirstPerm[numC, Nc, n ] ;
       While[PVc!= {},  
             Wck = Table[(PVcm1 = PVc ;
                          If[ PVcm1[[c]] > 0, 
                              PVcm1[[c]]-- ;
                              Dck[[c]] * ( 1 + 
                              If[OddQ[n], Lck1[PVcm1], Lck2[PVcm1]]), 
                              zVectorK]), 
                   { c, 1, numC} ] ;   

             Lambdac = Table[ ( Wtot = Tc[[c]] + Apply[Plus, Wck[[c]] ] ;
                           If[ Wtot > 0, PVc[[c]] / Wtot, 0 ] ), 
                  {c, 1, numC} ] ;

             Ltmp = Lambdac . Wck ; 
             If[OddQ[n], Lck2[PVc]=Ltmp, Lck1[PVc]=Ltmp ] ;
  
             PVc = NextPerm[numC, Nc, PVc ] ]; 
       If[OddQ[n], Clear[Lck1], Clear[Lck2]] 
   ] ;
    
   Wc = Apply[Plus, Wck, 1 ] ;
   Rhok = Lambdac . Dck ;
   Lk = Lambdac . Wck ;

   Return[{Wc,Lambdac,Lk,Rhok }] ;

] /; Length[Nc] == Length[Tc] == Length[Dck]



Approx[ Nc_?VectorQ, Tc_?VectorQ, Dck_?MatrixQ, epsilon_Real ] := 

Block[ { Flag, Wck, Lambdac, newL, Lk, Lck, Wc, Lk, Rhok,  
         numC = Length[Nc], numK = Dimensions[Dck][[2]] }, 
 
   Lck = N[Table[ Nc[[c]]/numK, {c, 1, numC}, {k, 1, numK} ] ] ;
   Flag = True ; 
   While[Flag==True,

       Lk = Apply[Plus, Lck ] ; 
       Wck = Table[ Dck[[c,k]] * 
             (1 + Lk[[k]] - Lck[[c,k]] + ((Nc[[c]]-1)/Nc[[c]]) Lck[[c,k]] ), 
             {c, 1, numC}, {k, 1, numK} ] ;
       Wc = Apply[Plus, Wck, 1 ] ; 
       Lambdac = Nc / ( Tc + Wc) ; 

       Flag = False ; 
       Lck = Table[(newL = Lambdac[[c]] Wck[[c,k]] ; 
                    If[ Abs[ Lck[[c,k]] - newL ] >= epsilon, Flag=True] ;
                    newL), {c, 1, numC}, {k, 1, numK} ] ;
   ] ;

(* Compute final results *)

   Rhok = Lambdac . Dck ;
   Lk = Lambdac . Wck ;
   Return[{Wc,Lambdac,Lk,Rhok }] ;

] /; Length[Nc] == Length[Tc] == Length[Dck]

End[] (* end 'private' context *)
EndPackage[]
